Spiral to stripe transition in the two-dimensional Hubbard model

We obtain an almost complete understanding of the mean-field phase diagram of the two-dimensional Hubbard model on a square lattice with a sizable next-nearest neighbor hopping and a moderate interaction strength. In particular, we clarify the nature of the transition region between the spiral and the stripe phase. Complementing previous [Phys. Rev. B 108, 035139 (2023)] real-space Hartree-Fock calculations on large finite lattices, we solve the mean-field equations for coplanar unidirectional magnetic order directly in the thermodynamic limit, and we determine the nature of the magnetic states right below the mean-field critical temperature $T^*$ by a Landau free energy analysis. While the magnetic order for filling factors $n \geq 1$ is always of N\'eel type, for $n \leq 1$ the following sequence of magnetic states is found as a function of increasing hole-doping: N\'eel, planar circular spiral, multi-spiral, and collinear spin-charge stripe states. Multi-spiral states are superpositions of several spirals with distinct wave vectors, and lead to concomitant charge order. We finally point out that nematic and charge orders inherited from the magnetic order can survive even in the presence of fluctuations, and we present a corresponding qualitative phase diagram.


I. INTRODUCTION
The two-dimensional Hubbard model on a square lattice plays a key role in the field of strongly correlated electron systems as a prototype model for competing and intertwined ordering tendencies.It captures the most salient features of electrons in the copper oxide planes of high-T c cuprates, namely antiferromagnetism and d-wave superconductivity [1].Thanks to remarkable advances in the development of computational techniques, fragments of the phase diagram of this important model have been established [2,3], but many regions in the large parameter space spanned by hopping amplitudes, interaction strength, electron filling, and temperature remain terra incognita.
A large variety of magnetic phases in the twodimensional Hubbard model emerges already in a conventional static mean-field approximation, that is, Hartree-Fock theory.While the regime of ordered states in the phase diagram is usually overestimated by mean-field theory, qualitative insights may serve as a guide for more sophisticated calculations, in particular, to interprete data from numerical simulations on finite lattices.Numerous Hartree-Fock studies of the two-dimensional Hubbard model have already been published.In many of them the magnetic order was restriced to certain pat-terns, such as ferromagnetic and Néel order [32][33][34] or, more generally, to spiral order with arbitrary wave vectors [11,12].Allowing for collinear magnetic order with arbitrary wave vectors or even for completely arbitrary spin configurations, spin-charge stripes have been discovered [17][18][19][20][21][22].
Mean-field theory yields magnetic order also at finite temperatures, below a transition temperature T * , violating thus the Mermin-Wagner theorem [35].However, magnetically ordered states at finite temperature become meaningful in theories of fluctuating magnetic order, where the electron is fractionalized into a fermion with a magnetically ordered pseudospin, and a fluctuating SU(2) rotation matrix which restores the SU(2) spin symmetry [36][37][38].
Recently, we have performed a comprehensive and unbiased mean-field analysis of magnetic and charge order in the Hubbard model with a moderate interaction strength on a square lattice, at both zero and finite temperatures [39].Completely unrestricted real-space Hartree-Fock calculations on large finite lattices were combined with a stability analysis of mean-field solutions restricted to Néel and spiral order in the thermodynamic limit.It turned out that in most parts of the phase diagram only three classes of magnetic states with a relatively simple structure are stabilized in the thermodynamic limit: Néel, circular spiral, and collinear stripe states.The stripes are usually unidirectional, but can also be bidirectional at very large hole doping in presence of a sizable next-nearest neighbor hopping.In spite of rather large lattices (up to 48 × 48) used in the realspace calculations, the analysis of the stripe states was still hampered by finite size effects, and the nature of the transition from the spiral to the stripe phase remained open.
In this paper we complete the mean-field analysis of our previous work [39] by performing several complementary calculations.First, we solve the mean-field equations for ground states with generic coplanar unidirectional order directly in the thermodynamic limit.This includes the Néel, spiral, and unidirectional stripe states found in the real-space calculations [39] as special cases.Second, we determine the magnetic ordering pattern right below T * from a Landau free energy analysis.Third, we clarify the nature of the instability of the spiral state upon increasing doping by analyzing its spin susceptibility, again directly in the thermodynamic limit.We find that the transition from the spiral to the stripe phase leads through a rather complex intermediate phase with a superposition of multiple spiral components with three or four distinct wave vectors.Finally, we present a qualitative discussion of fluctuation effects.Order parameter fluctuations restore the SU(2) spin symmetry at least at finite temperature, but nematic and charge orders found in the mean-field states may survive.
The remainder of our paper is structured as follows.In Sec.II we describe our three complementary methods used to compute the mean-field phase diagram and to clarify the nature of the various magnetic states.In Sec.III we present the corresponding results.In the Conclusion in Sec.IV we summarize and present a qualitative discussion of fluctuation effects.

II. MODEL AND METHOD
The Hubbard Hamiltonian for spin- 1  2 fermions with intersite hopping amplitudes t jj ′ and a local repulsive interaction U > 0 reads [2,3] where c jσ (c † jσ ) annihilates (creates) an electron on lattice site j with spin orientation σ ∈ {↑, ↓}, and n jσ = c † jσ c jσ .The hopping matrix t jj ′ depends only on the distance between the sites j and j ′ .We choose t jj ′ = −t if j and j ′ are nearest neighbor sites, t jj ′ = −t ′ if j and j ′ are next-to-nearest neighbors, and t jj ′ = 0 otherwise.We use the nearest neighbor hopping amplitude t as our energy unit.
In mean-field theory, the interaction term in (1) can be decoupled as [18,39] where c j is the two-component spinor composed of c j↑ and c j↓ , while σ 0 is the two-dimensional identity matrix and σ 1 , σ 2 , σ 3 are the Pauli matrices.The sign s a is plus one if a = 0, and minus one otherwise.The parameters ∆ a j are related to charge and spin expectation values as The mean-field decoupling in Eq. ( 2) captures both the Hartree (a = 0, 3) and the Fock (a = 1, 2) terms.
Previous unbiased and unrestricted real-space meanfield calculations on the Hubbard model [39] revealed that, except for very low electron densities, the solutions of the mean-field equations always converge to coplanar unidirectional phases.Thus, in this paper we focus on mean-field states characterized by one or more wave vectors of the form Q = (π − δ, π) or symmetry related (we call this property unidirectionality), and where all the spins lie in a common plane (coplanarity).This includes collinear spin states as special cases (with infinitly many common planes), and in particular the Néel state as the collinear state with Q = (π, π).
We analyze the different phases that one can obtain within mean-field theory in the Hubbard model by employing three distinct but mutually consistent techniques.
A) To find the magnetic ground state, we employ a Hartree-Fock ansatz that allows for a generic coplanar unidirectional state with an arbitrary integer periodicity P in x-direction (and antiferromagnetic order in ydirection).The contributing wave vectors have the form Q = (2πn/P, π), with integer numbers n ≥ 0. The meanfield equations are solved directly in the thermodynamic limit.The periodicity P determines the size of the realspace unit cell one has to deal with in the solution of the mean-field equations.We could reach unit cells as big as 220 sites along the x-axis, so that even incommensurate states without any translation symmetry (that is, P = ∞) are approximated very well.B) Close to the critical temperature T * at which meanfield magnetic order sets in, it is technically hard to obtain converged solutions of the mean-field equations.Therefore, to determine the pattern of the magnetic order setting in right below T = T * , we employ a Landau theory for mixed spin-charge order parameters, and microscopically compute its coefficients from the paramagnetic state at T = T * .Note that in the limit of a vanishing order parameter, Landau theory and Hartree-Fock theory yield the same type of order.
C) It has been previously observed [4][5][6][7][8][9][10][11][12][13][14][15][16]39] that a circular spiral magnetic state is favored for small hole doping (that is, slightly below half-filling) in presence of a finite t ′ < 0. We study the instabilities of the spiral state to other magnetic orders at larger hole dopings by computing the spin and charge susceptibilities in such a state within random phase approximation (RPA).This enables us to determine not only when the spiral state becomes unstable, but also the nature of the magnetic order emerging beyond the instability line.Note that the RPA is the unique conserving approximation for susceptibilities which is consistent with mean-field theory for the free energy, order parameters, and single-particle properties [40].
In the following, we provide a detailed description of the three methods mentioned above.
A. Mean-field theory for a generic coplanar unidirectional magnetic state We derive mean-field equations in momentum space which describe generic magnetic states characterized by the following three properties.
Coplanarity.The onsite magnetization should lie in a specific plane, which we choose, without loss of generality, to be the xy-plane.
Unidirectionality.Spins sitting on neighboring sites along the, say, y-direction are antiparallel and the charge densities are equal, while along the x-direction the magnetization amplitude and orientation, as well as the charge density, can be arbitrarily modulated.
Commensurability.Spin and charge orders display a periodicity with respect to translations along the x-axis with a finite integer period (denoted by P ) in units of the lattice spacing.This criterion implies a restriction to states with ordering wave vectors commensurate with the lattice.Incommensurate states can be approximated to any desired accuracy by choosing a sufficiently large value for P .
In a coplanar, unidirectional, commensurate state with periodicity P , the magnetization and charge density profiles can be expressed as follows where R j are the coordinates of lattice site j, and Q P = (2π/P, π).The first sum is running only over odd integers n because the spin order is antiferromagnetic in y-direction, while the second sum is restricted to even integers since ρ j is translation invariant in y-direction.
We define P S as the smallest positive integer satisfying P S Q P = (0, 0) modulo reciprocal lattice vectors.If P is even, one has P S = P , and P S = 2P if P is odd.Because nQ P is equivalent to (n + mP S )Q P with m, n ∈ Z, the summations in Eq. ( 4) run only over a finite number (P S /2) of terms.Moreover, since the spin and charge densities on the left hand side of Eqs.(4) are real, the coefficients M x n , M y n , and ϱ n must obey Spin and charge orderings break the translational symmetry of the original lattice, resulting in an enlarged unit cell containing P S inequivalent sites with distinct expectation values.Similarly, in momentum space, the size of the original Brillouin zone is reduced by a factor 1/P S .The new reciprocal lattice can be constructed by adding all vectors of the form nQ P with n = 0, ..., P S − 1 to the original reciprocal lattice vectors (see Fig. 1).The reduced Brillouin zone can then be defined as the set of all points that are closer to a given vector of the new reciprocal lattice than to any other (Wigner-Seitz construction).The reduced Brillouin zones for P = 6, 7, 8 are plotted in Fig. 1.
Inserting Eq. ( 4) into the Hamiltonian (1) with the mean-field decoupling (2), and Fourier transforming, we obtain the quadratic Hamiltonian where Fourier transform of the hopping parameters in Eq. ( 1).
We have also defined k is a shorthand for the integral k∈BZ Introducing a "Nambu spinor" with P S components with the convention ↑ = ↓ and ↓ = ↑, one can cast the Hamiltonian (6) in the form where we have dropped the constant term in Eq. ( 6) and defined k,σ has the form where we have defined p(ℓ) = + if ℓ is even and p(ℓ) = − if ℓ is odd, and Since ∆ 0 0 only shifts the chemical potential µ, in the following we redefine µ as µ − ∆ 0 0 and set ∆ 0 0 = 0.In Appendix A we report the explicit form of the matrix for the cases P = 3 and 6.
The parameters ∆ a n are self-consistently determined as where the matrices Γ a,n σ have been defined as k,↓ are related to each other by an inversion of the sign of ∆ 2 n , the expectation values on the right hand side of Eq. ( 12) take the same value for each of the two spin projections.For this reason, one can simplify Eq. ( 12) to where Ψ k = Ψ k,↑ and Γ a n = Γ a n,↑ .In other words, we can solve the mean-field equations using only the matrix H (P ) k,↑ .The right hand side of Eq. ( 14) is computed from Eq. ( 10) making an initial random assumption on the mean-field parameters ∆ a n , which are then updated using again Eq. ( 14).The procedure is repeated until convergence is reached.
To find the energetically best state, we converge Eq. ( 14) for different values of P and retain the state with the lowest mean-field free-energy.In practice, we discretize the original Brillouin zone BZ with N 2 k equally spaced points and, for a fixed N k , we only allow values of P that are divisors of N k .For every fixed set of parameters, we have offered the system over 90 integer values of P ranging from 2 to 220, each of them with a suitably adjusted N k such that N k /P ∈ N, with N k ranging from 116 to 220.
In the following, we discuss how several important familiar phases are captured as special cases within our general formalism.

Néel order
In the case of Néel order, one has where φ parametrizes the orientation of the spin order in the xy-plane, and M represents its amplitude.Néel order has the period P = 2, so that the matrix k,σ in Eq. ( 10) is two-dimensional and only the two parameters ∆ 1  1 and ∆ 2 1 contribute, where with ∆ = U M .

Circular spiral order
Circular spiral order has the form where, as in the case of Néel order, φ parametrizes the orientation of the spin order in the xy plane and M its (constant) amplitude.Q is a generic wave vector of the form (π − δ, π) with δ > 0. The "+" or "-" sign distinguishes between spirals rotating anti-clockwise and clockwise.This type of order emerges as a special case of our general formalism if ∆ a n is non-zero only for one mode n and its conjugate P S − n, with n odd and such that Q can be approximated by nQ P modulo a reciprocal lattice vector, with a suitably chosen P .Spiral order as in Eq. ( 17) is then described by with ∆ = U M ∈ R. All other ∆ a n are zero.The matrix H (P ) k,σ thus simplifies to a diagonal block matrix form with P S /2 matrices of size two on the diagonal.Indeed spiral order can be described by a simpler 2 × 2 meanfield Hamiltonian for each k-point, as previously used in mean-field calculations restriced to spiral states [6-9, 11, 13-16, 39].Note that the Eqs.(18) apply only if n ̸ = P S − n, which is fulfilled for any spiral state which is not a Néel state, that is, as long as Q ̸ = (π, π).

Stripe order
We define as stripe order any type of collinear order that differs from Néel antiferromagnetism.In this case, the magnetization and charge densities have the form where once again φ is an angle parameterizing the orientation of the spins in the xy-plane, while f S (R j ) and f ρ (R j ) are two functions defining the spatial modulation of the magnetization amplitude and charge density, respectively.They can be expressed in terms of their Fourier coefficients as M n e inQ P •Rj , (20a) Stripe order can be obtained as a particular case of our general formalism, with ∆ a n fulfilling Depending on the coefficients M n and ρ n , the spin and charge profiles can be sinusoidal or sharp (like domain walls) or anything in between.

B. Landau theory close to T *
To determine the type of magnetic order close to the critical temperature T * , we decouple the Hubbard interaction by introducing spin and charge order parameter fields via a Hubbard-Stratonovich transformation, and subsequently expand the resulting effective action in powers of the order parameters.

Derivation of the effective action
We write the Hubbard interaction as [41][42][43] where Ωj is an arbitrary site-and time-dependent unit vector.Intuitively, one can imagine Ωj as being the direction of the local (both in space and time) magnetization.Because Ωj can be arbitrarily chosen, in the path integral of the Hubbard model we take the average over all possible Ωj with a properly defined measure such that D Ω = 1.We cast the Hubbard interaction in the form (22), because this makes it compatible with our mean-field decoupling (see Eq. ( 2)).We perform a Hubbard-Stratonovich transformation to decouple each of the terms in Eq. ( 22) by means of two fields, ρ j and ρ S j , representing fluctuations of the charge and spin amplitude, respectively.Defining a spin field as ⃗ S j = ρ S j Ωj , we can represent the Hubbard interaction as where β = 1/T is the inverse temperature, and To keep the notation light, we have dropped the time dependence of the bosonic (ρ j , ⃗ S j ) and fermionic (c j , cj ) fields.
An effective action for the bosonic fields ρ j and ⃗ S j can be derived by integrating out the fermions.Except for a field-independent term, one obtains where G 0 is the Fourier transform to real space and imaginary time of the bare Matsubara Green's function , while ρ and ⃗ S are diagonal matrices in space and time defined as ).The trace is summing over space and time indices, G 0 A is the space-time matrix product j ′′ β 0 dτ ′′ G 0,jj ′′ (τ − τ ′′ ) A j ′′ j ′ (τ ′′ , τ ′ ), and 1 st = δ jj ′ δ(τ −τ ′ ) is the space-time unit matrix.

Taylor expansion of the effective action
We now expand the logarithm in Eq. ( 25) in powers of ρ j and ⃗ S j .Such an expansion is justified in the vicinity of the critical temperature T * , where the magnetic and charge order parameters are small.To this end we write tr ln [ with A = iρ+⃗ σ• ⃗ S. The trace does not depend on the representation.In the following we perform all calculations in momentum and frequency space.
The quadratic term in ⃗ S q , with ⃗ S q the (spatiotemporal) Fourier transform of ⃗ S j , takes the form k = T ν k is a shorthand for a sum over Matsubara frequencies and a momentum integration, and q = (q, Ω) is a collective variable comprising a lattice momentum and a bosonic Matsubara frequency.We define T * as the critical temperature where min q U −1 − Π 0 (q, 0) = 0, signaling an instability towards the formation of magnetic order.This condition is first met, in the most general case, at four symmetry related wave vectors in the Brillouin zone of the form q = ±Q x or q = ±Q y , with Q x = (π − δ, π) and Q y = (π, π − δ).This means that the magnetic order forming right below T * can be entirely characterized by these wave vectors.Note that, because of the opposite sign between the two terms on the right hand side of Eq. ( 22), the coefficient of the term quadratic in ρ q is U −1 + Π 0 (q), which is always positive.Thus, within mean-field theory, an instability towards charge order alone can never occur in the Hubbard model.
For a mean-field study of the Taylor-expanded effective action (25), we can therefore assume that ⃗ S q possesses solely modes at q = ±Q x and q = ±Q y close to T = T * , corresponding to the ansatz where ⃗ M x and ⃗ M y are constant complex vectors.We assume static fields, consistent with our mean-field treatment.
Third order terms involving only spin fields vanish due to time-reversal symmetry.The third order term involving two ⃗ S q fields and one ρ q field takes the form with a coupling function λ(q, q ′ ).Inserting Eq. ( 28) into this equation, we see that, within mean-field theory, the only charge modes that couple to ⃗ M x and ⃗ M y are those where q = 0, q = ±2Q x,y , and q = ±(Q x ± Q y ).Neglecting the q = 0 mode, which does not lead to any symmetry breaking, and higher order spin-charge interactions (this approximation will be justified below), we can write where ϕ x , ϕ y , and ϕ ± are complex constants, and Inserting Eq. ( 28) and (30) into Eq.( 25), and expanding up to quartic order in ⃗ The coefficients s, u 0 , u 1 , u 2 , u 3 , r 1 , r 2 , b 1 , and b 2 are determined by frequency and momentum integrals of products of bare propagators G 0 .The concrete expressions are presented in Appendix B.
The charge degrees of freedom can be eliminated from the theory by imposing ∂V /∂ϕ α = 0, for α = x, y, ±, which yields form (33) has previously been derived from general symmetry arguments [37,44,45].The form of the Landau theory restricted to the case a single mode ( ⃗ M x or ⃗ M y ) was derived earlier in Ref. [46].

C. Susceptiblities in the spiral state
With a proper redefinition of the local spin reference frame [47,48], spiral order as in Eq. ( 17) can be described in terms of a 2 × 2 Hamiltonian of the form Since the energy does not depend on the phase φ, we can choose, without loss of generality, φ = 0. Within the rotated reference frame, one can compute the charge and spin susceptibilities within random phase approximation (RPA) as where Γ 0 = 2U diag(−1, 1, 1, 1), and the bare susceptibility χ 0 (q, ω), as a function of the bosonic Matsubara frequency Ω, is given by Green's function.The real frequency susceptibility is obtained by substituting iΩ → ω +i0 + after performing the Matsubara sum.
To compute the susceptibilities in the physical (unrotated) spin reference frame, one must rotate Eq. ( 35) in the xy plane with a spatially dependent angle of Q • R j [47,48].Such a rotation will produce in general momentum off-diagonal components of the susceptibilities, as spiral order breaks translational invariance.However, for our purpose of a stability analysis of the spiral state it suffices to consider the susceptibilities in the rotated spin reference frame.

III. RESULTS
We now present the mean-field phase diagram of the two-dimensional Hubbard model as obtained from the three complementary methods described in the preceding section.We choose a sizable next-nearest neighbor hopping t ′ = −0.3t,as is frequently used to model the band structure and Fermi surface of the cuprate superconductor yttrium barium copper oxide (YBCO) [49].In the hole-doped region (n < 1) we obtain the same sequence of magnetic states also for other negative values of t ′ /t.For the interaction strength we choose U = 3t, which is strong enough to obtain magnetic order in spite of the magnetic frustration imposed by t ′ , but weak enough to obtain qualitatively plausible results from the Hartree-Fock approximation.In the cuprates the Hubbard interaction is much larger, but the effective interaction driving magnetic order and magnetic correlations is renormalized to smaller values by fluctuations.
Since the magnetic order at densities n ≥ 1 is generally of Néel type [39], we focus on the hole-doped regime n < 1, where an intriguing sequence of ordering patterns is found.

A. Ground state phase diagram
In Fig. 2 we show the ground state phase diagram as obtained from the mean-field solution described in Sec.II A. The Néel state at half-filling (n = 1) is immediately unstable toward a spiral state upon hole doping.At n ≃ 0.9, the spiral state becomes unstable, leading into a more complex phase at lower densities which is still coplanar and non-collinear, but with a modulated spin amplitude and charge density.We discuss this phase, which we call multi-spiral, in more detail in Sec.III C. Upon further increasing the hole-doping, a conventional stripe state with collinear spin order and charge density wave order is stabilized.The transitions from the Néel to the spiral state and from the spiral to the multi-spiral state are continuous, while the transition from the multispiral state to the stripe state might be first order.For a quantitative characterization of the various states, we define the average spin amplitude ⟨S⟩ and the average charge modulation ⟨δρ⟩ as where the lattice sum extends over one unit cell (with P S sites).The dominant wave vector is parametrized by the incommensurability with Q max = n max Q P , where n max is the index belonging to the largest magnetic gaps ∆ 1 n or ∆ 2 n .In Fig. 2 we see that, irrespective of the phase transitions occurring, both ⟨S⟩ and δ display a rather smooth and monotonic behavior.By contrast, ⟨δρ⟩ vanishes in the Néel and spiral phases, but then rises quickly in the multi-spiral regime, peaking at the transition to the stripe phase.It then slowly decays as the density is further decreased.

B. Phase diagram at T = T *
In Fig. 3 we show the various magnetic phases we obtain at T = T * by minimizing the effective potential (33).In a density regime near half-filling, at T = T * the maximum of Π 0 (q, 0) in Eq. (27b) occurs at q = (π, π), implying that the phase being realized right below the (mean- field) critical temperature is a Néel antiferromagnet.Reducing the density, Π 0 (q, 0) develops four identical maxima at q = ±Q x and ±Q y at T * .An analysis of the quartic terms in the effective potential (33) reveals that the Néel state is replaced by a spiral phase characterized by or by the same expression with ⃗ M x ↔ ⃗ M y .Here, ê1 and ê2 are two orthogonal real unit vectors.The spiral phase maintains a uniform charge density, but it breaks the C 4 rotational symmetry of the square lattice.
At larger hole dopings, the spiral phase is replaced by a stripe phase, such that with an arbitrary unit vector ê, or by the same expression with ⃗ M x and ⃗ M y interchanged.This phase displays collinear magnetic order, a modulation of the charge density ρ j − n ∝ cos(2Q x + 2φ), and it breaks the C 4 symmetry.
A smooth interpolation between (circular) spiral and stripe order is given by elliptical spiral order [46], with The parameter α allows for a smooth interpolation between a spiral (α = π/4) and a stripe (α = 0) phase.At the transition point between spiral and stripe order, the effective potential ( 33) is degenerate with respect to variations of α.For T < T * , this degeneracy is lifted by higher order terms (beyond quartic).
At even lower densities, we find a coplanar bidirectional stripe phase (CpBS), characterized by with orthogonal unit vectors ê1 and ê2 , and arbitrary phases φ 1 and φ 2 .The charge density is then modulated as It is therefore conceivable that, at a finite distance below the T = T * line, a new phase emerges between unidirectional stripe and CpBS orders, interpolating between the two.We have marked this possible intermediate phase with a question mark in Fig. 3.The white color in the low temperature regime of the CpBS phase indicates that we have not clarified the nature of this phase far below T * .More complex ordering patterns are possible there [39], but in this regime of very large hole doping any magnetic order is probably an artifact of mean-field theory, and thus of limited interest.

C. Instability of the spiral state
The spiral state is stable in a finite hole-doping range near half-filling.At larger hole-doping, collinear stripe states have the lowest energy.We now clarify the nature of the instability of the spiral state upon increasing hole-doping, and the transition to a stripe state.The instability of the spiral state can be detected by analyzing the static charge and spin susceptibilities.At zero frequency, the bare susceptibilities χab 0 in Eq. ( 36) vanish if a = 3 and b ̸ = 3 or a ̸ = 3 and b = 3.Hence, the a, b = 0, 1, 2 sector of the susceptibilities, corresponding to charge, spin amplitude, and in-plane spin orientation fluctuations, decouples from the a, b = 3 sector, which is associated with out-of-plane spin orientation fluctuations.An instability is signaled by a divergence and subsequent sign change of the susceptibilities in Eq. (35).Such a divergence must however be distinguished from divergences due to Goldstone modes.Within our conventions, the Goldstone modes of the spiral state manifest themselves as χ 22 (0, 0) = ∞ and χ 33 (±Q, 0) = ∞ in the rotated spin frame [47,48].The static out-of-plane spin susceptibility χ 33 (q, 0) remains always positive and finite for q ̸ = ±Q.
We therefore search for a diverging susceptibility in the a, b = 0, 1, 2 sector at q ̸ = 0, which is necessarily associated with an eigenvalue of the RPA denominator in Eq. ( 35) crossing zero.Hence, to determine the instability of the spiral state, and the nature of the magnetic order beyond the instability line, we study the eigenval-ues of the matrix with χ ab 0 (q) = χ ab 0 (q, 0).The spiral state is stable if the matrix has two positive and one negative eigenvalues for all q ̸ = 0, and viceversa when it is unstable.
We define Q ′ as the non-zero wave vector at which the absolute value of the second largest eigenvalue of D(q) has a global minimum.With some lengthy but straightforward algebra, one can prove that Q ′ is the momentum at which | χ 12 0 (q)| is maximal.In the ground state, Q ′ is entirely determined by the geometry of the Fermi surface in the spiral state, which, at least for small dopings, consists of two hole pockets centered at ( π+δ 2 , ± π 2 ) with δ > 0, see panel (b) of Fig. 4. In panel (a) of Fig. 4 we see that | χ 12 0 (q)| has pronounced peaks at crossing points of two ellipse-shaped lines in q space, on which | χ 12 0 (q)| exhibits a singularity.These lines are "2k F -lines" [50] corresponding to the set of wave vectors connecting points with parallel tangents on the Fermi surfaces of the hole pockets.They can be geometrically constructed by shifting the two hole pockets such that their centers coincide with the Γ point (0, 0), and rescaling them by a factor of two.The global maximum of | χ 12 0 (q)| occurs where the two 2k F -lines cross.As displayed in Fig. 4, there are two pairs of crossings, one occurring on the q x axis (q y = 0), and one on the q y axis (q x = 0).Since χ 12 0 (q) is identically zero along the q y axis [48], Q ′ and −Q ′ are determined as the points in momentum space where the two 2k F -lines cross on the q x axis.Using these prescriptions, an analytical expression for Q ′ = (q ′ , 0) can be derived: where μ = µ/(2t) and ∆ = ∆/(2t).
The eigenvector of D(q) corresponding to the smallest positive eigenvalue (in the regime of stability of spirals) or the largest negative one (in the regime of instability of spirals) can be shown to take the general form (ϱ 0 , 1, iγ), with ϱ 0 , γ ∈ R.This form can be deduced by using that χ 12 0 (q) = − χ 21 0 (q) and χ 02 0 (q) = − χ 20 0 (q) are purely imaginary, while all other entries of D(q) are purely real [48].The form of the eigenvector corresponding to the eigenvalue of D(q) that can cross zero enables us to derive the form of the magnetic and charge ordering occurring right beyond the instability line.In the rotated frame in which spiral order appears as ferromagnetic, the order parameters take the form where M ′ is an overall amplitude and φ ′ a phase.Assuming an anti-clockwise rotating spiral proportional to cos(Q • R j )ê x + sin(Q • R j )ê y , corresponding to γ > 0, rotating Eqs.(45) to the physical spin reference frame yields where The charge order parameter is left unchanged by the rotation.Eq. ( 46) describes a magnetic state with three overlapping spirals with distinct wave vectors, two of which propagate anti-clockwise (those with Q and Q ′ − ), and one clockwise (with Q ′ + ).Thus, one can label this state as 3Q spiral.Such a state is found also by our numerical calculations using the formalism in Sec.II A. In Table I, we report the ground state values of Q, Q ′ ± , Q ′ , γ, and ρ 0 as predicted by the analysis of the susceptibilities in the spiral state, and as computed by solving the mean-field equations from II A for n = 0.90, that is, right beyond the instability line (see Fig. 2 for comparison).The lowest energy state was found to have a period of P = 36.In Fig. 5, we show the spin and charge patterns for a 3Q spiral state.
If the strength of the spiral order is weakened by raising the temperature, χ 12 0 (q) approaches [Π 0 (q + Q, 0) − Π 0 (q − Q, 0)]/(2i), with Π 0 (q) the bare bubble defined as in Eq. (27b).At the onset of magnetic order Π 0 (q, 0) is peaked exactly at Q (and symmetry related), which implies that for very weak ∆ one has Q ′ = −2Q modulo a reciprocal lattice vector.Similarly, in this limit one observes that χ 01 0 (q), and χ 02 0 (q) become zero, and that χ 11 0 (q) and χ 22 0 (q) approach the same value.For this reason the eigenvector of D(q) that can cross zero takes the form (0, 1, i) in the limit of vanishing spiral order, that is, γ → 1.This, together with Q ′ → −2Q, implies that Eq. ( 46) takes the form of an elliptical spiral (see Eq. ( 41)).Thus, the multi-spiral phase smoothly turns into an elliptical spiral phase as one raises the temperature toward T * , as schematically shown in Fig. 6.
When moving away from the instability line of the spiral phase by increasing the doping, we find a fourth mode emerging, such that the spin order assumes the form where we have dropped possible phases in the sine and cosine functions, and We also observe, upon increasing doping, that Hence, S j in Eq. ( 47) could gradually turn into a collinear stripe order with two harmonics: At the lowest density n = 0.87 evaluated numerically in the multi-spiral regime, we find M 2 ≈ 0.5M 1 , while at n = 0.86 we already find a stripe phase with a single Q-vector.Hence, either M 2 grows to approach M 1 very quickly in a small density range, or the transition from multi-spiral to stripe is of first order.This is the reason why we have interrupted the lines that serve as a guide to the eye in Fig. 2 at the transition point.
, and M ′ (see text) as predicted from calculations in the spiral phase (first row) and as numerically obtained from the mean-field theory of Sec.II A (second row).6.Schematic phase diagram of the mean-field transition from spiral to stripe order.We call in general multispiral the intermediate phase between spiral and stripe order.Closer to the spiral phase, we find three contributing Q-vectors (Eq.( 46)), while close to the stripe phase, we observe four (Eq.( 47)).Close to the critical temperature, the multi-spiral phase asymptotically approaches elliptical spiral order with only one Q-vector.

IV. CONCLUSION
Complementing our previous real-space Hartree-Fock study [39] by various additional techniques, we have obtained an almost complete understanding of the meanfield phase diagram of the two-dimensional Hubbard model with a moderate interaction strength.A large variety of distinct magnetic states appears, some with and some without concomitant charge order.Since, in presence of a sizable next-nearest neighbor hopping, the magnetic states in the electron doped regime (filling n > 1) are always Néel ordered [39], we focused on the hole doped regime n < 1.
The analysis in Ref. [39] showed that the magnetic order of the Hubbard model is always coplanar and usually unidirectional, with wave vectors of the form (π − δ, π).Bidirectional order was found only at very small densities (large hole doping).Allowing for arbitrary coplanar and unidirectional order, we were able to solve the mean-field equations directly in the thermodynamic limit.In the ground state, we thereby confirmed the circular spiral order at low hole-doping and the stripe order at large hole-doping.In between, we discovered a new multispiral phase consisting of a superposition of various spirals with distinct but unidirectional wave vectors.Unlike the single-component spiral phase, the multi-spiral phase exhibits charge order similar to the stripe phase.Analyzing the spin-charge susceptibility of the spiral phase at its instability point, we found that the additional wave vectors contributing to the multi-spiral phase are related to nesting vectors of the hole-pockets in the simple spiral state.We complemented the ground state calculation by a Landau free-energy analysis of the magnetic states right below the mean-field transition temperature T * , where we found the following sequence of states as a function of increasing hole doping: Néel -circular spiral -unidirectional stripe -bidirectional stripe.The multispiral phase found in the ground state becomes narrower (in doping) upon increasing temperature, and collapses to a point at T * .Approaching that point from below (T < T * ), the multi-spiral state converges to an elliptical spiral with a single wave vector.
Our results are thus largely consistent with the previous real space Hartree-Fock calculation on large but finite lattices [39].Only the multi-spiral phase could not be identified in the real space calculation, since the superposition of three or more contributing wave vectors leads naturally to very large unit cells.
To keep our paper concise, we fixed the next-nearest neighbor hopping and the Hubbard interaction to one value, t ′ = −0.3tand U = 3t, respectively, in all numerical results.The nearest-neighbor hopping t sets the global energy scale.Qualitative changes of the phase diagram upon changing these parameters can be described as follows.Setting t ′ = 0, the phase diagram becomes electron-hole symmetric.The magnetic order of the ground states is either Néel (at half-filling), or stripe (away from half-filling).Spiral states near half- filling exist only at finite temperatures in this special case [39].From a continuity argument it is clear that for very small finite t ′ , spiral and stripe order will still be present also in the electron doped regime.However, already for t ′ /t = −0.15, the electron-doped regime is exclusively Néel ordered for U = 3t [39].Bidirectional stripe order at large hole doping appears only for a rather large t ′ .It is absent for t ′ /t = −0.15[39].Decreasing U obviously reduces the magnetically ordered regime, both in density and temperature.For sufficiently weak U and a negative t ′ /t, a Néel ground state can be stable even for (small) finite hole doping [6,13].For t ′ = 0 there is magnetic order at and near half-filling for any non-zero U , due to perfect nesting of the half-filled Fermi surface, while for t ′ ̸ = 0 a certain minimal interaction strength is required.
We finally discuss how order parameter fluctuations affect the phase diagram.The Mermin-Wagner theorem [35] dictates that the spin SU(2) symmetry cannot be broken at any finite temperature.Hence, magnetic long-range order appearing in mean-field theory is destroyed by order parameter fluctuations.However, some secondary order parameters emerging in the phases described above may survive in the form of vestigial order.Moreover, features of the spectral function for fermionic single-particle excitations in the magnetically ordered regime, such as the Fermi arcs in the Néel and spiral regimes [14,51,52], may also survive [36,38,53].The mean-field critical temperature T * then becomes a crossover temperature, below which the electronic spectral function develops gaps for momenta in the antinodal region and Fermi arcs near the nodal region.The magnetic phases obtained in mean-field theory can therefore be interpreted as pseudogap (PG) phases at T > 0 (see Fig. 7).Whether the ground state remains magnetically ordered depends on the strength of the quantum fluctuations.
The breaking of the discrete (not continuous) C 4 rotational symmetry in the spiral, multi-spiral, and stripe phases can survive fluctuations, leaving a nematic phase, whose (mean-field) transition temperature T nem is indicated in Fig. 7.At low hole doping nematic order sets in at T nem < T * , while at larger doping T nem and T * coincide (with T * of course not being sharply defined in the presence of fluctuations).Charge density wave (CDW) order, displayed by the multi-spiral and stripe phases, can also survive the presence of fluctuations.Fig. 2 indicates that the incommensurability δ is a continuous function of the electron density, at least on the mean-field level.This implies that, except for certain special fillings, CDW order is generically incommensurate.Such an order has an emergent U(1) symmetry [54] and can exist at finite temperature only in the form of topological order in a Berezinskii-Kosterlitz-Thouless (BKT) phase.In Fig. 7, we sketched the mean-field CDW transition temperature T CDW .The BKT transition temperature for CDW quasi long-range order (QLRO), T BKT CDW , is expected to be lower than T CDW .The temperature range T BKT CDW ≤ T ≤ T CDW is a regime of fluctuating short range CDW order.It is possible that at lower temperature CDW fluctuations will lock the period of the charge modulation to a commensurate value, realizing a phase with long range order via an incommensurate to commensurate transition [54].Spin fluctuations will render the multi-spiral and stripe phase qualitatively identical at finite temperature, as they both break the same symmetries.

Figure 1 .
Figure 1.Construction of the reduced Brillouin zone for P = 6 (a), P = 7 (b), and P = 8 (c).The blue dots represent the momenta nQP with n = 0, ..., PS − 1 modulo vectors of the original reciprocal lattice.The black lines separate distinct reduced Brillouin zones.The original Brillouin zone is divided into PS equivalent reduced Brillouin zones.

Figure 3 .
Figure 3. Mean field phase diagram of the Hubbard model for t ′ = −0.3tand U = 3t.The various types of magnetic order appearing at T = T * were determined from Landau theory, the states at T = 0 by solving the Hartree-Fock equations derived in Sec.II A. The phase transitions for 0 < T < T * are indicated only schematically by straight lines connecting the calculated transition points at T = T * and T = 0.

Figure 4 .
Figure 4. Panel (a): Momentum dependence of | χ 12 0 (q)|, displaying nonanalyticities on the two so-called 2kF -lines.The absolute maxima occur where these two lines cross on the qx axis.In the upper half of the Brillouin zone (qy > 0), the 2kF -lines are retraced by dashed white lines as a guide to the eye.Panel (b): Fermi surface in the spiral state, consisting of two hole pockets.The scattering processes with a momentum transfer of Q ′ , connecting opposite sides of the pockets with parallel tangents, are highlighted by red arrows.Parameters: T = 0, t ′ = −0.3t,n = 0.90, U = 3t, for which mean-field theory yields ∆ ≈ 0.637t and δ ≈ 0.167π.

Figure 5 .
Figure 5. Spin and charge order pattern for a 3Q spiral state.Bigger (smaller) red bubbles represent a higher (lower) local hole concentration.The black arrows represent the local magnetization vector ⃗ Sj.The parameters are the same as in Tab.I, except for M ′ , which we have enhanced for visualization purposes.

Figure
Figure 6.Schematic phase diagram of the mean-field transition from spiral to stripe order.We call in general multispiral the intermediate phase between spiral and stripe order.Closer to the spiral phase, we find three contributing Q-vectors (Eq.(46)), while close to the stripe phase, we observe four (Eq.(47)).Close to the critical temperature, the multi-spiral phase asymptotically approaches elliptical spiral order with only one Q-vector.

Figure 7 .
Figure 7. Schematic phase diagram in presence of magnetic order parameter fluctuations in the most relevant density regime, excluding very large hole doping (cf.Fig. 3 for the corresponding mean-field phase diagram).Various pseudogap phases with and without nematic (C4 broken) and/or charge density wave (CDW) order are obtained.The hatched red region derives from the multi-spiral phase in mean-field theory, which has the same symmetries as the stripe regime in the presence of thermal fluctuations.

Table I .
Comparison of the values of Q