Truncated-Unity Parquet Equations: Application to the Repulsive Hubbard Model
C. J. Eckhardt G. A. H. Schober J. Ehrlich C. Honerkamp


->-/.style=decoration=
  markings,
  mark=at position #1 with Latex[scale=1.6],postaction=decorate






















[]christian.eckhardt@rwth-aachen.de






^1 Institute for Theoretical Solid State Physics, RWTH Aachen University, 52074 Aachen,
Germany 

 ^2 Peter Grünberg Institute and Institute for Advanced Simulation, Forschungszentrum 52435 Jülich, Germany 

^3 JARA-FIT, Jülich Aachen Research Alliance - Fundamentals of Future Information Technology, Germany 















The parquet equations are a self-consistent set of equations for the effective two-particle vertex of an interacting many-fermion system. 
The application of these equations to bulk models is, however, demanding due to the complex emergent momentum and frequency structure of the vertex.
Here, we show how a channel-decomposition by means of truncated unities, which was developed in the context of the functional renormalization group to efficiently treat the momentum dependence, can be transferred to the parquet equations.
This leads to a significantly reduced complexity and memory consumption scaling only linearly with the number of discrete momenta.
We apply this technique to the half-filled repulsive Hubbard model on the square lattice and present approximate solutions for the channel-projected vertices and the full reducible vertex.








Truncated-Unity Parquet Equations: 
 Application to the Repulsive Hubbard Model
    C. Honerkamp^1,3
    December 30, 2023
===============================================================================







§ INTRODUCTION

The parquet approach was introduced more than five decades ago as a method to analyze interacting many-fermion systems <cit.>.
Since then, it has helped significantly in understanding the physics of magnetic impurities in metals <cit.> as well as the breakdown of Fermi liquid behavior in one-dimensional metals <cit.>.
For a recent application to Hubbard nanorings, see Ref. valli.

A main advantage of the parquet scheme is that it can be made self-consistent at the single-particle and two-particle level <cit.>.
It has been shown that the parquet approximation (to be introduced later) is a thermodynamically consistent and conserving, Φ-derivable approximation in Baym's sense at the one-particle level <cit.>.
From early on, it has been clear that the parquet approximation is closely related to perturbative renormalization group (RG) schemes.
At least regarding the modern fermionic functional renormalization group (fRG) flavors (for reviews, see e.g. Refs. metzner2012,platt13), it has been understood that the parquet approximation and the fRG in the usual truncations sum the same classes of diagrams <cit.>.
However, it has also been known that in the fRG, due to the unavoidable truncation of the hierarchy of flow equations, certain combinations of internal lines are suppressed compared to the contributions kept in the parquet approximation.
Parts of these missing contributions can be recollected by refined flow equations <cit.>.
Quite recently, a systematic multi-loop fRG scheme was proposed <cit.> that in simplified models reconciles the fRG results with those of the parquet approximation in a quantitative manner.
These attempts to lift the fRG on higher levels already indicate that fRG approaches may offer advantages despite the fact that they do not readily contain the full perturbative corrections of the parquet approximation.
Indeed, considering the vast more recent literature on standard zero- to two-dimensional correlated many-fermion lattice systems, the applications of RG schemes seem to outnumber clearly those of parquet schemes.

In order to make the comparison more specific, we mention two recent state-of-the-art parquet studies of the two-dimensional Hubbard model <cit.>. 
In these works, the finest momentum resolution reaches 6×6 due to memory constraints.
Hence, very few points are located in the vicinity of the Fermi surfaces, and many interesting questions like the generation of unconventional superconductivity, the opening of a pseudogap, or the tendency toward incommensurate or stripe ordering are hard to study. By contrast, in the fRG approaches, 𝒪(100) momentum-space patches were employed at an early stage without the use of parallel computers <cit.>, and 14×14 grids were also analyzed <cit.>.
It should be mentioned that these works did not consider the frequency dependence of the interactions, which is usually kept in the parquet studies. Yet, more recent fRG schemes are about to remedy this shortcoming while still reaching similar momentum resolutions.

A substantial progress in simplifying the description of the momentum structure has been provided by channel-decomposed fRG schemes <cit.>.
The main simplification there consists in expressing the effective two-particle interaction, which depends on three momenta (usually two incoming and one outgoing momentum, the fourth one being fixed by momentum conservation), by three interaction functions that describe the interaction between specific fermion bilinears in particle-particle and particle-hole channels. Each of these interaction functions depends strongly on one `bosonic' momentum but only weakly on two other wavevectors, hence the latter dependencies can be expanded in a suitable form-factor basis. In the simplest cases, this means that the interacting fermion bilinears live on nearby sites on the lattice, while longer bilinears are usually not relevant.
This then allows for a well-convergent and physically meaningful truncation of the form-factor expansion. Without frequency dependence, the channel-decomposed fRG schemes can be parallelized efficiently and pushed to very fine momentum resolutions with thousands of momenta in the Brillouin zone, in conjunction with convergence checks in the form-factor truncation <cit.>.
Currently, the inclusion of the frequency dependence and self-energy effects into these schemes is on the way (see e.g. Refs. giering,eberlein15,vilardi).
For the frequency dependence, a channel decomposition was also shown to yield meaningful results for impurity models <cit.>, but in general, the complex frequency structure requires a more sophisticated description <cit.>.







The main goal of the present article is to show how the advantages of the channel decomposition explored in fRG studies can be transferred to the parquet equations.
We show that without the frequency dependence, this readily gives meaningful results with a high momentum resolution. In particular, the channel decomposition of the two-particle interaction reduces the memory required for numerically evaluating the parquet equations from 𝒪(N^3) to 𝒪(N), where N is the number of momenta in the first Brillouin zone.
This makes is plausible that the channel decomposition will also be beneficial in cases where frequency-dependent interactions are considered.

The article is organized as follows. We begin in Sct. <ref> by briefly describing the parquet equations. After that, we introduce the projections onto the direct particle-hole, crossed particle-hole and particle-particle channel, and subsequently derive the truncated-unity (TU) parquet equations. In Sct. <ref>, we further derive the cross projections between the different channels, which are necessary for iteratively solving the TU parquet equations. In Sct. <ref>, we discuss the advantages of the channel decomposition with respect to computational complexity and memory cost. Following this, we provide more details of our numerical implementation in Sct. <ref>, and we present our results for the half-filled repulsive Hubbard model in Sct. <ref>. Finally, the appendix is concerned with the derivation of channel-decomposed parquet equations for general spin-SU(2)-symmetric systems.



§ CHANNEL DECOMPOSITION
 

We start from the parquet equations as described, for example, in Ref. held. 
These are formulated in terms of the full Green function G, the full (one-particle-irreducible) two-particle-reducible vertex F, 
the (in channel r two-particle-) reducible vertices ^r in the direct particle-hole channel (r = ph,d),
the crossed particle-hole channel (r = ph,c), and the particle-particle channel (r = pp),
as well as the corresponding (in channel r two-particle-) irreducible vertices ^r in each channel.
Explicitly, the parquet equations read as follows:[We mainly follow the conventions of Refs. Schober, Schober17, held; however, the two-point Green function used here contains an additional factor β = 1/(k_B T), and all four-point functions contain an additional factor 2 compared to those used in Refs. Schober, Schober17. Further note that the functions ^r on the left-hand side of the parquet equations are not identical to the corresponding functions used in the fRG framework.]


    _σ_1 σ_2 σ_3 σ_4^ph,d(p_1, p_2, p_3)
         = -k_B T/N∑_k, σ, σ' F_σ_1 σ' σσ_4(p_1, k + p_2 - p_3, k)   G(k)   G(k + p_2 - p_3)  ^ph,d_σσ_2 σ_3 σ'(k, p_2, p_3)  , 
    _σ_1 σ_2 σ_3 σ_4^ph,c(p_1, p_2, p_3)
         = k_B T/N∑_k, σ, σ' F_σ_1 σ' σ_3 σ(p_1, k + p_3 - p_1, p_3)   G(k)   G(k + p_3 - p_1)  ^ph,c_σσ_2 σ' σ_4(k, p_2, k + p_3 - p_1)  , 
    _σ_1 σ_2 σ_3 σ_4^pp(p_1, p_2, p_3)
         = -1/2k_B T/N∑_k, σ, σ' F_σ_1 σ_2 σσ'(p_1, p_2, k)   G(k)   G(p_1 + p_2 - k)  ^pp_σσ' σ_3 σ_4(k, p_1 + p_2 - k, p_3)  .


Here, the function arguments k ≡ (ω, k⃗) ≡ (k_0, k⃗) are multi-indices comprising fermionic Matsubara frequencies (at temperature T) and Bloch momenta (of which there are N in the first Brillouin zone).
Each vertex in Eqs. (<ref>)–(<ref>) depends on only three momenta/frequencies, since the fourth is always determined by momentum/energy conservation.
Furthermore, each vertex depends on four spin indices, while we have assumed the full Green function to be spin independent corresponding to the SU(2)-symmetric case. The different vertices appearing in Eqs. (<ref>)–(<ref>) are related through

    F =  + ^ph,d + ^ph,c + ^pp ,

where  denotes the fully (two-particle-) irreducible vertex,
and by

    ^r = F - ^r .

Throughout this article, we will use the parquet approxi­mation <cit.>, which identifies the fully irreducible vertex with the initial interaction given in the Hubbard model by Eq. (<ref>), hence
Λ≡ F^0.
Graphically, the parquet equations can be represented by means of Feynman diagrams as in Fig. <ref> (see also Ref. held).

Next, we define the bosonic Matsubara frequencies and corresponding 
transfer (t⃗, u⃗) or total (s⃗) momenta as

    t ≡ (t_0,t⃗)     = (ω_3-ω_2, p⃗_3-p⃗_2)  , 
    
    u ≡(u_0,u⃗)     = (ω_1-ω_3,p⃗_1-p⃗_3)  , 
    
    s ≡(s_0,s⃗)     = (ω_1+ω_2,p⃗_1+p⃗_2)  .

We introduce the following vertices, which differ from the original vertices merely in a relabeling of their arguments (for other conventions, see Refs. husemann, julian, Wang12):

    _k_1 k_2(t)     := (k_1,  k_2 - t,  k_2)  , 
    _k_1 k_2(u)     := (k_1,  k_2 - u,  k_1 - u)  , 
    _k_1 k_2(s)     := (k_1,  s - k_1,  s - k_2)  .

Note that the three functions on the left-hand side are actually different, although in our notation they are distinguished only by their respective argument (t, u, or s). With these definitions, the parquet equations can be rewritten compactly as follows (for deriving the last equation, one has to use that ^pp is antisymmetric with respect to its first two arguments, see the appendix:


    [^ph,d_σ_1 σ_2 σ_3 σ_4]_k_1 k_2(t)     = -k_B T/N∑_k, σ, σ'[F_σ_1 σ' σσ_4]_k_1 k(t)   G(k)   G(k - t)  [^ph,d_σσ_2 σ_3 σ']_k k_2(t)  , 
    [^ph,c_σ_1 σ_2 σ_3 σ_4]_k_1 k_2(u)     = k_B T/N∑_k, σ, σ'[F_σ_1 σ' σ_3 σ]_k_1 k(u)   G(k)   G(k - u)  [^ph,c_σσ_2 σ' σ_4]_k  k_2(u)  , 
    [^pp_σ_1 σ_2 σ_3 σ_4]_k_1 k_2(s)     = 1/2k_B T/N∑_k, σ, σ'[F_σ_1 σ_2 σσ']_k_1 k(s)   G(k)   G(s- k)  [^pp_σ'σσ_3σ_4]_k  k_2(s)  .



We generally expect that the three reducible vertices ^ph,d, ^ph,c, ^pp,
which together solve the parquet equations, have a strong dependence on their respective main transfer or total momentum (see e.g. Refs. husemann, husemann12, eberlein10, maier2013, wang, julian).
On the other hand, the dependencies on the remaining two fermionic momenta are expected to be rather weak. 
Therefore, we assume that these dependencies can be described to a sufficient accuracy by using only few smooth basis functions. 

Thus, we consider a set of form factors, i.e., functions in the first Brillouin zone (BZ) denoted by

    { f_ℓ⃗(k⃗); k⃗∈ 1^st BZ, ℓ⃗∈ℤ×ℤ} ,

which are assumed to be pairwise orthonormal and complete in the sense that

    1/N∑_k⃗f_ℓ⃗(k⃗)  f_ℓ⃗'^*(k⃗)     = δ_ℓ⃗ℓ⃗' , 
    ∑_ℓ⃗f_ℓ⃗(k⃗)  f_ℓ⃗^*(k⃗')     = N δ_k⃗k⃗'.

Later, the index ℓ⃗ will label the sites of the two-dimensional square lattice (see Eq. (<ref>)). For technical reasons, we further introduce the multi-indices ℓ = (ℓ_0, ℓ⃗) and the “frequency-dependent” form factors

    f_ℓ(k) := δ_ℓ_0,  k_0  f_ℓ⃗(k⃗)  ,

which are orthonormal and complete in the sense that

    1/N∑_k f_ℓ(k)  f_ℓ'^*(k)     = δ_ℓℓ' , 
    ∑_ℓ f_ℓ(k)  f_ℓ^*(k')     = N δ_k k' .

In fact, the following considerations would remain valid for more general, truly frequency-dependent form factors fulfilling Eqs. (<ref>)–(<ref>).
Finally, we define the projections of an arbitrary vertex  onto the direct particle-hole, crossed particle-hole and particle-particle channel as follows (cf. Ref. [Eqs. (18)-(20)]julian):

    D̂[]_ℓ_1 ℓ_2(t)     = 1/N^2∑_k_1,  k_2_k_1 k_2(t)   f_ℓ_1(k_1)   f_ℓ_2^*(k_2)  , 
    Ĉ[]_ℓ_1 ℓ_2(u)     = 1/N^2∑_k_1,  k_2_k_1 k_2(u)   f_ℓ_1(k_1)   f_ℓ_2^*(k_2)  , 
    P̂[]_ℓ_1 ℓ_2(s)     = 1/N^2∑_k_1,  k_2_k_1 k_2(s)   f_ℓ_1(k_1)   f_ℓ_2^*(k_2)  .

By combining these equations with Eqs. (<ref>)–(<ref>), we obtain the following formulae, by which the vertex  can be reconstructed from its respective projections (these equations will be used in Sct. <ref> to derive the “cross projections” between different channels):


    (p_1, p_2, p_3)     = ∑_ℓ_1, ℓ_2D̂[]_ℓ_1 ℓ_2(p_3 - p_2)   f_ℓ_1^*(p_1)   f_ℓ_2(p_3)  , 
    (p_1, p_2, p_3)     = ∑_ℓ_1, ℓ_2Ĉ[]_ℓ_1 ℓ_2(p_1 - p_3)   f_ℓ_1^*(p_1)   f_ℓ_2(p_1 + p_2 - p_3)  , 
    (p_1, p_2, p_3)     = ∑_ℓ_1, ℓ_2P̂[]_ℓ_1 ℓ_2(p_1 + p_2)   f_ℓ_1^*(p_1)   f_ℓ_2(p_1 + p_2 - p_3)  .


Using these relations, 
we will now derive parquet-type equations for the reducible vertices 
^ph,d, ^ph,c, ^pp, where each of these vertices is
projected onto its respective main momentum. 
For this purpose, we apply the mappings (<ref>)–(<ref>) to both sides of Eqs. (<ref>)–(<ref>)
and insert two partitions of unity of the form-factor basis (see Eq. (<ref>)) on both sides of the fermion loops.
This procedure is analogous to the derivation of flow equations
in the truncated-unity functional renormalization group (TUfRG) scheme <cit.>.
We thus arrive at the following self-consistent equations, which we call the truncated-unity (TU) parquet equations:


    D̂[^ph,d_σ_1 σ_2 σ_3 σ_4]_ℓ_1 ℓ_2(t)    =
    - ∑_ℓ, ℓ', σ, σ'D̂[F_σ_1 σ' σσ_4]_ℓ_1ℓ(t)   L^ph_ℓℓ'(t)  D̂[^ph,d_σσ_2 σ_3 σ']_ℓ' ℓ_2 (t)  ,
    
    Ĉ[^ph,c_σ_1 σ_2 σ_3 σ_4]_ℓ_1 ℓ_2(u)    =
    ∑_ℓ, ℓ', σ, σ'Ĉ[F_σ_1 σ' σ_3 σ]_ℓ_1ℓ(u)   L^ph_ℓℓ'(u)  Ĉ[^ph,c_σσ_2 σ' σ_4]_ℓ' ℓ_2 (u)  ,
    
    P̂[^pp_σ_1 σ_2 σ_3 σ_4]_ℓ_1 ℓ_2(s)    = 1/2 ∑_ℓ, ℓ', σ, σ'P̂[F_σ_1 σ_2 σσ']_ℓ_1 ℓ(s)   L^pp_ℓℓ'(s)  P̂[^pp_σ'σσ_3 σ_4]_ℓ' ℓ_2(s)  ,

where the particle-hole loop L^ph and the particle-particle loop L^pp are given by

    L^ph_ℓℓ'(t)     = k_B T/N∑_k G(k)   G(k - t)   f_ℓ(k)   f_ℓ'^*(k)  , 
    
     L^pp_ℓℓ'(s)     = k_B T/N∑_k G (k)   G(s - k)   f_ℓ(k)   f_ℓ'^*(k)  .


One main feature of Eqs. (<ref>)–(<ref>) is that they involve only matrix multiplications with respect to the internal summation indices ℓ and ℓ' of the form-factor basis (provided one uses a countable set of form factors).
This structure of the TU parquet equations is particularly advantageous for the numerical parallelization because
it allows for an independent evaluation for different values of t, u, and s, and thus for a distribution of the vertices over several compute nodes. We remark, however, that internode communication is still needed when invoking Eq. (<ref>). Furthermore, the above form of the parquet equations still requires the calculation of cross projections between different channels as we will explain in the following.



§ CROSS PROJECTIONS
 

The standard procedure for solving the self-consistent parquet equations is an iteration scheme, which takes the bare interaction as the initial vertex and in each step evaluates the parquet equations once to recalculate the vertices. To explain this in more detail for the TU parquet equations (<ref>)–(<ref>), let us assume that in one iteration step we have calculated the projected vertices
D̂[^ ph,d], Ĉ[^ph,c], and P̂[^pp].
Then, in the next step, we want to employ Eqs. (<ref>)–(<ref>) to recalculate these vertices. Consider, as an example, the projection D̂[F] of the total vertex, which appears on the right-hand side of Eq. (<ref>). This can be split into the initial interaction and the three channels (see Eq. (<ref>)), i.e.,

    D̂[F]_ℓ_1 ℓ(t)     = 
    D̂[F^0]_ℓ_1 ℓ + D̂[^ph,d]_ℓ_1 ℓ(t) 
           +
    D̂[^ph,c]_ℓ_1 ℓ(t) +
    D̂[^pp]_ℓ_1 ℓ(t) .

Similarly, D̂[^ph,d] can be calculated from Eq. (<ref>) as

    D̂[^ph,d]_ℓ' ℓ_2(t) = D̂[F]_ℓ' ℓ_2(t) - D̂[^ph,d]_ℓ' ℓ_2(t)  .

Now, the first term in Eq. (<ref>), i.e., the projection of the initial interaction, is known explicitly: combining Eqs. (<ref>) and (<ref>) gives

    D̂[F^0]_ℓ_1 ℓ    = 1/N^2∑_k,  k' f_ℓ_1(k)  f_ℓ_2^*(k')   U 
       ≡⟨ f_ℓ_1⟩⟨ f_ℓ^* ⟩  U  ,

where we have omitted the spin dependencies. The second term in Eq. (<ref>), namely
D̂[^ph,d], is directly available from the previous iteration step. By contrast, for calculating the projections
D̂[^ph,c] and D̂[^pp] from the previously obtained Ĉ[^ph,c] and P̂[^pp], it is necessary to invert the projections Ĉ and P̂. Hence, we can calculate D̂[F] by means of the formal identity


    D̂[F]_ℓ_1 ℓ(t)     = 
    D̂[F^0]_ℓ_1 ℓ + D̂[^ph,d]_ℓ_1 ℓ(t) 
           + D̂[Ĉ^-1[Ĉ[^ph,c]]]_ℓ_1 ℓ(t) 
           + D̂[P̂^-1[P̂[^pp]]]_ℓ_1 ℓ(t)  .



To calculate the “cross projections” between the different channels, we may use the definitions (<ref>)–(<ref>) as well as Eqs. (<ref>)–(<ref>). After some algebra, we thus obtain


    D̂[F]_ℓ_1 ℓ(t)     =
    D̂[F^0]_ℓ_1 ℓ + D̂[^ph,d]_ℓ_1 ℓ(t) 
           + 1/N^2∑_k'∑_ℓ_3, ℓ_4Ĉ[^ph,c]_ℓ_3 ℓ_4(k')  ∑_k_1 f_ℓ_3^*(k_1)   f_ℓ_4(k_1 - t)   f_ℓ_1(k_1)   f_ℓ^*(k_1 - k') 
           + 1/N^2∑_k'∑_ℓ_3, ℓ_4P̂[^pp]_ℓ_3 ℓ_4(k')  ∑_k_1 f_ℓ_3^*(k_1)   f_ℓ_4(k_1 - t)   f_ℓ_1(k_1)   f_ℓ^*(k' - k_1 + t)  .


Similarly, one can derive the cross projections between any two other channels, and this then allows one to iteratively solve the TU parquet equations.

We go on to describe some approximations which further simplify the iterative solution. First, we neglect the frequency dependencies of all vertices. The remaining frequency sums in the fermion loops (<ref>)–(<ref>) can be performed analytically <cit.>, giving 


    L^ph_ℓ⃗ℓ⃗'(t⃗)     := ∑_ℓ_0, ℓ'_0 L^ph_ℓℓ'(t⃗, t_0 = 0) 
     = 1/N∑_k⃗n_F(ε(k⃗)) - n_F(ε(k⃗ - t⃗))/ε(k⃗) - ε(k⃗ - t⃗)  f_ℓ⃗(k⃗)   f_ℓ⃗'^*(k⃗)  , 
    
    L^pp_ℓ⃗ℓ⃗'(s⃗)     := ∑_ℓ_0, ℓ'_0 L^pp_ℓℓ'(s⃗, s_0 = 0)  = 1/N∑_k⃗1 - n_F(ε(k⃗)) - n_F(ε(s⃗ - k⃗))/ε(k⃗) + ε(s⃗ - k⃗)  f_ℓ⃗(k⃗)   f_ℓ⃗'^*(k⃗)  .


Here, n_F(ε) = (1 + exp(βε))^-1 denotes the Fermi distribution function, which depends on the inverse temperature β = 1/k_B T.
Next, as has already been seen in (Ref. wang; see also Sct. 3.3 of Ref. Schober17), a further simplification can be achieved by using complex exponentials as form factors, i.e.,

    f_ℓ⃗(k⃗) =e^-i a ℓ⃗·k⃗ ,

where a denotes the lattice spacing. This set of functions naturally fullfills the requirements (<ref>)–(<ref>). Then, Eq. (<ref>) simplifies to

    D̂[F^0]_ℓ_1 ℓ = U  δ_ℓ⃗_1 , 0⃗ δ_ℓ⃗, 0⃗ .

Furthermore, by introducing Fourier-transformed projections D, C, P, such that

    D̂[F]_ℓ⃗_1 ℓ⃗(t⃗) = ∑_ℓ⃗'e^-i a ℓ⃗' ·t⃗D[F]_ℓ⃗_1 ℓ⃗(ℓ⃗')  ,

or conversely,

    D[F]_ℓ⃗_1 ℓ⃗(ℓ⃗') = 1/N∑_t⃗e^i a ℓ⃗' ·t⃗D̂[F]_ℓ⃗_1 ℓ⃗(t⃗)  ,

we can transform Eq. (<ref>) into

    D̂[F]_ℓ⃗_1 ℓ⃗(t⃗) = U  δ_ℓ⃗_1 , 0⃗ δ_ℓ⃗, 0⃗ + D̂[^ph,d]_ℓ⃗_1 ℓ⃗(t⃗) 
        + ∑_ℓ⃗'e^i a ℓ⃗' ·t⃗ C[^ph,c]_ℓ⃗' + ℓ⃗_1 - ℓ⃗,  ℓ⃗'(-ℓ⃗) 
        + ∑_ℓ⃗'e^i a ℓ⃗' ·t⃗ P[^pp]_ℓ⃗' + ℓ⃗_1,  ℓ⃗' - ℓ⃗(ℓ⃗)  .

We remark that instead of simple plane-wave functions given by Eq. (<ref>), one could also use form factors that explicitly respect the symmetry of the lattice, so-called lattice harmonics.
Fewer of these are required to reach the same accuracy as with simple plane-wave functions. For an introduction to lattice harmonics, we refer the interested reader to Sct. A3 of Ref. platt13.

Finally, the main approximation which the channel projetions aim at is to keep only a finite number of form factors f_ℓ⃗(k⃗), such that

    ℓ⃗∈ [-ℓ_cut, ℓ_cut]×[-ℓ_cut, ℓ_cut]  ,

with a cutoff parameter ℓ_cut∈ℕ.
The total number of form factors is then n_cut≡ 4 ℓ_cut^2.



To illustrate the implications of this approximation, we consider again Eq. (<ref>).
For each particular combination of arguments ℓ_1, ℓ, and t, the projected vertex D̂[F]_ℓ_1 ℓ(t) depends on Ĉ[^ph,c]_ℓ_3 ℓ_4(k') (and similarly on P̂[^pp]_ℓ_3,ℓ_4(k')) at all possible arguments ℓ_3, ℓ_4, and k'. Similar considerations hold for the Fourier-transformed vertices as in Eq. (<ref>).
In other words, it is in principle necessary to know the projected vertices at all possible arguments in one iteration step
before one can recalculate these vertices at any particular argument in the next iteration step.
Therefore, when keeping only a limited set of basis functions as specified by Eq. (<ref>),
the vertices Ĉ[^ph,c]_ℓ_3 ℓ_4 and P̂[^pp]_ℓ_3 ℓ_4 (or their Fourier transforms)
from one iteration step contribute to D̂[F]_ℓ_1 ℓ in the next iteration step only approximately (in the sense that also their strong dependence on the main momentum is not accounted for exactly).
This in turn reflects the fact that the projections D̂, Ĉ, and P̂ are actually not invertible if restricted to a limited set of form factors.
Nevertheless, the inverse mappings Ĉ^-1 and P̂^-1 as formally employed in Eq. (<ref>) can be approximated by means of Eq. (<ref>) or Eq. (<ref>), provided that a sufficiently large number of form factors is taken into account in these summations.
 
To summarize, it is not a priori clear why keeping only a limited number of form factors is a good approximation.
In particular, if the vertex has some sharp momentum structures in one channel (close to a phase transition at low temperatures, or caused by long-range initial interactions), then these sharp structures in one channel would in principle also lead to sharp structures in the other channels by means of the cross projections.
However, keeping only a few form factors would cause these sharp structures to be smeared out.
Notwithstanding this caveat, we expect that local contributions to the vertices are more important for determining critical temperatures.
Previous fRG studies (e.g. Ref. julian) have shown that for the t-t' Hubbard Model, including the Fermi-surface studied here
[The plot in the cited paper does not explicitly cover half filling.
The authors of the present manuscript have, however, convinced themselves that the same conclusion also holds for half-filling.
This has been done together with the authors of Ref. julian using the same code as in the cited work.]
, the inclusion of higher-order form factors does not change the critical scales considerably.




§ COMPLEXITY AND MEMORY COST
 

In this section, we briefly discuss the computational complexity and memory cost of our projection scheme
as compared to a direct implementation of the parquet equations.
We limit this discussion to the evaluation of the momentum dependencies of the vertices, since the other dependencies
(i.e., those on spin and frequency) are not affected by the projections (see Eqs. (<ref>) and (<ref>)–(<ref>)). 

First, consider a direct evaluation of the parquet equations (<ref>)–(<ref>).
For each combination of three “external” momentum arguments, a sum over one “internal” momentum has to be performed. Since each momentum ranges in the first BZ, the complexity of evaluating these equations is 𝒪(N^4),
where N is the number of discrete Bloch momenta.
This can possibly be improved by using matrix multiplications, for which efficient algorithms exist,
but a lower bound is always 𝒪(N^3) corresponding to the memory cost.
Furthermore, since vertices of three momentum arguments have to be stored, the memory consumption scales like 𝒪(N^3).

Next, consider the channel-decomposed parquet equations (<ref>)-(<ref>). For each argument of the vertices on the left-hand side, two sums over the form-factor basis have to be performed. Since the vertices have N × n_cut^2 arguments, this scales like 𝒪(N × n_cut^4). In order to evaluate the parquet equations, one further needs the projections of F onto each channel as well as the loop functions L^ph and L^pp.
Note that only the cross projection of F has to be computed, because the one of ^r can simply be obtained via Eq. (<ref>). For evaluating the cross projections we have presented two different schemes, namely Eq. (<ref>) and Eq. (<ref>). In the first case, one has to perform one momentum sum and two form-factor sums  for each argument of the projected vertex, which scales like 𝒪(N^2× n_cut^4).
The second sum in Eq. (<ref>) (over k_1) can be computed in advance and is therefore faster than the other calculations.
In the second case, i.e., when using Eq. (<ref>), only one sum over the form-factor basis has to be performed.
Again, this has to be done for each argument of the projected vertex, thus scaling like 𝒪(N × n_cut^3). Finally, one still needs to calculate the bubble functions via Eqs. (<ref>)–(<ref>). Here, one momentum sum is required for each argument of the bubbles, thus scaling like 𝒪(N^2× n_cut^2).

We conclude that when using Eq. (<ref>), the cross projections are actually the most expensive calculations such that the overall complexity is 𝒪(N^2× n_cut^4). On the other hand, when invoking Eq. (<ref>), the most expensive calculation is the evaluation of the bubble functions.
However, these bubble functions need to be calculated only once because they do not change during the iteration process (at least if no self-energy is calculated). Hence, in this second case the amortized costs are determined by the evaluation of Eqs. (<ref>)-(<ref>), and therefore the overall amortized costs scale like 𝒪(N × n_cut^4). In particular, the complexity then scales only linearly with the momentum resolution.














Concerning the memory costs, all vertices have to be stored as functions of their main momentum and two basis-function indices, which implies a scaling like 𝒪(N × n_cut^2).
This does not change when employing Eq. (<ref>): for this scheme, the Fourier-transformed vertices would have to be stored in addition to the original vertices. The Fourier-transformed vertices are not larger than the original vertices, but if the vertices are already highly memory-consuming, the additional storage space required might still not be negligible.

Thus, if memory is the main constraint, one can instead employ Eq. (<ref>) at the cost of a higher complexity.
In summary, if only a constant set of basis functions is kept, the memory consumption scales only linearly with the momentum resolution N.

Finally, we remark that previous implementations of the parquet equations <cit.> have shown that the “bottle neck” for evaluating them is in fact the memory consumption. In particular, when performing computations on heavily parallel machines, the internode communications usually limit the system sizes that can be treated. Thus, reducing the memory consumption to linear scaling in N as in the present TU parquet approach may indeed represent a major step forward.



§ NUMERICAL IMPLEMENTATION
 

In order to check the validity of our method, we have applied the TU parquet equations to the Hubbard model on the square lattice given by the Hamiltonian

    H = -t ∑_⟨ i,  j ⟩, σ a_i, σ^† a_j, σ + U ∑_i n_i, ↑ n_i, ↓ .

Here, a_i, σ and a_i, σ^† denote the annihilation and creation operators at site i with spin σ∈{↑, ↓}, and
n_i, σ = a_i, σ^† a_i, σ the number operator at site i. 
The sum in the first term in Eq. (<ref>) is only over nearest-neighbor sites, while the second term describes an onsite interaction. By diagonalizing the Hamiltonian, one obtains the energy dispersion

    ε(k⃗) = 2  t  (cos(k_x  a) + cos(k_y  a) )  .

The Hubbard model has been studied extensively in the literature (for some references, see e.g. Refs. maierDCA,gull,schaefer) as it is expected to be relevant for the high-temperature superconducting cuprates <cit.>. Therefore, this model can serve as a testing case for any newly developed quantum many-body method.

Our numerical solution of the TU parquet equations uses the approximations mentioned in Sct. <ref>, hence we neglect all frequency dependencies of the vertices, setting t_0 = u_0 = s_0 = 0.
As form factors we choose complex exponentials, which have already been given in Eq. (<ref>).
Furthermore, we neglect the self-energy and thereby replace the full Green function with the bare Green function: G = G_0.
Finally, as our aim here is a consistency check, we limit the discussion to the  contributions (thus considering only a single basis function) in the TU parquet equations.
It is known, however, that the critical scales in TUfRG for the same model at half filling change only insignificantly when higher form factors ℓ⃗≠0⃗ are taken into account <cit.>.
Of course, this changes when one moves away from half filling and non-local pairing becomes important. In any case, within the “onsite” ℓ⃗ = 0⃗ approximation, the evaluation of the parquet equations becomes particularly simple because then Eq. (<ref>) reduces to

    D̂[F]_0⃗0⃗(t⃗)     = U + D̂[^ph,d]_0⃗0⃗(t⃗) 
           + ⟨Ĉ[^ph,c]_0⃗0⃗⟩ + ⟨P̂[^pp]_0⃗0⃗⟩ .

Here, ⟨ . ⟩ denotes the mean of a vertex with respect to its main momentum argument (which coincides with the Fourier-transformed vertex evaluated at zero), i.e.,

    ⟨⟩ := 1/N∑_t⃗(t⃗)  .

Thus, the projection operation Ĉ[D̂^-1[D̂[  .  ]] simply reduces to taking the mean of D̂[  .  ] with respect to its main momentum argument, and the same applies to the other cross projections.


Before writing out the parquet equations within the above approximations, let us briefly discuss the spin dependence of the vertices.
In fact, since the Hubbard model has an SU(2)-symmetry (see e.g. Ref. SH00), it is possible to eliminate all spin indices and thereby to simplify the parquet equations.
We only state the results here and refer the interested reader to the appendix for a detailed derivation.
Generally, SU(2)-symmetry implies the following spin dependence of the vertices as a consequence of the so-called crossing relations <cit.>:


    ^ph,d_σ_1 σ_2 σ_3 σ_4(k_1,k_2,k_3,k_4)     = -V^ph,d(k_1,k_2,k_3,k_4)  δ_σ_1 σ_4δ_σ_2 σ_3 + V^ph,c(k_1,k_2,k_4,k_3)  δ_σ_1 σ_3δ_σ_2 σ_4 , 
    ^ph,c_σ_1 σ_2 σ_3 σ_4(k_1,k_2,k_3,k_4)     = -V^ph,c(k_1,k_2,k_3,k_4)  δ_σ_1 σ_4δ_σ_2 σ_3 + V^ph,d(k_1,k_2,k_4,k_3)  δ_σ_1 σ_3δ_σ_2 σ_4 , 
    ^pp_σ_1 σ_2 σ_3 σ_4(k_1,k_2,k_3,k_4)     = -V^pp(k_1,k_2,k_3,k_4)  δ_σ_1 σ_4δ_σ_2 σ_3 + V^pp(k_1,k_2,k_4,k_3)  δ_σ_1 σ_3δ_σ_2 σ_4 .


Conversely, the spin-independent V-functions can be obtained from the spin-dependent vertices by evaluating the latter at particular spin combinations, i.e.,
    V^ph,d(k_1, k_2, k_3, k_4) = -^ph,d_↑↓↓↑(k_1, k_2, k_3, k_4)  ,

and similarly for the other vertices. 

Furthermore, by evaluating the parquet equations (<ref>)–(<ref>) at spin arguments (σ_1 σ_2 σ_3 σ_4) = (↑↓↓↑) and using Eqs. (<ref>)–(<ref>), one can derive the corresponding parquet equations for the V-functions. These in turn can be transformed into a channel-decomposed version as shown in the appendix.


In summary, we have implemented approximate TU parquet equations for the 0⃗0⃗-components of the channel-projected V-functions, which we abbreviate as

    D(t⃗)     := D̂[V^ph, d]_0⃗0⃗(t⃗)  , 
    
     C(u⃗)     := Ĉ[V^ph, c]_0⃗0⃗(u⃗)  , 
    
     P(s⃗)     := P̂[V^pp]_0⃗0⃗(s⃗)  .












The equations which we have implemented read as follows (where averages are defined as in Eq. (<ref>), and where the 0⃗0⃗-components of the fermion loops can be read off from Eqs. (<ref>)–(<ref>)):


    D(t⃗)     = 2  (  U + D(t⃗) + ⟨ C ⟩ + ⟨ P ⟩)   L^ph_0⃗0⃗(t⃗)  (  U + ⟨ C ⟩ + ⟨ P ⟩) 
           - (  U + D(t⃗) + ⟨ C ⟩ + ⟨ P ⟩)   L^ph_0⃗0⃗(t⃗)  (  U + ⟨ D ⟩ + ⟨ P ⟩) 
           - (  U + ⟨ D ⟩ + C(t⃗) + ⟨ P ⟩)   L^ph_0⃗0⃗(t⃗)  (  U + ⟨ C ⟩ + ⟨ P ⟩)  , 
     
    
    C(u⃗)     = -  (  U + ⟨ D ⟩ + C(u⃗) + ⟨ P ⟩)   L^ph_0⃗0⃗(u⃗)  (  U + ⟨ D ⟩ + ⟨ P ⟩)  , 
    
    P(s⃗)     = -  (  U + ⟨ D ⟩ + ⟨ C ⟩ + P(s⃗) )   L^pp_0⃗0⃗(s⃗)  (  U + ⟨ D ⟩ + ⟨ C ⟩)  .




















Importantly, since the fermion loops as well as all mean values of the vertices can be computed in advance, the numerical evaluation of these equations scales only linearly with the momentum resolution.



§ RESULTS FOR THE HUBBARD MODEL
 











Finally, we present our numerical results for the Hubbard model at half filling and without next-nearest-neighbor hopping.
In this case, the non-interacting dispersion has the well-known fully nested Fermi surface, and the ground state at temperature T=0 and positive onsite interaction U should exhibit antiferromagnetic (AF) long-range order.
At nonzero temperatures, the Mermin-Wagner theorem prohibits long-range order, but as in the better-understood Heisenberg model there should still be longer-ranged AF correlations <cit.>.
A proper description of the state at T>0 with quantum many-body methods requires special care <cit.>, and many methods like DCA or fRG replace the low-T short-range ordered state with a long-range ordered state up to an artificial Néel temperature T_c.
The major goal of the present analysis is to show that our channel-decomposed parquet scheme can reproduce the approach to long-range order with leading AF correlations, with a smaller remnant T_c. 

The main observable that we study here is the full one-particle-irreducible vertex F.
Its momentum structure near the AF instability is well known from fRG studies and discussed e.g. in Sct. IIIB of Ref. metzner2012.
Very close to the instability, which occurs at some nonzero T_c in the usual fRG approximations, the vertex has the leading momentum dependence 

    V_crit (p⃗_1,p⃗_2,p⃗_3) = J/4 ( 2  δ_p⃗_1- p⃗_3, Q⃗ + δ_p⃗_3- p⃗_2, Q⃗)   ,

with J ∝ 1/|T-T_c|. Here, p⃗_1 and p⃗_2 are incoming momenta, and Q⃗ = (π,π). 
This effective interaction can be transformed onto the real lattice, leading to an infinitely-long-ranged 
AF spin-spin interaction

    J ∑_⟨ i,  j ⟩ e^i Q⃗· (R⃗_i-R⃗_j)S⃗_i ·S⃗_j  ,

with spin operators defined as

    S⃗_i = 1/2∑_α, βσ⃗_αβ  c^†_i,α c_i,β .

The same expression is also found in the random phase approximation (RPA) by using bare Green functions when ladder and bubble diagram chains are summed up (see Ref. scalapino93 Eq. (5) close to the divergence).
The fRG changes the RPA results by reducing T_c and adding more non-divergent structure to the vertex.








Figures <ref>–<ref> show the numerical results for the vertices D(t⃗), C(u⃗), P(s⃗) defined by Eqs. (<ref>)–(<ref>) as functions of their respective main momentum, which ranges in the first Brillouin zone. For these computations an initial interaction of U=2.0  t was used.
The temperature was set to a small value of T = 0.1  t, and the computations were performed on a 200 × 200 grid.  We further mention that in order to improve the convergence, we have transformed the recursive formulas according to <cit.>

    x_i+1 = f(x_i)   →   x_i+1 = α f(x_i) + (1-α)  x_i

with a constant α∈ (0,1], which in our implementation was set to 0.5. The evaluation of Eqs. (<ref>)–(<ref>) takes only a few minutes on a laptop.

One can see that the crossed particle-hole vertex C(u⃗) is always positive and strongly peaked at u⃗= (π,π).
Similarly, the direct particle-hole vertex D(t⃗) is peaked at t⃗=(π,π),
but it has also negative contributions and its absolute value stays below that of C(u⃗). The two strong peak features can be clearly associated with the two terms of Eq. (<ref>), where the t⃗ = (π,π) feature is represented by the first term and the 
u⃗= (π,π) peak by the second term.
By contrast, P(s⃗) peaks at s⃗=(0,0) and takes only negative values. This is also clearly understood and expected from the sign and momentum dependence of the particle-particle loop diagram. 
For completeness, we also show the full vertex

    V^F(p⃗_1, p⃗_2, p⃗_3) =    U + V^ph,d(p⃗_3 - p⃗_2)+
       V^ph,c(p⃗_1 - p⃗_3) + V^pp(p⃗_1 + p⃗_2)

in Fig. <ref>.
By setting p⃗_3 = (π/2, π/2), one clearly sees two peaked lines at p⃗_1 = (-π/2, -π/2) and at p⃗_2 = (-π/2, -π/2), which correspond to t⃗=(π, π) and u⃗=(-π,-π), respectively.
Furthermore, one observes a dip at p⃗_⃗1⃗=-p⃗_⃗2⃗, which results from the V^pp contribution.

When further decreasing the temperature, we find a divergence of C and D at (π,π) as shown in Fig. <ref>. The behavior close to the divergence is well approximated by

    C(T) ∼1/|T-T_c| ,

where T_c denotes the Néel critical temperature. In our implementation we find T_c≈ 0.035  t
(at this value, the vertex is more than three times larger than the bandwidth).
In addition, as shown in Fig. <ref>, we observe the relation C /D → 2 at (π,π) and close to the divergence at T_c, which is precisely the ratio between the two terms in Eq. (<ref>) thus reproducing the results from Ref. metzner2012.
This also matches the behavior of previous RPA studies close to the divergence as seen for instance in Ref. scalapino93. The parquet vertex in the TU approximation shows the same behavior of a nesting-driven AF ordering instability expected from RPA. 

The TU parquet results for the Hubbard model also agree with fRG treatments of the same model as described e.g. in Ref. metzner2012.
In fact, the full parquet vertex constructed from the three channels can be directly compared with fRG data for the same model parameters, as shown e.g. in Fig. 10 of Ref. metzner2012. There, F has been obtained from an N-patch fRG scheme in the `standard' level-two truncation <cit.>, also without self-energy corrections as in our work here. The comparison of our results with the fRG data obtained with the same code as in Ref. metzner2012 can be seen in Fig. <ref>. The results for the full vertex match qualitatively, in terms of the enhancement features and also where the vertex remains small. The qualitative agreement with these fRG works supports the basic validity of our method. It also confirms the sensibility of the zeroth-order truncation in the form factor expansion for this case, as the analysis in Ref. metzner2012 does not rely on form factors. The parquet data also agree with a recent preprint, Ref. tagliavini, that for the same parameters uses a form-factor expansion and at least partially includes non-local form factors. Hence, our method passes this qualitative sanity check. 


The observable differences in Fig. <ref> concern the resolution of the enhancement features. The parquet features are much narrower than the corresponding ones in fRG. One reason for this is that the fRG divergence scale T^fRG_c∼ 0.11 t for the approximation used is about three times higher than in parquet, and we had to choose T=0.4 t in order to find a similar magnitude of the largest couplings than in the parquet data shown here for T=0.1  t. Furthermore, the momentum resolution of the parquet result is given by the 200 × 200 grid points in the whole Brillouin zone, while the fRG works with just 96 points on the Fermi surface.  

The difference in the divergence scales should be discussed in more detail. The first point to keep in mind is that the true T_c for the AF ordering instability should be zero because of the Mermin-Wagner theorem. The self-consistent parquet approximation with the bare interaction as fully irreducible vertex, but including self-energy feedback on the internal lines, was argued<cit.> to correctly describe the finite-T correlations of Heisenberg-like order parameters, i.e., to fulfill the Mermin-Wagner theorem. In comparison with this, we remark that we do not include self-energies here because these would require a refined treatment of the frequency dependencies within our scheme (cf. Refs. held_new, held2). The same approximation holds for the N-patch data shown in Fig. <ref>. The parquet treatment includes more perturbative corrections than the fRG in the level-2 truncation, and indeed, the T_c found here is smaller, ∼ 0.035  t instead of T^fRG_c∼ 0.11t, as visible in Fig. <ref>.
This nourishes hope that a parquet approach with self-energy feedback could actually get close to fulfillment of the Mermin-Wagner constraints.
Furthermore, in the recent preprint tagliavini, the multiloop-fRG scheme including frequency dependence of the interaction and self-energy feedback was applied to the same situation in the Hubbard model. The multiloop corrections reconstruct parquet contributions that are missed in the level-2 truncation of the fRG and result also in a reduced divergence scale compared to the previous fRG results.
The upshot of this comparison is that the quantitative picture at low T depends on the further approximations used. The divergence scale becomes indeed smaller if a better approximation is used and should - theoretically, which may be hard in practice - reach zero if all parquet approximation terms and self-energies are included, and if sufficient momentum and frequency resolution is obtained. 
 Here, our numerically efficient parquet scheme may be a good starting point for further refinements.  











§ CONCLUSION AND OUTLOOK


We have presented a channel decomposition of the parquet equations, which effectively reduces the number of momentum arguments of the two-particle-reducible vertex functions. This method relies on introducing resolutions of unity in a form-factor basis, of which only a finite set of basis functions is kept. In this sense, the scheme can be called truncated-unity (TU) parquet approximation in analogy to the recently developed truncated-unity functional renormalization group (TUfRG) method <cit.>. In the TUfRG, the convergence in the number of form factors kept has turned out to be rather quick in most parameter regimes <cit.>.

In a numerical implementation of the parquet equations, the channel decomposition effectively reduces the memory consumption from 𝒪(N^3) to 𝒪(N), where N denotes the number of Bloch momenta taken into account. Furthermore, this method is particularly suitable for the parallelization on a large number of compute nodes. Since memory consumption is generally regarded as the “bottle neck” for implementing the parquet equations, the channel decomposition may allow one to reach a much higher precision in predicting ground-state orderings and critical scales in many-body models of condensed matter physics. 

To benchmark our method, we have implemented the TU parquet equations disregarding self-energy and frequency dependencies, and restricting attention to the lowest form factors. By means of the channel decomposition, we could study momentum resolutions with 𝒪(10^4) momenta in the first Brillouin zone with only a few minutes computing time on a standard laptop.
 Our results for the two-particle-reducible vertices in the fully nested Hubbard model on the square lattice qualitatively match those of previous fRG studies of the same model and get closer to a fulfillment of the Mermin-Wagner theorem than published fermionic fRG works. A more sophisticated implementation of the TU parquet equations which takes into account more form factors and/or the frequency dependencies of the vertices may improve this issue even more and is currently underway.













































































We thank K. Held, C. Hille, A. Kauch, F. Kugler, J. Lichtenstein, T. Reckling, M. Salmhofer, D. Sánchez de la Peña, A. Tagliavini, D. Vilardi, and J. von Delft for helpful discussions.
We further acknowledge the computing time granted on the RWTH Compute Cluster. This work was supported by the DFG grants HO 2422/10-1, 11-1, and 12-1, DFG RTG 1995, and by the Vienna Computational Materials Laboratory (FWF SFB ViCoM F41).





§ SU(2)-SYMMETRIC PARQUET EQUATIONS
 

In this appendix, we show how spin-SU(2)-symmetry can be used to facilitate the evaluation of the parquet equations.
In the SU(2)-symmetric case, all vertices effectively depend on only three spin arguments, the fourth being determined by spin conservation.
Correspondingly, in the parquet equations given in Ref. held,

the sum over internal spin indices can be limited to only a few spin configurations. 
Often, one also considers superpositions of such configurations to evaluate the parquet equations <cit.>. In this article, we use instead the following decomposition of vertices, which generally holds in the SU(2)-symmetric case (see e.g. Ref. SH00):

    _σ_1 σ_2σ_3σ_4(k_1,k_2,k_3,k_4) = -V^(k_1,k_2,k_3,k_4)  δ_σ_1 σ_4δ_σ_2 σ_3 + W^(k_1,k_2,k_4,k_3)  δ_σ_1 σ_3δ_σ_2 σ_4 .

For a general , the component functions V^ and W^ can be independent of each other. However, if we consider the full vertex, ≡ F, we can further employ its antisymmetry under the simultaneous interchange of two momentum and spin arguments, i.e., the so-called crossing relation <cit.>

    F_σ_1 σ_2σ_3σ_4(k_1,k_2,k_3,k_4) = -F_σ_1 σ_2σ_4σ_3(k_1,k_2,k_4,k_3)  .

From this equation, one can derive the following relation between the coefficient functions in Eq. (<ref>):

    V^F(k_1, k_2, k_3, k_4) = -F_↑↓↓↑(k_1, k_2, k_3, k_4) = F_↑↓↑↓(k_1, k_2, k_4, k_3) = W^F(k_1, k_2, k_3, k_4)  ,

hence V^F ≡ W^F. Similarly, the particle-particle vertex _pp is also antisymmetric under the exchange of two momentum and spin arguments, hence Eqs. (<ref>) and (<ref>) hold analogously for the particle-particle vertex. Another similar case is the bare interaction of the Hubbard model, which is given by

    F^0_σ_1 σ_2 σ_3 σ_4(k_1, k_2, k_3, k_4) = U ( δ_σ_1 σ_3δ_σ_2 σ_4 - δ_σ_1 σ_4δ_σ_2 σ_3)  .

On the other hand, for the particle-hole vertices the relations between the V- and W-functions are different. For these functions, we can use another crossing relation <cit.>, namely

    ^ph,d_σ_1σ_2σ_3σ_4(k_1,k_2,k_3,k_4) = -^ph,c_σ_1 σ_2 σ_4 σ_3(k_1,k_2,k_4,k_3)  .

With this, we obtain the following identities:

    V^ph,d≡ W^ph,c ,     and  
    V^ph,c≡ W^ph,d ,

where we have abbreviated  etc. Together, these relations imply Eqs. (<ref>)–(<ref>) in the main text.

Next, one can reformulate the parquet equations (<ref>)–(<ref>) in terms of the spin-independent V-functions. This calculation is analogous to the derivation of the SU(2)-symmetric RG equations in Ref. SH00, and hence we only state the result here (cf. also Ref. [Eqs. (2)–(4)]julian):

    V^ph,d(p_1, p_2, p_3)     = k_B T/N∑_k G(k)   G(k + p_2 - p_3) 
          ×(  2  V^F(p_1, k + p_2 - p_3, k)   [V^F - V^ph,d](k, p_2, p_3) 
        - V^F(p_1, k + p_2 - p_3, k)   [V^F - V^ph,c](k, p_2, k + p_2 - p_3) 
        - V^F(p_1, k + p_2 - p_3,  p_1 + p_2 - p_3)   [V^F - V^ph,d](k, p_2, p_3) )  , 
    
     V^ph,c(p_1, p_2, p_3)     = -k_B T/N∑_k G(k)   G(k + p_3 - p_1)   V^F(p_1, k + p_3 - p_1, p_3)   [ V^F - V^ph,c](k, p_2, k + p_3 - p_1)  , 
    
     V^pp(p_1, p_2, p_3)     = -k_B T/N∑_k G(k)   G(p_1 + p_2 - k)   V^F(p_1, p_2, k)   [ V^F - V^pp ](p_1 + p_2 - k, k, p_3)  .

The full vertex F is correspondingly given by

    V^F = U + V^ph,d + V^ph,c + V^pp .

Furthermore, these SU(2)-symmetric parquet equations can be projected onto the various channels defined in the main text. Performing the same steps as in the derivation of Eqs. (<ref>)–(<ref>) for the spin-dependent vertices, we arrive at the following TU parquet equations for the spin-independent vertices (cf. Ref. [Eqs. (22)–(24)]julian):

    D̂[V^ph,d]_ℓ_1 ℓ_2(t)     = ∑_ℓ, ℓ'( 2  D̂[V^F]_ℓ_1 ℓ(t)   L^ph_ℓℓ'(t)  D̂[V^F - V^ph,d]_ℓ' ℓ_2(t) 
        - D̂[V^F]_ℓ_1 ℓ(t)   L^ph_ℓℓ'(t)  Ĉ[V^F - V^ph,c]_ℓ' ℓ_2(t) 
        - Ĉ[V^F]_ℓ_1 ℓ(t)   L^ph_ℓℓ'(t)  D̂[V^F - V^ph,d]_ℓ' ℓ_2(t) )  , 
    Ĉ[V^ph,c]_ℓ_1 ℓ_2(u)     = -∑_ℓ, ℓ'Ĉ[V^F]_ℓ_1 ℓ(u)   L^ph_ℓℓ'(u)  Ĉ[V^F - V^ph,c]_ℓ' ℓ_2(u)  , 
    P̂[V^pp]_ℓ_1 ℓ_2(s)     = - ∑_ℓ, ℓ'P̂[V^F]_ℓ_1 ℓ(s)   L^pp_ℓℓ'(s)  P̂[V^F - V^pp]_ℓ' ℓ_2(s)  ,

where the loop terms are again given by Eqs. (<ref>)–(<ref>). Thus, we have shown that the crossing relations, i.e., the antisymmetry of _pp and the initial interaction, as well as Eq. (<ref>), allow one to express all spin-dependent vertices in terms of the spin-independent V-functions and thereby to reduce the number of functions one has to keep track off.
We note that other implementations of the parquet equations such as Refs. held_new,jarrell_sym do not exploit these crossing relations explicitly but
instead enforce them during the iteration process to improve the convergence.


























































