Low energy structure of spiral spin liquids

In this work we identify a previously unexplored type of topological defect in spiral spin liquids -- the momentum vortex -- and reveal its dominant role in shaping the low energy physics of such systems. Spiral spin liquids are a class of classical spin liquids featuring sub-extensively degenerate ground states. They are distinct from spin liquids on geometrically frustrated lattices, in which the ground state degeneracy is extensive and connected by local spin flips. Despite a handful of experimental realizations and many theoretical studies, a concrete physical picture of their spin liquidity has not been established so far. In this work, we study a 2D spiral spin liquid model to answer this question. We find that the local momentum vector field can carry topological defects in the form of vortices, which, however, have very different properties from the commonly known spin vortices. The fluctuations of such vortices lead the system into a liquid phase at intermediate temperatures. Furthermore, the effective low energy theory of such vortices indicates their equivalence to quadrupoles of fractons in a rank-2 U(1) gauge theory or, alternatively, to quadrupoles of disclinations in elasticity theory. At very low temperatures, the system freezes into a glassy state in which these vortices form a rigid network with straight-line domain walls. Our work sheds light on the nature of spiral spin liquids and also paves the way toward understanding their quantum limit.

A fascinating theme of condensed matter is emergence [1] -from the Anderson-Higgs mechanism in superconductors, to topological defects in superfluid 3 He, to gauge theories in the description of frustrated magnets -each condensed matter system is itself a different universe in a grain of sand. In particular, spin liquids, both classical and quantum, have been proven to be a fruitful field to search for exotic phases, (classical analogs of) fractionalized excitations and topological orders [2][3][4][5][6][7][8][9][10][11].
However, a clear physical picture of their spin liquidity has not been established to our knowledge. This is in contrast to spin liquids from geometric frustration (e.g. in pyrochlore spin ice systems), in which the ground state degeneracy is of local nature and, hence, the mechanism of fluctuations is much better understood: The system can visit different classical ground states by flipping only a few spins which can be easily accomplished by thermal fluctuations. Spiral spin liquids have a much smaller ground state degeneracy and, as a consequence, local manipulations are not sufficient to bring the system into different ground states. It is already known that in the thermodynamic limit, the spin configurations occupy all the ground state wave-vectors. Yet, since changing from one ground state spiral wave-vector to another is a global action on the spins, how exactly all the ground state wave-vectors can be populated remains unclear.
In this work, we thoroughly investigate the low energy behavior of 2D spiral spin liquids with XY spins to understand their spin liquidity. We find that, besides the commonly known spin vortices, there exists another type of topological defect that corresponds to vortices in the local momentum vector field on the coarse-grained lattice. These topological defects, dubbed local momentum vortices, have very different mathematical properties compared to spin vortices and play a crucial role in determining the low energy behavior of spiral spin liquids.
At intermediate temperatures, our effective continuum theory can be formulated in terms of the local momentum vector field, and is found to be similar to that of elasticity and scalar-charged rank-2 U(1) gauge theory [39][40][41][42][43][44][45][46]. The local momentum vector field plays the role of lattice distortions in elasticity, and their vortices are identified as vacancies/interstitials in elasticity or quadrupoles of fractons [47,48]. The liquid phase can then be understood as a state with mobile topological defects populating the system.
Proceeding to lower temperatures, the defects continuously lose their mobility, leading to a glassy state which we call a rigid network of momentum vortices. In this network state, the momentum vortices sit at the vertices, and the edges connecting them are narrow domain walls between regimes of different momenta. Particularly, due to the unusual low energy properties of momentum vortices, the domain walls must be straight lines, establishing the rigidity of the network. All these phenomena are thoroughly analysed and numerically demonstrated.
Our study of spiral spin liquids answers the key question of their low energy structure and reveals a previously unexplored type of topological defect with unexpected connections to fracton physics. It also paves the way for understanding the quantum limit of these models [14,25,49], which may give rise to quantum spin liquids. A particularly promising system to realize the phases studied here, but under the additional effects of quantum fluctuations, is the bilayer kagome material Ca 10 Cr 7 O 28 [33] which has recently been identified as a quantum spin liquid candidate.

II. THE SPIRAL SPIN LIQUID MODEL AND ITS PHASE DIAGRAM
In this work, our main focus is on 2D spiral spin liquids whose T = 0 ground state degeneracy is homotopic to a ring. Our showcase example is the classical square lattice XY spin model with couplings up to third nearest neighbours.
However, this is not the only model exhibiting spiral spin liquid physics. There are other models defined on the honeycomb and triangular lattices with further neighbour couplings [12,14,15,25] in 2D, and also on the diamond [28], face centered cubic [18,50,51], and body centered cubic [18] lattices in 3D. A more general overview on the construction of these models can be found in Refs. 24

and 25.
A. The square lattice XY spin model The Hamiltonian for the square lattice XY spin model is given by (1) Here, S i are normalized (|S i | = 1) classical XY (i.e. twocomponent) spins and J 1,2,3 are first, second, and third nearest neighbour couplings respectively, see Fig. 1(a). In the region of parameters the ground states are spin spirals of momentum q, S i = (cos(Φ(r i )), sin(Φ(r i )) = (cos(q · r i + φ), sin(q · r i + φ)) where r i is the position of site i, and φ corresponds to a global rotation of all spins [12]. Most importantly, the system exhibits a degenerate set of ground state spirals with momenta q satisfying the condition 2 cos 2 (q x ) + 2 cos 2 (q y ) + 4 cos(q x ) cos(q y ) = 1 The solutions q form a continuous 1D manifold isomorphic to a loop around the Brillouin zone center. The shape of the manifold depends on the value of J 2 , and is illustrated in Fig. 1(c). In the limit of small the spiral contour shrinks and becomes circular, obeying q = |q| = 4 √ δ, until for J 2 < 1/4 a simple ferromagnet is realized. An important property of XY models with this type of ground state degeneracy is that they exhibit two distinct types of U(1) symmetries. The first one is the standard global U(1) spin rotational symmetry which is generated by a simultaneous rotation of all spins, Φ(r i ) → Φ(r i ) + α. The second one is a U(1) symmetry in momentum space which changes the momentum q of a ground state spiral along the contour of solutions of Eq. (4). Note that this second symmetry is in principle only a property within the exact ground state manifold. However, in the limit of small δ and small momenta q where the system is approximately spherical, the effective low-energy U(1) symmetry still stands.
The consequences of U(1) symmetry in spin space are well understood. It can trigger a finite temperature Kosterlitz-Thouless transition associated with a proliferation of spin vortices [52,53]. On the other hand, the effective momentum U(1) symmetry is much less studied. While previous works, mostly focusing on numerics, have found that it gives rise to a spiral spin liquid phase where the system thermally fluctuates through the degenerate manifold of spiral states, the precise mechanism behind such fluctuations is poorly understood.
We emphasize that this mechanism causing the spin liquid physics is very different from other betterunderstood classical spin liquid models on frustrated lattices (e.g. pyrochlore spin ice [4,54]). Those models are endowed with local zero modes giving rise to an extensive ground state degeneracy. Hence, it is very intuitive that at small but finite temperatures, local spin flips enable the system to visit the degenerate manifold of ground states leading to a classical spin liquid phase. The situation in spiral spin liquids is very different. Even though the spiral ring degeneracy is sub-extensively large, a change of momentum q is still a global operation associated with the modification of a macroscopic number of spins. It is, hence, a priori unclear how thermal fluctuations may induce a liquid-like property in the system studied here.

B. Summary of main results
Before we study spiral spin liquids in detail in the next two sections, we give a brief outline of our main conclusions. Although we concentrate here on the square lattice XY model only, we expect our results to also apply to other models with a spiral degeneracy.
First, our core result is the identification of topological defects from the spiral ring U(1) symmetry (not the spin rotation symmetry). In these topological defects, the spin configuration Φ(r) (where Φ is the spins' inplane angle) varies smoothly without singularities. However, the corresponding coarse-grained momentum which takes its value only in the neighborhood of the spiral ring [ Fig. 1(c)] and varies in space, can have a non-trivial winding around a loop. As one would expect, these momentum vortices (see Fig. 2 for two vortex configuration examples) play a central role in shaping the low temperature physics. Here, we also note that the properties of such momentum vortices are very different from the commonly-known spin vortices, due to the fact that in the absence of spin vortices, the momentum field is subject to a curl-free condition. As we will elaborate later, a peculiar consequence of this restriction is that momentum vortices can only be realized with winding numbers n ≤ 1. Our numerical Monte Carlo results for small δ indicate four different temperature regimes, see phase diagram in Fig. 1(b) and exemplary spin configurations in Fig. 1(d)-(g). At large temperatures a trivial paramagnet is realized where the spins can be considered as uncorrelated [ Fig. 1(g)]. Upon cooling, the system first undergoes a crossover into a phase referred to as 'pancake liquid' in Ref. 23, see Fig. 1(f). While this regime already shows a certain degree of correlations between spins, the thermal fluctuations are too strong to restrict the momentum to the spiral contour. As a result, spiral configurations do not form, as is indicated by the absence of clear stripelike patterns in Fig. 1(f).
When decreasing the temperature further, the system undergoes a transition at T = T * into a regime with welldefined spin spirals which can be recognized as stripy configurations in Figs. 1(d) and 1(e). The investigation and characterization of the spiral regime at T < T * represents the main subject of our work. Despite the common spiral motif, the physical properties near the upper (T < ∼ T * ) and lower (T T * ) boundaries are distinctly different. We refer to the two regimes at T < T * as 'spiral spin liquid' (or alternatively as 'R2-U1/elasticity phase') and 'rigid vortex network' and discuss them separately in Secs. III and IV, respectively.
In the spiral spin liquid phase [Figs. 1(e)] the local momentum q is approximately confined along the spiral contour, however, the direction of q fluctuates strongly in real-space and Monte-Carlo time, establishing a liquidlike property. Most importantly, momentum (anti-) vortices with winding numbers n = ±1 are clearly discernible and represent the key source of fluctuations in this regime (for comparison, see the ideal examples of momentum vortices with n = ±1 in Fig. 2). The occurrence of momentum vortices prompts us to analytically investigate their precise nature and contrast them with spin vortices.
Another surprising property of the spiral spin liquid at T < ∼ T * is that its effective continuum theory for local momentum q can be mapped onto the elasticity theory of the shear modulus. In turn, elasticity theory is known to be dual to a rank-2 U(1) fracton tensor gauge theory [43,44] in which immobile and scalar fractons (i.e., charges of the gauge theory) correspond to lattice disclinations and subdimensional fracton dipoles are related to lattice dislocations. This justifies the name R2-U1/elasticity phase. We show, however, that the nature of the momentum degree of freedom q in our spiral model and, particularly, the integer-quantized spin vortices do not allow for the existence of isolated fractons and dipoles of fractons. Instead, the momentum vortices of winding number n = 1 correspond to fracton quadrupoles in R2-U1 theory, or vacancies or interstitials in elasticity. Hence, the spiral spin liquid phase can be effectively understood as a charge-and dipole-free rank-2 U(1) tensor gauge theory. Similarly, momentum vortices of winding number n < 0 can be understood as higher multipoles of vacancies and interstitials. The tensor gauge theory property of spiral spin liquids is also investigated numerically by demonstrating the occurrence of four-fold pinch points [55] in the electric-field correlator (Fig. 8).
The spiral spin liquid exhibits both spatial fluctuations in the momentum direction q/q and in the momentum amplitude q (where q remains in the vicinity of the spiral contour, q ≈ 4 √ δ). Since both variations in space cost energy, they smoothly freeze out as one further decreases the temperature. In our numerical results [see Figs. 1(d)] this freezing occurs via the formation of spiral domains each characterized by a well-defined momentum direction. In such states, the excitation energy from spatial variations of q is completely concentrated along narrow domain walls which form a rigid network spanning the entire system. Most importantly, momentum vortices correspond to the intersections of domains walls. The system also exhibits thermally excited spin vortices which are marked by circles in Fig. 1(d). However, they are associated with much larger excitation energies than momentum vortices. No indications for a Kosterlitz-Thouless transition and a binding into low-temperature vortex-antivortex pairs is found, neither in spin-nor in momentum space. We explain this by the special restrictions imposed on the effective degrees of freedom such as the curl-free condition for q.
Our numerical results at the lowest simulated temperatures also indicate that the domain wall network preferably forms rectangular patterns. This can partially be explained by the underlying square lattice nature of our model but is mainly due to the antivortices (n = −1) which realize the lowest excitation energy cost when four domain walls are radiating from the vortex core. Their characterization through the number and precise arrangement of radiating domain walls allows us to develop a general classification scheme for momentum vortices. Due to Mermin-Wagner theorem the formation of rectangular domains associated with the breaking of momentum symmetry [U (1) → Z 4 ] does not occur in a finite-temperature transition but rather in a smooth crossover separating the spiral spin liquid and the rigid vortex network regimes. Particularly, with decreasing temperature the average domain size continuously increases while thermalization and mobility of momentum vortices significantly slow down. Besides the shear modulus term in the system's continuum theory, further contributions not allowed in elasticity theory (such as a potential term for the momentum) become increasingly important in this low-temperature regime. As a result, a rank-2 U(1) gauge theory description is no longer appropriate which is also seen in the fading of the four-fold pinch points in the electric field correlator (Fig. 8).

III. THE R2-U1/ELASTICITY PHASE
In this section, we discuss the properties of the spiral spin liquid as well as its relation to a rank-2 U(1) gauge theory and elasticity theory. Since fluctuations of momentum vortices are the main driving force behind this phase, we start discussing their precise nature and the constraints imposed by the curl-free condition. At the end of the section, we confirm our findings with numerical results.

A. Vortices of local momentum vector field
Let us first qualitatively describe the arising of vortices from the local momentum vector field. We view the lattice in a coarse-grained way. Each coarse-grained block is over a few lattice sites so that we can define the local momentum vector q as the gradient of the spins' in-plane angle Φ, At the same time the block is not too big so that the momentum vector effectively does not vary within the block. In this way, we have defined a vector field q(r) of the local momentum on the coarse-grained lattice. At low temperatures, each momentum vector takes a value from the spiral momentum ring, or at least within a narrow region around it. Consequently, the configuration space of a momentum vector at each individual coarsegrained block is homotopic to a circle, so the winding number n of the momentum vector field is a well-defined topological quantity on an arbitrary loop on the lattice. Examples of winding number ±1 vortices are shown in Fig. 2.
Momentum vortices are distinctly different from spin vortices which are defined in terms of spin winding along closed loops. In fact, spin and momentum vortices are independent of each other. Furthermore, unlike normalized spins, there is no strict amplitude constraint on the momentum q. As a consequence, momentum vortices can be realized on continuously varying spin textures without any singularity. Another drastic and, at a first glance, unexpected difference concerns the possible winding numbers n of momentum vortices which we will elaborate on in the next subsection.

B. Curl-free constraint on momentum vortices
In the absence of spin vortices, Φ varies continuously on every point of the system. So q(r) must obey the curl-free condition by definition of q(r) [Eq. (7)]. This restriction on q(r) plays a central role in determining the low temperature properties of spiral spin liquids. If spin vortices are taken into consideration, then for any loop (not directly passing through the spin vortex core), one finds where n s is the spin winding number. This implies that a spin vortex at r is a point source of quantized curl for the momentum vector field, However, we are mostly concerned with the physics of spiral spin liquids in the absence of spin vortices. In fact, our numerical results demonstrate that spin vortices are associated with much higher excitation energies than momentum vortices. Hence, from now on, we assume the curl-free condition [Eq. (8)] to be a constraint on the spiral spin liquid system. The curl-free condition restricts the allowed configurations of the momentum vortices. The result and its mathematical deduction are first summarized below before elaboration: 1. As a starting point, the momentum vector field has to be curl-free [Eq. (8)].
2. In the neighbourhood of a singular point, the vector field cannot be of the types focus 3. Vortices of winding number n ≥ 2 are forbidden, because they require at least one elliptic sector. Vortices of winding number n ≤ 1 are allowed, but their configuration is still constrained by the condition above.
Two-dimensional vector fields and their vortex singularities are a well-studied topic in topology as tangent vector bundle sections in 2D, and also dynamical systems [56,57]. Around a vortex singularity (also known as critical point in mathematical literature), the vector field configurations are fully classified. One possibility is that it can be either a focus or a center [ Fig. 3(a,b)], surrounding the entire singularity. The other possibility is that the neighbourhood of the singularity is divided into a few sectors. Here, a sector is defined such that at its boundaries the vector field q(r) obeys q(r) r where the singularity is at r = 0. The vector field configuration in each sector is independent from the neighboring ones (as long as continuity on the boundaries is satisfied) For focus and center configurations the winding number is n = 1. For a singularity divided into several sectors, Poincaré-Bendixon theorem [56,57] states that the winding number is where e (h) is the number of elliptic (hyperbolic) sectors. This completes the classification of all vortex configurations with a point singularity. Now we can examine the curl for each class of vortices. First, it is straightforward to see that focus and center have non-zero curl at the singularity, and hence are excluded in our system. Second, the elliptic sector has nonzero curl in the entire region. This can be seen by following the integration path of the vector field [highlighted in red in Fig. 3(c)]. Since such a path forms a closed loop starting and ending at the singularity, and starting and ending vectors are rotated by a finite angle, we know the enclosed region must have non-zero curl. In particular, the singular point will have divergent curl. Hence these sectors are forbidden too. Third, the hyperbolic sector may have local non-zero curl but not necessarily, and is therefore allowed in our system. Finally, the parabolic sector is strictly curl-free and allowed. Now we can conclude on the vortex restrictions. First, since all elliptic sectors are forbidden, and Poincaré-Bendixon theorem states that vortices with winding numbers n ≥ 2 require at least two elliptic sectors, such vortices are not allowed (Fig. 4). Second, the only possible configuration for vortices of winding number n = 1 is parabolic in the entire neighborhood. Finally, vortices of winding number n ≤ −1 are generally allowed, but they cannot have elliptic sectors, and their hyperbolic sectors still need to be curl-free.
Another interesting property is that the absence of curl requires the vortices to carry a finite divergence distribution. Here, we briefly discuss the divergence of the simplest vortices with winding numbers ±1, since they are most common in our numerical simulations. We note that the winding number 1 vortices carry a net positive or negative, unquantized divergence, as shown in Fig. 5. The two types of winding number 1 vortices cannot be smoothly transformed into each other. Interestingly, the winding number −1 vortices carry a dominant quadrupole of the divergence distribution (Fig. 5). They can also carry a net scalar divergence and a dipole, but these two quantities can be smoothly tuned to zero.

C. Hamiltonian in the continuum limit
The next ingredient to understand the low-energy behavior of spiral spin liquids is the formulation of a continuum theory. Assuming that the spins, denoted by their angle Φ(r), are varying slowly on the lattice, we can rewrite the Hamiltonian of Eq. (1) in the continuum limit. In the parameter regime of Eq. (2), the Hamilto- nian becomes We then expand Φ(r + a) in a Taylor series, After inserting Eq. (13) into Eq. (12), we expand the cosine, keeping terms which contain up to four derivatives in total, which requires terms up to third orders in Eq. (13). The continuum limit is performed by replacing Using q = ∇Φ, and denoting the 2 × 2 Hessian matrix Q as where Greek indices denote Cartesian coordinates, µ, ν ∈ {x, y}, we obtain the Hamiltonian written in q only, with a potential term and a stiffness term: with Here, we have again defined δ = J 2 − 1/4. Our notation indicates that H sq-XY consists of two terms, a potential term H p which determines the energy cost of momentum q at each point. The other term H s determines the stiffness of q, i.e. the energy cost of spatially varying spiral configurations.
It is worth noting that in the limit of small δ and small momentum q, the model becomes rotationally symmetric: where we have only kept the leading order in δ for each power in q. This form makes it obvious that the potential term acquires a standard symmetric 'Mexican hat' shape.
A spin vortex corresponds to a quantized point singularity of the curl of q. Let us first assume that they do not appear in the system such that q is curl-free. Hence, we have and we can rewrite Q as We note that by identifying the momentum vector q in spiral spin liquids as the lattice distortion in elasticity, the matrix Q in Eq. (16) becomes the symmetric strain tensor U .
The term Q µν C µνρσ Q ρσ in Eq. (16) is the Hamiltonian of elasticity with the shear modulus only, and does not contain the compression modulus term [58]. After identifying the strain tensor, we can study the correspondence of topological defects in the two systems. In elasticity, the bond angle ∇ × u is a periodic quantity [44]. By definition it corresponds to ∇ × q in spiral spin liquids, which is the spin vortex density and is quantized as an integer instead of being periodic.
The fundamental topological defects in elasticity are disclinations, which correspond to a winding of the periodic bond angle [44]. Hence, its analog cannot appear in spiral spin liquids, since ∇ × q is not periodic. Similarly, dislocations b as dipoles of disclinations, defined as b µ = ρσ ∂ ρ ∂ σ u µ , do not appear in spiral spin liquids neither.
It turns out that the winding number 1 momentum vortices in spiral spin liquids correspond to quadrupoles of disclinations (or pairs of dislocations), which are vacancies or interstitials (i.e. additional atoms squeezed in the lattice) in elasticity. This is because inserting an additional atom in the lattice makes the other atoms move radially away from it, leading to a lattice distortion field u similar to our winding number 1 momentum vortices, see Fig. 5. These objects have vanishing disclination and dislocation densities. Instead they are manifested as nonzero divergence of u (q), which is also Since the winding number 1 momentum vortices cannot have any curl, the only possibility is to have one parabolic sector in the entire region, which has divergent divergence at the singularity. Examining the negative winding number vortices in the same way, we find them to carry multipoles of ∇ · q, i.e., they are multipoles of vacancies or interstitials, see Fig. 5.
Utilizing the duality between elasiticy and rank-2 U(1) gauge theory discovered by Pretko, Radzihovsky [43], and by Gromov [45], we can also establish a connection between spiral spin liquids and rank-2 U(1) theory. The latter theory is a generalization of Maxwell's electromagnetism, by upgrading the electric field to symmetric tensors, and modifying the definitions of charges, gauge fields, and magnetic fields accordingly. A more detailed analysis can be found in Refs. [39][40][41][42]. Here we only explain the part relevant to our work.
The mapping between the Hessian matrix in spiral spin liquids and the symmetric tensorial electric field is given by To see the existence of a Gauss's law, we notice that for smoothly varying q, the scalar charge (fracton) is always zero in the entire space, where we used ∂ µ ∂ ρ µν = 0. The fracton scalar charges correspond to disclinations in elasticity. There is also a conservation law for dipoles of fractons (i.e., dislocations in elasticity): Note that here, the first term vanishes because ∂ σ σβ q β = ∇ × q = 0, and the second term vanishes because Finally, the winding number 1 momentum vortices are quadrupoles of fractons, and negative winding number momentum vortices are higher multipoles. A quadrupole of fracton corresponds to vacancies or interstitials in elasticity, and manifests as non-zero trace of the electric field, or divergence of q (Fig. 5) We summarize all the relations between spiral spin liquids, elasticity theory and rank-2 U(1) theory in Table I.
We have now established the effective theory that sheds more light on the spiral spin liquid nature. It can be mapped onto classical elasticity (rank-2 U(1) electrostatics) where momentum vortices correspond to multipoles of vacancies or interstitials (fractons). These objects are free to move around in the lattice and lead to a spin liquid-like behavior. This is also consistent with the fact  16)] of homogeneous spiral states with momentum q, also referred to as potential term in the system's continuum theory. Shown is the Mexican hat shape of this potential for δ = 0.03 and two momentum directions, q = (qx, 0) or symmetry related directions (blue) and q = (qx, qx)/ √ 2 or symmetry related directions (red).
that fractons (disclinations) and dipoles of fractons (dislocations) have restricted mobility, while quadrupoles do not. The correspondence between the Hamiltonians is exact at the critical point but the spiral spin liquid effective theory gains additional terms (which are not of the shear modulus form Q µν C µνρσ Q ρσ ) when δ = 0, see Eq. (16).

E. Numerical results
We now demonstrate how our analytical results on momentum vortices in spiral spin liquids as obtained in the previous subsections manifest in actual numerical simulations of the J 1 -J 2 -J 3 square lattice XY model. Here, we focus on an intermediate temperature regime, particularly on T < ∼ T * where T is low enough to enable spiral formation but not too small such that spirals remain liquid due to thermal fluctuations. In the next section, we will investigate the low-temperature rigid vortex network regime.
Our classical Monte Carlo simulations are based on a standard Metropolis algorithm for a quadratic system with 400×400 = 160000 sites and periodic boundary conditions. To reduce the autocorrelation times two third of all spin updates are chosen as overrelaxation steps while the other third are regular Monte Carlo updates. One Monte Carlo run includes a total of 6 · 10 7 Monte Carlo steps during which the system is cooled down from T = 2 to T = 0.005 (in units of |J 1 | = 1) using an exponential cooling protocol. At various selected temperatures, measurements of the energy, heat capacity, spin structure factor, and momentum distribution are performed. Averages are taken over 10 independent Monte Carlo runs. We choose the parameter δ = J 2 − 1/4 = 0.03 in all our simulations below. This value turns out to be suitable for illustrating the physics of momentum vortices, since the potential term H p in Eq. (16) has an almost perfectly rotation symmetric Mexican hat shape (at least in the vicinity of the valley), as shown in Fig. 6. For larger δ, the valley of the Mexican hat potential would lose its circular shape, leading to an entropic selection of discrete momenta, which counteracts spiral spin liquid behavior. On the other hand, choosing this parameter too small, the energy gain from spiral formation relative to the ferromagnetic state would become negligible such that spiral states only appear at very low temperatures. For δ = 0.03, the classical ground states are spirals with q 0 ≈ 0.66 which corresponds to a spiral wave length of λ ≈ 9.5 lattice spacings.
Coming from the high-temperature side, the first pronounced feature in the heat capacity (see Fig. 7) is a broad shoulder at T ≈ 0.32 which marks the onset of magnetic correlations. Above this temperature, the system effectively behaves as a paramagnet without any noticeable features in the spin configurations [for an exemplary spin state see Fig. 1 For temperatures below the broad shoulder the system first enters a regime where correlated spin patterns become visible [see Fig. 1(f) obtained for T = 0.1], however, a clear spiral formation does not yet take place. In the spin structure factor S(k) defined as this manifests in a broad and featureless peak which roughly fills the area enclosed by the valley of the Mexican hat potential [ Fig. 8(b)], hence the name 'pancakeliquid' in Ref. 23. This property indicates that the system can access an extensive manifold of states in momentum space and that in this subspace the potential term H p (q) in Eq. (16) is largely irrelevant. Since the potential term constitutes a key difference to elasticity theory, its irrelevance implies that the system's behavior in this temperature regime is dictated by the analogy to elasticity, i.e., fracton gauge theory. We demonstrate this numerically by plotting the electric field correlation function As excepted for a system subject to a generalized Gauss's law ∂ µ ∂ ν E µν = 0 the electric field correlator needs to obey a specific projector form [55], so that shows a characteristic fourfold pinch point singularity [ Fig. 8(c), (d)]. The pinch point pattern at T = 0.1 does not exactly extend to the Γ point q = (0, 0) which is likely due to a thermal broadening.
Lowering the temperature further, the heat capacity shows a sharp peak at T = T * = 0.08 associated with the formation of spin spirals, see Fig. 1(e). Just below this peak a large density of momentum vortices is observed which are intimately connected with spatial fluctuations of spin spirals, establishing a liquid-like property. Particularly, the spin structure factor now has a ring-like shape indicating that amplitude fluctuations of the spiral momentum are suppressed [ Fig. 8(a)]. On the other hand, fluctuations of the direction of spiral momentum are still pronounced as evidenced by the relatively even distribution of signal along the ring. Note that the diameter of the ring is slightly smaller compared to the valley of the Mexican hat potential. This reduction of the momentum amplitude is characteristic for spin configurations with spatially varying direction of q which will be further studied in Sec. IV C in the context of domain walls. Since spatially varying momentum directions are the origin of momentum vortices, we interpret the reduction of ring size in Fig. 8(a) as indirect evidence for momentum vortices. The increasing effect of the potential term in the low-energy effective theory slightly reduces the intensity of the fourfold pinch points in Fig. 8(c) but the decreasing influence of thermal fluctuations sharpens the pattern at small k.
It is worth commenting on the nature of the apparent phase transition at T = T * = 0.08. Due to Mermin-Wagner theorem the continuous U (1) spin symmetry cannot be spontaneously broken at a finite temperature; it can only be broken algebraically which would correspond to a Kosterlitz-Thouless transition. However, as we will argue based on our low-temperature results in Sec. IV D, we do not observe any indications for such a transition. Spontaneous breaking of discrete time-reversal symmetry is also excluded for a 2D XY model since it acts as (S x , S y ) → (−S x , −S y ), i.e., it is identical to a π rotation in spin space which, however, is a subgroup of the continuous U (1) spin symmetry. In principle, U (1) momentum symmetry could be broken across T * which would correspond to a spiral selection in the degenerate manifold. Particularly, since this symmetry is not an exact continuous symmetry on all energy scales but only an approximate one at small energies and small δ, the breaking could even occur at a finite temperature. We studied the behavior of order parameters near T * which are sensitive with respect to a complete breaking of U (1) momentum symmetry as well as a partial breaking U (1) → Z 2 and U (1) → Z 4 . However, we do not observe any noticeable features of such order parameters at T * which also excludes these types of transitions.
We, therefore, conclude that the heat capacity peak at T * is not related to a real phase transition in the thermodynamic sense but only signals a change of magnetic correlations on finite length scales. Our spin configurations at T < ∼ T * [ Fig. 1(e)] indicate that in this temperature regime the system forms short range spirals whose momenta q are well defined on length scales of a few spiral wave lengths. Beyond this distance spiral momenta show pronounced real-space fluctuations without any long-range patterns, reminiscent of a liquid-like property. The system can gain much energy via this short range spiral formation, as the spin configurations are now pinned near the valley of the Mexican-hat potential. This sudden energy reduction manifests as a pronounced peak in the heat capacity.

IV. RIGID VORTEX NETWORK PHASE
As the temperature decreases, stricter energetic constraints are imposed on the system's spin configurations and the question about the precise nature of energetically optimized momentum vortices arises. We first analytically investigate such low-energy momentum vortices via a spiral domain construction where multiple domain walls are radiating from the vortex core. Thereafter, we demonstrate that such vortices are indeed found in lowtemperature numerical simulations.
The overall emerging picture can be summarized as follows: At small temperatures, the Mexican hat potential term H p becomes dominant, and q at every point in space essentially has to lie exactly on the spiral ring determined by Eq. (4). Furthermore, the curl-free condition [Eq. (8)] requires all domain walls to be straight lines. As a result, the system freezes into a rectangular network of winding number ±1 momentum vortices.

A. Vortices without momentum amplitude variation
Let us now consider the case where the local momentum vectors q(r) take values strictly on the spiral momentum ring [Eq. (4)], except of singular lines (domain walls) and points (vortices). For simplicity, we assume that the spiral momentum ring is an exact circle (as is the case for δ 1) so it can be parametrized by an angle θ , q(r) = q 0 (cos θ(r), sin θ(r)). Instead, the zero curl condition severely restricts the possible momentum vector field configurations. As we will show, the non-trivial, relevant solutions all have a singularity point (vortex core) and a few straight-line domain walls radiating from it. These constructions are also motivated by our numerical results which show a rather accurate segmentation of the system into spiral domains, see Fig. 1(d).

Single domain wall
Let us first consider a single domain wall that forms a straight line of canting angle φ 12 between two regions of constant momentum. In these two regions labeled as 1 and 2, we have respectively The curl-free law for the momentum field requires the projections of q 1 and q 2 in the direction of the domain wall to be the same, which is In terms of the spin textures, this is just requiring the spins on the two sides of the domain wall connect continuously. The non-trivial solution is For a given θ 1 , the momentum in region 2 has only one solution. The corresponding momentum distribution and spin texture is mirror-symmetric regarding the domain wall. An important corollary from this conclusion is that the domain wall between two regions of fixed momentum has to be straight.

Singularity with odd number of discrete domain walls
We now consider a singularity -which will turn out to be a vortex of the momentum vector field -with several domain wall branches radiating from it.
We start with the simplest case of a singularity with three branches, as illustrated in Fig. 9. The momentum vectors and slopes of domain walls are as labeled there.
The solutions are guaranteed to exist and are unique (up to reversing all momentum vectors), due to the fact that the m × m matrix always has a non-zero determinant det K m = 0, for odd m.
We now illustrate the momentum field solutions for different domain wall distributions. The cases of three domain walls at different angles are shown in Fig. 10. There, we plotted the momentum field solutions on top of heat maps of the local energy density, and also the spin textures in separate panels. In Fig. 10(a), the three domain walls are evenly distributed, and the momentum vector field forms a vortex of winding number 1 with vanishing curl everywhere. In terms of the spin texture, the contours of spins form triangle-shaped loops around the center.
As we squeeze the three domain walls to one side of the plane [ Fig. 10(b)], a transition point is reached where the winding number becomes ill-defined. Further narrowing down the angles of domain walls leads to a state illustrated in Fig. 10(c), in which the momentum vector field forms a vortex of winding number −1.
Another transition point is reached when the three domain walls are squeezed into a quarter of the plane [ Fig. 10(d)]. After the transition point, when all three domain walls are confined within an angle of π/2, the momentum vector field has zero winding number [ Fig. 10(e)].
This observation can be generalized to any odd m. When the domain walls are evenly distributed, the momentum vector field forms a vortex of winding number 1. As one "squeezes" the domain walls into a narrower angle, the momentum vector field can go through transitions into vortices of winding numbers −(m − 1)/2 to 0.

Singularity with even number of discrete domain walls
Next we consider a singularity with an even number of domain wall branches radiating from it. Again, determining the momentum vector field is equivalent to solving Eq. (44). However, the situation is different from the odd number case, because det K m = 0, for even m.
More specifically, the matrix K m has m − 1 linearly independent rows. The last row can be written as linear combination of the first m − 1 rows As a consequence, there could be infinitely many solutions if and no solution if this condition is not met. When there are solutions, one can treat θ m as a free parameter, and then solve θ 1 , . . . θ m−1 from the m − 1 linearly independent equations. The solutions then form a 1D manifold parametrized by θ m .
For a clearer physical picture, we illustrate the cases with four and six domain walls in Fig. 11 with varying parameter θ m . In the case of four domain walls, the momentum vector field can form a vortex of winding number 1 [ Fig. 11(a)], and after the transition point [ Fig. 11(b)] transforms into a vortex of winding number −1 [ Fig. 11(c)]. In the case of six domain walls, the momentum vector field can also form a vortex of winding number 1 [ Fig. 11(d)], and after the transition point [ Fig. 11(e)] transforms into a vortex of winding number −2 [ Fig. 11(f)].
The situation can be generalized to any even m. As the parameter θ m varies, the momentum vector field configuration can transit from winding number 1 to −(m − 2)/2.
To summarize, we have analytically constructed all excitations based on domain walls when q is restricted to be on the spiral ring. We have found that the momentum vector field can form vortices around singular points. A peculiar feature is that the vortices can have any negative winding number n < 0 and n = +1, but no winding number n > 1. This is a consequence of the curl-free condition [Eq. (8)], which we have already discussed in detail in Sec. III B.
Since domain walls are associated with an energy cost proportional to their length, all vortices constructed here have the property that their energy cost remains constant as a function of distance from the core. This is in stark contrast to usual spin vortices where the energy usually decays with 1/r 2 . As is intuitively clear but will be discussed in more detail in Sec. IV C, the excitation energy of a domain wall becomes larger with increasing momentum difference across the domain wall. Consequently, the energetically cheapest vortex with negative winding is the n = −1 vortex in Fig. 11(c) with four domain walls arranged in angles of π/2. In a real physical system, the momentum vector field q is allowed to slightly deviate away from the ring associated with a 'potential energy cost'. This lead to a softening of the sharp domain walls, and will be studied in detail in Sec. IV C. However, the key physics discovered here -that the low-energy excitations are vortices of winding number smaller than or equal to 1 -is preserved, and plays a major role in determining the low energy physics of the spiral spin liquid.

Smooth vortices
For the sake of completeness, here we treat the case where the vector field q is analytic everywhere except of the vortex core at r = 0. This will result in a smooth vortex with winding number n = 1 [see Fig. 2 (left)], which does not have the sharp domain walls discussed above. This case can also be understood as the limit of infinitely many domain walls, densely covering the entire plane. However, because all negative winding number vortices are made of a singularity and several straight line domain walls, the smooth winding number n = 1 vortex is actually not observed in numerical simulations due to its incompatibility with n ≤ 1 vortices.
Let us denote the momentum field as q(r) = q(θ(r)) where the polar angle θ = θ(r) defines the inplane orientation of q, see Eq. (33). The curl-free condition [Eq. (8)] then becomes Since we are constructing excited states, we first exclude the trivial case ∇θ = 0. Notice that ∂ θ q ⊥ q, therefore the condition in Eq. (50) implies that must hold. Furthermore, at each point r 0 = 0 the gradient ∇θ(r 0 ) is oriented perpendicular to the contour line of constant θ(r) running through r 0 . This means that the momentum field q(r) aligns in the direction of the contour line of constant θ(r). It follows immediately that contours of constant θ(r) must be straight lines; if they are bent, q would have different orientations on different points of the contour, against the property that the contour has constant θ(r).
In a field θ(r) with the property that lines of constant θ are straight, a singularity at r = 0 implies that these lines cross at r = 0. In other words, lines of constant θ(r) point radially away from r = 0 which means q(r) ∼ e r , where e r is the radial unit vector. Including the normalization of the momentum yields the two solutions q(r) = ±q 0 e r .
One solution is illustrated in Fig. 2. These states correspond to vortices in momentum space with a phase winding of n = +1, as one would intuitively expect. Our proof shows that they are the only viable states. By taking this configuration as a given input state and calculating its classical energy in a spiral spin liquid phase, we find that its energy decays as 1/r 2 , just like for a usual spin vortex. Note, however, that for a spin vortex, there is the freedom to perform a global U (1) rotation of all spins. This is different for a momentum vortex with n = +1 which only allows for the Z 2 transformation of globally inverting the momentum, q → −q, leading to the two signs in Eq. (52).
Finally, we mention briefly that proper mixtures of smooth and discrete domain wall configurations can also be constructed. A hybrid momentum vector field configuration constructed by gluing parts of Fig. 2(left) and Fig. 11(a) together.
cut and pasted together to form a new solution. An example is given in Fig. 12, where we take the right quarter of Fig. 11(a) and the rest of Fig. 2 (left). However, for the same reason of compatibility with negative winding number vortices, these hybrid vortices do not actually occur in numerical situations.

B. Formation of rigid vortex network
Let us now zoom out to the entire lattice, and argue the formation of a rigid vortex network.
First, vortices and straight line domain walls must form a network. Assuming periodic boundary conditions, the entire lattice will have vortices of both positive winding number, which can only be 1, and negative winding numbers. The negative winding number vortices do not have a smooth configuration -they must have several branches of straight, sharp domain walls radiating from the vortex cores. For this reason, the positive winding number vortices cannot have smooth configurations neither, as they are incompatible with domain walls from the negative winding number vortices. Hence, all vortices should have only straight-line domain walls which interconnect into a network. In this network, the edges (links) are the domain walls, and the nodes (vertices) are the vortices.
Furthermore, the configurations of vortices and domain walls that can actually appear in the system are selected by their energy cost, especially at very low temperatures. Numerically, we find that vortices with four perpendicular domain walls are the energetically most favorable ones. This property results from a balance between positive and negative winding number vortices: while winding number 1 vortices usually prefer more domain walls to lower their energy, this does not always apply to negative winding number vortices. Therefore, a rectangular network of vortices and domain walls is expected at low temperatures.
The vortex network is rigid. This can be deduced from the analytical solutions for a given configuration of domain walls [Eq. (44)]. Note that, around a vortex, the four domain wall angles must satisfy Eq. (49) in order for a momentum vector field solution to exist. If we move the positions of one or few vortices, then for some of the neighboring vortices, only one of the domain wall angles will change, which obviously will violate the condition in Eq. (49). The movement of domain walls or vortices must be of sub-system type at least. Hence, the vortex network becomes highly rigid at low temperatures, when the local momentum vectors are almost strictly on the spiral ring and Eq. (49) is enforced.
In a system with a frozen rigid vortex network, it becomes very difficult for the vortices to annihilate each other in order to further reduce the system's energy, since any movement of vortices requires a global, cooperative change of the spin configuration.
These arguments imply that Kosterlitz-Thouless-type behavior is not expected in our system. Fundamentally, this is because of the curl-free condition [Eq. (8)] and the appearance of vortices that break the momentum vector U(1) rotational symmetry. Note that, although the ground state manifold has such a degeneracy/symmetry, the excitations (momentum vortices) do not. More specifically, this leads to the formation of straight line domain walls, with an energy cost proportional to their length, in stark contrast to usual spin vortices. These differences significantly control the thermodynamic behavior of the rigid vortex network phase, and make the Kosterlitz-Thouless phase transition physics inapplicable.

C. Domain wall broadening
To make sure that our analysis in the previous section is valid, we still have to address the question whether our assumption that the momentum vector is strictly confined to the spiral ring actually applies. We will see in this section that, although sharp domain walls do not represent the exact energetically optimal spin configurations in a real system, reintroducing amplitude variations only broadens the domain walls to a finite extend, at a scale much smaller compared to the rigid vortex network. Consequently, the sharp domain wall assumption can effectively be considered as valid.
Let us now study the finite broadening of domain walls that happens in more realistic settings, when the momentum vector field q is allowed to fluctuate away from the spiral ring. Particularly, we will demonstrate that an energy-optimized broadened domain wall has a simple analytic description in the continuum limit. We will calculate the ideal form of the rounded domain wall and show that it comes with a low energy cost if the momentum difference across the domain wall is small. Note that this rounding is necessarily associated with amplitude variations of the momentum. We first specify the type of spin configuration that we will describe with the help of our continuum model. Let us start with a discontinuous (non-optimized) horizontal domain wall along the x axis which separates two homogeneous spin spirals with momenta q + and q − in the upper and lower half planes, respectively. Furthermore, we require that both spirals are ground states, i.e., |q + | = |q − | = q 0 , where the ground state momentum q 0 follows from minimizing the potential term H p in Eq. (16) in the limit δ 1, yielding q 0 = 4 √ δ. Since the spin configuration must be symmetric around the x axis due to the curl-free condition (see Fig. 13), the Cartesian components obey q + x = q − x and q + y = −q − y . In the following, we will characterize the domain wall by α ∈ [0, π/2] which is the angle enclosed by the x axis and q + (or q − ) such that Without loss of generality we have fixed the sign q + y > 0 in this equation. With these definitions, α describes the strength of the discontinuity where α = 0 stands for no domain wall and α = π/2 corresponds to the maximal momentum jump across the domain wall.
The above spin configuration is an excitation where all the energy cost is concentrated along the infinitely narrow domain wall. Upon optimization the energy will spread out into a strip of finite width, and contours of constant spin angle Φ become rounded (see Fig. 13). In the following, we will use the continuum model in Eq. (16) to calculate the ideal momentum distribution q(r) = (q x (r), q y (r)) of this state. Most importantly, the optimized momentum configuration is still translation invariant along the x axes, such that lines of constant Φ transform into each other by parallel shifts along the x axes. As a result, the functional dependencies of q(r) reduce to q(r) ≡ q(y) = (q x , q y (y)) , (54) which means that q x is constant across the entire spin configuration and q y is only a function of y such that the optimization becomes an effective 1D problem. This also guarantees that the curl-free condition [Eq. (8)] is not violated in the optimization process. Using this property to simplify the continuum model in Eq. (16) and exploiting Eq. (53) one finds that the energy E per length l of the domain wall in leading order in δ is given by Here, terms constant in y are neglected (these terms, however, may still have an α dependence). This functional needs to be minimized with respect to q y (y) where the boundary conditions follow from the fact that far away from the domain wall the spin configurations are given by the initial homogeneous spiral states with momenta q + and q − : Using Euler-Lagrange equation this leads to the differential equation The solution respecting the boundary conditions in Eq. (56) is found to have a simple form: q y (q) = q + y tanh(q + y y) .
It is worth highlighting two properties of this result. First, from the argument of the hyperbolic tangent, the width ∆y of the optimized domain wall is given by which indicates that domain walls with small α have a diverging width. On the other hand, amplitude variations of the momentum are small in this limit such that, in total, these domain walls are still energetically cheap. For a typical angle α = π/4 which is realized for the domain walls of momentum antivortices with the lowest possible energy one obtains ∆y = 1/ √ 8δ. Inserting δ = 0.03 as used in our numerics yields a rather small width of ∆y ≈ 2 lattice constants. [It is, of course, questionable whether our continuum model where Eq. (13) is truncated still approximates the system reasonably well when momentum variations occur on such short distances. Better estimates are expected for smaller δ.] Our numerics results below confirm that domain walls are narrow with a widths of only a few lattice spacings.
Second, it is instructive to calculate the maximal excitation energies in the center of a domain wall at y = 0. To this end, we compare the Lagrangian L in Eq. (55) (which provides the energy per unit area, or equivalently, per site) at y = 0 (where q y = 0 and ∂ y q y = q + y 2 ) and at y → ∞ (where q y = q + y and ∂ y q y = 0). The maximal excitation energy per site of a domain wall with angle α then reads as The scaling ∆E α ∼ α 4 , which holds for α 1, implies remarkably small excitation energies in the limit of vanishing α. One may conclude that the formation of domain walls with small α provides a natural source for low-energy thermal fluctuations.
For a typical domain wall with α = π/4 the excitation energy is given by ∆E π/4 = 8δ 2 . Comparing this value with the local excitation energy of a trivial ferromagnetic state, ∆E fm = 16δ 2 , which defines another characteristic energy scale of the system, one realizes that ∆E π/4 cannot be considered small. The fact that such domain walls are still observed at small temperatures (see below) is related to the fact that this energy is confined to a narrow strip.
One may conclude that this confinement of energy represents the key property of domain walls and implies that they remain well defined and narrow even far away from vortices.

D. Numerical results
Having discussed the physical properties of energyoptimized momentum vortices and their domain wall construction we now study their occurrence in our lowtemperature Monte-Carlo simulations at δ = 0.03.
Further decreasing the temperature in the spiral spin liquid regime investigated in Sec. III E, the heat capacity (Fig. 7) remains completely featureless below T * = 0.08. This indicates that all changes which a spiral spin liquid undergoes when cooling it down are smooth crossovers while sharp phase transitions are not observed. Note that in the low temperature limit, the heat capacity approaches C = 1/2, as is expected for a XY model with one quadratic mode per site.
The most obvious qualitative change between T = 0.05 and T = 0.005 in Figs. 1(d) and 1(e) is that the spiral domains become well-defined and the domain walls straighten, turning into a network of rigid lines running through the entire system. The intersections of domain walls define the locations of momentum vortices which have winding numbers n = 1 or n = −1. Note that in agreement with our analytical investigation of energy optimized spin configurations, all vortices with n = −1 have four domain walls radiating from the center with angles close to π/2 between them. As discussed above, n = 1 vortices could in principle be constructed in a smooth way without domain walls, see Fig. 2. However, due to their incompatibility with n = −1 vortices we observe that in most of our numerical outputs n = 1 vortices adapt the geometry of n = −1 vortices and, likewise, show four domain walls radiating from the core.
No binding effects between vortices and antivortices are observed; on the contrary, with decreasing temperature, their distance increases. We, hence, exclude the possibility of a Kosterlitz-Thouless transition resulting from U(1) momentum symmetry. The same applies to spin vortices which are occasionally seen in the numerical outputs, even at the lowest simulated temperatures, see encircled spin configurations in Fig. 1(d) (in the spin pattern of an ideal spiral state, a spin vortex has a similar shape as a dislocation in a regular crystal). None of our independent Monte Carlo runs reveals a binding into spin vortex-antivortex pairs.
In most of our low-temperature spin configurations from Monte Carlo, the energetic preference of n = −1 vortices with four domain walls leads to regular squareshaped patterns that spread over the entire lattice. This goes long with the selection of four momenta along the spiral ring, as is evident in the spin structure factor at T = 0.005, see Fig. 14(a). In the case of a perfect circular symmetric potential term H p , these four momenta could also occur at any rotated positions on the spiral ring. However, the small deviations from rotation symmetry at finite δ pin them along the cartesian q x and q y axes in Fig. 14(a). For this reason, a finite-temperature spontaneous breaking of U (1) momentum symmetry down to Z 4 cannot occur; rather the spiral selection is continuous as a function of temperature. Even for a perfectly rotation invariant potential H p (i.e. in the limit δ → 0) a finite temperature transition associated with spontaneous momentum symmetry breaking U (1) → Z 4 would be suppressed due to Mermin-Wagner theorem. Overall, the dominant effects of H p and the momentum pinning makes a description in terms of a rank-2 U(1) gauge theory inaccurate and, consequently, the four-fold pinch points in the correlator C EE (k) fade drastically, see Fig. 14(b).
To discuss the system's behavior below the spiral transition at T = 0.08 in more detail, we show in Fig. 15 the real-space momentum and energy distributions at T = 0.005 and T = 0.05. The formation of spiral patches and rigid networks of domain walls is clearly visible when comparing the local momentum directions in Fig. 15(a) and Fig. 15(b). Furthermore, since domain walls have reduced momenta, the network can be made visible by plotting the momentum amplitudes, see Fig. 15(c) and Fig. 15(d). While at T = 0.05 only faint indications of domain walls are visible, they become very pronounced at T = 0.005. The domain walls also show up as lines of enhanced energy, see Fig. 15(e), however, in a less pronounced way as in the momentum amplitude. [In order to make the network of domains walls visible in the energy a Gaussian smoothing of the data with a standard deviation of σ = 3 lattice spacings has been performed in Fig. 15(e) and 15(f).] The spin configuration in Fig. 15(f) also reveals that momentum antivortices (n = −1) cost more energy than momentum vortices (n = 1). This is because n = 1 vortices can reduce their energy by realizing spin patterns with an approximate circular symmetry near the vortex core. As an additional consequence of this freedom, we observe that n = 1 vortices typically show larger deviations from π/2 angles between their domain walls than n = −1 vortices. Clearly, however, the excitations with the largest energy in Fig. 15(f) are spin vortices which appear as narrow and high peaks in the energy landscape. This shows that topological defects in spin and momentum space occur on two different energy scales.
It should be emphasized, however, that at the lowest simulated temperatures (T ≈ 0.005) the typical time scales of thermal relaxation become exceedingly large. This is because of the rigidity of domain walls which cannot undergo any local moves but can only be modified when changing the spin state over the entire lattice. Therefore, it is possible (or even likely) that our numerically obtained low-temperature spin configurations do not represent thermal equilibrium but are rather metastable states.
In summary, our low-temperature Monte-Carlo results reveal a well-defined network of narrow domain walls with a tendency for rectangular patterns minimizing the energy of momentum antivortices. While the domain sizes grow with decreasing temperature a transition into a phase where U (1) momentum symmetry is spontaneously broken down to Z 4 (or lower symmetries) is not observed. In contrast to fluctuations in momentum space which take place on relatively small energy scales, spin vortices appear as massive and locally confined excitations.

V. DISCUSSION
In this work, we have investigated the low-energy behavior of spiral spin liquids and have identified various unexpected properties and connections to other fields in physics. The key reason for the non-trivial behavior of spiral spin liquids lies in its effective U(1) degree of free-dom which corresponds to the direction of the spiral momentum q on the degenerate ring-like ground state manifold. Unlike the elementary spin degree of freedom, the spiral momentum, which we define as a vector field q(r) on a coarse-grained lattice, is subject to a curl-free condition, ∇ × q(r) = 0. This has drastic consequences on the nature of excited spin configurations, particularly, we prove that momentum vortices can only have winding numbers equal or smaller than one while higher-windingnumber vortex types are strictly forbidden.
Interestingly, even though the momentum field q(r) is directly related to the spin degree of freedom via q(r) = ∇Φ(r) (where Φ is the spins' inplane angle) spin vortices and momentum vortices are independent excitations where the latter ones cost a much lower energy than the former. There is, hence, a temperature regime where thermal fluctuations mostly affect the momentum direction of spin spirals while the momentum amplitude is fixed near the ground state ring, leading to spin configurations with large densities of momentum vortices, see Fig. 1(e). This is the spiral spin liquid regime whose existence we have numerically demonstrated for the classical square lattice J 1 -J 2 -J 3 XY model.
We have further demonstrated that the precise mechanism behind the momentum fluctuations in the spiral spin liquid bears striking similarities with elasticity theory of crystals, particularly, our low-energy continuum spin model can be directly mapped onto the shear modulus term in elasticity. However, due to the constrained nature of the momentum field, the topological defects of crystals -disclinations and dislocations -do not have any analogues in a classical spiral spin liquid. Rather, a momentum vortex with positive winding can be considered either as a bound pair of dislocations or a quadrupole of disclinations. Given the mapping between elasticity theory and rank-2 U(1) gauge theory for fractons one may alternatively describe the low-energy behavior of spiral spin liquids by a tensor gauge theory electrostatics subject to a generalized Gauss law. However, since deconfined charges (i.e. fractons) and dipoles of charges are strictly absent in spiral spin liquids, the thermal fluctuations of momentum vortices follow the low-energy behavior of an electrostatics tensor gauge theory where only fracton quadrupoles (or higher multipoles) are present. We make this connection explicit by numerically resolving the four-fold pinch points in the electric field correlator which are characteristic for these fracton theories.
Cooling down the spiral spin liquid, the fluctuations in the momentum amplitude continuously freeze out. This freezing is again controlled by the curl-free condition for q which dictates that the low-energy momentum antivortices must be formed by the intersections of straight and narrow domain walls separating regions of different momenta q. We developed a classification scheme for momentum vortices based on the number of domain walls radiating from the vortex core and found that vortices with winding n = −1 have the smallest excitation energy when four domain walls with π/2 angles between them emanate from the center. Even though momentum vortices with n = 1 can in principle be constructed in a smooth way without any domain walls they are incompatible with n = −1 vortices when arranging both in the same system. As a result, our low-temperature numerical simulations show rigid networks of domain walls connecting momentum vortices and antivortices where angles of π/2 between momenta in adjacent domains are energetically preferred. Interestingly, in such a network state no local excitations can be made in a way that domain walls are only modified, inserted or deleted in a finite lattice region. The simplest low-energy process of changing the domain wall configuration consists of shifting a domain wall parallel to itself over its full length. As a result of this non-local dynamics, equilibration times are increasing significantly as the temperature is decreased such that the system gets easily trapped in metastable states. A typical observation in our numerical outputs is that the domain wall network becomes more wide-meshed with decreasing temperature.
Our work sheds light on the nature of 2D spiral spin liquids, and also opens gateways to analyzing a plethora of related problems. First, its generalization to 3D is highly non-trivial. The degenerate spiral ground state wave-vectors in 3D can form a 1D ring, a 2D sphere, a sphere with punctures [27][28][29], or other manifolds with or without boundaries [24]. In each of these cases, the topological defects of momentum vectors are different, and their classification will be an indispensable step in understanding 3D spiral spin liquids. Second, insights from this work may help us to understand the quantum version of spiral spin liquids. By identifying the classical spiral spin liquid as the electrostatics of a rank-2 U(1) theory, it is reasonable to speculate that the quantum model will at least carry some features of these effective theories. We note that our construction seems to have a natural extension to the quantum model very similar to the higher-rank deconfined quantum criticality studied by Ma and Pretko [59].
We also highlight that the local momentum vortices are a new type of topological defect with a rather exotic curl-free constraint. Their detailed properties await more in-depth study. For example, while spin vortices are only topological for XY spins, local momentum vortices are also well defined for Heisenberg spins. They have already been observed in studies of Heisenberg spiral spin liquid models on the honeycomb lattice by Shimokawa and Kawamura [23]. Closely related to this, the precise melting process of the rigid vortex network as temperature increases represents an interesting statistical physics problem.
Our theoretical study has direct applications in experiment. Among various spiral spin liquid materials McSc 2 S 4 [27][28][29], MgCr 2 O 4 [30], CoAl 2 O 4 [31], NiRh 2 O 4 [32], and Ca 10 Cr 7 O 28 [33][34][35][36][37][38], our theory is most relevant to the 2D quantum spin liquid material Ca 10 Cr 7 O 28 . It has been shown that in classical limit, this bilayer breathing kagome magnet can be mapped to a honeycomb lattice model at low temperatures, and exchange parameters place it very close to the spiral spin liquid phase on the honeycomb model [37,38]. This phase is essentially the same as the one studied in our work. The spiral ring in the spin structure factor [ Fig. 8(a)] has already been observed in neutron scattering [33,34]. More direct experimental tests of our theory would be the search for four-fold pinch points in the correlator defined in Eq. (31) [ Fig. 8(c)], or taking direct snapshot of spin configurations using cutting-edge technology like electron holography [60]. Despite the challenges of such studies, we strongly believe that the investigation of spiral spin liquids in real materials represents a fruitful research direction.