{'fraction_non_alphanumeric': 0.04623184323824086, 'fraction_numerical': 0.011346718705991228, 'mean_word_length': 4.203070481507328, 'pattern_counts': {'":': 0, '<': 8, '<?xml version=': 0, '>': 7, 'https://': 0, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 8, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
The low-energy phase-only action in a superconductor: a comparison with the XY model

The derivation of the effective theory for the phase degrees of freedom in a superconductor is still, to some extent, an open issue. It is commonly assumed that the classical XY model and its quantum generalizations can be exploited as effective phase-only models. In the quantum regime, however, this assumption leads to spurious results, such as the violation of the Galilean invariance in the continuum model. Starting from a general microscopic model, in this paper we explicitly derive the effective low-energy theory for the phase, up to fourth-order terms. This expansion allows us to properly take into account dynamic effects beyond the Gaussian level, both in the continuum and in the lattice model. After evaluating the one-loop correction to the superfluid density we critically discuss the qualitative and quantitative differences between the results obtained within the quantum XY model and within the correct low-energy theory, both in the case of s-wave and d-wave symmetry of the superconducting order parameter. Specifically, we find dynamic anharmonic vertices, which are absent in the quantum XY model, and are crucial to restore Galilean invariance in the continuum model. As far as the more realistic lattice model is concerned, in the weak-to-intermediate-coupling regime we find that the phase-fluctuation effects are quantitatively reduced with respect to the XY model. On the other hand, in the strong-coupling regime we show that the correspondence between the microscopically derived action and the quantum XY model is recovered, except for the low-density regime.


I. INTRODUCTION
The description of the low-energy dynamics of the phase θ of the superconducting order parameter ∆ = |∆|e iθ has recently attracted a renewed interest in connection to the phenomenology of high-T c superconducting cuprates. These materials are characterized by a strongly anisotropic (quasi-two-dimensional) crystal structure and by anomalous normal-state properties, which can be interpreted in terms of preformed Cooper pairs (non-zero amplitude |∆|) without phase coherence (vanishing superfluid density ρ s ). 1 Although in ordinary (weak-coupling) superconductors both |∆| and ρ s vanish at the critical temperature T c , due to thermally excited quasiparticles, it is possible that ρ s vanishes when |∆| is finite, due to phase fluctuations. If such is the case, phase fluctuations should play a role more relevant than expected in ordinary superconductors, both near T c and deeply in the superconducting state.
A crucial point in the analyses of phase-fluctuation effects is clearly the choice of the effective model, the most common being the classical XY model, as the paradigm for the universality class which is relevant to superconductivity. 2 The XY model, which has also been extensively used in contexts such as systems of resistively shunted Josephson junctions or granular superconductors, 3 assumes a phase field θ i defined on a coarse-grained lattice with a lattice constant determined by the coherence length ξ 0 , which sets the length scale above which the fluctuations of the amplitude |∆| become uncorrelated, and |∆| may be assumed to be fixed. Thus, at distances larger than ξ 0 , phase fluctuations only are relevant, and can be described in terms of a Josephson-like interaction between neighboring sites < i, j >, Here θ ij ≡ θ i − θ j ≃ | ∇θ|ξ 0 is the phase difference between nearest-neighboring sites, and D 0 is the bare coupling constant. Eq. (1) is written for a two-dimensional system (which is the case relevant for cuprate superconductors). In d = 2, D 0 is an energy which measures the superfluid stiffness, i.e., the energetic cost to produce phase variations in the system. In the continuum (Galilean invariant) model D 0 = ρ s (T = 0)/m = ρ/m, where m is the electron mass, ρ is the electron density, and we set = 1. In the case of a generic d-dimensional system, we still refer to D 0 as to the stiffness, even though the corresponding energy scale is D 0 ξ d−2

0
. Within the XY model, one can determine the depletion of the superfluid stiffness due to the collective-mode excitations, i.e., to the anharmonic terms in the gradient expansion of Eq. (1). At low temperature, the leading correction to the bare stiffness D 0 comes from the fourth-order term, and within perturbation theory with respect to   The specific issue of the temperature dependence of the superfluid density in high-T c superconductors has been analyzed in Ref. [5,6]. In general, according to Eq. (1), one would expect that Coulomb effects, lifting the phase mode to the plasma frequency (see Sec. II B), completely suppress the contribution of phase fluctuations to D(T ). However, in Ref. [6] it has been shown that the d-wave symmetry of the order parameter plays a crucial role in determining the behavior of the phase mode in the presence of both long-range Coulomb interactions and dissipation. Even though here we ignore dissipative effects, for the sake of completeness we still compare s-wave and d-wave superconductors. Once again, the formalism of the effective action makes this comparison, to some extent, particularly simple, and of transparent physical interpretation.
The plan of the paper is the following. In Sec. II A we derive the Gaussian phase-only action starting from the continuum model, in the absence of long-range Coulomb forces. The case of the charged system is dealt with in Sec. II B. Although most of the results presented in Sec. II are standard, our scope here is to introduce the formalism and a classification scheme for the interaction vertices, which turns out to be particularly useful in dealing with non-Gaussian terms. In Sec. III we apply the same procedure to a lattice model, which does not appear as a trivial extension of the continuum model. In Sec. IV we calculate the anharmonic terms, providing a classification of the diagrams (Sec. IV A), and determine the one-loop correction to the stiffness within the continuum model, both at T = 0 (Sec. IV B) and in the classical limit (Sec. IV C). The results are then compared with the results of the XY model. In Sec. V the same analysis is carried out in the lattice case. The strong-coupling regime of the lattice model is discussed in Sec. VI, and the results are summarized in Sec. VII. The details on the derivation of the phase-only action are reported in App. A, whereas in App. B we discuss the connection between the coefficients of the Gaussian action and the gauge-invariant electromagnetic response functions.

A. The neutral system
To introduce the general formalism adopted in this paper, and for the sake of completeness, we shortly discuss the standard formal steps to derive the phase-only action within the continuum BCS model H = H 0 + H I , with Most of our results are obtained for generic spatial dimension d of the system, unless explicitly indicated. In Eqs.
(3)-(4), c (+) σ is the annihilation (creation) operator of an electron of spin σ, Ω is the volume, and U > 0 is the pairing interaction strength. Hereafter we assume = k B = 1. The factor w(k) controls the symmetry of the Cooper-pair wave function, e.g., w(k) = 1 for s-wave superconductors, and w(k) = cos 2φ, with φ = atan(k y /k x ), for d-wave superconductors, in the continuum case and for d = 2. Note that, in the case of a d-wave symmetry, four nodes are present in the quasiparticle excitation gap, and the existence of gapless excitations is relevant in determining the low-temperature thermodynamics of the system.
We first discuss the case of a neutral system, while the Coulomb interaction will be introduced in Sec. II B, to deal with the charged case. As customary, the microscopic effective model for the collective modes is derived by considering the action corresponding to the Hamiltonian (3)-(4), within the finite-temperature Matsubara formalism, where τ is the imaginary time, β = 1/T , and ξ k = ǫ k − µ is the band dispersion with respect to the chemical potential µ, with ǫ k = k 2 /2m in the continuum case. To obtain the effective action in terms of the order-parameter collective degrees of freedom, the interaction is decoupled in the particle-particle channel by means of the Hubbard-Stratonovich transformation, introducing the auxiliary complex field ∆(x, τ ), as explicitly reported in App. A. It is then possible to make the dependence on the phase of ∆(q, τ ) explicit in the action, by performing the gauge transformation Eq. (A2) on the fermionic fields. Indeed the dependence on the phase is eliminated from S I , and made explicit in the free part of the action S 0 , which now readsS whereτ i are the Pauli matrices and the operator acts on the fermionic Nambu spinor φ(x) defined as the column vector (c ↑ (x), c + ↓ (x)). It is worth noting that the field θ appears in the action only through its time and spatial derivatives, as it is expected for a Goldstone (massless) field. Moreover we point out that, while in the s-wave case the transformation (A2) completely eliminates the phase θ from S I , in the d-wave case a residual dependence survives. Nonetheless, in App. A, where we report some more detail, we show that this residual dependence can be safely neglected in the hydrodynamic limit, which is the relevant regime when discussing the low-energy properties of the collective mode. Thus one can now integrate out the fermions, which appear quadratically inS 0 + S I , leading to the effective action (A7), written in terms of the collective variables only. This standard procedure allows one to analyze the starting problem (3)-(4) from a different point of view. Instead of perturbatively expanding the Hamiltonian in powers of the coupling U , the effective action (A7) is expanded by assuming that the fluctuations of the fields around their saddle-point values are small, and that the variations of the phase in space and time are slow. Moreover, if one is interested in the dynamics of the phase at low temperature, the fluctuations of the modulus |∆| can be neglected. As far as the specific problem of high-T c superconductors is concerned, this approximation seems particularly appropriate for underdoped cuprates, which exhibit an anomalously large binding energy when compared to the critical temperature. 1 On the other hand, one should expect that modulus fluctuations can still be relevant at distances smaller than the typical length scale over which |∆| fluctuates. It is then natural to introduce a spatial cut-off for phase fluctuations, of the order of the coherence length ξ 0 . At distances greater than ξ 0 the phase fluctuations are the only relevant degrees of freedom, and a description in terms of the phase-only action becomes meaningful. Accordingly, the cut-off ξ 0 for the phase-only action is defined as the characteristic length scale of the spatial decay of the correlation function for |∆|. This point is delicate, since the naive identification of the coherence length ξ 0 with the Cooper-pair size ξ pair is meaningless in d-wave superconductor, where ξ pair diverges at T = 0, due to the presence of the gapless quasiparticle excitations mentioned above, as it is discussed in Ref. [15].
Starting from Eq. (A7) of App. A, with |∆| = |∆|(0) + δ|∆|(q), and neglecting the fluctuations of the modulus δ|∆|(q), the sum U −1 q |∆|(q) 2 reduces to its q = 0 value. As a consequence, the self-energy (A8) does no longer depend on |∆|. The action (A7) is then decomposed as is the mean-field BCS action, and S ef f is the phase-only which takes into account phase fluctuations around the BCS saddle point. HereĜ 0 is the BCS matrix defined in Eq. (A9), and the trace Tr is explicitly given in k-space by where tr is the trace in the Nambu space only andΣ(k) is the Fourier transform of the expression (7). Since the self-energy (7) contains powers of the field θ up to the second, the expansion (8) is not directly an expansion in powers of θ and the terms with the same power of θ must be recollected, leading to the expansion which is completely equivalent to Eq. (8), and includes both the Gaussian (k = 2) and the anharmonic (k > 2) terms. Each term in the effective action (9) is given by a closed fermionic loop A k with k incoming lines corresponding to fluctuating θ fields. In forming the fermionic loops, according to Eq. (7), we have three kinds of fermion-θ vertices, depicted in Fig. 1: (a) The vertex with a single incoming dotted line corresponds to the insertion of the time derivative of the phase, i.e., the termΣ 1 B = (i∂ τ θ/2)τ 3 of Eq. (7). (b) The vertex with two incoming dashed lines corresponds to the insertion of the gradient squared of the phase, i.e., the termΣ 2 B = [( ∇θ) 2 /8m]τ 3 of Eq. (7). (c) The vertex with a single incoming dashed line corresponds to the insertion of the spatial derivative of the phase associated to the spatial derivative of the incoming and outgoing Green functions, i.e., the termΣ 1 The terms generated in the effective action by the self-energyΣ B =Σ 1 B +Σ 2 B , associated to theτ 3 matrix, are essentially the same which are met when deriving the phase-only action in a bosonic superfluid system. 17 We therefore refer to these terms as "bosonic" contributions. TheΣ 1 F term generates instead contributions which are strictly related to the existence of quasiparticle excitations. Therefore, we refer to them as "fermionic" contributions. Clearly, the expansion (8) also generates "mixed" contributions, obtained by combinations ofτ 0 andτ 3 terms of the self-energŷ Σ. Such a distinction will be relevant in the following, while discussing the derivation of the anharmonic terms of the phase-only action.
At Gaussian level (k = 2) the effective action describes the collective mode of the system. Beyond the Gaussian level (k > 2), interaction terms in the field θ appear, and we obtain the quantum generalization of the anharmonic terms of the XY model, Eq. (1). As in the classical case, the anharmonic terms are responsible for the reduction of the stiffness due to phase fluctuations. Before discussing this issue, we analyze the structure of the Gaussian phase propagator in the quantum regime. The Gaussian effective action is derived with long but straightforward calculations, and reads The coefficients Λ ρρ , Λ ρJ and Λ JJ , which are depicted in Fig. 2 and whose explicit expressions are reported in App. A, represent respectively the density-density, density-current, and current-current bubbles, evaluated with the BCS Green functions. Since the BCS approximation explicitly breaks the gauge invariance, 18 the Λ functions do not coincide with the physical correlation functions K 14 , which must obey the gauge invariance. However, as it is shown explicitly in App. B (where the general relationship between the Λ and the K bubbles is obtained), it is possible to The correspondence between different linestyles and the various θ insertions is the same as in Fig. 1. Notice that the full current-current correlation function Λ αβ J J = −ρMF /mδ αβ + Λ αβ jj is given by the sum of the first and last diagram.
relate the coefficients of Gaussian effective action with measurable quantities such as the bare compressibility κ(T ) and the stiffness D(T ) = ρ s (T )/m of the system, in the hydrodynamic limit. 19 Indeed, in this regime, we obtain the expression The Bogoljubov sound mode of the neutral system, with velocity c s = D(T )/κ 0 (T ), is thus recovered. It is worth noting that, since the phase and the density are conjugate variables, the collective mode which appears in the phasefield propagator is the same as the density mode (see, e.g., Ref. [20]), as it can be explicitly checked by deriving the physical density-density correlation function K ρρ according to Eq. (B1). By comparing Eq. (11) with the classical XY model, one can recognize the same term ρ s (∇θ) 2 which appears in Eq. (1). However, the quantum Gaussian propagator for phase fluctuations contains also the Ω 2 m term of Eq. (11) which takes into account the dynamic effects. Moreover, according to the results of App. B and to Eqs. (A13)-(A14), the temperature dependence of the bare stiffness D(T ) which appears in Eq. 11 is entirely due to the quasiparticle excitations. Indeed, the coefficient Λ ab jj in Eqs. (A13)-(A14) describes the normal-fluid component, whose contribution increases as the temperature increases In the absence of dissipation the coefficients Λ are not analytic at q = 0 and finite temperature, since in this case the static (Ω m = 0, q → 0) and the dynamic limit (q = 0, Ω m → 0) are different. This issue has been recently addressed in Ref. [13]. Since we are interested in the T = 0 correction to the superfluid density we shall not discuss further the consequences of such non analyticity, which would disappear in any case in the presence of dissipation.

B. The charged system
When dealing with a charged system, the long-range Coulomb interaction between the fermions, must be included in the Hamiltonian (3)-(4). V (q) is the Fourier transform of the three-dimensional Coulomb repulsion V (r) = e 2 /ε B r, projected onto a d dimensional system (e.g., d = 2 for a single layer). Thus, for generic d, V (q) = λe 2 /|q| d−1 , where λ is a constant which depends on the dimension d, e.g., λ = 4π/ε B for d = 3 (isotropic three-dimensional system) and λ = 2π/ε B for d = 2 (single layer). Here ε B is the dielectric constant of the ionic background. The background ensures overall charge neutrality and cancels out the apparent divergence of the Coulomb interaction V (q) at q = 0. We point out that Eq. (13) generically describes a density-density interaction mediated by the potential V (q). The particular choice of V (q) distinguishes the short-range from the long-range case. After including the term (13) in Eq. (5), we follow the standard procedure and introduce a Hubbard-Stratonovich field ρ HS associated to the electron density, thus decoupling the Coulomb interaction in the particle-hole channel. Since the short-range pairing in Eq. (4) is important for small center-of-mass momentum, while the Coulomb effects are important for small momentum transfer, one can reasonably assume that the breakup of the actual interaction in this manner is physically sensible and does not lead to any "overcounting". Note that in principle also the pairing term should be decoupled in the particle-hole channel, but these terms are safely negligible, at least for weak and intermediate U . Nevertheless, in Sec. VI we shall discuss the subtleties of such a procedure, in the case of strong pairing interaction.
As described in App. A, by separating ρ HS (q) = ρ HS (0)+δρ HS (q), and by neglecting the fluctuations of the modulus δ|∆|, the expansion of the phase-only action around the BCS saddle point leads to S charged = S MF + S ef f (θ, ρ HS ), where now andΣ The Hubbard-Stratonovich term δρ 2 HS (0)/V (0) does not appear in the mean-field action S MF as it is supposed to be canceled out, in order to preserve the charge neutrality of the system, by an equal (and opposite in sign) contribution coming from the interaction with the ionic background. This term survives instead in the case of shortrange interaction. Since the self-energyΣ c depends now also on the second Hubbard-Stratonovich field δρ HS , to obtain the phase-only action, one has to proceed in two steps, as is explicitly shown in App. A: first, the Gaussian action in both θ and δρ HS is derived and, then, after integrating out the field δρ HS , one obtains the Gaussian phase-only action which generalizes Eq. (10) for a charged system. The coefficients L, whose definitions are reported in App.A, are simply the coefficients Λ of the unscreened phase-only action (10) dressed by the density fluctuations within the random-phase approximation (RPA). This represents the only difference between the Gaussian actions (10) and (16), which have therefore the same formal structure. The inclusion of Coulomb interaction obviously reflects on the collective mode of the system. Indeed, the expression of the Gaussian action (16) in the hydrodynamic limit differs from Eq. (11) only in the coefficient of the Ω 2 m term, since both Λ a ρj and L a ρj vanish in the static limit, and L ab JJ behaves as Λ ab JJ , whereas L ρρ gives the RPA compressibility of the charged system which vanishes for q → 0. As a consequence, the collective mode described by the action is the plasma mode in d dimensions. For example, in the three-dimensional case Eq. (18) reads and the plasma mode is characterized by a finite energy at zero momentum, Ω P = 4πe 2 D(T ) = 4πe 2 ρ s (T )/m. At T = 0, where ρ s coincides with the particle density, Ω P is exactly the plasma frequency of a three-dimensional charged system.
The vanishing of the compressibility (17) in the presence of density-density interactions is peculiar of the Coulomb case. If one considers a short-range interaction V (q) =Ṽ , Eq. (17) simply gives the RPA dressing of the bare compressibility which inserted into Eq. (16) leads to the Gaussian action In the presence of short-range density-density interactions the collective mode is still a sound mode (known as

III. THE EFFECTIVE ACTION FOR THE LATTICE MODEL
We extend the standard procedures discussed in Sec. II, to deal with the more realistic case of a lattice model. This extension does not reduce to a trivial modifications of the continuum results, the most important differences being the appearance of vertices with more than two incoming dashed lines, and the impossibility to absorb the dressing of the fermionic lines in a chemical potential shift.
We first consider the neutral case, the generalization to the charged case being straightforward. The main difference between the continuum and the lattice model is in the kinetic term of the Hamiltonian. Here, we rewrite the free action (see Eqs. (3), (5)) introducing a hopping t between nearest-neighboring sites on a lattice of spacing a = 1, so that the free-electron dispersion now reads ξ k = −2t j cos k j − µ, with j = x 1 , . . . , x d . The different form of the kinetic term modifies the effect of the gauge transformation on the free action S 0 . As we showed in Sec. II, after the transformation (A2) the phase disappears from S I and modifies the free action S 0 according to Eq. (6). As a consequence, to obtain the self-energyΣ within the lattice model we analyze the effect of the gauge transformation on the free action (21), the subsequent steps being unchanged. We find two differences with respect to the continuum case: (i) The termsΣ 1 F andΣ 2 B of the self-energy (7) arise, after the gauge transformation, from the kinetic term ∇ 2 /2m of the free action. Since the analogous of the velocity k/m of the continuum case is ∂ξ k /∂k in the lattice case, it follows that for a lattice model each insertion of a spatial derivative of θ in the fermionic loops is associated to a same-order k-derivative of the band dispersion ξ k .
(ii) On the lattice the coupling of the fermions to θ induced by the gauge transformation of the kinetic term generates higher than second-order derivatives of θ in the self-energyΣ. In the language of Feynman diagrams this means that besides the vertices with one or two incoming θ-lines depicted in Fig. 1, we should add all the vertices with n ≥ 3 incoming θ-lines, each associated with a n-order k-derivative of ξ k , see Fig. 3. Observe that no mixing between different components (a = x, y, · · ·) of the k-derivative of ξ k is produced, at least in the case of nearest-neighbor hopping.
As a consequence of (i) and (ii), the self-energyΣ(k) appearing in Eq. (8) for S ef f (θ) reads noŵ where δ(q) is a delta-function andτ (n) corresponds toτ 0 for n odd (fermionic-like contributions) and toτ 3 for n even (bosonic-like contributions). The functions sin(q i /2) in Eq. (22) reduce to the factors q i /2 of the continuum case at FIG. 3: Vertices arising from the kinetic term of the lattice model after the gauge transformation. Each vertex with n θ-lines is associated to a same-order derivative of the band dispersion. The first diagram is present also in the continuum case (see diagram (b) in Fig. 1), with a factor 1/m. small momenta, while preserving the lattice periodicity. Since we are interested in deriving the hydrodynamic action, in the following they will be approximated by q i /2.
In this section we focus on the generalization of the Gaussian action derived for the continuum model to the lattice case. To this purpose, it is sufficient to take into account the modifications described in the item (i). Indeed, at Gaussian level, only fermionic loops with two θ lines contribute, i.e., the same diagrams depicted in Fig. 2, which acquire however a different meaning. In particular, the first fermionic loop of Fig. 2, which corresponds to the diamagnetic contribution to the current, does not correspond to the particle density divided by the mass, resulting instead (in the specific case of a square lattice with nearest-neighbor hopping) proportional to the average kinetic energy E kin of the particles where we indicate with E k the superconducting quasi-particle energy defined in App. A, and, for the sake of clarity, we made the dependence on the lattice spacing a explicit. As a consequence, the current-current correlation function Λ JJ becomes now Λ ab JJ = −T δ ab + Λ ab jj (q), where it is assumed that in the definitions (A12)-(A14) each factor k a /m is substituted with ∂ξ k /∂k a and the lattice dispersion ξ k appears now in the quasiparticle energy E k . The stiffness, defined as the static limit of the transverse part of Λ JJ , is the lattice analogous of Eq. (12), i.e.
At Gaussian level the only difference in the phase-only action between the continuum and the lattice case is in the definitions (12) or (24) of the superfluid stiffness, the remaining steps being the same. In the following, for the sake of simplicity, we write down the generic Gaussian actions as where χ indicates the static limit of the density-density correlation function χ(q) ≡ Λ ρρ (q), i.e., the compressibility, for neutral systems which has to be replaced by χ LR (q) ∼ 1/V (q) in the presence of Coulomb interaction. As a consequence, in the hydrodynamic regime, the Gaussian propagator for phase fluctuations will be generically given by The stiffness which appears in the Gaussian propagator is the bare one, which we indicated with a subscript D 0 in Sec. I. However, to simplify the notation, in the following we drop the subscript, while keeping in mind that the Gaussian value D is corrected by anharmonic phase fluctuation according to D → D +δD, where δD will be calculated within the one-loop approximation. When quantum corrections are finite, i.e. δD(T = 0) = 0, the zero-temperature stiffness does not coincide with the Gaussian value, defined by the T = 0 limit of the BCS formulas (12) and (24). The consequences of such an occurrence will be discussed in Sec. IV.

IV. ONE-LOOP CORRECTIONS TO D WITHIN THE CONTINUUM MODEL
Before deriving the anharmonic terms in the phase-only action, we discuss the phase-fluctuation correction to the stiffness at T = 0 in the quantum XY model, which was briefly anticipated in Sec. I. The quantum XY model includes quantum effects for the phase through a suitable modification of the Gaussian phase propagator only, while the interaction terms are powers of (θ i −θ j ) ∼ ξ 0 | ∇θ|, as in the expansion (1). Consequently δD = −(Dξ 2 0 /2d) ( ∇θ) 2 , which generalizes Eq. (2) in d dimensions. Since ( ∇θ) 2 = q 2 P (q), with the quantum propagator P (q) defined in Eq. (28), we find where 21 We considered the neutral case as an example, the generalization to the charged case being straightforward.
To compare this result with that expected for the microscopically derived phase-only action, a preliminary distinction arises between continuum and lattice models. The Gaussian phase propagator is given by Eq. (25) in both cases, and the temperature dependence of D(T ) is due to quasiparticles only. In spite of this analogy, the stiffness D has a different meaning in the continuum or lattice case. In the former D = ρ s /m, with the constraint ρ = ρ s + ρ n at any T , whereas in the latter D is related to the mean kinetic energy, see Eqs. (12) and (24), respectively. In particular, at T = 0, The equality ρ s = ρ is expected to hold at T = 0 in the continuum model, since Galilean invariance ensures that at T = 0, where quasiparticle excitations are absent, there is no preferred reference frame in the system, and all the electrons contribute to the superfluid component. 22 Thus, Eq. (30) states that the Gaussian approximations preserves that Galilean invariance of the microscopic model. Furthermore, despite the fact that, beyond the Gaussian level, anharmonic terms induce corrections to the bare Gaussian stiffness, such corrections are expected to vanish at T = 0 in the continuum model, to preserve the Galilean-invariant relation ρ s (0) = ρ. In the lattice case, instead, the presence of finite corrections to the Gaussian value D(0) in Eq. (31) is not forbidden. This is exactly what happens in the quantum XY model, which is a lattice model, according to Eq. (29). As a consequence, the quantum XY model cannot be a good approximation for the phase-only model of the continuum system, since it misses the cancellation of the T = 0 correction to D due to anharmonic phase fluctuations. This requirement is instead fulfilled by expanding the microscopically derived S ef f (θ) beyond the Gaussian level, and deriving one-loop corrections to D, as we shall show in the following. The failure of the quantum XY model for continuum systems is not a purely formal question, as it rises the issue of its quantitative validity both in the lattice case at T = 0, and in lattice and continuum systems in the classical regime, above some crossover temperature T cl (which depends on the parameter of the microscopic model). This issue is of particular relevance when high-T c superconductors are considered. Indeed, although in charged systems T cl can be a quite large energy scale, of the order of the plasma energy, dissipation in d-wave superconductors strongly reduces T cl , as discussed in Ref. [6], making the classical regime experimentally accessible at much lower temperature. This situation then requires a critical discussion of the corresponding classical model for phase fluctuations, with respect to the standard assumption of considering the XY model. Here we do not discuss neither dissipative effects, nor the precise definition of T cl , but we directly identify the classical regime by considering only the Ω m = 0 contribution in the Matsubara sums which define the phase-fluctuation corrections δD and δD XY . In the case of the XY model, by considering Eq. (29), we find that in the classical regime i.e., δD XY (T ) ∝ T , as already observed in Sec. I. We shall find below that also δD is linear in T , but with a much smaller coefficient in the weak-to-intermediate-coupling regime (see also Ref. [12]).

A. Classification of the diagrams
We start our analysis with the continuum neutral system, while the long-range Coulomb forces will be included later. To determine the one-loop correction δD to the stiffness induced by anharmonic terms in S ef f (θ), we derive the contributions up to S 4 in Eq. (9). In general, if we separate the effective action in the Gaussian (S G ) and anharmonic (S anh ) part, and make a Taylor expansion of the latter, by means of the Wick theorem we can evaluate the two-field propagator as the resummation of all the connected diagrams obtained by fixing two θ legs and contracting the others with the Gaussian propagator P (q), where · · · CON indicates connected diagrams only. All one-loop corrections to P (q), i.e., diagrams with only one integration on bosonic momenta, are obtained by approximating n (−1) n /n![S anh ] n ≃ S 4 + [S 3 ] 2 /2. Moreover, since we are specifically interested in the correction to D, which is the coefficient of the gradient term in the effective Gaussian action, we select the diagrams which have at least two ∇θ lines, as depicted in Figs. 4 and 5. To obtain one-loop corrections to D we consider all the possible contractions of two θ lines in one of the fourth-order vertices in Fig. 5, or in the combination of two third-order vertices in Fig. 4. The above classification in bosonic, fermionic and mixed contributions turns out to be particularly convenient while evaluating all the resulting one-loop corrections to D. We can indeed distinguish: • Fermionic corrections, which are given by the self-energy-and vertex-like diagrams obtained from the 4a loop in Fig. 5.
• Bosonic corrections, which arise from two loops of the 3b type and from the loop 4e of Fig. 5, see Fig. 6.
• Mixed corrections, which are the remaining combinations of third-and fourth-order diagrams.
The distinction between fermionic, bosonic, and mixed contributions to δD is motivated a posteriori by the fact that all the contributions to δD coming from the fermionic and mixed diagrams cancel out each other at T = 0, whatever is the symmetry of the gap and both in continuum and lattice models. This result is connected to the fact that the temperature dependence of the fermionic and mixed bubbles is controlled by the thermal excitation of the quasiparticles, which are suppressed at low temperature. However, while for s-wave superconductors this suppression is exponential in temperature due to the presence of a finite gap in the excitation spectrum, for a d-wave superconductor the presence of nodal quasiparticles obliges a careful evaluation of the fermionic contribution to δD. For example, one can see that the fermionic contributions generated by the vertex (4a) in Fig. 5 can be written as T Ω q,Ωm where v k = ∂ξ k /∂k and the hydrodynamic limit of the Q bubble is where f ′′′ (x) is the third-order derivative of the Fermi function. In the s-wave case the integral in Eq. (34) vanishes exponentially at T = 0, while in the d-wave case, by evaluating the contribution of nodal quasiparticles to Eq. (34), which is the relevant contribution at low temperature, one finds that it diverges as 1/T . This can be easily seen by considering that the density of states of nodal quasiparticles is linear in energy It should be stressed that the divergence in the d-wave case of all fermionic loops with N ≥ 4 Green functions, when evaluated at zero external momenta, is a characteristic feature of systems having Dirac quasiparticles. Nevertheless the one-loop expansion of D is well defined: indeed by retaining all the q, Ω m dependence on the fermionic bubble Q one has no longer to deal with any divergence of fermionic corrections to D, which result even vanishing in the zero-temperature limit. 24 Special attention must be devoted also to the bosonic (d1 and d2) and fermionic (d3) diagrams depicted in Fig.  7. Indeed, these diagrams appear as corrections to the bare Green function in the diagrams of Fig. 2 which define the the mean-field current-current correlation function Λ JJ in Eq. (A13). In particular, d1 and d2 in Fig. 7 are phase-fluctuation corrections to the diamagnetic term −ρ MF /m, and d3 in Fig. 7 is a phase-fluctuation correction to the bubble Λ jj , which gives the quasiparticle depletion of ρ s (T ) at finite temperature. One can easily see that these terms are due to the shift of ρ with respect to the mean-field value ρ MF , as induced by anharmonic phase fluctuations. On the other hand, working at fixed density, one assumes ρ MF ≡ ρ and takes into account these contributions by shifting the chemical potential µ → µ + δµ with respect to the mean-field value, so that where we use the definition χ = ∂ρ/∂µ for the compressibility, and the free-energy contribution ∆F due to phase fluctuations is evaluated within the Gaussian approximation, While evaluating one-loop corrections to D, one takes into account the chemical potential shift (35) by adding the contribution which cancels out exactly the contribution of the diagrams in Fig. (7). This can be easily understood, e.g., at T = 0. Because at T = 0 all the mixed diagrams vanish, only d1 and d2 contribute to δD. Let us consider the diagram d2. Its contribution to δD is Indeed, the fermionic loop of d2, given by the trace in the previous equation, corresponds to the hydrodynamic limit of the same density-density bubble χ(q) which appears as the coefficient of the Ω 2 m term in the Gaussian propagator P (q) (28), see Eqs.(26)- (27). 25 On the other hand, while evaluating ∂∆F/∂µ according to Eq. (36), one finds that the Dq 2 term of Eq. (36) contributes with where we used Eq. (28) and D = ρ s /m = [ρ − ρ n (T )]/m so that, at T = 0, D = ρ/m and ∂D/∂µ = χ/m. As a consequence, by means of Eqs. (35) and (37) we find that which is exactly −δD d2 . The same holds for δD d1 , δD d3 and the contribution arising by differentiating χΩ 2 m and ρ n (T ) (at finite temperature) in Eq. (36), leading to a cancellation of δD µ with the diagrams of Fig. 7. The cancellation of fermionic and mixed contributions as T → 0 can be demonstrated with lengthly but straightforward calculations, so we will not discuss it here, where instead we analyze the behavior of the remaining bosonic contributions of Fig. 6, and in particular the role of dynamics in restoring Galilean invariance in the continuum model. We note in passing that the cancellation of the fermionic corrections has been already observed in Ref. [12], where however the authors did not include both mixed terms and corrections to δD coming from the third-order loops in S 3 (θ). While neglecting mixed term at T = 0 is not problematic, since they cancel out at T = 0, neglecting the dynamic third-order vertex 3b of Fig. 4 leads to neglecting the bosonic diagrams B2 and B3 of Fig. 6, whose role at T = 0 is crucial.
Roughly speaking, the diagrams depicted in Fig. 6, arising from the bosonic self-energyΣ B only, correspond to where the angular brackets indicate the average with respect to the Gaussian action. As a consequence, the first diagram B1, arising from a vertex ∼ ( ∇θ) 4 in the effective action S ef f , is the analogous of the term d α=1 (∂ α θ) 4 in the quantum XY model (1). Indeed, the term B1 generates the correction where χ(q) is given by Eq. (26). By evaluating Eq. (38) in the hydrodynamic limit, in analogy with Eq. (29), we find that which is again a finite correction, as in the quantum XY model, see Eq. (29). By comparing Eq. (29) and (39), one finds that the strength of the term ( ∇θ) 4 in the continuum model is controlled by χ/m 2 , and not ξ 2 0 D, as in the XY model. Nevertheless, at T = 0, the result (39) cannot be conclusive, since it would lead to ρ s = ρ at zero temperature. Indeed, by adding the contribution of the remaining diagrams B2 and B3 we find that the bosonic one-loop corrections to D reduce to In writing Eq. (40) we are relying on the fact that both the third-and fourth-order vertices needed to calculate δD are expressed in terms of the fermionic bubble χ(q), Eq. (26). Indeed, both the staticΣ 2 B = (( ∇θ) 2 /8m)τ 3 and the dynamicΣ 1 B = (i∂ τ θ/3)t 3 part of the "bosonic" self-energyΣ B are tied to the Pauli matrixτ 3 , leading to the same fermionic bubble (26) as a coefficient of the phase-only action. As observed in Refs. [7,9], this is the relevant consequence of the Galilean-invariant form of the bosonic τ 3 term appearing in (7). Evaluating Eq. (40) in the hydrodynamic limit, we find where ω q = D/χ|q| is the sound mode, and b ′ (x) = −βe βx /(e βx − 1) 2 . As a consequence of the addition of the B2 and B3 contributions, the pole at iΩ m = −ω q < 0, responsible for a finite contribution at T = 0 in Eq. (39), is canceled in favor of a double pole at iΩ m = ω q > 0, leading to Eq. (41), where the standard Bogoljubov reduction of ρ s /m in a superfluid bosonic system is recognized 17 . Thus, for T → 0 the integral in Eq. (41) vanishes, i.e., by fully including the dynamic structure of the interaction for the phase, we obtain that in a Galilean-invariant system ρ s = ρ at T = 0. Notice that the result (41) arises from an exact compensation between static and dynamic bosonic vertices for phase fluctuations, both at Gaussian level, in determining the structure of the propagator P (q), and beyond the Gaussian level, in weighting the relative contributions of static and dynamic parts.
where k F is the Fermi wave-vector. Nevertheless, also ρ s , i.e., ρ (if we neglect quasiparticle effects 27 ), can be defined through k F , at least in the weak-to-intermediate-coupling regime (more precisely as long as T c < E F , i.e., U < t) It then follows that the ratio between the coefficient of the static interaction term ( ∇θ) 4 in the continuum or in the XY model is which means that, apart from a factor d/2, we have As a consequence, in the classical limit, the one-loop correction to D derived within the effective action is smaller than the corresponding first-order correction in the XY model as far as k −1 F < ξ 0 , i.e., in the weak-and intermediatecoupling regime for the pairing interaction.
Finally, let us consider the effect of the Coulomb interaction between the electrons. As we showed in Sec. II B, the presence of a density-density interaction in the microscopic Hamiltonian implies the dressing of the coefficients of the Gaussian action by the RPA series of the Coulomb potential V (q). As a consequence, in the hydrodynamic limit, the bare compressibility χ is replaced by the dressed bubble χ LR (q) ≈ 1/V (q), and the sound mode ω q of Eq. (41) is converted into the plasma mode ω P q of the d-dimensional system (see Eq. (19) for the three-dimensional case). In the same way, when deriving the anharmonic terms in S ef f , we must now include the RPA density fluctuations in all the vertices for S 3 and S 4 . The one-loop corrections to D are formally identical to Eqs. (40)-(41), with χ → χ LR . Thus at T = 0 we recover again the cancellations (41) of the bosonic diagrams, with ω q → ω P q . At the same time, since χ LR (q) vanishes as q → 0, the classical (Ω n = 0) term in Eq. (40), i.e., δD B1 given by Eq. (39), reads whereas the quantum XY model leads to the same result (32) of the neutral case. Thus we can roughly estimate where ε F = k 2 F /2m is the Fermi energy and ε C = λe 2 k F is a characteristic Coulomb energy scale. As a consequence, while within the XY model the Coulomb interaction modifies only the low-temperature behavior of ρ s , within the continuum model it affects also the high-temperature classical regime.

V. ONE-LOOP CORRECTIONS TO D WITHIN A LATTICE MODEL
To extend the previous results to a lattice model, we need to take into account the discussion of Sec. III, in which we derived the phase-only action starting from a nearest-neighbor tight-binding model for electrons in a lattice.
At Gaussian level, the item (i) of Sec. III allows us to define the superfluid stiffness of a lattice model according to Eq. (24). In particular, a factor 1 8 Λ a , with appears in the fermionic bubble which defines the diamagnetic term, i.e., the coefficient of the ( ∇θ) 2 term in the Gaussian action for phase fluctuations. As a consequence, D is controlled at T = 0 by the mean kinetic energy T , Eq. (23), instead of ρ/m of the continuum case. In deriving anharmonic terms in S ef f (θ), we must take into account item (i) in the proper definition of the fermionic bubbles of Figs. 4 and 5, which give the vertices for the interacting phase-only model. Moreover, we need to include also the fourth-order vertex depicted in Fig. 8, which is peculiar of the lattice case, and arises from the n = 4 term of Eq. (22) (see also Fig. 3). At this point the calculation of one-loop correction to D follows the same steps described in Sec. IV. The cancellation of fermionic and mixed one-loop corrections to D at T = 0 still holds, thus we can again focus on bosonic corrections only. As far as the corrections in Fig. 7 are concerned, they are no more canceled out by the contribution δD µ due to the chemical-potential shift beyond the Gaussian level. In particular, while discussing in detail in Sec. IV the case of the term d2, we stressed that in the continuum model such a cancellation holds because the fermionic bubble of diagram d2 is the same bubble which one obtains by differentiating the diamagnetic part of D with respect to µ, so that ∂D/∂µ = χ/m. Indeed, since in the continuum case the factor Λ α = 1/m associated to the ( ∇θ) 2 insertion is a constant, one obtains the same two-line fermionic loop both by expanding S ef f up to S 4 , leading to d2, and by differentiating the single-line diagram which defines ρ. However, in the lattice case the factor Λ α is not a constant. As a consequence the diagram obtained by differentiating the diamagnetic term T ≡ D(0) with respect to µ has one insertion of Λ α , while the the diagram d2, obtained by expanding S ef f according to item (i), has a Λ α factor for each θ-vertex (diagram L 2 in Fig. 9). The same holds for the d1 diagram.
As a consequence, in the lattice case, the contributions to D come from the bosonic diagrams depicted in Fig. 9, and from the correction δD µ induced by the chemical-potential shift, which do not cancel anymore the diagrams L 1 and L 2 in Fig. 9. While evaluating δD L2 we must consider that we have 2(d − 1) bubbles with Λ α = Λ β in the expansion of S anh in d dimensions, so that where η ρE (q) is the fermionic bubble with threeĜ 0 lines and one insertion of Λ α , and the χ ab (q) bubbles correspond to the insertion of one (χ ρE ) or two equal/different (χ EE /χ EE ) factors Λ α in the two-line bubble. As usual, when the dependence on q is not explicitly indicated, we are considering the q → 0, Ω m = 0 (static) limit. Notice that, since (∂ 4 ξ k /∂k 4 α ) = −a 2 (∂ 2 ξ k /∂k 2 α ), theT which appears in δD L3 is proportional to the mean kinetic energy defined by Eq. (23),T = a 2 T = a 4 E kin (the dependence on the lattice parameter a is made explicit, for the sake of clarity in the forthcoming discussion).
According to Eq. (35), the chemical-potential shift is where η(q) is the three-line bubble without factor Λ α , χ is the compressibility, which carries no Λ α insertion, and we used χ ρE = ∂D/∂µ. By means of Eq. (37), we obtain Finally, the overall one-loop correction to D coming from Eqs. (45)-(49) and (50) is To make a comparison with the continuum case, Eq. (41), we evaluate the above equation in the hydrodynamic limit. This corresponds to calculate the χ ab bubbles at zero incoming momentum and by carefully evaluating the q → 0 limit for the terms which involve η ρE (q) and η(q). Indeed, by simply calculating Ωm P (q)Ω 2 m η(q) → η Ωm P (q)Ω 2 m one finds that Ωm P (q)Ω 2 m → ∞. Instead, by retaining the Ω m dependence in the fermionic bubble one ensures the convergence of the Matsubara sum, calculating the q → 0 limit afterward. Thus at T = 0 limit we obtain that the overall correction to D is finite in the lattice case, The T = 0 value of χ is given by Eq. (27). The others bubbles are defined as Thus in the lattice case the T = 0 mean-field value of the superfluid density is corrected by one-loop phase fluctuations according to Eq. (52). It is interesting to observe that the continuum case is recovered from Eq. (52) by performing the limit a → 0, t → ∞, keeping 2a 2 t = 1/m finite. Indeed, since ξ k = −2tǫ k − µ, with ǫ k = a cos(ak a ), Λ α = (∂ 2 ξ k /∂k 2 a ) = 2ta 2 ǫ k ≈ 1/m in this limit, and according to Eqs. (53)-(57), one finds χ EE ,χ EE → χ/m 2 , χ ρE → χ/m,η ρE →η/m, whileT = a 2 D → 0. Thus in this limit, Eq. (52) correctly reduces to the continuum (Galilean-invariant) result δD = 0 of Eq. (41).
We now compare δD and δD XY in the lattice case. Differently from the continuum case, this comparison is now meaningful both in the classical limit and at T = 0, since now we find finite corrections to D. By comparing Eq. (52) and Eq. (29) at T = 0, we see that roughly speaking the strength ∼ ξ 2 0 D of the interaction in the quantum XY model is given by the quantity in square brackets in the first line of Eq. (52), i.e., it is of order ∼ χ EE ,χ EE , χ 2 ρE /χ. Following the same arguments which led us to Eq. (43), we can estimate where χ EE ,χ EE , χ 2 ρE /χ have the same role as χ/m 2 in Eq. (43), and k F ∼ 1/a in the lattice case. However, in the lattice case the dynamic contribution controlled by χ ρE /χ 2 does not cancel with the static contribution, and both contribute to δD at T = 0. Analogously, one finds that T = a 2 D ≈ D/k 2 Eq. (42), ∆ is the maximum gap value, and we used the fact that at weak coupling k F ξ 0 ≃ k F v F /∆ = µ/∆ ≫ 1. 15 Taking into account the prefactor g/ζχ in the second line of Eq. (52), we have whereη ρE ∼η/m. Thus, the contributions coming from the three-line fermionic bubbles are subleading, and at T = 0 Eq. (52) leads to an estimate of the same order of Eq. (44), within a numerical factor. The presence of Coulomb forces does not change qualitatively this conclusion, and introduces only minor quantitative corrections. This is due to the fact that the RPA expressions for χ EE ,χ EE have a finite limit for q → 0, contrary to the continuum case, When δD is evaluated with the long-range phase propagator P (q), one sees that the leading terms in powers of 1/(ξ 0 k F ) come from the finite q = 0 limit of χ EE ,χ EE in Eq. (60), and from theT term, all of order 1/(ξ 0 k F ) 2 with respect to the coefficient ξ 2 0 D of the XY model, whereas so that dynamical terms are subleading at T = 0 in the long-range case. This is clearly true as far as the constant in the parenthesis of Eq. (60) is a finite quantity. Indeed, when it vanishes by approaching the continuum limit, all (static and dynamics) terms are of the same order, canceling each other as discussed above. Finally, the same result (59) holds in the classical limit, both for the neutral and for the charged system. As we discussed above, at high temperature only the static correction to D contributes, with a coefficient controlled in the neutral case by χ EE ,χ EE ,T , or by their q = 0 finite limit in the charged case, which play the same role (and have the same estimate) as χ/m 2 in the continuum case.
As far as the effect of bosonic degrees of freedom is concerned, our result (59) suggests that the XY model can eventually lead to an overestimate of phase-fluctuation effects with respect to the result obtained by means of the microscopically derived effective action. Indeed, according to Eq. (59) the depletion of D due to phase fluctuations can be quite small at weak and intermediate coupling (particularly in the BCS limit) both in the quantum and in the classical regime.
On the other hand, as far as the temperature dependence of δD for temperatures 0 < T < T cl is considered, one should take into account also the effects of fermionic and mixed contributions. As we discussed above, quasiparticles contribute to the depletion of D at finite temperature already at Gaussian level, with a temperature dependence determined by the symmetry of the order parameter, which is linear in T in d-wave superconductors. One would expect that at T > 0 the one-loop corrections to D coming from fermionic and mixed diagrams correct quantitatively but not qualitatively the D(T ) dependence due to quasiparticles at Gaussian level. However, a detailed analysis of the temperature dependence of D(T ) as arising from bosonic, mixed and fermionic contribution at low temperature is a much harder task. From one side, one has to face the difficulty of evaluating the q = 0 limit of these expressions, as we discussed in Sec. IV A. On the other side, at finite temperature, non-local terms appear when deriving the phase-only action, making the hydrodynamic expansion not well-defined. These are the so-called Landau singularities, arising from the non analytic behavior at small q of the Λ JJ . 13 In any case the analysis of these issues requires separate investigation and is not addressed in the present work.

VI. THE STRONG-COUPLING REGIME
The results (59) for the lattice case was derived by assuming 1/(k F ξ 0 ) as a small expansion parameter. In particular, this assumption allowed us to estimate the fermionic bubble according to Eq. (42), i.e., χ ≃ N (ε F ), and to neglect the contribution to δD of theη ab bubbles, see Eq. (58). One should expect that such approximations are reliable in the weak-or intermediate-coupling regime, where the correlation length ξ 0 is larger than the lattice parameter a, which sets the scale of the Fermi wave-vector k F ∼ 1/a. However, in the strong-coupling regime for the pairing interaction several modifications should be included. The fermionic bubbles can no longer be estimated as the density of states at the Fermi level; the correlation length is generically expected to be of the same order of magnitude as the lattice parameter; one should include both the RPA fluctuations, also induced by U , in the particle-hole channel, and the fluctuation of |∆| ≃ (U/2) ρ(2 − ρ), which fluctuates because ρ fluctuates. 7 This last issue can be better understood by analogy with the superfluid bosonic systems, where the order parameter is directly the square-root of the density, and one expects that, in the strong-coupling regime, the fermionic system can be mapped onto an effective bosonic system, where Cooper pairs act as bosons with a weak residual repulsion. The above points are intimately connected, as we already discussed in Ref. 15. There, the behavior of the correlation length, properly defined as the characteristic length scale for the fluctuations of the modulus of the superconducting order parameter, was analyzed in the strong-coupling regime. The main finding is that ξ 0 attains a value of the order of the lattice parameter a at finite densities, and diverges as 1/ √ ρ in the low-density regime ρ → 0, recovering the well-known behavior of a bosonic system. 28 In the following we consider the s-wave case, which allows us for a transparent analytical treatment. We first address the neutral case, i.e. the model (3)-(4) with w(k) = 1 (negative-U Hubbard model), and then we extend our results to the charged case. In the presence of d-wave pairing the analysis becomes more difficult, even though the physics is not expected to be different from the s-wave case, as we briefly discuss at the end of this section.
As we said at the beginning of Sec. II A, the Hubbard-Stratonovich decoupling is a useful tool which allows us to substitute the microscopic model with an effective action written in terms of the relevant collective degrees of freedom. However, the choice of the relevant variables depends on the physical properties of interest and on the energetic scales involved in the problem. For example, the negative-U Hubbard model can be decoupled both in the particle-particle channel, as we did, and in the particle-hole channel. The latter decoupling would lead to an effective action depending on the particle density as well, S ef f = S ef f (θ, δ|∆|, δρ). 29 The resulting action is formally the same that we obtained by decoupling the additional interaction Hamiltonian (13) in the particle-hole channel, but with V (q) replaced by the same −U/2 which appears in the Hubbard term. The relevance of density fluctuations within the Hubbard model at strong coupling is addressed in Refs. [7,15]. In deriving S ef f up to Gaussian level a term proportional to the product of amplitude and density fluctuation, ∼ δ|∆|δρ, appears. The coefficient of such a term in the Gaussian action is vanishingly small at weak coupling. This means that density and modulus fluctuations decouple at weak and intermediate coupling, and one can safely neglect modulus fluctuations without loosing information on the contribution of the density fluctuations. As a consequence, one can derive the phase-only action according to the procedure described in the previous sections. However, at strong coupling the situation is more involved: as it is shown in Refs. [7,15], by increasing U the density-modulus coupling become sizable, and one finds that the density of particles and the modulus of the order parameter experience the same fluctuations. Since the phase mode has the same behavior of the density mode, in the strong-coupling regime we should also take into account amplitude fluctuations to recover a consistent description of the phase fluctuations.
In Sec.II B we found that in deriving the effective action for the phase of the order parameter the inclusion of density fluctuations reflects in the replacement of the "bare" coefficients Λ of the Gaussian action with their counterparts L evaluated within RPA resummation in the particle-hole channel. Beyond Gaussian level one expects the RPA dressing of the various fermionic loops which appear as coefficients of the effective action. As far as the δD correction in Eq. (52) is concerned, this leads, e.g., to substitute the χ ab fermionic bubbles with the RPA resummation with the potential −U/2 of the corresponding density or current irreducible correlation functions χ irr ab (q). We are then left with the problem of including modulus fluctuations. In principle, to take them into account one should derive the Gaussian effective action S G (θ, δ|∆|, δρ) by including all the fields and then integrate out modulus fluctuations as well, see, e.g., Refs. [7,15]. However, the contribution of amplitude fluctuations to the static long-wavelength limit of the fermionic bubbles χ irr ab (q → 0, ω n = 0) can be easily derived in a different way. Let us consider, e.g., the density-density correlation function χ irr . Since χ irr = ∂ρ/∂µ (see also Sec. IV), we can deduce it from the mean-field expression (A15) of the density of particles, by taking into account the dependence on µ of the gap amplitude |∆|(0) ≡ ∆ 0 itself, which appears in the definition of E k = ξ 2 k + ∆ 2 0 (see also Eqs. (A10)). At weak coupling ∂∆ 0 /∂µ can be neglected, and one finds the definition (27) of χ irr . At strong coupling, where modulus and density fluctuations are proportional, by working at fixed number of particles the chemical-potential variations needed to preserve the particle number reflect into amplitude variations, and ∂∆ 0 /∂µ gives a significant contribution to the irreducible bubbles. We then have (at T = 0): where ∂∆ 0 /∂µ can be obtained from the self-consistency saddle-point equation by differentiating both sides with respect to µ, , so that In the limit of strong pairing interaction, we can evaluate the previous bubble for t/U ≪ 1, by taking into account the strong-coupling solution of the coupled saddle-point equations for the amplitude ∆ 0 Eq. (61) and for the particle number Eq. (A15), where we introduced the short-hand notation δ ≡ 1 − ρ and the small parameter α ≡ 2t/U . We then obtain Thus, including also the RPA resummation of χ irr in the particle-hole channel, we finally get .
The expression (63) for the density-density electronic correlation function at q = 0 is the same which one would obtain by explicitly integrating out both the δρ and the δ|∆| field in the action S(θ, δ|∆|, δρ). The latter procedure was indeed adopted in Ref. [15] while deriving the propagator for the Hubbard-Stratonovich field associated to the density. The result (63) shows a large increase of χ at large U : this could be expected, since the electrons are strongly paired, naturally increasing the compressibility of the system. As far as the other fermionic bubbles are concerned, we follow the same steps. In analogy with the previous discussion where the stiffness is defined in Eq. (24), and at T = 0. By taking the derivative with respect to the hopping t one includes also ∂∆ 0 /∂t, as obtained again from Eq. (61). In the strong-coupling limit By introducing the dressed potentialŨ = (U/2)/[1 − (U/2)χ irr ] = U/4dα 2 , in analogy with Eq. (63), we find The inclusion of the modulus fluctuations always adds a second term to the weak-coupling definition of the irreducible bubbles, see, e.g., Eq. (62). This additional term is never of lower order in α with respect to the term obtained with the weak-coupling definitions (53)-(55). Since one expects that the same holds for the η bubbles, let us first estimate the RPA resummation of the weak-coupling definitions (56) and (57), which do not include modulus-fluctuation contribution. The resulting corrections to D, arising from the second line of Eq. (52), are subleading with respect to the contribution coming from the χ ab bubbles, first line of Eq. (52). Since modulus-fluctuation contributions do not change this result, we can neglect the correction δD coming from the η bubbles. By using the expressions (67)-(69), the strong-coupling expression and the fact thatT = a 2 D, we find that Eq. (52), at large U/t, reduces to where we reintroduced explicitly the dependence on the lattice parameter a, to allow one for a direct comparison with the result of the XY model, Eq. (29). We then obtain that By taking into account the previous results for ξ 0 , 15 we find that, since at finite densities ξ 0 ∼ a, the XY model appears as the proper "strong-coupling" effective quantum model for phase fluctuations in a lattice. A noticeable exception is found in the extreme low-density regime, where ξ 0 diverges as 1/ √ ρ. Moreover, since g ∼ 1/ξ d+1 0 , Eq. (70) shows that in the low-density limit the relative superfluid-stiffness correction δD/D vanishes as ρ d/2 . Observe that instead the XY model would lead to a relative correction δD XY /D ∼ ρ d/2−1 which does not vanish at ρ = 0 in d = 2.
In the classical regime, as we already discussed in the weak-coupling case, only static anharmonic terms survive, leading to Since the quantity in parenthesis behaves as the microscopically derived (ξ 0 /a) 2 at all densities, 15 up to a numerical constant of order one, we find that the in the classical regime the XY model turns out to be appropriate, regardless of the density. Observe that the numerical difference comes out from the fact that the XY model attributes a somewhat arbitrary coefficient to the anharmonic term ( ∇θ) 4 . If one assumed in Eq. (1) higher order harmonics as D n 4 <i,j> (1 − cos θ ij ) n , (n ≥ 2), which do not modify the Gaussian term, by adjusting the D n coefficients one could reproduce the correct ( ∇θ) 2n term as it arises in the microscopic expansion. In such a case the classical limit of this extended XY model and of the microscopic model would lead to exactly the same result. As it has been noted recently in Ref. 16 , the mapping of the classical phase-only action of a lattice systems on an extended XY model can also lead in two dimensions to an enhanced fluctiation region near the Kostelitz-Thouless transition. Let us consider now the case when also the Coulomb potential is present. The effect of modulus fluctuations on the irreducible fermionic bubbles is the same as already discussed above , so Eqs. (62) and (64)-(66) are still valid. As far as density fluctuations are concerned, we just have to take into account that the RPA resummation of the irreducible bubbles in the particle-hole channel must be performed with the full potential V (q) − U/2. As a consequence, even though we are considering a strong-coupling (i.e. large U ) limit, as q goes to zero the Coulomb interaction V (q) is always predominant. Thus in the hydrodynamic limit the contributions to the correction δD only come from the static interaction vertices, as already discussed in the weak-coupling limit, since the χ ρE bubble vanishes for small q, while the current-current correlation functions attain a finite value. More explicitly, by performing the RPA resummations (67)-(69) withŨ replaced by the full dressed potentialŨ →Ṽ (q) = −[V (q) + U/2]/[1 + (V (q) − U/2)χ irr ], it can be easily seen that As a consequence, in Eq. (51) for δD a coefficient (3χ irr EE +T )/8 = a 2 D/2 survives. This coefficient plays the same role as ξ 2 0 D/2 in the XY -model, see Eq. (29) (which however should be evaluated with the long-range propagator P (q)). Thus, Eq. (71) still holds, both at zero temperature and in the classical regime, where by definition only static anharmonic terms survive. Observe that this result is again not trivial, because only properly including modulus fluctuations one can find the equivalence between χ irr EE ,T and a 2 D in the strong-coupling regime for the pairing potential U .
Finally, we briefly comment on the extension of the previous results to the case of d-wave symmetry of the order parameter. As far as modulus fluctuations in the particle-hole channel are concerned, in the d-wave case this would reflect in a Hartree-Fock-like correction to the band dispersion. Since this makes the analytical treatment not viable, both for the coherence length and for δD itself we do not address this issue in the present work. However, at generic filling the physics is not expected to be different from the s-wave case, and the strong-coupling mapping of the effective action in the XY model should be preserved.

VII. CONCLUSIONS
In this paper we analyzed in detail the issue of the collective-mode description of a superconductor in the quantum regime. In particular, we addressed the issue from the point of view of the phase-fluctuation correction to the superfluid density, and we extended the previously known results for continuum microscopic models to the more realistic lattice case.
Usually, the evaluation of the phase-fluctuation contribution to the superfluid-density depletion is addressed within the quantum generalizations of the classical XY model, which however fails in providing a complete description of dynamical effects. Specifically, in the continuum system the finite correction to D down to T = 0 derived within the XY model implies the explicit violation of Galilean invariance. As we showed, the inadequacy of the quantum XY model can be understood as one derives the phase-only action starting from the microscopic model. This approach allows us to treat on the same footing both static and dynamic interaction terms in the phase field: besides the classical (static) anharmonic terms ( ∇θ) 4 , one has third-and fourth-order quantum (dynamic) interaction terms which contain the time derivative of θ. These quantum terms are absent in the quantum XY model, where the dynamics only appears at the Gaussian level, and induce a correction to D which cancels exactly, in the continuum case, the contribution due to the classical interaction, restoring the equality ρ s = ρ at T = 0. We point out that this evident failure of the quantum XY in the continuum case naturally poses the question of the quantitative validity of the quantum XY model estimates even in the case of a lattice microscopic model.
The detailed description of the main steps required to derive the phase-only action was motivated by the fact that this formalism also provides a more general framework for the description of different microscopic systems. As a consequence, we were able to investigate the phase mode both in the presence and in the absence of long-range Coulomb interactions, and for both s-wave and d-wave symmetry of the order parameter. As far as the lattice case is concerned, we also showed that the phase-only action acquires a more complex structure due to the fact that the minimal substitution couples the fermions to the electromagnetic field at all the orders. Moreover, since the stiffness at T = 0 is no longer related to the particle density but rather proportional to the average kinetic energy (for nearestneighbor hopping), no conservation law is violated by finding that classical and quantum interaction terms lead now to a finite one-loop correction δD. When compared with the result of the quantum XY model, one also sees that δD is of order 1/(k F ξ 0 ) 2 with respect to δD XY . A similar result also holds in the classical regime, where we found that the phase-fluctuation correction to D is smaller (in the weak-to-intermediate-coupling regime) than within the classical XY model by the same factor ∼ 1/(k F ξ 0 ) 2 , in both the continuum and the lattice model. The reduction of the phase-fluctuation effects for k F ξ 0 ≫ 1 is made even more pronounced in the continuum case by the inclusion of long-range Coulomb forces.
We also devoted particular emphasis to the discussion of the strong-coupling regime. As the interaction strength increases the fermionic model evolves toward a bosonic one, where the Cooper pairs act as individual bosons. Since in the superfluid the order parameter is connected to density, one expects that the description of the strongly-interacting electronic system requires the inclusion of both density and modulus fluctuations. We explicitly considered the negative-U Hubbard (s-wave) model. This model and its extensions are widely considered in the literature, as candidates to capture the main features of the intermediate-to-strong-coupling superconductors. This issue has attracted renewed interest due to its possible relevance to the description of the underdoped regime of the superconducting cuprates.
The inclusion of density fluctuations within the negative-U Hubbard model is performed by RPA resummation in the particle-hole channel of the fermionic bubbles which appear as coefficients of the phase-only action. As far as the modulus fluctuations are concerned, their inclusion in the fermionic bubbles is directly derived by the definition of the susceptibilities in terms of the physical quantities, although the formal treatment of modulus fluctuation at Gaussian level is possible. The discussion of the result strongly involves the knowledge of the behavior of the coherence length ξ 0 at strong coupling as a function of the fermion density, since it sets the spatial cut-off for phase fluctuations. This issue was addressed elsewhere, 15 and allows us to conclude that at strong coupling, both in the presence and in the absence of long-range Coulomb forces, the fermionic system maps onto the quantum XY model, except in the very low-density quantum regime.
Even though our attention was devoted to the correction to the superfluid density, our results for the structure of the phase-only action also suggest the possible outcomes of the analysis of phase-fluctuation effects on other quantities. The connection between the coefficients of the phase-only action and the various physical response functions (in specific regimes) allows us to relate the correction to the Gaussian phase-only action to the phase-fluctuation induced correction to the corresponding physical quantity. This has been done, e.g., for the optical conductivity, as discussed in Ref. [30] in connection to the physics of high-T c superconductors. A second possibility concerns the phase-fluctuation effect on the thermal conductivity, which is presently under investigation 31 .
where w(r) = k w(k)e ikr is the Fourier transform of w(k). After the transformation (A2), we obtain, e.g., for the first term linear in ∆ in Eq. (A3) If the interaction is local in real space, i.e., w(r) = δ(r), as in the (isotropic) s-wave case, the phase θ disappears from the exponential in Eq. (A4). In the d-wave case instead, one can expand e i[θ(R,τ )− θ(R− r 2 ,τ ) 2 − θ(R+ r 2 ,τ ) 2 ] ≃ 1 − ir α r β ∇ α ∇ β θ(R) + · · ·. The residual dependence on the phase affects the form of the phase-fluctuation propagator by introducing terms which are irrelevant in the hydrodynamical limit. As a consequence, in the following we completely discard the dependence on the phase field in S I , assuming that after the gauge transformation (A2) S I has the form (A3) with both ∆(x) and ∆ * (x) substituted by |∆(x)|.
As far as S 0 in Eq. (5) is concerned, the effect of the gauge transformation can be better expressed by first rewriting where we introduced the standard Nambu spinor φ + (r, τ ) = c + ↑ (r, τ ) c ↓ (r, τ ) . After the transformation (A2), S 0 gets modified intoS where the operator ↔ ∇≡ ( ← ∇ − → ∇) acts on the fermionic variables. Introducing the Nambu notation also for the interacting part of the action, after the gauge transformation we find where we put k = (k, ω n ) and q = (q, Ω m ), and the matrix A ij k ′ k can be deduced from Eqs. (A1) and (A5). As we observed in Sec. II A, the fermionic variables appear in the action (A6) only in the Gaussian term, whose coefficient A ij k ′ k contains the collective variables |∆| and θ. By integrating out the fermions, 32 and using the well-known identity: ln DetA ij k ′ k = Tr ln A ij k ′ k , where the trace Tr acts in the Nambu space and on the four-momenta k, k ′ , we finally get the effective action for the bosonic degrees of freedom The second term of Eq. (A7) acquires a more readable form by separating the part of the matrix A k,k ′ with the explicit structure of δ kk ′ from the rest. This leads to where the first, the third and the fourth terms were already present in Eq. (A1). The last term describes the interaction of the fermions with the ρ HS . This term is not affected by the gauge transformation (A2), since it is of the form c + c. Integrating out the fermions the effective action for the collective variables(A7) is now where the matrix elements A ij k ′ k include now the density field ρ HS , A k ′ k −→ A k ′ k + i T /Ωρ HS (k − k ′ )τ 3 . As a consequence, the density contributes to both the diagonal part of A ij k ′ k , whencê and to the self-energy matrixΣ →Σ + iρ HS (x)τ 3 . Again, we separate ρ HS (q) = ρ HS (0) + δρ HS (q). The value of ρ HS (0) is determined by the solution of the saddle-point equation for the action (A16), i.e.
− iρ HS (0) where ρ MF is the particle density evaluated at mean-field level according to Eq. (A15). Therefore, the correction (A17) to the Green function is the standard Hartree shift of the chemical potential µ → µ + V (0)ρ MF . 33 .

APPENDIX B: EFFECTIVE ACTION AND GAUGE INVARIANCE
In Sec. II A we observed that the BCS bubbles Λ(q) which appear as coefficients of the Gaussian effective action (10) break, as it is well known 18 , the gauge invariance of the theory. In this appendix we derive the general relationship between the coefficents Λ(q) and the corresponding gauge-invariant electromagnetic response functions K(q), which control the physically accessible quantities (see also Ref. [5,13]).