Multiloop functional renormalization group for the two-dimensional Hubbard model: Loop convergence of the response functions

We present a functional renormalization group (fRG) study of the two dimensional Hubbard model, performed with an algorithmic implementation which lifts some of the common approximations made in fRG calculations. In particular, in our fRG flow; (i) we take explicitly into account the momentum and the frequency dependence of the vertex functions; (ii) we include the feedback effect of the self-energy; (iii) we implement the recently introduced multiloop extension which allows us to sum up {\emph{all}} the diagrams of the parquet approximation with their exact weight. Due to its iterative structure based on successive one-loop computations, the loop convergence of the fRG results can be obtained with an affordable numerical effort. In particular, focusing on the analysis of the physical response functions, we show that the results become {\emph{independent}} from the chosen cutoff scheme and from the way the fRG susceptibilities are computed, i.e., either through flowing couplings to external fields, or through a"post-processing"contraction of the interaction vertex at the end of the flow. The presented substantial refinement of fRG-based computation schemes paves a promising route towards future quantitative fRG analyses of more challenging systems and/or parameter regimes.


Flow equations for the response functions 7
3 Numerical implementation 11 3.1 Full frequency and momentum parametrization 11 3.1.1 Truncated Unity fRG 12 3.1.2 Dynamical fRG 13 3.1.3 Flow equations for the TU-dynamical fRG 13 3.1.4 Calculation of the fermionic particle-hole and particle-particle excitation 15 3.1.5 Diagrammatic and lattice related symmetries

Introduction
Over the last two decades, functional renormalization group (fRG) methods have been broadly used for analyzing two-dimensional (2D) lattice electron systems (for reviews, see Refs. [1,2]). The main advantage of the fRG lies in the exploration of the leading low-energy correlations and instabilities towards long-range ordered states, similar to what has been investigated earlier for one-dimensional systems [3,4]. However, in one dimension, other methods like Bethe-Ansatz, bosonization [5,6] and DMRG [7] exist, which are for certain aspects more controlled. Hence, assessing the precision of RG methods in one-dimensional systems was not really in the foreground. The situation evidently changes for two-and three-dimensional systems, where the specific simplifications associated to the peculiar one-dimensional geometry are not applicable. At the same time, spatial correlations in 2D are strong enough to induce qualitative corrections [8,9] with respect to another class of rigorous many-body approaches, such as the Dynamical Mean-Field Theory (DMFT) [10][11][12] which allows one to include all purely local dynamical correlations. In fact, due to the intrinsic complexity of the many-electron problem in 2D, the development of unbiased quantitative methods applicable to a wide energy range from electronic structures on the scale of a few eV down to, e.g., ground state ordering in the (sub-)meV region, is still on the wishlist. This goal has motivated, in the last decade, the development of several algorithmic schemes for treating electronic correlations in 2D from different perspectives [1,13,14]. In this context, the fRG has already unveiled quite promising features: The fRG has the potential of resolving band structures and Fermi surface details and to treat competing orders on low energy scales in a rather unbiased way, since it does not require preliminary assumptions about dominating scattering channels. Recent applications range from studies of cuprate high-T c superconductors [15][16][17][18] over iron superconductors [2,19] to few-layer graphene systems [20,21], to cite a few.
We also note that, while the current applicability of the fRG is generally restricted to the weak to intermediate coupling regimes, its combination [22,23] with the DMFT might allow one, in the future, to access much more strongly correlated parameter regions, including the ones in proximity of the Mott-Hubbard metal-insulator transition. This is achieved by constructing a fRG flow starting from the DMFT solution of the considered lattice problem to the exact solution, i.e., in practice, using the DMFT to determine the initial conditions for the fRG flow [22]. Similarly to other diagrammatic extensions [14] of DMFT, such as the Dynamical Vertex Approximation (DΓA) [24] or the Dual Fermion [25] approach, one might work either with the physical degrees of freedom (as in the so-called DMF 2 RG [22]) or in the space of auxiliary (dual) fermions [26], introduced by means of a suitable [14,25] Hubbard-Stratonovich transformation.
Yet, what is hitherto missing is a thorough analysis of the quantitative reliability of the fRG for a well-defined test case. More precisely this would require to clarify how much the fRG results, going beyond the correct estimation of general physical trends, depend on the approximations inherent in the used fRG scheme. This study within the fRG would then also provide a solid basis for future comparisons with other numerical techniques.
The mentioned approximations can be grouped in three categories: (i) Momentum/frequency discretization: As the fRG algorithm typically exploits the flow of vertex functions that depend continuously on multiple momenta and frequencies, various approximations are performed to mitigate numerical and memory costs. Early on, N -patch discretizations of the momentum dependencies through the Brillouin zone were used. Later, it was noticed that channel-decompositions in conjunction with form factor expansions [27][28][29] lead to physically appealing approximations featuring advantageous momentum resolution and numerical performance [30]. Clever prescriptions for the treatment of the high-frequency tails of the vertex function have been devised [31][32][33] which are also used in this work.
(ii) Self-energy feedback: In many applications of the fRG the self-energy and its feedback on the flow of the n-particle (n > 1) vertex functions has not been accounted for. While there are arguments that the self-energy may be important mainly when the interactions are close to a flow to strong coupling (see Appendix in Ref. [34]), more quantitative results should overcome this deficit. In fact, neglecting the self-energy feedback was mainly motivated by the disregarded frequency dependence of the interactions in earlier fRG studies: Within a static treatment the self-energy lacks the effects of quasiparticle degradation, so that its inclusion became less important. Within the current frequency-dependent fRG treatments, the self-energy feedback can be included in a meaningful way. A number of works have already investigated the self-energy effects in the flows to strong coupling in Hubbard-type models [28,[35][36][37][38][39][40][41][42][43], mainly exploring the quantitative effects, besides signatures of pseudogap openings [38,39] and non-Fermi liquid behavior [28] in particular cases. (iii) Truncation of the flow equation hierarchy: Finally, one should also consider the truncation of the hierarchy of flow equations for the n-point one-particle irreducible (1PI) vertex functions. This is usually done at "level-II" as defined in Ref. [1], also referred to as one-loop (1 ) approximation, i.e., the 1PI six-point vertex is set to zero. Due to this truncation, the final result of an fRG flow might depend -to a certain degree-on the cutoff scheme adopted for the calculation.
In this perspective, it was noticed by Katanin [44] that replacing the so-called singlescale propagator in the loops on the r.h.s. of the flow equation for the four-point vertex by a scale-derivative of the full Green's function allows this scheme to become equivalent to oneparticle self-consistent (a.k.a. mean-field) theories in reduced models, and then to go beyond such self-consistent approximations in more general models. Another significant comparison can be made with the parquet-based approaches [45,46], such as the parquet approximation (PA) [32,33,[47][48][49]. The latter represents the "lowest order" solution of the parquet equations, where the two-particle irreducible vertex is approximated by the bare interaction. In fact, although the diagrams summed in the 1 truncation of the fRG are topologically the same as in the PA, the way the single contributions are generated during the flow leads in general to differences with respect to the PA [33,50]. This is due to some internal-line combinations, e.g., in particle-hole corrections to the particle-particle channel, which are suppressed by the cutoff functions attached to the propagators and not fully reconstructed during the flow because of the truncation. A quantitative analysis of this effect has been performed for the single impurity Anderson model in Ref. [33]. These differences are absent for single-channel summations (e.g. RPA), but could lead to more pronounced quantitative errors in presence of channel coupling, e.g., in the generation of superconducting pairing through spin fluctuations. Furthermore, while the Mermin-Wagner theorem is fulfilled within the PA [51], it is typically violated by 1 fRG calculations. First steps to remedy this shortcoming were undertaken in various works [42,52,53], but only recently a comprehensive path of how the PA contributions can be recovered in full extent was presented within the multiloop extension of the fRG (mfRG) [54][55][56]. The mfRG flow equations incorporate all contributions of the six-point vertex that complement the derivative of diagrams already part of the 1 flow, as organized by their loop structure. A key insight in this approach is that the higher-loop contributions can be generated by computing 1 flows for scale-differentiated vertices, with an effort growing only linearly with the loop order that is fully kept. The multiloop corrections stabilize the flow by enabling full screening of competing two-particle channels, ultimately recovering the selfconsistent structure of the PA. As the PA corresponds to a well-defined subset of diagrams, a converged mfRG flow able to reproduce the PA is by construction independent of the adopted cutoff.
In this paper, we present a fRG study of the 2D Hubbard model performed with an algorithm combining the most recent progress on all three approximation levels. We use (i) the so-called "truncated unity" fRG [30] (TUfRG) formalism to describe the momentum dependence of the vertex and, in addition, keep the full frequency dependence as a function of three independent frequencies. Differently from the approach adopted in Ref. [43], we employ a refined scheme to treat the high-frequency asymptotics [33] that allows us to reduce the numerical effort considerably. Within this scheme, we can consistently include (ii) the (frequency-dependent) self-energy feedback in our fRG flow equations. Finally, we present (iii) first data for the 2D Hubbard model computed with the multiloop extension proposed by Kugler and von Delft [54]. In this context, we have also generalized the multiloop formalism to compute the flow of the response functions, and illustrated the loop convergence of the fRG results for the 2D Hubbard model. In particular, we show that including up to 8 loops in the fRG flow yields a clear convergence of the data with the loop order and the final results are independent of the cutoff. This represents an important check and illustrates that fRG flows can be brought in quantitative control for 2D problems. Finally, our multiloop analysis of the response functions demonstrates that the two different ways to compute susceptibilities in the fRG, either by tracking the renormalization group flow of the couplings to external fields [1] or by contracting the final interaction vertex (see, e.g., Ref. 22), converge to the same value with increasing loop order. This confirms that the output of this improved fRG scheme can indeed be trusted on a quantitative level.
The paper is organized as follows: The formalism and theory of the linear response functions and their computation by mfRG flow equations are introduced in Section 2. In Section 3 we present the actual implementation scheme for the full momentum-and frequencydependent fRG. In Section 4 we show the results for the 2D Hubbard model, with a detailed analysis of the effects of the different approximation levels and in particular of the convergence with the loop order. A conclusion and outlook is provided in Section 5.

Theory and formalism 2.1 Definitions and formalism
In this section we provide the definitions of the linear response functions to an external field, before describing their computation with the fRG. We focus on correlation functions of fermionic bilinears. In particular, in a time-space translational-invariant system, we consider the charge (density) and spin (magnetic) bilinears, both charge invariant, and the non-charge invariant pairing (superconducting) bilinears where ψ andψ represent the Grassman variables and p (q) a fermionic (bosonic) quadrimomentum p = {iν o , p} (q = {iω l , q}). The integral includes a summation over the Matsubara frequencies (iν o ), normalized by the inverse temperature β, and an integral over the first Brillouine Zone normalized by its volume V BZ . The function f n (p, q), determines the momentum and frequency structure of the bilinears in the different physical channels. In the present case we restrict ourselves to a static external source field, such that the function f n (p, q) = f n (p) acquires only a momentum dependence, whose structure is specified by the subscript n and explicitly shown in Table 1 (in the present work we will mostly focus on the s-as well as d-wave momentum structure). Note that, when using a different frequency-momentum notation, centered in the center of mass of the scattering process (see "symmetrized" notation in Appendix A), one should account for an additional shift of the momentum dependence p by means of the momentum transfer q. After a reshift of the operators in Eq. (1) with respect to their average value ρ n d/m → ρ n d/m − ρ n d/m , we can now define the correlation functions of these bilinears in the three channels In linear response theory, these correlation functions correspond to the physical susceptibilities in the corresponding channels. Divergences in χ nn η (q), with η = {d, m, sc}, indicate spontaneous ordering tendencies or instabilities of the system. The above definition encodes not only the real-space pattern or wavevector with which the system starts ordering, but also the symmetry of the order parameter associated to the instability. In the 2D Hubbard model study presented here (see Section 4) we detect various response functions growing considerably towards low T , such as the spin-density wave (SDW) response, characterized by the isotropic s-wave magnetic susceptibility at q = (π, π), as well as s-and d-wave pairing response functions at q = (0, 0) and Pomeranchuk instabilities [57]. Inserting Eq. (1) or Eq. (2) into Eq. (3), the susceptibilities appear as two-particle Green's functions. In particular, they can be determined from the two-particle vertex γ 4 by where σ 0/3 represent the Pauli matrices (σ 0 = 1) and we made use of the spin conservation. Eqs. (4) can be considerably simplified by making use of the SU(2) symmetry. The "bare bubbles" Π η appearing in (4) read By exploiting the SU(2) symmetry, we can drop the spin dependencies for the bare bubbles. In presence of the above symmetries, we can introduce the following definitions for (spin-independent) channels of the two-particle vertex γ 4,d (q, p, p ) = 1 2 σ,σ γ 4;σσσ σ (p, p + q, p + q, p ) (7a) with the Levi-Civita symbol. The resulting spin-independent expression of the physical susceptibilities reads We conclude this section by recalling the definition of the so-called fermion-boson vertex [58], which, for the considered symmetries, reads Similarly to the susceptibility, one can rewrite Eqs. (9a) and (9b) in a form where the twoparticle vertex γ 4,η appears explicitly γ n 3,η (q, p) = f n (p) + dp f n (p )γ 4,η (q, p, p )Π η (q, p , p ) , where, because of the SU(2) symmetry, we dropped the spin dependence of the fermion-boson vertices γ n 3,sc;↓↑ = γ n 3,sc;↑↓ = γ n 3,sc . (11b)

Flow equations for the response functions
In this section we derive the mfRG [54] flow equations of the response functions and discuss the improvement with respect to the 1 version [1]. Note that one could analogously provide a formal analytical derivation [59] in the same spirit of the approach used in Ref. [56]. In the following we provide the main steps of the derivation in the 1PI formulation [1,60] (see also Ref. [57] for the Wick-ordered formulation), for the details we refer to Appendix B. Following the review of Metzner et al. [1], we introduce the coupling of the density operators in Eqs.
We note that, although J n η appears as a functional dependence in our derivation, it is not an integration variable since our system is fully fermionic (for an fRG formulation of coupled fermion-boson systems, see Refs. [1,[61][62][63]).
By expanding the scale-dependent effective action Γ Λ in powers of the fermionic fields, as well as of the external bosonic source field, we obtain Note that the index x = {σ, k} combines the spin index σ and the fermionic quadrivector k = (iν l , k) (here we disregard additional quantum dependencies, e.g., orbital), while y = {n, q} refers to the momentum structure of the coupling to the bilinears, n, and to the bosonic quadrivector q = {iω l , q} . In Eq. (13) the first term on the r.h.s. represents the expansion of the effective action in absence of external field (see Section 3), while the functional derivatives in the following terms represent the Λ-dependent susceptibility and the fermion-boson vertex in the different channels. Taking the derivative with respect to the scale parameter Λ (see Appendix B), yields the following flow equations for the susceptibility and fermion-boson vertex (assuming SU(2) symmetry and momentum-frequency as well as spin conservation) and respectively where represents the single-scale propagator. The functionγ 4 , differently from the (fermionic) twoparticle vertex γ 4 , represents a mixed bosonic-fermionic vertex, i.e., with two bosonic and two fermionic legs where we summed over its spin dependences while the spin-independent form for γ 5 used in Eqs. (15) reads The conventional approximations [1,57,60] disregard the first terms on the r.h.s. of Eqs. (14) and (15). This 1 approximation is consistent with the corresponding approximation of γ Λ 4 (see Appendix C) and justified in the weak-coupling regime. Using the notation of Refs. [44,54], one can rewrite the 1 approximation of Eqs. (14) and (15) in a more concise tensor-formχ We here introduced also the subscript ph and pp indicating the diagrammatic channels that will be referred to in Sec. 3. So far we pinpointed two possible ways to compute the susceptibility and fermion-boson vertex from an fRG calculation: (i) Solving Eqs. (14) and (15) alongside the ones for Σ and γ 4 (at the same level of approximation), and (ii) by means of Eqs. (8) and (10) at the end of the fRG flow, using Σ Λ final and γ Λ final 4 , later referred to as "post-processing". These two procedures are non-equivalent in the presence of approximations, e.g., if one restricts oneself to the 1 level. This leads to an ambiguity in practical implementations of the fRG. In fact, as shown in Appendix D, the two results deviate at O((γ Λ 4 ) 2 ) for the 1 case (for a larger number of loops the deviations occur at higher orders in the effective interaction γ Λ 4 ). In order to solve this ambiguity we note that the exact relations (8) and (10) are fulfilled in the PA. At the same time, the recently introduced multiloop extension allows one to sum up all parquet diagrams. Hence, generalizing the multiloop flow to the computation of the response functions recovers the equivalence of the two procedures.
In order to derive the mfRG equations for the response functions, we first recall the channel-decomposition of the two-particle vertex as known from the parquet formalism. The latter divides γ 4 in the two-particle reducible vertex φ (all diagrams that can be divided into two separate ones by removing two internal fermionic propagators) and the two-particle irreducible vertex I (which can be not be divided). Depending on the direction of the propagation lines the diagrams are reducible in either parallel, longitudinal antiparallel or transverse antiparallel, corresponding to the particle-particle, particle-hole, and particle-hole crossed channel, respectively. Besides this diagrammatic channel decomposition, there is also a distinct physical channel decomposition that identifies the components η = {d, m, sc} and which we will use in the following. Inserting this decomposition into the flow equation for the twoparticle vertex, we obtain While the usual diagrammatic channel decomposition [64] leads to simple expressions for the two-particle irreducible vertex I Λ η , the latter assumes a more complicated form in the physical channels where we approximated the fully two-particle irreducible vertex by its first-order contribution in the interaction ∼ U , which is known as PA.
We now derive the mfRG flow equations for the response functions, which mimic the effect of the mixed fermion-boson verticesγ Λ 4 and γ Λ 5 in the exact flow Eqs. (14) and (15). First, one performs the so-called Katanin substitution [44] in the 1 flow equations (19). One observes that all differentiated lines in these flow equations come fromΠ Λ η . Secondly, differentiated lines from the other channels are contained in the higher-loop terms of the expansion Using the channel decomposition (21), we can directly write down the 2 correction to the flow of the fermion-boson vertex, which accounts for the leading-order diagrams of the effective interaction and stem from γ Λ 5 in Eq. (15) On the three-and higher-loop level, we can now useİ in an analogous way. In addition, we have to consider the vertex corrections to the right of the differentiated lines, yieldinġ Considering the 1 flow equation of the susceptibility (19a), we see that the fermion-boson vertices provide vertex corrections on both sides of the differentiated lines inΠ Λ η . Hence, for all higher-loop corrections we can simply connectİ Λ( ) η to both fermion-boson vertices, thereby raising the loop order by two. We obtainχ Λ(2) η = 0, as well aṡ For a schematic representation of Eqs. (25) and (26)   two-particle vertex (see Section 3.2) allow us to sum up all differentiated parquet diagrams of γ Λ 3 and χ Λ . As a consequence, the aforementioned two ways of computing the response functions within the fRG become equivalent. We finally note that for a consistent fRG scheme, it is important to adopt the same level of approximation (truncating the sums in Eq. (23a) to a certain finite -loop level) for all flowing quantities.

Numerical implementation 3.1 Full frequency and momentum parametrization
In order to illustrate the fRG algorithm adopted in the present work, let us start from the flow equations for the 1PI fermionic vertex in the 1 fRG approximation. In the following, the SU(2), spin conserving symmetry will be always assumed. Exploiting this symmetry, the self-energy and two-particle fermionic vertices can be written as where the fourth argument of γ 4 is determined by k 4 = k 1 + k 3 − k 2 in a momentum and energy conserving system. The spin-independent flow equation for the self-energy readṡ where S Λ (p) represents the single-scale propagator specified in Eq. (16). We formulate the flow equation for γ 4 in the channel decomposed form suggested by Husemann and Salmhofer [27] γ where the diagrammatic channel index r = {pp, ph, ph} distinguishes between particle-particle, particle-hole and particle-hole exchange diagrams, and the first dependence of the functions T Λ r refers to the bosonic four-momentum transfer in the internal loop of their corresponding Each of the above equations depends, besides the aforementioned bosonic transfer dependence (k 1 + k 3 , k 2 − k 1 and k 3 − k 2 ), on two fermionic dependencies. Such mixed 'bosonic-fermonic' notation, referred to as 'non-symmetrized' notation, has been substituted in some work (e.g., in Ref. [30]) by a different notation where the dependencies of the four fermionic propagators involved in the scattering process have been chosen symmetrically with respect to the bosonic four-momentum transfer. This symmetrized notation simplifies the implementation of the symmetries exploited in the fRG code (see Appendix F and Ref. [30]) but leads to less compact flow equations. The equation (31) generates the two-particle reducible vertices T r =φ r of the diagrammatic parquet decomposition The two-particle fermionic vertex can be reconstructed by using Eq. (32). The use of a mixed 'bosonic-fermonic' notation allows us to identify the bosonic transfer four-momentum as the strongest dependence, while the two fermionic dependencies can be treated with controllable approximations. In the following we illustrated two efficient ways to simplify the treatment of both momentum and frequency dependencies.

Truncated Unity fRG
The approximation for the fermionic momentum dependencies in TUfRG [30] is done by the expansion of the fermionic momentum dependencies in form factors, illustrated here for the while the expansion of the φ ph and φ ph analogously defines D n,n (q) and C n,n (q). Following the conventions introduced in previous works [27,28,30,40,65,66], we choose the form factors such that they correspond to a specific shell of neighbors in the real space lattice. The unity inserted in the flow equations contains a complete basis set of form factors Converged results can be obtained already with a small set of form factors [30]. For a fast convergence it is convenient to include one shell after another, starting from the constant local form factor and increasing the distance of neighbors taken into account. The form factors used in this paper are listed in Table 1.
A major difficulty in this approach is the feedback of the different channels into each other. In addition to the dressing of the objects by the form factors, the translation of the notation in momentum and frequency from one to another channel has to be considered. Computationally time consuming integrations in momentum space can be avoided by Fourier transformation and evaluation in real space. Furthermore the expression of the projection in terms of a matrix multiplication allows for the precalculation of the projection matrices which can be found in the Appendix F.

Dynamical fRG
In frequency space, we adopt the simplifications proposed in Refs. [32,33]. For all systems with an instantaneous microscopic interaction one can use diagrammatic arguments to prove that, in the high-frequency regime, the fermionic two-particle vertex exhibits a simplified asymptotic structure. In this region one can reduce the three-dimensional frequency dependence of γ 4 using functions with a simplified parametric dependence. It is straightforward to see that, sending all three frequencies to infinity, γ 4 reduces to the instantaneous microscopic interaction, which in the present case is represented by the Hubbard on-site U . The contribution of the reducible vertices φ r to γ 4 becomes non-negligible if the bosonic frequency transfer is kept finite, while sending the two secondary fermionic frequencies to infinity. This contribution, depending on a single bosonic frequency transfer in a given channel r, is denoted by K 1,r (iω l , q). For models with an instantaneous and local microscopic interaction, one observes that the momentum dependencies disappear alongside the frequency dependencies when performing such limits. In the limit where just one fermionic frequency is sent to infinity, the vertex φ r can be parametrized by the function K 2,r (iω l , iν o , q, k) + K 1,r (iω l , q). By subtracting the asymptotic functions from the full object φ r we obtain the so-called [33],"rest-function" R(iω l , iν o , iν o , q, k, k ) which decays to zero within a small frequency box. The parametrization in terms of K 1/2 allows us to reduce the numerical cost of computing and storing the fermionic two-particle vertices. In fact, for any of the three channels, we calculate the fRG flow of the three-frequency dependent function R on a small low-frequency region and add the information on the high frequencies by computing the flow of the functions K 1 and K 2 which are numerically less demanding. The full two-particle reducible vertex φ r is then recovered by whereK 2,r can be obtained from K 2,r by exploiting the time reversal symmetry (see Appendix A.3).

Flow equations for the TU-dynamical fRG
Finally, applying the aforementioned projection on the form-factor basis we can write matrixlike 1 fRG flow equations for the self-energy, the two-particle vertex, the fermion-boson vertex and the susceptibility: where the multiplication of bold symbols has here to be understood as matrix multiplications with respect to the form factors. For a schematic visualization of the practical implementation of these equations, see Fig. 2. We note that, in order to derive Eqs. (36), we inserted the unity (34), truncated to a finite number of form factors, in Eqs. (19) as well as in (31). The full vertex γ 4,r , with r = {P, D, C} represents the fermionic two-particle vertex in the channel-specific mixed 'bosonic-fermionic' notations, while γ 4,η with η = {sc, d, m} is given by The 1 -fRG flow consists in integrating the coupled differential equations in (36) with the following initial conditions: Finally,Π S,η relative to the particle-hole η = {d/m(ph)} and to the particle-particle channels η = {sc(pp)}, are defined aṡ Self-energy e¡ort: Projection Figure 2: Schematic code structure specifying the array sizes and the numerical effort of the single steps. I f denotes the number of elements of the object f . N ν is the number of fermionic frequencies of the rest function, N q the number of bosonic momentum patches, N F F the number of form factors, N R FFT and N ν int the number of frequencies over which the internal fermionic bubble is integrated. The symmetries reduce the total number of elements I f to I n f independent elements which have to be calculated and to I s f which can be obtained by using symmetry relations. The arrows indicate the feedback of the different parts, namely the two-particle fermionic vertices (red), fermion-boson vertex (yellow), and the self-energy (green).
whereΠ Λ S,η (q, k) is defined in Eq. (20). In order to perform the momentum integration in Eqs. (39) we adopt a strategy which, exploiting the convolution theorem, represents a numerically convenient alternative to the use of adaptive integration algorithms. The latter is described in the following section.
3.1.4 Calculation of the fermionic particle-hole and particle-particle excitation We here present a numerically convenient way of calculating the fermionic particle-hole and particle-particle bubbles in the flow equations of the vertex (36), defined in Eqs. (39). Since the integral over momenta is very sensitive on the momentum mesh resolution near the Fermi surface and a refined adaptive integration is computationally time consuming, we rewrote the integrals in such a way to use the convolution theorem. The Green's function can then be transformed via the Fast-Fourier-Transform (FFT) to real space, where the real-space expression of the form factors is provided in Table 1. After some algebraic steps, we find an expression without momentum integratioṅ where F f (R) (k) is the Fourier transform which can be determined by using FFT-methods and the weight W n,n (R) is defined as The infinite sum of the lattice points in Eqs. (40b), (40a), and (41) is restricted by the finite range of the form factors for a specific truncation. For instance the sum in Eq. (41) is limited to the maximal shell taken into account by the form factors. Hence, the weight has a nonzero contribution only inside a shell twice as large the maximal shell of the form factors and therefore the sum in Eq. (40a) can be constrained to twice the distance of the maximal form factor shell. The momentum and real space grid for the Fourier transformations needed in the bubbles has to be chosen fine enough, especially at low temperatures. The convergence in terms of FFT-grid points N R FFT has to be checked separately from the bosonic momentum grid of the vertex. Recent works using the TUfRG [30,66] have demonstrated that, if needed, both low temperatures and high wavevector resolutions can be achieved by means of an adaptive integration scheme.

Diagrammatic and lattice related symmetries
Further numerical simplifications come from the extensive use of symmetries related to diagrammatic arguments and lattice-specific properties, which can be found in Appendix A.

The mfRG implementation
The mfRG flow introduced in Ref. [54] ameliorates the approximation induced by the truncation of the fRG hierarchy of flow equations as it incorporates all contributions from the six-point vertex γ 6 that can be computed at the same cost as the 1 flow considered so far. In fact, it includes all contributions coming from γ 6 that can be computed in an iterative 1 construction of four-point objects; hence, the numerical effort grows only linearly in the number of loops retained. It has been shown [54] that the multiloop prescription fully sums up all parquet diagrams. This gives rise to a number of advantageous properties, the most important of which are (i) that the multiloop corrections restore the independence on the choice of regulator, and (ii) that the multiloop flow fully accounts for the interplay between different two-particle channels and thus hampers spurious vertex divergences coming from ladder diagrams in the individual channels.
Let us briefly recall the multiloop vertex flow employing the same line of arguments as used for the flow of the response functions in Section 2.2. We consider the reducible vertices in the physical channels φ η={sc,d,m} . At first, one performs the Katanin substitution [44] and finds that, for every channel φ Λ η , all differentiated lines come fromΠ Λ η . Differentiated lines from the other channels are contained in higher-order terms of the loop expansion Using the channel decomposition (21), one has the two-loop correctioṅ where, according to Eq. (22),İ Λ,( ) η can be determined from theφ Λ,( ) η of the complementary channels η = η. All higher-loop terms are obtained in a similar fashion where one additionally accounts for vertex corrections to both sides ofİ In Ref. [54], it has further been pointed out that corrections to the self-energy flow (29) are necessary in order to generate all differentiated diagrams of the parquet self-energy. These corrections are included in the central part of the vertex flow withΣ given by Eq. (36a) and where the central part (see Eq. (45)) for the differentiated reducible verticesφ r={P,C,D} = {Ṗ,Ḋ,Ċ} is defined by and Table 1: Local and first nearest-neighbor form factors both in momentum and real space presentation. For each calculation we specify which form factors are used. A pure s-wave calculation restricts to the first line corresponding to the local form factor, the d-wave accounts for the first two nearest neighbors form factors, and a calculation with all nearest neighbors form factors includes all five form factors shown here.

Numerical results
In this section we show fRG numerical results obtained with the formalism and code described in the previous sections. After introducing our test system, namely the 2D Hubbard model at half-filling, we will test our full momentum-frequency resolved fRG implementation, together with the inclusion of the self-energy feedback, and study the effect of including multiloop corrections to the 1 approximated flow equations. If not specified differently, we will make use of a "smooth" frequency-dependent regulator throughout this work: where G 0 specifies the non-interacting Green's function of the 2D Hubbard model. The fRG scheme associated to such a regulator is referred to as Ω-flow [28].

2D Hubbard model at half filling as test system
As test model we consider the single-band two-dimensional (2D) Hubbard model on the square lattice. Its Hamiltonian readŝ whereĉ ( †) iσ annihilates (creates) an electron with spin σ at the lattice site R i (n iσ =ĉ † iσĉ iσ ), t is the hopping amplitude for electrons between neighboring sites, µ the chemical potential and U > 0 the repulsive on-site Coulomb interaction. In the present study, we consider U = 2t, µ = U/2 and different temperature regimes. Since the present model has been extensively studied in the theoretical literature (see, e.g., Refs. [10,13,14,45,[67][68][69][70]) as well as in fRG (for a review, see Ref. [1]), it constitutes a reference system to test our novel fRG implementation. Furthermore, the 2D Hubbard model constitutes a delicate case in the context of the Mermin-Wagner theorem [71], which prevents the onset of the antiferromagnetic ordering at finite temperature. Whereas the 1 fRG results exhibit a pseudocritical Néel temperature (T pc ), the inclusion of the multiloop corrections to the standard fRG flow should, from a theoretical perspective, recover the parquet solution, which is known to fulfill the Mermin-Wagner theorem [72]. Therefore, we expect T pc to be moved down to 0 in the (converged) multiloop fRG scheme. Despite the rich phase diagram of the 2D Hubbard model out of particle-hole symmetry, we restrict this study to the half-filled particle-hole symmetric case, in order to reduce the numerical efforts.
Let us stress that the bosonic momentum discretization of the first Brillouin zone (BZ) has been chosen such that one obtains a uniform grid along the x-and y-directions. This represents, though, not the unique choice of resolving the reciprocal space and one could adopt some sophisticated "patching" schemes [43], which should be accounted in future optimization of our code.

Convergence and stability study on the TUfRG-implementation
In the previous section 3.1, we presented an efficient parameterization of the vertex which combines the TUfRG scheme [30] for treating momenta and the dynamical fRG implementation proposed in Ref. [33]. In order to illustrate the efficiency of such merge, we have performed a convergence study of the (dominant) antiferromagnetic (AF) susceptibility, (χ AF = χ 00 m (iω l = 0, q = (π, π)), by means of Eq. (36f), as a function of the number of Matsubara frequencies, momenta and form factors, used in our algorithm. The convergence tests have been performed at temperatures T = 0.25t T pc and T = 0.125t ∼ T pc . Let us first consider the convergence in the number of fermionic frequencies N ν at which the low-frequency structure of the rest-function R is captured. For T = 0.25t in Fig. 3 (left panel) one observes that the susceptibility does not show significant changes as a function of N ν . In fact, it is known that, in weakly correlated electron systems, the frequency dependence of the vertex is less important because of power counting arguments [57,73] and as shown Figure 4: Inverse (1 ) AF susceptibility at q = (π, π) as a function of temperature, for U = 2t. Only local (s-wave) form factors are used, but including the nearest-neighbor form factors does not change the results within the accuracy. Besides the curve obtained using the full TU dynamical fRG scheme (i) (blue dots, "full fRG"), different approximations are shown: approximation (ii) (green diamonds, "no Σ"), (iii) (red squares, "no ω dep") and (iv) (yellow triangles, "no ω dep, no Σ").
by numerics for small numbers of fermionic Matsubara frequencies, e.g., in Ref. [41]. At T = 0.125t the convergence with respect to N ν is slower. According to our tolerance of 1% we get convergence at N ν = 8. Secondly, we analyze the dependence of the AF-susceptibility on the number of bosonic patching points, N q , as shown in Fig. 3. The data at T = 0.25t are already converged at N q = 64, while at T = 0.125t we need N q = 256. In the latter case, one sees that the convergence is more sensitive to N q than to N ν . This can be ascribed to the presence of a finite pseudocritical temperature since for T → T pc the AF fluctuations become long-ranged, requiring an increasingly finer momentum resolution. At the same time, the size of the objects to handle is only linear in N q while it is expected to scale up to cubic in N ν , depending on the quantity considered (see Fig. 2). Moreover, the number of independent momentum patching points can be substantially reduced by exploiting point-group symmetries of the lattice.
Last but not least, we have also verified that, for all T considered, the AF response function is fully converged with respect to the number of form factors (not shown).

Effects of different approximations
In our fRG scheme, we can choose different approximation levels regarding the treatment of the frequency dependence of the interactions and the self-energy. This allows us to gain a better understanding of the interplay of the different interaction channels and the role of the self-energy.
Here we define four approximation levels (i) to (iv) with decreasing rigor. Approximation (i) represents the fRG treatment described in Sec. (3) which merges the TUfRG scheme with an efficient inclusion of the vertex dynamics; (ii) denotes the flow with a frequency-dependent effective interaction but without the flow and feedback of the self-energy; (iii) is the frequencyindependent (static) approximation for the effective interaction and the self-energy, in which the fermion-fermion, fermion-boson and boson-boson vertices are approximated by their value at zero frequency; and (iv) combines the neglect of the self-energy feedback with a static approximation for the vertices.
Approximation (iv) has been the standard one adopted in many previous works, as those reviewed in Ref. [1]. Various other fRG works have already explored the changes occurring by using better approximations like (i) to (iii) introduced above. Earlier studies of the self-energy without explicit frequency dependence of the effective interaction pointed to the possibility of non-Fermi liquid behavior [38,74]. Later, channel-decomposed fRG [28,40] and N -patch fRG [41] were used to explore the effects of a frequency-dependent effective interaction and of the self-energy feedback. In the following, we rediscover some of their findings, with a more refined momentum-and frequency-dependent self-energy. Eberlein [42] used a channeldecomposed description of the interaction where each exchange propagator was allowed to depend on one bosonic frequency. He found that in the presence of antiferromagnetic hot spots on the Fermi surface, antiferromagnetic fluctuations lead to a flattening of the Fermi surface and increase the critical scales. Most recently, Vilardi et al. [43] presented a refined 1 study of the role of the various frequency structures in the interaction, parametrized by three frequencies, albeit with a reduced set of form factors. They argued that a one-frequency parametrization can in some cases lead to spurious instabilities. Our study differs from this work by the ability of taking into account more form factors, using a more economic description of the higher frequencies, and by implementing the multiloop corrections.
In Fig. 4 we show how differently the approximations affect the results for the AF susceptibility. More precisely, we plot the inverse AF susceptibility which decreases quite linearly, i.e., Curie-Weiss-like, upon lowering T . The intersection of the curve with the abscissa marks the pseudocritical temperature which, violating the Mermin-Wagner theorem, assumes a finite value in the 1 fRG scheme. One can observe that the full TU-dynamic fRG approach (i) leads to larger inverse AF susceptibilities, or smaller χ AF , than the other three approximations, shifting T pc to a smaller value.
Let us first compare the full calculation (i) with the calculation without self-energy but frequency-dependent interactions (ii). It is to be expected that the self-energy renormalizes the leading vertices and therefore also susceptibilities, as has also been observed in fRG studies [28,43]. This explains why the calculations without self-energy flow diverge at higher T pc with respect to scheme (i).
The flow variants with static interactions (iii) and (iv) differ only slightly. Compared to the fRG flow using scheme (ii), the AF tendencies in these static flows are somewhat weaker as their suppression by particle-particle processes increases when the pairing channel is approximated by its static part, for which it assumes the maximum value. The downwardshift in the inverse AF susceptibility from (iii) and (iv) to (ii) with the inclusion of the frequency dependence of the couplings is however overcompensated by the inclusion of the dynamical self-energy in (i).
Finally, we consider the pseudocritical temperature and the AF susceptibility for the combined approximation of no self-energy and no frequency dependence (iv). Without the screening effect of the self-energy, the pseudocritical temperature increases a little bit more with respect to the static approximation (iii). This has been already observed in Ref. [ 41]. The small difference may come from the real part of the self-energy that can be understood as upward-  renormalization of the hopping parameter, or equivalently a downward-renormalization of the density of states. This is consistent with the self-energy shown below in Fig. 6. For a detailed discussion on the pseudocritical temperatures on a wider range of parameters, we refer the reader to Ref. [40].

Computation of the self-energy
As already implied above, the implementation presented in Sec. 3.1 allows one to compute a frequency and momentum dependence of self-energy during the flow according to Eq. (36a). In Figs. 5 and 6, we present the results for the frequency-and momentum-dependence of the self-energy for different temperatures and momentum points. For the fermionic momentum patching we use the same momentum grid as for the bosonic transfer momentum of the vertex. In the results shown in Fig. 6 (left panel), we subtracted the Hartree contribution, which represents a rigid U/2 energy shift at half-filling. By looking at Fig. 5, we notice that the frequency dependence of the imaginary part of the self-energy is consistent with a Fermi-liquid, yet without any remarkable difference at different temperatures. As the slope of these curves determines the quasiparticle weight Z, we arrive at the conclusion that Z does not decrease steeply when we lower T towards the AF pseudocritical temperature, as Figure 7: AF susceptibility (upper panels) and fermion-boson vertex (lower panel) at (q = (π, π)) as a function of the number of loops, for U = 2t and T = 0.5t, 0.2t, 0.125t (from left to right). The susceptibility is evaluated at ω = 0 and the fermion-boson vertex at ω = 0 and ν = π/β. The blue line shows the behavior of the integrated Eq. (23a) up to = 8, while the green line the one obtained from the post-processed calculation by means of Eq. (8) for χ and of (10) for γ 3 ). The insets show the relative difference between the blue and the green lines, defined for the susceptibility as rel = (χ flow − χ post-proc )/χ =8 post-proc .
already observed in Refs. [41,42]. Figure 6 shows the momentum dependence of the real and imaginary part of the self-energy along a path in the first BZ defined by Γ = (0, 0), X = (0, π) and M = (π, π). The fermionic frequency is set to the first fermionic Matsubara frequency. The real part is positive at M and negative at Γ, while at X and Y it is zero. At lowest order, this momentum structure can be approximated by a positive nearest-neighbor hopping renormalization, which increases the bandwidth. The imaginary part of the self-energy shows two peaks around X and M/2= (π/2, π/2). This corresponds to a maximal scattering on the nested Fermi surface and minimal on the points Γ and M, which are at maximal distance from the Fermi surface. Note that this refers to the self-energy at small fixed imaginary frequency and not at real frequency equal to the excitation energy, i.e., this behavior does not contradict the typical behavior that the scattering rates for quasiparticles rise with distance from the Fermi surface.

Effect of the multiloop implementation
Let us now investigate the effect of including multiloop corrections to the flow equations of the susceptibility and the fermion-boson vertex as in Eq. (23a). As previously discussed, the inclusion of the multiloop corrections should allow us to recover the full derivative of Eq. (8) and (10) with respect to the scale parameter Λ. This means that the integration of the multiloop fRG flow equations should converge, by increasing the number of loops, to Eq. (8) and (10), as well as to the parquet equations for γ 4 and Σ as discussed in Ref. [54]. Although, in the half-filled case, the numerical effort is already reduced compared to the non-particle-hole symmetric situation, calculations for T < 0.5t are already quite demanding if a multiloop cycle is included. Therefore, the only calculations involving more than one form factor (i.e., s-wave) that will be presented here were performed at a rather high temperature of T = 0.5t. Despite this restriction, since the physics of the single band Hubbard model at half-filling is dominated by the AF fluctuations, the fRG results are already converged in the form factor numbers. Nevertheless, a meaningful part of the d-wave susceptibilities is still accessible, as it will be shown in the following, via the s-wave two-particle vertex.
In Fig. 7 we show the s-wave susceptibility χ (fermion-boson vertex γ 3 ) in the upper (lower) panels in the magnetic channel for iω l = 0 (iω l = 0 and iν o = π/β for γ 3 ) and q = (π, π) as a function of the number of loops considered in the mfRG calculation, for three selected temperatures T = {0.5t, 0.2t, 0.125t} (left to right). The blue lines show the value of χ and γ 3 calculated by the integration of Eq. (23a). On the other hand, the green lines show χ (γ 3 ) acquired at the end of the -loop fRG flow (Λ = Λ fin ) by means of Eq. (8) ((10)), where we inserted on the r.h.s γ Λ fin 4 and G Λ fin , referred to in Section 2.2 as "post-processed" method. In the present case, one sees how the convergence of the two lines is achieved after 8 for all temperatures presented. Thus, we have a dual convergence: as a function of the loop number and between two ways of computing the same quantity. Clearly, by decreasing the temperature and approaching the 1 fRG pseudocritical temperature (see Fig. 10), the antiferromagnetic (AF) susceptibility and γ 00 3,m (ω = 0, ν = π/β, q = (π, π)) = γ 3,AF increase and the green and blue lines for the two ways to compute the susceptibility exhibit the largest relative difference at = 1 of ∼ 25%. This difference decreases by increasing the loop number down to less then 1% for = 8.
It is interesting to see the main effect of the multiloop corrections occurs already at the 2 level, where the 1 results experience the strongest screening effect. Furthermore, as explicitly argued in Ref. [33] the inclusion of the two-loop corrections to the flow of the interaction allows to substantially enrich the virtual excitation content of the fRG equations. By looking at Fig. 7 one could deduce that, performing a post-processed evaluation of the susceptibility, as well as of the fermi-boson vertex, brings them closer to the converged values than the corresponding results coming from the fRG flow (blue curves). However, it has to be stressed that the convergence trend observed in the magnetic channel for the post-processed χ and γ 3 does not apply in general. Counterexamples can be observed, for instance, in the s-wave secondary channels (i.e., charge and superconducting), where the post-processed evaluation of the 1 susceptibility not only leads to an overscreening (i.e., an underestimation with respect to the converged result), but, e.g., in the charge channel, to even unphysical results, as can be observed in Fig. 8. Here, the s-wave susceptibility in the density channel is plotted at q = (0, 0) as a function of the bosonic Matsubara frequencies. One observes negative values of the post-processed susceptibility (green line) at finite bosonic frequencies, which are restored to positive values by the multiloop corrections (red line). An attempt to explain this different trend between the dominant (magnetic) and the secondary channels (density and superconducting) is extensively discussed in Appendix D and summarized in the following observations.
As explicitly derived in Appendix D, the Λ-derivative of the formal definition for the susceptibility reported in Eq. (8) (as well as Eq. (10) Fig. 14) contained in such terms act as a screening effect provided by the complementary channels (η ) to the one considered (η = η ). Because of the imbalance between the 1 approximation for the two-particle vertex γ 4 and Σ, and the 3 diagrams included in the modified "post-processed flow equation" for the susceptibility (see Appendix D), this screening effect ends up being overestimated. Nonetheless, it represents a minor effect on the dominant (magnetic) channel, where the imbalance effect is still governed by the large 1 antiferromagnetic contribution. It could however lead to major changes in the secondary channels, which are affected by the strong screening effect of the magnetic channel appearing on the 3 -like terms. The overscreening affects all frequencies, because of the internally summed diagrams. Therefore, it is particularly severe at nonzero frequencies where the susceptibility assumes small values. This explains the unphysical negative values of the density susceptibility in Fig. 8.
By applying different fRG cutoff schemes, we obtain further tests of the reconstruction of the full derivative of Eq. (8) provided by the multiloop approach. In Fig. 9 we compare the results shown already in Fig. 7 (central upper panel) for T = 0.2t using a frequencydependent regulator (Ω-flow) with the results for χ at the same temperature obtained by a trivial or flat regulator, also known as interaction or U -cutoff [75]. Differently from the Ω-flow, the U -flow just multiplies the bare propagator with a scale factor that is increased from 0 to 1. Hence, it does not provide any cutoff in energy during the fRG flow so that all energy scales are treated on an equal footing. The insertion of the multiloop corrections into the fRG flow equations, as already observed in a different system in Ref. [55], makes the mfRG calculation almost independent, at high enough loop-order, from the specific regulator considered. A more detailed analysis of our results revealed a persisting small discrepancy even for higher loops. Since it vanishes in absence of self-energy corrections, we attribute it to the truncation of the form factor basis in the vertex flow which prevents the reconstruction of the full derivative of the self-energy. The substantial reduction of the pseudocritical temperature (T pc ) provided by the multiloop corrections is shown in Fig. 10. Here, the inverse 1 fRG antiferromagnetic susceptibility (blue line) is plotted as a function of temperature and compared to the one computed with 8 mfRG calculation (green line). We note that, in principle, the formal equivalence between the mfRG and the parquet approximation should even guarantee the fulfillment of the Mermin-Wagner theorem [71] as this is fulfilled by the parquet approximation [64].
Hence, a frequency-momentum converged mfRG calculation should yield a complete suppression of the pseudocritical temperature down to zero. It is, however, very hard to prove this result numerically, due to the quasi-long-range nature of the spatial fluctuations responsible for the Mermin-Wagner theorem. In fact, the avoided onset of a true long-range antiferromagnetism at finite temperature T is associated with the appearance of antiferromagnetic fluctuations with an exponential growing correlation length (see, e.g, discussion in Ref. [72]). Their occurrence has been indeed explicitly verified in several many-body calculations [8,72,[76][77][78][79] compatible with the Mermin-Wagner theorem. While these low-temperature exponentially extended correlations make the overall physics of our system very similar to that of a true AF ordered phase [80], being associated with a rapid crossover towards a low-temperature insulating behavior, they also make it numerically impossible to access the T → 0 limit, because of the finiteness of any momentum grid discretization. In fact, in the temperature range where we could achieve a satisfactory momentum-convergence of our 8 results the antiferromagnetic susceptibility does not show yet any evidence of the exponential behavior expected in the low-temperature regime. This finding is, however, fully consistent with the most recent estimates of the temperature range, below which the exponential behavior of χ AF should become visible: According to the most recent DΓA and Dual Fermion studies [8,14,79,81] such a "crossover" temperature would be lower than the ordering temperature of DMFT. The latter, for U = 2 is T DM F T N ∼ 0.05 (β = 20), i.e., already twice smaller than the lowest temperature considered in the present work.
To conclude this section, we analyze the effect of the fRG multiloop corrections on some d-wave physical susceptibilities which, although suppressed in the particle-hole symmetric case, play an important role in describing the phase diagram of the 2D Hubbard model, most notably away from half filling [27,35,57,82,83]. In particular we analyze the static (ω = 0) d-wave susceptibility in the superconducting channel for q = (0, 0) (dSC), as well as the static d-wave susceptibility in the charge channel for a bosonic momentum transfer q = (0, 0) (dPom), which would become dominant in the case of the so-called "Pomeranchuk" instability. The staggered d-wave charge density wave (dCDW) susceptibility for q = (π, π) has not been shown because of its degeneracy with the correspondent d-wave superconducting one. In fact, one can formally demonstrate that in a SU(2) and particle-hole symmetric case, where the system becomes invariant under pseudospin rotation, the pairing susceptibility at q = (0, 0) associated to a specific symmetry of the order parameter is degenerate with the staggered (q = (π, π)) CDW associated to that specific symmetry. In Fig. 11 we display the result of a fRG calculation where, in addition to the s-wave form factor, the form factors indicated as 1 and 2 in Table 1 have been used. As in Fig. 7, the blue line indicates the fRG result obtained by the integration of Eq. (23a) up to a specific -loop order, alongside the corresponding ( -loop) mfRG equations for Σ and γ 4 . The green line represents the post-processed result for the d-wave susceptibilities calculated from a s+d-wave -loop order mfRG results for the self-energy (Σ) and the two-particle vertex (γ 4 ). The red line has been obtained, similarly to the green one, from a s-wave -order mfRG results for Σ and γ 4 . One notices that, differently to the antiferromagnetic case, the relative difference between blue and green lines with respect to the convergence value is, at the 1 -level, of the order of few percents and lowers even down ground at T = 0.5t, i.e., they would not flow strongly in their 'native' d-wave channels. Going to lower T and in particular out of half-filling, this will likely change, as the particle-particle diagrams will enhance any attractive pairing component. Therefore, it is a priori not clear if the d-wave susceptibilities computed at lower T by projecting the vertex made up from s-wave bilinears could provide satisfactory physical results. Nevertheless, we argue that they serve as useful theoretical test objects for the convergence in the order of the multiloop corrections. This is because the effective d-wave interactions captured this way can be understood as two-particle irreducible (2PI) interactions in the d-wave pairing or charge channels, generated purely by s-wave one-loop processes. These 2PI d-wave quantities are non-singular but zero at lowest order in U in typical cases. Hence they can be expected to be dominated by diagrams of finite order in U that should exhibit stronger multiloop effects. In contrast with these terms, the missing boosts in the respective native channels, e.g., in the pairing channel, would just be a higher-order ladder summation of, for T → 0, increasingly singular one-loop diagrams. Hence, if multiloop convergence is reached in the two-particle irreducible interactions, it is likely that the same degree of convergence would be found in the true susceptibilities. This idea leads us to consider the data shown in Fig. 12.
As already visible for T = 0.5t in Fig. 11, the post-processing calculations exhibit a weak dependence on the loop number (with a relative fluctuation less then 1 o / oo ). This is confirmed in Fig. 12 where the post-processed inverse d-wave susceptibilities in the aforementioned channels are calculated out of an s-wave 1 (blue and green lines) and 8 (red and yellow dashed lines) fRG flow. As it is apparent in the figure, the effects of the multiloop corrections are insignificant compared to the variation of the inverse susceptibilities in temperature.

Conclusions
We have presented a comprehensive study of forefront algorithmic implementations of the fRG for interacting fermions on 2D lattices. While we focused on the 2D Hubbard model, the methodological improvements discussed here can provide a useful guidance for the generalization to other systems. Our main goal is to illustrate the progress achieved when going beyond the approximations routinely made in most previous fRG computations. In particular, we have worked on the following aspects: (i) an accurate and converged treatment of both the momentum and frequency dependence of the vertex function together with its asymptotic structures; (ii) the inclusion of the self-energy and its feedback in the fRG flow; (iii) the implementation of the multiloop corrections beyond the standard 1 .
Regarding the first aspect (i), we have kept the more general dependence of the twoparticle vertex on all three Matsubara frequencies. We extend previous works [40,43,[83][84][85][86] by exploiting an "economic" description [33] provided by an efficient parametrization of the high frequency asymptotics [31]. We could show that this parametrization can be brought to convergence in the number of frequencies employed, i.e. it the results do not change if more frequencies are used. We combined this treatment of the frequency dependence with the truncated-unity technique for the momentum dependence, whose form-factor expansion was also shown to converge quickly for our test case [30].
With a frequency-dependent flowing interaction, we could also compute a momentum-and frequency-dependent self-energy, which has been fed back into the flow of the two-particle vertex. Through a systematic analysis of specific observables -in particular of the response functions -we could assess the effects of the improved algorithmic implementation with respect to previous results and demonstrate how, for the parameters studied, the fRG results can be converged in the number of considered frequencies. An analogous convergence could be also established for the 2D momentum dependence.
The major advancement achieved in this work is, however, the implementation of the multiloop corrections both for the flow of the two-particle vertex as well as for the flow of the coupling to external fields and the corresponding susceptibilities. The multiloop extension, so far only tested for a (prototypical) toy model [55], adds more virtual excitations to the flow of the two-particle vertex compared to the previously used 1 truncation. As it was diagrammatically shown [54,55], if truncated fRG results are converged with respect to the loop order, they exactly reproduce the parquet approximation (PA), not only concerning the topology of the summed diagrams, but also -quantitatively -their precise weight. This has been also recently confirmed by a formal analytical derivation of the multiloop fRG equations [56]. From this property, it follows that the results of a loop-converged fRG algorithm become completely independent from the employed cutoff scheme, at least if all modes are integrated out at sufficiently high temperature.
Previously, it was not clear how the contributions missing in the 1 truncation would influence the results quantitatively. On the numerical level, the effort for including the multiloop corrections to order only rises linearly in , i.e. the situation is far better than if one really had to compute all higher-loop diagrams. Our studies show that the multiloop corrections can be included also in 2D up to rather high orders of = 8. We find that the observables converge quite nicely when the multiloop order is increased. While it is not obvious that this quick convergence will hold for all model parameters and for all models of interest, our study shows that these checks can be performed with feasible numerical efforts. This adds an new important degree of quantitative control to the fRG, at least in the weak to intermediate coupling regime where the PA can be considered accurate. At stronger coupling, where low-frequency vertex corrections beyond the PA might appear [31,49,[87][88][89], the mfRG could provide a much better [14] setup for the proposed combination with the DMFT [22,26]. The loop convergence of our fRG results is also reflected in the progressive reduction of the dependence of our fRG results on the chosen cutoff scheme, which appears completely suppressed at the 8 level.
The incorporation of the multiloop contributions has also another rather appealing and quantitatively important aspect, giving rise to an additional very useful type of convergence. It has been known that response functions can be computed in two different ways in RG approaches and that the results differ due to the involved approximations. One way is to consider the flow of couplings of 'composite operator' bilinears in the primary degrees of freedom to external fields of appropriate type. Then the response function is obtained as renormalization of the propagator of the external field.
The other way, referred to as post-processing, is to compute the response functions by means of their diagrammatic expression, evaluated from the dressed bare fermion bubbles and the two-particle vertex at the end of the flow. In fact, in some cases arguments were made (see, e.g., Ref. 90 and references therein) that the external field methods should give more controlled results, i.e., that composite operators should be renormalized separately, because, at the level of the approximations made, the post-processed quantities, which involve the integration over all energies and momenta, are more afflicted by approximation errors. In our study, the multiloop extension of the response function flow allows us to show that also the flow of the response functions becomes an exact scale derivative of the post-processed response function. This establishes the formal equivalence of the two ways to compute response functions on the multiloop level. This formal equivalence is remarkably reflected by our numerical results, which exhibit a clear convergence of the two approaches: If the multiloop convergence is achieved, and frequency and momentum dependencies as well as the self-energy feedback are included appropriately, the fRG results for the response functions are unambiguous. The corresponding data can be used for quantitative studies and directly compared with other numerical techniques or with experiments, if the effective modelling of the problem is sufficiently realistic.
In summary, our study shows how the fRG algorithms for two-dimensional fermionic lattice models can be brought to a quantitatively reliable level at weak to moderate couplings, as long as the parquet approximation is appropriate. This goal has been reached by means of an economic, but accurate, treatment of the momentum and frequency dependencies which takes into account the asymptotic structure of the two-particle vertex and the self-energy during the fRG flow. This fRG framework has been supplemented with the implementation of the multiloop corrections to the 1 truncation scheme.
The algorithmic implementation presented in this paper can be applied to other model Hamiltonians, and -within the framework of the 2D Hubbard model -to broader parameter regimes in future works. For instance, if the Fermi surface displays a given curvature, due to the inclusion of, e.g., more hopping terms or changes of the band filling, the dominance of the AF channel will be weakened and the pseudo-critical scales will become smaller. For such cases the convergences of the different approximation might possibly vary. In particular, since the generation of d-wave pairing tendencies in third order of the bare coupling involves 2 diagrams that are only partially captured in the 1 truncation, we would expect the impact of the multiloop corrections to become more noticeable.

A Symmetries and symmetrized notation
Here we illustrate how diagrammatic and lattice related symmetries can be expressed in an easy way and how they are implemented in our code. Directly related to the symmetries is the question if one uses the symmetrized or the non-symmetrized notation for the momentum and frequency dependence of the channels. In Section 3.1 we argued that the non-symmetrized notation leads to more readable flow equations, bubbles and projection matrices. Therefore we adopted primarily this notation. The symmetries, however, are much easier to express in the symmetrized notation. While in the non-symmetrized notation, simple relations like the crossing relation involve multiple form factor combinations, in the symmetrized notation we find a one-to-one correspondence. Therefore we here use for both momentum and frequency the symmetrized notation (s), which is related to the non-symmetrized (ns) by (52c)

A.1 Lattice related symmetries
First we specify how lattice related symmetries are reflected in the form factor expansion of the channels in the symmetrized notation. The lattice symmetries always depend on the system and we here focus on the 2D Hubbard model on a square lattice, where we have for example the rotation of π/2 around the z-axis and the mirroring at the y-axis as independent symmetry operations. Under any of these operations, or combinations of them, applied simultaneously to all momentum dependencies, the expressions of the channels are invariant. This can be translated into the form factor expansion bŷ where F is any of the channels D, C or P . The frequency dependence is not affected and is therefore omitted. We here exploited the symmetry under consideration and introduced a variable change. If the form factors are chosen in such a way that under this symmetry operation any form factor is related to a linear combination of others, described by the matrix MR −1 (k), it holds in addition If moreover, the symmetry operation on every form factor yields a single other form factor expressed by the vector VR −1 , the above relation simplifies tô where the only difference is a possible sign change taken into account by S V R (n) . These assumptions hold for the form factors used in the present implementation (see Table 1), but are not necessarily valid for an arbitrary choice of form factors.

A.2 Diagrammatic symmetries
In addition to the lattice related symmetries, there are diagrammatic symmetries which are independent of the geometry of the system. Considering a two-particle fermionic vertex, we can apply the crossing symmetry simultaneously to the annihilation and the creation operators, recovering the following relations: time reversal and complex conjugation for which we refer to Ref. [91]. In the SU(2) symmetric case, by projecting the vertex φ to the form factor basis and adopting the symmetrized notation, one has that Eq. (56) gives where Π m is the parity associated to the momentum inversion of the form factor m defined as The time reversal symmetry reads and the complex conjugation A.3 Connection between K 2 andK 2 In Section 3.1.2 we argued thatK 2 can be obtained from K 2 by symmetry. For the pp and ph channel the time reversal symmetry exchanges the two fermionic dependencies while keeping the transfer frequency and momentum fixed. The same holds for the ph-channel by using the combination of the crossing and the time reversal symmetry. Taking the limit of large frequencies for the first and second fermionic frequency respectively, we obtain trivially B Formal derivation of the fRG flow equations for χ and γ 3 In this section we provide an explicit derivation of the flow equations for the response functions. As anticipated in Sec. 2.2, we start by coupling the fermionic bilinears to an external source field J, by adding the following scalar product where n indicates the momentum structure of the fermionic bilinears coupled to the field J n η . Since the density is in general not charge conserving, it is convenient to use the Nambu formalism that allows for a more concise derivation of the flow equations of the physical response functions. We rewrite Eqs. (1) and (2) in the Nambu basis [92,93] where s = ± represents the Nambu index and The matrices α η (with η = {d, m, sc}), which define the Nambu index structure in the different physical channels, are given by In order to derive the flow equations for the fermion-boson vertex of Eq. (15) and the susceptibility of Eq. (14) we start from the so-called Wetterich equation [94] ∂ where Γ Λ represents the scale-dependent effective action, which is a function of the functional variable J η and the Nambu field φ, Q Λ 0 is the inverse non-interacting Green's function and the dot denotes the derivative with respect to the flow parameter Λ. Further, the matrix and were we used, where ∂ and∂ applied to the effective action Γ Λ are a shorthand notation for the functional derivative of with respect to φ andφ, respectively. Following the derivation of Ref. [2], we introduce the matrix Thus, we can recast (Γ (2)Λ [J η , φ]) −1 = (1 − G Λ U Λ ) −1 G Λ and expand the inverse matrix in a geometric series We can now insert Eq. (70) in Eq. (67). Expanding up to second order yields represents the matrix diagonal form of the single scale propagator, and we exploited the cyclic property of the trace. After applying the trace to the matrices in the curly brackets, we can expand the effective action in powers of the fermionic Nambu fields and the external bosonic source field Note that the index x = {s, k} combines the Nambu index s and the fermionic quadrivector k = (ν, k) (here we disregard additional quantum numbers, as e.g., orbital), while y = {n, q} combines the momentum structure of the coupling to the bilinears, n, with the bosonic quadrivector q = (ω, q). Inserting this expansion in Eq. (71), we compare the expansion coefficient related to the same order on the fields on both sides of the equation. For n 1 = 0 we recover the standard fermionic hierarchy of flow equations [1,2]. For n 1 > 0 we can derive the flow equations for the fermion-boson vertex (n 1 = 1, m 1 = 1) as well as for the boson-boson vertices or susceptibilities (n 1 = 2, m 1 = 0). In Nambu notation, the flow equation for the susceptibility reads and the one for the fermion-boson vertex is In the second term on the r.h.s. of Eq. (73), γ Λ 2m 1 +n 1 =γ Λ 2+2 represents the functional derivative of the effective action with respect to two bosonic and two fermionic Nambu fields. Figure 13: Simplified diagrammatic representation of the flow equations for the susceptibily (first line) and the fermi-boson vertex (second line) illustrating the topological structure of the diagrams. The circle, triangle and the square represent the susceptibility χ, the fermi-boson vertex γ 3 , and the two-particle vertex γ 4 , respectively.
The two Eqs. (73) and (74) are schematically shown in Fig. 13. If one neglects the second term in both r.h.s., they correspond to the 1 fRG equations for the response functions. Both γ 3 and χ do not feed back into the flow equations for γ 4 and Σ.
C Connection between the vertex asymptotics and the response functions In this appendix we demonstrate that the integration of the fRG flow equations for the so-called kernel functions K 1 and K 2 mentioned in Section 3, coincide with the s-wave susceptibility and fermion-boson vertex resulting from the flow. Let us write explicitly the flow equation for the asymptotics K Λ 1,η andK Λ 2,η , with η = {sc,d,m}, obtained from Eq. (43) in the limit of infinite fermionic Matsubara frequencies ν and ν and whereφ Λ η is given byφ the bare vertex γ 0 4,η = ∓U corresponds to the Hubbard interaction (with the minus sign for η = sc, d and the plus sign for η = m), and the asymptotic vertex functionK Λ 2,η is related to K Λ 2,η by symmetry (see Appendix A). For local bare interactions, the only non-zero elements of the matricesK Λ 1,η and γ 0 4,η correspond to both form factors being equal to zero, and of K Λ 2,η (K Λ 2,η ) to a vanishing second (first) form factor. The connection between the vertex asymptotics and the response function is shown by induction using the assumption For the initial condition, it holds γ Λ init one has α = γ 0 4,η = ∓U δ n,0 δ n ,0 . Considering (γ 0 4,η + K Λ 1,η +K Λ 2,η ) for an arbitrary value of Λ, we can identify the flow equation of the asymptotics with the one of γ 3 , see Eq. (23a). Therefore Eq. (78) applies also for the following Λ step. As a consequence we can extract the fermion-boson vertex from the vertex asymptotics. Finally, inserting Eq. (78) into (75), we obtain the flow equation for the susceptibility (23a).
The s-wave fRG results for the susceptibility and the fermion-boson vertex can be extracted from the asymptotic vertex functions K Λ 1,η andK In this section we explicitly provide the scale derivative of Eqs. (8) and (10) for the case in which the Σ and γ 4 entering the r.h.s. are obtained from the integration of the corresponding 1 flow equations. We first consider Eq. (10) and, after introducing a Λ-dependence of the Green's functions and of γ 4 on the r.h.s., perform the full derivative with respect to Λ. For simplicity we here consider the magnetic vertex as example, which is directly related to the particle-hole crossed vertex by γ 4,m (q, k, k ) = γ 4,ph,↑↑ (q, k, k )−γ 4,ph,↑↓ (q, k, k ) = γ 4,ph,↑↓ (k −k, k, k+q) = −γ 4,ph,↑↓ (q, k, k ) , where the boxes indicate the conventional 1 approximation. The internal loops in red provide the particle-hole and particle-particle contributions respectively. The empty dot represents the bare AF fermion-boson vertex (γ 0 3,m ) n,m = δ n,m .
where we used the SU(2) and crossing symmetries [31]. The derivative of the fermion-boson vertex with respect to Λ, as obtained from Eq. (10), reads post-proc S,m +Π magnetic channel. Following the derivation of Eq. (80) one obtains with the diagrammatic representation is provided in Fig. 14 (second line). One observes the appearance of additional terms on the r.h.s. of the post-processing flow equations for γ 3 and χ with respect to their standard 1 equations, indicated by the boxes in Fig. 14. Besides the terms containing the Λ derivative of the self-energy (which are included in the Katanin corrections [44]), let us draw the attention to the last two diagrams appearing on the r.h.s. for both ∂ Λ (γ Λ 3 ) post-proc and ∂ Λ (χ Λ ) post-proc . The diagrammatic structure in terms of loops is of second order for γ 3 and of third for χ. The integration of these post-processed flow equations, along with the 1 flow equations for Σ and γ 4 , would generate the last two diagrams already at the first integration step Λ init + dΛ (with dΛ < 0 in the Ω-flow), providing the following The first term vanishes due to the Pauli principle (φ Λ init ph = 0, see Ref. [33]), and the last one provides a negative contribution which reduces the 1 term. In fact, the unscreened particleparticle bubble enteringĈ[φ Λ init pp ] n,m has the same sign of the unscreened (magnetic) S − G bubble. This overall suppression by the additional 3 -like terms is a general feature of the post-processed fRG scheme. The unbalance between the 1 γ 4 flow, which topologically cuts part of the parquet diagrams, and the additional 3 -like diagrams of the susceptibility flow, leads to an artificial overscreening of the conventional 1 calculation. Analogous conclusions can be drawn for the density and superconducting channels. Thus one expects a pronounced effect in the secondary channels because the dominant channel enters the internal loop of one of the two 3 -like additional diagrams, resulting in a reduction with respect to the converged data. In contrast, the dominant channel will not be affected that strongly, presenting only a slight overestimation of the post-processed susceptibility at the 1 level (see Fig. 7). Moreover, since this overscreening affects all frequencies, it may be responsible for the unphysical negative value of the density susceptibility observed at finite frequencies in Fig. 8. In particular, since the parquet diagrams disregarded in the 1 approximation depend on the cutoff, the detected unphysical results in the secondary channels were observed to be more severe for the interaction flow. We finally note that this opposite effect of the density and the superconducting channels with respect to the dominant magnetic channel has been observed also in Ref. [95] by analyzing the effect of the parquet decomposition of the vertex on the self-energy.  We here provide the derivation of the 2 corrections to the conventional 1 truncated flow equations. The derivation follows the scheme adopted for the flow equation of the twoparticle vertex as reported in Ref. [16]. Our goal is to include the feedback of γ Λ 5 onto the flow equation for γ Λ 3 , see Eq. (74), at the second order in the effective interaction. From the derivation provided in Appendix B, one sees that the differential equation for γ Λ 5 is given by the sum of all diagrams which have the topological structure depicted in Fig. 15. The first and the second diagrams on the r.h.s. are at least of third order in the effective interaction since γ Λ 7 (depicted by a heptagon) and γ 5 (depicted by a pentagon) are at least O((γ Λ 4 ) 3 ) and O((γ Λ 4 ) 2 ), respectively. Therefore, we can restrict ourselves to diagrams with a topological structure of the third one. Its contribution can be obtained by taking the following functional derivative evaluated at zero fields ∂ Λ γ Λ 5 (y, x 1 , x 2 , x 1 , x 2 ) = ∂ 5 ∂J η (y)∂φ(x 1 )∂φ(x 2 )∂φ(x 2 )∂φ(x 1 ) where x = {s, k}, y = {η, q} and ∂ Λ,S acts only on G Λ and returns the single-scale propagator S Λ . At this point we integrate the r.h.s. which is an easy operation once we take into account that i) one can replace S Λ = ∂ Λ,S G Λ by the full derivative ∂ Λ G Λ since their difference due the derivative of the self-energy is of higher order in the effective interaction γ Λ 4 , and ii) one can let the scale derivative act also on γ Λ 4 since its derivative is at least of order O((γ Λ 4 ) 2 ). According to these arguments, the r.h.s. of Eq. (84) can be approximated by the total derivative with respect to the Λ and integrated to γ Λ 5 (y, x 1 , x 2 , x 1 , x 2 ) = ∂ 5 ∂J η (y)∂φ(x 1 )∂φ(x 2 )∂φ(x 2 )∂φ(x 1 ) 1 3 tr The only terms surviving the functional derivative are all connected diagrams composed by two two-particle vertices γ Λ 4 and one fermion-boson vertex γ Λ 3 . What distinguishes the first and the second contributions of Eq. (85) is the position of γ Λ 3 which can be inserted at all∂∂ in the first line, while is restricted to a single∂∂ in the second one because of the conservation of Nambu particles. Moreover, the first term accounts for two-particle vertices whose external lines are always a particle and a hole, whereas in the second term they are attached to two particles ∂∂Γ Λ and two holes∂∂, respectively. The topological structure of these two contributions is schematically shown in Fig. 16. The last step consists in closing these diagrams in all possible ways by means of the single-scale propagator and adding them to the flow equation of γ Λ 3 . Hence, one obtains 2 approximated flow equations for γ Λ 3 which contain terms of the order O((γ Λ 4 ) 2 ) in the effective interaction. We can classify [16,39] the 2 corrections according to theit topological structure, with overlapping loops (Fig. 17 (b)) and non-overlapping loops ( Fig. 17 (a)). We observe that the latter can be included in the 1 equations by using the Katanin correction [44] where S Λ → S Λ + G ΛΣΛ G Λ . The remaining 2 corrections have as building block the 1 diagrams of the flow equation of γ Λ 4 . Translating our Nambu formalism to the physical fields, those corrections yield Eq. (24).

F Implementation details
Here we provide the explicit form of γ 4,{P,D,C} appearing on the r.h.s. of Eq. (36). By using the parquet decomposition in the diagrammatic channels (see Eq. (32)), the first contribution of the projections of the four-point vertex onto the different channels, is the projection of the fully two-particle irreducible vertex, approximated by its first order in the on-site Hubbard Secondly, every channel, written in its natural bosonic-fermionic notation on the l.h.s. of Eq. (36), need to be projected onto the complementary channels. The projection of one channel φ r to another leads to a linear combination of its frequency arguments (see Eq. (22) for the physical channels and Eq. (90a) to Eq. (90f) for the diagrammatic channels). In momentum space, the projection is more involved due to the form factor dependence. Following the procedure of Ref. [30], we identify the projection matrices which describe the momentum translation from one channel to another using a matrix multiplication where . . . stands for the channel specific translation of the frequency dependencies. We exemplify the projection for the channel D to P . In momentum space, it reads We now transform the form factors to real space and shift the momentum dependence in order to get the matrix form of (87) [P [φ ph ](iω l , iν o , iν o , q)] n,n = m,m π −π dK RR 1 R 2 The same procedure for every channel projection leads to the matrix equations with the following projection matrices for the non-symmetrized notation A P,D n,n ,m,m (l, q) = A P,C n,n ,m,m (l, q) = A D,P n,n ,m,m (l, q) = A D,C n,n ,m,m (l, q) = A C,P n,n ,m,m (l, q) = A C,D n,n ,m,m (l, q) =