Rico Pohle

Owen Benton

L D C Jaubert

Reentrance of disorder in the anisotropic shuriken Ising model

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

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.

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][2][3][4][5][6][7] or further nearest-neighbour couplings [8][9][10][11][12][13][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][16][17][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][24][25][26][27][28][29][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 longrange 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][25][26][27][28][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 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 JAA (red), while interactions between A-and B-sites (octagonal plaquettes) are described with JAB (black).

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 longrange 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 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 doublecrossover 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 Bsites on the octagonal plaquettes. The Hamiltonian of the model can be written as:

H = -J AA ij AA σ A i σ A j -J AB ij AB σ A i σ B j (1)
where we consider Ising spins σ i = ±1 with nearestneighbor 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 = 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 = J AB J AA , (2)
with ferro-and antiferromagnetic J AB .

II. PHASE DIAGRAM

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

σ A → -σ A J AB → -J AB (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 → ±∞), the shuriken lattice becomes a decorated square lattice, with A-sites sitting on the bonds between Bsites. Being bipartite, the decorated square lattice is not frustrated and orders via a phase transition of the 2D Ising Universality class 35 by spontaneous Z 2 symmetry breaking. Non-universal quantities such as the transition temperature can be exactly computed by using the decoration-iteration transformation [35][36][37] [see appendix D 2]

T c = 2J AB ln √ 2 + 1 + 2 + 2 √ 2 ≈ 1.30841J AB (4)
The low-temperature ordered phases, displayed in Fig. 2 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 38 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 39 , 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 squareplaquettes locally order in two different configurations equivalent to a superspin Ξ with Ising degree of freedom.

Ξ = σ A 1 -σ A 2 -σ A 3 + σ A 4 = ±4, (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 33 . At zero temperature, the frustration of the J AB bonds perfectly decouples the superspins Ξ 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 = 6 N uc be respectively the total number of unit cells and spins in the system, and X 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 BPM = k B ln 2 Nuc 2 2Nuc = N 2 k B ln 2 (6)
which turns out to be half the entropy of an Ising paramagnet. As for the magnetic susceptibility χ, it diverges as T → 0 + . But the reduced susceptibility χ T , which is nothing less than the normalized variance of the magnetization

χ T = 1 N   i,j σ i σ j -σ i σ j   , = 1 + 1 N i =j σ i σ j , (7)
converges to a finite value in the BPM

χ T | BPM = 1 3 . (8)

C. Classical spin liquid: |x| ∼ 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 40 , also known as classical spin liquids.

Due to the high entropy of these cooperative paramagnets, the SL 1,2 phases spread to the neighboring region of the phase diagram for |x| ∼ 1 and T > 0, continuously connected to the high-temperature paramagnet [see Fig. 2]. Hence, for |x| 1, the anisotropic shuriken model stabilizes a cooperative paramagnet above a nondegenerate 41 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 42 or pyrochlore [43][44][45] lattices, and possibly in the material of Er 2 Sn 2 O 7 19 . For |x| 1 on the other hand, multiple crossovers take place upon cooling which deserves 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 = ±1 and x = 0 respectively, and that all models for |x| 1 have extensively degenerate ground states. For x = ±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 |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. 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][47][48][49] , triangular 50 and kagome [50][51][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 nonmonotonic. χ 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. 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 multistep 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 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 1,2 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 = {σ A i =±1} {σ B i =±1} exp [-βH] (9)
where β = 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 Asites 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 = {σ B i =±1} Z ({σ B i }) ( 10
)
where the product is over all the square plaquettes of the lattice and Z ({σ B i }) 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 ⊠ on that lattice:

Z = {σ B i =±1} exp(-β ⊠ H ⊠ ) (11)
H ⊠ = -J 0 (T ) -J 1 (T ) ij σ B i σ B j + -J 2 (T ) ij σ B i σ B j -J ring (T ) i∈⊠ σ B i ( 12
)
where ⊠ is a sum over checkerboard plaquettes of Bspins. The effective Hamiltonian H ⊠ contains a constant term J 0 , a nearest neighbour interaction J 1 , a second nearest neighbour interaction J 2 , and a four-site ring interaction J ring . All couplings are functions of temperature J i = J i (T ) and are invariant under the transformation J AB -→ -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 ring vanish exponentially at low temperature T << 1. For intermediate temperatures T ∼ 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 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 βJi. Such an expansion is reasonable for |x| ≤ 1 since βJi << 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 ∼ 1, and short range correlations remain down to T = 0.

|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 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 ring ∼ T. (13)
The ratio of effective couplings to the temperature βJ i thus tends to a constant below T ∼ 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 β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 βJ i . For two B-spins chosen such that the shortest path between them is along nearest 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 ∼ 2.10 -5 is due to finite size effects. The y-axis is on a logarithmic scale.

neighbour J 1 bonds we obtain to leading order

σ B i σ B j = exp - r ij ξ BB (14)
ξ BB ≈ 1 √ 2 ln T J1(T ) (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 ∼ 1. On the other hand for x = -1, the correlation length enters a plateau for temperatures below T ∼ 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 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 = ±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 (qx, qy) = (60π, 0) or (0, 60π) 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.

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 ρ ) and Fourier (S q ) space. The function C ρ measures the correlation between a central spin σ 0 and all spins at distance ρ. 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 X ρ be the ensemble of sites at distance ρ from a given spin σ X 0 on the X = {A, B} sublattice. The correlation function is defined as

C X ρ = i∈D X ρ | σ X 0 σ i | i∈D X ρ ( 16
)
where the absolute value accounts for the antiferromagnetic correlations. As for the static structure factor S q , it is defined as

S q = σ q σ -q = 1 N uc i e -i q• ri σ i 2 . ( 17
)
C A ρ and C 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(ρ) ≈ 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 ξ reaches a maximum in the spin-liquid regime with ξ ≈ 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 ∼ 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 Asites 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][54][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 frustratedand 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][63][64][65][66][67][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 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: S = σ z and S ′ = σ ′ z. The interaction between them is dipolar of the form

D S • S ′ r 3 -3 ( S • r)( S ′ • r) r 5 = D r 3 σ σ ′ ( 18
)
where D is the strength of the dipolar interaction and 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 siteswould 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 = ±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 = ±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 χ 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. 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 D 1 and then give details of the calculation of the correlation length in Section D 3.

Derivation of the effective model on the checkerboard lattice

Consider the partition function of the anisotropic shuriken model

Z = {σ B i =±1} {σ A i =±1} exp [-β(H AA + H AB )] (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 plaque- ttes. Summing over configurations of A-spins, we obtain:

Z = {σ B i =±1} Z ({σ B i }) (D2)
where the product is over all the square plaquettes of the lattice and Z ({σ B i }) 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 Z :

Z ++++ = 2(2 + 4 cosh(4βJ AB ) + exp(-4βJ AA ) + exp(4βJ AA ) cosh(8βJ AB )) (D3) Z +++-= 2(3 + 3 cosh(4βJ AB ) + exp(-4βJ AA ) + exp(4βJ AA ) cosh(4βJ AB )) (D4) Z ++--= 4(1 + 2 cosh(4βJ AB ) + cosh(4βJ AA )) (D5) Z +-+-= 4(3 + cosh(4βJ AA )) (D6)
From these we can assign "free energies" F i = -T ln(Z i ) to each of the four possible inequivalent configurations of B-spins around a square plaquette, i.e.

F ++++ = -T ln(Z ++++ ) (D7) F +++-= -T ln(Z +++-) (D8) F ++--= -T ln(Z ++--) (D9) F +-+-= -T ln(Z +-+-) (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

Z = {σ B i =±1} exp -β ⊠ H ⊠ (D11)
The sum ⊠ is a sum over the elementary units of the checkerboard lattice. The function H ⊠ 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 ⊠ 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 J 0 , a nearest neighbour interaction J 1 , a second nearest neighbour interaction J 2 , and a four-site ring interaction J ring .

H ⊠ = -J 0 (T ) -J 1 (T ) ij σ B i σ B j -J 2 (T ) ij σ B i σ B j -J ring (T ) i∈⊠ σ B i (D12)
All couplings are functions of temperature J i = J i (T ).

The relationship between the temperature dependent couplings J i (T ) appearing in Eq. (D12) and the free en-ergies F j defined in Eqs. (D7)-(D10) is: The term 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 J 1 (T ). It is exactly soluble and the transition temperature T c = 1/β c is obtained by injecting Eq. (D17) into Onsager's solution of the Ising square lattice which gives the result of Eq. ( 4)

J 0 = - 1
T c = 2J AB ln √ 2 + 1 + 2 + 2 √ 2
≈ 1.30841 J AB (D20)

Correlation length

We observed in Section III B that for x ≤ 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 βJ i is thus justified. Where |x| < 1 this expansion is assymptotically exact in both high and low temperature regimes.

Here we show how to use this expansion to calculate the correlation function σ B 0 σ B m 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 FIG. 13. A path (in red) between two spins on the checkerboard lattice containing only nearest neighbour J1 bonds. The correlation function between two such spins in the disordered regime is calculated in section D 3.

σ B 0 σ B m = {σi±1} σ B 0 σ B m exp β ⊠ J 0 (T ) + J 1 (T ) ij σ B i σ B j + J 2 (T ) ij σ B i σ B j + J ring (T ) i∈⊠ σ B i {σi±1} exp β ⊠ J 0 (T ) + J 1 (T ) ij σ B i σ B j + J 2 (T ) ij σ B i σ B j + J ring (T ) i∈⊠ σ i = 1 N c {σi±1} σ B 0 σ B m ∞ n=0 1 n! β ⊠ J 1 (T ) ij σ B i σ B j + J 2 (T ) ij σ B i σ B j + J ring (T ) i∈⊠ σ i n 1 + 1 Nc {σi±1} ∞ n=1 1 n! β ⊠ J 1 (T ) ij σ B i σ B j + J 2 (T ) ij σ B i σ B j + J ring (T ) i∈⊠ σ B i n (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 σ B i and σ B j 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

σ B i σ B j ≈ (βJ 1 (T )) m = exp -m ln 1 βJ 1 (T ) (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 = m √ 2 (D23)
We therefore have a correlation length

ξ BB = 1 √ 2 ln 1 βJ1(T ) (D24)

FIG. 2 .

 2 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| 1, on cooling, the system undergoes a evolution from "gas crossover ------→ liquid transition ------→ solid". As for |x| 1, it provides a remarkable example of reentrance from "gas crossover ------→ liquid crossover ------→ gas".

FIG. 3 .

 3 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 χ T . The models correspond to a) x = ±1, b) x = ±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 = ±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 → +∞) = 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 kB is set to 1.

  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

FIG. 5 .

 5 FIG. 5.Behavior of the coupling constants of the 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 ∼ 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 Ji, βJi vanishes at high temperature and tends to a finite constant of magnitude |βJi(T )| << 1 at low temperature. The short range correlated, spin-liquid regime, thus extends all the way down to T = 0.

FIG. 9 .

 9 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 Tc = 2.788(5) for x = -3 and Tc = 0.0714(5) for x = -1.05

FIG. 11 .

 11 FIG. 11. Finite size scaling of the magnetization |M | (a,c) and susceptibility χ (b,d) for x = -1.05.

FIG. 12 .

 12 FIG.12. The checkerboard lattice formed by the set of Bspins on the shuriken lattice.

8 (FF

 8 ++++ + F +-+-+ 2F ++--+ 4F +++-) ++++ + F +-+--2F ++--) (D15) J ring = -1 8 (F ++++ + F +-+-+ 2F ++---4F +++-) (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 .J 2 =

 22 Transition temperature of the decorated square lattice In the limit x → +∞, 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 J ring = 0. (D18)

86

 86

β c J 1 (

 1

TABLE I .

 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 kB is set to 1.

	T → 0 +	Monte Carlo			Husimi tree	exact
	S(|x| = 1)	0.504(1)	1 6	ln	41 2	≈ 0.5034	n/a
	χ T (x = 1)	0.203(1)				0.2028	n/a
	χ T (x = -1) 1.766(1)				1.771	n/a
	S(|x| < 1)	0.347(1)	1 2	ln 2 ≈ 0.3466	1 2	ln 2
	χ T (|x| < 1)	0.333(1)					1 3	1 3

FIG.

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| ≤ 1 (total of 900 000 samples);

• system sizes varying from N = 2 400 to 15 000. Fig. 2 has been obtained for N = 2 400 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][80][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 ≈ 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), 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 ⇒ T c = 2.788 (5) (C1) x = -1.05 ⇒ T c = 0.0714 (5) (C2)

Based on these values of the transition temperature, we can define the reduced temperature ǫ = (T -T c )/T c .