Accessing the ordered phase of correlated Fermi systems: vertex bosonization and mean-field theory within the functional renormalization group

We present a consistent fusion of functional renormalization group and mean-field theory which explicitly introduces a bosonic field via a Hubbard-Stratonovich transformation at the critical scale, at which the order sets in. We show that a minimal truncation of the flow equations, that neglects order parameter fluctuations, is integrable and fulfills fundamental constraints as the Goldstone theorem and the Ward identity connected with the broken global symmetry. To introduce the bosonic field, we present a technique to factorize the most singular part of the vertex, even when the full dependence on all its arguments is retained. We test our method on the two-dimensional attractive Hubbard model at half-filling and calculate the superfluid gap as well as the Yukawa couplings and residual two fermion interactions in the ordered phase as functions of fermionic Matsubara frequencies. Our formalism constitutes a convenient starting point for the inclusion of order parameter fluctuations by keeping a full, non simplified, dependence on fermionic momenta and/or frequencies.


I. INTRODUCTION
When dealing with correlated Fermi systems, one has very frequently to face the breaking of one or more symmetries of the model through the development of some kind of order. Mean-field theory provides a relatively simple but often qualitatively correct description of ground state properties in the ordered phase. Remarkably enough, this is not limited to weak coupling calculations but it may survive at strong coupling [1,2]. On the other hand, fluctuations of the order parameter play a key role at finite temperature T [3] and in low dimensionalities. In particular, they are fundamental in twodimensional systems, where they prevent the formation of long-range order at any T = 0 [4,5], or in the specific case of the breaking of an abelian continuous symmetry group, specifically U(1), they are responsible for the formation of the Berezinskii-Kosterlitz-Thouless (BKT) phase, characterized by quasi long-range order [6,7].
The functional renormalization group (fRG) provides a framework to deal with interacting Fermi systems and ordering tendencies [8]. The inclusion of an infrared cutoff in the bare model allows for the treatment of different energy scales Λ in a unified approach. In the most typical cases of symmetry breaking, as those associated with the onset of magnetic or superfluid/superconducting orders, at high energies the system is in its symmetric phase, while by decreasing the scale Λ, the effective two fermion interaction grows until it reaches a divergence at a scale Λ c in one (or more) specific momentum channel [9][10][11]. In order to continue the flow into the low energy regime, Λ < Λ c , one has to explicitly introduce an order parameter taking into account spontaneous symmetry breaking. Various approaches are possible. One can, for example, decouple the bare interaction via a Hubbard-Stratonovich transformation and run a flow for a mixed boson-fermion system above and below the critical scale. In this way one is able to study fluctuation effects both in the symmetric and in the ordered phases [12][13][14][15][16]. More-over, the two fermion effective interaction generated by the flow can be re-bosonized scale by scale, with a technique called flowing bosonization, either decoupling the bare interaction from the beginning [17,18] or keeping it along the flow, reassigning to the bosonic sector only contributions that arise on top of it [19][20][21]. A different approach to fluctuation effects consists in including below the critical scale Λ c the anomalous terms arising from the breaking of the global symmetry, by keeping only fermionic degrees of freedom [22][23][24][25][26]. If one is not interested in the effects of bosonic fluctuations, as it could be for ground state calculations, a relatively simple truncation of flow equations can reproduce a mean-field (MF) like solution [22,[27][28][29].
Concerning the symmetric phase above the critical scale, recent developments have made the fRG a more reliable method for quantitative and/or strong coupling calculations. We refer, in particular, to the development of the multiloop fRG, that has been proven to be equivalent to the parquet approximation [30][31][32][33] and the fusion of the fRG with the dynamical mean-field theory (DMFT) [34] in the so called DMF 2 RG scheme [35,36]. Within these frameworks, the full dependence of the effective two fermion interaction on all three Matsubara frequencies is often kept.
On the other hand, many efforts have been made in order to reduce the computational complexity of the effective interaction with a full dependence on its fermionic arguments. This is mainly achieved by describing the fermion-fermion interaction process through the exchange of a small number of bosons. Many works treat this aspect not only within the fRG [20,21,37], but also within the DMFT, in the recently introduced single boson exchange approximation [38,39], its nonlocal extensions, the TRILEX approach for example [40], or the dual boson theory [41][42][43][44]. Describing the fermionic interactions in terms of exchanged bosons is important not only to reduce the computational complexity, but also to identify those collective fluctuations that play a fundamental role in the ordered phase.
In this paper, we present a truncation of the fRG flow equations, in which a bosonic field is explicitly introduced, and we prove it to be equivalent to the fusion of the fRG with MF theory introduced in Refs. [28,29]. These flow equations fulfill fundamental constraints as the Goldstone theorem and the global Ward identity connected with spontaneous symmetry breaking (SSB), and they can be integrated, simplifying the calculation of correlation functions in the ordered phase to a couple of self consistent equations, one for the bosonic field expectation value, and another one for the Yukawa coupling between a fermion and the Goldstone mode. In order to perform the Hubbard-Stratonovich transformation, we decompose the effective two fermion interaction in terms of an exchanged boson, which becomes massless at the critical scale, and a residual interaction, and we present a technique to factorize the fRG vertex when its full dependence on fermionic Matsubara frequencies is kept. We prove the feasibility and efficiency of our formalism by applying it to the two-dimensional half-filled attractive Hubbard model, calculating the superfluid gap, Yukawa couplings and residual two fermion interactions in the SSB phase. One notable aspect of our formalism is that the full dependence on fermionic momenta and/or frequencies can be retained. This makes it suitable for a combination with the newly developed methods within the fRG, to continue the flow with a simple truncation in those cases in which the effective two fermion interaction diverges. In the one loop truncation, both in plain fRG [45] and in the DMF 2 RG [36], these divergences are actually found at finite temperature, indicating the onset of spontaneous symmetry breaking. Our method can be also combined with the multiloop fRG, where no divergences are found at finite temperature in 2D [33], to study three-dimensional systems or zero temperature phases. Furthermore, the introduction of the bosonic field makes our method a convenient starting point for the inclusion of order parameter fluctuations on top of the MF, and paves the way for the study of the SSB phases with a full treatment of fermionic Matsubara frequency dependencies.
This paper is organized as follows. In Sec. II we give a short overview of the fRG and its application to correlated Fermi systems. In Sec. III we introduce the attractive Hubbard model, that will be the prototypical model for the application of our method. In Sec. IV we review the MF approximation within the fRG by making use only of fermionic degrees of freedom. In Sec. V we introduce our method by reformulating the fermionic MF approach with the introduction of a bosonic field and we prove the equivalence of the two methods. In Sec. VI we expose a strategy to extract a factorizable part from the effective two fermion interactions, necessary to implement the Hubbard-Stratonovich transformation. This strategy is suitable for the application to the most frequently used schemes within the fRG. In Sec. VII we present some exemplary results for the attractive Hubbard model. A conclusion in Sec. VIII closes the presen-tation.

II. FUNCTIONAL RENORMALIZATION GROUP
In this section we present a short review of the fRG applied to interacting Fermi systems and we refer to Ref. [8] for further details. Providing the bare fermionic action with a regulator R Λ , where the symbol (·, ·) indicates a sum over quantum numbers and fermionic Matsubara frequencies ν = (2j + 1)πT , with j ∈ Z, one can derive an exact differential equation for the effective action as a function of the scale Λ [46,47]: where Γ (2)Λ is the matrix of second derivatives of the effective action w.r.t. the fermionic fields, ∂ Λ is a derivative acting only on the explicit Λ-dependence of R Λ and the trace is intended to run over all the quantum numbers and Matsubara frequencies. In general, the regulator can be any generic function of the scale Λ and the fermionic "d+1 momentum" k = (k, ν) (with k being the spatial momentum), provided that R Λ→Λinit → ∞ and R Λ→Λ fin → 0. In this way, Eq. (2) can be complemented with the initial condition Eq. (2) is however very hard to tackle. A common procedure is to expand the effective action Γ Λ in polynomials of the fields up to a finite order, so that one is limited to work with a finite number of scale dependent couplings. Rather often, in the context of correlated Fermi systems, this truncation is restricted to a flow equation for the self-energy Σ Λ and a vertex V Λ , describing the two fermion effective interaction. The differential equations for these couplings can be inferred directly from Eq. (2). Furthermore, when working with systems that possess U(1) charge, SU(2) spin rotation and translational symmetries, the vertex V Λ as a function of the spin variables σ i and the four d + 1 momenta k i of the fermions (two incoming, two outgoing) can be written as where the fermions labeled as 1-2 are considered as incoming and the ones labeled as 3-4 as outgoing in the scattering process. Furthermore, thanks to translational invariance, the vertex is nonzero only when the total momentum is conserved, that is when k 1 + k 2 = k 3 + k 4 . So, one can shorten the momentum dependence to three momenta, the fourth being fixed by the conservation law. By exploiting the relation above, one is left with the calculation of a single coupling function V Λ that summarizes all possible spin combinations. Its flow equation reads (dropping momentum dependencies for the sake of compactness) where the last term contains the 3-fermion coupling Γ (6)Λ contracted with the single scale propagator ∂ Λ G Λ . This term is often neglected or treated in an approximate fashion in most applications. The remaining three terms can be expressed as loop integrals involving two fermionic propagators and two vertices V Λ . They are grouped in three channels, namely particle-particle (T Λ pp ), particlehole (T Λ ph ) and particle-hole-crossed (T Λ phx ), depending on which combination of momenta is transported by the loop. For the expressions of all the terms in Eq. (5) see Ref. [8].
In numerous applications of the fRG to various systems, the vertex function V Λ diverges before the numerical integration of Eq. (5) reaches the final scale Λ fin . This fact signals the tendency of the system to develop some kind of order by spontaneously breaking one (or more) of its symmetries. One can often trace back the nature of the order tendency by looking at which of the terms in Eq. (5) contributes the most to the flow of V Λ near the critical scale Λ c , where the divergence occurs.

III. MODEL
In this section we present the prototypical model that we use for the application of our method. This is the two-dimensional (2D) attractive Hubbard model, that exhibits an instability in the particle-particle channel, signaling the onset of spin-singlet superfluidity. Our formalism, however, can be extended to a wide class of models, including the 2D repulsive Hubbard model, to study the phases in which (generally incommensurate) antiferromagnetism and/or d-wave superconductivity appear. The bare action of the model describes spin-1 2 fermions on a 2D lattice experiencing an attractive on-site attraction where ν is a fermionic Matsubara frequency, ξ k is the bare band dispersion measured relative to the chemical potential µ, and U < 0 is the local interaction. The symbol k = T ν k d 2 k (2π) 2 (T being the temperature) denotes an integral over the Brillouin zone and a sum over Matsubara frequencies. In the next sections, we will assume that a fRG flow is run for this model, up to a stopping scale Λ s , very close to a critical scale Λ c where the vertex V Λ diverges due to a pairing tendency, but still in the symmetric regime. From now on, we will also assume an infrared regulator such that the scale Λ is lowered from Λ ini to Λ fin , so that the inequality Λ ini > Λ s Λ c > Λ fin holds.

IV. BROKEN SYMMETRY PHASE: FERMIONIC FORMALISM
In this section we will present a simple truncation of flow equations that allows to continue the flow beyond Λ s in the superfluid phase within a MF-like approximation, that neglects any kind of order parameter (thermal or quantum) fluctuations. This approximation can be formulated by working only with the physical fermionic degrees of freedom. In order to tackle the breaking of the global U(1) symmetry, we introduce the Nambu spinors A

. Flow equations and integration
In the SSB phase, the vertex function V acquires anomalous components due to the violation of particle number conservation. In particular, besides the normal vertex describing scattering processes with two incoming and two outgoing particles (V 2+2 ), in the superfluid phase also components with three (V 3+1 ) or four (V 4+0 ) incoming or outgoing particles can arise. We avoid to treat the 3+1 components, since they are related to the coupling of the order parameter to charge fluctuations [24], which do not play any role in a MF-like approximation for the superfluid state. It turns out to be useful to work with combinations that represent two fermion interactions in the longitudinal and transverse order parameter channels, respectively. They are related to the amplitude and phase fluctuations of the superfluid order parameter, respectively. In principle, a longitudinal-transverse mixed interaction can also appear, from the imaginary parts of the vertices in Eq. (8), but it has no effect in the present MF approximation because it vanishes at zero center of mass frequency [48]. Below the stopping scale, Λ < Λ s , we consider a trun-cation of the effective action of the form with the Nambu bilinears defined as where the Pauli matrices τ α are contracted with Nambu spinor indexes. The fermionic propagator G Λ (k) is given by the matrix where is the regulator, Σ Λ (k) is the normal self energy and ∆ Λ (k) is the superfluid gap. The initial conditions at the scale Λ = Λ s require ∆ Λs to be zero and both V Λs A and V Λs Φ to equal the vertex V Λs in the symmetric phase.
We are now going to introduce the MF approximation to the symmetry broken state, that means that we focus on the q = 0 component of V A and V Φ and neglect all the rest. So, from now on we drop all the q-dependencies. We neglect the flow of the normal self-energy below Λ s , that would require the inclusion of charge fluctuations in the SSB phase, which is beyond the MF approximation. In order to simplify the presentation, we introduce a matrixvector notation for the gaps and vertices. In particular, the functions V A and V Φ are matrices in the indices k and k , while the gap and the fermionic propagator behave as vectors. For example, in this notation an object of the type k V Λ A (k, k )∆ Λ (k ) can be viewed as a matrixvector product, V Λ A ∆ Λ . Within our MF approximation, we consider in our set of flow equations only the terms that involve only the q = 0 components of the functions V A and V Φ . This means that in a generalization of Eq. (5) to the SSB phase, we consider only the particle-particle contributions. In formulas we have: where we have defined the bubbles where δ k,k = (2π) 2 /T δ (2) (k−k )δ νν , and the trace runs over Nambu spin indexes. The last terms of Eqs. (12) and (13) involve the 6-particle interaction, which we treat here in the Katanin approximation, that allows us to replace the derivative acting on the regulator ∂ Λ of the bubbles with the full scale derivative ∂ Λ [49]. This approach is useful for it provides the exact solution of mean-field models, such as the reduced BCS, in which the bare interaction is restricted to the zero center of mass momentum channel [22]. In this way, the flow equation (12) for the vertex V A , together with the initial condition V Λs A = V Λs can be integrated analytically, giving where is the (normal) particle-particle bubble at zero center of mass momentum, is the fermionic normal propagator, and is the irreducible (normal) vertex in the particle-particle channel at the stopping scale. The flow equation for the transverse vertex V Φ exhibits a formal solution similar to the one in Eq. (15), but the matrix inside the square brackets is not invertible. We will come to this aspect later.

B. Gap equation
Similarly to the flow equations for vertices, in the flow equation of the superfluid gap we neglect the contributions involving the vertices at q = 0. We are then left with where is the anomalous fermionic propagator, with G defined as in Eq. (17), and with the normal self-energy kept fixed at its value at the stopping scale. By inserting Eq. (15) into Eq. (19) and using the initial condition ∆ Λs = 0, we can analytically integrate the latter, obtaining the gap equation [28] ∆ In the particular case in which the contributions to the vertex flow equation from other channels (different from the particle-particle) as well as the 3-fermion interaction and the normal self-energy are neglected also above the stopping scale, the irreducible vertex is nothing but −U , the (sign reversed) bare interaction, and Eq. (21) reduces to the standard Hartree-Fock approximation to the SSB state.

C. Goldstone Theorem
In this subsection we prove that the present truncation of flow equations fulfills the Goldstone theorem. We revert our attention on the transverse vertex V Φ . Its flow equation in Eq. (13) can be (formally) integrated too, together with the initial condition V Λs Φ = V Λs , giving However, by using the relation one can rewrite the matrix in angular brackets in the second line of Eq. (22) as Multiplying this expression by ∆ Λ (k ) and integrating over k , we see that it vanishes if the gap equation (21) is obeyed. Thus, the matrix in angular brackets in Eq. (22) has a zero eigenvalue with the superfluid gap as eigenvector. In matrix notation this fact can be expressed as Due to the presence of this zero eigenvalue, the above matrix is not invertible. This is nothing but a manifestation of the Goldstone theorem. Indeed, due to the breaking of the global U(1) symmetry, transverse fluctuations of the order parameter become massless at q = 0, leading to the divergence of the transverse two fermion interaction V Φ .

V. BROKEN SYMMETRY PHASE: BOSONIC FORMALISM
The SSB phase can be accessed also via the introduction of a bosonic field, describing the fluctuations of the order parameter, and whose finite expectation value is related to the formation of anomalous components in the fermionic propagator. In order to introduce this bosonic field, we express the vertex at the stopping scale in the following form: We assume from now on that the divergence of the vertex, due to the appearance of a massless mode, is absorbed into the first term, while the second one remains finite.
In other words, we assume that as the stopping scale Λ s approaches the critical scale Λ c at which the vertex is formally divergent, the (inverse) bosonic propagator m Λs (q) at zero frequency and momentum vanishes, while the the Yukawa coupling h Λs (k; q) and the residual two fermion interaction Q Λs (k, k ; q) remain finite.
In Sec. VI we will introduce a systematic scheme to extract the decomposition (26) from a given vertex at the stopping scale.

A. Hubbard-Stratonovich transformation and truncation
The bosonic field φ q can be introduced via a Hubbard-Stratonovich transformation on the factorized (and singular) part of the vertex. At the stopping scale Λ s the effective action will then read Note that we have avoided to introduce an interaction between equal spin fermions. Indeed, since we are focusing on a spin singlet superconducting order parameter, within the MF approximation this interaction has no contribution to the flow equations. We introduce Nambu spinors as in Eq. (7) and we decompose the bosonic field into its (flowing) expectation value plus longitudinal (σ) and transverse (π) fluctuations [16]: where we have chosen α Λ to be real. For the effective action at Λ < Λ s in the SSB phase, we use the following ansatz where the first three quadratic terms are given by and the fermion-boson interactions are with S α k,q as in Eq. (10). The residual two fermion interaction term is written as As in the fermionic formalism, in the truncation in Eq. (29) we have neglected any type of longitudinaltransverse fluctuation mixing in the Yukawa couplings, bosonic propagators and two fermion interactions because at q = 0 they are identically zero. In the bosonic formulation, as well as for the fermionic one, the MF approximation requires to focus on the q = 0 components of the various terms appearing in the effective action and neglect all the rest. So, from now on we drop all the qdependencies. We will make use of the matrix notation introduced in Sec. IV, for which the newly introduced Yukawa couplings behave as vectors and bosonic inverse propagators as scalars.

B. Flow equations and integration
Here we focus on the flow equations for two fermion interactions, Yukawa couplings and bosonic inverse propagators in the longitudinal and transverse channels within a MF approximation, that is we focus only on the Cooper channel (q = 0) and neglect all the diagrams containing internal bosonic lines or the couplings A, Φ at q = 0. Furthermore, we introduce a generalized Katanin approximation to account for higher order couplings in the flow equations. We refer to Appendix A for a derivation of the latter. We now show that our reduced set of flow equations for the various couplings can be integrated. We first focus on the longitudinal channel, while in the transverse one the flow equations possess the same structure.
The flow equation for the longitudinal bosonic mass (inverse propagator at q = 0) reads Similarly, the equation for the longitudinal Yukawa coupling is and the one for the residual two fermion longitudinal interaction is given by The above flow equations are pictorially shown in Fig. 1.
The initial conditions at Λ = Λ s read, for both channels, We start by integrating the equation for the residual two fermion longitudinal interaction A. Eq. (35) can be solved exactly as we have done in the fermionic formalism, obtaining for A where we have introduced a reduced residual two fermion interaction Q We are now in the position to employ this result and plug it in Eq. (34) for the Yukawa coupling. The latter can be integrated as well. Its solution reads where the introduction of a "reduced" Yukawa coupling is necessary. This Bethe-Salpeter-like equation for the Yukawa coupling is similar in structure to the parquetlike equations for the three-leg vertex derived in Ref. [39]. Finally, we can use the two results of Eqs. (37) and (39) and plug them in the equation for the bosonic mass, whose integration provides where, by following definitions introduced above, the "reduced" bosonic mass is given by In the transverse channel, the equations have the same structure and can be integrated in the same way. Their solutions read Eq. (45) provides the mass of the transverse mode, which, according to the Goldstone theorem, must be zero. We will show later that this is indeed fulfilled. It is worthwhile to point out that the combinations obey the same flow equations, Eqs. (12) and (13), as the vertices in the fermionic formalism and share the same initial conditions. Therefore the solutions for these quantities coincide with expressions (15) and (22), respectively. Within this equivalence, it is interesting to express the irreducible vertex V Λs of Eq. (18) in terms of the quantities, Q Λs , h Λs and m Λs , introduced in the factorization in Eq. (26): where Q Λs , h Λs and m Λs were defined in Eqs. (38), (40) and (42). For a proof see Appendix B. Relation (47) is of particular interest because it states that when the full vertex is expressed as in Eq. (26), then the irreducible one will obey a similar decomposition, where the bosonic propagator, Yukawa coupling and residual two fermion interaction are replaced by their "reduced" counterparts. This relation holds even for q = 0.

C. Ward identity for the gap
We now focus on the flow of the fermionic gap and the bosonic expectation value and express a relation that connects them. Their flow equations are given by (see Appendix A) and with F Λ given by Eq. (20). In Fig. 2 we show a pictorial representation. Eq. (48) can be integrated, with the help of the previously obtained results for A, h σ and m σ , yielding This equation is the Ward identity for the mixed bosonfermion system related to the global U(1) symmetry [16].
In Appendix C we propose a self consistent loop for the calculation of α, h π , through Eqs. 50 and 44, and subsequently the superfluid gap ∆.

D. Goldstone theorem
Let us now come back to the question of the Goldstone theorem. For the mass of the Goldstone boson to be zero, it is necessary for Eq. (45) to vanish. We show that this is indeed the case. With the help of Eq. (23), we can reformulate the equation for the transverse mass in the form where the Ward Identity ∆ = αh π was applied in the last line. We see that the expression for the Goldstone boson mass vanishes when α obeys its self consistent equation, Eq. (50). This proves that our truncation of flow equations fulfills the Goldstone theorem.

E. Equivalence of bosonic and fermionic formalisms
As we have proven in the previous sections, within the MF approximation the fully fermionic formalism of Sec. IV and the bosonized approach introduced in the present section provide the same results for the superfluid gap and for the effective two fermion interactions. Notwithstanding the formal equivalence, the bosonic formulation relies on a further requirement. In Eqs. (43) and (44) we assumed the matrix 1 − Q Λs Π Λ 22 to be invertible. This statement is exactly equivalent to assert that the two fermion residual interaction Φ remains finite. Otherwise the Goldstone mode would lie in this coupling and not (only) in the Hubbard-Stratonovich boson. This fact cannot occur if the flow is stopped at a scale Λ s coinciding with the critical scale Λ c at which the (normal) bosonic mass m Λ turns zero, but it could take place if one considers symmetry breaking in more than one channel. In particular, if one allows the system to develop two different orders and stops the flow when the mass of one of the two associated bosons becomes zero, it could happen that, within a MF approximation for both order types, the appearance of a finite gap in the first channel makes the two fermion transverse residual interaction in the other channel diverging. In that case one can apply the technique of the flowing bosonization [20,21], by reassigning to the bosonic sector the (most singular part of the) two fermion interactions that are generated during the flow. It can be proven that also this approach gives the same results for the gap and the effective fermionic interactions in Eq. (46) as the fully fermionic formalism.

VI. VERTEX BOSONIZATION
In this section we present a systematic procedure to extract the quantities in Eq. (26) from a given vertex, within an approximate framework. The full vertex in the symmetric phase can be written as [50,51] where V Λini is the vertex at the initial scale, and we call φ p pairing channel, φ m magnetic channel and φ c charge channel. Each of this functions depends on a bosonic and two fermionic variables. Within the so called 1-loop approximation, where one neglects the 3-fermion coupling in Eq. (5), in the Katanin scheme [49], or in more involved schemes, such as the 2-loop [52] or the multiloop [30,31], one is able to assign one or more of the terms of the flow equation (5) for V Λ to each of the channels, in a way that their last bosonic argument enters only parametrically in the formulas. This is the reason why the decomposition in Eq. (53) is useful. The vertex at the initial scale can be set equal to the bare (sign-reversed) Hubbard interaction −U in a weak-coupling approximation, or as in the recently introduced DMF 2 RG scheme, to the vertex computed via DMFT [35,36]. In order to simplify the treatment of the dependence on fermionic spatial momenta of the various channels, one often introduces a complete basis of Brillouin zone form factors {f k } and expands each channel in this basis [53] with X = p, m or c, and sgn p = −1, sgn c = sgn m = +1. For practical calculations the above sum is truncated to a finite number of form factors and often only diagonal terms, = , are considered. Within the form factor truncated expansion, one is left with the calculation of a finite number of channels that depend on a bosonic "d+1 momentum" q = (q, Ω) and two fermionic Matsubara frequencies ν and ν .
We will now show how to obtain the decomposition introduced in Eq. (26) within the form factor expansion. We focus on only one of the channels in Eq. (53), depending on the type of order we are interested in, and factorize its dependence on the two fermionic Matsubara frequencies. We introduce the so called channel asymptotics, that is the functions that describe the channels for large ν, ν . From now on we adopt the shorthand lim ν→∞ g(ν) = g(∞) for whatever g, function of ν. By considering only diagonal terms in the form factor expansion in Eq. (54), we can write the channels as [54]: with According to Ref. [54], these functions are related to physical quantities. K (1) X, turns out to be proportional to the relative susceptibility and the combination K (1) X, ) to the so called boson-fermion vertex, that describes both the response of the Green's function to an external field [55] and the coupling between a fermion and an effective boson. In principle one should be able to calculate the above quantities diagrammatically (see Ref. [54] for the details) without performing any limit. However, it is well known how fRG truncations, in particular the 1-loop approximation, do not properly weight all the Feynman diagrams contributing to the vertex, so that the diagrammatic calculation and the high frequency limit give two different results. To keep the property in the last line of Eq. (56), we choose to perform the limits. We rewrite Eq. (55) in the following way: where we have made the frequency and momentum dependencies explicit only in the second line and we have defined From the definitions given above, it is obvious that the rest function R X, decays to zero in all frequency directions.
Since the first term of Eq. (57) is separable by construction, we choose to identify this term with the first one of Eq. (26). Indeed, in many cases the vertex divergence is manifest already in the asymptotic K (1)Λ X, , that we recall to be proportional to the susceptibility of the channel. There are however situations in which the functions K (1) and K (2) are zero even close to an instability in the channel, an important example being the d-wave superconducting instability in the repulsive Hubbard model. In general, this occurs for those channels that, within a Feynman diagram expansion, cannot be constructed with a ladder resummation with the bare vertex. In the Hubbard model, due to the locality of the bare interaction, this happens for every = 0, that is for every term in the form factor expansion different than the s-wave contribution. In this case one should adopt a different approach and, for example, replace the limits to infinity in Eq. (57) by some given values of the Matsubara frequencies, ±πT for example.

VII. RESULTS FOR THE ATTRACTIVE HUBBARD MODEL AT HALF-FILLING
In this section we report some exemplary results of the equations derived within the bosonic formalism, for the attractive two-dimensional Hubbard model. We focus on the half-filled case. For pure nearest neighbors hopping with amplitude −t, the band dispersion ξ k is given by with µ = 0 at half-filling. We choose the onsite attraction and the temperature to be U = −4t and T = 0.1t. All results are presented in units of the hopping parameter t.

A. Symmetric phase
In the symmetric phase, in order to run a fRG flow, we introduce the Ω-regulator [50] so that the initial scale is Λ init = +∞ (fixed to a large number in the numerical calculation) and the final one Λ fin = 0. We choose a 1-loop truncation, that is we neglect the last term of Eq. (5), and use the decomposition in Eq. (53) with a form factor expansion. We truncate Eq. (54) only to the first term, that is we use only s-wave, f k ≡ 1, form factors. Within these approximations, the vertex reads where P, M, C, are referred as pairing, magnetic and charge channel, respectively. Furthermore, we focus only on the spin-singlet component of the pairing (the triplet one is very small in the present parameter region), so that we require the pairing channel to obey [56] P Λ νν (q) = P Λ Ω−ν,ν (q) = P Λ ν,Ω−ν (q), with q = (q, Ω). The initial condition for the vertex reads so that P Λinit = M Λinit = C Λinit = 0. Neglecting the fermionic self-energy, Σ Λ (k) ≡ 0, we run a flow for these three quantities until one (ore more) of them diverges. Under a technical point of view, each channel is computed by keeping 50 positive and 50 negative values for each of the three Matsubara frequencies (two fermionic, one bosonic) on which it depends. Frequency asymptotics are also taken into account, by following Ref. [54]. The momentum dependence of the channel is treated by discretizing with 38 patches the region B = {(k x , k y ) : 0 ≤ k y ≤ k x ≤ π} in the Brillouin zone and extending to the other regions by using lattice symmetries. The expressions of the flow equations are reported in Appendix D.
Due to particle-hole symmetry occurring at half-filling, pairing fluctuations at q = 0 combine with charge fluctuations at q = (π, π) to form an order parameter with SO(3) symmetry [57]. Indeed, with the help of a canonical particle-hole transformation, one can map the attractive half-filled Hubbard model onto the repulsive one. Within this duality, the SO(3)-symmetric magnetic order parameter is mapped onto the above mentioned combined charge-pairing order parameter and vice versa. This is the reason why we find a critical scale, Λ c , at which both C((π, π), 0) and P(0, 0) diverge, as shown in Fig. 3. On a practical level, we define the stopping scale Λ s as the scale at which one (or more, in this case) channel exceeds 10 3 t. With our choice of parameters, we find that at Λ s 0.378t both C and P cross our threshold. In the SSB phase, we choose to restrict the ordering to the pairing channel, thus excluding the formation of charge density waves. This choice is always possible because we have the freedom to choose the "direction" in which our order parameter points. Indeed, in the particle-hole dual repulsive model, our choice would be equivalent to choose the (antiferro-) magnetic order parameter to lie in the xy plane. This choice is implemented by selecting the particle-particle channel as the only one contributing to the flow in the SSB phase, as exposed in Secs. IV and V.
In order to access the SSB phase with our bosonic formalism, we need to perform the decomposition in Eq. (26) for our vertex at Λ s . Before proceeding, in order to be consistent with our form factor expansion in the SSB phase, we need to project V in Eq. (61) onto the swave form factors, because we want the quantities in the ordered phase to be functions of Matsubara frequencies only. Therefore we define the total vertex projected onto s-wave form factors Furthermore, since we are interested only in spin singlet pairing, we symmetrize it with respect to one of the two fermionic frequencies, so that in the end we are dealing  .
In order to extract the Yukawa coupling h Λs and bosonic propagator m Λs , we employ the strategy described in Sec. VI. Here, however, instead of factorizing the pairing channel P Λs alone, we subtract from it the bare interaction U . In principle, U can be assigned both to the pairing channel, to be factorized, or to the residual two fermion interaction, giving rise to the same results in the SSB phase. However, when in a real calculation the vertices are calculated on a finite frequency box, it is more convenient to have the residual two fermion interaction Q Λs as small as possible, in order to reduce finite size effects in the matrix inversions needed to extract the reduced couplings in Eqs. (38), (40) and (42), and in the calculation of h π , in Eq. (44). Furthermore, since it is always possible to rescale the bosonic propagators and Yukawa couplings by a constant such that the vertex constructed with them (Eq. (57)) is invariant, we impose the normalization condition h Λs (ν → ∞; q) = 1.
In formulas, we have m Λs (q) = 1 and h Λs (ν; q) = K (2)Λs p, =0 (ν; q) + K (1)Λs The limits are numerically performed by evaluating the pairing channel at large values of the fermionic frequencies. The extraction of the factorizable part from the pairing channel minus the bare interaction defines the rest function and the residual two fermion interaction Q We are now in the position to extract the reduced couplings, Q Λs , h Λs and m Λs , defined in Eqs. (38), (40), (42). This is achieved by numerically inverting the matrix (we drop the q-dependence from now on, assuming always q = 0) with and In Fig. 4 we show the results for the pairing channel minus the bare interaction, the rest function, the residual two fermion interaction Q and the reduced one Q at the stopping scale. One can see that in the present parameter region the pairing channel (minus U ) is highly factorizable. Indeed, despite the latter being very large because of the vicinity to the instability, the rest function R remains very small, a sign that the pairing channel is well described by the exchange of a single boson. Furthermore, thanks to our choice of assigning the bare interaction to the factorized part, as we see in Fig. 4, both Q and Q possess frequency structures that arise from a background that is zero. Finally, the full bosonic mass at the stopping scale is close to zero, m Λs 10 −3 , due to the vicinity to the instability, while the reduced one is finite, m Λs 0.237. Lower left: residual two fermion interaction. The choice of factorizing P Λs − U instead of P Λs alone makes the background of this quantity zero. Lower right: reduced residual two fermion interaction. As well as the full one, this coupling has a zero background value, making calculations of couplings in the SSB phase more precise by reducing finite size effects in the matrix inversions.

B. SSB Phase
In the SSB phase, instead of running the fRG flow, we employ the analytical integration of the flow equations described in Sec. V. On a practical level, we implement a solution of the loop described in Appendix C, that allows for the calculation of the bosonic expectation value α, the transverse Yukawa coupling h π and subsequently the fermionic gap ∆ through the Ward identity ∆ = αh π . In this section we drop the dependence on the scale, since we have reached the final scale Λ fin = 0. Note that, as exposed previously, in the half-filled attractive Hubbard model the superfluid phase sets in by breaking a SO(3) rather than a U(1) symmetry. This means that one should expect the appearance of two massless Goldstone modes. Indeed, besides the Goldstone boson present in the (transverse) particle-particle channel, another one appears in the particle-hole channel and it is related to the divergence of the charge channel at momentum (π, π). However, within our choice of considering only superfluid order and within the MF approximation, this mode is decoupled from our equations. Within our previously discussed choice of bosonizing P Λs − U instead of P Λs alone, the self consistent loop introduced in Appendix C converges extremely fast, 15 iterations for example are sufficient to reach a precision of 10 −7 in α. Once convergence is reached and the gap ∆(ν) obtained, we are in the position to evaluate the remaining couplings introduced in Sec. V through their integrated flow equations. In Fig. 5 we show the computed frequency dependence of the gap. It interpolates between ∆ 0 = ∆(ν → 0), its value at the Fermi level, and its asymptotic value, that equals the (signed reversed) bare interaction times the Cooper pair expectation value k ψ −k,↓ ψ k,↑ . ∆ 0 also represents the gap between the upper and lower Bogoliubov bands. Magnetic and charge fluctuations above the critical scale significantly renormalize the gap with respect to the Hartree-Fock calculation ( V = −U in Eq. (21)), that in the present case coincides with Bardeen-Cooper-Schrieffer (BCS) theory. Since ∆ is a spin singlet superfluid gap, and we have chosen α to be real, it obeys where the first equality comes from the spin singlet nature and the second one from the time reversal symmetry of the effective action. Therefore, the imaginary part of the gap is always zero. By contrast, a magnetic gap would gain, in general, a finite (and antisymmetric in frequency) imaginary part. In Fig. 6 we show the results for the residual two fermion interactions in the longitudinal and transverse channels, together with the total effective interaction in the longitudinal channel, defined as The analogue of Eq. (74) for the transverse channel cannot be computed, because the transverse mass m π is zero, in agreement with the Goldstone theorem. The key result is that the residual interactions A νν and Φ νν inherit the frequency structures of Q Λs νν and Q Λs νν , respectively, and they are also close to them in values (compare with Fig. 4). The same occurs for the Yukawa couplings, as shown in Fig. 7. Indeed, the calculated transverse coupling h π does not differ at all from the Yukawa coupling at the stopping scale h Λs . In other words, if instead of solving the self consistent equations, one runs a flow in the SSB phase, the transverse Yukawa coupling will stay the same from Λ s to Λ fin . Furthermore, the longitudinal coupling h σ develops a dependence on the frequency which does not differ significantly from the one of h Λs . This feature, at least for our choice of parameters, can lead to some simplifications in the flow equations of Sec. V. Indeed, when running a fRG flow in the SSB phase, one might let flow only the bosonic inverse propagators by keeping the Yukawa couplings and residual interactions fixed at their values, reduced or not, depending on the channel, at the stopping scale. This fact can be crucial to make computational costs lighter when including bosonic fluctuations of the order parameter, which, similarly, do not significantly renormalize Yukawa couplings in the SSB phase [16].

VIII. CONCLUSION
We have introduced a truncation of fRG flow equations which, with the introduction of a Hubbard-Stratonovich boson, has been proven to be equivalent to the MF equations obtained in Refs. [28,29]. These flow equations satisfy fundamental requirements such as the Goldstone theorem and the Ward identities associated with global symmetries, and can be integrated analytically, reducing the calculation of correlation functions in the SSB phase to a couple of self consistent equations for the bosonic expectation value α and the transverse Yukawa coupling h π . A necessary step to perform the Hubbard-Stratonovich transformation, on which our method relies, is to extract a factorizable dependence on fermionic variables k and k from the vertex at the critical scale. A strategy to accomplish this goal has been suggested for a vertex whose dependence on spatial momenta k and k is treated by a form factor expansion, making use of the vertex asymptotics introduced in Ref. [54]. Furthermore, we have tested the feasibility and efficiency of our method on a prototypical model, namely the half-filled attractive Hubbard model in two dimensions, focusing on frequency dependencies of the two fermion interactions, Yukawa couplings and fermionic gap. We have found a good convergence of the iterative scheme proposed. The remaining couplings introduced in our method have been computed after the loop convergence from their integrated flow equations.
Our method leaves room for applications and extensions. First, one can directly apply the MF method, as formulated in this paper, to access the SSB phase in those calculations for which the dependencies on fermionic momenta and/or frequencies cannot be neglected. Some ex-amples are the fRG calculations with a full treatment of fermionic frequencies, within a 1-loop truncation [45], in the recent implementations of multiloop fRG [32,33] or in the DMF 2 RG scheme [36]. These combinations can be applied to two-or three-dimensional systems. In the former case, even though in 2D order parameter fluctuations are expected to play a decisive role, our method can be useful to get a first, though approximate, picture of the phase diagram. Of particular relevance is the 2D repulsive Hubbard model, used in the context of high-Tc superconductors. An interesting system for the latter case, where bosonic fluctuations are expected to be less relevant, is the 3D attractive Hubbard model, which, thanks to modern techniques, can be experimentally realized in cold atoms setups.
Secondly, our method constitutes a convenient starting point for the inclusion of bosonic fluctuations of the order parameter, as done for example in Refs. [16,21], with the full dependence of the gap, Yukawa couplings and vertices on the fermionic momenta and/or frequencies being kept. In particular, by providing the Hubbard-Stratonovich boson with its own regulator, our MF-truncation of flow equations can be extended to include order parameter fluctuations, which in two spatial dimensions and at finite temperature restore the symmetric phase, in agreement with the Mermin-Wagner theorem. One may also adapt the bosonic field at every fRG step through the flowing bosonization [20,21]. This can be done by keeping the full frequency dependence of the vertex and Yukawa coupling, by applying the strategy discussed in Sec. VI to the flow equation for the vertex.
Finally, our MF method does not necessarily require a vertex coming from a fRG flow. In particular, one can employ the DMFT vertex, extract the pairing channel [56] (or any other channel in which symmetry breaking occurs) from it, and apply the same strategy as described in this paper to extract Yukawa and other couplings. This application can be useful to compute those transport quantities [58] and response functions in the SSB phase which, within the DMFT, require a calculation of vertex corrections [34]. The anomalous vertices can be computed also at finite q with a simple generalization of our formulas.

ACKNOWLEDGMENTS
We are grateful to W. Metzner and D. Vilardi for stimulating discussions and a critical reading of the manuscript. In this section we will derive the flow equations used in Sec. V.
We consider only those terms in which the dependence on the center of mass momentum q is fixed to zero by the topology of the relative diagram or that depend only parametrically on it. These diagrams are the only ones necessary to reproduce the MF approximation.
The flow equations will be derived directly from the Wetterich equation (2), with a slight modification, since we have to keep in mind that the bosonic field φ acquires a scale dependence due to the scale dependence of its expectation value. The flow equation reads (for real α Λ ): where Γ (2)Λ is the matrix of the second derivatives of the action with respect to the fields, the supertrace Str includes a minus sign when tracing over fermionic variables. The first equation we derive is the one for the flowing expectation value α Λ . This is obtained by requiring that the one-point function for σ q vanishes. Taking the σ q derivative in Eq. (A1) and setting the fields to zero, we have where we have defined From Eq. (A2) we get the flow equation for α Λ .
The MF flow equation for the fermionic gap reads tions for the transverse couplings:

(A15)
Appendix B: Calculation of the irreducible vertex in the bosonic formalism In this appendix we provide a proof of Eq. (47) by making use of matrix notation. If the full vertex can be decomposed as in Eq. (26) we can plug this relation into the definition of the irreducible vertex, Eq. (18). With some algebra we obtain where in the last line we have inserted a representation of the identity, in between the two matrices and we have made use of definitions (38) and (40). With a bit of simple algebra, we can analytically invert the matrix on the left in the last line of Eq. (B2), obtaining Where m is defined in Eq. (42). By plugging this result into Eq. (B2), we finally obtain that is the result of Eq. (47).
Appendix C: Algorithm for the calculation of the superfluid gap The formalism described in Sec. V allows us to formulate a minimal set of closed equations required for the calculation of the gap. We drop the Λ superscript, assuming that we have reached the final scale. The gap can be computed using the Ward identity, so we can reduce ourselves to a single self consistent equation for α, that is a single scalar quantity, and another one for h π , momentum dependent. The equation for α is Eq. (50). The transverse Yukawa coupling is calculated through Eq. (44). The equations are coupled since the superfluid gap ∆ = αh π appears in the r.h.s of both. We propose an iterative loop to solve the above mentioned equations. By starting with the initial conditions α (0) = 0 and h (0) π (k) = 0, we update the transverse Yukawa coupling at every loop iteration i according to Eq. (44), that can be reformulated in the following algorithmic form: with the matrix M (i) defined as This formulation of self consistent equations is not computationally lighter than the one in the fermionic formalism, but more easily controllable, as one can split the frequency and momentum dependence of the gap (through h π ) from the strength of the order (α). Moreover, thanks to the fact that h π is updated with an explicit expression, namely Eq. (C1), that is in general a well behaved function of k, the frequency and momentum dependence of the gap is assured to be under control.