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

import numpy as np
from triqs.gf import *

The Two-Particle Self Consistency (TPSC) method

The Two-Particle Self Consistent (TPSC) method is an approach to solving the single-band Hubbard Model first introduced by Vilk and Tremblay (https://arxiv.org/abs/cond-mat/9702188).

Theory

Among the simplest approximations to self-energy and susceptibilities are the Hartree-Fock (HF) and Random-Phase-Approximation (RPA), which, in the single-band Hubbard model are

\[\Sigma^\text{HF}_\sigma(i\omega_n, \mathbf{k}) = U \langle n_{-\sigma}\rangle\]

\[\chi^\text{RPA}_\text{sp}(i\omega_n, \mathbf{k}) = \frac{\chi_0(i\omega_n, \mathbf{k})}{1 - \frac{U}{2}\chi_0(i\omega_n, \mathbf{k})}, \qquad \chi^\text{RPA}_\text{ch}(i\omega_n, \mathbf{k}) = \frac{\chi_0(i\omega_n, \mathbf{k})}{1 + \frac{U}{2}\chi_0(i\omega_n, \mathbf{k})}\]

The HF-self-energy has the immediate disadvantage of being constant and thus absorbed in the chemical potential.

For the susceptibilities, there exist two sum rules

\[\begin{split}\begin{align} \frac{T}{N}\sum_{i\omega_n, \mathbf{k}}{\chi_\text{sp}(i\omega_n, \mathbf{k})} &= n - 2\langle n_\uparrow n_\downarrow\rangle \\ \frac{T}{N}\sum_{i\omega_n, \mathbf{k}}{\chi_\text{ch}(i\omega_n, \mathbf{k})} &= n + 2\langle n_\uparrow n_\downarrow\rangle - n^2 \end{align}\end{split}\]

which are violated by the RPA-approximations. Also, since the denominator of \(\chi_\text{sp}\) can vanish, RPA predicts a magnetically ordered phase, which, in two dimensions is prohibited by the Mermin-Wagner theorem.

To fix this, TPSC introduces RPA-like susceptibilities with screened interaction vertices \(U_\text{ch}\) and \(U_\text{sp}\)

\[\chi^\text{TPSC}_\text{sp}(i\omega_n, \mathbf{k}) = \frac{\chi_0(i\omega_n, \mathbf{k})}{1 - \frac{U_\text{sp}}{2}\chi_0(i\omega_n, \mathbf{k})}, \qquad \chi^\text{TPSC}_\text{ch}(i\omega_n, \mathbf{k}) = \frac{\chi_0(i\omega_n, \mathbf{k})}{1 + \frac{U_\text{ch}}{2}\chi_0(i\omega_n, \mathbf{k})}\]

\(U_\text{ch}\) and \(U_\text{sp}\) are chosen such that the sum rules are fulfilled; for this to work, we need to now the double occupancy \(\langle n_\uparrow n_\downarrow\rangle\) for which TPSC makes the Ansatz

\[U_\text{sp} \langle n_\uparrow\rangle\langle n_\downarrow\rangle = U \langle n_\uparrow n_\downarrow\rangle\]

Having obtained the renormalized vertices and susceptibilities, TPSC then approximates the self-energy as

\[\Sigma^\text{TPSC}_\sigma(i\omega_n, \mathbf{k}) = U \langle n_{-\sigma}\rangle + \frac{U}{8}\frac{T}{N}\sum_{i\nu_m, \mathbf{q}}{\left[ 3 U_\text{sp}\chi_\text{sp}(i\nu_m, \mathbf{q}) + U_\text{ch}\chi_\text{ch}(i\nu_m, \mathbf{q}) \right] G_{0\sigma}(i\nu_m+i\omega_n, \mathbf{q}+\mathbf{k})}\]

Triangular lattice

In the following, we apply the TPSC method to a triangular two-dimensional Hubbard model. The triangular lattice in 2 dimensions is spanned by the primitive lattice vectors

\[\mathbf{a}_1 = a\hat{x}, \qquad \mathbf{a}_2 = \frac{a}{2} \hat{x} + \frac{\sqrt{3}a}{2}\hat{y}\]

For simplicity, we set \(a = 1\); in TRIQS we can create such a lattice as follows:

from triqs.lattice.tight_binding import TBLattice

t = -1.0
mu = -2.0

units = [(1, 0, 0), (1/2, np.sqrt(3)/2, 0)]

hoppings = {(+1,+0) : [[t]],
            (-1,+0) : [[t]],
            (+0,+1) : [[t]],
            (+0,-1) : [[t]],
            (+1,-1) : [[t]],
            (-1,+1) : [[t]],
            (+0,+0) : [[mu]]}

H = TBLattice(units=units, hoppings=hoppings)

Warning: could not identify MPI environment!

Starting serial run at: 2025-04-09 11:27:10.617834

TRIQS also allows us to calculate the dispersion relation and store it in an object e_k. We can plot this dispersion along high symmetry lines:

kmesh = H.get_kmesh(n_k=(64, 64, 1))

e_k = H.fourier(kmesh)

# define the high-symmetry points of the triangular lattice
G = [0.0, 0.0, 0.0]
M = [1/2, 0.0, 0.0]
K = [1/3, -1/3, 0.0]

paths = [(G, M), (M, K), (K, G)]

from triqs_tprf.lattice_utils import k_space_path
k_vecs, k_plot, k_ticks = k_space_path(paths, bz=H.bz, relative_coordinates=False)
k_vecs_rel, k_plot, k_ticks = k_space_path(paths, bz=H.bz, relative_coordinates=True)

e_k_interp_ref = [ H.tb.fourier(k)[0,0].real for k in k_vecs_rel]
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.plot(k_plot, e_k_interp_ref, '--')
plt.xticks(k_ticks, [r'$\Gamma$',r'$M$',r'$K$',r'$\Gamma$'])
plt.ylabel(r'$\epsilon(\mathbf{k})$')
plt.grid(True)

Since it is a triangular lattice, the dispersion is not particle-hole symmetric. As we can see, this dispersion has a saddle point at the Fermi-level. In 2D, this will lead to a logarithmic divergence in the density of states, called a van Hove-singularity. The density of states can be calculated with the dos function:

from triqs.lattice.tight_binding import dos

DOS = dos(tight_binding=H.tb, n_kpts=2**10, n_eps=200, name='dos')[0]
plt.plot(DOS.eps, DOS.rho)
plt.xlabel('$\\varepsilon$')
plt.ylabel('$\\rho(\\varepsilon)$')

Text(0, 0.5, '$\\rho(\\varepsilon)$')

To get a better feeling for the dispersion, we can plot it as a function of \(\mathbf{k}\) within the first Brillouin-Zone. This also allows us to visualize the Fermi-Surface, which in this case is a hexagon (dashed lines):

e_k = H.fourier(H.get_kmesh(n_k=512))
k = np.linspace(-1.0, 1.0, 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()

g_abs = H.bz.units.T @ G
k_abs = H.bz.units.T @ K
m_abs = H.bz.units.T @ M

abs_path = [ g_abs, k_abs, m_abs, g_abs ]

K_X = [ p[0] for p in abs_path ]
K_Y = [ p[1] for p in abs_path ]

plt.plot(K_X, K_Y, '-', color='white')

plt.plot(g_abs[0], g_abs[1], 'co', label='G')
plt.plot(k_abs[0], k_abs[1], 'ro', label='K')
plt.plot(m_abs[0], m_abs[1], 'mo', label='M')

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], linestyles='dashed')
plt.xlabel(r'$k_x$'); plt.ylabel(r'$k_y$');
plt.legend(loc='upper left');

Here, the Fermi-Surface is indicated by dashed lines. In this case, it is hexagonal; this will lead to pronounced Fermi-Surface-Nesting, as described below.

Non-interacting susceptibility \(\chi_0\)

TRIQS TPRF offers functions to calculate non-interacting Green’s functions and susceptibilities from a known dispersion e_k. Here, we calculate the static susceptibility \(\chi_0(i\omega_n=0, \mathbf{k})\):

from triqs_tprf.lattice import lattice_dyson_g0_wk
from triqs_tprf.lattice_utils import imtime_bubble_chi0_wk

beta = 100.

wmesh = MeshDLRImFreq(
    beta=beta, statistic='Fermion',
    eps=1e-6, # Tune to desired accuracy
    w_max=10.0, # Tune to energy range in problem
    )

g0_wk = lattice_dyson_g0_wk(mu=0., e_k=e_k, mesh=wmesh)

chi0_wk = imtime_bubble_chi0_wk(g0_wk, nw=1, verbose=True)
iw0 = chi0_wk.mesh.components[0][0] # get zero Matsubara point to evaluate later

╔╦╗╦═╗╦╔═╗ ╔═╗  ┌┬┐┌─┐┬─┐┌─┐
 ║ ╠╦╝║║═╬╗╚═╗   │ ├─┘├┬┘├┤
 ╩ ╩╚═╩╚═╝╚╚═╝   ┴ ┴  ┴└─└
Two-Particle Response Function tool-box

beta  = 100.0
nk    = 262144
nw    = 32
norb  = 1

Approx. Memory Utilization: 0.38 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)

We can now plot the static susceptibility:

chi0_k_interp = np.vectorize(lambda kx, ky : chi0_wk(Idx(0), [kx, ky, 0])[0, 0].real)
chi0_k_interp = chi0_k_interp(Kx, Ky)

plt.figure()

extent = (k.min(), k.max(), k.min(), k.max())
plt.imshow(chi0_k_interp, cmap=plt.get_cmap('RdBu'),
           extent=extent, origin='lower')
plt.colorbar()
plt.contour(Kx, Ky, e_k_interp, levels=[0], linestyles='dashed')
plt.xlabel(r'$k_x$'); plt.ylabel(r'$k_y$');
plt.plot(np.pi, np.sqrt(3)*np.pi, 'o')
plt.plot(np.pi, -np.pi/np.sqrt(3), 'o')

[<matplotlib.lines.Line2D at 0x7f419c4cdb10>]

The analytic form of this quantity is known:

\[\chi_0^\text{stat}(\mathbf{k}) = -\frac{2}{N} \sum_{\mathbf{q}}^{\text{1.BZ}} { \frac{f_\mu(\varepsilon(\mathbf{q})) - f_\mu(\varepsilon(\mathbf{q} + \mathbf{k}))} {\varepsilon(\mathbf{q}) - \varepsilon(\mathbf{q} + \mathbf{k})}},\]

where \(f_\mu\) is the Fermi-Dirac-statistic. From this expression, we can see that the peaks of \(\chi_0\) correspond to all possible particle-hole excitations. Whenever there exists a wave vector \(\mathbf{k}\) which maps a large portion of the Fermi-surface to another part of the Fermi-surface, this will result in a pronounced peak at that wave vector. This is known as Fermi Surface Nesting (see also https://triqs.github.io/triqs/latest/userguide/python/tutorials/TwoParticleResponse/solutions/01s-Fermi_surface_nesting.html).

The peak in the upper right corner, for instance, can be understood as a result of fermi surface nesting with a wave-vector of \(\mathbf{k} = (\pi, \sqrt{3}\pi)^T\). Mapped back into the 1. BZ, this results in the peak at \((\pi, -\pi / \sqrt{3})^T\), and so forth.

The broadening of the peaks can be understood by noticing that if one slightly tilts its nesting vector, the resulting \(\mathbf{k}\)-vector will still map part of the fermi surface to another.

To get an even better feel for this, we can plot the static susceptibility along the high-symmetry paths in the 1. BZ:

chi0_stat_k = np.array([chi0_wk(Idx(0), k)[0,0,0,0].real for k in k_vecs])

fig, axes = plt.subplots(nrows=1, ncols=1)

axes.plot(k_plot, chi0_stat_k, label='$\\chi_0$')
axes.legend(bbox_to_anchor=(1,1))
axes.set_xticks(k_ticks, [r'$\Gamma$',r'$M$',r'$K$',r'$\Gamma$'])
axes.set_xlabel('$\\mathbf{k}$')
axes.set_ylabel('$\\chi_0(i\\omega_n=0, \\mathbf{k})$')
axes.grid(True)
axes.set_box_aspect(1)

Interacting charge (\(\chi_{ch}\)) and spin response (\(\chi_{sp}\))

TRIQS TPRF has a tpsc_solver-class that can be used to do TPSC-calculations for one-band Hubbard models. The solver has to be supplied with the target density, the strength of the Hubbard interaction, an imaginary frequency mesh and a dispersion relation. We will run the calculation at half-filling, \(n=1\), because there the Hubbard-Model becomes anti-ferromagnetic, which will be visible in the spin susceptibility.

from triqs_tprf.tpsc_solver import tpsc_solver

U = 2.0
S = tpsc_solver(n=1, U=U, wmesh=wmesh, e_k=e_k)

S.solve()

╔╦╗╦═╗╦╔═╗ ╔═╗  ┌┬┐┌─┐┌─┐┌─┐
 ║ ╠╦╝║║═╬╗╚═╗   │ ├─┘└─┐│
 ╩ ╩╚═╩╚═╝╚╚═╝   ┴ ┴  └─┘└─┘
Two-Particle Self-Consistent method

Using the TPSC-Ansatz for the double occupancy.

len(wmesh) = 32
nk = 262144
beta = 100.0
n = 1
U = 2.0

--> lattice_dyson_g0_wk

--> imtime_bubble_chi0_wk

--> get Usp

--> get Uch

--> get chisp_wk and chich_wk

--> get double occupancy

Summary first level approximation:
    Usp = 1.662732184900791, Uch = 2.9234160464925343
    <n_up*n_down> = 0.20784152311259887

FIXME need custom function to make chi DLR (target_shape!)

def make_chi_dlr(chi_wk):
    temp = Gf(mesh=chi_wk.mesh, target_shape=(1,1))
    temp.data[:,:,0,0] = chi_wk.data[:,:,0,0,0,0]
    return make_gf_dlr(temp)

We now plot the non-interacting, and the TPSC-charge- and spin-susceptibilities along high-symmetry paths in the 1. BZ. This allows us to easily compare them:

chi0_stat_k = make_chi_dlr(S.chi0_wk)(iw0, all)
chi0_stat_k = np.array([chi0_stat_k(k)[0,0].real for k in k_vecs_rel])

chich_stat_k = make_chi_dlr(S.chich_wk)(iw0, all)
chich_stat_k = np.array([chich_stat_k(k)[0,0].real for k in k_vecs_rel])

chisp_stat_k = make_chi_dlr(S.chisp_wk)(iw0, all)
chisp_stat_k = np.array([chisp_stat_k(k)[0,0].real for k in k_vecs_rel])

fig, axes = plt.subplots(nrows=1, ncols=2, figsize=(10, 10))

axes[0].plot(k_plot, chi0_stat_k, label='$\\chi_0$')
axes[0].plot(k_plot, chich_stat_k, label='$\\chi_{ch}$')
axes[0].plot(k_plot, chisp_stat_k, label='$\\chi_{sp}$')
axes[0].set_title('Absolute comparison')
axes[0].legend(bbox_to_anchor=(1,1))
axes[0].set_xticks(k_ticks, [r'$\Gamma$',r'$M$',r'$K$',r'$\Gamma$'])
axes[0].set_xlabel('$\\mathbf{k}$')
axes[0].grid(True)
axes[0].set_box_aspect(1)

axes[1].plot(k_plot, chi0_stat_k - chi0_stat_k[0], label='$\\chi_0$')
axes[1].plot(k_plot, chich_stat_k - chich_stat_k[0], label='$\\chi_{ch}$')
axes[1].plot(k_plot, chisp_stat_k - chisp_stat_k[0], label='$\\chi_{sp}$')
axes[1].set_title('Relative comparison')
axes[1].legend(bbox_to_anchor=(1,1))
axes[1].set_xticks(k_ticks, [r'$\Gamma$',r'$M$',r'$K$',r'$\Gamma$'])
axes[1].set_xlabel('$\\mathbf{k}$')
axes[1].grid(True)
axes[1].set_box_aspect(1)

plt.tight_layout()

In the left image, we have plotted the absolute values of the susceptibilities; this shows that, when compared to the bare bubble, the spin fluctuations are increased in the interacting system, whereas charge fluctuations are reduced.

In the right figure, we have plotted all three susceptibilities on the same scale; this shows as that the spin susceptibility has the strongest peak, whereas charge fluctuations are much more uniform across the 1. BZ. This is, of course, exactly what we should expect in an anti-ferromagnetic system.

At the same time, the pronounced peak at the M-point due to Fermi-surface nesting is still visible in all three susceptibilities. It is important to note that the TPSC method respects the Mermin-Wagner theorem which forbids magnetic order in 2D systems; we see that this is indeed the case by noticing that the TPSC susceptibilities do not introduce additional peaks when compared to the bare bubble.