Classification of Classical Spin Liquids: Detailed Formalism and Suite of Examples

The hallmark of highly frustrated systems is the presence of many states close in energy to the ground state. Fluctuations between these states can preclude the emergence of any form of order and lead to the appearance of spin liquids. Even on the classical level, spin liquids are not all alike: they may have algebraic or exponential correlation decay, and various forms of long wavelength description, including vector or tensor gauge theories. Here, we introduce a classification scheme, allowing us to fit the diversity of classical spin liquids (CSLs) into a general framework as well as predict and construct new kinds. CSLs with either algebraic or exponential correlation-decay can be classified via the properties of the bottom flat band(s) in their soft-spin Hamiltonians. The classification of the former is based on the algebraic structures of gapless points in the spectra, which relate directly to the emergent generalized Gauss's laws that control the low temperature physics. The second category of CSLs, meanwhile, are classified by the fragile topology of the gapped bottom band(s). Utilizing the classification scheme we construct new models realizing exotic CSLs, including one with anisotropic generalized Gauss's laws and charges with subdimensional mobility, one with a network of pinch-line singularities in its correlation functions, and a series of fragile topological CSLs connected by zero-temperature transitions.

The wide variety of observations and discovery has naturally motivated an attempt to classify as comprehensively as possible the various phases that are imaginable.This programme has made tremendous progress, e.g. for electronic band structures alone, we can now distinguish various kinds of stable [5,39,40] and fragile [41] topological insulators and multiple forms of topological semimetal [42][43][44].Indeed, such classification schemes have themselves taken on a wide variety of guises.They range from the rather compact description of the ten-fold way [5] to the classification of Z 2 quantum spin liquids, which has found a bewildering variety of possibilities [45], to a more general classification of two-dimensional [23] and three-dimensional [24] topological phases with internal or crystalline symmetries.
Spin liquids have been at the forefront of topological physics for quite some time, with the resonating valence bond liquid [46] proposed already in the early 1970s, but not discovered and identified as a topological phase until much later [47,48].
In a separate development, the search for disordered magnetic ground states was pursued in the context of spin glasses, where the role of magnetic frustration was identified as a crucial ingredient for destabilising conventional ordered states [49].Since the foundational works of Anderson and Villain, the field of frustrated magnets has grown into a huge field of its own, and proposals of spin liquids as well as candidate materials have become increasingly plentiful [50][51][52][53][54][55][56][57].
Despite their apparently simpler character, there exists, as yet, no similarly comprehensive formalism for CSLs.However, CSLs are of considerable interest in their own right.With their extensive degeneracies, they represent the extreme limit of the consequences of frustration.Even though CSLs are always found at fine-tuned points in parameter space at T = 0, their large entropy allows them to spread out in the surrounding phase diagram at finite T .They are thus relevant to the finite temperature behavior of real frustrated magnets.They may also serve as a starting point in discovering QSLs, with several of the most prominent QSL models having a classical counterpart with a CSL ground state [65,75,[81][82][83][84].
Amongst classical spin models with continuous spinsof which the Heisenberg model is the most familiar member -there exist a number of well-established CSLs (see Table III for a survey).The first was the Heisenberg antiferromagnet on the pyrochlore lattice, which exhibits an emergent U(1) gauge field in the low-energy description of its so-called Coulomb phase [58,59,61,62].
For a long time, it has seemed that the number of distinct CSLs, in the sense of a classification, is quite limited.However, recent work has begun to uncover a landscape of classical spin liquids beyond the "common" U(1) Coulomb liquids, both at the level of effective field theories and microscopic models.In this vein, there have been proposals of short-range correlated spin liquids [64], higher-rank Coulomb phases [68,69,85], and pinch-line liquids [66].This has brought the tantalising promise that there may be quite a large uncharted landscape of possibilities waiting to be discovered.
The present work is devoted to realising this promise.We provide a classification scheme for spin liquids occurring as ground state ensembles in classical continuousspin Hamiltonians (ref.Fig. 1) and apply it to a number of existing and new models (ref.Table II and Table III).This enables us to understand and distinguish different kinds of CSL in a way that goes beyond simply distinguishing algebraic from short-range correlations.We identify distinct kinds of algebraic and short-range correlated CSL and zero-temperature transitions between them, and uncover simple models exhibiting previously unseen forms of spin liquid.
In this article, we develop the classification theory in some detail, with numerous examples.A shorter companion paper [86], which illustrates the main ideas in the context of a single model on the kagome lattice, accompanies this article, and may be of use to any readers who do not require such a comprehensive exposition.
The example models we construct are themselves significant as they provide simple settings in which to realize novel physics.This includes a model realizing anisotropic Gauss's laws, in which derivatives with respect to different directions enter the Gauss's law with different powers [87,88] and concomitant subdimensional excitations; a spin-liquid with a network of line-like singularities (pinch lines) in the structure factor; and a series of topological CSLs connected by zero-temperature transitions.We therefore establish the utility of our classification scheme in the construction of new models realizing interesting phenomena.
Previous works have classified highly frustrated classical spin systems via constraint counting [59], via linearization around particular ground state configurations [89] and via supersymmetric connections between models [90].In recent work by two of the present authors, the possibility of distinct types of algebraic spin liquids distinguished by topological properties was explicitly demonstrated [69].Here we present a scheme which generalizes across different kinds of spin liquid and assists in the construction of new ones.It is based on the physics of the spin liquid as a whole, rather than individual spin configurations within it and, in the case of algebraic spin liquids, unveils the connection between the microscopic model and the Gauss's law which governs the long distance physics.
The classification scheme is based on a soft spin description of the CSL state.In such a description one neglects the spin-length constraints |S| i = 1, replacing it instead with an averaged constraint ⟨S i • S i ⟩ = 1.
The soft spin approximation is known to provide a good description of CSLs for many known examples [60,64,65,69,91,92].Nevertheless, classifying CSLs according to their properties within an approximate treatment such as the soft spin approximation, may seem unsatisfactory.It is, however, in keeping with the spirit of other classification schemes such as the use of PSGs to classify QSLs.The PSG analysis is based around the properties of a mean field description of a given QSL, but remains useful because the qualitative nature of the phase is more robust than the quantitative accuracy of the mean field theory.Here, similarly, we expect our classification to correctly distinguish between CSLs, the limitations of the soft-spin description notwithstanding.
If the Hamiltonian is bilinear in spins, then one may diagonalize it in momentum space, leading to a spectrum with a band structure that carries information about the low energy spin liquid state.Our classification is based on the algebraic and topological properties of this band structure, and is schematically illustrated in Fig. 1.
The common feature of the soft spin description of CSLs is the presence of one or more flat bands at the bottom of the spectrum (Fig. 1(a)).These flat bands correspond to the extensive number of degrees of freedom which remain free in the CSL ground state.The most basic distinction we can make between CSL soft spin band structures is whether or not the flat bands at the bottom of the spectrum are separated by a gap from the higher energy bands.Spectra with (without) a gap correspond to CSLs with short-ranged (algebraically-decaying) spin correlations.
We classify CSLs without a gap for the bottom band via the algebraic properties of their band structures around the gapless points in the Brillouin Zone (BZ).In particular, a Taylor expansion of the eigenvector(s) of the dispersive band(s) which come down to meet the low energy flat band(s) at the gapless point defines an effective Gauss's law which constrains the long wavelength fluctuations of the CSL (e.g.∇•E = 0 in the ordinary Coulomb phase).The form of this Gauss's law distinguishes different kinds of such CSLs.Table II lists representatives with different generalized Gauss's laws.We name such CSLs "algebraic CSLs", due to the fact that their correlation decays algebraically, and the emergent generalized Gauss's law depends on the algebraic structure of the gapless points (Fig. 1(b)).
For the short-range correlated CSLs , the classification is based on the topology of the soft spin band structure.Depending on the symmetries present, and the number of sites in the unit cell, the bands may possess topological invariants which are insensitive to small changes to the CSL ground state constraint.These topological invari-ants can be used to distinguish different classes of such CSLs.We find that the nontrivial topology of these CSLs is generically fragile, in the sense that it can be rendered trivial by adding additional spins in the unit cell.This motivates us to introduce the term "fragile topological CSL" (FT-CSL) as a descriptor of short-range correlated CSLs (Fig. 1(c)).
By tuning the Hamiltonian, it is possible to drive zero-temperature transitions between FT-CSLs.At these transitions, algebraic CSLs emerge.We hence arrive at a landscape of CSLs where the phases are occupied by the FT-CSLs, and the phase boundaries are algebraic CSLs, as shown in Fig. 1(d).
To illustrate all these ideas we introduce a number of new models which are of autonomous interest beyond the classification scheme, in that they represent hitherto unknown types of classical spin liquids, worthy of study in their own right.A summary of different algebraic and fragile topological CSLs known in literature is presented in Table III.
Our approach to the analysis of CSLs presents a comprehensive advancement in our understanding of these frustrated systems.It reveals a landscape of classical spin liquids as fragile topological CSLs separated by the algebraic CSLs, and encompasses all CSL models to the best of our knowledge (in the soft spin setting at least).While building the classification scheme, we have established close connections between CSLs and other fields of physics and mathematics including flat bands in electronic band theory, symmetry protection and fragile topology, and homotopy theory.Our classification scheme can also be easily reversely engineered to design new CSL models with desired properties.
The article is organised as follows.The next section (Sec.II) provides a non-technical overview of our central results.The main content starts from Sec. III, which reviews a few recently discovered new models of classical spin liquids, and motivates us to pose the question of classification.In Sec.IV, we formulate the problem on a more mathematical footing, to make it amenable to the algebraic and topological treatments later.Sec.V discusses the abstract aspect of one of the two main categories of CSLs: the algebraic CSLs, followed by Sec.VI which provides a handful of examples for concrete demonstration of the physics.Following a similar structure, Sec.VII discusses the abstract aspect of the fragile topological CSLs, followed by Sec.VIII which provides a concrete example.We then briefly discuss wider applications of our classification scheme in Sec.IX and show how previously established examples of CSLs fit into our scheme.Finally, we conclude with a summary and outlook of future directions and open issues in Sec.X.

II. SKETCH OF THE MAIN RESULTS
Here, we telegraphically list our main results to guide readers through the technical details later.A selfcontained, non-technical narrative can be found in the short sister paper Ref. 86.
We study spin models in the limit of a large number of spin components N .This is effectively a "soft spin" approach, where the spin length constraint is enforced 'on average' by the central limit theorem for N → ∞.This amounts to treating each spin components as a scalar, and this has been shown to be a good approximation for many, but not all, Heisenberg candidate CSLs.CSLs in such a description tend to have an extensive degeneracy of exact ground states.
Such CSL Hamiltonians can be generally be written in what we call the constrainer form: where for a given constrainer index I, C I (R) is the sum over a local cluster of spins around the unit cell located at R. The Hamiltonians we consider are the translationally invariant sums of such squared constrainers.
For simplicitly, we will mostly work with models with one constrainer (M = 1) and N sublattice sites per unit cell, and will note where deviations from this setup affect the classification.In this case, there are N − 1 bottom flat bands at zero energy that satisfy the constraint, and one higher dispersive band that violates it.The dispersive band's eigenvector, denoted T(q), can be algebraically determined by Fourier transforming the constrainer C(R).The dispersion of the higher band is exactly ω T (q) = |T(q)| 2 .
The overall spectra encode the information of the CSLs.They can be divided into two broad categories.

Algebraic CSL:
There is one or more gap-closing point between the bottom flat band and higher dispersive bands.In this case, the CSL is an algebraic CSL, i.e. the spin correlations decay algebraically.Furthermore, the ground states can be described by a charge-free Gauss's law, determined by the Taylor expansion of T(q), where q denotes the distance in momentum space from the band touching.More specifically, if the lowest order term in the Taylor expansion is the charge-free Gauss's law is then given by where we have defined a generalized differential operator D (ma) a of order m a ≥ 1 on sublattice a.A similar picture applies for models with multiple constraints per unit cell, and hence more than one T(q), with the subtlety that one must take care of the orthogonality between different T(q) around the band touching.

Fragile Topological CSL:
The bottom flat band is fully gapped from the higher dispersive band.In this case, T(q) is a non-zero and smoothly defined vector field in the target manifold CP N −1 (if it is complex) or RP N −1 (if it is real) over the entire BZ.It can then wind around the BZ (a d-torus, T d in d dimensions) in a nontrivial manner, captured by the homotopy class In the case where there is more than one constraint per unit cell, the target manifold may be something other than CP N −1 or RP N −1 .Adiabatic changes to the Hamiltonian which retain the constrainer form and do not close the gap between the bottom flat band and the upper bands cannot change the homotopy class.The homotopy class only changes when the gap closes.That is, at the boundaries of fragile topological CSLs are algebraic CSLs.
FT-CSLs can be rendered trivial by the addition of extra degrees of freedom to the unit cell, hence our use of the term 'fragile', in keeping with the notion of fragile topology in electronic band theory [41].
In the main text we provide numerous examples to show how the abstract theory above can be applied to concrete models.

III. MOTIVATING THE CLASSIFICATION PROBLEM
In this section we motivate the question of classification by reviewing two known examples of CSLs.The two models we call the honeycomb-snowflake model proposed in Ref. 69 and the Kagome-Hexagon model, Ref. 65.These are representatives for the two different categories of CSLs.The honeycomb-snowflake model with a varying parameter hosts several algebraic CSLs that realize different generalized U (1) Gauss's laws.Correspondingly, the spin correlations decay algebraically.The Kagome-Hexagon model hosts a qualitatively different CSL -the fragile topological CSL -that does not exhibit any U (1) Gauss's laws, and has exponentially decaying spin correlations.
After reviewing the two models, we summarize their common features to extract the most general set-up for the CSL models.At the end of this section, we will be ready to establish a classification scheme that, once a specific Hamiltonian is given, can mechanically analyze the CSL physics from that Hamiltonian.
A. Honeycomb-snowflake model FIG. 2. (The honeycomb-snowflake model [69], Eq. ( 5), exhibiting a series of distinct spin liquids as the model parameter γ is varied.(a) The honeycomb lattice, composed of two sublattices, colored red and blue respectively.(b) The constrainer defining the ground states of the model, applied to each hexagonal plaquette in the lattice.The sum of spins on each hexagon (1 to 6) plus the coefficient γ multiplied by the sum of spins linked to the exterior of the hexagon (1 ′ to 6 ′ ) must vanish on every hexagon (Eqs.( 6)-( 7)).
The honeycomb-snowflake model proposed in Ref. 69 serves to demonstrate how a series of distinct algebraic CSLs can be accessed by varying local constraints on a classical spin system.Its Hamiltonian is defined as a squared sum of Heisenberg spins around the hexagonal plaquettes on the Honeycomb lattice [Fig.2]: R∈all hexagons α=x,y,z The sum of R is taken over all hexagonal plaquettes of the lattice or, equivalently, over all unit cells.The sum of α = x, y, z is taken over all three spin components.The terms C γ HS,α (R) defined on the hexagons are weighted sums of spins around each "snowflake" shown in Fig. 2(b): The first sum in Eq. ( 6) is over spins on the hexagon at R (sites in Fig. 2(b) labelled 1 to 6) and the second is over neighboring spins connected to exterior of the hexagon (sites in Fig. 2(b) labelled 1 ′ to 6 ′ ).γ is a dimensionless parameter which we use to tune the model.Ground states of Eq. ( 5) satisfy the constraint: The case γ = 0 corresponds to the model of Ref. [64].
Let us now outline a description of the honeycombsnowflake model, equivalent to that in [69], based on the gap closings in the spectrum of the Fourier transformed Hamiltonian.First, we observe that the Hamiltonian is identical for the three components α = x, y, z.If we relax the spin norm constraint, |S i | = 1, and treat it only on average (⟨S i • S i ⟩ = 1), the spin components can be thought of as essentially indenpendent scalar variables.This step can be justified more formally taking the limit of a large number of spin components, N .The theory in which the spin norm is fixed only on average then corresponds to the leading order of a 1/N expansion.This approach has been for example successful, even quantitatively, in describing pyrochlore spin liquids with O(3) Heisenberg (N = 3), and even Ising spins (N = 1) [60], and has been widely used in the treatment of spin liquids since its introduction to the field in Ref. 63.In the remainder of the paper, we work within this large-N picture.This allows us to build our classification scheme, and in this sense we are working in the same spirit as other classification schemes in Condensed Matter Physics which are also derived from mean-field or non-interacting theories; with the expectation that the classification labels are robust even when the underlying approximate theory is not quantitatively accurate.Exceptions to this can, and do, occur however -such as in the case of the O(3) kagome Heisenberg model.While a large-N picture predicts a spin liquid, the order-by-disorder effects drive the O(3) system into an ordered phase at very low temperature [93][94][95][96][97].The approach we present here thus provides a tool for classifying CSLs but does not prove the stability of a CSL in any given hard-spin model, which is a task that generally requires simulations.
Working within the large-N theory, we can drop the component index label α and regard each spin S i as now a scalar instead of a vector.
Taking the Fourier transform of Eq. ( 5) results in a Hamiltonian written as a 2 × 2 interaction matrix J(q), Sa (−q)J γ ab (q) Sb (q) , where a, b index the two translationally inequivalent sublattices of the honeycomb lattice, and S(q) = ( S1 (q), S2 (q)) is the lattice Fourier transform of the spin fields on the sublattice sites 1, 2. The interaction matrix J γ (q) depends on the parameter γ, and can be computed straightforwardly.Its explicit form is lengthy and not of importance for now, but can be found in Eq. ( 108) in Sec.VI B when we revisit this model.Diagonalizing J γ (q) yields a 2-band spectrum, in which the lower band is flat at energy ω = 0, and the upper band is dispersive, its dispersion denoted as ω(q).The gap between the flat and dispersive bands closes at multiple points in the Brillouin zone [Fig.12].The Hamiltonian can then be represented as ω(q) Sa (−q) Ta (q) T * b (q) Sb (q) , ( where ω(q) is the eigenvalue of the dispersive (upper) band and T(q) is the corresponding normalized eigenvector of the top band.
The upper eigenvector can be used to give a momentum space description of the ground state constraints, Eq. (7).Any Fourier-transformed spin configuration on the two sublattices a = 1, 2 obeying the condition a=1,2 ω(q) T * a (q) Sa (q) ≡ a=1,2 T * a (q) Sa (q) = 0 ∀ q , (10) is a ground state.This is the momentum space representation of Eq. (7).
Eq. ( 10) can be seen as an orthogonality condition between the vector of sublattice Fourier transforms of the spin configuration S(q) and the upper band eigenvector T(q).(The upper band eigenvector is thus equivalent to the constraint vector L(q) introduced in Ref. 69).
The ground state phase diagram of the honeycombsnowflake model is shown in Fig. 3. Three distinct algebraic CSLs emerge as γ is varied (the CSLs at large negative γ and large positive γ are equivalent, as may be inferred from Eqs. ( 5)-( 6)).In Ref. 69, the distinction between these CSLs was understood in terms of topological defects in T(q).
It was found that the CSLs with a pinch point (singular pattern of the structure factor at the K point) [81,98,99] host a spin liquid described by the Gauss' law of a Maxwell U(1) gauge theory: Here, E = (E x , E y ) is an emergent vector electric field degree of freedom (DOF).At γ = 1/2, four of the pinch points merge at the K point, forming a 4-fold pinch point (4FPP) [100,101], and a more exotic Gauss's law describing the system in terms of a rank-2 electric field with a scalar charge [102] was found: We will come back to the emergence of different Gauss's laws in Sec.VI B. Finally, let us explain the plots in Fig. 3, and the similar plots which appear for other models throughout the paper.Fig. 3 shows the spectrum of J γ (q).Additionally, on each band i, we have also plotted the spin correlations defined as where v i (q) = (v 1 i (q), v 2 i (q)) is the eigenvector of the band i.
In the T = 0 limit of a large-N approximation the equal time structure factor is the sum of the structure factors S(ω, q) i over the flat bands only [61]: FIG. 3. Phase diagram of the honeycomb-snowflake model (Eq.( 5)) as a function of γ, showing a series of algebraically correlated CSLs.Transitions between spin liquids occur either by creation/annihilation of pinch points in the spectrum (γ = 1/3) or by merging of them (γ = 1/2) leading to a higher-rank Coulomb liquid with multi-fold pinch points.The phase diagram is based on Ref. 69. (a-e): band dispersion ω(q) (upper panels) and structure factor S(q) (lower panels) of the honeycomb-snowflake model with varying γ.There is always a gap closing at the K point of the Brillouin zone and additional ones are created or annihilated in pairs as γ is varied, giving rise to topological transitions between distinct CSLs.
In particular, the spin structure factor measured in inelastic neutron scattering contains valuable information about these pinch points and can be used to experimentally determine the nature of the CSL.The dispersion ω(q) vanishes at multiple points in the Brillouin zone (BZ).At these points, the upper band eigenvector T(q) (and hence Eq. ( 10)) is not uniquely defined and there are corresponding singularities in S(q).These singularities in the structure factor give rise to an algebraic form of the spin correlations when Fourier transformed back into real space, and also dictate the Gauss's law constraining the spin fluctuation of the ground states.The presence of non-trivial gap closings is, therefore, an essential part of the physics of these CSLs.

B. Kagome-Hexagon model
We now discuss the kagome-hexagon model [65] as an example of fragile topological CSL with short-ranged correlations at low temperature.Its Hamiltonian is defined on the kagome lattice: R∈all hexagons α=x,y,z where the sum of R runs over hexagonal plaquettes on the kagome lattice (indicated in Fig. 4(a,b)), or equivalently the centers of all unit cells.C α (R) is the sum of the six spins around each hexagon as labeled in Fig. 4(a,b): Ground states are hence defined by the constraint: on every hexagonal plaquette.Again, the Hamiltonian is the same for the three components α = x, y, z, and within the large-N approximation we treat this as three copies of a theory in which the spins are independent scalars.Multiplying out Eq. ( 15) we can rewrite the Hamiltonian in terms of bilinear exchange interactions.These interactions couple first, second and third neighbor spins across the hexagon with equal strength.15)-( 17)).(c) Spectrum ω(q) that arises from diagonalizing the Hamiltonian (Eq.( 15)) in momentum space.There are two degenerate flat bands at the bottom of the spectrum and a dispersive upper band with no band touchings between the upper and lower bands.(d) Spin structure factor showing an absence of singularities.
Taking the lattice Fourier transform of the interactions results in a 3 × 3 interaction matrix J(q), since there are three sites per unit cell.Diagonalizing J(q) yields a spectrum with three bands, of which the lowest two are flat and degenerate, with the upper band being dispersive (Fig. 4

(c)).
There are no band touchings between the upper band and the two flat bands at any point in the Brillouin zone.Accordingly, the real space correlations remain short ranged with a correlation length on the order of the nearest-neighbor distance at T = 0. Also, the ground state fluctuations are not described by any effective Gauss's law.
The CSL state of this model seems to be qualitatively distinct from a trivial paramagnet, as evidenced by the fractionalization of "orphan" spins around a cluster of introduced vacancies [65].But this raises questions: are there different types of non-trivial, short range correlated CSLs?If so, how do we distinguish these and are they separated by sharp transitions?These are some of the core questions this work will address.

C. The question of classification
Having examined some sample models, we can now sharpen the question of classification.The common feature of the CSLs is that they are described by the type of Hamiltonian H = R [C(R)] 2 , where C is a sum over a local cluster of spins.The Hamiltonian can be diagonlized in momentum space, characterized by a matrix J(q).Its spectrum has one or several flat band(s) at the bottom, below one or several bands that are generally dispersive.
However, depending on the structure of the spectrum, different CSLs can have very distinct properties.It is thus important to understand the mechanism that leads to such distinction, and provide a classification scheme that puts all CSLs in their place.
The first fundamental difference between CSL mod-els can be seen by comparing the honeycomb snowflake model with the Kagome-Hexagon model.Some CSLs, such as in the former model, have gap closing points between the bottom flat band(s) and the higher band(s), while others do not.This leads to the most crucial differences between the two categories of CSLs.The former has an emergent generalized Gauss's law describing the ground state fluctuations, and algebraicallydecaying spin correlations.The latter exhibits no emergent Gauss's law, and its spin correlations decay exponentially.
Within each category, we still need to make finer distinctions between CSLs.For the first category, which we call algebraic CSLs, the main question is how many and what kind of generalized Gauss's laws describe the ground state fluctuations, and how a particular type of Gauss's law appears.We will show that the number and structure at the gap-closing points determines this, and also explains the transitions between different algebraic CSLs.
For the second category, which we call fragile topological CSLs, there is no generalized Gauss's law.It is then natural to ask what can distinguish the different members of this group.We will show that the homotopy class of the eigenvector, defined as a map over the BZ, see Eq. ( 4), is a good topological quantity for the classification of the fragile topological CSLs.

IV. SETTING UP THE CLASSIFICATION FORMALISM A. Constrainer Hamiltonian and its spectrum
Let us now define the CSLs in sufficiently general and accurate terms for the development of a robust classification scheme.
First, we work with spins in the large N limit, or equivalently with soft spins.That is, we treat every spin com-ponent S α i as a real number, S α i ∈ R. We ignore nonlinear constraints that may apply to a real system.The common types of non-linear constraints include those on Ising and Potts variables that take a finite range of discrete values; or on classical Heisenberg spins that are three-dimensional vectors of unit length (S 2 = 1).Very often, the soft spin treatment provides a good approximation to real spins, and can correctly capture the physics of the actual CSLs.Exceptions do exist, and we will discuss this point in later sections of this Article.
For Heisenberg models, each vector spin has three DOFs S x , S y , S z .However, because they decouple from each other in the soft-spin limit, and are described by the same Hamiltonian individually, we can just analyze one copy of them.From now on, we therefore treat spin components as independent scalars, and collectively denote them as S a (R) where a = 1, . . ., N are the number of DOFs in a unit cell labelled by its position R.
Second, we work with bilinear Hamiltonians with finite range of interactions, which is natural for most physical systems.We specifically investigate the CSLs where the dimension of the ground state manifold grows linearly with the system volume.We thereby exclude spiral spin liquids [103][104][105][106][107], which have subextensive degeneracies.An equivalent statement is that we study the systems where the spectrum of the Hamiltonian has one or more flat bands at its bottom.
Such CSLs can be written in forms of what we call constrainer Hamiltonians.Such Hamiltonians are written in real space as where C I (R) is a linear combination of spins around the unit cell centered at R (not necessarily restricted to nearest neighbors).Different spins can have different realvalued coefficients (weights) in this sum.The index I in the summation runs over all constrainers (there could be more than one) in a unit cell.Here we denote the number of constraints per unit cell by M , and the number of spin sites in the unit cell by N .In real space, the ground states of a classical spin liquid are the spin configurations s.t.C I (R) = 0 at all unit cells and for all I's, hence the name constrainer.Given N DOFs in a unit cell and N > M , then generally such ground states exist, at least within the large-N approximation, because there are more DOFs than constraints.
The constrainer Hamiltonian formalism includes all the canonical classical Heisenberg spin liquids.For such models, one can always add a term i (S i ) 2 or i (S z i ) 2 with the correct coefficient to turn the Hamiltonians into the constrainer formalism.This added term does not affect the physics because the spin length is fixed in the hard-spin Heisenberg model.
Let us now write down the constrainer Hamiltonian in a more explicit form by specifying C I (R).First, it is convenient to encode a given constrainer C I (0) for the unit cell at the origin (R = 0) in a N −component vector encoding the information of how different sublattice sites are summed in the constrainer.It is written as Here, r is a variable that we use to visit all sites on the lattice to see if a spin at the location r is involved in C I (0) (it is if and only if a b,j appears in C I (0, r)).The first component [C I (0, r)] 1 records information of all the first sub-lattice sites in different unit cells that are involved in C I (0).Their locations are at a 1,j 's relative to the center at 0, where a 1,j , pointing to nearby unit cells, is always an integer multiples of lattice vectors plus a constant shift to the center at 0. The coefficients for different spins summed in C I (0) are c 1,j 's.Similarly, the b th component [C I (0, r)] b records the information of how the b th sublattice sites are summed in C I (0).Hence, given C I (0, r), we have the complete information of how the constrainer C I (0) is defined.
For the constrainer in a unit cell at a general location R, we need to perform a translation on C I (0, r) to get The real space Hamiltonian is written explicitly as Here, S(r) = (S 1 , . . ., S N )(r) is the vector array formed of the N sublattice sites.For example, S b (r) is the b-th sublattice site at location r.The term r S(r) is the explicit form of the constrainer C I (R) as shown in Eqs.(6,16).With C I (R, r) given, we now do not need to rely on pictorial description of the constrainers.Instead, we now have their algebraic description ready for mathematical treatment in what follows.Now let us diagonalize the Hamiltonian in momentum space.A billinear Hamiltonian can be diagonalized in momentum space as Sa (−q)J ab (q) Sb (q) Here, Sa is the Fourier-transformed spin field S a , and a = 1, . . ., N labels the sub-lattice sites.J is the N -by-N matrix of the interactions.For constrainer Hamiltonians, there is a simple expression for J based on C I .Each constrainer C I can be Fourier transformed into momentum space as the FT-constrainer T I (q), J ab (q) explicitly reads Note that using the constrainer at either 0 or a general unit cell position R to define FT-constrainer T I (q) does not affect J(q), since it only adds an overall phase to T I (q) that is cancelled in J(q).
In momentum space, we can examine the spectrum of J(q) (we will now slightly abuse the notation and refer to both H and J(q) as the Hamiltonian).Given M constrainers in the Hamiltonian, there will generally be M upper bands and N − M bottom flat bands.The upper bands may touch the bottom flat ones at some special points (or in some cases along special lines or planes).The higher bands' eigenvectors are those in the space spanned by all T I (q)'s, but not necessarily T I (q)'s themselves: note that two different constrainers T I (q) and T J (q) are not required to be orthogonal to each other.The bottom degenerate flat bands' eigenvectors are those orthogonal to all T I (q)'s.
The information of the ground states of the CSL is encoded in the bottom bands and their eigenvectors.Equivalently, one can also access such information from the higher bands and their eigenvectors T I (q), since the two sets of eigenvectors are orthogonal to each other and span the full N -dimensional vector space.It is often easier to look at the higher bands since all T I (q)'s are known explicitly from the definition of the constrainer C I 's.
Let us now analyze the structure of the ground states.First, we note that they span a linear subspace in the space of all spin configurations, and the ground state fluctuations span an isomorphic linear space.Starting from a ground state that satisfies C I (R) = 0 for all R's and I's, we can then consider a fluctuation that keeps C I (R) = 0. Note that the C I 's are linear in the spin variables, thus, at the level of the soft spin approximation, any ground state and any such fluctuation can be added linearly with the system remaining in a ground state.Mathematically speaking, all the ground states span a linear (vector) space, so the ground states manifold and the manifold of fluctuations between ground states are isomorphic.In more physical terms, we can start with any initial ground state, and then every other ground state is bijectively mapped to a fluctuation from the initial ground state to it.
Just like the constrainers describe energetically costly spin configurations in real space, and their Fourier transforms describe the higher bands in the spectrum, their counterparts describe the ground state fluctuations.Let us consider the local fluctuations that satisfies the C I (R) = 0 condition for all R's and I's.Since the bottom band is (N −M )-fold degenerate, we know that there should be (N −M ) such linearly-independent local fluctuations.We name these fluctuators and abstractly denote them as F I (R), where I = 1, . . ., N − M .We express each fluctuator as an N −component operator acting linearly in the spin vector space, just as we did with the constrainers C I 's in Eq. ( 19), and denote them as F I (R, r): The components in the fluctuator F I (R, r) describe quantitatively the local spin fluctuations that keep all constrainers zero, i.e., the flucturator is a zero eigenmode of the Hamiltonian.
Fluctuators and constrainers are orthogonal: The FT-fluctuator, defined as the Fourier transform to the momentum space is then orthogonal to all the FT-constrainers T J (q).The FT-fluctuators are exactly the eigenvectors spanning the (N − M ) degenerate bottom flat bands.
The sample models studied in this paper have only one constrainer (M = 1), so we can drop the index I: In this way, the physics can be clearly demonstrated without too much notational complication.Correspondingly, the Hamiltonian matrix has N −1 flat bands with eigenvalue zero and 1 dispersive higher band.T(q) as the only FT-constrainer is also the unnormalized eigenvector of the higher band.The higher band dispersion is The top band may or may not touch the bottom bands, depending of the specific form of the constrainer C(R, r) and its Fourier transform T(q).Since there is only one top band but several bottom bands, it is easier to analyze the top band rather than the bottom ones.The physics is easily generalizable to the cases with multiple higher bands.
Depending on whether the top band touches the bottom bands, the CSL falls into one of two broad categories.
The algebraic CSLs have band-touching points, and are controlled by the physics around those points; they have algebraically-decaying correlations described by emergent, generalized U(1) Gauss's laws.
The fragile topological CSLs have no band-touching points, and the correlations in the bulk are short-ranged.However, as we shall demonstrate below, they have quantized topological properties that cannot be changed without closing the gap between the top and bottom flat bands.All the topological information is encoded in the FT-constrainer T(q).After introducing several mathematical tools, we will demonstrate in detail how to extract the information about the algebraic and fragile topological CSLs from the Hamiltonian.

B. Tools from flat band theory
Since our analysis focuses on flat bands, known results from flat band theory (for fermionic/bosonic hopping models) can be applied here.In this subsection we review these results, with a view to applying them later.
The key properties for the CSLs are encoded in the flat bands at the bottom of the spectrum of the Hamiltonian matrix J(q).In the context of classical spin systems, the bottom bands are related to the fluctuations between ground states, as discussed in Sec.IV A .The real space local fluctuators F I (R, r), or equivalently the momentum space FT-fluctuators B I (q), describe these fluctuations.
The key to the physics of a flat band in a free hopping model is that a flat band in momentum space corresponds to a compact local state (CLS, not to be confused with classical spin liquid, CSL) in real space.The compact local state is an exact eigenstate of the Hamiltonian, and is only supported on a finite, local region of the lattice.Their existence is proven in Appendix A of Ref. 108.Such a locally supported state usually does not exist for a dispersive band.Compact local states in real space, and the flat band in momentum space, are two facets of the same physics.For a rigorous proof of this statement, see Sec.II.A of Ref. 108.
The connection to CSLs is the following: the compact local state in the hopping model corresponds to the fluctuator in CSL.Let us use nearest-neighbor-hopping kagome model, as an example.Its CSL version is the nearest neighbor AFM kagome model, which is a classical spin liquid in the large-N description (although order-by-disorder at very low temperatures cuts off the spin liquid behaviour for O(3) Heisenberg spins [93][94][95][96][97]).A more detailed analysis of the kagome model will be presented in Sec.VI A 2.
Here we state a few basic facts of it.Its Hamiltonian is Given the hopping Hamiltonian Eq. ( 31), one can find by inspection the compact local state of the model.The wavefunction of the compact local state at location R can be generically encoded in an N -component fluctuator vector F r,R via the relation Here, the a th component of F R,r encodes the information of the a th sub-lattice site's contribution to the compact local state.And |r, a⟩ denotes an electron occupying the sublattice site a at unit cell r.In the case of kagome model (Eq.( 36)), N = 3, and F 0,r is The corresponding compact local state wavefunction is illustrated in Fig. 5.One can apply the hopping Hamiltonian to it and find that the hopping amplitude of the compact local state to any other site is exactly zero.
We can now illustrate the connection between the compact local state in Eq. ( 33) and the bottom band eigenvector of the CSL model in Eq. (32).Indeed, the compact local state in Fig. 5 has the property, when reformulated in the language of spin components, that the sum of spins on each triangle remains 0, as expected from Eq. ( 32).More formally, Fourier transforming the fluctuator (34) into momentum space yields On the other hand, diagonalizing both Hamiltonians in momentum space, we obtain (up to adding an additional TABLE I. Connection between the physics of flat band theory and classical spin liquid, using the language of compact local states (CLS) and non-local loop states (NLS) [108].
, and a 3 = −a 1 −a 2 encode the lattice geometry.We can directly confirm the flat band is at ω = 0 with eigenvector where is the nomalization factor.This is exactly the normalized FT-fluctuator in Eq. (35).
We have thus established that the compact local state formulation of the flat band hopping Hamiltonian is related, via the Fourier tranform, to the momentum-space eigenvector (a.k.a. the fluctuator) of the ground state spin configuration in a CSL.Having established this connection, we can now translate known properties of compact local states into the language of spin liquids.
As mentioned before, the leading-order criterion for classification is whether there is band touching between the flat bands and upper bands or not.This is also one of the main topics in the study of compact local states.Depending on the hopping model, there can be three scenarios: no band touching, non-singular band touching, and singular band touching.Each of these cases is reflected in the structure of the corresponding compact local states.

Non-singular band touching
Roughly speaking, a non-singular band touching is an "accidental" band touching that does not qualitatively affect the physics of the flat band.More precisely, it can be defined in terms of the CLSs being linearly independent of the eigenvectors of the non-flat band.The  31)), which can also be interpreted as the local and loop fluctuators in the classical spin liquid model (Eq.( 32)).One can check that the hopping amplitude from these states to any other site is zero.The CSLs are not linearly independant: adding all of them on the entire lattice yields zero.(b) Spectrum of the Hamiltonian in Eq. (36).
simplest example of this is two completely decoupled systems I and II, each with its own bands.Obviously, if a band in system I touches the flat band in system II, there is nothing special happening at the band-touching point, and the band touching can be lifted trivially.Another example is that the vicinity of the band touching in a two band system can be written as p(k x , k y )(σ z + σ 0 ), where p(k x , k y ) vanishes at the band touching.Since the matrices σ z and σ 0 commute, the two modes can be trivially separated by shifting the dispersive band upwards via addition of a term E 0 σ z .Such non-singular band touching can thus be smoothly deformed to a gapped spectrum.
Let us first discuss the physics of flat band with no band touching or with non-singular band touchings only.In this case, the eigenvector B(q) of the bottom band is well-defined globally so that a vector bundle associated with the flat band exists globally -this is known as a trivial vector bundle.(The reader may be familiar with nontrivial (complex) vector bundles, which can possess nontrivial Chern number.The perfectly flat bands resulting from the constrainer formulation of the CSL can be shown to have zero Chern number if they are separated by the gap, see Section VII below and Ref. 118.)To make B(q) well-defined without any singularity requires |B(q)| > 0 in the entire BZ.This is exactly the condition that the bottom band is separated by a gap from the dispersive higher bands.
In real space, that means the L x L y compact local states generated by applying lattice translations to a single CLS (assuming the lattice has L x L y unit cells) are all linearly independent, so they span the entire flat band [108], encoding the L x L y states on this band exactly.The same applies to the CSL models: if the total L x L y fluctuators F R,r on different unit cells are linearly independent, then the corresponding FT-fluctuator B(q) is non-vanishing everywhere in the BZ, and there is no band-touching between the bottom bands and the top ones.The Kagome-Hexagon model (Eq.( 15)) is an example of such a system (with a slight complication that the flat bands are twofold degenerate).

Singular band touching
Let us move on to the case of singular band-touching between the bottom flat bands and higher bands.In this case, the bottom band eigenvector B(q) vanishes at certain q's, which are the band touching points.A single band accounts for L x L y states of the Hamiltonian, so a flat band with band touching points should account for L x L y + n states, where the additional n states come from the degeneracy at the band touching points.The exact value of n depends on the type of band touching ponts.Therefore, the L x L y compact local states (related by spatial translations) are not enough to account for all the L x L y + n states on the flat band.Moreover, in the presence of singular band touchings, it can be shown (see e.g.Ref. 108) that the L x L y compact local states are not linearly independent.
Where are the missing states?It turns out that there are new, non-local loop (or other topological) states (NLS), which are eigenstates of the Hamiltonian.They are new in the sense of being linearly independent from the compact local states.They account for the states on the flat band which are missing due to the linear dependence of the compact local states [109], as well as additional n states from the band touching points.Singular band touching, linear-dependence of compact local states, and the existence of nontrivial loop states are different facets of the same physics.
The physics can be translated to CSLs too.In this context, the local fluctuators F(R, r) are not linearly independent, and B(q) becomes zero at the singular band touching points.There are loop fluctuators F loop accounting for the degeneracy at the band touching points.The consequence of them -emergence of the generalized U(1) structure -will be analyzed in detail in the next section.
There is one band touching per BZ, and the total number of zero energy states is L x L y + 1.We can see that the eigenvector B(q) (Eq.( 35)) becomes zero at q = 0, where a singularity exists, i.e., B(q) is not smooth there.This is in contrast to the non-singular band touching point, where B(q) can be written down smoothly.In real space, that means an equal weighted sum (i.e.phase distribution of q = 0) of all the L x L y compact local states (Fig. 5) vanishes, meaning they are linearly dependent [119,120].Removing any one of them results in L x L y − 1 linearly independent states.In addition to the compact local states, there are two non-local loop states supported on winding loops on the lattice (Fig. 5) that are also eigenstates of the Hamiltonian.They cannot be constructed from the compact local states, so we have in total L x L y + 1 states at the energy of the flat band.They account for all the states on the flat band and at the point of band touching.
Finally, we comment that a complete set of compact local states accounting for all states on the flat band can always be found in 1D systems.Therefore all flat bands in 1D system have no band touching, or at most nonsingular ones [121,122].We will therefore concentrate on 2D and 3D examples in this article.

V. ALGEBRAIC CSL CLASSIFICATION: EMERGENT GAUSS'S LAWS
The common feature of algebraic CSLs is that the gap between bottom flat band(s) and higher band(s) closes at some points in the momentum space in a singular manner, or, in the words of flat band theory, these singular band touching points determine the class of the algebraic CSL.By examining the eigenvector configuration or, equivalently, the effective Hamiltonian near a band crossing point, one can derive the generalized U(1) Gauss's law emerging there.The ground state fluctuations are essentially effective electric field fluctuations that obey a charge-free condition, in which the charge is defined via the generalized U(1) Gauss's law.The statements above are already well-understood for conventional U(1) spin liquids like the pyrochlore (and kagome N ≥ 4) Heisenberg models.In this paper we generalize it to other types of U(1) spin liquids with a simple algorithm to identify the Gauss's law.

A. U(1) structure of the ground state manifold
We first show that with singular band touching points, the linear space of ground states has a U(1) structure.
As we have established in previous section, in this case, there are local fluctuators encoded in F I (R), and also loop fluctuators we denote abstractly as F loop 1 , F loop 2 etc that are linearly independent from the local ones.Together they account for all states on the bottom flat bands and the additional states at the band touching points.
Hence the ground states can be divided into equivalence classes in the following sense: two ground states are equivalent if and only if there are some local fluctuators that take one to the other.Then, applying a loop fluctuator to a ground state takes it to another equivalence class.Note that, each loop fluctuator comes with a real coefficient c.The equivalence classes hence have an uncompactified U (1) (or R) structure.The loop fluctuators play the role similar to logical operator in topological orders, taking the ground state from one equivalence class (a.k.a.superselection sector) to another.This is schematically shown in Fig. 6.Now the question is: how to describe the U (1) structure?It turns out that if the band-touching manifold is a point (or a few points), then associated with each point one can derive a generalized Gauss's law by examining the eigenvector structure around the point.This will be the central result for the algebraic CSL classification.
In more exotic cases, the band-touching manifold is not a point (or a few points) but a higher dimensional object (curves, membranes etc.).Then it is no longer possible to write down the long-wavelength physics as an expansion around a point and obtain a Gauss's law to capture all the physics -because there are infinitely many gapless points elsewhere.In what follows we will mostly focus on the former case of isolated touching point(s).

B. Generalized Gauss's laws and their physics
While the Maxwell U(1) gauge theory and its reincarnation in classical spin ice are well known [61,62,75,[123][124][125], the concept and consequences of generalized U(1) gauge theories may be an unfamiliar topic to some readers.In this section we introduce the electrostatics of these new theories, since many of the algebraic CSLs are described in this language.Given our focus on classical spin liquids, we focus here on the classical electrostatic sector of the generalized Maxwell theory, without the magnetic fields which would introduce quantum dynamics.The Gauss's law of Maxwell U(1) theory is written as The spin liquid ground states are described by an electrostatic theory requiring the charge-free condition to be satisfied everywhere on the lattice.As the simplest Lorentz invariant gauge theory, the Maxwell U(1) gauge theory describes one of the fundamental forces of the universe as well as the emergent behavior of various many-body systems.Obviously, electric field fluctuations obeying the charge-free condition preserve the net Noether charge of the system A difference between condensed matter systems and the universe is that the Lorentz symmetry, including continuous rotational symmetry of space, can be broken in the former cases.This means the emergent theories describing solid-state systems need not have Lorentz invariance.Instead, only a lower set of symmetries (e.g.discrete rotational symmetry of the lattice or even less) need to be satisfied.Applying this principle to the CSLs, means that one can write down generalized U(1) gauge theories and their Gauss's laws that do not necessarily respect Lorentz/rotational symmetry.Some of the preeminent examples in recent years are the rank-2 symmetric U(1) gauge theories [102,126].
Here we briefly review the so-called scalar charged case [102].The theory respects rotational symmetry of space but not the Lorentz symmetry.Its electric field is a rank-2 symmetric tensor E αβ , which can be chosen to be traceless or not.Its (scalar) charge is defined as One exotic consequence is the conservation of charge dipole or higher multipoples.For the example given above, the total electric dipole in the γ spatial direction where dΣ α denotes an integral over the boundary surface normal to the component α.This implies that the total dipole moment is entirely determined by the value of the fields at the boundary of the system, which further implies that it cannot be changed by any local rearrangement of the electric field in the bulk.Thus any local dynamics must conserve the electric dipole, with the consequence that isolated charges cannot move, in contrast to Maxwell U(1) gauge theory.Such immobile charges are dubbed fractons, which have received much theoretical attention in the past decade (see e.g.Ref. 127 for review and references therein).We can take a further step in the generalization [128].We need two pieces of data to define a generalized electromagnetism: (1) the electric field and (2) the Gauss's laws that define the charges.The electric field does not need to be in the form of vector or tensor, since we do not enforce the rotational symmetry in the first place.Instead, we just label different components of the electric field as E i , where i = 1, 2, . . ., n E .Correspondingly, the charges do not need to be a scalar, vector, or tensor.Instead it can have several components labeled as ρ j where j = 1, 2, . . ., n c .Each component is defined via the Gauss's law as Here, D j i 's are linear differential operators.In the case of Maxwell electromagnetism, Gauss's law is explicitly written in Eq. (39), and in the case of rank-2 U(1) symmetry gauge theory it is written in Eq. (42).One can also write down any other choice of D j i to define a new U(1) electromagnetism.
For a generalized gauge theory, the conserved quantities are for any set of functions {f 1 , f 2 , . . ., f nc } that satisfy Here, Dj i is a linear differential operator, related to D j i by multiplying every term in D j i that has n derivatives by (−1) n .It is obvious that total charge conservation, i.e. f j = constant holds for any generalized Gauss's law.But depending on the form of D j i , there can be other sets of {f j } that satisfy Eq. (46).For instance, choosing f j = r j would correspond to the dipole moment conservation in Eq. (42).The above generalization encapsulates new conservation laws in the form of charge dipoles, multipoles, or combinations thereof.Like the rank-2 symmetric U(1) gauge theory, such multipole conservation laws lead to immobility of isolated charge excitations, which are fractons.
Eqs. ( 44)-( 45) complete the definition of electrostatics (i.e. the classical sector) of the generalized U(1) gauge theory.We will show that the algebraic CSLs are described by the low energy effective theory, written here in the Hamiltonian form where E i emerges from the spin degrees of freedom (see section V C for the detailed derivation).The ground state fluctuations are then described by a generalized Gauss's law and the requirement that all charges vanish Given the definition of electric field and charge (Eq.( 44)), it is also straightforward to write down the gauge transformations (more accurately speaking, gauge redundancy) and construct the magnetic field as objects invariant under these gauge transformations.The synthetic magnetic field encodes the fluctuations within the classical manifold of degenerate states and is necessary to describe the quantum spin liquid that originate from its 'parent' CSL, see e.g. the well-known U(1) description of quantum spin ices [75].This completes the construction of electromagnetism of the generalized U(1) gauge theory; interested readers can refer to Refs.[128,129] for more details.

C. Extracting Gauss's laws: one constrainer models
The generalized Gauss's laws introduced above provide a description of the ground state fluctuations in terms of the generalized charge-free condition in the corresponding U(1) theory.Hence, the Gauss's law distinguishes different algebraic CSLs.
We will describe the general mathematical recipe to determine the Gauss's law in this section, and then apply it to concrete examples in Sec.VI.Since the only terms in the Hamiltonian are the constrainers, they must dictate the emergent Gauss's law.In momentum space, FT-constrainers (i.e. the eigenvectors of the higher bands) describe the energetically costly spin configurations.Upon the inverse Fourier transform into real space, these become the (generalized) derivatives D j i E i (see Eq. ( 44)) in the long-wave length limit, which turn out to be precisely the formulation of Gauss's law.
In real space, the Hamiltonian is given by the constrainer form Eq. (18).To lighten the notation, we assume one constrainer in what follows: In Sec.IV A, we have analyzed the mathematical detail of this type of Hamiltonians.It has one dispersive top band and N − 1 bottom flat bands, where N is the number of sub-lattice sites in a unit cell.The Fourier transformed constrainer (FT-constrainer) T(q) has N components.
The Hamiltonian in momentum space is then represented by an N × N matrix in Eq. ( 29) The eigenvector of the top band is T(q), and its eigenvalue (dispersion) is ω top (q) = |T(q)| 2 .The N − 1 bottom bands are at energy 0, whose eigenvectors are those orthogonal to T(q) ≡ T(q)/|T(q)|.
Since we are studying the cases in which singular bandtouching happens, there must be one (or more) wavevector q 0 where the dispersive band has zero eigenvalue: ω top (q 0 ) = |T(q 0 )| 2 = 0.At this point, all components of T(q 0 ) are identically zero.This is reflecting the singular nature of the band-touching point: due to the nonsmoothness of the eigenvector configuration around the singular gap-closing point, the only way to write it down continuously is to have T(q 0 ) = 0.If the band-touching point is non-singular, then such a requirement does not apply, and one can choose T(q) in such a way as to be smooth and non-vanishing in the neighborhood of q 0 .Expanding T around q 0 for small k = q − q 0 , we get Note that by construction of the FT-constrainer T(q) (Eq.( 23)), q x , q y always appear in exponential forms as exp (iq • a i,j ), we can then expand each component Ta (k) as a polynomial of ik x , ik y which satisfies Ta (0, 0) = 0, for a = 1, . . ., N .That is, there is no constant term in the polynomial, so the leading term must have finite powers of k x , k y .The emergent Gauss's law is encoded in the algebraic form of the FT-constrainer T(k x , k y ).Note that T(k x , k y ) lives on the top band, so it describes the spin configurations that cost energy.That is, it encodes the generalized electric charge in terms of the spins S 1 , S 2 , . . ., S N .
Before describing the most general scenario, let us look at a simple example.Consider a system with N = 2 degrees of freedom per unit cell, and Then the bottom band eigenvector ( S1 , S2 ) satisfies Identifying the Fourier modes of the emergent electric field with the spins: Ẽ ≡ ( Ẽ1 , Ẽ2 ) = ( S1 , S2 ), this condition ik • Ẽ(k) = 0 is exactly the Fourier transformed conventional U (1) charge-free constraint in real space using ik x → ∂ x , ik y → ∂ y .The long-wavelengh effective Hamiltonian is then formulated as in real space.This imposes exactly the two dimensional electrostatics of the Maxwell U(1) gauge theory, i.e., the electric field configuration has to obey charge-free condition at low energy.Now let us formulate the general description.For each polynomial Ta (k x , k y ), we only need to keep the leadingorder terms in ik x , ik y , since higher-order terms become negligibly small for sufficiently small k x , k y .Suppose for a component T * a (k), the leading order term is of power m a ≥ 1, then it takes the general form The emergent Gauss's law in momentum space is then written as If the expansion is around a general wavevector point q 0 , then the c ij 's can be complex.The Fourier mode of spin field S i is also complex.It is reconciled with the fact that the spins are real scalars by the constraint that This guarantees the Fourier mode expansion of the real scalar field is also real after taking into consideration of both q 0 and −q 0 .This also means we also have to take into account of what happens at −q 0 .We have so that imposes the complex conjugated version of Eq. ( 57).We then have a complex Gauss's law whose charge-free condition around q 0 is The Gauss's law at −q 0 is the complex conjugate of it, so we only need to consider one copy of them.
Let us elaborate on the meaning of the Gauss's law appearing at a general wavevector q 0 in real space.We first define the "phase-shifted derivative" ∂ (q0) α .For derivative in a general direction a, we define For example, for q 0 = (π, π) on a square lattice of lattice constant 1, we have which agrees with how we extract the soft mode from an anti-ferromagnetic background [130].More generally, S(r − a) does not have to be on the lattice site if we take a proper coarse-graining procedure, and ∂ (q0) a S(r) is complex.
The phase-shifted derivative ∂ (q0) α is the correct spatial derivative from the expansion around general wavevector q 0 .When it acts on S(r), it yields the correct Gauss's law in momentum space.For example, This again confirms the relation of ik α ↔ ∂ α (omitting some factors from lattice constants).We see that here, although S(r) is real, its phase-shifted derivative can be complex.So indeed the emergent Gauss's law (Eq.( 61)) is defined over complex fields.However, we did not double the number of DOFs or the constraints.This is because we have so the other copy of Gauss's law at −q 0 , which contains shifted derivatives of the form ∂ (−q0) a S(r), is automatically obeyed when the original Gauss's law is.Therefore, nothing gets doubled.Another equivalent point of view is that the DOFs and constraints around q 0 and −q 0 combine together to form the complex-valued field that obeys the complex Gauss's law.Because the complex Gauss's law has two constraints (one on the real component and one on the imaginary one), the counting of DOFs and constraints remain correctly unchanged.
Finally, once Eq. ( 61) is written down, we can separate its real and imaginary components to form two copies of a real Gauss's law.
A special situation -which actually happens often -is when the FT-contrainer is purely real, i.e. we have the condition T(k) = T * (−k).This happens if q 0 is some high symmetry point so that q 0 and −q 0 are identified.For example, if q 0 = 0, or their difference q 0 − (−q 0 ) is a reciprocal lattice vector (q 0 is often on the BZ boundary in this case).Then we have all c ij real, and Eq. ( 61) (or equivalently, its charge-conjugate) has the real space interpretation as the charge-free condition for a generalized Gauss's law where we have defined a generalized differential operator D (ma) a of order m a ≥ 1 on site a.The effective longwavelength Hamiltonian is then (67) in real space.Note that the number of sublattice sites in a unit cell N is not necessarily the number of components of the electric field.The equation (66) needs to be regrouped in terms of different D (ma) a 's.We will see plenty of examples later.

D. Extracting the Gauss's laws: multiple constrainer models
We now discuss the physics when there are multiple constrainers per unit cell.In this case, the Hamiltonian is in its most general form (repeating Eq. ( 21)) There are M FT-constrainers T 1 , T 2 , . . ., T M .At a general momentum q, these FT-constrainers span the space for eigenvectors for the higher dispersive bands.However, different FT-constrainers are not necessarily orthogonal to each other, and each FT-constrainer is not necessarily the eigenvector of a certain band.In this case, there are two possible ways to close the gap.The first way is the same as the single constrainer case, i.e., one (or several) of the FT-constrainers vanishes at q 0 .The second is when a subset of the FT-constrainers become linearly dependent, so that the dimension of the linear space they span (i.e. the number of the non-flat higher bands) decreases.
To extract the Gauss's law, the core idea is the same as before: we would like to know the eigenvector configuration on the higher dispersive band in the vicinity of the q-points where it becomes gapless.However, more care is needed since the FT-constrainers themselves are not necessarily the eigenvectors we look for.To find the eigenvector, one has to make sure that the orthogonality condition is satisfied.This is just an exercise in linear algebra.
Let us use the case of two FT constrainers T 1,2 (q) as an example.In the first case, when one of the constrainers Breathing pyrochlore model [68] (when generalized to 3D) Equation ( 125) anisotropic U (1): Anisotropic honeycomb-snowflake Model, Sec.VI C vanishes, let us assume T 1 (0) = 0 without loss of generality.We then have T 1 (k) as a vector polynomial Taylor expansion in powers of k α , and we keep only the leading order term in each of its components.The Gauss's law should be extracted using Here, the second term on the right hand side is to project out the part of T 1 that is along the direction of T 2 , so that the rest, T, is orthogonal to T 2 .Since T is still in the space spanned by the FT-constrainers, it is then guaranteed to be the eigenvector of the band that becomes gapless at 0. We can use T 2 (0) instead of T 2 (k) because only the leading order term needs to be kept.
In the second case mentioned above, the FTconstrainers T 1 and T 2 become linearly dependent at q = 0. Let us separate T 1 via So we know and its Taylor expansion is some polynomial of k α for each of its components.The Gauss's law can then be extracted via The above considerations can be generalized to the case of more constrainers.In each case, suppose we need to do Taylor expansion on T 1 or δT 1 , then we should first find an orthognal basis of the linear space spanned by T 2 , . . ., T M .Let us denote the unit vectors of this basis by T ′ 2 , . . ., T ′ M , then Eq. ( 68) should be replaced by and Eq. ( 71) should be replaced by

E. Transitions between different algebraic CSLs
We can classify different algebraic CSLs by examining their gap-closing points.Specifically, two algebraic CSLs belong to the same class if one can smoothly transform the constrainer Hamiltonian and the Gauss's law of one CSL into that of the other, without encountering singular processes that involve merging, splitting, or lifting any of these points.On the other hand, two algebraic CSLs are considered distinct if they have a different number of gapclosing points or if their associated Gauss's laws involve a different number of effective electric field degrees of freedom or a different order of ∂ x and ∂ y .It is impossible to make these gap-closing points identical without going through certain singular transitions.
By identifying the emergent Gauss's law with the structure of the gap-closing point, we can also study the transition between different algebraic CSLs as merging/splitting of the gapless points on the bottom flat band.
The simplest structure of the band-touching point is the one associated with the (complex) Maxwell Gauss's law, shown in the first row of Table .II.Let us call it the basic band-touching point.Other band-touching points corresponding to more exotic Gauss's laws can often be obtained by merging some of the basic band touching point.
Often, the scenario is the following (see e.g.[69]).We start with an algebraic CSL with only basic bandtouching points in momentum space.By tuning some parameters of the Hamiltonian, the positions of the basic band-touching points can be changed, or new basic bandtouching points can emerge when higher bands come down to zero energy.When the parameters are tuned to certain critical values, several basic band-touching points can merge to form a new band-touching point.The new band-touching point is then described by a different generalized Gauss's law.
For readers familiar with topological band theory, this scenario is very similar to the knowledge that the Weyl point is the "basic" gap-closing point containing divergent Berry curvature at the singularity, and merging a few Weyl points together generates other types of gapclosings.In fact, the basic band-touching point in the spectrum of the CSL is exactly equivalent to two merged Weyl cones.
From the perspective of the effective theory, this tells us that by taking a few copies of Maxwell electrostatics and tuning them to a critical point, one can obtain more general forms of U(1) electrostatics.
A summary of the important points in the classification of algebraic CSL is shown in Table .I. In Sec.VI B, we will see concretely how transitions between algebraic CSLs happen in the case of the honeycomb-snowflake model.

VI. ALGEBRAIC CSL MODELS
In this section, we will analyze many old and new examples of algebraic CSLs using our classification scheme as well as tools from flat band theory introduced in Sec.IV B. A survey of various CSL models and prior studies, all fitting within the present classification, can be found Table.III.

A. Checkerboard, kagome, and pyrochlore AFM
To understand how the classification scheme works on concrete examples, let us first apply it to the checkerboard [58], kagome [58,63], and pyrochlore [58,59] antiferromagnets in the large-N limit.These models, due to their geometric frustration, were the first ones discovered to host spin liquids, and are perhaps the most familiar to readers.

Emergent Gauss's law from the checkerboard AFM
Let us demonstrate our classification scheme with the checkerboard lattice model.The model is illustrated in Fig. 7(a).The spins sit on the edges of the square lattice, and the constrainer Hamiltonian is (74) Note that there are N = 2 inequivalent sites in the periodic unit cell.Without loss of generality, we can take the spin on sites 1, 4 in Fig. 7(b) to be the first and second sublattice DOFs in one unit cell, respectively.In this convention, the spins on site 2 and 3 are related by lattice translation to the other two sites: the spin on site 2 is a second sublattice DOF in the unit cell to the left, and the spin on site 3 is a first sublattice DOF in the unit cell below.Therefore, the constrainer is (see Fig. 8 on how each spin maps to each term in the constrainer) The FT-constrainer is then The Hamiltonian in momentum space is Its spectrum is illustrated in Fig. 7(c).We see that it has gapless points at q = (±π, ±π).We can expand the FT-constrainer around q = (π, π) to get (upon adding an overall factor −i) This gives us ground state constraints which is exactly the expected Maxwell U(1) Gauss's law upon identifying the spin sites with the components of the electric field: E x ≡ S 2 , E y ≡ S 1 .The charge-free Gauss's law is shown as pinch poings around these gapless points in the equal-time spin structure factor (Fig. 7(d)).
Finally we note that when writing down C CB (R, r), we made the "gauge choice" equivalent to treating two sublattice sites to be at their physical locations in the unit cell.One can also use other gauge choice (for example, assuming they are at the same position in the unit cell) as long as the complex phase factor is taken care of.
We also note that, on the checkerboard lattice, if the constrainer is symmetric regarding inversion about the center of the constrainer (the vertex of the lattice), the  74)-( 75)).(c) Spectrum ω(q) that arises from diagonalizing the Hamiltonian Eq. ( 77).There is one flat band at the bottom of the spectrum and a dispersive upper band with gap-closing points between them.(d) Spin structure factor showing pinch points at the position of gap-closing points.FIG. 8. How to write down the vector form constrainer C CB (R, r) (Eq.( 75)) from its real space image (Fig. 7(b)).
spectrum is guaranteed to be gapless at (±π, ±π).Such constrainers include the one we used above, and also more generalized ones containing spins on sites farther from the vertex.
The argument, which works for the checkerboard lattice (but not all other lattices), is the following.For the first sublattice sites, if the constrainer involves a spin at site a 1,1 = (r ′ x , r ′ y ) relative to its center set at R = 0 (the vertex) with coefficient c 1,1 , then it also involves a spin at site −a 1,1 , with the same coefficient for the second spin.So the first element of the constrainer must have a pair of terms in the form of Note that, due to the symmetry, any term in the constrainer appears in the form above. Hence, we know the first component of the FT-constrainer must look like [T CB (q)] 1 = 2c 1,1 cos(q • a 1,1 ) + 2c 1,2 cos(q • a 1,2 ) + . . . .(81) Since the vector a 1,1 , pointing from the lattice vertex to the sublattice site on the checkerboard lattice, must be where n x,y are integers, the term cos(q • a 1,1 ) is guaranteed to vanish for q = (±π, ±π) . ( This applies to any other terms in T CB (q), so the FTconstrainer must vanish at q = (±π, ±π), at which point the spectral gap between the dispersive top band and the bottom flat band closes.We hence conclude that given the checkerboard lattice crystalline symmetry, and the properly-chosen action of the constrainer under the crystalline symmetry, existence of gapless points in the spectrum is guaranteed, i.e., the algebraic CSL is protected by symmetry.
Such analysis can be generalized to all crystalline symmetries and their associated constrainer behaviors.Given the proper combination of them, the band touching points are protected and the CSL has to be an algebraic CSL.A systematic examination of all crystalline symmetries and constrainer behaviors is achievable, but lies beyond the scope of this work.

Emergent Gauss's law from the kagome AFM
Next we discuss the kagome lattice model with AFM interactions (Fig. 9 84)-( 86)).(c) Spectrum ω(q) that arises from diagonalizing the Hamiltonian Eq. ( 84).There are one flat band at the bottom of the spectrum and two dispersive upper bands with gap-closing points between them.
The two constrainers, written in vector form, are The FT-constrainers are Since there are two constrainers, there is one flat bottom band and two upper dispersive bands at a general momentum q.However, at q = 0, the two constrainers become linearly dependent, which means a gap closing happens there, as shown in the spectrum in Fig. 5. Hence we have to expand T KGM1 around q = 0, and take its component perpendicular to T 0 , which is Expanding T KGM2 and extracting its perpendicular component gets the same result.Fourier-transforming the constrainer to the real space as in Section V C, we obtain Gauss's law in the form of Maxwell's U(1) theory: Note that the number of sublattice sites does not necessarily need to equal to the number of components of the electric field.Here, the DOF (S 1 + S 2 + S 3 )/ √ 3 is not involved in the low energy physics.It is instead relevant to the third band on top whose eigenvector is T 0 .
The same physics can also be obtained by analyzing the bottom band eigenvector and fluctuator, which has been discussed in Sec.IV B. However, in general it is easier to use the higher dispersive bands because their eigenvectors can be obtained analytically as shown here.

Emergent Gauss's law from the pyrochlore AFM model
The third model we review is the Pyrochlore AFM Model.The lattice is a network of tetrahedra, shown in Fig. 11(a).Its Hamiltonian also contains two constrainers (Fig. 11(b)), written as The treatment is very similar to that of the kagome AFM model.For completeness, let us write down all the steps again.
The two constrainers, written in the vector form, are where a i 's are along the edges of the tetrahedron: The FT-constrainers are Again, since there are two constrainers, there are two flat bottom bands and two higher dispersive bands at a general momentum q.However, at q = 0, the two constrainers become linearly dependent, which means a gap closing happens there (ref.spectrum in Fig. 11(c)).Thus we expand T PC1 around q = 0, and take its component perpendicular to T 0 , which is This yields the Gauss's law of 3D Maxwell U(1): The gapless points also appears at q = (± √ 3π/2, 0, 0) and its cubic rotations.At these points, the equaltime spin correlation, shown in Fig. 11(d), exhibits the two-fold pinch points (2FPP), which is the canonical hallmark of the emergent U(1) electrostatics physics.We note that, although the spectrum is also gapless at q = (0, 0, 0), the equal-time spin correlation does not observe a pinch point there.But this is merely due to the cancellation of intensity when all spin correlation channels are summed together.

Emergent Gauss's laws from the honeycomb-snowflake model
Now let us apply the classification algorithm to the honeycomb-snowflake model [69], which we introduced in Sec.III A.
The model is defined on the honeycomb lattice, with a spin on each site, which we treat as a scalar within the large-N approximation.The Hamiltonian is The sum of R is taken over all unit cells, which is best visualized as hexagonal plaquettes.The constrainer C γ (R) defined on the hexagons contains weighted sums of spins around each hexagon shown in Fig. 12: The constrainer reads  105)).This figure is a replication of Fig. 2, reproduced here for convenience.
Here, r j is the vector from the center of the snowflake to the corresponding site j labeled in Fig. 12(b).Figure 13 shows how the first element of the constrainer is constructed by going over all the first sub-lattice sites in the adjacent unit cells.
The FT-constrainer is then obtained by Fourier transforming C γ HS (R, r), And the Hamiltonian in momentum space is The spectrum structure is plotted in Fig. 3 for different values of γ.It has two bands.The top band always undergoes a gap closing at wavevector Let us now examine the physics for small k = q − q 0 for two cases: γ = 0 and γ = 1/2.
When γ = 0, we have, at leading order, the FTconstrainer So the spin fluctuations around the ground state satisfy the constraint (111) Reorganizing the DOFs, we have which upon Fourier tranformation into the real space yields a Maxwell Gauss's law which acts on a complex electric field E = (i( S1 + S2 ), S1 − S2 ).Note here that because of the phase shift arising from expanding around finite momentum q 0 , Sa are themselves complex in real space and therefore E x is not purely imaginary and E y is not pureley real.
If we separate Eq. ( 113) into real and imaginary parts we may consider it as two real Gauss's laws.States satisfying these two Gauss's laws are also guaranteed to satisy the Gauss's laws that would be obtained from an expansion around the singular band touching at q = −q 0 , due to the property Sa (q) * = Sa (−q).We thus have two real Gauss's laws in total, which is to be expected, as there are two singular band touchings per BZ.
When γ = 1/2, on the other hand, we obtain the FT-constrainer T1/2 We then have the emergent Gauss's law in the form which we rewrite as ) Again this is a complex Gauss's law.If we identify a traceless, symmetric complex matrix to be then the Gauss's law becomes which is a (complex) realization of the electrostatics for a symmetric rank-2 U(1) gauge theory.
Breaking the complex Gauss's law into real and imaginary parts, we obtain two real Gauss's laws.As before, these also take care of the constraints arising from the band touching at −q 0 , and the presence of two Gauss's laws agrees with the presence of two band touchings in the BZ.

Transition between algebraic CSLs
Let us now study the transition between different algebraic CSLs in the honeycomb-snowflake model.We will study the transitions near two critical points: γ = 1/3 and γ = 1/2.
For the critical point γ = 1/3, a new band touching point emerges at the mid point of the BZ boundary, as illustrated in Fig. 14(b).This happens, as the top band gradually moves down, when the top band touches the bottom band, as γ → 1/3 − .The new band touching then splits into two band touching points as γ increases above 1/3, see Fig. 14(c).Each single band touching point is associated with a Maxwell's U(1) Gauss's law.
For 1/3 < γ < 1/2, there are band touching points on the BZ corner and boundary.Each single band touching point is associated with a Maxwell U(1) Gauss's law, as just shown.This can also be seen from Fig. 3, where the structure factor (defined in Eq. ( 14)) on each band touching point exhibits the characteristic two-fold pinch point.
As γ increases and approaches 1/2, three band touching points on the BZ boundary move toward the fourth one on the BZ corner, as illustrated in Fig. 15.At the critical point γ = 1/2, the four points merge together (see Fig. 15c), and form a new band touching point with a different structure: one associated with the rank-2 U(1) Gauss's law shown in Eq. (118).
The lesson we learn here is that the transition between different algebraic CSLs can be understood as the emergence/disappearance and merging/splitting of the band touching points in their spectrum.Mathematically, such transitions are described in the same way as in topological band theory, and much prior knowledge can be borrowed to understand transitions of algebraic CSLs.This will be a topic for future study.

Symmetry and topological protection of the gapless points
In the honeycomb-snowflake model these gap closings are symmetry-and topologically protected, provided that the Hamiltonian respects inversion symmetry.This is because inversion symmetry requires the two components of the constrainer to obey Hence, the components in the FT-constrainer T(q) are related by Combining this with the normalization of the eigenvector implies T(q) can always be written in the form and can thus be represented by a point on the unit circle ϕ(q).We can thereby define a winding number of the vector field around closed paths in reciprocal space: The topologically stable gap closing points correspond to vortices of v(q) with a finite integer winding number for closed paths encircling them.The gapless points are thus topologically stable and cannot be removed by small changes to the ground state constraint, provided that inversion symmetry is maintained.
Let us now revisit the ground state phase diagram of the honeycomb-snowflake model shown in Fig. 3. Transitions between distinct CSLs occur as γ is varied, via pair creation/annihilation of vortices, or by coalescence of vortices with like winding number.The various CSLs have a distinct arrangement of singularities (known as pinch points) in their spin correlation functions, affirming their distinctive nature.At γ = 0 the model is an algebraic CSL with gap closings and corresponding pinch points at the Brillouin zone corners (K points) (Fig. 3  (a)).On increasing γ this remains the case until γ = 1/3, at which point pairs of oppositely charged vortices nucleate at the M points of the BZ (Fig. 3 (b) and Fig. 14(bc)), in addition to the existing pinch points at the K wavevector.This leads to a new CSL with 8 pinch points per BZ instead of only 2, with all pinch points on the zone boundaries (the points shared by several adjacent Brillouin zones are included only once in this count).As γ is further increased the vortices formed at the M points migrate towards the K points, such that three vortices of one charge converge on one of the opposite charge (Fig. 3 (c) and Fig. 15(c)).This leads to the formation of vortices with winding number +/ − 2 at the K points when γ = 1/2, and four-fold pinch points in the spin structure factor (Fig. 3(d) and Fig. 15(c)).This is indicative of a spin liquid described by a higher-rank U(1) gauge theory [68,100] in Eq. ( 118), as explained earlier.On increasing γ further the vortices at the zone corners separate again and the system enters a new CSL, with 8 pinch points per BZ but now with 6 of them in the interior of the BZ rather than on the boundary (Fig. 3 (l,p)).
For negative γ, the story is similar, and the readers can refer to the original paper [69] for more detail.

C. Anisotropic U(1) CSL
The honeycomb-snowflake model with γ = 1/2 provides a simple example of a classical spin liquid with an isotropic Gauss's law, as shown in Ref. 69.Here, we propose a simple model exhibiting a spin liquid described by an anisotropic Gauss's law, and demonstrate its nature using the algebraic classification from Section V.The model can be considered as a generalisation of the honeycomb-snowflake model with explicit lattice symmetry breaking.Specifically, we take the honeycombsnowflake model with γ = 0 such that the constraint on each hexagon only involves spins belonging to that hexagon.We then adjust the contribution of each spin to the constraint according to a new parameter β, such that spins at the top and bottom of the hexagon contribute to the ground state constraint with weight β and the others with weight 1: with the sites numbered around each hexagon as shown in Fig. 16(b).The case β = 1 corresponds to the isotropic honeycomb model from Ref. 64.
Upon increasing β from β = 1, gap closing points migrate along the Brillouin zone boundaries normal to the q y axis.At β = 2 they merge at the M point of the Brillouin zone (q = q M ).At this merging point, we can expand the FT-constrainer T(q) around q = q M .Here, similarly to Section VI A, we use a gauge in which we reference the spins to the position of the centre of their unit cell, rather than their physical position on the lattice.This leads us to: The dispersion ω(q) is anisotropic around the band touching, having a form: To obtain the Gauss law we use T * a (q) • Sa (q) = 0.By adding a phase S′ 2 = S2 e −i2π/3 , we have and therefore a real space Gauss' law: where we have identified the electric field components with the suitable combination of the fluctuating spin variables Si .
We found the conserved quantities to be the following.First, the obvious one is the net charge conservation: We can also look for other conservation laws defined by a suitably chosen function f (x, y) in the integrand: after integration by parts.To make sure this is zero, we need to choose f such that the following two conditions are simultaneously satisfied: The solutions are for any choice of real numbers a i .Hence we deduce that the second conserved quantity is the charge dipole in the x direction: The charge therefore has reduced mobility as in fracton theories [127]: it is immobile in x direction but can move in the transverse (here, y) direction.This model can actually be viewed as a 2D cut of the generalized U(1) gauge theory for Haah's code (before Higgsing).Its structure factor is essentially identical to that proposed in Ref. 88, featuring pinch points with parabolic contours.The evolution of the spin correlations on tuning through the critical point is shown in Fig. 17.
It is worth noting that despite the anisotropic Gauss's law, this spin liquid is not a "Type II" fracton phase [88], which requires an infinite number of conservation laws.Here we have only a finite number.The parabolic pinch points should therefore be understood as a signature of an anisotropic Gauss's law, and not necessarily as a signature of "Type II" fracton phases, as proposed in Ref. 88.This spin liquid occurs at the special point of parameter space β = 2.For β > 2, there are no band touchings in the Brillouin Zone and the momentum space correlations are smooth.The model for β > 2 connects smoothly to the β → ∞ limit, which is a trivial paramagnet in which a pair of spins is coupled within each unit cell, but there is no inter-unit-cell coupling.Hence, the anisotropic spin liquid occurs at the transition point between a Coulomb phase and a short-range correlated trivial paramagnet.D. Higher-dimensional gapless manifolds: pinch lines etc.
The discovery of various forms of algebraic spin liquid based on gapless points of the Hamiltonian J(q) in Eq. ( 29) leads to a natural question: are there spin liquids associated to nodal lines of J(q) and, if so, what are their properties?The connection between gapless points in J(q) and pinch points in the spin structure factor S(q) (see Eq. 14) is suggestive of a generalization to nodal lines, i.e. the extension of pinch point singularities along lines of reciprocal space, namely pinch lines.
Such features have previously been found in a classical spin liquid based on an anisotropic spin Hamiltonian [66], and indeed a soft spin treatment of this model finds nodal lines in the dispersion attached to the flat bands at the bottom of the spectrum.The spin liquid in [66] thus establishes one example of a nodal line spin liquid.
Here we present a new, simple model, of a nodal line spin liquid with isotropic spin interactions, based on the concept of the symmetry protected topological phases, which we have previously applied to other algebraic CSL.
To motivate the construction, we consider once more the honeycomb-snowflake model (see Section III A and earlier subsection VI B for the definition of this model).The Fourier transformed constrainer T(q) has two components, listed in Eqs. ( 107) and (120), corresponding to the two sites per unit cell and obeys the relation T 1 (q) = T 2 (q) * due to inversion symmetry.When normalised, as done in Eq. (120), this means that T(q) lies on the unit circle and its evolution in reciprocal space can support stable vortices corresponding to the nontrivial homotopy classes π 1 (S 1 ) of the phase ϕ(q).The singularities in the center of such momentum-space vortices correspond to band touchings in J(q) and pinch points in the equal-time spin structure factor S(q) (see Eq. ( 14)).These singularities are then protected in the sense that they cannot be removed by small changes to the ground state constraint which respect inversion symmetry.The presence of these vortices arose directly from a two site unit cell and a symmetry constraining T(q) onto the unit circle.
In two dimensions vortices are point like, but in three dimensions they are line like.The above considerations lead us to expect that a classical spin liquid with two sites and one constraint per unit cell, and inversion symmetry should support pinch lines.
One such example is found on the lattice shown in Fig. 18(a).This lattice is formed from octahedral units which share edges in the xy plane and join at vertices in the z direction.There are two sites per unit cell, indicated in red and blue in Fig. 18(a).We write down a Hamiltonian on this lattice, as a sum over octahedra: The resulting soft-spin dispersion ω(q) has two bands, a lower flat band and a dispersive upper band.The upper band meets the flat band along the edges of the q z = ±π faces of the Brillouin zone, i.e. along q = (q x , π, π), q = (π, q y , π) and equivalent directions.The location of the nodal lines is illustrated in Fig. 18(b).
The structure factor S(q) for the moodel is depicted in Fig. 19, as a series of cuts at fixed values of q y .These cuts intersect the pinch lines along the lines q = (±π, q y , ±π), and thus four pinch points are visible in each panel where the plane cuts the pinch line.
We thus establish a simple model for a spin liquid with pinch line singularities.Based on the topological considerations outlined above, pinch lines should be common to inversion symmetric three-dimensional classical spin liquids with two sites and one constraint per unit cell.There is a two site unit cell, with inequivalent sites here indicated in red and blue.Defining a local constraint on the octahedra leads to a classical spin liquid with nodal lines in J(q) and hence pinch lines in S(q).(b) Location of nodal lines for the model defined in Eq. ( 135).The nodal lines appear at wavevectors q = (qx, π, π), q = (π, qy, π) and equivalent, creating a network of nodal lines along the edges of the qz = ±π faces of the Brillouin zone.

VII. FRAGILE TOPOLOGICAL CSL CLASSIFICATION: EIGENVECTOR HOMOTOPY A. The topological classification
Next we discuss the other category of classical spin liquids: the fragile topological CSLs with short-range spin correlations (the meaning of the qualifier 'fragile' will be explained later in section VII B 2).A fundamental difference between this category and the algebraic one is that in fragile topological CSLs there is no band-touching between the higher bands the bottom flat one(s).In-stead, the bottom flat bands are gapped from all other bands in the spectrum.In real space, this means that all L x L y local fluctuators (due to translation symmetry on an L x × L y lattice) are linearly independent and form a complete basis, thus accounting for all the ground states in the flat band.The absence of band-touching points means there are no emergent Gauss's laws describing the CSLs.For the same reason, the spin correlations decay exponentially instead of algebraically with distance.
In this category, we can still ask the question about the classification of fragile topological CSL models.More precisely, we consider two CSL systems A, B that have the same number of DOFs per unit cell, and same number of gapped flat bands, and ask if it is possible to adiabatically tune CSL A into B, while keeping the system in a CSL state (i.e., maintaining the flatness of bottom bands)?In terms of constrainers, this is to ask if we can smoothly change the L x L y constrainers C A (R)'s into C B (R)'s without making them linearly dependent at some point (for simplicity we use the one-constrainer Hamiltonian but its generalization is straightforward).Although all of C A (R)'s are linearly independent (i.e. the corresponding T(q) never vanishes and there is no band touching point), and so are the C B (R)'s, in the process of tuning we may have to go through a boundary point in parameter space whose constrainer C X (R)'s are not linearly independent anymore.In the spectrum of the Hamiltonian, this would manifest itself in a gap closing.If such an intermediate gapless point is unavoidable, then we say that the two CSLs A and B belong to distinct equivalence classes.If, on the contrary, an adiabatic tuning of the constrainers from C A (R) into C B (R) is possible without closing the spectral gap, we identify the two CSLs as belonging to the same equivalence class.
FIG. 19.Spectrum and S(q) for the octahedral nodal line model (Eq.( 135)), taking cross sections at a series of fixed values of qy and qz, cutting through the nodal lines at different points (see Fig. 18).A pinch point is present at q = (±π, qy, ±π) for all values of qy, thus forming an extended, line-like singularity: a pinch line.
The reason we make this distinction is because given the short-ranged spin correlations, one may naively expect all CSLs in this category to be equivalent to a trivial paramagnet.The trivial paramagnet is defined as systems where spins only interact within a unit cell, and there is no inter-unit cell couplings.Using the breathing kagome lattice as an example, a trivial paramagnet is given by a Hamiltonian which constrains the spins only on the up-pointing triangles As we see here, it is a model with two DOFs in the unit cell freely fluctuating, while the other DOF completely frozen to be zero.More importantly, there is no interunit cell coupling, so the higher dispersive band has a constant eigenvector T (q) = (1, 1, 1)/ √ 3. Given another breathing kagome model with one constrainer per unit cell one can ask the following question: if we keep the one-constrainer form of the Hamiltonian but change the constrainer smoothly to tune the model from the trivial model, Eq. ( 136) to the new model, can this procedure happen without closing the gap in the spectrum at any step?We will see later an example of a FT-CSL CSL which can be shown by such an argument to not be equivalent to a trivial paramagnet (Section VIII).
As one may expect, if adiabatic transitions between two CSLs are obstructed, there must be some mathematical quantities that distinguishes them.The idea is very similar to the notion of Chern insulators in band theory, wherein two theories with different Chern numbers cannot be adiabatically transformed into each other by tuning the Hamiltonian, without the gap closing.The classification can be further enriched by symmetry: while there are paths to deform C A to C B without closing gaps, it is only possible to do so when the path breaks a symmetry.In such a symmetry-enforced scenario, the two states A and B are still considered to be different.
We will show that up to the equivalence of adiabatic connection, the CSLs can then be divided into different topological classes.What is the topological quantity that distinguishes the different topological classes?Since the bottom band eigenvector is globally well-defined (in mathematical terms, it is a section of a trivial vector bundle), the band always has zero Chern number, so that is not the quantity we look for.Instead, we found that the fragile topological CSLs are classified by the homotopy class of the bottom band eigenvector subspace configuration on the torus of the BZ.When there is only one bottom band, the eigenvector subspace is simply the Ncomponent unit vector modulo an overall phase.When there is more then one bottom band, the eigenvector subspace is a higher dimensional subspace of the total space of all the eigenvectors.We now consider these two cases in more detail.

One bottom band
Let us consider the simplest case first: a 2D model with N spins per unit cell, and one bottom flat band in the spectrum of its Hamiltonian (generalizations to 3D exist, but we shall primarily focus on 2D models in what follows).The flat band has normalized eigenvector configuration, which we refer to as the fluctuator, see Eq. ( 23): The fact that the bottom band is gapped from the other bands means B(q) is well-defined and non-vanishing everywhere in the momentum space.
At a fixed wavevector q, B(q) and e iθ B(q) correspond to the same physical spin fluctuation.Therefore, the physical configuration space for B(q) is the complex projective space CP N −1 .Often, we have additional inversion or time reversal symmetries that constrain B(q) to be real, in which case the physical configuration space is the real projective space RP N −1 .From now on, we take B(q) to denote a ray in the target space CP N −1 or RP N −1 .Now B(q) defines a map from the torus of the (twodimensional) BZ to the space of CP N −1 or RP N −1 : B(q) : T 2 → CP N −1 (or RP N −1 ); q → B(q) .(138) The equivalence classes of such maps are classified by the relative homotopy group The homotopy classes are the topological quantities that distinguish different fragile topological CSLs.Without closing the gap, i.e., having B(q) vanishing and B(q) ill-defined at some momentum point, the homotopy class cannot be changed.Hence two fragile topological CSLs of different homotopy classes cannot be adiabatically turned into each other without closing the gap.Obviously, the comparison of homotopy classes is only sensible when B(q)'s have the same number of components.This indicates that the topology classification is a fragile concept (hence our use of the term fragile topological CSL), which we will explain in more detail below in subsection VII B 2.
In general, the homotopy group [T 2 , X] is not easy to compute.However, if the target manifold X is simply connected (i.e., π 1 (X) = 0 and path-connected), then we have the homotopy group isomorphic to the second homotopy group of X: Since π 1 (CP N −1 ) = 0 for any N − 1 ≥ 1, we have This is the homotopy class for complex eigenvector B(q).That is, in general, the homotopy classes are labeled by an integer number in Z.
However, in some scenarios, we can consistently assign directions to the RP N −1 eigenvectors smoothly over momentum space, without encountering any inconsistencies with the BZ periodic boundary condition.Then, we can treat the eigenvectors as unit vectors pointing on the S N −1 sphere.For N − 1 ≥ 2, π 1 (S N −1 ) = 0 so S N −1 is simply connected.In this case we have We see that here the only non-trivial case is when the model has N = 3 degrees of freedom.The integer homotopy invariant is then nothing but the skyrmion number on the 2-torus.Given the eigenvector configuration, this skyrmion number n sk can be computed by We mention in passing that here the skyrmion number should not be confused with the winding of effective magnetic field for a two level Hamiltonian J(q) = B(q)• σ.In that case, the skyrmion configuration makes it impossible to smoothly define the phase of a band, which is equivalent to the statement of a nontrivial Chern number of the bottom band.By contrast, the skyrmion characterizing the band eigenvector in a CSL has a fundamentally different physical meaning.In fact, when a band's eigenvector is well defined in the BZ (so we can talk about its skyrmion number to begin with), there is no problem in smoothly defining the phase of the band at all wavevectors -it is then a section of a trivial vector bundle, with zero Chern number.In particular, exact flat bands with finite range of interactions have been shown to be alway have zero Chern number [118].

N − 1 bottom bands
Another equally simple case is when we have a single constrainer in the Hamiltonian Eq. ( 1), resulting in N −1 bottom bands and one disperstive top band separated by a gap.In this case we can examine the homotopy for eigenvector T(q) of the top band instead.All the analysis in the previous sub-sub-section carries over by replacing B(q) with T(q).Equations ( 140)-( 141) applied to a single bottom band or a single top band explicitly tell us the possible homotopy classes of the corresponding cases.In Sec.VIII, we will present a concrete microscopic spin model which exhibits several topological classes [T 2 , S 2 ] = Z and transitions between them, as the Hamiltonian is tuned.

Other cases
The more complicated situation is to have N − M degenerate bottom bands where 1 < M < N − 1.In this case, the target space is not a ray in CP N −1 or RP N −1 , but a projective plane (for two bottom bands), or generally the projective (N − M )-dimensional subspace in CP N −1 (or RP N −1 ) generated by the N − M eigenvectors B1 (q), . . ., BN−M (q).These homotopy classes are in principle calculable, though we are not aware of a simple closed-form expression.

B. Properties of the fragile topological CSLs
We now discuss the general properties of the fragile topological CSLs.For concreteness, we use the one bottom flat band (or equivalently one top band) cases for demonstration.Our discussion is straightforward to generalize to multi-band cases.

Transition between homotopy classes
The homotopy equivalence class remains unchanged upon adiabatically tuning the Hamiltonian while keeping the bottom bands flat and gapped.The changing of topological class thus requires the gap between the higher bands and bottom bands to close, so that the bottom band eigenvector configuration can go through singular changes at the gap closing point.
From this point of view, while all CSLs' Hamiltonians are fine-tuned, the fragile topological CSLs are the more common ones over the entire parameter space.The algebraic CSLs, with the spectral gap closing, require additional tuning and denote the critical boundaries between different fragile topological CSLs, or higher order critical points where the critical boundaries intersect.In this sense, the algebraic spin liquids are more fine-tuned than the topological CSLs.A schematic phase diagram indicating both types of spin liquids is shown in Fig. 1.

Eigenvector homotopy is fragile
The homotopy class is a "fragile" topological quantity (see e.g.discussion in Refs.[5,135,136]), in that upon adding a new spin DOF per unit cell to the model, the previous non-trivial homotopy class of the N -component eigenvector configuration may become trivial as a (N +1)component eigenvector configuration.
Let us demonstrate this using the following model.Consider the original model with N sublattice sites and only one constrainer C: Such a model has N −1 degenerate flat bands and 1 higher dispersive band.The top band has eigenvector T 0 (q) obtained by Fourier transforming the vector C 0 (r, R).
We now add a new DOF S N +1 to the system, and introduce a parameter γ to tune the interactions.The new Hamiltonian is which has one higher dispersive band and N bottom flat bands.The corresponding constrainer in its vector form is then and the FT-transformed constrainer is Now, we adiabatically tune the Hamiltonian parametrized by γ going from 0 to 1.
Note that the norm square of T(γ, q) is always positive everywhere, so the gap between the bottom band and the higher bands never closes.However, at the end of this adiabatic tuning, the eigenvector becomes T(1, q) = (0, . . ., 0, 1) , which belongs to the trivial homotopy class for the N + 1 band model.By the above argument, the N band model of any homotopy class can be adiabatically changed to the trivial homotopy class of the N +1 band model.We can also join two such processes together to adiabatically change between the two different homotopy classes of the N band model without closing any gap: , q) = (0, . . ., 0, 1) We would like to stress that such a construction is not possible without introducing the new DOF.That is to say, the homotopy class of the eigenvector configuration is a fragile topological quantity only when allowing for arbitrary 'padding' the unit cell with new DOFs.However, when restricted to the original degrees of freedom, the homotopy class is well-defined and can only be changed via spectral gap closing.

Absence of algebraic boundary correlations
One may naturally wonder whether fragile topological CSLs have a notion of bulk-boundary correspondence, in analogy with topological insulators.Specifically, given the association between gapless points in the band structure of J(q) and algebraic correlations of a CSLs, one may have hoped that topological CSLs would host gapless edge states of J(q), and therefore algebraic correlations at their boundaries.
However, as we have argued above that topological CSLs are generically fragile in nature; and since fragile topology does not guarantee gapless boundary modes [41], the scenario of algebraic boundary correlations is not realized.Fragile topological CSLs with open boundary conditions will generally have short ranged correlations on the edge, as well as in the bulk.
It is known that fragile topology can have an associated bulk-boundary correspondence in the presence of specially chosen twisted boundary conditions [137].However, the naturalness of such twisted boundary conditions in a CSL is doubtful, so we do not pursue this topic further in this work, and leave its investigation for the future studies.How the fragile topological CSLs manifest its non-trivial topology in experimentally measurable quantities is still an open and important question.In this section we introduce a generalization of the kagome-hexagon model [65] introduced earlier in Eq. ( 15), which is used to demonstrate the application of our scheme of fragile topological CSL and establishes the possibility of transitions between distinct fragile topological CSLs.We refer to this generalized model as the kagome-star model.The Hamiltonian is: where ζ is a dimensionless (real) tuning parameter and The two contributions to the constrainer C KSα in Eq. ( 151) are illustrated in Fig. 20.The first sum over i is a sum over spins 1, . . ., 6 belonging to the interior of the hexagon centred at R. The second sum over j is a sum over spins 1 ′ , . . ., 6 ′ connected to the exterior of the hexagon, forming the points of a six-pointed star.The ground states are those which satisfy the constraints Since the three components α = x, y, z are identical and decouple from each other, we can focus on one copy of them now and drop the index α.There is one star motif -and hence one constrainer -per unit cell.The Fourier transformed constrainer is: For ζ = 0, this model reduces to the Kagome-Hexagon model [65], in which T 0 (q) is well-defined and non-zero for all q.Correspondingly, the soft spin band structure is gapped everywhere, with two flat bands at the bottom of the spectrum separated from one top band.
The top band is topologically non-trivial, as can be seen through calculating the momentum space skyrmion number associated to T ζ (q): Q sk (ζ) = 1 4π EBZ d 2 q Tζ (q) • Tζ (q) ∂q x × Tζ (q) ∂q y , (154) where T(q) = T(q)/|T(q)| is the normalized constrainer.The integral is taken over the extended Brillouin zone (EBZ), corresponding to the periodicity of T(q), which due to the relative phase of different sites within the unit cell has double the period of the primitive Brillouin zone.
Provided that inversion symmetry is maintained and no further sites are added to the unit cell, Q sk ∈ Z takes integer values which can only be changed by tuning the model through a gapless point, as we discussed in detail in Sec.VII.The evolution of Q sk with increasing ζ is shown in Fig. 21.
The Skyrmion number jumps discontinously at ζ = 1/2 and ζ = 1.These changes in Q sk indicate zero temperature transitions between distinct CSLs, all with short range correlations, but distinguished by the homotopy of the momentum space constrainer.At the boundaries between the fragile topological CSLs in parameter space, the soft spin dispersion has gapless points.Based on our discussion above, this implies the emergence of algebraic CSLs at the boundaries.Indeed, calculating the equaltime spin structure factor S(q), defined in Eq. ( 14  As discussed in VII B 3, the fragile topological nature of the CSLs does not guarantee bulk-boundary correspondence, meaning that there no additional gapless points arising at the edge with open boundary conditions, and that the correlations remain short ranged up to the edge of the lattice.This is demonstrated in Fig. 23 where we plot the soft spin dispersions with open boundary conditions in one direction and with fully periodic boundary conditions, for ζ = 0.There is no additional gap closing at the boundary and hence no algeberaic boundary correlations, underscoring the fragility of the topology underlying the short-range correlated CSLs.
We thus establish the kagome-star model as exemplifying a series of distinct fragile topological CSLs, with algebraic CSLs at the boundary between them.
Another example based on the modified kagome-star model can be found in the companion article Ref. 86, which provides a richer phase diagram of different FT-CSL phases and their algebraic CSL boundaries.

IX. ASPECTS OF APPLICATION
A. Self-consistent Gaussian approximation and lattices with symmetry-inequivalent sites Our classification formalism employs a scheme for soft spins in constrainer Hamiltonians, i.e., each spin component is a real scalar free of any non-linear constraints, and the Hamiltonian is written as a sum of squared linear constrainers.Adding to the discussion when introducing this formalism, we comment on aspects of what happens in the siutation when the spins have hard constraints (particularly Heisenberg spin with the hard constraint S 2 = 1), and the bare Hamiltonian does not take the constrainer form.First, note that the constrainer Hamiltonian of soft spins is largely equivalent to the Self-Consistent Gaussian Approximation (SCGA) [63], the applicability of which for hard spins can be justified by Luttinger-Tisza ideas [140,141].In the SCGA, the hardspin constraint is enforced only on average, by using a Lagrange multiplier that can also viewed as a chemical potential term.Thus, up to this chemical potential shift, the classification scheme we presented extends to the SCGA scheme.
In fact, the constrainer Hamiltonian form of Heisenberg spins has been used for analyzing frustrated magnets already.Perhaps the most instructive example is that of the pyrochlore spin models [60,142], where the Hamiltonian is naturally written as a sum over constrainers on all tetrahedra: all tet.i∈tet.
≡ H constrainer + const. ( The Hamiltonian H constrainer is actually what SCGA yields at the limit T → 0. A subtlety to be aware of is when the sublattice sites are not equivalent to each other via space group symmetries.In such cases, as in e.g. the centered pyrochlore lattice, the bottom flat band from diagonalizing the bare Hamiltonian H bare may not satisfy Eq. ( 155) for all sublattice sites.In such cases, a generalized Luttinger-Tisza method proposed by Lyon and Kaplan [143] and improved by Schmidt and Richter [144] can be used to properly derive the physical ground states by "renormalizing" the spins.Nevertheless, in the end, the renormalized Hamiltonian still hosts bottom flat bands, and can be analyzed with our classification scheme.A detailed application of this formalism to the centered pyrochlore lattice model [138] can be found in Ref. 139.This may also be able to explain the observation of disappearance of pinch points for the square-kagome model with varying exchange parameters in Ref. 71.

B. Survey of known CSLs in literature
The classification scheme we propose provides a comprehensive view of both FT-CSLs and algebraic CSLs, and encompasses the majority of classical spin liquid models known in literautre.In this section, we provide Table .III of CSL models found in the literature, as well as those constructed in this work, along with brief comments on their place within our classification scheme.In this table, we see a variety of models that realize different types of algebraic CSLs and fragile topological CSLs, as well as demonstrate the transitions between different CSLs.They all fit snugly into the landscape of CSLs we propose in Fig. 1.

X. SUMMARY AND OUTLOOK
In this work, we have presented a classification scheme for classical spin liquids.The scheme includes two main categories, namely algebraic CSLs and fragile topological CSLs, with a finer classification within each category based on the emergent Gauss's law and homotopy of eigenvectors.Along with placing known examples from the literature into the landscape of CSLs, we introduce new models to illustrate the major aspects of the classification scheme.We also make connections to flat band theory to analyze the structure of the ground state degeneracy in real space.
The classification scheme is a useful tool for under-standing known and new CSL models and constructing new ones using the constrainer Hamiltonian formalism.
We do note that the large-N /SCGA treatment may fail for classical magnets, as it treats the non-linear hard spin constraints only 'on average'.There are known examples where hard spin constraints lead to different -and often interesting -new physics.One such case is the kagome Heisenberg antiferromagnet, which is found to be magnetically ordered, contrary to the large-N prediction [93][94][95][96][97], although our analysis applies for N ≥ 4-component spins on that lattice [145].Another type of cases is when the spins are discrete (Ising spins for example), so that e.g.satisfying the constrainers on triangles with their odd number of sites becomes impossible.The discreteness of the spins is a fundamental obstacle here, and the inter- esting physics of triangular and kagome Ising models lies beyond our present analysis.Therefore, case-by-case analytical and numerical studies are still necessary for specific models of interest.While we believe our scheme to be comprehensive for those spin liquids where a soft spin approximation is appropriate, there remains the interesting possibility of spin liquids outside our scheme in which the non-linear constraints are a crucial element of the effective description.The identification and classification of such cases remains an open question.
It is also interesting to speculate on the fate of CSLs in the presence of quantum dynamics.Although quantum models are usually not solvable, the algebraic CSLs do provide numerous realizations of the electric sector of the generalized rank-2 U(1) electrodynamics, serving as an interesting starting point for constructing quantum spin liquid models that host exotic emergent particles with reduced mobility.
While the constrainer Hamiltonian approach is mathematically convenient for analyzing classical spin liquids, it relies on fine-tuned interactions between nearest and farther neighbor spins for realising interesting new spin liquids in the T → 0 limit.But even when this fine-tuning is not precisely met, there is a good chance that signatures of the spin liquid under consideration will be present at moderately low temperatures.This holds the promise that future developments, both in the realm of magnetic materials but also in cold atomic systems, will provide a realization of some of these models.

FIG. 4 .
FIG. 4. (a) The kagome lattice.It has three sites in the unit cell, forming three sublattices indicated here in red, blue and green.(b) Constrainer of the Kagome-Hexagon model.Classical spins are arranged on a kagome lattice, with ground states defined by the constraint that the sum of spins on each heaxagonal plaquette must vanish (Eqs.(15)-(17)).(c) Spectrum ω(q) that arises from diagonalizing the Hamiltonian (Eq.(15)) in momentum space.There are two degenerate flat bands at the bottom of the spectrum and a dispersive upper band with no band touchings between the upper and lower bands.(d) Spin structure factor showing an absence of singularities.
Flat band theory Classical spin liquid CLS : local eigenstate of the flat band local spin fluctuation within ground states NLSs : non-local eigenstate of flat band non-local spin fluctuation within ground states a singular band touching point effective Hamiltonian indicates the Gauss's law multiple singular band touching points coexistence of different Gauss's laws merging/splitting of singular band touching points transition between different algebraic CSLs no band touching on the flat bands fragile topological CSLs c1 3×3 to shift all bands by a constant)

FIG. 5 .
FIG. 5. (a) Compact local states (CLS) and non-local loop states (NLS) of the kagome model (Eq.(31)), which can also be interpreted as the local and loop fluctuators in the classical spin liquid model (Eq.(32)).One can check that the hopping amplitude from these states to any other site is zero.The CSLs are not linearly independant: adding all of them on the entire lattice yields zero.(b) Spectrum of the Hamiltonian in Eq.(36).

FIG. 6 .
FIG. 6. U (1) Structure of the ground states.States connected by local fluctuators F local are in the same equivalence class, and different equivalence classes are connected via the non-local F loop fluctuators.

FIG. 7 .
FIG. 7. (a) Checkerboard lattice.(b) Constrainer of the checkerboard model.Classical spins are arranged on the edges of a square lattice, with ground states defined by the constraint that the sum of spins on each vertex must vanish (Eqs.(74)-(75)).(c) Spectrum ω(q) that arises from diagonalizing the Hamiltonian Eq. (77).There is one flat band at the bottom of the spectrum and a dispersive upper band with gap-closing points between them.(d) Spin structure factor showing pinch points at the position of gap-closing points.
(a)), which we have already introduced in Sec.IV B in the context of the flat band theory.The Hamiltonian contains two constrainers, as shown in

FIG. 9 .
FIG. 9. (a) Kagome lattice.(b) Two constrainers of the kagome model shown in shaded regions.Classical spins are arranged on the edges of a square lattice, with ground states defined by the constraint that the sum of spins on each vertex must vanish (Eqs.(84)-(86)).(c) Spectrum ω(q) that arises from diagonalizing the Hamiltonian Eq. (84).There are one flat band at the bottom of the spectrum and two dispersive upper bands with gap-closing points between them.(d) Spin structure factor showing pinch points at the position of gap-closing points.
FIG. 9. (a) Kagome lattice.(b) Two constrainers of the kagome model shown in shaded regions.Classical spins are arranged on the edges of a square lattice, with ground states defined by the constraint that the sum of spins on each vertex must vanish (Eqs.(84)-(86)).(c) Spectrum ω(q) that arises from diagonalizing the Hamiltonian Eq. (84).There are one flat band at the bottom of the spectrum and two dispersive upper bands with gap-closing points between them.(d) Spin structure factor showing pinch points at the position of gap-closing points.

FIG. 14 .
FIG. 14.The transition between different algebraic CSLs as the emergence and splitting of the gap-closing points.This figures shows one such transition in the honeycomb-snowflake model around γ = 1/3.The three plots are the zoomed-in view of the spectrum at the center of the BZ edge.The insets on the top left corner show the position of gap-closing points in the BZ (actual distance is exaggerated for better visibility).(a) At γ = 1/3 − 0.03, there is no gap-closing there, but a higher dispersive band moves down to approach the bottom flat band.(b) At γ = 1/3, a new gap-closing point appears as the higher dispersive band touches the bottom flat band.(c) As γ increases to γ = 1/3 + 0.03, the gap-closing point split into two and moves toward the conner of BZ.

FIG. 17 .
FIG. 17. Structure factor S(q) for the anisotropic honeycomb model defined by the constrainer Eq. (124).(a) At β = 1, the quadratic band touching is visible as standard pinch points in the structure factor.(b) As β is increased, two gap-closing points migrate along the Brillouin zone boundary toward the M point (edge center of the Brillouin zone).(c) At β = 2, two gap-closing point merge at the M point of the Brillouin zone, creating a parabolic pinch point singularity.(d) For β > 2 the system enters a trivial paramagnetic phase with smooth spin correlations throughout the Brillouin zone.

FIG. 18 .
FIG. 18.(a) Frustrated lattice composed of octahedra which share edges in the xy plane and connect via vertices in the z direction.There is a two site unit cell, with inequivalent sites here indicated in red and blue.Defining a local constraint on the octahedra leads to a classical spin liquid with nodal lines in J(q) and hence pinch lines in S(q).(b) Location of nodal lines for the model defined in Eq. (135).The nodal lines appear at wavevectors q = (qx, π, π), q = (π, qy, π) and equivalent, creating a network of nodal lines along the edges of the qz = ±π faces of the Brillouin zone.
FIG. 20.(a) The kagome lattice.(b) The constrainer of the kagome-star model.The constraint is defined on each hexagon of the lattice, with the spins on the interior of the hexagon (1 − 6) contributing to the constraint with coefficient 1, and the spins connected to the exterior (1 ′ − 6 ′ ) contributing with coefficient ζ (Eq.(151)-(152)).
) as a function of ζ reveals that pinch point singularities appear precisely at ζ = 1/2 and ζ = 1 (Fig.22) demonstrating the emergence of algebraic classical spin liquids where fragile topological classical spin liquids meet.

FIG. 21 .
FIG. 21.Evolution of momentum space Skyrmion number, Q sk , as a function of the tuning parameter ζ in The kagomestar model.Jumps in Q sk at ζ = 1/2 and ζ = 1 reveal zero temperature transitions between distinct, short range correlated, CSLs.Algebraic CSLs emerge at the boundaries.

FIG. 22 .
FIG.22.Evolution of spectrum and structure factor S(q), as a function of the tuning parameter ζ in the kagome-star model (Eq.(150)).S(q) is smooth throughout momentum space for generic values of ζ (panels (a), (b), (d)), indicating the short range correlations of the fragile topological CSLs.At the boundaries between these fragile topological CSLs algebraic CSLs appear, with gap-closing points in the spectrum and pinch point singularities in S(q) (panels (c), (e)).

TABLE II .
Some common algebraic classes of CLS

TABLE III .
Survey of known CSL models.See respective references for detailed definition of the models.