# Reentrance of disorder in the anisotropic shuriken Ising model

Rico Pohle, Owen Benton, L. D. C. Jaubert

► To cite this version:

Rico Pohle, Owen Benton, L. D. C. Jaubert. Reentrance of disorder in the anisotropic shuriken Ising model. Physical Review B: Condensed Matter and Materials Physics (1998-2015), 2016, 94, pp.014429. 10.1103/PhysRevB.94.014429. hal-01542103

HAL Id: hal-01542103

https://hal.science/hal-01542103

Submitted on 19 Jun 2017

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

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

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

Distributed under a Creative Commons Attribution 4.0 International License
# Reentrance of disorder in the anisotropic shuriken Ising model

Rico Pohle,¹ Owen Benton,¹ and L.D.C. Jaubert¹

¹Okinawa Institute of Science and Technology Graduate University, Onna-son, Okinawa 904-0495, Japan

(Dated: February 22, 2016)

For a material to order upon cooling is common sense. What is more seldom is for disorder to reappear at lower temperature, which is known as reentrant behavior. Such resurgence of disorder has been observed in a variety of systems, ranging from Rochelle salts to nematic phases in liquid crystals. Frustration is often a key ingredient for reentrance mechanisms. Here we shall study a frustrated model, namely the anisotropic shuriken lattice, which offers a natural setting to explore an extension of the notion of reentrance between magnetic *disordered* phases. By tuning the anisotropy of the lattice, we open a window in the phase diagram where magnetic disorder prevails down to zero temperature. In this region, the competition between multiple disordered ground states gives rise to a double crossover where both the low- and high-temperature regimes are less correlated than the intervening classical spin liquid. This reentrance of disorder is characterized by an entropy plateau, a multi-step Curie law crossover and a rather complex diffuse scattering in the static structure factor. Those results are confirmed by complementary numerical and analytical methods: Monte Carlo simulations, Husimi-tree calculations and an exact decoration-iteration transformation.

PACS numbers: 75.10.Hk,75.30.Kz,75.10.Kt

Recent progress in frustrated magnetism has delivered entire maps of long-range ordered and disordered phases, obtained for example via the variation of bond anisotropy$^{1-7}$ or further nearest-neighbour couplings$^{8-14}$. Such phase diagrams have allowed to put a series of frustrated materials onto a global and connected map, that can be experimentally explored via physical or chemical pressure$^{15-18}$. On such phase diagrams, when two ordered phases meet, an enhancement of the classical ground-state degeneracy takes place$^{19}$. This degeneracy can either be lifted by thermal fluctuations, giving rise to multiple phase transitions$^{20,21}$, or may destroy any kind of order down to (theoretically) zero temperature. This is where spin liquids appear. But this picture is less clear at the frontier between ordered and (possibly multiple) disordered ground states. In particular how do disordered phases compete with each other at finite temperature ?

The frustrated shuriken lattice$^{22}$ – also known as square-kagome$^{23-30}$, squagome$^{31,32}$, squa-kagome$^{33}$ or L4-L8$^{33}$ lattice – provides an interesting model-example for such competition. Being made of corner-sharing triangles, it is locally similar to the famous kagome lattice, but with the important difference that the shuriken lattice is composed of two inequivalent sublattices [see Fig. 1]. Such asymmetry offers a natural setup for lattice anisotropy. In the asymptotic limits of this anisotropy, a promising zero-temperature phase diagram has emerged for quantum spin-1/2, ranging from a bipartite long-range ordered phase to a highly degenerate ground state made of tetramer clusters of spins$^{33}$. However, while the quantum ground states$^{22,26,30,33}$ and the influence of a magnetic field$^{22,24-29,34}$ have been studied to some extent, little is known about the finite-temperature properties in zero field$^{23,31}$.

In this paper, our goal is to develop a comprehensive and precise understanding of the frustrated phase

<figure><img src="image_2.png" /><figcaption>FIG. 1. **The shuriken lattice** as seen in real space (**left**) and Fourier space (**right**). There are 6 sites per unit-cell with two sublattices **A** and **B**. Interactions between A-sites (square plaquettes) are described with coupling constant $J_{AA}$ (red), while interactions between A- and B-sites (octagonal plaquettes) are described with $J_{AB}$ (black).</figcaption></figure>

diagram of the Ising model on the anisotropic shuriken lattice, relying on a combination of numerical and analytical methods (Monte Carlo simulations, Husimi tree calculations and decoration-iteration transformation). Using the lattice anisotropy as a tuning parameter, we find that this model supports two long-range ordered phases (ferromagnet (FM) and staggered ferromagnet (SFM)), two classical spin liquids (SL$_{1,2}$) characterized by complex static structure factors, and a *zero-temperature* paramagnet composed of two kinds of isolated (super)spins with strictly zero correlations between them. We shall refer to this latter phase as a *binary paramagnet*. Over an extended region of the phase diagram, there is a double crossover from the high-temperature paramagnet to the spin liquids and finally into the low-temperature binary paramagnet. This double crossover gives rise to a non-monotonic behavior of the correlation length, which can be seen as an analogue of reentrant behavior between disordered
<figure><img src="image_3.png" /><figcaption>FIG. 2. Phase Diagram of the Ising model on the anisotropic shuriken lattice. (a) The circles (triangles) correspond to phase transitions (crossovers), obtained by Monte Carlo simulations (Husimi-tree calculations) [see appendix A for further details]. As a function of the coupling ratio $x = J_{AB}/J_{AA}$, the model supports a long-range ordered ferromagnet (FM) [see panel (b)], a long-range ordered staggered ferromagnet (SFM) [see panel (c)], a binary paramagnet (BPM) [see panel (d)], and two classical spin liquids (SL1,2). The BPM illustrated in panel (d) is made of antiferromagnetically ordered square plaquettes, decoupled from each other and from the intermediate spins sitting on the B sublattice. For $|x| \gtrsim 1$, on cooling, the system undergoes a evolution from “gas crossover$\rightarrow$ liquid transition$\rightarrow$ solid”. As for $|x| \lesssim 1$, it provides a remarkable example of reentrance from “gas crossover$\rightarrow$ liquid crossover$\rightarrow$ gas”.</figcaption></figure>

phases. As a by-product, we notice an essentially perfect agreement between Husimi-tree analytics and Monte Carlo simulations in the disordered regimes.

The paper is divided as follows. The model is introduced in section I, followed by its phase diagram in section II. In section III, we analyze in details the double-crossover region between disordered regimes. We conclude the paper by discussing possible experimental realizations of the shuriken lattice and summarizing our results in sections IV and V respectively. Most technical details are given in the appendices.

## I. ANISOTROPIC SHURIKEN MODEL

The shuriken lattice is made of corner-sharing triangles with 6 sites per unit cell [see Fig. 1]. As opposed to its kagome parent where all spins belong to hexagonal loops, the shuriken lattice forms two kinds of loops made of either 4 or 8 sites. As a consequence, 2/3 of the spins in the system belong to the A-sublattice, while the remaining 1/3 of the spins form the B-sublattice. Let us respectively define $J_{AA}$ and $J_{AB}$ as the coupling constants between A-sites on the square plaquettes, and between A- and B-sites on the octagonal plaquettes. The Hamiltonian of the model can be written as:

$$H = -J_{AA} \sum_{\langle ij \rangle_{AA}} \sigma_i^A \sigma_j^A - J_{AB} \sum_{\langle ij \rangle_{AB}} \sigma_i^A \sigma_j^B \quad (1)$$

where we consider Ising spins $\sigma_i = \pm 1$ with nearest-neighbor coupling.

There is no frustration for ferromagnetic $J_{AA} = +1$ where the system undergoes a phase transition with spontaneous $Z_2$ symmetry breaking for $J_{AB} \neq 0$. We shall thus focus on antiferromagnetic $J_{AA} = -1$, which will be our energy and temperature scale of reference. The thermodynamics will be discussed as a function of the coupling ratio [25,33,30]

$$x = \frac{J_{AB}}{J_{AA}}, \quad (2)$$

with ferro- and antiferromagnetic $J_{AB}$.

## II. PHASE DIAGRAM

The Hamiltonian of equation (1) is invariant under the transformation

$$\begin{aligned} \sigma^A &\rightarrow -\sigma^A \\ J_{AB} &\rightarrow -J_{AB} \end{aligned} \quad (3)$$

All quantities derived from the energy, and especially the specific heat $C_h$ and entropy $S$, are thus the same for $x$ and $-x$. Their respective magnetic phases are related by reversing all spins of the A-sublattices.
A. Long-range order: $|x| > 1$

When the octagonal plaquettes are dominating ($x \to \pm\infty$), the shuriken lattice becomes a decorated square lattice, with A-sites sitting on the bonds between B-sites. Being bipartite, the decorated square lattice is not frustrated and orders via a phase transition of the 2D Ising Universality class³⁵ by spontaneous $Z_2$ symmetry breaking. Non-universal quantities such as the transition temperature can be exactly computed by using the decoration-iteration transformation³⁵⁻³⁷ [see appendix D 2]

$$T_c = \frac{2J_{AB}}{\ln(\sqrt{2} + 1 + \sqrt{2} + 2\sqrt{2})} \approx 1.30841J_{AB} \quad (4)$$

The low-temperature ordered phases, displayed in Fig. 2.(b) and 2.(c), remain the ground states of the anisotropic shuriken model for $x < -1$ and $x > 1$ respectively. The persistence of the 2D Ising Universality class down to $|x| = 1^+$ is not necessarily obvious, but is confirmed by finite-size scaling from Monte Carlo simulations [see appendix C].

These two ordered phases are respectively ferromagnetic (FM, $x < -1$) and staggered ferromagnetic (SFM, $x > 1$) [see Fig. 2.(b,c)]. The staggering of the latter comes from all spins on square plaquettes pointing in one direction, while the remaining ones point the other way. This leads to the rather uncommon consequence that fully antiferromagnetic couplings – both $J_{AA}$ and $J_{AB}$ are negative for $x > 1$ – induce long-range ordered (staggered) ferromagnetism, reminiscent of Lieb ferrimagnetism³⁸ as pointed out in Ref. [33] for quantum spins. The existence of ferromagnetic states among the set of ground states of Ising antiferromagnets is not rare, with the triangular and kagome lattices being two famous examples. But such ferromagnetic states are usually part of a degenerate ensemble where no magnetic order prevails on average. Here the lattice anisotropy is able to induce ferromagnetic order in an antiferromagnetic model by lifting its ground-state degeneracy at $|x| = 1$ (see below). This is interestingly quite the opposite of what happens in the spin-ice model³⁹, where frustration prevents magnetic order in a ferromagnetic model by stabilizing a highly degenerate ground state.

B. Binary paramagnet: $|x| < 1$

The central part of the phase diagram is dominated by the square plaquettes. The ground states are the same for all $|x| < 1$. A sample configuration of these ground states is given in Fig. 2.(d), where antiferromagnetically ordered square-plaquettes are separated from each other via spins on sublattice B. The antiferromagnetic square-plaquettes locally order in two different configurations equivalent to a superspin $\Xi$ with Ising degree of freedom.

$$\Xi = \sigma_1^A - \sigma_2^A - \sigma_3^A + \sigma_4^A = \pm 4, \quad (5)$$

where the site indices are given in Fig. 1. These superspins are the classical analogue of the tetramer objects observed in the spin-1/2 model³³. At zero temperature, the frustration of the $J_{AB}$ bonds perfectly decouples the superspins $\Xi$ from the B-sites. The system can then be seen as two interpenetrating square lattices: one made of superspins, the other one of B-sites. We shall refer to this phase as a binary paramagnet (BPM).

The perfect absence of correlations beyond square plaquettes at $T = 0$ allows for a simple determination of the thermodynamics. Let $N_{uc}$ and $N = 6N_{uc}$ be respectively the total number of unit cells and spins in the system, and $\langle X \rangle$ be the statistical average of $X$. There are $N_{uc}$ square plaquettes and $2N_{uc}$ B-sites, giving rise to an extensive ground-state entropy

$$S_{\text{BPM}} = k_B \ln (2^{N_{uc}} 2^{2N_{uc}}) = \frac{N}{2} k_B \ln 2 \quad (6)$$

which turns out to be half the entropy of an Ising paramagnet. As for the magnetic susceptibility $\chi$, it diverges as $T \to 0^+$. But the reduced susceptibility $\chi T$, which is nothing less than the normalized variance of the magnetization

$$\begin{aligned} \chi T &= \frac{1}{N} \left( \sum_{i,j} \langle \sigma_i \sigma_j \rangle - \langle \sigma_i \rangle \langle \sigma_j \rangle \right), \\ &= 1 + \frac{1}{N} \sum_{i \neq j} \langle \sigma_i \sigma_j \rangle, \end{aligned} \quad (7)$$

converges to a finite value in the BPM

$$\chi T|_{\text{BPM}} = \frac{1}{3}. \quad (8)$$

C. Classical spin liquid: $|x| \sim 1$

There is a sharp increase of the ground-state degeneracy at $|x| = 1$, when the binary paramagnet and the (staggered) ferromagnet meet. As is common for isotropic triangle-based Ising antiferromagnets, 6 out of 8 possible configurations per triangle minimize the energy of the system. As opposed to the BPM one does not expect a cutoff of the correlations [see section III C], making these phases cooperative paramagnets⁴⁰, also known as classical spin liquids.

Due to the high entropy of these cooperative paramagnets, the SL₁,₂ phases spread to the neighboring region of the phase diagram for $|x| \sim 1$ and $T > 0$, continuously connected to the high-temperature paramagnet [see Fig. 2]. Hence, for $|x| \gtrsim 1$, the anisotropic shuriken model stabilizes a cooperative paramagnet above a non-degenerate⁴¹ long-range ordered phase. This is a general property of classical spin liquids when adiabatically tuned away from their high-degeneracy point, as observed for example in Heisenberg antiferromagnets on the kagome⁴² or pyrochlore⁴³⁻⁴⁵ lattices, and possibly in the material of Er₂Sn₂O₇¹⁹. For $|x| \lesssim 1$ on the other hand,
<figure><img src="image_4.png" /></figure>

FIG. 3. Multiple crossovers between the paramagnetic, spin-liquids and binary regimes as observed in the specific heat $C_h$, entropy $S$ and reduced magnetic susceptibility $\chi T$. The models correspond to a) $x = \pm 1$, b) $x = \pm 0.9$ and c) $x = 0$. There is no phase transition for this set of parameters, which is why the Husimi tree calculations (lines) perfectly match the Monte Carlo simulations (circles) for all temperatures. The double crossover is present for $x = \pm 0.9$, with the low-temperature regime being the same as for $x = 0$, as confirmed by its entropy and susceptibility. The entropy is obtained by integration of $C_h/T$, setting $S(T \to +\infty) = \ln 2$. The vertical dashed lines represent estimates of the crossover temperatures determined by the local specific-heat maxima. The temperature axis is on a logarithmic scale. All quantities are given per number of spins and the Boltzmann constant $k_B$ is set to 1.

multiple crossovers take place upon cooling which de-
serves a dedicated discussion in the following section III.

### III. REENTRANCE OF DISORDER

#### A. Double crossover

First of all, panels (a) and (c) of Fig. 3 confirm that the classical spin liquids and binary paramagnet persist down to zero temperature for $x = \pm 1$ and $x = 0$ respectively, and that all models for $|x| \le 1$ have extensively degenerate ground states. For $x = \pm 0.9$ there is a double crossover indicated by the double peaks in the specific heat $C_h$ of Fig. 3.(b). These peaks are not due to phase transitions since they do not diverge with system size. The double crossover persists for $0.5 \lesssim |x| < 1$. Upon cooling, the system first evolves from the standard paramagnet to a spin liquid before entering the binary

paramagnet. The intervening spin liquid takes the form of an entropy plateau for $|x| = 0.9$ [see Fig. 3.(b)], at the same value as the low-temperature regime for $|x| = 1$ [see Fig. 3.(a)]. All relevant thermodynamic quantities are summarized in Table I.

While the mapping of equation (3) ensures the invariance of the energy, specific heat and entropy upon reversing $x$ to $-x$, it does not protect the magnetic susceptibility. The build up of correlations in classical spin liquids is known to give rise to a Curie-law crossover$^{46}$ between two $1/T$ asymptotic regimes of the susceptibility, as observed in pyrochlore$^{46-49}$, triangular$^{50}$ and kagome$^{50-52}$ systems. This is also what is observed here on the anisotropic shuriken lattice for $x = \{-1, 0, 1\}$ [see Fig. 4]. But for intermediate models with $x = \{-0.99, -0.9, 0.9, 0.99\}$, the double crossover makes the reduced susceptibility non-monotonic. $\chi T$ first evolves towards the values of the spin liquids $SL_1$ (resp. $SL_2$) for $x < 0$ (resp. $x > 0$) before converging to $1/3$ in the binary paramagnet.
<table><caption>TABLE I. Entropies *S* and reduced susceptibilities *χT* as *T* → 0⁺ for the anisotropic shuriken lattice with coupling ratios *|x|* ≪ 1. The results are obtained from Monte Carlo simulations, Husimi tree analytics and the exact solution for the binary paramagnet. All quantities are given per number of spins and the Boltzmann constant *k*<sub>B</sub> is set to 1.</caption><thead><tr><th>$T \to 0^+$</th><th>Monte Carlo</th><th>Husimi tree</th><th>exact</th></tr></thead><tbody><tr><td>$S(|x| = 1)$</td><td>0.504(1)</td><td>$\frac{1}{6} \ln \frac{41}{2} \approx 0.5034$</td><td>n/a</td></tr><tr><td>$\chi T(x = 1)$</td><td>0.203(1)</td><td>0.2028</td><td>n/a</td></tr><tr><td>$\chi T(x = -1)$</td><td>1.766(1)</td><td>1.771</td><td>n/a</td></tr><tr><td>$S(|x| < 1)$</td><td>0.347(1)</td><td>$\frac{1}{2} \ln 2 \approx 0.3466$</td><td>$\frac{1}{2} \ln 2$</td></tr><tr><td>$\chi T(|x| < 1)$</td><td>0.333(1)</td><td>$\frac{1}{3}$</td><td>$\frac{1}{3}$</td></tr></tbody></table>

<figure><img src="image_5.png" /><figcaption>FIG. 4. Reduced susceptibility χT with coupling ratios of *x* = ±1, ±0.99, ±0.9 and 0, obtained from Husimi-tree calculations (solid lines) and Monte Carlo simulations (circles). The Curie-law crossover of classical spin liquids is standard, i.e. χT is monotonic, for *x* = ±1 and 0, and takes a multi-step behavior for intermediate values of *x*, due to the double crossover. The characteristic values of the entropy and reduced susceptibility are given in Table I. The temperature axis is on a logarithmic scale</figcaption></figure>

Beyond the present problem on the shuriken lattice, this multi-step Curie-law crossover underlines the usefulness of the reduced susceptibility to spot intermediate regimes, and thus the proximity of different phases. From the point of view of renormalization group theory, the (*x*, *T*) = (±1, 0) coordinates of the phase diagram are fixed points which deform the renormalization flows passing in the vicinity.

## B. Decoration-iteration transformation

The phase diagram of the anisotropic shuriken model and, in particular, the double crossover observed for |*x*| < 1 [see Fig. 2] can be further understood using an exact mapping to an effective model on the checkerboard lattice, a method known as decoration-iteration transformation [see Ref. [37] for a review]. In short, by summing over the degrees of freedom of the A-spins, one can arrive at an effective Hamiltonian involving only the B-spins, which form a checkerboard lattice. The coupling constants of the effective Hamiltonian are functions of the temperature *T* and for |*x*| < 1 they vanish at both high and low temperatures, but are finite for an intermediate regime. This intermediate regime may be identified as the SL<sub>1,2</sub> cooperative paramagnets of Fig. 2, whereas the low-temperature region of vanishing effective interaction corresponds to the binary paramagnet (BPM). This mapping is able to predict a non-monotonic behavior of the correlation length.

In this section we give a brief sketch of the derivation of the effective model, before turning to its results. Details of the effective model are given in Appendix D.

To begin, consider the partition function for the system, with the Hamiltonian given by Eq. (1)

$$Z = \sum_{\{\sigma_i^A = \pm 1\}} \sum_{\{\sigma_i^B = \pm 1\}} \exp[-\beta H] \quad (9)$$

where $\beta = 1/T$ is the inverse temperature and the sums are over all possible spin configurations. Since in the Hamiltonian of Eq. (1) the square plaquettes of the A-sites are only connected to each other via their interaction with the intervening B-sites, it is possible to directly take the sum over configurations of A-spins in Eq. (9) for a fixed (but completely general) configuration of B-spins. Doing so, we arrive at

$$Z = \sum_{\{\sigma_i^B = \pm 1\}} \prod_{\Box} \mathcal{Z}_{\Box}(\{\sigma_i^B\}) \quad (10)$$

where the product is over all the square plaquettes of the lattice and $\mathcal{Z}_{\Box}(\{\sigma_i^B\})$ is a function of the four B-spins immediately neighbouring a given square plaquette. The B-spins form a checkerboard lattice, and Eq. (10) can be exactly rewritten in terms of an effective Hamiltonian $H_{\Box}$ on that lattice:

$$Z = \sum_{\{\sigma_i^B = \pm 1\}} \exp(-\beta \sum_{\Box} H_{\Box}) \quad (11)$$

$$H_{\Box} = -\mathcal{J}_0(T) - \mathcal{J}_1(T) \sum_{\langle ij \rangle} \sigma_i^B \sigma_j^B + \\ - \mathcal{J}_2(T) \sum_{\langle\langle ij \rangle\rangle} \sigma_i^B \sigma_j^B - \mathcal{J}_{\text{ring}}(T) \prod_{i \in \Box} \sigma_i^B \quad (12)$$

where $\sum_{\Box}$ is a sum over checkerboard plaquettes of B-spins. The effective Hamiltonian $H_{\Box}$ contains a constant term $\mathcal{J}_0$, a nearest neighbour interaction $\mathcal{J}_1$, a second
<figure><img src="image_6.png" /><figcaption>FIG. 5. Behavior of the coupling constants of the effective checkerboard lattice model as a function of temperature [Eq. (12)] for $x = -0.9$ (upper panel) and $x = -1$ (lower panel). Upper panel: All couplings vanish at both high and low temperatures with an intermediate regime at $T \sim 1$ where the effective interactions are stronger. The intermediate regime corresponds to the spin liquid region of the phase diagram Fig. 2, with the high- and low-temperature regimes corresponding to the paramagnet and binary paramagnet respectively. Lower panel: For all couplings $J_i$, $\beta J_i$ vanishes at high temperature and tends to a finite constant of magnitude $|\beta J_i(T)| \ll 1$ at low temperature. The short range correlated, spin-liquid regime, thus extends all the way down to $T = 0$.</figcaption></figure>

nearest neighbour interaction $J_2$, and a four-site ring interaction $J_{\text{ring}}$. All couplings are functions of temperature $J_i = J_i(T)$ and are invariant under the transformation $J_{AB} \mapsto -J_{AB}$ because the degrees of freedom of the A-sites have been integrated out. Expressions for the dependence of the couplings on temperature are given in Appendix D.

The temperature dependence of the effective couplings $J_i = J_i(T)$ can itself give rather a lot of information about the behavior of the shuriken model.

First we consider the case $|x| < 1$. In this regime of parameter space, all effective interactions $J_1, J_2, J_{\text{ring}}$ vanish exponentially at low temperature $T \ll 1$. For intermediate temperatures $T \sim 1$ the effective interactions in Eq. (12) become appreciable before vanishing once more at high temperatures. This is illustrated for the case $x = -0.9$ in the upper panel of Fig. 5. Seeing the problem in terms of these effective couplings gives some intuition into the double crossover observed in simulations. As the temperature is decreased the effective couplings

<figure><img src="image_7.png" /><figcaption>FIG. 6. Correlation lengths in the effective checkerboard model, calculated from Eq. (15), for $x = -0.9$ and $x = -1$. The correlation length is calculated to leading order in a perturbative expansion of the effective model in powers of $\beta J_i$. Such an expansion is reasonable for $|x| \le 1$ since $\beta J_i \ll 1$ for all $T$ (see Fig. 5). For $x = -0.9$ the behavior of the correlation length is non-monotonic. The correlation length is maximal in the spin liquid regime but correlations remain short ranged at all temperatures. In the binary paramagnet regime, the correlation length vanishes linearly at low temperature. For $x = -1$, the correlation length enters a plateau at $T \sim 1$, and short range correlations remain down to $T = 0$.</figcaption></figure>

$|J_i|$ increase in absolute value and the system enters a short range correlated regime. However, as the temperature decreases further, the antiferromagnetic correlations on the square plaquettes of A-spins become close to perfect, and act to screen the effective interaction between B-spins. This is reflected in the exponential suppression of the couplings $J_1, J_2, J_{\text{ring}}$.

In the case $|x| = 1$, the effective interactions $J_i$ no longer vanish exponentially at low temperature, but instead vanish linearly

$$J_1, J_2, J_{\text{ring}} \sim T. \qquad (13)$$

The ratio of effective couplings to the temperature $\beta J_i$ thus tends to a constant below $T \sim 1$, as shown in the lower panel of Fig. 5. Thus, the zero temperature limit of the shuriken model can be mapped to a finite temperature model on the checkerboard lattice for $|x| = 1$ and to an infinite temperature model for $|x| < 1$.

The behavior of the spin correlations in the shuriken model can be captured by calculating the correlation length between B-spins in the checkerboard model. Since $\beta J_i$ is small for all of the interactions $J_i$, at all temperatures $T$ (see Fig. 5), this can be estimated using a perturbative expansion in $\beta J_i$. For two B-spins chosen such that the shortest path between them is along nearest
<figure><img src="image_8.png" /><figcaption>FIG. 7. Spin-spin correlations in the vicinity of the spin liquid phases for $x = -0.9$ (a,b) $-1$ (c,d) and $-1.05$ (e,f), obtained from Monte Carlo simulations. The temperatures considered are $T = 0.01$ (♦), $1$ (■) and $891.25$ (●). Because of the anisotropy of the lattice, we want to separate correlation functions which start on A-sites (a,c,e) and B-sites (b,d,f). The radial distance is given in units of the unit-cell length. The agglomeration of data points around $C \sim 2.10^{-5}$ is due to finite size effects. The y-axis is on a logarithmic scale.</figcaption></figure>

neighbour $J_1$ bonds we obtain to leading order

$$ \langle \sigma_i^B \sigma_j^B \rangle = \exp \left( - \frac{r_{ij}}{\xi_{BB}} \right) \quad (14) $$

$$ \xi_{BB} \approx \frac{1}{\sqrt{2 \ln \left( \frac{T}{J_1(T)} \right)}} \quad (15) $$

where we choose units of length such that the linear size of a unit cell is equal to 1. Details of the calculation are given in Appendix D.

The correlation length between B-spins, calculated from Eq. (15), is shown for the cases $x = -0.9$ and $x = -1$ in Fig. 6. For $x = -0.9$ the correlation length

shows a non-monotonic behavior, vanishing at both high and low temperature with a maximum at $T \sim 1$. On the other hand for $x = -1$, the correlation length enters a plateau for temperatures below $T \sim 1$ and the system remains in a short range correlated regime down to $T = 0$. The extent of this plateau agrees with the low-temperature plateau of the reduced susceptibility in Fig. 4
<figure><img src="image_9.png" /><figcaption>FIG. 8. Static structure factors of the anisotropic shuriken lattice for (a) $x = -1$, (b) $x = 0$ and (c) $x = 1$ at zero temperature, obtained from Monte Carlo simulations. For $x = \pm 1$, the scatterings are strongly inhomogeneous (as opposed to a standard paramagnet) and non-divergent (i.e. without long-range order), confirming the spin liquid nature of these phases. The structure factors of the $x = +1$ and $x = -1$ models are similar by a $(q_x, q_y) = (60\pi, 0)$ or $(0, 60\pi)$ translation. The patterns are related to the 6-site unit cell of the shuriken lattice, as visible from Fig. 1. On the other hand for (b) $x = 0$, the black background underlines the absence of correlations in the binary paramagnet beyond the size of the superspins (square plaquettes), which is responsible for the finite extension of the dots of scattering. In order to restore ergodicity, a local update flipping the four spins of square plaquettes was used in the simulations. A video showing the temperature dependence of the static structure factor for $x = 0.9$ is available in the Supplementary Materials.</figcaption></figure>

### C. Correlations and Structure factors

The non-monotonic behavior of the correlation length estimated in the previous section III B can be measured by Monte Carlo simulations. Let us consider the microscopic correlations both in real ($C_\rho$) and Fourier ($S_q$) space. The function $C_\rho$ measures the correlation between a central spin $\sigma_0$ and all spins at distance $\rho$. Because of the nature of the binary paramagnet, one needs to make a distinction between central spins on the A and B sublattices. Let $D_\rho^X$ be the ensemble of sites at distance $\rho$ from a given spin $\sigma_0^X$ on the $X = \{A, B\}$ sublattice. The correlation function is defined as

$$C_{\rho}^{X} = \frac{\sum_{i \in D_{\rho}^{X}} |\langle \sigma_{0}^{X} \sigma_{i} \rangle|}{\sum_{i \in D_{\rho}^{X}}} \qquad (16)$$

where the absolute value accounts for the antiferromagnetic correlations. As for the static structure factor $S_q$, it is defined as

$$S_q = \langle \sigma_{\vec{q}} \sigma_{-\vec{q}} \rangle = \left\langle \left| \frac{1}{N_{uc}} \sum_i e^{-i\vec{q} \cdot \vec{r}_i} \sigma_i \right|^2 \right\rangle. \qquad (17)$$

$C_\rho^A$ and $C_\rho^B$ are respectively plotted on the left and right of Fig. 7. Let us first consider what happens in absence of reentrant behavior. For $x = -1.05$ [see panels (a,b)], the system is ferromagnetic at low temperature with $C(\rho) \approx 1$ over long length scales. Above the phase transition, the correlations are exponentially decaying.

When $x = 1$ [see panels (c,d)], the correlations remain exponentially decaying down to zero temperature. The correlation length $\xi$ reaches a maximum in the spin-liquid regime with $\xi \approx 0.3$. The quantitative superimposition of data for $T = 0.01$ and $T = 1$ is in agreement with the low-temperature plateau of the correlation length in Fig. 6. The spin liquid remains essentially unchanged all the way up to $T \sim 1$, when defects are thermally excited. However even if the correlations are exponential, they should not be confused with paramagnetic ones, as illustrated by their strongly inhomogeneous structure factors [see Fig. 8 and Supplementary Materials].

Once one enters the double-crossover region [see Fig. 7.(e,f) for $x = -0.9$], the correlation function becomes non-monotonic with temperature, as predicted from the analytics of Fig. 6. In the binary paramagnet, the B-sites are perfectly uncorrelated, while the A-sites have a finite cutoff of the correlation that is the size of the square plaquettes (superspins). This is why $S_q$ takes the form of an array of dots of scattering, whose width is inversely proportional to the size of the superspins [see Fig. 8]. Please note that the dip of correlations for the nearest-neighbors in Fig. 7.(e) is because half of the nearest-neighbors of any A-site are on the uncorrelated B sublattice.

The intervening presence of the spin liquids between the two crossovers is conceptually reminiscent of reentrant behavior$^{53-55}$. Not in the usual sense though, since reentrance is usually considered to be a feature of ordered phases surrounded by disordered ones. But the present scenario is a direct extension of the concept of reentrance applied to disordered regimes. This reentrance is quan-
titatively characterized at the macroscopic level by the double-peak in the specific heat, the entropy plateau and the multi-step Curie-law crossover of Fig. 3.(b), and microscopically by the non-monotonic evolution of the correlations [see Figs. 6, 7 and 8]. As such, it provides an interesting mechanism to stabilize a gas-like phase “below” a spin liquid, where (a fraction of) the spins form fully correlated clusters which i) can then fluctuate independently of the other degrees-of-freedom while ii) lowering the entropy of the gas-like phase below the one of the spin liquids.

## IV. THE SHURIKEN LATTICE IN EXPERIMENTS ?

Finally, we would like to briefly address the experimental situation. Unfortunately we are not aware of an experimental realization of the present model, but several directions are possible, each of them with their advantages and drawbacks.

The shuriken topology has been observed, albeit quite hidden, in the dysprosium aluminium garnet (DAG)$^{56,57}$ [see Ref. [58] for a recent review]. The DAG material has attracted its share of attention in the 1970's, but its microscopic Hamiltonian does not respect the geometry of the shuriken lattice – it is actually not frustrated – and is thus quite different from the model presented in equation (1). However it shows that the shuriken topology can exist in solid state physics.

Cold atoms might offer an alternative. Indeed, the necessary experimental setup for an optical shuriken lattice has been proposed in Ref. [32]. The idea was developed in the context of spin-ice physics, i.e. assuming an emergent Coulomb gauge theory whose intrinsic Ising degrees of freedom are somewhat different from the present model. Nonetheless, optical lattices are promising, especially if one considers that the inclusion of “proper” Ising spins might be available thanks to artificial gauge fields$^{59}$.

But the most promising possibility might be artificial frustrated lattices, where ferromagnetic nano-islands effectively behave like Ising degrees-of-freedom. Since the early days of artificial spin ice$^{60}$, many technological and fundamental advances have been made$^{61}$. In particular, while the thermalization of the Ising-like nano-islands had been a long-standing issue, this problem is now on the way to be solved$^{62-68}$. Furthermore, since the geometry of the nano-array can be engineered lithographically, a rich diversity of lattices is available, and the shuriken geometry should not be an issue. Concerning the Ising nature of the degrees-of-freedom, nano-islands have recently been grown with a magnetization axis $\vec{z}$ perpendicular to the lattice$^{69,70}$.

To compute their interaction$^{69,70}$, let us define the

Ising magnetic moment of two different nano-islands: $\vec{S} = \sigma\vec{z}$ and $\vec{S}' = \sigma'\vec{z}$. The interaction between them is dipolar of the form

$$D \left( \frac{\vec{S} \cdot \vec{S}'}{r^3} - 3 \frac{(\vec{S} \cdot \vec{r})(\vec{S}' \cdot \vec{r})}{r^5} \right) = \frac{D}{r^3} \sigma \sigma' \quad (18)$$

where D is the strength of the dipolar interaction and $\vec{r}$ is the vector separating the two moments. The resulting coupling is thus antiferromagnetic and quickly decays with distance. Hence, at the nearest-neighbour level, a physical distortion of the shuriken geometry – by elongating or shortening the distance between A and B sites – would precisely reproduce the anisotropy of equation (1) for $x > 0$. However, the influence of interactions beyond nearest-neighbours has successively been found to be experimentally negligible$^{69}$ and relevant$^{70}$ on the kagome geometry. Thus the phase diagram of Fig. 2.(a) could possibly be observed at finite temperature, but would likely be influenced by longer-range interactions at relatively low temperature.

## V. CONCLUSION

The anisotropic shuriken lattice with classical Ising spins supports a variety of different phases as a function of the anisotropy parameter $x = J_{AB}/J_{AA}$: two long-range ordered ones for $|x| > 1$ (ferromagnet and staggered ferromagnet) and three disordered ones [see Fig. 2]. Among the latter ones, we make the distinction, at zero temperature, between two cooperative paramagnets $SL_{1,2}$ for $x = \pm 1$, and a phase that we name a binary paramagnet (BPM) for $|x| < 1$. The BPM is composed of locally ordered square plaquettes separated by completely uncorrelated single spins on the B-sublattice [see Fig. 2.(d)].

At finite temperature, the classical spin liquids $SL_{1,2}$ spread beyond the singular points $x = \pm 1$, giving rise to a double crossover from paramagnet to spin liquid to binary paramagnet, which can be considered as a reentrant behavior between disordered regimes. This competition is quantitatively defined by a double-peak feature in the specific heat, an entropy plateau, a multi-step Curie-law crossover and a non-monotonic evolution of the spin-spin correlation, illustrated by an inhomogeneous structure factor [see Figs. 3, 4,6, 7 and 8]. The reentrance can also be precisely defined by the resurgence of the couplings in the effective checkerboard model [see Fig. 5].

Beyond the physics of the shuriken lattice, the present work, and especially Fig. 3, confirms the Husimi-tree approach as a versatile analytical method to investigate disordered phases such as spin liquids. Regarding classical spin liquids, Fig. 4 illustrates the usefulness of the reduced susceptibility $\chi T$ [46], whose temperature evolution quantitatively describes the successive crossovers between disordered regimes. Last but not least, we
hope to bring to light an interesting facet of distorted frustrated magnets, where extended regions of magnetic disorders can be stabilized by anisotropy, such as on the Cairo$^{71,72}$, kagome$^{51,73}$ and pyrochlore$^{74}$ lattices. Such connection is particularly promising since it expands the possibilities of experimental realizations, for example in Volborthite kagome$^{75}$ or breathing pyrochlores$^{76,77}$.

Possible extensions of the present work can take different directions. Motivated by the counter-intuitive emergence of valence-bond-crystals made of resonating loops of size 6 [30], the combined influence of quantum dynamics, lattice anisotropy $x^{30,33}$ and entropy selection presented here should give rise to a plethora of new phases and reentrant phenomena. As an intermediary step, classical Heisenberg spins also present an extensive degeneracy at $x = 1^{26,33}$, where thermal order-by-disorder is expected to play an important role in a similar way as for the parent kagome lattice, especially when tuned by anisotropy $x$. The addition of an external magnetic field$^{25,29}$ would provide a direct tool to break the invariance by transformation of equation (3), making the phase diagram of Fig. 2.(a) asymmetric. Furthermore, the diversity of spin textures presented here offers a promising framework to be probed by itinerant electrons coupled to localized spins via double-exchange.

## ACKNOWLEDGMENTS

We are thankful to John Chalker, Arnaud Ralko, Nic Shannon and Mathieu Taillefumier for fruitful discussions and suggestions. This work was supported by the Theory of Quantum Matter Unit of the Okinawa Institute of Science and Technology Graduate University.

## Appendix A: Methods

Classical Monte Carlo simulations have been performed based on the single-spin-flip algorithm. Let a Monte Carlo step (MCs) be the standard Monte Carlo unit of time made of $N$ attempts to flip a spin chosen at random. Typical simulations in this paper consist of

* $10^7$ MCs, including $10^6$ MCs for equilibration;

* 1 measurement every 50 MCs for $|x| > 1$ (total of 180 000 samples);

* 1 measurement every 10 MCs for $|x| \le 1$ (total of 900 000 samples);

* system sizes varying from $N = 2400$ to 15 000. Fig. 2 has been obtained for $N = 2400$ sites.

In order to improve the statistics, a large number of temperatures were simulated, and data were averaged over 4 neighboring temperatures.

To avoid any potential problems of ergodicity breaking in Monte Carlo simulations, we combined this numerical approach to analytical calculations on a Husimi tree$^{78}$, a method that has already demonstrated success in frustrated magnets$^{46,79–81}$. In a nutshell, the Husimi tree is a recursive approach on a Bethe lattice where all vertices are replaced by a cluster of spins. The clusters are connected to each other via their external corners, without making any closed loops. This allows to correctly take into account the interactions within each cluster, where frustration can be encoded.

Because the shuriken lattice is made of corner-sharing triangles, a natural choice would have been to consider triangles as building blocks of the Husimi-tree recursion. However a single triangle does not properly include the geometry of the anisotropy presented in Fig. 1. This is why, in the same way as for the 16-vertex model$^{82,83}$, we chose a larger building block made of four triangles forming a “shuriken” [see Fig. 1], which includes the anisotropy between A- and B-sites.

On the other hand, it neglects correlations on the length scale of the octagonal plaquettes and beyond. As such, the Husimi tree remains a mean field approximation which can only be qualitative in the vicinity of a critical point below its upper critical dimension. Since the 2D Ising Universality class is obviously not mean field, the Husimi tree underestimates the transition temperatures for $|x| > 1$ by a factor of $\approx 0.7$. This is why the boundaries of the FM and SFM phases have been determined with Monte Carlo simulations [open circles in Fig. 2.(a)].

But as far as disordered phases are concerned, the Husimi tree is quantitatively correct, as shown by Fig. 3 and table I. Being analytical, it provides an accurate way to determine the local maxima of the specific heat during crossovers [open triangles in Fig. 2.(a)].

## Appendix B: Pauling entropy of the spin liquids

In the isotropic case ($x = 1$), and by symmetry for $x = -1$ as well [see equation (3)], a simple Pauling argument is possible for the calculation of the entropy$^{84}$. If $N$ is the number of Ising spins, then there are $2N/3$ triangles in the system. Out of the $2^N$ possible configurations, the Pauling argument states that approximately $(6/8)^{2N/3}$ are allowed in the ground state, giving a total number of ground states in the spin liquids $SL_{1,2}$

$$ \Omega_{\text{SL-Pauling}} = 2^N \left(\frac{6}{8}\right)^{2N/3} = \left(\frac{9}{2}\right)^{N/3} \quad (\text{B1}) $$

giving an entropy

$$ \begin{aligned} S_{\text{SL-Pauling}} &= \frac{N}{3} k_B \ln \frac{9}{2} \\ &= \frac{N}{6} k_B \ln \frac{40.5}{2} \\ &\approx N k_B 0.50136 \end{aligned} \quad (\text{B2}) $$
<figure><img src="image_10.png" /><figcaption>FIG. 9. Finite size effects on the specific heat for $x = -3$ (left) and $x = -1.05$ (right). Insets: The transition temperature is scaled as a function of $1/L$ where $L$ is the linear system size. In the thermodynamic limit, we find $T_c = 2.788(5)$ for $x = -3$ and $T_c = 0.0714(5)$ for $x = -1.05$.</figcaption></figure>

<figure><img src="image_11.png" /><figcaption>FIG. 10. Finite size scaling of the magnetization $|M|$ (a,c) and susceptibility $\chi$ (b,d) for $x = -3$.</figcaption></figure>

The small difference between the Pauling estimate of equation (B2) and Monte Carlo results (0.50366) is mostly corrected by considering shurikens as building blocks in the Husimi-tree calculations (0.50340) [see Table I].

## Appendix C: 2D Ising Universality class

For $|x| > 1$, the anisotropic shuriken model orders at low temperature via a spontaneous $Z_2$ symmetry breaking [see Fig. 2]. We know it is a critical point of the 2D Ising Universality class for large $|x|$ [see section II A]. In

this appendix, our goal is to confirm numerically that it remains in the same Universality class as $|x| \to 1^+$, by considering two different values of the coupling ratio: $x = -3$ and $x = -1.05$. By symmetry of equation (3), the results directly apply to $x > 1$ also.

In Fig. 9, we analyze the specific heat $C_h$ for four different system sizes $N = \{2400, 5400, 9600, 15000\}$. The transition temperature scales like $1/N^{1/3}$ to its thermodynamic limit found at

$$x = -3 \Rightarrow T_c = 2.788(5) \qquad (C1)$$

$$x = -1.05 \Rightarrow T_c = 0.0714(5) \qquad (C2)$$

Based on these values of the transition temperature, we can define the reduced temperature $\epsilon = (T - T_c)/T_c$.
<figure><img src="image_12.png" /><figcaption>FIG. 11. Finite size scaling of the magnetization $|M|$ (a,c) and susceptibility $\chi$ (b,d) for $x = -1.05$.</figcaption></figure>

Following standard finite size scaling$^{85}$, we confirm in Figs. 10 and 11 that the nature of the phase transition is consistent with the 2D Ising Universality class with critical exponents

$$ \beta = 0.125, \quad \gamma = 1.75, \quad \nu = 1 \qquad (C3) $$

## Appendix D: Details of the decoration-iteration transformation

In this Appendix we give the details of the map to the effective model on the checkerboard lattice derived in Section III B. We give the derivation in Section D1 and then give details of the calculation of the correlation length in Section D3.

### 1. Derivation of the effective model on the checkerboard lattice

Consider the partition function of the anisotropic shuriken model

$$ Z = \sum_{\{\sigma_i^B = \pm 1\}} \sum_{\{\sigma_i^A = \pm 1\}} \exp[-\beta(H_{AA} + H_{AB})] \quad (D1) $$

where $H_{AA}$ and $H_{AB}$ are respectively the Hamiltonian of the square plaquettes of A- spins and the Hamiltonian coupling the intermediate B-spins to the square plaquette.

<figure><img src="image_13.png" /><figcaption>FIG. 12. The checkerboard lattice formed by the set of B-spins on the shuriken lattice.</figcaption></figure>

Summing over configurations of A- spins, we obtain:

$$ Z = \sum_{\{\sigma_i^B = \pm 1\}} \prod_{\square} \mathcal{Z}_{\square}(\{\sigma_i^B\}) \qquad (D2) $$

where the product is over all the square plaquettes of the lattice and $\mathcal{Z}_{\square}(\{\sigma_i^B\})$ depends on the configuration of the four B-spins immediately neighbouring a given square plaquette.

There are sixteen possible arrangements of the four spins B-spins surrounding a square plaquette of which
only four are inequivalent from the point of view of sym-

metry. These give rise to four possible values for $\mathcal{Z}^{\square}$:

$$ \mathcal{Z}_{++++} = 2(2 + 4 \cosh(4\beta J_{AB}) + \exp(-4\beta J_{AA}) + \exp(4\beta J_{AA}) \cosh(8\beta J_{AB})) \quad (D3) $$

$$ \mathcal{Z}_{++++} = 2(3 + 3 \cosh(4\beta J_{AB}) + \exp(-4\beta J_{AA}) + \exp(4\beta J_{AA}) \cosh(4\beta J_{AB})) \quad (D4) $$

$$ \mathcal{Z}_{+---} = 4(1 + 2 \cosh(4\beta J_{AB}) + \cosh(4\beta J_{AA})) \quad (D5) $$

$$ \mathcal{Z}_{+--} = 4(3 + \cosh(4\beta J_{AA})) \quad (D6) $$

From these we can assign “free energies” $F_i = -T \ln(\mathcal{Z}_i)$ to each of the four possible inequivalent configurations of B- spins around a square plaquette, i.e.

$$ F_{++++} = -T \ln(\mathcal{Z}_{++++}) \quad (D7) $$

$$ F_{+---} = -T \ln(\mathcal{Z}_{+---}) \quad (D8) $$

$$ F_{+---} = -T \ln(\mathcal{Z}_{+---}) \quad (D9) $$

$$ F_{+--} = -T \ln(\mathcal{Z}_{+--}) \quad (D10) $$

The B-spins form a checkerboard lattice as illustrated in Fig. 12. Using Eqs. (D3)-(D10) we can rewrite Eq. (D2) in terms of an effective Hamiltonian on the checkerboard lattice

$$ \mathcal{Z} = \sum_{\{\sigma_i^B = \pm 1\}} \exp \left[ -\beta \sum_{\Box} H_{\Box} \right] \quad (D11) $$

The sum $\sum_{\Box}$ is a sum over the elementary units of the checkerboard lattice. The function $H_{\Box}$ is a function only of the four B-spins around a checkerboard unit and returns one of the four $F_i$ defined in Eqs. (D7)-(D10) as appropriate to the configuration of those four spins.

We can rewrite $H_{\Box}$ explicitly in terms of interactions between the spins on the checkerboard lattice. The resultant effective Hamiltonian for the spins on the checkerboard lattice contains a constant term $\mathcal{J}_0$, a nearest neighbour interaction $\mathcal{J}_1$, a second nearest neighbour interaction $\mathcal{J}_2$, and a four-site ring interaction $\mathcal{J}_{\text{ring}}$.

$$ H_{\Box} = -\mathcal{J}_0(T) - \mathcal{J}_1(T) \sum_{\langle ij \rangle} \sigma_i^B \sigma_j^B - \mathcal{J}_2(T) \sum_{\langle\langle ij \rangle\rangle} \sigma_i^B \sigma_j^B - \mathcal{J}_{\text{ring}}(T) \prod_{i \in \Box} \sigma_i^B \quad (D12) $$

All couplings are functions of temperature $\mathcal{J}_i = \mathcal{J}_i(T)$.

The relationship between the temperature dependent couplings $\mathcal{J}_i(T)$ appearing in Eq. (D12) and the free en-

nergies $F_j$ defined in Eqs. (D7)-(D10) is:

$$ \mathcal{J}_0 = \frac{-1}{8} (F_{++++} + F_{+--} + 2F_{+---} + 4F_{+--+}) \quad (D13) $$

$$ \mathcal{J}_1 = \frac{-1}{8} (F_{++++} - F_{+--}) \quad (D14) $$

$$ \mathcal{J}_2 = \frac{-1}{8} (F_{++++} + F_{+--} - 2F_{+--+}) \quad (D15) $$

$$ \mathcal{J}_{\text{ring}} = \frac{-1}{8} (F_{++++} + F_{+--} + 2F_{+---} - 4F_{+--+}) \quad (D16) $$

We have thus succeeded in mapping the original model on the shuriken lattice, onto an effective model on the checkerboard lattice [Eq. (D12)].

## 2. Transition temperature of the decorated square lattice

In the limit $x \to +\infty$, one obtains the decorated square lattice. Applying $J_{AA} = 0$ to Eqs. (D3)-(D10) and then injecting the results into Eqs. (D14)-(D16), one obtains

$$ \mathcal{J}_1 = \frac{1}{2\beta} \ln(\cosh(2\beta J_{AB})), \quad (D17) $$

$$ \mathcal{J}_2 = \mathcal{J}_{\text{ring}} = 0. \quad (D18) $$

The term $\mathcal{J}_0$ does not cancel, but it only appears as a prefactor in the partition function of Eq. (D11) and thus does not influence the critical point.

Our effective model thereby becomes a square lattice with a temperature dependent nearest-neighbour coupling $\mathcal{J}_1(T)$. It is exactly soluble and the transition temperature $T_c = 1/\beta_c$ is obtained by injecting Eq. (D17) into Onsager’s solution of the Ising square lattice$^{86}$

$$ \begin{aligned} \beta_c \mathcal{J}_1(T_c) &= \frac{1}{2} \ln(\cosh(2\beta_c J_{AB})) \\ &= \frac{1}{2} \ln(\sqrt{2} + 1) \quad (\text{Onsager}) \quad (D19) \end{aligned} $$

which gives the result of Eq. (4)

$$ \begin{aligned} T_c &= \frac{2J_{AB}}{\ln(\sqrt{2} + 1 + \sqrt{2 + 2\sqrt{2}})} \\ &\approx 1.30841 J_{AB} \quad (D20) \end{aligned} $$
### 3. Correlation length

We observed in Section IIIB that for $x \le 1$ the couplings of the effective model are small compared to the temperature, for all values of temperature.

An expansion of the partition function of the effective model in powers of $\beta J_i$ is thus justified. Where $|x| < 1$ this expansion is asymptotically exact in both high and low temperature regimes.

Here we show how to use this expansion to calculate the correlation function $\langle \sigma_0^B \sigma_m^B \rangle$ for a pair of B-spins. For simplicity and concreteness we will do the calculation for a pair separated by a path such as that in Fig. 13, where the shortest route between them traverses only $J_1$ bonds and contains $m$ such bonds. However, there is no difficulty in making the calculation for other cases.

We have

<figure><img src="image_14.png" /><figcaption>FIG. 13. A path (in red) between two spins on the checkerboard lattice containing only nearest neighbour $J_1$ bonds. The correlation function between two such spins in the disordered regime is calculated in section D3.</figcaption></figure>

$$ \begin{aligned} \langle \sigma_0^B \sigma_m^B \rangle &= \frac{\sum_{\{\sigma_i \pm 1\}} \sigma_0^B \sigma_m^B \exp \left[ \beta \sum_{\Box} \mathcal{J}_0(T) + \mathcal{J}_1(T) \sum_{\langle ij \rangle} \sigma_i^B \sigma_j^B + \mathcal{J}_2(T) \sum_{\langle\langle ij \rangle\rangle} \sigma_i^B \sigma_j^B + \mathcal{J}_{\text{ring}}(T) \prod_{i \in \Box} \sigma_i^B \right]}{\sum_{\{\sigma_i \pm 1\}} \exp \left[ \beta \sum_{\Box} \mathcal{J}_0(T) + \mathcal{J}_1(T) \sum_{\langle ij \rangle} \sigma_i^B \sigma_j^B + \mathcal{J}_2(T) \sum_{\langle\langle ij \rangle\rangle} \sigma_i^B \sigma_j^B + \mathcal{J}_{\text{ring}}(T) \prod_{i \in \Box} \sigma_i \right]} \\ &= \frac{1}{N_c} \frac{\sum_{\{\sigma_i \pm 1\}} \sigma_0^B \sigma_m^B \sum_{n=0}^{\infty} \frac{1}{n!} \left[ \beta \sum_{\Box} \mathcal{J}_1(T) \sum_{\langle ij \rangle} \sigma_i^B \sigma_j^B + \mathcal{J}_2(T) \sum_{\langle\langle ij \rangle\rangle} \sigma_i^B \sigma_j^B + \mathcal{J}_{\text{ring}}(T) \prod_{i \in \Box} \sigma_i \right]^n}{1 + \frac{1}{N_c} \sum_{\{\sigma_i \pm 1\}} \sum_{n=1}^{\infty} \frac{1}{n!} \left[ \beta \sum_{\Box} \mathcal{J}_1(T) \sum_{\langle ij \rangle} \sigma_i^B \sigma_j^B + \mathcal{J}_2(T) \sum_{\langle\langle ij \rangle\rangle} \sigma_i^B \sigma_j^B + \mathcal{J}_{\text{ring}}(T) \prod_{i \in \Box} \sigma_i \right]^n} \end{aligned} \quad (D21) $$

where $N_c$ is the total number of spin configurations of the checkerboard model.

The leading non-zero term in Eq. (D21) comes from the $n=m$ part of the sum in the numerator, and corresponds to covering the shortest path between $\sigma_i^B$ and $\sigma_j^B$ with $J_1$ interactions. There are $m!$ ways of ordering the product of terms, which cancels the $1/n!$ occurring in the denominator. We thus obtain

$$ \langle \sigma_i^B \sigma_j^B \rangle \approx (\beta J_1(T))^m = \exp \left[ -m \ln \left( \frac{1}{\beta J_1(T)} \right) \right] \quad (D22) $$

In our choice of units of length, made such that the linear size of a unit cell equals 1, the distance between the spins is

$$ r = \frac{m}{\sqrt{2}} \quad (D23) $$

We therefore have a correlation length

$$ \xi_{BB} = \frac{1}{\sqrt{2} \ln \left( \frac{1}{\beta J_1(T)} \right)} \quad (D24) $$

1. L. Savary and L. Balents, Phys. Rev. Lett. 108, 037202 (2012).

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

3. Z. Hao, A. G. R. Day, and M. J. P. Gingras, Phys. Rev. B 90, 214430 (2014).

4. A. L. Chernyshev and M. E. Zhitomirsky, Phys. Rev. Lett. 113, 237202 (2014).

5. O. Götze and J. Richter, Phys. Rev. B 91, 104402 (2015).

6. J. Oitmaa and R. R. P. Singh,

Phys. Rev. B 93, 014424 (2016).

7. K. Essafi, O. Benton, and L. D. C. Jaubert, Nature Communications 7, 10297 (2016).

8. L. Messio, B. Bernu, and C. Lhuillier, Phys. Rev. Lett. 108, 207204 (2012).

9. Y.-C. He and Y. Chen, Phys. Rev. Lett. 114, 037201 (2015).

10. Y.-C. He, S. Bhattacharjee, F. Pollmann, and R. Moessner, Phys. Rev. Lett. 115, 267209 (2015).

11. S. Bieri, L. Messio, B. Bernu, and C. Lhuillier, Phys. Rev.
B **92**, 060407 (2015).

12 Y. Iqbal, H. O. Jeschke, J. Reuther, R. Valentí, I. I. Mazin, M. Greiter, and R. Thomale, Phys. Rev. B **92**, 220404 (2015).

13 P. A. McClarty, O. Sikora, R. Moessner, K. Penc, F. Pollmann, and N. Shannon, Phys. Rev. B **92**, 094418 (2015).

14 P. Henelius, T. Lin, M. Enjalran, Z. Hao, J. G. Rau, J. Altosaar, F. Flicker, T. Yavors'kii, and M. J. P. Gingras, Phys. Rev. B **93**, 024402 (2016).

15 H. D. Zhou, J. G. Cheng, A. M. Hallas, C. R. Wiebe, G. Li, L. Balicas, J. S. Zhou, J. B. Goodenough, J. S. Gardner, and E. S. Choi, Phys. Rev. Lett. **108**, 207206 (2012).

16 Z. L. Dun, M. Lee, E. S. Choi, A. M. Hallas, C. R. Wiebe, J. S. Gardner, E. Arrighi, R. S. Freitas, A. M. Arevalo-Lopez, J. P. Attfield, H. D. Zhou, and J. G. Cheng, Phys. Rev. B **89**, 064401 (2014).

17 C. R. Wiebe and A. M. Hallas, APL Materials **3**, 041519 (2015).

18 J. G. Rau and M. J. P. Gingras, Phys. Rev. B **92**, 144417 (2015).

19 H. Yan, O. Benton, L. D. C. Jaubert, and N. Shannon, arXiv:1311.3501 (2013).

20 L. D. C. Jaubert, O. Benton, J. G. Rau, J. Oitmaa, R. R. P. Singh, N. Shannon, and M. J. P. Gingras, Phys. Rev. Lett. **115**, 267208 (2015).

21 J. Robert, E. Lhotel, G. Remenyi, S. Sahling, I. Mirebeau, C. Decorse, B. Canals, and S. Petit, Phys. Rev. B **92**, 064425 (2015).

22 H. Nakano and T. Sakai, Journal of the Physical Society of Japan **82**, 083709 (2013).

23 R. Siddharthan and A. Georges, Physical Review B **65**, 014417 (2001).

24 J. Richter, O. Derzhko, and J. Schulenburg, Physical Review Letters **93**, 107206 (2004).

25 O. Derzhko and J. Richter, The European Physical Journal B **52**, 23 (2006).

26 J. Richter, J. Schulenburg, P. Tomczak, and D. Schmalfuss, Condensed Matter Physics **12**, 507 (2009).

27 O. Derzhko, J. Richter, O. Krupnitska, and T. Krokhmalskii, Phys. Rev. B **88**, 094426 (2013).

28 H. Nakano, T. Sakai, and Y. Hasegawa, Journal of the Physical Society of Japan **83**, 084709 (2014).

29 H. Nakano, Y. Hasegawa, and T. Sakai, Journal of the Physical Society of Japan **84**, 114703 (2015).

30 A. Ralko and I. Rousochatzakis, Physical Review Letters **115**, 167202 (2015).

31 P. Tomczak and J. Richter, Journal of Physics A: Mathematical and General **36**, 5399 (2003).

32 A. W. Glaetzle, M. Dalmonte, R. Nath, I. Rousochatzakis, R. Moessner, and P. Zoller, Phys. Rev. X **4**, 041037 (2014).

33 I. Rousochatzakis, R. Moessner, and J. van den Brink, Physical Review B **88**, 195109 (2013).

34 O. Derzhko, J. Richter, O. Krupnitska, and T. Krokhmalskii, Low Temperature Physics **40**, 513 (2014), 1312.1111.

35 S. Chun-Feng, K. Xiang-Mu, and Y. Xun-Chang, Communications in Theoretical Physics **45**, 555 (2006).

36 M. E. Fisher, Phys. Rev. **113**, 969 (1959).

37 J. Strecka and M. Jascur, Acta Physica Slovaca **65**, 235–367 (2015).

38 E. H. Lieb, Phys. Rev. Lett. **62**, 1201 (1989).

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

40 J. Villain, Zeitschrift Für Physik B-Condensed Matter **33**, 31 (1979).

41 Besides the trivial time-reversal symmetry.

42 M. Elhajal, B. Canals, and C. Lacroix, Phys. Rev. B **66**, 014422 (2002).

43 B. Canals, M. Elhajal, and C. Lacroix, Phys. Rev. B **78**, 214431 (2008).

44 G.-W. Chern, arXiv (2010), arXiv:1008.3038 [cond-mat.str-el].

45 P. A. McClarty, P. Stasiak, and M. J. P. Gingras, Phys. Rev. B **89**, 024425 (2014).

46 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).

47 S. V. Isakov, K. S. Raman, R. Moessner, and S. L. Sondhi, Phys. Rev. B **70**, 104418 (2004).

48 I. A. Ryzhkin, Journal of Experimental and Theoretical Physics **101**, 481 (2005).

49 P. H. Conlon and J. T. Chalker, Phys. Rev. B **81** (2010).

50 M. Isoda, Journal of Physics: Condensed Matter **20**, 315202 (2008).

51 W. Li, S.-S. Gong, Y. Zhao, S.-J. Ran, S. Gao, and G. Su, Phys. Rev. B **82**, 134434 (2010).

52 A. J. Macdonald, P. C. W. Holdsworth, and R. G. Melko, Journal of Physics-Condensed Matter **23**, 164208 (2011).

53 J. Hablützel, Helv. Phys. Acta **12**, 489 (1939).

54 V. Vaks, A. Larkin, and Y. N. Ovchinnikov, Soviet Physics JETP **22**, 820 (1966).

55 P. Cladis, Molecular Crystals and Liquid Crystals **165**, 85 (1988).

56 D. P. Landau, B. E. Keen, B. Schneider, and W. P. Wolf, Phys. Rev. B **3**, 2310 (1971).

57 W. P. Wolf, B. Schneider, D. P. Landau, and B. E. Keen, Phys. Rev. B **5**, 4472 (1972).

58 W. P. Wolf, Brazilian Journal of Physics **30**, 794 (2000).

59 J. Struck, M. Weinberg, C. Ölschläger, P. Windpassinger, J. Simonet, K. Sengstock, R. Höppner, P. Hauke, A. Eckardt, M. Lewenstein, and L. Mathey, Nature Physics **9**, 738 (2013).

60 R. F. Wang, C. Nisoli, R. S. Freitas, J. Li, W. McConville, B. J. Cooley, M. S. Lund, N. Samarth, C. Leighton, V. H. Crespi, and P. Schiffer, Nature **439**, 303 (2006).

61 C. Nisoli, R. Moessner, and P. Schiffer, Rev. Mod. Phys. **85**, 1473 (2013).

62 V. Kapaklis, U. B. Arnalds, A. Harman-Clarke, E. T. Papaioannou, M. Karimipour, P. Korelis, A. Taroni, P. C. W. Holdsworth, S. T. Bramwell, and B. Hjorvarsson, New Journal of Physics **14** (2012), 10.1088/1367-2630/14/3/03500

63 A. Farhan, P. M. Derlet, A. Kleibert, A. Balan, R. V. Chopdekar, M. Wyss, L. Anghinolfi, F. Nolting, and L. J. Heyderman, Nature Physics **9**, 1 (2013).

64 A. Farhan, P. M. Derlet, A. Kleibert, A. Balan, R. V. Chopdekar, M. Wyss, J. Perron, A. Scholl, F. Nolting, and L. J. Heyderman, Phys. Rev. Lett. **111**, 057204 (2013).

65 J. P. Morgan, J. Akerman, A. Stein, C. Phatak, R. M. L. Evans, S. Langridge, and C. H. Marrows, Phys. Rev. B **87**, 024405 (2013).

66 C. Marrows, Nature Physics **9**, 324 (2013).

67 L. Anghinolfi, H. Luetkens, J. Perron, M. G. Flokstra, O. Sendetskyi, A. Suter, T. Prokscha, P. M. Derlet, S. L. Lee, and L. J. Heyderman, Nature Communications **6**, 8278 (2015).

68 U. B. Arnalds, J. Chico, H. Stopfel, V. Kapaklis, O. Brenbold, M. A. Verschuuren, U. Wolff, V. Neu, A. Bergman, and B. Hjorvarsson,
New Journal of Physics **18**, 023008 (2016).

69 S. Zhang, J. Li, I. Gilbert, J. Bartell, M. J. Erickson, Y. Pan, P. E. Lammert, C. Nisoli, K. K. Kohli, R. Misra, V. H. Crespi, N. Samarth, C. Leighton, and P. Schiffer, Phys. Rev. Lett. **109**, 087201 (2012).

70 I. A. Chioar, N. Rougemaille, A. Grimm, O. Fruchart, E. Wagner, M. Hehn, D. Lacour, F. Montaigne, and B. Canals, Phys. Rev. B **90**, 064411 (2014).

71 I. Rousochatzakis, A. M. Läuchli, and R. Moessner, Phys. Rev. B **85**, 104415 (2012).

72 M. Rojas, O. Rojas, and S. M. de Souza, Phys. Rev. E **86**, 051116 (2012).

73 W. Apel and H.-U. Everts, Journal of Statistical Mechanics: Theory and Experiment **2011**, P09002 (2011).

74 O. Benton and N. Shannon, Journal of the Physical Society of Japan **84**, 104710 (2015).

75 Z. Hiroi, M. Hanawa, N. Kobayashi, M. NHara, H. Takagi, Y. Kato, and M. Takigawa, Journal of the Physical Society of Japan **70**, 3377 (2001).

76 Y. Okamoto, G. J. Nilsen, J. P. Attfield, and Z. Hiroi, Phys. Rev. Lett. **110**, 097203 (2013).

77 K. Kimura, S. Nakatsuji, and T. Kimura, Phys. Rev. B **90**, 060414 (2014).

78 K. Husimi, Journal of Chemical Physics **18**, 682 (1950).

79 P. Chandra and B. Doucot, Journal of Physics A-Mathematical and General **27**, 1541 (1994).

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

81 J. Strečka and C. Ekiz, Phys. Rev. E **91**, 052143 (2015).

82 L. Foini, D. Levis, M. Tarzia, and L. Cugliandolo, Journal of Statistical Mechanics, P02026 (2013).

83 D. Levis, L. Cugliandolo, L. Foini, and M. Tarzia, Phys. Rev. Lett. **110**, 207206 (2013).

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

85 D. P. Landau and K. Binder, *A guide to Monte Carlo simulation in statistical physics* (Cambridge University Press, Berlin, 2009).

86 L. Onsager, Physical Review **65**, 117 (1944).