Second-order perturbation

Second-order perturbation#

Theory#

We consider the Hubbard model on a square lattice at half filling, whose Hamiltonian is given by

H = - t \sum_{⟨ i, j ⟩} c_{i σ}^{†} c_{j σ} + U \sum_{i} n_{i ↑} n_{i ↓} - μ \sum_{i} (n_{i ↑} + n_{i ↓}),

where $μ = U / 2$ , $t (= 1)$ a hopping integral, $c_{i σ}^{†} (c_{i σ})$ a creation (annihilation) operator for an electron with spin $σ$ at site $i$ , $μ$ chemical potential. The non-interacting band dispersion is given by

ϵ (k) = - 2 (\cos k_{x} + \cos k_{y}),

where $k = (k_{x}, k_{y})$ . Hereafter, we ignore the spin dependence.

The corresponding non-interacting Green’s function is given by

G (i ν, k) = \frac{1}{i ν - ϵ (k) - μ},

where $ν$ is a fermionic frequency.

In the imaginary-time domain, the self-energy reads

Σ (τ, r) = U^{2} G^{2} (τ, r) G (β - τ, r),

where $r$ represents a position in real space. Note that the first-order term of the self-energy, i.e., a Hartree term, is absorbed into the chemical potential. Thus, we take $μ = 0$ for half filling. The Fourier transform of Green’s functions and self-energies between the momentum and real-space domains is defined as follows:

A (i ν, r) = \frac{1}{N} \sum_{k} e^{- i k \cdot r} A (i ν, k),

A (i ν, k) = \sum_{r} e^{i k \cdot r} A (i ν, r),

where $A$ is either a Green’s function or a self-energy.

We use regular 2D meshes for the real-space and momen spaces. More specifically, a point in the real space $r$ is represented by three integers as $r = (r_{1}, r_{2})$ , where $0 \leq r_{α} \leq L - 1$ . A point in the momentum space $k$ is represented by three integers as $k = (2 π k_{1} / L, 2 π k_{2} / L)$ , where $0 \leq k_{α} \leq L - 1$ .

G_{m, r_{1}, r_{2}} = \frac{1}{N} F (G_{m, k_{1}, k_{2}}),

G_{m, k_{1}, k_{2}} = N F^{- 1} (G_{m, r_{1}, r_{2}}) .

Here, $F$ denotes two-dimentional discrete Fourier transform (DFT) and $F^{- 1}$ its inverse transform, which are defined as

x_{n_{1}, n_{2}} = \sum_{m_{1}, m_{2} = 0}^{L - 1} x_{m_{1}, m_{2}} e^{- \frac{2 π \sum_{α} n_{α} m_{α}}{L} i} = F (x_{m_{1}, m_{2}}),

y_{m_{1}, m_{2}} = \frac{1}{L^{2}} \sum_{n_{1}, n_{2} = 0}^{L - 1} x_{n_{1}, n_{2}} e^{\frac{2 π \sum_{α} n_{α} m_{α}}{L} i} = F^{- 1} (x_{n_{1}, n_{2}}) .

The DFT can be performed efficiently by the fast Fourier transform algorithm, which is implemented in many numerical libraries such as FFTW and FFTPACK. The computational complexity of FFT scales as $O (L^{2} \log L)$ , which is considerably smaller than that of the naive implementation [ $O (L^{3})$ ] for large $L$ .

Implementation#

import numpy as np
%matplotlib inline
import matplotlib.pyplot as plt
plt.rcParams['font.size'] = 15
import sparse_ir
from numpy.fft import fftn, ifftn

Step 1: Generate IR basis and associated sampling points#

The following code works only for a single-orbital problem. numpy.linalg.inv can be used instead for computing matrix-valued Green’s functons, i.e., for multi-orbital systems. The condition numbers are around 180 and 250 for imaginary-time and Matsubara-frequency samplings, repspectively. This means that we may lose two significant digits when fitting numerical data on the sampling points. Thus, the results of fits are expected to have an accuracy of five significant digits ( $10^{- 7} \times 10^{2} = 10^{- 5}$ ). If you want more significant digits, you can decrease the cutoff.

lambda_ = 1e+5
beta = 1e+3
eps = 1e-7
# Number of divisions along each reciprocal lattice vector
# Note: For a smaller nk (e.g., 64), an unphysical small structures appear in the self-energy at low frequencies.
#nk_lin = 256
nk_lin = 256
U = 2.0 # Onsite repulsion

wmax = lambda_/beta
basis = sparse_ir.FiniteTempBasis("F", beta, wmax, eps=eps)
L = basis.size

# Sparse sampling in tau
smpl_tau = sparse_ir.TauSampling(basis)
ntau = smpl_tau.sampling_points.size
print("cond (tau): ", smpl_tau.cond)

# Sparse sampling in Matsubara frequencies
smpl_matsu = sparse_ir.MatsubaraSampling(basis)
nw = smpl_matsu.sampling_points.size
print("cond (matsu): ", smpl_matsu.cond)

cond (tau):  162.99257542854215
cond (matsu):  465.80046864182265

Step 2#

Compute the non-interacting Green’s function on a mesh. The Green’s function computed at $k = Γ = (0, 0)$ is plotted below.

kps = (nk_lin, nk_lin)
nk = np.prod(kps)
nr = nk
kgrid = [2*np.pi*np.arange(kp)/kp for kp in kps]
k1, k2 = np.meshgrid(*kgrid, indexing="ij")
ek = -2*(np.cos(k1) + np.cos(k2))
print(k1.shape, k2.shape, ek.shape)

iv = 1j*np.pi*smpl_matsu.sampling_points/beta
gkf = 1.0 / (iv[:,None] - ek.ravel()[None,:])

plt.plot(smpl_matsu.sampling_points*np.pi/beta, gkf[:,0].imag, label=r"$k=\Gamma$", marker="x")
plt.xlabel(r"$\nu$")
plt.ylabel(r"$\mathrm{Im}~G(\mathrm{i}\nu)$")
plt.legend()
plt.tight_layout()
plt.show()

(256, 256) (256, 256) (256, 256)

../_images/1c01f8084475effe70b5fcd9e9a9fe9320b786c63c306070de301fc791b90e68.png

Step 3#

Transform the Green’s function to sampling times.

# G(l, k): (L, nk)
gkl = smpl_matsu.fit(gkf, axis=0)

assert gkl.shape == (L, nk)

plt.semilogy(np.abs(gkl[:,0]))
plt.semilogy(basis.s, ls="--")
plt.xlabel(r"$l$")
plt.ylabel(r"$|G(l, \Gamma)|$")
plt.show()

../_images/c9ffeead806a88319fda0dba95158c908ff6c1c6919ec9502b2a8ccedf9c0adf.png

gkt = smpl_tau.evaluate(gkl)

assert gkt.shape == (ntau, nk)

plt.plot(smpl_tau.sampling_points, gkt[:,0].real, label=r'$\Gamma$', marker="x")
plt.plot(smpl_tau.sampling_points, 
    gkt.reshape(-1,*kps)[:,nk_lin//2,nk_lin//2].real, label=r'$M$', marker="+")
plt.xlabel(r"$\tau$")
plt.ylabel(r"$\mathrm{Re}~G(\tau, k)$")
plt.legend()
plt.tight_layout()
plt.show()

../_images/7f0bcce84e2e3247041ca09fbbe2c42ceb3810ebd8cd6efe5553bb55a116e6ea.png

Step 4#

Transform the Green’s function to the real space and evaluate the self-energy

# Compute G(tau, r): (ntau, nk)
#  (1) Reshape gkt into shape of (ntau, nk_lin, nk_lin).
#  (2) Apply FFT to the axes 1, 2.
#  (3) Reshape the result to (ntau, nk)
# G(tau, k): (ntau, nk)
grt = fftn(gkt.reshape(ntau, *kps), axes=(1,2)).reshape(ntau, nk)/nk

plt.plot(smpl_tau.sampling_points, grt[:,0].real, label='$r=(0,0)$')
plt.xlabel(r"$\tau$")
plt.ylabel(r"$\mathrm{Re}~G(\tau, r)$")
plt.legend()
plt.tight_layout()
plt.show()

../_images/5ee0c03faf6ae2384efddaf904e3a6bba7fa5c1271a025cfca30121cfe0d405c.png

Compute the second-order term of the self-energy on the sampling points

# Sigma(tau, r): (ntau, nr)
srt = U*U*grt*grt*grt[::-1,:]

plt.plot(smpl_tau.sampling_points, srt[:,0].real, label='$r=(0,0)$', marker="x")
plt.xlabel(r"$\tau$")
plt.ylabel(r"$\mathrm{Re}~\Sigma(\tau, r)$")
plt.legend()
plt.tight_layout()
plt.show()

../_images/21b5ac0d6fad78a1e5305b682fc3f7040c0335ae56a4f839e176796092c4a48c.png

Step 5#

Transform the self-energy to the IR basis and then transform it to the momentum space

# Sigma(l, r): (L, nr)
srl = smpl_tau.fit(srt)
assert srl.shape == (L, nr)

plt.semilogy(np.abs(srl[:,0]), label='$r=(0,0)$')
plt.xlabel(r"$l$")
plt.ylabel(r"$|\Sigma(l, r)|$")
plt.ylim([1e-8,1])
plt.legend()
plt.show()

../_images/3947b1d5faed634399eedeb4c67a5a6586febbdc73e974bf95717f8e118e8fcc.png

# Sigma(l, k): (L, nk)
srl = srl.reshape(L, *kps)
skl = ifftn(srl, axes=(1,2)) * nk
skl = skl.reshape(L, *kps)

#plt.semilogy(np.max(np.abs(skl),axis=(1,2)), label="max, abs")
plt.semilogy(np.abs(skl[:,0,0]), label="$k=\Gamma$")
plt.semilogy(basis.s, label="Singular values", ls="--")
plt.xlabel(r"$l$")
plt.ylabel(r"$|\Sigma(l, k)|$")
plt.legend()
plt.tight_layout()
plt.show()

../_images/6379f3d236ff7e7ab0e5fa13a3a7df54beec4650af02547216259ad5ab1c8d9b.png

Step 6#

Evaluate the self-energy on sampling frequencies

sigma_iv = smpl_matsu.evaluate(skl, axis=0)
plt.plot(smpl_matsu.sampling_points*np.pi/beta, sigma_iv[:,0,0].imag, marker='x', label=r"$k=\Gamma$")
plt.xlabel(r"$\nu$")
plt.ylabel(r"$\mathrm{Im}~\Sigma(\mathrm{i}\nu)$")
plt.legend(frameon=False)
plt.tight_layout()
plt.show()

../_images/a974a6d967fceabe044780f1d140358a73ae9d0960ec1ad92bb1ba2026ad4751.png

You can evaluate the self-energy on arbitrary frequencies.

my_freqs = 2*np.arange(-10000, 10000, 20) + 1 # odd
smpl_matsu2 = sparse_ir.MatsubaraSampling(basis, my_freqs)
res = smpl_matsu2.evaluate(skl, axis=0)
plt.plot(my_freqs*np.pi/beta, res[:,0,0].imag, marker='x', label=r"$k=\Gamma$")
plt.xlabel(r"$\nu$")
plt.ylabel(r"$\mathrm{Im}~\Sigma(\mathrm{i}\nu)$")
plt.legend(frameon=False)
plt.tight_layout()
plt.show()

/home/jovyan/.venv/default/lib/python3.11/site-packages/sparse_ir/sampling.py:129: ConditioningWarning: Sampling matrix is poorly conditioned (kappa = 5.5e+14)
  warn(f"Sampling matrix is poorly conditioned "

../_images/d8b7949c9bb139f2aea821415cb7dc488fbe901b0b9a7fec0034aadf2ef373b5.png