Dynamical functional renormalization group computation of order parameters and critical temperatures in the two-dimensional Hubbard model

We analyze the interplay of antiferromagnetism and pairing in the two dimensional Hubbard model with a moderate repulsive interaction. Coupled charge, magnetic and pairing fluctuations above the energy scale of spontaneous symmetry breaking are treated by a functional renormalization group flow, while the formation of gaps and order below that scale is treated in mean-field theory. The full frequency dependences of interaction vertices and gap functions is taken into account. We compute the magnetic and pairing gap functions as a function of doping $p$ and compare with results from a static approximation. In spite of strong frequency dependences of the effective interactions and of the pairing gap, important physical results from previous static functional renormalization group calculations are confirmed. In particular, there is a sizable doping regime with robust pairing coexisting with N\'eel or incommensurate antiferromagnetism. The critical temperature for magnetic order is interpreted as pseudogap crossover temperature. Computing the Kosterlitz-Thouless temperature from the superfluid phase stiffness, we obtain a superconducting dome in the $(p,T)$ phase diagram centered around 15 percent hole doping.


I. INTRODUCTION
Shortly after the discovery of high-temperature superconductivity in cuprates, Anderson 1 proposed the twodimensional Hubbard model to describe the behavior of the valence electrons in the copper-oxygen planes. Indeed, the model captures the most prominent ordered phases observed in high-T c cuprates: antiferromagnetism and d-wave superconductivity. 2 Antiferromagnetism in the Hubbard model is not always of Néel type, that is, with antiparallel spin orientation between adjacent lattice sites. Magnetic order with (generally incommensurate) wave vectors away from the Néel point (π, π) has been found away from half-filling in several mean-field studies, [3][4][5][6][7] and also, including fluctuations, by expansions in the hole-density. [8][9][10][11] Unbiased evidence for superconductivity with a sizable energy scale already at moderate interaction strengths has been established from functional renormalization group (fRG) calculations, [12][13][14][15][16][17] and from quantum cluster methods at intermediate to strong coupling. [18][19][20][21][22][23][24][25] Recently, a fRG flow starting from the dynamical mean-field solution (instead of the bare action) has confirmed robust pairing with d-wave symmetry at strong coupling. 26 The fRG flow is defined by a successive, scaledependent integration of the fermionic fields in a path integral representation of the effective action. 17 Spontaneous symmetry breaking is signaled by a divergence of effective interactions at a critical energy scale Λ c . Accessing the ordered phase by continuing the flow beyond the critical scale is possible, 27 but rather complicated due to a proliferation of anomalous interaction terms. 28,29 One option is to decouple the fermionic interaction by introducing a bosonic order parameter field via a Hubbard-Stratonovich transformation, and to study a coupled flow involving fermions and bosons. 30 Alternatively, as a "poor man's" approximation, one may restrict the effective interactions below the critical scale Λ c to those (reduced) interactions which drive the symmetry breaking. 31,32 The flow of these interactions is governed only by a single channel (one for each order parameter), and the final order parameters at the end of the flow are given by simple gap equations, which can be derived from the fermionic flow equation 32 or by introducing a bosonic order parameter field at the critical scale. 33 Consistently formulated, this procedure is equivalent to a mean-field theory (MFT) for the degrees of freedom below the critical scale. In a regime where pairing is the only instability, it turned out that the ground state pairing gap of the Hubbard model obtained from this fRG + MFT approach agrees well with results from the much more involved flow equations with coupled interaction channels. The method was then extended to parameter regions where pairing coexists with Néel 32 or spiral 34 antiferromagnetism. Coexistence of antiferromagnetism and d-wave pairing was also found in quantum cluster calculations at stronger interactions. 19,[21][22][23]25 So far, fRG calculations of order parameters for the repulsive two-dimensional Hubbard model have relied on a static approximation, that is, the frequency dependences of the effective interaction and gap functions have been neglected. 35 However, dynamical fRG flows in the symmetric regime (before reaching the critical scale Λ c ) have indicated a rather strong impact of the frequency dependence already at moderate interaction strengths. 36,37 The effective two-particle interactions generally develop a strong dependence on all three Matsubara frequencies, and the critical scale is enhanced compared to the static approximation. More importantly, it turned out that the frequency dependence leads to an expansion of the parameter regime where antiferromagnetism is the first instability (at the critical scale Λ c ), with rather weak pairing interactions at that scale. This cast some doubt on the robust pairing tendencies obtained in the static fRG.
In the present work we address this issue by using a dynamical extension of the fRG + MFT method with full arXiv:2010.06267v1 [cond-mat.str-el] 13 Oct 2020 frequency dependence to compute magnetic and pairing gap functions for the repulsive two-dimensional Hubbard model with a moderate interaction strength. While magnetism is indeed the leading instability in a broad doping range, we find that robust pairing with a sizable pairing gap emerges in coexistence with antiferromagnetism at energy scales below Λ c around optimal doping. The size of the gap is not reduced compared to the one obtained in a static fRG + MFT approximation. We also compute the superfluid phase stiffness which enables us to estimate the Kosterlitz-Thouless transition temperature for superconductivity as a function of doping.
Our paper is structured as follows. In Sec. II we define the model and we present the fRG flow equations. Results for the gap functions as a function of doping and frequency, for the phase stiffness, and the Kosterlitz-Thouless temperature are presented and discussed in Sec. III. A summary and conclusion in Sec. IV closes the article.

A. Model
The Hubbard model 38 describes spin-1 2 lattice fermions with quantum mechanical hopping amplitudes and a local interaction. Its Hamiltonian has the form where c † i,σ (c i,σ ) creates (annihilates) a fermion on site i with spin orientation σ (↑ or ↓). We consider the twodimensional case on a square lattice with a repulsive interaction U > 0 at finite temperatures T . We restrict the hopping matrix t ij to nearest and next-to-nearest neighbors, with amplitudes −t and −t , respectively. Fourier transforming t ij yields the dispersion relation k = −2t (cos k x + cos k y ) − 4t cos k x cos k y . (2)

B. Functional renormalization group
We compute the gap functions for magnetism and superconductivity from a truncated fRG flow. The flow is defined by a successive integration of the fermionic fields, which is implemented via a flowing cutoff applied to the bare propagator. 17 Here we choose a smooth frequency cutoff of the form 39 where G 0 (k, ν) = [iν − ( k − µ)] −1 is the bare fermion propagator as a function of the crystal momentum k and the fermionic Matsubara frequency ν (odd integer multiples of πT ). The flow parameter Λ is reduced continuously from infinity to zero. At low temperatures the flow needs to be divided in two qualitatively distinct regimes. For Λ > Λ c all symmetries of the bare Hamiltonian are conserved, while for Λ < Λ c the SU(2) spin symmetry, the U(1) charge symmetry, or both are spontaneously broken. The instabilities are signalled by divergencies of effective interactions at the critical scale Λ c . The latter is non-zero for temperatures below a pseudocritical temperature T * , which depends on the hopping amplitudes, the interaction strength, and the band filling.
In the symmetric regime we approximate the flow by a second order (one-loop) truncation for the effective two-particle interaction, discarding self-energy feedback and contributions from the three-particle interaction. All fluctuation channels (charge, magnetic, pairing) and the coupling between these channels are taken into account on equal footing. Fluctuation driven instabilities such as d-wave pairing from magnetic fluctuations are captured by this weak-coupling truncation. In the symmetrybroken regime (Λ < Λ c ) we further simplify the flow by keeping only those interactions which generate the order parameters. The flow of these reduced interactions is decoupled, that is, each of them is determined by a single channel only. The flowing reduced interactions determine the flow of the gap functions, which appear as anomalous self-energy contributions. Their feedback on the flow of the interactions is crucial in the symmetrybroken regime 27 and is therefore taken into account. The fusion of a complete one-loop flow for Λ > Λ c with a single-channel truncation for Λ < Λ c corresponds to a mean-field approximation with effective interactions extracted from the flow at the critical scale Λ c . 32

C. Symmetric regime
The two-particle vertex V Λ σ1σ2σ3σ4 (k 1 , k 2 , k 3 , k 4 ) with k i = (k i , ν i ) is generally a function of four momentum, frequency, and spin variables, where the labels 1 and 2 correspond to ingoing, the labels 3 and 4 to outgoing particles. Translation invariance in space and time implies momentum and frequency conservation, k 1 +k 2 = k 3 +k 4 . We can thus drop the redundant variable k 4 in V Λ . In case of SU(2) spin rotation invariance, all spin components of the vertex can be expressed by a single momentum and frequency dependent function To parametrize the momentum dependence of the vertex, we use the channel decomposition introduced by Husemann and Salmhofer, 39 where the vertex is expressed as a sum of the bare interaction and fluctuation induced effective interactions in the pairing, magnetic, and charge channels, The dependence on the two fermionic momenta in each channel is more regular than the dependence on the bosonic linear combination in the last argument. We will neglect it completely in the magnetic and charge channel, where only relatively weak momentum dependences are expected, and approximate it by a constant and a d-wave term in the pairing channel, that is, with q = (q, ω) and d k = cos k x −cos k y . The dependence on the remaining bosonic momentum q is kept, and the dependence on all three Matsubara frequencies ν 1 , ν 2 , and ω is fully taken into account. The parametrization of V Λ (k 1 , k 2 , k 3 ) via Eqs. (5) -(8) has already been used in Ref. 37. The flow equations for the four functions S Λ , D Λ , M Λ , and C Λ can be found there.

D. Symmetry-broken regime
In the regime Λ < Λ c at least one symmetry of the Hubbard Hamiltonian is spontaneously broken. The flow of the effective interactions 17 and most studies by other methods 2 indicate antiferromagnetism (Néel, stripes, spiral etc.) and d-wave pairing as the key instabilities. We restrict the zoo of possible magnetic order patterns to spiral order with a single wave vector Q, which includes Néel order as the special case where Q = (π, π). Spiral order is planar and can be oriented in any plane. The most convenient choice is a plane perpendicular to the spin-quantization axis, that is, the xy-plane for the standard spin basis. Spiral magnetic order in the xy-plane is associated with anomalous expection values ψ ↑ (k)ψ * ↓ (k + Q) , where Q = (Q, 0), and singlet pairing with anomalous expectation values of the form ψ ↑ (k)ψ ↓ (−k) . Here and in the following ψ * σ (k) and ψ σ (k) are Grassmann fields corresponding to fermion creation and annihilation operators in momentum representation, respectively. Symmetry breaking leads to anomalous terms in the effective action Γ Λ [ψ, ψ * ]. Spiral order leads to spin-flips combined with a momentum shift Q, and pairing leads to pair creation and annihilation terms. The quadratic part of the effective action thus has the general form where k * = (k, −ν), (2π) 2 is a shorthand notation for momentum integrals and Matsubara frequency sums. Σ Λ (k) is the normal self-energy, which we neglect in this work. ∆ Λ m (k) and ∆ Λ p (k) are (generally) complex functions, which we refer to as magnetic gap function and pairing gap function, respectively. The frequency dependence of the (spin-singlet) pairing gap is symmetric, that is, with a 4-component Nambu spinor and a 4 × 4 Nambu propagator The prime at the integral indicates that the momentum integration is restricted to a reduced magnetic Brillouin zone. Its shape depends on Q. For example, for Q = (π − 2πη, π), a suitable reduced Brillouin zone is given by {k; |k x | ≤ π, |k y | ≤ π/2}. The matrix elements of the inverse Nambu propagator are determined by Eq. (9) as Our central approximation in the regime Λ < Λ c is that we keep only those effective interactions which contribute directly to the flow of the gap functions. The corresponding reduced quartic part of the effective action has the form The coupling functions V Λ m (k, k ) and V Λ p (k, k ) parametrize reduced normal interactions.
For example, V Λ p (k, k ) corresponds to the coupling function of the reduced BCS model, which is restricted to the Cooper channel (with vanishing total momentum of ingoing and outgoing particles). W Λ m (k, k ) and W Λ p (k, k ) parametrize anomalous interaction terms which are generated in the symmetry broken regime. 41 At the critical scale Λ c (and above it), the anomalous terms vanish, while the normal reduced coupling functions can be obtained from the full two-particle vertex is real, since the imaginary parts cancel in the spin-singlet component.
Due to the restrictions of momenta in the reduced interactions, the flows of the magnetic and pairing coupling functions are decoupled from each other, and each of them is governed by one channel only. For example, the flow of the pairing coupling is determined by the particleparticle channel. The flow of the gap functions is entirely determined by the amplitude coupling 27,28,42 with X = m, p, so that we consider only this linear combination in the following. The transverse coupling is related to a Ward identity and the Goldstone theorem, which are fulfilled in the fRG + MFT approach. 33 In line with our parametrization in the symmetric regime, we simplify the momentum dependencies of the gap functions and the coupling functions by using a small set of form factors. In the magnetic channel, the dependences on k and k are weak and will be neglected.
In the pairing channel we keep only the d-wave components, since there is no pairing instability with any other symmetry. Hence, we approximate the momentum dependences of the gap and coupling functions by a simple ansatz, namely for the gap functions, and for the coupling functions. The dependences on the Matsubara frequencies ν and ν are fully taken into account. At the critical scale Λ c , we have ∆ Λc m (ν) = ∆ Λc p (ν) = 0, and These are the initial conditions for the flow in the symmetry-broken regime Λ < Λ c . In our numerical solution we will switch from the full one-loop flow in the symmetric regime to the reduced single-channel flow slightly above the critical scale Λ c , and we insert tiny gap values as initial condition for the gap functions to get the symmetry breaking started.
Since only a single channel contributes to the flow of each coupling function A Λ X , the right hand side of the flow equation is a quadratic form in A Λ X . It is convenient to view A Λ X as a matrix with matrix elements A Λ X (ν, ν ). The flow equation for A Λ X can then be written in matrix form, The propagators in Eqs. (24) and (25) are defined by the expectation values where G Λ αα (k) are matrix elements of the Nambu propagator defined in Eq. (12).
The flow equation (23) can be formally integrated. With the initial condition A Λ X = A Λc X = V Λc X for Λ = Λ c , one obtains the solution where [. . . ] −1 denotes a matrix inversion, and V Λc X and V Λc X are related by the Bethe-Salpeter equation Hence,Ṽ Λc X is the two-particle irreducible part of V Λc X . The flow equations for the gap functions are obtained from the first equation in the flow equation hierarchy, which relates the flow of the self-energy to the twoparticle vertex. 17 Inserting the two-particle vertex as described above, one finds where F Λ m (k) and F Λ p (k) are the anomalous propagators corresponding to magnetic order and pairing, respectively. The derivative denoted by ∂ Λ acts only on the scale dependence directly introduced by the cutoff function, not on the scale dependence from the flowing gap functions appearing in the expressions for F Λ m (k) and F Λ p (k). By contrast, the scale derivative ∂ Λ in Eq. (23) is a total derivative, which includes self-energy feedback terms generated from tadpole contractions of threeparticle vertices. 27,44 The flow equations (31) and (32) can be formally integrated to 32 whereṼ Λc X is the irreducible part of V Λc X . These are nonlinear integral equations for the gap functions ∆ Λ X (ν). Finding a self-consistent solution by iteration is difficult. We found that it is much easier to compute the gap functions from a numerical integration of the flow equations (31) and (32).
The pairing gap function is symmetric in frequency, that is, ∆ p (ν) = ∆ p (−ν), and can be chosen real for all frequencies, sinceṼ Λc p (ν, ν ) is real for all frequencies. The magnetic gap function can be chosen such that ∆ m (−ν) = ∆ * m (ν). Its imaginary part cannot be removed by choosing a suitable global phase, since the effective interactionṼ Λc m (ν, ν ) is complex.

III. RESULTS
Before presenting our results, we mention here a few technical details regarding the numerical solution of the equations. In the symmetric regime, we solve the flow equation of the competing channels by taking into account about 90 Matsubara frequencies for each frequency argument and about 320 Brillouin zone patches for the transfer momentum. When required, we extend the frequency range with a numerical projection. 45 We have checked that the frequency boxes are large enough to obtain converged results. The symmetric one-loop flow is stopped when the effective interaction in one of the channels reaches the value 400t, that is, at a scale Λ very close to the critical scale Λ c . In the symmetry-broken regime, the two order parameters are initialized with a small value of the order 10 −3 t. The irreducible vertices, Eq. (30), and the flowing vertices, Eq. (29), in the symmetry-broken regime are calculated with a matrix inversion in Matsubara frequency space. The momentum Q characterizing the spiral order is determined from the maximum of the magnetic interaction upon approaching the critical scale Λ c in the symmetric regime. It has the general form Q = (π − 2πη, π), or symmetry-related. We have checked that the peak momentum of the magnetic interaction at Λ c is very close to the momentum Q which minimizes the free energy at the end of the flow.
We choose a fixed hopping amplitude ratio t /t = −0.16 and a moderate interaction strength U = 3t in all calculations. We use natural units such that = k B = 1. The lattice constant is also set to one. Quantities with dimension energy (= temperature) are presented in units of t.

A. Order parameters
In Fig. 1, we show the amplitudes (maxima in momentum space) of the magnetic and the pairing gaps as a function of the hole doping p = 1 − n at a fixed temperature T = 0.027t. Results from the dynamical fRG with full frequency dependence as described above are compared to results from a static approximation as employed by Wang et al. 32 and Yamase et al., 34 where all frequency dependences are neglected. The gap amplitudes obtained from the dynamical fRG are shown at the lowest positive Matsubara frequency, that is, ν 0 = πT . The pairing gap amplitude is 2∆ p (ν), since the maximum value of the dwave form factor d k = cos k x − cos k y is two. Within our ansatz, the magnetic gap is momentum independent. Its "amplitude" is thus simply ∆ m (ν). It is generally complex, but the imaginary part is small for ν = ν 0 = πT . The imaginary part of the analytic continuation of ∆ m (ν) to the entire complex frequency plane vanishes for ν → 0. In the figure we only show the real part of ∆ m (ν 0 ). In the static approximation the magnetic gap is frequency independent and can be chosen real.
Antiferromagnetic order extends from half-filling to a doping value of around 20%, with Néel order up to about 13%, and incommensurate spiral order beyond. The transition between Néel order and incommnsurate spiral order is discontinuous, with a pronounced jump of the incommensurability η. In the range between 8% and 20% doping, a sizable pairing order parameter appears. Hence, our calculation confirms the coexistence of antiferromagnetism and d-wave pairing in the two-dimensional Hubbard model, as previously obtained from static fRG calculations. 31,32,34 . The doping dependence of the incommensurability η also agrees well with previous static fRG calculations. 34 The pairing mechanism is magnetic, that is, the attraction in the d-wave pairing channel is predominantly generated by antiferromagnetic fluctuations. The suppression of the pairing amplitude close to half-filling is caused by the magnetic gap leading to a truncation of the Fermi surface to small hole pockets. The gap amplitudes obtained from the static fRG are somewhat smaller than those obtained from the dynamical fRG, but they exhibit a similar qualitative doping dependence. The magnetic order in the static approximation is weaker for two reasons, which have been revealed already previously. 37 First, there is a minimum of the magnetic effective interaction M Λ Q,ω=0 (ν 1 , ν 2 ) at the lowest fermionic Matsubara frequencies |ν i | = πT . In the static approximation this minimum value is practically extended to all frequencies ν 1 and ν 2 . Second, in the static approximation the suppression of magnetic interactions from other fluctuation channels is overestimated. 37 The pairing gap is also reduced in the static approximation. Since the pairing mechanism is mostly magnetic, one reason for this reduction is certainly the weaker magnetic interaction. On the other hand, the decay of the effective pairing interaction at large frequencies is neglected in the static approximation, leading to an enhancement of pairing tendencies. 36,37 The net effect The influence of the frequency dependence of the magnetic interaction vertex on pairing was previously analyzed by Kitatani et al. 46 at an intermediate coupling strength (U = 6t) within the dynamical vertex approximation. 47 In agreement with our results, they found a minimum at low frequencies. They concluded that this minimum leads to a significant reduction of the energy scale for pairing. This is not in conflict with our results, because they compared their result to that from a simple random phase approximation for the magnetic interaction, which grossly overestimates its strength in the relevant frequency range.
Yamase et al. 34 observed a pronounced dip of the magnetic gap at van Hove filling, where the pairing gap suppresses the magnetic order completely. This feature is not visible in Fig. 1, neither in the static nor in the dynamical approximation. The discrepency is probably due to the finite temperature in our calculation. At present, we can access lower temperatures only by neglecting frequency dependences. In Fig. 2 we show the doping dependence of the gap amplitudes at T = 0.006t as obtained from the static fRG. Here, a dip of the magnetic gap at van Hove filling (p vH ≈ 14% for our parameters) is clearly visible. We expect that the dip will become even more pronounced upon further lowering the temperature.

B. Flow and frequency dependence
In this section we present some details on the flow and on the frequency dependences for a fixed doping p = 0.12 and a fixed temperature T = 0.027t. In Fig. 3 we show the flow of the magnetic gap ∆ Λ m (ν 0 ) at the lowest posi- tive Matsubara frequency ν 0 = πT together with the flow of the magnetic amplitude coupling A Λ m (ν 0 , ν 0 ). The flow looks qualitatively the same as for the BCS model. 27 The effective interaction increases rapidly upon approaching the critical scale Λ c from above and decreases again below it, approaching eventually a moderate finite value. The peak at Λ c is regularized by the small initial gap term inserted by hand. The gap increases monotonically from its tiny initial value at Λ c to a much larger value for Λ → 0. The flows of ∆ Λ m (ν) and A Λ m (ν, ν ) have the same qualitative behavior for all choices of the Matsubara frequencies ν and ν . The flow of the pairing gap and pairing coupling looks similar, too, but the onset of the gap and the peak in the coupling is situated at a lower scale.
In Fig. 4 we show the frequency dependence of the gap amplitudes at the end of the flow. The pairing gap can be chosen real for all frequencies, while the magnetic gap is necessarily complex. It is obvious that the frequency dependence of Re∆ m (ν) is very weak, while ∆ p (ν) decays rather quickly over the first few Matsubara frequencies.
Hence, an accurate continuation of the pairing gap to zero frequency is difficult.
The frequency dependence of effective magnetic and pairing interactions is shown in Fig. 5. There is a pronounced structure on the diagonal ν = ν in the magnetic channel, and for ν = ±ν in the pairing channel. Note that these features are not present in the coupling functions M Λ q,ω (ν, ν ) and D Λ q,ω (ν, ν ), respectively, but are rather due to peaks of the other coupling functions as a function of the bosonic frequency ω at ω = 0. For the magnetic interaction, a qualitatively similar frequency structure is also obtained from a dynamical mean-field approximation, which neglects all correlations except the local ones. 48 The frequency dependence is obviously generated in the symmetric regime, so that it is fully developed already at the critical scale Λ c . The amplitude coulings A Λ X (ν, ν ) at the end of the flow (Λ → 0) exhibit the same frequency dependence as the two-particle irreducible interaction partsṼ Λc X (ν, ν ) at the critical scale Λ c . Hence, the frequency structure is not significantly changed in the symmetry-broken regime. In the pairing channel there is even quantitative agreement betweeñ V Λc p (ν, ν ) and A Λ p (ν, ν ). This indicates thatṼ Λc p Π Λ p becomes small for Λ → 0, see Eq. (29), which must be due to a strong suppression of Π Λ p by the gap formation. For large frequencies and away from the special lines ν = ±ν, the reduced magnetic interactionṼ Λc m (ν, ν ) tends to the bare Hubbard coupling U , while the reduced pairing interactionṼ Λc p (ν, ν ) decays to zero. The latter behavior is the reason for the decay of the pairing gap at large frequencies described above.

C. Superfluid stiffness and phase diagram
We finally compute the superfluid phase stiffness, which allows us to estimate the Kosterlitz-Thouless temperature T KT for the onset of superconductivity. Together with the temperature T * for the onset of antiferromagnetism, we can thus draw a phase diagram in the plane spanned by doping and temperature.
A general expression for the phase stiffness in a meanfield state with coexisting spin-singlet superconductivity and antiferromagnetism (Néel or spiral) has been derived in a recent work by Yamase and one of us. 49 The phase stiffness is fully determined by the bare dispersion relation and the magnetic and pairing gaps. The gaps have been assumed to be frequency independent in the derivation. We therefore neglect the frequency dependence of the gaps, and insert the gap at the lowest Matsubara frequency ν 0 = πT . In a spiral state with an ordering vector Q = (π − 2πη, π) with η > 0 the phase stiffness in x-and y-direction is slightly different. In Fig. 6 we plot the phase stiffnesses J x and J y as a function of doping at the fixed temperature T = 0.027t. The pairing gap amplitude from Fig. 1 is also reproduced for direct comparison. In the Néel state the stiffness is isotropic, J x = J y , while in the spiral regime J y is slightly smaller than J x . The stiffness and the gap amplitude have a comparable size in the regime where both are larger than the temperature, but the stiffness decreases much faster at low doping, since it is more strongly suppressed by thermal excitations than the gap. By contrast, a static fRG calculation at zero temperature indicated that in the ground state the stiffness decreases slightly more slowly than the gap amplitude upon approaching half-filling from the holedoped side. 49 Note that, in general, there is no direct relation between the size of the gap and the size of the phase stiffness in a superconductor.
In a two-dimensional system, the thermal phase transition between the superfluid and the normal phase is a Kosterlitz-Thouless transition associated with topological excitations (vortices). 50 Magnetic order or (noncritical) magnetic fluctuations do not affect the universal properties of this transition. In an isotropic system, the transition temperature T KT is related to the phase stiffness J by the universal relation T KT = π 2 J(T KT ). 50 In the spiral state the phase stiffness is slightly anisotropic, that is, we have J x = J y , while there is still a unique transition temperature. Generalizing the relation between T KT and J to anisotropic systems by a simple rescaling of the length scales in the phase action, we find Using this relation we are able to compute the Kosterlitz-Thouless temperature from the stiffnesses J α (T ), as long as T KT is higher than the lowest temperature T = 0.027t we can access. In Fig. 7 we plot the resulting Kosterlitz-Thouless temperature as a function of doping, together with the critical temperature T * for the onset of antiferromagnetism. The latter is determined as the lowest temperature at which the fRG flow does not encounter any magnetic instability down to Λ = 0. We also show the pairing temperature T p at which the pairing gap ∆ p vanishes. One can see that T p is much higher than T KT , especially at lower doping. Hence, a sizable temperature window with a pairing gap and superconducting fluctuations opens between T KT and T p . Our mean-field approximation yields magnetic longrange order for temperatures below T * . General arguments and numerical studies show that fluctuations destroy this long-range order, giving rise to a pseudogap state with strong short-ranged magnetic correlations. Hence, T * should be interpreted as the onset temperature for pseudogap behavior. Implementing the fluctuations that turn the magnetically ordered state into a pseudogap state will be an interesting extension of our present theory. Since we cannot access temperatures below T = 0.027t with our present dynamical fRG code, we cannot calculate the Kosterlitz-Thouless temperature below 10 percent doping. Low temperature results from a static fRG flow indicate that T KT vanishes linearly in doping upon approaching half-filling. 49

IV. CONCLUSION
We have performed a dynamical fRG analysis of magnetic order and superconductivity in the two-dimensional repulsive Hubbard model at a moderate interaction strength U = 3t. A one-loop flow with coupled charge, magnetic and pairing interaction channels in the symmetric regime above the critical energy scale Λ c was combined with a mean-field approximation with decoupled reduced interactions in the symmetry-broken regime below Λ c . The full frequency dependences of interaction vertices and gap functions were taken into account. The momentum dependences were approximated by suitable form factors.
While magnetism appears as the leading instability at the critical scale Λ c in a broad doping range from halffilling to about 20 percent, robust pairing with a sizable pairing gap emerges between 10 and 20 percent doping, in coexistence with antiferromagnetism. The size of the pairing gap is slightly enhanced compared to results from a static fRG, probably as a consequence of the enhanced magnetic interactions. The effective interactions exhibit strong frequency dependences, as already observed in previous dynamical fRG calculation, which were however limited to the symmetric regime. 36,37 The magnetic gap depends only weakly on frequency, while the pairing gap is peaked at the lowest Matsubara frequency and decays rapidly at higher frequencies.
We have also computed the superfluid phase stiffness and the Kosterlitz-Thouless transition temperature T KT as a function of doping p. Combining T KT (p) with the onset temperature T * (p) for magnetism yields a phase diagram with a dome-shaped superconducting regime under the quasi-parabolic roof defined by T * (p). Including magnetic order parameter fluctuations would replace the ordered antiferromagnet in our theory by a pseudogap state with strong magnetic correlations at any finite temperature. The shape of T * (p) and T KT (p) agrees qualitatively with the pseudogap temperature and the superconducting transition temperature, respectively, in high-T c cuprates.
Our work can be naturally extended in two directions. First, the fRG + MFT approach can be extended to the strongly interacting regime by using the dynamical meanfield solution 51,52 of the Hubbard model as a starting point for the fRG flow. The combination of dynamical mean-field theory and fRG was proposed some time ago, 53 and recently applied to the two-dimensional Hubbard model at strong coupling, albeit only in the symmetric regime up to the point where effective interactions diverge. 26 Magnetic and pairing gap functions could now be computed by continuing the flow with reduced but dynamical interactions into the symmetry broken regime.
Second, the mean-field solution for magnetic order in the symmetry-broken regime could be improved by implementing thermal and quantum fluctuations of the spin orientation. A most promising route to do this is to treat the magnetic regime with an SU(2) gauge theory as recently developed by Scheurer et al. 54 The critical temperature T * for magnetism obtained from the fRG flow assumes the role of the pseudogap temperature in that improved theory.