SU(2) gauge theory of the pseudogap phase in the two-dimensional Hubbard model

We present a SU(2) gauge theory of fluctuating magnetic order in the two-dimensional Hubbard model. The theory is based on a fractionalization of electrons in fermionic chargons and bosonic spinons. The chargons undergo N\'eel or spiral magnetic order below a density dependent transition temperature $T^*$. Fluctuations of the spin orientation are described by a non-linear sigma model obtained from a gradient expansion of the spinon action. The spin stiffnesses are computed from a renormalization group improved random phase approximation. Our approximations are applicable for a weak or moderate Hubbard interaction. The spinon fluctuations prevent magnetic long-range order of the electrons at any finite temperature. The phase with magnetic chargon order exhibits many features characterizing the pseudogap regime in high-$T_c$ cuprates: a strong reduction of charge carrier density, a spin gap, Fermi arcs, and electronic nematicity.


I. INTRODUCTION
Besides their exceptionally high transition temperatures for superconductivity, a peculiar and fairly universal feature of hole-doped cuprate superconductors is their pseudogap behavior for temperatures above T c , observed in a broad doping range from the underdoped into the optimally doped regime [1,2]. The pseudogap behavior sets in at a temperature T * which is much higher than T c in the underdoped regime, and merges with T c near optimal doping. It is characterized by a spin gap, a reduction of the charge carrier concentration, a suppression of the electronic density of states, a gap for single-particle excitations in the antinodal region of the Brillouin zone, and a reconstructed Fermi surface which in photoemission looks like Fermi arcs. It is also associated with a tendency to electronic nematicy, where the electronic state breaks the tetragonal symmetry of the crystal. Suppressing superconductivity by high magnetic fields, the pseudogap regime extends into the region which in the absence of magnetic fields is superconducting -for doping concentrations as high as about 20 percent [2].
There is convincing numerical evidence for pseudogap behavior in the strongly interacting two-dimensional Hubbard model, in particular from quantum cluster calculations [3]. An unbiased "fluctuation diagnostics" [4] of the contributions to the self-energy has revealed that the pseudogap is generated predominantly by antiferromagnetic fluctuations, that is, by spin fluctuations with wave vectors at or near (π, π).
While the guidance provided by numerical results is clearly very valuable, a deeper understanding of the pseudogap phenomenon and data with a higher momentum resolution remain desirable. The momentum resolution of the self-energy and of other momentum dependent quantities is limited in all cluster methods, since the computational effort grows exponentially with the cluster size. Also long-range correlations (in real space) beyond the cluster size cannot be captured. The same limitations hold of course for direct numerical simulations of finite systems.
In this situation approximate analytic theories can provide further insights, especially concerning longrange correlations and the fine-structure in momentum space. Early theories of the pseudogap phenomenon were based on weak-coupling diagrammatic perturbation expansions, most notably Moriya's renormalized theory [5] and the two-particle self-consistent theory by Vilk and Tremblay [6]. The Mermin-Wagner theorem on the absence of spin symmetry breaking at finite temperatures is respected in these theories, but the pseudogap seems to develop only for fairly large magnetic correlation lengths, while the numerical data show that strong short-ranged correlations are sufficient.
More recently it was shown by Sachdev, Scheurer, and coworkers that many features of the pseudogap behavior observed in cuprates can be captured by a SU(2) gauge theory [7][8][9][10][11]. This approach is based on a fractionalization of the electron into a fermionic "chargon" and a charge neutral "spinon". The latter is a SU(2) matrix providing a space and time dependent local reference frame [12]. The local spin rotations can be parametrized by a SU(2) gauge field, and the fractionalization leads to a gauge redundancy. One can then consider states where the chargons exhibit some sort of magnetic order (for example, Néel or spiral), while the spinon fluctuations prevent symmetry breaking and magnetic longrange order of the physical spin-carrying electrons [13][14][15]. Quantities involving only charge degrees of freedom behave essentially as in a conventional magnetically ordered state, and the Fermi surface gets correspondingly reconstructed. While long-range order is absent, at least at finite temperatures, the electrons are subject to a "topological" order in the sense that smoothly varying local spin rotations can map the fluctuating spin configurations to an ordered pattern, that is, there is no proliferation of topological defects [11].
In this paper we formulate a SU(2) gauge theory for the fluctuating antiferromagnet in a way that allows us to compute, in a decent approximation for weak or moderate Hubbard interactions, effective low-energy parameters and physical quantities as a function of the microscopic model parameters. The chargon order parameter arXiv:2207.00829v1 [cond-mat.str-el] 2 Jul 2022 is computed from a renormalized mean-field theory [16] which takes high energy (above T * ) spin, charge, and pairing fluctuations into account on equal footing. We allow for Néel and planar spiral order with generally incommensurate ordering wave vectors. The spinon dynamics is described by a non-linear sigma model (NLσM). The parameters of the NLσM, that is, the spin stiffnesses, are computed from a renormalized RPA for the SU(2) gauge field response of the chargons [17], and the ultraviolet cutoff is estimated via the magnetic coherence length. The NLσM is evaluated in a large N expansion. Applying the general theory to the Hubbard model with next and next-nearest neighbor hopping at a moderate interaction strength (about half band width), we obtain a broad finite temperature pseudogap regime on the hole doped side and a narrower pseudogap region for electron doping. Nematic order is present at sufficiently low temperatures for hole doping, but not for electron doping. There is no magnetic long-range order at T > 0, in agreement with the Mermin-Wagner theorem, and the spin excitations are gapped. The spinon quantum fluctuations are not strong enough to destroy magnetic long-range order in the ground state, except possibly near the edge of the pseudogap regime at large hole doping. In the hole doped pseudogap regime, the Fermi surfaces extracted from the single-particle spectral function have the form of hole pocket boundaries with a truncated back side. Their topology is thus the same as for the experimentally observed Fermi arcs.
The paper is structured as follows. In Sec. II we derive the general structure of the SU(2) gauge theory for the pseudogap phase with Néel or spiral order in the chargon sector. In Sec. III we describe how we compute the parameters of the gauge theory, in particular the spin stiffnesses, from the underlying microscopic model. Sec. IV deals with the solution of the nonlinear sigma model for the spinon fluctuations in a large N expansion. Results for the two-dimensional Hubbard model are presented in Sec. V. We conclude with a summary and a final discussion of our theory in Sec. VI.

A. Fractionalizing the electron field
We consider the Hubbard model on a square lattice with units of length such that the lattice spacing is one. The Hubbard action in imaginary time reads where c jσ = c jσ (τ ) and c * jσ = c * jσ (τ ) are Grassmann fields corresponding to the annihilation and creation, respectively, of an electron with spin orientation σ at site j, and n jσ = c * jσ c jσ . The chemical potential is denoted by µ, and U > 0 is the strength of the (repulsive) Hubbard interaction. To simplify the notation, we write the dependence of the fields on the imaginary time τ only if needed for clarity.
The action in (1) is invariant under global SU(2) rotations acting on the Grassmann fields as where c j and c * j are two-component spinors composed from c jσ and c * jσ , respectively, while U is a SU(2) matrix acting in spin space.
To separate collective spin fluctuations from the charge degrees of freedom, we fractionalize the electronic fields as [12][13][14][15] where R j ∈ SU(2), to which we refer as "spinon", is composed of bosonic fields, and the components ψ js of the "chargon" spinor ψ j are fermionic. According to (2) and (3) with V j ∈ SU(2). Hence, the components ψ js of ψ j carry a SU(2) gauge index s, while the components R j,σs of R j have two indices, the first one (σ) corresponding to the global SU(2) symmetry, and the second one (s) to SU(2) gauge transformations. We now rewrite the Hubbard action in terms of the spinon and chargon fields. The quadratic part of (1) can be expressed as [14] where we have introduced a SU(2) gauge field, defined as with ∂ µ = (i∂ τ , ∇). Here, the nabla operator ∇ is defined as generator of translations on the lattice, that is, e −r jj ·∇ with r jj = r j − r j is the translation operator from site j to site j .
To rewrite the interacting part in (1), we use the decomposition [12,14,18] where n j = n j,↑ + n j,↓ is the charge density operator, σ = (σ 1 , σ 2 , σ 3 ) are the Pauli matrices, andΩ j is an arbitrary time-and site-dependent unit vector. Inserting the decomposition (3), the interaction term of the Hubbard action can therefore be written as where n ψ j = ψ * j ψ j is the chargon density operator, S ψ j = 1 2 ψ * j σψ j is the chargon spin operator, andΩ R j is a unit vector obtained by rotatingΩ j as Using (7) again, we obtain with n ψ js = ψ * js ψ js . Therefore, the final form of the action S = S 0 + S int is nothing but the Hubbard model action where the physical electrons have been replaced by chargons coupled to a SU(2) gauge field.
Since the chargons do not carry the physical spin degree of freedom, a global breaking of their SU(2) gauge symmetry ( S ψ j = 0) does not necessarily imply long range order for the electrons. The matrices R j describe directional fluctuations of the order parameter S j , where, at low temperatures, the most important ones vary slowly in space and time.

B. Nonlinear sigma model
We now derive a low energy effective action for the spinon fields R j by integrating out the chargons, Since the action S is quartic in the fermionic fields, the functional integral must be carried out by means of an approximate method. In previous works [12][13][14] a Hubbard-Stratonovich transformation has been applied to decouple the chargon interaction, together with a saddle point approximation on the auxiliary bosonic (Higgs) field. We will employ an improved approximation based on the functional renormalization group [19], which we describe in Sec. III. The effective action for the spinons can be obtained by computing the response functions of the chargons to a fictitious SU(2) gauge field. Since we assign only low energy long wave length fluctuations to the spinons in the decomposition (3), the spinon field R j is slowly varying in space and time. Hence, we can perform a gradient expansion. To second order in the gradient ∂ µ R j , the effective action S eff [R] has the general form where x = (τ, r) combines imaginary time and space coordinates, T = [0, β] × R 2 is the integration region, and repeated indices are summed. We have expanded the gauge field A µ in terms of the SU(2) generators, with a running from 1 to 3. In line with the gradient expansion, the gauge field is now defined over a continuous space-time. The coefficients in (12) do not depend on x and are given by where • ( • c ) denotes the (connected) average with respect to the chargon Hubbard action. The first and second order current vertices have been defined as where x jj and y jj are the x and y components, respectively of r jj = r j − r j .
In Appendix A we will see that the linear term in (12) vanishes. We therefore consider only the quadratic contribution to the effective action. Defining the adjoint representation R of the SU(2) rotation R via we obtain the non-linear sigma model (NLσM) action for the directional fluctuations (see Appendix B) where P µν = 1 2 Tr[J µν ]1 − J µν . The structure of the matrices J µν and P µν depends on the magnetically ordered chargon state. In the trivial case S ψ j = 0 all the stiffnesses vanish and no meaningful low energy theory for R can be derived. A well-defined low-energy theory emerges, for example, when Néel antiferromagnetic order is realized in the chargon sector, that is, whereû is an arbitrary fixed unit vector. Choosingû = e 1 = (1, 0, 0), the spin stiffness matrix in the Néel state has the form with (J µν ) = diag(−Z, J, J). In this case the effective theory reduces to the well-known O(3)/O(2) S 2 nonlinear sigma model [20,21] whereΩ a = R a1 , and |Ω| 2 = 1. Another possibility is planar spiral magnetic ordering of the chargons, where Q is a fixed wave vector as obtained by minimizing the chargon free energy, whileû 1 andû 2 are two arbitrary mutually orthogonal unit vectors. The special case Q = (π, π) corresponds to the Néel state. Fixingû 1 toê 1 and u 2 toê 2 ≡ (0, 1, 0), the spin stiffness matrix assumes the form where for a ∈ {⊥, 2}. In this case, the effective action maintains its general form (18) and it describes the O(3)×O(2)/O(2) symmetric NLσM, which has been previously studied in the context of geometrically frustrated antiferromagnets [22][23][24][25][26]. This theory has three independent degrees of freedom, corresponding to one in-plane and two out-of-plane Goldstone modes. In the following we will restrict the magnetic ordering pattern of the chargons to Néel or planar spiral order. Néel or spiral antiferromagnetism has been found in the two-dimensional Hubbard model over broad regions of the parameter space by several approximate methods, such as Hartree-Fock [27], slave boson mean-field theory [28], expansion in the hole density [29], moderate coupling fRG [30], and dynamical mean-field theory [31,32]. In our theory the mean-field order applies only to the chargons, while the physical electrons are subject to order parameter fluctuations.

III. COMPUTATION OF PARAMETERS
In this section, we describe how we evaluate the chargon integral in Eq. (11) to compute the magnetic order parameter and the stiffness matrix J µν . The advantage of the way we formulated our theory in Sec. II is that it allows arbitrary approximations on the chargon action. One can employ various techniques to obtain the order parameter and the spin stiffnesses in the magnetically ordered phase. We use a renormalized mean-field (MF) approach with effective interactions obtained from a functional renormalization group (fRG) flow. In the following we briefly describe our approximation of the (exact) fRG flow, and we refer to Refs. [19,33,34] for the fRG, and to Refs. [16,30,35,36] for the fRG+MF method.

A. Symmetric regime
We evaluate the chargon functional integral by using an fRG flow equation [19,33,34], choosing the temperature T as flow parameter [37]. Temperature can be used as a flow parameter after rescaling the chargon fields as ψ j → T 3 4 ψ j , and defining a rescaled bare Green's function, , where ν n = (2n + 1)πT is the fermionic Matsubara frequency, and k is the Fourier transform of the hopping matrix in (1).
We approximate the exact fRG flow by a second order (one-loop) flow of the two-particle vertex V T , discarding self-energy feedback and contributions from the threeparticle vertex [19]. In a SU(2) invariant system the twoparticle vertex has the spin structure where k α = (k α , iν αn ) are combined momentum and frequency variables. Translation invariance imposes momentum conservation so that k 1 + k 2 = k 3 + k 4 . We perform a static approximation, that is, we neglect the frequency dependency of the vertex. To parametrize the momentum dependence, we use the channel decomposition [38][39][40][41] where the functions φ p,T , φ m,T , and φ c,T capture fluctuations in the pairing, magnetic, and charge channel, respectively. The dependences of these functions on the linear combination of momenta in the brackets are typically much stronger than those in the subscripts. Hence, we expand the latter dependencies in form factors [38,42], keeping only the lowest order s-wave, extended s-wave, p-wave and d-wave contributions. We run the fRG flow from the initial temperature T ini = ∞, at which V Tini = U , down to a critical temperature T * at V T diverges, signaling the onset of spontaneous symmetry breaking (SSB). If the divergence of the vertex is due to φ m,T , the chargons develop some kind of magnetic order.

B. Order parameter
In the magnetic phase, that is, for T < T * , we assume an order parameter of the form ψ * k,↑ ψ k+Q,↓ , which corresponds to Néel antiferromagnetism if Q = (π, π), and to spiral order otherwise.
For T < T * we simplify the flow equations by decoupling the three channels φ p,T , φ m,T , and φ c,T . The flow equations can then be formally integrated, and the formation of an order parameter can be easily taken into account [16]. We focus on magnetic order and ignore the pairing instability to analyze the non-superconducting "normal" state. In the magnetic channel one thus obtains the magnetic gap equation [30] where f (x) = (e x/T + 1) −1 is the Fermi function, k is a shorthand notation for d 2 k (2π) 2 , and E ± k are the quasiparticle dispersions The effective couplingV m k,k (Q) is the particle-hole irreducible part of V T * in the magnetic channel, which can be obtained by inverting a Bethe-Salpeter equation at the critical scale, where V m,T k,k (q) = V T (k−q/2, k +q/2, k −q/2, k+q/2), and the particle-hole bubble is given by Although V m,T * k,k (q) diverges at certain wave vectors q = Q c , the irreducible couplingV m k,k (q) is finite for all q. The dependence ofV m k,k (q) on k and k is rather weak and of no qualitative importance. Hence, to simplify the calculations, we discard the k and k dependencies of the effective coupling by taking the momentum averagē V m (q) = k,k V m k,k (q). The magnetic gap then becomes momentum independent, that is, ∆ k = ∆. While the full vertex V m,T k,k (q) depends very strongly on q, the dependence of its irreducible partV k,k (q) on q is rather weak. The calculation of the stiffnesses in the subsequent section is considerably simplified approximatinḡ V m (q) by a momentum independent effective interaction U m eff =V m (Q c ). The optimal ordering wave vector Q is found by minimizing the mean-field free energy of the system where the chemical potential µ is determined by keeping the density n = k =± f (E k ) fixed. The optimal wave vectors Q at temperatures T < T * generally differ from the wave vectors Q c at which V T * k,k (q) diverges. Eq. (26) has the form of a mean-field gap equation with a renormalized interaction that is reduced compared to the bare Hubbard interaction U by fluctuations in the pairing and charge channels. This reduces the critical doping beyond which magnetic order disappears, compared to the unrealistically large values obtained already for weak bare interactions in pure Hartree-Fock theory (see e.g. Ref. [27]).

C. Spin stiffnesses
The NLσM parameters, that is, the spin stiffnesses J ab µν , are obtained by evaluating Eq. (15). These expressions can be viewed as the reponse of the chargon system to an external SU(2) gauge field in the low energy and long wavelength limit, and they are equivalent to the stiffnesses defined by an expansion of the inverse susceptibilities to quadratic order in momentum and frequency around the Goldstone poles [17,43]. The following evaluation is obtained as a simple generalization of the RPA formula derived in Ref. [17] to a renormalized RPA with effective interactions U m eff and U c eff . Since in the magnetic state the translational symmetry is broken, the Fourier transforms of the response functions depend on two distinct momenta q and q , where q can assume the values q, q ± Q, and q ± 2Q. However, to compute J ab µν , we only need to deal with the limit q, q → 0.
The temporal components of the stiffness matrix, that is, J ab 00 , are given by the uniform spin susceptibility in the dynamical limit [17] J ab where χ ab (q, q , ω) is the Fourier transform of Note that, in a metallic system, the static uniform susceptibility obtained from q, q → 0 after setting ω = 0 differs from the quantity defined in Eq. (31).
The spin susceptibility can be most conveniently computed in a rotating spin frame defined by the transformation [43,44] since in the rotated frame the magnetically ordered system appears translation invariant. Hence, the rotated susceptibility is diagonal in momentum space and can therefore be written as χ ab (q, ω), with a single momentum variable q.
Consistently with the mean-field theory for the magnetic order parameter, we compute the susceptibilities in the magnetic state via a random phase approximation (RPA) with renormalized interactions as obtained from the fRG. In a spiral state with a generic wave vector Q, the spin susceptibility is coupled to the charge susceptibility [44]. Hence, we extend the definition of the spin susceptibility in Eq. (32) to a combined charge-spin susceptibility by including the value 0 for the indices a and b, in addition to the values 1, 2, 3, and defining σ 0 as the two-dimensional unit matrix. The prefactor 1 4 in Eq. (32) implies that χ 00 is actually a quarter of the conventional charge susceptibility.
In RPA, the rotated susceptibility χ can be written as where q = (q, ω), and Γ ab (q) is the RPA effective interaction in the rotated spin frame. The "bare" susceptibility χ ab 0 (q) is given by the particle-hole bubble [43] with G(k, iν n ) the mean-field chargon Green's function in the rotated basis The RPA effective interaction is obtained from a ladder sum, leading to the linear matrix equation where with V c,T k,k (q) = 2V T (k−q/2, k +q/2, k+q/2, k −q/2) − V T (k−q/2, k +q/2, k −q/2, k+q/2).
Here we keep the dependence on q since it does not complicate the calculations. The off-diagonal (a = b) elements of χ ab (0, 0, ω) with a, b = 1, 2, 3 vanish for ω → 0 both in the spiral and in the Néel state, so that we need to deal only with the diagonal spin susceptibility components χ aa (q, q, ω).
In the limits ω → 0 and q → 0 or Q, several offdiagonal matrix elements of the RPA effective interaction Γ ab (q, ω) vanish [43]. The expressions for Z a can therefore be simplified to [17] and where superscripts + and − attached to χ 0 indicate that the susceptibilities are formed with the ladder operators S ± = 1 2 S 1 ± S 2 .
In the Néel state, terms which contribute to χ 11 (q, q , ω) and χ 22 (q, q , ω) with q = q ± 2Q for Q = (π, π), contribute to the momentum diagonal susceptibilities since 2Q ≡ 0 for Q = (π, π). The transformation of the susceptibilities from the rotated to the unrotated basis then reads [43] Since χ 11 (Q, ω) = 0 and χ 22 (q + Q) = χ 33 (q) for Q = (π, π), we obtain with Z = χ 22 (0, 0, ω → 0) in the Néel state, which can be evaluated by the same expression as the one for Z 2 in Eq. (42). The spatial components of the stiffness matrix J ab αβ with α, β = 1, 2 are obtained from the spatial components of the uniform gauge field kernel K ab αβ in the static limit [17] where K ab αβ (q, q , ω) = K p,ab αβ (q, q , ω) + δ ab K d αβ is the Fourier transform of K ab αβ,jl (τ ) = K p,ab αβ,jl (τ ) + K d αβ δ ab δ jl δ(τ ), with the paramagnetic and diamagnetic contributions (cf. Eq. (15)) respectively. The diamagnetic contribution is translation invariant (in a spiral or Néel state). Fourier transforming, and evaluating the expectation value in Eq. (48) by the renormalized RPA, the paramagnetic part of the spin stiffness can be written as where Γ ab (q, q , ω) is the RPA effective interaction (37) in the original (non-rotated) spin basis. The bare paramagnetic response kernel is given by where γ (1) (k) = (1, ∂ kx k , ∂ ky k ) is the Fourier transform of γ (1) (j, j ) in Eq. (16a). The chargon Green's function in the original spin basis reads The diamagnetic part of the spin stiffness can be obtained from the Green's function as where γ (2) αβ (k) = ∂ kα ∂ k β k is the Fourier transform of the second order current vertex γ (2) αβ (j, j ) in Eq. (16b). The various contributions to the response kernel are represented diagrammatically in Fig. 1.
In the Néel state one finds, in close analogy to the temporal components J ab 00 of the stiffness matrix, J 11 αβ = 0 and J 22 αβ = J 33 αβ = Jδ αβ , which can be most easily computed from the right hand side of Eq. (55).
In our low energy theory of the spinons we have ignored possible imaginary contributions from Landau damping of the Goldstone modes. In a Néel state, they are of the same order in the gradient expansion as the (real) temporal and spatial stiffness terms [45]. The same is true for the Landau damping of the in-plane mode in a spiral state, but the damping of the out-of-plane mode is of higher order [43]. Moreover, it requires the existence of hot spots (connected by Q) of the reconstructed Fermi surface. In our large N evaluation of the NLσM, the in-plane modes of the spiral state do not contribute. Hence, for the spiral state, Landau damping is irrelevant for our theory, while in the Néel state their might be a quantitative (not qualitative) modification of our results.
We conclude this section by comparing our theory to the SU(2) gauge theory of the half-filled Hubbard model derived by Borejsza and Dupuis [14]. They used the same fractionalization of the electron in chargons and spinons, and the chargon order was treated in (plain) mean-field theory. Our expressions for the spin stiffnesses agree with theirs (at half-filling) if we replace our renormalized interaction U m eff by the bare Hubbard interaction U , although their derivation differs from ours. Following earlier work by Haldane [20,21] for the Heisenberg model, Borejsza and Dupuis obtained their expressions for the spin stiffnesses by splitting the spinon into a "Néel field" and a "canting field" describing ferromagnetic fluctuations. Integrating out the fermions and the canting field they obtained a NLσM for the Néel field, where the stiffnesses are given by the RPA. We obtain the same expressions (with a renormalized coupling) more directly from the RPA evaluation of the gauge field response, without introducing the canting field.

IV. EVALUATION OF SIGMA MODEL
To solve the NLσM, we resort to a saddle point approximation in the CP N −1 representation, which is exact in the large N limit [46,47].

A. CP 1 representation
The matrix R can be expressed as a triad of orthonormal unit vectors: whereΩ i ·Ω j = δ ij . We represent these vectors in terms of two complex Schwinger bosons z ↑ and z ↓ [45] with z = (z ↑ , z ↓ ) andΩ ± =Ω 1 ∓ iΩ 2 . The Schwinger bosons obey the nonlinear constraint The parametrization (59) is equivalent to Inserting the expressions (58) and (59) into Eq. (18) and assuming a stiffness matrix J µν of the form (23), we obtain the CP 1 action for fluctuating spiral order with sum convention for repeated greek indices and the current operator For the Néel case, the CP 1 action is given by the same expression with J 2 µν = 0. We recall that x = (τ, r) comprises the imaginary time and space variables, and T = [0, β] × R 2 .

B. Large N expansion
The current-current interaction in Eq. (62) can be decoupled by a Hubbard-Stratonovich transformation, introducing a U(1) gauge field A µ , and implementing the constraint (60) by means of a Lagrange multiplier λ. The resulting form of the action describes the so-called massive CP 1 model [48] where D µ = ∂ µ − iA µ is the covariant derivative. The numbers M µν are the matrix elements of the mass tensor of the U(1) gauge field, where J 2 and J ⊥ are the stiffness tensors built from the matrix elements J 2 µν and J ⊥ µν , respectively. To perform a large N expansion, we extend the twocomponent field z = (z ↑ , z ↓ ) to an N -component field z = (z 1 , . . . , z N ), and rescale it by a factor N/2 so that it now satisfies the constraint To obtain a nontrivial limit N → ∞, we rescale the stiffnesses J ⊥ µν and J 2 µν by a factor 2/N , yielding the action This action describes the massive CP N −1 model [49], which in d > 2 dimensions displays two distinct critical points [47,48,50]. The first one belongs to the pure CP N −1 class, where M µν → 0 (J 2 µν = 0), which applies, for example, in the case of Néel ordering of the chargons, and the U(1) gauge invariance is preserved. The second is in the O(2N) class, where M µν → ∞ (J ⊥ µν = J 2 µν ) and the gauge field does not propagate. At the leading order in N −1 , the saddle point equations are the same for both fixed points, so that we can ignore this distinction in the following.
At finite temperatures T > 0 the non-linear sigma model does not allow for any long-range magnetic order, in agreement with the Mermin-Wagner theorem. The spin correlations decay exponentially and the spin excitations are bounded from below by a spin gap m s = i λ /Z ⊥ . In the large N limit, the spin gap m s is related to the spin stiffness by the following equation (see Appendix C for a derivation) where Λ uv is an ultraviolet momentum cutoff. The constant J is an "average" spin stiffness given by and c s = J/Z ⊥ is the corresponding average spin wave velocity. In Sec. IV C, we shall discuss how to choose the value of Λ uv . For m s c s Λ uv , and T c s Λ uv , the magnetic correlation length ξ s = 1 2 c s /m s , behaves as with the critical stiffness The correlation length is finite at each T > 0. For J > J c , ξ s diverges exponentially for T → 0, while for J < J c it remains finite in the zero temperature limit. At T = 0, the bosons may condense and the saddle point condition yields where n 0 = | z 1 | 2 is the fraction of condensed bosons. Eq. (72) can be easily solved, yielding (if m s Λ uv ) The Mermin-Wagner theorem is thus respected already in the saddle-point approximation to the CP N −1 representation of the nonlinear sigma model, that is, there is no long-range order at T > 0. In the ground state, longrange order (corresponding to a z boson condensation) is obtained for a sufficiently large spin stiffness, while for J < J c magnetic order is destroyed by quantum fluctuations even at T = 0, giving rise to a paramagnetic state with a spin gap.

C. Choice of ultraviolet cutoff
The impact of spin fluctuations described by the nonlinear sigma model depends strongly on the ultraviolet cutoff Λ uv . In particular, the critical stiffness J c separating a ground state with magnetic long-range order from a disordered ground state is directly proportional to Λ uv . The need for a regularization of the theory by an ultraviolet cutoff is a consequence of the gradient expansion. While the expansion coefficients (the stiffnesses) are determined by the microscopic model, there is no systematic way of computing Λ uv .
A pragmatic choice for the cutoff is given by the ansatz where C is a dimensionless number, and ξ A is the magnetic coherence length, which is the characteristic length scale of spin amplitude correlations. This choice may be motivated by the observation that local moments with a well defined spin amplitude are not defined at length scales below ξ A [14]. The constant C can be fixed by matching results from the nonlinear sigma model to results from a microscopic calculation in a suitable special case (see below). The coherence length ξ A can be obtained from the connected spin amplitude correlation function χ A (r j , r j ) = (n j · S ψ j )(n j · S ψ j ) c , wheren j = S ψ j /| S ψ j |. At long distances between r j and r j this function decays exponentially with an exponential dependence e −r/ξ A of the distance r. Fourier transforming and using the rotated spin frame introduced in Sec. III C, the long distance behavior of χ A (r j , r j ) can be related to the momentum dependence of the static correlation function χ ab (q, 0) in the amplitude channel a = b = 1 for small q, which has the general form The magnetic coherence length is then given by where . The constant C in Eq. (74) can be estimated by considering the Hubbard model with pure nearest neighbor hopping (with amplitude −t) at half-filling. At strong coupling (large U ) the spin degrees of freedom are then described by the antiferromagnetic Heisenberg model, which exhibits a Néel ordered ground state with a magnetization reduced by a factor n 0 ≈ 0.6 compared to the mean-field value [51]. On the other hand, evaluating the RPA expressions for the Hubbard model in the strong coupling limit, one recovers the mean-field results for the spin stiffness and spin wave velocity of the Heisenberg model with an exchange coupling J H = 4t 2 /U , namely J = J H /4 and c s = √ 2J H . Evaluating the RPA spin amplitude correlation function yields ξ A = 1/ √ 8 in this limit. With the ansatz (74), one then obtains n 0 = 1 − 4C/π. Matching this with the numerical result n 0 ≈ 0.6 yields C ≈ 0.3 and Λ uv ≈ 0.9.
We finally note that we are not overcounting any fluctuations in our theory. In general, the electron fractionalization in Eq. (3) introduces redundant degrees of freedom associated with the gauge symmetry, Eq. (4). We have not explicitly fixed a gauge but, due to our (renormalized) mean-field treatment of the chargons, fluctuations of the magnetic order parameter are captured exclusively by the spinons.

V. RESULTS
In this section we present and discuss results obtained from our theory for the two-dimensional Hubbard model, both in the hole-(n < 1) and electron-doped (n > 1) regime. We allow for next and second nearest neighbor hopping with amplitudes −t and −t , respectively, and we fix the ratio of the hopping amplitudes as t /t = −0.2, and we choose a moderate interaction strength U = 4t. The energy unit is t in all plots.

A. Chargon mean-field phase diagram
The critical temperatures T * m and T * p at which the vertex V T (k 1 , k 2 , k 3 , k 4 ) diverges are shown in Fig. 2. In a density range from n ≈ 0.83 to n ≈ 1.08, the divergence of the vertex is due to a magnetic instability. Beyond the edges of this density interval, the leading instability occurs in the d-wave pairing channel. Pairing extends into the magnetic regime at lower temperatures (below T * m ) as a secondary instability. Vice versa, magnetic order is possible at temperatures below T * p in the regime where pairing fluctuations dominate [16,30]. In Fig. 2, we also show the irreducible effective magnetic interaction U m eff defined in Sec. III B. The effective interaction U m eff is strongly reduced from its bare value (U = 4t) by the non-magnetic channels in the fRG flow, while its density dependence is not very strong.
From now on we ignore the pairing instability and focus on the magnetic order of the chargons. We compute the magnetic order parameter ∆ together with the optimal wave vector Q as described in Sec. III B. In Fig. 3, we show results for ∆ in the ground state (T = 0) as a function of the filling. We find a stable magnetic solution extending deep into the hole doped regime down to n ≈ 0.73. On the electron doped side magnetic order terminates abruptly already at n ≈ 1.08. This pronounced electron-hole asymmetry and the discontinuous transition on the electron doped side has already been observed in previous fRG+MF calculations for a slightly weaker interaction U = 3t [30].
The onset temperature T * for magnetic order of the chargons as obtained from the renormalized mean-field theory is shown in Fig. 2. At densities where magnetic interactions dominate, it coincides with the temperature T * m at which V T diverges. At densities where the interaction diverges in the pairing channel, T * is lying only slightly below T * p on the hole doped side, while it vanishes on the electron doped side. While the magnetic gap in the ground state reaches its peak at n = 1, as expected, the pseudocritical temperature T * and the irreducible effective interaction U m eff exhibit their maximum in the hole doped regime slightly away from half-filling.
The magnetic states are either Néel type or spiral with a wave vector of the form Q = (π − 2πη, π), or symmetry related, with an "incommensurability" η > 0. In Fig. 3 results for η in the ground state are shown as a function of the density. At half-filling and in the electron doped region only Néel order is found, as expected and in agreement with previous fRG+MF studies [30]. Hole doping instead immediately leads to a spiral ground state with η > 0. Whether the Néel state persists at small hole doping depends on the hopping parameters and the interaction strength. Its instability toward a spiral state is favored by a larger interaction strength [29]. Indeed, in a previous fRG+MF calculation at weaker coupling the Néel state was found to survive up to about 10 percent hole doping [30].
At low and moderate hole doping, there is a transition between a Néel state at high temperatures and a spiral state at low temperatures. Since the spiral state breaks the tetragonal symmetry of the square lattice, spiral order entails electronic nematicity. In Fig. 2 we show the corresponding nematic transition temperature T nem as a function of density. T nem merges with T * at n ≈ 0.88. For lower densities the magnetic order is spiral with η > 0 at any temperature below the magnetic transition temperature. Within the spiral regime there is a topological transition of the quasiparticle Fermi surface (indicated by the black dashed line in Fig. 2), where hole pockets merge. The Fermi surface extracted from the singleparticle spectral function develops Fermi arcs on the right hand side of this transition, while it resembles the large bare Fermi surface on the left (see Sec. V C).

B. Spinon fluctuations
Once the magnetic order parameter ∆ of the chargons and the wave vector Q have been computed, we are in the position to calculate the NLσM parameters from the expressions presented in Sec. III C.
In Fig. 4, we plot results for the spatial and temporal spin stiffnesses J a αα and Z a in the ground state. In the spiral state (for n < 1) out-of-plane and in-plane stiffnesses are distinct, while in the Néel state (for n ≥ 1) they coincide. Actually the order parameter defines an axis, not a plane, in the latter case. All the quantities except Z 2 exhibit pronounced jumps between half-filling and infinitesimal hole-doping. These discontinuities are due to the appearance of hole pockets around the points ( π 2 , π 2 ) in the Brillouin zone [43]. The spatial stiffnesses are almost constant over a broad range of hole-doping, with a small spatial anisotropy J a xx = J a yy . The temporal stiffnesses Z a exhibit a stronger doping dependence. The peak of Z ⊥ at n ≈ 0.79 is associated with a van Hove singularity of the quasiparticle dispersion [43]. On the electron doped side all stiffnesses decrease almost linearly with the electron filling. The off-diagonal spin stiffnesses J a xy and J a yx vanish both in the Néel state and in the spiral state with Q = (π − 2πη, π) and symmetry related.
In Fig. 5 we show the magnetic coherence length ξ A and the average spin wave velocity c s in the ground state. The coherence length is rather short and only weakly doping dependent from half-filling up to 15 percent holedoping, while it increases strongly toward the spiral-toparamagnet transition on the hole-doped side. On the electron-doped side it almost doubles from half-filling to infinitesimal electron doping. This jump is due to the for-  mation of electron pockets upon electron doping. Note that ξ A does not diverge at the transition to the paramagnetic state on the electron doped side, as this transition is first order. The average spin wave velocity exhibits a pronounced jump at half-filling, which is inherited from the jumps of J ⊥ αα and Z ⊥ . Besides this discontinuity it does not vary much as a function of density.
We now investigate whether the magnetic order in the ground state is destroyed by quantum fluctuations or not. To this end we compute the boson condensation fraction n 0 as obtained from the large-N expansion of the NLσM. This quantity depends on the ultraviolet cutoff Λ uv . As a reference point, we may use the half-filled Hubbard model at strong coupling, as discussed in Sec. IV C, which yields Λ uv ≈ 0.9, and the constant in the ansatz Eq. (74) is thereby fixed to C ≈ 0.3.
In Fig. 6 we show the condensate fraction n 0 com- puted with two distinct choices of the ultraviolet cutoff: Λ uv = Λ uv (n) = C/ξ A (n) and Λ uv = C/ξ A (n = 1). For the former choice the cutoff vanishes at the edge of the magnetic region on the hole-doped side, where ξ A diverges. One can see that n 0 remains finite for both choices of the cutoff in nearly the entire density range where the chargons order. Only near the hole-doped edge of the magnetic regime, n 0 vanishes slightly above the mean-field transition point, if the ultraviolet cutoff is chosen as density independent. The discontinuous drop of n 0 upon infinitesimal hole doping is due to the corresponding drop of the out-of-plane stiffness, while the discontinuous increase of n 0 upon infinitesimal electron doping, for the density dependent cutoff choice Λ uv (n) = C/ξ A (n), is due to the discontinuity of ξ A (n). In the weakly hole-doped region there is a substantial reduction of n 0 below one, for both choices of the cutoff. Except for the edge of the magnetic region on the hole-doped side, the choice of the cutoff has only a mild influence on the results, and the condensate fraction remains well above zero. Hence, we can conclude that the ground state of the Hubbard model with a moderate coupling U = 4t is magnetically ordered over wide density range. The spin stiffness is sufficiently large to protect the magnetic order against quantum fluctuations of the order parameter.

C. Electron spectral function
Fractionalizing the electron operators as in Eq. (3), the electron Green's function assumes the form To simplify this expression, we decouple the average RR * ψψ * as RR * ψψ * , yielding [9,10,14] [ The spinon Green's function can be computed from the NLσM in the continuum limit. Using the Schwinger boson parametrization (61), we obtain, in the large N limit, The boson propagator D(r, τ ) is the Fourier transform of with the bosonic Matsubara frequency ω n = 2πnT . Fourier transforming Eq. (78), the electron Green's function is obtained in momentum representation as where G(k, k , ν n ) is the chargon Green's function. We see that when n 0 = 0, the electron Green's function is diagonal in momentum, that is, it is translational invariant, as the diagonal components of the chargon Green's function entering the trace are nonzero only for k = k . Furthermore, in this case G e is proportional to the unity matrix in spin space, since there is no spin SU(2) symmetry breaking, and is thus given by a single normal state Green's function G e (k, ν n ). Performing the Matsubara sum in Eq. (81) and continuing to real frequencies, we get where and n B (x) = (e x/T − 1) −1 is the Bose distribution function. In the right column of Fig. 7 we show the spectral function obtained as the imaginary part of the retarded electron Green's function at zero frequency as a function of momentum for various electron densities in the hole doped regime. The temperature T = 0.05t is below the chargon ordering temperature in all cases. The Fermi surface topology is the same as the one obtained from a mean field approximation of spiral spin density wave order [52]. At low hole doping it originates from a superposition of hole pockets (see left column of Fig. 7), where the spectral weight on the back sides is drastically suppressed by coherence factors, so that only the front sides are visible. The spinon fluctuations lead to a broadening of the spectral function, so that the Fermi surface is smeared out. Since the spinon propagator does not depend on the fermionic momentum, the broadening occurs uniformly in the entire Brillouin zone. Hence, the backbending at the edges of the "arcs" obtained in our theory for n = 0.9 is more pronounced than experimentally observed in cuprates. The backbending edges could be further suppressed by including a momentum dependent self-energy which has a larger imaginary part in the antinodal region [53].

VI. CONCLUSIONS
We have presented a SU(2) gauge theory of fluctuating magnetic order in the two-dimensional Hubbard model. The theory is based on a fractionalization of the electron operators in chargons and spinons [12][13][14][15]. The chargons are treated in a renormalized mean-field theory with effective interactions obtained from a functional renormalization group flow. They undergo Néel or spiral magnetic order in a broad density range around half-filling below a density dependent temperature T * . Fluctuations of the spin orientation are described by a non-linear sigma model obtained from a gradient expansion of the spinon degrees of freedom. The parameters of the sigma model, the spin stiffnesses, have been computed from a renormalized RPA. Our approximations are applicable for a weak or moderate Hubbard interaction U . While magnetic long-range order of the electrons is still possible in the ground state, at any finite temperature the spinon fluctuations prevent long-range order -in agreement with the Mermin-Wagner theorem. We expect that at strong coupling even the ground state becomes disordered already at relatively low hole-doping, since fluctuations are then enhanced due to the shorter magnetic coherence length.
In spite of the moderate interaction strength chosen in our explicit calculations, the phase with magnetic chargon order below T * exhibits all important features characterizing the pseudogap regime in high-T c cuprates. The Fermi surface reconstruction yields a reduction of the electronic density of states. At low hole doping the Fermi surface obtained from the spectral function for singleparticle excitations looks like Fermi arcs. The spinon fluctuations generate a spin gap at any finite temperature. The spinon fluctuations do not contribute to quantities involving only charge degrees of freedom, such as the longitudinal or Hall conductivities. It was already shown previously that Néel or spiral order of the chargons can explain the drastic charge carrier drop observed at the onset of the pseudogap regime in hole-doped cuprates [32,52,[54][55][56][57][58].
In the Néel regime, the structure of our theory is very similar to the SU(2) gauge theory of the pseudogap phase derived by Sachdev and coworkers [7][8][9][10][11]15]. Besides our extension to spiral states, the major new aspect of our work is that we compute the magnetic order parameter and spin stiffnesses instead of fitting the parameters of the theory. This computation revealed in particular an important particle-hole asymmetry of the stiffnesses. Spiral order of the chargons entails nematic order of the electrons. At low hole doping, the chargons form a Néel state at T * , and a spiral state at a lower temperature T nem . The electrons thus undergo a nematic phase transition at a critical temperature below the pseudogap temperature. Evidence for a nematic transition at a temperature T nem < T * has been found recently in slightly underdoped YBCO [59]. For large hole doping instead, the nematic transition occurs right at T * , while nematic order is completely absent for electron doping, that is, above half-filling.
In the ground state of the two-dimensional Hubbard model there is a whole zoo of possible magnetic ordering patterns, and away from half filling Néel or spiral order do not always minimize the ground state energy. The most important competitor is stripe order, that is, collinear spin order associated with charge order, where holes accumulate in one-dimensional lines [3]. Stripe order in the ground state has been established rather convincingly for special cases, such as pure nearest neighbor hopping and doping concentration 1/8 [60]. The energy difference between distinct order patterns can be very small. At finite temperatures, the issue of the proper choice of the magnetic order reappears for the chargons. A classification of the numerous possibilities has been provided recently by Sachdev et al. [61]. We have focused on Néel and spiral states because any other state leads to a fractionalization of the Fermi surface into numerous tiny pieces (infinitly many for incommensurate wave vectors), which is in conflict with the experimental observation of only four arcs in the pseudogap phase of cuprates. Moreover, it is hard to explain the sharp carrier drop observed at the edge of pseudogap regime in high magnetic fields via collinear magnetic order [62]. Hence, to us Néel or spiral order of the chargons seems the most promising starting point to understand the universal features of the pseudogap phase. Refinements are required to capture also secondary instabilities, that is, charge order and superconductivity.
At finite temperatures, we obtain a "pseudogap" phase with a reconstructed Fermi surface and a spin gap also for the electron doped Hubbard model. In contrast, in electron doped cuprates one observes a comparatively broad (in doping) Néel phase, and no or only a very narrow pseudogap regime. Néel order at finite temperature is possible due to the interlayer coupling in cuprates. On the hole doped side, interlayer coupling stabilizes the Néel state only in a very narrow regime near half-filling. This electron-hole asymmetry can be explained by the asymmetry of the spin stiffnesses, which are much smaller on the hole doped side (see Fig. 4), enhancing thus the impact of spin fluctuations.