The Two-Particle Self Consistent method (TPSC)

The Two-Particle Self Consistent method (TPSC) is a non-perturbative modification of the first order interaction expansion for single-particle and two-particle quantities of the Hubbard model introduced by Vilk and Tremblay, see ArXiV:9702188. Compared to simpler approximations it has the essential feature of obeying the Mermin-Wagner theorem forbidding symmetry breaking at finite temperatures in two dimensions.

The triangular lattice Hubbard model

In this tutorial we use the TPSC solver triqs_trpf.tpsc_solver to determine the properties of the Hubbard model on the triangular lattice

\[H = H_{t} + H_U = -t \sum_{\langle i, j \rangle, \sigma} c^\dagger_{i \sigma} c_{j\sigma} + U \sum_{i} c^\dagger_{i\uparrow} c_{i\uparrow} c^\dagger_{i \downarrow} c_{i\downarrow}\]

where \(\langle i, j \rangle\) denote nearest neighbous, \(t\) is the nearest neighbour hopping, \(U\) the local Hubbard interaction, and \(\sigma \in \{\uparrow, \downarrow\}\).

The hopping \(t\) will be used as the unit of energy and set to unity (\(t = 1\)) and the chemical potential \(\mu = 2t\) will be used.

Background

First order Hartree-Fock (HF) and Random-Phase Approximation (RPA)

The first order approximation for single-particle Green’s function is the Hartree-Fock (HF) approximation, which in the Hubbard model gives a static and momentum independent single-particle self-energy

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

The two-particle response in terms of the charge- (\(\chi_{\text{ch}}\)) and spin-susceptibility (\(\chi_{\text{sp}}\)) are to first order in the interaction \(U\) given by the Random-Phase-Approximation (RPA)

\[\begin{split}\chi_\text{sp}(i\omega_n, \mathbf{k}) \approx \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})} \, \\ \chi_\text{ch}(i\omega_n, \mathbf{k}) \approx \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})} \, .\end{split}\]

The differing signs in the denominators cause the charge- and spin-response to behave very differently. Due to the negaitve sign in the denominator of the spin-response the RPA response diverges at some interaction, in violation of the Mermin-Wagner theorem.

The problem with unpysical ordering in RPA can be connected to the violation of the 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}\]

the charge- and spin-susceptibility.

The Two-Particle Self Consistent method (TPSC)

The central idea with the Two-Particle Self Consistent method (TPSC) is to retain the simple mathematical form of RPA

\[\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})} \, ,\]

while fulfilling the two sum rules using by introducing screened static effective interactions \(U_{\text{sp}}\) and \(U_{\text{ch}}\). The screened interactions are in turn self-consistently tuned as to fulfill the two sum rules.

It is the resulting screening of the spin-vertex \(U_{\text{sp}}\) that reinstates the the Mermin-Wagner theorem, reducing \(U_{\text{sp}}\) so that the denominator in \(\chi^{\text{TPSC}}_\text{sp}\) never goes to zero (at finite temperature).

To get a closed set of equations we require the expectation value of the double occupancy \(\langle n_\uparrow n_\downarrow\rangle\), and in is original form this is obtained in TPSC using the Ansatz

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

However, there are extensions of TPSC that use higher order methods (like DMFT) to get a better approximation for the double occupancy.

Once the TPSC spin- and charge-response is known it is also possible to compute an improved approximation of the single-particle self energy using first order shell-diagrams with the effective spin- and charge-interaction

\[\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})} \, .\]

The resulting self-energy is momentum dependent and in this form TPSC was the first approximation that captures the cold- and hot-spots along the Fermi surface in the half-filled Hubbard model on the square lattice [REF?].

Calculation using Triqs and TPRF

Using Triqs, TPRF and the triqs_trpf.tpsc_solver we can now setup the triangular lattice Hubbard model Hamiltonian and obtain the self-consistent TPSC solution.

Tight binding Hamiltonian

The triangular lattice is two-dimensional and 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}\]

where \(a\) is the lattice spacing which we set to unity (\(a=1\)), for simplicity. In TRIQS the hopping part \(H_t\) of the Hamiltonian can be represented by a tight-binding lattice instance (TBLattice).

The Fourier transform of the hopping Hamiltonian \(H_t\) to momentum-space gives the diagonal representation

\[H_t = \sum_{\mathbf{k}\sigma} \epsilon(\mathbf{k}) \, c^\dagger_{\mathbf{k} \sigma} c_{\mathbf{k} \sigma}\]

where \(\epsilon(\mathbf{k})\) is the so called dispersion relation. In TRIQS the dispersion relation can be obtained from the tight binding Hamiltonian as a Green’s function object e_k with a single k-space mesh. The dispersion \(\epsilon(\mathbf{k})\) can be visualized by plotting its value along the high symmetry path \(\Gamma\) - \(M\) - \(K\) - \(\Gamma\) in the Brillouin zone.

t = 1.0 # Nearest-neighbour hopping strength
mu = 2.0 # Chemical potential

from triqs.lattice.tight_binding import TBLattice

H_t = TBLattice(units=[(1, 0, 0), (1/2, np.sqrt(3)/2, 0)],
    hoppings = { coord : [[-t]] for coord in [(+1,+0), (-1,+0), (+0,+1), (+0,-1), (+1,-1), (-1,+1)]})

e_k = H_t.fourier(H_t.get_kmesh(n_k=96))

Warning: could not identify MPI environment!

Starting serial run at: 2025-06-01 10:33:03.239174

from triqs_tprf.lattice_utils import k_space_path

# High-symmetry points of the triangular lattice
G, M, K = [0.0, 0.0, 0.0], [1/2, 0.0, 0.0], [1/3, -1/3, 0.0]
k_vecs, k_plot, k_ticks = k_space_path([(G, M), (M, K), (K, G)], bz=H_t.bz, relative_coordinates=False)

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, mu + 0*k_plot, '--', lw=2, color='black', label=r'Chemical potential $\mu$')
plt.xticks(k_ticks, [r'$\Gamma$',r'$M$',r'$K$',r'$\Gamma$']); plt.legend(loc='best'); plt.grid(True);
plt.xlabel(r'Momentum high-symmetry path'); plt.ylabel(r'Dispersion relation $\epsilon(\mathbf{k})$');

In contrast to the square lattice the triangular lattice is not particle-hole symmetric around \(\epsilon(\mathbf{k}) = 0\). However, it still has a saddle point at the \(M\) point which generates a van Hove singularity at the energy \(\epsilon = 2\). (In anticipation of this we perviously selected a chemical potential \(\mu = 2\) tuned to this very van Hove singularity.)

The singularity is more prominently seen as a logarithmic divergence in the density of states (DOS) \(\rho(\epsilon)\) as a function of energy \(\epsilon\). The density of states is defined as

\[\rho(\epsilon) = \frac{1}{N_k} \sum_{\mathbf{k}} \delta( \epsilon(\mathbf{k}) - \epsilon )\]

where \(\delta(\cdot)\) is the Dirac delta function. For a TRIQS tight binding object the DOS can be computed using the dos function:

from triqs.lattice.tight_binding import dos

DOS = dos(tight_binding=H_t.tb, n_kpts=2**10, n_eps=200, name='dos')[0]

plt.plot([mu]*2, [0, 0.6], '--', color='black', lw=2, label='Chemical potential $\mu$')
plt.plot(DOS.eps, DOS.rho, lw=2)
plt.fill_between(DOS.eps, DOS.rho, 0, alpha=0.25)

plt.xlim([-7, 4]); plt.ylim(bottom=0, top=0.6); plt.grid(True); plt.legend(loc='best')
plt.xlabel('Energy $\\varepsilon$'); plt.ylabel('Density of states $\\rho(\\varepsilon)$');

Brillouin-zone and high symmetry path

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

# Two dimensional interpolation of the dispersion

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)(Kx, Ky)

plt.figure(figsize=(12,5))

plt.subplot(1, 2, 1)
plot_brillouin_zone()
plot_high_symmetry_path(linecolor='black')
plt.legend(loc='best'); plt.axis('square');
plt.xlim([-2*np.pi, +2*np.pi]); plt.ylim([-2*np.pi, +2*np.pi])
plt.xlabel(r'$k_x$'); plt.ylabel(r'$k_y$');

plt.subplot(1, 2, 2)
plt.contour(Kx, Ky, e_k_interp, levels=[mu], linestyles='dashed')
plt.plot([], [], color='black', linestyle='dashed', label='Fermi\nsurface')
plt.legend(loc='upper left');
plot_brillouin_zone()
plot_high_symmetry_path()
plt.imshow(e_k_interp - mu, cmap=plt.get_cmap('RdBu'),
           extent=(k.min(), k.max(), k.min(), k.max()), origin='lower', vmin=-6, vmax=6)
plt.colorbar(label=r'$\epsilon(\mathbf{k}) - \mu$')
plt.axis('square');
plt.xlim([-2*np.pi, +2*np.pi]); plt.ylim([-2*np.pi, +2*np.pi])
plt.xlabel(r'$k_x$'); plt.ylabel(r'$k_y$');

plt.tight_layout()

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 = 30.
#beta = 20.

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=mu, e_k=e_k, mesh=wmesh)

g0_w = Gf(mesh=g0_wk.mesh[0], target_shape=[1, 1])
g0_w.data[:] = np.sum(g0_wk.data, axis=1) / len(g0_wk.mesh[1])
n = density(g0_w)[0, 0].real
print(f'n = {n}')

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

n = 0.7500925150418974

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

beta  = 30.0
nk    = 9216
nw    = 27
norb  = 1

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)
MatsubaraFreq(n: 0, beta: 30.0, statistic: Boson)

Similarly as the dispersion we can now plot the static susceptibility in momentum space:

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('YlOrRd'),
           extent=extent, origin='lower', vmin=0)
plt.colorbar()

plt.contour(Kx, Ky, e_k_interp, levels=[mu], linestyles='dashed', colors='black')

plt.plot(np.pi, np.sqrt(3)*np.pi, 'o')
plt.plot(np.pi, -np.pi/np.sqrt(3), 'o')

plt.xlabel(r'$k_x$'); plt.ylabel(r'$k_y$');

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 the Triqs tutorial.

For instance, the peak in the upper right corner, 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 first Brillouin zone, 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.

While the two dimensional color plot of the susceptibility \(\chi_0\) is visually pleasing, to see quantitative features of the response it is more informative to plot \(\chi_0(\mathbf{k})\) along the high-symmetry path in the first Brillouin zone:

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_ylim(bottom=0)
axes.set_box_aspect(1)

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

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 \(\langle n \rangle\), the strength of the Hubbard interaction \(U\), an imaginary frequency mesh wmesh and a dispersion relation \(\epsilon(\mathbf{k})\).

We will run the calculation with the Fermi level at the van Hoove singularity \(\mu=2\) which corresponds to \(n=1.5\), because there the triangular lattice has perfect nesting.

from triqs_tprf.tpsc_solver import tpsc_solver

U = 2.0
S = tpsc_solver(n=2*n, 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) = 27
nk = 9216
beta = 30.0
n = 1.500185030083795
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.8102117702601144, Uch = 6.046972734724986
    <n_up*n_down> = 0.509247671995815

We now plot the bare (\(\chi_0\)), and the charge- (\(\chi_{ch}\)) and spin-susceptibilities (\(\chi_{sp}\)) along high-symmetry path:

chi0_k = get_zero_frequency(S.chi0_wk)
chi0_k_interp = [chi0_k(k).flatten().real for k in k_vecs]

chi_ch_k = get_zero_frequency(S.chich_wk)
chi_ch_k_interp = [chi_ch_k(k).flatten().real for k in k_vecs]

chi_sp_k = get_zero_frequency(S.chisp_wk)
chi_sp_k_interp = [chi_sp_k(k).flatten().real for k in k_vecs]

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

axes[0].plot(k_plot, chi0_k_interp, label='$\\chi_0$')
axes[0].plot(k_plot, chi_ch_k_interp, label='$\\chi_{ch}$')
axes[0].plot(k_plot, chi_sp_k_interp, label='$\\chi_{sp}$')
#axes[0].semilogy([], [])
axes[0].set_ylim(bottom=0)
axes[0].set_title('Linear scale')
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].set_ylabel('Re[$\chi(i\omega = 0, \\mathbf{k}$]')
axes[0].grid(True)
axes[0].set_box_aspect(1)

axes[1].plot(k_plot, chi0_k_interp, label='$\\chi_0$')
axes[1].plot(k_plot, chi_ch_k_interp, label='$\\chi_{ch}$')
axes[1].plot(k_plot, chi_sp_k_interp, label='$\\chi_{sp}$')
axes[1].semilogy([], [])
axes[1].set_ylim([1e-1, 1e2])
axes[1].set_title('Semilog y-axis')
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].set_ylabel('Re[$\chi(i\omega = 0, \\mathbf{k}$]')
axes[1].grid(True)
axes[1].set_box_aspect(1)

plt.tight_layout()

In the left figure, we have plotted the susceptibilities on a linear scale, and to the right on a logarithmic y-axis. When compared to the bare bubble \(\chi_0\), one sees that the spin fluctuations are enhanced \(\chi_{sp} > \chi_0\) in the interacting system, whereas charge fluctuations are suppressed \(\chi_{ch} < \chi_0\).

The pronounced peak at the M-point, due to Fermi-surface nesting, is 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 at finite temperature. This effect is only seen when reducing temperature, where the tendency of divergence in the spin response is suppressed by the reduction in in the screened interaction \(U_{sp}\) (not shown).