Local Ward identities for collective excitations in fermionic systems with spontaneously broken symmetries

We derive Ward identities for fermionic systems exhibiting a gauge symmetry that gets globally broken. In particular, we focus on the relation that connects the gauge field response functions to the transverse susceptibilities of the order parameter. We find that the long-wavelength and zero energy limit of the former are related to the coefficients of a low-energy expansion of the latter. We examine three physical cases: the superconductor, the N\'eel antiferromagnet and the spiral magnet. In the case of a metallic spiral magnet that completely breaks the SU(2) spin symmetry we explicitly show that the Ward identities are fulfilled within the random phase approximation. We subsequently derive microscopic expressions for the spin stiffnesses and spin wave velocities, which can be plugged into low energy models to study the effect of long-wavelength bosonic fluctuations on top of mean-field solutions.


I. INTRODUCTION
Symmetries play a prominent role in several branches of physics.In condensed matter systems, they often allow for making rigorous statements about the nature of a phase transition, or provide insightful information about entire phases.A renowned example is the Goldstone theorem [1], which predicts the number of gapless modes when a global continuous symmetry is spontaneously broken, irrespective of the microscopic details of the system.In some systems, a global invariance is promoted to a local one with the introduction of a gauge field.
In this case, constraints for the correlation functions can be derived, corresponding to fundamental conservation laws.These exact relations, going under the name of Ward (or Ward-Takahashi) identities [2,3], play a key role in the development of approximations to tackle the many body problem [4][5][6][7], as it is desirable for a method to satisfy them.Many works have been produced with this aim, focusing, among others, on the consistency of many-body approximations [8][9][10][11][12], or on the gauge invariance of response functions [13][14][15].Ward identities have also been analyzed in systems with broken symmetries, such as ferromagnets [9,16] or superconductors [17][18][19][20].Some of the most considered Ward identities are those which impose constraints relating single particle properties (as the self-energy) to the two-particle correlators (the vertex functions), providing a benchmark on the consistency of the approximations employed to compute these quantities In this paper, we derive Ward identities descending from a local symmetry that gets globally broken in a fermionic system.We consider the gauge field as an external perturbation and probe the linear response of the system to it, thus preserving the massless character of the Goldstone modes.We derive a functional identity from which a hierarchy of exact constraints for correlation functions can be extracted.We focus on the Ward identity connecting the linear response to the external gauge field and the transverse susceptibilities of the or-der parameter, containing information on the Goldstone gapless excitations.In particular, we show that the low frequency and momentum expansion of the inverse of the susceptibilities is related to precise limits of the gauge field response kernel.We discuss three examples: a superconductor breaking the U(1) charge symmetry, a Néel antiferromagnet breaking the SU(2) spin symmetry leaving a residual U (1), and a spiral magnet, completely breaking SU (2).Similar relations can be found in the literature for special cases (see for example Ref. [21]).Here, we provide a formal derivation within a formalism that allows to extend them to any system and any possible kind of long-range order that breaks a continuous symmetry.These identities are of fundamental relevance, as from the sole knowledge of the local symmetry group, one can infer the form of the transverse susceptibilities for small energies and long wavelengths.This is reminiscent of and intimately connected to hydrodynamic theories for systems exhibiting spontaneous symmetry breaking, such as superfluid helium [22] or magnets [23,24].Moreover, when the gauge field is of "pure gauge", that is, it can be removed by means of a local transformation, the Ward identities are useful to derive a low-energy theory for the slow directional fluctuations of the order parameter in the fashion of nonlinear sigma models [25].
We focus particularly on the case of doped antiferromagnets, that host low-lying fermionic excitations, forming small Fermi surfaces.We show that from the Ward identities one can compute a dynamical transverse susceptibility χ ⊥ dyn ≡ lim ω→0 χ ⊥ (q = 0, ω) and a spin stiffness J and extract the spin wave velocity as c s = J/χ ⊥ dyn .This is in contrast with the hydrodynamic approach of Refs.[23,24], which predicts c s = J/χ ⊥ , with χ ⊥ ≡ lim q→0 χ ⊥ (q, ω = 0) the uniform transverse susceptibility.As it was noticed in Ref. [26], the equality χ ⊥ = χ ⊥ dyn holds only for insulating antiferromagnets at low temperatures, as the spin systems considered in Refs.[23,24].Furthermore, as a consequence of the doping, Néel antiferromagnetism is often destroyed, making way to a coplanar spiral magnet.In this state, the spins are no longer antiferromagnetically ordered but the magnetization rotates in a plane, completely breaking the SU(2) rotational symmetry and giving rise to three rather than two Goldstone modes [27][28][29][30].Spiral states have been found to emerge in the two dimensional Hubbard model and in the t-J model at moderate hole doping away from half filling, and they are possible candidates for the incommensurate magnetic correlations observed in the cuprate superconductors [31][32][33][34][35][36][37][38][39][40].They have also been proposed in the context of frustrated antiferromagnets [41][42][43][44][45][46][47].We present a direct evaluation of the Ward identities for a metallic spiral magnet via the random phase approximation, and derive expressions for the spin stiffness and dynamical susceptibility of each gapless mode.These formulas generalize to finite temperature and doping and to spiral magnetic states the previously derived expressions for spin stiffnesses in the ground state at half filling [21,48].This formalism represents a convenient starting point to study longwavelength fluctuation effects on top of mean-field solutions.
The paper is organized as follows.In Sec.II we derive the anomalous Ward identities descending from a gauge symmetry that gets globally broken in a fermionic system.We consider three different systems: the superconductor, the antiferromagnet, and the spiral magnet.In Sec.III we present an explicit and detailed evaluation of the Ward identities for a spiral magnet and derive approximate expressions for the spin stiffnesses and dynamical susceptibilities.A conclusion in Sec.IV closes the presentation.

II. WARD IDENTITIES
In this section we derive and discuss the Ward identities connected with a specific gauge symmetry which gets globally broken due to the onset of long-range order in the fermionic system.We focus on two specific symmetry groups: the (abelian) U(1) charge symmetry and the (nonabelian) SU(2) spin symmetry.All over the paper we employ Einstein's notation, that is, a sum over repeated indices is implicit.

A. U(1) symmetry
We now analyze the Ward identities in a superconductor.According to the Elitzur's theorem [49], the breaking of the U(1) gauge symmetry is not possible, leading to the interpretation of the superconductor as a topologically ordered state [50], due to a subtle interplay of the dynamical gauge field and the order parameter [51].In this work, we treat the gauge field as a classical external perturbation and measure the system's response to it, so that a discussion on the topological nature of the order is not required.
We consider the generating functional of susceptibilities of the superconducting order parameter and gauge kernels, defined as: where ψ = (ψ ↑ , ψ ↓ ) (ψ = (ψ ↑ , ψ ↓ )) are Grassmann spinor fields corresponding to the annihilation (creation) of a fermion, A µ is the electromagnetic field, J (J * ) is a source field that couples to the superconducting order parameter ψ ↑ ψ ↓ (ψ ↓ ψ ↑ ), and S[ψ, ψ, A µ ] is the action of the system.The index µ = 0, 1, . . ., d, with d the system dimensionality, runs over temporal (µ = 0) and spatial (µ = 1, . . ., d) components.In the above equation and from now on, the expression (A, B) has to be intended as x A(x)B(x), where x is a collective variable consisting of a spatial coordinate x (possibly discrete, for lattice systems), and an imaginary time coordinate τ , and x is a shorthand for d d x β 0 dτ , with β the inverse temperature.Even in the case of a lattice system, we define the gauge field over a continuous space time, so that expressions involving its gradients are well defined.
We let the global U(1) charge symmetry be broken by an order parameter that, to make the treatment simpler, we assume to be local (s-wave) where the average is computed at zero source and gauge fields, and, without loss of generality, we choose φ 0 ∈ R.
A generalization to systems with nonlocal order parameters, such as d-wave superconductors, is straightforward.The functional (1) has been defined such that its second derivative with respect to J and J * at zero J, J * and A µ gives (minus) the susceptibility of the order parameter χ(x, x ′ ), while (minus) the gauge kernel K µν (x, x ′ ) can be extracted differentiating twice with respect to the gauge field.In formulas Let us now consider the constraints that the U(1) gauge invariance imposes on the functional G.Its action on the fermionic fields is with θ(x) a generic function.Similarly, the external fields transform as where ∂ µ = (i∂ τ , ∇).In Eqs. ( 4) and ( 5) the spatial coordinate x of the spinors ψ and ψ, as well as the sources J and J * may be a lattice one, while the gauge field A µ and the parameter θ are always defined over a continuous space.To keep the notation lighter, we always indicate the space-time coordinate as x, keeping in mind that its spatial component could have a different meaning depending on the field it refers to.For G to be invariant under a U(1) gauge transformation, it must not to depend on θ(x): Considering an infinitesimal transformation, that is |θ(x)| ≪ 1, from Eqs. ( 5) and ( 6), we obtain We now consider the change of variables ) is a source field coupling to longitudinal (transverse) fluctuations of the order parameter, and the functional Γ, defined as the Legendre transform of G, where . The gauge kernel can be computed from Γ as well: because, thanks to the Legendre transform properties, δΓ/δA µ (x) = δG/δA µ (x).Differently, differentiating Γ twice with respect to the fields ϕ a returns the inverse correlator which obeys a reciprocity relation [52] x ′′ with the generalized susceptibility χ ab (x, x ′ ), defined as Eq. ( 7) can be expressed in terms of Γ as Eq. ( 14) is an identity for the generating functional Γ stemming from U(1) gauge invariance of the theory.Taking derivatives with respect to the fields, one can derive an infinite set of Ward identities.We are interested in the relation between the gauge kernel and the transverse inverse susceptibility C 22 (x, x ′ ).For this purpose, we differentiate Eq. ( 14) once with respect to ϕ 2 (x ′ ) and once with respect to A ν (x ′ ), and then set the fields to zero.We obtain the set of equations where ⟩, and we have defined the quantity Combining (15a) and (15b), we obtain Fourier transforming Eq. ( 17) and rotating to real frequencies, we have with q = (q, ω) a collective variable combining momentum and real frequency.We now define the superfluid stiffness J αβ and the uniform density-density susceptibility χ n [53] as where the minus sign in (19a) has been introduced so that J αβ is positive definite.Notice that, even though the limits q → 0 and ω → 0 in Eq. (19b) have been taken in the opposite order compared to what is conventionally done, they commute in a s-wave superconductor because of the absence of gapless fermionic excitations.
In the above equation and from now on, we employ the convention that the indices labeled as µ, ν include temporal and spatial components, whereas α and β only the latter.Taking the second derivative with respect to q on both sides of (18), we obtain where ∂ 2 qαq β and ∂ 2 ω are shorthands for ∂ 2 ∂qαq β and ∂ 2 ∂ω 2 , respectively.Moreover, we have made use of the Goldstone theorem, reading C 22 (0, 0) = 0. To derive Eq. ( 20) from (18) we have exploited the finiteness of the gauge kernel K µν (q) in the q → 0 and ω → 0 limits.Eq. (20) states that the superfluid stiffness and the uniform densitydensity correlation function are not only the zero momentum and frequency limit of the gauge kernel, but also the coefficients of the inverse transverse susceptibility when expanded for small q and ω, respectively.Inverting Eq. ( 12), C 22 (q) can be expressed in terms of χ ab (q) as In the limit q → 0 = (0, 0), χ 22 (q) diverges for the Goldstone theorem, while the second term in the denominator vanishes like some power of q.This implies that, for small q, From this consideration, together with (20), we can deduce that the transverse susceptibility can be written as for small q and ω.
The above form of the χ 22 (q) can be also deduced from a low energy theory for the phase fluctuations of the superconducting order parameter.Setting J and J * to zero in (1), and integrating out the Grassmann fields, one obtains an effective action for the gauge fields.The quadratic contribution in A µ is where q is a shorthand for dω 2π d d q (2π) d .Since we are focusing only on slow and long-wavelength fluctuations of A µ , we replace K µν (q) with K µν (0).Considering a pure gauge field, A µ (x) = −∂ µ θ(x), where θ(x) is (half) the phase of the superconducting order parameter (ϕ(x) = φ 0 e −2iθ(x) ), we obtain with θ(x) ∈ [0, 2π] a periodic field.The above action is well known to display a Berezinskii-Kosterlitz-Thouless (BKT) transition [54,55] for d = 1 (at T = 0) and d = 2 (at T > 0), while for d = 3 (T ≥ 0) or d = 2 (T = 0), it describes a gapless phase mode known as Anderson-Bogoliubov phonon [19].

B. SU(2) symmetry
In this Section, we repeat the same procedure we have applied in the previous one to derive the Ward identities connected to a SU(2) gauge invariant system.We consider the functional where A µ (x) = A a µ (x) σ a 2 is a SU(2) gauge field, ⃗ σ are the Pauli matrices, and ⃗ J(x) is a source field coupled to the fermion spin operator 1  2 ψ(x)⃗ σψ(x).Similarly to the previous section, derivatives of G with respect to A µ and ⃗ J at zero external fields give minus the gauge kernels and spin susceptibilities, respectively.In formulas, We let the SU(2) symmetry be broken by a (local) order parameter of the form with v(x) a position-dependent unit vector pointing along the local direction of the magnetization [56].
A SU(2) gauge transformation on the fermionic fields reads where R(x) ∈ SU( 2) is a matrix acting on the spin indices of ψ and ψ.The external fields transform as where R(x) is the adjoint representation of R(x) The SU(2) gauge invariance of G can be expressed as Writing R(x) = e iθa(x) σ a 2 , R † (x) = e −iθa(x) σ a 2 , and considering an infinitesimal transformation |θ a (x)| ≪ 1, we obtain the functional identity where ε abc is the Levi-Civita tensor.Γ[A µ , ⃗ ϕ] is the Legendre transform of G, defined as with δJa(x) .The inverse susceptibilities C ab (x, x ′ ), defined as, obey a reciprocity relation with the spin susceptibilities χ ab (x, x ′ ) similar to (12).Defining the quantities we obtain from (35) the set of equations where (39a), (39b), have been obtained differentiating (35) with respect to ϕ b (x ′ ) and A ν (x ′ ), respectively, and setting the fields to zero.Eq. (39c) simply comes from (35) computed at zero A µ , ϕ a .According to Eq. ( 30), the expectation value of ⃗ ϕ(x) takes the form ⟨ ⃗ ϕ(x)⟩ = φ 0 v(x).Combining (39a), (39b), and (39c), we obtain the Ward identity which connects the gauge kernels with the inverse susceptibilities.
In the following, we analyze two concrete examples where the above identity applies, namely the Néel antiferromagnet and the spiral magnet.We do not consider ferromagnets or, in general, systems with a net average magnetization, as in this case the divergence of the transverse components of the kernel K ab 00 (q) for q → 0 leads to changes in the form of the Ward identities.In this case, one can talk of type-II Goldstone bosons [57], characterized by a nonlinear dispersion.

Néel order
We now consider the particular case of antiferromagnetic (or Néel) ordering for a system on a d-dimensional bipartite lattice.In this case v(x) takes the form (−1) x v, with (−1) x being 1 (−1) on the sites of sublattice A (B), and v a constant unit vector.In the following, without loss of generality, we consider v = (1, 0, 0).Considering only the diagonal (a = b) components of (40), we have Despite Néel antiferromagnetism breaking the lattice translational symmetry, the components of the gauge Kernel considered above depend only on the difference of their arguments x − x ′ , and thus have a well-defined Fourier transform.Eq. (41a) implies q µ q ν K µν 11 (q, ω) = 0, as expected due to the residual U(1) gauge invariance in the Néel state.In particular, one obtains lim q→0 K 11 αβ (q, 0) = 0, and lim ω→0 K 11 00 (0, ω) = 0. Eqs.(41b) and (41c) are the same equation as we have ), again because of the residual symmetry.If we rotate them onto the real time axis and perform the Fourier transform, we get where J αβ is the spin stiffness, Q = (π/a 0 , . . ., π/a 0 ), with a 0 the lattice spacing, and we name χ ⊥ dyn as transverse dynamical susceptibility [53].In the above equations we have made use of the Goldstone theorem, which in the present case reads Furthermore, to derive Eq. ( 42) from (41b), we have used the finiteness of the q → 0 and ω → 0 limits of the gauge kernels.Following the argument given in the previous section, for q = (q, ω) close to Q = (Q, 0), we can replace C 33 (q) by 1/χ 33 (q) in ( 42), implying Notice that in Eq. ( 44) we have neglected the imaginary parts of the susceptibilities, that, for doped antiferromagnets, can lead to Landau damping of the Goldstone modes [58].Also for Néel ordering, form (44) of the transverse susceptibilities can be deduced from a low energy theory for the gauge field A µ (x), that is, Considering a pure gauge field with R(x) a SU(2) matrix, we obtain the action where n(x) = (−1) x R(x)v(x), with R(x) defined as in Eq. ( 33), and |n(x)| 2 = 1.Eq. ( 47) is the well-known O(3)/O(2) nonlinear sigma model (NLσM) action, describing low-energy properties of quantum antiferromagnets [25,59].
Writing R(x) = e iθa(x) σ a 2 , and expanding to first order in θ a (x), n(x) becomes n(x) ≃ (1, θ 2 (x), −θ 3 (x)).Considering the expression ⃗ ϕ(x) = (−1) x φ 0 n(x) for the order parameter field, we see that small fluctuations in n(x) only affect the 2-and 3-components of ⃗ ϕ(x).The transverse susceptibilities can be therefore written as which is the result of Eq. ( 44).In Eq. ( 48) we have made use of the propagator of the n-field dictated by the action of Eq. ( 47), that is, Eq. ( 48) predicts two degenerate magnon branches with linear dispersion for small q − Q.In the case of an isotropic antiferromagnet (J αβ = Jδ αβ ), we have ω q = c s |q|, with the spin wave velocity given by c s = J/χ ⊥ dyn .

Spiral magnetic order
We now turn our attention to the case of spin spiral ordering, described by the magnetization direction where at least one component of Q is neither 0 nor π/a 0 .
In this case, it is convenient to rotate the field ⃗ ϕ(x) to a basis in which v(x) is uniform.This is achieved by the transformation [30] with In this way, the inverse susceptibilities are transformed into If we now apply the Ward identity (40), we obtain with We remark that an order parameter of the type (50) completely breaks the SU(2) spin symmetry, which is why none of the right hand sides of the equations above vanishes.We have considered a coplanar spiral magnetic order, that is, we have assumed all the spins to lie in the same plane, so that out of the three Goldstone modes, two are degenerate and correspond to out-of-plane fluctuations, and one to in-plane fluctuations of the spins.Furthermore, translational invariance is broken, so the Fourier transforms of the gauge kernels K ab µν (q, q ′ , ω) and inverse susceptibilities C ab (q, q ′ , ω) are nonzero not only for q − q ′ = 0 but also for q − q ′ = ±Q or ±2Q.Time translation invariance is preserved, and the gauge kernels and the inverse susceptibilities depend on one single frequency.However, in the basis obtained with transformation (52), translational invariance is restored, so that the Fourier transform of C ab (x, x ′ ) only depends on one spatial momentum.With this in mind, we can extract expressions for the spin stiffnesses and dynamical susceptibilities from (54).After rotating to real frequencies, and using the property that for a spiral magnet the gauge kernels are finite in the limits q = q ′ → 0 and ω → 0, we and where the labels ⊥ and □ denote out-of-plane and inplane quantities, respectively.In the equations above, we have defined K ab µν (q, ω) as the prefactors of the components of the gauge kernels K ab µν (q, q ′ , ω) which are proportional to (2π) d δ(q − q ′ ).From Eqs. (56) and (57) it immediately follows that J ⊥,1 αβ = J ⊥,2 αβ ≡ J ⊥ αβ , and χ ⊥,1 dyn = χ ⊥,2 dyn ≡ χ ⊥ dyn , as expected in the case of coplanar order [41].To derive the equations above, we have made use of the Goldstone theorem, which for spiral ordering reads (see for example Refs.[30,34]) Notice that the above relations can be also derived from a functional identity similar to (35) but descending from the global SU(2) symmetry.Moreover, close to their respective Goldstone points ((0, 0) for C 22 , and (±Q, 0) for C 33 ), C 22 (q) can be replaced with 1/ χ 22 (q), and C 33 (q) with 1/ χ 33 (q), with the rotated susceptibilities defined analogously to (53).If the spin spiral state occurs on a lattice that preserves parity, we have C aa (q, ω) = C aa (−q, ω), from which we obtain Neglecting the imaginary parts of the susceptibilities, giving rise to dampings of the Goldstone modes [58], from Eq. ( 59) we can obtain expressions for the susceptibilities near their Goldstone points Expressions ( 60) can be deduced from a low energy model also in the case of spin spiral ordering.Similarly to what we have done for the Néel case, we consider a pure gauge field, giving the nonlinear sigma model action where R(x) ∈ SO(3) is defined as in Eq. (33), and now ∂ µ denotes (−∂ t , ⃗ ∇).The matrix P µν is given by with for a ∈ {□, ⊥}.Action ( 61) is a NLσM describing low energy fluctuations around a spiral magnetic ordered state.It has been introduced and studied in the context of frustrated antiferromagnets [41][42][43][44].
We now write the field ⃗ we get with ê1 = (1, 0, 0), and ⃗ θ ′ (x) = M(x) ⃗ θ(x).Inserting (64) into (61), we obtain We are finally in the position to extract the form of the susceptibilities for small fluctuations which is the result of Eq. ( 60).In the above equations we have used the correlators of the θ field descending from action (66).Form (60) of the susceptibilities predicts three linearly dispersing Goldstone modes, two of which (the out-of-plane ones) are degenerate and propagate with velocities c ⊥ are the eigenvalues of J ⊥ αβ and n = 1, . . ., d.Similarly, the in-plane mode velocity is given by c

III. EXPLICIT CALCULATION FOR A SPIRAL MAGNET
In this section, we present an explicit calculation for a spiral magnet, emerging from a Hubbard model, of the left and right hand sides of ( 56) and (57), and show that the Ward identities are fulfilled.We compute the spin susceptibilities and the gauge kernels within the random phase approximation (RPA), which is a conserving approximation in the sense of Baym and Kadanoff [4], and it is therefore expected to fulfill the Ward identities.

A. Mean-field approximation
In a fermionic lattice system, the amplitude φ 0 of spin spiral order parameter can be expressed as where k is a shorthand for an integral over the full d-dimensional Brillouin zone, that is k∈BZ 69), we deduce that spiral ordering couples only the electron states (k, ↑) and (k + Q, ↓), and the meanfield fermion Green's function can be expressed as a 2 × 2 matrix: where ν n = (2n + 1)πT (n ∈ Z) is a fermionic Matsubara frequency, ξ k = ϵ k − µ, with ϵ k the single particle dispersion and µ the chemical potential, and ∆ is the spiral order parameter.Diagonalizing (70), one obtains the quasiparticle dispersions with g k = (ξ k + ξ k+Q )/2, and The Green's function can be conveniently written as a linear combination of quasiparticle poles where the coefficients u ℓ k are given by with σ 0 = 1, and e k = h 2 k + ∆ 2 .We assume the spiral states to emerge from a lattice model with onsite repulsive interactions (Hubbard model), with imaginary time action where t jj ′ describes the hopping amplitude between the lattice sites labeled by j and j ′ and U is the Hubbard interaction.The Hartree-Fock or mean-field (MF) gap equation at temperature T reads where νn denotes a sum over the fermionic Matsubara frequencies, and f (x) = 1/(e x/T + 1) is the Fermi function.The bosonic order parameter φ 0 is related to ∆ as ∆ = U φ 0 .Finally, the optimal Q-vector is obtained minimizing the mean-field free energy where n is the fermion density.

B. Random phase approximation
In this Section, we summarize the RPA for a spiral magnet, closely following Ref.[58].Further details can be also found in Ref. [30].
The charge and spin susceptibilities are coupled in a spin spiral state.It is therefore convenient to consider the retarded real-time correlators [58] where j and j ′ label two lattice sites, t is the time, [•, •] denotes the commutator, Θ(t) is the Heaviside function, and a and b run from 0 to 3. The generalized spin operators are defined in terms of the fermionic fields as The a = 1, 2, 3 components correspond to the usual spin operator, while a = 0 gives half the fermion density n j = ψ j,↑ ψ j,↑ + ψ j,↓ ψ j,↓ .
It is convenient to work in the rotated basis [30] introduced in Eqs. ( 51) and ( 52), the susceptibilities transform as with where r j represents the coordinates of the lattice site j.
As previously mentioned, the χ jj ′ (t) are translationally invariant, so that their Fourier transform depends only on a single spatial momentum q.Within RPA, the susceptibilities are given by where χ 0 (q) is a matrix in the a, b = 0, 1, 2, 3 indices, 1 is the 4 × 4 identity matrix, and Γ 0 = 2U diag(−1, 1, 1, 1).The bare susceptibilities on the real frequency axis are given by where Ω m = 2mπT (m ∈ Z) denotes a bosonic Matsubara frequency and the substitution iΩ m → ω + i0 +  I. Symmetries of the bare susceptibilities.The first two signs in each cell represent the sign change of χ ab 0 (q) by flipping the sign of q and ω, respectively.The third sign corresponds to the sign change of χ ab 0 (q) when a and b are interchanged.
has been performed.Using (72), one can perform the frequency sum, obtaining where we have defined with P the principal value, and the coherence factors Notice that in Eq. ( 84) we have neglected the imaginary part of 1/(ω , which is proportional to a δ function, and is responsible for the Landau damping of the collective spin fluctuations [58]. The bare susceptibilities have a well defined sign under inversion of q and ω, and under the exchange of the indices a and b.These symmetries can be derived from Eqs. (83), (84), and (85) (see, for example, Ref. [58]) and are listed in Table I.Since neglecting the imaginary part of F ℓℓ ′ (k, q, ω) makes the matrix χ ab 0 (q) hermitian, we deduce that those elements which are even (odd) under the exchange a ↔ b are purely real (imaginary).At least within the RPA, the symmetry properties derived for the bare χ 0 (q) hold also for the full susceptibilities χ(q).
Goldstone modes associated with spontaneous symmetry breaking of the SU(2) spin symmetry appear in the susceptibilities.In particular, considering the limit q → 0, all the off-diagonal elements of the type χ 2a 0 (q) and χ a2 0 (q) vanish as they are odd in q or ω, while χ 22 0 (q) takes the simple form The full susceptibility at q = 0 therefore reads We immediately see that when the gap equation ( 75) is fulfilled, the denominator of (87) vanishes, signaling a gapless mode.This Goldstone mode represents fluctuations of the spins within the plane in which the magnetization lies.We therefore refer to it as in-plane mode in the following.With a similar reasoning, one proves that in the limit q → ±Q (with Q = (Q, 0)), all the off diagonal elements of the type χ 3a 0 (q) and χ a3 0 (q) vanish, leaving Noticing that χ 33 0 (±Q) = χ 22 0 (0) [58], one sees that the denominator of (88) is vanishing.The susceptibility χ 33 (q) therefore hosts two poles at q = ±Q, which are due to gapless excitations of the spins outside the plane in which the magnetization lies.We refer to these modes as out-of-plane modes.As already stated, in the spiral state three Goldstone modes appear, as all three generators of the symmetry group SU(2) are broken.However, the two out-of-plane modes have degenerate dispersions, as it is expected for coplanar order.

C. Small q and ω expansion of the susceptibilities
We are now in the position to compute the right hand sides of Eqs. ( 56) and (57).The small q and ω expansion of the susceptibilities in the spiral state has already carried out in Ref. [58], but for completeness we revisit it here in a slightly different form.

In-plane mode
Using (81), the in-plane susceptibility can be conveniently written as with where Γ ab 0,2 and χ ab 0,2 (q) are matrices obtained from Γ ab 0 and χ ab 0 (q) removing the components where a = 2 and/or b = 2, and 1 3 denotes the 3 × 3 identity matrix.For later convenience, we notice that for q = 0, all the off-diagonal elements χ 2a 0 (q) and χ a2 0 (q) vanish, so that Γ 2 (0) can be obtained from the full expression selecting only the components in which the indices take the values 0,1, or 3.

D. Gauge kernels
To calculate the gauge kernels, that is, the left hand sides of Eqs. ( 56) and ( 57), we couple our system to a SU(2) gauge field via a Peierls substitution in the quadratic part of action (74): where e −r jj ′ •∇ is the translation operator from site j to site j ′ , with r jj ′ = r j − r j ′ .Notice that under the transformation ψ j → R j ψ j , with R j ∈ SU(2), the interacting part of the action (74) is left unchanged, while the gauge field transforms according to (32b).Since the gauge kernels correspond to correlators of two gauge fields, we expand (101) to second order in A µ .After a Fourier transformation one obtains where the first order coupling is given by γ µ k = (1, ∇ k ϵ k ), and the second order one is γ αβ k = ∂ 2 kαk β ϵ k .Analyzing the coupling of the temporal component of the gauge field to the fermions in (101) and (102), we notice that the temporal components of the gauge kernel are nothing but the susceptibilities in the original (unrotated) spin basis K ab 00 (q, q ′ , ω) = χ ab (q, q ′ , ω), where ω is a real frequency.The spatial components of the gauge kernel can be expressed in the general form (see Fig. 1) K ac para,α0 (q, q ′′ , ω)× ×Γ cd (q ′′ , q ′′′ , ω)K db para,0β (q ′′′ , q ′ , ω), where Γ(q ′ , q ′′ , ω) is the effective interaction (92) expressed in the unrotated basis.Within the RPA, the paramagnetic terms are given by The Green's function in the unrotated basis takes the form where δ k,k ′ is a shorthand for (2π) d δ d (k − k ′ ), and The diamagnetic term does not depend on q, q ′ and ω, and is proportional to the unit matrix in the gauge indices.It evaluates to We can now compute the spin stiffnesses and dynamical susceptibilities from the gauge kernels.

In-plane mode
The in-plane spin stiffness is defined as  2) gauge field, the solid lines the electronic Green's functions, the black triangles the paramagnetic vertex γ µ k σ a , the black circle the diamagnetic one γ αβ k σ 0 , and the dashed line the effective interaction Γ(q, q ′ , ω).
The dynamical susceptibility is defined as From Eq. ( 80) we deduce that Remarking that for ω = 0 all the off-diagonal elements of the bare susceptibilities with one (and only one) of the two indices equal to 3 vanish, we obtain the RPA expression for χ □ dyn

Out-of-plane modes
To compute the the out-of-plane stiffness, that is, we find that all the paramagnetic contributions to the gauge kernel that mix temporal and spatial components vanish in the ω → 0 and q = q ′ → 0 [60] limits.Moreover, the q → 0 limit of the momentum diagonal paramagnetic contribution can be written as where we have defined σ ± = (σ 1 ± iσ 2 )/2.The out-ofplane spin stiffness is thus given by Finally, we evaluate the dynamical susceptibility of the out-of-plane modes.This is defined as Applying transformation (79), we can express the momentum-diagonal component of χ 22 (q, q ′ , ω) in terms of the susceptibilities in the rotated basis as ζ=± χ 11 (q + ζQ, ω) where we have used (see Sec. III B) χ 12 (q) = − χ 21 (q).Sending q to 0 in (120), and using the symmetry properties of the susceptibilities for q → −q (see Table I), we obtain with χ −+ (q) = ⟨S − (−q)S + (q)⟩, and S ± (q) = (S 1 (q) ± iS 2 (q))/2.It is convenient to express χ −+ (Q, ω) as where we have defined In the limit ω → 0, χ −3 0 (Q, ω) and χ 3+ 0 (Q, ω) vanish as they are odd in frequency (see Table I).We can now cast the dynamical susceptibility in the form or, equivalently, We remark that in the formulas above we have not specified in which order the limits q → ±Q and ω → 0 have to be taken as they commute.

E. Equivalence of RPA and gauge theory approaches
In this Section, we finally prove that the expressions for the spin stiffnesses and dynamical susceptibilities obtained in Sec.III C coincide with those of Sec.III D via a direct evaluation.

In-plane mode
We start by computing the first term in Eq. ( 93).The second derivative of the 22-component of the bare susceptibility can be expressed as where f ′ (x) = df /dx is the derivative of the Fermi function.On the other hand, the bare contribution to J □ αβ (Eq.( 110)) reads The second (diamagnetic) term can be integrated by parts, giving where we have used the properties Summing up both terms, we obtain Performing the Matsubara sum, we arrive at which is the same result as in (126).Furthermore, one can show that Inserting results (126), (130), and (132) into ( 93) and (112), we prove that these two expressions give the same result for the in-plane stiffness.Explicit expressions for κ 30 α (0) and κ 31 α (0) are given in Appendix A. If we now consider the dynamical susceptibility, it is straightforward to see that which, if inserted into Eqs.( 94) and (115), proves that the calculations of χ □ dyn via gauge kernels and via the low-energy expansion of the susceptibilities provide the same result.

Out-of-plane modes
With the help of some lengthy algebra, one can compute the second momentum derivative of the bare susceptibility χ 33 0 (q), obtaining Similarly to what we have done for the in-plane mode, we integrate by parts the diamagnetic contribution to the gauge kernel.Its sum with the paramagnetic one gives Performing the Matsubara sums, one can prove the equivalence of the RPA and gauge theory approach for the calculation of J ⊥ αβ .Similarly, we obtain for the second frequency derivative of the bubble χ 33 0 (q) Furthermore, one can prove that for a = 0, 1, 2. Inserting results (136) and (137) into Eqs.( 89) and (124), one sees that the RPA and gauge theory approaches are equivalent for the calculation of χ ⊥ dyn .In Appendix B we provide explicit expressions for the off diagonal bare susceptibilities χ −a 0 (Q).

Remarks on more general models
We remark that in the more general case of an interaction of the type producing, in general, a k-dependent gap, the identities we have proven above do not hold anymore within the RPA, as additional terms in the derivative of the inverse susceptibilities emerge, containing expressions involving first and second derivatives of the gap with respect to the spatial momentum and/or frequency.In fact, in the case of nonlocal interactions, gauge invariance requires additional couplings to the gauge field in S int , complicating our expressions for the gauge kernels.Similarly, even for action (74), approximations beyond the RPA produce in general a k-dependent ∆, and vertex corrections in the kernels are required to obtain the same result as the one obtained expanding the susceptibilities.
From these considerations, we obtain for the spin stiffness which implies that J αβ is given by Eq. ( 131).If the underlying lattice is C 4 -symmetric, the spin stiffness is isotropic in the Néel state, that is, J αβ = Jδ αβ .Similarly, for the dynamical susceptibility, we have which, combined with (133), implies with χ 33 0 (0, ω → 0) given by Eq. ( 133).We notice that the dynamical susceptibility is obtained from the susceptibility by letting q → 0 before ω → 0. This order of the limits removes the intraband terms (that is, the ℓ = ℓ ′ terms in Eq. ( 83)), which instead would yield a finite contribution to the uniform transverse susceptibility χ ⊥ ≡ lim q→0 χ 22 (q, 0).In the special case of an insulator at low temperature T ≪ ∆, the intraband contributions vanish and one has the identity χ ⊥ dyn = χ ⊥ , leading to the hydrodynamic relation for the spin wave velocity [23] c s = J/χ ⊥ (in an isotropic antiferromagnet).As noticed in Ref. [26], in a doped antiferromagnet this hydrodynamic expression does not hold anymore, and one has to replace the uniform transverse susceptibility with the dynamical susceptibility.Since J = 0 and χ ⊥ dyn = 0 in the symmetric phase due to SU(2) gauge invariance, the expression c s = J/χ ⊥ dyn yields a finite value c s at the critical point ∆ → 0, provided that J and χ ⊥ dyn scale to zero with the same power of ∆, as it happens within mean-field theory.Note that in the symmetric phase SU(2) gauge invariance does not pose any constraint on χ ⊥ , which is generally finite.
In the simpler case of perfect nesting, that is, when ξ k = −ξ k+Q , corresponding to the half-filled particle-hole symmetric Hubbard model, and at zero temperature, expressions for J and χ ⊥ have been derived in Refs.[21,48] for two spatial dimensions, and it is straightforward to check that our results reduce to these in this limit.Moreover, Eqs.31-34 in Ref. [21] are similar to our Ward identities but no derivation is provided.

IV. CONCLUSION
In conclusion, we have derived Ward identities for fermionic systems in which a gauge symmetry is globally broken.In particular, we have shown that the zeroenergy and long-wavelength components of the gauge kernels are connected to the transverse susceptibilities of the order parameter by exact relations.We have analyzed several examples, namely a superconductor, a Néel antiferromagnet, and a spiral magnet.In the latter case, we have performed an explicit calculation of the transverse susceptibilities and of the gauge kernels within the random phase approximation and verified that the Ward identities are indeed fulfilled.Furthermore, we have considered the Néel limit Q → (π/a 0 , . . ., π/a 0 ) and found that our RPA expressions for the spin stiffnesses and susceptibilities reduce to those previously obtained in Refs.[21,48].We have also shown that the hydrodynamic expression for the magnon velocity c s = J/χ ⊥ does not hold in presence of gapless fermionic excitations, and it must be replaced by c s = J/χ ⊥ dyn .While χ ⊥ is computed by taking the ω → 0 limit before letting q → 0 in the transverse susceptibility, χ ⊥ dyn is obtained reversing the order of the limits.The equality χ ⊥ dyn = χ ⊥ holds only for insulating antiferromagnets at low temperatures T ≪ ∆, as it is the case for the spin systems considered in Refs.[23,24].
Our findings find immediate application in the study of large distance properties of interacting fermion systems.In fact, via a microscopic calculation, one can evaluate the coefficients of the effective actions listed in Sec.II, and study phenomena which are purely fluctuationsdriven.Examples are the Berezinskii-Kosterlitz-Thouless transition in two-dimensional superconductors (see, for example, the microscopic calculation of Ref. [61]) or strong magnetic fluctuations leading to the formation of a pseudogap in low-dimensional electron systems [48,62].Special care, however, must be taken when the fermions form Fermi surfaces, as the low-energy modes of the order parameter can decay in particle-hole pairs, producing a Landau damping [58] that can affect the macroscopic and critical properties of the system [26,63].An interesting extension of the present work would be to derive Ward identities for the damping of the Goldstone modes.
Appendix B: Expressions for χ −a 0 (Q) We report here the RPA expressions for the off-diagonal bare susceptibilities χ −a 0 (Q), with a = 0, 1, 2. They can all be obtained by computing the trace and the Matsubara summation in Eq. (B).We obtain with F ℓℓ ′ (k, q, ω) defined as in Eq. (84).
We found a mistake in the original paper, namely the statement that the second derivative of the Γ-functional with respect to the gauge field returns the gauge kernel (Eq.( 10) of the original paper).The Ward identities (WIs) presented in the original paper are nonetheless correct, provided that one does not interpret the object labeled as K µν as the gauge kernel but as the second derivative with respect to the gauge field of the generating functional Γ computed at zero fields.The physical meaning, the usefulness, and the computation of such object are left for future research.
In the following, we present a revised derivation of the WIs for the actual gauge kernel K µν , both for the case of a U(1) and of an SU(2) symmetry.

I. U(1) SYMMETRY
We consider a fermionic theory, described by the Grassmann fields ψ and ψ and coupled to a U(1) gauge field A µ , and governed by the action S[ψ, ψ, A µ ].We define the U(1) gauge kernel as with j µ (x) ≡ δS δAµ(x) | Aµ=0 the paramagnetic current operator, P µν (x, x ′ ) ≡ − δ 2 S δAµ(x)δAν (x ′ ) | Aµ=0 the diamagnetic current operator, and where the (connected) average is taken by means of the path integral defined by S. The gauge kernel can also be defined as the second derivative at zero source fields of the G-generating functional of the original paper (Eq.( 1)) The U(1) gauge kernel obeys the relation even in presence of spontaneous symmetry breaking.This can be readily derived from Eq. ( 7) of the original paper.The above relation is somewhat obvious as it implies charge conservation, which must be fulfilled even in the superconducting state [17].Fourier transforming Eq. ( 3), one obtains where α is only a spatial index and a sum over it is implied.Here and henceforth we use Einstein's convention for repeated indices, unless otherwise specified.Note that for the derivative in Eq. ( 3) to make sense, A µ needs to be defined on a continuous space-time, which implies, in case of a lattice system, that Eq. ( 4) is valid only for those values of |q| that are much smaller then the inverse lattice spacing.Setting ν = 0, q = 0, and ν = β (with β also a spatial index), ω = 0, respectively, one gets where J αβ (q) a function approaching the superfluid stiffness for q → 0 in the superconducting state.Eq. (5b) is the most general form that K αβ (q, 0) can take to obey q α K αβ (q, 0) = 0.The above relations clearly show that Eqs. ( 19) of the original paper cannot be correct if the correlator K appearing there is the gauge kernel.While δG/δA µ = δΓ/δA µ descends from the properties of the Legendre transform, a similar relation does not apply to the second derivative of the G-and Γ-functionals with respect to the gauge field.The correct form of Eq. ( 10) is with A proof of this relation is given in Appendix A. This proves that the K-correlator in Eq. ( 10) of the original paper is not the gauge kernel.Thus, in Eqs.(15)(16)(17)(18)(19)(20) of the original paper one should replace K µν (q, ω) with (minus) the right hand side of Eq. ( 6).Eqs.(20) of the original paper, despite being correct if one properly replaces K µν (q, ω), are not particularly useful as it is not clear how to compute the second derivative of the Γ-functional with respect to the gauge field within a diagrammatic approach.In the following, we present a revised and more useful derivation of the Ward identities making use of the sole G-functional.
We consider Eq. ( 7) of the original paper using the parametrization in Eqs.(8).We write J 1 (x) = h + δJ 1 (x), where h is a uniform symmetry breaking field.Taking the derivative of Eq. ( 7) of the original paper once with respect to A µ and once with respect to J 2 and setting A µ = J 2 = δJ 1 = 0, we obtain where , and L h µ,a (x, x ′ ) are the gauge kernel, susceptibility, condensate fraction, and L-correlator, respectively (see Eq. ( 7)), computed in presence of a uniform explicit symmetry breaking field h.Combining the two equations above, we obtain where ∂ ′ ν is a derivative over x ′ .Considering the global U(1) symmetry of the problem, one can obtain a relation similar to Eq. ( 7) of the original paper: where x is an integral over the spatio-temporal coordinates.From the above relation, one readily obtains This, combined with the properties χ h 22 (−q, ω) = χ h 22 (q, −ω) = χ h 22 (q, ω), implies that the Fourier transform of the right hand side of Eq. ( 9) is at least quadratic in frequency and/or momentum.We can thus assume the following small-q and small-ω expansion of χ h 22 (q, ω) Taking the Fourier transform of Eq. ( 9), one gets Finally, taking the h → 0 limit, we obtain where χ n and J αβ are the coefficients of the expansion (12) obtained in absence of a symmetry breaking field.In particular, J αβ is the superfluid stiffness.Also, φ 0 and χ 22 (q, ω) are the condensate fraction and susceptibility obtained for zero h.
The equations above are a more useful version of the Ward identities presented in the original paper.The calculation of lim h→0 K h µν can be performed in the same way as one computes the gauge kernel, introducing a symmetry breaking field, and sending it to zero at the end of the calculation.The equations above also clearly prove that the dynamical (Eq.(14a)) or static (Eq.(14b)) limits do not commute with the h → 0 limit in the gauge kernel.Inverting their order, one would obtain Eqs.(5) instead of those above.

II. SU(2) SYMMETRY
Most of the considerations above also apply to the case of a spontaneously broken SU(2) symmetry, with minor modifications.The SU(2) gauge kernel has a similar definition to its U(1) counterpart.Also in this case, the SU(2) gauge kernel can only be obtained by taking the derivative of the G-functional (and not from the Γ functional).
Relations similar to Eqs. ( 6), (7), can be obtained for a system with SU(2) symmetry, making, also in this case, Eq. ( 40) of the original paper not useful for practical calculations.In the following, we derive WIs for the actual SU(2) gauge kernel.

III. EXPLICIT CALCULATION FOR A SPIRAL MAGNET
Sec. III of the original paper is fully correct, as the microscopic expressions for the spin stiffnesses and dynamical susceptibilities correspond to those derived in presence of a small symmetry breaking field that is sent to zero at the end of the calculation.
A misprint is present in Eq. (A4) of the original paper.Its correct form is that is, compared to the original paper, the prefactor is ∆ and not ∆ 2 .

IV. NOTE
In Ref. [65] the expressions for the spin stiffnesses and dynamical susceptibilities derived in the original manuscript have been used.Even though it was not explicitly stated, they have been derived applying a small symmetry breaking field to the system and sending it to zero after performing the dynamical or static limit.The expressions in Ref. [65] are therefore correct within the (renormalization group improved) random phase approximation employed in the paper.
During the review process of this Erratum, Ref. [66] appeared, which presents a similar derivation of the Ward identities above, and, as also discussed here, identifies the presence of an infinitesimal symmetry breaking field as crucial to obtain the correct formulas for the spin stiffnesses.

Figure 1 .
Figure 1.Diagrams contributing to the spin stiffnesses.The wavy line represents the external SU(2) gauge field, the solid lines the electronic Green's functions, the black triangles the paramagnetic vertex γ µ k σ a , the black circle the diamagnetic one γ αβ k σ 0 , and the dashed line the effective interaction Γ(q, q ′ , ω).