{'fraction_non_alphanumeric': 0.04752063730300756, 'fraction_numerical': 0.014200773538070774, 'mean_word_length': 4.293082375947201, 'pattern_counts': {'":': 0, '<': 2, '<?xml version=': 0, '>': 3, 'https://': 0, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 4, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
Fulfillment of sum rules and Ward identities in the multiloop functional renormalization group solution of the Anderson impurity model

We investigate several fundamental characteristics of the multiloop functional renormalization group (mfRG) flow by hands of its application to a prototypical many-electron system: the Anderson impurity model (AIM). We first analyze the convergence of the algorithm in the different parameter regions of the AIM. As no additional approximation is made, the multiloop series for the local self-energy and response functions converge perfectly to the corresponding results of the parquet approximation (PA) in the weak- to intermediate-coupling regime. Small oscillations of the mfRG solution as a function of the loop order gradually increase with the interaction, hindering a full convergence to the PA in the strong-coupling regime, where perturbative resummation schemes are no longer reliable. By exploiting the converged results, we inspect the fulfillment of (i) sum rules associated to the Pauli principle and (ii) Ward identities related to conservation laws. For the Pauli principle, we observe a systematic improvement by increasing the loop order and including the multiloop corrections to the self-energy. This is consistent with the preservation of crossing symmetries and two-particle self-consistency in the PA. For the Ward identities, we numerically confirm a visible improvement by means of the Katanin substitution. At weak coupling, violations of the Ward identity are further reduced by increasing the loop order in mfRG. In this regime, we also determine the precise scaling of the deviations of the Ward identity as a function of the electronic interaction. For larger interaction values, the overall behavior becomes more complex, and the benefits of the higher-loop terms are mostly present in the contributions at large frequencies.


I. INTRODUCTION
The many-electron problem poses a formidable challenge to modern solid-state physics, involving a large number of degrees of freedom at different energy scales. In general, although the exact solution cannot be computed, some of its fundamental properties are known a priori. Specifically, the exact solution is guaranteed to obey the Pauli principle, which manifests itself in sum rules and the crossing symmetry of fourpoint correlators. At the same time, it also fulfills Ward identities (WIs) related to thermodynamic and quantum-mechanical principles. This knowledge usually provides an important "compass" for constructing suitable approximation schemes. For a given approximation, however, the preservation of all fundamental features of the exact solution cannot be guaranteed [1]. For instance, approximate schemes constructed from the Luttinger-Ward functional, so-called conserving approximations, maintain WIs-an important aspect when comparing with spectroscopic experiments-but violate the crossing symmetries. On the other hand, it is known that approximate approaches specifically designed to guarantee the crossing symmetries, such as the parquet approximation (PA) [1][2][3][4][5][6], violate the WIs to a certain degree [7][8][9]. Hence, investigating how this trade-off actually manifests itself in advanced quantum many-body methods will provide significant theoretical insight.
In this work, we analyze these issues within the functional renormalization group (fRG) for interacting Fermi systems [10][11][12][13][14], which can be used as a framework for introducing powerful new approximation schemes. Specifically, we consider the recent multiloop extension (mfRG) of fermionic fRG in the vertex expansion [9,15,16] and apply it to the Anderson impurity model (AIM), a paradigmatic model of many-body physics. Reasons for focusing on this particular model are given below.
Computation schemes based on the fermionic fRG can be designed to treat the characteristic scale-dependent behavior of correlated electrons in a flexible and unbiased way. The most commonly used implementations employ the oneloop (1 ) truncation of the exact hierarchy of flow equations. There, one neglects three-particle and higher vertices, which can be justified, e.g., from a perturbative perspective. Several studies have discussed the nonconserving nature of 1 fRGbased schemes, and possible routes for mitigating the violation of the associated WIs have been proposed [10,[17][18][19][20]. An important example is the widely used Katanin substitution [17]. In its most common form, it incorporates some contributions from the three-particle vertex as two-loop contributions to the flow of the two-particle vertex via (one-particle) selfenergy corrections. The Katanin substitution was designed to better fulfill WIs-we here present first numerical results to quantitatively assess this aspect. Conversely, WIs have also been used to propose new truncation schemes [21][22][23].
The multiloop extension of the fRG approach, mfRG, includes all contributions of the three-particle vertex to the flow of the two-particle vertex and self-energy that can be computed with numerical costs proportional to the 1 flow. In doing so, it sums up all parquet diagrams, formally reconstructing the PA if loop convergence is achieved [9,15,16]. This ensures self-consistency at the one-and two-particle level, in that the PA is a solution of the self-consistent parquet equations [1]. It also ensures the validity of one-particle conservation laws, but not of two-particle ones [9].
Whether or not mfRG yields quantitative improvements arXiv:2110.07455v2 [cond-mat.str-el] 21 Apr 2022 over the 1 truncation depends on the context. For a zerodimensional model with a logarithmically divergent perturbation theory, it was recently shown [24] that the leading logarithms can be obtained in the 1 truncation, in which case the higher-loop contributions incorporated via mfRG thus are subleading. Similarly, 1 fRG treatments of the interacting resonant level model [25][26][27][28] as well as inhomogeneous Tomonaga-Luttinger liquids [29][30][31][32][33][34][35] should yield a proper summation of the leading logs governing the infrared behavior of these systems. By contrast, a quantitatively precise description of the weakly interacting two-dimensional Hubbard model could only be achieved with a full multiloop computation [36][37][38]. It is thus of interest to analyze the multiloop series for a model whose perturbation series lacks a leadinglog classification, but which is less complex than the Hubbard model. This criterion is satisfied by the AIM. We study it here at finite temperature in the imaginary-frequency Matsubara formalism. A Matsubara treatment of the AIM suits our purpose for two further reasons. First, a numerically exact solution is available as a benchmark via Quantum Monte Carlo (QMC) methods [39]. Second, recent algorithmic and methodological advances [5,6] make it possible to track the full frequency dependence of the two-particle vertex functions of the AIM [40][41][42], including their non-trivial asymptotic structure [6]. The numerical (m)fRG equations can be then solved to great accuracy and without any further approximations. This sets our study apart from recent mfRG applications [36,37] to more complex systems (where additional approximations for the momentum dependence [43,44] of the problem were necessary) and builds upon previous frequency-dependent fRG studies of the AIM [45][46][47][48][49][50], paving the way for a systematic inspection of sum rules and WIs in mfRG and parquet approaches.
From a more general perspective, we note that flows of the truncated fermionic fRG or mfRG can a priori be expected to be reliable for weak to intermediate interaction strengths only. However, nonperturbative [51] parameter regimes of, e.g., the Hubbard model can be accessed [52,53] via fRG by proceeding as follows: first evoke dynamical mean-field theory (DMFT) [54] to solve a self-consistent AIM (by non-fRG means, e.g. QMC or the numerical renormalization group [55]); then use fRG to systematically include nonlocal correlations missed by DMFT [42]. This procedure defines the so-called DMF 2 RG scheme [56,57]. So far, it has been implemented in the 1 truncation, but multiloop extensions are conceivable, too. Our careful investigation of the mfRG solution of the AIM may also provide valuable methodological information for future multiloop DMF 2 RG developments. For example, analogous vertex frequency parametrizations are needed for a mfRG treatment of the AIM and for the mfRG part of DMF 2 RG computations. Furthermore, the study of the mfRG convergence properties as well as of crossing symmetries and WIs for different coupling strengths will represent an important guidance for DMF 2 RG calculations relying on multiloop resummations.
The structure of our paper reflects the main scientific questions raised above. After introducing the required formalism in Sec. II, we present a detailed analysis of the mfRG solution of the AIM in Sec. III. We illustrate how the convergence to the corresponding results of the PA is perfectly achieved in the weak-to intermediate-coupling regime and also discuss the appearance of increasing multiloop oscillations in the strongcoupling regime. Having defined the parameter region of convergence for mfRG applied to the AIM, we analyze in Sec. IV the fulfillment of the sum rules associated to the Pauli principle as well as of the WIs related to conservation laws. We discuss the systematic effects observed as a function of loop order, and separately consider the low-and high-frequency parts of the WIs. Throughout, we also include results obtained via the Katanin substitution, allowing its merits to be compared to those of 1 or higher-loop schemes. We summarize our conclusions and perspectives for future developments in Sec. V, and discuss additional technical aspects relevant to a more specialized readership in the Appendices.

II. FORMALISM
In this section, we concisely introduce the methods and concepts underlying the calculations presented in the following sections. For brevity, we reduce the formal derivations to a minimum, referring to prior works for more explicit discussions. In Sec. IV, we extend the formalism where needed for the analysis of the Pauli principle and WIs.
A. Anderson impurity model Throughout this paper, we consider the AIM [58] close to the wide-band limit [59]. The Hamiltonian is given bŷ whered † σ (d σ ) is the creation (annihilation) operator of electrons localized on the impurity site, andĉ † k,σ ,ĉ k,σ are the corresponding operators for the bath electrons. The energy on the impurity site is denoted by d and the dispersion relation in the bath by k . The first term in the second line represents the local interaction, where U is the interaction strength and n σ =d † σdσ with σ ∈ {↑, ↓}. The second term accounts for the hopping onto/off the impurity site, where we consider a k-independent hybridization strength V k = V . We set V = 2, thereby measuring energy in units of V /2. We use a boxshaped DOS for the bath electrons, ρ( ) = 1/(2D)Θ(D−| |), with half-bandwidth D = 10. Further, we consider half filling, where d = −U/2 is exactly canceled by the Hartree self-energy. Thus, the (bare) propagator is purely imaginary: For |ν| D, we find ∆ ν ≈ −isgn(ν)∆ 0 with the characteristic hybridization strength ∆ 0 = πV 2 /(2D) = π/5 0.63.
For prior works using this specific AIM, see Ref. [60] and especially Ref. [51], where the physical regimes relevant for this paper are also discussed. For a more general introduction of the physics of the AIM, we refer to Refs. [59,61]. The values U = 1, 1.5, 2, 3 and 4 studied below correspond to U/∆ 0 1.59, 2.39, 3.18, 4.77 and 6.37, respectively. Throughout we fix the inverse temperature to β = 10.
B. Numerical approaches fRG flows and PA-We briefly discuss here the structure of the flow equations, both on the one-and multiloop level, as well as the PA, for the one-particle self-energy Σ, the twoparticle vertex F , and the susceptibilities χ of the AIM.
The fRG flow describes the evolution of Σ, F , χ upon tuning the scale or flow parameter Λ from an initial value Λ i to a final value Λ f . The flow parameter Λ is introduced in the quadratic part of the action, i.e., the bare (one-particle) propagator [Eq. (2)]. We consider two cutoff functions: the frequency flow (or Ω-flow for short), and the interaction flow (or U -flow) [62], With G Λi 0 = 0, the initial values of Σ and F are Σ Λi = 0, where the Hartree term is absorbed in G 0 , and F Λi = F 0 , the bare vertex of magnitude U (in our convention F σσ 0 = −U δσ σ , where↑ =↓ and vice versa). The fRG flow of Σ is determined by the two-particle vertex F contracted with the single-scale propagator S = −G(∂ Λ G −1 0 )G, which is related to the differentiated propagatorĠ ≡ ∂ Λ G byĠ = S + GΣG. For simplicity, we omit the superscript Λ here and henceforth. The flow equation for F further involves the three-particle vertex Γ (6) . If Γ (6) was known at all scales, the flow of Σ and F would be exact. This would imply, in particular, that every specific Λ dependence or cutoff choice, as in Eqs. (3) or (4), yields the same result at the end of the flow. In practice, however, Γ (6) can hardly be treated numerically and its effect on the flow of Σ and F can only be accounted for approximately. As a consequence, the results of such truncated fRG flows will generically depend on the choice of the cutoff.
The most widely used fRG implementations neglect Γ (6) , yielding approximate 1 flow equations for Σ and F . The contributions of Γ (6) that amount to self-energy derivatives can be added to the vertex flow by substituting S →Ġ. This "Katanin substitution" [17], labeled by 1 K in the following, was argued to yield a better fulfillment of WIs. A further refinement, which effectively incorporates the three-particle vertex to third order in the renormalized interaction, is obtained by the two-loop (2 ) vertex corrections [17,63,64].
Subsequently, the mfRG extension [9,15,16] was introduced to incorporate all those contributions of Γ (6) to the flow of Σ and F ensuring that their right-hand sides are total scale derivatives-which is not the case for the 1 flow-thus guaranteeing by construction that the final results are independent of the choice of cutoff. In fact, the corresponding higherloop terms of the mfRG represent the minimal additions to the conventional 1 flow required to obtain cutoff-independent results. They also provide the maximal amount of diagrammatic contributions that can be added in a numerically feasible manner. Indeed, due to the iterative structure based on successive 1 computations, these higher-loop contributions can be computed very efficiently [36]. Besides the two-dimensional Hubbard model [37], recent applications also include spin models [65,66] within the pseudofermion representation.
The mfRG was shown to formally reproduce the diagrammatic resummation of the PA. We use the vertex decomposition where γ r are the two-particle reducible vertices in channel r ∈ {a, p, t} and R 2PI the fully two-particle irreducible (2PI) vertex (notation as in Ref. [9] [67]). The PA then corresponds to the approximation R 2PI = F 0 . An analogous approximation is performed in truncated fRG flows: Neglecting Γ (6) , the vertex flow equation is of the formḞ = rγ r . Thus, only the reducible parts are renormalized, while the fully irreducible part does not flow and remains at its initial value R 2PI = F 0 . In the following, we recall the flow equations for the selfenergy, two-particle vertex, and susceptibilities, as well as the parquet equations. We will use the compact symbolic notation introduced in Ref. [9]; the explicit dependence on spin and frequencies will be given where needed.
One-loop flow-The 'standard' fRG self-energy flow [10] isΣ as illustrated in Fig. 1(a). Figure 1(b) shows an exemplary depiction of the 1 flow of the vertices, given bẏ andḞ = rγ r .Π r,S corresponds to the differentiated twoparticle propagator in channel r, with S used instead ofĠ. The flow equation for the susceptibilities χ r can be derived from the corresponding reducible vertex in the limit of large fermionic frequencies, i.e., the so-called K 1 contribution [6], where λ r are the renormalized three-point vertices (for further details and the 1 flow equation of the latter, see Ref. [10]). The 1 K flow with Katanin substitution is obtained by replacing S →Ġ, i.e.,Π r,S →Π r , in Eqs. (7) and (8). Since it includes self-energy (and not vertex) corrections from Γ (6) , we will display the 1 K results between those for = 1 and = 2.
Multiloop flow-The multiloop flow further includes the contributions from Γ (6) which are generated by vertex corrections. These can be ordered by loops, leading to the expansionγ r = ≥1γ ( ) r [15,16]. Here,γ (1) already includes the Katanin substitution to account for the self-energy corrections as above. The higher-loop terms, > 1, are determined bẏ where γr = r =r γ r . Equation (9a) with = 2 corresponds to the 2 flow, while the so-called center partγ In order to fully generate all parquet diagrams, the selfenergy flow also acquires a multiloop correction [15], whereγt ,C = ≥3γ ( ) t,C (see above). While not relevant for our AIM study, we note that additional approximations, such as the low-order expansion in form factors for the momentumdependence of the vertex functions, useful for reducing the numerical effort in treating lattice problems, require extra adaptations of the flow equations for the mfRG solution to converge to the PA [37,68].
The multiloop flow equation for the susceptibilities readṡ with the scale derivative of the two-particle irreducible verteẋ Fig. 2 for an exemplary diagrammatic representation. For more details and the equations for λ r , we refer to Refs. [9,36].
PA-In parquet approaches [1], a set of self-consistent equations for the self-energy and vertex is solved by iteration. First, Σ is related to F by the Schwinger-Dyson equation (SDE) Second, in the decomposition (5), the two-particle reducible vertices γ r are related to two-particle irreducible vertices I r by the Bethe-Salpeter equations (BSEs) In the PA, R 2PI = F 0 . Finally, the susceptibilities χ r can be directly deduced from F (and Σ via the propagators) by Here, λ r,0 are the bare three-point vertices encoding the relation of the composite bosonic degrees of freedom of χ r to the original fermionic ones.
In the parquet context, Eqs. (12)- (14) do not involve a scale parameter Λ. However, as they hold for any underlying bare propagator, they can also be applied when the bare propagator is G Λ 0 . These relations can then be used to derive the multiloop flow equations [9], and thus Eqs. (12)- (14) are fulfilled exactly in mfRG [16,36]. In other truncated schemes, they can be exploited as additional post-processing (PP) relations for computing (i) the self-energy from the SDE (12), (ii) the reducible vertices from the BSEs (13), and (iii) the susceptibilities using Eq. (14), instead of using the corresponding results of the flow. We recall that, for a generic truncated fRG scheme (including the standard 1 truncation), the PP values of Σ, γ r , and χ r differ from their counterparts obtained directly from the flow. In fact, the equivalence between the flowing and PP results for Σ, γ r , and χ r (upon convergence) represents, besides the independence from the choice of the cutoff function, a hallmark of the mfRG [36]. For this reason, we will also compute the PP results for Σ and χ r , and analyze their evolution with loop order.
QMC-Next to the fRG and PA described above, we use a state-of-the-art Quantum Monte Carlo [39] (QMC) solver to obtain numerically exact benchmark results of the AIM. We employ continuous-time QMC in the hybridization expansion (CT-HYB) [39] provided by the open-access W2DYNAMICS [69] package. Further details on the calculations are provided in Appendix B 3.

III. mfRG SOLUTION OF THE AIM
We now apply the mfRG, briefly summarized in Sec. II, to the half-filled AIM at the inverse temperature β = 10 and discuss the results. For details on the implementation, we refer to Refs. [6,36]. We just note here that, for the reducible vertices, we adopt the parametrization γ r = K 1r +K 2r +K 2 r +K 3r proposed in Ref. [6]. The K 1r and K 2 ( ) r functions with one and two frequency arguments, respectively, describe the highfrequency asymptotics, while the remaining full dependence at low frequencies is contained in K 3r . This reduces the numerical cost, allowing for the calculation of the vertices on a larger Matsubara frequency range (see Appendix B 1 for computational details). The (flowing) susceptibilities are conveniently extracted through χ r = −K 1r /U 2 .
We start the presentation of our numerical results by showcasing the central quantities of our study of the AIM, i.e., the self-energy Σ, the magnetic susceptibility χ m (= −χ ↑↓ a ), and the reducible vertex γ m of the impurity site in the magnetic channel, computed in the weak-coupling regime (U = 1) by means of all the approaches mentioned in Sec. II. Figure 3  (bosonic) Matsubara frequencies. The corresponding numerical data would also allow one to estimate important physical quantities (e.g., the quasiparticle mass renormalization and life time) relevant for the description of the Fermi-liquid state of the impurity problem [70,71] as well as to quantify the temporal fluctuations of the local magnetic moment on the impurity site [72][73][74].
Consistent with the small U value of these illustrative calculations, all approaches yield qualitatively the same behavior and deviations to numerically exact QMC data are hardly visible. In particular, we note that the converged mfRG solution (orange squares), perfectly matches the PA (dashed black line) for all quantities, ImΣ, χ m , and γ m (not shown). The results at the highlighted Matsubara frequencies are then used in the following Sec. III A for a quantitative study of the mfRG convergence as a function of loop order . There, we also showcase two hallmark qualities of the converged mfRG solution: (i) It is cutoff-independent, reflecting the fact that it reproduces the PA solution, which, as a self-consistent diagrammatic resummation, by construction is defined without reference to any cutoff. (ii) For quantities that can be computed either via their own RG flow equations or via PP relations, the results agree. (If the susceptibility flow is computed separately, and not via that of the K 1r part of the vertex, this requires to further use multiloop flow equations for the susceptibilities and the three-point vertices [9,36].) In Sec. III B, we extend this analysis to larger values of U .

A. Multiloop convergence to PA
In Fig. 4, we analyze in detail the loop convergence of the mfRG flow for U = 1. The four panels display both the flowing and PP results for ImΣ(ν = πT ) and χ m (ω = 0) as well as the flowing results for K 2m (ω = 0, ν = πT ) and K 3m (ω = 0, ν = ν = πT ), as a function of loop order obtained with the two cutoffs, i.e., the Ω-flow (red circles) and the U -flow (blue triangles). For comparison, the PA (black dashed line) and QMC (green solid line) solutions are also reported. One readily notices that the mfRG solution for both cutoffs converges to the PA for all considered quantities. Throughout the paper, the label '∞' refers to the infinite loop-order mfRG solution (see Appendix B 2 for its numerical definition). The high quality of the mfRG convergence can be appreciated by looking at the corresponding insets, showing the data restricted to higher loop orders. While the gray area in the main panels marks 1% deviation with respect to the PA, the blue area in the insets corresponds to 0.1%.
It is worth stressing that for some quantities and specific values of , the mfRG and PA solution may be accidentally close, e.g. the 3 Ω-flow result for χ m (ω = 0) or the 3 Uflow result for K 3m (ω = 0, ν = ν = πT ). Of course, this does not mean that the mfRG procedure has already converged at 3 : Full convergence implies the equivalence of mfRG and PA for all quantities and both cutoffs up to differences smaller than a given , e.g., here 0.1%. For the U = 1 calculations, this is clearly achieved for ≥ 8. Looking at the insets, the Ω-flow appears to converge systematically faster than the Uflow. We note that all U -flow results shown in the paper are obtained via a frequency extrapolation (see Appendix B 1), which is required to achieve the highly precise convergence to PA demonstrated in the inset.
Another important property of the converged mfRG solution is the equivalence of the flowing and PP results, shown both for ImΣ(ν = πT ) and χ m (ω = 0) in the upper panels of Fig. 4. Except for the 1 and 1 K results for the self-energy, the PP data (dotted lines with '×' or '+' symbols) are always found to be closer to the PA than the flowing data (for the susceptibility, this trend was previously reported in Ref. [37]). For both cutoffs, flowing and PP results agree with the PA Insets show a zoom for ≥ 4. The gray areas mark 1% deviation from the PA result, the blue ones in the insets 0.1%. The label '∞' represents the fully converged mfRG result. In this and similar figures below, the data points plotted between those at = 1 and = 2 represent the 1 K results (Katanin substitution).
for ≥ 8, highlighting the perfect convergence of the mfRG scheme in this parameter regime. The loop convergence can also be seen from calculations with a single cutoff, as there are no more changes larger than a small in all quantities when going from to + 1, and flowing and PP results agree with one another. Finally, let us note that adopting the PP procedure has also important implications for the fulfillment of sum rules, which are studied in Sec. IV A.
For values of U slightly larger than U = 1, the convergence behavior is qualitatively the same (see Fig. 17 in Appendix A for U = 1.5), albeit with increasing interaction, as expected, more loop orders are required to reach convergence.
For U = 2, the dependence on loop order is shown in Fig. 5. While the mfRG solution quickly approaches the PA for low , the path towards full convergence for higher becomes visibly slower as the curves describing the loop dependence of the mfRG calculations keep oscillating around the PA solution. The Ω-flow results are generally found to be more accurate than the U -flow data (note that for the Uflow at = 11, 12, no solution could be obtained; see also Appendix B 2). Yet, even with the Ω-flow, we did not reach perfect convergence up to = 40. Different from the situation at U = 1 and U = 1.5, the results obtained by PP do not show a clear improvement. Instead, they seem to follow a slightly different oscillation pattern, somewhat shifted  from the flowing data (see insets of Fig. 5). Further insight on the oscillations characterizing the mfRG convergence with increasing interaction can be gained from Fig. 6. Here, we show the relative difference between the mfRG results and the corresponding PA solutions for different values of U . By comparing the (flowing) results of the Ω-flow for different interaction strengths U = 1, 1.5 and 2, one notices the presence of "nodes" in the multiloop oscillations, i.e., of loop orders at which mfRG and PA yield numerically very similar results for the quantity under consideration. The location of these nodes, however, depends on the observable. [While, e.g., = 7 for χ m (ω = 0) is close to the PA for all values of U , for ImΣ(ν = πT ) this is not the case.] For larger interactions, the oscillations become stronger. Already for U = 2, the amplitude of the self-energy oscillations hardly decreases with increasing loop order, making a full convergence numerically challenging as discussed above. (The U -flow shows similar behavior, see Fig. 18 in Appendix A.) This effect gets even more pronounced for U = 3 displayed in the insets, together with U = 2 for comparison. There, higher loop orders, especially for χ m (ω = 0), yield a progressively enhanced deviation from the PA for increasing loop order. Therefore, we conclude that, within our current implementation and the given settings of the AIM, the mfRG loop resummation ceases to converge for U = 3. Such a lack of loop convergence serves as a built-in red-flag indicator that a parameter regime lies outside the zone of safe applicability of the approach. This outcome, however, is not entirely unexpected since, for the specific AIM considered, the interaction strength U = 3 already corresponds to the strong-coupling regime, where nonperturbative [51,[75][76][77]] divergences of two-particle irreducible vertices [78][79][80][81][82][83][84][85], which are-per construction-beyond the PA, were detected by means of QMC calculations [51,60]. More speculatively, one might then suppose a relation between the breakdown of the mfRG convergence and the entrance into the nonperturbative parameter regime, where the PA itself yields results significantly different from the exact solution [51]. In this respect, the oscillations of increasing size could be seen as a precursor of the breakdown of perturbative resummation schemes. We finally compare the results of low loop orders to the PA and the exact solution, as a function of U . For very low values of U , the deviations of mfRG and PA schemes from QMC can be qualitatively understood from general perturbation-theory considerations. Already for U > 1, however, the interpretation becomes more complicated, and the accuracy of the different schemes depends on the observable considered. Among the Ω-flow results up to U = 4 in Fig. 7, the plain 1 flow performs worst for all quantities. Comparing 1 K and the PA to the exact QMC for large U ≥ 2, we find the best results for ImΣ(ν = πT ) with 1 K , similar deviations for Im Σ(ν = πT ) − Im Σ(ν = 3πT ) with 1 K and the PA, and the best results for χ m (ω = 0) with the PA. For the physical interpretation of the strong-coupling regime, we refer to Ref. [51] and the corresponding supplemental material. There it was shown that both the PA and fRG schemes yield a qualitatively correct description of the magnetic channel; in particular, the proper behavior of χ m (ω = 0) as a function of T is found, reflecting the formation of a local magnetic moment and its screening. However, both methods fail in describing the associated suppressed fluctuations in the charge sector, which are heavily affected by the emergence of the local magnetic moment. Hence, at strong-coupling, the truncated fRG, mfRG or PA resummations of diagrams describe the formation of a local moment without the intrinsic physical implications onto the charge channel. This can be regarded [51,86] as an insufficient transfer of information between the magnetic and the charge sector, formally corresponding to the impossibility of generating the irreducible vertex divergences in these approximate methods.
On a more general perspective, we note that the loop convergence of the mfRG procedure is mostly controlled by the ratio between the local interaction U and other relevant energy scales of the system under consideration (e.g., in the case of the AIM: π∆ or the temperature T ) rather than by the ratio between the temperature and the Kondo temperature [51]. In future dedicated studies, it may be interesting to verify to what extent the grade of the loop convergence itself might be regarded as an additional independent marker of central physical aspects of the underlying exact solution of the problem.

IV. PAULI PRINCIPLE AND WARD IDENTITY
Both the Pauli principle and the WIs are fundamental features of the many-electron physics. They are deeply rooted in quantum mechanics and pose important constraints on manybody correlation functions. An exact solution must evidently obey all such constraints. In approximate treatments, however, their fulfillment is not guaranteed a priori. As mentioned in the Introduction, it is commonly reckoned [1] that approximate many-body approaches either obey sum rules imposed by the Pauli principle or satisfy WIs. Hence, fulfilling both the Pauli principle and the WIs would represent a specific hallmark of the exact solution. On a more formal level, a pertinent example of such a trade-off in the context of parquet-based approximations can be obtained by exploiting explicit relations between the self-energy and four-point vertices [7][8][9] in the parquet formalism.
In the following, we utilize our converged numerical results for the AIM to analyze, on a quantitative level, to what extent the Pauli principle and WIs are fulfilled for the important class of approximate many-body approaches ranging from the conventional fRG to the mfRG and PA. such a constraint affects the many-body correlation functions in several ways, e.g., through sum rules they must obey.
In this context, a relevant correlation function for the physics of the AIM is the equal-spin density-density susceptibility, Here, n σ = n σ , and T τ denotes (imaginary) time ordering (for brevity, we omit here the particle-hole channel label). This susceptibility is directly affected by the Pauli principle through the operator identityn 2 σ =n σ . Indeed, an evaluation at τ = 0 yields a value, which is fully determined by the single-particle expectation value n σ . Furthermore, as the equal-time correlator χ σσ (τ = 0) is identical to the sum over all its Fourier components χ σσ ω , the following sum rule [87] must hold: At SU(2) spin symmetry and half filling, the result is 1/4. For the purposes of the subsequent discussions, it is useful to elaborate on the quantum-field-theoretical relations which underlie Eq. (17). To this end, we recall that the Pauli principle can be translated from an operator identity ( to the crossing symmetry of four-point correlators. For illustration, let us briefly use a compact notation where all arguments of an electronic operator are summarized in a single index. Then, for G 1,2;2 ,1 . Furthermore, the susceptibility can be represented through (full) propagators G and the (full) two-particle vertex F by as illustrated in Fig. 8(a). The first term of Eq. (18) summed over all frequencies ω, i.e., taken at τ = 0, gives which yields already the entire sum rule [Eq. (17)]. Consequently, the vertex contributions must vanish when summed over all frequencies ω. This is indeed guaranteed by the crossing symmetry, as we show below. Consider the summed vertex contribution of Eq. (18), For F σσ , the vertex with equal spins on all legs, the crossing symmetry simply gives F σσ ν 1 ,ν 2 ;ν1,ν2 = −F σσ ν 1 ,ν 2 ;ν2,ν1 . After inserting this into Eq. (21), we relabel the summation indices according toω = ν − ν,ν = ν + ω: This reproduces the original expression for the summed vertex correction 1 Sum rule of χ σσ : Numerical results-As mentioned in Sec. II, there are two ways [10,36] of computing susceptibilities in fRG: (i) one can use Eq. (18) to obtain χ from Σ and F in a PP fashion, or (ii) one can deduce χ from its own flow equation. In the former approach the sum rule of χ σσ is fulfilled per construction, as long as the vertex used in the computation obeys the crossing symmetry, see Eqs. (20) and (23), while, in the latter scheme, this property is not guaranteed.
Not surprisingly, strategies (i) and (ii) then yield different results within 1 fRG (see Figs. 4 and 5), suggesting that the susceptibility computed from a 1 flow does not fulfill the sum rule. Indeed, one can easily convince oneself that the multiloop vertex corrections to the flow of χ σσ do not vanish when summing over all frequencies, cf. Fig. 9(a). On the other hand, we already noted that, for a converged mfRG calculation, both schemes of computing susceptibilities become equivalent [9,36]. Therefore, the sum rule of χ σσ will be consistently fulfilled, no matter the strategy employed.
On the basis of these considerations, we now turn to our numerical mfRG data. In Fig. 10, we show the loop dependence of 1 β ω χ σσ ω for the flowing susceptibility (obtained in the Ωflow) for different values of U . With increasing loop order, the fulfillment of the sum rule [Eq. (17)], indicated by a dashed black line, is approached. Altogether, we observe a similar behavior as in Sec. III: While, at low interaction values, the exact value is quickly reached, multiloop oscillations characterize the behavior at larger interaction (U = 2). Nevertheless, even for large U , the results at large are much closer to the fulfillment of the sum rule than the ones at low loop order. As for the PP susceptibility (not shown), we confirmed numerically that it fulfills the sum rule for all , consistent with the above explanations.
High-frequency asymptote of Σ: Formal aspects-Beside its natural link to the density susceptibility, the Pauli principle also affects the self-energy, albeit more indirectly. From the moments of the single-particle spectral function, known through expectation values of operators, one can determine the high-frequency expansion of the propagator G, and thereby of the self-energy Σ [87]. One finds Next to the constant Hartree shift U nσ, the 1/ν coefficient coincides with the r.h.s. of the sum rule for χσσ [Eq. (17)].
Indeed, Eq. (24) can be equivalently rewritten [88] as More insight about the quantum-field-theoretical relations underlying the asymptotic behavior of Σ can be gained from the SDE, see Fig. 8(b). To this end, let us replace the vertex by its bare contribution, F σσ 0 = −U , and use the first propagator in Eq. (26), G σ ν+ω , to factor out the dominant contribution for large ν ω, G σ ν+ω ∼ 1/(iν). The remainder is a GG bubble summed over both frequencies ω and ν . Hence, we find that the second-order contribution, already provides the correct asymptotic behavior (24). This is similar to the sum rule of χ σσ , where Eqs. (19)-(20) give the entire result, while the summed vertex corrections vanish [Eq. (23)]. Via Eq. (25), the same cancellation of vertex corrections occurs for the self-energy asymptote, as we explicitly show in Appendix C 1.
Within an fRG treatment, the standard flow equation for the self-energyΣ std in terms of the vertex F is in principle exact, as long as the exact vertex F is available. As this is almost never the case, the flowΣ std must be considered approximate. In mfRG, the multiloop corrections to the selfenergy flow [cf. Eq. (10)] effectively generate contributions tȯ Σ std which would require-when using the termΣ std onlyvertex diagrams beyond the PA (and thus beyond 1 fRG). Indeed, one can generally show that vertex diagrams beyond the PA (and thus beyond 1 fRG), such as the envelope diagram, do contribute toΣ std to order 1/ν in the large-frequency limit [cf. Fig. 9(b)]. Therefore, the Σ asymptote [Eq. (24)] is violated when using a 1 or multiloop vertex flow while keeping the standard self-energy flow. This problem is circumvented by including the multiloop corrections to the self-energy flow [16], which guarantee a perfect equivalence to the SDE and, thereby, that the correct asymptote will be restored.
High-frequency asymptote of Σ: Numerical results-In Fig. 11, we show (flowing) results for the asymptotic behavior of Σ as obtained from Ω-flow calculations for U = 1. The left panel displays νImΣ ν /U as a function of for a fixed, large value of ν l ≈ 19.16. At this frequency, νImΣ ν is expected to be slightly lower (in absolute value) than the corresponding asymptotic value of −1/4 for ν → ∞. The correct asymptotic description of the mfRG results (red circles) for large is demonstrated by their perfect match with the corresponding PA results, as the latter yield the correct high-frequency asymptotic by construction. As explained above, this would have not been the case without multiloop corrections to the self-energy flow. In fact, the gold pentagon line shows results which are obtained byΣ std without multiloop additions to the self-energy flow (these start at = 3) and notably deviate from the correct value.
The right panel shows the frequency dependence of νImΣ ν in a frequency window around ν l (ν l is represented by the vertical blue line). For fRG results at lower loop order, the high-frequency asymptote is incorrect, reflecting the fact that the SDE relation is not fulfilled. For the same reason, all approaches satisfying the SDE lie on top of each other, i.e., the PA (black dashed line), mfRG (orange solid line), and QMC (green line) [89] yield the correct high-frequency behavior of Σ. While the improvement of the high-frequency results is not monotonous for the lowest loop orders, we observe that rather accurate results are obtained already at the 3 level, where the first multiloop corrections to the self-energy flow appear. In this respect, it is also interesting to note that the standard self-energy flowΣ std provides a large-frequency asymptote in agreement with Eq. (25), but with χ σσ obtained by a one-loop flow and thus violating the sum rule [Eq. (17)]. We derive this result in Appendix C 1 and show explicitly which multiloop additions toΣ std contribute to the asymptote.

B. Ward identities
Formal aspects-The WIs play an essential role in the many-electron theory as they define how the information encoded in the continuity equations at a microscopical level is reflected onto response functions and macroscopic quantities. More specifically, a continuity equation is an operator relation of the form ∂ τρ = −[ρ,Ĥ]. Ifρ is a symmetry of the Hamiltonian, [ρ,Ĥ] = 0, thenρ is a conserved quantity, ∂ τρ = 0. In this case, the continuity equation describes a conservation law. However, even if this is not the case, continuity relations can be used for deriving relevant WIs, in particular when [ρ,Ĥ]albeit nonzero-yields a simple expression.
In practice, WIs can be derived for n-point correlation func-tions of arbitrary n. Ifρ and [ρ,Ĥ] involve n 1 and n 2 = n 1 + δn fermionic operators, respectively, then relates an n to an (n + δn)-point function. Typically, one mostly considers the WI connecting two-and four-point functions (i.e., the WIs ensuring the physical consistency between the one-and the two-particle description) and restricts oneself to the (local or global) charge or spin operators, substituting them forρ. A recent derivation, applicable to lattice and impurity systems, as well as references to prior work can be found in Refs. [90,91]. Here, we consider explicitly the (local) charge,ρ = σn σ , as done in several preceding works [17,92]. The resulting WI for the AIM, formulated in a way that allows for an optional SU(2) spin symmetry breaking (e.g. by a Zeeman field), reads We introduce the short-hand W L (ν, ω) for the left and W R (ν, ω) for the right side of the above equation, which is illustrated diagrammatically in Fig. 8(c). There, we use For our numerical results we exploit the SU(2) spin symmetry, which-together with the crossing symmetry-entails I ↑↑ t;ν 1 ,ν 2 ;ν1,ν2 = I ↑↓ t;ν 1 ,ν 2 ;ν1,ν2 − I ↑↓ a;ν 1 ,ν 2 ;ν2,ν1 .
Eventually, we briefly recall that one often refers to functional WIs, such as δΣ δG = −I t . These are a cornerstone of Φderivable approaches [93], where δΦ δG = Σ, and δ 2 Φ δG 2 = −I t . Since the functional derivative cannot be evaluated numerically, it mostly serves as a formal tool. However, by choosing a specific variation δG in the functional WI, one can derive more practical relations (as necessary but not sufficient conditions of the functional WIs). For instance, one can easily deduce Eq. (29) in the limit ω → 0 by varying G w.r.t. frequency (see Ref. [94] for a related treatment). Moreover, one can derive the standard fRG self-energyΣ std by varying G through the scale parameter [9]. Numerical results-Since the Pauli principle is preserved in the PA as well as (loop-converged) mfRG, one expects-on general grounds-these approximate schemes to violate the WIs to a certain extent. Arguably, the size of such violation should increase for increasing interaction strength, driven by the leading terms of the exact solution (where all fundamental relations are fulfilled) which are neglected in either approximate approach. Furthermore, it is known [10,18] that the 1 truncation leads to violations of the WIs. Katanin [17] proposed schemes to mitigate this deficiency. In particular, the 1 K flow is widely used and often argued to better fulfill WIs. However, no explicit numerical studies were presented thus far. Here, we intend to fill this gap and investigate quantitatively the fulfillment of WIs in fRG using our numerical results for the AIM. We focus on flowing (m)fRG results obtained with the Ω-flow, in order to avoid the frequency extrapolation required for the U -flow (see Appendix B 1).
We start with Fig. 12, where the top row shows W L (ν, ω) (squares) and W R (ν, ω) (hexagons) for ω = 2πT as a function of ν for U = 1, as obtained from the flow. We find that the  To better quantify the deviations between both sides of the WI, we focus on the quantity δW(ν, ω) = W L (ν, ω) − W R (ν, ω) at ω = 2πnT (n ∈ N) for two different choices for ν: In the first case, we fix ν to ν s = − n/2 2πT + πT , which gives the fermionic frequency closest to the symmetry axis ν = −ω/2, where the largest absolute deviations are found (e.g. ν s = −πT for ω = 2πT in Fig. 12, see also Fig. 13 discussed below). In the second case, we sum |δW| for ν in a finite frequency box. Specifically, we sum over 11 frequencies to the left and 11 frequencies to the right of the symmetry axis, adding also the contribution right at ν = −ω/2 if n is odd. In this way, we incorporate the behavior at larger frequencies, while avoiding numerical inaccuracies from the finite-frequency box effect of the high-frequency parametrization in our implementation [6] (see Appendix B 1). When comparing results for different transfer frequency ω = 2πnT , we divide by n to obtain more comparable results. The bottom row of Fig. 12 shows δW for the n = 1 data reported at the top. The plot confirms that, at weak-coupling, already the first multiloop corrections strongly improve the fulfillment of the WI. In particular, the minimal value for δW at ν = −πT (left panel) is found at = 2 and for |δW| summed over ν (right panel) at = 3. Hence, our U = 1 calculations show that the finite deviation from the exact fulfillment of the WI expected to occur in the loop-converged mfRG/PA results is notably smaller in comparison to 1 or 1 K , and that it quantitatively represents a marginal effect in the weak-coupling regime. This trend is also confirmed regarding relative deviations |δ r W| = |δW/W L |, as we explicitly show in Fig. 19 in Appendix A.
Next, we extend the analysis to larger values of ω = 2πnT and show in Fig. 13 loop-converged mfRG results for 1 ≤ n ≤ 5. The plot demonstrates that the mfRG data provide satisfactory agreement between W L (empty squares) and W R (filled symbols) for all values of ω and ν, and that the largest absolute deviation indeed occurs for ν around ν s , i.e., the frequency closest to the symmetry axis ν = ω/2 (see above). Figure 14 presents δW as a function of for n up to 40. Again, the fulfillment of the WI is slightly improved when going from 1 to 1 K and strongly improved starting from 2 , for all values of ω (confirmed also by Fig. 19 in Appendix A). However, the details in the change from = 2 to ∞ depend on ω. In general, we observe that the WI is better fulfilled for larger values of ω. In fact, a perfect match is given for ω → ∞ and → ∞, since the WI reproduces the SDE for ω → ∞ (see Appendix C 2), which is exactly fulfilled in mfRG and the PA. This can be clearly seen in both insets of Fig. 14. The inset of the right panel uses a logarithmic scale, where one can also spot the onset of oscillations in the multiloop convergence, in spite of their small amplitude.
Finally, we analyze the effects of the interaction strength, by progressively increasing its value up to U = 4. In Fig. 15, we examine δW for ω = 2πnT at n = 1 and n = 11, comparing results of (m)fRG flows at low loop order with the PA. At large interaction, the pure 1 flow is evidently unreliable, violating the WI with very large values of δW. The situation visibly improves in 1 K , 2 , and PA. In particular, for U ≤ 2, 1 K is farther off than 2 and PA. Interestingly, however, the 1 K deviations display a highly non-trivial behavior with increasing U -they are non-monotonous in the top left panel and have a decreasing slope in the other panels-and thereby yield comparatively small values of δW at larger U . By contrast, for the PA results, |δW| starts rather small but increases monotonously with increasing U . Overall, for intermediate to large values of U , it seems that 1 K provides the most accurate description of the WI at small frequencies (n = 1), while mfRG and the PA lead to a smaller violation of the WI for larger frequencies (here n = 11). Further details on the individual deviations of W L and W R are given in Appendix A.
As a last step, we compare the numerical deviations δW as   Fig. 15 for n = 1, but including data for very small interaction values U < 1, using a log-log scale. The straight lines correspond to f (x) = αx + d fits using the first data points between U = 0.1 and U = 0.5, yielding the exponents of the ∼ U α behavior.
a function of U focusing on small interaction values U < 1. Figure 16 shows δW, similarly as in Fig. 15, but on a log-log scale. Using a f (x) = α x + d fit, we extract the exponents of the deviations of the WI, δW ∼ U α , for the (m)fRG flow and PA scheme. Our analysis shows perfect agreement with the theoretical predictions of Ref. [17]: the 1 scheme displays deviations that grow with the third power of U (α ≈ 3, solid lines), and the 2 results are in agreement with a U 4 growth (α ≈ 4, dotted lines). The 1 K results at small U also manifest O(U 3 ) deviations. This is in agreement with the an-alytic arguments of Ref. [17] since, for the commonly used 1 K scheme, only part of the 2 corrections are included by substituting S →Ġ (as described in Sec. II). Hence, some terms violating the WI at O(U 3 ) remain, as seen in our numerical data in Fig. 16 (α ≈ 3, dashed-dotted lines). Note that the behavior at larger interaction values, as discussed above, is beyond the reach of the present analysis applicable at small values of U . Further, concerning the loop-converged mfRG/PA results, we find deviations of the WI, which behave as O(U 5 ) (dashed lines). In general, one expects the PA/mfRG schemes to deviate from the exact solution as O(U 4 ). However, at half filling, the combination of the particle-hole symmetry and spin symmetry of our problem causes the contributions to the WI from the forth-order "envelope" diagrams to exactly cancel, as we show explicitly in Appendix C 3. For completeness, we also note that the same behavior as in Fig. 16 is found for other frequency choices as well (e.g. for n = 11 used in the lower panel of Fig. 15).

V. CONCLUSION AND OUTLOOK
We investigated several essential features of the recently introduced mfRG approach by performing a quantitative study of the particle-hole symmetric AIM for different coupling strengths. As the numerical implementation of the mfRG applied to the AIM does not require additional algorithmic approximations (such as the form factor expansion used for the Hubbard model [36,37]), we were able to demonstrate how the precise convergence of the mfRG series to the corresponding PA results is readily obtained in the entire weak-to intermediate-coupling regime. A thorough inspection further confirmed the pivotal features of a converged mfRG solution, i.e., its independence of the specific RG cutoff adopted as well as the equivalence between flowing and post-processed results. Hence, in the parameter regimes where a fast loop convergence of the mfRG is found, the application of this method offers potential advantages w.r.t. to the full iterative solution of the PA through the intrinsic flexibility of the underlying fRG framework.
By increasing the value of the electronic interaction, we studied the oscillatory behavior emerging in the loop dependence of the mfRG series, which eventually hinders the convergence to the PA solution in the strong-coupling regime. Interestingly, the parameter region where a multiloop convergence could not be achieved appears roughly to match the one in which previous Quantum Monte Carlo studies [51,60] have shown an explicit breakdown of perturbative resummations to occur at the two-particle level. In this respect, the strong oscillatory behavior of the non-converging mfRG series could be plausibly regarded as a further hallmark of the nonperturbative [51,75,76,78] parameter regime, where significant physical differences between the PA and the exact solution of the AIM are found [51].
The numerical data obtained in the region of proper convergence of the mfRG algorithm were then used for a quantitative investigation of the fulfillment of fundamental features of the many-electron problem, namely those linked to (i) the Pauli principle and (ii) the WIs. For (i) the Pauli principle, we observed a sizable violation of sum rules in the conventional 1 fRG results, which gets systematically reduced by increasing the loop order. This is consistent with the fact that mfRG converges to the PA solution, and that the PA obeys the Pauli principle by construction, realized through the crossing symmetry and two-particle self-consistency. We also note that the indirect effects of the Pauli principle on the high-frequency asymptotic behavior of one-particle quantities are only recovered by including the multiloop additions to the self-energy flow, which start from the third loop onwards. For (ii) the WIs, these are generally neither fulfilled in fRG nor in the PA. For weak to intermediate coupling, our results demonstrated that adding higher-loop terms systematically reduces the overall violation of WIs. In particular, while a first improvement can be already observed by including the one-loop Katanin (1 K ) substitution, higher loop orders and the PA yield quantitatively much smaller deviations. By increasing the interaction, however, the situation becomes more complex. Going beyond the 1 fRG level, whose description of the WIs is largely unreliable, we find that 1 K mitigates most efficiently the WI violations at low frequencies, while higher-loop mfRG and the PA yield better results for large frequencies. This is consistent with our observation that the WI reproduces the SDE for ω → ∞. Additionally, we confirmed the predictions of Ref. [17] for the asymptotic weak-coupling behavior of the WI deviations as a function of U for the 1 and 2 scheme. Our numerical results for the mfRG/PA scheme revealed a O(U 5 ) deviation, smaller than the expected O(U 4 ), which we showed to be related to the particle-hole and spin symmetry used in our computations.
The insights gained in our study, which might be extended in the future to other regimes (e.g., out of half filling, and/or in the presence of a magnetic field) and more complex systems, are important for several reasons. On the one hand, they improve the understanding of the convergence of the mfRG procedure, whose relevance extends to more complex contexts than the basic AIM considered here. Such insights may be particularly important if the mfRG is used to include nonlocal correlations on top of the DMFT solution of strongly correlated lattice problems, thus extending the DMF 2 RG algorithms beyond the conventional (1 ) fRG used so far [52,53]. In that context, the mfRG might offer important advantages over corresponding parquet-based implementations. In contrast to the latter, the mfRG flow does not rely on the numerical manipulation of two-particle irreducible vertex functions, which display multiple divergences in the intermediateto-strong coupling regime of different many-electron models [51,60,[78][79][80][81][82][83][84][85]. This should allow the circumvention of several of the problems faced by parquet-based DMFT extensions [42] constructed upon such potentially diverging irreducible vertices, such as parquet DΓA [4,95] or QUADRILEX [96].
On the other hand, the possible relation of the loop convergence properties in mfRG with the breakdown of the perturbation expansion might have interesting theoretical and algorithmic implications, calling for an extension of our study to more complex physical situations than those considered here.
Together with our precise analysis of the fulfillment or violation of sum rules and WIs, this might shed new light on fundamental aspects of the many-electron theory and help to further develop refined calculation strategies for the most challenging parameter regimes.

VI. ACKNOWLEDGMENTS
The authors thank C. Eckhardt

APPENDIX
In the Appendix, we provide additional results, details on the numerical treatment as well as diagrammatic derivations, in order to specify our approach and further support the messages of the main part. The additional results are in Appendix A, mainly focused on the U -flow and the fulfillment of the WI. Details on our numerical approach, especially the dependence of different quantities on the number of Matsubara frequencies included in the computations, are discussed in Appendix B. Finally, we give the diagrammatic derivations of several relations used in Sec. IV in Appendix C.

Appendix A: Additional results
In Fig. 17, we report the results for U = 1.5 (β = 10, half filling), which were anticipated in Sec. III B. For this parameter set, too, the mfRG scheme converges perfectly in loop order. For ≥ 15, both regulators lead to identical results for all quantities, and the PP (dotted lines with '×' or '+' symbols) and flowing data coincide. As stated in the main text, no qualitative difference in the convergence behavior is observed, apart from the fact that, for U = 1.5, more loop orders are necessary to reach it.
In Fig. 18, the relative comparison between U -flow results and the PA for U = 1, 1.5, 2 is shown in the same fashion as in Fig. 6 for the Ω-flow. While there is no qualitative difference, quantitatively the U -flow shows larger relative differences with respect to the PA. Note that we were unable to converge the U -flow calculation for = 11, 12; see also Appendix B 2.
Finally, we add further analyses on the fulfillment of the WI, namely (i) on the relative deviations for the cases discussed in the main text, and (ii) more details on the devia-tions as a function of U in the different approaches. Concerning (i), Fig. 19 is a combined plot of Figs. 12 and 14 of the main text, but instead of δW, we show |δ r W(ν, ω)| = |W L (ν, ω) − W R (ν, ω)/W L (ν, ω)|. In the top row, |δ r W| is shown for ω = 2πT , similarly as in the top row of Fig. 12. Note that the y-axis is cut at |δ r W| = 100% in order to present the behavior of |δ r W(ν, ω = 2πT )| for the various values of ν with sufficient resolution. The reason for the peak of |δ r W| at one specific Matsubara frequency is the sign change (and hence the closeness to zero) of |W L |. The bottom panels and the corresponding insets show the relative deviation |δ r W(ν, ω = n2πT )| for ν = ν s (left) as well as for an averaged sum over a finite frequency box (see main text for both). Due to the averaging effect of the factor 1/N ν in 1 Nν ν |δ r W(ν, ω = n 2πT )|, where N ν is the number of elements summed over, the factor 1/n used in Fig. 14 is omitted. In general, Fig. 19 confirms the trend described in the main text. One notices how the increase of the loop order leads to a reduction of the relative deviations for all frequencies ω and ν. As pointed out in Sec. IV, the WI is exactly fulfilled for the mfRG/PA solution at n → ∞, which is also confirmed in Fig. 19 (see insets). An important difference to Fig. 14 is that for the 1 , 1 K and 2 scheme, δ r W(ν = ν s , ω) is roughly constant, or even grows as n is increased. This reflects the fact that these approaches do not respect the SDE, and hence do not fulfill the WI exactly for n → ∞.
Regarding (ii), we use the numerically exact QMC solution (fulfilling the WI) as a reference and compare W L and W R obtained by fRG/PA for ν = ν s (see main text) individually with the QMC result.  Fig. 20: on the one hand, there are the deviations of the fRG/PA results from the numerically exact QMC results, on the other hand, the fact that the fRG/PA results do not fulfill the WIs, and are hence not conserving. As discussed in the main part in Fig. 7, the deviations between PA/fRG calculations and the QMC results grow with U , which can also be seen in Fig. 20. The solution of a conserving approximation would show this deviation, but would not show a difference between the left and right side, i.e., the full and the dashed lines would coincide. Hence, it is not the value on the y-axis itself, but the difference in the deviation of ∆W L and ∆W R , which turns out to be instructive. As can be seen in Fig. 20, for most cases, it is the right side of the WI that deviates more from the QMC solution, the Ω-flow 1 K -results for n = 11 represent the extreme case. While in the PA, the 1 and 2 results show a steadily growing difference between the solid and the dashed line, the situation is less monotonous for the 1 K approach. From its data for n = 1 (top), one clearly notices the change in behavior as ∆W R changes sign, leading to the sign change of δW 1 K L (ν = −πT, ω = 2πT ) seen in Fig. 15.  Appendix B: Details on the numerical approach

fRG and mfRG calculations
Our fRG, mfRG, and PA computations for the AIM are based on the implementation used in Refs. [6,36]. As stated in the main text, we employ the following parametrization of the reducible vertex functions [6] γ r = K 1r +K 2r +K 2 r +K 3r . The high-frequency asymptotics are included in the K 1r and K 2 ( ) r functions with one and two frequency arguments, respectively. The remaining full frequency dependence, which has a relevant contribution at low Matsubara frequencies, is contained in K 3r . These contributions increase with increasing interaction values, and it is hence necessary to extend the frequency box, i.e., the number of frequencies where the full frequency dependence of K 3r is taken into account. In Ta-ble I, we provide the number of positive fermionic frequencies of K 3r , N f+ , for different approaches and values of U . The parameter N f+ also dictates all other frequency ranges in the same way as detailed in Ref. [6]. Outside the finite frequency box, the K 3r functions are set to zero, which is the core of the high-frequency asymptotics approximation. While this affects all quantities calculated with the different approaches, the difference in the results observed by comparing computations with different box sizes is negligible for the Ω-flow and PA. By contrast, for the U -flow, an extrapolation in N f+ is necessary, as detailed in the following subsection.  Fig. 12, as a function of ν. The y-axis is cut at 100% to provide enough resolution for |δrW(ν, ω = 2πT )| at the various values of ν. Bottom: As Fig. 14 but showing |δrW| instead. In the right panel, the normalizing factor 1/n of the main text is replaced by 1/Nν , where Nν is the number of frequencies summed over (see text). All quantities are given in percent [%]. In order to achieve agreement between the ∞-loop mfRG solution using the U -flow and the corresponding PA result to the precision chosen in the main part of the paper (0.1% in the insets), it is necessary to perform a frequency extrapola-tion. To this end, several calculations for the same parameter set are performed with different sizes of N f+ (see Table I). In Fig. 21, we showcase this for U = 1 and β = 10, i.e., the case discussed in Sec. III A. The open blue symbols represent the results for ImΣ(ν = πT ) as obtained by different ∞loop U -flow calculations, plotted as a function of 1/N f+ . For comparison, the results of corresponding PA calculations with different box sizes are shown as open black symbols, which hardly display any dependence on N f+ at this scale. Using a f (x) = A + B/x fit (blue line), we obtain the extrapolated value (filled blue triangle), which lies on-top of the PA result for N f+ = 32 (dashed black line). For comparison, we also plot the result of an ∞-loop Ω-flow calculation (open red circle) using N f+ = 32, which highlights that, for the Ω-flow, no frequency extrapolation is necessary to reach agreement with PA at this precision, as stated above.  results shown), where a A + B/x + C/x 2 fit turned out to be the best choice.

mfRG calculations
In this part of the Appendix, we provide further details on our multiloop calculations. In particular, we specify how the ∞-loop mfRG solution is obtained and concisely discuss the iteration of Σ.

a. ∞-loop mfRG solution
At each step of the fRG flow, the changes in all quantities for all Matsubara frequencies when going from to + 1 are measured. As soon as the relative (absolute) changes are lower than a given , in our case 10 −5 (10 −7 ), the calculation of higher loop orders is stopped. This speeds up the computation especially at the beginning of the flow, where usually a low loop order is sufficient; for more details on this, see Ref. [97]. While for obtaining the solution of loop order , the multi-  loop calculation is stopped at this specific , it is continued until the changes are smaller than to calculate the ∞-loop order solution. In Table II, we provide the actual number of loops needed ( max ) to obtain the ∞-loop order solution for the different flows and parameter sets.
b. Iteration of Σ Part of the mfRG scheme is also the iteration of Σ at each step of the flow [9,15,16]. The effect of these self-energy it-erations was analyzed in great detail in Ref. [37]. Throughout our calculations, their impact proved to be small, e.g., comparing χ m (ω = 0) with and without the iteration of Σ for the Ω-flow at U = 1 leads to a difference of O(10 −4 ). In Table II, we provide the necessary number of iterations for the ∞-loop mfRG solution (N Σ-iter ) to arrive at differences smaller than (given above, see Appendix B 2 a) when comparing iteration i with i+1. For all other loop orders, the same condition was used. As it turns out, the number of necessary iterations proved to be very similar, except for the = 11, 12 U -flow calculations for U = 2. There, the number of required Σ iterations increased considerably, preventing our numerical calculation from converging in a reasonable amount of time.
Lowering the maximum number of Σ iterations did not allow for obtaining a converged result, as the adaptive solver used for our computations did no longer converge in this case.
For completeness, Table II also lists the number of Runge-Kutta integration steps in Λ during the fRG flow for both regulators (N step ), as well as the number of PA iterations (N PA-iter ). Note that, since the calculations for U = 3 were numerically very costly, as they required a large frequency box for K 3r , we restricted the number of iterations for the Ω-flow computations shown in the main part to 3.

QMC calculations
As stated in the main part, we employed the W2DYNAMICS [69] package (version 1.0.0) as a continuous-time QMC [39] solver. We used the default sampling method for all calculations shown apart from the data for Fig. 11. There, we performed Worm sampling [98,99] computations with symmetric improved estimators [100] instead, which reduces the high-frequency noise. While we used about 2000 CPU hours for the former computations, the Worm sampling calculations were done using up to 25000 CPU hours. In this section, we will derive relations between the highfrequency asymptote of Σ and the sum rule of χ σσ . First, we will show that the two are directly related through the SDE in parquet-type approaches. Then, we move on to fRG flows. We will show that the standard self-energy flow also relates the Σ asymptote to the susceptibility sum rule, with χ σσ given by its one-loop flow. Since the latter does not fulfill the sum rule, the former violates the exact asymptote. Both the sum rule and the asymptote are fulfilled in multiloop fRG. We will show which terms of the multiloop corrections toΣ complete the relation, so that the Σ asymptote is determined by χ σσ obtained in a multiloop flow. The entire derivation will proceed diagrammatically.

Connection through the SDE
In Fig. 22(a) we recall Eq. (18), which expresses χ σσ through a GG bubble and corrections in terms of the full fourpoint vertex F . The vertex F is contracted by pairs of propagators on both sides. Therefore, one can also express χ σσ through a (full) three-point vertex λ on either the left or the right side, as illustrated in Fig. 22(b). The vertex λ is particularly useful when considering F in the limit of large fermionic frequencies.
Indeed, to find the self-energy asymptote, we will consider a large fermionic frequency ν. In Figs. 22(c,d), we show which diagrams of the vertex F , carrying ν on the external legs marked in red, remain nonzero in the limit ν → ∞, i.e., which diagrams are independent of ν. We use the symbol '=' for that purpose, signifying equality up to O(1/ν). To have nonzero contributions when ν → ∞, the red (amputated) external legs must directly meet at a bare interaction vertex. This is clearly fulfilled for F = F 0 , but there can also be arbitrary vertex corrections after the two red legs have met. If ν is on the lower two legs [ Fig. 22(c)], such corrections are a subset of the vertex γ t , reducible in transverse (vertical) particle-hole lines. If it is on the left two legs [ Fig. 22(d)], the corrections belong to γ a , reducible in antiparallel (horizontal) lines. The bare vertex and and the corrections are summarized by the three-point vertex λ. To see this, one may insert the BSEs for γ t/a , connecting the irreducible vertices I t/a to the full vertex F . Since I t/a are irreducible in their respective channels, they collapse to F 0 in the limit ν → ∞, and one obtains λ similarly as in going from Fig. 22(a) to Fig. 22(b). Now, by means of the SDE (26), the self-energy (minus its static Hartree part) is determined by the vertex F connected to three propagators, as we recall in Fig. 22(e). To find Σ to first order in 1/ν, we need F to zeroth order. We choose to transport ν through the propagator at the bottom. Then, we can directly use the relation in Fig. 22(c) to replace F by λ up to corrections O(1/ν 2 ) (signified by the symbol '='). Using Fig. 22(a), we obtain χ σσ through λ. The last step is similar to Eq. (27): Take the red propagator as Gσ ν+ω . For ν ω, we can replace Gσ ν+ω by 1/(iν) up to corrections O(1/ν 2 ). This leaves χ σσ ω summed over all ω, and, with a prefactor U 2 from the two interaction lines, we obtain Eq. (25) for Σσ.

Standard self-energy flow
Next, we turn to fRG flows. The standard self-energy flow is given byΣ σ std = −F σσ ·S σ , where '·' denotes the contraction of the top two vertex legs by the following propagator, S is the single-scale propagator, and a sum over σ is understood. For formal derivations, it is helpful to analyzeΣ σ std by means of its equivalent skeleton version [9],Σ σ 1 = −I σσ t ·Ġ σ , illustrated in Fig. 23(a). As before, a line with a doubled orthogonal slash denotesĠ, and dashed dark and light colors indicate a summation over spin.
As mentioned previously,Σ 1 is exact only for an exact vertex, which is not available in practice. Instead, we will consider the much more relevant case of a vertex obtained in the PA or, equivalently, a multiloop flow. In this case,Σ 1 is approximate. We will show that it generates a high-frequency asymptote of similar type as the exact relation Fig. 22(f), but with χ σσ obtained by its (approximate) one-loop flow. The connection from the general, Λ-independent statement Fig. 22(f) to an fRG flow is made by taking the scale derivative ∂ Λ on the entire equation. In this way, ∂ Λ is subsequently applied to the trivial Hartree part, to the red propagator alongside χ σσ , and finally to χ σσ itself. Indeed, we will precisely find such a structure, where the derivative ∂ Λ χ σσ is approximated byχ σσ 1 , see Fig. 23(b). The one-loop flowχ σσ 1 is given by the first summand of Fig. 23(c). (The long double slash denotes a differentiated two-particle propagator,Π =ĠG+GĠ.) The multiloop corrections toχ σσ , which are compactly encoded in the second summand of Fig. 23(c) and will be considered more closely in the next part, restore equivalence to the general susceptibility-vertex relation shown in Fig. 22(a).
To derive Fig. 23(b), we start from I t = F 0 + γ a + γ p in the PA. The bare vertex F 0 immediately gives the differentiated Hartree part as the first summand of Fig. 23(b). From Fig. 23(d) onward, we analyze the effect of γ a/p using their BSEs. The analysis is slightly more complicated than in Fig. 22(e): There, we had just a single vertical interaction line; now, we have two spin-dependent vertices, where same-spin propagators can meet both vertically and horizontally.
In Figs. 23(d,e), we insert the BSE of γ σσ a , with a summation on the spin carried byĠ σ . This gives three terms: (i) γ σσ a , where the antiparallel two-particle propagator Π a necessarily has two opposite spins; γ σσ a , where Π a has (ii) both spins equal to σ and (iii) both spins equal toσ. Cases (i) and (ii) are contained in Fig. 23(d), with a spin sum encoded in the dashed colors. Regarding ∼ ν −1 contributions, Fig. 23(d) contains all diagrams where the lower two legs of I a and F directly meet at vertical interaction lines and the large frequency is transported through the bottom propagator. Since both I a   FIG. 23. (a) The standard self-energy flowΣ std ≡Σ1 in skeleton form [9]. (b) In the PA, it yields asymptotic contributions of the same structure as Fig. 22(f), whereχ σσ is approximated by its oneloop flow. (c) One-loop ( = 1) and multiloop terms ( > 1) forχ σσ [9]. To derive the relation in (b), we insert the BSEs for γa, γp as part of It; (d,e) concerns γa, while (f,g) concerns γp. and F contain F 0 +γ t , their ∼ ν 0 contributions are expressed through λ according to Fig. 22(c). Proceeding with case (iii), Fig. 23(e) contains all diagrams where the left (right) legs of I a (F ) directly meet at horizontal interaction lines and the large frequency is transported through the top propagator. While F contains both F 0 and γ a , I a contains only F 0 . Hence, their ∼ ν 0 contributions are expressed through λ and a bare interaction line, respectively, according to Fig. 22(d).
We continue with γ p and insert in Fig. 23(f) the BSE of γ σσ p (Ġ is in light color), where the parallel two-particle propagator Π p is summed over both spins (and thus the typical prefactor 1/2 is kept). Since both I p and F are crossing symmetric, contributions stemming from vertical and horizontal interaction lines enter equivalently. Indeed, in the first (second) summand of Fig. 23(f), the red propagator passes by vertical (horizontal) interaction lines. Since both I p and F contain F 0 +γ a +γ t , we replace their ∼ ν 0 contributions by λ using Fig. 22(c,d), and we end up with two equivalent terms. In Fig. 23(g), we insert the BSE of γ σσ p , where Π p must also carry spins σ (the prefactor 1/2 remains). Again, the red propagator can pass by vertical and horizontal interaction lines, and we get two equivalent terms expressed through λ.
Finally, we see that Fig. 23(e) gives the second summand of Fig. 23(b) [by means of Fig. 22(b)], and the sum of Figs. 23(d,f,g) reproduces the first summand of Fig. 23(c). This yields the last part of Fig. 23(b), thus concluding the derivation.  Fig. 23(b), the asymptote ofΣ2 must be related to the multiloop corrections ofχ σσ . To show this, we split the contraction of γt ,C into four summands: (c) γ σσ a,C · G σ already yields the desired expression; (d) γ σσ p,C ·G σ vanishes up to corrections O(1/ν 2 ); (e) (γ σσ a,C +γ σσ p,C )·Gσ cancel to that order, as can be seen after factoring out Gσ ∼ 1/(iν) for the first and Gσ ∼ 1/(−iν) for the second summand.
It is interesting to analyze how Fig. 24(b) comes about. Through the spin sum and the composite nature ofγt ,C ,Σ 2 has four contributions, stemming fromγ σσ a,C ,γ σσ p,C ,γ σσ a,C , anḋ γ σσ p,C . We will show that the first term already gives the desired result in Fig. 24(b). Up to corrections O(1/ν 2 ), the second term vanishes while the last two terms cancel.
Insertingγ σσ a,C , the only way to get ∼ ν −1 contributions is to transport the large frequency through the loop propagator at the top, marked red in Fig. 24(c) (note that both two-particle propagators Π a are summed over spin). Further, all red lines must directly meet at (horizontal) interaction lines. Hence, the four-point vertices thus yielding Fig. 24(b).
For the remaining terms, one immediately sees in Fig. 24(d) thatγ σσ p,C has no ∼ ν −1 contribution: The external legs and the loop propagator would need to directly meet as two out-going (in-going) lines at a bare interaction line of the left (right) vertex. However, they all have the same spin, and the bare interaction requires in-and out-going lines to have opposite spin. Next, the opposite-spin contributionγ σσ a,C +γ σσ p,C is shown in Fig. 24(e). By choosing fixed spin labels for the two Π p enteringγ p,C , we eliminate the typical prefactor (1/2) 2 . The upper loop propagator carrying the ν dependence goes in opposite directions for the first compared to second summand. Hence, after factoring out the dominant 1/(iν), we get opposite signs for the ∼ ν −1 contributions between the a and p channel. The remaining part for both is summed over all internal frequencies, including ω, as indicated by the closed wiggly line. Their sum cancels, as can be checked explicitly at low orders. Note that, for this to work, one needs the same number of diagrams in γ σσ a and γ σσ p at each interaction order, as is indeed the case [101].

Deriving the SDE from the WI
The WI (29) relates a difference of self-energies, Σ σ ν+ω − Σ σ ν , to a vertex contracted by a combination of propagators. For infinitely large ω, while ν remains finite, the first selfenergy simplifies to its static value, Σ σ ν+ω → U nσ, and we thus obtain a relation for Σ σ ν alone. This relation is precisely the SDE (26), as we show now.
We start by restating Eq. (29) in the form Here, we labeled I t by only three frequencies, chosen in the natural parametrization of the t channel, with the bosonic frequency ω as a superscript and the two fermionic frequencies ν and ν as subscripts. We also introduce a diagrammatic representation of the WI that is slightly different from Fig. 8(c): In Fig. 25(a), we have the difference in self-energies on the left and a difference of the vertices, each contracted by a different propagator on the right. Indeed, each vertex is contracted by only the propagator corresponding to the long line. All the short, external legs are amputated; they do not contribute to the diagram. In particular, the short wavy line only serves to ensure energy conservation for each vertex; it does not enter the equation itself. We recall that dark and light colors distinguish the two spin species; dashed lines with dark and light colors symbolize a sum over spin.
If we take the limit ω → ∞ in Eq. (C1) or Fig. 25(a), the l.h.s. simplifies to −Σ, whereΣ σ ν = Σ σ ν − U nσ is the self-energy without its static Hartree part. In this limit, the r.h.s. simplifies as well. First, we express I t , the vertex irreducible in the t channel, as a sum of the fully irreducible vertex R 2PI and the vertices reducible in the complementary channels, γ a and γ p . Fully irreducible vertex diagrams beyond the bare vertex, F 0 , decay in all frequency arguments; therefore, lim ω→∞ R = F 0 . However, F 0 makes no contribution to Eq. (C1), as it is frequency independent and thus leads to cancellation in the ν sum. In contrast to R 2PI , the reducible vertices γ r have specific contributions that are independent of certain (fermionic) frequencies. By substituting γ a +γ p for I t in Eq. (C1), we get (γ σσ ;ν −ν a;ν+ω,ν +γ σσ ;ν+ν +ω p;ν+ω,ν )(G σ ν +ω −G σ ν ).
where we relabeledν = ν +ω in the p channel. This relation is the first equality in Fig. 25(b). The symbol '=' here means that both sides agree up to O(1/ω), i.e., they are equal in the limit ω → ∞.
For the rest of the analysis, we refrain from spelling out the equations and proceed diagrammatically. In Fig. 25(b), lines in red colors (dark and light for the two spins) carry the large frequency ω. All of these are amputated external legs, for, otherwise, the result would vanish in the limit ω → ∞. This means that only those diagrams of γ p and γ a contribute where the red legs directly meet at the same bare vertex. Thereby, ω is transferred without entering an actual propagator, and the result is completely independent of ω. We can gather all those diagrams by inserting the BSEs for the reducible vertices. This is done in the second equality of Fig. 25(b). The first two and last two summands per row differ by the choice of spin in the propagator loop on top of γ a/p . In the BSE for γ σσ p , we fixed the spin σ at the bottom propagator, thus eliminating the prefactor 1/2.
By virtue of the BSEs, we have two out-going red legs attached to I p , and an in-and an out-going red leg attached to the left of I a . Since I p and I a are irreducible in parallel and antiparallel lines, respectively, the only diagram for each that allows the red legs to meet directly is the bare vertex F 0 . Furthermore, as F 0 is only nonzero between different spins, the result collapses to the three contributions (without any spin summation) shown in Fig. 25(c). At this point, the red lines meet at a bare vertex, and the ω dependence (and thus the wiggly line) can be simply removed. We see that the first two terms in Fig. 25(c) cancel. It remains to use the crossing symmetry of F to transform the last summand of Fig. 25(c) into the expression of Fig. 25(d). The latter is precisely the SDE (the Hartree term is absorbed inΣ) in the form known from Eq. (26) and Fig. 8(b).
It is no coincidence that the first two summands of Fig. 8(c) canceled, and the nonzero contribution toΣ is the one from the equal-spin vertex γ σσ a coming from I σσ t . In fact, the WI (C1) also holds without spin sum, 2PI does not contribute to the WI. This is because pairs of diagrams cancel in the difference of Eq. (29), as apparent after a suitable transformation ν → ν in, say, the second term. (d) Cancellation of diagrams (i) and (ii) from (a) inserted in the first and second term of the WI. (e) Same for diagrams (iii) and (iv) from (b).
For convenience, we check this explicitly at second-order in U in Fig. 26. Equation (C3) can be found by deriving the WI not only using the local charge operator,ρ 1 = σn σ , but also the local spin operatorρ 2 = σ τ z σσnσ , where τ z is the third Pauli matrix. For the latter, the WI reads

Envelope diagrams in the WI
The vertex in the PA deviates from the exact vertex starting at fourth order in the bare interaction U . Through the SDE, the self-energy is exact up to order U 4 , while errors start at order U 5 . A priori, one thus expects the PA to violate the WI to order U 4 , as induced by the missing U 4 vertex diagramsthe so-called envelope diagrams. However, it could also happen that this specific class of vertex diagrams does not con-tribute to the WI, i.e., that the envelope diagrams cancel out in Eq. (29). Indeed, numerically, we found the PA to violate the WI to order U 5 instead of U 4 . In this section, we show analytically that in the special case of particle-hole and spin symmetry, the envelope diagrams cancel in the WI. Figures 27(a) and 27(b) show the envelope diagrams with different and equal spins (distinguished by light and dark colors) on the external legs, respectively. (For brevity, we collapsed the interaction line to a dot.) These diagrams can also be found in Figs. 14 and 15 of Ref. [40]. Together, they form the fully irreducible vertex R 2PI at fourth order in U ; R (4)σσ 2PI has two and R (4)σσ 2PI four diagrams, as enumerated by Roman numbers. Figure 27(c) states that neither R (4)σσ 2PI nor R (4)σσ 2PI contribute to the WI given particle-hole and spin symmetry.
Indeed, using G σ ν = −G σ −ν and G ↑ ν = G ↓ ν , one can always find pairs of envelope diagrams that cancel in the subtraction inherent to the WI. This cancellation becomes apparent after suitably transforming the summation frequency in, say, the second term of the WI (thus changing ν to ν ). In Figs. 27(d) and 27(e), we establish the cancellation by explicitly writing frequency labels on all internal lines. with diagram (i) and (ii) in the first and second term of the WI, respectively. The same set of frequency labels occurs in both terms, differing at most by minus signs. Both terms have a total of three global minus signs in their frequency labels; using G σ ν = −G σ −ν , these minus signs can be pulled out of the equation. One obtains a mathematically identical expression for both terms and thus a vanishing difference. The case of R (4)σσ 2PI with diagram (ii) in the first and diagram (i) in the second term proceeds analogously. Indeed, one can transform one case into the other by flipping the arrows on the two horizontal and the two diagonal lines (accordingly changing ν → −ν on their frequency labels) and suitably changing the spin labels. The cancellation works just as before. Further, the cancellation also works for diagrams (i) and (ii) of R (4)σσ 2PI (again for both orders). In this case, one must further invoke spin symmetry, G ↑ ν = G ↓ ν . Finally, Fig. 27(e) treats the case of R (4)σσ 2PI with diagram (iii) and (iv) in the first and second term of the WI, respectively. The argumentation is the same as before. Thanks to the transformation ν → ν in the second term, one has a mathematically identical expression for both terms and thus a vanishing difference. Again, interchanging the role of diagrams (iii) and (iv) in the two terms of the WI merely amounts to flipping the arrows; here, this affects the two vertical and the two diagonal lines, all of which have spinσ. As apparent from Fig. 27(e), no minus signs are involved, and matching pairs of frequency labels also have the same spins. The cancellation of diagrams (iii) and (iv) thus works also in the absence of particle-hole or spin symmetry.