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

Abstract The parquet approach to vertex corrections is unbiased but computationally demanding. Most applications are therefore restricted to small cluster sizes or rely on various simplifying approximations. We have recently shown that the bosonization of the parquet diagrams provides interpretative and algorithmic advantages over the original purely fermionic formulation. Here, we present first results of the numerical implementation of this method by applying it to the half-filled Hubbard model on the square lattice at weak coupling. The improved algorithmic performance allows us to evaluate the parquet approximation for a \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$16\times 16$$\end{document}16×16 lattice, retaining the full momentum and frequency structure of the various vertex functions. We discuss their symmetries and consider parametrizations of their momentum dependence using the truncated-unity approximation. Graphical abstract


Introduction
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 [1], who expressed the latter in terms of the Green's function (G), the screened interaction (W ), and a vertex correction (γ). 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. [2] and references therein.
It is hard to go beyond the GW approximation, although it is desirable in cases of gross quantitative discrepancies to experiment [2] or, for example, at strong coupling where vertex corrections may alter the interaction between fermions and bosons qualitatively [3]. 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 [4]; a positive semi- definite real-axis spectrum may also serve as a stringent criterion [5,6]; and/or one may want to include strong correlation effects [7,8]. Here, on yet a different note, we are interested in an unbiased approach to the vertex correction γ as it is provided, for example, by the parquet approach [9,10] or by the functional renormalization group (fRG, [11,12]), 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 [13]; or, vice versa, one may say that the parquet approach was recast exactly into the GW γ form. As such, it requires as an input the fully irreducible vertex Λ of the parquet formalism, where fully irreducible implies that it cannot be cut into two parts by removing two Green's function lines (GGirreducible, [14]). The quantities that appear in Hedin's equations are, however, irreducible with respect to the bare interaction (U -irreducible, [13]). Therefore, the reformulated parquet equations actually useΛ = Λ − U as a fundamental building block, where U is the Hubbard interaction. In the application presented in this work, we consider the parquet approximation, wherẽ Λ 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 [15] 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. [16], where it is also noted that keeping track of yet another quantity, the U -irreducible Hedin vertex γ 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. [13] 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. [17], are representable in terms of the bold objects γ and W (see Fig. 1). The remaining Uirreducible diagrams, to which we refer as multi-boson exchange (M ), do not permit a representation in terms of γ, W alone, but instead capture repeated exchange of bosons. The resulting picture of bosons mediating effective interactions [18] is physically appealing and remains valid even at strong coupling [3,19].
Furthermore, 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 × 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 [20][21][22][23]. We note in passing that the performance may be also improved through a nonlocal formulation of the parquet approach [24,25]. However, 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 [26]. 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 [9,10], progressing further along the path of pioneering applications to the Anderson impurity model [27] and small Hubbard clusters [28,29].
Recently, the parquet approach has also been unified with the multi-loop functional renormalization group (mfRG, [30][31][32]). By extension, the latter can be recast in terms of boson exchange as well, a corresponding theory is presented in Ref. [33]. Further efforts aim at unbiased extensions of the dynamical mean-field theory (DMFT, [34]) to reach the strongcoupling regime [24,[35][36][37]. 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 [38,39]. 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 [23] and the vertex asymptotics [26]. In the latter case, our presentation extends to nonlocal correlations the investigation of Ref. [19], 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 Sect. 2. The screened interaction and Yukawa couplings are presented in Sect. 3, various four-point vertex functions are examined in Sect. 4. The convergence of the truncated unity is benchmarked in Sect. 5. We conclude in Sect. 6.
where t ij denotes the hopping between nearest neighbors i and j, its absolute value t = 1 sets the unit of energy. c, c † are the annihilation and creation operators with the spin index σ =↑, ↓. We denote the Hubbard repulsion between the densities n σ = c † σ c σ as U ; we consider the weak-coupling regime, 2 ≤ U/t ≤ 4. The lattice size is fixed to 16 × 16. The temperature is We solve the Hubbard model (1) using the parquet approximation [9,10]. 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. [13]. The notation used in this work is fully equivalent to Ref. [13], 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. [18,33], 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 Here, Λ is the fully GG-irreducible vertex as explained in the introduction. The Φ's denote the vertices GGreducible in the horizontal (ph), vertical (ph), and particle-particle (pp) channel. Each vertex, e.g., Φ ph,α (k, k , q) carries a flavor label, in the particle-hole channel α = ch/sp corresponds to charge or spin, and k = (k, ν), q = (q, ω) denote fermionic, bosonic momentum, and Matsubara frequency, respectively. The parquet decomposition is shown at the bottom of Fig. 1.
On the other hand, Refs. [13,24] introduced a bosonized parquet formalism where vertex diagrams are further decomposed, namely, the full vertex is expressed through the SBE decomposition [17] The Δ's represent the U -reducible diagrams which can be cut in two parts by removing a bare interaction [13,33]. They are given in terms of the Yukawa coupling (Hedin vertex) and the screened interaction; for example The bare interaction arises as the leading order of all the Δ's, and it is therefore subtracted twice in Eq. (3) to avoid overcounting. Notice that in Eq. (3), it also carries a flavor label, U ch/sp = ±U .
In turn, Λ Uirr is the fully U -irreducible vertex given through a parquet-like decomposition whereΛ = Λ−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 which connects the traditional and the bosonized parquet quantities on the left-and right-hand side, respectively. Finally, we introduce the parquet approximatioñ 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 (Σ = GW γ) and one defines a bosonic self-energy (Π = GGγ), which determines the screened interaction via another Dyson-like equation (W = U + UΠW ). However, to keep the presentation concise and general, we refer to Refs. [13,18,33] for detailed information including a calculation cycle or (m)fRG flow equations, respectively.
In our numerical application, we use N γ ν = 32 fermionic and N γ ω = 32 bosonic Matsubara frequencies for the Yukawa couplings γ. The M 's are evaluated on a smaller fermionic frequency grid, using N M ν = 16 and N M ω = 32 for bosonic frequencies. Even though frequency summations like ν M (ν, ν , ω)G(ν )G(ν + ω) decay by one power of ν faster compared to a summation over the corresponding Φ's of the traditional parquet approach, a cut-off error arises in the γ's for ν ≈ (2N M ν + 1)πT . Using smaller momentum grids, we checked that our results for small ν presented in the following are not affected qualitatively by the frequency cut-off error. Quantitative convergence analysis for the 16 × 16 grid is however beyond computational capability of the current implementation. As γ determines the key observables G and W , it is desirable to achieve convergence in γ with respect to frequencies which would correspond to a very high standard of convergence for the parquet approach. An asymptotic treatment of γ (beyond the leading constant γ = 1 + · · · ) lies, however, outside the established theory of vertex asymptotics [26]; this problem may be considered elsewhere in the future.

Screened interaction and Yukawa coupling
Fermionic properties of the Hubbard model at weak coupling, in particular the formation of a pseudogap in the spectral function due to long-ranged spin fluctuations, have been discussed in great detail in the recent literature; see, for example, Refs. [39,40]. However, electronic correlations renormalize also the Yukawa coupling between fermions and bosons, an effect which has received much less attention [41]. 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 [39]. Nevertheless, we still observe an interesting evolution of γ as antiferromagnetic fluctuations begin to build up.

Screened interaction
The left panels of Fig. 2 show the static screened interactions W ch/sp (q, ω = 0) along the high-symmetry path. For comparison, we normalize it with the absolute value U of the bare interaction. The sign of the different curves, therefore, signals repulsion (W/U > 0) or attraction (W/U < 0) and the amplitude indicates whether the interaction U ch/sp = ±U in the respective channel is screened (|W/U | < 1) or enhanced (|W/U | > 1). As expected, with increasing U , a strong attractive interaction develops in the spin channel along the q = (π, π) direction.
Overall, γ ch/sp depend most strongly on q, less strongly on ν, and the least on k (dependence on ω will be considered elsewhere). However, this can not be generalized, as γ ch (k, ν = πT, q, ω = 0) shows a sizable k-dependence for U/t = 3, whereas γ sp is largely independent of k for the same set of parameters and labels. In the non-interacting system, the Yukawa coupling is unity; Fig. 2 shows that a weak interaction leads to screening (γ ch/sp < 1). Notice that γ determines both the fermionic (Σ = GW γ), as well as the bosonic self-energy (Π = GGγ), which also enters Σ via W . Close to an instability an increase of γ, even by a few percent, can drastically enhance W . Indeed, we showed recently that even for the harmless parameters U/t = 2, T/t = 0.2, the screening of γ sp is indispensable to obtain a reasonable approximation for Σ [13].
Furthermore, as the system is driven to the antiferromagnetic instability, fermions decouple from the soft bosons (γ sp → 0 as W sp → −∞), since the Goldstone excitations of the ordered phase are pro- Fig. 3 Numerical validation of Eq. (9) for the static charge Yukawa coupling γ ch (k, q, ν, ω = 0): shifting k by −q is the same as going from k to −k tected (Adler principle, [42,43]). Indeed, comparing U/t = 2 and U/t = 3 in Fig. 2, we see that γ sp is much more strongly screened around q = (π, π) for the larger interaction, which corresponds to a longer correlation length (see also Sect. 5). On the other hand, we found in recent investigations that, as soon as fermionic states are destroyed due to the feedback from the spin fluctuations, this requirement is lifted and γ sp rises again for those k where a pseudogap opens, resulting in a nodal/antinodal dichotomy of γ sp with respect to k [3,41]. There hence exists a subtle interplay between bosonic fluctuations, Fermi surface features, and the Yukawa couplings, which needs to be considered when dependencies of the latter are neglected or parametrized.

Symmetries
We also discuss symmetries of the Yukawa couplings, see Refs. [44,45]. First, we note that inversion symmetry of the lattice, as well as time-reversal and SU(2) symmetry are required for the derivation in Ref. [13] and by our implementation. This set of symmetries allows to interchange the fermionic labels of the full vertex function F (k, k , q) = F (k , k, q); see also Refs. [14,45]. Since the Yukawa coupling is just a four-point vertex with tapered Green's function legs on one side (plus 1) [13], the symmetry of the full vertex implies that it does not matter on which side the legs are attached. As a result, the left-and right-handed Yukawa couplings shown in Fig. 1 are identical. It is important to keep in mind, however, that in a more general setting, our formalism needs to be re-derived using left-and righthanded Yukawa couplings [18,33].
A symmetry valid by definition is due to complex conjugation, γ * (k, q) = γ(−k, −q). On the other hand, the γ's are invariant under symmetry operations of the point group of the lattice [46]. For example, inversion symmetry implies γ(k, q, ν, ω) = γ(−k, −q, ν, ω). Since the symmetry operations needs to be applied to both momenta at the same time, in a practical implementation only one of the momenta can be mapped to the irreducible wedge of the lattice. Hence, for the 16 × 16 square lattice, each Yukawa coupling requires ∼ 45 256 (256) 2 N γ ν N γ ω complex numbers. Inversion combined with complex conjugation further implies γ(k, q, ν, ω) = γ * (k, q, −ν, −ω), and hence Re γ(k, q, ν, ω = 0) =Re γ(k, q, −ν, ω = 0). (8) In the considered weak-coupling regime, we find that the imaginary part Im γ is negligibly small. Moreover, it vanishes for momenta k and k + q on the Fermi surface [e.g., for any k on the Fermi surface and q = (π, π)]. For these reasons, we can ignore the imaginary part in the following.
Finally, we verify numerically that a nontrivial symmetry of the γ's is respected by our implementation. Namely, the full four-point vertex satisfies by definition the "swapping symmetry" F kk q = F k +q,k+q,−q [47]. Together with F kk q = F k kq , it follows [48] that γ ch/sp (k−q, q) = γ ch/sp (k, −q) 1 . We set q = (q, ω = 0), resulting in In the last line, we applied the inversion symmetry. Equation (9) implies for γ ch/sp (k, q, ν, ω = 0) that shifting k → k−q has the same effect as k → −k. That this is indeed the case in our implementation can be seen in Fig. 3 which shows Re γ ch for, e.g., q = (π/2, π/4) and q = (3π/4, 0). We chose here γ ch for U/t = 3, as it depends strongly on k (see Fig. 2), and incommensurate q for a generic result. Symmetries put strong conditions on the γ's which are useful to verify code during debugging, or to save memory space.

Single-and multi-boson exchange
We analyze the quantities Φ, Δ, and M in Eq. (6). These are four-point vertex functions depending on three momenta k, k , q, and three frequencies ν, ν , ω.
To get a grasp of these quantities, we focus on fermionic momenta k F , k 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 ν = ν = πT or ν = −ν = −πT . We focus on the static limit ω = 0 and first set the bosonic transfer momentum to q = (π, π), which always guides scattered quasiparticles to final states on the Fermi surface. In this manner, we plot the real part of M sp (k F , k F , q =  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 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, ν → −πT . According to Eq. (8) in the previous section and for vanishing Im γ sp , γ sp (ω = 0) is symmetric with respect to ν. Since the frequency dependence of the Δ's stems from the γ's, Δ(ω = 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 sp whose momentum structure is completely overturned under the sign flip of ν. It was observed already in Refs. [17,19] that the fully U -irreducible vertex changes drastically when going from the sectors sgn(ν) = sgn(ν ) to sgn(ν) = −sgn(ν ). Apparently, in case of nonlocal correlations, this is intertwined with its dependence on the fermionic momenta.
The patterns visible in Φ sp arise from the superposition of those in M sp with the more symmetric ones in Δ sp , with an optically astounding result. Notice, however, that the color plot overemphasizes small variations in these quantities. It is |Δ sp − U sp | |M sp |, because the former inherits a large absolute value from W sp (q = (π, π), ω = 0), and a weak k dependence from γ 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 ch | is comparable to |Δ ch − U ch | at small frequencies and depends more strongly on k F and k F ; see Fig. 6. The resulting patterns in Φ ch are therefore dominated by M ch . Again, Δ ch is symmetric with respect to momenta  and under a sign flip of ν, whereas M 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, q = (π/2, π/4), in Fig. 7. Although Δ ch retains some regularity compared to M ch , it loses much of its symmetry with respect to momenta, but remains symmetric under under a sign flip of ν.

Convergence of the truncated unity
While in this work, we kept the full momentum dependence of the various vertex functions, this is in general undesirable beyond applications to simple model systems. It is therefore, on one hand, a question of practical interest to parametrize the momentum dependencies in a memory-efficient way. On the other hand, the formal construction of the theory should also work toward this goal. Here, for example, the single-boson exchange Δ is by construction parametrized through W and γ. However, if a simplified parquet or fRG scheme keeps also the multi-boson exchange M , the question arises whether the bosonized theory offers any advantages over a traditional fermionic formulation using the Φ's. Moreover, the vertex asymptotics [26] is often used to parametrize the Φ's at high frequencies. Since the vertex asymptote corresponds itself to high-frequency limits of the Δ's [17,19], the bosonized theory may only offer advantages in the low-frequency regime.
In this regard, Ref. [19] recently demonstrated that the Δ's capture resonant low-frequency features of the local full vertex function F loc of the Anderson impurity model (AIM). Even though other low-frequency features reside in the M 's, the two-particle quantities of the DMFT approximation are recovered to good accuracy using only the Δ's (cf. Fig. 1; F loc was approximated by neglecting all of the M 's). If however lowfrequency information in the Δ's is also neglected, the parametrization of F loc fails at strong coupling [19]. Concretely, we find the difference between Δ sp and its asymptotic expression as follows: and Ref. [19] showed for the AIM (k, k , q → ν, ν , ω) that an approximation for F loc should keep the term Δ sp R , which vanishes asymptotically for |ν| → ∞ and/or |ν | → ∞.
Here, we draw an analogy to the present investigation: while the effective AIM of the DMFT approximation exhibits strong local spin fluctuations at strong coupling, here the Hubbard model at weak coupling  Fig. 6 for incommensurate bosonic momentum q = (π/2, π/4) develops long-ranged spin-density wave fluctuations. Physically, these scenarios are of course quite different; for example, in the AIM γ sp seems to diverge for small ν and low temperature [19], while Fig. 2 shows that in the Hubbard model, 0 < γ sp < 1 is screened. However, a similarity is that the screened interaction W sp is large, either due to the local moment in the AIM, or, here, because of the growing antiferromagnetic correlation length ξ. In the latter case, it is therefore plausible that the term Δ sp R in Eq. (10) grows with ξ, and at the same time also develops a strong dependence on the bosonic momentum q. In this case, it could be advantageous to keep Δ sp R parametrized as a part of Δ sp , rather than to assign it to a memory-intensive four-point vertex. This is what we show in the following.
To this end, we expand the q-dependence (see remarks, Sect. 5.2) of various vertices in the form-factor basis [23] and observe the convergence with respect to the number of expansion coefficients; see also Ref. [13] where this was done for U/t = 2 and T /t = 0.2. To keep the maximal number of form factors f ( , q) small, we use results for an 8 × 8 lattice. We transform, for example, Φ sp to the form-factor basis and back into q-space, after discarding all but N form factors (11) where we set ν = ν = πT, ω = 0, k = k = ( π 2 , π 2 ) fixed. The complete q dependence is thus recovered for N = 64, but the series may be truncated at a smaller N if the expanded vertex is sufficiently short-ranged in real space (truncated unity). Blue lines in Fig. 8 show for q = (π, π) the thus expanded Φ sp , the reducible vertex of the traditional parquet formalism, for U/t = 2, 3, 4. Notice that in the considered regime, the antiferromagnetic correlation length ξ increases monotonously with U . Namely, we find for U/t = 2 and 3 that ξ ≈ 1.5 and 2, respectively, which are consistent with our calculations for the 16 × 16 lattice. For U/t = 4, we expect a sizable finite-size effect for the 8 × 8 lattice [49], which arises for ξ on the order of half the linear lattice size or larger.
Since the form-factor expansion of Φ sp with respect to q converges only slowly, Ref. [24] introduced the idea, within the bosonized parquet approach, to expand only the multi-boson exchange M sp in form factors, while the full momentum dependence of Δ sp was kept. Using Eq. (6), this corresponds to the approximation Φ sp (q) ≈ M sp (q, N ) + Δ sp (q) − U sp . Red lines in Fig. 8 show this result again for q = (π, π). For U/t = 2, 3, 4, this approximation lies close to the fully converged Φ sp even for N = 1. This indicates, remarkably, that the relative importance of M sp compared to Δ sp does not increase with ξ at all (even if the correlations described by M sp grow in range as ξ increases, they do not grow faster than it is the case for Δ sp ).
On the other hand, we show now that the relative importance of the term Δ sp R compared to Δ sp does increase with the correlation length. To this end, we expand this term together with M sp , such an approximation may be written as This corresponds to a parametrization of Φ sp where its highfrequency limits are given through the vertex asymptote, Δ sp − Δ sp R , retaining full momentum dependence, while the rest function M sp + Δ sp R is expanded in form factors. The convergence of this parametrization can be observed in the green lines drawn in Fig. 8. As expected, the convergence with form factors worsens considerably as the correlation length increases at larger U/t; in fact, for U/t = 4, it becomes comparable to the slow convergence of Φ sp . We conclude that it is advantageous to keep Δ R parametrized through Δ, rather than to combine it with M sp in a rest function.

Remarks
Several remarks are in order to put the result reported in Fig. 8 into perspective. First, we recall that the truncated unity is intended foremost to parametrize the dependence on fermionic momenta k, k , which is often much weaker than the q dependence. However, an unbiased approach to two-particle correlations, such as parquet or fRG schemes, requires channel projections which map the q dependence in one channel to the k, k dependence in another. It was therefore noted in Ref. [23] that the truncated-unity cut-off unfortunately also appears in bosonic arguments. This explains the fast convergence of the truncated unity in Refs. [24,41], where it was only applied to the M 's. In this respect, it is also encouraging that the relative importance of Fig. 8 Convergence of the truncated unity applied in three different ways: blue lines show the direct application to Φ sp , cf. Ref. [23]; red lines the application only to M sp , cf. Refs. [24,41]. Green lines indicate application to M sp + Δ sp R , the rest function of the vertex asymptote [26] M sp compared to Δ sp appears to be independent of the correlation length (Fig. 8), so that the quality of a fixed truncated-unity cut-off N does not deteriorate with growing ξ. Hence, while a purely fermionic mfRG scheme using the truncated unity approximation is quantitatively accurate for small ξ [38], the fast convergence of the truncated unity in our bosonized formalism becomes crucial for larger ξ [24]. We expect this benefit to translate directly to the corresponding bosonized (m)fRG schemes of Refs. [18,33].
Compared to the traditional parquet formalism, the improved performance of our implementation, and the generally weaker momentum dependence of the quantities calculated in it, are reminiscent of similar observations in the context of vertex-corrected GW approaches [2,50].
On the other hand, one has to keep in mind that the practical advantage of the bosonized formalism depends on the physical regime and the correlation functions of interest. For example, we find in the half-filled Hubbard model at weak coupling that |Δ sp − U sp | is much larger than |M sp |; however, in the charge channel, we find the opposite in the low-frequency regime. In particular in applications to pseudogaps induced by spin-density wave fluctuations the charge sector is of a lesser interest; however, it remains to be seen how much improvement the bosonized formalism offers in other physical settings. One may hope that in a regime which exhibits strong charge fluctuations, the importance of Δ ch may be enhanced over M ch .
However, a case where a breakdown of the fast convergence of the truncated unity can be expected is, for example, a regime of long-ranged d-wave singlet fluctuations. They are captured by the corresponding M s of the particle-particle channel [18]. How much the results suffer from this may depend on the importance of the feedback of the d-wave fluctuations on other channels, which requires a projection operation, as discussed above. In this regard, it is intriguing to consider a re-bosonization and suitable parametrization (through new γ's and W 's) of the corresponding strongly fluctuating channel captured by the M 's. As the example of the d-wave shows, the bosonized formalism does not come with an autopilot for improved performance. However, in any case, the interpretative advantages of the bosonization remain, and there are, to our knowledge, no disadvantages associated with it.

Conclusions
We applied the parquet approximation to the Hubbard model on a 16 × 16 lattice and presented two-particle correlation functions corresponding to the bosonized parquet formalism introduced in Refs. [13,24].
The vertex functions reveal intriguing patterns as functions 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 [24,51], 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. [19] our analysis reveals that, in the considered setting, the formalism extends the asymptotic parametrization of the vertex functions [26] in a practically useful way to low frequencies. In particular, it facilitates fast convergence of the truncated unity approximation even in the presence of long-ranged antiferromagnetic correlations.
Our implementation can be used to investigate the properties of parquet-based approximations in their pure form for reasonably large lattice sizes, such as the fulfillment of Ward identities [52,53] or nontrivial sum rules for the vertex functions [54], without any additional approximations.