# Magnetic properties of the spiral spin liquid and surrounding phases in the
square lattice XY model

Matías G. Gonzalez  Helmholtz-Zentrum Berlin für Materialien und Energie,
Hahn-Meitner-Platz 1, 14109 Berlin, Germany  Dahlem Center for Complex Quantum
Systems and Fachbereich Physik, Freie Universität Berlin, 14195 Berlin,
Germany  Anna Fancelli  Helmholtz-Zentrum Berlin für Materialien und Energie,
Hahn-Meitner-Platz 1, 14109 Berlin, Germany  Dahlem Center for Complex Quantum
Systems and Fachbereich Physik, Freie Universität Berlin, 14195 Berlin,
Germany  Han Yan  Institute for Solid State Physics, The University of Tokyo,
Kashiwa, Chiba 277-8581, Japan  Department of Physics and Astronomy, Rice
University, Houston, TX 77005, USA  Smalley-Curl Institute, Rice University,
Houston, TX 77005, USA  Johannes Reuther  Helmholtz-Zentrum Berlin für
Materialien und Energie, Hahn-Meitner-Platz 1, 14109 Berlin, Germany  Dahlem
Center for Complex Quantum Systems and Fachbereich Physik, Freie Universität
Berlin, 14195 Berlin, Germany

###### Abstract

Spiral spin liquids possess a subextensively degenerate ground-state manifold,
represented by a continuum of energy minima in reciprocal space. Since a small
change of the spiral state wavevector requires a global change of the spin
configuration in real space, it is a priori unclear how such systems can
fluctuate within the degenerate ground state manifold. Only recently it was
proposed that momentum vortices are responsible for the liquidity of the
spiral phase and that these systems are closely related to an emergent rank-2
U(1) gauge theory [H. Yan and J. Reuther, Phys. Rev. Research 4, 023175
(2022)]. As a consequence of this gauge structure, four-fold pinch-point
singularities were found in a generalized spin correlator. In this article, we
use classical Monte Carlo and molecular dynamics calculations to embed the
previously studied spiral spin liquid into a broader phase diagram of the
square lattice XY model. We find a multitude of unusual phases and phase
transitions surrounding the spiral spin liquid such as an effective four-state
Potts transition into a colinear double-striped phase resulting from the
spontaneous breaking of two coupled $\mathds{Z}_{2}$ symmetries. Since this
phase is stabilized by entropic effects selecting the momenta away from the
spiral manifold, it undergoes a re-entrance phenomenon at low temperatures
into a nematic spiral phase. We also observe a region of parameters where the
phase transition into the spiral spin liquid does not break any symmetries and
where the critical exponents do not match those of standard universality
classes. We study the importance of momentum vortices in driving this phase
transition and discuss the possibility of a Kosterlitz-Thouless transition of
momentum vortices. Finally, we explore the regime where the rank-2 U(1) gauge
theory is valid by investigating the four-fold pinch point singularities
across the phase diagram.

##  I Introduction

Spin liquids (SLs) are one of the most sought-after states in the field of
modern condensed matter physics [1, 2, 3]. This is because these disordered
yet highly correlated spin states offer an ideal platform to study a wide
range of captivating emergent phenomena [4, 5, 6, 7, 8]. Examples include
hidden topological order, long-range entanglement, fractional excitations,
fracton physics, and many more. In the classical limit, where spins behave as
unit vectors, much attention has been focused on lattices of corner-sharing
frustrated motifs such as triangles or tetrahedra [9, 10, 11, 12, 13]. These
lattices allow for a re-writing of a nearest neighbor Heisenberg Hamiltonian
into a sum of complete squares over the triangular or tetrahedral building
blocks. Hence, all ground states are determined by local spin constraints such
as the condition of a vanishing total spin on each of these motifs [14, 15,
16, 17]. Prominent examples occur on kagome and pyrochlore lattices, where
extensively degenerate ground state manifolds are connected by local spin
flips and the spin structure factors exhibit characteristic two-fold pinch-
point singularities. Upon adding quantum fluctuations these highly degenerate
ground states comprise a well-known route toward inducing quantum spin liquid
behavior [18, 19, 20, 21, 22]. Particularly fascinating is the possibility of
realizing emergent gauge theories via these local spin constraints, which in
pyrochlore quantum spin ice systems gives rise to an effective U(1) quantum
electrodynamics [23, 24].

Another kind of spin liquid is the spiral spin liquid induced by competition
between the nearest and further neighbor couplings. Contrary to the
aforementioned spin liquids which result from local spin constraints, this
kind of system exhibits only a subextensively degenerate classical ground-
state manifold given by a continuum of minima in reciprocal space, homotopic
to a ring in two dimensions or sphere (or other 2D surfaces) in three
dimensions [25, 26, 27, 28, 29, 30, 31, 32, 33, 34]. In these systems,
however, a small change in the spiral wave vector implies a global,
coordinated change of all spins in real space, contrary to the local spin
flips of kagome or pyrochlore magnets. Thus, it is not a priori clear how
spiral spin liquids can explore the whole ground-state manifold. A solution to
this issue was proposed very recently, identifying the momentum vortices as
the effective local degrees of freedom responsible for the liquidity of the
phase [35]. The proliferation of these topological defects allows each small
patch of spins in the system to visit different wave vectors along the spiral
surface at very little energy cost.

Furthermore, a low-energy effective theory for a classical spiral spin liquid
was derived, which possesses remarkable similarities with a rank-2 U(1) gauge
theory. The family of such theories constitutes higher-rank generalizations of
U(1) electromagnetism known for their unusual kinetically constrained fracton
quasiparticles [36, 37, 35]. As a result of the tensor structure of the
associated Gauss’ law, a rank-2 U(1) gauge theory exhibits four-fold pinch-
point singularities in the electric-field correlator [38] instead of the usual
two-fold pinch points of the electric-field correlator in standard U(1)
electromagnetism. The effective theory also has a close relation to the
fracton theory of smectic matter [39, 40, 41, 42, 43]. In particular, the
topological defects known as momentum vortices have the same exotic feature as
those in smectic matter, namely the winding number cannot be larger than 1,
but can be any arbitrary negative number. The mathematics behind this has been
discussed in detail in Refs. [35, 43].

In the spiral system of Ref. [35] – a classical square lattice XY model with
up to third nearest neighbor interactions – the emergent electric field
translates into a second-order derivative of the spin angle field and the
corresponding correlation function of these objects has indeed been
numerically confirmed to show four-fold pinch points. Nevertheless, the range
of validity of the rank-2 U(1) gauge theory emerging from a spiral spin liquid
remains uncharted since the exact mapping requires the radius of the spiral
contour in momentum space to be nearly vanishing. A further question left
unanswered in Ref. [35] is the nature of the thermal phase transition observed
as a sharp heat capacity peak when the system enters the spiral spin liquid
regime.

Motivated by the previous results but also by remaining open questions, in
this article, we revisit the spiral phase of the classical XY model on the
square lattice using classical Monte Carlo and molecular dynamics
calculations. Compared to Ref. [35] where the system was only studied for a
single parameter set for a fixed, small, and almost circular spiral ring,
here, we explore a larger region of the phase diagram with varying radii of
the spiral contour, revealing the broader interplay of phases in which the
spiral spin liquid is embedded. Particularly, moving away from the circular
spiral ring we study the range of validity of the correspondence between the
spiral spin liquid and the rank-2 U(1) gauge theory. We also uncover a rich
phase diagram in temperature and parameter space that contains a variety of
interesting magnetic phases beyond the previously studied spiral spin liquid,
pancake liquids, and rigid vortex networks. This includes double stripe states
stabilized by two intertwined entropic selection mechanisms that lead to a
$\mathds{Z}_{2}\times\mathds{Z}_{2}$ symmetry breaking of lattice translation
and rotation symmetries. However, since these selection mechanisms have a
finite energy cost, at lower temperatures the system regains one of the broken
$\mathds{Z}_{2}$ symmetries, showing a re-entrant behavior into a nematic
spiral phase. The different broken symmetries lead to phase transitions
belonging to the 2D Ising and four-state Potts universality classes.
Furthermore, between the pancake and spiral spin liquids, we observe a phase
transition with a logarithmically divergent specific heat that is not related
to any symmetry breaking. We investigate the behavior of momentum vortices at
this transition and discuss the possibility of a Kosterlitz-Thouless
transition of momentum vortices.

The remainder of the article is organized as follows: Section II introduces
the XY model on the square lattice and the key properties of the spiral
contour. In Sec. III we show the finite-temperature phase diagram as obtained
in our calculations, with different subsections devoted to the individual
phases. In Sec. IV we investigate the dynamical spin structure factor and
discuss the characteristic spectral features of the different phases. In Sec.
V we discuss the roles of spin and momentum vortices, especially at the phase
transition into the spiral spin liquid regime. Finally, in Sec. VI we study
the validity of the rank-2 U(1) gauge theory throughout the phase diagram.
Sec. VII contains the conclusions of our work.

##  II spiral spin model

Figure 1: Square lattice with up to third-nearest-neighbor exchange
interactions shown in orange, purple, and blue, respectively, emanating from
an arbitrary site (black).

The Hamiltonian for the XY model on the square lattice considered in this
article is given by

$\mathcal{H}=J_{1}\sum_{\langle
ij\rangle_{1}}\mathbf{S}_{i}\cdot\mathbf{S}_{j}+J_{2}\sum_{\langle
ij\rangle_{2}}\mathbf{S}_{i}\cdot\mathbf{S}_{j}+J_{3}\sum_{\langle
ij\rangle_{3}}\mathbf{S}_{i}\cdot\mathbf{S}_{j},$ (1)

where $\langle ij\rangle_{1}$, $\langle ij\rangle_{2}$, and $\langle
ij\rangle_{3}$ indicate first-, second-, and third-nearest neighbor sites,
respectively (see Fig. 1). The XY spins $\mathbf{S}_{i}=(S_{ix},S_{iy})$ are
normalized as $|\mathbf{S}_{i}|=1$ and represented by an angle
$\phi\in\left[0,2\pi\right)$ via $\mathbf{S}_{i}=(\cos\phi_{i},\sin\phi_{i})$.
Throughout this work we constrain the coupling parameters as follows: We fix
the nearest-neighbor interaction to be ferromagnetic, $J_{1}=-1$, and consider
antiferromagnetic $J_{2}$ and $J_{3}$. Furthermore, we fix $J_{2}>0.25$ and
$J_{3}=J_{2}/2$, as this gives rise to the spiral ground states of interest
here [35]. It is thus convenient to use the parametrization
$\delta=J_{2}-0.25$, such that the transition between a pure ferromagnet and
spin spirals occurs at $\delta=0$. The degenerate momenta
$\mathbf{q}=(q_{x},q_{y})$ of the spin spirals lie on a ring in reciprocal
space, determined by the solutions of the equation

$\cos(q_{x})+\cos(q_{y})=\frac{2}{4\delta+1}.$ (2)

In Fig. 2 we show the continuous manifolds of ground state wave-vectors in
reciprocal space for varying values of $\delta$ from $\delta=0$ to $0.25$.
When $\delta=0.25$ the spiral ring reaches the special wave-vectors
$(\pm\frac{\pi}{2},0)$ and $(0,\pm\frac{\pi}{2})$.

From Eq. (2) it is clear that $\delta=0$ only allows for the ferromagnetic
solution $\mathbf{q}=(q_{x},q_{y})=(0,0)$. When $\delta$ is finite but
$\delta\ll 1$, solutions in the vicinity of $\mathbf{q}=(0,0)$ can be found by
approximating $\cos(q)=1-\frac{q^{2}}{2}+\mathcal{O}(q^{4})$. Inserting this
into Eq. (2) results in circular rings given by

$\left(\frac{q_{x}}{4\sqrt{\delta}}\right)^{2}+\left(\frac{q_{y}}{4\sqrt{\delta}}\right)^{2}=1,$
(3)

where the radius is $4\sqrt{\delta}$ [35]. This is evidenced in Fig. 2, where
contour lines with small $\delta$ close to the origin become more circular. It
is in this limit that the spiral spin liquid physics is governed by momentum
vortices whose effective theory corresponds to a rank-2 U(1) gauge theory
[35].

Figure 2: Continuous ground-state manifolds of wave-vectors as a function of
$\delta$ corresponding to the spiral solutions of Eq. (2). Contour lines are
shown for $\delta=0.001$, 0.01, 0.03, 0.05, 0.10, 0.15, and 0.20 from smallest
to largest.

##  III Phase Diagram

We perform classical Monte Carlo (cMC) calculations in a range of the phase
diagram given by $\delta\in[0,0.25]$, from the ferromagnetic limit to the
point where the spiral ring touches the wave vectors $(\pm\frac{\pi}{2},0)$
and $(0,\pm\frac{\pi}{2})$. Note that the last wave vectors correspond to
spirals with a periodicity of four sites along one lattice direction. The
behavior in this case is expected to be substantially different from the
spiral liquid studied in Ref. [35]. Our core calculations are performed with
the following specifications. We consider square systems with periodic
boundary conditions, with up to $N=200\times 200=40000$ sites. A logarithmic
cool-down protocol is implemented with 120 temperature steps from $T=2$ down
to $T=0.02$ (where the units of energy are eliminated via $|J_{1}|=1$). Each
temperature step consists of $5\times 10^{5}$ Monte Carlo steps composed of
$N$ Metropolis trials and $N$ over-relaxation steps each. The acceptance ratio
remains close to $50\%$ thanks to the adaptive Gaussian step [8]. Data for $e$
and $c_{v}$ is collected after performing half of the cMC steps at a given
temperature. Results are averaged over 10 independent runs. Apart from these
core calculations, we perform additional calculations, e.g. for the spin
structure factors, starting from stored thermalized states.

Figure 3: Phase diagram of the Hamiltonian in Eq. (1) obtained by cMC. The
color indicates the size of $c_{v}$. Open symbols and lines indicate the
positions of the peaks or local maxima in $c_{v}$ (see main text), except for
the orange symbols that are calculated from an energy crossover. The type of
phase transition (for example, the universality class) is indicated in the
legend (SB = symmetry breaking), and the phases are indicated in the
corresponding regions (SL = spiral liquid, VN = vortex network). The red star
corresponds to the transition into the spiral liquid as obtained in Ref. [35].

The energy and specific heat calculations exhibit finite-temperature phase
transitions in the full range of investigated $\delta>0$. The phase diagram is
presented in Fig. 3, where the background color indicates the magnitude of the
specific heat $c_{v}(T,\delta)$. The colors of the lines and symbols indicate
the types of phase transitions manifested by peaks in the specific heat, which
we will explain in detail in the following subsections. Overall, the phase
diagram reveals many interesting features such as the existence of two phase
transitions (red and green lines) at $\delta\geq 0.16$, where the lower
critical temperature $T_{c}$ goes to zero as $\delta\to 0.25$. These
transitions surround a double-striped phase which we denote $\langle
2\rangle$. The central part of the phase diagram is dominated by
incommensurate lattice nematic spirals which can be understood as stripy spin
states in one of the two Cartesian directions. At small $\delta\lesssim 0.05$,
the system shows a crossover at high temperatures from the paramagnetic regime
into a pancake liquid phase [44]. Upon lowering the temperature further, the
system undergoes a transition into a spiral spin liquid phase (brown line)
without displaying any signs of symmetry breaking. The red star on this line
highlights the phase transition found in Ref. [35]. At even lower
temperatures, this phase turns into a rigid vortex network.

These phase transitions and crossovers depict a rich phase diagram including
the spiral spin liquid phase at finite temperatures, which we will
characterize in the following. Before we continue describing the individual
phases in detail, it is important to note that the Mermin-Wagner theorem
forbids the spontaneous breaking of the continuous U(1) spin rotation symmetry
at any finite temperature. Therefore, $\langle\mathbf{S}_{i}\rangle=0$ holds
for any site $i$ at $T>0$. This means that whenever we speak of ‘symmetry
breaking order’ we refer to orders that manifest in spin-rotation invariant
order parameters (such as spin correlations
$\langle\mathbf{S}_{i}\mathbf{S}_{j}\rangle$) and that break discrete lattice
symmetries.

###  III.1 Double colinear stripes

Let us start with the parameter regime corresponding to the $\langle 2\rangle$
phase (2-up-2-down phase), which consists of nematic two-site-wide stripes of
antiparallel spins [see Fig. 4 (a)]. This phase breaks both lattice
translation and rotation symmetries and is typically observed in systems with
short-ranged ferromagnetic and long-ranged antiferromagnetic interactions [45,
46, 47]. In the spin structure factor,

$S(\mathbf{q})=\frac{1}{N}\sum_{i,j}\
\langle\mathbf{S}_{i}\cdot\mathbf{S}_{j}\rangle\
e^{i\mathbf{q}\cdot\mathbf{r}_{ij}}$ (4)

where $\mathbf{r}_{ij}=\mathbf{r}_{j}-\mathbf{r}_{i}$ is the distance between
spins $\mathbf{S}_{i}$ and $\mathbf{S}_{j}$, the $\langle 2\rangle$ phase is
characterized by peaks at ${\mathbf{q}}=(\pm\pi/2,0)$ or $(0,\pm\pi/2)$ [see
Fig. 4 (b)]. It is important to note that these points only belong to the
spiral ground-state manifold at exactly $\delta=0.25$. For $\delta<0.25$ there
is another transition at lower temperatures that takes the system to the
ground-state manifold.

Figure 4: Characterization of phase $\langle 2\rangle$: (a) Snapshot of a
real-space spin configuration from cMC for $L=200$ at $\delta=0.17$ and
$T=0.119$ (the inset shows a zoom-in), where colors indicate the angle
$\phi_{i}$ at each site $i$. (b) Spin structure factor $S(\mathbf{q})$ (on a
logarithmic scale) for $\delta=0.17$ and $T=0.119$. (c) Order parameter
$\mathcal{O}_{1}$ to test for lattice translation symmetry breaking [see main
text and Eq. (6)].

The selection of q-points on the Cartesian axes in reciprocal space is
entropic and happens due to the discrete nature of the lattice. Naively,
breaking the discrete $\mathds{Z}_{2}$ lattice rotation symmetry that
transforms the wave vectors ${\mathbf{q}}=(\pm\pi/2,0)$ and $(0,\pm\pi/2)$
into each other would lead to an Ising transition, as in the well-known
$J_{1}$-$J_{2}$ square lattice Heisenberg model for $J_{2}>0.5~{}J_{1}$, where
simple antiferromagnetic stripes are selected [48, 49]. However, this is not
the case here since there is an extra accidental degeneracy that is
responsible for the formation of the $\langle 2\rangle$ phase. Specifically,
the Hamiltonian in Eq. (1) implies that all coplanar spin stripes with the
periodic four-site pattern
$[\phi_{A}\phi_{B}\overline{\phi}_{A}\overline{\phi}_{B}]$ along one of the
two Cartesian axis have the same energy regardless of the values of $\phi_{A}$
and $\phi_{B}$ ($\phi_{A}$ and $\phi_{B}$ indicate two different angles and
$\overline{\phi}_{A}=\phi_{A}+\pi$ is the opposite direction of $\phi_{A}$).
Within this manifold there are spin spirals with a period of four lattice
sites and spatially homogeneous $\pi/2$ rotation angles between neighboring
spins (i.e., $\phi_{A}$ and $\phi_{B}$ differ by $\pi/2$) as well as $\langle
2\rangle$ phases when $\phi_{A}=\phi_{B}$ or $\phi_{A}=\overline{\phi}_{B}$.
The entropic term of the free energy selects the $\langle 2\rangle$ states,
indicating a preference for colinear states. The two possibilities
$\uparrow\downarrow\downarrow\uparrow\uparrow\downarrow\downarrow\uparrow\cdots$
or
$\downarrow\downarrow\uparrow\uparrow\downarrow\downarrow\uparrow\uparrow\cdots$
of breaking the lattice translation symmetry represent an extra
$\mathds{Z}_{2}$ broken symmetry. Typically, there is no reason why the two
$\mathds{Z}_{2}$ symmetries should be broken at the same temperature. However,
we observe only one phase transition belonging to the 4-state Potts
universality class [50], confirmed by finite-size scaling from cMC
calculations (see Appendix B). This indicates a strong order parameter
coupling, leading to a broken
$\mathds{Z}_{2}\times\mathds{Z}_{2}=\mathds{Z}_{4}$ symmetry. Such an emergent
4-state Potts transition through order-by-disorder was already observed in the
$J_{1}$-$J_{3}$ kagome lattice Heisenberg model but as a result of a true
four-fold symmetry and not a combination of two coupled two-fold symmetries
[51]. One possible reason for the merging of the two types of symmetry
breaking is that the two order parameters are fundamentally dependent on each
other: Translation symmetry breaking along a lattice direction can only occur
if nematic symmetry breaking has selected this lattice direction.

As already mentioned, the phase $\langle 2\rangle$ breaks translation symmetry
by one lattice spacing. As we will see later, other phases of our system (such
as the nematic spiral) also break translation symmetry but remain invariant
under combined translation and spin rotation. What is special about the
$\langle 2\rangle$ phase is that it also breaks all possible combinations of
translation symmetry by one lattice spacing and global spin rotations. To
explicitly prove this it is useful to consider the following order parameter

$\mathcal{O}_{1}^{x}=\frac{1}{N}\left|\sum_{x,y}(-1)^{x}\langle\mathbf{S}_{x,y}\cdot\mathbf{S}_{x+1,y}\rangle\right|,$
(5)

where we have slightly changed the notation in the sense that the subscripts
$x$, $y$ in $\mathbf{S}_{x,y}$ now indicate the Cartesian coordinates of the
square lattice sites $i$. This order parameter takes into account all nearest
neighbor correlations in the $x$ direction with alternating signs.
Analogously, $\mathcal{O}_{1}^{y}$ can be defined that detects translation
symmetry breaking by one lattice spacing in the $y$-direction. To take into
account both types of translation symmetry breaking along the $x$ and $y$ axes
we consider the order parameter

$\mathcal{O}_{1}=\mathcal{O}_{1}^{x}+\mathcal{O}_{1}^{y}.$ (6)

It is easy to see that for an ideal $\langle 2\rangle$ state (with either
stripes along the $x$ or $y$ directions), $\mathcal{O}_{1}=1$. In the same
manner, one can check that $\mathcal{O}_{1}=0$ for ferromagnetic,
antiferromagnetic Néel, and spiral states with homogeneous rotation angles.

The results for $\mathcal{O}_{1}$ as a function of $T$ and $\delta$ are shown
in Fig. 4 (c). The order parameter $\mathcal{O}_{1}$ takes a finite value only
in the region encapsulated by the green and red curves at high values of
$\delta$, that has been denoted $\langle 2\rangle$ in Fig. 3. Even though
favored by thermal fluctuations, the $\langle 2\rangle$ phase (and more
generally, all $[\phi_{A}\phi_{B}\overline{\phi}_{A}\overline{\phi}_{B}]$
states) is not a part of the ground-state manifold for $\delta<0.25$.
Therefore, we cannot strictly speak of order by disorder due to thermal
fluctuations. As a consequence, another phase transition must occur at lower
temperatures for $\delta<0.25$, taking the system to the ground-state
manifold. This is confirmed by the recovery of the translation symmetry at
lower temperatures (below the green line), where $\mathcal{O}_{1}$ takes again
small values, evidencing a re-entrance phenomenon.

We note that there is still some reminiscence of the $\langle 2\rangle$ phase
after the system passes through it at lower temperatures, manifested by the
non-zero value of $\mathcal{O}_{1}$ below the green line at high values of
$\delta$. The boundary conditions are responsible for this spurious effect, as
the system tries to evolve from a 4-spin periodic structure ($\langle
2\rangle$ phase) onto an incommensurate spiral phase. The frustration induced
by the boundary conditions may forbid the system from fully re-arranging into
a nematic spiral. Another consequence of this is the larger energy difference
between these states when we reach $T=0$ and the exact ground-state energy
(see Appendix A).

###  III.2 Nematic coplanar spirals

Figure 5: Characterization of the nematic spiral phase: (a) Real-space
configuration for $L=200$ at $\delta=0.14$ and $T=0.0001$ (the inset shows a
zoom-in), where colors indicate the angle $\phi_{i}$ at each site $i$. (b)
Spin structure factor (in logarithmic scale) for $\delta=0.14$ and $T=0.02$.
(c) Order parameter $\mathcal{O}_{2}$ to test for lattice rotation symmetry
breaking (see main text).

The nematic spiral phase exists below the green curve in the phase diagram of
Fig. 3 and only breaks the lattice rotation symmetry by choosing two
equivalent points on the spiral ring [see Fig. 5 (a) and (b)]. The spiral
states in this phase are characterized by spatially homogeneous rotation
angles between spins on neighboring sites. For $\delta\geq 0.16$, the nematic
spiral phase is reached through the recovery of the translation symmetry that
is broken in the $\langle 2\rangle$ phase. The corresponding transition is
found to belong to the Ising universality class (see Appendix B) which becomes
evident when approaching the phase boundary from the low-temperature side.
Below the phase transition, the system only breaks $\pi/2$ lattice rotation
symmetry while translation symmetry (in combination with a properly chosen
spin rotation) is intact. On the other hand, above the phase transition (i.e.
in the $\langle 2\rangle$ phase), an extra $\mathds{Z}_{2}$ translation
symmetry is broken, explaining the Ising nature of the transition.

For $0.05\lesssim\delta\leq 0.15$, the nematic spiral phase is directly
reached from the paramagnetic regime, populating a large region of the phase
diagram. In this $\delta$ region, the transition is also in the Ising
universality class due to the broken $\mathds{Z}_{2}$ lattice rotation
symmetry. For most values of $\delta$, these spirals are incommensurate since
the wavelength changes smoothly with $\delta$. For example, commensurate
spirals with wavelengths of 4, 5, and 6 lattice spacings are realized for
values $\delta$ given by 0.25, $\sim 0.131$, and $\sim 0.083$, respectively
[see Eq. (2)]. As a consequence, finite lattices have an extra degree of
frustration due to the boundary conditions which lead to deviations from the
ideal wave vectors $(q_{x},0)$ and $(0,q_{y})$. These effects, however, vanish
in the thermodynamic limit where order-by-disorder effects can select the
exact wave vectors $(q_{x},0)$ and $(0,q_{y})$ on the spiral ring.

An order parameter that detects the lattice rotation symmetry breaking can be
constructed in terms of the local momentum (i.e. local wave vector) of the
spiral. At each site, we can define the local momentum via
$\mathbf{q}=\nabla\phi$ [35] (where $\nabla$ is implemented as a discrete
derivative on the square lattice) which is a two-component vector in the
$xy$-plane. Then, our order parameter for lattice rotation symmetry breaking
is defined by

$\mathcal{O}_{2}=\frac{1}{N}\left|\sum_{i\in
N}\langle\left|\hat{\mathbf{q}}_{i}^{x}\right|\rangle-\langle\left|\hat{\mathbf{q}}_{i}^{y}\right|\rangle+\langle\left|\hat{\mathbf{q}}_{i}^{x+y}\right|\rangle-\langle\left|\hat{\mathbf{q}}_{i}^{x-y}\right|\rangle\right|,$
(7)

where $\hat{\mathbf{q}}_{i}=\mathbf{q}_{i}/|\mathbf{q}_{i}|$ is the normalized
unit vector of $\mathbf{q}_{i}$ at site $i$ of the square lattice.
Furthermore, the superscripts indicate the components of
$\hat{\mathbf{q}}_{i}$ where $x+y$ and $x-y$ are symbolic notations for the
two diagonal lattice directions. Even though these two directions are not
favored by entropy, we include them to account for canted stripes that appear
due to finite-size effects. In practice, this makes the results slightly
smoother at small $\delta$. The $\mathcal{O}_{2}$ order parameter vanishes
when there is no imbalance between the total momentum of the spirals in the
$x$ and $y$ directions (or the $x+y$ and $x-y$ directions). On the other hand,
$\mathcal{O}_{2}=1$ if a spiral is aligned along a Cartesian or a diagonal
direction.

The results for $\mathcal{O}_{2}$ are shown in Fig. 5 (c). We find that
$\mathcal{O}_{2}$ is always sizeable below the green curve. This confirms that
the peak in the heat capacity at $0.05\lesssim\delta\lesssim 0.15$, identified
as an Ising transition, is due to lattice rotation symmetry breaking. The
small but non-zero signal observed at $\delta\lesssim 0.05$ can be attributed
to the phase’s slow dynamics (or bad thermalization) and a small number of
independent runs trapped in a nematic state near the orange boundary. The
$\mathcal{O}_{2}$ order parameter provides evidence of the existence of the
phase transition shown in orange, which is not captured by the specific heat
and which we will discuss in detail later. It is also important to note that
the poor signal obtained in the $\langle 2\rangle$ region originates from the
discrete definition of the momentum. Due to the colinear double-stripe nature
of the phase, the relative angles $\phi_{i}-\phi_{j}$ between neighboring
sites $i$ and $j$ vary vastly on the scale of one lattice distance, making
discrete derivatives ill-defined. Still, the non-zero values of
$\mathcal{O}_{2}$ confirm that the $\langle 2\rangle$ phase also breaks the
lattice rotation symmetry.

###  III.3 Pancake and spiral liquids

Figure 6: Pancake liquid phase: (a) Real-space configuration for $L=200$ at
$\delta=0.04$ and $T=0.119$ (the inset shows a zoom-in), where colors indicate
the angle $\phi_{i}$ at each site $i$, and (b) the corresponding spin
structure factor (in logarithmic scale). Spiral liquid phase: (c) Real-space
configuration for $\delta=0.04$ and $T=0.081$ where one momentum vortex and
antivortex are indicated by an up and down arrow, respectively, and (d)
corresponding spin structure factor.

Finally, for $\delta\lesssim 0.06$, the heat capacity peak indicating a
finite-temperature phase transition splits into two, a wide shoulder at higher
$T$ signaling a crossover and a sharp peak at lower $T$ that gets smaller and
shifts to $T=0$ as $\delta$ decreases (brown curve in Fig. 3). The broad
shoulder, shown by the light-blue curve in Fig. 3, is only discernible for
$\delta\leq 0.05$ and does not exhibit critical behavior (no scaling with
system size) indicating a crossover. Generally, this part of the phase diagram
is well represented by the results for $\delta=0.03$ in Ref. [35]. The onset
of spin correlations at the broad high-temperature peak leads to a phase known
as pancake liquid [44], where the spin states have contributions from spiral
wave vectors ranging from $\mathbf{q}=0$ up to the spiral ring [see Fig. 6 (a)
and (b)]. This property leads to a spin structure factor featuring an almost
uniform disk-like signal (pancake shape) within the spiral surface.

Further decreasing the temperature leads to a phase transition into the spiral
spin liquid phase [see Fig. 6 (c) and (d)], which does not show any symmetry
breaking. Contrary to the pancake liquid, this phase presents a nearly
isotropic signal in the spin structure factor only along a circle inside the
spiral ring [highlighted by dashed white lines in Fig. 6 (d)]. This phase
shows liquid-like fluctuations between spin spirals with all possible
directions of spiral wave vectors, justifying the name spiral spin liquid.
These fluctuations give rise to well-defined momentum vortices and
antivortices representing local defects in the spiral configurations [35] [see
up and down arrows in Fig. 6 (c)]. The circle formed in the spin structure
factor slightly enlarges as the temperature decreases until it reaches the
spiral surface corresponding to the ground-state manifold. At the same time,
the dynamics of the phase are slowed down and thermalization becomes difficult
in the cMC runs. The transition separating the pancake and spiral spin liquids
has a logarithmically divergent specific heat (see Appendix B), corresponding
to a critical exponent $\alpha=0$. Typically, $\alpha=0$ could be associated
with an Ising transition, as it is found for $\delta>0.05$. However, in this
case, the critical exponent for the correlation length $\nu$ is not consistent
with an Ising universality class, and the spiral spin liquid does not exhibit
any signs of symmetry breaking. This evidences that the transition falls out
of the standard paradigm of phase transitions, and we will argue in Sec. V
that it is possibly a Kosterlitz-Thouless transition driven by momentum
vortices.

Figure 7: Energy-level crossings between the nematic spiral and spiral
liquid/vortex network phases. Dashed lines indicate decreasing $\delta$-sweeps
starting from $\delta=0.07$. Full lines indicate increasing $\delta$-sweeps
starting from $\delta=0.03$. Black dots indicate the crossing points.

As mentioned before, two different phases are observed at low temperatures
around $\delta\approx 0.05$. On one hand, we find nematic spirals for
$\delta\gtrsim 0.05$, which break the lattice rotation symmetry (green curve
in Fig. 3). On the other hand, we find a spiral liquid for $\delta\lesssim
0.05$, which does not break any symmetries (brown curve in Fig. 3). Thus, a
phase transition must exist between these two phases. However, no such phase
transition is captured by the specific heat, even when very small $\delta$
steps are considered (see Appendix C). Furthermore, performing $\delta$-sweeps
at constant $T$, coming from either of the two phases, we could not tune the
system from one phase into the other. This can be attributed to the large
change of spin configuration needed to go from one phase to the other, the
slow dynamics of spin fluctuations at these temperatures, and/or the
frustration generated by the boundary conditions. Other update methods for cMC
could provide a solution to this challenging problem which, however, we leave
for future work.

Despite the difficulties in evolving the two phases into each other, we can
infer the existence of a phase transition from energy considerations.
Specifically, Fig. 7 shows an energy-level crossing in our $\delta$-sweeps at
constant $T$, when either starting in the nematic spiral phase and decreasing
$\delta$ (dashed lines) or starting in the spiral liquid/vortex network phase
and increasing $\delta$ (full lines). The crossing points are denoted by black
dots and are shown by orange symbols in Fig. 3. This line roughly agrees with
the region in the phase diagram where the order parameter $\mathcal{O}_{2}$,
detecting lattice rotation symmetry breaking, vanishes. Nonetheless, it is
worth noticing that this method of locating the phase transition is only
strictly valid at $T=0$ or when the two phases have the same entropy (which is
probably not the case here). Therefore, the orange line in Fig. 3 should only
be understood as a rough estimate for the transition. Since we could not
observe the system’s evolution from one phase to the other, we could also not
detect any hysteresis behavior.

###  III.4 Vortex network or ripple state

As mentioned in the previous section, the ring-like signal in the spin
structure factor characterizing the spiral liquid (see Fig. 6) grows as the
temperature is lowered until it reaches the spiral ring of the exact ground
state manifold. Around this temperature, a crossover to a rigid vortex network
state occurs, where the spirals are well-defined and correspond to four
different wave vectors $\mathbf{q}=(\pm q,0)$ and $(0,\pm q)$ [35] selected by
an order-by-disorder mechanism. A typical spin configuration is shown in Fig.
8 (a), with the corresponding spin structure factor presented in Fig. 8 (b).
This rigid vortex network does not show any indications of symmetry breaking
and realizes an approximate square arrangement of momentum vortices and
antivortices (therefore the name vortex network) connected by straight domain
walls [35]. Furthermore, these structures are rigid, in the sense that they
evolve very slowly in Monte Carlo time. This makes the system very prone to
getting stuck in a metastable spin configuration. For example, while the wave
vectors $\mathbf{q}=(\pm q,0)$ and $(0,\pm q)$ are favored by entropic effects
at finite temperature, we still find configurations with $\mathbf{q}=\pm(q,q)$
and $\pm(q,-q)$, indicating that the cMC runs are affected by thermalization
issues, as well as by finite sizes and boundary conditions.

Figure 8: Vortex network phase: (a) Real-space configuration for $L=200$ at
$\delta=0.03$ and $T=0.02$, and (b) the corresponding spin structure factor.
(c) Real-space configuration for $\delta=0.03$ and $T=0.0001$, where the black
arrow shows an elongated momentum antivortex. (d) Ripple state obtained for
$\delta=0.03$ and $T=0.02$ with open boundary conditions.

Other studies of a closely related model on the honeycomb lattice found that
the low-temperature phase is a ripple state for systems with open boundary
conditions [44]. The ripple state consists of a single momentum vortex from
which spin spirals extend in all directions. Interestingly, we also find such
a state in the vortex network regime when imposing open boundary conditions
[see Fig. 8 (d)]. The discrepancy between the two boundary conditions
indicates that they play an essential role in determining the bulk spin
configurations even for large lattice sizes. This applies in particular to the
limit $\delta\to 0$ where the spiral wavelength diverges such that it is not
possible to make a sensible extrapolation to the thermodynamic limit.
Nonetheless, at finite temperatures, a ripple state in the thermodynamic limit
seems peculiar because this would imply that spins far away from the ripple
center are affected by it, regardless of the distance. In what follows, we
provide a simple understanding of why the single vortex state is favored by
the open boundary condition.

Figure 9: Bottom part of a lattice with an open boundary and the remaining
couplings to the system of a spin at the edge (black circle).

As an example, we examine the bottom boundary of a square-shaped system (i.e.,
a boundary oriented along the $x$-axis) and investigate the Hamiltonian at the
boundary. The other three boundaries of the system can be treated analogously.
Assuming that the spins are in a spiral state with wave vector
$\mathbf{q}=(q_{x},q_{y})$, the energy of a spin closest to the boundary can
be obtained by computing the Hamiltonian in Eq. (1) for the spiral state but
omitting the contributions from the bonds that are cut by the boundary. This
leads to an energy contribution from such a spin given by (see Fig. 9)

$\begin{split}2e_{\text{bdry}}=&2J_{1}\cos{(q_{x}a)}+J_{1}\cos{(q_{y}a)}\\\
&+J_{2}\cos{(q_{x}a+q_{y}a)}+J_{2}\cos{(-q_{x}a+q_{y}a)}\\\
&+2J_{3}\cos{(2q_{x}a)}+J_{3}\cos{(2q_{y}a)},\end{split}$ (8)

where we recall that $J_{1}=-1$, $J_{2}=1/4+\delta$, $J_{3}=1/8+\delta/2$ and
$a=1$ is the nearest-neighbor distance. Minimizing $e_{\text{bdry}}$ as a
function of $q_{x}$ and $q_{y}$ leads to a wave vector
$\mathbf{q}^{\prime}=(0,\pm q_{y}^{\prime})$ with $q_{y}^{\prime}$ given by

$\cos(q_{y}^{\prime}a)=\frac{1-4\delta}{1+4\delta}.$ (9)

Note that since $\cos(0)+\cos(q_{y}^{\prime}a)=\frac{2}{1+4\delta}$, the wave
vector $(0,\pm q_{y}^{\prime})$ is exactly on the spiral ring [see Eq. (2)].

This analysis shows that open boundaries energetically prefer spirals whose
stripe-like spin patterns are parallel to the boundary, as is (approximately)
realized for the ripple state in Fig. 8 (d). While the bulk energy is
minimized by a domain with a single wave vector $\mathbf{q}$ on the spiral
ring, such a state costs much energy on the boundary (growing linearly with
the linear system size $L$) since a single domain cannot simultaneously
satisfy the energetic preferences at all boundaries. Hence, the system forms a
single momentum vortex in the bulk which reduces the energy costs at the
boundary (note that the bulk energy of a single momentum vortex only scales
logarithmically in the system size $L$ [35]). This mechanism stabilizes the
single vortex ripple state in Fig. 8 (d) at low temperatures when open
boundary conditions are imposed.

Upon further decreasing the temperature, systems with periodical boundary
conditions show a re-organization of momentum vortices associated with a
bending of the spiral wavefronts. This leads to spin configurations shown in
Fig. 8 (c) obtained at $T=0.0001$. The bending of the wavefronts gives rise to
momentum antivortices with a very elongated shape along one Cartesian
direction [shown by the black arrow in Fig. 8 (c)]. Selecting one of the two
lattice directions corresponds to a $\mathds{Z}_{2}$ symmetry breaking and
implies a certain degree of lattice nematicity in the system. This spontaneous
breaking of a $\mathds{Z}_{2}$ symmetry is expected at low but finite
temperatures. However, due to the aforementioned thermalization issues in the
vortex network phase, most independent cMC runs do not present any signatures
of this transition in the specific heat such that the numerical evidence for
the spontaneous symmetry breaking remains rather weak. On the other hand, the
decreasing density of momentum vortices at low temperatures is a very robust
observation. When the distance between momentum vortices becomes comparable to
the system size, the results are largely affected by the boundary conditions
and therefore cannot be trusted. We thus conclude that the vortex network
regime is the numerically most challenging part of the phase diagram and a
thorough understanding of this phase is still lacking.

##  IV Dynamics

To study the signatures of the different phases in the dynamical spin
structure factor which is experimentally accessible through inelastic neutron
scattering, we perform molecular dynamics (MD) calculations. These simulations
start from XY spin configurations obtained from cMC for a set of coupling
parameters at a given temperature and then calculate the system’s time
evolution by solving the Landau-Lifshitz equations (without a damping term)
[12, 52]. The dynamical spin structure factor is then calculated as

$S(\mathbf{q},\omega)=\frac{1}{N\sqrt{N_{t}}}\sum_{n_{t}=0}^{N_{t}}\sum_{i,j}^{N}\
\langle\mathbf{S}_{i}(t)\cdot\mathbf{S}_{j}(0)\rangle\
e^{i(\mathbf{q}\cdot\mathbf{r}_{ij}-\omega t)},$ (10)

where the time is given by $t=n_{t}\delta t$, for a given time step $\delta t$
[53, 54]. It is important to note that the Landau-Lifshitz dynamics are only
well-defined for three-component spins. Hence, for the MD calculations, we use
the XY Hamiltonian from Eq. (1) but let it act on three-component spins [55],
i.e., the spin’s $z$-component (which vanishes in the initial spin state from
cMC) does not contribute to the energy. During the MD runs, the energy per
spin is conserved up to the fifth decimal digit, indicating that the time
evolution is well performed. The dynamical spin structure factors in Figs. 10
and 11 are obtained by averaging over 120 independent initial configurations.
We point out that the vertical high-intensity signals extending from
$\omega=0$ to the upper boundary of the plotted regions are artifacts
generated by the Fourier transform in time, which can be explained as follows.
The mismatch between the period of spin oscillations and the simulation time
generates a finite offset in the Fourier-transformed signal. Although this
offset is much smaller than the height of the peak, it is visible for the
momenta where the intensity is highest. The MD results are compared with the
zero-temperature spin wave bands from linear spin-wave theory (LSWT),
calculated with the software SpinW [56].

Figure 10: Spin structure factor $S(\mathbf{q},\omega)$ from MD calculations
at $\delta=0.18$ for two different temperatures, corresponding to the colinear
$\langle 2\rangle$ phase (top panel) and the coplanar nematic spiral phase
(bottom panel). The dashed white lines are the spin wave bands obtained by
LSWT at $T=0$, and the dashed blue vertical lines indicate momenta on the
spiral ring. The inset shows the plotted path in momentum-space (red) and the
spiral ring corresponding to the ground-state manifold (blue). The white
dashed boxes indicate spurious ghost signals (see main text).

In Fig. 10 we show results for $\delta=0.18$ at two different temperatures,
above and below the low-temperature phase transition between the $\langle
2\rangle$ phase and the nematic spiral phase. The bottom panel shows results
for the nematic spiral phase, where the magnetic Bragg peak occurs at a wave
vector on the spiral ring close to $\mathbf{q}=(0,\pi/2)$. Momenta on the
spiral ring that lie on the plotted path are indicated by dashed blue lines.
As a reference, we also plot the spin wave bands obtained with LSWT, shown as
dashed white lines in Fig. 10. Overall, the spectral weight obtained within MD
closely follows the dispersive bands from LSWT. Another relevant feature
revealed by both methods is that, away from the Bragg peak the spectrum is
always gapped, even for momenta $\mathbf{q}$ on the spiral ring, which is a
consequence of the XY-anisotropy. A band minimum with a small but finite gap
is found close to $\mathbf{q}=(0,\pi/2)$ and is related to the proximity to
the $\langle 2\rangle$ phase at $\delta=0.25$. Finally, it should be noted
that some lower energy “ghost” excitations are also observed, highlighted by
dashed white rectangles. These are spurious signals originating from the
passing through the colinear $\langle 2\rangle$ phase, as was already found in
the translation symmetry-breaking order parameter $\mathcal{O}_{1}$ (see Fig.
4). These features should be disregarded as remnants from the $\langle
2\rangle$ phase.

The top panel of Fig. 10 corresponds to the higher-temperature colinear
$\langle 2\rangle$ phase, where the magnetic Bragg peak is located at
$\mathbf{q}=(0,\pi/2)$. Note that in this case, LSWT cannot be applied because
the order is unstable at $T=0$ (it does not belong to the ground-state
manifold for $\delta\neq 0.25$). In addition to the gapless excitations at the
Bragg peak, there are further low-energy band minima at the wave vectors where
the plotted path crosses the spiral ring (indicated by dashed blue lines). In
contrast to the coplanar nematic spiral phase, the spectrum is more dispersed
and the weight is not concentrated along well-defined narrow excitation modes.
Overall, apart from the low-energy signal around the Bragg peak, the intensity
distributions in the colinear $\langle 2\rangle$ and the coplanar spiral
phases are rather different. This is a remarkable observation given that the
wave vectors of the two corresponding orders are very similar.

Figure 11: Spin structure factor $S(\mathbf{q},\omega)$ from MD calculations
at $\delta=0.03$ for two different temperatures, corresponding to the pancake
liquid (top panel) and the spiral spin liquid (bottom panel). The white and
green dashed lines correspond to the spin wave bands obtained by LSWT at $T=0$
for ground state nematic spirals with wave vectors $\mathbf{q}=(q,0)$ and
$\mathbf{q}=(0,q)$ on the spiral ring, respectively. The dashed blue vertical
lines indicate the momenta on the ground state spiral ring. The plotted path
in $\mathbf{q}$-space is shown in the inset (red), together with the spiral
ring (blue). Note that the color scale is logarithmic.

We now turn to the other end of the phase diagram at $\delta=0.03$, where the
pancake and spiral spin liquids are found. Figure 11 shows the MD results for
the dynamical spin structure factor for two different temperatures
corresponding to both phases. The pancake liquid phase receives its name due
to the nearly homogeneous contributions from all wave vectors inside the
spiral ring. This property is also observed in the dynamical spin structure
factor of the pancake liquid, displayed in the top panel of Fig. 11. The
intensity distribution shows a strong signal at low energies for all wave
vectors enclosed by the spiral ring (i.e., between the dashed blue lines)
without any particular Bragg peak. Here, the white and green dashed lines
represent the LSWT bands at $T=0$ for wave vectors $\mathbf{q}=(q,0)$ and
$\mathbf{q}=(0,q)$ on the spiral ring, respectively. In Fig. 11 one can
appreciate that the spectral weight of the pancake liquid follows these bands
while displaying a continuum of excitations in the region enclosed by them.

When decreasing the temperature at $\delta=0.03$, the system enters the spiral
spin liquid phase where only wave vectors close to the spiral ring coexist.
This implies that spin spirals are well-defined but their momentum direction
can vary in real space. The dynamical spin structure factor of the spiral spin
liquid in the bottom panel of Fig. 11 displays this property where the low
energy signal is now more concentrated around the spiral ring (blue dashed
lines). Also at higher energies the spectral weight is sizeable only close to
the $T=0$ LSWT bands, while the region enclosed by them remains empty,
representing a significant difference from the pancake liquid. Overall these
features help to distinguish the two types of spin liquids by their excitation
spectrum.

##  V Effects of spin and momentum vortices

As mentioned in Sec. III.3, the phase transition between the pancake liquid
and the spiral spin liquid shows unusual properties. It is not associated with
any spontaneous symmetry breaking and its scaling behavior falls out of the
standard classification of second-order phase transitions. Specifically, the
specific heat is logarithmically divergent (see Appendix B) but the critical
exponent $\nu$ is not consistent with the known universality classes. Here, we
revisit the nature of this phase transition and discuss whether it could be
driven by vortices, either formed by spin or momentum. Since both degrees of
freedom are planar quantities, a Kosterlitz-Thouless transition seems to be a
natural possibility [57, 58].

Figure 12: Top panel: Spin vortex density $n_{v}$ (defined by the number of
spin vortices per number of sites) as a function of temperature $T$ and
$\delta$. Bottom panel: Derivative of $n_{v}$ with respect to the temperature.
The lines and circles indicate the phase boundaries obtained from peaks in the
specific heat.

To obtain insights into whether spin or momentum vortices show any special
behavior across the phase transition, we first examine the density of vortices
within cMC. We start investigating spin vortices which are identified by the
total spin rotation angle along elementary four-site loops of the square
lattice and whose density we denote $n_{v}$. The results for $n_{v}$ are shown
in the top panel of Fig. 12 while the bottom panel displays its derivative
with respect to temperature, $\partial n_{v}/\partial T$. In the small
$\delta$-regime the spin-vortex density $n_{v}$ shows a pronounced drop at
rather large temperatures $T\sim 0.4$. As a result of this pronounced
decrease, at the low-temperature phase transition between the pancake liquid
and the spiral spin liquid (brown line) $n_{v}$ is negligibly small. This
makes it evident that the transition is not related to spin vortices and we,
therefore, rule out a Kosterlitz-Thouless transition from the spin degree of
freedom. Furthermore, $n_{v}$ increases with $\delta$, and close to the Ising
transition indicated by the green line the system is still populated by a
considerable number of spin vortices.

It is worth emphasizing that our investigation of the vortex density is only
to check whether there is a general connection between vortices and the spiral
liquid transition. The vortex density is expected to be smooth across a
Kosterlitz-Thouless transition [59, 60], however, the derivative has been
shown to exhibit a peak at the same position as that of the specific heat
[59]. Despite coinciding, these peaks occur at temperatures slightly above the
Kosterlitz-Thouless transition and do not diverge on the square lattice XY
model [61]. Thus, $n_{v}$ is no good quantity to precisely locate a
Kosterlitz-Thouless transition. We, nevertheless, study the vortex density
here since more standard quantities for identifying a Kosterlitz-Thouless
transition, such as the correlation length, are difficult to calculate
accurately due to the spiral nature of the magnetic parent state (which leads
to oscillating correlation functions).

Next, we investigate the density of momentum vortices $n_{qv}$ in the
$\delta$-$T$ phase diagram. If $\phi(\mathbf{r})$ evolves smoothly throughout
the lattice, we can define the momentum field
$\mathbf{q}(\mathbf{r})=\nabla\phi(\mathbf{q})$ locally as the discrete
derivative of $\phi(\mathbf{r})$ and identify momentum vortices in the local
momentum texture $\mathbf{q}(\mathbf{r})$. The results are shown in the top
panel of Fig. 13 while the derivative $\partial n_{qv}/\partial T$ is
presented in the bottom panel. One sees that for all values of $\delta$, the
momentum vortices have a strong presence at temperatures above the phase
transitions. Furthermore, $n_{qv}$ drops abruptly at the transitions as
evidenced by its derivative, which is peaked along the brown line. At
$\delta=0$, where there is no finite-temperature phase transition, the
momentum vortex density remains high down to the lowest temperatures. This
observation indicates that momentum vortices play a key role in driving the
phase transition between the pancake liquid and the spiral spin liquid.

Figure 13: Top panel: Momentum vortex density $n_{qv}$ (defined by the number
of momentum vortices per number of sites) as a function of temperature $T$ and
$\delta$. Bottom panel: Derivative of $n_{qv}$ with respect to the
temperature. The lines and circles indicate the phase boundaries from the
specific heat.

The close connection between the spiral liquid transition and momentum vortex
proliferation motivates us to analyze more closely whether this transition
could be a Kosterlitz-Thouless transition of momentum vortices. To do so, we
investigate the behavior of the momentum-momentum correlation length
$\xi_{qq}$, which is expected to diverge at a Kosterlitz-Thouless transition
$T_{\text{KT}}$ according to [58, 61]

$\xi_{qq}(T)\sim e^{{b\sqrt{\frac{T_{\text{KT}}}{T-T_{\text{KT}}}}}}$ (11)

when approaching $T_{\text{KT}}$ from above. Here, $b$ is a non-universal,
dimensionless number. To obtain $\xi_{qq}(T)$ we calculate the momentum-
momentum correlation function
$\langle\mathbf{q}(\mathbf{r})\cdot\mathbf{q}(\mathbf{r}^{\prime})\rangle$ for
distances $\mathbf{r}-\mathbf{r}^{\prime}$ along the $x$ and $y$ Cartesian
directions. We average the correlation function over the whole lattice and
over 10 independent cMC runs. The results are shown by circles in the top
panel of Fig. 14 for $L=200$ and $\delta=0.03$, and the colors indicate the
temperature ranging from 0.409 (red) to 0.029 (blue). The black circles
indicate the temperature of the peak in the specific heat for this lattice
size. The lines correspond to exponential fits, from which we extract the
correlation length $\xi_{qq}$ whose temperature dependence is shown in the
inset. It becomes evident that there is a sudden growth of the correlation
length close to the critical temperature (indicated by the dashed black line),
which can be interpreted as the start of a divergence. The expected divergence
in the correlation length for the Kosterlitz-Thouless transition corresponds
to a power-law decay in the entire temperature region $T<T_{\text{KT}}$.
However, in our case, a clean divergence $\xi_{qq}\rightarrow\infty$ cannot be
detected due to the finite size of our system. The pink curve in the inset of
Fig. 14 shows a fit of $\xi_{qq}$ in the range from $T=0.204$ to $0.084$ above
the critical temperature to the aforementioned functional dependence of
$\xi_{qq}(T)$ according to the Kosterlitz-Thouless theory. Apart from the
high-temperature regime and the temperature region near the observed
transition, the fit shows good agreement with our data. Finally, the bottom
panel of Fig. 14 displays the behavior of $\xi_{qq}$ across the whole phase
diagram, showing the sudden increase of $\xi_{qq}$ at the phase transition in
a wider $\delta$ region.

Figure 14: Top panel: Decay of momentum-momentum correlations
$\langle\mathbf{q}(\mathbf{r})\cdot\mathbf{q}(\mathbf{r}^{\prime})\rangle$ as
a function of the distance between sites $|\mathbf{r}-\mathbf{r}^{\prime}|$ at
$\delta=0.03$. The temperature ranges from $T=0.409$ (red) to 0.029 (blue),
and lines correspond to exponential fits (black lines and points indicate data
near the critical temperature). The extracted correlation lengths $\xi_{qq}$
are shown in the inset and fitted by the pink curve between $T=0.204$ and
$0.084$. Bottom panel: Momentum-momentum correlation length $\xi_{qq}$ as a
function of temperature $T$ and $\delta$.

Overall, our numerical results are consistent with a Kosterlitz-Thouless
transition associated with momentum vortices, but cannot ultimately resolve
the nature of this transition. If present, such a transition would constitute
a rather uncommon and previously unexplored occurrence of Kosterlitz-Thouless
physics. In fact, momentum degrees of freedom are quite different from the
usual microscopic conditions of a Kosterlitz-Thouless transition, which, at
first sight, rather speaks against a Kosterlitz-Thouless transition from
momentum vortices. First, the system does not have an exact global U(1)
rotation symmetry of the momentum $\mathbf{q}$. Rather, an exact energy-
conserving U(1) transformation of spiral momenta is only possible in the
ground state manifold. Second, from
$\mathbf{q}(\mathbf{r})=\nabla\phi(\mathbf{q})$ it follows that the momentum
field $\mathbf{q}(\mathbf{r})$ is curl-free,
$\nabla\times\mathbf{q}(\mathbf{r})=0$, a condition that does typically not
exist for more standard planar degrees of freedom such as XY spins. Third, as
a result of the last property, it was shown in Ref. [35] that momentum
antivortices have a larger excitation energy than momentum vortices giving
rise to a peculiar disparity between both vortex types. These arguments led
the authors of Ref. [35] to be reluctant to the possibility of a Kosterlitz-
Thouless transition from momentum vortices. If our present numerical
indications in favor of such a transition are true, this would indicate a
striking robustness of Kosterlitz-Thouless physics. However, to eventually
resolve this question, further studies are necessary.

##  VI Effective rank-2 electrodynamics and pinch-points

An effective continuum theory for a spiral spin liquid was recently derived in
the small $\delta$ limit where the spiral ring is approximately circular [35].
One of the assumptions of this theory is that the spin texture
$\phi(\mathbf{r})$ is a smooth function in real space so that the momentum
field $\mathbf{q}(\mathbf{r})=\nabla\phi(\mathbf{r})$ is curl-free,
$\nabla\times\mathbf{q}(\mathbf{r})=0$. Note that this excludes the
possibility of spin vortices which represent a local source of curl,
$\nabla\times\mathbf{q}(\mathbf{r})=\delta(\mathbf{r})$. The Hamiltonian of
the continuum theory for small $\delta$ and $q=|\mathbf{q}|$ reads as

$\mathcal{H}=\int d^{2}\mathbf{r}\left(\frac{q^{4}}{16}-2\delta
q^{2}\right)+\int
d^{2}\mathbf{r}\left(\mathcal{Q}_{\mu\nu}C_{\mu\nu\rho\sigma}\mathcal{Q}_{\rho\sigma}\right),$
(12)

where
$C_{\mu\nu\rho\sigma}=\frac{1}{16}\left(\delta_{\mu\rho}\delta_{\nu\sigma}+\delta_{\mu\sigma}\delta_{\nu\rho}-\delta_{\mu\nu}\delta_{\rho\sigma}\right)$
is a combination of Kronecker deltas and
$\mathcal{Q}_{\mu\nu}=\partial_{\mu}\partial_{\nu}\phi$ is the Hessian matrix
of $\phi(\mathbf{r})$ [35]. The Hamiltonian in Eq. (12) has two terms [35];
the first term is an effective potential for spin spirals that governs the
energy cost of a homogeneous spiral with momentum $\mathbf{q}$. The second
term describes a spiral stiffness that captures the energy cost of deforming
spiral configurations in real space. Importantly, within this effective theory
a matrix-valued field $E_{\mu\nu}$ can be defined via

$E_{\mu\nu}=\epsilon_{\mu\rho}\epsilon_{\nu\sigma}\mathcal{Q}_{\rho\sigma},$
(13)

that is subject to the constraint

$\partial_{\mu}\partial_{\nu}E_{\mu\nu}=0,$ (14)

again under the assumption that spin vortices are absent. In Eq. (13),
$\epsilon_{\mu\rho}$ denotes the Levi-Civita symbol. The condition in Eq. (14)
can be interpreted as a generalized Gauss law in a charge-free rank-2
electrodynamics theory, where the emergent electric field $E_{\mu\nu}$ is a
rank-2 tensor as opposed to the vector-valued electric field in conventional
electrodynamics. Inspired by the unusual kinetic properties of the associated
charged matter fields (in this case called fractons) which retain their
mobility only within subdimensional manifolds [62, 36, 37], such higher-rank
versions of electrodynamics have become a topical research field. For further
details about the mapping between a spiral spin liquid and a rank-2
electrodynamics theory, we refer the interested reader to Ref. [35].

A characteristic feature of a rank-2 electrodynamics theory is a four-fold
pinch point (4FPP) in the electric field correlation function defined by

$\mathcal{C}_{EE}(\mathbf{q})=\frac{1}{N}\sum_{\mathbf{r},\mathbf{r^{\prime}}}\langle
E_{xx}(\mathbf{r})E_{yy}(\mathbf{r^{\prime}})\rangle\
e^{i\mathbf{q}\cdot(\mathbf{r^{\prime}}-\mathbf{r})}.$ (15)

For an exact fulfillment of Eq. (14), the 4FPP in the electric field
correlator has the form

$\mathcal{C}_{EE}(\mathbf{q})\propto\frac{q_{x}^{2}q_{y}^{2}}{|\mathbf{q}|^{4}}\propto\sin^{2}(2\theta)$
(16)

at small $q$. In the right-most expression of Eq. (16) we have used a polar
representation of the momentum, i.e., $q_{x}=q\cos(\theta)$ and
$q_{y}=q\sin(\theta)$. This expression makes it obvious that an ideal 4FFP
does not have a dependence on $q$. Importantly, the sharpness of the non-
analyticity of Eq. (16) at $\mathbf{q}=0$ serves as a useful measure for the
fulfillment of the generalized Gauss law in Eq. (14). In Ref. [35], well-
defined 4FPP in $\mathcal{C}_{EE}(\mathbf{q})$ have already been observed in
the spiral spin liquid phase at $\delta=0.03$ using numerically-obtained spin
configurations.

The mapping between a spiral spin liquid and a rank-2 electrodynamics theory,
however, depends on various assumptions such as small $\delta\ll 1$ and the
absence of spin vortices, both of which are not fulfilled in large parts of
the phase diagram. Furthermore, the exact 4FPP shape of the electric field
correlator in Eq. (16) does not take into account possible thermalization
problems, which, however, are known to occur in our spiral model and which
limit the configuration space the system can explore. Therefore, it is a
priori unclear to what extent the analogy between the spiral spin liquid and
the rank-2 electrodynamics theory is valid across the phase diagram. To answer
this question, we calculate the correlation function
$\mathcal{C}_{EE}(\mathbf{q})$ in Eq. (15) from spin configurations obtained
by cMC in the entire $\delta$-$T$ phase diagram and investigate the intactness
of 4FPPs.

Figure 15: Electric field correlator $\mathcal{C}_{EE}(\mathbf{q})$ [see Eq.
(15)]: (a) for $\delta=0.03$, $T=0.081$ and (b) for $\delta=0.10$, $T=0.119$.
(c) and (d) show the corresponding signal along circular paths
$\theta\in[0,2\pi]$ with radii $q$ ranging from 0.1 to 1.0 at 0.1 intervals
[indicated by color lines in the first quadrants of (a) and (b)].

Note that even for an exact lattice realization of a rank-2 electrodynamics
theory Eq. (16) is only valid at small $q$ while near the boundary of the
Brillouin zone, lattice effects play a role that leads to deviations from the
exact 4FPP shape. For a large part of the phase diagram, we observe four
symmetric lobes in $\mathcal{C}_{EE}(\mathbf{q})$ resembling the shape of
4FPPs at intermediate values of $q$. The decisive property, however, is the
behavior of these lobes at $q\rightarrow 0$.

We observe distinctly different behaviors in this limit across the
$\delta$-$T$ phase diagram, for which we show two representative examples in
Fig. 15. Panels (a) and (c) display $\mathcal{C}_{EE}(\mathbf{q})$ from cMC
results at $\delta=0.03$ and $T=0.081$ (close to the phase transition between
the pancake and spiral spin liquids). Ideal pinch-point singularities have an
intensity that only depends on the angle $\theta$ in momentum space and not on
the distance $q$ from the origin. To check this, we plot in Fig. 15 (c) the
signal along the circular paths $\theta\in[0,2\pi]$ with constant $q$ as
indicated in Fig. 15 (a). We observe that all curves overlap for different
radii $q$ down to $q=0.1$, very close to the origin, indicating an intact
4FPP. On the other hand, we show in Fig. 15 (b) and (d) results for
$\delta=0.10$ and $T=0.119$ (close to the phase transition between the
paramagnet and the nematic spiral). In this case, we see that the intensity
vanishes as the radius $q$ decreases, signaling the absence of a singularity.
Such a smearing of pinch points may indicate that emergent charge fluctuations
are soft and become thermally activated [6].

To investigate the behavior of 4FPP at small momenta more systematically and
over the whole phase diagram, we define a quality measure of pinch points in
the following way. For the 10 different values of the radii $q_{i}$ ranging
from 0.1 to 1.0 in steps of 0.1, we fit the data to
$A(q_{i})\sin^{2}(2\theta)$ and extract $A(q_{i})$. We then define our quality
measure as $Q=b^{2}/|a|$ where $a$ and $b$ are obtained from a fit to a linear
function $A(q)=aq+b$. The choice of defining $Q$ with $b^{2}$ in the numerator
as opposed to just $b$ is to suppress the influence of a small fluctuating
background signal for $a,b\sim 0$. Then, $Q$ takes large values when $A(q)$ is
approximately constant ($a\to 0$), signaling a sharp singularity. On the other
hand, $Q$ takes negligible values if the slope $a$ is finite and/or $b$
vanishes (featureless electric field correlator), indicating deviations from
the exact pinch point shape in Eq. (16). The results for $Q$ as a function of
$\delta$ and $T$ are shown in Fig. 16, where it becomes clear that the quality
of 4FPPs is highest around and below the phase transition into the spiral spin
liquid (brown line). This demonstrates that the fluctuations in the spiral
spin liquid can indeed be described by an effective Gauss law associated with
a higher-rank generalization of electromagnetism. In the nematic spiral phase
(below the green line), four-fold symmetric lobes with
$\mathcal{C}_{EE}(\mathbf{q})\sim A\sin^{2}(2\theta)$ are also observed, as
shown in Fig. 15. However, in this case, the magnitude $A=A(q)$ depends on the
distance $q$ from the origin.

Figure 16: Pinch-point singularity quality $Q$ as a function of $T$ and
$\delta$ (see main text). Lines and symbols correspond to the phase
transitions indicated by the heat capacity.

##  VII Conclusions

We have used classical Monte Carlo and molecular dynamics simulations to study
the finite-temperature phase diagram of a prototypical spiral spin model on
the square lattice with XY spins. We have identified a variety of interesting
phases and emergent phenomena in this model. First investigating the regime of
large $\delta$, at $\delta=0.25$ we found that the ground-state manifold
contains spiral stripes with a wavelength of four sites, as well as colinear
stripes with two-up-two-down magnetic unit cells, called the $\langle
2\rangle$ order. At finite temperatures, entropic effects select ordering wave
vectors along the $x$ or $y$ directions corresponding to a discrete
$\mathds{Z}_{2}$ symmetry breaking of the $C_{4}$ lattice rotation symmetry
down to $C_{2}$. On top of this, entropy selects the colinear $\langle
2\rangle$ states over the homogeneous spin spirals which constitutes an
additional $\mathds{Z}_{2}$ symmetry breaking related to lattice translation.
We found that the two $\mathds{Z}_{2}$ symmetries are broken at the same
temperature due to order-by-disorder effects, leading to a
$\mathds{Z}_{2}\times\mathds{Z}_{2}=\mathds{Z}_{4}$ symmetry breaking into the
$\langle 2\rangle$ phase. We verified that the critical exponents agree with a
second-order Ashkin-Teller (or four-state Potts) phase transition.

We found that the phase transition into the $\langle 2\rangle$ phase persists
for $0.16\leq\delta\leq 0.25$, even though the $\langle 2\rangle$ order is no
longer part of the ground-state manifold. As a consequence, a second phase
transition into a nematic spiral state arises at lower temperatures. This is
an example of a re-entrance phenomenon since it implies a restoration of the
broken lattice translation symmetry (while lattice rotation symmetry remains
broken). We verified that this second phase transition belongs to the Ising
universality class. For $0.05<\delta\leq 0.15$, the system shows only one
finite-temperature phase transition into the nematic spiral state which breaks
lattice rotation symmetry. This transition is again associated with a broken
$\mathds{Z}_{2}$ symmetry and the transition belongs to the Ising universality
class.

For $\delta<0.05$, the peak in the specific heat splits into two peaks. A wide
shoulder shifts to higher temperatures as $\delta$ decreases, while a sharp
low-temperature peak decreases in intensity and shifts towards $T\rightarrow
0$ as $\delta\rightarrow 0$. The high-temperature feature shows an off-
critical behavior, indicating a crossover into a pancake spin liquid, where
spirals with all wave vectors inside the spiral ring coexist. On the other
hand, the low-temperature transition leads into a spiral spin liquid phase and
shows a logarithmic divergence of the specific heat in system size, consistent
with a critical exponent $\alpha=0$ and reminiscent of an Ising transition.
However, no spontaneous symmetry breaking is observed across this transition
and the critical exponent $\nu$ is inconsistent with an Ising transition,
excluding a standard second-order transition. Instead, we found indications
that the transition is driven by momentum vortices. While the density of
momentum vortices shows a sudden decrease at the transition, the momentum-
momentum correlation length increases sharply, pointing to the possibility of
an unusual and hitherto unexplored Kosterlitz-Thouless transition from
momentum degrees of freedom.

We also performed molecular dynamics calculations to characterize the
dynamical spin structure factor of the different phases. The colinear $\langle
2\rangle$ phase shows low-energy spin-wave bands at all momenta on the spiral
ring, contrary to the nematic spiral phase where all modes are gapped, except
for the Goldstone mode. On the other hand, the pancake liquid shows abundant
low-energy excitations for all wave vectors inside the spiral ring with a
continuum of excitations extending to higher energies. For the spiral spin
liquid, low-energy modes are only observed in the vicinity of the spiral ring
while at higher energies a fading of the continuum and spectral weight mostly
concentrated along well-defined spin wave bands is found.

Finally, we verified the claim of Ref. [35] that the fluctuations in the
spiral spin liquid can be captured by an effective rank-2 Gauss law for an
emergent matrix-valued electric field as it appears in a so-called rank-2 U(1)
electrodynamics theory. We took the sharpness of the characteristic 4FPP
singularities in the electric-field correlator as a measure to determine the
regions in the phase diagram where the mapping onto a rank-2 electrodynamics
theory holds. We, indeed, verified that the sharpest 4FPP are observed in the
spiral spin liquid regime, demonstrating the close connection between this
phase and an emergent rank-2 electrodynamics.

In total, our work presents a coherent picture of a spin model that, even on
the classical level, displays a wealth of fascinating features ranging from
spiral spin liquids, re-entrance phenomena, and unusual Kosterlitz-Thouless
transitions to emergent higher-rank gauge theories. While our work resolves
many unanswered questions from previous studies, it also points to several
aspects that deserve further investigation. For example, developing a deeper
understanding and verification of the proposed Kosterlitz-Thouless transition
from momentum degrees of freedom would be a worthwhile future research effort.
Furthermore, the inclusion of quantum fluctuations and their effects on the
rank-2 electrodynamics theory may give rise to even more fascinating emergent
phenomena in this spin system.

## Acknowledgements

The authors would like to thank Johannes Knolle, Piet Brouwer, Achim Rosch,
Matthias Vojta, Randall Samien and Leo Radzihovsky for inspiring discussions.
M.G.G. and A.F. would like to thank the CURTA cluster of ZEDAT, Freie
Universität Berlin, and the JUWELS cluster at the Forschungszentrum Jülich for
computing time.

##  Appendix A Exact vs. cMC ground-state energy

The ground-state energy can be obtained via cMC as a continuation of the
cooldown process to $T=0$. In Fig. 17 we compare these results for $L=200$ and
periodical boundary conditions with the exact values of the ground-state
energy for the spiral solutions. At low values of $\delta$, we can see that
the energy of the ferromagnetic state (green line) lies very close above the
spiral state. This can explain why at finite temperatures, when these
ferromagnetic states can be realized, a pancake liquid emerges that displays
strong ferromagnetic and spiral correlations up to the edge of the spiral ring
in reciprocal space.

On the other hand, due to the frustrating boundary conditions for
incommensurate spiral phases, it is expected that the exact ground-state
energy differs from the cMC result (bottom panel). However, we can see that
the difference is small. Two parts show larger differences: One is the region
where the system passes through the $\langle 2\rangle$ phase. As explained in
the main text, there are reminiscences of the passage through this phase in
the form of small but finite $\mathcal{O}_{1}$ order parameter values [see Eq.
(6)] that lead to a higher ground-state energy in cMC. On the other hand, also
when the system goes through the vortex network, the cMC ground-state energy
differs from the exact one, showing that the dynamics of these structures are
slow and that the systems need more space to accommodate properly.

Figure 17: Top panel: Ground-state energy per site $e_{0}$ obtained from cMC
(blue circles) compared to the exact energies for the spiral ground state
(orange line) and the ferromagnetic state (green line). Bottom panel: Energy
difference $e_{0}^{\text{MC}}-e_{0}$ between the cMC result and the exact
value.

##  Appendix B Finite-size scaling and universality classes of phase
transitions

To determine the order of the different phase transitions and the
corresponding universality classes, we perform cMC calculations on a large
range of lattice sizes starting from $L=4$ up to $L=80$. In some cases, the
peaks in the specific heat $c_{v}(T)$ get sharp (thin and high) very fast with
increasing system size. In such cases, it becomes difficult to obtain an
accurate value of the peak, $c_{v}^{\text{max}}$, while $T_{c}$ is well-
defined.

The correlation length of the system near the critical point diverges like
$\xi\propto|T-T_{c}|^{-\nu}$, where $\nu$ is the critical exponent. However,
in finite systems, the correlation length cannot be larger than $\sim L$.
Therefore, we can assume that $|T_{c}(\infty)-T_{c}(L)|^{-\nu}\propto L$ where
$T_{c}(L)$ is the position of the peak in the specific heat for a given linear
size $L$ and $T_{c}(\infty)$ is the value in the thermodynamic limit
$L\to\infty$. Then we obtain the finite-size scaling law [63, 64, 65]

$T_{c}(L)=T_{c}(\infty)+a\ L^{-1/\nu}.$ (17)

Using this, we can fit the positions of the peak in $c_{v}(T)$ as a function
of $1/L$ and obtain $T_{c}(\infty)$, $a$, and $\nu$ from the fit. While
typically $T_{c}(\infty)$ is well defined by these fits, $\nu$ tends to
present larger uncertainties because the curvature of $T_{c}(L)$ is affected
by the values on small lattices, where the effect of the boundary conditions
on the state of the systems is more visible. On the other hand, the specific
heat behaves like $c_{v}\propto|T-T_{c}|^{-\alpha}$ around the critical point,
such that it is straightforward to see that [65]

$c_{v}^{\text{max}}(L)\propto L^{\alpha/\nu}.$ (18)

In this case, a log-log plot allows the extraction of the quotient
$\alpha/\nu$ as the slope of a linear fit, whereas in the case of Eq. (17) an
unknown constant $T_{c}(\infty)$ needs to be removed before using the same
method. For this reason, the determination of $\alpha/\nu$ is more precise
than that of $\nu$ from $T_{c}(L)$. However, obtaining a precise value for
$c_{v}^{\text{max}}(L)$ is complicated for large lattices because the peak
becomes too thin and high, leading to large differences in
$c_{v}^{\text{max}}(L)$ for small deviations in $T_{c}(L)$.

In practice, there are two different routes for calculating the specific heat
and obtaining $c_{v}^{\text{max}}(L)$ and $T_{c}(L)$. One consists of
calculating $c_{v}(T)=N(\langle e(T)^{2}\rangle-\langle
e(T)\rangle^{2})/T^{2}$, where $e(T)$ is the internal energy per site and
$\langle.\rangle$ is the cMC average over different Monte Carlo steps at a
given temperature. Then, $c_{v}(T)$ is averaged over several independent runs.
The second route consists of taking the discrete derivative of $e(T)$ in the
temperature grid available. This is done after averaging $e(T)$ over
independent runs. Ideally, both approaches should lead to the same result, but
in practice, we find that the second option leads to less noisy results.
Still, the specific heat calculations are not always smooth for small
temperature steps $\Delta T$ because small uncertainties in the energy $e$
lead to large ones when calculating $\Delta e/\Delta T$. To overcome this, we
use a Savitzky–Golay filter to smooth the curves and improve the finite-size
scaling [66].

Figure 18: Finite-size scaling of $c_{v}^{\text{max}}$ and $T_{c}$ for the
phase transition at $\delta=0.25$ between the paramagnetic phase and the
$\langle 2\rangle$-ordered phase. Points correspond to the cMC calculations,
while fits to extract the critical exponents are shown with lines (see main
text).

In Fig. 18, we show the results for the high-temperature phase transition at
$\delta=0.25$, where the system goes from a paramagnetic phase to a phase with
$\langle 2\rangle$ order. We obtain $\alpha/\nu\approx 1.0$, which is
consistent with a 4-state Potts (or Ashkin-Teller) transition in two
dimensions (top panel). These four possible states are realized by
combinations of horizontal or vertical stripe configurations and the two
possibilities for the two patterns of stripes
($\phi\phi\overline{\phi}\overline{\phi}$ and
$\phi\overline{\phi}\overline{\phi}\phi$) related by a lattice translation of
one nearest neighbor distance. In the bottom panel of Fig. 18, we see that
$\nu$ takes values close to the expected $\nu=2/3$ for the 4-state Potts
transition.

Figure 19: Finite-size scaling of $c_{v}^{\text{max}}$ and $T_{c}$ for the
phase transition at $\delta=0.17$ between the $\langle 2\rangle$ phase and the
nematic spiral phase. Note that the top panel uses a logarithmic scale only on
the $x$-axis but not on the $y$-axis because $\alpha=0$ for an Ising
transition and the leading term becomes $c_{v}^{\text{max}}(L)\propto\log(L)$.
The parameter that indicates the goodness of the fit, $R^{2}$, is shown in the
legend and is close to the optimal value of 1.

In Fig. 19, we show the results for the low-temperature phase transition
between the $\langle 2\rangle$ phase and the nematic spiral phase at
$\delta=0.17$. In this case, the peaks in the specific heat become very thin
as the system size increases, and it becomes difficult to determine the height
accurately (top panel). However, the critical temperature shown in the bottom
panel is very well fitted by $\nu\approx 1.1$. This is consistent with $\nu=1$
expected for an Ising transition in two dimensions. We can think of this
transition coming from the nematic spiral phase at lower temperatures. Then,
when the temperature increases and the system evolves to the $\langle
2\rangle$ phase, there are two possible choices of double stripes
($\phi\phi\overline{\phi}\overline{\phi}$ and
$\phi\overline{\phi}\overline{\phi}\phi$). An Ising transition is
characterized by $\alpha=0$ and the next relevant term in Eq. (18) is
$\propto\log(L)$. The top panel of Fig. 19 shows that the numerical result for
$c_{v}^{\text{max}}(L)$ as a function of $\log(L)$ is indeed well fitted by a
linear function.

Figure 20: Finite-size scaling of $c_{v}^{\text{max}}$ and $T_{c}$ for
$\delta=0.10$, extracted from the transition between the disordered and the
nematic spiral state. Note that the top panel does not show a log-log plot but
a lin-log plot to obtain the leading term
$c_{v}^{\text{max}}(L)\propto\log(L)$ for an Ising transition where
$\alpha=0$.

In Fig. 20 we show the results for $\delta=0.10$, where there is a transition
between a disordered and a nematic spiral state. As in the previous case, the
transition is expected to be of Ising type because of the selection of stripes
along the $x$ or $y$ directions. This is confirmed by the linear growth of the
peak in the specific heat $c_{v}^{\text{max}}(L)$ as a function of $\log(L)$
(top panel). On the other hand, $T_{c}(L)$ shows an approximate linear
behavior as a function of $1/L$, indicating $\nu=1$, as expected for an Ising
transition. The large deviation in the values $T_{c}(L)$ for small systems can
be associated with the frustration induced by the boundary conditions, which
is stronger for incommensurate spirals in the case of small lattices.

Figure 21: Finite-size scaling of $c_{v}^{\text{max}}$ and $T_{c}$ for the
low-temperature peak at $\delta=0.03$, where the system evolves from the
pancake to the spiral spin liquid. The orange square corresponds to the result
of Ref. [35] for $L=400$.

Finally, in Fig. 21 we show the finite-size scaling for the low-temperature
peak at $\delta=0.03$, where the system goes from the pancake liquid to the
spiral spin liquid. In this case, we observe that $c_{v}^{\text{max}}(L)$
diverges logarithmically with system size, indicating that the corresponding
critical exponent vanishes, $\alpha=0$. However, $T_{c}(L)$ does not evolve
linearly with $1/L$, as expected for an Ising transition ($\alpha=0$).
Instead, the critical exponent is closer to $\nu=1/2$, which corresponds to a
linear behavior in $L^{-1/\nu}=N^{-1}$. As explained in the main text, this
transition is connected to the proliferation of momentum vortices and
antivortices and may be a Kosterlitz-Thouless transition of momentum vortices.
However, a more thorough investigation is needed to confirm this scenario.

##  Appendix C Transition between the spiral spin liquid and nematic spiral
phase

Figure 22: Specific heat $c_{v}$ as a function of $\delta$ and $T$ in a
smaller region of the $\delta$-$T$ space than in Fig. 3, for lattice size
$L=200$ and periodical boundary conditions.

As discussed in the previous section and in the main text, the phase
transition for $\delta>0.05$ belongs to the Ising universality class and leads
to a nematic stripe spiral order. On the other hand, for $\delta<0.05$ the
transition does not fit in any universality class and presents no symmetry
breaking. Since the two low-temperature phases are fundamentally different in
terms of broken symmetries, a phase transition has to exist between them.
However, as shown in Fig. 3 and in Fig. 22, such a phase transition is not
observed in the specific heat. For the enlarged view in Fig. 22, we performed
independent cMC runs for $\delta$ between 0.03 and 0.07 in 0.001 steps. At
$\delta=0.07$, there is a clear transition from a disordered state into a
nematic spiral state. As $\delta$ is lowered, the transition is weakened,
i.e., the peak decreases and shifts to lower temperatures. Furthermore, a wide
shoulder emerges at higher temperatures (as evidenced by the contour lines).
The latter indicates the onset of correlations leading to the pancake liquid
phase. The low-energy transition indicates the passage into the spiral spin
liquid without symmetry breaking. As shown in the previous section, all along
the low-temperature phase transition, the specific heat diverges
logarithmically. The fact that the specific heat calculations cannot resolve
this transition even with this small $\delta$ step size indicates that the
single-spin update algorithm might not be sufficient to capture it and further
studies are needed to confirm its presence. Nonetheless, as discussed in the
main text, the transition can be detected via the calculation of appropriate
order parameters and evidenced by energy-level crossings in $\delta$-sweeps.

## References

  * Balents [2010] L. Balents, Spin liquids in frustrated magnets, Nature 464, 199 (2010).
  * Savary and Balents [2016] L. Savary and L. Balents, Quantum spin liquids: a review, Reports on Progress in Physics 80, 016502 (2016).
  * Broholm  _et al._ [2020] C. Broholm, R. J. Cava, S. A. Kivelson, D. G. Nocera, M. R. Norman, and T. Senthil, Quantum spin liquids, Science 367, eaay0668 (2020).
  * Bramwell and Gingras [2001] S. T. Bramwell and M. J. P. Gingras, Spin Ice State in Frustrated Magnetic Pyrochlore Materials, Science 294, 1495 (2001).
  * Castelnovo  _et al._ [2008] C. Castelnovo, R. Moessner, and S. L. Sondhi, Magnetic monopoles in spin ice, Nature 451, 42 (2008).
  * Fennell  _et al._ [2009] T. Fennell, P. P. Deen, A. R. Wildes, K. Schmalzl, D. Prabhakaran, A. T. Boothroyd, R. J. Aldus, D. F. McMorrow, and S. T. Bramwell, Magnetic Coulomb Phase in the Spin Ice Ho2Ti2O7, Science 326, 415 (2009).
  * Morris  _et al._ [2009] D. J. P. Morris, D. A. Tennant, S. A. Grigera, B. Klemke, C. Castelnovo, R. Moessner, C. Czternasty, M. Meissner, K. C. Rule, J.-U. Hoffmann, K. Kiefer, S. Gerischer, D. Slobinsky, and R. S. Perry, Dirac Strings and Magnetic Monopoles in the Spin Ice Dy2Ti2O7, Science 326, 411 (2009).
  * Alzate-Cardona  _et al._ [2019] J. D. Alzate-Cardona, D. Sabogal-Suárez, R. F. L. Evans, and E. Restrepo-Parra, Optimal phase space sampling for Monte Carlo simulations of Heisenberg spin systems, J. Phys. Condens. Matter 31, 095802 (2019).
  * Chalker  _et al._ [1992] J. T. Chalker, P. C. W. Holdsworth, and E. F. Shender, Hidden order in a frustrated system: Properties of the Heisenberg Kagomé antiferromagnet, Phys. Rev. Lett. 68, 855 (1992).
  * Huse and Rutenberg [1992] D. A. Huse and A. D. Rutenberg, Classical antiferromagnets on the Kagomé lattice, Phys. Rev. B 45, 7536 (1992).
  * Isakov  _et al._ [2004] S. V. Isakov, K. Gregor, R. Moessner, and S. L. Sondhi, Dipolar Spin Correlations in Classical Pyrochlore Magnets, Phys. Rev. Lett. 93, 167204 (2004).
  * Moessner and Chalker [1998a] R. Moessner and J. T. Chalker, Properties of a Classical Spin Liquid: The Heisenberg Pyrochlore Antiferromagnet, Phys. Rev. Lett. 80, 2929 (1998a).
  * Hopkinson  _et al._ [2007] J. M. Hopkinson, S. V. Isakov, H.-Y. Kee, and Y. B. Kim, Classical Antiferromagnet on a Hyperkagome Lattice, Phys. Rev. Lett. 99, 037201 (2007).
  * Davier  _et al._ [2023] N. Davier, F. A. Gómez Albarracín, H. D. Rosales, and P. Pujol, Combined approach to analyze and classify families of classical spin liquids, Phys. Rev. B 108, 054408 (2023).
  * Yan  _et al._ [2023a] H. Yan, O. Benton, R. Moessner, and A. H. Nevidomskyy, Classification of Classical Spin Liquids: Typology and Resulting Landscape (2023a), arXiv:2305.00155 [cond-mat.str-el] .
  * Yan  _et al._ [2023b] H. Yan, O. Benton, A. H. Nevidomskyy, and R. Moessner, Classification of Classical Spin Liquids: Detailed Formalism and Suite of Examples (2023b), arXiv:2305.19189 [cond-mat.str-el] .
  * Fang  _et al._ [2023] Y. Fang, J. Cano, A. H. Nevidomskyy, and H. Yan, Classification of Classical Spin Liquids: Topological Quantum Chemistry and Crystalline Symmetry (2023), arXiv:2309.12652 [cond-mat.str-el] .
  * Sachdev [1992] S. Sachdev, Kagomé- and triangular-lattice Heisenberg antiferromagnets: Ordering from quantum fluctuations and quantum-disordered ground states with unconfined bosonic spinons, Phys. Rev. B 45, 12377 (1992).
  * Lecheminant  _et al._ [1997] P. Lecheminant, B. Bernu, C. Lhuillier, L. Pierre, and P. Sindzingre, Order versus disorder in the quantum Heisenberg antiferromagnet on the kagomé lattice using exact spectra analysis, Phys. Rev. B 56, 2521 (1997).
  * Yan  _et al._ [2011] S. Yan, D. A. Huse, and S. R. White, Spin-Liquid Ground State of the S = 1/2 Kagome Heisenberg Antiferromagnet, Science 332, 1173 (2011).
  * Messio  _et al._ [2012] L. Messio, B. Bernu, and C. Lhuillier, Kagome Antiferromagnet: A Chiral Topological Spin Liquid?, Phys. Rev. Lett. 108, 207204 (2012).
  * Gingras and McClarty [2014] M. J. P. Gingras and P. A. McClarty, Quantum spin ice: a search for gapless quantum spin liquids in pyrochlore magnets, Reports on Progress in Physics 77, 056501 (2014).
  * Huse  _et al._ [2003] D. A. Huse, W. Krauth, R. Moessner, and S. L. Sondhi, Coulomb and Liquid Dimer Models in Three Dimensions, Phys. Rev. Lett. 91, 167004 (2003).
  * Hermele  _et al._ [2004] M. Hermele, M. P. A. Fisher, and L. Balents, Pyrochlore photons: The $U(1)$ spin liquid in a $S=\frac{1}{2}$ three-dimensional frustrated magnet, Phys. Rev. B 69, 064404 (2004).
  * Bergman  _et al._ [2007] D. Bergman, J. Alicea, E. Gull, S. Trebst, and L. Balents, Order-by-disorder and spiral spin-liquid in frustrated diamond-lattice antiferromagnets, Nature Physics 3, 487 (2007).
  * Mulder  _et al._ [2010] A. Mulder, R. Ganesh, L. Capriotti, and A. Paramekanti, Spiral order by disorder and lattice nematic order in a frustrated Heisenberg antiferromagnet on the honeycomb lattice, Phys. Rev. B 81, 214419 (2010).
  * Gao  _et al._ [2017] S. Gao, O. Zaharko, V. Tsurkan, Y. Su, J. S. White, G. Tucker, B. Roessli, F. Bourdarot, R. Sibille, D. Chernyshov, T. Fennell, A. Loidl, and C. Rüegg, Spiral spin-liquid and the emergence of a vortex-like state in MnSc2S4, Nature Physics 13, 157 (2017).
  * Iqbal  _et al._ [2018] Y. Iqbal, T. Müller, H. O. Jeschke, R. Thomale, and J. Reuther, Stability of the spiral spin liquid in ${\mathrm{MnSc}}_{2}{\mathrm{S}}_{4}$, Phys. Rev. B 98, 064427 (2018).
  * Buessen  _et al._ [2018] F. L. Buessen, M. Hering, J. Reuther, and S. Trebst, Quantum Spin Liquids in Frustrated Spin-1 Diamond Antiferromagnets, Phys. Rev. Lett. 120, 057201 (2018).
  * Ghosh  _et al._ [2019] P. Ghosh, Y. Iqbal, T. Müller, R. T. Ponnaganti, R. Thomale, R. Narayanan, J. Reuther, M. J. P. Gingras, and H. O. Jeschke, Breathing chromium spinels: a showcase for a variety of pyrochlore heisenberg hamiltonians, npj Quantum Materials 4, 63 (2019).
  * Niggemann  _et al._ [2019] N. Niggemann, M. Hering, and J. Reuther, Classical spiral spin liquids as a possible route to quantum spin liquids, Journal of Physics: Condensed Matter 32, 024001 (2019).
  * Yao  _et al._ [2021] X.-P. Yao, J. Q. Liu, C.-J. Huang, X. Wang, and G. Chen, Generic spiral spin liquids, Frontiers of Physics 16, 53303 (2021).
  * Gao  _et al._ [2022] S. Gao, M. A. McGuire, Y. Liu, D. L. Abernathy, C. d. Cruz, M. Frontzek, M. B. Stone, and A. D. Christianson, Spiral Spin Liquid on a Honeycomb Lattice, Phys. Rev. Lett. 128, 227201 (2022).
  * Graham  _et al._ [2023] J. N. Graham, N. Qureshi, C. Ritter, P. Manuel, A. R. Wildes, and L. Clark, Experimental Evidence for the Spiral Spin Liquid in ${\mathrm{LiYbO}}_{2}$, Phys. Rev. Lett. 130, 166703 (2023).
  * Yan and Reuther [2022] H. Yan and J. Reuther, Low-energy structure of spiral spin liquids, Phys. Rev. Res. 4, 023175 (2022).
  * Pretko [2017a] M. Pretko, Generalized electromagnetism of subdimensional particles: A spin liquid story, Phys. Rev. B 96, 035119 (2017a).
  * Pretko [2017b] M. Pretko, Subdimensional particle structure of higher rank $U(1)$ spin liquids, Phys. Rev. B 95, 115139 (2017b).
  * Prem  _et al._ [2018] A. Prem, S. Vijay, Y.-Z. Chou, M. Pretko, and R. M. Nandkishore, Pinch point singularities of tensor spin liquids, Phys. Rev. B 98, 165140 (2018).
  * Hsieh  _et al._ [2022] T.-C. Hsieh, H. Ma, and L. Radzihovsky, Helical superfluid in a frustrated honeycomb Bose-Hubbard model, Phys. Rev. A 106, 023321 (2022).
  * Zhai and Radzihovsky [2021] Z. Zhai and L. Radzihovsky, Fractonic gauge theory of smectics, Annals of Physics 435, 168509 (2021), special issue on Philip W. Anderson.
  * Radzihovsky [2020] L. Radzihovsky, Quantum Smectic Gauge Theory, Phys. Rev. Lett. 125, 267601 (2020).
  * Radzihovsky [2023] L. Radzihovsky, Critical Matter (2023), arXiv:2306.03142 [cond-mat.stat-mech] .
  * Machon  _et al._ [2019] T. Machon, H. Aharoni, Y. Hu, and R. D. Kamien, Aspects of Defect Topology in Smectic Liquid Crystals, Communications in Mathematical Physics 372, 525–542 (2019).
  * Shimokawa and Kawamura [2019] T. Shimokawa and H. Kawamura, Ripple State in the Frustrated Honeycomb-Lattice Antiferromagnet, Phys. Rev. Lett. 123, 057202 (2019).
  * Fisher and Selke [1980] M. E. Fisher and W. Selke, Infinitely Many Commensurate Phases in a Simple Ising Model, Phys. Rev. Lett. 44, 1502 (1980).
  * Seifert and Vojta [2019] U. F. P. Seifert and M. Vojta, Theory of partial quantum disorder in the stuffed honeycomb Heisenberg antiferromagnet, Phys. Rev. B 99, 155156 (2019).
  * Blesio  _et al._ [2023] G. G. Blesio, F. T. Lisandrini, and M. G. Gonzalez, Partially disordered Heisenberg antiferromagnet with short-range stripe correlations, Phys. Rev. B 107, 134418 (2023).
  * Chandra  _et al._ [1990] P. Chandra, P. Coleman, and A. I. Larkin, Ising transition in frustrated Heisenberg models, Phys. Rev. Lett. 64, 88 (1990).
  * Weber  _et al._ [2003] C. Weber, L. Capriotti, G. Misguich, F. Becca, M. Elhajal, and F. Mila, Ising Transition Driven by Frustration in a 2D Classical Model with Continuous Symmetry, Phys. Rev. Lett. 91, 177202 (2003).
  * Wu [1982] F. Y. Wu, The Potts model, Rev. Mod. Phys. 54, 235 (1982).
  * Grison  _et al._ [2020] V. Grison, P. Viot, B. Bernu, and L. Messio, Emergent Potts order in the kagome ${J}_{1}-{J}_{3}$ Heisenberg model, Phys. Rev. B 102, 214424 (2020).
  * Moessner and Chalker [1998b] R. Moessner and J. T. Chalker, Low-temperature properties of classical geometrically frustrated antiferromagnets, Phys. Rev. B 58, 12049 (1998b).
  * Conlon and Chalker [2009] P. H. Conlon and J. T. Chalker, Spin Dynamics in Pyrochlore Heisenberg Antiferromagnets, Phys. Rev. Lett. 102, 237206 (2009).
  * Zhang  _et al._ [2019] S. Zhang, H. J. Changlani, K. W. Plumb, O. Tchernyshyov, and R. Moessner, Dynamical Structure Factor of the Three-Dimensional Quantum Spin Liquid Candidate ${\mathrm{NaCaNi}}_{2}{\mathrm{F}}_{7}$, Phys. Rev. Lett. 122, 167203 (2019).
  * Nho and Landau [2002] K. Nho and D. P. Landau, Spin-dynamics simulations of the triangular antiferromagnetic XY model, Phys. Rev. B 66, 174403 (2002).
  * Toth and Lake [2015] S. Toth and B. Lake, Linear spin wave theory for single-Q incommensurate magnetic structures, Journal of Physics: Condensed Matter 27, 166002 (2015).
  * Kosterlitz and Thouless [1973] J. M. Kosterlitz and D. J. Thouless, Ordering, metastability and phase transitions in two-dimensional systems, Journal of Physics C: Solid State Physics 6, 1181 (1973).
  * Kosterlitz [1974] J. M. Kosterlitz, The critical properties of the two-dimensional xy model, Journal of Physics C: Solid State Physics 7, 1046 (1974).
  * Ota  _et al._ [1992] S. Ota, S. B. Ota, and M. Fahnle, Microcanonical Monte Carlo simulations for the two-dimensional XY model, Journal of Physics: Condensed Matter 4, 5411 (1992).
  * Lee  _et al._ [2005] K. W. Lee, C. E. Lee, and I. mook Kim, Helicity modulus and vortex density in the two-dimensional easy-plane Heisenberg model on a square lattice, Solid State Communications 135, 95 (2005).
  * Nguyen and Boninsegni [2021] P. H. Nguyen and M. Boninsegni, Superfluid Transition and Specific Heat of the 2D x-y Model: Monte Carlo Simulation, Applied Sciences 11, 4931 (2021).
  * Ma and Pretko [2018] H. Ma and M. Pretko, Higher-rank deconfined quantum criticality at the lifshitz transition and the exciton bose condensate, Phys. Rev. B 98, 125105 (2018).
  * Ferdinand and Fisher [1969] A. E. Ferdinand and M. E. Fisher, Bounded and Inhomogeneous Ising Models. I. Specific-Heat Anomaly of a Finite Lattice, Phys. Rev. 185, 832 (1969).
  * Landau [1976] D. P. Landau, Finite-size behavior of the Ising square lattice, Phys. Rev. B 13, 2997 (1976).
  * Privman [1990] V. Privman,  _Finite Size Scaling and Numerical Simulation of Statistical Systems_ (WORLD SCIENTIFIC, 1990).
  * Savitzky and Golay [1964] A. Savitzky and M. J. E. Golay, Smoothing and Differentiation of Data by Simplified Least Squares Procedures., Analytical Chemistry 36, 1627 (1964).