Quantitative functional renormalization-group description of the two-dimensional Hubbard model

Using a forefront algorithmic implementation of the functional renormalization group (fRG) for interacting fermions on two-dimensional lattices, we provide a detailed analysis of its quantitative reliability for the Hubbard model. In particular, we show that the recently introduced multiloop extension of the fRG flow equations for the self-energy and two-particle vertex allows for a precise match with the parquet approximation (PA) {\sl also} for two-dimensional lattice problems. The refinement with respect to previous fRG-based computation schemes relies on an accurate treatment of the frequency and momentum dependences of the two-particle vertex, which combines a proper inclusion of the high-frequency asymptotics with the so-called"truncated unity"fRG for the momentum dependence. The adoption of the latter scheme requires, as an essential step, a consistent modification of the flow equation of the self-energy. We quantitatively compare our fRG results for the self-energy and momentum dependent susceptibilities and the corresponding solution of the parquet approximation to determinant Quantum Monte Carlo data, demonstrating that fRG is remarkably accurate up to moderate interaction strengths. The presented methodological improvements illustrate how fRG flows can be brought to a quantitative level for two-dimensional problems, providing a solid basis for the application to more general systems.


I. INTRODUCTION
Renormalization group (RG) methods have a long history in theoretical physics, ranging from a way to treat divergences in quantum field theories [1], critical phenomena [2,3] and quantum impurity problems [4,5] to current attempts to elucidate deep learning algorithms by physics [6].In general, RG methods connect specific quantities of a theory, like coupling constants or correlation functions, at a given scale with those at another scale via differential equations.This leads to a flow of these quantities which under appropriate circumstances distills out the dominating, to some degree universal, properties of the system.
The origins of the RG for electron lattice models date back to the second decade of the cuprate hightemperature superconductors more than 20 years ago [7][8][9][10][11].Here, the RG was utilized as a tool to deal with competing ordering tendencies in the Hubbard model (which are also seen in the cuprates) and as a method to understand in principle the stability of the Landau-Fermi liquid state [12].While the aptness of RG schemes in competing-order situations had already been known from impurity [4] and one-dimensional models [13], the systematic and versatile functional RG schemes [12,[14][15][16] turned out to be advantageous in the study of twodimensional (2D) lattice models such as the Hubbard model, beyond more general considerations [17,18].
The name functional RG (fRG) can be understood as having a twofold reason: on the one hand, the RG flow is derived from an exact flow equation for a generating functional of the theory when a (suitably chosen) parameter Λ in the free quadratic part of the action is changed.On the other hand, one usually investigates the flow of continuous functions of variables k such as wavevectors and frequencies with Λ, i.e., one deals with differential equations in Λ for functions f Λ (k) = f (Λ, k) of k and Λ.This marks a difference w.r.t. the conventional RG, where only a small finite number of constants is flowing.The main objects of interest in the fRG flows we discuss here are the electron single-particle selfenergy Σ Λ (k 1 , k 2 ) and the two-particle (interaction) vertex V Λ (k 1 , k 2 , k 3 , k 4 ).While the self-energy determines the changes of the single-particle excitations due to the mutual interactions, the two-particle correlation functions and collective properties are mostly controlled by the two-particle scattering processes, the so-called vertex correction terms.A dominant role played by the latter clearly emerged from realistic calculations of correlated metals [19][20][21].
Numerically, the costly parts are the evaluation of the fRG differential flow equations, usually determined by Feynman-type loop diagrams, as well as the high number of components of the flowing functions such as Initially, fRG in the form used here was developed arXiv:2002.02733v1[cond-mat.str-el]7 Feb 2020 and employed in 2D Hubbard models in view of high-T c cuprates.First works using so-called N -patch techniques to resolve the emerging wavevector dependence of the flowing interaction focused on the exploration of the leading ordering tendencies of the model [7][8][9][10] and understanding similarities to one-dimensional models.This included the question whether the pseudogap of the cuprates is foreshadowed in the fRG flow of weak to moderately coupled Hubbard models [11].Parallel developments were conducted in the field of inhomogeneous one-dimensional systems [22,23].Already at early stages, they were able to perform some quantitative benchmarking with exact numerical techniques (DMRG), and to include partial self-energy corrections into fRG.For 2D models, incorporating the fRG flow of the self-energy was technically more difficult for a number of reasons [11] but works indicating the opening of pseudogaps in the 2D flows came out a few years later [24,25].Yet, at that time, little attempts were made to assess the quantitative precision of the method in comparison with other theoretical methods.Rather, the Npatch fRG techniques were used in numerous applications to broad classes of experimentally relevant material systems, like the iron superconductors [26,27], graphene [28][29][30], and in the search for interaction-induced topological states [31,32].Here fRG was mainly used as a flexible tool to explore ground state phase diagrams over a wide range of parameters.
Besides these applications, the formal development saw a steady evolution.Recently, the long-standing challenge of relating fRG schemes to the parquet approach has been addressed [33,34].Within the so-called multiloop extension of fRG, all higher-loop contributions to the flow of parquet diagrams are accounted for with their exact weight.Conveniently, the effort only grows linearly with the loop order.The equivalence between the parquet approximation (PA) and multiloop fRG has been rigorously established at the level of spinless [33] and spinful [35] impurity models.Its numerical verification in the 2D Hubbard model poses additional challenges, due to the necessity of also treating the 2D-momentum variables.Specifically, we will demonstrate how it is possible to render the numerical effort for a quantitive treatment of the fRG flow manageable in 2D, by exploiting the multichannel decomposition [36] of the interaction vertex in combination with form factor expansions [30,37].This frees resources to incorporate frequency dependencies and selfenergy corrections in various cases [38,39], and to conduct flows into symmetry-broken states [40,41].Strongly inspired by earlier channel-decomposition schemes [37] and the singular-mode fRG by Wang et al. [30], the socalled truncated unity fRG was set up [42].This formalism combines various technical improvements and allows for a developement of a highly parallelizable and fast-performing code, mainly involving 2D-(or 3D-) integrations and matrix multiplications [43].The name 'truncated unity' comes from the insertion of unity into loop diagrams.These unities are sums over form factors, which are truncated by only considering the relevant form factors.As the form factors can be related to fermion bilinears on the real-space lattice, one obtains a physically appealing understanding what the given truncation captures and what is left out.Furthermore, one can check the convergence of the truncation by varying the number of form factors kept [42].Notably, the truncated unity fRG methodology is useful in another diagrammatic approach, the parquet scheme, whose performance can be boosted significantly in the truncated unity PA [44,45].A major step forward was achieved by including the frequency dependence and self-energy corrections into the truncated unity fRG framework [46] and related schemes [47,48].Supplemented by multiloop corrections [46], this paved the way for a) internal convergence of fRG as a function of the number of loops, and b) internal consistency in fRG, as different ways to compute response functions (by the flow of external-field couplings or by post-processing of the final interaction vertex) and flows with different cutoff schemes are found to agree.
In this paper we now aim at demonstrating the accuracy of our approach to describe one of the most challenging quantum many-body model of correlated electrons: the 2D Hubbard model.We do this by comparing data obtained by the truncated unity fRG in the multiloop extension (TU-mfRG, hereafter denoted by fRG* with the adapted self-energy flow and by fRG with the conventional self-energy flow) with other numerical methods.Going beyond the preliminary results of Ref. [46], restricted to the calculation of static response functions, we address here the computation of the self-energy in detail.In particular, investigating the form of the multiloop flow equation for the self-energy, we find that a form factor expansion of the two-particle vertex prevents the full reconstruction of the Schwinger-Dyson equation [49,50] (SDE) at loop convergence.The deviations can be traced back to the approximations introduced by the form factor projections in the different channels, and -more importantly -can be cured by using the direct derivative of the SDE instead.Including this methodological improvement, we provide a detailed analysis of the quantitative reliability of fRG for the 2D Hubbard model by verifying the agreement with the solution of the PA.Further, by comparing to Quantum Monte Carlo data, we show that fRG is remarkably accurate up to moderate interaction strengths.This demonstration may have considerable impact as it shows that the fRG can be pushed to become a quantitative quantum many-body method.Among its benefits for condensed matter material research are the unbiased treatment of various fluctuation channels down to low energy scales, including the indication of ordering transitions and (pseudo-)gap openings.This important feature is accompanied by the conceptual transparency of the approach that can be exploited to pin down the processes responsible for emerging physical effects.On a longer term perspective, a quantitatively reliable implementation of the fRG approach provides the most suited set-up for its combination [51,52] with complementary quantum many-body approaches, such as the dynamical mean-field theory (DMFT).The improvements described in this paper could represent, thus, an essential step for an accurate description of two-dimensional electron systems even in the intermediate-to-strong-coupling regime.
The paper is organized as follows: In Section II we introduce the Hubbard model and the main observables.In Section III we discuss the multiloop flow equation of the self-energy and present its extension in the truncated unity fRG framework, showing that it correctly accounts for the form factor projections in the different channels.Section IV contains a brief description of our benchmarking methods, the PA and determinant Quantum Monte Carlo (dQMC).In Sections V and VI, we illustrate our results for the self-energy and the susceptibilities for the 2D Hubbard model at half filling and out of it, provide numerical evidence for the convergence to the PA, and perform benchmarks with the dQMC data.We conclude with a summary and an outlook in Section VII.

A. 2D Hubbard model
We consider the single-band Hubbard model in 2D, where ĉ( †) iσ annihilates (creates) an electron with spin σ at the lattice site i (n iσ = ĉ † iσ ĉiσ ), t ij = −t is the hopping between neighboring and t ij = −t between next-nearest neighboring sites, µ the chemical potential, and U the on-site Coulomb interaction.The bare propagator is In the following we use t ≡ 1 as the energy unit.

B. Susceptibilities and self-energy
We compute different susceptibilities describing the linear response of a system to a weak external perturbation, as obtained from fRG*, PA and dQMC.In Matsubara frequency space, the susceptibilities are defined via the Fourier transform with respect to the imaginary time τ where η = M/D/SC indicates the magnetic, density and superconducting (s-and d-wave) channel, respectively.
Furthermore, we compute the self-energy

III. MULTILOOP FRG A. Channel decomposition of the vertex
Next to the self-energy, a central object of fRG and parquet algorithms is the one-particle-irreducible twoparticle vertex F .Using the SU (2) spin symmetry [53], we can restrict ourselves to one spin component, V = F ↑↑↓↓ .From the vertex, the susceptibilities can be computed either via the flow of the response vertices or by the contraction of the two-particle vertex at the end of the flow ("post-processed") [46].
In the channel or parquet decomposition of the vertex we can identify the two-particle reducible contributions Φ pp/ph/ph in the particle-particle, particle-hole and crossed (or transverse) particle-hole channel, respectively.We have where the reducible vertices on the r.h.s.have been parametrized according to a single generalized transfer momentum and two fermionic ones.In the parquet approximation, the two-particle irreducible vertex is approximated by Λ 2PI = U .Fully accounting for the interplay between different two-particle channels is a central motivation for the multiloop extension of fRG, described next.

B. Multiloop extension of fRG: a brief overview
The Wetterich equation is an exact one-loop (1 ) flow equation for the generating functional of one-particle irreducible vertices [15].Expanding in the vertices leads to an infinite hierarchy of one-loop flow equations for the vertices.However, objects like the three-particle vertex are intractable for numerical treatments.The fundamental approximation in many fRG flows is therefore the truncation in the hierarchy of flow equations [54].Setting the three-particle vertex to zero yields an approximate 1 flow equation for the self-energy and two-particle vertex.A common way to reintroduce some of the lost contributions is to reuse the scale derivative of the self-energy Σ in the flow of the two-particle vertex.This Katanin substitution [55] already leads to significantly improved results, labelled by 1 K in the following.A further refinement, which effectively incorporates the three-particle vertex to third order in the renormalized interaction, is realized by reusing the 1 results in a 2 addition to the vertex flow [40,55,56].
The multiloop fRG [33,34] extends these schemes to arbitrary loop order.One finds that the multiloop additions complete the scale derivative of all the diagrams which are only partly generated in the 1 flow.Thus, it lift the dependence of the results from the particular choice of the regulator-the quantitative reliability of these results is the subject of the present paper.Due to the iterative structure of the multiloop corrections, the numerical effort grows linearly with the loop order.Starting from an efficient algorithm for the 1 flow, the implementation of the multiloop equations is straightforward.
An alternative derivation of the multiloop flow equations [57] highlights the close connection between fRG and the parquet approach, see also Section III C. Starting from the parquet equations, a scale dependence of propagators and vertices can be introduced by making the bare propagator scale dependent.Taking the scale derivative of the self-consistent parquet equations then leads to the multiloop flow equations.The equivalence is exact, provided that the self-energy and vertices are treated without further approximation of their momentum and frequency dependence.
A crucial ingredient for such an equivalence and overall quantitative accuracy is a good resolution of the twoparticle vertex and the self-energy in momentum and frequency space.Since the numerical effort grows extremely fast with the number of momentum patches and the size of the frequency window, an efficient vertex parametrization of the two-particle vertex is very important.We use the truncated unity fRG [42,58] for the momentum dependence, together with an accurate treatment of the frequency dependence which includes the high-frequency asymptotics [59].For the detailed description and their implementation in the multiloop extension we refer to Ref. [46], where we also provide the expression of the employed smooth frequency cutoff.Here we further use a refined momentum grid to resolve the peak at q = (π, π) in the AF susceptibility at half filling.This can be easily accounted for in the truncated-unity formulation with precalculated projection matrices in real space [46].The accurate description of the long-range AF fluctuations is also of major importance in order to fulfill the Mermin-Wagner theorem [60].
The effect of the multiloop extension on the frequency structure of the two-particle vertex V is exemplified in Fig. 1.Here, we compare 1 K with fully converged fRG* results, where the latter are obtained by using the revised multiloop flow equation of the self-energy presented in Section III D. While Φ pp is enhanced in fRG*, Φ ph and Φ ph are screened and the absolute values are smaller.The screening of Φ ph with more loops is also reflected in the AF susceptibility, shown in Fig. 2, which is dominated by the contribution of the magnetic channel Φ M = −Φ ph .The 1 -scheme strongly overestimates the peak at momentum transfer q = (π, π) leading to an AF ordering at finite interaction strength in contrast with the Mermin-Wagner theorem [60].With increasing loop order, the AF peak is reduced.At an inverse temperature of 1/T = 5 and on the scale of this plot, the differences in the susceptibility between the 2 result and fRG* are hardly visible.We note that for U = 3 fRG* is not fully converged w.r.t.loop order and for this reason no results for larger values of U are displayed.The convergence threshold we use is of 1% for χ AF , Im Σ(k = (π, 0), iν 0/1 ) and Im Σ(k = (π/2, π/2), iν 0/1 ).For completeness, we report the parameters used for the benchmark analysis in Sections V and VI in Appendix A.

C. Parquet approximation and post-processing of the self-energy
In many-body theory, there are various exact relations between one-and two-particle quantities, such as the Bethe-Salpeter equations [61,62] connecting different parts of the two-particle vertex and the SDE relating self-energy and vertex.While these equations are used in parquet approaches to iteratively find a self-consistent solution on the one-and two-particle level, applying them to the final self-energy and vertex of any method is a way to check its consistency.
In our previous work [46], we focused on two-particle quantities like the susceptibilities and compared their outcome directly from the fRG flow with their postprocessed result using the final self-energy and vertex.Here, we use the SDE to analogously determine the selfenergy.The SDE involves the self-energy itself through the full propagator as well as the vertex and reads FIG. 2. Antiferromagnetic susceptibility χAF(iω = 0) defined in Eq. 5 as a function of the bare interaction U , for 1/T = 5.
Here the sum over Matsubara frequencies includes the normalization factor of T ; the fourth dependence of the vertex can be recovered by momentum and frequency conservation.Its diagrammatic representation is shown on the left-hand side of Fig. 3.Note that we take the Hartree part implicitly into account by shifting the chemical potential by U nσ ; half filling then corresponds to µ = 0.
In the parquet decomposition (8), the vertex is split into the fully irreducible part and the three two-particle channels.As depicted on the right-hand side of Fig. 3, the SDE is then determined by four parts

Σ(k
where f n (k) are the form factors in TUfRG.The first diagram can be calculated in real space using the convolution theorem twice.The remaining parts in the ph-, ph-, and pp-channel are computed by using where we used the same conventions as in Ref. [46].Applying the above Eq.( 10) as post-processing procedure to compute the self-energy at the end of the fRG flow, we find deviations of up to 20% in both the real and the imaginary part with respect to the solution of the conventional flow equation where the first lines corresponds to the 1 flow equation.
While at 1 , a difference between the flowing and the post-processed result is not surprising, we expect these differences to vanish in a (loop) converged multiloop fRG solution.As we will show in the following, the remaining discrepancies originate from the truncated form factor expansion, for which the flow of the self-energy has to be replaced by the direct derivative of the SDE [63].

D. Self-energy flow in a form factor expansion
The different fRG results for the self-energy obtained from the (multiloop) flow and the post-processing via the SDE can be traced back to the truncated unity treatment with a reduced number of form factors [64].In this case, some of the identities used in the general derivation of the self-energy flow [57] do not hold any more.Consequently, the approximation of the vertex in terms of form factors destroys the equivalence of the flow equation and the SDE.However, this problem can be overcome by using the direct derivative of the SDE instead.
In more detail, each summand of the SDE contains two vertices and three propagators.In Ref. [57], multiple transformations which interchange these propagators have been used, amounting to a translation between the different two-particle channel descriptions.However, in the truncated unity parametrization, channel transformations are only information-loss free in the infinite form-factor limit.With a finite number of form factors, the invariance of the SDE under the exchange of propagators does not hold anymore.
Here we propose the fRG* extension which exactly reproduces the SDE equation for the self-energy in a form factor expansion of the two-particle vertex.Technically, the flow of the self-energy is replaced by the direct derivative of the SDE (for details we refer to the Appendix B), which can be found by introducing a scale dependence in the SDE and taking the derivative with respect to it [57].An example illustrating the advantage of the new fRG* self-energy flow compared to the conventional flow is shown in Fig. 4. We focus on two specific differentiated diagrams contributing to the flow of the self-energy.In fRG*, each summand of the SDE is directly differentiated w.r.t.Λ.For concreteness, we consider the last summand of Fig. 3, insert the lowest-order diagram for Φ ph , and let the Λ derivative act on the two spin-down propagators (dashed lines).Thereby, we obtain the two self-energy diagrams on the r.h.s. of Fig. 4 (upper panel), where the propagators with a diagonal dash symbolize differentiated propagators.The same contributions should be part of the standard fRG self-energy flow, shown in the lower panel of Fig. 4. Indeed, the first diagram on the r.h.s.simply follows from a second-order Φ ph diagram with two ph bubbles, and the second diagram originates from the Φ pp part of F , where the ph bubble (orange) is inserted into another pp bubble (blue).The crucial point is that the r.h.s. of both panels are formally equivalent, but the form-factor truncation applies in a less favorable way on the right diagram in fRG (lower panel): The two diagrams in the upper panel and the first in the lower panel are exactly described by only s-wave form factors.However, in the lower right diagram, the ph contribution is averaged in the process of translating it to the pp channel due to the s-wave form factor truncation.
In Fig. 5 we present the real and imaginary part of the self-energy Σ(iν = iπT ) as a function of momentum, for U = 2 and 1/T = 5, both for the conventional fRG (blue) and fRG* (red).We show that unlike fRG, fRG* yields an excellent agreement between the flowing (solid lines) and the post-processed self-energy (dashed lines), which is determined by Eq. ( 10) with the final vertices and selfenergies at the end of the flow.The non-improved fRG results obtained from the flow exhibit pronounced deviations with respect to the fRG* self-energy.For the post-processed ones, these deviations are significantly reduced since in this case the self-energy is updated by the SDE at the end of the flow.

E. Self-energy iterations
We now analyze the effect of the self-energy iterations in a complete multiloop flow, in fRG as well as in fRG*.As described in detail in the Appendix B, these are needed because the Katanin substitution ) in the vertex flow depends on the total self-energy flow-including the multiloop corrections to the self-energy flow [34] or parts of the differentiated SDE which in turn involve the vertex flow.In order to study their effect, we compare the fully converged results to the ones obtained by solely inserting the 1 self-energy flow in the Katanin substitution.In Fig. 6 we show the self-energy results with (solid lines) and without iterations (dotted lines), as obtained by fRG (red) and fRG* (blue).We note that the slight kink between the 4th and 5th frequency and also between the 8th and 9th frequency corresponds to the crossing of the low frequency tensor range and the high frequency asymtotics of the two-particle vertex.This ef- fect is more pronounced in fRG than in fRG* since the conventional self-energy flow involves a direct contraction of the vertex as opposed to a nested integration in the SDE.In the conventional fRG flow, the effect of the self-energy iterations amounts to 5 − 10% except for the first Matsubara frequency which appears to be well described without any additional iterations.In contrast, in fRG* which accounts for the form factor projections in the different channels, the self-energy iterations lead to much smaller overall corrections but are of relevance for the lowest frequencies (see also Fig. 17 in Appendix C for the effect on the momentum dependence).
The self-energy iterations also affect the AF susceptibility χ AF displayed in Fig. 7.For U = 2 and 1/T = 5, a sizeable effect is only observed for zero Matsubara frequency, highlighted in the inset.Neglecting the selfenergy iterations in fRG* (red dotted lines) overestimates the AF peak by 2%.In fRG, their impact on χ AF is smaller, according to a small effect on the first Matsubara frequency of the self-energy observed in Fig. 6.Anticipating the comparison of fRG* (red) with PA (grey) in Fig. 8, we conclude that while the conventional scheme (blue) leads to a deviation of 6% for zero Matsubara frequency, fRG* (red) shows a perfect agreement.This is confirmed also by the post-processed result (red dashed lines).

A. Parquet approximation
The PA results where obtained with the truncated unity implementation of the parquet equations [45].The parquet equations are solved by iterating the Bethe-Salpeter equations, the parquet equation ( 8) and the SDE (10) until self-consistency is reached.The momentum dependence of the vertices is parametrized using the form factor expansion [44], formally identical to the scheme used in the truncated unity fRG.Although the equivalence of the PA and the multiloop fRG has been formally shown [57], the actual PA calculations substantially differ from the ones in fRG, since no flow parameter or cutoff dependence is introduced.In the PA no differential equations are solved and the convergence to a fixed point is achieved iteratively, starting from an initial guess for the vertices (in our case given by the lowest order diagrams).
In order to account for the finite frequency box, we use the asympotics as introduced in Ref. [65] and also used in Ref. [66].The treatment of the frequency asymptotics thus differs from the one used in the fRG calculations.All further computational details of the truncated unity implementation of the parquet equations can be found in Ref. [45], as well as a detailed analysis of the convergence in the number of form factors, showing that the single form factor approximation employed here at half filling is justified.
All results are converged in the number of discrete lattice momenta N q and positive fermionic Matsubara frequencies N f + .Specifically, we use a uniform momentum grid with N q = 32 × 32, and a frequency box with N f + = 32 positive Matsubara frequencies.

B. Determinant quantum Monte Carlo
The dQMC algorithm, proposed by Blankenbecler, Scalapino, and Sugar [67], is a state-of-the-art, numerically exact method and is commonly applied for finitetemperature [67,68] calculations of interacting fermion systems.The basic idea of the dQMC algorithm is to decouple the two-body interaction into noninteracting fermions coupled with auxiliary fields, and to compute the fermionic observables via importance sampling of the fields.To achieve that, a Trotter decomposition and a Hubbard-Stratonovich (HS) transformation are successively used, after the discretization for the inverse temperature as β = M ∆τ .The systematic error from finite ∆τ can be removed by extrapolations with several different ∆τ values.For further details, we refer to reviews [69,70].In this work, we choose ∆τ t = 0.02 which has been tested to safely reach the ∆τ → 0 limit.In this work, we have also implemented our most recent improvements [71,72] of the dQMC algorithm.For the computation of dynamical quantities, we first measure the imaginary-time correlation functions, and then obtain the imaginary-frequency observables via Fourier transformation.Specifically, for the self-energy we implemented the Legendre polynomial representation [73] for the imaginary-time single-particle Green's function G(k, τ ) to compute G(k, iν), and subsequently applied the Dyson equation.This yields smooth self-energy re-sults even for high frequencies.All dQMC data presented here are found to converge to the thermodynamic limit for a linear system size of L = 28 (with the number of lattice sites being N = L 2 ) for half filling and L = 24 away from it.As for statistics, we typically use in total 10 5 measurement samples after the Markov Chain equilibrium process.The error bars are significantly smaller than the corresponding symbol and thus neglected in the plots.

V. RESULTS AT HALF FILLING
We now compare different physically relevant quantities as obtained from fRG*, the PA, and the numerically exact dQMC.In particular, we first focus on the various susceptibilities in Section V A and then present the results for the self-energy in Section V B.

A. Susceptibilities
We first present the results for the leading AF susceptibility in the half-filled 2D Hubbard model χ AF as a function of the bare interaction strength U .In Fig. 8 we report fRG* (red), PA (grey), and dQMC data (black), together with the relative difference of fRG* with respect to PA and dQMC shown in the inset.Up to U = 2, fRG* and PA coincide with a relative difference of ≤ 1% (indicated by the green shaded area).For larger values of U , the convergence of fRG* in frequencies and also in loop numbers becomes numerically challenging and is not reached yet [74].This leads to the observed deviations from the PA solution.The differences between the PA and the numerically exact dQMC data are essentially due to the fully two-particle irreducible diagrams not included in the PA.These diagrams contribute to fourth order in U ; the corresponding relative difference amounts to ∆ rel 0.003 U 4 .A second source of the differences between the PA solution and dQMC is given by the form factor expansion of the two-particle vertex which accounts only for the local s-wave part.Due to perfect nesting, the physics at half filling is dominated by magnetic fluctuations peaked at q = (π, π) and at small coupling, there is only a small quantitative correction due to the form factor truncation [45].Away from half filling, where we expect superconducting d-wave components to become relevant, we will extend the form factor truncation to include them in Section VI.
Figure 9 shows the magnetic susceptibility at zero frequency χ M (q, iω = 0), for U = 2 (and 1/T = 5).The results of fRG*, the PA, and dQMC exhibit the same momentum dependence, with an excellent quantitative agreement.The largest deviation is found at M = (π, π) which corresponds to the AF susceptibility shown in Fig. 8 for different values of U .We note that for all other frequencies iω = 0 the AF susceptibility of fRG* perfectly agrees with the one of the PA.While the AF peak height obtained from fRG* does not converge perfectly to the PA for U > 2, the correlation length ξ extracted from its width shows a very good agreement between the different methods, see inset in Fig. 9.The correlation length is fitted to all points of χ M (q, iω = 0) within a distance of 0.3π from M , through which reduces to the Ornstein-Zernike form for small momentum differences q x − π and q y − π (cf.Ref. [75 and 76]).The number of momenta taken into account for the fit are between 33 and 45 in fRG*, 69 and 161 in PA, and 57 in dQMC.The maximal standard deviation error is 0.027 in fRG*, 0.025 in PA and 0.028 in dQMC.
In Fig. 10 we show different subleading susceptibilities as a function of U , for 1/T = 5.In fRG*, the compressibility κ and the (s-wave) charge density wave χ CDW susceptibilities are obtained from the flow.Note that the equivalence χ SC,s = χ CDW holds in general for SU (2) spin and charge (particle-hole) symmetry, see Appendix D for the proof.In contrast, the d-wave superconducting susceptibility χ SC,d is determined through a post-processing of the vertex at the end of the flow.The quantitative agreement of fRG*, PA and dQMC results for the subleading susceptibilities is very good, with a relative difference of fRG* w.r.t. the PA at U = 3 being 13% for κ, 10% for χ CDW , and 1.5% for χ SC,d , and w.r.t.dQMC 18% for κ, 12% for χ CDW , and 2.6% for χ SC,d .Note that the compressibility κ is also consistent with Ref. [77].The good agreement between fRG* and PA, both affected by the form-factor truncation, and the exact dQMC justifies a computation with only the local s-wave form factor.
All subleading susceptibilities κ, χ CDW , χ SC,s and χ SC,d decrease with U , in contrast to χ AF .We ana-FIG.9. Magnetic susceptibility χM(q, iω = 0) as obtained by fRG* (red), the PA (grey), and dQMC (black), for U = 2 and 1/T = 5.The inset shows the correlation length ξ extracted from χM(q, iω = 0) as a function of U , for 1/T = 5 (see text for the details of the fitting procedure).lyze this behavior by presenting a more detailed (postprocessing) analysis of the different contributions to the susceptibility and in particular of the importance of the vertex corrections, see Fig. 11.Comparing the postprocessed susceptibilities to the ones obtained from the flow of the response functions, provides an indication of the fRG* convergence with respect to momenta, frequencies, and loop number.The agreement is within numerical accuracy except for χ AF , where for U > 2 we have difficulties to converge the fRG* calculations in frequencies and loop numbers.We furthermore show the 'uncorrelated' susceptibilities in terms of dressed Green's functions (without vertex corrections) where the form-factor index 0 stands for s-wave and 1 for d-wave.They all decrease with U , as a consequence of self-energy screening effects.The vertex contributions, given by the difference to the susceptibilities, exhibit a richer physical behavior: they lead to a reduction or screening of the bare κ and χ SC,s susceptibilites, whereas χ SC,d and most prominently χ AF are enhanced with respect to their bare values.For χ SC,d the vertex corrections are not strong enough to induce an overall increasing susceptibility.This occurs only for χ AF , where the vertex corrections are dominant.

B. Self-energy
We now discuss the frequency and momentum dependence of the self-energy and their comparison as obtained from the different methods.In Fig. 12 we show the imaginary part at the nodal k = (π/2, π/2) and antinodal k = (π, 0) point as a function of U , for 1/T = 5.The agreement between fRG* and the PA is almost perfect for small values of U , with increasing deviations for larger U .In particular, for the first Matsubara frequency (upper panels) the relative difference between fRG* and PA amounts to only 2% for U = 3, while for the second Matsubara frequency (lower panels) the deviations are more pronounced, reaching 6%.In contrast, the dQMC results for the first Matsubara frequency differ considerably for U = 3, with an important physical implication: the onset of the pseudogap opening [78][79][80] in dQMC resulting from quasi equal values for the first two Matsubara frequencies at k = (π, 0) is not observed in fRG* and the PA.There, the absolute value for the first Matsubara frequency is 11% smaller and therefore the gap opening will set in only at larger interactions.
We also compare the behavior of the self-energy as a function of frequency in Fig. 13, for a representative value of U = 2 (and 1/T = 5).At small frequencies the selfenergy shows typical Fermi-liquid behavior, both at the nodal and the antinodal point.The antinodal point is affected more strongly by correlation effects, with an increased absolute value for the lowest Matsubara frequencies.In general, both the upper and lower panel indicate that for these parameters the resulting self-energy does not develop a momentum selective gap.In Fig. 13 fRG* and the PA exhibit larger deviations in the intermediate to high frequency range, despite the excellent agreement with the corresponding post-processed results (see also Fig. 18 in Appendix C).These differences are due to the specific implementations of the high frequency asymptotics of the the two-particle vertex [59]: in fRG* the asymptotic functions are calculated and stored explicitly [46], retaining a smaller tensor for the low frequency range with respect to the PA, where a large tensor over many fermionic and bosonic frequencies is used and the values outside are constructed from the ones at the edges [65].The former is numerically more efficient, but has the drawback of kinks arising at the transition between the different tensors, see Fig. 13.On a quantitative level, in the full Green's function G(k, iν) the differences between the frequency dependence of fRG* and the PA self-energy are almost negligible w.r.t. the large iν contribution of the bare Green's function G 0 (k, iν), see Eq. ( 2).We checked that for smaller interactions excellent convergence in frequencies, momenta and loops can be achieved (compare Fig. 19 for U = 0.5 in Appendix C).
Finally, a comparison of the self-energy as a function of momentum for the same representative parameters of U = 2 and 1/T = 5 is performed in Fig. 14.Concerning the agreement of fRG* with the PA (grey), we observe small differences with respect to fRG* only in the imaginary part of the self-energy, for momenta which are far away from the Fermi surface (see lower panel).Since it is the Green's function and not self-energy that directly enters the calculation of observables, these momenta have little influence.Moreover, the relative difference between the methods for the Green's function at the Γ and M points is less than 3‰ [81].Also for the self-energy as a function of momentum, the differences disappear for smaller interactions and we find that fRG* perfectly reproduces the PA solution, see Fig. 20 for U = 0.5 in Appendix C.

VI. RESULTS AWAY FROM HALF FILLING
In presence of a finite doping and an additional nextnearest neighbor hopping t , the physical behavior is much richer and not exclusively driven by AF fluctuations any more.Since we expect the superconducting d-wave component of the two-particle vertex to become more important here, we include also the d-wave form factor.
In the following we present fRG* results for the evolution of the different susceptibilities away from half filling, together with their comparison to PA and dQMC data.Specifically, we consider the parameters t = −0.2 for the next-nearest neighbor hopping and µ = −0.35,−0.7, −1.4,−2 for the chemical potential.Due to the self-energy flow in fRG* the initial chemical potential is renormalized leading to a different filling at the end of the flow.This effect is very small close to half filling and increases with the doping δ = 1− n , see Fig. 15 where the structure of the magnetic susceptibility χ M (q, iω = 0) in momentum space is shown for U = 2 and 1/T = 5.In Fig. 16 we report the compressibility κ, the charge density wave χ CDW , and superconducting χ SC (s-and dwave) susceptibility as a function of doping, for the same parameters.We note that here χ CDW and χ SC,s are not equivalent any more.
The magnetic susceptibility dominates for all dopings and is maximal at the commensurate antiferromagnetic wave vector (π, π) for µ = −0.35 and µ = −0.7 and at incommensurate wave vectors for larger values of the doping, consistent with previous fRG findings [48,82].In particular, we did not find a pairing instability for 1/T = 5 at any doping.The subleading susceptibilities are still reduced with respect to the (in-)commensurate peak in χ AF , but for increasing δ their difference decreases and for µ = −2 the maximum in Fig. 15 is only twice as large as χ SC,d .While the AF susceptibility gradually evolves from the beginning of the flow [9,37,48,83], the superconducting d-wave susceptibility emerges only  in proximity of the critical scale.This indicates that the AF fluctuations are responsible for the d-wave pairing.The parameter regime presented here is far away from any instability.Hence, for a finite doping we expect the d-wave pairing susceptibility to increase only at lower temperatures.Due to the high computational cost of low-T calculations (specifically to be able to accurately parametrize the frequency dependence of the two-particle vertex), we cannot access the superconducting transition temperature at the moment.For the temperatures considered here, the onset of large d-wave pairing interaction is likely a high-temperature precursor of a superconducting phase at lower temperature: as the temperature is further decreased, the relevance of the d-wave pairing should increase until, eventually, the d-wave pairing becomes larger than the magnetic one.
The agreement with the PA and the numerically exact dQMC data is very good also away from half filling.We refrain here from providing relative differences because the data correspond to different fillings.Moreover, due to the high numerical cost the present calculation including s-and d-wave form factors is not fully converged in frequencies.This hardly affects the susceptibilities, while the quantitative accuracy of the self-energy appears to be more sensitive.

VII. CONCLUSIONS AND OUTLOOK
In our work, we illustrated how it is possible to achieve, by means of fRG, a quantitatively accurate description of correlated electrons on two dimensional lattices, by implementing proper enhancements to the conventional algorithms.Our starting point was the significant progress recently obtained [46], which combined the truncated unity fRG [37], a clever frequency representation [53,59] and the multiloop extension of the approach described in Refs.[33,34].While the latter advances suffice for computing impurity models [33,35], the missing piece for a quantitatively accurate description of the electronic correlations in 2D is to make the self-energy flow consistent with the truncated-unity scheme.Specifically, we show that replacing the corresponding flow equation by the direct derivative of the Schwinger-Dyson equation (SDE) allows us to sum up the contributions of the different channels in the correct proportion.The post-processed computation from the propagator and interaction vertex at the end of the flow exactly fulfills the SDE.As a consequence, the resulting self-energy is independent of the chosen cutoff scheme.
This methodological improvement is needed for converging-at a high degree of numerical accuracy-the multiloop fRG results to the PA ones and for obtaining quantitative reliable fRG data for the 2D Hubbard model at half-filling as well as upon hole doping.In particular, by comparing the converged fRG data to the PA and dQMC, a satisfactory agreement between the corresponding values of the self-energies and physical response functions could be established up to intermediate interaction strength.We stress that such quantitative agreement for the 2D Hubbard model cannot be obtained by exploiting conventional (e.g., 1-loop-based) truncations of the fRG flow.Minor deviations between fRG and PA (on the one side) and dQMC (on the other side) are instead observed, as expected, by increasing the interaction values.Further optimization and parallelization of the code [43] will allow us to overcome the restriction to essentially a single s-wave form factor.This is necessary to explore broader parameter regions.
The presented advancements of fRG-based computation schemes constitute the basis for its extensions to more general systems and for its combination [51,52] with non-perturbative many-body methods, such as DMFT.Note that also for ab-initio investigations, the consistent summation of different scattering channels is of clear importance [84], and fRG might prove useful in this regard as well.Hence, our study paves a promising route towards quantitative fRG analyses of electronic phase diagrams at all coupling strengths and for an fRGbased investigation of emerging energy scales, competing instabilities and response functions in wider classes of quantum materials of high relevance for the cutting-edge theory in condensed matter research.Calculations have been done in part on the Vienna Scientific Cluster (VSC).The authors also gratefully acknowledge the computing time granted through JARA on the supercomputer JURECA at Forschungszentrum Jülich [85].We report the technical parameters for the results of Sections V and VI, in Table I  As the calculations with an additional d-wave form factor are numerically challenging, the frequency range was fixed to N f + = 2.
Appendix B: Implementation of the self-energy flow In the fRG* approach we replace the conventional flow equation of the self-energy together with its multiloop corrections by where the Hartree part is implicitly taken into account by shifting the chemical potential by U nσ (with µ = 0 for half filling).The first part does not depend on the effective vertex and reproduces second order perturbation theory with renormalized propagators.Applying the convolution theorem twice, it can be calculated efficiently by using Fast-Fourier-Transform routines F in real space Each of the vertex-dependent contributions to the selfenergy flow (B1), Σph (k, iν), Σph (k, iν) and Σpp (k, iν) in the respective channels, include similar contributions regarding the Λ derivative: Applying the product rule, we obtain the derivative of the two-particle reducible vertex, of the excitation bubble, and of the propagator closing the external loop.As the full expression of the self-energy flow equation can be easily derived from the SDE (10), here we report only the ph-contribution in its explicit form with V Λ=0,s s = 4π 2 U due to the normalization of the form factors [46].The three parts are similar in structure and can be calculated using the same procedure just exchanging the vertex or the propagator at a specific scale with its derivative at this scale.In practice, the bubble Π Λ ph,m n (q, iω, iν) and its derivative are already calculated for the conventional 1 self-energy and the 2 vertex flow.The multiplications inside the square brackets are of the same form as the ones on the r.h.s. of the flow equations for the vertex, where to the right of the excitation bubble only the bare vertex is inserted.The resulting effective vertices are projected into the purely fermionic notation by the form-factors f * m (k) and f 0 (k) and then contracted with a propagator (or scale derivative of a propagator) as in the conventional 1 -flow equation for the self-energy.Hence, only the computation of ΣGGG (k, iν) has to be implemented separately.
Both the fRG* as well as the conventional multiloop corrections of the self-energy flow can only be calculated after the r.h.s. of the vertex, as they depend on the scale derivative of the vertex.The resulting derivative of the self-energy which affects the Katanin substitution of the single-scale propagator has then to be included iteratively until the change in Σ is negligible.In contrast to the conventional 1 flow equation for the self-energy, the selfenergy flow equation in 'Schwinger-Dyson form' depends explicitly on Ġ and therefore also on Σ.Therefore, we use the conventional 1 flow as first estimate of Σ from which to start the iteration.Using the 1 flow in the Katanin substitution yields the results of Figs. 6, 7 and 17, as compared to the full iteration.
Appendix C: Supplementary self-energy results

Self-energy in the different fRG schemes
In addition to the results presented in Section III D, we illustrate in Fig. 17 the effect of the self-energy iterations on the momentum dependence of the self-energy at iν = iπT , both for fRG and fRG*.This supplements the results for the frequency dependence for U = 2 and 1/T = 5 shown in Fig. 6.Consistent with the observation there, the self-energy iterations lead to small corrections in fRG*, while in the conventional fRG flow the effect on the lowest Matsubara frequency is almost negligible.
Furthermore, in Fig. 18 we present the post-processed results for the self-energy as a function of frequency, in analogy to the ones for the momentum dependence shown in Fig. 5. Also here, fRG* which accounts for the form factor projections in the different channels yields perfect agreement between the flowing and the post-processed results for the whole frequency range, unlike the conventional fRG scheme.
2. Comparison between fRG* and PA at U = 0.5 The differences in the self-energy as obtained from fRG* and the PA disappear for smaller interactions and also by increasing the size of the low-frequency tensor, as can be seen in Figs.19 and 20.For U = 0.5 we find that fRG* perfectly reproduces the PA solution for the self-energy.These calculations require less momentum patches and therefore a larger low-frequency box with 8 (instead of 4) positive fermionic frequencies can be used, suppressing the box effects present at U = 2, where 6 were used.Appendix D: Equivalence χSC(q = 0) = χD(q = π) at half filling The proof of equivalence between χ SC and χ CDW at half filling traced in Ref. [86]  (D1) This charge symmetry implies the degeneracy between the local s-wave spin-singlet pairing and the charge density wave channel χ SC (q = 0) = χ D (q = π), which we observe also in the numerical results presented in Section V.
In order to prove the charge SU (2) symmetry for the Hamiltonian (D1), we introduce the following spinor operators, ĝi = ĉi↑ (−1) i ĉ † where Pi and P † i are given by Pi = 1 2 (−1) i ĝ † i (iτ 2 )(ĝ † i ) T = ĉ † i↑ ĉi↓ (D4a) with τ 2 the second Pauli matrix.Performing a global SU (2) transformation for the spinor operators ĝi → U † ĝi and ĝ † i → ĝ † i U , where U ∈ SU (2) is a 2 × 2 matrix which does not depend on the lattice site, we can show that the Hamiltonian (D3) is invariant under this transformation.For the first hopping term, it is clear that (ĝ † i ĝj + h.c.) is invariant under this transformation.For the interaction term, we consider the operator Pi under the above SU (2) transformation Pi and applying the property det(U ) = 1 of the 2 × 2 SU (2) matrix, we can actually prove that U (iτ 2 )U T = iτ 2 .This means that under the above transformation, the operator Pi is invariant, and so is the P † i operator.Thus, the total Hamiltonian (D3) is invariant under the SU (2) transformation.In the main text, the density susceptibility χ D (q) of Eq. ( 6) is defined by using (−1) i Di and the superconducting one χ SC (q) of Eq. ( 7) by using Re ∆i .Thus, as a result of the charge SU (2) symmetry as described above, the susceptibilities are related by χ SC (q = 0) = χ D (q = π).
Both finite next-nearest-neighbor hopping amplitudes and doping away from half filling break the above charge SU (2) symmetry of the Hamiltonian (1), lifting the degeneracy between the on-site s-wave spin-singlet pairing and charge density wave channel, as also demonstrated by the results of Section VI.

FIG. 1 .
FIG. 1. Two-particle vertex decomposed as in Eq. 8 in the different s-wave channel contributions at zero bosonic frequency as a function of fermionic frequencies, 1 including the Katanin correction (upper panel) vs. multiloop data (lower panel), for U = 2 and 1/T = 5.

FIG. 3 .
FIG. 3. R.h.s. of the Schwinger-Dyson equation for the self-energy (10), illustrating the contributions of the different channels, where solid (dashed) lines carry spin up (down).The first diagram on the right-hand side can be calculated using the convolution theorem and Fast-Fourier transform algorithms.The other contributions can be determined by first combining the two-particle reducible vertex (Φ) and the bare vertex (dot) through the propagator pair of the corresponding channel (red).Finally, each diagram is closed by a propagator (black) through the direct summation over frequency and momentum.

FIG. 4 .
FIG.4.Illustration of the self-energy flows in fRG*, restricted to the part with Φ ph (upper panel), and fRG (lower panel).The r.h.s.shows two exemplary differentiated diagrams contributing to the self-energy flow at third order, where solid (dashed) lines carry spin up (down), and the diagonal dash symbolizes a scale-differentiated bare propagator.In fRG*, two bubbles from the same channel (colored) are combined and then closed with the black line.By contrast, in fRG, the second diagram requires to insert a ph (orange) into a pp (blue) bubble, before closing with the differentiated propagator (black).

FIG. 5 .
FIG. 5. Real and imaginary part of the self-energy as obtained by conventional fRG (blue) and the fRG* flow (red), together with the respective post-processed results (dashed lines), for U = 2 and 1/T = 5.Within fRG*, the postprocessed (dashed red) lie exactly on top of the fRG* flow results (red solid).

FIG. 6 .
FIG. 6. Imaginary part of the self-energy at the nodal and antinodal point, as obtained by fRG (blue) and fRG* (red), with and without self-energy iterations, for U = 2 and 1/T = 5.Inset: relative difference.

FIG. 12 .
FIG.12.Imaginary part of the self-energy at the nodal and antinodal point as a function of U , as obtained by fRG* (red), the PA (grey), and dQMC (black), for 1/T = 5.
Appendix A: Tables with technical parameters of fRG*

FIG. 17 .
FIG. 17.Real and imaginary part of the self-energy as obtained by the conventional fRG (blue) and fRG* (red), with and without self-energy iteration respectively, for U = 2 and 1/T = 5.For comparison, also the result of the PA is shown (grey).

FIG. 18 .
FIG. 18. Imaginary part of the self-energy at the nodal and antinodal point, as obtained by fRG (blue) and fRG* flow (red), together with the respective post-processed results (dashed lines), for U = 2 and 1/T = 5.

2 i
FIG. 19.Imaginary part of the self-energy at the nodal and antinodal point, as obtained by fRG* (red) and the PA (grey), for U = 0.5 and 1/T = 5.The relative difference is always below 1%.

FIG. 20 . 1 2
FIG.20.Real and imaginary part of the self-energy as obtained by fRG* (red) and the PA (grey), for U = 0.5 and 1/T = 5.The relative difference is always below 1%.

TABLE I .
and II respectively.fRG* parameters used in Section V.The additional 24 bosonic patching points in Nq are distributed around k = (π, π).N k is the number of points in the momentum integration of the fermionic bubble.The frequency ranges of the vertex and vertex asymptotics are proportional to the number of positive fermionic frequencies N f + of the low-frequency object with three dependencies.Due to computational limits, the calculations for U > 2 are not converged with respect to the number of loops N and self-energy iterations NΣ-iter.

TABLE II .
fRG* parameters used in Section VI.The calculations are converged with respect to Nq, N k , N and NΣ-iter.