Square lattice susceptibility and the Random Phase Approximation (RPA)
A nice example of two-particle response functions is the two-dimensional square lattice with nearest neighbour hopping \(t\). Here we will use TRIQS and TPRF routines to compute the bare generalized susceptibility
and the generalized susceptibility in the random phase approximation (RPA)
In order to do this we first have to setup the tight binding Hamiltonian \(H\), the single particle dispersion \(\epsilon(\mathbf{k})\), and the single-particle Green’s function \(G_0(i\omega_n, \mathbf{k})\) of the system.
Tight binding Hamiltonian
The tight binding Hamiltonian of the square lattice with nearest neigbhour hopping \(t\) takes the form
where \(c^\dagger_{\sigma i}\) creates a fermion with spin \(\sigma \in \{\uparrow, \downarrow\}\) on lattice site \(i\), and \(\langle i, j \rangle\) denotes summation over nearest neighbours \(i\) and \(j\).
A representation of \(H\) can be constructed using the TBLattice class in triqs_tprf.tight_binding where the hoppings are given as a dictionary with relative coordinate vectors as keys and hopping matrices as values. The unit vectors of the lattice, the position of the site local orbitals and names also needs to be setup, see below.
[2]:
from triqs_tprf.tight_binding import TBLattice
t = 1.0
H = TBLattice(
units = [(1, 0, 0), (0, 1, 0)],
hopping = {
# nearest neighbour hopping -t
( 0,+1): -t * np.eye(2),
( 0,-1): -t * np.eye(2),
(+1, 0): -t * np.eye(2),
(-1, 0): -t * np.eye(2),
},
orbital_positions = [(0,0,0)]*2,
orbital_names = ['up', 'do'],
)
Starting on 1 Nodes at : 2019-06-05 17:03:01.636154
Dispersion \(\epsilon(\mathbf{k})\) and Fermi surface
The Fourier transform of the real space hopping to momentum space \(\mathbf{k}\) gives the momentum space dispersion \(\epsilon(\mathbf{k})\) and the alternative representation of \(H\) as
To visualize \(\epsilon(\mathbf{k})\), we can plot it along the high-symmetry path \(\Gamma - X - M - \Gamma\) in the Brillouin zone. We use the helper function k_space_path to build the \(\mathbf{k}\)-space path.
[3]:
kmesh = H.get_kmesh(n_k=(32, 32, 1))
e_k = H.fourier(kmesh)
G = [0.0, 0.0, 0.0]
X = [0.5, 0.0, 0.0]
M = [0.5, 0.5, 0.0]
paths = [(G, X), (X, M), (M, G)]
from triqs_tprf.lattice_utils import k_space_path
k_vecs, k_plot, k_ticks = k_space_path(paths, bz=H.bz)
e_k_interp = np.vectorize(lambda k : e_k(k)[0, 0].real, signature='(n)->()')
plt.plot(k_plot, e_k_interp(k_vecs))
plt.xticks(k_ticks, [r'$\Gamma$',r'$X$',r'$M$',r'$\Gamma$'])
plt.ylabel(r'$\epsilon(\mathbf{k})$')
plt.grid(True)
From this we can see that the bandwidth is \(4t\) (with \(t=1\)) and that we have a Van Hove singularity at the \(M\) point.To visualize the dispersion in the \((k_x, k_y)\) plane and the half-filled fermi surface where \(\epsilon(\mathbf{k}) = 0\) we interpolate e_k on a two-dimensional grid and use matplotlib.pyplot.contour to plot the zero contour.
[4]:
e_k = H.fourier(H.get_kmesh(n_k=(32, 32, 1)))
k = np.linspace(-0.5, 0.5, num=200) * 2. * np.pi
Kx, Ky = np.meshgrid(k, k)
e_k_interp = np.vectorize(lambda kx, ky : e_k([kx, ky, 0])[0, 0].real)
e_k_interp = e_k_interp(Kx, Ky)
plt.figure()
extent = (k.min(), k.max(), k.min(), k.max())
plt.imshow(e_k_interp, cmap=plt.get_cmap('RdBu'),
extent=extent, origin='lower')
plt.colorbar()
plt.contour(Kx, Ky, e_k_interp, levels=[0])
plt.xlabel(r'$k_x$'); plt.ylabel(r'$k_y$');
I.e. the Fermi surface has a diamond shape with perfect nesting of the fermi surface for the momentum transfer vectors \(\mathbf{k}=\pm(\pi, \pi)\) and \(\mathbf{k}=\pm(\pi, -\pi)\).
Bare generalized susceptibility \(\chi_0\)
The bare generalized (particle-hole) susceptibility \(\chi_0\) is defined as the bubble diagram of a propagating particle hole excitation, i.e.
where \(G\) is the single-particle Green’s function.
For a non-interacting system the single particle Green’s function \(G_0\) is known analytically
and the product of two Green’s functions can be performed analytically using the Lindhard formula
where \(f(\cdot)\) is the Fermi distribution function.
For reference this formula is implemented in TPRF, for details on the generalization to multiorbital and accidentally degenerate dispersion see the documentation for lindhard_chi00.
[5]:
from triqs_tprf.lattice import lindhard_chi00
from triqs.gf import MeshImFreq, Idx
wmesh_lind = MeshImFreq(beta=5.0, S='Boson', n_max=1)
chi00_wk_analytic = lindhard_chi00(e_k=e_k, mesh=wmesh_lind, mu=0.)
A computationally more efficient way of computing \(\chi_0\) is to perform the product of two single-particle Green’s functions in imaginary time and real space and then fourier transform back to frequency and momentum space, i.e.
To do this we first compute the single particle Green’s function of the lattice using the lattice Dyson equation
[6]:
from triqs_tprf.lattice import lattice_dyson_g0_wk
wmesh = MeshImFreq(beta=5.0, S='Fermion', n_max=30)
g0_wk = lattice_dyson_g0_wk(mu=0., e_k=e_k, mesh=wmesh)
And then perform the bubble product using the dedictated function
[7]:
from triqs_tprf.lattice_utils import imtime_bubble_chi0_wk
chi00_wk = imtime_bubble_chi0_wk(g0_wk, nw=1)
╔╦╗╦═╗╦╔═╗ ╔═╗ ┌┬┐┌─┐┬─┐┌─┐
║ ╠╦╝║║═╬╗╚═╗ │ ├─┘├┬┘├┤
╩ ╩╚═╩╚═╝╚╚═╝ ┴ ┴ ┴└─└
Two-Particle Response Function tool-box
beta = 5.0
nk = 1024
nw = 60
norb = 2
Approx. Memory Utilization: 0.01 GB
--> fourier_wk_to_wr
--> fourier_wr_to_tr
--> chi0_w0r_from_grt_PH (bubble in tau & r)
--> chi_wk_from_chi_wr (r->k)
Physical response function \(\chi_{S_z, S_z}\)
From the generalized susceptibility \(\chi_{abcd}(i\omega_n, \mathbf{k})\) we can obtain any two particle response function by contracting the rank 4 tensor indices with two quadratic operators. E.g. the spin-spin response along the z-axis is given by
[8]:
def chi_SzSz_contraction(chi):
""" Computes the trace Tr[ S_z \chi S_z ]"""
Sz = np.diag([+0.5, -0.5])
chi_SzSz = chi[0, 0, 0, 0].copy()
chi_SzSz.data[:] = np.einsum('wqabcd,ab,cd->wq', chi.data, Sz, Sz)[:, :]
chi_SzSz = chi_SzSz[Idx(0), :]
return chi_SzSz
To look at the structure of \(\chi_0\) we can plot it along the high symmetry path in the Brillouin zone (using the call operator of a triqs Green’s function on MeshBrZone).
[9]:
def interpolate_chi(chi, k_vecs):
assert( k_vecs.shape[1] == 3 )
chi_interp = np.zeros(
[k_vecs.shape[0]] + list(chi.target_shape), dtype=complex)
for kidx, (kx, ky, kz) in enumerate(k_vecs):
chi_interp[kidx] = chi((kx, ky, kz))
return chi_interp
With some additional plot helper functions we can now look at the spin-spin response
[10]:
def plot_chi_1D(chi, label=None):
chi_SzSz = chi_SzSz_contraction(chi)
chi_interp = interpolate_chi(chi_SzSz, 2*np.pi*k_vecs)
plt.plot(ks, chi_interp.real, label=label)
plt.grid()
plt.xticks(ticks=eps_k['K'], labels=[r'$\Gamma$',r'$X$',r'$M$',r'$\Gamma$'])
plt.title(r'Spin-response $\chi_{S_z S_z}(\mathbf{q}, \omega=0)$')
[11]:
plot_chi_1D(chi00_wk, label=r'$U=0$')
[12]:
def plot_chi(chi, label=None):
k = np.linspace(-0.75, 0.75, num=100) * 2.*np.pi
Kx, Ky = np.meshgrid(k, k)
k_vecs = np.vstack((Kx.flatten(), Ky.flatten(), 0*Kx.flatten())).T
chi_SzSz = chi_SzSz_contraction(chi)
chi_interp = interpolate_chi(chi_SzSz, k_vecs)
chi_interp = chi_interp.real.reshape(Kx.shape)
plt.imshow(
chi_interp, cmap=plt.get_cmap('magma'),
extent=(k.min(), k.max(), k.min(), k.max()),
origin='lower', vmin=0, vmax=chi_interp.max())
plt.title(label); plt.colorbar()
plt.xlabel(r'$k_x$'); plt.ylabel(r'$k_y$')
[13]:
label = r'Lindhardt response $\chi^{00}_{S_z S_z}(\mathbf{q}, \omega=0)$'
plot_chi(chi00_wk, label=label)
plt.tight_layout()
Random phase approximation (RPA)
The full Bethe-Salpeter equation in the particle-hole channel has the form
where all components are momentum and (three) frequency dependent quantities. In the random phase approximation (RPA) the vertex is apprximated by the (frequency independent) local interaction \(\Gamma^{(PH)} \approx V\).
The locality and frequency independence simplifies the BSE which formally is solved as
.
In TPRF, there is the dedicated routine for RPA (triqs_tprf.lattice.solve_rpa_PH) that can treat the general multiorbital case.
For generality, we take the nessecary steps to obtain the four index tensor \(V\) from the operator expression of the Hubbard interaction.
[14]:
from triqs.operators import n, c
H_int = n(0, 0) * n(0, 1)
fundamental_operators = [c(0, 0), c(0, 1)]
from triqs_tprf.OperatorUtils import quartic_tensor_from_operator
from triqs_tprf.OperatorUtils import quartic_permutation_symmetrize
V_int_abcd = quartic_tensor_from_operator(H_int, fundamental_operators)
V_int_abcd = quartic_permutation_symmetrize(V_int_abcd)
V_abcd = np.zeros_like(V_int_abcd)
from itertools import product
for a, b, c, d in product(list(range(V_abcd.shape[0])), repeat=4):
V_abcd[a, b, c, d] = V_int_abcd[b, d, a, c]
Now we can tune the interaction strength and compute the corresponding RPA generalized susceptibility.
[15]:
from triqs_tprf.lattice import solve_rpa_PH
chi_wk_vec = []
U_vec = np.arange(1.0, 5.0, 1.0)
for u in U_vec:
chi_wk = solve_rpa_PH(chi00_wk, u * V_abcd)
chi_wk_vec.append(chi_wk)
[16]:
plot_chi_1D(chi00_wk, label=r'$U=0$')
for U, chi_wk in zip(U_vec, chi_wk_vec):
plot_chi_1D(chi_wk, label=r'$U=%2.2f$' % U)
plt.legend();
We see that the system has an instability at the \(M\) point
[17]:
plt.figure(figsize=(7, 3))
subp = [1, 2, 1]
plt.subplot(*subp); subp[-1] += 1
label = r'Lindhardt response $\chi^{00}_{S_z S_z}(\mathbf{q}, \omega=0)$'
plot_chi(chi00_wk, label=label)
plt.subplot(*subp); subp[-1] += 1
label=r'RPA U=%1.0f $\chi^{00}_{S_z S_z}(\mathbf{q}, \omega=0)$' % U_vec[-1]
plot_chi(chi_wk_vec[-1], label=label)
plt.tight_layout()
