The plain and simple parquet approximation: single- and multi-boson exchange in the two-dimensional Hubbard model

Methods of quantum field theory represent a cornerstone of many-body physics. In their most general form they require the computation of multi-point correlation functions, whose dependence on several momentum and frequency labels lies beyond any practical implementation in most cases. An elegant formalism for the derivation of computationally feasible approximations for the electronic self-energy was introduced by Hedin Hedin65 , who expressed the latter in terms of the Green’s function ($G$), the screened interaction ($W$), and a vertex correction ($\gamma$). The simplest, so-called $GW$ approximation already includes the feedback of collective excitations on fermions and has become a standard tool of electronic structure theory, see, for example, Ref. Kutepov16 and references therein.

It is hard to go beyond the $GW$ approximation, although it is desirable in cases of gross quantitative discrepancies to experiment Kutepov16 or, for example, at strong coupling where vertex corrections may alter the interaction between fermions and bosons qualitatively Krien21 . However, it is not a trivial task to even define proper strategies to extend the $GW$ approximation: as usual, derivability from a potential leads to approximations that respect conservation laws Almbladh99 ; on the other hand one may also prefer, for example, a positive semi-definite real-axis spectrum as a stringent criterion Leeuwen14 ; Leeuwen16 ; or aim at including strong correlation effects Biermann03 ; Ayral15 . Here, on yet a different note, we are interested in an unbiased approach to the vertex correction $\gamma$ as it is provided, for example, by the parquet approach Diatlov57 ; Dominicis64-2 or by the functional renormalization group (fRG, Metzner12 ; Dupuis21 ), which respect the crossing symmetry of two-particle correlation functions.

Following this path, we recently introduced a variation of Hedin’s equations which is equivalent to the parquet approach Krien21-2 ; or, vice versa, one may say that the parquet approach was recast exactly into the $GW\gamma$ form. As such, it requires as an input the fully irreducible vertex $\Lambda$ of the parquet formalism, where fully irreducible implies that it can not be cut into two parts by removing two Green’s function lines ($GG$-irreducible, Rohringer12 ). The quantities that appear in Hedin’s equations are, however, irreducible with respect to the bare interaction ($U$-irreducible, Krien21-2 ). Therefore, the reformulated parquet equations actually use $\tilde{\Lambda}=\Lambda-U$ as a fundamental building block, where $U$ is the Hubbard interaction. In the application presented in this work we consider the parquet approximation, where $\tilde{\Lambda}$ vanishes, leading nevertheless to a highly nontrivial approximation for the Hubbard model.

To put the unification of Hedin’s formalism with the parquet approach into perspective, we recall that a key technique of quantum field theory is the boldification of Feynman diagrams: summarily denoting a partial series of diagrams by an effective quantity, as for example the self-energy, reduces the number of Feynman diagrams that need to be evaluated, at the expense of keeping track of the effective quantity which has to be computed self-consistently. In Hedin’s original approach further diagrams are summarized in the screened interaction $W$ and in the polarization, representing, respectively, a boson and a bosonic self-energy. The corresponding reduction in the number of Feynman diagrams is concisely put on display in Ref. Molinari06 , where it is also noted that keeping track of yet another quantity, the $U$-irreducible Hedin vertex $\gamma$ which mediates a Yukawa-like coupling between fermions and bosons, can be used to boldify diagrams even further. In this spirit, the key theoretical step taken in Ref. Krien21-2 is to boldify a subset of diagrams arising from the Bethe-Salpeter equations, which are of a simple structure. Namely, the $U$-reducible diagrams, coined single-boson exchange (SBE) in Ref. Krien19-2 , are representable in terms of the bold objects $\gamma$ and $W$ (see Fig. 1). The remaining $U$-irreducible diagrams, to which we refer as multi-boson exchange ($M$), do not permit a representation in terms of $\gamma,W$ alone, but instead capture repeated exchange of bosons. The resulting picture of bosons mediating effective interactions Bonetti22 is physically appealing and remains valid even at strong coupling Harkov21 ; Krien21 .

Further, it is plausible that fewer Feynman diagrams correspond to reduced computational cost in practical applications. Indeed, using the bosonized parquet approach, we are in a position to evaluate the parquet approximation for the Hubbard model on a $16\times 16$ lattice, which is, to our knowledge, the hitherto largest cluster size reached before any approximate parametrization of the momentum-dependent vertex functions as, for example, the truncated unity (TU) approximation Husemann09 ; Wang12 ; Lichtenstein17 ; Eckhardt20 . We note in passing that the performance may be also improved through a nonlocal formulation of the parquet approach Krien20 ; Astretsov19 . But we refrain from applying any further approximations (besides the parquet approximation itself). Therefore, the computational cost is reduced here only by the asymptotic decay of the vertex functions $M$ after the SBE diagrams are treated separately, because the latter determine the parquet vertices asymptotically Wentzell20 . As a result, frequency summations involving the vertex functions $M$ decay by one power faster compared to diagrams arising in the traditional parquet approach and, hence, the number of Matsubara frequencies can be reduced and the momentum grid refined. We thus arrive at the full-fledged parquet approximation for the Hubbard model, as envisioned in the seminal papers Diatlov57 ; Dominicis64-2 , progressing further along the path of pioneering applications to the Anderson impurity model Chen92 and small Hubbard clusters Yang09 ; Tam13 .

Recently, the parquet approach has also been unified with the multi-loop functional renormalization group (mfRG, Kugler18-3 ). By extension the latter can be recast in terms of boson exchange as well, a corresponding theory is presented in Ref. Walter22 . Further efforts aim at unbiased extensions of the dynamical mean-field theory (DMFT, Georges96 ) in order to reach the strong-coupling regime Toschi07 ; Taranto14 ; Rohringer18 ; Krien20 . Implementation details of the different methods vary widely and often additional approximations need to be applied. Therefore, we put here a spotlight on the plain parquet approximation as a (comparatively) simple reference case, which nevertheless provides a quantitative description of the Hubbard model in the weak-coupling limit Hille20 ; Schaefer21 . The aim of this paper is therefore twofold: On one hand, we discuss the qualitative behavior of various correlation functions, evaluated within the parquet approximation for the half-filled square lattice at weak coupling. On the other hand, with the full momentum dependence of the vertex functions readily available, we put two important tools to the test, namely, the TU approximation Eckhardt20 and the vertex asymptotics Wentzell20 . In the latter case our presentation extends to nonlocal correlations the investigation of Ref. Harkov21 , which compared the SBE diagrams to the vertex asymptotics for the Anderson impurity model.

The paper is structured as follows. We recollect definitions of the bosonized parquet formalism in Sec. 2. The screened interaction and Yukawa couplings are presented in Sec. 3, various four-point vertex functions are examined in Sec. 4. The convergence of the truncated unity is benchmarked in Sec. 5. We conclude in Sec. 6.

We consider the paramagnetic Hubbard model on the square lattice at half-filling,

$\displaystyle H=$

$\displaystyle-\sum_{\langle ij\rangle\sigma}{t}_{ij}c^{\dagger}_{i\sigma}c_{j\sigma}+U\sum_{i}n_{i\uparrow}n_{i\downarrow},$

(1)

where $t_{ij}$ denotes the hopping between nearest neighbors i and j, its absolute value $t=1$ sets the unit of energy. $c,c^{\dagger}$ are the annihilation and creation operators with the spin index $\sigma=\uparrow,\downarrow$. We denote the Hubbard repulsion between the densities $n_{\sigma}=c^{\dagger}_{\sigma}c_{\sigma}$ as $U$; we consider the weak-coupling regime, $2\leq U/t\leq 4$. The lattice size is fixed to $16\times 16$. The temperature is $T/t=0.2$.

We solve the Hubbard model (1) using the parquet approximation Diatlov57 ; Dominicis64-2 . In the following, we recollect only the most essential definitions. Readers with a background in parquet theory find a complete set of definitions, derivations, and the calculation cycle of our implementation in Ref. Krien21-2 . The notation used in this work is fully equivalent to Ref. Krien21-2 , it corresponds to a compromise between notations frequently used in the parquet and $GW$ literature. On the other hand, readers more familiar with the fRG find the corresponding definitions in Refs. Bonetti22 ; Walter22 , which use a notation more consistent with the fRG literature.

In the traditional parquet formalism the full vertex function is given in terms of the parquet decomposition,

$\displaystyle F=\Lambda+\Phi^{{ph}}+\Phi^{\overline{ph}}+\Phi^{{pp}}.$

(2)

Here, $\Lambda$ is the fully $GG$-irreducible vertex as explained in the introduction. The $\Phi$’s denote the vertices $GG$-reducible in the horizontal (${ph}$), vertical ($\overline{ph}$), and particle-particle (${pp}$) channel. Each vertex, e.g., $\Phi^{{ph},\alpha}(k,k^{\prime},q)$ carries a flavor label, in the particle-hole channel $\alpha=\text{ch}/\text{sp}$ corresponds charge or spin, and $k=(\mathbf{k},\nu)$, $q=(\mathbf{q},\omega)$ denote fermionic, bosonic momentum and Matsubara frequency, respectively. The parquet decomposition is shown at the bottom of Fig. 1.

On the other hand, Refs. Krien20 ; Krien21-2 introduced a bosonized parquet formalism where vertex diagrams are further decomposed, namely, the full vertex is expressed through the SBE decomposition Krien19-2 ,

$\displaystyle F=\Lambda^{\text{Uirr}}+\Delta^{{ph}}+\Delta^{\overline{ph}}+\Delta^{{pp}}-2U.$

(3)

The $\Delta$’s represent the $U$-reducible diagrams which can be cut in two parts by removing a bare interaction Krien21-2 ; Walter22 . They are given in terms of the Yukawa coupling (Hedin vertex) and the screened interaction, for example,

$\displaystyle\Delta^{{ph},\alpha}(k,k^{\prime},q)=\gamma^{\alpha}(k,q)W^{\alpha}(q)\gamma^{\alpha}(k^{\prime},q).$

(4)

The bare interaction arises as the leading order of all the $\Delta$’s, it is therefore subtracted twice in Eq. (3) to avoid overcounting. Notice that in Eq. (3) it also carries a flavor label, $U^{\text{ch}/\text{sp}}=\pm U$.

In turn, $\Lambda^{\text{Uirr}}$ is the fully $U$-irreducible vertex given through a parquet-like decomposition,

$\displaystyle\Lambda^{\text{Uirr}}=\tilde{\Lambda}+M^{{ph}}+M^{\overline{ph}}+M^{{pp}},$

(5)

where $\tilde{\Lambda}=\Lambda-U$ is the fully $GG$-irreducible vertex with the bare interaction removed. The $M$’s represent the multi-boson exchange, they are $GG$-reducible but fully $U$-irreducible vertices, whose momentum-energy dependence does not dissociate in the manner of Eq. (4).

Inserting Eq. (5) into Eq. (3), the resulting vertex decomposition of the bosonized parquet approach is shown in Fig. 1, above the horizontal line. It is convenient to add and subtract the bare interaction, represented by a black dot, so that the diagrams above the horizontal line are arranged consistently: summing the diagrams in a column yields the corresponding vertex of the traditional parquet formalism drawn below the horizontal line, for example,

$\displaystyle\Phi^{{ph},\alpha}=\Delta^{{ph},\alpha}-U^{\alpha}+M^{{ph},\alpha},\;\alpha=\text{ch},\text{sp},$

(6)

which connects the traditional and the bosonized parquet quantities on the left- and right-hand-side, respectively. Lastly, we introduce the parquet approximation:

$\displaystyle\tilde{\Lambda}\equiv 0.$

(7)

As a matter of fact, this is a rich approximation with nontrivial properties, as our results exemplify.

In this work we consider only the particle-hole quantities in Eqs. (4) and (6). In the traditional parquet formalism the set of equations is closed via the Bethe-Salpeter, Dyson, and Schwinger-Dyson equations. In the bosonized formalism the latter is replaced with the Hedin equation ($\Sigma=GW\gamma$) and one defines a bosonic self-energy ($\Pi=GG\gamma$), which determines the screened interaction via another Dyson-like equation ($W=U+U\Pi W$). However, to keep the presentation concise and general, we refer to Refs. Krien21-2 ; Bonetti22 ; Walter22 for detailed information including a calculation cycle or (m)fRG flow equations, respectively.

In our numerical application we use $N^{\gamma}_{\nu}=32$ fermionic and $N^{\gamma}_{\omega}=32$ bosonic Matsubara frequencies for the Yukawa couplings $\gamma$. The $M$’s are evaluated on a smaller fermionic frequency grid, using $N^{M}_{\nu}=16$ and $N^{M}_{\omega}=32$ for bosonic frequencies. Even though frequency summations like $\sum_{\nu^{\prime}}M(\nu,\nu^{\prime},\omega)G(\nu^{\prime})G(\nu^{\prime}+\omega)$ decay by one power of $\nu^{\prime}$ faster compared to a summation over the corresponding $\Phi$’s of the traditional parquet approach, a cutoff error arises in the $\gamma$’s for $\nu\approx(2N^{M}_{\nu}+1)\pi T$. Using smaller momentum grids we checked that our results for small $\nu$ presented in the following are not affected qualitatively by the frequency cutoff error. Quantitative convergence analysis for the $16\times 16$ grid is however beyond computational capability of the current implementation. As $\gamma$ determines the key observables $G$ and $W$, it is desirable to achieve convergence in $\gamma$ with respect to frequencies which would correspond to a very high standard of convergence for the parquet approach. An asymptotic treatment of $\gamma$ goes however beyond the established theory of vertex asymptotics Wentzell20 ; this problem may be considered elsewhere in the future.

Fermionic properties of the Hubbard model at weak coupling, in particular the formation of a pseudogap due to long-ranged spin fluctuations, have been discussed in great detail in the recent literature, see, for example, Refs. Hille20-2 ; Schaefer21 . However, electronic correlations renormalize also the Yukawa coupling between fermions and bosons, an effect which has received much less attention Krien20-2 . The parquet approach respects the crossing symmetry and hence provides us by construction with the full dependence of the Yukawa couplings on fermionic and bosonic momentum. Notice that we do not enter the pseudogap regime, which requires roughly 1000 lattice sites to avoid a finite-size effect Schaefer21 . However, we still observe an interesting evolution of $\gamma$ as antiferromagnetic fluctuations begin to build up.

We analyze the quantities $\Phi,\Delta,$ and $M$ in Eq. (6). These are four-point vertex functions depending on three momenta $\mathbf{k},\mathbf{k}^{\prime},\mathbf{q}$, and three frequencies $\nu,\nu^{\prime},\omega$. To get a grasp of these quantities, we focus on fermionic momenta $\mathbf{k}_{F},\mathbf{k}^{\prime}_{F}$ on the Fermi surface which traverse the path shown in Fig. 4, thereby passing through all four antinodal points. The fermionic frequencies are set to $\nu=\nu^{\prime}=\pi T$ or $\nu=-\nu^{\prime}=-\pi T$. We focus on the static limit $\omega=0$ and first set the bosonic transfer momentum to $\mathbf{q}=(\pi,\pi)$, which always guides scattered quasiparticles to final states on the Fermi surface.

In this manner we plot $M^{\text{sp}}(\mathbf{k}_{F},\mathbf{k}_{F}^{\prime},\mathbf{q}=(\pi,\pi),\nu=\pi T,\nu^{\prime}=\pi T,\omega=0)$ for $U/t=2$ in the top left panel of Fig. 5. Comparison with $\Delta^{\text{sp}}$ with the same labels, drawn in the center, shows that the latter exhibits a higher symmetry with respect to the fermionic momenta. Finally, $\Phi^{\text{sp}}$ on the right is obtained as the sum of $M^{\text{sp}}$ and $\Delta^{\text{sp}}$, with the bare interaction $U^{\text{sp}}=-U$ subtracted [cf. Eq. (6) and compare the magnitude of the color bars]. The high symmetry of $\Delta^{\text{sp}}$, which repeats along each of the four edges of the Fermi surface (cf. Fig. 4), implies that in a scattering event of two quasiparticles, mediated by this vertex, it is irrelevant to which of the four edges their initial momenta belong. In contrast, the lower symmetry of $M^{\text{sp}}$ implies that it mediates scattering events where it does matter whether the respective scattering partner lives on the same, an adjacent, or on the opposite edge of the Fermi surface.

Let us now consider the effect of flipping the sign of one fermionic frequency, $\nu\rightarrow-\pi T$. According to Eq. (8) in the previous section, $\gamma^{\text{sp}}(\omega=0)$ is symmetric with respect to $\nu$. Since the frequency dependence of the $\Delta$’s stems from the $\gamma$’s, $\Delta(\omega=0)$ is also invariant under the sign flip, which can be seen in the bottom center panel of Fig. 5. The situation is again quite different for $M^{\text{sp}}$ whose momentum structure is completely overturned under the sign flip of $\nu$. It was observed already in Refs. Krien19-2 ; Harkov21 that the fully $U$-irreducible vertex changes drastically when going from the sectors $\text{sgn}(\nu)=\text{sgn}(\nu^{\prime})$ to $\text{sgn}(\nu)=-\text{sgn}(\nu^{\prime})$. Apparently, in case of nonlocal correlations this is intertwined with its dependence on the fermionic momenta.

The patterns visible in $\Phi^{\text{sp}}$ arise from the superposition of those in $M^{\text{sp}}$ with the more symmetric ones in $\Delta^{\text{sp}}$, with an optically astounding result. Notice however that the color plot overemphasizes small variations in these quantities. It is $|\Delta^{\text{sp}}|\gg|M^{\text{sp}}|$, because the former inherits a large absolute value from $W^{\text{sp}}(\mathbf{q}=(\pi,\pi),\omega=0)$, and a weak $\mathbf{k}$ dependence from $\gamma^{\text{sp}}$ (cf. Fig. 2). We find that for larger interaction the difference in magnitude is even more enhanced and a discussion of the tiny variations is moot.

However, in the charge channel we find that $M^{\text{ch}}$ is larger than $\Delta^{\text{ch}}$ at small frequencies, see Fig. 6. The resulting $\Phi^{\text{ch}}$ is thus dominated by $M^{\text{ch}}$. Again $\Delta^{\text{ch}}$ is symmetric with respect to momenta and under a sign flip of $\nu$, whereas $M^{\text{ch}}$ not only changes its asymmetric momentum structure completely under the sign flip, but also its magnitude by a factor $4$ to $8$. Finally, we also present the charge quantities for an incommensurate bosonic momentum, $\mathbf{q}=(\pi/2,\pi/4)$, in Fig. 7. Although $\Delta^{\text{ch}}$ retains some regularity compared to $M^{\text{ch}}$, it loses much of its symmetry with respect to momenta, but remains symmetric under under a sign flip of $\nu$.

We applied the parquet approximation to the Hubbard model on a $16\times 16$ lattice and presented two-particle correlation functions corresponding to the bosonized parquet formalism introduced in Refs. Krien20 ; Krien21-2 .

The vertex functions reveal intriguing patterns as a function of the momenta, and the few shown examples scratch only the surface of the diverse variations that we observed in our calculations. It is an exciting outlook to consider the effects of next-nearest neighbor hopping, doping, larger interaction Krien20 ; Chalupa21 , and other modifications, where one or the other of the patterns may emerge as a physically important one.

We applied the truncated unity to quantities defined in the bosonized parquet formalism and benchmarked its convergence with the number of form factors. Similar to Ref. Harkov21 our analysis reveals that, in the considered setting, the formalism extends the asymptotic parametrization of the vertex functions Wentzell20 in a practically useful way to low frequencies. In particular, it facilitates fast convergence of the truncated unity approximation even in presence of long-ranged antiferromagnetic correlations.

Our implementation can be used to investigate properties of parquet-based approximations in their pure form for reasonably large lattice sizes, such as the fulfillment of Ward identities Janis17 ; Chalupa21-2 or nontrivial sum rules for the vertex functions Mermin67 , without any additional approximations.

Acknowledgments We acknowledge financial support from the Austrian Science Fund (FWF) through Projects No. P32044 and No. P30997.

A.K. implemented the algorithm of Ref. Krien21-2 . Both authors analyzed the results and jointly prepared the text.