Topological Route to New and Unusual Coulomb Spin Liquids

Coulomb spin liquids are topological magnetic states obeying an emergent Gauss' law. Little distinction has been made between different kinds of Coulomb liquids. Here we show how a series of distinct Coulomb liquids can be generated straightforwardly by varying the constraints on a classical spin system. This leads to pair creation, and coalescence, of topological defects of an underlying vector field. The latter makes higher-rank spin liquids, of recent interest in the context of fracton theories, with attendant multi-fold pinch points in the structure factor, appear naturally. New Coulomb liquids with an abundance of pinch points also arise. We thus establish a new and general route to uncovering exotic Coulomb liquids, via the manipulation of topological defects in momentum space.

One of the most appealing aspects of many-body systems is their ability to produce large scale behavior described by emergent degrees of freedom very different from the system's microscopic constituents. The emergence of a gauge field at low energy is taken as a defining feature of a topological state of matter. A prominent example occurs in the Coulomb phase of spin ice [1][2][3][4][5][6][7], in which the strongly interacting spins are described by the fluctuations of a (otherwise free) field constrained by Gauss' law: Analogous descriptions exist for other frustrated models [8][9][10][11] in which the low energy configurations obey a constraint which can be encoded as a Gauss' law for an emergent field or fields. The resulting state is called a Coulomb spin liquid, in analogy to the Coulomb phase of a U(1) gauge theory.
Recent work has highlighted extensions in which the emergent field is a higher-rank symmetric tensor E µν , and obeys a generalisation of Gauss' law [12][13][14][15][16], e.g.: Higher rank Coulomb liquids are of particular interest due to a close connection with fracton theories [17], sought after for their exotic nature and their potential usefulness in quantum information contexts [18].
There are only few examples of simple spin models exhibiting Coulomb liquids, with higher-rank cases rarer still. Also, we lack a systematic understanding for how to distinguish, and account for, various different types of Coulomb liquid. Further, there has been little exploration of possible transitions between types, or of how lower-rank Coulomb liquids relate to higher rank ones.
Here we introduce an approach which, besides filling in those gaps, enables the discovery, and provides a description of, new Coulomb liquids in classical spin systems almost mechanically, and shows how higher rank Coulomb liquids arise at transitions between distinct lower-rank ones. In so doing, we establish a new understanding of the topological nature of these spin liquids.
In the remainder of this paper we first introduce the general framework, which we then illustrate using two classical spin models, one on the honeycomb lattice [ Fig.  1-2] and one on the octochlore lattice [ Fig. 3-4] (corner sharing octahedra) which between them realise a multitude of new and unusual Coulomb liquids.
General approach: Our approach is based on relating the ground state constraints in real space (Eq. (8) in the honeycomb example below) imposed by the spin liquid's Hamiltonian (Eq. (7)) to the topological properties of a vector function in momentum space L(q). For Heisenberg spins, the real-space constraint quite generally takes the form where c are some real space clusters (e.g. triangles in the kagome lattice, tetrahedra in the pyrochlore lattice) and η i are some real coefficients [19].
If the system has translational symmetry, then Eq. (3) can be rewritten in momentum space: where there are n u sites per unit cell, and the sum in Eq. (4) runs over the n u translationally inequivalent sublattices. S m (q) is the lattice Fourier transform of the spin configuration on the m th sublattice. Eq. (5) defines a set of n u -component vectors, with members of the set being indexed by (p). The number of vectors L (p) is equal to the number of real space constraints in one unit cell. The sum in Eq. (5) runs over spins belonging to sublattice m within a single constrained cluster. r i −r cp is the position of spin i relative to the center of the cluster. This representation is inspired by Henley [9], who related spin correlation functions to a projection matrix P(q), enforcing Eq. (4) by projecting out all L (p) (q):  m and n are sublattice indices and κ is a normalisation constant to enforce the sum rule on the correlation function [20]. Eq. (6) is an approximation: specifically it may be considered as a zero-temperature limit of the leading order of a 1/N expansion, with N the number of spin components. However, we emphasize that this approximation is not necessary to define L (p) (q), although it is helpful in establishing the relationship between L (p) (q) and the spin correlation functions.
Points in q-space where one of the L (p) (q) vanishes, or where one L (p) (q) becomes linearly dependent on the others, have special significance since at those momenta Eq. (4) is satisfied trivially for at least one value of p. At such points the projection matrix P mn (q) becomes singular, and there will be a corresponding singularity in the structure factor (a pinch point in most known examples).
Crucially, depending on n u , and the symmetries of the system, the zeros of L (p) (q) can carry a topological charge. They will therefore be robust against small modifications to the constraint (3), provided that the symmetries necessary to define the charge are maintained. Larger modifications can cause pair-creation or annihilation of topological defects, and thereby a transition to a different spin liquid with more/fewer singularities in the q-space correlations. By manipulating Eq. (3) one can also coalesce topological defects into higher-charged objects.
We find that manipulations of this kind can generate new Coulombic spin liquids and transitions between them. Further, the coalescence of q space defects into objects with higher topological charge is associated to the formation of multi-fold pinch points in the structure factor, (i.e. pinch points with more than 2 intense lobes). Such singularities are signatures of higher-rank spin liquids described by tensor gauge theories [15,16], and we therefore establish a simple mechanism for the construction of higher-rank Coulomb spin liquids.
Honeycomb Model: We consider a model on the honeycomb lattice [ Fig. 1], with the Hamiltonian: The sum in Eq. (7) is over hexagonal plaquettes of the lattice. The first sum in Eq. (8) is over spins on the plaquette (filled black circles in Fig. 1(a)) and the second is over spins connected to the outside of the plaquette (red open circles in Fig. 1(a)). γ is a dimensionless parameter with which we tune the model. Ground states of Eq. (7) everywhere satisfy the constraint For γ = 0 this model reduces to the one studied in [11]. We this constraint in Fourier space, using Eqs. (4)-(5). With only one constraint per unit cell, i.e. only one vector L(q), we can drop the index (p). There are two sites per unit cell, so L(q) has two, complex, components: L 1 (q), L 2 (q), which can be calculated from Eq. (5). Since Eq. (4) is invariant under rescaling L(q), we normalizeL(q) = L(q)/|L(q)|.
Inversion symmetry imposes thatL 1 (q) =L 2 (q) * . L(q) can thus be written as Thus,L(q) for the honeycomb model [Eq. (7)] can always be mapped to a position on the unit circle, and can host stable vortices. At the position of these vortices L(q) must vanish, and there is a corresponding singularity in the spin correlation function [Eq. (6)]. The vortices have an integer winding number defined for closed paths C in momentum space: The total winding number is conserved under smooth changes to the constraint Eq. (3) (e.g. varying γ). New vortices can only be created in oppositely charged pairs. At γ = 0, calculation ofL(q) reveals vortices at the corners of the Brillouin Zone (BZ) and their positions coincide as expected with the location of pinch point singularities in the structure factor [11] Evolution of S(q) across distinct honeycomb spin liquids [ Fig. 1]. For small values of γ there are vortices in L(q) and corresponding pinch points in S(q) at the Brillouin Zone corners. At γ = 1/3 there is a nucleation of pairs of oppositely charged topological defects at the zone boundary which then migrate towards the zone corners. At γ = 1/2 the coalescence of these defects into objects with charge Q = ±2 manifests in the appearance of 4-fold pinch points in the structure factor, associated to the emergence of a Rank-2 Coulomb liquid. On increasing γ further these defects separate and the system enters a new spin liquid with 6 pinch points in the interior of the Brillouin zone, as well as at the zone corners. The left half of each panel is the result of a Monte Carlo simulation of N = 1920 spins at T = 0.002J and the right half is a calculation using the projection approach [Eq. (6)].
The evolution of S(q) at finite γ is shown in Fig. 2. The left half of each panel shows S(q) calculated from a Monte Carlo simulation [21], the right half a calculation using the projection method [Eq. (6)] [9]. The good agreement between the two indicates that conclusions drawn from analysis of L(q), the basis of the projection calculation, are robust.
For γ < 1/3 the only vortices are those at the BZ corners, and the corresponding pinch points in S(q) remain robust. At γ = 1/3 pairs of vortices with opposite winding numbers nucleate at the zone boundaries. For γ > 1/3 these vortex pairs separate along the zone boundary. There are now 8 rather than 2 vortices per BZ, and the same number of pinch points in S(q). The appearance of these new singularities in the structure factor demonstrates that the system is in a qualitatively distinct Coulomb liquid [ Fig. 1(b)].
As γ increases further, the vortices migrate along the zone boundary. At γ = 1/2, three positive (negative) vortices converge on a single negative (positive) vortex at the zone corner, to make vortices with winding number Q = ±2. This has a striking consequence for S(q): the appearance of four-fold pinch point singularities. These are known as signatures of higher rank Coulomb spin liquids with tensor electromagnetic fields [15,16]. Indeed, on coarse graining the model we find that the spin liquid at γ = 1/2 can be described in terms of fluctuations of a traceless, symmetric tensor m µν subject to the constraint: The relationship between m µν and the microscopic spins is given in the Supplemental Material [21]. The components of m µν are formed from the local order parameter for antiferromagnetic order with wavevector at the BZ corners. Assuming slow variation of m µν in real space (or, equivalently, expanding the q space constraints around the BZ corners), one can use a gradient expansion to turn the microscopic contraints on the spin configuration into constraints on the spatial variation of m µν . At γ = 1/2 the first derivative term in this expansion vanishes, leaving a second derivative term given by Eq. (13), which is a generalized Gauss' law for a higher rank Coulomb liquid [14]. It is thus apparent that the merging of topological defects in the BZ generates a higher rank spin liquid. The existence of this higher rank liquid is confirmed by our Monte Carlo simulations of S(q [ Fig. 2], which exhibit four-fold pinch points at γ = 1/2 [15].
However, this higher rank spin liquid is unstable against varying γ. For γ > 1/2 the topological defects separate once again, and the four fold pinch point splits up into four two-fold pinch points. Two of these per BZ remain at the BZ corners, while the others migrate into the interior of the BZ. The process of appearance, merging and separation of the topological defects is similar to the behavior of Dirac points in models of graphene with further neighbor hopping [22,23].
In summary, the honeycomb model [Eq. (7)] exhibits three distinct classical spin liquids with zero temperature transitions between them. These are associated in one case to topological defect creation/annihilation and in the other to a coalescence of defects leading to a higherrank spin liquid.
Octochlore model. We now turn to our second example: a model on the three dimensional octochlore lattice [ Fig.  3]. The Hamiltonian is M oct,α,β = i∈oct S i + α i∈ oct The sum in Eq. (14) is over octahedra. The first sum in Eq. (15) is over spins on the octahedron (red sites in  Table I. Fig. 3(b)), the second is over spins nearest to the outside of the octahedron (blue in Fig. 3(b)) and the third over the next nearest spins to the octahedron (yellow in Fig.  3(b)). α and β are tuning parameters. Ground states obey the constraint M oct,α,β = 0 everywhere.
There are three sites, and one constraint, per unit cell so there is one three-component constraint vector L(q). Inversion symmetry forces L(q) to be real. In three dimensions, this allows topological defects in L(q), with integer topological charge given by the skyrmion winding number of the normalized vectorL(q) around closed 2D surfaces.
Varying α, β enables the generation of several liquids, distinguished through number and arrangement of topological defects in reciprocal space. The phase diagram is shown in Figure 3(c), with the definition of the different liquids in terms of the arrangement of defects in the BZ in Table I. Defects always appear at the zone corners q = (±π, ±π, ±π), and can additionally appear at momenta displaced from the zone corner along high (π, π, π) (π, π, π)+ (π, π, π)+ (π, π, π)+ (q, q, q) (q, q, 0) (q, 0, 0)  Fig. 3(c), distinguished by the number and arrangement of topological defects in the Brillouin Zone. Defects always appear at the Brillouin zone corners (π, π, π) or positions displaced from the zone corner along high symmetry directions q-space. Entries in the table refer to the number and charge of such defects in each CSL.  Fig. 3(c). L(q) exhibits topological defects at (π, π, π) and at wavevectors removed along the (1, 0, 0), (1, 0, 0) and (1, 1, 1) directions from that point. The structure factor exhibits pinch point singularities at the corresponding positions. (b) S(q) at the boundary between regions II and VI of Fig. 3(c). There are Q = 7 topological defects at the zone corners which manifest as multi-fold pinch points, associated to a higher rank Coulomb liquid. The left half of each panel is the result of a Monte Carlo simulation of N = 31944 spins at T = 0.01J and the right half is a calculation using the projection approach.
symmetry directions. Once again, defects in L(q) correspond to singularities in S(q). As seen in Table I, some liquids feature many defects in the BZ, including up to 27 in the case of spin liquid IV. This results in an abundance of pinch points in S(q) as shown in Fig. 4(a).
Multifold pinch points indicating higher rank spin liquids emerge at certain boundaries, for example the boundary between Liquids II and VI. In this case, starting from II, 8 Q = −1 defects placed at (π ±δ, π ±δ, π ±δ) converge on a Q = 1 defect at q = (π, π, π). The resulting Q = −7 defect at the transition results in a complicated structure in S(q), which has the appearance of a 6-fold pinch point when cut through the (h, h, l) plane [ Fig. 4(b)].
Summary and Outlook -We have demonstrated a new approach to the discovery of models exhibiting exotic Coulomb liquids, including higher rank Coulomb liquids. This approach is built on relating the ground state constraints which define classical spin liquids to the topological properties of a vector in reciprocal space. We discover several new Coulomb liquids, all corresponding to simple Hamiltonians on local clusters, some with large numbers of pinch points per Brillouin Zone in S(q), and topological transitions between them, some of which exhibit higher rank spin liquids connected to the potential realization of fractons.
The ability to distinguish classical spin liquids using topological properties of the constraint vector L(q) suggests the possibility of a comprehensive classification of classical spin liquids, in the spirit of the established classification of topological insulators and Weyl semi-metals [24][25][26][27]. We will explore this classification in forthcoming work.
1 Max Planck Institute for the Physics of Complex Systems, Nöthnitzer Str. 38, Dresden 01187, Germany

DETAILS OF MONTE CARLO SIMULATIONS
Here we give some details about the Monte Carlo simulation used to calculate S(q) for the honeycomb and octochlore models in Fig. 3  For simulating S(q) we take a cluster of N sites, and seek to calculate S(q) at a temperature T 0 . The simulations are started at a higher temperature T = T start and gradually cooled down to the base temperature with N eq MCS at each of N T logarithmically spaced intermediate temperatures. After a further N eq MCS for equilibriation at T = T 0 , S(q) is averaged over N measure measurements, each separated by N sep MCS. The same procedure occurs over N par parallel runs of the simulation, with the final result reported being the average of those runs.
For the honeycomb simulations we calculated S(q) at T 0 = 0.002J on a cluster with N = 1920 spins. The other settings of the simulation were T start = 10J, N T = 60, N eq = 10000, N sep = 1000, N measure = 6000, N par = 25.
For the octochlore simulations we calculated S(q) at T 0 = 0.01J on a cubic cluster with N = 31944 spins. The other settings of the simulation were T start = 10J, N T = 50, N eq = 5000, N sep = 2000, N measure = 1250, N par = 200.

LONG WAVELENGTH THEORY OF THE HONEYCOMB MODEL AT γ = 1/2
Here we give a derivation of Eq. (13) of the main text, in which constraints on the honeycomb spin liquid at γ = 1/2 are, after coarse graining, expressed as a constraint on the second derivatives of a traceless tensor field m µν .
∂ µ ∂ ν m µν = 0 (S1) To obtain this result, we consider the constraint on a single hexagon h, centred on position r h of the lattice [ Fig.  S1].
The ground state constraint on h is where the first sum runs over spins belonging to h and the second sum runs over spins connected to the exterior of h. σ is a spin component index. Since different spin components are not coupled by the constraint we drop this index from now on. The only difference this makes is that for O(N ) spins we should remember that there are in principle N copies of each of the fields defined below, one for each spin component. We introduce some complex variables S K,j , defined as where r j is the position of site i and with a the honeycomb bond length. We then further define where ζ j = ±1 is a sign factor which is +1 for sites on the 'A' sublattice and -1 for on the 'B' sublattice (see Fig. S1), and t j , v j are real. In terms of t j and v j the constraint is then To coarse grain, we now introduce continuous fieldst(r),ṽ(r), defined such that when they are evaluated at the lattice sites r = r j , they return the values of the t j , v j on those sites. Assuming smooth variation oft(r),ṽ(r), we can Taylor expand around the center of the hexagon r = r h to obtain Inserting this into Eq. (S6) we obtain the coarse grained constraint ont(r),ṽ(r) 3ai(2γ − 1)(∂ x t + ∂ y v) + 3a 2 4 (1 + 4γ)(−∂ 2 x t + ∂ 2 y t + 2∂ x ∂ y v) = 0.
When γ = 1/2, the first derivative term vanishes, and we are left with If we define a traceless matrix: Eq. (S10) then takes exactly the form of Eq. (13) of the main text.