Superconductivity in the attractive Hubbard model: functional renormalization group analysis

We present a functional renormalization group analysis of superconductivity in the ground state of the attractive Hubbard model on a square lattice. Spontaneous symmetry breaking is treated in a purely fermionic setting via anomalous propagators and anomalous effective interactions. In addition to the anomalous interactions arising already in the reduced BCS model, effective interactions with three incoming legs and one outgoing leg (and vice versa) occur. We accomplish their integration into the usual diagrammatic formalism by introducing a Nambu matrix for the effective interactions. From a random-phase approximation generalized through use of this matrix we conclude that the impact of the 3+1 effective interactions is limited, especially considering the effective interactions important for the determination of the order parameter. The exact hierarchy of flow equations for one-particle irreducible vertex functions is truncated on the two-particle level, with higher-order self-energy corrections included in a scheme proposed by Katanin. Using a parametrization of effective interactions by patches in momentum space, the flow equations can be integrated numerically to the lowest scales without encountering divergences. Momentum-shell as well as interaction-flow cutoff functions are used, including a small external field or a large external field and a counterterm, respectively. Both approaches produce momentum-resolved order parameter values directly from the microscopic model. The size of the superconducting gap is in reasonable agreement with expectations from other studies.


Introduction
Renormalization group (RG) techniques facilitate the systematic investigation of interacting many-particle systems even in situations where perturbation theory breaks down. A central idea underlying many renormalization group methods, pioneered by Wilson [1], is to take into account the degrees of freedom successively, ordered by their kinetic energy. Formally, this is accomplished by integrating a set of differential equations called flow equations, where some energy scale (usually a cutoff) plays the role of the flow parameter. The RG flow generates a family of effective actions which interpolates smoothly between the bare action of the system and the final effective action where all degrees of freedom have been integrated out.
The first significant applications of RG methods to interacting Fermi systems in more than one dimension were undertaken by mathematicians aiming at a rigorous control of certain systems and their properties [2,3,4]. The usefulness of RG concepts for Fermi systems was soon recognized also by physicists [5,6]. RG calculations for Fermi systems are often termed "functional" as they involve the flow of functions instead of just a few running couplings. An early computational success of the functional RG (fRG) was unambiguous evidence for d-wave superconductivity in the repulsive twodimensional Hubbard model [7,8,9].
Of particular interest in interacting Fermi systems is the spontaneous breaking of symmetries. Formally, such a symmetry breaking occurs if a product of two fermionic fields -a bilinear -which does not appear in the original Hamiltonian acquires a nonzero expectation value. This bilinear can be replaced by a bosonic field through a Hubbard-Stratonovich transformation, which simplifies the treatment of the order parameter and its fluctuations [10]. However, these simplifications come at the price of introducing a bias towards breaking the symmetry in the channel fixed by the Hubbard-Stratonovich field. Working with the original fermionic degrees of freedom this bias can be avoided. However, the symmetry-broken phase is not easily accessible by a purely fermionic RG flow. Calculations starting from the symmetric phase run into divergences at a finite energy scale, that is before all degrees of freedom have been integrated out.
Most previous functional RG studies of interacting Fermi systems are based on a flow equation for the generator of the one-particle irreducible (1PI) vertex functions [11]. This equation is equivalent to an infinite hierarchy of flow equations for the 1PI vertex functions. Truncating the hierarchy in a weak coupling expansion, the RG flow has been computed for various low-dimensional Fermi systems [8,12,13,14,15,16]. Instead of the common choice of a momentum cutoff, it was found to be advantageous in certain situations to use as flow parameter the physical temperature [18], the strength of the interaction [19], or a cutoff in frequency space [20,21].
Spontaneous symmetry breaking within the functional RG framework has frequently been studied by introducing bosonic fields for the order parameter, leading to a coupled flow of bosons and fermions [22,23,24]. A simple way to access spontaneous symmetry breaking in a purely fermionic setting is to combine the fRG with a mean-field calculation for the low-energy degrees of freedom [17]. However, order parameter fluctuations are thereby discarded. To continue the fermionic functional RG flow into the symmetry-broken phase, an improved truncation of the hierarchy with an efficient treatment of self-energy corrections, first formulated by Katanin [25], is very promising. Employing this improvement, the Bardeen-Cooper-Schrieffer (BCS) mean-field model can be solved exactly already by a low-order truncation [26], if a small external field is used to trigger the symmetry-breaking. The vaunted divergence is then regularized, but there is still a flow regime with strong coupling. The method was extended to the breaking of a discrete symmetry and finite temperatures [27]. Employing the interaction flow [19], the external field can be compensated by a counterterm, thus enabling the study of first-order phase transitions which were previously out of the method's reach [28]. Furthermore, this approach yields results at zero external field without having to pass through a strong coupling regime in the case of discrete-symmetry breaking. Strong attractors for the flows of various counterterm strengths were found.
In this work, we apply the fermionic fRG to the attractive Hubbard model, using both the external-field as well as the counterterm approach. The repulsive Hubbard model [29,30,31] was originally introduced to study ferromagnetism of itinerant electrons, but has since been employed to elucidate many aspects of strong correlations in solid-state physics. The attractive version exhibits s-wave superconductivity under certain conditions, including the case we study here, namely for a square lattice at zero temperature and weak or moderate coupling [32]. It is a good testbed because the correct order parameter is affected by fluctuations even in the weak coupling limit [33], while the model is fairly well-studied, allowing comparisons of the results obtained using the new method. Lattice fermions with attractive interactions can nowadays be realized by cold atoms in optical lattices, and s-wave superfluidity in such systems has been observed very recently [34].
Our goals in this work are threefold. First, we would like to study the role of anomalous effective interactions with three ingoing particles and one outgoing particle (or vice versa), which arise in systems with both Cooper and forward scattering. We find that these interactions have a limited impact. Second, we would like to determine whether the application of the Katanin-truncated fRG to the Hubbard model works as well as for mean-field models. For the momentum-shell approach with a small external field, we obtain reasonable flows resulting in order parameter values comparable to the literature. It is possible to choose the strength of the external field two orders of magnitude below the final order parameter value without running into artificial divergences. For the interaction-flow approach, we find remnants of the strong-attractor behavior known from the mean-field models. Third, we compute the superconducting order parameter in the ground state of the attractive Hubbard model as a function of coupling and filling. This paper is structured as follows. In section 2, the attractive Hubbard model is introduced formally. The formalism permitting the study of symmetry breaking in general and the 3 + 1 anomalous effective interactions in particular is detailed in section 3. It is employed in the context of a generalized random-phase approximation (RPA) of the model's behavior in section 4. Section 5 presents the result of the fRG studies, containing momentum-shell and interaction flows as well as a study of the order parameter strength and a comparison to earlier calculations.

Attractive Hubbard model
We discuss the attractive Hubbard model on a two-dimensional square lattice specified by the Hamiltonian N denotes the number of lattice points, U 0 > 0 is the onsite attraction, and k takes all values in the first Brillouin zone of the square lattice. The dispersion relation contains a probability amplitude t for nearest-neighbor hopping, a probability amplitude t ′ for next-nearest-neighbor hopping, and a chemical potential µ. An external field is added to H to trigger the instability towards superconductivity. This field is either kept small or turned off gradually in the fRG flow to approximate the situation with spontaneous symmetry breaking.

Functional integral and field substitution
We study the partition function Z in a functional integral representation within the Matsubara formalism, employing the grand canonical action where ψ andψ are the Grassmann fields matching c and c † , respectively. Furthermore, k = (iω n , k) and volume as well as temperature factors are implicitly contained in the summation symbol , a convention we also adopt in the following. We employ the substitutionφ

Propagator and effective interaction
In (7), the external field appears on the off-diagonal of a matrix specifying the action's quadratic part. According to the Dyson equation, inverting this matrix after subtracting the self-energy leads to the propagator or one-particle Green's function where E(k) = (ξ k − Σ(k)) 2 + |∆| 2 and Σ(k) is the normal self-energy. ∆ denotes the superconducting order parameter and contains the anomalous self-energy as well as the external field. To simplify the diagrammatics, we antisymmetrize the interaction part where i, if used as an index, is short for k i N i and 2 and 4 are exchanged in comparison with (7). If (N 1 N 3 ) labels the rows and (N 2 N 4 ) the columns of a matrix to the Nambu base {(−−), (−+), (+−), (++)}, the antisymmetrized bare interaction reads By employing the Nambu field substitution (6), a formalism containing four potentially different normal (particle-number conserving) and twelve potentially different anomalous (particle-number non-conserving) effective interactions arises. Here, as we concentrate on conventional s-wave pairing, we require time reversal invariance to hold. This implies the following behavior of an effective interaction component under flipping of all its Nambu indices: if an odd number of the Nambu indices are (+), the sign of the effective interaction's real part is flipped; otherwise, the sign of the imaginary part is flipped.
Exploiting this symmetry, we assign symbols to the various effective interactions: The Ω i are collectively referred to as 3 + 1 effective interactions as they have three incoming and one outgoing leg or vice versa if drawn according to the ψ-field representation of (5). They are generated in the symmetry-broken phase by diagrams of the type in figure 1. W is a 4 + 0 anomalous effective interaction arising only in the symmetry-broken phase. V , X, and U are normal 2 + 2 effective interactions. Only V and U have non-zero bare values.

Mean-field and generalized RPA
Here we employ a mean-field approximation of the self-energy and a generalized randomphase approximation (RPA) of the effective interactions. This represents a useful first step toward the fRG treatment, as it allows us to assess the basic behavior and the frequency-and momentum-dependence of normal and anomalous vertices. In particular, we find that anomalous 3+1 vertices have a sizable feedback on certain effective interactions.

Gap equation
The system's self-energy can be calculated approximatively by resumming a subset of perturbation theory diagrams as in figure 2. This subset contains all diagrams corresponding to bubbles, sunrise diagrams, and trees of either or both. It is equivalent to self-consistent Hartree-Fock theory. If we were to restrict the Hubbard Hamiltonian to only include the Cooper channel and forward scattering processes with zero momentum transfer, the resulting diagrammatic self-consistency equation would be exact. Analogous statements can be made about the one-loop fRG treatment for such a restricted Hamiltonian [35]. The self-energy equation can be expanded into two distinct parts, The Matsubara sums are performed analytically. The result is In this approximation, both Σ and ∆ are constant in momentum space if the external field is. The Hartree contribution Σ can be absorbed into the chemical potential and is subsequently ignored.

Bethe-Salpeter equation
A subset of all diagrams of the perturbation expansion for the effective interaction can be resummed into a Bethe-Salpeter equation as shown in figure 3. This corresponds to a generalization of the random-phase approximation in the sense that particleparticle ladders in the ψ-representation are included in addition to the usual particle-hole bubbles. Symbolically, this Bethe-Salpeter equation reads where, as a parameter, i = 1 . . . 4 is short for p i N i . By the introduction of the loop matrix and by exploiting momentum and energy conservation p 1 + p 2 = p 3 + p 4 while setting p = (p, iν m ) := p 1 − p 3 , (15) can be rewritten as where the Nambu indices of U and U 0 are organized into pairs as in (11). U(p 1 . . . p 4 ) only depends on p in this approximation, as can be seen by solving for U, which yields

Numerical setup
The Matsubara sums in (18) can be performed analytically. The integration over the Brillouin zone and the matrix inversion on the right-hand side of (18) is accomplished numerically. In the following, a system with t ′ = t/6, µ = −1.5t, U 0 = 2t, and T = 0.035t is assumed. This choice of parameters corresponds to an approximately quarter-filled system far below its mean-field critical temperature. ∆ is chosen to be real. We mostly employ an external field ∆ ext = 10 −4 t. This keeps the numerics wellbehaved while having no appreciable influence on U(p) except near p = 0.

3 + 1 effective interactions
The Ω i are generated by diagrams as in figure 1. In the reduced mean-field models studied previously these vertices remain zero due to the restricted momentum-structure of the bare interaction where only particle pairs with total momentum equal to zero interact. However, in the general case with scattering of particle pairs with arbitrary total momentum they have to be considered. Due to the appearance of an anomalous propagator figure 1, these diagrams and hence the Ω i are non-zero only in the symmetrybroken phase. To obtain an overview of the momentum and frequency dependence of the 3 + 1 effective interactions, we study the graphs in figure 4. The right-hand graphs show that the 3 + 1 effective interactions exhibit no divergence in the p 1 − p 3 = 0 channel.
The Ω i at larger frequencies are suppressed by decreasing the external field. From the figure's left-hand diagram, we learn that the Ω i 's moduli are large around p = 0 and around p = (π, π) for ν = 0.

Large interaction vertices and Goldstone boson
The instability toward superconductivity gives rise to a growth of certain combinations of the interaction vertices. From previous fRG studies on reduced mean-field models [26,27] and from the solution of the Bethe-Salpeter equation, the mean-field picture of the leading components for zero total incoming momentum and frequency has been obtained as follows. There is one strongly growing combination of effective interactions that occurs in the flow equation for the anomalous self-energy. It typically reaches its maximum at the critical scale or temperature T c and does not diverge below.  In the BCS mean-field model, this combination is given by V + W for real ∆ [26]. Another combination of the effective interactions, V − W , diverges also below T c . This divergence is associated with the massless phase mode dubbed Goldstone boson present if a continuous symmetry, in our case a global U(1)-symmetry, is broken. For cases where no continuous symmetry is broken as for the half-filled charge-density-wave mean-field model with a two-fold degenerate ground state [27], the cutoff-dependence of the effective interaction is comparable to the cutoff dependence of the combination V + W of normal and anomalous effective interactions of a BCS mean-field model. For the Hubbard interaction studied in this paper, the decay of these large vertices away from zero-total-momentum or -frequency becomes important, as other channels not present in the mean-field models might pick up large contributions from the symmetrybreaking channel. In this section, we investigate V + W and V − W with regard to their momentum and energy dependence as well as to the impact of the 3 + 1 effective interactions.

Momentum dependence
We analyze the momentum dependence of V − W and V + W at ν m = 0. We fit 1 to data obtained from numerical calculations for small up to intermediate p, see figure 5. For the non-divergent part V +W , the fit reproduces the numerically-calculated behavior. For small momenta, the agreement is better, while for larger momenta, deviations become larger, but vanish along the axes. This is due to an anisotropy in the numerical data which is not captured by the approximation (19). For the divergent part V − W , we note that A as obtained from the fitting procedure is proportional to ∆ ext . This is confirmed by calculations for ∆ ext = 10 −3 t where deviations of about 10% are common. The peak at p = 0 is well-captured by the approximation (19), which permits us to confirm the dispersion relation for the low-lying Goldstone mode in the next section.
We consider the impact of the 3+1 effective interactions. By setting loops combining an anomalous and a normal propagator to zero in (15), the 3 + 1 effective interactions are set to zero. The modulus of the relative change of a given effective interaction combination Γ upon switching on the 3 + 1 effective interactions is plotted for Γ = V + W and Γ = V − W in figure 6. We

Energy dependence
The energy dependence of V + W and V − W is to be included in the approximate formula (19) by adding β ν ν 2 m to the denominator, i.e. by considering the approximation function 1 We restrict ourselves to the real part of the effective interactions for brevity.
Contemplating the left-hand graph of figure 7, we become aware of a shoulder close  The blue lines are least-squares fits to (21). The parameters used in the fits are, from left to right, 33t}. Both diagrams: U 0 = 2t, µ = −1.5t, t ′ = t/6, T << T c , ∆ ext = 10 −4 t to zero frequency, which is not captured by (21). This implies a small non-linearity of the dispersion relation of the amplitude mode. The ansatz (21) reproduces the numerically-obtained results for V − W upon properly adjusting the parameters, as the right-hand part of figure 7 illustrates. The agreement is good for small Matsubara frequencies. In the inset, we see that the largefrequency behavior is not captured optimally, which is due to our choice of γ, having taken it from the fit in figure 5. However, the exponent 2 for ν m is confirmed by the fit. This together with the momentum-dependence fit in section 4.5.1 confirms the linearity of the dispersion relation of the Goldstone mode connected to V − W .
To understand the influence of the Ω i (i = 1 . . . 4) on the energy dependence of V + W and V − W , we plot the relative change at p = 0 as in (20), see figure 8. We glean that for V + W , the only appreciable change occurs at ν m = 0 and is quite small in the limit ∆ ext → 0. For V − W , the converse is true: Only at ν m = 0 is there no change, while even for small ν m > 0, a sizable change is observed which is only suppressed at large ν m . This implies that it would not be prudent to neglect the Ω i , while the smallness of their influence on the ν m = p = 0 part of V + W suggests that their quantitative impact on the superconducting gap obtained by an fRG calculation is likely small. The latter supposition is motivated by the observation that V + W drives the flow of the gap, as already shown for the BCS mean-field model in [26].

Functional renormalization group
In this section, we describe the fRG treatment of the attractive Hubbard model on the two-dimensional square lattice at zero temperature. Our main goals are to study the flows produced as well as to obtain reliable quantitative results for the order parameter. Compared to the partial resummation in the last section, the fRG takes into account all one-loop diagrams. The consequential increase in numerical complexity is balanced by adopting a cruder discretization of momentum and frequency space. By taking into account all one-loop diagrams, the pairing channel is renormalized by other processes with non-zero total momentum such as spin-and charge fluctuations. We will find that this causes a reduction of the superconducting order parameter with respect to the mean-field result. Furthermore, the fRG treatment allows for a k-dependence of the self-energy around the Fermi surface.

Momentum-shell and interaction-flow procedures
We introduce a cutoff Λ into our action by substituting We call momentum-shell cutoff function [1] and interaction-flow cutoff function [19]. The former name is reflected in the expansion of the system by an infinitesimal shell in momentum space if Λ → Λ+dΛ. The latter name can be understood by scaling the fields ψ in (5) to 4 √ Λψ, causing the non-interacting part of the action to be cutoff-independent and the interacting, quartic part to be proportional to Λ. For the momentum-shell flow, setting Λ = Λ i > max k |ξ k | provides a solvable, non-interacting starting point. For Λ = Λ f = 0, the original, fully interacting system is recovered. For the interaction flow, Λ = Λ i = 0 suppresses the interaction according to the rescaling of the fields outlined above and provides a solvable starting point. Λ = Λ f = 1 recovers the original, fully interacting system.
Employing the momentum-shell cutoff, we introduce a small external field ∆ ext to select the symmetry-broken phase as the endpoint of the fRG calculation. This means our results in this case apply to a situation with spontaneously-broken symmetry only in the limit ∆ ext → 0. For the interaction-flow cutoff, we add zero in the form ∆ ext − ∆ ext to the quadratic part of the action, including ∆ ext in the non-interacting part and −∆ ext in the interacting part of the action. This approach has been shown to yield symmetrybroken results at zero external field [28] for Λ = Λ f .

fRG equations in the Katanin truncation
The flow equations which we employ are shown in figure 9, where the single-scale propagator S = −G ∂ Λ G −1 0 G is used. The name can be understood by noting that because of the differentiation acting on Θ, S lives only on the energy scale Λ for the momentum-shell flow. These flow equations can be obtained by truncating according to Katanin [25] the infinite hierarchy of one-particle irreducible (1PI) flow equations of [11]. Solving them yields the exact order parameter and effective interaction in the thermodynamic limit for any model with only two-particle interactions which is meanfield exact [27,36,26].
Symbolically, the flow equation for U readṡ (25) where p = p 1 + p 2 , q = p 3 −p 1 , and q ′ = p 3 −p 2 . Note that the last two terms correspond to the second diagram pictured above, as the antisymmetrization is denoted explicitly in the symbolic equation. Prefactors arising from Fourier transformations are absorbed into the summation symbols.
The symbolic form of the flow equation for the self-energy readṡ Here, the single-scale propagator is rewritten into a form better suited for numerical evaluation [27].

Numerical results
In the following, we simplify (25) and (26) and subsequently present results stemming from a numerical quadrature of the simplified flow equations. While we focus on the momentum and frequency structure of the effective interaction in section 4, in this section we are mainly interested in determining the order parameter beyond mean-field theory, taking into account all one-loop diagrams in the fRG. However, this comes at the price of an increase in numerical complexity, which we balance by adopting a cruder discretization of the effective interaction. We first drop the frequency-dependence of U. Consequently, no energy dependence of Σ is generated in the flow according to (26). We can therefore do all Matsubara sums analytically. We furthermore drop the normal part of the self-energy, thus neglecting Fermi surface shifts. All quantities can be considered real in our approximation. A parametrization taking into account partially the momentum dependence of the effective interaction is the so-called "patching" [9,37]. The Brillouin zone is split into patches as shown in figure 10. The functional value of the effective interaction at the mid-angle of the patch and on the Fermi surface is assumed to be also valid for all other momenta in the patch. Thereby, only the momentum dependence along the Fermi surface is taken into account. This is motivated by the finding that the dependence perpendicular to the Fermi surface is irrelevant in the RG sense for the momentum-shell flow without symmetry-breaking ( [5,6]). We specify the momentum arguments of U as triples of patch indices, so V (1, 2, 3) refers to any V (p 1 , p 2 , p 3 ) with p 1 in patch 1, p 2 in patch 2, and p 3 in patch 3. Note that to obtain momenta in the ψ representation from these patch numbers, signs must be adjusted according to (6).

Momentum-shell flows
We replace the step function from (23) by This simplifies the numerical calculation significantly, although one could suspect that the sharp cutoff might eliminate an integration. However, in the Katanin truncation, terms proportional to the derivative of the order parameter arise in addition to the terms proportional to the derivative of the cutoff function. Note that our results do not depend on α in a reasonable range, which adjusts the sharpness of the cutoff.
In figure 11, we plot flows of the order parameter in the first quadrant of the Brillouin zone. We note that the flow saturates for the two external field strengths employed. However, if we employ too small an external field, the effective interaction flows diverge and the method breaks down. Fortunately, the external field can be chosen at least a hundred times smaller than the final value of the superconducting gap. In calculations for mean-field models [27], the modification of the order parameter due to the external field is less than 20% in this case.
From the resummation studies in section 4, we have learned that the Goldstone phase mode V − W dominates for zero momentum transfer in the Nambu formalism, which is equivalent to zero total momentum if we transform back to the original ψ fields according to (6). In figure 12, we plot several flows from this channel. They dominate in magnitude all other effective interaction flows, including those not shown in the figure, and exhibit a spread of ∼ 10% under variation of the momentum combination. They are quite similar to the flows found in [26]. Note that the maximal final value modulus increases stronger than proportionally to ∆ −1 ext . This is a precursor of the divergence to arise for even smaller ∆ ext .
Complementing the phase mode is the effective interaction combination V + W driving the order parameter. Flows of V + W are plotted in figure 13. The spread of the final values increases with decreasing ∆ ext , while the average remains approximately constant. The graphs resemble those found when a discrete symmetry is broken [27]. However, the monotonicity below the critical scale which had been established analytically for the mean-field model is violated by one of the flows in the right-hand figure. Overall, all momentum-shell flows exhibited resemble the flows for the mean-field models.

Interaction flows
In this section, we employ the interaction flow cutoff (24) with an initial anomalous self-energy ∆ s.e.,i = ∆ ext and a counterterm ∆ c = ∆ ext which we include in the bare propagator. In figure 14 superconducting order parameter ∆ = ∆ c − χ∆ s.e. for various counterterms. We note that if we could treat the flow equations without truncation and further approximations, all counterterms should produce the same final self-energy and vertices. Indeed, the counterterm-fRG for a half-filled charge-density-wave mean-field model exhibited a single strong attractor for the flow of the self-energy [28].
For the more general model studied here, the truncation error becomes noticeable, and for certain values of the pairing counterterm, the interaction vertices still diverge. The flows are terminated if the maximal effective interaction exceeds four times the bandwidth, i.e., if the flow leaves the weak-coupling regime, or if Λ reaches Λ f = 1. We see in the left diagram that the flows no longer converge onto a single strong attractor as for the mean-field case [28]. Studying the final 5% of the flow shown in the right plot of figure 14 reveals that there exists a set of small counterterms for which the flows cross and leave the weak-coupling regime, and a minimal counterterm for which the effective interactions remain within the above limit up to Λ f = 1. The flows for counterterms whose strength is just beyond this minimal counterterm terminate closer to each other than do the flows for larger counterterms. The minimal counterterm is approximately twice as strong as the final order parameter value obtained in the minimal counterterm flow. The behavior described above is interpreted as being analogous to the strong-attractor behavior in the mean-field case. The overestimation of the final order parameter value for larger counterterms is attributed to the overestimation of the order parameter in the flow due to the neglect of the order parameter's energy dependence and momentum dependence perpendicular to the Fermi surface. Consequently, we choose the minimal order parameter obtainable by terminating the flow at Λ f in the weak-coupling regime as our approximation for the physical order parameter. Another way to argue for the so chosen counterterm is that the continuous symmetry breaking necessitates a divergent vertex exactly at the end of the flow in order to render the Goldstone mode correctly. This infinite value is suppressed by our approximation, but the flow for the chosen counterterm yields the maximal final value for the Goldstone mode within the above limit. Finally, note that the order parameter obtained by these rules is approximately 20% smaller than the value obtained with the momentum-shell method in section 5.3.1, which is enhanced by a finite, albeit small, external field. figure 15 shows the magnitude of the superconducting order U 0 1.0t 1.5t 2.0t 2.5t 3.0t Figure 15. Magnitude of the superconducting order parameter ∆ as calculated via the interaction flow procedure for varying interaction strength U 0 and chemical potential µ, t ′ = −0.1t parameter ∆ versus the chemical potential as calculated using the interaction flow method including a counterterm as described in the previous section. For small bare coupling U 0 , ∆ is maximal if the van Hove points lie on the bare Fermi surface, which is the case for µ = 4t ′ . For larger bare coupling, the maximum gradually shifts to smaller µ. We compare the results from the fRG to results from mean-field theory in figure 16. We note that the suppression of the mean-field values by the fluctuations taken into account by the fRG varies between a factor 1.6 and 2.0. This is quite similar to values found by second-order perturbation theory [33]. Table 1 compares our results to results from a combination of fRG and mean-field theory [38] for finite t ′ . While the values are of the same order of magnitude, the results obtained by employing only the fRG are larger by between 20 − 30%.

Conclusions
We have applied the fermionic one-particle irreducible (1PI) version of the functional renormalization group (fRG) in the Katanin truncation to the symmetry-broken phase of a two-dimensional, non-mean-field model exhibiting superconductivity at zero temperature. The attractive Hubbard model has been chosen as the prototype system because it is sufficiently well-studied to provide a good testbed for new many-body methods. In this case, particle-hole as well as particle-particle loops are important for the symmetry-breaking properties of the system. The fRG takes both into account on an equal footing. Our analysis has proceeded in three steps. First, the usual Nambu formalism has been extended to cover 3 + 1 anomalous effective interactions featuring an odd number of incoming legs, whose behavior and importance was a focus of our interest. We have found that for the case of superconductivity, the effective interaction had to be extended into a 4 × 4 matrix containing certain redundancies due to symmetries. This extension followed from a substitution of the fields in which the model is originally formulated by Nambu fields which are chosen specifically to produce an action of canonical shape while encompassing singlet Cooper pairs in the quadratic part. 3+1 effective interactions were found to arise naturally within this formalism.
In the second step, a study of the system has been performed by means of a resummation of a subset of all perturbation theory diagrams. The strength of the symmetry breaking was determined by a mean-field gap equation. A class of diagrams containing bubble chains, ladders, and all chains comprising both was resummed into a Bethe-Salpeter equation. In the context of this equation, the impact of the 3 + 1 effective interactions was studied by comparing effective interaction values calculated with and without the 3 + 1 effective interactions. This impact was found to be zero for the zero-frequency part of the phase mode, and to be dependent on the external field for non-zero frequency for both phase and amplitude mode.
In the third step, the fRG has been numerically applied to the attractive Hubbard model. Two distinct variants of the fRG were employed. First, a momentum-shell cutoff was used together with an external field providing a finite initial order parameter. The Goldstone mode divergence is regularized by the field. Despite the approximate nature of the flow and the explosive character of the fRG differential equation, flows reaching down to zero scale were found employing external field strengths two orders of magnitude smaller than the final order parameter values. The relative error of the results for the order parameter due to such an external field was estimated to be at most 20%, judging from earlier calculations for mean-field models. For even smaller external field values, a regime of diverging flows was found. Second, we have calculated interaction flows, balancing an external field with a counterterm increasing continuously with the flow. We have noted that vestiges of strong-attractor behavior found previously for a mean-field model are also present for the attractive Hubbard model flows. They permit the determination of order parameter values which are comparable to results from the literature. While these values proved to be smaller than values calculated by mean-field theory alone, as expected, they are generally slightly larger than results obtained by combining symmetric-phase Wick-ordered fRG calculations and mean-field theory. This deviation is likely due to small differences between the 1PI and Wickordered fRG schemes, since symmetric-phase 1PI calculations already yield larger critical scales than the corresponding Wick-ordered calculations. In summary, the 1PI fRG in the Katanin truncation was found to yield fairly accurate quantitative results for weak and moderate interaction strengths (up to U 0 /t = 3). We reckon that the diverging flows present for very small external fields in our calculations can be pushed back by improving the resolution of the discretization in both frequency and momentum space. This would furthermore improve the accuracy of the results. The method may also be employed to study symmetry breaking in interacting-fermion models for which a stronger momentum dependence of the order parameter is expected, for example to study the repulsive Hubbard model. The method can also be applied to situations with competing order. The interaction-flow scheme with possible symmetry-breaking in two different channels has recently been tested in a one-dimensional model with competing long-range correlations. There it allowed for a correct determination of a quantum critical point between a charge-gapped phase and a compressible phase with dominant superconducting correlations [39]. Finally, it will also be interesting to analyze how non-zero temperature affects the renormalization of the mean-field-like instability encountered in the flows.