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

HAL
open science

# Curie-law crossover in spin liquids

Rico Pohle, Ludovic Jaubert

► To cite this version:

Rico Pohle, Ludovic Jaubert. Curie-law crossover in spin liquids. Physical Review B, 2023, 108, pp.024411. 10.1103/PhysRevB.108.024411. hal-04079565

HAL Id: hal-04079565

https://hal.science/hal-04079565

Submitted on 24 Apr 2023

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_2.png" /></figure>

Distributed under a Creative Commons Attribution 4.0 International License
Curie-law crossover in spin liquids

Rico Pohle¹,² and Ludovic D. C. Jaubert³

¹Theory of Quantum Matter Unit, Okinawa Institute of Science and Technology Graduate University, Onna-son, Okinawa 904-0412, Japan

²Department of Applied Physics, University of Tokyo, Hongo, Bunkyo-ku, Tokyo 113-8656, Japan

³CNRS, University of Bordeaux, LOMA, UMR 5798, F-33400 Talence, France

(Dated: April 24, 2023)

The Curie-Weiss law is widely used to estimate the strength of frustration in frustrated magnets. However, the Curie-Weiss law was originally derived as an estimate of magnetic correlations close to a mean-field phase transition, which – by definition – is absent in spin liquids. Instead, the susceptibility of spin liquids is known to undergo a Curie-law crossover between two magnetically disordered regimes. Here, we study the generic aspect of the Curie-law crossover by comparing a variety of frustrated spin models in two and three dimensions, using both classical Monte Carlo simulations and analytical Husimi tree calculations. Husimi tree calculations fit remarkably well the simulations for all temperatures and almost all lattices. We also propose a Husimi Ansatz for the reduced susceptibility $\chi T$, to be used in complement to the traditional Curie-Weiss fit in order to estimate the Curie-Weiss temperature $\theta_{cw}$. Applications to materials are discussed.

# I. INTRODUCTION

The Curie-Weiss law is a simple and useful tool to estimate the behaviour of the susceptibility $\chi$ for conventional magnets at high temperatures¹⁻⁵

$$ \chi = \frac{\mathcal{C}}{T - \theta_{cw}}, \qquad (1) $$

with $\mathcal{C}$ the Curie constant, and $\theta_{cw}$ the Curie-Weiss temperature. In a Landau mean-field treatment⁶, $|\theta_{cw}|$ represents the transition temperature. The sign of $\theta_{cw}$ indicates dominant ferromagnetic ($\theta_{cw} > 0$) or antiferromagnetic ($\theta_{cw} < 0$) interactions, while the limit $\theta_{cw} \to 0$ represents the susceptibility of a paramagnet, given by the Curie law, $\chi = \mathcal{C}/T$. For more details about the application of the Curie-Weiss law in susceptibility measurements, we refer the reader to the recent tutorial by Mugiraneza & Hallas [5].

In frustrated magnets, the Curie-Weiss temperature is often used to measure the “frustration index”⁷

$$ f = \frac{|\theta_{cw}|}{T^*}, \qquad (2) $$

by comparing the transition, or freezing, temperature of a material, $T^*$, to its mean-field expectation, $|\theta_{cw}|$, for an unfrustrated system. Large values of $f$ account for strong frustration in the system. For a spin liquid where $T^* \to 0^+$ theoretically, the frustration index diverges. Being a priori readily accessible to experiments, this quantity $f$ has become a convenient tool to gauge how frustrated a system is.

But as many successful, broadly used indicators, a few shortcomings are inevitable. Deviations from the standard Curie-Weiss law have been studied in a variety of magnetic systems, such as spin glasses⁸,⁹, the pyrochlore molybdate $Y_2Mo_2O_7$¹⁰, the valence bond glass $Ba_2YMoO_6$¹¹, or Kitaev materials with strong spin-orbit coupling¹², to cite but a few. For example in anisotropic

lattices, the high-temperature Curie constant and low-temperature transition temperature may be set by different energy scales, giving rise to an artificially large parameter $f$ even when the system is barely frustrated¹³,¹⁴.

In spin liquids, this deviation has been rationalized as the onset of a Curie-law crossover¹⁵,¹⁶ between the universal high-temperature Curie law and a low-temperature, model-specific, spin-liquid Curie law¹⁶⁻¹⁹. The problem is that fitting the susceptibility of spin liquids with a Curie-Weiss law always gives an answer, but not necessarily the right one, as illustrated for the Ising kagome antiferromagnet in Fig. 1. Beyond the traditional difficulties to measure the Curie-Weiss temperature⁵,¹², frustration precisely prevents the phase transition in spin liquids that would justify the Curie-Weiss fit as a mean-field approximation of a scaling law with critical exponent $\gamma = 1$. Eq. (1) is only valid at high temperature as a first order perturbation of the Curie law. And whether this high-temperature regime is accessible to experiments then becomes an important question⁵,¹². Internal energy scales such as a single-ion crystal field, a band gap, the structural distortion of the lattice or even the melting of the materials might prevent access to the necessary high temperatures. In that case, the values of the Curie constant and Curie-Weiss temperature strongly depend on the temperature range of the fitting procedure¹⁶,²⁰. The latter can even change sign when the exchange coupling is particularly small (see e.g. Refs. [21 and 22] for $Dy_2Ti_2O_7$). And as a high-temperature expansion of the susceptibility, the Curie-Weiss fit is not designed to capture the spin-liquid behaviour at low temperatures.

To summarise the issue, applying the Curie-Weiss fit to frustrated magnets means applying a method that has been derived around a mean-field critical point, to a class of systems where this critical point is absent by definition.

In this paper we want to rationalise this conceptual divergence of viewpoints. Is it possible to quantify how the magnetic susceptibility deviates from the Curie-Weiss law, not just for a specific model but for frustrated mag-
nets in general? In particular, can we identify generic features? Practically, understanding the limits of the Curie-Weiss fit will help estimate the appropriate temperature window to measure the Curie-Weiss temperature, and what to do when this window is not experimentally available.

For the sake of generality, we will at first focus on a variety of traditional frustrated lattices in two and three dimensions [Fig. 2], made either of triangular or tetrahedral frustrated unit cells, with Ising spins. Our motivation here is not to study these models individually. That has already been done extensively in the literature; see e.g. the following references for the two dimensional triangular$^{18,23,24}$, kagome$^{19}$, square-kagome$^{13,25,26}$, checkerboard$^{27}$ and ruby$^{28}$ lattices, and for the three dimensional trillium$^{29}$, hyperkagome$^{30}$, and pyrochlore$^{16,17,31,32}$ lattices. Instead, we will compare these models together, understand why similarities appear between some of them, and build an overall intuition for the phenomenon of the Curie-law crossover in spin liquids.

On the theoretical front, comparing unbiased classical Monte Carlo simulations to the analytical Husimi-Tree approximation shows that thermodynamic quantities are, to a large extent, independent of the lattice dimension, and even of the structure of the lattice beyond the minimal frustrated unit cells [Fig. 5]. What essentially matters is simply the type of frustrated unit cell (triangle, tetrahedron, ...) and the local connectivity between them.

On the experimental front, one of our take-home messages is that the reduced susceptibility $\chi T$ (that is frequently used in chemistry) is a very useful complement to the inverse susceptibility $1/\chi$ for frustrated magnets. The Curie-law crossover is especially transparent in this quantity, between two horizontal asymptotic lines. $\chi T$ thus immediately tells us (i) how far we are from the high-temperature Curie law, and (ii) the presence or absence of a low-temperature region where the onset of a spin-liquid regime can be expected above a potential ordering or freezing temperature $T^*$. In order to describe the Curie-law crossover in its entirety, we introduce the following fitting Ansatz [Fig. 1(b)]:

$$ \chi T^{\text{fit}} = \frac{1 + b_1 \exp[c_1/T]}{a + b_2 \exp[c_2/T]}, \quad (3) $$

inspired by the above analogy between disparate models and Husimi-tree calculations. This empirical Ansatz provides a complementary estimate of the Curie constant and Curie-Weiss temperature,

$$ \mathcal{C} = \frac{1 + b_1}{a + b_2} \quad \& \quad \theta_{\text{cw}} = \frac{b_1 c_1}{1 + b_1} - \frac{b_2 c_2}{a + b_2}, \quad (4) $$

that is not based on a high-temperature expansion. Hence, if Eq. (4) agrees with values obtained from a Curie-Weiss fit, then it is reasonable to consider them as an accurate description of the material. On the

<figure><img src="image_3.png" /><figcaption>FIG. 1. Curie-law crossover in spin liquids. Both panels compare susceptibility results from Monte Carlo (MC) simulations of the Ising antiferromagnet on the kagome lattice [Eq. (5) and Fig. 2(b)] (open black circles) to their corresponding results on the Husimi tree “HT(3,2)” [Eq. (C19) and Fig. 3(a)] (solid red line). (a) Inverse susceptibility $1/\chi$ on a linear temperature scale. The Curie-Weiss fit has been obtained from fitting data for $2 < T/J < 10$ (blue dashed line), giving $\theta_{\text{cw}}^{\text{fit}} \approx -5.6J$ different from the known exact value of $-4J$ (red dashed line, obtained from a Curie-Weiss fit of the HT(3,2) curve). (b) Same results are plotted for the reduced susceptibility $\chi T$ on a semi-logarithmic plot. The Husimi tree “HT(3,2)” result matches quantitatively with MC simulations, and shows the crossover between two different Curie constants at high-T ($\mathcal{C}_{\infty} = 1$ in paramagnetic phase) and low-T ($\mathcal{C}_0 = 0.2$ in spin liquid phase), corresponding to two different Curie laws. If the fit is done in the intermediate crossover region $(2 - 10|J|)$, which is typically the region accessible to experiments (see Section V), the resulting Curie-Weiss law quickly deviates from simulations.</figcaption></figure>

other hand if there is a noticeable mismatch, then it is likely that experimental data have not reached the high-temperature regime where the Curie-Weiss law is valid.

The remainder of this Article is structured as follows. In Sec. II, we introduce the models of classical spin liquids, defined on a variety of frustrated lattices in two and three dimensions [Fig. 2]. These models will be solved numerically with classical Monte Carlo simulations and analytically on their corresponding Husimi trees [Fig. 3].

In Sec. III, we present thermodynamic quantities for all spin liquids considered in this article and discuss analogies and signatures of their Curie-law crossover. In particular, we discuss the reason for the very good match
<figure><img src="image_4.png" /><figcaption>FIG. 2. Corner-sharing lattices in two and three dimensions involve different lengths of minimal loops $\mathcal{L}$ between frustrated cells. (a) triangular lattice ($\mathcal{L} = 3$), (b) kagome lattice ($\mathcal{L} = 6$), (c) square-kagome lattice ($\mathcal{L} = 4$), (d) checkerboard lattice ($\mathcal{L} = 4$), (e) ruby lattice ($\mathcal{L} = 3$), (f) trillium lattice ($\mathcal{L} = 5$), (g) hyperkagome lattice ($\mathcal{L} = 10$), and (f) pyrochlore lattice ($\mathcal{L} = 6$). Thermodynamic observables for each lattice (see Fig. 5) have been obtained numerically with classical Monte Carlo simulations of Hamiltonian $\mathcal{H}$ [Eq. (5)], as described in Appendix B. While commonly referred to “edge sharing” in the literature, we describe the triangular lattice as corner sharing to emphasise its analogy with the trillium lattice in three dimensions and the corresponding Husimi tree HT(3,3) in Fig. 3(c). Numbers on lattice sites indicate the Manhattan distance $\ell$, used in Fig. 6.</figcaption></figure>

between Monte-Carlo simulations and Husimi-tree calculations, despite the different lattice structure.

In Sec. IV, we discuss the limitations of the conventional Curie-Weiss fit, showing the advantage to use the reduced susceptibility $\chi T$. We introduce and benchmark the Husimi Ansatz [Eq. (3)] to numerical simulations.

In Sec. V, we apply this Ansatz to experimental data for the pyrochlore NaCaNi$_2$F$_7$ [33], the square-kagome KCu$_6$AlBiO$_4$(SO$_4$)$_5$Cl [34] and the spiral spin liquid FeCl$_3$ [35].

In Sec. VI, we conclude with a brief summary and implications for future experiments on spin liquid materials.

Details on the lattice geometries, Monte Carlo simulations, Husimi tree calculations, connection to Coulomb gauge field theory, and structure factors are given in Appendices A, B, C, D and E respectively. In particular, we refer the reader interested in Husimi trees to Appendix C 4 where a couple of non-trivial exact results are derived in presence of local easy-axes anisotropy.

## II. MODELS AND METHODS

### A. The Ising model

In Sec. II and III, we focus on thermodynamic properties of minimal spin-liquid models,

$$ \mathcal{H} = J \sum_{\langle ij \rangle} \vec{S}_i \cdot \vec{S}_j, \qquad (5) $$

for Ising spins $\vec{S}_i = \sigma_i \vec{e}_i$, with $\sigma_i = \pm 1$, and nearest-neighbour coupling $J$, applied to a variety of lattices, as shown in Fig. 2. We shall consider two types of Ising spins, either collinear along the same global z-axis $\vec{e}_z$, or oriented along their local easy axis $\vec{e}_i$ attached to the sublattice of site $i$. The latter is motivated from single-ion anisotropy, as found, for example, in kagome materials like Dy$_3$Mg$_2$Sb$_3$O$_{14}$$^{36}$ and spin ices like Dy$_2$Ti$_2$O$_7$ and Ho$_2$Ti$_2$O$_7$ on the pyrochlore lattice$^{37,38}$, and EuPtSi$^{29,39,40}$ on the trillium lattice. We shall refer to each system as “global-axis” and “local-axis” models, as illustrated in Fig. 4. All local easy axes relevant for this work are defined in Appendix A. Global-axis and local-axis models are equivalent, up to a simple rescaling of the coupling constant $J^{41,42}$

$$ J_{\text{local}} = J_{\text{global}} (\vec{e}_i \cdot \vec{e}_j), \qquad (6) $$
<figure><img src="image_5.png" /><figcaption>FIG. 3. Husimi Trees for various corner sharing lattices. Frustrated cells from real lattices [Fig. 2] are arranged on the Husimi tree (HT), with the central cell in red, the 1st shell in blue and the 2nd shell in green. (a) HT(3,2): for the kagome and hyperkagome lattice, with corner sharing triangular plaquettes. (b) HTS: for the square-kagome lattice. (c) HT(3,3): for the triangular and trillium lattice, were three triangular plaquettes share one corner. (d) HT(4,2): for the checkerboard, ruby and pyrochlore lattice, which is made of corner-sharing square/tetrahedron plaquettes.</figcaption></figure>

where *i* and *j* are two neighbouring sites. For lattices considered here, the scalar product $(\vec{e}_i \cdot \vec{e}_j)$ is the same for all neighbouring pairs, which means that the energy, specific heat and entropy of the two models are the same up to rescaling (6). However, magnetic quantities such as the susceptibility differ. In this work, the exchange coupling is always antiferromagnetic $J_{global} > 0$ (ferromagnetic $J_{local} < 0$) for global-axis (local-axis) models, in order to stabilise a spin-liquid ground state. From now on, all energies and temperatures are given in units of $J_{global} = 1$, understanding that the rescaling of Eq. (6) is always applied for local-axis models.

<figure><img src="image_6.png" /><figcaption>FIG. 4. While Ising models usually consider collinear spins (a),(c), the crystal field in materials may impose a local easy axis (b), (d) respecting the symmetry of the magnetic-ion environment, as illustrated here for the kagome and pyrochlore lattice. All local easy axes are defined in Appendix A.</figcaption></figure>

## B. Spin liquids on the Husimi tree

The frustrated Ising model [Eq. (5)] on corner-sharing lattices [see Fig. 2] can be solved, regardless of its physical dimension, by numerical methods such as classical Monte Carlo simulations [see Appendix B]. On the analytical front, however, the question is more delicate. Since correlations play a major role, one needs a method

beyond standard mean-field theory, but nonetheless valid for frustrated models across different dimensions. The Husimi-Tree (HT) calculation precisely fits our needs, by incorporating the local frustrated constraints of the lattice, irrespectively of its dimension. HT recursively extends from a central frustrated cell – e.g. a triangle or a tetrahedron – into a non-reciprocal lattice, without any internal loop beyond the frustrated cell [Fig. 3]. As a consequence, its boundary is of comparable size to the volume of the bulk$^{51,52}$ and the HT remains a mean-field approach. It is thus inaccurate at critical points, except above their upper critical dimensions$^{53–55}$. But since we explicitly study models away from phase transitions, we expect pertinent analytical insights from the HT, spurred on by encouraging results on frustrated systems in the literature$^{16,53–62}$. Technical aspects of the HT method are explained in Appendix C.

We will compare a variety of physical lattices, with different numbers of internal loops and frustrated unit cells [Fig. 2], to their pseudo-lattice counterparts on the Husimi tree, which do not have any internally closed loops [Fig. 3]. Let us define $\mathcal{L}$ as the smallest internal loop formed by frustrated cells on the physical lattice. We relate all physical lattices, as introduced in Fig. 2, to their corresponding HT trees, according to the number of sites per frustrated cell and their connectivity:

(i) HT(3,2) [Fig. 3(a)] contains three sites in the frustrated cell, where each site is connected between two cells. It relates to the kagome ($\mathcal{L} = 6$), square-kagome ($\mathcal{L} = 4$) and hyperkagome ($\mathcal{L} = 10$) lattice. Considering the complexity of the frustrated cell in the square-kagome lattice, we also included the Husimi tree HTS [Fig. 3(b)] to improve the mean-field approximation.

(ii) HT(3,3) [Fig. 3(c)] contains three sites in the frustrated cell, where each site is connected between three cells, and relates to the triangular ($\mathcal{L} = 3$) and trillium ($\mathcal{L} = 5$) lattice.
<table><caption>TABLE I. Residual entropy S |<sub>T→0</sub> and spin liquid Curie constant C₀, obtained from Monte Carlo (MC) simulations and Husimi tree (HT) calculations. Table cells including two lines for C₀ display results for global-axis (first line) and local-axis models (second line). HT calculations of C₀ are detailed in Appendix C4. The column “other methods” compares the HT estimate, also known as Pauling entropy, to exact results when available (except for the 3D pyrochlore lattice obtained from series expansion). As a side remark, one should be aware that the Pauling entropy is not always a lower bound of the residual entropy on the corresponding real lattices (see Ref. [50] and Appendix C5).</caption><thead><tr><th rowspan="2">Lattice</th><th colspan="3">S |<sub>T→0</sub></th><th colspan="2">C<sub>0</sub> ≡ χT |<sub>T→0</sub></th></tr><tr><th>Monte Carlo</th><th>Husimi Tree</th><th>other methods</th><th>Monte Carlo</th><th>Husimi Tree</th></tr></thead><tbody><tr><td>kagome</td><td>0.502(1)</td><td><sup>1</sup>&frasl;<sub>3</sub> ln <sup>9</sup>&frasl;<sub>2</sub> ≈ 0.5014<sup>43</sup></td><td>0.50183<sup>44</sup></td><td>0.201(1)<br>1.988(1)</td><td>1/5<br>2</td></tr><tr><td>hyperkagome</td><td>0.502(1)</td><td><sup>1</sup>&frasl;<sub>3</sub> ln <sup>9</sup>&frasl;<sub>2</sub> ≈ 0.5014<sup>43</sup></td><td>n/a</td><td>0.200(1)<br>1.500(1)</td><td>1/5<br>3/2</td></tr><tr><td>square-kagome</td><td>0.504(1)<sup>13</sup></td><td><sup>1</sup>&frasl;<sub>6</sub> ln <sup>41</sup>&frasl;<sub>2</sub> ≈ 0.5034<sup>13</sup></td><td>n/a</td><td>0.203(1)</td><td>0.2028</td></tr><tr><td>triangular</td><td>0.323(2)</td><td>ln <sup>3</sup>&frasl;<sub>2</sub> ≈ 0.4055<sup>29</sup></td><td>0.323066<sup>45,46</sup></td><td>0.162(8)</td><td>1/7</td></tr><tr><td>trillium</td><td>0.392(1)<sup>29</sup></td><td>ln <sup>3</sup>&frasl;<sub>2</sub> ≈ 0.4055<sup>29</sup></td><td>n/a</td><td>0.135(1)<br>0.969(1)</td><td>1/7<br>1</td></tr><tr><td>ruby</td><td>0.194(1)</td><td><sup>1</sup>&frasl;<sub>2</sub> ln <sup>3</sup>&frasl;<sub>2</sub> ≈ 0.2027<sup>47</sup></td><td>n/a</td><td>0.0</td><td>0</td></tr><tr><td>checkerboard</td><td>0.216(1)</td><td><sup>1</sup>&frasl;<sub>2</sub> ln <sup>3</sup>&frasl;<sub>2</sub> ≈ 0.2027<sup>47</sup></td><td><sup>3</sup>&frasl;<sub>4</sub> ln <sup>4</sup>&frasl;<sub>3</sub> ≈ 0.2158<sup>48</sup></td><td>0.0</td><td>0</td></tr><tr><td>pyrochlore</td><td>0.206(1)</td><td><sup>1</sup>&frasl;<sub>2</sub> ln <sup>3</sup>&frasl;<sub>2</sub> ≈ 0.2027<sup>47</sup></td><td>0.205006(9)<sup>49</sup></td><td>0.0<br>2.002(1)<sup>17</sup></td><td>0<br>2</td></tr></tbody></table>

(iii) HT(4,2) [Fig. 3(d)] contains four sites in the frustrated cell, where each site is connected between two cells, and relates to the checkerboard ($\mathcal{L} = 4$), ruby ($\mathcal{L} = 3$) and pyrochlore ($\mathcal{L} = 6$) lattice.

The similarity between a given lattice and its Husimi tree, taken individually, makes sense – except maybe for the triangular lattice, which will be discussed separately in Section III D. In this set up we shall investigate the Curie-law crossover by comparing thermodynamic quantities between the physical lattice (as obtained by classical Monte Carlo simulations) and their corresponding pseudo lattice on the Husimi tree in the next section.

### III. THE CURIE-LAW CROSSOVER

The Curie-law crossover describes the evolution of the magnetic susceptibility between two different Curie laws¹⁶, whose origin becomes obvious when considering the reduced susceptibility $\chi T$:

$$ \begin{aligned} \chi T &= \frac{1}{N} \sum_{i,j} [\langle \vec{S}_i \cdot \vec{S}_j \rangle - \langle \vec{S}_i \rangle \langle \vec{S}_j \rangle] \\ &= 1 + \frac{1}{N} \sum_{i \neq j} \langle \vec{S}_i \cdot \vec{S}_j \rangle. \end{aligned} \quad (7) $$

In magnetically disordered systems, as studied here, $\langle \vec{S}_i \rangle = 0$ for all temperatures, while translational invariance implies additionally that

$$ \chi T = 1 + \sum_{i \neq 0} \langle \vec{S}_0 \cdot \vec{S}_i \rangle = \sum_i \langle \vec{S}_0 \cdot \vec{S}_i \rangle, \quad (8) $$

where $\vec{S}_0$ is an arbitrary “central” spin. In a paramagnet with uncorrelated spins, Eq. (8) gives the Curie constant

$$ \mathcal{C}_{\infty} \equiv \chi T |_{T \to \infty} = 1. \quad (9) $$

At zero temperature, the Curie constant is renormalised by the correlations of the spin liquid

$$ \mathcal{C}_0 \equiv \chi T |_{T \to 0} = \sum_i \langle \vec{S}_0 \cdot \vec{S}_i \rangle_{T \to 0}. \quad (10) $$

In fact, C₀ is nothing less than the integration of spin correlations in real space, with C₀ smaller (greater) than 1 for dominating antiferromagnetic (ferromagnetic) correlations.

#### A. Thermodynamics

Fig. 5 displays thermodynamic observables: energy E, specific heat C, entropy S and reduced susceptibility $\chi T$, obtained by simulating the Hamiltonian $\mathcal{H}$ [Eq. (5)] with classical Monte Carlo simulations for the physical lattices [Fig. 2] and analytical calculations on their corresponding Husimi trees [Fig. 3]. As explained in the introduction, these systems have often been studied in the literature; see e.g. Refs. 13,16–19,23–32,43–49. Our interest here is not to study them individually, but to see how their thermodynamic properties compare to each other. In particular, classical spin liquids are known for their residual entropy as $T \to 0^+$, that measures the degeneracy of the spin-liquid ground state. It can be categorized into three groups43–49 (corresponding to the three columns in
<figure><img src="image_7.png" /></figure>

FIG. 5. Thermodynamic signatures of the Curie-law crossover in spin liquids. Comparison of the energy *E*, specific heat *C*, entropy *S* and reduced susceptibility χ*T* per spin for results obtained from classical Monte Carlo simulations (open symbols) on the physical lattices [Fig. 2] and analytical calculations (solid black lines) on their corresponding Husimi trees [Fig. 3]. Observables are shown on a semi-logarithmic plot. (a) Lattices with triangular cells, where each site belongs to two frustrated cells, HT(3,2): kagome, square-kagome and hyperkagome. (b) Lattices with triangular cells, where each site belongs to three frustrated cells, HT(3,3): trillium and triangular. (c) Lattices with tetrahedral cells, HT(4,2): checkerboard, ruby and pyrochlore. Results are given for global-axis and local-axis models, respectively labeled “global” and “local”, as explained in Sec. II. All systems perform a crossover from a high-temperature paramagnetic regime into a low-temperature classical spin-liquid regime. This is seen by two different Curie laws at high and low temperatures in the reduced susceptibility χ*T* [Eq. (8)]. Technical details on simulations and analytics are given in Appendix B and C, respectively.

Fig. 5), (i) kagome, square-kagome and hyperkagome lattices with $S(T \to 0) \approx 0.5$, (ii) triangular and trillium lattice with $S(T \to 0) \sim 0.3-0.4$, and (iii) ruby, checkerboard and pyrochlore lattices with $S(T \to 0) \approx 0.2$. The HT estimate of the residual entropy is also known as Pauling entropy, which, as a side-note, is not always a lower bound [see Appendix C 5].

The behavior of the entropy is accompanied by a change in magnetic correlations from a high-temperature regime with $\mathcal{C}_\infty = 1$ to a model-dependent value $\mathcal{C}_0$ at low temperatures [see also Table I]. The low-temperature Curie constant $\mathcal{C}_0$ is not universal, making its value a characteristic property of the underlying spin liquid.

In some models, the value of $\mathcal{C}_0$ is easy to understand. For the ruby, checkerboard and pyrochlore lattice with

global axis spins, $\mathcal{C}_0$ is zero [Fig. 5 (c)]. This is because their ground state respects the so-called ice rules$^{41,42}$ with two up spins and two down spins per frustrated cell. The magnetization $M = |\sum_i \vec{e}_z \sigma_i|$ is thus not only globally zero on average, $\langle M \rangle = 0$, but also locally zero for all frustrated units. No fluctuations of the magnetization are allowed in the spin liquid, resulting in $\langle M^2 \rangle = 0$ and thus a vanishing reduced susceptibility. In other words, we get $\mathcal{C}_0 = 0$ as can be expected for any system with a zero-magnetization plateau. For triangular frustrated units, the opposite reasoning applies because the magnetization cannot be canceled with three collinear Ising spins. Magnetic fluctuations persist down to zero temperature, and $\langle M^2 \rangle$ and $\chi T$ remain finite.
<figure><img src="image_8.png" /><figcaption>FIG. 6. Absolute value of the real-space spin-spin correlation length $|D(\ell)|$ at low temperature, deep in the spin liquid regime of $\mathcal{H}$ [Eq. (5)], as obtained from classical Monte Carlo simulations for global-axis Ising spins on their physical lattices [Fig. 2], and their corresponding Husimi tree (HT). (a) Exponential decay $D(\ell) \sim e^{-\ell/\xi}$ on the kagome ($\xi = 1.10(2)$), hyperkagome ($\xi = 1.03(2)$), square-kagome ($\xi = 1.18(2)$) and ruby lattice ($\xi = 0.76(3)$) compare semi-quantitatively well with the exponential decay on the Husimi tree ($\xi_{HT} = 0.91$) [Eq. (15)]. (b) Algebraic decay $D(\ell) \sim 1/\ell^\alpha$ on the triangular ($\alpha = 0.476(10)$), checkerboard ($\alpha = 2$), and pyrochlore ($\alpha = 3$), lattice. On the checkerboard and pyrochlore lattices, correlations are known to decay algebraically (see Section III C) but follow the HT exponential decay for the first three or four nearest neighbours, i.e. over a distance larger than $\xi_{HT}$. Correlations on the trillium ($\xi = 1.55(2)$) and triangular lattice deviate more strongly from HT expectations, which we believe causes the small, but visible, mismatch of thermodynamic quantities in Fig. 5(b). Note that for a relevant comparison between physical lattices and HT, we used the Manhattan distance $\ell$, defined on paths for each lattice as shown in Fig. 2. For lattices in panel (b), the Manhattan distance is also the Euclidian distance. Panels (a) and (b) are respectively on a semi-log and log-log plot.</figcaption></figure>

Remarkably, thermodynamic observables match well within each group of lattices, despite their different physical dimensions. It was already known$^{43–49}$ that some models had very similar residual entropies as $T \to 0^+$. Here this similarity is further illustrated with the value of the spin-liquid Curie law $C_0$ [see Table I]. But more importantly, thermodynamic quantities are essentially the same for all temperatures within each group of models. For example, the 2D square-kagome model compares well with the 2D kagome, as recently noticed for quantum spins–1/2 [26], but also the 3D hyperkagome, while the 2D ruby model matches with 3D pyrochlore for all temperatures. Furthermore, thermodynamic observables for each group are well reproduced by their corresponding HT, suggesting that correlations barely depend on the physical dimension of the lattice. In the following we will try to understand why.

## B. Husimi tree sets the correlation length

As seen in Eq. (10), $C_0$ corresponds to the integration of real-space correlations in the spin liquid [see also Table I]. Let us consider HT(3,2) whose $C_0 = 1/5$. This value deviates from Monte-Carlo results on the kagome and hyperkagome lattice within less than 1 %. For the square-kagome lattice, the mismatch drops from 2% to 0.1% by including a more evolved version of the Husimi

tree (see HTS in Fig. 3(b)), which contains a larger frustrated unit cell and includes internal loop lengths of 4 sites. Such a trend suggests the presence of a particularly small correlation length $\xi$ in these systems.

To confirm our suggestion, let us define spin-spin correlations on the HT:

$$D(\ell) = \langle \vec{S}_0 \cdot \vec{S}_\ell \rangle = \langle \sigma_0 \sigma_\ell \rangle, \quad (11)$$

assuming that all spins are collinear along a global axis $\vec{S}_i = \sigma_i \vec{e}_z$. The fact that (i) there is no closed loop in the HT (beyond the size of the frustrated unit), (ii) the Hamiltonian is invariant under time-reversal symmetry, and (iii) the HT is by definition locally the same at each vertex, allows us to formulate an exact expression for the spin-spin correlations

$$ \begin{aligned} D(\ell) &= \langle \sigma_0 \sigma_\ell \rangle = \langle \sigma_0 \sigma_1^2 \sigma_2^2 \dots \sigma_{\ell-1}^2 \sigma_\ell \rangle \\ &= \langle \sigma_0 \sigma_1 \rangle \langle \sigma_1 \sigma_2 \rangle \dots \langle \sigma_{\ell-1} \sigma_\ell \rangle \\ &= \langle \sigma_0 \sigma_1 \rangle^\ell. \end{aligned} \quad (12) $$

The nearest-neighbor spin-spin correlation averaged over the ensemble of ground states can be easily calculated. And it turns out to be the same result for the three kinds of Husimi trees, HT(3,2), HT(3,3) and HT(4,2):

$$\langle \sigma_0 \sigma_1 \rangle = -1/3. \quad (13)$$

This means that correlations decay exponentially on
Husimi trees, following the formula

$$D(\ell) = \left(-\frac{1}{3}\right)^{\ell} = (-1)^{\ell} e^{-\ell \ln 3}, \quad (14)$$

for all Husimi trees considered here, giving a correlation length

$$\xi_{HT} = (\ln 3)^{-1} \approx 0.91. \quad (15)$$

More generally, for a Husimi tree whose frustrated units are made of $N_u$ Ising spins fully connected between each other via antiferromagnetic couplings, the correlation length in the degenerate ground state is

$$\xi_{HT} = \frac{1}{\ln N_u}, \quad \text{if } N_u \text{ is odd,} \quad (16)$$

$$\xi_{HT} = \frac{1}{\ln(N_u - 1)}, \quad \text{if } N_u \text{ is even.} \quad (17)$$

Eq. (15) means that correlations decay typically over the nearest-neighbor distance in Husimi trees. This length scale is smaller than any loop in the real lattice, suggesting that correlations in real lattices may decay in a similar way at short distances. Monte Carlo simulations confirm this assumption on the kagome, hyperkagome³⁰, square-kagome¹³ and ruby²⁸ lattice at low temperatures [Fig. 6(a)], whose correlation lengths are roughly the same as $\xi_{HT}$. Since the correlation length is expected to decrease upon heating, this short correlation length is consistent with the success of the HT approximation over the whole temperature range for global- and local-axis models alike.

For the sake of clarity, we should point out that the value of $\mathcal{C}_0$ is not coming from a cutoff of the correlations beyond $\xi_{HT}$. Indeed, it would be tempted to see classical spin liquids as an ensemble of independent clusters of superspins (on each triangle or tetrahedron), and the spin-liquid Curie law as a form of superparamagnetism, as observed with ferromagnetic nanoparticles⁶³. However, we cannot recover the value $\mathcal{C}_0 = 0.2$ for kagome-type systems from such an argument. Appendix C4 shows that the resulting error scales like $-\frac{6}{5} \left(-\frac{2}{3}\right)^{L+1}$ on a Husimi tree of L layers. Even if correlations ultimately vanish at long distance, the cutoff necessary to approximately recover the value of $\mathcal{C}_0$ is much larger than $\xi_{HT}$. To understand the similarity between simulations and analytics, it would be more accurate to see the paths connecting the central spin to the many spins on layer L on the infinite-dimensional Husimi tree of Fig. 3 as virtual paths of correlations connecting a pair of Lth nearest neighbors on the corresponding real lattice of Fig. 2. This picture is nearly exact up to the nth nearest neighbor before closing the minimal loop of size L on the real lattice ($n = L/2 - 1$), which is why deviations between Monte Carlo and Husimi tree grow inversely with L in Fig. 6: first hyperkagome, then kagome and finally square-kagome.

### C. Coulomb field theory and flat bands

On the other hand, correlations on the checkerboard and pyrochlore lattices are algebraic at low temperature, scaling like $1/r^d$ [64], with d the physical dimension of the lattice. Their angular dependence is dipolar though, which means that the integration of these correlations over the entire system in Eq. (10) does not diverge, and $\mathcal{C}_0$ is well defined. The dipolar nature of these correlations comes from the fact that their ground states are ice models, described by an emergent Coulomb field theory²⁷. With respect to the exponential decay of the HT, these algebraic correlations only differ beyond the third or fourth neighbour; see the comparison to the black curve on Fig. 6(b). In that sense the correlation length $\xi_{HT}$ remains effectively relevant at short distances. That being said, one would have been forgiven to expect larger corrections to $\mathcal{C}_0$ coming from the long-range algebraic tail. Here again we are left with the question: why are these corrections so small ?

For models with a global axis, $\mathcal{C}_0$ is known to be exactly zero (see discussion in Section III A); the ice rules prevent magnetic fluctuations for all tetrahedra, and thus conveniently prevent any corrections. But this does not explain the match of Fig. 5(c) for the local-axis pyrochlore model, a.k.a. spin ice, where $\mathcal{C}_0 \approx 2.0$. Magnetic fluctuations are allowed in the spin-ice ground state. Additionally, since the spin-ice model is ferromagnetic, the sum of Eq. (8) contains mostly positive terms, as opposed to the alternating series encountered for the integration of correlations in antiferromagnets [Appendix C4]. For the latter, potential corrections coming from algebraic correlations would partially cancel out; while they would a priori add up in the ferromagnetic model. This suggests that an alternative point of view is necessary.

Let us temporarily step away from Husimi trees and consider the other facet of spin ice, as a U(1) Coulomb gauge field. As mentioned previously this gauge-field texture comes from the ice rules, that can be rewritten as a divergence-free constraint on the magnetisation field $\vec{M}$ [27]. But spin ice is not the only model supporting this type of texture. The ground state of the pyrochlore antiferromagnet with classical Heisenberg spins is a U(1)×U(1)×U(1) Coulomb gauge field that has often been described as three copies of spin ice²⁷,⁶⁴. The susceptibility of these divergence-free fields is readily available using the Self-Consistent Gaussian Approximation (SCGA). It means that with the proper normalisation, SCGA offers an alternative approach to compute $\mathcal{C}_0$ and $\mathcal{C}_\infty$ [see Appendix D]. In particular it tells us that the ratio $\mathcal{C}_\infty/\mathcal{C}_0$ is due to the topology of the magnetic band structure of the pyrochlore lattice⁶⁵⁻⁶⁸; the ground state is composed of two degenerate flat bands, accounting for half ($\mathcal{C}_\infty/\mathcal{C}_0 = 1/2$) of the total number of bands.

To summarise, since $\mathcal{C}_0$ comes from the integration of correlations [Eq. (10)], it is remarkable that algebraic correlations in real lattices give almost the same
result as exponential correlations in Husimi trees [see pyrochlore and checkerboard results in Fig. 5]. This is because $C_0 = 0$ is protected by the absence of local fluctuations for global-axis models, while $C_0 \approx 2.0$ is a direct consequence of the topology of the band structure for local-axis models.

<figure><img src="image_9.png" /><figcaption>FIG. 7. Signatures of the Curie-law crossover in coulombic spin liquids. (a)–(c) equal-time structure factor $S(q)$ [Eq. (E1)] of $H$ [Eq. (5)] for Ising spins in their global axis on the checkerboard lattice, obtained from classical Monte Carlo simulations. The pinch-point gets broader upon heating. (d) Temperature-dependent full width at half maximum (FWHM) of pinch-points. FWHM has been obtained from a Lorentzian fit for line cuts of the pinch point along its singular $q_x$ direction (see inset of (d)). The FWHM illustrates the Curie-law crossover in a similar way as the reduced susceptibility $\chi T$.</figcaption></figure>

Before closing our discussion on the checkerboard and pyrochlore lattices, let us take advantage of these dipolar correlations, whose signatures in the equal-time structure factor [see Appendix E] present sharp, singular features known as pinch points$^{27,69,70}$. Upon heating, these singular features broaden as topological-charge excitations disrupt the Coulomb field [Fig. 7 (a)–(c)]$^{71}$. By measuring their breadth, pinch points offer a quantitative way to measure the establishment of the spin liquid. Fig. 7(d) shows the full width at half maximum (FWHM) of the pinch point on the checkerboard lattice as a function of temperature. Our point is that the Curie-law crossover, as seen in $\chi T$, is able to grasp the evolution of FWHM, i.e. the build up of the spin liquid. And while only a fraction of spin liquids have characteristic, singular, patterns such as pinch points [see e.g. Appendix E for the hy-

perkagome and trillium lattices], the Curie-law crossover is a generic property of all spin liquids. This vindicates the Curie-law crossover as a useful signature of the onset of a spin liquid, and the reduced susceptibility $\chi T$ as a suitable observable to measure it.

### D. The triangular and trillium lattice

Let us now consider two systems with noticeably different geometries; the triangular and trillium lattice. While the latter is three dimensional and made of cornersharing triangles, the former is two dimensional and usually seen as made of edge-sharing triangles. From the view point of Husimi trees, $HT(3,3)$ is clearly a reasonable approximation for the trillium lattice, with each spin belonging to three triangles. But, even if less conventional, it can also be used for the triangular lattice$^{62,72}$, since each spin can similarly be seen as shared by three triangles (see colored lattice in Fig. 2(a)). The obvious caveat of this choice of Husimi tree (made of 3 spins) is that loops that are ignored, are of the same size than the frustrated triangular unit cell itself. However, by direct comparison between MC and $HT(3,3)$ results in Fig. 5(b), the reduced susceptibility, $\chi T$, of the two antiferromagnetic models overlap with a quantitative difference appearing only below $T \lesssim 1$ [Fig. 8].

<figure><img src="image_10.png" /><figcaption>FIG. 8. Reduced susceptibility $\chi T$ of the Ising antiferromagnet on the triangular and trillium lattice, emphasising the difference between the Husimi tree (HT) and Monte Carlo (MC) simulations. The global-axis triangular antiferromagnet possesses a non-monotonic reduced susceptibility $\chi T$, with a small but distinct minimum at $T = 0.9$.</figcaption></figure>

The excellent match above $T \gtrsim 1$ is in part due to the fact that the nearest-neighbour correlations in the degenerate ground state is $\langle\sigma_0\sigma_1\rangle = -1/3$ [23], for triangular and trillium systems in accordance with their corresponding Husimi tree [see Eq. (13)]. Indeed, the ground state energy is $E_{gs} = -N_{bond}J/3 = N_{bond}\langle\sigma_0\sigma_1\rangle$, where $N_{bond}$ is the number of nearest-neighbour bonds. For $T \lesssim 1$, correlations beyond nearest neighbours apparently start to play a role on the real lattices. From Fig. 8, the deviation from the HT curve indicates a dominant antiferromagnetic (resp. ferromagnetic) contribution for the trillium (resp. triangular) lattice [Eq. (8)]. In the triangular
case, the third nearest-neighbour correlations are known to be strongly ferromagnetic$^{23,24}$, with $\langle\sigma_0\sigma_3\rangle > |\langle\sigma_0\sigma_1\rangle|$, as $T \to 0^+$. It is likely that this increase of ferromagnetic correlations in the ground state causes an upturn of the reduced susceptibility [Fig. 8]. Accordingly, integrated correlations in the triangular Ising antiferromagnet are more antiferromagnetic at finite temperature, for $T \approx 0.9$, than in the spin-liquid ground state. Such a non-monotonic behaviour of the reduced susceptibility $\chi T$ is unusual, but not rare.

It is even more pronounced for the trillium lattice with easy axes. The reduced susceptibility $\chi T$ of easy-axes models necessarily increase upon cooling from high temperature, because nearest neighbor correlations are always ferromagnetic (the scalar product in Eq. (6) is always negative). For the trillium lattice, however, one can show that $\mathcal{C}_0 \approx \mathcal{C}_\infty = 1$ [see Appendix C 4 g]. It means that $\chi T$ has to decrease at low temperature.

The phenomenon of reentrance with bond-dependent interaction anisotropy is yet another example of non-monotonic $\chi T$, and discussed in detail elsewhere [13,14].

# IV. HUSIMI ANSATZ FOR THE CURIE-LAW CROSSOVER

## A. Limitation of the Curie-Weiss fit

As mentioned in the introduction, the Curie-Weiss temperature $\theta_{cw} = -zJ$ is a mean-field estimate of the transition temperature $T_c$ for a system with connectivity $z$, where the Curie-Weiss law is a consequence of critical scaling invariance with critical exponent $\gamma = 1$. Even though the concept of conventional order does not apply to spin liquids, $\theta_{cw}$ does represent a meaningful quantity, as a measure of interaction strength. The practical question is, how accurately can this quantity be measured in experiments ?

Best estimates can only be made at high temperatures, since $\theta_{cw}$ is the first-order correction to the Curie law

$$ \frac{1}{\chi} = \frac{T}{\mathcal{C}} \left[ 1 - \frac{\theta_{cw}}{T} + \mathcal{O}\left(\frac{1}{T^2}\right) \right]. \quad (18) $$

And here is the main issue with the Curie-Weiss temperature $\theta_{cw}$. In magnets, the high-temperature regime is frequently not accessible, since it is two or three orders of magnitude larger than the characteristic exchange coupling $J$. For example in magnets with 3d valence electrons, $J$ is often of the order of $\sim 100$ K and the high-temperature regime is inaccessible because it lies above the melting point of the crystal. On the other hand for magnets with 4f valence electrons, $J$ is much smaller, of the order of $\sim 1$ K. But 4f ions have a large single-ion degeneracy, lifted by the local crystal field. This crystal field introduces a single-ion anisotropy with an associated energy scale, which varies a lot from one material to another, but the lowest single-ion excitation is usually of the order of 10 - 100 K. The high-temperature

region is thus difficult to access because the nature of magnetic moments changes with temperature$^{12}$. We refer the reader to the useful tutorial written by Mugiraneza & Hallas [5] for a practical, step-by-step, application of the Curie-Weiss fit.

The susceptibility measures the evolution of the spin-spin correlations [Eq. (8)]. And the problem is that, as we have seen throughout this paper, this evolution from paramagnetism to spin liquid takes place over several orders of magnitude in temperatures. It is thus naturally best seen on a logarithmic scale. Applying the Curie-Weiss law, which is a linear fit, can be dangerous. What appears to be a reasonable temperature window on a linear scale might actually only measure a small evolution of the spin-spin correlations. The Curie-Weiss fit will always give a result of course, but the outcome will depend on the window of measurement [Fig. 1]. And if the high-temperature regime is not available, then it is not possible to check if the value is correct or not, causing a potentially (largely) inaccurate estimate of $\theta_{cw}$.

## B. The Husimi Ansatz

To measure the Curie-Weiss temperature in spin liquids, a complementary approach, relying on data points within an experimentally accessible temperature region, would be welcome.

While very high-temperatures are often physically not accessible, very low-temperatures are also not ideal. Irrespectively of the possible difficulty to thermalise the sample, perturbations beyond the spin-liquid Hamiltonian usually set a temperature scale $T^*$ below which the physics of the spin liquid is lost; the system may order or fall out-of-equilibrium. The most appropriate window in experiments is thus at intermediate temperatures, precisely where the crossover between the two Curie laws takes place. And while low- and high-temperature expansions are the least accurate in this regime, Section III A has shown that HT calculations are quantitatively reliable over the entire temperature region for corner-sharing lattices.

Appendix C 2 gives the analytical formula of the susceptibility for different geometries of the Husimi tree. We notice that the reduced susceptibility is always of the form

$$ \chi T|_{HT} = \frac{\sum_i \alpha_i e^{\kappa_i/T}}{\sum_i \alpha'_i e^{\kappa'_i/T}}. \quad (19) $$

This expression is sufficiently generic that it should be able to fit any form of $\chi T$. But as it is written, Eq. (19) is unpractical. Fortunately, it turns out that only a few terms are usually necessary. The simplest pertinent form of Eq. (19) is

$$ \chi T|_{HA} = \frac{1 + b_1 \exp[c_1/T]}{a + b_2 \exp[c_2/T]} . \quad (20) $$
<figure><img src="image_11.png" /><figcaption>FIG. 9. Empirical fit [Eq.(20)] of the reduced susceptibility $\chi T$ on the pyrochlore lattice [Fig.2(h)], as obtained from classical Monte Carlo simulations of $\mathcal{H}$ in Eq. (5) and $\mathcal{H}_{6\times2}$ in Eq. (26). (a) pyrochlore Ising, global axis and (b) pyrochlore Heisenberg antiferromagnet show dominant antiferromagnetic correlations $\chi T|_{T \to 0} = 0$, while (c) pyrochlore Ising, local axis (spin ice) and (d) the pyrochlore pseudo-Heisenberg antiferromagnet show dominating ferromagnetic correlations $\chi T|_{T \to 0} = 2$. The fitting windows are given by the coloured bars at the bottom of each figure. Examples for different fitting windows are shown for high-temperature (1st row) and low-temperature regions (2nd row). The last row shows for each model a minimal fitting window, which is sufficient to reproduce $\chi T$ over the full range of temperatures. Simulations were done for system sizes of linear dimensions $L = 16$, i.e. $N = 65536$ spins.</figcaption></figure>

We shall refer to Eq. (20) as the Husimi Ansatz. In this form the Curie constant and Curie-Weiss temperature can be directly extracted from the fitting parameters:

$$ \mathcal{E}_{\infty}^{\text{HA}} = \frac{1 + b_1}{a + b_2}, \qquad (21) $$

$$ \theta_{\text{cw}}^{\text{HA}} = \frac{b_1 c_1}{1 + b_1} - \frac{b_2 c_2}{a + b_2}. \qquad (22) $$

Eq. (20) will be our primary phenomenological Ansatz for the rest of this paper. Intuitively, we understand that the $c_1$ and $c_2$ parameters correspond to effective energy scales in the Boltzmann factor. However, two energy scales might be too minimal to describe the physics of some models, especially if different types of couplings are involved. This is why we will also consider an extended

Ansatz to fit $\chi T$

$$ \chi T|_{\text{HA}} = \frac{1 + b_1 \exp[c_1/T]}{a + b_2 \exp[c_2/T] + b_3 \exp[c_3/T]}, \qquad (23) $$

where the Curie constant and Curie-Weiss temperature become

$$ \mathcal{E}_{\infty}^{\text{HA}} = \frac{1 + b_1}{a + b_2 + b_3}, \qquad (24) $$

$$ \theta_{\text{cw}}^{\text{HA}} = \frac{b_1 c_1}{1 + b_1} - \frac{b_2 c_2 + b_3 c_3}{a + b_2 + b_3}. \qquad (25) $$

### C. Benchmark of the Husimi Ansatz

The purpose of this section is to benchmark the Husimi Ansatz of Eq. (20) in a controlled way on various model
Hamiltonians. In Fig. 9 we fit the Curie-law crossover with Eq. (20) for pyrochlore models with global-axis and local-axis Ising spins. In order to test the Ansatz on a general framework, beyond the Ising models used to build our Husimi-based intuition, we also consider continuous spins on the Heisenberg antiferromagnet (HAF)⁷³⁻⁷⁵, and pseudo-Heisenberg antiferromagnet (pHAF)⁷⁶⁻⁷⁹. The pHAF is defined on the XXZ model as follows:

$$ \mathcal{H}_{\text{XXZ}} = \sum_{\langle ij \rangle} [J_{zz}S_i^z S_j^z - J_{\pm}(S_i^{\pm} S_j^{-} + S_i^{-} S_j^{\pm})], \quad (26) $$

with $S_i^z$ along the local [111] easy-axis, as defined in Tab. IV, for parameters $J_{\pm}/J_{zz} = -0.5^{79}$. This model is thermodynamically equivalent to the HAF, but with different magnetic correlations, and thus a distinct evolution of the Curie-law crossover.

Fig. 9(a) and (b) show vanishing $\mathcal{C}_0 = 0$, induced by the zero-divergence constraint on the ground state manifold, imposing zero magnetisation in all tetrahedra (see discussion in Section III A). Fig. 9(c) and (d) show $\mathcal{C}_0 = 2$, as a result of dominant ferromagnetic correlations. Entering the spin-liquid regime at low $T$ for (a) and (c) for models with Ising degrees of freedom shows a rather sharp kink below $T/|J| \lesssim 1$, while on the opposite, models with continuous degrees of freedom in (b) and (d) enter the low-$T$ regime rather smoothly.

Results were obtained from classical MC simulations (black circles) and have been fitted with the Husimi Anstaz (solid lines) from Eq. (20) over different temperature windows. Examples of three different fitting windows are shown for high-temperature (1st row), and low-temperature (2nd row) fits. The range of fitting windows are indicated by blue, yellow and green bars on the bottom of each plot, and allow to judge their reliability in comparison to MC data. It becomes clear that fitting windows, which include only one Curie-law regime (either at low or high temperature), do not accurately reproduce the Curie-law crossover. This is especially important for Ising models, because of the relatively sharp kink when entering the spin-liquid regime.

On the other hand, fitting windows including the intermediate temperature window, with only the onset of high- and low-temperature regimes quantitatively reproduce $\chi T$ over the full temperature range. The 3rd row of panels shows the “minimal” fitting window. By using Eqs. (21, 22) we can precisely extract the Curie constant $\mathcal{C}_{\infty}$ and Curie-Weiss temperature $\theta_{cw}$ from those fits. Fitted and exact solutions match perfectly within error bars (see Tab. II).

This benchmark shows that the Husimi Ansatz correctly reproduces the Curie-law crossover over the full range of temperatures for several distinct models with Ising and continuous spins. It requires a fitting window spanning typically 1 or 2 orders of magnitude in temperature, in the intermediate regime that is usually accessible to experiments [see the bottom row of Fig. 9]. This is a useful theoretical proof of concept, that now needs to be applied to experiments.

<table><caption>TABLE II. Curie constant $\mathcal{C}_{\infty}$ [Eq. (21)] and Curie temperature $\theta_{cw}$ [Eq. (22)], obtained for the fit of the reduced susceptibility $\chi T$ with minimally sufficient fitting window, as shown in the 3rd row of Fig. 9.</caption><thead><tr><th rowspan="2">model</th><th colspan="2">C<sub>∞</sub></th><th colspan="2">θ<sub>cw</sub></th></tr><tr><th>fit</th><th>exact</th><th>fit</th><th>exact</th></tr></thead><tbody><tr><td>Ising global</td><td>1.00(1)</td><td>1.0</td><td>-6.0(1)</td><td>-6</td></tr><tr><td>HAF</td><td>1.00(2)</td><td>1.0</td><td>-2.03(5)</td><td>-2</td></tr><tr><td>Ising local</td><td>1.00(2)</td><td>1.0</td><td>2.0(1)</td><td>2</td></tr><tr><td>pHAF</td><td>1.00(1)</td><td>1.0</td><td>0.65(2)</td><td>2/3</td></tr></tbody></table>

# V. THE HUSIMI ANSTAZ IN EXPERIMENTS

## A. NaCaNi₂F₇

First, let us consider a material where the Ansatz gives similar results to the Curie-Weiss fit. To do so, let us consider one of the closest materials to the canonical HAF.

NaCaNi₂F₇ is a spin-1 pyrochlore material, well described by a weakly perturbed nearest-neighbour Heisenberg Hamiltonian⁸⁰,⁸¹. It shows spin freezing at $T_f \approx 3.6$ K, which has been assumed to originate from Na¹⁺/Ca²⁺ charge disorder, however, no long-range magnetic order has been observed³³.

In Fig. 10, we plot the magnetic susceptibility of NaCaNi₂F₇, extracted from Ref. [33], on a semilogarithmic scale for $\chi T$ and on a linear scale for $1/\chi$. The data points are well fitted by the Husimi Anstaz of Eq. (20) over the whole range of accessible temperatures.

We fit the Husimi Ansatz within physically relevant temperatures $\Delta T = [3.6$ K, ..., $300$ K], above the freezing transition up to the maximally available datapoints, and obtain a Curie-Weiss temperature $\theta_{cw}^{HA} = -122(1)$ K, and a Curie constant $\mathcal{C}_{\infty}^{HA} = 1.67(1)$ (emu K)/(Oe mol-Ni), which gives an effective magnetic moment of $\mu_{eff}^{HA} = 3.65(1)\mu_B/Ni$. All these quantities are in good agreement with a standard Curie-Weiss fit over a temperature window $\Delta T = [150$ K, ..., $300$ K], which reveals $\theta_{cw} = -129(1)$ K with $\mu_{eff} = 3.6(1)\mu_B/Ni$. This strongly suggest, as also qualitatively visible from the straight behaviour of $1/\chi$ in Fig. 10(b), that experimentally measured data points reach the high-temperature regime where a standard Curie-Weiss fit becomes a reliable estimate.

## B. KCu₆AlBiO₄(SO₄)₅Cl

$KCu₆AlBiO₄(SO₄)₅Cl$ is a promising $S = 1/2$ quantum spin liquid candidate on the distorted square-kagome lattice, as reported by M. Fujihala et al. in Ref. [34]. Measurement of specific heat and susceptibility did not
<figure><img src="image_12.png" /><figcaption>FIG. 10. Fit of the experimental susceptibility for the pyrochlore material NaCaNi₂F₇ (black circles) plotted on (a) a semi-logarithmic plot for $\chi T$ and (b) a linear scale for $1/\chi$. Fitting the Husimi Ansatz [Eq. (20)] within a temperature window $\Delta T = [3.6 \text{ K}, \dots, 300 \text{ K}]$ (red straight line), gives an estimate of $\theta_{cw} = -122(1) \text{ K}$. Our result is in good agreement with a standard Curie-Weiss fit (blue dashed line). Experimental data were extracted from Ref. [33].</figcaption></figure>

find any signatures of long-range order down to 1.8 K, while $\mu$SR confirmed the absence of spin order and spin freezing down to 58 mK.

In Fig. 11, we plot the magnetic susceptibility of KCu₆AlBiO₄(SO₄)₅Cl, as kindly provided by M. Fujihala [34], on a semi-logarithmic scale for $\chi T$ and on a linear scale for $1/\chi$. In comparison to NaCaNi₂F₇ in Fig. 10(b), it becomes evident that $1/\chi$ for KCu₆AlBiO₄(SO₄)₅Cl shows a rather strong deviation from a linear behaviour over nearly the whole range of experimentally accessible temperatures. A Curie-Weiss fit for the high-temperature tail within $\Delta T = [200 \text{ K}, \dots, 300 \text{ K}]$ gives $\theta_{cw} = -237(2) \text{ K}$ and $\mu_{eff} = 1.96 \mu_B/Cu$ with a Landé factor $g = 2.25^{34}$. The Husimi Ansatz from Eq. (20) gives a noticeably different outcome though. We find $\theta_{cw}^{HA} = -154 \pm 28 \text{ K}$. The large error bar comes from the choice of the fitting temperature window [$T_{min}, T_{max}$] (see the spread of the red curve in Fig. 11), where we fix $T_{max} = 300 \text{ K}$ at the highest available temperature, and vary $T_{min}$ between 10 and 30 K. The non-linearity of $1/\chi$ and spread of the Husimi estimate suggest that 300 K is too far from the high-temperature limit for a conclusive estimate of $\theta_{cw}$. The noticeable difference between the outcomes of the Curie-Weiss fit and Husimi Ansatz, however, makes us wonder which of the two estimates is

<figure><img src="image_13.png" /><figcaption>FIG. 11. Fit of the experimental susceptibility for the $S = 1/2$ square-kagome quantum spin liquid candidate KCu₆AlBiO₄(SO₄)₅Cl (black circles), plotted on (a) a semi-logarithmic plot for $\chi T$ and (b) a linear scale for $1/\chi$. Experimental data are plotted together with exact diagonalization (ED) data (18 sites) of an effective $J_1-J_2-J_3$ Heisenberg model (green dots), as proposed in [34]. We fit experimental data with the Husimi Ansatz [Eq. (20)] for $T_{min} \le T \le T_{max}$, where we set $T_{max} = 300 \text{ K}$ to the highest available temperature from experiment, and vary $T_{min}$ between 10 and 30 K (shaded respectively from light red to red). Experimental data and ED results were kindly provided by M. Fujihala [34].</figcaption></figure>

more reliable.

From a microscopic analysis in [34] we understand that KCu₆AlBiO₄(SO₄)₅Cl is not an ideal square-kagome lattice; the three bonds of a triangle in Fig. 2(c) are inequivalent. All triangles are distorted in the same way and form three distinct “nearest-neighbour” couplings, $J_1, J_2, J_3$, on each triangle. M. Fujihala et al.³⁴ built a microscopic Hamiltonian which describes its magnetic susceptibility at high temperature, using exact diagonalization (ED) and finite-temperature Lanczos methods, as shown on Fig. 11. ED results fit the experimental data very well down to $T \approx 40 \text{ K}$, below which finite-size effects make further estimates difficult. M. Fujihala et al. obtained

$$J_1 = -135 \text{ K}, J_2 = -162 \text{ K}, J_3 = -115 \text{ K}, \quad (27)$$

with a Landé factor $g = 2.11$. This high-temperature analysis cannot rule out low-energy perturbations, but it establishes the model of M. Fujihala et al as a solid parametrisation of KCu₆AlBiO₄(SO₄)₅Cl in the temperature regime relevant to $\theta_{cw}$, which is straightforward to
estimate from Eq. (27)

$$ \theta_{\text{cw}}^{\text{ED}} = \frac{S(S+1)}{3} 4 \left( \frac{J_1 + J_2 + J_3}{3} \right) \quad (28) $$

$$ = \left( \frac{J_1 + J_2 + J_3}{3} \right) \approx -137 \text{ K.} \quad (29) $$

Eq. (29) leads to a couple of remarks. Firstly, the ED results are in better agreement with the Husimi Ansatz than the Curie-Weiss fit, which a posteriori validates the former. Secondly, $\theta_{cw}$ here simply corresponds to the average value of the three inequivalent exchange couplings. $J_1, J_2, J_3$ fit within the energy window set by $\theta_{cw} \pm \delta J$, thus defining the anisotropic energy scale $\delta J = 25$ K. Using $T_{min} = \delta J$ as a lower bound of our fitting temperature window, we obtain from the Husimi Ansatz $\theta_{cw}^{HA} = -136$ K with a Landé factor $g = 2.1$, which is essentially the same result as from ED$^{82}$. This suggests that the main difficulty to estimate $\theta_{cw}$ comes from the lattice anisotropy of KCu$_6$AlBiO$_4$(SO$_4$)$_5$Cl. And while the Curie-Weiss law is not adapted to account for multiple energy scales in this intermediate regime, the Husimi Ansatz has been designed to be a flexible fitting function for the crossover that happens in this intermediate regime. We believe it is the reason why the Husimi Ansatz, albeit its large error bar, gives a better result than the Curie-Weiss fit.

### C. FeCl$_3$

#### 1. Experiments

As seen from the two previous materials with negative Curie-Weiss temperatures, spin liquids usually show dominant antiferromagnetic couplings. However, there also exist frustrated magnets where the interplay between ferro- and antiferro-magnetism can lead to multiple Curie-law crossovers$^{13,14}$. An important example relevant to materials are spiral spin liquids. They form a class of classical spin liquids where spiral states compete and form a sub-extensive ground state manifold with characteristic ring features in momentum space$^{83–87}$.

The Van der Waals magnet FeCl$_3$ is a prototype of a spiral spin liquid. At first, investigated as a member of crystallized anhydrous ferric chlorides$^{88}$, the history of FeCl$_3$ dates far back into the 1930's, where it already attracted attention due to its unusual magnetic properties at low temperature. Susceptibility measurements reported a Curie-Weiss temperature of $\theta_{cw} \approx -12$ K, however, noticing already at that time a deviation from the conventional Curie-Weiss law$^{89}$. Furthermore, inelastic neutron-scattering (INS) measurements$^{90}$, magnetic susceptibility$^{91}$, Möessbauer effect$^{92}$, magnetic field$^{93}$, and NMR measurements$^{94}$ confirmed a phase transition into an unusual spiral ground state at about $T_N \approx 10$ K. But it was only recently, with the work of S. Gao et al. [35], that continuous ring features around the $\Gamma$-point

could be observed in INS experiments; a clear evidence of spiral spin liquid physics in FeCl$_3$.

In Fig. 12 we show the magnetic susceptibility of FeCl$_3$, as kindly provided by M. McGuire [35], on a semi-logarithmic scale for $\chi T$ and on a linear scale for $1/\chi$. In contrast to the materials above (see Fig. 10 and Fig. 11), it seems that $\chi T$ reaches the plateau of the high-temperature Curie-Weiss regime already at about 300 K. For the traditional $1/\chi$ vs $T$ plot [Fig. 12(b)], the Curie-Weiss law shows a good fit over the temperature window $\Delta T = [100 \text{ K}, \dots, 350 \text{ K}]$, which gives $\theta_{cw} = -11(1)$ K, in agreement with previous measurements$^{89}$. However, when plotting the reduced susceptibility $\chi T$ instead [Fig. 12(a)], the Curie-Weiss law is seen to noticeably deviate from experimental data below 50 K. In fact, after careful consideration, experimental data do not seem to reach a plateau, but rather show a broad maximum at about $T \approx 300$ K, suggesting that the reduced susceptibility $\chi T$ is not monotonic. It means that FeCl$_3$ might be an experimental realisation of a multi-step Curie-law crossover. This motivates us to fit the available experimental data with the extended Husimi Ansatz of Eq. (23) which allows for non-monotonic behaviour. It fits the experimental data well over the whole temperature range and indeed presents a slight downturn at high temperatures above $T > 500$ K. Unfortunately, there are not enough data points after the downturn of $\chi T$ to extract a reliable estimate of $\theta_{cw}$; and since susceptibility measurements are naturally more noisy at high temperature, one should remain cautious. That being said, the Husimi Ansatz suggests a positive Curie-Weiss temperature in FeCl$_3$ – as opposed to previous measurements$^{88,89,91}$ – and thus a multi-step Curie-law crossover with dominant ferromagnetic interactions.

Interestingly, this outcome is consistent with INS experiments$^{35}$. The spiral spin liquidity of FeCl$_3$ comes from the first- and second-neighbor couplings. When they are dominantly ferromagnetic (resp. antiferromagnetic), they form ring features in the structure factor around the $\Gamma$ (resp. K) points. In FeCl$_3$, S. Gao et al measured clear, round, circles around the $\Gamma$ points in INS experiments$^{35}$. Further neighbor exchanges may ultimately invert the sign of $\theta_{cw}$, but INS measurements are consistent with a positive $\theta_{cw}$.

#### 2. Simulations

Since FeCl$_3$ is structurally unstable at higher temperatures, it is difficult to push the experimental analysis any further. Therefore, to conclude this discussion on the multi-step Curie-law crossover, we shall turn to classical Monte Carlo (MC) simulations. Magnetic Fe$^{3+}$ ($S = 5/2$) ions cover honeycomb layers, which are stacked in an ABC arrangement along the c axis. By comparing LSW theory and SCGA results to INS data, S. Gao et al. [35] proposed a series of models with up to nine coupling parameters. For the sake of simplicity, we will focus on
<figure><img src="image_14.png" /><figcaption>FIG. 12. Fit of the experimental susceptibility for the S = 5/2 magnet FeCl$_{3}$ (black circles) plotted on (a) a semi-logarithmic plot for χT and (b) a linear scale for 1/χ. Fitting the extended Husimi Ansatz [Eq. (23)] for all available data points (red straight line) reveals a multistep (non-monotonic) Curie-law crossover with a slight downturn at high temperatures, and hence a positive θ$^{HA}$. A standard Curie-Weiss fit (blue dashed line) gives a very different result of θ$_{cw}$ = −11(1) K, while showing a strong deviation from experimental data points below 50 K. Experimental data were kindly provided by M. McGuire [35].</figcaption></figure>

<figure><img src="image_15.png" /><figcaption>FIG. 13. Fit of the numerical susceptibility (black circles) plotted on (a) a semi-logarithmic plot for χT and (b) a linear scale for 1/χ. Numerical data were obtained from classical Monte Carlo simulations for the J$_{1}$-J$_{2}$-J$_{c1}$ Heisenberg model [Eq. (30)] on the ABC stacked honeycomb lattice for model parameters as proposed for FeCl$_{3}$ 35 [Eq. (31)]. Fitting the extended Husimi Ansatz [Eq. (23)] for all available data points (red straight line) clearly reveals a multistep (non-monotonic) Curie-law crossover with a significant downturn at high temperatures, and hence a positive θ HA cw = 8.7(2) K.</figcaption></figure>

their minimal model, which is able to reproduce the ring features of a spiral spin liquid; the J$_{1}$-J$_{2}$-J$_{c1}$ Heisenberg model.

$$
\mathcal{H}_{J_1 J_2 J_{c1}} = J_1 \sum_{\langle ij \rangle, \text{nn}}^{\text{intra}} \mathbf{S}_i \mathbf{S}_j + J_2 \sum_{\langle ij \rangle, \text{nnn}}^{\text{intra}} \mathbf{S}_i \mathbf{S}_j + J_{c1} \sum_{\langle ij \rangle, \text{nn}}^{\text{inter}} \mathbf{S}_i \mathbf{S}_j, \quad (30)
$$

where

$$
J_1 = -0.3 \text{ meV}, \quad J_2 = 0.075 \text{ meV}, \quad J_{c1} = 0.15 \text{ meV}(31)
$$

The couplings J$_{1}$ and J$_{2}$ respectively account for nearest-
neighbor and next nearest-neighbor interactions within
individual honeycomb layers, while J$_{c1}$ is the nearest-
neighbor antiferromagnetic interlayer coupling.

In Fig. 13 we show the susceptibility, measured from
MC simulations of $\mathcal{H}_{J_1 J_2 J_{c1}}$ [Eq. (30)]. Now the multistep
Curie-law crossover becomes evident and the extended
Husimi Ansatz from Eq. (23) gives $\theta_{cw}^{HA} = +8.7(2)$ K,
with a Curie constant $\mathcal{E}_{\infty}^{HA} = 3.6(1)$.

D. Summary about experimental comparison

In this section, we analyzed magnetic properties for three spin liquid candidates, namely, the *S* = 1 pyrochlore fluoride NaCaNi₂F₇, the *S* = 1/2 square-kagome material KCu₆AlBiO₄(SO₄)₅Cl, and the *S* = 5/2 spiral spin liquid on the honeycomb lattice FeCl₃. All three materials showed a Curie-law crossover over a wide temperature range, from ~ 1 K up to ~ 10⁴ K. Considering those examples, it becomes clear that a conventional Curie-Weiss fit applied to spin liquids can be reliable, but does not always have to. Depending on the microscopic model parameters and the nature of the underlying spin liquid, the high-temperature Curie-Weiss regime might not be practically accessible. We showed, that the comparison between the conventional Curie-Weiss fit and the Husimi Ansatz, as introduced in Sec. IV B, allows us to quantify whether the high-temperature regime of a material is reached or not.

NaCaNi₂F₇ is an example where experiments could
reach to the high-temperature regime, and results from
Husimi Ansatz and Curie-Weiss fit gave nearly the same
results. On the other hand, KCu₆AlBiO₄(SO₄)₅Cl,
shows a rather nonlinear behavior of 1/χ in Fig. 11(b)
for the available temperatures in experiment, which re-
sults in a mismatch between standard the Curie-Weiss fit and the Husimi Ansatz. The latter, however, agrees with independent ED results. Last but not least, FeCl₃ is probbaly an example of a multi-step Curie-law crossover. Such non-monotonic behaviour of magnetic correlations cannot be described by a conventional Curie-Weiss law, and therefore requires extra caution. By comparison to a minimal Heisenberg model we showed that an extended Husimi Ansatz [Eq. (23)] is able to capture such a non-monotonic Curie-law crossover, predicting a Curie-Weiss temperature which is noticeably different compared to the one obtained from a standard Curie-Weiss fit.

## VI. CONCLUSIONS

The Curie-Weiss temperature $\theta_{cw}$ is a useful quantity to estimate the strength of frustration in frustrated magnets [Eq. (2)]. However, the Curie-Weiss law is an estimate of the magnetic susceptibility close to a mean-field critical point, which – by definition – is absent in frustrated magnets. In this Article, we show that the concept of a Curie-law crossover¹⁶ is a generic feature for spin liquids and a more accurate description of their thermodynamic properties, that can be used to partially classify them. We systematically study the Curie-law crossover among a variety of frustrated Ising models in two and three dimensions [Fig. 2], and motivate its relevance to thermodynamic signatures, as seen in experiments. Comparing unbiased Monte Carlo simulations with the analytical Husimi-Tree approximation shows that the Curie-law crossover is determined by the type of frustrated unit cell (triangle, tetrahedron, ...) and the connectivity between them, rather than the physical dimension of the lattice. As a side note, the Husimi-Tree approximation proves to be quantitatively accurate for all temperatures and for many spin-liquid models [Fig. 5].

As a consequence of the Curie-law crossover, we recommend using the reduced susceptibility $\chi T$, complementary to the usual $1/\chi$ plot, when studying a potential spin liquid. It is often difficult to estimate whether $1/\chi$ has reached the asymptotic linear behavior, while $\chi T$ quickly indicates how far we are from the high-temperature Curie law.

Based on the success of the Husimi-Tree approximation, we propose an empirical Ansatz [Eq. (3)] as a useful complement to the Curie-Weiss law. The Husimi Ansatz is easy to use and designed to be a flexible fitting function for the crossover in $\chi T$ that takes place in the temperature regime which is typically accessible to experiments. This means that the Husimi Ansatz can be used on a broader temperature window than the Curie-Weiss fit, which is necessarily limited to the region where $1/\chi$ is linear in $T$. In its extended form [Eq. (23)], the Husimi Ansatz can also take into account the competition between ferro- and antiferromagnetic couplings in multi-step Curie-law crossovers.

It should be noted that the approach developed here

works for frustrated magnets in general. Frustration doesn't need to be geometric in origin, it may come from further neighbor or anisotropic spin exchange, as present in Kitaev materials. In addition, some materials might support the physics of a spin liquid at low but finite temperature, before ordering (or spin-freezing) at ultra low temperatures. In this case, the low-temperature kink of the reduced susceptibility $\chi T$ would indicate if the “hidden” spin-liquid regime has been reached before any potential phase transition. And even if the simulations and calculations were based on classical spins in this paper, the Husimi Ansatz can be applied to quantum materials in the crossover regime, as done in Section V.

## ACKNOWLEDGMENTS

The authors thank Harald Jeschke, Elsa Lhotel, Rodolphe Clérac, Claire Lhuillier for fruitful discussions, and Shang Gao, Jason Krizan and Yukitoshi Motome for critically reading the manuscript. This work was supported by the Theory of Quantum Matter Unit, OIST and “Quantum Liquid Crystals” JSPS KAKENHI Grant No. JP19H05825. L.D.C.J. acknowledges financial support from CNRS (PICS France-Japan MEFLS) and from the French "Agence Nationale de la Recherche" under Grant No. ANR-18-CE30-0011-01. Numerical calculations were carried out using HPC Facilities provided by OIST, and the Supercomputer Center of the Institute for Solid State Physics, the University of Tokyo.

## Appendix A: Definition of local easy axes

We provide positions and definitions for the local easy axes $\vec{e}_i$ of Ising spins [see Eq. (6)] for the kagome (Tab. III), pyrochlore (Tab. IV), hyperkagome (Tab. V) and trillium lattice (Tab. VI). Models with global and local axes are equivalent up to a rescaling in the exchange coupling $J$ given in each table caption.

<table><caption>TABLE III. The three sublattices of the kagome lattice with their local easy axes $\vec{e}_i$ and their positions. The corresponding lattice vectors are $\vec{a} = (1, 0)$, $\vec{b} = \frac{1}{2}(-1, \sqrt{3})$. The rescaling of exchange coupling between local and global axes is $J_{\text{local}}^{\text{kagome}} = -2J_{\text{global}}^{\text{kagome}} < 0$.</caption><thead><tr><th>site index <i>i</i></th><th><i>e</i><sub><i>i</i></sub></th><th>position</th></tr></thead><tbody><tr><td>1</td><td>(0, 1)</td><td>(<sup>1</sup>&frasl;<sub>2</sub>, 0)</td></tr><tr><td>2</td><td><sup>1</sup>&frasl;<sub>2</sub>(&minus;&radic;3, &minus;1)</td><td><sup>1</sup>&frasl;<sub>4</sub>(3, &radic;3)</td></tr><tr><td>3</td><td><sup>1</sup>&frasl;<sub>2</sub>(&radic;3, &minus;1)</td><td><sup>1</sup>&frasl;<sub>4</sub>(1, &radic;3)</td></tr></tbody></table>
<table><caption>TABLE IV. The 16 sublattices in the cubic unit cell of the pyrochlore lattice with their local easy axes $\vec{e}_i$ and their positions. The corresponding lattice vectors respect the cubic symmetry of the lattice $\vec{a} = (1, 0, 0)$, $\vec{b} = (0, 1, 0)$, $\vec{c} = (0, 0, 1)$. The rescaling of exchange coupling between local and global axes is $J_{\text{local}}^{\text{pyrochlore}} = -3J_{\text{global}}^{\text{pyrochlore}} < 0$.</caption><thead><tr><th>site index <i>i</i></th><th><b>S&#773;<sub><i>i</i></sub></b></th><th>position</th></tr></thead><tbody><tr><td>1</td><td><sup>1</sup>&frasl;<sub>&radic;3</sub>(+1, +1, +1)</td><td><sup>1</sup>&frasl;<sub>8</sub>(-3, -3, 1)</td></tr><tr><td>2</td><td><sup>1</sup>&frasl;<sub>&radic;3</sub>(-1, +1, -1)</td><td><sup>1</sup>&frasl;<sub>8</sub>(-1, -3, 3)</td></tr><tr><td>3</td><td><sup>1</sup>&frasl;<sub>&radic;3</sub>(+1, -1, -1)</td><td><sup>1</sup>&frasl;<sub>8</sub>(-3, -1, 3)</td></tr><tr><td>4</td><td><sup>1</sup>&frasl;<sub>&radic;3</sub>(-1, -1, +1)</td><td><sup>1</sup>&frasl;<sub>8</sub>(-1, -1, 1)</td></tr><tr><td>5</td><td><sup>1</sup>&frasl;<sub>&radic;3</sub>(+1, +1, +1)</td><td><sup>1</sup>&frasl;<sub>8</sub>(1, -3, -3)</td></tr><tr><td>6</td><td><sup>1</sup>&frasl;<sub>&radic;3</sub>(-1, +1, -1)</td><td><sup>1</sup>&frasl;<sub>8</sub>(3, -3, -1)</td></tr><tr><td>7</td><td><sup>1</sup>&frasl;<sub>&radic;3</sub>(+1, -1, -1)</td><td><sup>1</sup>&frasl;<sub>8</sub>(1, -1, -1)</td></tr><tr><td>8</td><td><sup>1</sup>&frasl;<sub>&radic;3</sub>(-1, -1, +1)</td><td><sup>1</sup>&frasl;<sub>8</sub>(3, -1, -3)</td></tr><tr><td>9</td><td><sup>1</sup>&frasl;<sub>&radic;3</sub>(+1, +1, +1)</td><td><sup>1</sup>&frasl;<sub>8</sub>(-3, 1, -3)</td></tr><tr><td>10</td><td><sup>1</sup>&frasl;<sub>&radic;3</sub>(-1, +1, -1)</td><td><sup>1</sup>&frasl;<sub>8</sub>(-1, 1, -1)</td></tr><tr><td>11</td><td><sup>1</sup>&frasl;<sub>&radic;3</sub>(+1, -1, -1)</td><td><sup>1</sup>&frasl;<sub>8</sub>(-3, 3, -1)</td></tr><tr><td>12</td><td><sup>1</sup>&frasl;<sub>&radic;3</sub>(-1, -1, +1)</td><td><sup>1</sup>&frasl;<sub>8</sub>(-1, 3, -3)</td></tr><tr><td>13</td><td><sup>1</sup>&frasl;<sub>&radic;3</sub>(+1, +1, +1)</td><td><sup>1</sup>&frasl;<sub>8</sub>(1, 1, 1)</td></tr><tr><td>14</td><td><sup>1</sup>&frasl;<sub>&radic;3</sub>(-1, +1, -1)</td><td><sup>1</sup>&frasl;<sub>8</sub>(3, 1, 3)</td></tr><tr><td>15</td><td><sup>1</sup>&frasl;<sub>&radic;3</sub>(+1, -1, -1)</td><td><sup>1</sup>&frasl;<sub>8</sub>(1, 3, 3)</td></tr><tr><td>16</td><td><sup>1</sup>&frasl;<sub>&radic;3</sub>(-1, -1, +1)</td><td><sup>1</sup>&frasl;<sub>8</sub>(3, 3, 1)</td></tr></tbody></table>

<table><caption>TABLE V. The twelve sublattices of the hyperkagome lattice with their local easy axes $\vec{e}_i$ and their positions. The corresponding lattice vectors respect the cubic symmetry of the lattice $\vec{a} = (1, 0, 0)$, $\vec{b} = (0, 1, 0)$, $\vec{c} = (0, 0, 1)$. The rescaling of exchange coupling between local and global axes is $J_{\text{local}}^{\text{hyperK}} = -3J_{\text{global}}^{\text{hyperK}} < 0$.</caption><thead><tr><th>site index <i>i</i></th><th><b>e&#773;<sub><i>i</i></sub></b></th><th>position</th></tr></thead><tbody><tr><td>1</td><td><sup>1</sup>&frasl;<sub>&radic;3</sub>(+1, +1, +1)</td><td><sup>1</sup>&frasl;<sub>8</sub>(-3, -3, 1)</td></tr><tr><td>2</td><td><sup>1</sup>&frasl;<sub>&radic;3</sub>(-1, +1, -1)</td><td><sup>1</sup>&frasl;<sub>8</sub>(-1, -3, 3)</td></tr><tr><td>3</td><td><sup>1</sup>&frasl;<sub>&radic;3</sub>(-1, -1, +1)</td><td><sup>1</sup>&frasl;<sub>8</sub>(-1, -1, 1)</td></tr><tr><td>4</td><td><sup>1</sup>&frasl;<sub>&radic;3</sub>(+1, +1, +1)</td><td><sup>1</sup>&frasl;<sub>8</sub>(1, -3, -3)</td></tr><tr><td>5</td><td><sup>1</sup>&frasl;<sub>&radic;3</sub>(+1, -1, -1)</td><td><sup>1</sup>&frasl;<sub>8</sub>(1, -1, -1)</td></tr><tr><td>6</td><td><sup>1</sup>&frasl;<sub>&radic;3</sub>(-1, -1, +1)</td><td><sup>1</sup>&frasl;<sub>8</sub>(3, -1, -3)</td></tr><tr><td>7</td><td><sup>1</sup>&frasl;<sub>&radic;3</sub>(+1, +1, +1)</td><td><sup>1</sup>&frasl;<sub>8</sub>(-3, 1, -3)</td></tr><tr><td>8</td><td><sup>1</sup>&frasl;<sub>&radic;3</sub>(-1, +1, -1)</td><td><sup>1</sup>&frasl;<sub>8</sub>(-1, 1, -1)</td></tr><tr><td>9</td><td><sup>1</sup>&frasl;<sub>&radic;3</sub>(+1, -1, -1)</td><td><sup>1</sup>&frasl;<sub>8</sub>(-3, 3, -1)</td></tr><tr><td>10</td><td><sup>1</sup>&frasl;<sub>&radic;3</sub>(-1, +1, -1)</td><td><sup>1</sup>&frasl;<sub>8</sub>(3, 1, 3)</td></tr><tr><td>11</td><td><sup>1</sup>&frasl;<sub>&radic;3</sub>(+1, -1, -1)</td><td><sup>1</sup>&frasl;<sub>8</sub>(1, 3, 3)</td></tr><tr><td>12</td><td><sup>1</sup>&frasl;<sub>&radic;3</sub>(-1, -1, +1)</td><td><sup>1</sup>&frasl;<sub>8</sub>(3, 3, 1)</td></tr></tbody></table>

<table><caption>TABLE VI. The four sublattices of the trillium lattice with their local easy axes $\vec{e}_i$ and their positions. The corresponding lattice vectors respect the cubic symmetry of the lattice $\vec{a} = (1, 0, 0)$, $\vec{b} = (0, 1, 0)$, $\vec{c} = (0, 0, 1)$. The rescaling of exchange coupling between local and global axes is $J_{\text{local}}^{\text{trillium}} = -3J_{\text{global}}^{\text{trillium}} < 0$. The explicit position of each site within the unit cell is given by the crystal parameter $u = 0, 138$ in order to compare to previous work$^{29,40,95}$.</caption><thead><tr><th>site index <i>i</i></th><th><b>e<sub><i>i</i></sub></b></th><th>position</th></tr></thead><tbody><tr><td>1</td><td><sup>1</sup>&frasl;<sub>&radic;3</sub>(+1, +1, +1)</td><td>(<i>u</i>, <i>u</i>, <i>u</i>)</td></tr><tr><td>2</td><td><sup>1</sup>&frasl;<sub>&radic;3</sub>(+1, -1, -1)</td><td>(<sup>1</sup>&frasl;<sub>2</sub> + <i>u</i>, <sup>1</sup>&frasl;<sub>2</sub> - <i>u</i>, 1 - <i>u</i>)</td></tr><tr><td>3</td><td><sup>1</sup>&frasl;<sub>&radic;3</sub>(-1, +1, -1)</td><td>(1 - <i>u</i>, <sup>1</sup>&frasl;<sub>2</sub> + <i>u</i>, <sup>1</sup>&frasl;<sub>2</sub> - <i>u</i>)</td></tr><tr><td>4</td><td><sup>1</sup>&frasl;<sub>&radic;3</sub>(-1, -1, +1)</td><td>(<sup>1</sup>&frasl;<sub>2</sub> - <i>u</i>, 1 - <i>u</i>, <sup>1</sup>&frasl;<sub>2</sub> + <i>u</i>)</td></tr></tbody></table>

## Appendix B: Monte Carlo simulations

Numerical Monte Carlo (MC) simulations of the Hamiltonian $\mathcal{H}$ [Eq. (5)] for Ising spins (Ising model) were performed by updating the site dependent Ising variable $\sigma = \pm 1$ for systems larger than $N = 10000$ spins. To account for statistically independent samples at very low temperatures a local single-spin flip Metropolis update algorithm as been used in combination with parallel tempering$^{96,97}$, and a worm-update algorithm$^{98,99}$ in the case of the checkerboard, pyrochlore and ruby lattice. A single MC step consists of $N$ local single spin-flip updates on randomly chosen sites, and 5 worm updates (checkerboard, pyrochlore and ruby lattice), performed in parallel for replicas at 100 to 200 different temperatures, with replica-exchange initiated by the parallel tempering algorithm every $10^2$ MC step.

MC simulations of the Hamiltonian $\mathcal{H}$ [Eq. (5)] for Heisenberg spins (Heisenberg model) were performed by using a local heat-bath algorithm$^{100,101}$, in combination with parallel tempering$^{96,97}$, and over-relaxation$^{102}$. Here, a single MC step consists of $N$ local heat-bath updates on randomly chosen sites, with $N$ over-relaxation steps, flipping the spin direction about their local exchange field, and replica-exchange every $10^2$ MC step.

In both cases, simulations for Ising and Heisenberg models, thermodynamic quantities were averaged over $10^6$ statistically independent samples, after $10^6$ MC steps for simulated annealing and $10^6$ MC steps for thermalization.

## Appendix C: Husimi Tree

### 1. Explicit calculations for the kagome Ising antiferromagnet

In this section, the Husimi tree calculation shall be explained on the example of HT(3,2) [see Figs.3(a) and 14(a)]. Branches of nonintersecting triangular cells
<figure><img src="image_16.png" /><figcaption>FIG. 14. The Husimi tree HT(3,2) for the kagome lattice. (a) Triangular frustrated cells arranged in shells, where shell 0 (red) corresponds to the central unit. The Husimi tree is equivalent to the real kagome lattice up to its 2nd shell (green). (b) All possible spin configurations for an isolated triangular cell with corresponding Boltzmann weights $g_0$ and $g_1$ for global axis Ising spins [see Eq. (C4)].</figcaption></figure>

spread out from the central unit (shell 0, drawn in red). Let us consider the Hamiltonian Eq. (5) for Ising spins $\sigma_i$ on sites $i$ with an additional external magnetic field $h$

$$ \mathcal{H} = J \sum_{\langle ij \rangle} \sigma_i \sigma_j - h \sum_i \sigma_i . \qquad (C1) $$

At the end of the calculations, the magnetic field will be taken vanishingly small in order to obtain the susceptibility $\chi$. The magnetisation on one of the central site (chosen arbitrarily) is

$$ \langle \sigma_1 \rangle = \frac{1}{Z_0} \sum_{\{\sigma_1, \sigma_2, \sigma_3\}} \sigma_1 \left( \prod_{\langle ij \rangle} g_{ij} \right) \left( \prod_{i=1}^{3} \alpha_i \right) \cdot Z_1(\sigma_1) Z_1(\sigma_2) Z_1(\sigma_3), \qquad (C2) $$

$$ \text{with } Z_0 = \sum_{\{\sigma_1, \sigma_2, \sigma_3\}} \left( \prod_{\langle ij \rangle} g_{ij} \right) \left( \prod_{i=1}^{3} \alpha_i \right) \cdot Z_1(\sigma_1) Z_1(\sigma_2) Z_1(\sigma_3), \qquad (C3) $$

being the total partition function. $\prod_{\langle ij \rangle}$ denotes the product over all nearest-neighbour pairs within the central triangular plaquette. $Z_1(\sigma_i)$ is the partition function of the branch of the Husimi tree moving outwards and starting from the central spin $i$ with orientation $\sigma_i$. Let us label $Z_n(\sigma_j)$ the partition function starting on site $j$ belonging to the $n^{th}$ layer of the tree. The Boltzmann weights are

$$ g_{ij} = e^{-\beta J \sigma_i \sigma_j}, \qquad (C4) $$

$$ \alpha_i = e^{\beta h \sigma_i}, \qquad (C5) $$

taking the values $g_0 = e^{\beta J}$, and $g_1 = e^{-3\beta J}$ here [Fig. 14(b)]. Eq. (C2) then gives explicitly

$$ \langle \sigma_1 \rangle = \frac{g_0(Y_1 - Y_1^2) + g_1(1 - Y_1^3)}{3g_0(Y_1 + Y_1^2) + g_1(1 + Y_1^3)}, \qquad (C6) $$

where we introduced the ratio between partition functions of a spin on shell $n$, pointing $\uparrow$ ($\sigma = 1$) and a spin pointing down $\downarrow$ ($\sigma = -1$)$^{52}$ as

$$ \frac{Z_n(\downarrow)}{Z_n(\uparrow)} = Y_n e^{2\beta h}, \qquad (C7) $$

and where successive layers of the Husimi tree are related recursively

$$ Y_n = \frac{g_0 Y_{n+1}^2 + g_1 (1 + 2Y_{n+1})}{g_0 + g_1 Y_{n+1} (2 + Y_{n+1})} \qquad (C8) $$

To solve the Husimi tree, we calculate the limit $Y_n \xrightarrow{n \to +\infty} Y$ and replace it in Eq. (C6), $Y_1 = Y^{103}$. In absence of an external magnetic field $Y = 1$, since the disordered system does not prefer any spin direction. This gives $\langle \sigma_1 \rangle = 0$ as trivially expected. But other observables such as the energy $E$, specific heat $C$ and entropy $S$ can be derived analytically from the partition function $Z_0$. These calculations are relatively straightforward and explicit results for the different Husimi trees are given in Appendix C 2.

In this section, we will further show the calculation of the susceptibility. An external magnetic field $h$ causes a perturbation $\epsilon$ away from the trivial value

$$ Y = 1 - \epsilon, \qquad (C9) $$

which can be used together with Eqs. (C7), (C8) to obtain $\epsilon$ in first order of $h$

$$ \epsilon = 2\beta h \frac{3g_0 + g_1}{5g_0 - g_1}. \qquad (C10) $$

The first-order expansion in $h$ is sufficient to compute the magnetic susceptibility, since higher-order terms vanish as $h \to 0$. Introducing Eqs. (C7)–(C10) into Eq. (C6) gives the temperature-dependent magnetisation

$$ \langle \sigma_1 \rangle = \beta h \frac{g_0 + 3g_1}{5g_0 - g_1}. \qquad (C11) $$
and the reduced susceptibility

$$ \chi T \equiv \frac{1}{\beta} \left. \frac{\partial \langle \sigma_1 \rangle}{\partial h} \right|_{h \to 0^+} = \frac{g_0 + 3g_1}{5g_0 - g_1}. \quad (\text{C12}) $$

## 2. Analytic Equations

Next to the magnetization and reduced susceptibility [Eq. (C12)], thermodynamic observables like energy E, specific heat C and entropy S are directly obtained from the partition function of the Husimi tree [Eq. (C3)]¹⁵.

$$ E = -\frac{1}{Z_0} \left. \left( \frac{\partial Z_0}{\partial \beta} \right) \right|_{\{h, \epsilon\} \to 0}, \quad (\text{C13}) $$

$$ C = -\beta^2 \left. \left( \frac{\partial E}{\partial \beta} \right) \right|_{\{h, \epsilon\} \to 0}, \quad (\text{C14}) $$

$$ S = \beta E + \log \frac{Z_0}{Z_1(\uparrow)} A_0, \quad (\text{C15}) $$

where $A_0$ is fitted such that $S|_{T \to \infty} = \log(2)$.

Here we show analytic expressions for thermodynamic observables, as obtained by HT calculations. All thermodynamic observables are plotted in Fig. 5, for $J = 1$, inducing antiferromagnetic correlations between spins.

Husimi tree HT(3,2) corresponding to the kagome and hyperkagome lattices:

$$ E = 2J \frac{-g_0 + g_1}{3g_0 + g_1}, \quad (\text{C16}) $$

$$ C = 32J^2 \beta^2 \frac{g_0 g_1}{(3g_0 + g_1)^2}, \quad (\text{C17}) $$

$$ S = 2J\beta \frac{-g_0 + g_1}{3g_0 + g_1} + \frac{2}{3} \log \left[ \frac{1}{\sqrt{2}}(3g_0 + g_1) \right], \quad (\text{C18}) $$

$$ \chi_{\text{glob}} T = \frac{g_0 + 3g_1}{5g_0 - g_1}, \quad (\text{C19}) $$

where $g_0 = e^{\beta J}$, and $g_1 = e^{-3\beta J}$.

Husimi tree HTS corresponding to the square-kagome lattice:

$$ E = \frac{2}{3} J \frac{-41g_0 + 30g_2 + 8g_3 + 3g_4}{41g_0 + 52g_1 + 30g_2 + 4g_3 + g_4}, \quad (\text{C20}) $$

$$ C = \frac{16}{3} J^2 \beta^2 \frac{41g_0(26g_1 + 60g_2 + 18g_3 + 11g_4) + 30g_2(2g_3 - g_4 + 26g_1) + 26g_1(16g_3 + 9g_4) - g_4(g_4 + 22g_3)}{(41g_0 + 52g_1 + 30g_2 + 4g_3 + g_4)^2} \quad (\text{C21}) $$

$$ S = \frac{2}{3} J \beta \frac{-41g_0 + 30g_2 + 8g_3 + 3g_4}{41g_0 + 52g_1 + 30g_2 + 4g_3 + g_4} + \frac{1}{6} \log \left[ \frac{1}{2}(41g_0 + 52g_1 + 30g_2 + 4g_3 + g_4) \right] \quad (\text{C22}) $$

$$ \chi T = \frac{2}{3} \frac{2g_0 + 7g_1 + 15g_2 + 5g_3 + 3g_4}{17g_0 + 30g_1 + 16g_2 + 2g_3 - g_4} + \frac{4}{3} \frac{(65g_0 + 381g_1 + 605g_2 + 601g_3) + g_4(275g_1 + 103g_2 + 15g_3 + 3g_4)}{(41g_0 + 52g_1 + 30g_2 + 4g_3 + g_4)(17g_1 + 13g_2 + 3g_3 - g_4)} \quad (\text{C23}) $$

where $g_0 = e^{4\beta J}$, $g_1 = 1$, $g_2 = e^{-4\beta J}$, $g_3 = e^{-8\beta J}$, and $g_4 = e^{-12\beta J}$.

Husimi tree HT(3,3) corresponding to the triangular and trillium lattice:

$$ E = 3J \frac{-g_0 + g_1}{3g_0 + g_1}, \quad (\text{C24}) $$

$$ C = 48J^2 \beta^2 \frac{g_0 g_1}{(3g_0 + g_1)^2}, \quad (\text{C25}) $$

$$ S = 3J\beta \frac{-g_0 + g_1}{3g_0 + g_1} + \log \left[ \frac{1}{2}(3g_0 + g_1) \right], \quad (\text{C26}) $$

$$ \chi T = \frac{g_0 + 3g_1}{7g_0 - 3g_1}, \quad (\text{C27}) $$

where $g_0 = e^{\beta J}$, and $g_1 = e^{-3\beta J}$.

Husimi tree HT(4,2) corresponding to the checker-

board, ruby and pyrochlore lattice:

$$ E = -3J \frac{g_0 - g_2}{3g_0 + 4g_1 + g_2}, \quad (\text{C28}) $$

$$ C = 24J^2 \beta^2 \frac{g_0 g_1 + 4g_0 g_2 + 3g_1 g_2}{(3g_0 + 4g_1 + g_2)^2}, \quad (\text{C29}) $$

$$ S = -3J \beta \frac{g_0 - g_2}{3g_0 + 4g_1 + g_2} + \frac{1}{2} \log \left[ \frac{1}{2}(3g_0 + 4g_1 + g_2) \right], \quad (\text{C30}) $$

$$ \chi T = 2 \frac{g_1 + g_2}{3g_0 + 2g_1 - g_2}, \quad (\text{C31}) $$

where $g_0 = e^{2\beta J}$, $g_1 = 1$ and $g_2 = e^{-6\beta J}$.
### 3. High-temperature expansion of the susceptibility

As discussed in Section IV A, $\theta_{cw}$ contributes to the first order correction of the Curie law:

$$ \frac{1}{\chi} = \frac{T}{\mathcal{C}} \left[ 1 - \frac{\theta_{cw}}{T} (1 + \Delta(T)) \right] . \quad (C32) $$

The same high-temperature expansion can be done for the results from Husimi tree calculations, where second and higher-order terms will account for the deviation from the Curie-Weiss law. Curie-constant $\mathcal{C}$, Curie temperature $\theta_{cw}$ and higher-order corrections $\Delta(T)$, extracted from the inverse susceptibility $1/\chi$ for global axes spins are summarised as follows:

HT(3,2):

$$ \begin{aligned} \mathcal{C} &= 1, \\ \theta_{cw} &= -4J, \\ \Delta(T) &= \frac{J}{T} - \frac{J^2}{3T^2} - \frac{5J^3}{3T^3} + \dots. \end{aligned} \quad (C33) $$

HTS:

$$ \begin{aligned} \mathcal{C} &= 1, \\ \theta_{cw} &= -4J, \\ \Delta(T) &= \frac{J}{T} - \frac{J^2}{3T^2} - \frac{4J^3}{3T^3} + \dots. \end{aligned} \quad (C34) $$

HT(3,3):

$$ \begin{aligned} \mathcal{C} &= 1, \\ \theta_{cw} &= -6J, \\ \Delta(T) &= \frac{J}{T} - \frac{J^2}{3T^2} - \frac{5J^3}{3T^3} + \dots. \end{aligned} \quad (C35) $$

HT(4,2):

$$ \begin{aligned} \mathcal{C} &= 1, \\ \theta_{cw} &= -6J, \\ \Delta(T) &= \frac{J}{T} - \frac{4J^2}{3T^2} - \frac{5J^3}{3T^3} + \dots. \end{aligned} \quad (C36) $$

Since $J=1$, all models show negative values for $\theta_{cw}$, indicating dominating antiferromagnetic interactions. Furthermore, their absolute values correspond to the number of nearest neighbor sites, and measures the local exchange field (Weiss field or molecular field) acting on every individual spin. The deviation $\Delta(T)$ of $\theta_{cw}$ scales independently of the type of the Husimi tree with $1/T$ in leading order. However, the deviation in second-order terms of $1/T^2$ differs between Husimi trees, made of triangular plaquettes and square plaquettes. And from this comparison it becomes evident that HTS shows only a small difference of 2% compared to HT(3,2) [see Tab. I], since their deviation happens from third-order $1/T^3$.

### 4. An alternative way to compute $\mathcal{C}_0$

In Appendix C1, the susceptibility was calculated as the linear response to an external magnetic field $h$, when $h \to 0$. At zero temperature it is also possible to calculate it as the sum of spin-spin correlations, following Eq. (8). When applied to the ground-state ensemble, this method allows to extract the value of the spin-liquid Curie constant $\mathcal{C}_0$ as has been done for spin-ice related models$^{15,16,19,104,105}$. For ease of calculations, let us consider that the Husimi tree is made of $L$ layers of spins, centred around a central site instead of a central frustrated unit. It is then common practice to consider this central spin $\vec{S}_0$ as the spin representative of the bulk of the real lattice. This is because $\vec{S}_0$ is the furthest away from the boundary of the tree, and thus less sensitive to surface effects. For a HT of $L$ layers, the spin-liquid Curie constant of Eq. (10) becomes

$$ \mathcal{C}_0(L) = 1 + \sum_{\ell=1}^{L} g_{\ell} \langle \vec{S}_0 \cdot \vec{S}_{\ell} \rangle, \quad (C37) $$

where $\langle \vec{S}_0 \cdot \vec{S}_\ell \rangle$ is the correlation between the central spin and one of the spins on layer $\ell \in [1: L]$, in the ground state. $g_\ell$ is the number of spins in this layer.

#### a. Kagome-type Husimi tree with global axis

For HT(3,2), the number of sites per layer is $g_\ell = 2 \times 2^\ell$. Using Eq. (C37) and Eq. (14) with global Ising axis, one gets

$$ \begin{aligned} \mathcal{C}_0(L) &= 1 + \sum_{\ell=1}^{L} 2 \cdot 2^\ell \left(-\frac{1}{3}\right)^\ell \\ &= 1 + 2 \left( \frac{1 - (-2/3)^{L+1}}{1 + 2/3} - 1 \right) \\ &= 0.2 - \frac{6}{5} \left(-\frac{2}{3}\right)^{L+1} \underset{L \to \infty}{\longrightarrow} 0.2. \end{aligned} \quad (C38) $$

The value of 0.2 is recovered in the thermodynamic limit of the alternating (antiferromagnetic) series of spin-spin correlations.

#### b. Trillium-type Husimi tree with global axis

For HT(3,3), the number of sites per layer is $g_\ell = (3/2) \times 4^\ell$. As a consequence, the series of Eq. (C37)
becomes alternating divergent, because of the boundary

$$
\begin{align}
\mathcal{C}_0(L) &= 1 + \sum_{\ell=1}^{L} \frac{3}{2} \cdot 4^\ell \left(-\frac{1}{3}\right)^\ell \nonumber \\
&= 1 + \frac{3}{2} \left( \frac{1 - (-4/3)^{L+1}}{1 + 4/3} - 1 \right) \nonumber \\
&= \frac{1}{7} - \frac{9}{14} \left(-\frac{4}{3}\right)^{L+1} . \tag{C39}
\end{align}
$$

If the size of the boundary grows faster than the cor-
lations decay, then the series diverges. That being said,
even if the calculation is mathematically ill posed, it is
interesting to notice that the constant term, 1/7, is the
same as the one obtained from the complete Husimi-tree
calculation [see Eq. (C27) in the limit $\beta \to +\infty$ and Ta-
ble I].

c. Pyrochlore-type Husimi tree with global axis

For HT(4,2), the number of sites per layer is $g_\ell = 2 \times 3^\ell$. As a consequence, the sum of Eq. (C37) becomes alternating,

$$
\begin{equation}
\begin{split}
\mathcal{C}_0(L) = 1 &+ \sum_{\ell=1}^{L} 2 \cdot 3^\ell \left(-\frac{1}{3}\right)^\ell \\
& \xrightarrow[L \to \infty]{} 1 + 2 \left(\frac{1}{1+1} - 1\right) = 0
\end{split}
\tag{C40}
\end{equation}
$$

As was the case for HT(3,3), even if the calculation is
mathematically ill posed, it is interesting to notice that
the outcome is the exact result [Table I].

d. Kagome Husimi tree with local easy axes

Considering local axes makes the calculation a bit more
complex, because spins are not collinear anymore. For
the kagome lattice, the local easy axes are given in Table
III, giving $\vec{e}_i \cdot \vec{e}_j = -1/2$ for spins on different sublattices.
Eq. (C37) then becomes

$$
\mathcal{C}_0(L) = 1 + \sum_{\ell=1}^{L} u_{\ell} \left(-\frac{1}{3}\right)^{\ell} + \sum_{\ell=1}^{L} v_{\ell} \left(-\frac{1}{2}\right)^{\ell} \left(-\frac{1}{3}\right)^{\ell} . \quad (C41)
$$

From now on, $u_\ell$ (resp. $v_\ell$) are the number of spins on
layer $\ell$ belonging to the same (resp. a different) sublattice
as the central spin of reference, $\vec{S}_0$. By definition, we have
$u_\ell + v_\ell = g_\ell = 2 \times 2^\ell$ for HT(3,2). It is not difficult to
see that these sequences are related by recursion

$$
\begin{align}
u_{\ell+1} &= v_{\ell} \notag \\
v_{\ell+1} &= v_{\ell} + 2 u_{\ell} \tag{C42}
\end{align}
$$

which gives

$$
\begin{align}
u_{\ell} &= \frac{2}{3} 2^{\ell} + \frac{4}{3} (-1)^{\ell} \nonumber \\
v_{\ell} &= \frac{4}{3} 2^{\ell} - \frac{4}{3} (-1)^{\ell} \tag{C43}
\end{align}
$$

Injecting Eq. (C43) into Eq. (C41), and taking the limit
$L \to +\infty$, finally gives $\mathcal{C}_0 = 2$ for the kagome lattice with
local easy axes.

e. Spin-ice Husimi tree with local easy axes

For 3D spin ice on the pyrochlore lattice [Table IV], the calculation is very similar. The scalar product between spins on different sublattices is now $\vec{e}_i \cdot \vec{e}_j = -1/3$, and the number of spins belonging to the same, $u_\ell$, and different, $v_\ell$, sublattices are

$$
\begin{align}
u_{\ell} &= \frac{1}{2} 3^{\ell} + \frac{3}{2} (-1)^{\ell} \nonumber \\
v_{\ell} &= \frac{3}{2} 3^{\ell} - \frac{3}{2} (-1)^{\ell} \tag{C44}
\end{align}
$$

which gives $\mathcal{C}_0 = 2$ for the pyrochlore lattice with local easy axes. Please note this is the same value, up to a normalisation, as the one calculated for the dielectric constant of cubic ice$^{104,105}$.

f. Hyperkagome Husimi tree with local easy axes

There are four different types of spin orientations in
the hyperkagome lattice [see Table V], labelled 1, 2, 3,
4. Let us assume that the central spin of reference has
orientation 1, at no cost in generality. When posing the
problem, one quickly sees that the number of spins with
orientation 1 in layer $\ell$ is not obvious to calculate, be-
cause there are four types of triangles, with orientations
{1, 2, 3}, {1, 2, 4}, {1, 3, 4}, {2, 3, 4}. Among the $u_\ell$ sites
with orientation 1 on layer $\ell$, we need to make a distinc-
tion between:

* the $a_\ell$ spins that have a site with orientation 1 as second neighbor in the internal layers ($n < \ell$),

* the $b_\ell$ spins that do not have a site with orientation 1 as second neighbor in the internal layers.

We have $u_\ell = a_\ell + b_\ell$ and $u_\ell + v_\ell = g_\ell = 2 \times 2^\ell$ sites on layer $\ell$ for HT(3,2). If we impose the local geometry of the hyperkagome lattice on HT(3,2), one gets the following recursion relations

$$
\begin{align}
a_{\ell+2} &= b_{\ell+1} + a_{\ell} + b_{\ell} \nonumber \\
b_{\ell+3} &= 2(a_{\ell+1} + a_{\ell} + b_{\ell}) \tag{C45}
\end{align}
$$

Imposing the appropriate initial conditions, one gets
$$a_{\ell} = 2^{\ell-2} + (-1)^{\ell} + 2^{\ell/2-2} \left( \frac{3}{\sqrt{7}} \sin \left( \ell \left( \pi - \tan^{-1} \left( \sqrt{7} \right) \right) \right) - \cos \left( \ell \left( \pi - \tan^{-1} \left( \sqrt{7} \right) \right) \right) \right) \quad (C46)$$

$$b_{\ell} = 2^{\ell-2} - 2^{\ell/2-2} \left( \frac{1}{\sqrt{7}} \sin \left( \ell \left( \pi - \tan^{-1} \left( \sqrt{7} \right) \right) \right) - 3 \cos \left( \ell \left( \pi - \tan^{-1} \left( \sqrt{7} \right) \right) \right) \right) \quad (C47)$$

whose sum can be simplified into

$$u_{\ell} = 2^{\ell-1} + (-1)^{\ell} + \frac{(-1)^{\ell+1} 2^{\ell/2+\frac{1}{2}} \sin \left( (\ell-1) \tan^{-1} \left( \sqrt{7} \right) \right)}{\sqrt{7}} \quad (C48)$$

Since the easy axes of the hyperkagome lattice give $\vec{e}_i \cdot \vec{e}_j = -1/3$ for spins with different orientations, we get

$$\begin{aligned}
\mathcal{C}_0 &= 1 + \sum_{\ell=1}^{+\infty} u_{\ell} \left(-\frac{1}{3}\right)^{\ell} + \sum_{\ell=1}^{+\infty} v_{\ell} \left(-\frac{1}{3}\right) \left(-\frac{1}{3}\right)^{\ell} \\
&= 1 + \frac{2}{3} \sum_{\ell=1}^{+\infty} \left(-\frac{1}{3}\right)^{\ell} \left(2u_{\ell} - 2^{\ell}\right) \\
&= \frac{5}{3} - \frac{4\sqrt{2}}{3\sqrt{7}} \sum_{\ell=1}^{+\infty} \left(\frac{\sqrt{2}}{3}\right)^{\ell} \sin\left((\ell-1)\tan^{-1}\left(\sqrt{7}\right)\right) \\
&= \frac{5}{3} - \frac{4\sqrt{2}}{3\sqrt{7}} \times \frac{1}{8}\sqrt{\frac{7}{2}} = \frac{3}{2} \quad (C49)
\end{aligned}$$

for the hyperkagome lattice with local easy axes.

### g. Trillium Husimi tree with local easy axes

There are four sublattices in the minimal unit cell of the trillium lattice, labelled 1, 2, 3, 4. Let us assume that the central spin of reference is on sublattice 1, at no cost in generality. As for the hyperkagome case in Appendix

C4f, there are four types of triangles, with sublattices {1, 2, 3}, {1, 2, 4}, {1, 3, 4}, {2, 3, 4}. Among the $v_\ell$ sites that do not belong to sublattice 1 on layer $\ell$, we need to make a distinction between:

* the $c_\ell$ spins that were in a triangle with a sublattice-1 site in layer $\ell-1$,

* the $d_\ell$ spins that were not in a triangle with a sublattice-1 site in layer $\ell-1$.

We have $v_\ell = c_\ell + d_\ell$ and $u_\ell + v_\ell = g_\ell = \frac{3}{2} \times 4^\ell$ sites on layer $\ell$ for HT(3,3). If we impose the local geometry of the trillium lattice on HT(3,3), one gets the following recursion relations

$$\left\{
\begin{array}{l}
u_{\ell+1} = c_{\ell} + 2d_{\ell} \\
c_{\ell+1} = 4u_{\ell} + c_{\ell} + 2d_{\ell} \\
d_{\ell+1} = 2c_{\ell}
\end{array}
\right.
\quad (C50)$$

which gives a self-consistent recursion relation for the number of sites in sublattice 1

$$u_{\ell+3} = u_{\ell+2} + 8 u_{\ell+1} + 16 u_{\ell} \quad (C51)$$

whose solution is

$$u_{\ell} = 3 \left( 2^{2\ell-3} + 2^{\ell-3} \left( \frac{1}{\sqrt{7}} \sin \left( \ell \left( \pi - \tan^{-1} \left( \frac{\sqrt{7}}{3} \right) \right) \right) \right) + 3 \cos \left( \ell \left( \pi - \tan^{-1} \left( \frac{\sqrt{7}}{3} \right) \right) \right) \right) \quad (C52)$$

Since the easy axes of the trillium lattice give $\vec{e}_i \cdot \vec{e}_j = -1/3$ for spins on different sublattices, we get

$$\begin{aligned}
\mathcal{C}_0(L) &= 1 + \sum_{\ell=1}^{L} \left(-\frac{1}{3}\right)^{\ell} \left(u_{\ell} - \frac{v_{\ell}}{3}\right) = 1 + \sum_{\ell=1}^{+\infty} \left(-\frac{1}{3}\right)^{\ell} \left(\frac{4}{3}u_{\ell} - \frac{4^{\ell}}{2}\right) \\
&= 1 + \frac{1}{2} \sum_{\ell=1}^{L} \left(-\frac{2}{3}\right)^{\ell} \left(\frac{1}{\sqrt{7}} \sin\left(\ell \left(\pi - \tan^{-1}\left(\frac{\sqrt{7}}{3}\right)\right)\right) + 3 \cos\left(\ell \left(\pi - \tan^{-1}\left(\frac{\sqrt{7}}{3}\right)\right)\right)\right) \quad (C53)
\end{aligned}$$

The sum of Eq. (C53) converges to zero for $L \to +\infty$, which is why the Husimi tree for the trillium lattice with easy axes gives $\mathcal{C}_0 = 1$ [see Table I]. However, the first

term of the sum is positive (it is 2/3 for $L=1$), which means that the build up of correlations at short distance is primarily ferromagnetic. This is consistent with the
increase of the reduced susceptibility $\chi T$ when cooling from high temperature in Fig. 5(b).

## 5. Comment on the Pauling entropy

For ice problems, the Pauling entropy provides a lower bound on the exact value of the entropy$^{50}$. Ice problems are defined as systems of connected vertices, where each link between two vertices has a direction (the spin), and each vertex possesses as many inward as outward links – the so-called ice rules. The ground state of the checkerboard and pyrochlore lattices are ice problems, and their Pauling entropy are indeed lower than their exact residual entropy [Table I]. The ground state of the ruby lattice is, however, not an ice problem$^{28}$, even if it is also made of corner-sharing tetrahedra with two spins up and two spins down. This is because the centre of the tetrahedra – the above-mentioned vertices – form a kagome lattice, which is not bipartite but tripartite. There are three kinds of tetrahedra on the ruby lattice, labeled for convenience red, green and blue. If an up spin is mapped to an outward (inward) link in a red (green) tetrahedron, what happens in the blue tetrahedra? It is easy to show that all-in/all-out states then appear in the blue tetrahedra, and the ground-state ensemble is thus not an ice problem. The ground state of the Ising ruby antiferromagnet is actually a $\mathbb{Z}_2$ spin liquid, as opposed to the U(1) gauge structure on pyrochlore$^{28}$. Nevertheless, despite these fundamental differences, the thermodynamic quantities of these three models (ruby, checkerboard and pyrochlore) are semi-quantitatively the same for all temperatures, including their residual entropy.

## Appendix D: $\mathcal{C}_0$ appearing in Coulomb gauge field theory

The spin-ice ground state is famously known as a U(1) Coulomb gauge field$^{27}$. This gauge-field texture comes from the ice rules (“2 in - 2 out”), that can be rewritten as a divergence-free constraint on the magnetisation field $\vec{M}(\vec{r})$ at position $\vec{r}$. At lowest order, the probability distribution of $\vec{M}(\vec{r})$ is$^{75}$

$$ \mathcal{P} \propto \exp \left( -\frac{\kappa_0}{2v_{\text{cell}}} \int d\vec{r}' |\vec{M}|^2(\vec{r}') \right), \quad (\text{D1}) $$

where $v_{\text{cell}}$ is the volume of the primitive unit cell. From Eq. (D1), the entropic stiffness $\kappa_0$ is also the inverse of the variance of the magnetisation in the spin-ice ground

state (up to a prefactor), i.e.

$$ \kappa_0 \propto \frac{1}{\mathcal{C}_0}, \quad (\text{D2}) $$

It means that $\mathcal{C}_0$ is a measure of the (inverse of) the strength of entropic interactions between topological-charge excitations in spin ice$^{106}$. To conclude, the stiffness is also the Lagrange multiplier appearing in the Self-Consistent Gaussian Approximation (SCGA) that ensures the spin-length constraint on average$^{67}$. For many models with continuous spins, this Lagrange multiplier can be computed analytically in the limit of zero and infinite temperatures, and thus offers an alternative way to compute the ratio $\mathcal{C}_0/\mathcal{C}_\infty$ and to connect it to the number of flat bands in the system (see discussion in Section III C).

## Appendix E: S(q) – equal-time structure factor

In this Appendix we present the equal-time (energy-integrated) structure factor,

$$ S(\vec{q}) = \sum_{\alpha, \beta=x,y,z} \left( \delta_{\alpha\beta} - \frac{q_{\alpha}q_{\beta}}{q^2} \right) S^{\alpha\beta}(\vec{q}), \quad (\text{E1}) $$

$$ S^{\alpha\beta}(\vec{q}) = \frac{1}{N} \sum_{i,j}^{N} e^{i\vec{q}(\vec{r}_i - \vec{r}_j)} \langle S_i^{\alpha} \cdot S_j^{\beta} \rangle, \quad (\text{E2}) $$

as observed in neutron scattering experiments where only correlations perpendicular to wave-vector $\vec{q}$ are measured.

Fig. 15 shows the structure factor for the three dimensional pyrochlore, hyperkagome, and trillium lattice, for Ising spins in their local easy axis plane [see definitions in Appendix A], as obtained from classical Monte Carlo simulations$^{29,30,107}$. Hereby, the pyrochlore lattice shows sharp pinch-point singularities at the Brillouin zone center characteristic of the two-in/two-out ice rules and emergent Coulomb field$^{27}$. Similar pinch points are also seen on the checkerboard lattice with global Ising spins (same ice problem projected onto 2D plane, Fig. 7). The hyperkagome lattice is a depleted pyrochlore lattice, with the same local easy axes but only 3 spins per tetrahedron$^{30,108-110}$. Its structure factor is reminiscent to the one of spin ice, but with broadened pinch points because the missing 4th spin per tetrahedron prevents the ice rules and subsequent emergent Coulomb field. The structure factor for the trillium lattice does not show any pinch-points.

<footnote n="1">¹ P. Curie, *Propriétés magnétiques des corps à diverses températures*, Ph.D. thesis, Faculty of Sciences, University of Paris, Paris: Gautheir-Villars (1895).</footnote>

<footnote n="2">² P. Weiss, *J. Phys. Theor. Appl.* **6**, 661 (1907).</footnote>

<footnote n="3">³ C. Kittel, *Introduction to Solid State Physics* (Wiley, 2004).</footnote>
<figure><img src="image_17.png" /><figcaption>FIG. 15. Equal-time structure factor $S(q)$ [Eq. (E1)], as observed from inelastic neutron scattering experiments of the Hamiltonian Eq. (5) for Ising spins along the local easy axis for the 3-dimensional pyrochlore (spin ice), hyperkagome and trillium lattice [see Fig 2]. Results were obtained from classical Monte Carlo simulations and are plotted for the [h,h,l] and [h,k,0] planes at $T/J = 0.01$, deep within the spin liquid regime. All lattices show diffuse scattering patterns, indicating the absence of magnetic order. (a)–(b) the pyrochlore lattice (spin ice) shows pinch-point singularities at the Brillouin zone center, a signature of the two-in/two-out ground state manifold (ice rules)$^{27,69}$. (c)–(d) the hyperkagome lattice, a depleted pyrochlore lattice, with only 3 spins per tetrahedron, shows similar features to spin ice, however with absence of pinch-points due to the violation of the ice rules. (e)–(f) the trillium lattice shows an absence of pinch-points, indicating a different nature of the spin liquid ground state, compared to spin ice.</figcaption></figure>

<footnote n="4">⁴ N. Ashcroft and N. Mermin, *Solid State Physics*, HRW international editions (Holt, Rinehart and Winston, 1976).</footnote>

<footnote n="5">⁵ S. Mugiraneza and A. M. Hallas, *Communications Physics 5*, 1 (2022).</footnote>

<footnote n="6">⁶ L. D. Landau, Zh. Eksp. Teor. Fiz. **7**, 19 (1937).</footnote>

<footnote n="7">⁷ A. P. Ramirez, *Annu. Rev. Mater. Sci.* **24**, 453 (1994).</footnote>

<footnote n="8">⁸ H. Nagasawa, *Physics Letters A* **25**, 475 (1967).</footnote>

<footnote n="9">⁹ A. F. J. Morgownik and J. A. Mydosh, *Phys. Rev. B* **24**, 5277 (1981).</footnote>

<footnote n="10">¹⁰ H. J. Silverstein, K. Fritsch, F. Flicker, A. M. Hallas, J. S. Gardner, Y. Qiu, G. Ehlers, A. T. Savici, Z. Yamani, K. A. Ross, B. D. Gaulin, M. J. P. Gingras, J. A. M. Paddison, K. Foyevtsova, R. Valenti, F. Hawthorne, C. R. Wiebe, and H. D. Zhou, *Phys. Rev. B* **89**, 054433 (2014).</footnote>

<footnote n="11">¹¹ M. A. de Vries, A. C. Mclaughlin, and J.-W. G. Bos, *Phys. Rev. Lett.* **104**, 177202 (2010).</footnote>

<footnote n="12">¹² Y. Li, S. M. Winter, D. A. S. Kaib, K. Riedl, and R. Valenti, *Phys. Rev. B* **103**, L220408 (2021).</footnote>

<footnote n="13">¹³ R. Pohle, O. Benton, and L. D. C. Jaubert, *Phys. Rev. B* **94**, 014429 (2016).</footnote>

<footnote n="14">¹⁴ B. Schmidt and P. Thalmeier, *Phys. Rev. B* **96**, 214443 (2017).</footnote>

<footnote n="15">¹⁵ L. D. C. Jaubert, *Topological Constraints and Defects in Spin Ice*, Ph.D. thesis, Ecole Normale Supérieure de Lyon (2009).</footnote>

<footnote n="16">¹⁶ L. D. C. Jaubert, M. J. Harris, T. Fennell, R. G. Melko, S. T. Bramwell, and P. C. W. Holdsworth, *Phys. Rev. X* **3**, 011014 (2013).</footnote>

<footnote n="17">¹⁷ S. V. Isakov, K. S. Raman, R. Moessner, and S. L. Sondhi, *Phys. Rev. B* **70**, 104418 (2004).</footnote>

<footnote n="18">¹⁸ M. Isoda, *Journal of Physics: Condensed Matter* **20**, 315202 (2008).</footnote>

<footnote n="19">¹⁹ A. J. Macdonald, P. C. W. Holdsworth, and R. G. Melko, *Journal of Physics-Condensed Matter* **23**, 164208 (2011).</footnote>

<footnote n="20">²⁰ A. Nag and S. Ray, *Journal of Magnetism and Magnetic Materials* **424**, 93 (2017).</footnote>

<footnote n="21">²¹ S. T. Bramwell, M. N. Field, M. J. Harris, and I. P. Parkin, *Journal of Physics: Condensed Matter* **12**, 483 (2000).</footnote>

<footnote n="22">²² T. T. A. Lummen, I. P. Handayani, M. C. Donker,</footnote>
D. Fausti, G. Dhalenne, P. Berthet, A. Revcolevschi, and P. H. M. van Loosdrecht, Phys. Rev. B **77**, 214310 (2008).

23 J. Stephenson, *Journal of Mathematical Physics* **5**, 1009 (1964).

24 E. Rastelli, A. Tassi, and L. Reatto, *Il Nuovo Cimento B* (1971-1996) **42**, 120 (1977).

25 E. Jurčišinová and M. Jurčišin, Physica A: Statistical Mechanics and its Applications **492**, 1798 (2018).

26 J. Richter, O. Derzhko, and J. Schnack, Physical Review B **105**, 144427 (2022).

27 C. L. Henley, *Annual Review of Condensed Matter Physics* **1**, 179 (2010).

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

29 T. E. Redpath and J. M. Hopkinson, Phys. Rev. B **82**, 014410 (2010).

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

31 L. Bovo, L. D. C. Jaubert, P. C. W. Holdsworth, and S. T. Bramwell, *Journal of Physics: Condensed Matter* **25**, 386002 (2013).

32 L. Bovo, M. Twengström, O. A. Petrenko, T. Fennell, M. J. P. Gingras, S. T. Bramwell, and P. Henelius, Nature Communications **9**, 1999 (2018), arXiv:1805.09332 [cond-mat.str-el].

33 J. W. Krizan and R. J. Cava, Phys. Rev. B **92**, 014406 (2015).

34 M. Fujihala, K. Morita, R. Mole, S. Mitsuda, T. Tohyama, S.-i. Yano, D. Yu, S. Sota, T. Kuwai, A. Koda, H. Okabe, H. Lee, S. Itoh, T. Hawai, T. Masuda, H. Sagayama, A. Matsuo, K. Kindo, S. Ohira-Kawamura, and K. Nakajima, Nature Communications **11**, 3429 (2020).

35 S. Gao, M. A. McGuire, Y. Liu, D. L. Abernathy, C. d. Cruz, M. Frontzek, M. B. Stone, and A. D. Christianson, Phys. Rev. Lett. **128**, 227201 (2022).

36 J. A. M. Paddison, H. S. Ong, J. O. Hamp, P. Mukherjee, X. Bai, M. G. Tucker, N. P. Butch, C. Castelnovo, M. Mourigal, and S. E. Dutton, Nature Communication **7**, 13842 (2016).

37 S. T. Bramwell and M. J. Harris, *Journal of Physics: Condensed Matter* **32**, 374010 (2020).

38 M. Udagawa and L. Jaubert, *Spin Ice*, Springer Series in Solid-State Sciences (Springer International Publishing, 2021).

39 D. Adroja, B. Padalia, S. Malik, R. Nagarajan, and R. Vijayaraghavan, *Journal of Magnetism and Magnetic Materials* **89**, 375 (1990).

40 J. M. Hopkinson and H.-Y. Kee, Phys. Rev. B **74**, 224441 (2006).

41 S. T. Bramwell and M. J. Harris, *Journal of Physics-Condensed Matter* **10**, L215 (1998).

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

43 A. S. Wills, R. Ballou, and C. Lacroix, Phys. Rev. B **66**, 144407 (2002).

44 K. Kanô and S. Naya, *Progress of Theoretical Physics* **10**, 158 (1953).

45 G. H. Wannier, Phys. Rev. **79**, 357 (1950).

46 G. H. Wannier, Phys. Rev. B **7**, 5017 (1973).

47 L. Pauling, *Journal of the American Chemical Society* **57**, 2680 (1935).

48 E. H. Lieb, *Physical Review* **162**, 162 (1967).

49 J. F. Nagle, *Journal of Mathematical Physics* **7**, 1484 (1966).

50 E. H. Lieb and F. Y. Wu, “Two-dimensional ferroelectric models,” in *Phase Transitions and Critical Phenomena*, Vol. 1, edited by C. Domb and M. S. Green (Academic, New York, 1972) pp. 331–490.

51 H. A. Bethe, *Proceedings of the Royal Society of London A: Mathematical, Physical and Engineering Sciences* **150**, 552 (1935).

52 R. J. Baxter, *Exactly Solved Models in Statistical Mechanics* (Dover books on physics, Courier Corporation, 1982).

53 K. Wada and Y. Ogawa, *Journal of the Physical Society of Japan* **67**, 112 (1998).

54 L. D. C. Jaubert, J. T. Chalker, P. C. W. Holdsworth, and R. Moessner, Phys. Rev. Lett. **100**, 067207 (2008).

55 L. D. C. Jaubert, J. T. Chalker, P. C. W. Holdsworth, and R. Moessner, Phys. Rev. Lett. **105**, 087201 (2010).

56 J. L. Monroe, *Physica A: Statistical Mechanics and its Applications* **256**, 217 (1998).

57 S. ichi Yoshida, K. Nemoto, and K. Wada, *Journal of the Physical Society of Japan* **71**, 948 (2002).

58 L. D. C. Jaubert, J. T. Chalker, P. C. W. Holdsworth, and R. Moessner, *Journal of Physics: Conference Series* **145**, 012024 (2009).

59 M. Udagawa, H. Ishizuka, and Y. Motome, Phys. Rev. Lett. **104**, 226405 (2010).

60 E. Jurčišinová and M. Jurčišin, Phys. Rev. E **89**, 032123 (2014).

61 E. Jurčišinová, M. Jurčišin, and A. Bobák, *Journal of Statistical Physics* **154**, 1096 (2014).

62 E. Jurčišinová, M. Jurčišin, and A. Bobák, *Physics Letters A* **378**, 1448 (2014).

63 L. Néel, *Comptes Rendus Hebdomadaires Des Seances De L Académie Des Sciences* **224**, 1488 (1947).

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

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

66 P. H. Conlon and J. T. Chalker, Phys. Rev. B **81** (2010), DOI 10.1103/PhysRevB.81.224413.

67 P. Conlon, *Aspects of Frustrated Magnetism*, Ph.D. thesis, University of Oxford (2010).

68 O. Benton, *Classical and Quantum Spin Liquids on the Pyrochlore lattice*, Ph.D. thesis, University of Bristol (2014).

69 R. W. Youngblood and J. D. Axe, Phys. Rev. B **23**, 232 (1981).

70 M. J. Harris, S. T. Bramwell, D. F. McMorrow, T. Zeiske, and K. W. Godfrey, Phys. Rev. Lett. **79**, 2554 (1997).

71 T. Fennell, P. P. Deen, A. R. Wildes, K. Schmalzl, D. Prabhakaran, A. T. Boothroyd, R. J. Aldus, D. F. McMorrow, and S. T. Bramwell, Science **326**, 415 (2009).

72 J. L. Monroe, *Physica A: Statistical Mechanics and its Applications* **256**, 217 (1998).

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

74 B. Canals and D. A. Garanin, *Canadian Journal of Physics* **79**, 1323 (2001).

75 C. L. Henley, Phys. Rev. B **71**, 014424 (2005).

76 K. A. Ross, L. Savary, B. D. Gaulin, and L. Balents, Phys. Rev. X **1**, 021002 (2011).

77 S. Onoda and Y. Tanaka, Phys. Rev. B **83**, 094411 (2011).

78 S. Lee, S. Onoda, and L. Balents, Phys. Rev. B **86**, 104412 (2012).

79 M. Taillefumier, O. Benton, H. Yan, L. D. C. Jaubert, and N. Shannon, Phys. Rev. X **7**, 041057 (2017).
80 K. W. Plumb, H. J. Changlani, A. Scheie, S. Zhang, J. W. Krizan, J. A. Rodriguez-Rivera, Y. Qiu, B. Winn, R. J. Cava, and C. L. Broholm, *Nature Physics* **15**, 54 (2019).

81 S. Zhang, H. J. Changlani, K. W. Plumb, O. Tchernyshyov, and R. Moessner, *Physical Review Letters* **122**, 167203 (2019).

82 It is worth noting that a Landé factor of 2.1 is a more realistic value for Cu(II) ions than the value of 2.25 found with the Curie-Weiss fit.

83 N. Niggemann, M. Hering, and J. Reuther, *Journal of Physics: Condensed Matter* **32**, 024001 (2019).

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

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

86 C.-J. Huang, J. Q. Liu, and G. Chen, *Phys. Rev. Res.* **4**, 013121 (2022).

87 H. Yan and J. Reuther, *Physical Review Research* **4**, 023175 (2022).

88 Lallemand, A., *Ann. Phys.* **11**, 97 (1935).

89 C. Starr, F. Bitter, and A. R. Kaufmann, *Phys. Rev.* **58**, 977 (1940).

90 J. W. Cable, M. K. Wilkinson, E. O. Wollan, and W. C. Koehler, *Phys. Rev.* **127**, 714 (1962).

91 E. R. Jones, O. B. Morton, L. Cathey, T. Auel, and E. L. Amma, *The Journal of Chemical Physics* **50**, 4755 (1969).

92 J. P. Stampfel, W. T. Oosterhuis, B. Window, and F. d. Barros, *Phys. Rev. B* **8**, 4371 (1973).

93 P. B. Johnson, S. A. Friedberg, and J. A. Rayne, *Journal of Applied Physics* **52**, 1932 (1981).

94 B. Kang, C. Kim, E. Jo, S. Kwon, and S. Lee, *Journal of Magnetism and Magnetic Materials* **360**, 1 (2014).

95 S. V. Isakov, J. M. Hopkinson, and H.-Y. Kee, *Phys. Rev. B* **78**, 014404 (2008).

96 R. H. Swendsen and J.-S. Wang, *Phys. Rev. Lett.* **57**, 2607 (1986).

97 D. J. Earl and M. W. Deem, *Phys. Chem. Chem. Phys.* **7**, 3910 (2005).

98 G. T. Barkema and M. E. J. Newman, *Physical Review* **E** **57**, 1155 (1998).

99 N. Prokof'ev and B. Svistunov, *Phys. Rev. Lett.* **87**, 160601 (2001).

100 J. A. Olive, A. P. Young, and D. Sherrington, *Phys. Rev. B* **34**, 6341 (1986).

101 Y. Miyatake, M. Yamamoto, J. J. Kim, M. Toyonaga, and O. Nagai, *Journal of Physics C: Solid State Physics* **19**, 2539 (1986).

102 M. Creutz, *Phys. Rev. D* **36**, 515 (1987).

103 Please note that we use this numbering of layers – increasing outwards – essentially for pedagogical reasons. Mathematically, the problem is well posed if we consider a change of variable: layer 0 would be on the outer boundary, infinitely far away from the central layer which would be labeled *n*. Following this notation, this is why the limit *Y* can be used in Eq. (C6).

104 W. Gobush and C. A. J. Hoeve, *The Journal of Chemical Physics* **57**, 3416 (1972).

105 A. Yanagawa and J. F. Nagle, *Chemical Physics* **43**, 329 (1979).

106 C. Castelnovo, R. Moessner, and S. L. Sondhi, *Phys. Rev. B* **84**, 144435 (2011).

107 S. T. Bramwell and M. J. P. Gingras, *Science* **294**, 1495 (2001).

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

109 C. L. Henley, *Journal of Physics: Conference Series* **145**, 012022 (2009).

110 A. Mizuno, Y. Shuku, M. M. Matsushita, M. Tsuchiizu, Y. Hara, N. Wada, Y. Shimizu, and K. Awaga, *Phys. Rev. Lett.* **119**, 057201 (2017).