Fermions on the square lattice & perfect nesting

This tutorial is the first in a series of tutorials on two-particle response were we will use TRIQS and the Two-Particle Response-Function toolbox (TPRF) to compute;

1. the non-interacting Green’s function of fermions on the square lattice with nearest-neighbour hopping and study the Fermi surface [01]

2. the non-interacting two-particle response, also called the susceptibility [03]

3. the Random-Phase Approximation (RPA) susceptibility for weak interactions, studying the anti-ferromagnetic divergence at (\(\pi,\pi)\) [05]

4. the Two-Particle Self Conistent (TPSC) susceptibility and show that it satisfies the Pauli principle, while RPA does not [07]

5. and show that TPSC obeys the Mermin-Wagner theorem, since its spin susceptibility does not diverge at finite temperature [09]

%matplotlib inline
%config InlineBackend.figure_format = 'svg'
from triqs.plot.mpl_interface import plt

import numpy as np

Square lattice with nearest-neighbour hopping

The square lattice with nearest-neighbour hopping \(t\) appeared in earier tutorials, were the dispersion relation

:nbsphinx-math:`begin{equation}

epsilon(mathbf{k})=-2t(cos{k_x}+cos{k_y}),

end{equation}`

was computed using TRIQS in more than one way.

However, in TRIQS and TPRF there are a number of helper routines for lattice models that simplifies the study of general tight-binding models. Here we will therefore use these standard routines.

Insead of constructing \(\epsilon(\mathbf{k})\) directly in momentum space we construct a real-space tight binding lattice Hamitonian \(H(\mathbf{r})\) using the ``TBLattice` class <https://triqs.github.io/triqs/latest/documentation/python_api/triqs.lattice.tight_binding.TBLattice.html?highlight=tblattice#triqs.lattice.tight_binding.TBLattice>`__ corresponding to the square lattice with nearest neighbour hopping \(t=1\).

from triqs.lattice.tight_binding import TBLattice

t = 1.0 # nearest-neigbhour hopping amplitude

H_r = TBLattice(
    units=[
        (1,0,0), # basis vector in the x-direction
        (0,1,0), # basis vector in the y-direction
    ],
    hoppings={
        (+1,0) : [[-t]], # hopping in the +x direction
        (-1,0) : [[-t]], # hopping in the -x direction
        (0,+1) : [[-t]], # hopping in the +y direction
        (0,-1) : [[-t]], # hopping in the -y direction
    })

print(H_r)

Warning: could not identify MPI environment!
Tight Binding Hamiltonian on Bravais Lattice with dimension 2, units
[[1,0,0]
 [0,1,0]
 [0,0,1]], n_orbitals 1
with hoppings [
   [1,0] :
[[(-1,0)]]
   [-1,0] :
[[(-1,0)]]
   [0,1] :
[[(-1,0)]]
   [0,-1] :
[[(-1,0)]] ]

Starting serial run at: 2023-08-29 11:09:00.866052

The real-space Hamiltonian \(H(\mathbf{r})\) can both construct a discretized momentum mesh and evaluate the dispersion \(\epsilon(\mathbf{k})\) on a given momentum mesh, by

n_k = 128
kmesh = H_r.get_kmesh(n_k=n_k)
e_k = H_r.fourier(kmesh)
print(e_k)

Greens Function  with mesh Brillouin Zone Mesh with linear dimensions (128 128 1)
 -- units =
[[0.0490874,0,0]
 [0,0.0490874,0]
 [0,0,6.28319]]
 -- brillouin_zone: Brillouin Zone with 2 dimensions and reciprocal matrix
[[6.28319,0,0]
 [0,6.28319,0]
 [0,0,6.28319]] and target_shape (1, 1):

Since the square lattice is two-dimensional it is possible to visualize the dispersion using a color plot.

Here is an example that plots \(\epsilon(\mathbf{k})\) in the entire Brillouin zone as well as the shape of the Fermi surface at \(\omega = 0\) (dotted line).

k = np.linspace(-np.pi, np.pi, num=100)
kx, ky = np.meshgrid(k, k)

e_k_interp = np.vectorize(lambda kx, ky : e_k((kx, ky, 0)).real)(kx, ky)

plt.figure()

plt.pcolormesh(kx, ky, e_k_interp, rasterized=True, cmap='RdBu')
plt.colorbar().ax.set_ylabel(r'$\epsilon(\mathbf{k})$');

plt.contour(kx, ky, e_k_interp, levels=[0], linestyles='dotted')

plt.xlabel(r'$k_x$'); plt.ylabel(r'$k_y$');
k_ticks, k_labels = [-np.pi, 0, np.pi], [r"$-\pi$", r"0", r"$\pi$"]
plt.xticks(k_ticks, k_labels); plt.yticks(k_ticks, k_labels);
plt.axis('square');

# -- High-symmetry path G-X-M-G

pts = [
    (0., 0., r'$\Gamma$', 'w', 'bottom', 'right'),
    (np.pi, 0., r'$X$', 'k', 'center', 'left'),
    (np.pi, np.pi, r'$M$', 'r', 'bottom', 'left'),
    ]
for x, y, label, color, va, ha in pts:
    plt.plot(x, y, 'o', color=color, clip_on=False, zorder=110)
    plt.text(x, y, label, color=color, fontsize=18, va=va, ha=ha)

X, Y, _, _, _, _ = zip(*(pts+[pts[0]]))
plt.plot(X, Y, '-m', zorder=100, lw=4);

Momentum dependent quantities can also be visualized along high-symmetry paths in the Brillouin zone, see above for the high-symmetry points \(\Gamma\), \(X\) and \(M\) of the square lattice.

Here is an example that plots the dispersion \(\epsilon(\mathbf{k})\) along thepath \(\Gamma - X - M - \Gamma\) in k-space using the ``triqs_tprf.lattice_utils.k_space_path` function <https://triqs.github.io/tprf/unstable/reference/python_reference.html#triqs_tprf.lattice_utils.k_space_path>`__.

G = [0.0, 0.0, 0.0]
X = [0.5, 0.0, 0.0]
M = [0.5, 0.5, 0.0]

path = [(G, X), (X, M), (M, G)]

from triqs_tprf.lattice_utils import k_space_path

k_vecs, k_plot, k_ticks = k_space_path(path, num=32, bz=H_r.bz)

e_k_interp = np.vectorize(lambda k : e_k(k).real, signature='(n)->()')

plt.plot(k_plot, e_k_interp(k_vecs))
plt.xticks(k_ticks, labels=[r'$\Gamma$', '$X$', '$M$', r'$\Gamma$'])
plt.ylabel(r'$\epsilon(\mathbf{k})$')
plt.grid(True)

In the following we will re-purpose these visualization scripts to study the one-particle and two-particle Green’s functions of the square lattice model.

Non-interacting lattice Green’s function

Given the dispersion \(\epsilon(\mathbf{k})\) the non-interacting Green’s function \(G_0(i\omega_n, \mathbf{k})\) is given by

:nbsphinx-math:`begin{equation}

G_0(iomega_n, mathbf{k}) = frac{1}{iomega_n - epsilon(mathbf{k})} , .

end{equation}`

As shown in the Basic Tutorial it is of course possible to compute \(G_0\) using a loop over frequency and momentum:

from triqs.gf import Gf, MeshImFreq, MeshProduct

wmesh = MeshImFreq(beta=2.5, S='Fermion', n_max=128)
wkmesh = MeshProduct(wmesh, kmesh)
g0_wk = Gf(mesh=wkmesh, target_shape=e_k.target_shape)

for w, k in wkmesh:
    g0_wk[w, k] = 1/(w - e_k[k])

However, TPRF has Dyson equation solvers that are OpenMP+MPI parallell and all implemented in C++, see the TPRF documentation. Here we will use these fast routines!

Exercise 1:

---

Use `triqs_tprf.lattice.lattice_dyson_g0_wk <https://triqs.github.io/tprf/latest/reference/python_reference.html#triqs_tprf.lattice.lattice_dyson_g0_wk>`__ to compute \(G_0(i\omega_n, \mathbf{k})\) at inverse temperature \(\beta = 2.5\) using a fermionic ``MeshImFreq` frequency mesh <https://triqs.github.io/triqs/latest/documentation/python_api/triqs.gf.meshes.MeshImFreq.html?highlight=meshimfreq#triqs.gf.meshes.MeshImFreq>`__ with 128 Matsubara frequencies and name the resulting Green’s function g0_wk.

Check the properties of g0_wk by printing it, i.e.

print(g0_wk)

# Write your code here

from triqs.gf import MeshImFreq
wmesh = MeshImFreq(beta=2.5, S='Fermion', n_iw=32)

from triqs_tprf.lattice import lattice_dyson_g0_wk
g0_wk = lattice_dyson_g0_wk(mu=0., e_k=e_k, mesh=wmesh)

print(g0_wk)

Greens Function  with mesh Imaginary Freq Mesh with beta = 2.5, statistic = Fermion, n_iw = 32, positive_only = false, Brillouin Zone Mesh with linear dimensions (128 128 1)
 -- units =
[[0.0490874,0,0]
 [0,0.0490874,0]
 [0,0,6.28319]]
 -- brillouin_zone: Brillouin Zone with 2 dimensions and reciprocal matrix
[[6.28319,0,0]
 [0,6.28319,0]
 [0,0,6.28319]] and target_shape (1, 1):

Questions

• How many meshes does the Green’s function have?

• What is the k-space discretization?

• How is the reciprocal basis vectors of the Brillouin zone related to the lattice (units) vectors of the tight binding lattice H_r above?

Exercise 2:

---

Save the Green’s function g0_wk into a hdf5 file named g0_wk.h5 using `h5.HDFArchive <https://triqs.github.io/triqs/latest/documentation/manual/triqs/hdf5/ref.html>`__, so that we can read g0_wk from file in the following tutorials. Do the same with the dispersion e_k and save it to a file with the name e_k.h5.

# Write your code here

from h5 import HDFArchive

with HDFArchive("g0_wk.h5", "w") as R:
    R['g0_wk'] = g0_wk

with HDFArchive("e_k.h5", "w") as R:
    R['e_k'] = e_k

Fermi surface nesting

We will now study the Fermi surface of Fermions on the square lattice with nearest neighbour hopping, which has a special property called perfect nesting. A Fermi surface is said to be nested if parts of the Fermi surface map to each other by a single momentum vector \(\mathbf{Q}\), called the nesting vector.

For a non-interacting system the Fermi surface is the surface in k-space defined by

\[\epsilon(\mathbf{k}) - \mu = 0 \, ,\]

where \(\mu\) is the chemical potential. In terms of the spectral function \(A(\omega, \mathbf{k})\) this corresponds to large values of \(A\) at \(\omega=0\)

\[ A(\omega = 0, \mathbf{k}) = \frac{1}{\pi} \text{Im} \left[ \frac{1}{ 0 - \epsilon(\mathbf{k}) + \mu - i\delta } \right] \gg 1 \, ,\]

which also generalizes to interacting systems.

Exercise 3:

---

Make a color plot of the zero-frequency spectral function \(A(k, \omega=0)\) over the Brillouin zone, using the approximation

\[A(k, \omega=0) \approx -\frac{1}{\pi} \text{Im}[ G_0(\mathbf{k}, i\omega_0) ] \, ,\]

where we neglect the fact that the first fermionic Matsubara frequency \(i\omega_0\) is not exactly \(0\).

The right hand side can be evaluated using the Triqs Green’s function g0_wk and the interpolation feature:

n = 0 # Matsubara frequency index
kx, ky, kz = 0., 0., 0.
k_vec = (kx, ky, kz)
g0_wk(n, k_vec)

# Write your code here

kgrid1d = np.linspace(-np.pi, np.pi, n_k + 1, endpoint=True)
kx, ky = np.meshgrid(kgrid1d, kgrid1d)

A_k = np.vectorize(lambda kx, ky: -g0_wk( 0, (kx, ky, 0) ).imag / np.pi)
A_inv_k = np.vectorize(lambda kx, ky: (1 / g0_wk( 0, (kx, ky, 0) )).real)

plt.contour(kx, ky, A_inv_k(kx, ky), levels=[0], colors='white')
plt.pcolormesh(kx, ky, A_k(kx, ky), rasterized=True)

plt.colorbar().ax.set_ylabel(r"$G_0(i\omega_0, \mathbf{k})$")
plt.xticks([-np.pi, 0, np.pi],[r"$-\pi$", r"0", r"$\pi$"])
plt.yticks([-np.pi, 0, np.pi],[r"$-\pi$", r"0", r"$\pi$"])
plt.axis('square'); plt.xlabel(r"$k_x$"); plt.ylabel(r"$k_y$");

Hint: Re-purpose the code for the color plot of \(\epsilon(\mathbf{k})\) above.

Questions

• How can we see from the plot that the Fermi surface is nested?

• What is the nesting vector?

• Actually the Fermi surface is perfectly nested. What do you think is the difference between nesting and perfect nesting?

Exercise 4:

---

For non-interacting systems the Fermi surface can also be observed in the Fermionic density distribution in k-space \(\rho(\mathbf{k})\), which can be computed from the Matsubara Green’s function \(G_0(i\omega_n, \mathbf{k})\).

In TRIQS this is available through the .density() method of Matsubara frequency Green’s functions. For the lattice Green’s function it can be evaluated for a given k-vector using

kx, ky, kz = 0., 0., 0.
k_vec = (kx, ky, kz)
rho_k = g0_wk(all, k_vec).density().real

Plot the momentum distribution \(\rho(\mathbf{k})\) along the Brillouin zone high-symmetry path \(\Gamma - X - M - \Gamma\).

rho_k_interp = np.vectorize(lambda k: g0_wk(all, k).density(), signature='(n)->()')

plt.plot(k_plot, rho_k_interp(k_vecs).real)
plt.xticks(k_ticks, labels=[r'$\Gamma$', '$X$', '$M$', r'$\Gamma$'])
plt.ylabel(r'$\rho(\mathbf{k})$')
plt.grid(True)

Hint: Re-purpose the code above plotting \(\epsilon(\mathbf{k})\) along the same path.

Questions

• What is the value of \(\rho(\mathbf{k})\) at the Fermi surface?

• What is the sign of \(\epsilon(\mathbf{k})\) in the regions of k-space where \(\rho(\mathbf{k}) \approx 1\) and \(\approx 0\), respectively?

• How is this related to the Fermi distribution function

  \[f(\epsilon) = \frac{1}{1 + e^{\beta \epsilon}}\]

  plotted below?

beta = 2.5
e = np.linspace(-4., 4.)
f = lambda e : 1/(1 + np.exp(beta * e))
plt.plot(e, f(e))
plt.xlabel(r'$\epsilon$'); plt.ylabel(r'$f(\epsilon)$'); plt.grid(True);