# Interacting-cluster spin liquids with robust flat bands evolving into higher-rank half-moon phases and topological Lifshitz transitions

Naïmo Davier, Ludovic D. C. Jaubert

► To cite this version:

Naïmo Davier, Ludovic D. C. Jaubert. Interacting-cluster spin liquids with robust flat bands evolving into higher-rank half-moon phases and topological Lifshitz transitions. Physical Review B, 2026, 113 (5), pp.054428. 10.1103/np33-ngyr. hal-05515885

HAL Id: hal-05515885

https://hal.science/hal-05515885v1

Submitted on 17 Feb 2026

HAL is a multi-disciplinary open access
archive for the deposit and dissemination of sci-
entific research documents, whether they are pub-
lished or not. The documents may come from
teaching and research institutions in France or
abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire HAL, est
destinée au dépôt et à la diffusion de documents
scientifiques de niveau recherche, publiés ou non,
émanant des établissements d'enseignement et de
recherche français ou étrangers, des laboratoires
publics ou privés.

<figure><img src="image_1.png" /></figure>

Distributed under a Creative Commons CC BY 4.0 - Attribution - International License
# Interacting-cluster spin liquids with robust flat bands evolving into higher-rank half-moon phases and topological Lifshitz transitions

Naïmo Davier\* and Ludovic D.C. Jaubert†

CNRS, Université de Bordeaux, LOMA, UMR 5798, 33400 Talence, France

(Dated: February 17, 2026)

Classical spin liquids are disordered magnetic phases, governed by local constraints that often give rise to flat-band ground states. When constraints take the form of a zero-divergence field within a cluster of spins, the spin liquid is often described by an emergent Coulomb gauge theory. Here we introduce an interaction $\eta$ between these clusters of spins which compete with the zero-divergence field. Using a framework embracing both the connectivity matrices of graph theory and the topology of band structures, we develop a generic theory of interacting-cluster Hamiltonians. We show how flat bands remain at zero energy up to finite interaction $\eta$, until a dispersive band becomes negative, stabilizing a spiral spin liquid with a hypersurface of ground-state manifold in reciprocal space. This hypersurface can be interpreted as an effective Fermi surface in the spectrum of the parent system, acting as a tunable energy selector despite the absence of particle filling. This effective Fermi surface serves as a mold for the apparition of the half-moon patterns in the equal-time structure factor. Our generic approach enables to extend the notion of half moons to the perturbation of higher-rank Coulomb fields and pinch-line spin liquids. In particular, multi-fold half moons appear when unconventional gauge charges, such as potential fractons, are stabilized in the ground state. Finally, half-moon phases can be tuned across the equivalent of a Lifshitz transition, when the hypersurface manifold changes topology.

## I. INTRODUCTION

Classical spin liquids emerge in geometrically frustrated magnets where conventional magnetic ordering is suppressed, even at zero temperature. Instead of long-range order, these systems exhibit highly degenerate ground-state manifolds governed by local constraints, which often take the form of flat bands in reciprocal space. When these constraints are written as a divergence-free condition – analogous to a Gauss law – they give rise to an emergent gauge structure reminiscent of classical electromagnetism [1, 2]. This leads to distinctive features such as algebraically decaying spin-spin correlations and characteristic pinch points in the structure factor, which are directly observable in neutron scattering experiments [3, 4]. Paradigmatic examples include the nearest-neighbor Heisenberg antiferromagnet on the pyrochlore and kagome lattices [5, 6].

Pinch points in the equal-time structure factor are often accompanied by half-moon patterns at finite energy coming from the lowest dispersive band [7–12]. Flat-band engineering [13–19] can modify the band spectrum in order to destabilize the spin liquid and make the half moons a signature of the ground state in the equal-time structure factor on kagome and pyrochlore systems [10, 13].

Using tools recently developed for the classification of classical spin liquids [20–24] together with graph-theory methods based on the connectivity matrix [10, 13, 25, 26], we build a generic theory of interacting-cluster systems for a seemingly infinite variety of lattices and clusters

in two and three dimensions (or higher). This interaction $\eta$ creates an energetic competition with the zero-divergence field required by the Coulomb gauge theory, and will be our flat-band engineering parameter. We quantify the stability of the flat-band spin liquid until $\eta$ reaches a critical value, where the energetic competition ultimately brings the lowest dispersive band below the flat bands. This crossing selects a hypersurface of parent-system states that are promoted to the ground state of the interacting system, defining an interaction-tunable effective Fermi level despite the absence of particle filling in spin systems. At this point, pinch points transform into half moons, the signature of a spiral spin liquid where the magnetic texture of gauge charges is stabilized in the ground state. Then we extend the notion of half moons to higher-rank U(1) gauge fields, including models with pinch lines. Since multifold half moons come from the dispersive band of a parent rank-$n$ Coulomb field, they represent the signature of rank-$n$ gauge charges, such as fractons, in the ground state. Finally, we show how half-moon phases can be adiabatically tuned across the equivalent of a Lifshitz transition separating distinct spiral spin liquids.

This paper is organized as follows. Section II defines cluster systems and characterizes them using the constraint-vector (IIB) and connectivity-matrix (IIC) formalisms. We then outline the Coulomb phase that emerges in these systems (IID) and discuss a simple criterion for when this description is expected to break down (IIE). Section III introduces the interacting-cluster Hamiltonian, presenting its general properties (III A) and the band structure of a subclass of interacting cluster systems (IIIB). In Section IV, we show how this band structure produces half-moon patterns in the structure factor when the cluster-cluster coupling $\eta$ is sufficiently

<footnote>\* naimo.davier@u-bordeaux.fr</footnote>

<footnote>† ludovic.jaubert@u-bordeaux.fr</footnote>
strong (IV A). This leads to the concepts of high-rank half-moons and half-moon surfaces associated with high-rank classical spin liquids (IV B), which we illustrate through two examples (IV C, IV D). Finally, we conclude with a description of the topological Lifshitz transitions that occur beyond the half-moon phase (IV E).

## II. CLUSTER HAMILTONIANS

### A. Definitions

Cluster Hamiltonians can be decomposed as a sum of energy terms over clusters $n$ [20–23],

$$ \mathcal{H} = \sum_n |\mathcal{C}_n|^2 \qquad (1) $$

where $\mathcal{C}_n$ is called a *constrainer*. This object can be seen as a weighted local magnetization for cluster $n$,

$$ \mathcal{C}_n = \sum_{i \in n} \gamma_i^n \mathbf{S}_i, \qquad (2) $$

where $\mathbf{S}_i$ are classical Heisenberg spins of unit length, $|\mathbf{S}_i| = 1$. All spins interact with each other in a given cluster via the coupling $2\gamma_i^n\gamma_j^n$, including a self-interacting term $(\gamma_i^n)^2$ that amounts to an irrelevant constant energy shift. There is a considerable variety of models, as there is a large freedom in choosing the geometry and size of each cluster, which can be connected to each other via a vertex, an edge, or a face, as illustrated in Table I.

The ground state of Hamiltonian (1) is obtained by minimizing the value of $|\mathcal{C}_n|$ for all clusters $n$ in the system. In this paper, we will consider models where it is possible to impose $\mathcal{C}_n = 0, \forall n$. Depending on the weights $\gamma_i^n$, this constraint is not always permitted, but such counter-examples usually lead to ordered ground states (see e.g. [36]).

Before introducing the interacting-cluster Hamiltonian and presenting our results, we shall now provide a concise review of the literature on the classification of classical spin liquids [20–24] and the connectivity matrix framework in spin liquids [10, 13, 25, 26].

### B. Flat-band ground state

Let us define the Fourier transform of the spin component $\alpha \in \{x, y, z\}$ on sublattice $\mu \in \{1, .., n_s\}$ in unit cell $i \in \{1, .., N_{u.c}\}$,

$$ \begin{aligned} S_{\mu}^{\alpha}(\mathbf{q}) &= \sum_{i} S_{\mu,i}^{\alpha} e^{-i(\mathbf{R}_{i}+\mathbf{r}_{\mu})\cdot\mathbf{q}}, \\ S_{\mu,i}^{\alpha} &= \frac{1}{N_{u.c}} \sum_{\mathbf{q}} S_{\mu}^{\alpha}(\mathbf{q}) e^{i(\mathbf{R}_{i}+\mathbf{r}_{\mu})\cdot\mathbf{q}}, \end{aligned} \qquad (3) $$

where $\mathbf{R}_i$ is the position of the unit cell and $\mathbf{r}_\mu$ is the relative position of the sublattice within the unit cell. Hamiltonian (1) is then rewritten in Fourier space [20]

$$ \mathcal{H} = \frac{1}{N_{u.c}} \sum_{\mathbf{q}} \sum_{X,\alpha} |\mathbf{L}_{X}^{\alpha}(\mathbf{q}) \cdot \mathbf{S}^{\alpha}(\mathbf{q})|^2 \qquad (4) $$

where $X \in \{1, .., n_c\}$ refers to the sum over the $n_c$ different types of clusters in the system; e.g. for the kagome lattice [Table. I], there are two types of clusters corresponding to the two types of triangles $\Delta, \nabla$. The constraint vector $\mathbf{L}_X^\alpha$, as introduced by Benton & Moessner [20], encodes the geometric structure of the clusters,

$$ (L_X^\alpha)_\mu(\mathbf{q}) = \sum_{j \in \mu \cap X} \gamma_j^X e^{i \mathbf{r}_j \cdot \mathbf{q}} \qquad (5) $$

where $\mathbf{r}_j$ denotes the vector position of site $j$ relative to the center of the cluster $X$, and the sum runs over all sites within cluster $X$ that belong to the sublattice $\mu$. Note that the dimension of vectors $\mathbf{S}^\alpha(\mathbf{q})$ and $\mathbf{L}_X^\alpha(\mathbf{q})$ is $n_s$, where each component corresponds to a different sublattice $\mu \in \{1, .., n_s\}$.

From now on, we shall work within the Luttinger-Tisza approximation (LTA) [37–39], where the spin length constraint is only enforced on average over the entire system, $\sum_i |\mathbf{S}_i| = N$ with $N$ the total number of spins, which are known to be suitable for cluster-system classical spin liquids, see Appendix C. A known limitation of the LTA is that it cannot always properly account for thermal order by disorder [27, 40], i.e. when thermal fluctuations entropically lift the classical ground-state degeneracy and favor magnetic order at very low temperatures. It is, however, possible to estimate when order by disorder is expected or not, i.e. when flat bands are expected to lead to a magnetically disordered spin-liquid ground state (see section II E). And even if the system ultimately orders at very low temperature $T_c$, the results of this paper remain valid at temperatures above $T_c$. Indeed, even if the LTA is an approximation, this method is particularly well suited to the study of classical spin liquids whose disordered magnetic texture accommodates well the averaged spin length constraint, as confirmed by many comparison between theory and simulations [10, 20, 22, 23, 41, 42]. And since we shall only consider *isotropic* exchange coupling here, the spin component $\alpha \in \{x, y, z\}$ is actually irrelevant and will be omitted from now on.

Once the Hamiltonian is written under the form (4), it becomes clear that all spin modes $\mathbf{S}_\perp(\mathbf{q})$ orthogonal to all constraint vectors $\mathbf{L}_X(\mathbf{q})$,

$$ \mathbf{L}_X(\mathbf{q}) \cdot \mathbf{S}_\perp(\mathbf{q}) = 0, \qquad X = 1, \dots, n_c \qquad (6) $$

are zero energy modes. When Eq. (6) is valid for all wavevectors $\mathbf{q}$, then the ground state is a flat band. The number of flat bands $n_{f.b}$ is thus equal to the total number of bands – i.e. the number of sublattices $n_s$ – minus
<table><thead><tr><th>2D Lattice</th><th>Kagome</th><th>Checkerboard</th><th>Hexagonal</th><th>Kagome<br>hexagonal</th><th>Square<br>octagon</th><th>Ruby</th><th>Square kagome</th><th>Decorated<br>square kagome</th><th>Octagonal<br>kagome</th></tr></thead><tbody><tr><th>Lattice<br>Scheme</th><td></td><td></td><td></td><td></td><td></td><td></td><td></td><td></td><td></td></tr><tr><th>n<sub>s</sub></th><td>3</td><td>2</td><td>2</td><td>3</td><td>4</td><td>6</td><td>6</td><td>6</td><td>14</td></tr><tr><th>n<sub>c</sub></th><td>2</td><td>1</td><td>1</td><td>1</td><td>2</td><td>3</td><td>4</td><td>5</td><td>9</td></tr><tr><th>n<sub>f,b</sub></th><td>1</td><td>1</td><td>1</td><td>2</td><td>2</td><td>3</td><td>2</td><td>1</td><td>5</td></tr><tr><th>F</th><td>0</td><td>1</td><td>1</td><td>3</td><td>2</td><td>3</td><td>0</td><td>-3</td><td>1</td></tr><tr><th>Premedial<br>Lattice<br>Scheme</th><td></td><td></td><td></td><td></td><td></td><td></td><td></td><td></td><td></td></tr><tr><th>N<sub>c,l</sub></th><td>3</td><td>2</td><td>3</td><td>3</td><td>6</td><td>6</td><td>6</td><td>10</td><td>18</td></tr><tr><th>References</th><td>[22, 23, 27]</td><td>[22, 23]</td><td>[20, 22]</td><td>[22, 23]</td><td>[28]</td><td>[29, 30]</td><td>[31]</td><td>[32]</td><td></td></tr></tbody></table>

<table><caption>Table I. Examples of possible cluster systems in two and three dimensions. The $n_c$ different types of clusters among a single lattice are depicted with different colors. The number of flat bands $n_{f,b}$ can be computed explicitly from the number of sublattices $n_s$ and the number of cluster type $n_c$ as $n_{f,b} = n_s - n_c$. The number of zero modes per unit cell $F$ can also be derived from these two numbers for Heisenberg spins as $F = 2n_s - 3n_c$. Each cluster type corresponds to a different sublattice of the premedial lattice [2], as depicted for the two-dimensional examples. There are $N_{c,l}$ inequivalent links in the premedial lattice, each one depicted with a different color (when their number allows for it). For the octagonal kagome lattice, the unit cell is composed of nine clusters and is highlighted with a dashed contour. In 3D, the kagome bipyramidal lattice is composed of by piramidal clusters with six faces, while the quadrupahedral lattice is composed of triangular-based pyramids touching through a face.</caption><thead><tr><th>3D Lattice</th><th>Pyrochlore</th><th>Octahedral</th><th>Kagome<br>bipyramidal</th><th>Quadrupahedral</th><th>Hyperkagome</th></tr></thead><tbody><tr><th>Lattice<br>Scheme</th><td></td><td></td><td></td><td></td><td></td></tr><tr><th>n<sub>s</sub></th><td>4</td><td>3</td><td>5</td><td>5</td><td>12</td></tr><tr><th>n<sub>c</sub></th><td>2</td><td>1</td><td>2</td><td>4</td><td>6</td></tr><tr><th>n<sub>f,b</sub></th><td>2</td><td>2</td><td>3</td><td>1</td><td>6</td></tr><tr><th>F</th><td>2</td><td>3</td><td>4</td><td>-2</td><td>6</td></tr><tr><th>N<sub>c,l</sub></th><td>4</td><td>3</td><td>5</td><td>16</td><td>9</td></tr><tr><th>References</th><td>[1, 22]</td><td>[20, 33]</td><td></td><td></td><td>[34, 35]</td></tr></tbody></table>

the number of independent constraints (6) per unit cell, which is given by the number of cluster types $n_c$, leading to $n_{f,b} = n_s - n_c$ [22, 23, 43]. The expected number of flat bands for various cluster systems is listed in Table I.

These flat bands solely depend on the geometry of the clusters – namely, their type, arrangement, and how they tile the lattice – but not on the specific values of the spin-spin interaction coefficients $\gamma_i^n$. These coefficients $\gamma_i^n$ influence the internal structure of the constraint vectors and thereby determine the detailed form of the dispersive bands, but do not generally affect these geometric flat bands. For finely tuned values, these coefficients $\gamma_i^n$ can accidentally flatten a dispersive band; a trivial example is when $\gamma_i^n = 0$ for all clusters $n$ belonging to a certain type of clusters $X'$ [26].

Since the number of constraints – and thus the flat-band structure – is governed by the way spins are grouped into interacting clusters, it is natural to expect that this information can also be encoded in the connectivity of the model itself.

### C. The connectivity matrix

The cluster Hamiltonian (1) can be expressed in terms of connectivity matrices [10, 25, 26], expressing it no more
<figure><img src="image_2.png" /><figcaption>Figure 1. Illustration of the mapping between the cluster Hamiltonian (1) and the connectivity matrices of Eq. (8). First, we consider a cluster at random, that is arbitrarily depicted as an hexagon here. Let us place a virtual site at its center; these central sites form the premedial lattice. Second, let us link all cluster sites i to the central site n (represented as curved oriented lines here). This is precisely the definition of the connectivity matrix $h^{v \leftarrow c}$ with coefficients $h_{in}^{v \leftarrow c} = \gamma_i^n$. Finally, let us link the center site n with all cluster sites i via the transpose of the connectivity matrix $h^{c \leftarrow v} = (h^{v \leftarrow c})^T$. The resulting matrix $H^{(1)} = h^{v \leftarrow c} h^{c \leftarrow v}$ thus pairs, within each cluster, all cluster sites two-by-two, and with the appropriate coupling constant as defined in Hamiltonian (1).</figcaption></figure>

as a sum over clusters, but rather as a sum over all lattice sites

$$ \mathcal{H} = \sum_{i,j} H_{i,j}^{(1)} \mathbf{S}_i \cdot \mathbf{S}_j. \qquad (7) $$

Since each spin interact with all other spins within the cluster, the interaction matrix $H^{(1)}$ can be simply expressed as

$$ H^{(1)} = h^{v \leftarrow c} h^{c \leftarrow v} \qquad (8) $$

where $h^{v \leftarrow c} = (h^{c \leftarrow v})^t$ denotes the $N \times N_c$ connectivity matrix that connects the $N$ lattice sites to the $N_c$ virtual cluster central sites that form the premedial lattice [2]. Its coefficients $h_{in}^{v \leftarrow c}$ are thus equal to $\gamma_i^n$ (see Eq. (2)) if the spin $i$ belongs to the cluster $n$, and 0 otherwise. The mapping of Eq. (8) is illustrated in Fig.1.

While there is a priori no theoretical limit on the number of clusters that can be considered for a given lattice, it might be somewhat artificial to have more clusters than spins in the system. One shall thus focus on models where $N > N_c$. Under this hypothesis, and since the interaction matrix $H^{(1)}$ is the product of two $N \times N_c$ rectangular matrices, the dimension of the kernel of $H^{(1)}$ cannot be smaller than $N - N_c = N_{u.c}(n_s - n_c)$ [26]. Hence, the connectivity matrix directly provides the same information about the number $(n_s - n_c)$ of flat bands in the

system than in the previous section (since each flat band contains $N_{u.c}$ points in reciprocal space). In other words, there are less constraints than total degrees of freedom, leaving some freedom for an extensive ground state manifold. As a side note, we see that the ground-state manifold of models with more clusters than spins, $N_c > N$, should a priori be over-constrained and are thus expected to be ordered.

Using the translation invariance of the Hamiltonian matrix (8), $H_{\mu,\nu}(\mathbf{R}_i, \mathbf{R}_j) = H_{\mu,\nu}(\mathbf{R}_j - \mathbf{R}_i)$, one gets the Fourier transform of Hamiltonian (7) [10]

$$ \mathcal{H} = \frac{1}{N_{u.c}} \sum_{\mathbf{q}} \sum_{\mu, \nu} H_{\mu\nu}^{(1)}(\mathbf{q}) \mathbf{S}_{\mu}(\mathbf{q}) \cdot \mathbf{S}_{\nu}(-\mathbf{q}), \qquad (9) $$

where

$$ H_{\mu\nu}^{(1)}(\mathbf{q}) = \sum_{j} H_{(0,\mu),(j,\nu)}^{(1)} e^{i(-\mathbf{R}_{j}+\mathbf{r}_{\mu}-\mathbf{r}_{\nu})\cdot\mathbf{q}} \qquad (10) $$

is the Fourier transform of the $n_s \times n_s$ Hamiltonian incident matrix for unit cell $j$ and between $\mu$ and $\nu$ sublattices. Using (8) to decompose

$$ H^{(1)}(\mathbf{q}) = h^{v \leftarrow c}(\mathbf{q}) h^{c \leftarrow v}(\mathbf{q}) \qquad (11) $$

where

$$ \begin{aligned} h_{\mu,X}^{v \leftarrow c}(\mathbf{q}) &= \sum_{j} h_{(0,\mu),(j,X)}^{v \leftarrow c} e^{i(-\mathbf{R}_{j}+\mathbf{r}_{\mu}-\mathbf{r}_{X}^{c})\cdot\mathbf{q}} \\ &= \sum_{k} h_{(k,\mu),(0,X)}^{v \leftarrow c} e^{i(\mathbf{R}_{k}+\mathbf{r}_{\mu}-\mathbf{r}_{X}^{c})\cdot\mathbf{q}}, \qquad (12) \\ &= [h^{c \leftarrow v}(\mathbf{q})]_{\mu,X}^{\dagger} \end{aligned} $$

is the $n_s \times n_c$ connectivity matrix in reciprocal space, with $\mathbf{r}_X^c$ denoting the position of the cluster center of type X among unit cell $j$. Here we recover the previous results. For $n_s > n_c$, the rectangular nature of the connectivity matrix implies at least $n_s - n_c$ flat bands associated with a zero energy. Furthermore, as the non-zero eigenvalues $\lambda_\mu(\mathbf{q})$ can be expressed as the square modulus of the singular values of $h_{\mu,X}^{v \leftarrow c}(\mathbf{q})$, i.e. $\lambda_\mu(\mathbf{q}) = |a_\mu(\mathbf{q})|^2$, all the dispersive bands are positive and thus located above the flat bands for cluster Hamiltonians. While the number of flat bands relies on the number of constraint – i.e. the number of cluster $n_c$ – per unit cell, the structure of the dispersive bands depends on the interaction strength associated with each bond, ruled by the coefficients $\gamma_i^n$ present in the connectivity matrix $h^{v \leftarrow c}$ and the constraint vectors, which are thus the proper tools to characterize the band structure of the system.

Finally, since $h_{(k,\mu),(0,X)}^{v \leftarrow c} \equiv \gamma_{k,\mu}^X$ in Eq. (12), one automatically obtains that the columns of the connectivity matrix in Fourier space are the constraint vectors $L_X(\mathbf{q})$

$$ h_{\mu,X}^{v \leftarrow c}(\mathbf{q}) = (L_X)_\mu(\mathbf{q}). \qquad (13) $$

A similar result has been obtained in [24] where the column matrix $[\mathbf{L}_1(\mathbf{q}), \mathbf{L}_2(\mathbf{q}), \dots, \mathbf{L}_X(\mathbf{q}), \dots]$ represents the hopping term between the two sublattices of a virtual bipartite construction, similar to the connectivity matrix.
## D. Band structure and Coulomb physics

It has been established [20–24] that in cluster systems exhibiting flat bands at the bottom of their spectrum, the constrainers and associated constraint vectors play a crucial role in determining the nature of the emergent spin liquid. Specifically, if the dimension of the vector space spanned by the constraint vectors –that is the one associated with dispersive bands– decreases at some momentum point corresponding to a dispersive band touching the set of flat bands, then it implies the emergence of a Coulomb phase. This phase manifests as a pinch point in the structure factor, reflecting long-range spin-spin correlations. The Coulomb phase is governed by a Gauss law acting on an effective electric field that emerges from the spin components. The nature of this emergent field is determined by the structure of the band-touching point. If the local dispersion near the contact point is of order $2n$, the system obeys a Gauss law of order $n$ acting on a rank-$n$ tensor field. Such higher-order laws are typically associated with $n$-fold symmetric pinch points in the structure factor [20–23, 44, 45].

While a full discussion of the various Coulomb phases is not the purpose of this work, we introduce here the main features associated with the band structure in order to be able to discuss the impact of cluster-cluster interactions on the emerging Coulomb phase. We refer the interested reader to [22–24] for a detailed presentation. The key point is that the nature of the Coulomb phase is entirely determined by the structure of the vector space associated with the dispersive bands, which itself is generated by the constraint vectors introduced in Eq. (4). These constraints are encoded in the connectivity matrix $h^{v \leftarrow c}(\mathbf{q})$, whose columns are the constraint vectors. Since ground state configurations must be orthogonal to all constraint vectors, they lie in the vector space orthogonal to the one spanned by the constraints—i.e., the space associated with the dispersive bands. Consequently, the spin-spin correlation functions in reciprocal space, encoded in the static structure factor $S(\mathbf{q})$, must be proportional [1] to the projection operator $\Pi$ onto this orthogonal space

$$S(\mathbf{q}) = \sum_{\mu, \nu} \langle \mathbf{S}_{\mu}(\mathbf{q}) \cdot \mathbf{S}_{\nu}(\mathbf{q}) \rangle \propto \sum_{\mu, \nu} \Pi_{\mu\nu}. \quad (14)$$

Because the Fourier-transformed connectivity matrix $h^{v \leftarrow c}(\mathbf{q})$ contains the constraint vectors as columns, the projector $\Pi$ can be constructed as [1, 23]

$$\Pi = I - h^{v \leftarrow c} \left( [h^{v \leftarrow c}]^{\dagger} h^{v \leftarrow c} \right)^{-1} [h^{v \leftarrow c}]^{\dagger}. \quad (15)$$

This projector becomes singular when the matrix $M \equiv [h^{v \leftarrow c}]^{\dagger} h^{v \leftarrow c}$ is non-invertible—i.e., when its determinant vanishes. Since this determinant is the product of the dispersive band dispersions, a zero determinant indicates that a dispersive band is touching a flat band. Thus, any such band touching point necessarily results in a pinch

point in the structure factor [21]. This occurs when one of the constraint vectors either vanishes or becomes linearly dependent on the others; see Appendix A for further details. At such touching point $\mathbf{q}^*$, there exists a critical vector $\mathbf{L}_c$ that can be expressed as a linear combination of the constraint vectors

$$\mathbf{L}_c = \sum_{X=1}^{n_c} \alpha_X \mathbf{L}_X, \quad \|\alpha\|^2 = 1 \quad (16)$$

satisfying $\mathbf{L}_c(\mathbf{q}^*) = 0$. Since all ground state configurations must satisfy the constraint equations

$$\mathbf{L}_X(\mathbf{q}) \cdot \mathbf{S}(\mathbf{q}) = \sum_{\mu=1}^{n_s} L_X^{\mu}(\mathbf{q}) \mathbf{S}_{\mu}(\mathbf{q}) = 0 \quad (17)$$

for every constraint vector $\mathbf{L}_X(\mathbf{q})$, the same must hold for the critical vector

$$\mathbf{L}_c(\mathbf{q}) \cdot \mathbf{S}(\mathbf{q}) = 0. \quad (18)$$

However, at the contact point $\mathbf{q}^*$, this last constraint becomes trivial since $\mathbf{L}_c(\mathbf{q}^*) = 0$. In the vicinity of this point, the critical vector can be expanded in powers of $\delta\mathbf{q} = \mathbf{q} - \mathbf{q}^*$, yielding the relation

$$\sum_{\mu} \sum_{j=0}^{n} a_{\mu,j} (\delta q_1)^j (\delta q_2)^{n-j} \mathbf{S}_{\mu}(\mathbf{q}_{\parallel}) = 0 \quad (19)$$

in planes $\delta\mathbf{q}_{\parallel} = (\delta q_1, \delta q_2)$ that contain the pinch point, where $a_{\mu,j}$ are the Taylor coefficients of $\mathbf{L}_c$ around $\mathbf{q}^*$. These equations correspond in real space to Gauss laws of order $n$

$$\sum_{\mu} \sum_{j=0}^{n} a_{\mu,j} (\partial_1)^j (\partial_2)^{n-j} \mathbf{S}_{\mu}(\mathbf{R}_{\parallel}) = 0, \quad (20)$$

applying on the coarse-grained spin fields, denoted by $\mathbf{S}_{\mu}(\mathbf{R}_i) = \mathbf{S}_{\mu, i}$. The order $n$ of the lowest nontrivial term in the expansion of the critical vector $\mathbf{L}_c$ determines the local dispersion of the touching band, which must scale as $|\delta\mathbf{q}|^{2n}$ (see Appendix A 3). This implies that the order of the Gauss law is dictated by the local structure of the band-touching, which in turn controls the $2n$-fold symmetry of the corresponding pinch point in the structure factor [22]. The constraint vectors thus provide a natural and effective framework for analyzing the emergent physics of cluster systems.

## E. When flat bands lead to spin liquids

Our discussion so far relies on the Luttinger–Tisza approximation (LTA). In order for the results based on the flat band states obtained within this framework to remain valid, these modes must be appropriately combined to construct real-space ground state configurations that satisfy the spin length constraint $|\mathbf{S}_i| = 1$ for each
spin. This construction requires a sufficient number of flat bands to ensure enough freedom when combining flat band modes into normalized real-space configurations. A natural question in this context is: how can one estimate – at least qualitatively – the minimal number of flat bands needed to preserve the extensive ground-state degeneracy and the corresponding Coulomb phase description ? This question can be addressed by evaluating, in real space, the effective number of zero modes per cluster or unit cell, this time incorporating the strict constraint of fixed spin length [23, 43].

The effective number of zero modes per unit cell F is the number of effective degrees of freedom per unit cell that remain unconstrained in the ground state. Each unit cell contains an effective number of sites that is equal to the number of sublattices $n_s$, together with a number of clusters $n_c$. On each of these sites there is a spin of dimension $n_d$ ($n_d = 3$ for Heisenberg spins), which possesses $n_d$ degrees of freedom but has a constrained length $|S| = 1$, and is thus associated with $n_d - 1$ effective degrees of freedom. As each cluster undergoes a constraint $C = 0$ of dimension $n_d$, this implies that there are $n_c n_d$ constraints per unit cells. The effective number of free degrees of freedom (or zero modes) per unit cell can thus simply be expressed as

$$F = n_s(n_d - 1) - n_c n_d. \qquad (21)$$

If the number of zero modes is strictly positive this means that there exists some free degrees of freedom within each clusters, allowing for an extensive degeneracy of the ground state manifold, in good agreement with the existence of flat bands in the spectrum. In fact this guarantees the number of flat bands is big enough to be able to build real space ground states from flat bands states that are satisfying the spin length constraint $|S| = 1$.

If the number of zero modes per unit cell is zero, as it is the case for the nearest-neighbor kagome antiferromagnet [43], this does not necessarily mean that there exists no zero modes, but that such modes do not rely on a single unit cell. For the kagome lattice for example, there exists zero modes called weather-vane modes that are relying on a group of six triangular clusters [27, 46]. In this marginal case, there exist some possibilities to combine flat band states to produce an extensive number of real space ground states satisfying the spin fixed length constraint.

If the number of zero modes obtained using Eq. (21) is negative, it means that recombining flat bands states to produce real space spins configurations while enforcing the spin length constraint amounts to losing the extensive degeneracy. In other words, restraining the LTA ground state manifold by imposing the spin length constraint reduces this manifold to a sub-manifold whose dimension does not scale with the system size; it likely leads to magnetic order at very low temperature. Note however that these systems can remain good candidates for quantum spin liquids [47, 48], as relaxing the spin length constraint can be seen as a way to incorporate quantum fluctuations

of the spins[49].

The good condition to look for isotropic cluster systems hosting classical spin liquid is thus to ask for $F \ge 0$, that is equivalent to ask the number of flat bands to satisfy

$$n_{f,b} = n_s - n_c \ge \frac{n_s}{3} \qquad (22)$$

for Heisenberg spins, and where $n_s$ is to be seen here as the total number of bands. For Heisenberg spins, this corresponds to systems where at least one third of the band are flat.

These discussions have highlighted why cluster-based systems are excellent candidates for hosting classical spin liquid phases. This naturally raises the question: to what extent can these systems be modified while preserving the intrinsic properties responsible for the emergence of flat bands and Coulomb phases? In the following section, we propose an extension of the cluster Hamiltonian that modifies both the interaction matrix while preserving the band structure that underlies the emergent classical spin liquid phase of the parent system.

### III. INTERACTING CLUSTERS HAMILTONIAN

Let us consider a generalization of the cluster Hamiltonian (1) with an additional term mediating an interaction between different clusters,

$$\mathcal{H} = \alpha \sum_n |\mathcal{C}_n|^2 + 2\eta \sum_{\langle m,n \rangle} \mathcal{C}_m \cdot \mathcal{C}_n. \qquad (23)$$

where the constraints $\mathcal{C}_n$ are defined as before by Eq. (2). For $\eta = 0$, one recovers the parent cluster Hamiltonian with $\alpha$ the intra-cluster interaction strength. The second sum runs over pairs of neighboring clusters with $\eta$ a real coefficient fixing the interaction strength of neighboring clusters. Specific examples of Hamiltonian (23) have been studied on the checkerboard, pyrochlore and kagome lattices [10, 50–54] where the first term was the standard nearest-neighbor antiferromagnetic term on tetrahedral and triangular units respectively, while the second term was due to second- and third-neighbor couplings between spins (with $J_2 = J_3$). It means that despite its apparent complexity due to its generic nature, relatively simple and realistic models can be described by Hamiltonian (23).

#### A. Generic properties

##### 1. Constraint vector framework

Let us define $z_c$ the coordination number of a spin cluster, i.e. the number of clusters linked to a given cluster. The interacting cluster Hamiltonian (23) can then
be rewritten as

$$ \mathcal{H} = (\alpha - z_c \eta) \sum_n |\mathcal{C}_n|^2 + \eta \sum_{\langle m, n \rangle} (\mathcal{C}_m + \mathcal{C}_n)^2 . \quad (24) $$

Note that if this coordination number is unique for most cluster lattices, it is not necessarily always true. One could take for example the square octagon lattice of Table I where octagonal cluster have coordination number $z_{c,o} = 8$ while square ones have $z_{c,s} = 4$. In this case the correct Hamiltonian would be

$$ \mathcal{H} = \sum_n (\alpha - z_{c,n}\eta) |\mathcal{C}_n|^2 + \eta \sum_{\langle m, n \rangle} (\mathcal{C}_m + \mathcal{C}_n)^2 . \quad (25) $$

This does not change the outcome of our discussion though, and we shall consider a unique $z_c$ for pedagogical reasons. In particular, Hamiltonian (24) makes it clear that configurations satisfying $\mathcal{C}_n = 0$ for all clusters $n$ in the system have zero energy, and that for $\alpha - z_c \eta > 0$ these are the ground state configurations. In reciprocal space we get

$$ \mathcal{H} = \frac{\alpha - z_c \eta}{N_{u.c}} \sum_{\mathbf{q}} \sum_{X} |\mathbf{L}_X(\mathbf{q}) \cdot \mathbf{S}(\mathbf{q})|^2 \\ + \frac{\eta}{N_{u.c}} \sum_{\mathbf{q}} \sum_{(XY)} |\mathcal{L}_{XY}(\mathbf{q}) \cdot \mathbf{S}(\mathbf{q})|^2 \quad (26) $$

where the sum $\sum_{(XY)}$ runs over the $N_{c,l}$ types of cluster to cluster links; see Table I for the value of $N_{c,l}$ for different lattices together with illustrations of the premedial-lattice sites depicting these different types of links. Cluster-pair constraint vectors $\mathcal{L}_{XY}(\mathbf{q})$ can be defined from cluster constraint vectors $\mathbf{L}_X$ as

$$ \mathcal{L}_{XY}(\mathbf{q}) = e^{-i\mathbf{r}_{XY} \cdot \mathbf{q}} \mathbf{L}_X(\mathbf{q}) + e^{i\mathbf{r}_{XY} \cdot \mathbf{q}} \mathbf{L}_Y(\mathbf{q}) \quad (27) $$

where the vector $\mathbf{r}_{XY}$ is half the vector that links the center of a cluster of type $X$ to a neighboring cluster of type $Y$. From the expression (26) of the Hamiltonian it appears that zero energy modes forming flat bands are spin modes that satisfy both

$$ \mathbf{L}_X(\mathbf{q}) \cdot \mathbf{S}(\mathbf{q}) = 0 \quad \& \quad \mathcal{L}_{XY}(\mathbf{q}) \cdot \mathbf{S}(\mathbf{q}) = 0. \quad (28) $$

Since the cluster-pair constraint vectors $\mathcal{L}_{XY}$ are simply linear combinations of cluster constraint vectors $\mathbf{L}_X$, this flat band manifold turns out to be identical to the one of the parent cluster Hamiltonian when $\eta = 0$, i.e. the vector space orthogonal to all cluster constraint vectors $\mathbf{L}_X$. This implies that the number of flat bands associated with zero energy, as well as their attached eigenstates, will remain the same when switching on the inter cluster interaction $\eta$. Note that this result remains valid even if we attach a specific coefficient $\eta_{XY} = \eta_{YX}$ for each type of cluster pairs $X-Y$ among the system.

As introduced in section IID, ground-state flat bands can support pinch points in the structure factor at wavevectors $\mathbf{q}^*$ where a dispersive band becomes gapless [21]. Since the dispersive bands are entirely defined

<figure><img src="image_3.png" /><figcaption>Figure 2. Schematic representation of Eq. (29). The scalar product of two constrainers associated with two neighboring clusters can be translated in term of connectivity matrices by introducing two different types of virtual bonds. The first types are the one already used for cluster Hamiltonians. These virtual bonds that link all vertices sites from a cluster to its virtual central site, are here depicted as black lines, and are associated with the connectivity matrix $h^{v \leftarrow c}$ as already discussed. The second type of virtual links are bonding the central sites from neighboring clusters together, they are depicted in red and can be encoded in the incident matrix $h^{c \leftarrow c}$ which coefficients $h_{mn}^{c \leftarrow c}$ are unity if the central sites $m$ and $n$ belong to two interacting clusters and zero otherwise.</figcaption></figure>

by the cluster constraint vectors, Eqs. (27) and (28) imply not only the persistence of flat bands up to $\eta < \alpha/z_c$, but also that the structure and number of pinch points, if any, are expected to be conserved. We shall now make this result more transparent using connectivity matrices. We shall focus on $\eta \ge 0$ even if many of our results also apply to the negative case.

## 2. Connectivity-matrix formalism

The connectivity-matrix framework can also be directly applied to analyze the cluster interacting Hamiltonian (23) as the second part of this Hamiltonian, proportional to $\eta$, naturally translates in this formalism as

$$ 2 \sum_{\langle m, n \rangle} \mathcal{C}_m \cdot \mathcal{C}_n \rightarrow H^{(2)} = h^{v \leftarrow c} h^{c \leftarrow c} h^{c \leftarrow v} \quad (29) $$

where the matrix $h^{c \leftarrow c}$ is a $N_c \times N_c$ connectivity matrix linking cluster centers to centers of neighboring clusters. Its coefficients $h_{mn}^{c \leftarrow c}$ are thus equal to 1 if the clusters $m$ and $n$ are neighboring clusters, and zero otherwise. The factor 2 in front of the left term comes from double counting when the connectivity matrix $H^{(2)}$ is inserted into Hamiltonian (7). The interacting-cluster Hamiltonian thus writes as

$$ H = h^{v \leftarrow c} (\alpha I_{N_c} + \eta h^{c \leftarrow c}) h^{c \leftarrow v} = h^{v \leftarrow c} g^{c \leftarrow c} h^{c \leftarrow v} \quad (30) $$

with $I_{N_c}$ the identity matrix of dimension $N_c$, and where $g^{c \leftarrow c}$ is a square matrix of dimension $N_c \times N_c$. Once expressed in Fourier space this Hamiltonian becomes

$$ H(\mathbf{q}) = h^{v \leftarrow c}(\mathbf{q}) g^{c \leftarrow c}(\mathbf{q}) h^{c \leftarrow v}(\mathbf{q}) \quad (31) $$
where

$$g_{X,Y}^{c \leftrightarrow c}(\mathbf{q}) = \sum_{j} g_{(0,X),(j,Y)}^{c \leftrightarrow c} e^{i(-\mathbf{R}_{j} + \mathbf{r}_{X}^{c} - \mathbf{r}_{Y}^{c}) \cdot \mathbf{q}}$$

is a $n_c \times n_c$ matrix. In this situation the general rank relation

$$\mathrm{Rank}(AB) \leq \mathrm{Min}(\mathrm{Rank}(A), \mathrm{Rank}(B)) \quad (32)$$

can be used to justify that, because $g_{X,Y}^{c \leftrightarrow c}(\mathbf{q})$ has a rank bounded by its dimension $n_c$, and as already stated $h^{v \leftrightarrow c}(\mathbf{q})$ has a rank inferior or equal to its smallest dimension $n_c$, the product (31) must also have a rank smaller or equal to $n_c$. This imposes that the Hamiltonian of an interacting cluster system must in general be associated with a minimal number $n_{f,b} = n_s - n_c$ of flat bands, as for the parent cluster system; a useful property for flat-band engineering of tight-binding models [13–16, 18, 19].

## B. Band dispersion of uniform interacting-cluster models

For a family of interacting-cluster models respecting the following conditions, the bands dispersions can be derived directly from the ones of the parent cluster system

(i) There is a unique type of cluster with $\xi$ spins per cluster and weighting coefficients $\{\gamma_i\}_{i=1,...,\xi}$.

(ii) Only nearest-neighbor clusters share spins, whose number $\omega$ is fixed: all clusters are either corner-sharing ($\omega = 1$), or bond-sharing ($\omega = 2$ a priori), or face-sharing ($\omega \ge 3$).

(iii) For any neighboring clusters $n$ and $n'$, the quantity

$$\Omega = \sum_{\substack{i=1 \\ i \in n \cap n'}}^{\omega} \gamma_i^n \gamma_{i'}^{n'}, \quad (33)$$

summed over all spins $i$ belonging to both clusters, is a constant.

Even if condition (iii) looks somewhat complex, it is actually a natural property for a cluster of spins, as long as equivalent spins within a cluster have the same weight $\gamma_i$. Our point here is that we do not need this full equivalence within a cluster, but simply a weaker version in the form of condition (iii). For convenience, we shall refer to a Hamiltonian respecting these three conditions as a uniform interacting-cluster model.

For such models, squaring the cluster Hamiltonian connectivity matrix leads to

$$\begin{align} (H^{(1)})^2 &= h^{v \leftrightarrow c} (h^{c \leftrightarrow v} h^{v \leftrightarrow c}) h^{c \leftrightarrow v} \nonumber \\ &= h^{v \leftrightarrow c} (\Omega h^{c \leftrightarrow c} + \Xi I_{N_c}) h^{c \leftrightarrow v} \tag{34} \\ &= \Omega H^{(2)} + \Xi H^{(1)} \nonumber \end{align}$$

as illustrated in Fig. 2. The identity contribution $\Xi I_{N_c}$ counts all two-step paths that start at a cluster center, visit a single vertex, and return to the same center. Since the path $n \to i \to n$ carries weight $\gamma_i^n$ on both hops, the coefficient is

$$\Xi = \sum_{i \in n} (\gamma_i^n)^2. \quad (35)$$

If condition (i) was not satisfied, this sum would depend on the cluster type, yielding a coefficient $\Xi_X$ for each type of cluster $X$, and the corresponding term in Eq. (34) would no longer be proportional to the identity.

The term $\Omega h^{c \leftrightarrow c}$ counts two-step paths connecting neighboring cluster centers via spins shared by both clusters. If conditions (ii) and (iii) were not satisfied, this coefficient would depend on the pair, $\Omega_{nn'}$, and the corresponding contribution in Eq. (34) would no longer be proportional to $h^{c \leftrightarrow c}$.

Using Eq. (34) to express $H^{(2)}$ in terms of $H^{(1)}$ and substituting into the Hamiltonian (23) yields

$$H = \left( \alpha - \eta \frac{\Xi}{\Omega} \right) H^{(1)} + \frac{\eta}{\Omega} \left( H^{(1)} \right)^2. \quad (36)$$

This connectivity matrix is a polynomial of the matrix $H^{(1)}$ and thus possesses the same eigenbasis as each eigenvector of $H^{(1)}$ is automatically an eigenvector of $H$. Denoting the eigenvalues of the parent cluster Hamiltonian $H^{(1)}$ as $\lambda_\mu(\mathbf{q})$ in reciprocal space, the energy spectrum of the interacting cluster Hamiltonian is

$$\Lambda_{\mu}(\mathbf{q}) = \left( \alpha - \eta \frac{\Xi}{\Omega} \right) \lambda_{\mu}(\mathbf{q}) + \frac{\eta}{\Omega} \lambda_{\mu}(\mathbf{q})^2. \quad (37)$$

This dispersion imposes the conservation of the flat bands manifold as $\lambda_\mu(\mathbf{q}) = 0$ naturally imposes $\Lambda_\mu(\mathbf{q}) = 0$, meaning the number of flat bands and their eigenvectors are conserved. This means that for this class of systems, not only the eigenbasis of the dispersive bands is conserved, but each of the Hamiltonian eigenspaces are individually conserved.

This derivation can be straightforwardly generalized to cluster systems composed of inequivalent cluster types and involving multiple kinds of geometrical intersections between clusters, such as the two variants of the square kagome lattice shown in Table I. Details of this generalization are presented in Appendix D.

The definition of Eq. (33) leaves the sign of $\Omega$ undefined a priori. But there is a hidden gauge degree of freedom in Hamiltonian (36) since $\Omega$ always appears as part of the ratio $\eta/\Omega$. We shall from now on consider that

$$\eta/\Omega > 0, \quad (38)$$

which corresponds to the interesting regime where the first term of Hamiltonian (36) can change sign as a function of $\eta$. Since we chose to study the regime $\eta \ge 0$, it means that we restrict our analysis to positive $\Omega$.
The dispersion relation of Eq. (37) reveals that in interacting-cluster systems, the lowest dispersive band becomes ground state when

$$ \Lambda_{ld}(\mathbf{q}) \le 0 \Leftrightarrow \lambda_{ld}(\mathbf{q}) \le \Xi \left(1 - \frac{\alpha}{\eta} \frac{\Omega}{\Xi}\right), \quad (39) $$

where $ld$ denotes the lowest dispersive band of the interacting-cluster model ($\Lambda_{ld}$) and of the parent model ($\lambda_{ld}$). The relevant control parameter is the ratio

$$ \zeta \equiv -\frac{\alpha}{\eta}, \quad \zeta_c = -\frac{\Xi}{\Omega} < 0, \quad (40) $$

so the lowest dispersive band becomes the ground state when $\zeta > \zeta_c$. If the parent model has a gap $\Delta$ between the flat band(s) and $\lambda_{ld}$, then the interacting system retains a flat-band ground state up to

$$ \zeta < \zeta_c(\Delta) = \left(1 - \frac{\Delta}{\Xi}\right) \zeta_c, \quad (41) $$

with $\zeta_c(0) = \zeta_c$. Thus a finite gap requires stronger inter-cluster coupling to destabilize the flat ground state, as detailed in Sec. IV E 2.

If the parent model is gapless ($\Delta = 0$), at the threshold $\zeta = \zeta_c$, the local band dispersion around the band-touching point $\mathbf{q}^*$ is squared compared to the rank order $n$ of the parent cluster Hamiltonian

$$ \Lambda_{ld}(\mathbf{q}^* + \delta\mathbf{q}) = \frac{\eta_c}{\Omega} \lambda_{\mu}(\mathbf{q}^* + \delta\mathbf{q})^2 \propto \frac{\eta_c}{\Omega} (\delta\mathbf{q})^{2n}, \quad (42) $$

at lowest order in $\delta\mathbf{q}$ and with $\eta_c = -\alpha/\zeta_c$. Naively, this observation would suggests that the critical point $\zeta = \zeta_c$ may be generically associated with the emergence of an exotic Coulomb phase, characterized by higher-rank Gauss laws and multifold pinch points (see explanations in section IID). However, the correspondence between the dispersion order $2n$ near a contact point and the order $n$ of the associated Gauss law crucially depends on the Hamiltonian being a Gram matrix constructed from the constraint vectors; we refer the interested reader to Appendix A for a technical explanation. This condition no longer holds in the case of the interacting cluster Hamiltonian, which instead takes the form of a second-order polynomial of such a Gram matrix. It means that the nature of the Coulomb phase remains unchanged for any $\zeta \le \zeta_c$, including $\zeta_c$, as long as the flat bands remain the lowest-energy states and the eigenvector space of the dispersive bands is preserved.

From Eqs. (24) and (26), one might naively conjecture that the critical value of $\zeta$ takes the form $\zeta_c = -z_c$, which in turn would suggest

$$ z_c = \frac{\Xi}{\Omega} = \frac{\sum_{i \in n} (\gamma_i^n)^2}{\sum_{i \in n \cap n'} \gamma_i^n \gamma_{i'}^{n'}}. \quad (43) $$

This identity can be enforced by fine-tuning the interaction weights $\gamma_i^n$, but it does not hold generically.

Nevertheless, Eq. (43) is valid under certain geometric conditions. Assuming $\gamma_i^n = \text{cst}$ for all sites $i$ and clusters $n$, the weights factor out of Eq. (43), yielding $\Xi/\Omega = \xi/\omega$. If, in addition, each site belongs to exactly two neighboring clusters, then $\xi = \omega z_c$ follows automatically, and hence $\zeta_c = -z_c$. This condition is naturally satisfied by all corner-sharing lattices (e.g., kagome and pyrochlore, as derived in [10], as well as trillium, hyperkagome, checkerboard, ruby); it also holds for bond-sharing lattices provided the two sites of each bond belong to exactly two clusters (e.g., the square-octagon lattice in Table I when only the octagonal clusters are considered).

Beyond these specific models, the failure of Eq. (43) in generic interacting-cluster systems indicates that the minima of the second terms in Eqs. (24) and (26) are not generally zero. Consequently, within the constraint-vector formalism one must explicitly evaluate

$$ \sum_{(XY)} |\mathcal{L}_{XY}(\mathbf{q}) \cdot \mathbf{S}(\mathbf{q})|^2 \quad (44) $$

to determine the threshold $\zeta_c$ above which the flat bands of the interacting-cluster system cease to be the ground state. This is illustrated by systems with a single type of cluster—and thus a single constraint vector—for which the interaction-matrix formalism allows an explicit computation of the sum. In such cases the cluster-pair constraint vectors are proportional to the unique constraint vector, implying

$$ \sum_{(XY)} |\mathcal{L}_{XY}(\mathbf{q}) \cdot \mathbf{S}(\mathbf{q})|^2 = b(\mathbf{q}) |\mathbf{L}(\mathbf{q}) \cdot \mathbf{S}(\mathbf{q})|^2, \quad (45) $$

with a real, non-negative coefficient $b(\mathbf{q})$ encoding the cluster-cluster connectivity. Combining this with Eq. (26) gives the dispersion

$$ \Lambda_{\mu}(\mathbf{q}) = (\alpha - z_c \eta) |\mathbf{L}(\mathbf{q})|^2 + \eta b(\mathbf{q}) |\mathbf{L}(\mathbf{q})|^2. \quad (46) $$

Since $\lambda(\mathbf{q}) = |\mathbf{L}(\mathbf{q})|^2$, a direct comparison between Eqs. (37) and (46) yields

$$ \sum_{(XY)} |\mathcal{L}_{XY} \cdot \mathbf{S}|^2 = \left( z_c - \frac{\Xi}{\Omega} + \frac{\lambda(\mathbf{q})}{\Omega} \right) |\mathbf{L} \cdot \mathbf{S}|^2. \quad (47) $$

Thus, this term does not, in general, have zero as its minimum value. Hence, one cannot infer the value of $\zeta_c$ from Eq. (24) alone. For uniform interacting-cluster systems where the dispersion relation (37) applies, the interaction-matrix formalism is therefore the appropriate framework to determine the critical parameter $\zeta_c$ governing the transition from pinch points to, as we shall now see, half-moons.

## IV. HALF MOONS

From now on, we shall focus on gapless cluster models where $\Delta = 0$ (with the exception of section IV E 2),
whose structure factor bears pinch points at $\alpha = 1, \eta = 0$ corresponding to the limit $\zeta \to -\infty$. What the previous discussion tells us is that these pinch points persist for $\zeta \le \zeta_c$, and something changes as $\zeta > \zeta_c$. The first term of Hamiltonian (24) is ultimately expected to break the $C_n = 0$ constraint as $\zeta$ increases, implying that the flat bands cannot be ground state anymore. In parallel, minimizing $|C_n + C_m|$ induces strong correlations between neighboring clusters that are known to be responsible for half-moon patterns in the equal-time structure factor of kagome and pyrochlore lattices with further-neighbor spin couplings [9, 10, 50-52]. Here, it is this formation of zero-energy half moons that we shall generalize beyond the canonical kagome and pyrochlore lattices and, noticeably, beyond regular two-fold pinch points.

But to describe half moons, we need to know what an equal-time structure factor is. We refer the reader to Appendix B for its definition and for a presentation of the Self-Consistent Gaussian Approximation [6, 55] (SCGA) used to compute it. For the sake of the paper, one simply needs to understand that, in the SCGA formalism, the static structure factor takes the form

$$S(\mathbf{q}) = \sum_{\kappa=1}^{n_s} \frac{W_{\kappa}(\mathbf{q})}{\lambda + \beta\varepsilon_{\kappa}(\mathbf{q})}, \quad (48)$$

where

$$W_{\kappa}(\mathbf{q}) = \sum_{\mu, \nu} [\psi_{\kappa}^{*}(\mathbf{q})]_{\mu} [\psi_{\kappa}(\mathbf{q})]_{\nu} \quad (49)$$

represent the weight of the eigenstates $\psi_\kappa$, and thus obey the sum rule over all bands $\kappa$

$$\sum_{\kappa=1}^{n_s} W_{\kappa}(\mathbf{q}) = n_s. \quad (50)$$

### A. General theory of half moon patterns on uniform interacting-cluster systems

It is possible to derive a theory for the apparition of half-moon patterns in the structure factor of uniform interacting-cluster models thanks to the known form of the dispersive band given by Eq. (37). For $\zeta > \zeta_c = -\frac{\Xi}{\Omega}$, the energy minima respects the equation

$$\nabla_{\mathbf{q}} \Lambda_{ld} = \left[ \alpha + \frac{\eta}{\Omega} (-\Xi + 2\lambda_{ld}) \right] \nabla_{\mathbf{q}} \lambda_{ld} = 0. \quad (51)$$

The condition $\nabla_q \lambda_{ld} = 0$ is the same as for the parent cluster Hamiltonian ($\alpha = 1, \eta = 0$) where flat bands were ground state; its physics thus does not apply here where $\zeta > \zeta_c$. On the other hand, canceling the prefactor of Eq. (51) defines the condition on $\lambda_{ld}$ for the energy minimum of the interacting-cluster system. In reciprocal

space, the position of this energy minimum is given by the manifold $Q$

$$\{ \mathbf{q} \in Q \mid \lambda_{ld}(\mathbf{q}) = \frac{\Xi}{2} \left( 1 - \frac{\zeta}{\zeta_c} \right) \equiv \lambda_c \}. \quad (52)$$

Since we are working on a gapless parent cluster model ($\Delta = 0$), the eigenvalue $\lambda_{ld}$ is greater or equal than 0, and the manifold $Q$ is thus properly defined for $\zeta > \zeta_c$. On this manifold, the minimum energy of the interacting-cluster model is always negative, i.e. below the flat band,

$$\Lambda_{ld}(\mathbf{q} \in Q) = -\frac{\alpha \Xi \zeta_c}{4 \zeta} \left( \frac{\zeta}{\zeta_c} - 1 \right)^2 < 0. \quad (53)$$

Since $\lambda_{ld}(\mathbf{q}^*) = 0$ by definition of the band-touching point $\mathbf{q}^*$, $\mathbf{q}^*$ remains a band-touching point for all values of $\zeta$, even when the dispersive band $\Lambda_{ld}$ becomes negative. By analytical continuity of Eq. (37), the manifold $Q$ typically forms a closed contour in 2D or a closed surface in 3D, in other words a hypersurface, around $\mathbf{q}^*$ in reciprocal space. The evolution of the band structure as $\zeta$ increases can be understood visually by noting that the dispersion relation in Eq. (37) is given by a polynomial $\Lambda(\lambda)$ acting on the band structure of the parent system, as illustrated in Fig. 3. The spectrum of the interacting system is therefore entirely determined by the behavior of $\Lambda(\lambda)$ over the spectrum of the parent system. In particular, the minimum $\lambda_c$ of the polynomial $\Lambda(\lambda)$ plays the role of an effective Fermi level. Indeed, in spin systems there is no notion of filling, and the low-temperature behavior is therefore fully determined by the ground state. However, the introduction of cluster-cluster interactions allows parent-system states at energy $\lambda_c$ to become ground states of the interacting system. This mechanism enables one to probe the parent-system energy states located around $\lambda_c$. In this sense, $\lambda_c$ acts as an effective Fermi level that can be continuously tuned by varying the strength of the cluster-cluster interactions controlled by $\zeta$.

At low temperature, the eigenstates attached to $\Lambda_{ld}(\mathbf{q} \in Q)$ dominate the spin-spin correlations responsible for the structure factor. If the density of these eigenstates were uniform along $Q$, one would obtain an approximately circular and spherical contour in 2D and 3D respectively. However, the sum rule imposes that the density of states must be uniform once summed up over all bands. Hence, since the highest intensity region of the flat band(s) is the bow-tie shape of the pinch point, the highest intensity region of the dispersive band(s) must be the complementary region of the bow tie in reciprocal space [9, 10]. Ultimately, it is the intersection between the complementary region of the bow tie and the hypersurface $Q$ that forms the half-moon patterns of the structure factor.

In absence of extensively degenerate flat bands, this half-moon phase is not a traditional spin liquid. But since the manifold $Q$ is sub-extensive and isomorphic to
<figure><img src="image_4.png" /><figcaption>Figure 3. Effective Fermi-level illustration. The band structure of an interacting-cluster system can be written as a polynomial function $\Lambda(\lambda)$ of the band structure of the parent cluster system (see Eq. (37)). (a) Since the parent spectrum contains only positive energies, a positive value of the linear coefficient $\alpha - \eta\Xi/\Omega$ results in a simple, non-uniform rescaling of the parent band structure. (b-c) When the linear coefficient $\alpha - \eta\Xi/\Omega$ becomes negative, the polynomial ceases to be strictly increasing, and the minima of the polynomial select the intersections with the parent band structure as the new ground-state manifold. (b) Just beyond the transition, the new ground-state manifold $Q$ consists of closed contours encircling the contact points of the parent system. (c) When the polynomial minimum $\lambda_c$ crosses a saddle point of the parent-system band structure, the topology of the ground-state manifold changes, corresponding to a Lifshitz transition. Note that the flat bands of the parent system (depicted in red), which are located at zero energy, remain pinned at zero energy in the interacting-cluster spectrum, since the polynomial $\Lambda(\lambda)$ has no constant term.</figcaption></figure>

an $(d-1)$-sphere within a given Brillouin zone, the half-moon phase corresponds to a so-called spiral spin liquid [56–61] with unusual, non-local, spin dynamics [62]. At very low temperatures, beyond the Luttinger–Tisza approximation, thermal order-by-disorder (ObD) may entropically lift the degeneracy of the manifold $Q$, particularly in $d=3$ dimensions. In $d=2$ dimensions, however, the Hohenberg–Mermin–Wagner theorem forbids the emergence of long-range magnetic order, although the breaking of discrete lattice rotational symmetries remains possible [63]. Such symmetry breaking, potentially associated with a Berezinskii–Kosterlitz–Thouless–type transition in the case of XY spins, can only occur at very low temperatures. At intermediate temperatures, thermal fluctuations are instead expected to restore the effective degeneracy of $Q$, even slightly away from fine tuning. As a consequence, the half-moon features in the equal-time structure factor are expected to remain robust above the temperature regime where ObD takes place.

In any case, for $\zeta > \zeta_c$, the energy minimum of the interacting-cluster Hamiltonian ($\Lambda_{ld}$) comes from the lowest excited dispersive band of the parent cluster Hamiltonian ($\lambda_{ld}$). Since the flat bands support an emergent Coulomb gauge theory, the lowest dispersive band supports a gapless form of the gauge-charge excitations of the Coulomb theory [9, 10]. *From a real-space perspective*, the appearance of half-moon features therefore signals that these gauge charges are no longer merely thermal excitations but become an integral part of the ground-state manifold. For two-fold half moons, these gauge charges take the form of magnetic monopoles, cor-

responding to flux defects that locally violate the flux-conservation constraint imposed by $C=0$. Note that for spins of dimension $n_s$, there are $n_s$ equivalent copies of such scalar gauge charges. In simple systems where the constrainer corresponds to the local cluster magnetization, combinations of the $n_s$ types of gauge charges can be interpreted as effective magnetic moments localized on individual clusters. For multifold half moons, the situation is more subtle. In this case, the emergent gauge field is no longer vectorial, and the corresponding gauge charges cannot be interpreted as simple magnetic dipoles. Nevertheless, they remain associated with local violations of a tensorial Gauss law, and thus correspond to higher-rank, spatially structured charge defects in real space. These objects represent generalized gauge charges whose correlations define the characteristic magnetic textures underlying multifold half-moon patterns. The type of flat-band engineering we have introduced via the cluster interaction $\eta$ thus stabilizes magnetic textures where gauge charges have become part of the ground state. At low temperature, fluctuations are restricted to perturbations close to the hypersurface manifold. But at intermediate temperatures, when the flat bands start to be populated, there is a temperature window where a dense phase of gauge charges co-exist with the vacuum of the flat bands. Although this regime may appear reminiscent of a conventional Coulomb phase, the cluster-cluster interaction introduces an additional dominant interaction between gauge charges. This interaction sets a characteristic length scale for their spatial organization, directly related to the inverse radius of the ground-state
hypersurface. Because spins are shared between neighboring clusters, it is generally impossible to simultaneously maximize the norms of oppositely oriented neighboring constrainers. This geometric frustration leads to nontrivial charge arrangements in real space, where a maximally charged cluster is typically surrounded by oppositely charged clusters carrying smaller individual charges, as previously observed in Ref. [10]. By varying the parameter $\eta$, one can include more or less of the dispersive band below the flat bands, i.e. qualitatively more or less gauge charges in the low-energy regime.

For models where the thermal order-by-disorder selects a long-range magnetic order at very low temperatures, the ObD transition would correspond to the crystallization of these gauge charges [10, 64].

## B. How to extend the theory to higher-rank half moons

As presented in section II D, higher-rank U(1) spin liquids support emergent Coulomb gauge fields of rank $-n$ where the emergent electric field can be a vector ($n=1$), a matrix ($n=2$), or a tensor of higher rank [65]. Their higher-rank structure confers additional conserved quantities on the spin texture, where excitations may naturally take the form of fractons [66]. The most famous signature of higher-rank U(1) spin liquids is probably the $2n$-fold pinch points in the structure factor [44, 45]. At the band-touching point $\mathbf{q}^*$, the bow-tie shape of regular pinch points splits into $2n$ branches.

To the best of our knowledge, half-moon patterns have only been studied out of regular (rank $-1$) Coulomb gauge theories [9, 10]. It is thus tempting to see if our approach could account for multifold half moons, and the physics of higher-rank gauge charges stabilized at zero temperature.

There is, however, an issue. In absence of anisotropic exchange [45, 67-69], higher-rank spin liquids are only known on extended cluster systems [20-23], i.e. where there is a large overlap between adjacent clusters, and spins belong to more than 2 clusters. This is because the higher-rank nature requires to cancel the first ($n-1$) order terms in the expansion of each constraint vector component $\mathbf{L}_\mu(\mathbf{q})$ (see Eq. (19)). The only degrees of freedom available to enable this cancellation are the couple of parameters ($\gamma_{i,\mu}$, $\mathbf{r}_{i,\mu}$) associated with each site belonging both to the cluster and to sublattice $\mu$. To be able to enforce this cancellation for $n > 1$ thus requires to have multiple sites from a same sublattice encapsulated in a same cluster, corresponding to the case of an extended cluster.

The inconvenience is that extended-cluster systems do not satisfy condition (ii) of section III B. Consequently, their lowest dispersive band $\Lambda_{ld}$ does not, in general, take the convenient polynomial form of Eq. (37). We usually

have

$$ \Lambda_{ld}(\mathbf{q}) = \alpha\lambda_{ld}(\mathbf{q}) + \eta g(\mathbf{q}) \quad (54) $$

with a given function $g(\mathbf{q})$. The minimum of $\Lambda_{ld}(\mathbf{q})$ has no reason to form a closed contour in reciprocal space, nor to be related to the band-touching point $\mathbf{q}^*$ of the pinch point. Instead its gradient a priori cancels on a discrete set of $\mathbf{q}$ points only, since $\nabla_{\mathbf{q}}\Lambda_{ld} = 0$ corresponds to a set of three equations with three unknown variables ($q_x, q_y, q_z$). The intersection between the bow tie of a pinch point and the manifold of minimum energy is not expected to produce any remarkable feature. A priori, this would suggest to preclude observing higher-rank half moons within the simple interacting-cluster models governed by Hamiltonian (23).

Fortunately, there is a way to circumvent the problem. The strategy is to engineer a model in which the lowest dispersive eigenvalue of the interacting system, $\Lambda_{ld}$, becomes a simple function of the parent eigenvalue $\lambda_{ld}$. For simplicity, let us consider parent clusters that satisfy condition (i), i.e. that admit a single constraint vector $\mathbf{L}$ and a single set of weights $\gamma_i^n = \gamma_i \forall n$. We then seek an interacting model with a dispersive band $\Lambda$ that is a quadratic polynomial of $\lambda \equiv |\mathbf{L}|^2$,

$$ \Lambda(\mathbf{q}) = [\alpha + \eta (a + b\lambda(\mathbf{q}))] \lambda(\mathbf{q}) \quad (55) $$

with $a$ and $b$ some real scalars coefficients. In other words, the function $g(\mathbf{q})$ of Eq. (54) would be a linear function of $\lambda(\mathbf{q})$, and the Fourier transform of the Hamiltonian would take the form

$$ \mathcal{H} = \frac{1}{N_{u.c}} \sum_{\mathbf{q}} \left[ \alpha + \eta (a + b|\mathbf{L}|^2) \right] |\mathbf{L}|^2 |\mathbf{1} \cdot \mathbf{S}|^2, \quad (56) $$

where $\mathbf{l}(\mathbf{q}) = \mathbf{L}(\mathbf{q})/|\mathbf{L}(\mathbf{q})|$ is the dispersive-band eigenvector. We shall soon explain why this type of Hamiltonian supports higher-rank half moons in the equal-time structure factor. But before that, let us justify our strategy by explaining how to design a model with Hamiltonian (56). To this end, we generalize the model with cluster interactions beyond nearest neighbors, that can vary between different pairs

$$ \mathcal{H} = \alpha \sum_n |\mathbf{C}_n|^2 + \eta \sum_{n,p} K_{np} \mathbf{C}_n \cdot \mathbf{C}_p, \quad (57) $$

with translation invariance of the scalar kernel $K_{np} = K_{\mathbf{R}_{pn}}$. $\mathbf{R}_{pn}$ is the vector connecting the centers of clusters $n$ and $p$. Since we only consider one type of cluster, there must be a cluster associated with each unit cell (by definition of a unit cell). $\mathbf{R}_{pn}$ thus also connects the centers of unit cells. In Fourier space, Hamiltonian (57) becomes

$$ \mathcal{H} = \frac{1}{N_{u.c}} \sum_{\mathbf{q}} \left[ \alpha + \eta K(\mathbf{q}) \right] |\mathbf{L}(\mathbf{q}) \cdot \mathbf{S}(\mathbf{q})|^2, \quad (58) $$
with

$$K(\mathbf{q}) = \sum_{p} K_{\mathbf{R}_{cp}} e^{-i\mathbf{q} \cdot \mathbf{R}_{cp}} \qquad (59)$$

summed over all clusters p with respect to a reference cluster c. Since the system is translation invariant, $K(\mathbf{q})$ is independent of the reference cluster. The simplest way to reach the target form of Eq. (56) is to choose $K_{\mathbf{R}}$ such that

$$K(\mathbf{q}) \equiv |\mathbf{L}(\mathbf{q})|^2, \qquad (60)$$

where $a=0$ and $b=1$. This is achieved by

$$K_{\mathbf{R}} = \sum_{\mu} \sum_{i \in c \cap \mu} \sum_{j \in c \cap \mu} \gamma_{i, \mu} \gamma_{j, \mu} \delta_{\mathbf{R}, \mathbf{r}_{i, \mu} - \mathbf{r}_{j, \mu}} \quad (61)$$

for which

$$\begin{aligned} K(\mathbf{q}) &= \sum_{\mu} \sum_{i \in c \cap \mu} \sum_{j \in c \cap \mu} \gamma_{i, \mu} \gamma_{j, \mu} e^{i(\mathbf{r}_{i, \mu} - \mathbf{r}_{j, \mu}) \cdot \mathbf{q}} \\ &= \sum_{\mu} \left( \sum_{i \in c \cap \mu} \gamma_{i, \mu} e^{i\mathbf{r}_{i, \mu} \cdot \mathbf{q}} \right) \left( \sum_{j \in c \cap \mu} \gamma_{j, \mu} e^{-i\mathbf{r}_{j, \mu} \cdot \mathbf{q}} \right) \\ &= |\mathbf{L}(\mathbf{q})|^2. \end{aligned}$$

This construction relies on the requirement that, within a given sublattice, all relative position vectors $\mathbf{r}_{i, \mu} - \mathbf{r}_{j, \mu}$ are lattice translations $\mathbf{R}_n$, which is naturally satisfied when clusters encapsulate only vertex spins. The resulting interactions between pairs of clusters are illustrated in Fig. 4(b) for the example of the generalized checkerboard model.

The resulting dispersive band of Hamiltonian (57) is given by Eq. (55): $\Lambda(\mathbf{q}) = [\alpha + \eta\lambda(\mathbf{q})] \lambda(\mathbf{q})$, whose minimum is either the same as $\lambda(\mathbf{q})$ (i.e. at the band-touching point $\mathbf{q}^*)$ or for $\lambda(\mathbf{q}) = \zeta/2$. Since $\lambda(\mathbf{q})$ is always positive (or zero), the second solution only exists for $\zeta > 0 \Rightarrow \alpha < 0$. While negative values of $\alpha$ will be considered in a following section, let us push our reasoning a little further to see if we can find half moons while preserving the positivity of $\alpha$.

The problem is that condition $K(\mathbf{q}) = |\mathbf{L}(\mathbf{q})|^2$ is too simple. The Ansatz of Eq. (61) includes a troublesome, and somewhat unphysical, self-interaction which cancels out the $\eta$-contribution to the prefactor of $\sum_n |C_n|^2$ in the Hamiltonian. Going from Hamiltonian (23) to its rewritten form (24), the prefactor of $\sum_n |C_n|^2$ changes from $\alpha$ to $(\alpha - z_c\eta)$. Thanks to this negative term, $-z_c\eta$, the $C_n = 0$ condition ultimately breaks for large enough $\eta$, which is what makes the dispersive band ground state. Since the Ansatz of Eq. (61) does not “renormalize” $\alpha$, the flat band always remains ground state. This problem is easily solved as it is enough to remove the on-site ($\mathbf{R} = 0$) term from the Ansatz.

$$\begin{aligned} K_{\mathbf{R}} &= \frac{1}{\omega} \left( \sum_{\mu} \sum_{i,j \in c \cap \mu} \gamma_{i,\mu} \gamma_{j,\mu} (\delta_{\mathbf{R}, \mathbf{r}_{i,\mu} - \mathbf{r}_{j,\mu}} - \delta_{\mathbf{R},0}) \right) \\ &= \sum_{\mu} \sum_{i,j \in c \cap \mu} \frac{\gamma_{i,\mu} \gamma_{j,\mu}}{\omega} \delta_{\mathbf{R}, \mathbf{r}_{i,\mu} - \mathbf{r}_{j,\mu}} - \frac{\Xi}{\omega} \delta_{\mathbf{R},0}. \end{aligned} \qquad (62)$$

The division by $\omega$ in the above Ansatz is not necessary, but it makes for easier comparison with previous results on models with nearest-neighbor cluster interaction only, as it avoids multiple countings of the same nearest-neighbor link. Here $\omega$ is the number of spins shared between adjacent unit cells. In Fourier space we get

$$K(\mathbf{q}) = \frac{|\mathbf{L}(\mathbf{q})|^2}{\omega} - \frac{\Xi}{\omega}. \qquad (63)$$

and

$$\mathcal{H} = \frac{1}{N_{u.c}} \sum_{\mathbf{q}} \left[ \alpha + \frac{\eta}{\omega} (-\Xi + |\mathbf{L}|^2) \right] |\mathbf{L}|^2 |\mathbf{l} \cdot \mathbf{S}|^2, \quad (64)$$

which matches Eq. (47) for simple systems. This generalized interacting-cluster model thus has a dispersive band

$$\Lambda(\mathbf{q}) = \left[ \alpha + \frac{\eta}{\omega} (-\Xi + \lambda(\mathbf{q})) \right] \lambda(\mathbf{q}). \qquad (65)$$

If we redefine the critical parameter as

$$\zeta_c = -\frac{\Xi}{\omega}, \qquad (66)$$

then, using the same reasoning as for the simple interacting-cluster model, the ground-state manifold for $\zeta > \zeta_c$ again forms a hypersurface,

$$\left\{ \mathbf{q} \in \mathcal{Q} \mid \lambda(\mathbf{q}) = \lambda_c = \frac{\Xi}{2} \left( 1 - \frac{\zeta}{\zeta_c} \right) \right\}, \qquad (67)$$

surrounding the band-touching point $\mathbf{q}^*$ of the parent model. We obtain a spiral spin liquid, and since this strategy is valid for parent models made of extended clusters that can support multifold pinch points at $\mathbf{q}^*$, the generalized interacting-cluster Hamiltonian can support multifold half moons, as illustrated in the next section.

Half moons can thus take exotic, multifold, shapes indicating the proximity of higher-rank U(1) gauge fields in parameter space. The lines of minimum intensity in these multifold half moons always mark the position of the pinch-point branches in the related tensorial spin liquid, which has an exciting consequence. Multifold half moons belong to the lowest excited dispersive band of the parent cluster Hamiltonian ($\lambda_{ld}$). Since the flat bands support a higher-rank Coulomb gauge theory, the lowest dispersive band supports a gapless form of the higher-rank gauge-charge excitations, which, in some models, are expected to take the form of fractons. In that sense, multifold half moons are the smoking gun for fractons in the ground state, albeit confined to the Q-manifold vicinity. As discussed in the previous section for regular half moons, one can expect here a regime of higher-rank gauge charges of tunable density, co-existing with the vacuum of the flat bands at intermediate temperatures. And if order by disorder takes place, which is probable in three dimensions, then one could observe the crystallization of fractonic matter at very low temperature.
<figure><img src="image_5.png" /><figcaption>Figure 4. (a) Cluster structure of the generalized checkerboard model. The plain and crossed circles depict the two sublattices, while blue, red and yellow circles respectively appear with coefficient 1, $\gamma_1$ and $\gamma_2$ in the constrainer $C \square$ of Eq. (68). (b) Effective strength of the cluster interactions associated with the kernel $K_R$ (62) for the generalized checkerboard, assuming $\gamma_2 = 0$ for clarity of the figure. The vector between two clusters is defined as $R_{np}$. The strength of each cluster interaction $(n, p)$ is determined by all the combinations of pairs of sites $i$ and $j$ belonging to sublattice $\mu$, whose positions in the cluster respect $r_{i,\mu} - r_{j,\mu} = R_{np}$. Each of these combinations is associated with the products of the weights $\gamma_{i,\mu}\gamma_{j,\mu}/\omega$, with $\omega = 1$ here. Since the two sublattices are equivalent, let us apply our reasoning on one arbitrary sublattice. One has a weight 1 associated with vectors $\pm e_1$, and $\gamma_1$ associated with vectors $\pm u_1, \pm u_2$. The weight of a cluster-to-cluster link is thus always a combination of terms proportional to either $1^2, 1 \times \gamma_1$ or $(\gamma_1)^2$ as depicted in blue. Because of the two possible orientations of the bond vector $R_{np}$ each cluster-to-cluster link weight must be proportional to two. Note that some links weights can get contributions from two sublattices, as it is the case for the bonds linking first or second neighbors clusters.</figcaption></figure>

<figure><img src="image_6.png" /><figcaption>Figure 5. Generalized checkerboard model Figures (a,b) are obtained with constrainer parameters $\gamma_1 = -1/2, \gamma_2 = 0, \alpha = 1$, which fixes $\Xi = 6$ and thus $\zeta_c = -6$. The parameter $\zeta$ is thus set to $\zeta = -4 > \zeta_c$. The figure (a) presents the structure factor, showing fourfold half moons, in good agreement with figure (b) representing the weight $W_2(q)$ of the dispersive band and the minimum of this dispersive band depicted as a red line. Figures (c-d) are obtained with $\gamma_1 = 1, \gamma_2 = 1/3$, imposing $\Xi = 112/9$ and $\zeta_c = -112/9 \simeq -12.4$, suggesting the choice $\zeta = -10$. The structure factor (c) presents sixfold half moons, in good agreement with the weight $W_2(q)$ depicted on panel (d) that shows sixfold pinch points. For both cases calculations of the structure factors are performed via SCGA (see Appendix B) with inverse temperature $\beta = 20$ and with the $q_i$ in $a^{-1}$ units.</figcaption></figure>

### C. Example: higher-rank half moons on the generalized checkerboard model

The generalized checkerboard model supports several distinct spin liquids [23], with the following constrainer

$$C \square = \sum_{i \in \square} S_i + \gamma_1 \sum_{i \in \langle \square \rangle} S_i + \gamma_2 \sum_{i \in \langle \square \rangle} S_i. \quad (68)$$

as illustrated in Fig. 4. This model hosts fourfold pinch points in the structure factor when $\gamma_2 = -1 - 2\gamma_1$, and

sixfold pinch points for $\gamma_2 = (1 - 2\gamma_1)/3$. It is therefore well suited to study how higher-rank pinch points evolve into half-moon patterns for $\zeta > \zeta_c$.

The extended-cluster geometry imposes $\Xi = 4(1 + 2\gamma_1^2 + \gamma_2^2)$ and $\omega = 1$. At low temperature, for $\zeta > \zeta_c$ with $\gamma_1 = -1/2$ and $\gamma_2 = -1 - 2\gamma_1 = 0$, the structure factor exhibits fourfold half-moons [Fig. 5(a)]. The corresponding dispersive-band weight $W_2(q)$ [Fig. 5(b)] shows fourfold pinch points that are the complement of the flat-band weight $W_1(q)$ (see the dark regions in Fig. 5(b)). Because the dispersive band attains local
<figure><img src="image_7.png" /><figcaption>Figure 6. Generalized octochore model: The generalized interacting cluster Hamiltonian is set with parameters, $γ_1 = -1/2$, $γ_2 = 1$ to sit at the parameter point with a pinch-line singularity when $η = 0$. In this case $Ξ = 18$ and thus $ζ_c = -18$ suggesting to use $ζ = -17$ to be slightly above the transition. The first row presents the structure factor, in good agreement with the second row representing the weight $W_2(q)$ of the dispersive band whose minimum is given by the red lines. Each column corresponds to a different planar cut orthogonal to the [111] axis in reciprocal space: $q_1 = q \cdot (1, -1, 0)/\sqrt{2}$, $q_2 = q \cdot (1, 1, -2)/\sqrt{6}$, $q_{111} = q \cdot (1, 1, 1)$. Calculations of the structure factors are performed using SCGA (see Appendix B) with $α = 1$ and inverse temperature $β = 300$ and with the $q_i$ in $a^{-1}$ units.</figcaption></figure>

<figure><img src="image_8.png" /><figcaption>Figure 7. Cluster structure of the generalized octochore model. The red dots indicate site counted with a coefficient 1 in the constrainer definition, while green and magenta dots are encapsulated with coefficients $γ_1$ and $γ_2$ (see Eq. (69)).</figcaption></figure>

minima along the $Q$ manifold depicted as red contour in Fig. 5(b), those fourfold singularities lie at finite energy and are thus absent from the low-temperature structure factor. It is only the intersection between the intensity of $W_2(q)$ and the $Q$ manifold that appears at low energy, and is thus responsible for a half-moon pattern with four-fold symmetry. Note that we also have regular twofold pinch points in this system, giving rise to regular half moons.

Tuning to $γ_1 = 1$ and $γ_2 = (1 - 2γ_1)/3 = -1/3$ produces the expected sixfold pinch points in $W_2(q)$ [Fig. 5(d)], which, upon intersecting with $Q$, generate sixfold half-moons in the structure factor [Fig. 5(c)].

### D. Example: from pinch line to half-moon surfaces on the generalized octochore model

But not all higher-rank spin liquids take the form of multifold pinch points. In three dimensions, the singularity may extend into a one-dimensional line in reciprocal space, forming a so-called pinch line [22, 24, 33, 70, 71]. How does this string singularity evolve when the interaction between clusters is strong enough to transform the flat-band gauge field fluctuations into excitations ?

Here we shall consider the generalized octochore model [20], composed of corner-sharing octahedra, where the cluster is extended to encapsulate the sites of the six neighboring octahedra (see Fig. 7(a)). Its constrainer is defined as

$$ \mathbf{c}_{\otimes} = \sum_{i \in \otimes} \mathbf{S}_i + \gamma_1 \sum_{i \in \langle \otimes \rangle} \mathbf{S}_i + \gamma_2 \sum_{i \in \langle \langle \otimes \rangle \rangle} \mathbf{S}_i. \quad (69) $$

as illustrated in Fig. 7. This system supports a pinch line for $γ_1 = -1/2$, $γ_2 = 1$ [33]. This parent cluster system can be turned into a generalized interacting cluster model following the procedure of Sec. IV B. Using parameters $γ_1 = -1/2$, $γ_2 = 1$ imposes $ζ_c = -Ξ/ω = -18$. Choosing $ζ_c = -17$ should therefore allow to observe the half
<figure><img src="image_9.png" /><figcaption>Figure 8. **Generalized checkerboard model:** (a) The interacting cluster Hamiltonian is set with parameters $\gamma_1 = -1/2$, $\gamma_2 = 0$, for varying values of $\zeta$ while fixing $\eta = 1$. The right panel represents the weight $W_2(\mathbf{q})$ of the dispersive band with contours corresponding to the energy minimum of this dispersive band for $\zeta = 0$ (blue), $\zeta = \zeta_s \approx 3.48$ (red) and $\zeta = 8$ (green). The respective structure factors are plotted on the three left panels using SCGA (see Appendix B) for inverse temperature $\beta = 20$ and with the $q_i$ in $a^{-1}$ units. (b) 3D representation of $\lambda(\mathbf{q})$, with consecutive intersections with planes located at height $\lambda_c(\zeta)$. For $\zeta < \zeta_s$, the energy minima form closed line surrounding the position of pinch points located at $\Gamma$ and $M$ special points. When $\zeta$ reaches the critical value $\zeta_s$ that is such that $\lambda_c$ reaches the saddle energy, these closed lines enter in contact to form, for $\zeta > \zeta_s$, another set of lines that are centered around $\lambda$ maxima.</figcaption></figure>

moon pattern associated with the presence of a pinch line in the dispersive band weight $W_2(\mathbf{q})$. The pinch lines on the generalized octochlore model are along the [111] axes (and equivalent ones). It means that in order to visualize the lines, one needs to look at multiple cuts orthogonal to a given [111] axis in reciprocal space. This is visible in the second row of Fig. 6 which presents the weight $W_2(\mathbf{q})$ of the dispersive band. Keeping in mind that these figures are approximately the “photographic negative” of the flat-band weight $W_1(\mathbf{q})$, the persistence of pinch points along [111] is a signature of the (negative of the) pinch line. In particular, the plane at $q_{111} = \pi/3$ is a high-symmetry plane where multiple pinch lines cross, giving rise to apparent multifold pinch points [69, 71]. Since the second line of Fig. 6 presents the weight of the dispersive band, we can delimit the $\mathcal{Q}$ manifold of lowest energy in reciprocal space (see red contours), which corresponds to the region of high-intensity in the structure factor displayed in the first row.

When the singularity in $W_2(\mathbf{q})$ is well isolated, one sees clear half moons in the structure factor. Keeping in mind that the half moon propagates along the [111] direction, the resulting pattern is actually an elongated surface in the 3D reciprocal space, similar to a long cylinder cut in two along its length. In order to emphasize the link between these patterns and the original half moons, we shall call them “half-moon surfaces”. However, since pinch lines can cross at specific $\mathbf{q}$ points, it means that some of the singularities are actually close to each other. At the high-symmetry point $q_{111} = \pi/3$ where they cross, we get multifold half moons, that deform continuously as one moves along the [111] axis. Looking at those patterns at $q_{111} = 5\pi/12$ and $\pi/2$, one could be forgiven not to recognize any half-moon pattern; those patterns are nonetheless the signature of the proximity, in parameter space, of extended one-dimensional singularities, and thus of the proximity of a particular form of higher-rank

spin liquid.

### E. Topological Lifshitz transition at the edge of half-moon phases

For interacting-cluster systems, the post-threshold behavior for $\zeta > \zeta_c$ depends crucially on whether the parent system is gapless or gapped. We therefore discuss these two cases separately.

#### 1. Gapless systems

For gapless parent models with pinch points, the system first transitions from a Coulomb spin liquid to a half-moon phase when $\zeta$ crosses $\zeta_c$ (see section III B). Then, a second, topological, transition may follow. On kagome and pyrochlore, it had been noticed that a change of the dispersive band topology could transform half-moon contours into so-called “star” patterns in the equal-time structure factor [10, 54, 72, 73]. Now we can rationalize this change in the context of the present generic theory, based on the physics of the parent Hamiltonian.

For models with a single cluster type (e.g. checkerboard, octochlore), one has $\lambda(\mathbf{q}) = |\mathbf{L}(\mathbf{q})|^2$. If the parent system dispersive band $\lambda(\mathbf{q})$ possesses a saddle point, then, as the threshold $\lambda_c(\zeta)$ rises past that saddle energy $E_s$ (see Eqs. (52,67)), the topology of the ground-state manifold $\mathcal{Q}$ changes, as illustrated in Figs. 3(b,c) and 8(b) for the generalized checkerboard lattice. The key point to understand is that the manifold $\mathcal{Q}$ varies continuously with $\zeta$. For small $\zeta$, the manifold $\mathcal{Q}$ encircles the pinch points as expected for half moons; see the blue contours on the weight $W_2(\mathbf{q})$ of the dispersive band (right panel of Fig. 8(a)) and the half moons in the structure
factor (left panel).

As $\zeta$ increases, the contour of $Q$ gets bigger and bigger, until eventually neighboring contours touch each other at a critical value $\zeta_s$

$$ \lambda_c(\zeta_s) = \frac{\Xi}{2} \left(1 - \frac{\zeta_s}{\zeta_c}\right) = E_s \Leftrightarrow \zeta_s = \zeta_c \left(1 - \frac{2E_s}{\Xi}\right) \quad (70) $$

Since we find that $E_s > \Xi/2$ in our model, we have $\zeta_s > 0$; more precisely, $\zeta_s \approx 3.48$ [74]. Contours touching at the saddle point are depicted as red lines in the right panel of Fig. 8(a) and will give rise to stronger soft-mode excitations at the saddle-point wavevectors. If not protected by the Hohenberg-Mermin-Wagner theorem, one can expect order by disorder at those wavevectors for $\zeta = \zeta_s$.

For $\zeta > \zeta_s$, the contours disconnect from each other, surrounding a different region in reciprocal space, away from the pinch-point singularity; see the green contours in the right panel of Fig. 8(a) and the corresponding structure factor for $\zeta = 8$. One gets a distinct spiral spin liquid without half moons, but with closed rings. Since the manifold $Q$ does not encircle a pinch point anymore, it does not intersect with its branches, and the weight $W_2(\mathbf{q})$ of the dispersive band has thus no reason to cancel along $Q$.

Making an analogy between the contour formed by manifold $Q$ and the Fermi surface in the electronic structure of a material, we have here the equivalent of a topological Lifshitz transition [58, 75]. At the level of the band structure, there is no symmetry breaking between the two phases, but simply a change of the topology of the contour. This Lifshitz transition thus delimits two distinct regions of the phase diagram that can be distinguished from a clear change in the structure factor. If the parent $\lambda(\mathbf{q})$ hosts multiple minima and saddle points, multiple Lifshitz transitions can occur as $\zeta$ is increased.

## 2. Gapped systems

While pinch points and half moons are coming from gapless band spectrum, this is not necessary for Lifshitz transitions. What is necessary though, is that the minimum energy manifold $Q$ forms a hypersurface in reciprocal space; this way, this hypersurface can change its topology (closing around distinct $\mathbf{q}$ points) across the phase diagram. As discussed in section IV B, not all interacting-cluster models respect that condition. But many do, such as the family of generalized interacting-cluster systems (see section IV B) which, when restricted to nearest-neighbor cluster interaction, includes the uniform cluster models of section III B.

Let us consider a parent cluster system with a finite energy gap $\Delta > 0$ at $\mathbf{q}_\Delta$; hence, $\lambda_{ld}(\mathbf{q}_\Delta) = \Delta$. The gap requires a stronger interaction $\eta$, i.e. a higher value of $\zeta = \zeta_c(\Delta) > \zeta_c$, for the dispersive band to make contact

with the flat band (see Eq. (41)). But the main difference with gapless parent systems is that, by definition, there is no fixed point $\mathbf{q}^*$ where $\lambda_{ld}(\mathbf{q}^*) = 0$ would have imposed $\Lambda_{ld}(\mathbf{q}^*) = 0$ (see Eqs. (37) and (65)). As $\zeta$ increases above $\zeta_c(\Delta)$, the dispersive band $\Lambda_{ld}(\mathbf{q})$ simply goes through the flat band and preserves its minimum energy at $\mathbf{q}_\Delta$ since

$$ \nabla_q \lambda_{ld}(\mathbf{q}_\Delta) = 0 \Rightarrow \nabla_q \Lambda_{ld}(\mathbf{q}_\Delta) = 0. \quad (71) $$

The ground-state “manifold” is therefore a finite set of points in reciprocal space; we expect the system to long-range order at wave-vector $\mathbf{q}_\Delta$.

But Eq. (51) is still valid and the manifold $Q$ – as defined by Eqs. (52) and (67) – remains an extremum of the dispersive band $\Lambda_{ld}(\mathbf{q})$. The gap simply shifted its domain of existence to larger values of $\zeta$

$$ \begin{align} \lambda_c &= \frac{\Xi}{2} \left(1 - \frac{\zeta}{\zeta_c}\right) \ge \Delta \nonumber \\ \Leftrightarrow \zeta &\ge \zeta_\Delta = \zeta_c \left(1 - \frac{2\Delta}{\Xi}\right) = \zeta_c(\Delta) - \frac{\Delta}{\Xi} \tag{72} \\ \Leftrightarrow \zeta &\ge \zeta_\Delta > \zeta_c(\Delta) \nonumber \end{align} $$

When $\zeta = \zeta_\Delta$, we have $\lambda_c = \Delta$ and the manifold $Q$ includes all $\mathbf{q}_\Delta$ points. By analytic continuation, as $\zeta$ increases, manifold $Q$ will form a closed hypersurface around $\mathbf{q}_\Delta$ points. This is the same reasoning as before for the apparition of half moons, except that the $\mathbf{q}_\Delta$ points are not attached to pinch points, and the resulting patterns in equal-time structure factor will be model dependent.

Remarkably, as $\zeta$ increases from $-\infty$ (parent cluster model) to positive values, the ground state of the model starts as a fragile topological spin liquid (using the terminology of Ref. [22]), goes through a nonsingular band touching at $\zeta_c(\Delta)$, then orders magnetically at wavevector $\mathbf{q}_\Delta$ up to $\zeta_\Delta$ where it becomes a spiral spin liquid. And if the lowest dispersive band of the parent system $\lambda_{ld}$ has one (or more) saddle point(s), a Lifshitz transition(s) shall occur when $\lambda_c(\zeta)$ exceeds a saddle energy $E_s$, separating distinct spiral spin liquids.

As an illustration, let us consider the hexagonal kagome model from Table I, known to be gapped [22, 23]. Choosing the simplest constrainer

$$ \mathcal{C}_\odot = \sum_{i \in \mathcal{O}} \mathbf{S}_i, \quad (73) $$

places the model in the class satisfying conditions (i)–(iii), with $\Xi = \xi = 6$, and, since the lattice is corner-sharing, $\Omega = \omega = 1$. For this system $\lambda_c = 3(1 + \zeta/6)$ and the three relevant thresholds for $\zeta$ are thus

$$ \zeta_c = -6, \quad \zeta_\Delta = 0, \quad \zeta_s = 2, \quad (74) $$

as $\Delta = 3$ and $E_s = 4$. The three successive regimes are clearly visible in the structure factors of Fig. 9. For $\zeta < \zeta_c$
<figure><img src="image_10.png" /></figure>

**Figure 9.** **Hexagonal checkerboard model:** The top line contains structure factors computed for different values of ζ using the SCGA (see Appendix B) with η = 1 and inverse temperature β = 100. The bottom line shows the dispersive band weight W₃(q) together with the ground state manifold Q depicted as green dots or lines. In both cases the plots are obtained with qₓ, qᵧ ∈ [-2π, 2π] in inverse units of a cluster edge. The structure factors show four different phases as increasing ζ. For ζ < ζ_c the structure factor is similar to the one of the non interacting system, presenting no pinch point as the band structure is gapped. When ζ becomes larger than ζ_c the flat bands cease to be the ground state, which is now constituted of the points qΔ for which λ(qΔ) = Δ. The structure factor thus depicts a selection of the dispersive band weight around these minima, producing bright regions increasing when ζ approaches ζΔ. At this point λ_c(ζ) reaches the value Δ, and the ground state manifold becomes an hypersurface, producing bright lines in the structure factor. When ζ finally reaches ζ_s the structure of the ground state manifold changes, corresponding to a Lifshitz transition, that appears in the structure factor as the bright lines changing of pattern.

the structure factor closely matches that of the non-
interacting parent model (η = 0), showing only smooth
features as the ground state is the flat-band manifold
that contains no singularities. For ζ_c < ζ < ζΔ, bright
spots appear and broaden as ζ → ζΔ, reflecting that
the ground state now consists only of the points qΔ lo-
cated at the *K* high-symmetry points in this model. At
ζ = ζΔ, the ground-state set becomes a family of tri-
angular contours surrounding the *K* points, producing
triangular bright lines in the structure factor. These con-
tours grow with ζ until they touch at the *M* points when
ζ = ζ_s, after which they reconnect into circular contours
encircling the maxima of λ at the Γ points for ζ > ζ_s.
These circular contours are expected to shrink into bright
regions, as when ζ becomes large enough so that λ_c(ζ)
equates the maxima of the dispersive band λ, the ground
state manifold *Q* becomes the set of points q_m such that
λ(q_m) = Max[λ(q)].

V. DISCUSSION

Somewhat counter-intuitively, the tools recently devel-
oped for the classification of classical spin liquids [20–
24] prove useful to explore the physics of systems be-
yond spin liquids. In order to destabilize the spin-liquid
ground state, we used methods of flat-band engineering
coming from graph theory, in particular connectivity ma-
trices that are well suited to account for interaction be-
tween spin clusters [10, 13, 25, 26]. The two frameworks,

connectivity matrix and constrainers of spin liquids, are
actually two faces of the same coin (see Eq. (13) and
Ref. [24]).

Here we have built a generic theory for multifold half
moons, encompassing the evolution of pinch lines into ex-
tended surfaces as signatures of spiral spin liquids whose
magnetic texture hosts higher-rank gauge charges in the
ground state. The manifold of minimum energy, forming
a closed contour in reciprocal space, behaves as an ef-
fective Fermi surface of the parent noninteracting system
and undergoes topological changes at Lifshitz transitions,
when the dispersive band of the parent cluster Hamilto-
nian possesses saddle points. In most, if not all, three-
dimensional models, thermal order by disorder will take
place at very low temperature and long-range order the
system. As was the case on pyrochlore [10], this mag-
netic ordering can be seen as the crystallization of the
gauge charges that populate the dispersive band. In the
case of higher rank spin liquids, one might thus be able
to observe a crystallization of fractons. When the parent
cluster Hamiltonian is gapped, half moons disappear but
spiral spin liquids and Lifshitz transitions remain.

The present work is primarily aimed at establishing a
generic and unifying theoretical framework for interact-
ing cluster systems, rather than proposing specific micro-
scopic Hamiltonians tailored for immediate experimental
realization. In particular, the cluster-cluster interactions
required for the emergence of half-moon patterns may
effectively generate couplings between sites at relatively
long distances, which render the models discussed here
challenging to realize in conventional magnetic materials. Nevertheless, cluster-based descriptions are a well-established and physically motivated route for stabilizing classical spin liquids, and the mechanisms identified in this work are expected to apply broadly to any system whose low-energy physics admits an effective cluster formulation.

The present framework offers an infinite variety of theoretical models supporting a dense phase of gapless gauge charges at intermediate temperature, stabilized by cluster interaction. Their non-local dynamics is an open question, and their connection with local momentum vortices found in spiral spin liquids [62] – equivalent to quadrupoles of fractons – is an appealing direction to investigate. The investigation of half moons in quantum spins $S = 1/2$ has given interesting physics on kagome [54, 72, 73] and can now be applied to higher-rank and pinch-line spin liquids.

*Acknowledgements:* The authors thank Paul McClarty and Jérôme Cayssol for useful discussions. They acknowledge financial support from Grant No. ANR-23-CE30-0038-01.

## Appendix A: Band structure and structure factor relation

Both the Hamiltonian matrix $H^{(1)} = h^{v \leftarrow c} [h^{v \leftarrow c}]^\dagger$ and the matrix $M = [h^{v \leftarrow c}]^\dagger h^{v \leftarrow c}$ ruling the projector into the flat bands manifold are built from the same connectivity matrix $h^{v \leftarrow c}$. This implies a deep connection between the band structure and the structure factor that we detail in this Appendix.

### 1. Equivalence between contact point and pinch point

Here we show that for cluster Hamiltonians there is an equivalence between the presence of a band touching between dispersive and flat bands, and the presence of pinch point in the static structure factor. The rectangular matrix $h^{v \leftarrow c}$ admits a singular value decomposition

$$h^{v \leftarrow c} = U \Sigma V^{\dagger}. \qquad (A1)$$

where $U$ is an $n_s \times n_s$ unitary matrix, $V$ is an $n_c \times n_c$ unitary matrix and $\Sigma$ is an $n_s \times n_c$ diagonal matrix of singular values $\sigma_\nu$,

$$\Sigma = \begin{bmatrix} \sigma_1 & 0 & \dots & 0 \\ 0 & \sigma_2 & \dots & 0 \\ \vdots & \vdots & \ddots & \vdots \\ 0 & 0 & \dots & \sigma_{n_c} \\ 0 & 0 & \dots & 0 \\ \vdots & \vdots & \ddots & \vdots \\ 0 & 0 & \dots & 0 \end{bmatrix}. \qquad (A2)$$

This expression allows to get that

$$H^{(1)} = U \Sigma \Sigma^\dagger U^\dagger, \quad M = V \Sigma^\dagger \Sigma V \qquad (A3)$$

where $\Sigma \Sigma^\dagger$ is the diagonalized Hamiltonian, composed of $n_c$ trivial eigenvalues and $n_s - n_c$ dispersive eigenvalues $\lambda_i \equiv (\sigma_i)^2$, and $\Sigma^\dagger \Sigma$ is the diagonal matrix with diagonal elements $\lambda_i$. This imposes that the determinant of the matrix $M$, that is equal to the product of its eigenvalues, is equal to the product

$$\det M(\mathbf{q}) = \prod_{i=1}^{n_c} \lambda_i(\mathbf{q}) \qquad (A4)$$

of the dispersive bands dispersions $\lambda_i(\mathbf{q})$. Hence, at a contact point $\mathbf{q}^*$ where the lowest dispersive band vanishes, $\det M(\mathbf{q}^*) = 0$ and the flat-sector projector

$$\Pi = I - h^{v \leftarrow c} M^{-1} [h^{v \leftarrow c}]^\dagger \qquad (A5)$$

becomes singular. This nonanalyticity appears in the equal-time structure factor as a pinch-type singularity. Thus in cluster systems any contact point in the structure factor induces a pinch point in the static structure factor.

### 2. Contact points in cluster models are necessarily singular

The matrix $M(\mathbf{q}) = h^\dagger h$ is the Gram matrix of the constraint vectors $\{\mathbf{L}_X(\mathbf{q})\}_{X=1}^{n_c}$ (columns of $h$):

$$M_{XY}(\mathbf{q}) = \mathbf{L}_X(\mathbf{q})^* \cdot \mathbf{L}_Y(\mathbf{q}). \qquad (A6)$$

A basic property of Gram matrices is that $\det M(\mathbf{q}) = 0$ if and only if the family $\{\mathbf{L}_X(\mathbf{q})\}$ loses linear independence, i.e. the span of dispersive bands eigenvectors drops rank. Consequently, any gap closing (contact point) requires a rank drop, which can occur only if either

(i) a constraint vector vanishes: $\mathbf{L}_X(\mathbf{q}^*) = 0$, or

(ii) a nontrivial linear combination vanishes: there exists a linear combination $\mathbf{L}_c(\mathbf{q})$ of the constraint vectors $\mathbf{L}_X(\mathbf{q})$ such that $\mathbf{L}_c(\mathbf{q}^*) = 0$.

Both mechanisms are singular in the sense that some dispersive bands eigenvectors (or linear combination thereof) vanishes at $\mathbf{q}^*$. Therefore, in the cluster class, every contact point is singular.

The same mechanism can produce singular gap closings along higher-dimensional manifolds in momentum space (pinch lines/planes), corresponding to rank loss on sets of codimension 1 or higher. Note that in these cases this does not correspond to the emergence of a pinch point in the structure factor if the dimension of the contact manifold is the system dimension minus one[23], but that it does always imply the existence of non localized flat band states[22, 23].
### 3. Order of the Gauss law and contact point dispersion

Assume first that a single constrainer vanishes at $\mathbf{q}^*$: $\mathbf{L}_c(\mathbf{q}^*) = 0$. Let $V(\mathbf{q})$ be the span of the remaining constraint vectors at $\mathbf{q}$, and $P_V$ the orthogonal projector onto $V(\mathbf{q})$. The Gram determinant factorizes as

$$ \det M(\mathbf{q}) = \|\mathbf{L}_c(\mathbf{q}) - P_V \mathbf{L}_c(\mathbf{q})\|^2 \det M_V(\mathbf{q}), \quad (A7) $$

where $M_V$ is the Gram matrix on $V(\mathbf{q})$. If the first non-vanishing term in the multivariate Taylor expansion of $\mathbf{L}_c(\mathbf{q})$ about $\mathbf{q}^*$ has total degree $n$,

$$ \mathbf{L}_c(\mathbf{q}) = \mathbf{\Phi}(\delta\mathbf{q})^{(n)} + \mathcal{O}(\|\delta\mathbf{q}\|^{n+1}), \quad \delta\mathbf{q} = \mathbf{q} - \mathbf{q}^*, $$

then $\|\mathbf{L}_c(\mathbf{q}) - P_V \mathbf{L}_c(\mathbf{q})\|^2 = \mathcal{O}(\|\delta\mathbf{q}\|^{2n})$ while $\det M_V(\mathbf{q}) = \mathcal{O}(1)$. Thus

$$ \det M(\mathbf{q}) = \mathcal{O}(\|\delta\mathbf{q}\|^{2n}). \quad (A8) $$

More directly, because $M(\mathbf{q})$ is positive semidefinite, the Courant-Fischer theorem gives

$$ \lambda_1(\mathbf{q}) = \min_{\|\mathbf{x}\|=1} \mathbf{x}^\dagger M(\mathbf{q}) \mathbf{x}. \quad (A9) $$

At $\mathbf{q}^*$ one has $\alpha^\dagger M(\mathbf{q}^*)\alpha = 0$ for $\alpha = e_c$ (the unit vector selecting $\mathbf{L}_c$), so near $\mathbf{q}^*$ the minimizer remains $\alpha$ to leading order, and

$$ \lambda_1(\mathbf{q}) = \alpha^\dagger M(\mathbf{q})\alpha = \|\mathbf{L}_c(\mathbf{q})\|^2 = \mathcal{O}(\|\delta\mathbf{q}\|^{2n}). \quad (A10) $$

Hence the lowest dispersive band has local dispersion of order $2n$, which matches a Gauss-law of differential order $n$ (Maxwell for $n=1$, higher-rank for $n > 1$).

Now consider the linear-dependence mechanism (ii). There exists a unit vector $\alpha \in \mathbb{C}^{n_c}$ such that

$$ \mathbf{L}_c(\mathbf{q}) \equiv h^{v \leftarrow c}(\mathbf{q})\alpha = \mathbf{\Phi}(\delta\mathbf{q})^{(n)} + \mathcal{O}(\|\delta\mathbf{q}\|^{n+1}). \quad (A11) $$

Since $M = h^\dagger h$,

$$ \lambda_1(\mathbf{q}) \le \alpha^\dagger M(\mathbf{q})\alpha = \|\mathbf{L}_c(\mathbf{q})\|^2 = \mathcal{O}(\|\delta\mathbf{q}\|^{2n}), \quad (A12) $$

and positivity implies this is the actual leading scaling. Thus, again, the lowest dispersive band softens with order $2n$, and the emergent Gauss law is of order $n$.

**Degenerate contact point.** If the rank drops by $r > 1$ at $\mathbf{q}^*$, then $r$ independent combinations $\mathbf{L}_c^{(a)}$ vanish simultaneously, and the $r$ lowest eigenvalues soften. The product of these $r$ eigenvalues scales as $\mathcal{O}(\|\delta\mathbf{q}\|^{2(n_1+\cdots+n_r)})$, with each $n_a$ the lowest total degree for the corresponding critical combination. The single band touching case above is recovered for $r=1$.

**Visibility in $S(\mathbf{q})$.** While form-factor weights can redistribute intensity and obscure a pinch point at special symmetry points, the singular character (rank loss) of the projector is unaffected. Hence the contact point $\Leftrightarrow$ pinch-type singularity equivalence remains valid.

## Appendix B: Self Consistent Gaussian Approximation

Consider any isotropic Heisenberg Hamiltonian that can be expressed as

$$ \begin{aligned} \mathcal{H} &= \frac{1}{2} \sum_{i,j} H_{i,j} \mathbf{S}_i \cdot \mathbf{S}_j \\ &= \frac{1}{2N} \sum_{\mathbf{q}} \sum_{\mu,\nu} H_{\mu\nu}(\mathbf{q}) \mathbf{S}_{\mu}(\mathbf{q}) \cdot \mathbf{S}_{\nu}(-\mathbf{q}) \end{aligned} \quad (B1) $$

with $N$ the number of unit cells. Further consider there are $n_s$ sublattices, and denote the $n_s$ eigenvalues and eigenvectors from the Hamiltonian $H(\mathbf{q})$ as

$$ \varepsilon_{\kappa}(\mathbf{q}), \quad \psi_{\kappa}(\mathbf{q}). \quad (B2) $$

The self consistent Gaussian approximation[55] corresponds to take into account the spin length constraint $|\mathbf{S}_i| = 1$ only in average, enforcing only the constraint $\langle|\mathbf{S}_i|\rangle = 1$. Using that the three spin components are equivalent, the constraint that is really enforced in practice is

$$ \langle (S_i^\alpha)^2 \rangle = \frac{1}{3}. \quad (B3) $$

This can be done using a single Lagrange multiplier $\lambda$ to enforce the constraint in the Hamiltonian, while enforcing the hard constraint for every spin would amount to introduce $Nn_s$ Lagrange multipliers. Also shifting the energy origin by defining $\varepsilon_0 = \text{Min}(\varepsilon_\mu)$ and considering the positive semi definite matrix $K(\mathbf{q}) = H(\mathbf{q}) - \varepsilon_0 I$, the Hamiltonian can finally be expressed as

$$ \mathcal{H} = \frac{1}{2N} \sum_{\mathbf{q}} \left( K_{\mu\nu}(\mathbf{q}) + \frac{\lambda}{\beta} \delta_{\mu\nu} \right) \mathbf{S}_{\mu}(\mathbf{q}) \cdot \mathbf{S}_{\nu}(\mathbf{q}). \quad (B4) $$

The eigenvalues of $K$ are equal to the ones of $H$ but shifted of $\varepsilon_0$ and are therefore positive, while its eigenvectors are simply identical. For a semi-positive definite matrix $A$ the general formula

$$ \frac{\int \prod_i dx_i x_\mu x_\nu e^{-\frac{1}{2} \mathbf{x}^\dagger A \mathbf{x}}}{\int \prod_i dx_i e^{-\frac{1}{2} \mathbf{x}^\dagger A \mathbf{x}}} = (A^{-1})_{\mu\nu} \quad (B5) $$

can be used to write that, defining $M = \beta K + \lambda I$,

$$ \begin{aligned} \langle S_{\mu}^{\alpha}(-\mathbf{q}) S_{\nu}^{\alpha}(\mathbf{q}) \rangle &= \frac{\int \prod_i dS_i^{\alpha} S_{\mu}^{\alpha} S_{\nu}^{\alpha} e^{-\frac{1}{2N} \sum_{\mathbf{q}'} \sum_{\eta, \xi} M_{\eta\xi} S_{\eta}^{\alpha} S_{\xi}^{\alpha}}}{\int \prod_i dS_i^{\alpha} e^{-\frac{1}{2N} \sum_{\mathbf{q}'} \sum_{\eta, \xi} M_{\eta\xi} S_{\eta}^{\alpha} S_{\xi}^{\alpha}}} \\ &= N (M^{-1})_{\mu\nu} = N [\lambda I + \beta K(\mathbf{q})]_{\mu\nu}^{-1}. \end{aligned} $$
This allows to self consistently fix the Lagrange multiplier $\lambda$ by enforcing the constraint

$$ \begin{aligned} \langle (S_i^\alpha)^2 \rangle &= \langle S_i^\alpha S_i^\alpha \rangle = \frac{1}{N} \sum_i \langle S_i^\alpha S_i^\alpha \rangle \\ &= \frac{1}{N^2} \sum_{\mathbf{q}} \sum_{\mu=1}^{n_s} \langle S_\mu^\alpha (-\mathbf{q}) S_\mu^\alpha (\mathbf{q}) \rangle \quad (B6) \\ &= \frac{1}{N} \sum_{\mathbf{q}} \sum_{\mu=1}^{n_s} [\lambda I + \beta K(\mathbf{q})]_{\mu\mu}^{-1} = \frac{1}{3}, \end{aligned} $$

where we used that all spins in the system are equivalent by symmetry to write the first line. This also allows to express the static structure factor as

$$ \begin{aligned} S(\mathbf{q}) &= \frac{1}{3N} \sum_{\mu, \nu} \langle \mathbf{S}_{\mu}(-\mathbf{q}) \cdot \mathbf{S}_{\nu}(\mathbf{q}) \rangle = \frac{1}{N} \sum_{\mu, \nu} \langle S_{\mu}^{\alpha}(-\mathbf{q}) \cdot S_{\nu}^{\alpha}(\mathbf{q}) \rangle \\ &= \sum_{\mu, \nu} [\lambda I + \beta K(\mathbf{q})]_{\mu\nu}^{-1} = \sum_{\kappa=1}^{n_s} \sum_{\mu, \nu} \frac{[\psi_{\kappa}^{*}(\mathbf{q})]_{\mu} [\psi_{\kappa}(\mathbf{q})]_{\nu}}{\lambda + \beta \varepsilon_{\kappa}(\mathbf{q})} \end{aligned} $$

because the three spin components are equivalent, and then the spin components can be treated as three independent variables.

In the limit of zero temperature $\beta \to \infty$, only flat bands contributions are selected. In this case the static structure factor becomes simply proportional to the sum of the components of the projector into flat band manifold

$$ \begin{aligned} \Pi_{\mu\nu} &= \sum_{\kappa=1}^{n_{\text{f,b}}} [\psi_{\kappa}^{*}(\mathbf{q})]_{\mu} [\psi_{\kappa}(\mathbf{q})]_{\nu} \\ &= I_{\mu\nu} - \sum_{\kappa=n_{\text{f,b}}+1}^{n_s} [\psi_{\kappa}^{*}(\mathbf{q})]_{\mu} [\psi_{\kappa}(\mathbf{q})]_{\nu}. \end{aligned} \quad (B7) $$

As the manifold spanned by the dispersive bands eigenvectors $\psi_\kappa$, $\kappa > n_{\text{f,b}}$ is identical to the one spanned by the constraint vectors, this projector can be re-expressed as in Eq. (15).

## Appendix C: Relevance of the GLT method for clusters systems

For cluster systems, being able to express the Hamiltonian as a sum of squares

$$ \begin{aligned} \mathcal{H} &= \sum_{i,j} H_{ij} \mathbf{S}_i \cdot \mathbf{S}_j \\ &= \sum_{n,X} |\mathcal{C}_{n,X}|^2 - N_{u.c} \sum_{X} \sum_{i \in X} (\gamma_i^X)^2 \end{aligned} \quad (C1) $$

corresponds in fact to identify that the Langrange multipliers allowing to find a real space physical ground state are $\gamma_i^X$, as this allows to express the Hamiltonian as

a positive definite term which is minimized by configurations satisfying the ground state constraint $\mathcal{C}_{n,X} = 0 \forall n, X$. Fourier transforming this Hamiltonian leads to

$$ \begin{aligned} \mathcal{H} = & \sum_{\mathbf{q}} \sum_{X} \sum_{\mu, \nu} [h_{\mu}^{X}(\mathbf{q}) \mathbf{S}_{\mu}(\mathbf{q})] \cdot [h_{\nu}^{X}(-\mathbf{q}) \mathbf{S}_{\nu}(-\mathbf{q})] \\ & - \lambda N_{u.c} \sum_{\mu} \alpha_{\mu}^{2} \end{aligned} \quad (C2) $$

where we defined the parameters $\lambda, \alpha_\mu$ such that

$$ \lambda \sum_{\mu} \alpha_{\mu}^{2} = \sum_{X} \sum_{i \in X} (\gamma_{i}^{X})^{2} \quad (C3) $$

where $h_\mu^X(\mathbf{q})$ is the $\mu$ component of Fourier transform of the constrainer $\mathcal{C}_{0,X}$ as explained in the main text. Under this form the Hamiltonian corresponds to the Generalized Luttinger Tisza (GLT) method [39]. This procedure encodes a soft spin length constraint

$$ \sum_{n, \mu} \alpha_{\mu}^{2} |\mathbf{S}_{n, \mu}|^{2} = \sum_{\mathbf{q}} \sum_{\mu} \alpha_{\mu}^{2} |\mathbf{S}_{\mu}(\mathbf{q})|^{2} = N_{u.c} \sum_{\mu} \alpha_{\mu}^{2}, \quad (C4) $$

thanks to a unique Lagrange multiplier $\lambda$, associated with a set of parameters $\alpha_\mu$, specific to the different sublattices, that can be cleverly chosen by analyzing the symmetries of the problem. For a generic system, encoding this soft constraint would lead to a Hamiltonian

$$ \begin{aligned} \mathcal{H} = & \sum_{\mathbf{q}} \sum_{\mu, \nu} (H_{\mu\nu}(\mathbf{q}) - \alpha_{\mu}\alpha_{\nu}\delta_{\mu\nu}) \mathbf{S}_{\mu}(\mathbf{q}) \cdot \mathbf{S}_{\nu}(-\mathbf{q}) \\ & + \lambda N_{u.c} \sum_{\mu} \alpha_{\mu}^{2}. \end{aligned} \quad (C5) $$

Therefore, solving the ground state would amount to diagonalize the matrix

$$ F_{\mu\nu} = \frac{1}{\alpha_{\mu}\alpha_{\nu}} H_{\mu\nu}(\mathbf{q}) \quad (C6) $$

as the configuration associated with its lowest eigenvalue will naturally minimize the Hamiltonian (C5). Finally, one would need to check the resulting configuration fulfills the strong condition $|\mathbf{S}_i| = 1$ to ensure the built state is indeed a physical state.

For cluster systems, however, the story is different. Indeed, because the real-space structure of the Hamiltonian imposes a macroscopic number of ground-state configurations, one should not expect a unique configuration associated with the lowest eigenvalue of $F$, but rather a macroscopic number, since the lowest eigenvalues of $F$ form a flat band in reciprocal space. This macroscopic degeneracy originates from local constraints encoded in the cluster structure. Consequently, for these systems, the appropriate choice of parameters $\alpha_\mu$ is given by Eq. (C3) and corresponds to selecting parameters such that the matrix $F$ admits flat bands at the bottom of its spectrum. In this setting, checking that a single obtained configuration satisfies the strong constraint
$|S_i| = 1$ should be replaced by demonstrating that there exists a macroscopic family of superpositions of flat-band states that fulfill this constraint. Since flat bands are generically associated with the existence of compact localized states [14, 22], these can typically be superposed to form physical ground-state configurations, even in systems with inequivalent sublattices. This scenario has been confirmed, for example, by Monte Carlo simulations for the centered pyrochlore lattice [36].

The GLT method aims to determine the ground state of spin systems, but it does not address the excited energy bands. When discussing the upper part of the spectrum of the Hamiltonian matrix, one should therefore refer to the large-$\mathcal{N}$ approximation, in which the strong constraint $|S_i| = 1$ is relaxed and the spin degrees of freedom are treated as scalar fields. While, for the reasons discussed above in the context of the GLT, this approximation provides excellent results for studies focusing on the properties of the ground-state manifold, there is no a priori reason to consider it reliable for higher, dispersive energy bands. However, in the particular case of spiral spin liquids discussed in the present manuscript, the ground-state manifold obtained within the large-$\mathcal{N}$ approximation also displays a large degeneracy, suggesting again the existence of superpositions of these states that realize physical ground states. This explains the excellent agreement between Monte Carlo simulations and large-$\mathcal{N}$ predictions in previous observations of half-moon patterns for interacting cluster versions of the kagome and pyrochlore lattices [10]. Since the ground-state manifold realizes a cut of the non-interacting band-system band structure, such an agreement indicates that the band structure obtained within the large-$\mathcal{N}$ approximation is reliable, at least for systems with equivalent sublattices. A complementary way to understand the success of the large-$\mathcal{N}$ approximation for cluster spin liquid systems is that the extensive entropy associated with the massively degenerate ground-state manifold allows for strong thermal fluctuations even at low temperature, thereby washing out small discrepancies between the true ground state and the states predicted within the approximation when considering observables such as the structure factor.

## Appendix D: General case of non equivalent interacting clusters

The main text derivation of Sec. III B can in fact be generalized to the case of cluster systems composed of

inequivalent types of clusters and mixing different type of geometrical intersections between clusters, as for example the two versions of the square kagome lattice depicted in Table. I. This requires to generalize the Hamiltonian (24) to allow for more complex cluster-cluster connectivity, setting it as

$$ \mathcal{H} = \alpha \sum_n |\mathcal{C}_n|^2 + 2\eta \sum_{\langle m,n \rangle} \Gamma_{mn} \mathcal{C}_m \cdot \mathcal{C}_n. \quad (D1) $$

This Hamiltonian can also be well expressed in the framework of connectivity matrices, defining the connectivity matrix $h^{c\leftarrow c}$ as a $N_c \times N_c$ connectivity matrix linking cluster centers to centers of neighboring clusters. Its coefficients $h_{mn}^{c\leftarrow c}$ are equal to $\Gamma_{mn}$ if the clusters $m$ and $n$ are neighboring clusters, and zero otherwise. Now one can choose the cluster-cluster links weights $\Gamma_{mn}$ such that

$$ h^{c\leftarrow v} h^{v\leftarrow c} = h^{c\leftarrow c} + \Xi I_{N_c}. \quad (D2) $$

This corresponds in fact to represent the link between a pair of clusters $(m, n)$ as a sum of the different two move paths going from the cluster center-$m$ to a vertex shared by the clusters $m$ and $n$ and then to the center of the cluster $n$. With this construction the Hamiltonian matrix expresses as

$$ H = (\alpha - \eta \Xi) H^{(1)} + \eta (H^{(1)})^2 \quad (D3) $$

and the band structure of the interaction cluster model can therefore be expressed as

$$ \Lambda_{\mu}(\mathbf{q}) = (\alpha - \eta \Xi) \lambda_{\mu}(\mathbf{q}) + \eta \lambda_{\mu}(\mathbf{q})^2. \quad (D4) $$

For systems satisfying the set of conditions (i) to (iii) the present choice cluster-cluster interaction would correspond to encapsulate the coefficient $\Omega$ directly within $\eta$.

[1] C. L. Henley, Phys. Rev. B 71, 014424 (2005).

[2] C. L. Henley, Annu. Rev. Condens. Matter Phys. 1, 179 (2010).

[3] C. Castelnovo, R. Moessner, and S. L. Sondhi, Nature 451, 42 (2008).

[4] S. V. Isakov, K. Gregor, R. Moessner, and S. L. Sondhi, Phys. Rev. Lett. 93, 167204 (2004).

[5] P. W. Anderson, Phys. Rev. 102, 1008 (1956).

[6] D. A. Garanin and B. Canals, Phys. Rev. B 59, 443 (1999).
[7] J. Robert, B. Canals, V. Simonet, and R. Ballou, Phys. Rev. Lett. **101**, 117207 (2008).

[8] S. Guitteny, J. Robert, P. Bonville, J. Ollivier, C. Decorse, P. Steffens, M. Boehm, H. Mutka, I. Mirebeau, and S. Petit, Phys. Rev. Lett. **111**, 087201 (2013).

[9] H. Yan, R. Pohle, and N. Shannon, Phys. Rev. B **98**, 140402 (2018).

[10] T. Mizoguchi, L. D. Jaubert, R. Moessner, and M. Udagawa, Physical Review B **98**, 144446 (2018).

[11] S. Zhang, H. J. Changlani, K. W. Plumb, O. Tchernyshyov, and R. Moessner, Phys. Rev. Lett. **122**, 167203 (2019).

[12] H. Yan, J. Romhányi, A. Thomasen, and N. Shannon, Phys. Rev. B **110**, 195117 (2024).

[13] T. Mizoguchi and M. Udagawa, Physical Review B **99**, 235118 (2019).

[14] J.-W. Rhim and B.-J. Yang, Phys. Rev. B **99**, 045107 (2019).

[15] A. Graf and F. Piéchon, Phys. Rev. B **104**, 195128 (2021).

[16] W. Maimaiti, A. Andreanov, and S. Flach, Phys. Rev. B **103**, 165116 (2021).

[17] K. Roychowdhury, J. Attig, S. Trebst, and M. J. Lawler, Phys. Rev. Res. **6**, 043273 (2024).

[18] M. Udagawa, H. Nakai, and C. Hotta, Phys. Rev. B **110**, L201109 (2024).

[19] H. Nakai, M. Udagawa, and C. Hotta, Phys. Rev. B **111**, 235154 (2025).

[20] O. Benton and R. Moessner, Phys. Rev. Lett. **127**, 107202 (2021).

[21] H. Yan, O. Benton, R. Moessner, and A. H. Nevidomskyy, Phys. Rev. B **110**, L020402 (2024).

[22] H. Yan, O. Benton, A. H. Nevidomskyy, and R. Moessner, Phys. Rev. B **109**, 174421 (2024).

[23] N. Davier, F. A. Gómez Albarracín, H. D. Rosales, and P. Pujol, Phys. Rev. B **108**, 054408 (2023).

[24] Y. Fang, J. Cano, A. H. Nevidomskyy, and H. Yan, Phys. Rev. B **110**, 054421 (2024).

[25] H. Katsura, I. Maruyama, A. Tanaka, and H. Tasaki, Europhysics Letters **91**, 57007 (2010).

[26] K. Essafi, L. D. C. Jaubert, and M. Udagawa, Journal of Physics: Condensed Matter **29**, 315802 (2017).

[27] J. T. Chalker, P. C. W. Holdsworth, and E. F. Shender, Phys. Rev. Lett. **68**, 855 (1992).

[28] A. Codello, Journal of Physics A: Mathematical and Theoretical **43**, 385002 (2010).

[29] J. Rehn, A. Sen, and R. Moessner, Phys. Rev. Lett. **118**, 047201 (2017).

[30] R. Verresen, M. D. Lukin, and A. Vishwanath, Phys. Rev. X **11**, 031005 (2021).

[31] R. Siddharthan and A. Georges, Phys. Rev. B **65**, 014417 (2001).

[32] M. G. Gonzalez, Y. Iqbal, J. Reuther, and H. O. Jeschke, Communications Materials **6**, 96 (2025), arXiv:2410.10385 [cond-mat.str-el].

[33] N. Niggemann, Y. Iqbal, and J. Reuther, Phys. Rev. Lett. **130**, 196601 (2023).

[34] J. M. Hopkinson, S. V. Isakov, H.-Y. Kee, and Y. B. Kim, Phys. Rev. Lett. **99**, 037201 (2007).

[35] Y. Okamoto, M. Nohara, H. Aruga-Katori, and H. Takagi, Phys. Rev. Lett. **99**, 137207 (2007).

[36] R. P. Nutakki, L. D. C. Jaubert, and L. Pollet, SciPost Physics **15**, 040 (2023).

[37] J. M. Luttinger and L. Tisza, Phys. Rev. **70**, 954 (1946).

[38] J. M. Luttinger, Phys. Rev. **81**, 1015 (1951).

[39] T. Kaplan and N. Menyuk, Philosophical Magazine **87**, 3711 (2007).

[40] J. Villain, R. Bidaux, J.-P. Carton, and R. Conte, Journal de Physique **41**, 1263 (1980).

[41] J. N. Reimers, A. J. Berlinsky, and A. C. Shi, Phys. Rev. B **43**, 865 (1991).

[42] P. H. Conlon and J. T. Chalker, Phys. Rev. B **81**, 224413 (2010).

[43] R. Moessner and J. T. Chalker, Phys. Rev. B **58**, 12049 (1998).

[44] A. Prem, S. Vijay, Y.-Z. Chou, M. Pretko, and R. M. Nandkishore, Phys. Rev. B **98**, 165140 (2018).

[45] H. Yan, O. Benton, L. D. C. Jaubert, and N. Shannon, Phys. Rev. Lett. **124**, 127203 (2020).

[46] M. E. Zhitomirsky, Phys. Rev. B **78**, 094423 (2008).

[47] K. Morita and N. Shibata, Journal of the Physical Society of Japan **85**, 033705 (2016), https://doi.org/10.7566/JPSJ.85.033705.

[48] N. Caci, K. Karlová, T. Verkholyak, J. Strecka, S. Wessel, and A. Honecker, Phys. Rev. B **107**, 115143 (2023).

[49] I. Kimchi and A. Vishwanath, Phys. Rev. B **89**, 014414 (2014).

[50] J. G. Rau and M. J. P. Gingras, Nature Communications **7**, ncomms12234 (2016).

[51] M. Udagawa, L. D. C. Jaubert, C. Castelnovo, and R. Moessner, Phys. Rev. B **94**, 104416 (2016).

[52] T. Mizoguchi, L. D. C. Jaubert, and M. Udagawa, Phys. Rev. Lett. **119**, 077207 (2017).

[53] K. Tokushuku, T. Mizoguchi, and M. Udagawa, Phys. Rev. B **100**, 134415 (2019).

[54] T. Lugan, L. D. C. Jaubert, M. Udagawa, and A. Ralko, Phys. Rev. B **106**, L140404 (2022).

[55] D. A. Garanin, Phys. Rev. B **53**, 11593 (1996).

[56] E. Rastelli, A. Tassi, and L. Reatto, Physica B+C **97**, 1 (1979).

[57] D. Bergman, J. Alicea, E. Gull, S. Trebst, and L. Balents, Nature Physics **3**, 487 (2007).

[58] P. Balla, Y. Iqbal, and K. Penc, Phys. Rev. B **100**, 140402 (2019).

[59] R. Pohle, H. Yan, and N. Shannon, Phys. Rev. B **104**, 024426 (2021).

[60] X.-P. Yao, J. Q. Liu, C.-J. Huang, X. Wang, and G. Chen, Frontiers of Physics **16**, 53303 (2021).

[61] S. Gao, G. Pokharel, A. F. May, J. A. M. Paddison, C. Pasco, Y. Liu, K. M. Taddei, S. Calder, D. G. Mandrus, M. B. Stone, and A. D. Christianson, Phys. Rev. Lett. **129**, 237202 (2022).

[62] H. Yan and J. Reuther, Phys. Rev. Res. **4**, 023175 (2022).

[63] A. Mulder, R. Ganesh, L. Capriotti, and A. Paramekanti, Phys. Rev. B **81**, 214419 (2010).

[64] M. E. Brooks-Bartlett, S. T. Banks, L. D. C. Jaubert, A. Harman-Clarke, and P. C. W. Holdsworth, Phys. Rev. X **4**, 011007 (2014).

[65] M. Pretko, Phys. Rev. B **96**, 035119 (2017).

[66] M. Pretko, Phys. Rev. B **95**, 115139 (2017).

[67] D. Lozano-Gómez, V. Noculak, J. Oitmaa, R. R. P. Singh, Y. Iqbal, J. Reuther, and M. J. P. Gingras, Proceedings of the National Academy of Sciences **121**, e2403487121 (2024), https://www.pnas.org/doi/pdf/10.1073/pnas.2403487121.

[68] D. Lozano-Gómez, O. Benton, M. J. P. Gingras, and H. Yan, An atlas of classical pyrochlore spin liquids (2024), arXiv:2411.03547 [cond-mat.str-el].
[69] K. T. K. Chung, Mapping the phase diagram of a frustrated magnet: Degeneracies, flat bands, and canting cycles on the pyrochlore lattice (2024), arXiv:2411.03429 [cond-mat.str-el].

[70] O. Benton, L. D. C. Jaubert, H. Yan, and N. Shannon, Nature Communications **7**, 11572 (2016).

[71] N. Davier, F. A. Gómez Albarracín, H. D. Rosales, P. Puljol, and L. D. C. Jaubert, Phys. Rev. B **((in press))**).

[72] F. L. Buessen and S. Trebst, Phys. Rev. B **94**, 235138 (2016).

[73] D. Kiese, F. Ferrari, N. Astrakhantsev, N. Nigge- mann, P. Ghosh, T. Müller, R. Thomale, T. Neupert, J. Reuther, M. J. P. Gingras, S. Trebst, and Y. Iqbal, Phys. Rev. Res. **5**, L012025 (2023).

[74] As a side note, the positivity of $\zeta_s$ was found on all models we have tried (not shown here), which suggests a possible relation between the saddle-point energy $E_s$ and $\Xi$.

[75] I. M. Lifshitz, Zh. Eksp. Teor. Fiz. **38**, 1569 (1960), [Sov. Phys. JETP 11, 1130 (1960)].