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

Pietro M. Bonetti  Max Planck Institute for Solid State Research,
Heisenbergstrasse 1, D-70569 Stuttgart, Germany  Walter Metzner  Max Planck
Institute for Solid State Research, Heisenbergstrasse 1, D-70569 Stuttgart,
Germany

(August 27, 2024)

###### Abstract

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éel
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.

###### pacs:

##  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 $(\pi,\pi)$.

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 long-range 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 long-range 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 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$\sigma$M). The parameters of the
NL$\sigma$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$\sigma$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.

##  II SU(2) gauge theory

###  II.1 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

$\displaystyle\mathcal{S}[c,c^{*}]$ $\displaystyle=$
$\displaystyle\int_{0}^{\beta}\\!d\tau\bigg{\\{}\sum_{j,j^{\prime},\sigma}c^{*}_{j\sigma}\left[\left(\partial_{\tau}-\mu\right)\delta_{jj^{\prime}}+t_{jj^{\prime}}\right]c_{j^{\prime}\sigma}$
(1) $\displaystyle+\;U\sum_{j}n_{j\uparrow}n_{j\downarrow}\bigg{\\}},$

where $c_{j\sigma}=c_{j\sigma}(\tau)$ and
$c^{*}_{j\sigma}=c^{*}_{j\sigma}(\tau)$ are Grassmann fields corresponding to
the annihilation and creation, respectively, of an electron with spin
orientation $\sigma$ at site $j$, and
$n_{j\sigma}=c^{*}_{j\sigma}c_{j\sigma}$. The chemical potential is denoted by
$\mu$, 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 $\tau$ only if needed for clarity.

The action in (1) is invariant under _global_ SU(2) rotations acting on the
Grassmann fields as

$c_{j}\to\mathcal{U}c_{j},\quad c^{*}_{j}\to
c^{*}_{j}\,\mathcal{U}^{\dagger},$ (2)

where $c_{j}$ and $c^{*}_{j}$ are two-component spinors composed from
$c_{j\sigma}$ and $c^{*}_{j\sigma}$, respectively, while $\mathcal{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]

$c_{j}=R_{j}\,\psi_{j},\quad c^{*}_{j}=\psi^{*}_{j}\,R^{\dagger}_{j},$ (3)

where $R_{j}\in\mbox{SU(2)}$, to which we refer as “spinon”, is composed of
bosonic fields, and the components $\psi_{js}$ of the “chargon” spinor
$\psi_{j}$ are fermionic. According to (2) and (3) the spinons transform under
the global SU(2) spin rotation by a _left_ matrix multiplication, while the
chargons are left invariant. Conversely, a U(1) charge transformation acts
only on $\psi_{j}$, leaving $R_{j}$ unaffected. We have therefore separated
the spin degrees of freedom of the physical electrons, now encoded in the
spinons, from their charge, carried by the chargons. The transformation in Eq.
(3) introduces a redundant SU(2) gauge symmetry, acting as

$\displaystyle\psi_{j}\to\mathcal{V}_{j}\,\psi_{j},\quad\psi^{*}_{j}\to\psi^{*}_{j}\,\mathcal{V}^{\dagger}_{j},$
(4a) $\displaystyle R_{j}\to R_{j}\,\mathcal{V}_{j}^{\dagger},\quad
R^{\dagger}_{j}\to\mathcal{V}_{j}\,R^{\dagger}_{j},$ (4b)

with $\mathcal{V}_{j}\in\mbox{SU(2)}$. Hence, the components $\psi_{js}$ of
$\psi_{j}$ carry a SU(2) gauge index $s$, while the components $R_{j,\sigma
s}$ of $R_{j}$ have two indices, the first one ($\sigma$) 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]

$\displaystyle\mathcal{S}_{0}[\psi,\psi^{*},R]$ $\displaystyle=$
$\displaystyle\int_{0}^{\beta}\\!d\tau\bigg{\\{}\sum_{j}\psi^{*}_{j}\left[\partial_{\tau}-\mu-
A_{0,j}\right]\psi_{j}$ (5)
$\displaystyle+\,\sum_{j,j^{\prime}}t_{jj^{\prime}}\,\psi^{*}_{j}\,e^{-\mathbf{r}_{jj^{\prime}}\cdot\left(\bm{\nabla}-i\mathbf{A}_{j}\right)}\,\psi_{j}\bigg{\\}},$

where we have introduced a SU(2) gauge field, defined as

$A_{\mu,j}=(A_{0,j},\mathbf{A}_{j})=iR^{\dagger}_{j}\partial_{\mu}R_{j},$ (6)

with $\partial_{\mu}=(i\partial_{\tau},\bm{\nabla})$. Here, the nabla operator
$\bm{\nabla}$ is defined as generator of translations on the lattice, that is,
$e^{-\mathbf{r}_{jj^{\prime}}\cdot\bm{\nabla}}$ with
$\mathbf{r}_{jj^{\prime}}=\mathbf{r}_{j}-\mathbf{r}_{j^{\prime}}$ is the
translation operator from site $j$ to site $j^{\prime}$.

To rewrite the interacting part in (1), we use the decomposition [18, 12, 14]

$n_{j\uparrow}n_{j\downarrow}=\frac{1}{4}(n_{j})^{2}-\frac{1}{4}(c^{*}_{j}\,\vec{\sigma}\cdot\hat{\Omega}_{j}\,c_{j})^{2},$
(7)

where $n_{j}=n_{j,\uparrow}+n_{j,\downarrow}$ is the charge density operator,
$\vec{\sigma}=(\sigma^{1},\sigma^{2},\sigma^{3})$ are the Pauli matrices, and
$\hat{\Omega}_{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

$\mathcal{S}_{\mathrm{int}}[\psi,\psi^{*},R]=\int_{0}^{\beta}\\!d\tau\,U\sum_{j}\left[\frac{1}{4}(n_{j}^{\psi})^{2}-\frac{1}{4}(\vec{S}^{\psi}_{j}\cdot\hat{\Omega}^{R}_{j})^{2}\right],$
(8)

where $n^{\psi}_{j}=\psi^{*}_{j}\psi_{j}$ is the chargon density operator,
$\vec{S}^{\psi}_{j}=\frac{1}{2}\psi^{*}_{j}\vec{\sigma}\psi_{j}$ is the
chargon spin operator, and $\hat{\Omega}^{R}_{j}$ is a unit vector obtained by
rotating $\hat{\Omega}_{j}$ as

$\vec{\sigma}\cdot\hat{\Omega}^{R}_{j}=R^{\dagger}_{j}\,\vec{\sigma}\cdot\hat{\Omega}_{j}\,R_{j}.$
(9)

Using (7) again, we obtain

$\mathcal{S}_{\mathrm{int}}[\psi,\psi^{*},R]=\int_{0}^{\beta}\\!d\tau\,U\sum_{j}n^{\psi}_{j\uparrow}n^{\psi}_{j\downarrow},$
(10)

with $n^{\psi}_{js}=\psi^{*}_{js}\psi_{js}$. Therefore, the final form of the
action $\mathcal{S}=\mathcal{S}_{0}+\mathcal{S}_{\mathrm{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 ($\langle\vec{S}^{\psi}_{j}\rangle\neq
0$) does not necessarily imply long range order for the electrons. The
matrices $R_{j}$ describe directional fluctuations of the order parameter
$\langle\vec{S}_{j}\rangle$, where, at low temperatures, the most important
ones vary slowly in space and time.

###  II.2 Nonlinear sigma model

We now derive a low energy effective action for the spinon fields $R_{j}$ by
integrating out the chargons,

$e^{-\mathcal{S}_{\mathrm{eff}}[R]}=\int\\!\mathcal{D}\psi\mathcal{D}\psi^{*}\,e^{-\mathcal{S}[\psi,\psi^{*},R]}.$
(11)

Since the action $\mathcal{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 $\partial_{\mu}R_{j}$, the effective action
$\mathcal{S}_{\mathrm{eff}}[R]$ has the general form

$\mathcal{S}_{\mathrm{eff}}[R]=\int_{\mathcal{T}}\\!dx\Big{[}\mathcal{B}^{a}_{\mu}A_{\mu}^{a}(x)+\textstyle{\frac{1}{2}}\mathcal{J}^{ab}_{\mu\nu}A_{\mu}^{a}(x)A_{\nu}^{b}(x)\Big{]},$
(12)

where $x=(\tau,\mathbf{r})$ combines imaginary time and space coordinates,
$\mathcal{T}=[0,\beta]\times\mathbb{R}^{2}$ is the integration region, and
repeated indices are summed. We have expanded the gauge field $A_{\mu}$ in
terms of the SU(2) generators,

$A_{\mu}(x)=A_{\mu}^{a}(x)\,\sigma^{a}/2,$ (13)

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

$\displaystyle\mathcal{B}^{a}_{\mu}$ $\displaystyle=$
$\displaystyle\frac{1}{2}\sum_{j,j^{\prime}}\gamma^{(1)}_{\mu}(j,j^{\prime})\langle\psi^{*}_{j}(0)\sigma^{a}\psi_{j^{\prime}}(0)\rangle,$
(14) $\displaystyle\mathcal{J}_{\mu\nu}^{ab}$ $\displaystyle=$
$\displaystyle-\frac{1}{4}\sum_{j,j^{\prime}}\sum_{l,l^{\prime}}\gamma^{(1)}_{\mu}(j,j^{\prime})\gamma^{(1)}_{\nu}(l,l^{\prime})$
(15)
$\displaystyle\times\int_{0}^{\beta}d\tau\,\big{\langle}\left(\psi^{*}_{j}(\tau)\sigma^{a}\psi_{j^{\prime}}(\tau)\right)\left(\psi^{*}_{l}(0)\sigma^{b}\psi_{l^{\prime}}(0)\right)\big{\rangle}_{c}$
$\displaystyle+$
$\displaystyle\frac{1}{4}\sum_{j,j^{\prime}}\gamma^{(2)}_{\mu\nu}(j,j^{\prime})\langle\psi^{*}_{j}(0)\psi_{j^{\prime}}(0)\rangle\,\delta_{ab},$

where $\langle\bullet\rangle$ ($\langle\bullet\rangle_{c}$) denotes the
(connected) average with respect to the chargon Hubbard action. The first and
second order current vertices have been defined as

$\displaystyle\gamma^{(1)}(j,j^{\prime})=$
$\displaystyle\left(\delta_{jj^{\prime}},i\,x_{jj^{\prime}}\,t_{jj^{\prime}},i\,y_{jj^{\prime}}\,t_{jj^{\prime}}\right),\hskip
28.45274pt$ (16a) $\displaystyle\gamma^{(2)}(j,j^{\prime})=$
$\displaystyle-\left(\begin{array}[]{ccc}0&0&0\\\
0&x_{jj^{\prime}}x_{jj^{\prime}}\,t_{jj^{\prime}}&x_{jj^{\prime}}y_{jj^{\prime}}\,t_{jj^{\prime}}\\\
0&y_{jj^{\prime}}x_{jj^{\prime}}\,t_{jj^{\prime}}&y_{jj^{\prime}}y_{jj^{\prime}}\,t_{jj^{\prime}}\\\
\end{array}\right),\hskip-14.22636pt$ (16e)

where $x_{jj^{\prime}}$ and $y_{jj^{\prime}}$ are the $x$ and $y$ components,
respectively of
$\mathbf{r}_{jj^{\prime}}=\mathbf{r}_{j}-\mathbf{r}_{j^{\prime}}$.

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 $\mathcal{R}$ of the SU(2) rotation $R$ via

$R^{\dagger}\,\sigma^{a}R=\mathcal{R}^{ab}\sigma^{b},$ (17)

we obtain the non-linear sigma model (NL$\sigma$M) action for the directional
fluctuations (see Appendix B)

$\mathcal{S}_{\mathrm{NL\sigma
M}}[\mathcal{R}]=\int_{\mathcal{T}}\\!dx\,\frac{1}{2}{\rm
Tr}\left[\mathcal{P}_{\mu\nu}(\partial_{\mu}\mathcal{R}^{T})(\partial_{\nu}\mathcal{R})\right],$
(18)

where $\mathcal{P}_{\mu\nu}=\frac{1}{2}{\rm
Tr}[\mathcal{J}_{\mu\nu}]\mathbb{1}-\mathcal{J}_{\mu\nu}$.

The structure of the matrices $\mathcal{J}_{\mu\nu}$ and
$\mathcal{P}_{\mu\nu}$ depends on the magnetically ordered chargon state. In
the trivial case $\langle\vec{S}^{\psi}_{j}\rangle=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,

$\langle\vec{S}^{\psi}_{j}\rangle\propto(-1)^{\bm{r}_{j}}\hat{u},$ (19)

where $\hat{u}$ is an arbitrary fixed unit vector. Choosing
$\hat{u}=\hat{e}_{1}=(1,0,0)$, the spin stiffness matrix in the Néel state has
the form

$\mathcal{J}_{\mu\nu}=\left(\begin{array}[]{ccc}0&0&0\\\ 0&J_{\mu\nu}&0\\\
0&0&J_{\mu\nu}\end{array}\right),$ (20)

with $(J_{\mu\nu})={\rm diag}(-Z,J,J)$. In this case the effective theory
reduces to the well-known ${\rm O(3)/O(2)}\simeq S_{2}$ non-linear sigma model
[20, 21]

$\mathcal{S}_{\mathrm{NL\sigma
M}}=\frac{1}{2}\int_{\mathcal{T}}dx\,\left(Z|\partial_{\tau}\hat{\Omega}|^{2}+J|\vec{\nabla}\hat{\Omega}|^{2}\right),$
(21)

where $\hat{\Omega}^{a}=\mathcal{R}^{a1}$, and $|\hat{\Omega}|^{2}=1$.

Another possibility is planar _spiral_ magnetic ordering of the chargons,

$\langle\vec{S}^{\psi}_{j}\rangle\propto\cos(\mathbf{Q}\cdot\mathbf{r}_{j})\hat{u}_{1}+\sin(\mathbf{Q}\cdot\mathbf{r}_{j})\hat{u}_{2},$
(22)

where $\mathbf{Q}$ is a fixed wave vector as obtained by minimizing the
chargon free energy, while $\hat{u}_{1}$ and $\hat{u}_{2}$ are two arbitrary
mutually orthogonal unit vectors. The special case $\mathbf{Q}=(\pi,\pi)$
corresponds to the Néel state. Fixing $\hat{u}_{1}$ to $\hat{e}_{1}$ and
$\hat{u}_{2}$ to $\hat{e}_{2}\equiv(0,1,0)$, the spin stiffness matrix assumes
the form

$\mathcal{J}_{\mu\nu}=\left(\begin{array}[]{ccc}J_{\mu\nu}^{\perp}&0&0\\\
0&J_{\mu\nu}^{\perp}&0\\\ 0&0&J_{\mu\nu}^{\Box}\end{array}\right),$ (23)

where

$(J_{\mu\nu}^{a})=\left(\begin{array}[]{ccc}-Z^{a}&0&0\\\
0&J_{xx}^{a}&J_{xy}^{a}\\\ 0&J_{yx}^{a}&J_{yy}^{a}\end{array}\right).$ (24)

for $a\in\\{\perp,\Box\\}$. In this case, the effective action maintains its
general form (18) and it describes the O(3)$\times$O(2)/O(2) symmetric
NL$\sigma$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
$\mathcal{J}_{\mu\nu}$. 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. [33, 19, 34] for
the fRG, and to Refs. [16, 30, 35, 36] for the fRG+MF method.

###  III.1 Symmetric regime

We evaluate the chargon functional integral by using an fRG flow equation [33,
19, 34], choosing the temperature $T$ as flow parameter [37]. Temperature can
be used as a flow parameter after rescaling the chargon fields as $\psi_{j}\to
T^{\frac{3}{4}}\psi_{j}$, and defining a rescaled bare Green’s function,
$G_{0}^{T}(\mathbf{k},i\nu_{n})=T^{\frac{1}{2}}/(i\nu_{n}-\epsilon_{\mathbf{k}}+\mu)$,
where $\nu_{n}=(2n+1)\pi T$ is the fermionic Matsubara frequency, and
$\epsilon_{\mathbf{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 three-particle vertex [19]. In a SU(2) invariant system the two-
particle vertex has the spin structure

$\begin{split}V^{T}_{\sigma_{1}\sigma_{2}\sigma_{3}\sigma_{4}}(k_{1},k_{2},k_{3},k_{4})&=V^{T}(k_{1},k_{2},k_{3},k_{4})\,\delta_{\sigma_{1}\sigma_{3}}\,\delta_{\sigma_{2}\sigma_{4}}\\\
&-V^{T}(k_{2},k_{1},k_{3},k_{4})\,\delta_{\sigma_{1}\sigma_{4}}\,\delta_{\sigma_{2}\sigma_{3}},\end{split}$

where $k_{\alpha}=(\mathbf{k}_{\alpha},i\nu_{\alpha 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]

$\displaystyle
V^{T}(\mathbf{k}_{1},\mathbf{k}_{2},\mathbf{k}_{3},\mathbf{k}_{4})=U-\phi^{p,T}_{\frac{\mathbf{k}_{1}-\mathbf{k}_{2}}{2},\frac{\mathbf{k}_{3}-\mathbf{k}_{4}}{2}}(\mathbf{k}_{1}+\mathbf{k}_{2})$
$\displaystyle+\,\phi^{m,T}_{\frac{\mathbf{k}_{1}+\mathbf{k}_{4}}{2},\frac{\mathbf{k}_{2}+\mathbf{k}_{3}}{2}}(\mathbf{k}_{2}-\mathbf{k}_{3})+\frac{1}{2}\phi^{m,T}_{\frac{\mathbf{k}_{1}+\mathbf{k}_{3}}{2},\frac{\mathbf{k}_{2}+\mathbf{k}_{4}}{2}}(\mathbf{k}_{3}-\mathbf{k}_{1})$
$\displaystyle-\frac{1}{2}\phi^{c,T}_{\frac{\mathbf{k}_{1}+\mathbf{k}_{3}}{2},\frac{\mathbf{k}_{2}+\mathbf{k}_{4}}{2}}(\mathbf{k}_{3}-\mathbf{k}_{1}),$
(25)

where the functions $\phi^{p,T}$, $\phi^{m,T}$, and $\phi^{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_{\mathrm{ini}}=\infty$, at
which $V^{T_{\mathrm{ini}}}=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 $\phi^{m,T}$, the chargons develop
some kind of magnetic order.

###  III.2 Order parameter

In the magnetic phase, that is, for $T<T^{*}$, we assume an order parameter of
the form
$\langle\psi^{*}_{\mathbf{k},\uparrow}\psi_{\mathbf{k}+\mathbf{Q},\downarrow}\rangle$,
which corresponds to Néel antiferromagnetism if $\mathbf{Q}=(\pi,\pi)$, and to
spiral order otherwise.

For $T<T^{*}$ we simplify the flow equations by decoupling the three channels
$\phi^{p,T}$, $\phi^{m,T}$, and $\phi^{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]

$\Delta_{\mathbf{k}}=\int_{\mathbf{k}^{\prime}}\bar{V}^{m}_{\mathbf{k},\mathbf{k}^{\prime}}(\mathbf{Q})\,\frac{f(E^{-}_{\mathbf{k}^{\prime}})-f(E^{+}_{\mathbf{k}^{\prime}})}{E^{+}_{\mathbf{k}^{\prime}}-E^{-}_{\mathbf{k}^{\prime}}}\,\Delta_{\mathbf{k}^{\prime}},$
(26)

where $f(x)=(e^{x/T}+1)^{-1}$ is the Fermi function, $\int_{\mathbf{k}}$ is a
shorthand notation for $\int\\!\frac{d^{2}\mathbf{k}}{(2\pi)^{2}}$, and
$E^{\pm}_{\mathbf{k}}$ are the quasiparticle dispersions

$E^{\pm}_{\mathbf{k}}=\frac{\epsilon_{\mathbf{k}}+\epsilon_{\mathbf{k}+\mathbf{Q}}}{2}\pm\sqrt{\frac{1}{4}\left(\epsilon_{\mathbf{k}}-\epsilon_{\mathbf{k}+\mathbf{Q}}\right)^{2}+\Delta_{\mathbf{k}}^{2}}\,-\mu.$
(27)

The effective coupling
$\bar{V}^{m}_{\mathbf{k},\mathbf{k}^{\prime}}(\mathbf{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,

$\begin{split}V^{m,T^{*}}_{\mathbf{k},\mathbf{k}^{\prime}}(\mathbf{q})&=\bar{V}^{m}_{\mathbf{k},\mathbf{k}^{\prime}}(\mathbf{q})\\\
&-\int_{\mathbf{k}^{\prime\prime}}\bar{V}^{m}_{\mathbf{k},\mathbf{k}^{\prime\prime}}(\mathbf{q})\,\Pi^{T^{*}}_{\mathbf{k}^{\prime\prime}}(\mathbf{q})\,V^{m,T^{*}}_{\mathbf{k}^{\prime\prime},\mathbf{k}^{\prime}}(\mathbf{q}),\end{split}$
(28)

where
$V^{m,T}_{\mathbf{k},\mathbf{k}^{\prime}}(\mathbf{q})=V^{T}(\mathbf{k}-\mathbf{q}/2,\mathbf{k}^{\prime}+\mathbf{q}/2,\mathbf{k}^{\prime}-\mathbf{q}/2,\mathbf{k}+\mathbf{q}/2)$,
and the particle-hole bubble is given by

$\Pi^{T}_{\mathbf{k}}(\mathbf{q})=\sum_{\nu_{n}}G_{0}^{T}\left(\mathbf{k}-\mathbf{q}/2,i\nu_{n}\right)G_{0}^{T}\left(\mathbf{k}+\mathbf{q}/2,i\nu_{n}\right).$
(29)

Although $V^{m,T^{*}}_{\mathbf{k},\mathbf{k}^{\prime}}(\mathbf{q})$ diverges
at certain wave vectors $\mathbf{q}=\mathbf{Q}_{c}$, the irreducible coupling
$\bar{V}^{m}_{\mathbf{k},\mathbf{k}^{\prime}}(\mathbf{q})$ is finite for all
$\mathbf{q}$.

The dependence of $\bar{V}^{m}_{\mathbf{k},\mathbf{k}^{\prime}}(\mathbf{q})$
on $\mathbf{k}$ and $\mathbf{k}^{\prime}$ is rather weak and of no qualitative
importance. Hence, to simplify the calculations, we discard the $\mathbf{k}$
and $\mathbf{k}^{\prime}$ dependencies of the effective coupling by taking the
momentum average
$\bar{V}^{m}(\mathbf{q})=\int_{\mathbf{k},\mathbf{k}^{\prime}}\bar{V}^{m}_{\mathbf{k},\mathbf{k}^{\prime}}(\mathbf{q})$.
The magnetic gap then becomes momentum independent, that is,
$\Delta_{\mathbf{k}}=\Delta$. While the full vertex
$V^{m,T}_{\mathbf{k},\mathbf{k}^{\prime}}(\mathbf{q})$ depends very strongly
on $\mathbf{q}$, the dependence of its irreducible part
$\bar{V}_{\mathbf{k},\mathbf{k}^{\prime}}(\mathbf{q})$ on $\mathbf{q}$ is
rather weak. The calculation of the stiffnesses in the subsequent section is
considerably simplified approximating $\bar{V}^{m}(\mathbf{q})$ by a momentum
independent effective interaction $U_{\rm
eff}^{m}=\bar{V}^{m}(\mathbf{Q}_{c})$.

The optimal ordering wave vector $\mathbf{Q}$ is found by minimizing the mean-
field free energy of the system

$F(\mathbf{Q})=-T\int_{\mathbf{k}}\sum_{\ell=\pm}\ln\left(1+e^{-E^{\ell}_{\mathbf{k}}(\mathbf{Q})/T}\right)+\frac{\Delta^{2}}{2U_{\rm
eff}^{m}}+\mu n,$ (30)

where the chemical potential $\mu$ is determined by keeping the density
$n=\int_{\mathbf{k}}\sum_{\ell=\pm}f(E^{\ell}_{\mathbf{k}})$ fixed. The
optimal wave vectors $\mathbf{Q}$ at temperatures $T<T^{*}$ generally differ
from the wave vectors $\mathbf{Q}_{c}$ at which
$V^{T^{*}}_{\mathbf{k},\mathbf{k}^{\prime}}(\mathbf{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]).

###  III.3 Spin stiffnesses

The NL$\sigma$M parameters, that is, the spin stiffnesses
$\mathcal{J}_{\mu\nu}^{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 [43, 17]. 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_{\rm eff}^{m}$ and $U_{\rm eff}^{c}$. Since in
the magnetic state the translational symmetry is broken, the Fourier
transforms of the response functions depend on two distinct momenta
$\mathbf{q}$ and $\mathbf{q}^{\prime}$, where $\mathbf{q}^{\prime}$ can assume
the values $\mathbf{q}$, $\mathbf{q}\pm\mathbf{Q}$, and $\mathbf{q}\pm
2\mathbf{Q}$. However, to compute $\mathcal{J}_{\mu\nu}^{ab}$, we only need to
deal with the limit $\mathbf{q},\mathbf{q}^{\prime}\to\mathbf{0}$.

The temporal components of the stiffness matrix, that is,
$\mathcal{J}_{00}^{ab}$, are given by the uniform spin susceptibility in the
dynamical limit [17]

$\mathcal{J}_{00}^{ab}=-\chi_{\rm dyn}^{ab}=-\lim_{\omega\to
0}\chi^{ab}(\mathbf{0},\mathbf{0},\omega),$ (31)

where $\chi^{ab}(\mathbf{q},\mathbf{q}^{\prime},\omega)$ is the Fourier
transform of

$\chi_{jl}^{ab}(\tau)=\frac{1}{4}\big{\langle}\left(\psi^{*}_{j}(\tau)\sigma^{a}\psi_{j}(\tau)\right)\left(\psi^{*}_{l}(0)\sigma^{b}\psi_{l}(0)\right)\big{\rangle}_{c}.$
(32)

Note that, in a metallic system, the static uniform susceptibility obtained
from $\mathbf{q},\mathbf{q}^{\prime}\to\mathbf{0}$ after setting $\omega=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 [44, 43]

$\widetilde{\psi}_{j}=e^{-i\frac{\mathbf{Q}}{2}\cdot\mathbf{r}_{j}}e^{i\sigma^{3}\frac{\mathbf{Q}}{2}\cdot\mathbf{r}_{j}}\psi_{j},\quad\widetilde{\psi}^{*}_{j}=\psi^{*}_{j}\,e^{-i\sigma^{3}\frac{\mathbf{Q}}{2}\cdot\mathbf{r}_{j}}e^{i\frac{\mathbf{Q}}{2}\cdot\mathbf{r}_{j}},$
(33)

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 $\widetilde{\chi}^{ab}(\mathbf{q},\omega)$, with a
single momentum variable $\mathbf{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 $\mathbf{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 $\sigma^{0}$ as the two-
dimensional unit matrix. The prefactor $\frac{1}{4}$ in Eq. (32) implies that
$\chi^{00}$ is actually a quarter of the conventional charge susceptibility.

In RPA, the rotated susceptibility $\widetilde{\chi}$ can be written as

$\widetilde{\chi}^{ab}(q)=\widetilde{\chi}_{0}^{ab}(q)+\sum_{a^{\prime},b^{\prime}=0}^{3}\widetilde{\chi}_{0}^{aa^{\prime}}(q)\widetilde{\Gamma}^{a^{\prime}b^{\prime}}(q)\widetilde{\chi}_{0}^{b^{\prime}b}(q),$
(34)

where $q=(\mathbf{q},\omega)$, and $\widetilde{\Gamma}^{ab}(q)$ is the RPA
effective interaction in the rotated spin frame. The “bare” susceptibility
$\widetilde{\chi}_{0}^{ab}(q)$ is given by the particle-hole bubble [43]

$\displaystyle\widetilde{\chi}_{0}^{ab}(\mathbf{q},\omega)$ $\displaystyle=$
$\displaystyle-\frac{1}{4}\int_{\mathbf{k}}T\sum_{\nu_{n}}{\rm
Tr}\big{[}\sigma^{a}\widetilde{\mathcal{G}}(\mathbf{k}+\mathbf{q},i\nu_{n}+\omega+i0^{+})$
(35)
$\displaystyle\times\,\sigma^{b}\widetilde{\mathcal{G}}(\mathbf{k},i\nu_{n})\big{]}\,,$

with $\widetilde{\mathcal{G}}(\mathbf{k},i\nu_{n})$ the mean-field chargon
Green’s function in the rotated basis

$\widetilde{\mathcal{G}}(\mathbf{k},i\nu_{n})=\left(\begin{array}[]{cc}i\nu_{n}-\epsilon_{\mathbf{k}}+\mu&-\Delta\\\
-\Delta&i\nu_{n}-\epsilon_{\mathbf{k}+\mathbf{Q}}+\mu\end{array}\right)^{-1}.$
(36)

The RPA effective interaction is obtained from a ladder sum, leading to the
linear matrix equation

$\widetilde{\Gamma}^{ab}(q)=\widetilde{\Gamma}^{ab}_{0}(\mathbf{q})+\sum_{a^{\prime},b^{\prime}=0}^{3}\widetilde{\Gamma}^{aa^{\prime}}_{0}(\mathbf{q})\,\widetilde{\chi}^{a^{\prime}b^{\prime}}_{0}(q)\,\widetilde{\Gamma}^{b^{\prime}b}(q),$
(37)

where

$\widetilde{\Gamma}^{ab}_{0}(\mathbf{q})=\Gamma^{ab}_{0}(\mathbf{q})=2\,{\rm
diag}\left[-U_{\rm eff}^{c}(\mathbf{q}),U_{\rm eff}^{m},U_{\rm eff}^{m},U_{\rm
eff}^{m}\right].$ (38)

The effective charge interaction is given by $U_{\rm
eff}^{c}(\mathbf{q})=\int_{\mathbf{k},\mathbf{k}^{\prime}}\bar{V}^{c}_{\mathbf{k},\mathbf{k}^{\prime}}(\mathbf{q})$,
where the irreducible coupling
$\bar{V}^{c}_{\mathbf{k},\mathbf{k}^{\prime}}(\mathbf{q})$ is obtained by
inverting a Bethe-Salpeter equation similar to Eq. (28),

$\begin{split}V^{c,T^{*}}_{\mathbf{k},\mathbf{k}^{\prime}}(\mathbf{q})&=\bar{V}^{c}_{\mathbf{k},\mathbf{k}^{\prime}}(\mathbf{q})\\\
&+\int_{\mathbf{k}^{\prime\prime}}\bar{V}^{c}_{\mathbf{k},\mathbf{k}^{\prime\prime}}(\mathbf{q})\,\Pi^{T^{*}}_{\mathbf{k}^{\prime\prime}}(\mathbf{q})\,V^{c,T^{*}}_{\mathbf{k}^{\prime\prime},\mathbf{k}^{\prime}}(\mathbf{q}),\end{split}$
(39)

with

$\displaystyle V^{c,T}_{\mathbf{k},\mathbf{k}^{\prime}}(\mathbf{q})$
$\displaystyle=$ $\displaystyle
2V^{T}(\mathbf{k}\\!-\\!\mathbf{q}/2,\mathbf{k}^{\prime}\\!+\\!\mathbf{q}/2,\mathbf{k}\\!+\\!\mathbf{q}/2,\mathbf{k}^{\prime}\\!-\\!\mathbf{q}/2)$
$\displaystyle-$ $\displaystyle
V^{T}(\mathbf{k}\\!-\\!\mathbf{q}/2,\mathbf{k}^{\prime}\\!+\\!\mathbf{q}/2,\mathbf{k}^{\prime}\\!-\\!\mathbf{q}/2,\mathbf{k}\\!+\\!\mathbf{q}/2).$

Here we keep the dependence on $\mathbf{q}$ since it does not complicate the
calculations. The off-diagonal ($a\neq b$) elements of
$\chi^{ab}(\mathbf{0},\mathbf{0},\omega)$ with $a,b=1,2,3$ vanish for
$\omega\to 0$ both in the spiral and in the Néel state, so that we need to
deal only with the diagonal spin susceptibility components
$\chi^{aa}(\mathbf{q},\mathbf{q},\omega)$.

In a spiral state with $\mathbf{Q}\neq(\pi,\pi)$, the diagonal (both in
momentum and spin indices) spin susceptibility components are related to the
susceptibility components in the rotated basis as [43]

$\displaystyle\chi^{11}(\mathbf{q},\mathbf{q},\omega)$
$\displaystyle=\chi^{22}(\mathbf{q},\mathbf{q},\omega)$
$\displaystyle=\frac{1}{4}\sum_{s=\pm}\big{[}\widetilde{\chi}^{11}(q+sQ)+\widetilde{\chi}^{22}(q+sQ)\big{]}$
$\displaystyle\quad+2i\widetilde{\chi}^{12}(q+Q)+2i\widetilde{\chi}^{21}(q-Q),$
(40a) $\displaystyle\chi^{33}(\mathbf{q},\mathbf{q},\omega)$
$\displaystyle=\widetilde{\chi}^{33}(q),$ (40b)

where $Q=(\mathbf{Q},0)$. The momentum diagonal components of
$\chi^{11}(\mathbf{q},\mathbf{q}^{\prime},\omega)$ and
$\chi^{22}(\mathbf{q},\mathbf{q}^{\prime},\omega)$ are equal, and the limit
$\omega\to 0$ in (31) is nonzero for all diagonal components, yielding

$\mathcal{J}_{00}=\left(\begin{array}[]{ccc}-Z^{\perp}&0&0\\\
0&-Z^{\perp}&0\\\ 0&0&-Z^{\Box}\end{array}\right),$ (41)

with $Z^{\perp}=\chi^{22}(\mathbf{0},\mathbf{0},\omega\to 0)$ and
$Z^{\Box}=\chi^{33}(\mathbf{0},\mathbf{0},\omega\to 0)$. The quantities
$Z^{\Box}$ and $Z^{\perp}$ parametrize the low frequency dependence of the in-
plane and out-of-plane spin susceptibility, respectively, near the Goldstone
poles [43, 17].

In the limits $\omega\to 0$ and $\mathbf{q}\to\mathbf{0}$ or $\mathbf{Q}$,
several off-diagonal matrix elements of the RPA effective interaction
$\widetilde{\Gamma}^{ab}(\mathbf{q},\omega)$ vanish [43]. The expressions for
$Z^{a}$ can therefore be simplified to [17]

$Z^{\Box}=\frac{\widetilde{\chi}_{0}^{33}(\mathbf{0},\omega\to 0)}{1-2U_{\rm
eff}^{m}\,\widetilde{\chi}_{0}^{33}(\mathbf{0},\omega\to 0)}\,,$ (42)

and

$\displaystyle Z^{\perp}$ $\displaystyle=$ $\displaystyle
2\widetilde{\chi}_{0}^{-+}(\mathbf{Q},0)$ (43) $\displaystyle+$ $\displaystyle
2\sum_{a,b=0,1,2}\widetilde{\chi}_{0}^{-a}(\mathbf{Q},0)\widetilde{\Gamma}^{ab}(\mathbf{Q},0)\widetilde{\chi}_{0}^{b+}(\mathbf{Q},0)\,,\hskip
14.22636pt$

where superscripts $+$ and $-$ attached to
$\widetilde{\chi}_{0}^{\phantom{a}}$ indicate that the susceptibilities are
formed with the ladder operators $S^{\pm}=\frac{1}{2}\left(S^{1}\pm
S^{2}\right)$.

In the Néel state, terms which contribute to
$\chi^{11}(\mathbf{q},\mathbf{q}^{\prime},\omega)$ and
$\chi^{22}(\mathbf{q},\mathbf{q}^{\prime},\omega)$ with
$\mathbf{q}^{\prime}=\mathbf{q}\pm 2\mathbf{Q}$ for $\mathbf{Q}\neq(\pi,\pi)$,
contribute to the momentum diagonal susceptibilities since
$2\mathbf{Q}\equiv\mathbf{0}$ for $\mathbf{Q}=(\pi,\pi)$. The transformation
of the susceptibilities from the rotated to the unrotated basis then reads
[43]

$\displaystyle\chi^{11}(\mathbf{q},\mathbf{q},\omega)=\widetilde{\chi}^{11}(q+Q),$
(44a)
$\displaystyle\chi^{22}(\mathbf{q},\mathbf{q},\omega)=\widetilde{\chi}^{22}(q+Q),$
(44b)
$\displaystyle\chi^{33}(\mathbf{q},\mathbf{q},\omega)=\widetilde{\chi}^{33}(q).$
(44c)

Since $\widetilde{\chi}^{11}(\mathbf{Q},\omega)=0$ and
$\widetilde{\chi}^{22}(q+Q)=\widetilde{\chi}^{33}(q)$ for
$\mathbf{Q}=(\pi,\pi)$, we obtain

$\mathcal{J}_{00}=\left(\begin{array}[]{ccc}0&0&0\\\ 0&-Z&0\\\
0&0&-Z\end{array}\right),$ (45)

with $Z=\chi^{22}(\mathbf{0},\mathbf{0},\omega\to 0)$ in the Néel state, which
can be evaluated by the same expression as the one for $Z^{\Box}$ in Eq. (42).

Figure 1: Diagrams contributing to the spin stiffness. The black triangles and
circles represent the first and second order current vertices
$\gamma_{\alpha}^{(1)}(\mathbf{k})$ and
$\gamma_{\alpha\beta}^{(2)}(\mathbf{k})$, respectively, while the dashed line
represents the effective interaction
$\Gamma^{ab}(\mathbf{q},\mathbf{q}^{\prime},\omega)$.

The spatial components of the stiffness matrix
$\mathcal{J}_{\alpha\beta}^{ab}$ with $\alpha,\beta=1,2$ are obtained from the
spatial components of the uniform gauge field kernel
$\mathcal{K}_{\alpha\beta}^{ab}$ in the static limit [17]

$\mathcal{J}_{\alpha\beta}^{ab}=-\lim_{\mathbf{q}\to\mathbf{0}}\mathcal{K}_{\alpha\beta}^{ab}(\mathbf{q},\mathbf{q},0)\,,$
(46)

where
$\mathcal{K}_{\alpha\beta}^{ab}(\mathbf{q},\mathbf{q}^{\prime},\omega)=\mathcal{K}_{\alpha\beta}^{{\rm
p},ab}(\mathbf{q},\mathbf{q}^{\prime},\omega)+\delta_{ab}\mathcal{K}_{\alpha\beta}^{\rm
d}$ is the Fourier transform of

$\mathcal{K}_{\alpha\beta,jl}^{ab}(\tau)=\mathcal{K}_{\alpha\beta,jl}^{{\rm
p},ab}(\tau)+\mathcal{K}_{\alpha\beta}^{\rm
d}\,\delta_{ab}\delta_{jl}\delta(\tau),$ (47)

with the paramagnetic and diamagnetic contributions (cf. Eq. (15))

$\displaystyle\mathcal{K}_{\alpha\beta,jl}^{{\rm p},ab}(\tau)$
$\displaystyle=$
$\displaystyle\frac{1}{4}\sum_{j^{\prime},l^{\prime}}\gamma_{\alpha}^{(1)}(j,j^{\prime})\,\gamma_{\beta}^{(1)}(l,l^{\prime})$
(48) $\displaystyle\times$
$\displaystyle\big{\langle}\left(\psi^{*}_{j}(\tau)\sigma^{a}\psi_{j^{\prime}}(\tau)\right)\left(\psi^{*}_{l}(0)\sigma^{b}\psi_{l^{\prime}}(0)\right)\big{\rangle}_{c}\,,\hskip
17.07164pt$ $\displaystyle\mathcal{K}_{\alpha\beta}^{\rm d}$ $\displaystyle=$
$\displaystyle-\frac{1}{4}\sum_{j^{\prime}}\gamma_{\alpha\beta}^{(2)}(j,j^{\prime})\,\big{\langle}\psi^{*}_{j}(0)\psi_{j^{\prime}}(0)\big{\rangle}\,,$
(49)

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

$\displaystyle\mathcal{J}_{\alpha\beta}^{{\rm p},ab}$ $\displaystyle=$
$\displaystyle-\lim_{\mathbf{q}\to\mathbf{0}}\Big{[}\mathcal{K}_{0,\alpha\beta}^{{\rm
p},ab}(\mathbf{q},\mathbf{q},0)$ (50) $\displaystyle+$
$\displaystyle\sum_{a^{\prime},b^{\prime}=0}^{3}\int_{\mathbf{q}^{\prime},\mathbf{q}^{\prime\prime}}\mathcal{K}_{0,\alpha
0}^{{\rm p},aa^{\prime}}(\mathbf{q},\mathbf{q}^{\prime},0)$
$\displaystyle\times\,\Gamma^{a^{\prime}b^{\prime}}(\mathbf{q}^{\prime},\mathbf{q}^{\prime\prime},0)\,\mathcal{K}_{0,0\beta}^{{\rm
p},b^{\prime}b}(\mathbf{q}^{\prime\prime},\mathbf{q},0)\Big{]},$

where $\Gamma^{ab}(\mathbf{q},\mathbf{q}^{\prime},\omega)$ is the RPA
effective interaction (37) in the original (non-rotated) spin basis. The bare
paramagnetic response kernel is given by

$\displaystyle\mathcal{K}_{0,\mu\nu}^{{\rm
p},ab}(\mathbf{q},\mathbf{q}^{\prime},0)$ $\displaystyle=$
$\displaystyle-\frac{1}{4}\int_{\mathbf{k},\mathbf{k}^{\prime}}T\sum_{\nu_{n}}\gamma_{\mu}^{(1)}(\mathbf{k}+\mathbf{q})\,\gamma_{\nu}^{(1)}(\mathbf{k}^{\prime})$
(51) $\displaystyle\hskip-56.9055pt\times\,{\rm
Tr}\left[\sigma^{a}\mathcal{G}(\mathbf{k}+\mathbf{q},\mathbf{k}^{\prime}+\mathbf{q}^{\prime},i\nu_{n})\sigma^{b}\mathcal{G}(\mathbf{k}^{\prime},\mathbf{k},i\nu_{n})\right],\hskip
19.91692pt$

where
$\gamma^{(1)}(\mathbf{k})=(1,\partial_{k_{x}}\epsilon_{\mathbf{k}},\partial_{k_{y}}\epsilon_{\mathbf{k}})$
is the Fourier transform of $\gamma^{(1)}(j,j^{\prime})$ in Eq. (16a). The
chargon Green’s function in the original spin basis reads

$\mathcal{G}(\mathbf{k},\mathbf{k}^{\prime},i\nu_{n})=\left(\begin{array}[]{cc}G_{+\mathbf{Q}}(k)\,\delta_{\mathbf{k},\mathbf{k}^{\prime}}&F_{+\mathbf{Q}}(k)\,\delta_{\mathbf{k}+\mathbf{Q},\mathbf{k}^{\prime}}\\\
F_{-\mathbf{Q}}(k)\,\delta_{\mathbf{k}-\mathbf{Q},\mathbf{k}^{\prime}}&G_{-\mathbf{Q}}(k)\,\delta_{\mathbf{k},\mathbf{k}^{\prime}}\end{array}\right),$
(52)

with

$\displaystyle
G_{\pm\mathbf{Q}}(k)=\frac{i\nu_{n}-\epsilon_{\mathbf{k}\pm\mathbf{Q}}+\mu}{(i\nu_{n}\\!-\\!\epsilon_{\mathbf{k}}\\!+\\!\mu)(i\nu_{n}\\!-\\!\epsilon_{\mathbf{k}\pm\mathbf{Q}}\\!+\\!\mu)-\Delta^{2}},$
(53a) $\displaystyle
F_{\pm\mathbf{Q}}(k)=\frac{\Delta}{(i\nu_{n}\\!-\\!\epsilon_{\mathbf{k}}\\!+\\!\mu)(i\nu_{n}\\!-\\!\epsilon_{\mathbf{k}\pm\mathbf{Q}}\\!+\\!\mu)-\Delta^{2}}.$
(53b)

The diamagnetic part of the spin stiffness can be obtained from the Green’s
function as

$\mathcal{J}_{\alpha\beta}^{{\rm
d},ab}=\frac{\delta_{ab}}{4}\int_{\mathbf{k},\mathbf{k}^{\prime}}T\sum_{\nu_{n}}\gamma_{\alpha\beta}^{(2)}(\mathbf{k})\,{\rm
tr}\left[\mathcal{G}(\mathbf{k},\mathbf{k}^{\prime},i\nu_{n})\right]\,$ (54)

where
$\gamma_{\alpha\beta}^{(2)}(\mathbf{k})=\partial_{k_{\alpha}}\partial_{k_{\beta}}\epsilon_{\mathbf{k}}$
is the Fourier transform of the second order current vertex
$\gamma_{\alpha\beta}^{(2)}(j,j^{\prime})$ in Eq. (16e). The various
contributions to the response kernel are represented diagrammatically in Fig.
1.

The off-diagonal ($a\neq b$) components of
$\mathcal{K}_{\alpha\beta}^{ab}(\mathbf{q},\mathbf{q},0)$ vanish for
$\mathbf{q}\to\mathbf{0}$ in the spiral and in the Néel state, so that we need
to consider only the diagonal components. For $a=b=1,2$ only the first (bare)
term in Eq. (50) contributes to the stiffness [17]. In a spiral state with
$\mathbf{Q}\neq(\pi,\pi)$, one thus obtains the out-of-plane stiffness as

$J_{\alpha\beta}^{\perp}=\mathcal{J}_{\alpha\beta}^{11}=\mathcal{J}_{\alpha\beta}^{22}=-\lim_{\mathbf{q}\to\mathbf{0}}\mathcal{K}_{0,\alpha\beta}^{{\rm
p},22}(\mathbf{q},\mathbf{q},0)-\mathcal{K}_{\alpha\beta}^{\rm d}\,.$ (55)

For $a=b=3$, there are non vanishing components of the kernel that mix
temporal and spatial indices, namely

$\displaystyle\mathcal{K}_{0,\alpha 0}^{{\rm
p},30}(\mathbf{q},\mathbf{q}^{\prime},0)$
$\displaystyle=\mathcal{K}_{0,0\alpha}^{{\rm
p},03}(\mathbf{q},\mathbf{q}^{\prime},0)=\mathcal{K}_{0,\alpha 0}^{{\rm
p},30}(\mathbf{q},0)\,\delta_{\mathbf{q},\mathbf{q}^{\prime}}\,,$ (56a)
$\displaystyle\mathcal{K}_{0,\alpha 0}^{{\rm
p},31}(\mathbf{q},\mathbf{q}^{\prime},0)$
$\displaystyle=\mathcal{K}_{0,0\alpha}^{{\rm
p},13}(\mathbf{q},\mathbf{q}^{\prime},0)$ $\displaystyle=\mathcal{K}_{0,\alpha
0}^{{\rm
p},31}(\mathbf{q},0)\,\frac{\delta_{\mathbf{q}+\mathbf{Q},\mathbf{q}^{\prime}}+\delta_{\mathbf{q}-\mathbf{Q},\mathbf{q}^{\prime}}}{2}\,,$
(56b) $\displaystyle\mathcal{K}_{0,\alpha 0}^{{\rm
p},32}(\mathbf{q},\mathbf{q}^{\prime},0)$
$\displaystyle=\mathcal{K}_{0,0\alpha}^{{\rm
p},23}(\mathbf{q},\mathbf{q}^{\prime},0)$ $\displaystyle=\mathcal{K}_{0,\alpha
0}^{{\rm
p},32}(\mathbf{q},0)\,\frac{\delta_{\mathbf{q}+\mathbf{Q},\mathbf{q}^{\prime}}-\delta_{\mathbf{q}-\mathbf{Q},\mathbf{q}^{\prime}}}{2i}\,.$
(56c)

Using $\mathcal{K}_{0,\alpha 0}^{{\rm
p},31}(\mathbf{q}\to\mathbf{0},0)=\mathcal{K}_{0,\alpha 0}^{{\rm
p},32}(\mathbf{q}\to\mathbf{0},0)$, one thus obtains the in-plane stiffness in
the form [17]

$\displaystyle J_{\alpha\beta}^{\Box}$ $\displaystyle=$
$\displaystyle\mathcal{J}_{\alpha\beta}^{33}=-\lim_{\mathbf{q}\to\mathbf{0}}\mathcal{K}_{0,\alpha\beta}^{{\rm
p},33}(\mathbf{q},\mathbf{q},0)-\mathcal{K}_{\alpha\beta}^{\rm d}$ (57)
$\displaystyle-$ $\displaystyle\lim_{\mathbf{q}\to
0}\sum_{a,b=0,1}\mathcal{K}_{0,\alpha 0}^{{\rm
p},3a}(\mathbf{q},0)\,\widetilde{\Gamma}^{ab}(\mathbf{q},0)\,\mathcal{K}_{0,\beta
0}^{{\rm p},3b}(\mathbf{q},0)\,,\hskip 19.91692pt$

where $\widetilde{\Gamma}^{ab}(q)$ is the effective interaction in the spin
rotated basis, see Eq. (37).

In the Néel state one finds, in close analogy to the temporal components
$\mathcal{J}_{00}^{ab}$ of the stiffness matrix,
$\mathcal{J}_{\alpha\beta}^{11}=0$ and
$\mathcal{J}_{\alpha\beta}^{22}=\mathcal{J}_{\alpha\beta}^{33}=J\delta_{\alpha\beta}$,
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 $\mathbf{Q}$) of the reconstructed Fermi surface. In our large
$N$ evaluation of the NL$\sigma$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_{\rm eff}^{m}$ 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$\sigma$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$\sigma$M, we resort to a saddle point approximation in the
$\text{CP}^{N-1}$ representation, which is exact in the large $N$ limit [46,
47].

###  IV.1 CP1 representation

The matrix $\mathcal{R}$ can be expressed as a triad of orthonormal unit
vectors:

$\mathcal{R}=\big{(}\hat{\Omega}_{1},\hat{\Omega}_{2},\hat{\Omega}_{3}\big{)},$
(58)

where $\hat{\Omega}_{i}\cdot\hat{\Omega}_{j}=\delta_{ij}$. We represent these
vectors in terms of two complex Schwinger bosons $z_{\uparrow}$ and
$z_{\downarrow}$ [45]

$\displaystyle\hat{\Omega}_{-}=z(i\sigma^{2}\vec{\sigma})z,$ (59a)
$\displaystyle\hat{\Omega}_{+}=z^{*}(i\sigma^{2}\vec{\sigma})^{\dagger}z^{*},$
(59b) $\displaystyle\hat{\Omega}_{3}=z^{*}\vec{\sigma}z,$ (59c)

with $z=(z_{\uparrow},z_{\downarrow})$ and
$\hat{\Omega}_{\pm}=\hat{\Omega}_{1}\mp i\hat{\Omega}_{2}$. The Schwinger
bosons obey the nonlinear constraint

$z^{*}_{\uparrow}z_{\uparrow}+z^{*}_{\downarrow}z_{\downarrow}=1\,.$ (60)

The parametrization (59) is equivalent to

$R=\left(\begin{array}[]{cc}z_{\uparrow}&-z_{\downarrow}^{*}\\\
z_{\downarrow}&\phantom{-}z_{\uparrow}^{*}\end{array}\right).$ (61)

Inserting the expressions (58) and (59) into Eq. (18) and assuming a stiffness
matrix $\mathcal{J}_{\mu\nu}$ of the form (23), we obtain the $\rm CP^{1}$
action for fluctuating spiral order

$\displaystyle\mathcal{S}_{\text{CP}^{1}}[z,z^{*}]$ $\displaystyle=$
$\displaystyle\int_{\mathcal{T}}dx\,\Big{[}2J^{\perp}_{\mu\nu}(\partial_{\mu}z^{*})(\partial_{\nu}z)$
(62)
$\displaystyle-\,2(J^{\perp}_{\mu\nu}-J^{\Box}_{\mu\nu})j_{\mu}j_{\nu}\Big{]}\,,$

with sum convention for repeated greek indices and the current operator

$j_{\mu}=\frac{i}{2}\left[z^{*}(\partial_{\mu}z)-(\partial_{\mu}z^{*})z\right]\,.$
(63)

For the Néel case, the $\rm CP^{1}$ action is given by the same expression
with $J^{\Box}_{\mu\nu}=0$. We recall that $x=(\tau,\mathbf{r})$ comprises the
imaginary time and space variables, and
$\mathcal{T}=[0,\beta]\times\mathbb{R}^{2}$.

###  IV.2 Large N expansion

The current-current interaction in Eq. (62) can be decoupled by a Hubbard-
Stratonovich transformation, introducing a U(1) gauge field
$\mathcal{A}_{\mu}$, and implementing the constraint (60) by means of a
Lagrange multiplier $\lambda$. The resulting form of the action describes the
so-called massive $\rm CP^{1}$ model [48]

$\displaystyle\mathcal{S}_{\text{CP}^{1}}[z,z^{*},\mathcal{A}_{\mu},\lambda]$
$\displaystyle=$
$\displaystyle\int_{\mathcal{T}}dx\Big{[}2J^{\perp}_{\mu\nu}(D_{\mu}z)^{*}(D_{\nu}z)$
(64) $\displaystyle+$
$\displaystyle\frac{1}{2}M_{\mu\nu}\mathcal{A}_{\mu}\mathcal{A}_{\nu}+i\lambda(z^{*}z-1)\Big{]}\,,\hskip
14.22636pt$

where $D_{\mu}=\partial_{\mu}-i\mathcal{A}_{\mu}$ is the covariant derivative.
The numbers $M_{\mu\nu}$ are the matrix elements of the mass tensor of the
U(1) gauge field,

${\rm M}=4\big{[}1-{\rm J}^{\Box}({\rm J}^{\perp})^{-1}\big{]}^{-1}{\rm
J}^{\Box}\,,$ (65)

where ${\rm J}^{\Box}$ and ${\rm J}^{\perp}$ are the stiffness tensors built
from the matrix elements $J_{\mu\nu}^{\Box}$ and $J_{\mu\nu}^{\perp}$,
respectively.

To perform a large $N$ expansion, we extend the two-component field
$z=(z_{\uparrow},z_{\downarrow})$ to an $N$-component field
$z=(z_{1},\dots,z_{N})$, and rescale it by a factor $\sqrt{N/2}$ so that it
now satisfies the constraint

$z^{*}z=\sum_{\alpha=1}^{N}z^{*}_{\alpha}z_{\alpha}=\frac{N}{2}\,.$ (66)

To obtain a nontrivial limit $N\to\infty$, we rescale the stiffnesses
$J^{\perp}_{\mu\nu}$ and $J^{\Box}_{\mu\nu}$ by a factor $2/N$, yielding the
action

$\displaystyle\mathcal{S}_{\text{CP}^{N-1}}[z,z^{*},\mathcal{A}_{\mu},\lambda]$
$\displaystyle=$
$\displaystyle\int_{\mathcal{T}}dx\Big{[}2J^{\perp}_{\mu\nu}(D_{\mu}z)^{*}(D_{\nu}z)$
(67) $\displaystyle+$
$\displaystyle\\!\frac{N}{4}M_{\mu\nu}\mathcal{A}_{\mu}\mathcal{A}_{\nu}+i\lambda\Big{(}z^{*}z-\frac{N}{2}\Big{)}\Big{]}.\hskip
19.91692pt$

This action describes the massive ${\rm CP}^{N-1}$ model [49], which in $d>2$
dimensions displays two distinct critical points [48, 47, 50]. The first one
belongs to the pure ${\rm CP}^{N-1}$ class, where $M_{\mu\nu}\to 0$
($J^{\Box}_{\mu\nu}=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_{\mu\nu}\to\infty$
($J^{\perp}_{\mu\nu}=J^{\Box}_{\mu\nu}$) 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}=\sqrt{i\langle\lambda\rangle/Z^{\perp}}$. 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)

$\frac{1}{4\pi
J}\int_{0}^{c_{s}{\Lambda_{\mathrm{uv}}}}\\!\frac{\epsilon\,d\epsilon}{\sqrt{\epsilon^{2}+m_{s}^{2}}}\,\mathrm{coth}\left(\frac{\sqrt{\epsilon^{2}+m_{s}^{2}}}{2T}\right)=1\,,$
(68)

where ${\Lambda_{\mathrm{uv}}}$ is an ultraviolet momentum cutoff. The
constant $J$ is an “average” spin stiffness given by

$J=\sqrt{\mathrm{det}\left(\begin{array}[]{cc}J^{\perp}_{xx}&J^{\perp}_{xy}\\\
J^{\perp}_{yx}&J^{\perp}_{yy}\end{array}\right)}\,,$ (69)

and $c_{s}=\sqrt{J/Z^{\perp}}$ is the corresponding average spin wave
velocity. In Sec. IV.3, we shall discuss how to choose the value of
${\Lambda_{\mathrm{uv}}}$. For $m_{s}\ll c_{s}{\Lambda_{\mathrm{uv}}}$, and
$T\ll c_{s}{\Lambda_{\mathrm{uv}}}$, the magnetic correlation length
$\xi_{s}=\frac{1}{2}c_{s}/m_{s}$, behaves as

$\xi_{s}=\frac{c_{s}}{4T\,\sinh^{-1}\\!\left[\frac{1}{2}e^{-\frac{2\pi}{T}(J-J_{c})}\right]}\,,$
(70)

with the critical stiffness

$J_{c}=\frac{c_{s}{\Lambda_{\mathrm{uv}}}}{4\pi}\,.$ (71)

The correlation length is finite at each $T>0$. For $J>J_{c}$, $\xi_{s}$
diverges exponentially for $T\to 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

$n_{0}+\frac{1}{4\pi
J}\int_{0}^{c_{s}{\Lambda_{\mathrm{uv}}}}\\!\frac{\epsilon\,d\epsilon}{\sqrt{\epsilon^{2}+m_{s}^{2}}}=1\,,$
(72)

where $n_{0}=|\langle z_{1}\rangle|^{2}$ is the fraction of condensed bosons.
Eq. (72) can be easily solved, yielding (if $m_{s}\ll{\Lambda_{\mathrm{uv}}}$)

$\displaystyle\begin{cases}&m_{s}=0\\\
&n_{0}=1-\frac{J_{c}}{J}\end{cases}\hskip 14.22636pt\text{for }J>J_{c}\,,$
(73a) $\displaystyle\begin{cases}&n_{0}=0\\\ &m_{s}=2\pi
J\left[\left(J_{c}/J\right)^{2}-1\right]\end{cases}\hskip 14.22636pt\text{for
}J<J_{c}.$ (73b)

The Mermin-Wagner theorem is thus respected already in the saddle-point
approximation to the ${\rm CP}^{N-1}$ representation of the nonlinear sigma
model, that is, there is no long-range order at $T>0$. In the ground state,
long-range 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.

###  IV.3 Choice of ultraviolet cutoff

The impact of spin fluctuations described by the nonlinear sigma model depends
strongly on the ultraviolet cutoff ${\Lambda_{\mathrm{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
${\Lambda_{\mathrm{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 ${\Lambda_{\mathrm{uv}}}$.

A pragmatic choice for the cutoff is given by the ansatz

${\Lambda_{\mathrm{uv}}}=C/\xi_{A}\,,$ (74)

where $C$ is a dimensionless number, and $\xi_{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 $\xi_{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 $\xi_{A}$ can be obtained from the connected spin
amplitude correlation function
$\chi_{A}(\mathbf{r}_{j},\mathbf{r}_{j^{\prime}})=\big{\langle}(\hat{n}_{j}\cdot\vec{S}^{\psi}_{j})(\hat{n}_{j^{\prime}}\cdot\vec{S}^{\psi}_{j^{\prime}})\big{\rangle}_{c}$,
where
$\hat{n}_{j}=\langle\vec{S}^{\psi}_{j}\rangle/|\langle\vec{S}^{\psi}_{j}\rangle|$.
At long distances between $\mathbf{r}_{j}$ and $\mathbf{r}_{j^{\prime}}$ this
function decays exponentially with an exponential dependence $e^{-r/\xi_{A}}$
of the distance $r$. Fourier transforming and using the rotated spin frame
introduced in Sec. III.3, the long distance behavior of
$\chi_{A}(\mathbf{r}_{j},\mathbf{r}_{j^{\prime}})$ can be related to the
momentum dependence of the static correlation function
$\widetilde{\chi}^{ab}(\mathbf{q},0)$ in the amplitude channel $a=b=1$ for
small $\mathbf{q}$, which has the general form

$\widetilde{\chi}^{11}(\mathbf{q},0)\propto\frac{1}{J^{A}_{\alpha\beta}q_{\alpha}q_{\beta}+m_{A}^{2}}\,.$
(75)

The magnetic coherence length is then given by

$\xi_{A}=\sqrt{J_{A}}/(2m_{A})\,,$ (76)

where
$J_{A}=\left(J_{xx}^{A}J_{yy}^{A}-J_{xy}^{A}J_{yx}^{A}\right)^{\frac{1}{2}}$.

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}\approx 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}=\sqrt{2}J_{H}$. Evaluating the RPA spin amplitude correlation function
yields $\xi_{A}=1/\sqrt{8}$ in this limit. With the ansatz (74), one then
obtains $n_{0}=1-4C/\pi$. Matching this with the numerical result
$n_{0}\approx 0.6$ yields $C\approx 0.3$ and ${\Lambda_{\mathrm{uv}}}\approx
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^{\prime}$, respectively, and we fix the ratio of
the hopping amplitudes as $t^{\prime}/t=-0.2$, and we choose a moderate
interaction strength $U=4t$. The energy unit is $t$ in all plots.

###  V.1 Chargon mean-field phase diagram

Figure 2: Pseudocritical temperatures $T^{*}$ and nematic transition
temperature $T_{\rm nem}$ as a function of density $n$. The symbols labeled by
$T^{*}_{m}$ and $T^{*}_{p}$ indicate the temperatures at which the effective
interaction diverges, in the magnetic or in the pairing channel, respectively.
The black solid line labeled by $T^{*}$ indicates the onset of magnetic order
of the chargons and thus the boundary of the pseudogap regime in the absence
of superconductivity. $T^{*}$ coincides with the divergence temperature
$T_{m}^{*}$ for densities where the vertex diverges in the magnetic channel,
and in the hole doped regime it is only slightly lower than $T_{p}^{*}$ when
the leading divergence occurs in the pairing channel. The labels “Néel” and
“Spiral” refer to the type of chargon order. The dashed black line indicates a
topological transition of the quasiparticle Fermi surface within the spiral
regime. The inset shows the irreducible magnetic effective interaction $U_{\rm
eff}^{m}$ as a function of density.

The critical temperatures $T_{m}^{*}$ and $T_{p}^{*}$ at which the vertex
$V^{T}(\mathbf{k}_{1},\mathbf{k}_{2},\mathbf{k}_{3},\mathbf{k}_{4})$ diverges
are shown in Fig. 2. In a density range from $n\approx 0.83$ to $n\approx
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_{\mathrm{eff}}^{m}$ defined in Sec. III.2. The effective interaction
$U_{\mathrm{eff}}^{m}$ 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.

Figure 3: Magnetic gap $\Delta$ (left axis) and incommensurability $\eta$
(right axis) at $T=0$ as functions of the density.

From now on we ignore the pairing instability and focus on the magnetic order
of the chargons. We compute the magnetic order parameter $\Delta$ together
with the optimal wave vector $\mathbf{Q}$ as described in Sec. III.2. In Fig.
3, we show results for $\Delta$ 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\approx 0.73$. On the electron doped side magnetic
order terminates abruptly already at $n\approx 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_{\rm eff}^{m}$ 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 $\mathbf{Q}=(\pi-2\pi\eta,\pi)$, or symmetry related, with an
“incommensurability” $\eta>0$. In Fig. 3 results for $\eta$ 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 $\eta>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_{\rm nem}$ as a function of density. $T_{\rm nem}$
merges with $T^{*}$ at $n\approx 0.88$. For lower densities the magnetic order
is spiral with $\eta>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 single-particle
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.3).

###  V.2 Spinon fluctuations

Figure 4: Out-of-plane (left panel) and in-plane (right panel) spatial ($J$)
and temporal ($Z$) spin stiffnesses in the ground state ($T=0$) as functions
of the filling $n$. In the Néel state (for $n\geq 1$) out-of-plane and in-
plane stiffnesses coincide.

Once the magnetic order parameter $\Delta$ of the chargons and the wave vector
$\mathbf{Q}$ have been computed, we are in the position to calculate the
NL$\sigma$M parameters from the expressions presented in Sec. III.3.

In Fig. 4, we plot results for the spatial and temporal spin stiffnesses
$J^{a}_{\alpha\alpha}$ 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\geq 1$) they coincide. Actually the order parameter defines
an axis, not a plane, in the latter case. All the quantities except $Z^{\Box}$
exhibit pronounced jumps between half-filling and infinitesimal hole-doping.
These discontinuities are due to the appearance of hole pockets around the
points $(\frac{\pi}{2},\frac{\pi}{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}\neq J^{a}_{yy}$. The temporal stiffnesses
$Z^{a}$ exhibit a stronger doping dependence. The peak of $Z^{\perp}$ at
$n\approx 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 $\mathbf{Q}=(\pi-2\pi\eta,\pi)$ and symmetry related.

Figure 5: Magnetic coherence length $\xi_{A}$ (left axis) and average spin
wave velocity $c_{s}$ in the ground state as functions of the filling $n$.

In Fig. 5 we show the magnetic coherence length $\xi_{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
hole-doping, while it increases strongly toward the spiral-to-paramagnet
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 formation of electron pockets upon electron doping. Note that $\xi_{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_{\alpha\alpha}^{\perp}$ and $Z^{\perp}$. 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$\sigma$M.
This quantity depends on the ultraviolet cutoff ${\Lambda_{\mathrm{uv}}}$. As
a reference point, we may use the half-filled Hubbard model at strong
coupling, as discussed in Sec. IV.3, which yields
${\Lambda_{\mathrm{uv}}}\approx 0.9$, and the constant in the ansatz Eq. (74)
is thereby fixed to $C\approx 0.3$.

Figure 6: Fraction of condensed $z$-bosons $n_{0}$ at $T=0$ for two distinct
choices of the ultraviolet cutoff ${\Lambda_{\mathrm{uv}}}$ as a function of
the filling.

In Fig. 6 we show the condensate fraction $n_{0}$ computed with two distinct
choices of the ultraviolet cutoff:
${\Lambda_{\mathrm{uv}}}={\Lambda_{\mathrm{uv}}}(n)=C/\xi_{A}(n)$ and
${\Lambda_{\mathrm{uv}}}=C/\xi_{A}(n=1)$. For the former choice the cutoff
vanishes at the edge of the magnetic region on the hole-doped side, where
$\xi_{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
${\Lambda_{\mathrm{uv}}}(n)=C/\xi_{A}(n)$, is due to the discontinuity of
$\xi_{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.

###  V.3 Electron spectral function

Fractionalizing the electron operators as in Eq. (3), the electron Green’s
function assumes the form

$\displaystyle[\mathcal{G}^{e}_{jj^{\prime}}(\tau)]_{\sigma\sigma^{\prime}}$
$\displaystyle=$ $\displaystyle-\langle
c_{j^{\prime}\sigma^{\prime}}(\tau)c^{*}_{j\sigma}(0)\rangle$ $\displaystyle=$
$\displaystyle-\langle[R_{j^{\prime}}(\tau)]_{\sigma^{\prime}s^{\prime}}[R_{j}^{*}(0)]_{\sigma
s}\,\psi_{j^{\prime}s^{\prime}}(\tau)\psi^{*}_{js}(0)\rangle\,.$

To simplify this expression, we decouple the average $\langle
RR^{*}\psi\psi^{*}\rangle$ as $\langle
RR^{*}\rangle\langle\psi\psi^{*}\rangle$, yielding [14, 9, 10]

$[\mathcal{G}^{e}_{jj^{\prime}}(\tau)]_{\sigma\sigma^{\prime}}=-\langle[R_{j^{\prime}}(\tau)]_{\sigma^{\prime}s^{\prime}}[R_{j}^{*}(0)]_{\sigma
s}\rangle\langle\psi_{j^{\prime}s^{\prime}}(\tau)\psi_{js}^{*}(0)\rangle\,.$
(78)

The spinon Green’s function can be computed from the NL$\sigma$M in the
continuum limit. Using the Schwinger boson parametrization (61), we obtain, in
the large $N$ limit,

$\displaystyle\langle[R(\mathbf{r}_{j^{\prime}},\tau)]_{\sigma^{\prime}s^{\prime}}[R^{*}(\mathbf{r}_{j},0)]_{\sigma
s}\rangle$ $\displaystyle=$
$\displaystyle-D(\mathbf{r}_{j}\\!-\\!\mathbf{r}_{j^{\prime}},\tau)\,\delta_{\sigma\sigma^{\prime}}\delta_{ss^{\prime}}$
(79) $\displaystyle+\,n_{0}\,\delta_{\sigma
s}\delta_{\sigma^{\prime}s^{\prime}}\,.$

The boson propagator $D(\mathbf{r},\tau)$ is the Fourier transform of

$D(\mathbf{q},\omega_{n})=\frac{1}{Z^{\perp}\omega_{n}^{2}+J_{\alpha\beta}^{\perp}q_{\alpha}q_{\beta}+Z^{\perp}m_{s}^{2}}\,,$
(80)

with the bosonic Matsubara frequency $\omega_{n}=2\pi nT$. Fourier
transforming Eq. (78), the electron Green’s function is obtained in momentum
representation as

$\displaystyle\mathcal{G}^{e}(\mathbf{k},\mathbf{k}^{\prime},\nu_{n})$
$\displaystyle=$ $\displaystyle-T\sum_{\omega_{m}}\int_{\mathbf{q}}{\rm
tr}\left[\mathcal{G}(\mathbf{k}-\mathbf{q},\mathbf{k}^{\prime}-\mathbf{q},\nu_{n}-\omega_{m})\right]$
(81) $\displaystyle\times
D(\mathbf{q},\omega_{m})\,\mathbb{1}+n_{0}\,\mathcal{G}(\mathbf{k},\mathbf{k}^{\prime},\nu_{n})\,,$

where $\mathcal{G}(\mathbf{k},\mathbf{k}^{\prime},\nu_{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
$\mathbf{k}=\mathbf{k}^{\prime}$. Furthermore, in this case $\mathcal{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}(\mathbf{k},\nu_{n})$. Performing the Matsubara sum in Eq. (81)
and continuing to real frequencies, we get

$\displaystyle G^{e}(\mathbf{k},\omega)$ $\displaystyle=$
$\displaystyle\sum_{\ell=\pm}\sum_{p=\pm}\int_{|\mathbf{q}|\leq{\Lambda_{\mathrm{uv}}}}\frac{1}{4Z^{\perp}\omega^{\mathrm{sp}}_{\mathbf{q}}}\left(1+\ell\frac{h_{\mathbf{k}-\mathbf{q}}}{e_{\mathbf{k}-\mathbf{q}}}\right)$
(82) $\displaystyle\times$
$\displaystyle\frac{f(pE^{\ell}_{\mathbf{k}-\mathbf{q}})+n_{B}(\omega^{\mathrm{sp}}_{\mathbf{q}})}{\omega+i0^{+}-E^{\ell}_{\mathbf{k}-\mathbf{q}}+p\,\omega^{\mathrm{sp}}_{\mathbf{q}}}+\\{\mathbf{k}\to-\mathbf{k}\\}\,,\hskip
19.91692pt$

where

$\omega^{\mathrm{sp}}_{\mathbf{q}}=\sqrt{(J^{\perp}_{\alpha\beta}q_{\alpha}q_{\beta})/Z^{\perp}+m_{s}^{2}}\,,$
(83)

and $n_{B}(x)=(e^{x/T}-1)^{-1}$ is the Bose distribution function.

Figure 7: Quasiparticle Fermi surfaces defined as zeros of the chargon
quasiparticle energies $E_{\mathbf{k}}^{\pm}$ (left column) and momentum
dependence of electron spectral function at zero frequency (right column) for
various electron densities. The temperature is $T=0.05t$.

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 single-particle 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 [54, 55, 52, 56, 57, 58, 32].

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 [15, 7,
8, 9, 10, 11]. 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_{\rm 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_{\rm 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.

## Acknowledgements

We are very grateful to Andres Greco, Elio König and Demetrio Vilardi for
valuable discussions.

##  Appendix A Linear term in the gauge field

In this Appendix we show that the linear term in Eq. (12) vanishes. Fourier
transforming the vertex and the expectation value, the coefficient
$\mathcal{B}_{\mu}^{a}$ can be written as

$\mathcal{B}_{\mu}^{a}=\frac{1}{2}\int_{\mathbf{k}}T\sum_{\nu_{n}}\gamma_{\mu}^{(1)}(\mathbf{k}){\rm
Tr}\left[\sigma^{a}\mathcal{G}(\mathbf{k},\mathbf{k},\nu_{n})\right]\,.$ (84)

Inserting $\mathcal{G}$ from Eq. (52) one immediately sees that
$\mathcal{B}_{\mu}^{1}=\mathcal{B}_{\mu}^{2}=0$ for $\mu=0,1,2$, and
$\mathcal{B}_{0}^{3}=0$, too. Performing the Matsubara sum for
$\mathcal{B}^{3}_{\alpha}$ with $\alpha=1,2$, we obtain

$\mathcal{B}^{3}_{\alpha}=\frac{1}{2}\int_{\mathbf{k}}\sum_{\ell=\pm}\left[(\partial_{k_{\alpha}}\epsilon_{\mathbf{k}})u^{\ell}_{\mathbf{k}}f(E^{\ell}_{\mathbf{k}})+(\partial_{k_{\alpha}}\epsilon_{\mathbf{k}+\mathbf{Q}})u^{-\ell}_{\mathbf{k}}f(E^{\ell}_{\mathbf{k}})\right]\,,$
(85)

where $u^{\ell}_{\mathbf{k}}=\frac{1}{2}\big{(}1+\ell
h_{\mathbf{k}}/\sqrt{h_{\mathbf{k}}^{2}+\Delta^{2}}\big{)}$ with
$h_{\mathbf{k}}=\frac{1}{2}(\epsilon_{\mathbf{k}}-\epsilon_{\mathbf{k}+\mathbf{Q}})$.
One can see by direct calculation that this term vanishes if $\partial
F(\mathbf{Q})/\partial\mathbf{Q}$ with $F(\mathbf{Q})$ given by Eq. (30)
vanishes. Hence, $\mathcal{B}^{3}_{\alpha}$ vanishes if $\mathbf{Q}$ minimizes
the free energy. A similar result has been obtained in Ref. [26].

##  Appendix B Derivation of the NL$\sigma$M

Here we derive the NL$\sigma$M action (18) from Eq. (12). We first prove the
identity

$\partial_{\mu}\mathcal{R}=-i\mathcal{R}\,\Sigma^{a}A_{\mu}^{a},$ (86)

where $\mathcal{R}$ is defined by Eq. (17), and $\Sigma^{a}$ are the
generators of the SU(2) in the adjoint representation,

$\Sigma^{a}_{bc}=-i\varepsilon^{abc},$ (87)

with $\varepsilon^{abc}$ the Levi-Civita tensor. Rewriting Eq. (17) as

$\mathcal{R}^{ab}=\frac{1}{2}{\rm
Tr}\left[R^{\dagger}\sigma^{a}R^{\phantom{\dagger}}\sigma^{b}\right]\,,$ (88)

we obtain the derivative of $\mathcal{R}$ in the form,

$\displaystyle\partial_{\mu}\mathcal{R}^{ab}$ $\displaystyle=$
$\displaystyle{\rm
Tr}\left[R^{\dagger}\sigma^{a}\,(\partial_{\mu}R)\sigma^{b}\right]={\rm
Tr}\left[R^{\dagger}\sigma^{a}RR^{\dagger}(\partial_{\mu}R)\sigma^{b}\right]$
(89) $\displaystyle=$
$\displaystyle-i\mathcal{R}^{ac}\Sigma^{d}_{cb}A_{\mu}^{d},$

which is the identity in (86).

We now aim to express the object
$\frac{1}{2}\mathcal{J}^{ab}_{\mu\nu}A_{\mu}^{a}A^{b}_{\nu}$ in terms of the
matrix field $\mathcal{R}$. We write the stiffness matrix in terms of a new
matrix $\mathcal{P}_{\mu\nu}$ via

$\mathcal{J}^{ab}_{\mu\nu}={\rm
Tr}[\mathcal{P}_{\mu\nu}]\delta_{ab}-\mathcal{P}^{ab}_{\mu\nu}={\rm
Tr}\left[\mathcal{P}_{\mu\nu}\Sigma^{a}\Sigma^{b}\right].$ (90)

Using $\mathcal{R}^{T}\mathcal{R}=\mathbb{1}$, we obtain

$\displaystyle\frac{1}{2}\mathcal{J}^{ab}_{\mu\nu}A_{\mu}^{a}A^{b}_{\nu}$
$\displaystyle=$ $\displaystyle\frac{1}{2}{\rm
Tr}\left[\mathcal{P}_{\mu\nu}\,\Sigma^{a}\,\mathcal{R}^{T}\mathcal{R}\,\Sigma^{b}\right]A_{\mu}^{a}A^{b}_{\nu}$
(91) $\displaystyle=$ $\displaystyle\frac{1}{2}{\rm
Tr}\left[\mathcal{P}_{\mu\nu}(\partial_{\mu}\mathcal{R}^{T})(\partial_{\nu}\mathcal{R})\right],$

where we have used Eq. (86) in the last line. The above equation yields Eq.
(18). Relation (90) can be easily inverted using ${\rm
Tr}[\mathcal{J}_{\mu\nu}]=2{\rm Tr}[\mathcal{P}_{\mu\nu}]$.

##  Appendix C Details on the large-$N$ expansion

In this Appendix, we describe some details regarding the saddle point
equations of the CPN-1 action. Integrating out the $z$-bosons from Eq. (67),
we obtain the effective action [46]

$\displaystyle\mathcal{S}[\mathcal{A}_{\mu},\lambda]$ $\displaystyle=$
$\displaystyle
N\int_{\mathcal{T}}dx\Big{[}\ln\left(-2J^{\perp}_{\mu\nu}D_{\mu}D_{\nu}+i\lambda\right)-\frac{i}{2}\lambda$
(92)
$\displaystyle+\,\frac{1}{4}M_{\mu\nu}\mathcal{A}_{\mu}\mathcal{A}_{\nu}\Big{]}\,.$

In the large $N$ limit the functional integral for its partition function is
dominated by its saddle point, which is determined by the stationarity
equations

$\frac{\delta\mathcal{S}}{\delta\mathcal{A}_{\mu}}=\frac{\delta\mathcal{S}}{\delta\lambda}=0\,.$
(93)

The first condition implies $\mathcal{A}_{\mu}=0$, that is, in the large-$N$
limit the U(1) gauge field fluctuations are totally suppressed. The variation
with respect to $\lambda$ gives, assuming a spatially uniform average value
for $\lambda$,

$n_{0}+T\sum_{\omega_{n}}\int_{\mathbf{q}}\frac{1}{Z^{\perp}\omega_{n}^{2}+J^{\perp}_{\alpha\beta}q_{\alpha}q_{\beta}+i\langle\lambda\rangle}=1\,,$
(94)

where $n_{0}$ is the fraction of condensed bosons, which can be nonzero at
$T=0$. Performing the sum over the bosonic Matsubara frequencies
$\omega_{n}=2n\pi T$, inserting the identity

$1=\int_{0}^{\infty}\\!d\epsilon\,\delta\Big{(}\epsilon-\sqrt{J^{\perp}_{\alpha\beta}q_{\alpha}q_{\beta}/Z^{\perp}}\,\Big{)},$
(95)

and performing the $\mathbf{q}$-integral, we obtain Eq. (68) at $T>0$ and Eq.
(72) at $T=0$.

## References

  * Keimer  _et al._ [2015] B. Keimer, S. A. Kivelson, M. R. Norman, S. Uchida, and J. Zaanen, From quantum matter to high-temperature superconductivity in copper oxides, Nature 518, 179 (2015).
  * Proust and Taillefer [2019] C. Proust and L. Taillefer, The Remarkable Underlying Ground States of Cuprate Superconductors, Annu. Rev. Condens. Matter Phys. 10, 409 (2019).
  * Qin  _et al._ [2022] M. Qin, T. Schäfer, S. Andergassen, P. Corboz, and E. Gull, The Hubbard Model: A Computational Perspective, Annu. Rev. Condens. Matter Phys. 13, 275 (2022).
  * Gunnarsson  _et al._ [2015] O. Gunnarsson, T. Schäfer, J. P. F. LeBlanc, E. Gull, J. Merino, G. Sangiovanni, G. Rohringer, and A. Toschi, Fluctuation Diagnostics of the Electron Self-Energy: Origin of the Pseudogap Physics, Phys. Rev. Lett. 114, 236402 (2015).
  * Moriya and Ueda [2000] T. Moriya and K. Ueda, Spin fluctuations and high temperature superconductivity, Adv. Phys. 49, 555 (2000).
  * Vilk and Tremblay [1996] Y. M. Vilk and A.-M. S. Tremblay, Destruction of Fermi-liquid quasiparticles in two dimensions by critical fluctuations, Europhys. Lett. 33, 159 (1996).
  * Sachdev and Chowdhury [2016] S. Sachdev and D. Chowdhury, The novel metallic states of the cuprates: Topological Fermi liquids and strange metals, Prog. Theor. Exp. Phys. 2016, 12C102 (2016).
  * Chatterjee  _et al._ [2017a] S. Chatterjee, S. Sachdev, and M. S. Scheurer, Intertwining Topological Order and Broken Symmetry in a Theory of Fluctuating Spin-Density Waves, Phys. Rev. Lett. 119, 227002 (2017a).
  * Scheurer  _et al._ [2018] M. S. Scheurer, S. Chatterjee, W. Wu, M. Ferrero, A. Georges, and S. Sachdev, Topological order in the pseudogap metal, Proc. Natl. Acad. Sci. USA 115, E3665 (2018).
  * Wu  _et al._ [2018] W. Wu, M. S. Scheurer, S. Chatterjee, S. Sachdev, A. Georges, and M. Ferrero, Pseudogap and Fermi-Surface Topology in the Two-Dimensional Hubbard Model, Phys. Rev. X 8, 021048 (2018).
  * Sachdev [2019] S. Sachdev, Topological order, emergent gauge fields, and Fermi surface reconstruction, Rep. Prog. Phys. 82, 014001 (2019).
  * Schulz [1995] H. Schulz, Functional Integrals for Correlated Electrons, in  _The Hubbard Model_, edited by D. Baeriswyl (Plenum, New York, 1995).
  * Dupuis [2002] N. Dupuis, Spin fluctuations and pseudogap in the two-dimensional half-filled Hubbard model at weak coupling, Phys. Rev. B 65, 245118 (2002).
  * Borejsza and Dupuis [2004] K. Borejsza and N. Dupuis, Antiferromagnetism and single-particle properties in the two-dimensional half-filled Hubbard model: A nonlinear sigma model approach, Phys. Rev. B 69, 085119 (2004).
  * Sachdev  _et al._ [2009] S. Sachdev, M. A. Metlitski, Y. Qi, and C. Xu, Fluctuating spin density waves in metals, Phys. Rev. B 80, 155129 (2009).
  * Wang  _et al._ [2014] J. Wang, A. Eberlein, and W. Metzner, Competing order in correlated electron systems made simple: Consistent fusion of functional renormalization and mean-field theory, Phys. Rev. B 89, 121116(R) (2014).
  * Bonetti [2022] P. M. Bonetti, Local Ward identities for collective excitations in fermionic systems with spontaneously broken symmetries (2022), arXiv:2204.04132 .
  * Weng  _et al._ [1991] Z. Y. Weng, C. S. Ting, and T. K. Lee, Path-integral approach to the Hubbard model, Phys. Rev. B 43, 3790 (1991).
  * Metzner  _et al._ [2012] W. Metzner, M. Salmhofer, C. Honerkamp, V. Meden, and K. Schönhammer, Functional renormalization group approach to correlated fermion systems, Rev. Mod. Phys. 84, 299 (2012).
  * Haldane [1983a] F. D. M. Haldane, Nonlinear Field Theory of Large-Spin Heisenberg Antiferromagnets: Semiclassically Quantized Solitons of the One-Dimensional Easy-Axis Néel State, Phys. Rev. Lett. 50, 1153 (1983a).
  * Haldane [1983b] F. Haldane, Continuum dynamics of the 1-D Heisenberg antiferromagnet: Identification with the O(3) nonlinear sigma model, Physics Letters A 93, 464 (1983b).
  * Azaria  _et al._ [1990] P. Azaria, B. Delamotte, and T. Jolicoeur, Nonuniversality in helical and canted-spin systems, Phys. Rev. Lett. 64, 3175 (1990).
  * Azaria  _et al._ [1992] P. Azaria, B. Delamotte, and D. Mouhanna, Low-temperature properties of two-dimensional frustrated quantum antiferromagnets, Phys. Rev. Lett. 68, 1762 (1992).
  * Azaria  _et al._ [1993a] P. Azaria, B. Delamotte, and D. Mouhanna, Spontaneous symmetry breaking in quantum frustrated antiferromagnets, Phys. Rev. Lett. 70, 2483 (1993a).
  * Azaria  _et al._ [1993b] P. Azaria, B. Delamotte, F. Delduc, and T. Jolicoeur, A renormalization-group study of helimagnets in $D=2+\varepsilon$ dimensions, Nuclear Physics B 408, 485 (1993b).
  * Klee and Muramatsu [1996] S. Klee and A. Muramatsu, SO(3) nonlinear $\sigma$ model for a doped quantum helimagnet, Nucl. Phys. B 473, 539 (1996).
  * Igoshev  _et al._ [2010] P. A. Igoshev, M. A. Timirgazin, A. A. Katanin, A. K. Arzhnikov, and V. Y. Irkhin, Incommensurate magnetic order and phase separation in the two-dimensional Hubbard model with nearest- and next-nearest-neighbor hopping, Phys. Rev. B 81, 094407 (2010).
  * Frésard  _et al._ [1991] R. Frésard, M. Dzierzawa, and P. Wölfle, Slave-Boson Approach to Spiral Magnetic Order in the Hubbard Model, Europhys. Lett. 15, 325 (1991).
  * Chubukov and Musaelian [1995] A. V. Chubukov and K. A. Musaelian, Magnetic phases of the two-dimensional Hubbard model at low doping, Phys. Rev. B 51, 12605 (1995).
  * Yamase  _et al._ [2016] H. Yamase, A. Eberlein, and W. Metzner, Coexistence of Incommensurate Magnetism and Superconductivity in the Two-Dimensional Hubbard Model, Phys. Rev. Lett. 116, 096402 (2016).
  * Vilardi  _et al._ [2018] D. Vilardi, C. Taranto, and W. Metzner, Dynamically enhanced magnetic incommensurability: Effects of local dynamics on nonlocal spin correlations in a strongly correlated metal, Phys. Rev. B 97, 235110 (2018).
  * Bonetti  _et al._ [2020] P. M. Bonetti, J. Mitscherling, D. Vilardi, and W. Metzner, Charge carrier drop at the onset of pseudogap behavior in the two-dimensional Hubbard model, Phys. Rev. B 101, 165142 (2020).
  * Berges  _et al._ [2002] J. Berges, N. Tetradis, and C. Wetterich, Non-perturbative renormalization flow in quantum field theory and statistical physics, Physics Reports 363, 223 (2002).
  * Dupuis  _et al._ [2021] N. Dupuis, L. Canet, A. Eichhorn, W. Metzner, J. Pawlowski, M. Tissier, and N. Wschebor, The nonperturbative functional renormalization group and its applications, Physics Reports 910, 1–114 (2021).
  * Bonetti [2020] P. M. Bonetti, Accessing the ordered phase of correlated Fermi systems: Vertex bosonization and mean-field theory within the functional renormalization group, Phys. Rev. B 102, 235160 (2020).
  * Vilardi  _et al._ [2020] D. Vilardi, P. M. Bonetti, and W. Metzner, Dynamical functional renormalization group computation of order parameters and critical temperatures in the two-dimensional Hubbard model, Phys. Rev. B 102, 245128 (2020).
  * Honerkamp and Salmhofer [2001] C. Honerkamp and M. Salmhofer, Temperature-flow renormalization group and the competition between superconductivity and ferromagnetism, Phys. Rev. B 64, 184516 (2001).
  * Husemann and Salmhofer [2009] C. Husemann and M. Salmhofer, Efficient parametrization of the vertex function, $\Omega$ scheme, and the $t,{t}^{{}^{\prime}}$ Hubbard model at van Hove filling, Phys. Rev. B 79, 195125 (2009).
  * Husemann  _et al._ [2012] C. Husemann, K.-U. Giering, and M. Salmhofer, Frequency-dependent vertex functions of the ($t,{t}^{{}^{\prime}}$) Hubbard model at weak coupling, Phys. Rev. B 85, 075121 (2012).
  * Vilardi  _et al._ [2017] D. Vilardi, C. Taranto, and W. Metzner, Nonseparable frequency dependence of the two-particle vertex in interacting fermion systems, Phys. Rev. B 96, 235110 (2017).
  * Vilardi  _et al._ [2019] D. Vilardi, C. Taranto, and W. Metzner, Antiferromagnetic and $d$-wave pairing correlations in the strongly interacting two-dimensional Hubbard model from the functional renormalization group, Phys. Rev. B 99, 104501 (2019).
  * Lichtenstein  _et al._ [2017] J. Lichtenstein, D. Sánchez de la Pẽna, D. Rohe, E. Di Napoli, C. Honerkamp, and S. A. Maier, High-performance functional renormalization group calculations for interacting fermions, Comput. Phys. Commun. 213, 100 (2017).
  * Bonetti and Metzner [2022] P. M. Bonetti and W. Metzner, Spin stiffness, spectral weight, and Landau damping of magnons in metallic spiral magnets, Phys. Rev. B 105, 134426 (2022).
  * Kampf [1996] A. P. Kampf, Collective excitations in itinerant spiral magnets, Phys. Rev. B 53, 747 (1996).
  * Sachdev  _et al._ [1995] S. Sachdev, A. V. Chubukov, and A. Sokol, Crossover and scaling in a nearly antiferromagnetic Fermi liquid in two dimensions, Phys. Rev. B 51, 14874 (1995).
  * Auerbach [1994] A. Auerbach,  _Interacting Electrons and Quantum Magnetism_ (Springer-Verlag, New York, 1994).
  * Chubukov  _et al._ [1994a] A. V. Chubukov, S. Sachdev, and T. Senthil, Quantum phase transitions in frustrated quantum antiferromagnets, Nuclear Physics B 426, 601 (1994a).
  * Azaria  _et al._ [1995] P. Azaria, P. Lecheminant, and D. Mouhanna, The massive CPN-1 model for frustrated spin systems, Nuclear Physics B 455, 648 (1995).
  * Campostrini and Rossi [1993] M. Campostrini and P. Rossi, The $1/N$ expansion of two-dimensional spin models, Riv. Nuovo Cim. 16, 1 (1993).
  * Chubukov  _et al._ [1994b] A. V. Chubukov, T. Senthil, and S. Sachdev, Universal magnetic properties of frustrated quantum antiferromagnets in two dimensions, Phys. Rev. Lett. 72, 2089 (1994b).
  * Manousakis [1991] E. Manousakis, The spin-$\frac{1}{2}$ Heisenberg antiferromagnet on a square lattice and its application to the cuprous oxides, Rev. Mod. Phys. 63, 1 (1991).
  * Eberlein  _et al._ [2016] A. Eberlein, W. Metzner, S. Sachdev, and H. Yamase, Fermi Surface Reconstruction and Drop in the Hall Number due to Spiral Antiferromagnetism in High-$T_{c}$ Cuprates, Phys. Rev. Lett. 117, 187001 (2016).
  * Mitscherling and Metzner [2021] J. Mitscherling and W. Metzner, Non-Hermitian band topology from momentum-dependent relaxation in two-dimensional metals with spiral magnetism, Phys. Rev. B 104, L201107 (2021).
  * Storey [2016] J. G. Storey, Hall effect and Fermi surface reconstruction via electron pockets in the high $T_{c}$ cuprates, Europhys. Lett. 113, 2700 (2016).
  * Storey [2017] J. G. Storey, Simultaneous drop in mean free path and carrier density at the pseudogap onset in high-$T_{c}$ cuprates, Supercond. Sci. Technol. 30, 104008 (2017).
  * Chatterjee  _et al._ [2017b] S. Chatterjee, S. Sachdev, and A. Eberlein, Thermal and electrical transport in metals and superconductors across antiferromagnetic and topological quantum transitions, Phys. Rev. B 96, 075103 (2017b).
  * Verret  _et al._ [2017] S. Verret, O. Simard, M. Charlebois, D. Sénéchal, and A.-M. S. Tremblay, Phenomenological theories of the low-temperature pseudogap: Hall number, specific heat, and Seebeck coefficient, Phys. Rev. B 96, 125139 (2017).
  * Mitscherling and Metzner [2018] J. Mitscherling and W. Metzner, Longitudinal conductivity and Hall coefficient in two-dimensional metals with spiral magnetic order, Phys. Rev. B 98, 195126 (2018).
  * Grissonnanche  _et al._ [2022] G. Grissonnanche, O. Cyr-Choinière, J. Day, R. Liang, D. A. Bonn, W. N. Hardy, N. Doiron-Leyraud, and L. Taillefer, No nematicity at the onset temperature of the pseudogap phase in the cuprate superconductor $\rm YBa_{2}Cu_{3}O_{y}$ (2022), arXiv:2205.05233 .
  * Zheng  _et al._ [2017] B.-X. Zheng, C.-M. Chung, P. Corboz, G. Ehlers, M.-P. Qin, R. M. Noack, H. Shi, W. S. R., S. Zhang, and G. K.-L. Chan, Stripe order in the underdoped region of the two-dimensional Hubbard model, Science 358, 1155 (2017).
  * Sachdev  _et al._ [2019] S. Sachdev, H. D. Scammell, M. S. Scheurer, and G. Tarnopolsky, Gauge theory for the cuprates near optimal doping, Phys. Rev. B 99, 054516 (2019).
  * Charlebois  _et al._ [2017] M. Charlebois, S. Verret, A. Foley, O. Simard, D. Sénéchal, and A.-M. S. Tremblay, Hall effect in cuprates with an incommensurate collinear spin-density wave, Phys. Rev. B 96, 205132 (2017).