Single boson exchange representation of the functional renormalization group for strongly interacting many-electron systems

We present a reformulation of the functional renormalization group (fRG) for many-electron systems, which relies on the recently introduced single boson exchange (SBE) representation of the parquet equations [Phys. Rev. B 100, 155149 (2019)]. The latter exploits a diagrammatic decomposition, which classifies the contributions to the full scattering amplitude in terms of their reducibility with respect to cutting one interaction line, naturally distinguishing the processes mediated by the exchange of a single boson in the different channels. We apply this idea to the fRG by splitting the one-loop fRG flow equations for the vertex function into SBE contributions and a residual four-point fermionic vertex. Similarly as in the case of parquet solvers, recasting the fRG algorithm in the SBE representation offers both computational and interpretative advantages: the SBE decomposition not only significantly reduces the numerical effort of treating the high-frequency asymptotics of the flowing vertices, but it also allows for a clear physical identification of the collective degrees of freedom at play. We illustrate the advantages of an SBE formulation of fRG-based schemes, by computing through the merger of dynamical mean-field theory and fRG the susceptibilities and the Yukawa couplings of the two-dimensional Hubbard model from weak to strong coupling, for which we also present an intuitive physical explanation of the results. The SBE formulation of the one-loop flow equations paves a promising route for future multiboson and multiloop extensions of fRG-based algorithms.


I. INTRODUCTION
A major challenge for the theory of many-electron systems is to correctly describe the competing microscopical processes occurring on very different length and time scales. This task becomes particularly hard in the nonperturbative regime of intermediate to strong interactions, where several fascinating phenomena of condensed matter physics take place.
The recently developed [1][2][3] merger of the dynamical mean field theory (DMFT) [4,5] and the functional renormalization group (fRG) [6][7][8][9][10], coined as DMF 2 RG [1], can be regarded as one of the diagrammatic extensions [11] of DMFT designed to be applied in the most challenging parameter regimes of quantum many-body Hamiltonians. By combining the unbiased [9] diagrammatic structure of the fRG to the intrinsic nonperturbative content [12] of DMFT, the DMF 2 RG offers a particularly promising route to tackle crucial nonperturbative features of the many-electron physics. However, the impact of DMF 2 RG calculations will eventually depend on the numerical performance of its actual implementations, where the evident bottleneck is posed by the treatment of the two-particle vertex functions [13,14], due to the large number of momentum and frequency variables needed for their definition.
On a general level, we recall that applying the SBE formalism corresponds to recasting the two-particle diagrams in terms of their reducibility/irreducibility with respect to the cutting of an interaction line. Due to the two-particle nature of the electronic (Coulomb) interaction, the SBE classification shares important qualitative features with the parquet formalism as, for example, the high-frequency asymptotic properties [13,14,26] of the corresponding irreducibile diagrams. At the same time, the SBE classification of diagrams circumvents important problems (such as the multiple divergences of the irreducible vertices [12,[27][28][29][30][31][32][33][34]) which affect parquet-based approaches in the nonperturbative regime.
The specific application of the SBE formalism to the fRG presented here relies on the partial bosonization of the vertex function [35][36][37], similar to the channel decomposition [3,[38][39][40][41][42] already adopted in the context of fRG and parquet solvers [24,43] (for recent developments in this direction see also Refs. [44][45][46][47]). In addition to the screened interaction, a fermion-boson Yukawa coupling [48,49] (or Hedin vertex [50]) is determined from the vertex asymptotics, similarly to the construction of the kernel functions describing the high-frequency asymptotics, see Ref. [14]. Their relation allows to re-arXiv:2105.11749v4 [cond-mat.str-el] 5 Jun 2023 cover the flow equations of the screened interaction and Yukawa coupling in the SBE representation. We note that the obtained structure is apparently the same as the one reported in Ref. [51], where instead of the highfrequency limit the zero frequency value has been used.
Several bosonization procedures have been already developed for the weak coupling fRG, by applying, for instance, the Hubbard-Stratonovich transformation on the bare action [52][53][54][55][56], or by means of the dynamical bosonization [36,57,58]. Instead, the description of the DMF 2 RG vertex in terms of exchanged bosons is a highly nontrivial task. Indeed, the complex frequency structure of the initial DMFT vertex function [13] prevents a straightforward application of the Gaussian integration in the Hubbard-Stratonovic procedure or of the dynamical bosonization. In this perspective, the SBE decomposition offers a relatively simple way to circumvent this problem, expressing the DMFT vertex in terms of bosonic propagators and Yukawa couplings, which can then be used as initial conditions of a mixed boson-fermion flow. In this respect, it is worth stressing that the fermion-boson approach [59,60] has proven to be useful also in the context of diagrammatic extensions of DMFT [61] not based on the fRG formalisms.
The purpose of the present manuscript is twofold. On the one hand, we show how to apply the recently introduced SBE decomposition to both the fRG and the DMF 2 RG treatment of the 2D Hubbard model [62]. On the other hand, we analyze the quality of the approximation resulting from the SBE decomposition on the fRG/DMF 2 RG flows from weak to strong coupling, at half filling as well as away from it. In this respect, we also exploit the transparent nature of the SBE formulation to investigate the physical mechanisms responsible for the enhanced d-wave pairing fluctuations found in the DMF 2 RG calculations.
The paper is organized as follows: In Sec. II we introduce the Hubbard model and present the SBE decomposition together with its implementation in both the fRG and DMF 2 RG flow. In Secs. III and IV we discuss the results for the susceptibilities and Yukawa couplings at half filling and out of it, providing also a comparison with the fermionic formalism. In order to identify the mechanisms responsible for strong d-wave pairing correlations, we perform a diagnostics [28,[63][64][65][66] of the corresponding fluctuations. We finally conclude with a summary and an outlook in Sec. V.

A. The Hubbard model
We consider the single-band Hubbard model in two dimensions, Here we choose the particle-particle (pp) channel as an example. (a) two-particle-and U -reducible diagram in the pp channel. Since this diagram has two interaction vertices at its both ends, it contributes to the screened interaction. (b) Two-particle-and U -reducible diagram, contributing to the Yukawa coupling, as there is only one bare interaction vertex which can be removed to disconnect the diagram. (c) Diagram which is two-particle reducible but U -irreducible, therefore contributing to the SBE rest function.
where c iσ (c † iσ ) annihilates (creates) an electron with spin σ at the lattice site i (n iσ = c † iσ c iσ ), t ij = −t is the hopping between nearest-neighbor sites, t ij = −t the hopping between next-nearest-neighbor sites (the Fourier transform of t ij gives the bare dispersion k ), and U the on-site Coulomb interaction. The filling is fixed by adjusting the chemical potential µ.
In the following we use t ≡ 1 as the energy unit.

B. SBE decomposition
In this section we present the SBE decomposition as introduced in Ref. [15]. All diagrams contributing to the two-particle vertex can be divided into U -reducible and U -irreducible, depending on whether the removal of a bare interaction vertex cuts the diagram into two disconnected parts or not, respectively. Moreover, among the U -reducible diagrams, we can identify three different channels, depending on how the fermionic legs are connected to the removed interaction U . In particular, we find U -particle-particle (U -pp), U -particle-hole (U -ph) and U -particle-hole-crossed (U -ph) reducible diagrams [15]. This classification of diagrams is alternative, and not equivalent, to the more common one, based on the notion of two-particle reducibility (for the relation to the latter we refer to the Appendix A). Indeed, there are some diagrams which are two-particle reducible in a given channel, but U -irreducible. These diagrams are the so-called box diagrams (see Fig. 1). In general, a diagram which is U -reducible in a given channel is also two-particle reducible in the same channel, but not vice versa, the only exception being the diagram composed by a single interaction vertex. Within the SBE decomposition, it is quite natural to interpret the collection of all U -reducible diagrams in a given channel as an effective interaction between two fermions, mediated by the exchange of a boson, representing a collective fluctuation. As mentioned in Refs. [15,24,67] and illustrated explicitly in the following, the SBE decomposition not only significantly reduces the numerical complexity of the vertex functions, but it also allows for a clear physical identi-fication of the collective degrees of freedom arising in a strongly correlated electron system.
In a system with the U(1)-charge and SU(2)-spin symmetries, along with translational invariance, the twoparticle vertex can be expressed as [13]: where σ i represents the spin quantum number, and k i = (k i , ν i ) is a collective variable including the crystal momentum and a fermionic Matsubara frequency. The fourth fermionic variable, k 4 , on which the vertex depends, is fixed by momentum and energy conservation. According to the SBE decomposition, we represent the function V = V ↑↓↑↓ as where k 4 = k 1 + k 2 − k 3 , S, (M + C)/2, and M represent the sum of all U -pp, U -ph, and U -ph reducible diagrams, respectively, the function Λ U irr accounts for all fully Uirreducible diagrams, and a term 2U has been subtracted in order to avoid double counting of the bare interaction, which is already included in the U -reducible channels. The functions S, M, and C, corresponding to the pairing, charge and magnetic channel respectively, depend on two fermionic variables, and a bosonic one, indicated in brackets in Eq. (3). Since each of them is U -reducible in a given channel, their dependencies on the fermionic arguments can be factorized. We can therefore express them as where we name h (h) as left-sided (right-sided) Yukawa coupling and D as screened interaction, which plays the role of an effective bosonic propagator. The right-sided Yukawa couplingsh can be related to their respective left-sided ones, h through the relations By considering the symmetries where Eq. (6a) corresponds to the simultaneous exchange of the ingoing and outgoing variables, and Eq. (6b) to the simultaneous exchange of the two ingoing variables with the two outcoming ones, one can easily prove that h X k+q (−q) = h X k (q), with X = m, c. Therefore the relation holds for all channels. For this reason, from now on we label by h X both the left-and the right-sided Yukawa couplings. It is worthwhile to stress that the equivalence h =h holds because of the choice of notation (3), different choices may lead to to more complicated relations between h andh. The screened interactions D are related to the physical susceptibilities [15,67], that is where χ m , χ c , and χ s are the magnetic, charge, and pairing susceptibilities of the system, respectively. The Yukawa couplings are connected to the so-called threelegged correlators via where the symbol · · · 1PI indicates the connected average with amputated external propagators, and the fermionic bilinears are defined aŝ Here, and from now on, the symbol k = T ν B.Z.
denotes a sum over fermionic Matsubara frequencies and a momentum integration over the Brillouin zone. Notice that in Eq. (9) the division by a term 1 ± U χ is necessary to avoid double counting as it removes from the threelegged correlators all those U -reducible diagrams which are already included in the screened interaction D.
It is interesting to consider the asymptotic high frequency behavior of the Yukawa couplings and bosonic propagators: The limits of large bosonic frequency Ω are quite trivial to prove, as in this case both the susceptibilities and the three-legged correlators are zero. More interesting is the large fermionic frequency limit for the Yukawa coupling. In this limit, one can show by diagrammatic arguments [14] that the three-legged correlator approaches ±U χ (the sign depending on the channel), which makes the Yukawa coupling approaching 1.

C. Flow equations
In this section we derive the one-loop (1 ) fRG flow equations within the SBE decomposition. The dependencies of the various functions on the RG scale Λ are implicit to keep the notation lighter.
We start by focusing on the flow equations within the channel decomposition of Husemann and Salmhofer [39,51], who divided the different contributions to the vertex according to the notion of two-particle reducibility. The flow equation for a given two-particle reducible channel φ reads with X = m, c or s. The bubbles are defined as with G(k) = [(Θ Λ (k)/(iν − k +µ)) −1 −Σ(k)] −1 the propagator, and Σ the self-energy. The symbol ∂ Λ denotes a derivative only on the explicit RG scale dependence of the propagator introduced by the cutoff function Θ Λ (k). The L-functions are given by with V defined in Eqs. (2) and (3). Each two-particle reducible channel consists of two contributions, one which is also U -reducible and one which is not. We can therefore express φ as where we identify the U -irreducible contribution R X as rest function. It obeys the asymptotic relations [14] lim With the help of Eqs. (3) and (11), we further notice that By inserting the limits (11) and (16) into Eq. (12), we hence obtain the flow equations for the screened interactions, Yukawa couplings and rest functions: where . Alternatively, the above equations can be derived through the explicit introduction of three bosonic fields via as many Hubbard-Stratonovich transformations (see Appendix B).
We note that a similar decomposition has been used in Ref. [51], where the flow equations for the screened interactions and Yukawa couplings exhibit the same structure. However, instead of using the vertex asymptotics, they are determined by its value at the lowest Matsubara frequency.
By applying the same line of reasoning, one can generalize the flow Eq. (18), derived within the 1 truncation, to the multiloop extension introduced in Refs. [42,44,68,69]. In the following we discuss its application to the conventional fRG as well as its merger with the DMFT in the DMF 2 RG.

D. fRG
Within the conventional fRG, the fully U -irreducible vertex function Λ Uirr corresponds to the sum of the rest functions of the three channels, and hence does not include any contributions from diagrams that are fully twoparticle irreducible.

Cutoff scheme
In order to run a fRG flow, we regularize the bare Green's function by making use of the so-called Ωcutoff [39,51], that is where the cutoff function is given by

Initial conditions
The fRG initial conditions at Λ = Λ ini read as that, by comparing with Eq. (3), is equivalent to impos- The DMF 2 RG [1, 3] flow differs from the fRG one in its initial conditions, which are defined by the DMFT self-consistent solution of the same problem (see Fig. 2), as well as in the cutoff scheme used. Consistent to previous DMF 2 RG studies [1,3], our DMFT calculations have been performed using the DOS corresponding to the same 2D square lattice dispersion used for the subsequent fRG flow. It should be further noticed that, in order determine the initial conditions within the SBE implementation, one needs to apply the decomposition in Eq. (3) also to the DMFT vertex. The local bosonic propagators and Yukawa couplings can be computed directly from the Anderson impurity model (AIM), from Eqs. (8)- (9), while the fully U -irreducible vertex can be extracted by subtraction from Eq. (3).

Cutoff scheme
When introducing a scale dependence on the bare fermionic propagator through a cutoff, one has to keep in mind that the regularized bare propagator has to smoothly interpolate between the bare propagator of the self-consistent AIM and the bare lattice one, specifically: where ∆ hyb (ν) is the hybridization function of the AIM. Furthermore, we require the DMFT solution for the local propagator to be conserved at each step of the flow, that is [3] B.Z.
the full Green's function of the AIM, and Σ dmft the DMFT self-energy. We therefore make the following choice for the regularized Green's function where we choose the function Θ Λ (ν) to be a smooth frequency cutoff, namely the same as in Eq. (20), which imposes the boundary values for the RG scale Λ ini = +∞, and Λ fin = 0. The function Ξ Λ (ν) is determined by the DMFT self-consistency relation (23). With the above definitions it is straightforward to check that the boundary conditions (22) are consistently fulfilled. This choice for the DMF 2 RG cutoff (24) has a direct intuitive physical meaning: at a given scale Λ, indeed, the fermionic modes with energies |ν| Λ are nonlocal and their contributions on top of the DMFT solution have been included, while the ones with energies |ν| Λ still belong to the AIM and do not yet contribute to the generation of nonlocal correlations.

Initial conditions
The DMF 2 RG employs the DMFT solution as a correlated starting point for the RG flow. The initial conditions for the screened interactions and Yukawa couplings therefore read where D loc and h loc are the bosonic propagator and Yukawa couplings of the self-consistent AIM. Concerning the rest functions, we set their initial value to zero, so that they will represent only the nonlocal contributions, the local ones being already included in the U -irreducible vertex function of the AIM, Λ loc U irr . The flow equations for the screened interactions, Yukawa couplings, and rest functions are the same as in the plain fRG.

F. Form factor decomposition
In order to simplify the set of flow equation described above, throughout this work we project the Yukawa coupling dependence onto the secondary spatial momentum k onto s-wave form factors, f s k ≡ 1, that is we approximate with k a shorthand for B.Z.
We finally note that if the residual vertex develops a strong momentum dependence, considering a limited number of form factors is not justified a priori and has to be verified.

G. d-wave pairing channel
In the doped regime we include d-wave pairing fluctuations into our parametrization of the two-particle vertex by where M νν (q), C νν (q), and S νν (q) are the the s-wave projections of the couplings of Eq. (4), and, in order to deal with a positive quantity, we have chosen D with a minus sign in front of it. The function Λ loc U irr (ν 1 , ν 2 , ν 3 ) represents the sum of U -irreducible diagrams at the DMFT level and therefore does not flow, that is in the doped regime we neglect the flow of rest functions in the magnetic, charge, and s-wave pairing channels. The d-wave pairing channel D kk (q) is given by with the d-wave form factor f d k = cos k x − cos k y . In essence, the function D νν (q), for transfer momentum q = 0, represents the sum of all diagrams that are reducible in the two-particle-pp channel and at the same time exhibit a d-wave symmetry in the dependency on secondary momenta k and k . Due to the locality of the bare Hubbard interaction U , all the above mentioned diagrams are U -irreducible, that is the d-wave pairing channel will consist exclusively of its rest function. Its flow equation reads and where k = (k, ν), and k = (k , ν ). This corresponds to restrict ourselves to the pure d-wave first harmonic contribution. The flow equations for the Yukawa couplings and bosonic propagators of the magnetic, charge, and s-wave pairing channels are the same as in the previous section, with the difference that the feedback of the d-wave pairing channel on them is included through Eq. (29).

III. RESULTS AT HALF FILLING
To illustrate the potential of our SBE implementation of the DMF 2 RG scheme (see Appendix C for the details on the numerical implementation), we first consider the testbed case of half filling. The expected physical behavior, encoded in the response functions of the different channels, is quite clear in this situation, due to the underlying particle-hole symmetry of the problem.
In the spirit of the SBE, for most of the calculations presented in this section we have neglected the flow of the nonlocal corrections to the rest functions (R X = 0). Such an SBE-based simplification allows for a substantial gain in memory cost, compared to the conventional formalism, as the Yukawa couplings and the screened interactions depend on less arguments than the full vertex.
This turns also into a gain in computational time, which is however not trivial to estimate. At the same time, it is important to underline that, within the DMF 2 RG framework, this approximation corresponds to neglect only the flow describing the nonlocal corrections to the rest function, as its purely local part is per construction included nonperturbatively via the DMFT initial condition.
The quality of our SBE-based approximation to the DMF 2 RG flow has been eventually tested by performing specific calculations with and without the inclusion of the flow for the nonlocal corrections to R X and by comparing the obtained results. The quite positive outcome of this test has been explicitly discussed in Sec. III B.
A. Susceptibilities

Weak-coupling fRG and DMF 2 RG results
We start by comparing different methods in the weakcoupling regime, where also conventional perturbative approximations can be applied as a benchmark. In particular, for the calculations of physical susceptibilities presented in Fig. 3, we have exploited: i) the random phase approximation (RPA); ii) the one-loop fRG, where we have also neglected all the R X functions [70]; iii) the multiloop extension of the fRG [42], which amounts to the resummation of all diagrams of the parquet approximation (PA), thus fulfilling the Mermin-Wagner theorem in 2D [18]; iv) the DMFT; and v) the DMF 2 RG, where the nonlocal R X functions are neglected. For the DMF 2 RG, similarly as for the fRG, neglecting the nonlocal rest functions does not produce sizable deviations in this first parameter regime considered (see App. A). In order to allow for a quantitative comparison of all the approaches considered on an equal footing, these specific weak-coupling fRG and DMF 2 RG calculations have been performed including the feedback of the self-energy in the RG flow.
In the upper panels of Fig. 3, we show the screened interactions D X (q, Ω) for zero transfer frequency Ω = 0 and a specific momentum transfer. The lower panels show the corresponding susceptibilities that are connected with the screened interactions via Eq. (8). As expected, for the lowest values of U , the different calculations display a very good agreement, which is gradually reduced for larger values of the coupling. Consistent to the physics expected in the half-filled Hubbard model, all methods outline predominant antiferromagnetic (AF) correlations for the coupling values considered, though they differ on the quantitative level. In particular, the RPA yields, for all values of U , the largest magnetic effects. These get suppressed, to different amounts, in the results of the other approaches, because they all capture, differently from RPA, the competition with the other channels. In particular, a sizable suppression of the AF susceptibility for U > t is clearly observed, within fRG, already at the level of the 1 calculations. Not surprisingly, the suppression due to the channel interference becomes stronger in the PA results, reflecting the inclusion of the corresponding fluctuation effects at all loop orders. Due to the fully two-particle irreducible diagrams which are not included in the PA, it appears to systematically underestimate χ m as compared to numerically exact determinant quantum Monte Carlo data [44,71] for larger values of U . At the given temperature T = 0.2t, suppression effects are observed also in the DMFT data. The reduction of the magnetic fluctuations with respect to the 1 fRG is ascribable to self-energy effects, which reduce more quickly the electronic coherence in DMFT than in fRG. This is consistent with the observation that DMF 2 RG results closely resemble the ones of DMFT in this parameter regime. Indeed, the nonlocal correlations included at the 1 level on top of the DMFT results, induce, as expected, a further reduction of the magnetic correlations, but their quantitative impact remains marginal for U < 3t. We note that nonlocal multiloop corrections, which recover the PA result, have a stronger impact on the suppression than the local DMFT ones at 1 level.
We turn now to discuss the other sectors, which are secondary in the physics of the half-filled model. The corresponding results are reported in the central and right panels of Fig. 3, showing that the correlation effects beyond RPA are overall weaker and significantly less dependent by the particular approach chosen than in the magnetic sector. In particular, the slight suppression of the uniform charge and s-wave static susceptibilities induced by the interplay with their complementary channels get reflected in a slight increase of corresponding screened interaction with respect to the RPA values.
We note that including the flow of the fermionic rest functions R X leads only to negligible corrections of the results shown in Fig. 3, see Appendix A. In the weakcoupling regime, the screened interactions and susceptibilities obtained without R X correctly describe the physical behavior, justifying the application of the SBE approximation.

Intermediate-to strong-coupling DMF 2 RG results
After this preliminary comparison, we move towards the more challenging parameter regime of intermediate to strong couplings. Hence, from now on, we will focus on DMF 2 RG results for the physical susceptibilities of the different sectors, whereas, for simplicity, we turn off the flow of the nonlocal self-energy corrections. Especially in the strong-coupling regime, this simplification is not expected [3,11,72,73] to affect the final results for the response functions in a significant way. We start our analysis of the strong-coupling regime, by presenting in Fig. 4 our DMF 2 RG results for the whole momentum dependence of the static susceptibilities (at zero bosonic frequency) in the magnetic and in the charge sectors for the highest coupling considered U = 16t, which is significantly beyond the critical interaction value of the Mott-Hubbard metal to insulator transition of DMFT (U MIT (T = 0) ∼ 12t), and for a temperature slightly above the DMFT Néel temperature (see leftmost panel of Fig. 8 for the precise location in the DMFT phasediagram). Note that at half filling particle-hole symmetry implies χ s (q, ω) = χ c (q + Q) for the pairing susceptibility, with Q = (π, π). The magnetic susceptibility exhibits a very pronounced peak around momentum q = (π, π), a hint of strong AF fluctuations. This indicates that, similarly as for weak-coupling data, the nonlocal correlations captured by the 1 DMF 2 RG do not suppress the ordering temperature of DMFT in a significant way. At the same time, the static response in the charge (and pairing) sector bears very clear hallmarks of the strong-coupling physics. Except for a residual momentum modulation, these response functions appear massively suppressed (note the different order of magnitude of the scales in the two panels), reflecting the almost insulating nature of the Mott-Hubbard physics at finite T .
Following the analysis of Ref. [12], it is instructive to relate the results presented so far to the behavior of the corresponding generalized susceptibilities, which describe the underlying fermionic scattering processes. We recall that, in general [74], the susceptibilities χ X (q) describing the physical response of the system can be obtained by the generalized ones χ X νν (q) by summing over all the fermionic Matsubara frequencies ν, ν . In our notation, the explicit definition reads, also referred to as "postprocessing": where From a numerical side, we note that a larg number of Matsubara frequencies are required for the internal summation in order to account correctly for the highfrequency asymptotics. A restricted frequency summation yields unphysical negative values in the charge response function due to the formation of the local moment [12]. For illustrative reasons, we here explicitly discuss the Ω = 0 case of the generalized charge and magnetic susceptibilities for two rather different interaction values, namely U = 4t and U = 16t, which correspond, in DMFT, to a Fermi-liquid paramagnetic metallic (PM) and Mott-Hubbard paramagnetic insulating (PI) regime, respectively. The corresponding data are shown in Figs. 5 and 6, where we report the on-site [75] generalized magnetic (left panels) and charge (right panels) DMFT solutions used as input of our DMF 2 RG calculations in the upper panels and the corresponding DMF 2 RG results for the uniform (q = 0) generalized susceptibility in the lower ones. As discussed in Ref. [12], the frequency dependence of generalized charge susceptibilities represents a sensitive compass for identifying the underlying physics in fundamental many electron models, since it directly describes the impact of correlations on the electronic mobility and unveils its link to the changes of magnetic response. A first glance to the data of Figs. 5-6 shows that the qualitative modifications in the frequency structures of the generalized susceptibilities occurs indeed in the charge sector. In particular, sign changes in relevant frequency structures of the generalized charge susceptibility take place from weak (U = 4t) to strong coupling (U = 16t), reflecting important differences of the dominating physical mechanisms at play in the two regimes.
For U = 4t shown in Fig. 5, the leading frequency dependence in the generalized charge susceptibilities of both the AIM (DMFT) and DMF 2 RG results appears in the main diagonal structure. This assumes positive values, larger at low frequencies ν = ν and slowly decaying for larger ν = ν values. This feature is a typical hallmark of the metallic physics in the perturbative regime, as it arises from the bubble term contribution χ 0 νν = −βΠ c νν (Ω = 0) = −βG(ν)G(ν )δ νν > 0, built upon a metallic G(ν). The role of vertex corrections, for U = 4t, appears merely quantitative, yielding an overall slight suppression (enhancement) of all entries of χ c νν (χ m νν ), responsible for the emergence of small negative (positive) off-diagonal contributions (faint bluish/reddish color for small ν = ν in the right/left upper panels of Fig. 5). The moderate size of the vertex corrections is highlighted by the comparison of the two channels, whose results are both dominated by a positive diagonal frequency structure. Due to the proximity to an AF instability, the inclusion of nonlocal correlations in DMF 2 RG strongly affects the generalized susceptibility of the predominant magnetic channel, which gets significantly enhanced with respect to the corresponding AIM data. At the same time, the impact of nonlocal correlations appears marginal in the charge sector.
The data computed at U = 16t (Fig. 6) display relevant differences, which are induced by significant vertex correction effects. In the magnetic channels these drive an overall enhancement of the generalized susceptibility (diffuse reddish zone clearly visible in the left panels), which gradually overcomes the residual diagonal contribution of the bubble term. Such generalized enhancement of χ m νν encodes the typical Curie behavior of the local magnetic response in this regime. At the same time, the corresponding vertex corrections in the charge sector strongly suppress the physical response, by flipping the sign of the diagonal entries of χ c νν up to an energy scale of order U . Specifically, the largely negative low-frequency diagonal structure visible in the right panels of Fig. 6 is directly responsible [12,28] for freezing the local density fluctuations in the Mott PI regime [76]. Physically, it can be interpreted [12] as the charge counterpart of the local moment formation in the magnetic sector.
We note that the strong differentiation of the charge and magnetic response, which is a typical hallmark the Mott PM phase, has a clear nonperturbative origin: The negative diagonal structure emerging at low-frequency is associated to negative eigenvalues of χ c νν . This implies that sign flips with respect to the (positive) eigenvalues of the weak-coupling regime must occur by in-creasing U . By any vanishing eigenvalue, however, the matrix χ c νν becomes singular leading to multiple divergences [27,28,30,[32][33][34]77] of the irreducible vertex function and to the corresponding [31] problem of the multivaluedness [78][79][80][81][82] of the Luttinger-Ward functionals, which have been extensively discussed in the recent literature. As these singularities cannot be captured by any perturbative approach (including truncated fRG and PA) [12], for a proper description of the intermediate to strong-coupling regime it is pivotal to include the associated nonperturbative physics (such as local moment formation and its Kondo screening [12]) through a DMFT starting point, emphasizing the necessity of a DMF 2 RG treatment. In this specific case, the DMF 2 RG data (lower panels of Fig. 6) display a further increase of the generalized magnetic susceptibility at q =(π, π), due to the proximity to the AF instability in the phase diagram, as well as a further suppression of the generalized charge susceptibility at q = 0 with respect to the AIM results, consistent with the incompressible nature of the Mott insulating ground state.

B. Comparison with conventional implementations
All calculations presented above have been performed neglecting the flow of the fermionic rest functions R X , consistently within the SBE framework. While this approximation allows for a significant reduction of the numerical effort, it is important to verify to what extent its application is justified in the different coupling regimes.
We briefly recall here that, when fully considering the flow of the rest functions R X , the SBE-based implementation becomes formally equivalent to the one used in the conventional fRG based on the 1PI vertex function, see Ref. [3] for the specific case of DMF 2 RG [1]. We first analyze the frequency dependence of the (previously neglected) R X functions, evaluated for the same bosonic variables (q, Ω) of the susceptibilities shown in the preceding subsection, at U = 4t and U = 16t. The corresponding results are reported in Fig. 7. This highlights a key feature of the rest function, namely its characteristic decay to zero at large frequencies. Qualitatively, this reflects a general feature of the SBE diagrammatics [15]: the high-frequency asymptotics of the two-particle vertex functions is fully captured by the effective interactions, D X , and by the Yukawa couplings, h X . From a more quantitative perspective, the high-frequency decay appears particularly pronounced at large couplings: In contrast to U = 4t, for U = 16t the rest function exhibits a significant frequency dependence only for the lowest Matsubara frequencies, which, at this temperature, represents the leading frequency dependence of the whole 1PI vertex function. Despite the large values of the rest functions at strong coupling, the single-boson contributions h X D X h X are of the same order or even larger over an infinitely broad frequency range. Furthermore, we checked that the R X have a marginal feedback effect on the flow of the Yukawa couplings due to cancellations with the fermionic bubbles in the insulating phase.
After examining the high-frequency decay properties of R X , which highlight the overall convenience of an SBE-based formalism, we turn now to the analysis of the specific impact that neglecting them has on the final results of our DMF 2 RG calculations. In Fig. 8, we report the susceptibilities of the predominant magnetic channels, computed with and without the inclusion of R X , for three different values of U depicted on the phase diagram in the inset. The parameter set chosen is particularly relevant due to its proximity to the AF instability of the DMFT calculations, and the expected relevant contribution of large magnetic fluctuations. We note immediately that at half filling the magnetic susceptibilities obtained without R X correctly describe the physical properties expected in the whole parameter region considered, including the strong-coupling Mott-insulating regime at U = 16t. Furthermore, the corrections to the magnetic response induced by the inclusion of R X in our DMF 2 RG calculations appear to be overall rather marginal (see insets of Fig. 8) and further decreasing at larger interaction values. As a consequence of this, the inclusion of the rest function has also a minor impact on the determination of the temperature of the AF instability (Néel temperature T N ), which can be finite in all approaches where the Mermin-Wagner theorem is violated. In Fig. 9 the inverse magnetic susceptibility is shown as a function of the temperature, with and without the inclusion of R X . Both data-sets clearly display a linear mean-field like critical behavior, resulting in T N = 0.4042t and T N = 0.3986t respectively. These values, both slightly lower than the corresponding DMFT results, are quite close.
Hence, the overall effect of neglecting R X appears to be marginal, even quantitatively, within the 1 DMF 2 RG scheme, in particular for the dominant magnetic channel. Not surprisingly, the corrections induced by R X are even lower in the other channels: charge and pairing susceptibilities determined with and without the rest function display almost no difference (not shown).

C. Yukawa couplings
In this subsection, we focus on the second main ingredient of the SBE formalism, that is the (fermion-boson) Yukawa couplings, s. also Refs. [47,67,83,84]. The corresponding results for the magnetic and the charge sector are shown in Fig. 10, while the s-wave pairing one can be obtained from the relation h s ν (q, Ω) = h c ν (q + Q, Ω), with Q = (π, π), valid when the particle-hole symmetry is realized.
At weak coupling (U = 4t) we observe for both channels a moderate dependence of the static (Ω = 0) Yukawa couplings on the fermionic frequency ν. This results in an overall slight suppression of the asymptotic value of the coupling (h X ∼ 1) and encodes the metallic screening effects [85,86] occurring in the proximity of the Fermi level. The inclusion of nonlocal effects in DMF 2 RG induces mild corrections with respect to the DMFT (AIM) results (black dashed line). In particular, it can be noted that, in both sectors, the fluctuations associated to a transfer momentum of q = (π, π) tend to magnify the frequency dependence displayed by the local DMFT (AIM) calculations, while the uniform ones q = (0, 0) have the At strong-coupling, the situation appears essentially reversed. The DMFT (AIM) calculations performed at U = 16t show a much stronger dependence of the couplings on the fermionic Matsubara frequency and a dycothomic behavior in the two channels -a typical hallmark [12,87] of large vertex corrections. In particular, the magnetic Yukawa coupling gets significantly enhanced at low-frequencies with respect to the h m = 1 value, whereas the largest value is found in DMF 2 RG for the AF ordering wavevector. The charge channel, instead, exhibits a strong suppression, that even makes h c negative for the smallest frequencies. Here, the strongest suppression, according to our DMF 2 RG results occurs for q = (0, 0), consistent to a vanishing isothermal compressibility in the Mott-Hubbard ground state [5,87,88]. Physically, these opposite behaviors should be regarded as the nonperturbative fingerprints [12] of the local moment in the Yukawa couplings, since its formation enhances the system's response to a magnetic perturbation and, at the same time, freezes the fluctuations of the electronic density.

IV. RESULTS AT FINITE DOPING
To showcase the validity of our method in a physically more relevant parameter regime, we have also performed DMF 2 RG calculations for the non particle-hole symmetric case. In particular, we consider an intermediatecoupling of U = 8t and a fairly low temperature T = 0.044t, doping the system to the filling n = 0.82, and frustrating the magnetic correlations with a next to nearest neighbor hopping t = −0.2t. As the Mermin-Wagner theorem is not fulfilled at the level of the 1 truncation of the DMF 2 RG [1, 3, 42, 68], we observe a magnetic instability in the considered parameter region. Hence, in order to highlight the relation between superconductivity and magnetism, we stop the flow before reaching the final scale, as customarily done in the fRG treatment of pseudocritical transitions. In particular, the divergence of the magnetic propagator has been prevented, by stopping the flow when the magnetic propagator D m (q) exceeds a value of 8 × 10 3 t, which, for U = 8t, corresponds to a susceptibility of ∼ 120/t, resulting in a stopping scale Λ cr 0.067t. Although the criterion for stopping the flow is arbitrary, close to an instability the screened interactions diverge quickly, implying that even large changes in the stopping criterion translate into tiny differences in the stopping scales. This makes the results weakly dependent on the specific choice of the condition to interrupt the flow.
Eventually, as the PA fulfills the Mermin-Wagner theorem in 2D [18], we expect the future inclusion of multiloop corrections to suppress the magnetic screened interaction, allowing to continue the flow down to lower scales (to Λ = 0 at loop convergence).
To account for d-wave pairing correlations, we have included the flow of the d-wave pairing channel (consisting exclusively of its rest function), as explained in Sec. II G. We neglect the rest functions of the magnetic, charge and s-wave pairing channels since these appear to yield only minor corrections, together with nonlocal corrections to the self-energy.

A. Susceptibilities
Here, we will first discuss the results obtained for the magnetic χ m , charge χ c , and d-wave pairing χ d susceptibilities computed at the stopping scale Λ cr . While the first two can be directly extracted from the respective propagators, see Eqs. (8a) and (8b), the d-wave pairing one is computed through where we have neglected terms of the type T 2 νν Π sd ν (q)L s νν (q)Π sd ν (q), with Π sd ν (q) the mixed sd-wave pairing bubble, as it is nonzero only for q = 0 and generally rather small [89]. In Fig. 11, we show the static susceptibilities in the whole Brillouin zone. Specifically, the magnetic susceptibility is close to a divergence at momenta (π, π ± 2πη) (as well as, by symmetry, at (π ±2πη, π)) with η 0.08, which indicates the tendency towards an incommensurate magnetic instability [9,[90][91][92]. We recall that magnetic long-range orders with incommensurate wave vectors have been found in several mean-field studies [93][94][95][96] of the Hubbard model at finite doping. Similar conclusions have been also drawn when including fluctuations beyond the static mean-field, for example, by means of expansions in the hole-density [97][98][99][100], or exploiting extensions of DMFT [101,102]. DMFT calculations suggest that the ordering wave vector is related to the Fermi surface geometry not only at weak, but also at strong coupling [103].
Differently, the charge susceptibility exhibits a rather weak dependence on q, that is, only moderate deviations are observed from the local description of DMFT in this sector. We expect this feature not to depend on the fact that we have a small finite Λ cr . At a closer inspection, anyway, one notes two peaks located at (π, 0) and (0, π). These signal the presence of mild charge-stripes correlations.
Eventually, we focus on the d-wave pairing susceptibility. Although its absolute values are not excessively large, it presents a pronounced q-dependence with a welldefined peak at q = 0. Hence, one can reasonably expect that, if we were able to continue the flow below Λ cr , χ d would further increase, possibly even diverging at some finite scale. This heuristic expectation is supported by the analysis presented in the following sections.

B. Yukawa couplings
In this section, we analyze the real part of the magnetic, charge, and s-wave pairing Yukawa couplings at the stopping scale Λ cr shown in Fig. 12. The magnetic Yukawa coupling exhibits a behavior that qualitatively resembles the one in the weak-coupling regime of the halffilled case (compare with upper left panel of Fig. 10) despite the relatively large value of U = 8t. This happens because the local moment and, thus, its fingerprints in the Yukawa coupling, gets weakened by hole-doping: its low-frequency part is suppressed with respect to the highfrequency (h m = 1) value, due to the electronic screening.
Finally, the s-wave pairing Yukawa coupling appears overall suppressed, in spite a relatively small upturn at the lowest frequencies. Not surprisingly, it displays a rather marginal momentum dependence and, thus, small deviations from the AIM result.
C. d-wave pairing correlations

Diagnostics of the correlations
The aim of this section is to thoroughly inspect all terms contributing to the sizable enhancement of the susceptibility χ d , which we briefly discussed above. In the spirit of the post-processing procedures [28,[63][64][65] recently applied [15,[64][65][66]104] to several quantum manybody approaches to the Hubbard model, we decompose Eqs. (33) and (36), in order to distinguish the different contributions to the d-wave susceptibility. In particular, by combining Eq. (29) with (33), we get where we have introduced the functions with X = M or C, depending only on frequencies. We can therefore split the function L d in three distinct contributions: and D νν (q). By inserting this decomposition into the expression for the susceptibility (36), we obtain Here, χ d(0) is the bare bubble term, while χ d(m) and χ d(c) represent the Maki-Thompson contributions to the susceptibility, as diagrammatically shown in Fig. 13. Finally, since the d-wave pairing channel entails all diagrams which are two-particle (pp) reducible but Uirreducible, the remaining χ d(N b ) term can be interpreted, as the sum of N -boson processes. The physics encoded in the latter processes will be analyzed separately in the final part of this section. Here, we only Figure 12. Yukawa couplings in the magnetic, charge, and s-wave pairing channel (from left to right) at zero bosonic frequency, as a function of the fermionic Matsubara frequency ν, for the same parameters as in Fig. 11 and for various choices of the spatial momentum q.
remark that this representation becomes exact only in the multiloop extension. Before commenting on the results, we note that, from an fRG perspective, the AF susceptibility gradually evolves from the beginning of the flow [3,39,86,90], while the superconducting d-wave susceptibility emerges only in proximity of the critical scale [91,105].
In Fig. 14, we plot the different contributions to the d-wave static susceptibility defined above, as obtained at the stopping scale Λ cr . We observe that in the considered parameter region the magnetic Maki-Thompson processes yield the most relevant contribution. Moreover, we find a fairly large contribution of the multiboson terms. Indeed, we expect that, by approaching a d-wave pairing instability, the multiboson term would dominate over the Maki-Thompson ones, eventually driving the divergence of the associated susceptibility at the thermodynamic instability. More in general, if one were able to separate the N <N (χ d (N b <N ) ) from the N ≥N (χ d(N b ≥N ) ) bo- son processes in the susceptibility, arbitrarily close to the critical point one would detect for every finiteN . This happens because the divergence of the susceptibility is due to a term (1 − λ d ) −1 , with λ d , the maximum eigenvalue of the matrix product between the d-wave bubble and the two-particle irreducible vertex in the d-wave pairing channel, approaching 1. Since the term (1 − λ d ) −1 results from a resummation of infinite order diagrams, it can only be encoded in the N ≥N term, with the N <N one scaling asN λN d ∼N close to the instability. In general, a measure of the maximum N for which Eq. (42) is fulfilled, might be exploited to quantify the actual proximity to a thermodynamic (dwave pairing) instability.
As anticipated above, we eventually draw our attention on the last term in Eq. (40) identified by the presence of D νν (q), which entails all the reducible scattering processes in the d-wave pairing channel. Due to the d-wave symmetry in momentum space, D νν (q) arises exclusively from scattering events involving the exchange of two or more bosons. A refined inspection of these scattering processes is then possible via a corresponding decomposition of D νν (q).
According to Eqs. (31) and (39), we can classify the different contributions to D νν (q) directly from the corresponding flow equation: The diagrammatic representation of the terms D mm , D cc , and D mc is shown in Fig. 15. These two-boson processes can be associated to the so-called Aslamazov-Larkin diagrams. As one might expect, a closer inspection of Eqs. (43a)-(43c) shows that they are not fully reconstructed during the flow because the functions L d(m) , and L d(c) (as well as the nonlocal self-energy if included) also depend on the fRG scale. This feature is a typical artifact of the 1 truncation, which can be fully resolved  in the framework of future multiloop extensions of the approach. With this caveat, it is nonetheless possible to interpret D mm , D cc , and D mc as the two-boson contributions to the d-wave pairing channel, and the remainder D N b ≥3 in terms of higher order processes in the number of exchanged bosons. In Fig. 16 we plot the various terms contributing to D νν (q) as a function of the fRG scale Λ. We observe that the multiboson (N b ≥ 3) term develops at significantly lower scales as compared to the twoboson ones. At the same time, it increases considerably by approaching the stopping scale. This behavior is consistent to the following general consideration: If the system is close enough (in temperature, doping, or other parameters) to a thermodynamic instability, at some point D N b ≥3 will overtake the other terms in Eq. (40), and eventually diverge at the transition itself. In fact, similar as discussed for the susceptibility, sufficiently close to the transition the N ≥N term is going to exceed the sum of the N <N terms for every finiteN . Consistent with these general consideration, in the framework of our 1 DMF 2 RG, we indeed observe that the multiboson term D N b ≥3 becomes dominating at a finite scale just before the end of the flow. Hence, for the selected parameter choice, an important precondition for the onset of a dwave pairing instability has been already realized. This represents a promising hint that a true superconducting transition may be unveiled at lower temperatures by means of higher loop-order calculations.
In Fig. 17 we show the frequency structure of the twoboson contributions to the d-wave pairing channel as well as D itself, at q = (0, 0) as functions of the fermionic Matsubara frequencies ν, ν . Being associated only to U -irreducible diagrams, the d-wave pairing channel is a rapidly decaying function of the Matsubara frequencies. It exhibits a structure centered around the fre- quencies ν = ±ν 0 , ν = ±ν 0 (ν 0 = πT ), where it assumes fairly large values. About 50% of this structure is generated by two magnetic boson processes and (most of) the rest by multiboson ones, as the D cc and D mc terms play a very marginal role in the formation of dwave pairing correlations. It appears hence natural to conclude that among all multiboson terms, the one consisting only of multiple magnetic boson processes will have the largest weight. This observation suggests that d-wave pairing correlations in the normal phase of the Hubbard model are mostly generated by (incommensurate or not) AF fluctuations, consistent with previous findings in fRG [3,39,86,91] and DMF 2 RG [1,3], as well as in recent numerical analyses [63,65,106].

Remarks on the bosonization of the d-wave pairing channel
In view of further reductions of the numerical complexity in future applications, it might be helpful to describe also the d-wave pairing channel D νν (q) in terms of single boson processes. Since the bare interaction U has s-wave symmetry, the diagrammatic argument of the SBE decomposition does not hold in this case. Hence, a proper decomposition of the d-wave pairing channel that factorizes the dependence on the fermionic frequencies is needed.
Inspired by earlier fRG works where the Yukawa cou-pling has been calculated from the channel functions at given values of the fermionic frequencies [51], we can define the d-wave screened interaction and boson-fermion coupling as whereν 0 is a fixed fermionic frequency (eventually qdependent). We notice that in Eq. (44) the symmetrization over ±ν 0 is necessary to guarantee the correct symmetries for D d and h d . In this way the d-wave pairing channel can be expressed in a similar form as the other channels: A key point to keep in mind when choosingν 0 is that, once a d-wave pairing instability occurs, the divergence of the effective interaction must be reabsorbed into the bosonic propagator, while the Yukawa coupling and the rest function should remain finite, similarly to what happens in the other competing channels.

V. CONCLUSIONS
We have applied the recently introduced SBE representation to the fRG and DMF 2 RG, which relies on a diagrammatic decomposition in contributions mediated by the exchange of a single boson in the different channels. Specifically, the (1 ) flow equations for the two-particle vertex are recast into SBE contributions and a residual four-point fermion vertex.
The SBE-based formulation leads to a substantial reduction of the numerical effort, since the corresponding rest function is significantly localized in frequency space, especially in the strong-coupling regime. This justifies the approximation to significantly restrict the total number of frequencies taken into account in the RG flow or even to fully neglect the nonlocal corrections to the DMFT rest function. The reduced numerical effort facilitates the applicability of the fRG and DMF 2 RG to the most interesting regime of intermediate to strong correlations and/or low T . The advantage of this implementation is well illustrated by hands of our DMF 2 RG calculations of the 2D Hubbard model performed up to very strong interaction (U = 16t) at and out of half-filling. In this case, we specifically analyze the impact of including/excluding the rest function flow and observed a marginal effect from weak to strong coupling. Moreover, the SBE decomposition naturally allows for a clear physical identification of the relevant degrees of freedoms. As pertinent example, we exploited this specific feature of our approach to diagnose the tendency towards a d-wave superconducting instability of the doped Hubbard model in terms of magnetic and charge driven processes.
The derivation of the SBE-based RG flow within the 1 truncation represents the natural starting point for future multiloop extensions [42,44,68,69], as well as for the systematic inclusion of multiboson contributions. Within these methodological extensions, fRG-and DMF 2 RGbased computation schemes can be brought to a quantitative level for all coupling strengths. Physically, we expect that nonlocal correlation effects, associated to higher loop order processes, might become progressively more pronounced, especially in the most challenging lowtemperature regime.
Finally, the present formalism offers the possibility to explicitly introduce bosonic fields and study the flow of a mixed boson-fermion system. This would allow, for example, to study the effect of bosonic fluctuations on top of mean-field solutions [92,[107][108][109] below the (pseudo) critical scale, where symmetry breaking occurs.

ACKNOWLEDGMENTS
The authors thank P. Chalupa Figure 18. Comparison between the SBE and the conventional fermionic formalism, both for the fRG (upper panels) and DMF 2 RG (lower panels), for the same parameters as in Fig. 3. Main panels: magnetic screened interactions, computed with and without the inclusion of the rest functions. Insets: difference between the two above mentioned screened interactions.
Yukawa coupling is related to the K (2)X asymptotic function by with At this stage the SBE decomposition seems to offer no substantial computational gain as compared to the channel asymptotics, the only exception being in the vicinity of a critical point, where in the asymptotic formalism both K (1)X and K (2)X acquire large values, while in the SBE one this occurs only for D X , the Yukawa coupling being always finite. However, the most important difference is visible in the rest functions. Indeed, the SBE and the asymptotic rest functions are related via (see also Ref. [109]) (A6) From Eq. (A5) we notice two important facts: First, in the asymptotics formalism some rest functions can diverge at a given critical point as they contain D X . Therefore, in this regime it is not safe to neglect them. Second, R SBE,X contains less diagrams than R asymX , as the latter still includes some U reducible contribution. An approximation which neglects R SBE,X can be always justified as the selection of a well-defined class of diagrams.
For completeness we report the fRG and DMF 2 RG results of Fig. 3 obtained by including the flow of the rest functions R X , where the SBE-based implementation is equivalent to the fRG and DMF 2 RG respectively. In Fig. 18 we show the impact of the inclusion/neglection of the rest function in the magnetic screened interaction, which displays the largest difference, for both conventional fRG and DMF 2 RG. We notice that, as discussed in the main text for stronger coupling, sizable differences between the two approaches (with and without R X ) arise only in the vicinity of the Néel transition. We then conclude that, away from the AF transition point, the tiny differences fully justify the approximation to neglect the rest function, while close to the critical transition the inclusion/neglection of the rest functions will result in a (generally small) change of the values for the critical temperature or coupling.
(B5) Its value at some initial scale Λ ini depends on the formalism used. For instance, in the plain fRG we impose R Λ→Λini (k) → ∞, such that, from the saddle point equation of the functional integral, the initial conditions are determines by the bare ones.
Differently, in the DMF 2 RG, the cutoff must fulfill so that the rescaled action at Λ ini is where S AIM ψ,ψ is the action of the self-consistent AIM. Here the same procedure involving HST for the Hubbard interaction at the impurity has been performed. In this way, the initial condition for the resulting effective action reads Γ Λini ψ,ψ, Φ = Γ AIM ψ,ψ, Φ , where Γ AIM ψ,ψ, Φ is the effective action of the selfconsistent AIM. By expanding it in terms of 1PI functions, one recovers the initial conditions given in the main text, where the screened interactions D X and the Yukawa couplings h X play the role of bosonic propagators and fermion-boson interactions respectively.
Note that we do not explicitly introduce a regulator for bosonic fluctuations. Indeed, the bosonic path integral at this stage has to be performed in some way to recover the fermionic formalism. For this scope, a second scale Λ b can be introduced, responsible for the integration over Φ and which is associated to a second cutoff function R Λ b b regularizing the bare bosonic propagator. In this context, the bosonic flow is shown to be ineffective when integrated before the fermionic flow integration [54].

Appendix C: Numerical implementation
In this Appendix we discuss a few numerical details. Regarding the DMFT calculation, we solve the AIM with the exact diagonalization (ED) by discretizing the bath into 4 sites. The local magnetic, charge and pairing screened interactions D X and Yukawa couplings h X are computed by Lehmann representation as in Ref. [42] in a relatively large frequency range as well as the two-particle Green's function needed to extract the U -irreducible vertex Λ loc U irr . Regarding the flow equations of D X and h X , we use a frequency domain ranging from 30 to 64 positive values for both the bosonic and the fermionic Matsubara frequencies, depending on the convergence of the calcu-lation. The treatment of the frequency asymptotics here simplifies since the screened interactions D X and Yukawa couplings h X tend to 1 and U , respectively, at large frequencies, see Eqs. (21) or (25).
Regarding the Matsubara summation, in general we select a range from 100 to 256 positive frequencies, depending on the physical regime. Contrary to the fRG, in the DMF 2 RG the single-scale propagator decays faster at large frequency, simplifying the convergence in the Matsubara summation. Moreover, it is noteworthy to mention that, as stated in Sec. III A, the computation of the charge susceptibility in the Mott phase requires a larger number of frequencies in the Matsubara summation, which must be then extended, despite the relatively high temperature range, to 1000 to recover its physical value [12].
While part of the momentum dependence is projected onto form-factors as explained in the main text, the transfer momentum dependence has been patched similarly to Ref. [108], retaining 38 patches in the reduced Brillouin zone B red = {(k x , k y ) : 0 ≤ k y ≤ k x ≤ π}.
Finally, the flow equations have been solved using the adaptive Runge-Kutta Cash-Karp 54 method and the momentum integration over the Brillouin zone is carried out via an adaptive cubature technique.