Fluctuation-induced spin nematic order in magnetic charge-ice

Disorder in materials may be used to tune their functionalities, but much more strikingly, its presence can entail entirely new behavior. This happens in charge-ice where structural disorder is not weak and local, but strong and long-range correlated. Here, two cations of different charge occupy a pyrochlore lattice, arranging themselves such that all tetrahedra host two cations of each type. The ensuing correlated disorder is characterized by randomly packed loops of a single cation-type. If the cations are magnetic and interact antiferromagnetically, a new type of magnet with strong interactions along the loops, but frustrated interactions between loops, emerges. This results in an ensemble of intertwined Heisenberg spin chains that form an algebraic spin liquid at intermediate temperatures. At lower temperatures, we find these non-local degrees of freedom undergo a discontinuous transition to a spin nematic. While this phase does not break time reversal symmetry, its spin symmetry is reduced resulting in a dramatically slower spin relaxation. The transition is sensitive to the statistics of the cation loops, providing both a direct thermodynamic signature of otherwise elusive structural information and a structural route to engineering nematic phase stability.

The discrete translational symmetry of the crystalline state underlies the utility of many functional materials.Introducing random disorder into crystalline materials can play a crucial role in modifying their static and dynamical properties to obtain new or improved functionality, for example by producing pinning centers in superconductors 1 or tuning transition temperatures in multiferroics 2,3 .Recently it has been suggested that a kind of correlated disorder based on tiling high symmetry lattices with low symmetry motifs may be a route to novel functionalities via the interplay of the disorder with crystal properties such as lattice dynamics or electronic conduction 4 .Magnetism is another material property that may be controlled by disorder.Usually uncorrelated variations of exchange strength or coordination (via uncorrelated doping of magnetic ions) are expected to produce spin glasses 5 .Here, we show that more correlated types of structural disorder may result in distinct equilibrium and out-of-equilibrium properties.
An ice rule is a well known way to introduce such correlated disorder.Here, a simple constraint on the local configuration of binary degrees of freedom allows the construction of an extensively degenerate manifold of states, in which the correlation among local degrees of freedom decays not exponentially, but by a 1/r 3 (dipolar) powerlaw 8 .Such rules have become common currency for describing spin configurations in geometrically frustrated magnetic materials and arrays of nano-magnets, notably spin ice 9 , quantum spin ice 10 , and artificial spin ice 11,12 .The equivalence of spins and charges on this lattice was first noted by Anderson in an investigation of the Vewey transition in magnetite 13 , where it was pointed out that cations would obey the charge-ice rule which requires that each tetrahedron is occupied by two cations of each type.If the different cations carry magnetic moments one obtains a model of magnetic charge-ice, which is directly relevant to pyrochlores of the type AMM'F 6 (for example CsCrNiF 6 , see Ref. 14 ).In pyrochlores such as R 2 MM'O 7 and AA'M 2 F 7 , see respective Refs. 15,16, the charge-ice formed amongst the non-magnetic spectator ions may introduce more subtly correlated bond disorder amongst the magnetic atoms.More generally, geometric frustration of charge order causes correlated distributions of species and thus of the interactions among their associated degrees of freedom, resulting in specific material properties that reflect the correlated nature of the underlying disordered structure.
Strong and correlated magnetic exchange disorder via a charge-ice.Here we consider how an exchange network templated by a charge-ice cation configuration affects the low temperature properties of a classical Heisenberg spin system.In a simple model of magnetic charge-ice 7 , two types of magnetic atoms a and b populate the pyrochlore lattice according to the 2:2 charge-ice rule resulting in correlated site disorder characterised by a randomly packed set of single cationtype nearest-neighbour connected loops of even length (Fig. 1a).Fig. 1b shows the loop size distribution in which we distinguish four classes of loops: those that are non-winding or winding with respect to the periodic boundaries of the system, and, additionally, the largest and second-largest winding loop.For a given charge-ice realisation these latter two will be of different chemical type and we refer to them as giant loops.The fraction of sites occupied by the four loop classes tends (with increased sampling and system size) to f nw = 0.06 for nonwinding, f w = 0.22 for winding, f 2 = 0.31 and f 1 = 0.41 for the second and largest (giant) loops respectively, in agreement with Ref. 6 .
We describe the magnetic structure by unit-length arXiv:2311.05004v1[cond-mat.stat-mech]8 Nov 2023 FIG. 1. Charge-ice correlated exchange disorder and its low temperature magnetic properties.a) The charge-ice structure of two types of cations, a and b, distributed on the pyrochlore lattice is characterized by an ensemble of closed loops of the same atom species.Two loops of different species (blue and red) may share a number of tetrahedra.b) Sample-averaged normalized histogram of loop lengths in which non-winding and winding are distinguished, as well as the two largest (giant) loops.The power-law exponents are well understood from the perspective of diffusion 6 .The normalized distributions of loop length l are plotted as a function of l/N where N = 16L 3 with L = 20.c) Magnetic charge-ice is characterized by the nearest neighbour exchange constants Jaa, J bb , and J ab between moments associated with species a and b.The ground-state phase diagram 7 is found to be rich in structure with regions I-III hosting long-range ordered phases.The boundary to region IV (JaaJ bb = J 2 ab ) hosts a classical Heisenberg pyrochlore AFM at (−1, −1).Within IV the ground state is a (less degenerate) classical spin liquid with perfect AFM order on each loop but no correlations between loops.d) At low temperature, the heat capacity of region IV reveals a transition toward a phase in which the giant and other large loops are collinear.This nematic transition is reflected in the structure factor (upper inset) which has a smooth pinch point structure above the transition and becomes patchy below it.
classical Heisenberg spins on the sites, that are connected by the nearest neighbour exchange constants J aa , J bb , and J ab , into which we absorb the size of the different cation magnetic moments.The resulting Heisenberg Hamiltonian displays correlated bond disorder that derives from the spatial structure of the cation loops.Banks and Bramwell 7 identified four regions of the ground-state phase diagram for this model, as shown in Fig. 1c.We focus on region IV, where J aa J bb > J 2 ab , with J aa and J bb both promoting intra-species antiferromagnetic (AFM) alignment, so that the zero temperature ground states have perfect AFM arrangements on each loop, but are degenerate with respect to the orientation of the AFM alignment axis (the Néel vector) of any loop due to the inter-loop couplings J ab being perfectly frustrated.In the work of Ref. 7 , Monte Carlo simulations at T /|J| < 0.12 for regime IV (J aa = J bb = J with J/|J ab | = −1.2) revealed a pinch-point-like structure factor and a vanishing Edwards-Anderson parameter down to T /|J| = 0.012 (indicating no spin freezing), suggesting an algebraic spin liquid.
Low temperature spin nematic order.Focusing on the case J aa = J bb = J < 0 with J/|J ab | = −2 (region IV), our Monte Carlo heat bath algorithm 18 simulations FIG. 2. A first order nematic phase transition.a) The one dimensional AFM loop structure factor progressively increases as the temperature is reduced.Above T0 all loops exhibit 1D domain wall activity and thus short range order at a finite temperature characteristic of finite Heisenberg spin-chains (thin lines represent the J ab = 0 case of non-interacting loops).However, around T0, the giant loops show a discontinuous jump in their AFM order as their Néel vectors align.Such alignment is diagnosed by the bulk quadrupolar order parameter Q, b), whose average magnitude, ⟨|Q|⟩/N (|Q| 2 ≡ Tr[Q 2 ]), jumps at T0 concomitant with strong fluctuations c) defined through the generalized susceptibility χ |Q| = (⟨|Q 2 |⟩ − ⟨|Q| 2 ⟩)/(T N ).The first order nature of this phase transition is revealed by the generalized Binder cumulant 17 which becomes negative close to T0 indicating a bi-modal order parameter distribution (histograms of |Q| below, close-to, and above T0 are shown in the inset).In all panels, an estimate of the infinite size T0/J equal to 0.0103 is indicated by the vertical dashed lines derived from the Binder cumulant crossing seen in panel d).
show that on cooling (0.125 ≲ T /|J| ≲ 1.5) the system evolves from the paramagnetic state into a low temperature state with an energy per site approaching that of the expected ground state value, a strongly suppressed magnetization, and a well developed plateau in the heat capacity reflecting the low temperature behaviour of a classical Heisenberg system (Fig. 1d).At our lowest temperatures, the heat capacity c V has a value just below unity, indicating a significantly more constrained system than the pyrochlore Heisenberg AFM (PHAFM), which asymptotes to c V = 3/4 (see Ref. 19 and Fig. 3a).A small peak in c V at T 0 /|J| ≈ 0.01 (Fig. 1d), which sharpens with increasing system size, suggests a previously unnoticed phase transition, whose nature we now elucidate.
The static structure factor (Fig. 1d, inset panels), taken as and thermally averaged over statistically independent spin configurations for a single charge-ice realisation, shows the distinctive diffuse scattering and pinch-points associated with dipolar spin correlations on the pyrochlore lattice above T 0 .Below T 0 this pinch-point structure becomes patchy, like that of similarly sized individual charge-ice ground states 7 .The pinch-point structure arises because all spins that share the same loop are AFM correlated, which implies power-law spin correlations.This is not unlike the case of the PHAFM where the ground state manifold consists of all possible AFM close-packed loop realizations combining to give a smooth diffuse scattering profile at these system sizes.The patchiness is therefore due to a restricted sub-set of the full PHAFM manifold, which does not self-average at our finite system size.It can be removed by averaging over many charge-ice realizations or by considering larger system sizes.
To reveal the structure of the low T phase we investigate the 1D AFM structure factor of the ith loop: where l i is its length.
When S loop AFM,i = 1, the loop has complete AFM order, while the sign of its Néel vector may still fluctuate.Fig. 2a plots S loop AFM,l for loops of various sizes as a function of temperature.Generally, they develop smoothly as the temperature is reduced, but for the two giant loops, the structure factor jumps up abruptly at a temperature T 0 .As a reference, data are also shown for the case of non-interacting loops (J ab = 0).Above T 0 , loops of all sizes in the full system behave similarly as noninteracting loops, for which the thermal properties are known analytically 20,21 .Indeed, rescaling the temperature axis of the J ab = 0 data by the factor 0.8 results in almost perfect overlap with the J ab ̸ = 0 data for T > T 0 , suggesting the full system is well described in this temperature regime by an ensemble of non-interacting spinchains with the renormalized coupling ∼ 0.8J.
Inspection of the low temperature spin configurations reveals that for T < T 0 the Néel vectors of the two giant loops align collinearly, motivating the use of the bulk quadrupolar or nematic order parameter 22 where Q is a traceless symmetric tensor with components Q µν i = s µ i s ν i − δ µν /3.A non-zero ⟨Q⟩ signals breaking of rotational symmetry, but not necessarily of time reversal symmetry, being invariant under spin reversal (of entire loops).Figs.2b and c plot the average of the magnitude of the quadrupolar order parameter and a measure of its fluctuations near the transition, indicating a rapid turn-on of quadrupolar order that sharpens with increasing system size.Since Q and −Q describe qualitatively different spin structures the Landau free energy does not need to be invariant under a sign change of Q and will generally contain a cubic term Tr(Q 3 ), ruling out a continuous phase transition.This is confirmed by the generalised Binder cumulant 17 for Q, which becomes increasingly negative just above T 0 with increasing system size (Fig. 2d), due to a bi-modal distribution of the order-parameter magnitude reflecting phase coexistence at T 0 (see inset in Fig. 2d).
Role of loop lengths and loop-loop coupling.These results might suggest the giant loops are the essential ingredient for the transition to occur.This is not the case, since breaking up the giant loops through a modified charge-ice algorithm (Fig. SM1) or using open boundary conditions (not shown) has little effect on the transition as long as sufficiently large loops remain present.The observation that for a given J/|J ab |, sufficiently small loops do not order, suggests that an ordered charge-ice structure, consisting of 4L 2 linear loops of length 4L (system I, see Methods) should not nematically order for a small enough L. Indeed, for L = 8, order is absent for moderate J/|J ab | = −2, and only sets in for J/|J ab | ≳ −4/3 (Fig. SM3).
Order-by-disorder and symmetry reduction from Heisenberg to Ising loops.Insight into why the nematic structure is selected can be gained from the low temperature thermal properties of a single tetrahedron.Expanding the corresponding Hamiltonian to quadratic order with respect to transverse fluctuations around a ground state configuration defined by the angle ϕ between the Néel vectors of the two species (Sec.SM1 C 1), yields a fluctuational entropy, ∆S (2) ab /J 2 cos 2 ϕ, favoring collinear alignment of the Néel vectors.This indicates the observed first-order transition is driven by an order-by-disorder mechanism 19,[23][24][25] at the tetrahedral level.A more accurate estimate of the entropic advantage of the nematic state is obtained via a similar quadratic calculation for a magnetic charge-ice ground state, in which the free energy due to transverse fluctuations of AFM correlated loops with aligned Néel vectors is compared to that of loops with randomly oriented Néel vectors.The entropy of nematic order with a finite Q z = ⟨cos 2 ϕ⟩−1/3 exceeds that of randomly aligned configurations by 1 2 Q z J 2 ab /J 2 per spin (Sec.SM1 C 2), which is a number comparable to the entropy gain of a single tetrahedron.
Using this entropic interaction, we first investigate the possibility of loop alignment at temperatures T sufficiently low such that the correlation length L T = J/T of a Heisenberg chain exceeds the length of a loop i.In this regime, it is well characterised by its Néel vector n i .Two such loops, i and j, therefore experience the entropic interaction −T N ij ∆S (2) [(n i •n j ) 2 ], where N ij is the number of shared tetrahedra, and whose leading order term is Here (J ab /J) 2 can be viewed as the entropic coupling parameter.For the simple loop connectivity of the ordered charge-ice structure, a meanfield description may be developed (Sec.SM1 D 1) in which each equivalent loop is embedded in a symmetry breaking quadrupolar field, Q z = ⟨n 2 z ⟩ − 1/3.Requiring self-consistency then gives the temperature independent condition on the loop length, l > l 0 ≈ 13.5(J/J ab ) 2 for Q z to become non-zero.For longer loops the entropy gain l∆S (2) ∆Q z (∆Q z ≈ 0.2) from partial alignment overcompensates the entropy lost, O(1k B ), from constraining the fluctuations of the loop Néel vectors with the nematic order emerging discontinuously at loop lengths l = l 0 .Thus, the nematic phase disappears in charge-ice structures with too short loops and/or too weak entropic couplings.This also suggests that in an ordered structure of alternating short loops and winding loops of the size of the system, the nematic transition is strongly suppressed due to the strong fluctuations of the small loops.This is indeed the case for the ordered charge-ice system II (see methods and Fig. SM2b) where no signature of the nematic transition is seen (Fig. SM3).The presence of the nematic phase and the value of T c are thus both sensitive to the distribution of loop length and their intertwining (connectivity).The fact that short loops tend to fluctuate strongly also rationalizes the ordering tendency in a general charge-ice structure, in particular the results of Fig. 3, where progressively smaller loops align as J/|J ab | increases towards -1.
The density of normal mode frequencies of the linearized charge-ice Hamiltonian shows the 1D spinchain asymptotic form ρ(ε) ∼ 1/ √ ε for small ε (Fig. SM4), motivating a 1D Heisenberg spin chain Hamiltonian in the presence of a symmetry breaking mean-field anisotropy term, −Q z J 2 ab /J 2 (s 2 zi − 1/3) (see Sec. SM1 D 2).In the long-wavelength continuum limit, this is solved exactly via numerical transfer-integral methods 26 , predicting a first order nematic transition at T 0 /J ≈ 0.05cJ 2 ab /J 2 where c ≲ 1 represents the overall fraction of tetrahedra touched by two loops long enough to undergo nematic alignment.Note that at T 0 , the correlation length L T (T 0 ) becomes of order l 0 , such that entropy gained from nematic alignment compensates the entropy lost by constraining the fluctuations of correlated loop segments.This rationalizes the observed temperature scale of T 0 and its decrease as J/|J ab | becomes more negative.Indeed, for J/|J ab | = −2 the mean-field prediction gives T 0 = 0.0125c, which is remarkably close to that of simulation (Fig. 2) for c ≈ 0.95 which, for charge-ice, is the approximate fraction of sites involved in loops longer than l 0 .Mean-field theory also explains the anomalous temperature dependence of the heat capacity in the nematic phase, giving c V ≈ 1−3/4 T J 2 ab /J 3 (Figs.1d and  3a) and tracing it to the quenching of the entropy of the softest spin waves due to the increasingly strong entropic interaction (Sec.SM1 D 3).
An s 2 z anisotropy does not break time-reversal symmetry and no long range spin order is expected within the nematic phase.However, such a spontaneously emerging anisotropy reduces the O(3) global symmetry of the Heisenberg Hamiltonian to an Ising Z 2 symmetry, which remains unbroken on the chains in accord with the Mermin-Wagner theorem 27 .Above T 0 , loops fluctuate and equilibrate rapidly due to long wavelength spin waves, whereas below T 0 the reduced spin symmetry entails a many orders of magnitude larger spin relaxation time due to the tiny Gibbs factor exp(−J/T ) associated with the nucleation and separation of an Ising domainwall.With such kinetics nearly frozen out, the loops maintain their nearly perfect AFM order for very long times, with spin relaxation times of order τ 1 ∼ exp(J/T ).
Conclusions and outlook.While discontinuous transitions were found in related frustrated systems upon perturbing homogeneously the interactions and thereby lifting the ground state degeneracy [28][29][30][31] , those are driven by the essentially local competition between energy and entropy.In contrast, charge-ice establishes a complex connectivity among strongly correlated non-local cluster degrees of freedom, which reflects the precise realisation of the correlated disorder -and it is with respect to these degrees of freedom that the first order transition takes place.The predicted spin nematic breaks spin rotation symmetry, but preserves (statistical) lattice symmetries.It is thus quite distinct from lattice nematics, that break lattice rotation invariance at the level of the spinspin correlation function 32,33 .The continuous rotational symmetry of the Heisenberg Hamiltonian is reduced to a discrete Ising symmetry, entailing an emergent slow dynamics and a new type of sudden spin-liquid freezing, in which sufficiently large loops fall out of equilibrium and become AFM ordered on mesoscopic timescales.This differs strongly from the effect of random couplings, which may induce glassy spin freezing 34 , with slow dynamics deriving from a complex energy landscape, but occurring at temperatures far below the dominant exchange energy scale.
Our work shows that correlated structural disorder can produce non-trivial behavior due to the emergence of non-local degrees of freedom tied to lower-dimensional clusters (loops/strings).Solids in which similarly correlated disorder is known (or expected) to exist are numerous 35 , with corner-sharing tetrahedra being only one example of a more general class of materials whose corner or edge-sharing plaquettes may show qualitatively different magnetic behavior 25,36,37 .Moreover, transferring the paradigm of interacting non-local intertwined magnetic degrees of freedom that arise from correlated disorder to the realm of continuous phase transitions might offer the possibility of entirely new universality classes 38 .
Quantitatively understanding the relation between such correlated structural disorder and emergent collective degrees of freedom and their thermodynamic signatures is a formidable but not intractable problem.Indeed, experimentally observing the predicted nematic phase transition through magnetic birefringence would give indirect evidence for the existence of large loops and the presence of correlated disorder.Moreover, if it is possible to vary the exchange constants, either chemically or through a global distortion, and monitor the transition temperature and the order parameter magnitude, one might extract additional information on the distribution of loop lengths, establishing an experimental link between correlated disorder and the thermodynamics it entails.
Acknowledgments The authors wish to thank Sam Garratt, Afonso Dos Santos Rufino, and Hugo Bocquet for helpful discussions.We also thank Christian Rüegg for doctoral supervision of AH.The work is supported by the European Union Horizon 2020 research and innovation program under the Marie Skodowska-Curie Grant agreement No. 884104 (PSI-FELLOW-III-3i) and the Swiss National Science Foundation (grant number 200020 182536).
Author Contributions PMD and TF instigated the project; AH, PMD, KE and MT performed the simulations; AH, KE, TF, and PMD carried out the analyses; PMD and MM made the theoretical calculations; PMD, TF and MM wrote the paper with input from the other authors.
Methods: Monte Carlo A single-site Monte Carlo approach was found to be sufficient for the present work.Since a wide range of temperature scales are probed, the Monte Carlo heat bath algorithm was found to be most suitable.Here, an MC move entails randomly selecting a site and calculating exactly the probability density function for that spin with all other spins fixed.This distribution is then sampled to find the new state of the chosen spin.Whilst there is a computational cost in sampling this distribution, it has the advantage of all moves being accepted and of automatically reducing the scale of variations in spin as the temperature is decreased.For more details see, for example, Ref. 18 .
Methods: Sample Creation To produce a pyrochlore sample satisfying the charge-ice constraint on each tetrahedron, the pyrochlore lattice of size L (containing 16L 3 atoms) is constructed and initially populated with a and b sites according to an ordered structure consisting of [110] and [1 10] chains of sites of one or the other type of cation, respectively.Under periodic boundary conditions, this may be seen as a regular array of 4L 2 winding loops of length 4L.The connectivity of such a structure is characterized by any two loops sharing either zero or one tetrahedron.This initial structure will be referred to as an ordered charge-ice system I.To disorder it, a loop consisting of alternating site types is identified via a worm algorithm and all site types are interchanged, preserving the charge-ice structure.This procedure is repeated until variations in loop structure satisfy the known statistical properties of the loops as detailed in Figs.1b and in Ref. 6 .These samples will be referred to as a charge-ice system.An alternative ordered charge-ice system may be constructed consisting of (100) planes of [110] chains of sites separated by regions fully populated by hexagonal loops of length l = 6.This is referred to as the ordered charge-ice system II and contains L 2 loops of length 4L and 2L 3 of length 6. See Fig. SM2 which visualizes both charge ordered systems.

SM1. SUPPLEMENTARY MATERIAL A. Breaking up the giant loops
To investigate the robustness of the observed phase transition with respect to the size of the giant loops, we perform Monte Carlo simulations on an L = 8 system for which the loop structure generation was biased to generating smaller loops.This bias was achieved by only allowing changes in the structure which reduced the sum of the square of loop lengths.In particular, this procedure was applied to the L = 8 sample used in the main text, resulting in a sample (referred to as the "small loop" sample) with over 54 loops, the largest ten of which had lengths 1386, 1376, 1236, 1218, 1058, 826, 426, 140, 86, and 30.This should be compared with the original sample which had 35 loops, the largest four of which are 3926, 3908, 60 and 26 in length.Fig. SM1a displays the resulting heat capacity compared to the original L = 8 charge-ice system showing little change in the transition temperature T 0 .Fig. SM1b displays the loop-loop orientation correlation ⟨Tr( Ql1 Ql2 )⟩/(⟨|Q l1 | 2 ⟩⟨|Q l2 | 2 ⟩) 1/2 below the critical temperature, demonstrating that the growth in the bulk quadrupolar order parameter is due to the alignment of these larger non-giant loops.Fig. SM1c visualizes the 8 largest loops of the sample, showing that all but the eighth largest loop are winding.

B. Ordered charge-ice
Fig. SM2 displays the two ordered charge ice systems I and II for the case of L = 8.Both structures are investigated for inter-chain couplings J/|J ab | equal to -2, -4/3, -8/7, -16/15.Fig. SM3 displays a) the heat capacity and b) the average magnitude of the bulk quadrupolar order parameter for system I showing that for J/|J ab | = −2 the first order phase transition is entirely suppressed, as predicted by mean field theory (Sec.SM1 D 1).However as J/|J ab | becomes less negative, the transition appears, with T 0 again increasing as J/|J ab | approaches -1, as for the case of the general charge-ice structure (Figs.1-3 in the main text).Fig. SM3 also shows similar data for the ordered charge-ice system II where due to the linear loops not sharing any tetrahedra, interacting only via hexagonal loops of length six that do not align, the nematic transition is generally absent for all choices of J/|J ab | except very close to the PHAFM case where the hexagons also begin to align.For both ordered charge-ice structures, when no transition is observed, the heat capacity plateaus to a value equal to 1 − N l /(16L 2 ), where N l is the number of loops.This originates from the zero modes within the system 4 , which for charge ice is equal to twice the number of loops.For ordered charge ice I, N l = 4L 2 and for ordered charge ice II N l = L 2 + 2L 3 giving the respective heat capacity plateaus of 0.97 and 0.87.

C. Harmonic transverse spin fluctuations
The classical Heisenberg spin Hamiltonian may be written as for which the local field at each site i is To investigate the transverse spin fluctuations, s ⊥,i with respect to a given spin configuration ŝ0,i each spin is written as where If the magnetic configuration ŝ0,i is at a local energy minimum then all ŝ0,i will be parallel to their local fields, b 0,i .Then to quadratic order in the transverse components, the Hamiltonian may be written as H = H (0) + ∆H (2) where

FIG. SM2.
Ordered charge ice structures.The two ordered charge ice structures, a) system I and b) system II.The loops can be identified via similarly coloured bonds with the colour also reflecting the cation type.In system I, each linear loop will share either zero or at most one tetrahedron with any other linear loop.In system II, the linear loops do not share tetrahedra between themselves.
In the above, the off-diagonal term J ij is the full 3D Hessian whereas the second diagonal term is a correction to the 3D Hessian which projects the taken derivatives onto the tangent space of each spin.
Representing the 2D tangent space of spin i as e 1,i and e 2,i with e 1,i × e 2,i = ŝ0,i , the ith spin may be written as where χ α,i are real numbers.The choice of e α,i is not unique and we follow Ref. 1 .Together the above yields a symmetric matrix M iα,jβ of rank 2N , represented as an N × N matrix of 2 × 2 block elements, whose (i, j)th block element is Λ ij e α,i • e β,j .Solving the corresponding eigen-problem yields the normal modes of Eqn.SM4 that govern the fluctuations of this quadratic Hamiltonian.It is noted that for disordered/frustrated systems, ŝ0,i and thus the local tangent space, defined via e α,i , will be different for each spin.Thus the normal modes presently calculated are non-trivially related to the corresponding spin-wave modes which arise from a linearisation of the Landau-Lifschitz equation.
At the level of the quadratic approximation to the Hamiltonian, the resulting partition function becomes a simple Gaussian integral, evaluating to where the λ n are the non-zero eigenvalues of the fluctuation matrix M , from which the free energy may be calculated as F = −T log Z giving In the thermodynamic limit this can be evaluated as an integral ∞ 0 + dλ ρ(λ) log λ using the density of eigenvalues (or density of states DOS), ρ(λ), which is normalized to 2N .

Single tetrahedron
The Hamiltonian for a single tetrahedron satisfying the charge-ice rule is given by where we recall that we focus on the parameter regime where the J a and J b couplings are negative (AFM).For a ground state configuration of region IV, we have the AFM configurations between spins of the same type: ŝa,1 = −ŝ a,2 and ŝb,1 = −ŝ b,2 , and an angle ϕ between the alignment axis.This gives the ground state energy −J aa − J bb independent of ϕ.Using the formalism of the previous section, the quadratic Hamiltonian is represented as a matrix of rank 8: whose four non-zero eigenvalues give the free energy contribution The above can be conveniently written as ∆F (2) [T, ϕ] = ∆F (2) [T ] − T ∆S (2) [cos 2 ϕ] (SM10) where (with k B = 1) is the (temperature independent) fluctuational entropy evaluated for a given angle ϕ between the orientations of the two equal species pairs.Thus the angle-constrained free energy has minima at ϕ = 0, π and maxima at ϕ = ±π/2.Alignment or anti-alignment thus results in maximal fluctuational entropy, where ∆S (2) ≡ ∆S (2) (ϕ = 0) − ⟨∆S (2) (ϕ)⟩ ϕ ≈ 0.1 for J ab / √ J aa J bb = 1/2.

Full system
A similar harmonic analysis may be carried out for the full charge-ice system, where now the eigen-system of the Hessian M (calculated via Eqn.SM4) must be solved numerically for a particular choice of the T = 0 reference configuration.Fig. SM4 displays the normal mode density of states (DOS) for two ground-state configurations: one with nematic order, in which all loop Néel vectors are aligned; and one where they are randomly orientated with respect to each other (random loop AFM or RLA).These states are both members of the manifold of ground states identified by Banks and Bramwell 2 , and are indistinguishable in terms of their internal energies.The DOS of the RLA depends somewhat on the particular ground state configuration, but for sufficiently large samples self-averaging reduces such differences.For smaller samples an average over many choices of random alignment results in a converged DOS.Both RLA and nematic order reveal 2N loop zero-modes reflecting the individual O(3) symmetry of each AFM loop, and whilst there are differences between the nematic and random ground state configurations (e.g. the enhanced density of low frequency states and more discrete structure at higher frequencies in the nematic ground state), the similarities at low frequency are more revealing.In particular, a log-log plot (inset of Fig. SM4) reveals the asymptotic form ρ(λ) ∼ 1/ √ λ for small λ -a hallmark signature of the fluctuation spectrum of AFM-ordered reference configurations of 1D spin chains.For comparison the DOS derived from the harmonic Hessian with J ab = 0 is also shown, which consists of N loop non-interacting finite 1D antiferromagnetic spin chains.This is in agreement with the known analytical form (apart from finite size corrections due to a small fraction of short loops) ρ(λ) ∝ λ(4J − λ).
Within the harmonic approximation, the difference in fluctuational entropy between nematic and RLA states is given by − ∞ 0 + dλ (ρ nematic (λ) − ρ RLA (λ)) ln(λ) which we find to be a positive quantity.This originates from the coupling-induced softening of low frequency normal modes, which enhances the 1/ √ λ tail.This softening effect is strongest for the nematically aligned configuration (hence the enhancement at low frequency regime compared to a randomly aligned ground state configuration), which is thus entropically favored.

Nematic alignment
In what follows, we consider loops of characteristic length l at low temperatures T ≪ |J|/l, such that the persistence (or correlation) length of an infinite Heisenberg chain L T = |J|/T ≫ l.In this limit, each loop i can be characterised by a single Néel vector, ni (which may still flip slowly), and exhibits fast but small transverse fluctuations around it.From Sec.SM1 C 1 the total free energy of an ensemble of such loops is where N ij is the number of shared tetrahedra between loops i and j.This gives the temperature independent partition function, As a specific example we consider a periodic chargeordered system of size L, which contains N l = 4L 2 loops of length l = 4L, where each loop shares either zero or one tetrahedron with any other loop.Thus N ij = 1 and where the j summation spans the l loops that share a tetrahedron with the ith loop.A mean-field construction is performed by replacing the summand by an average with respect to nj , yielding an effective single loop weight: where in the last equality we have expanded ∆S (2) (n i • nj ) 2 to leading order in J ab , assuming J aa = J bb = J.
We now determine whether it is consistent to assume that dyads of ni acquire a finite expectation value.Choosing the z axis as the symmetry breaking axis we assume Q z = ⟨n 2 z ⟩ − 1/3, where a non-zero value would spontaneously break the rotational invariance.Assuming rotational symmetry around the z-axis, ⟨n 2 x ⟩ = ⟨n 2 y ⟩ = (1 − ⟨n 2 z ⟩)/2, and ⟨n x ⟩ = ⟨n y ⟩ = 0, the above average evaluates to Self consistency of the mean field now requires that ⟨n 2 z ⟩− 1/3 computed with the effective single loop weight equal Q z .Calling λ = 3(J ab /J) 2 /4, we thus seek the stable solution of the mean field equation: For large λ, Q z tends to 2/3.This symmetry breaking solution disappears at lλ ≈ 10.1 = (lλ) 0 where the order parameter discontinuously drops from Q z ≈ 0.205 at (lλ) 0 to zero, signaling a first order transition.It is noted that a spinodal instability of the disordered phase exists at (lλ) sp = 45/4 = 11.25, which is however preempted by the above first order transition -as required for a nematic transition, which cannot be continuous.Low temperature order is thus predicted to exist only for loops larger than l 0 = (lλ) 0 /λ = 4(lλ) 0 (J/J ab ) 2 /3.For the case of our ordered charge-ice where l = 4L, the nematic phase transition will only occur for periodic samples of size L, when J ab /|J| > J 0 ab /|J| = (4(lλ) 0 /3l) 1/2 = ((lλ) 0 /3L) 1/2 .For the case of L = 8 this requires J ab /J > 0.65 or J/|J ab | < −1.53.
A fully self-consistent mean-field theory with respect to normalized loop distributions for type a and b, P a/b (l), now follows by writing the average quadrupolar field component for sites of type a/b as , where loop l 1 is a giant loop whose quadrupolar field is oriented along ẑ with magnitude Q a/b z and l 2 is a small loop of type b/a.Fig. SM5b compares this to the data of Fig. 3c showing very good agreement, and quantitatively confirming the initial assumption entailed in Eqn.SM12 and the general mean-field approach.For this system Q a z /(2/3) = 0.983 and Q b z /(2/3) = 0.971, giving the effective fraction of tetrahredra participating in the nematic alignment.These numbers are comparable to the value c ≈ 0.975 obtained when c is given by the fraction of tetrahedra that touch those loops participating in the nematic alignment (see main text).

Finite temperature nematic phase transition
We now consider the full loop-resolved effective spin Hamiltonian H = l H l , where for each loop l we have (SM21) Here nn(i) are the ith site nearest neighbours of opposite type.In the above we assume that temperatures are sufficiently low that AFM order exists at the length-scale of the tetrahedron, so that the spin directions define the local orientation of the AFM structure.The factor of 1/4 takes into account that the tetrahedron free energy as calculated in Sec.SM1 C 1 involves four J ab bonds.Expanding the orientational entropy with respect to J ab results in the leading order term We assume quadrupolar order to set in, and choosing the polarization axis along z gives (such that the trace equals 1), while rotational invariance implies ⟨s xi ⟩ = ⟨s yi ⟩ = ⟨s zi ⟩ = 0. Substitution of the above into Eqn.SM24 with the first term dropped, finally gives, and the mean field loop Hamiltonian: When performing the simple gauge transformation from AFM to FM, s i → (−1) i s i , in Eqn.SM24, the corresponding free energy f HB per site of long loops (without the last mean field term) can be found within the continuum approximation 3 , valid in the limit T ≪ J, with respect to the order parameter Q z .Rewriting T ≡ A|J|τ , such that D|J|/T 2 = AQ z |J|/T = Q z /τ , we have Since ϕ[ρ] ∼ ρ 2 for small ρ, a local minimum will exist with Q z = 0.It may be shown that this minimum eventually becomes unstable at increasing temperature, however before this happens, a second minimum at finite Q z gives f MF = 0 indicating a first order transition.This happens if, for a positive Q z , one finds simultaneous solutions of f MF = 0 and df MF /dQ z = 0, or Multiplying the second equation by Q z /2 we find for ρ = Q z /τ the equation From its solution, ρ * , one obtains the order parameter at the first order transition, Figs. 3a-b display the heat-capacity and quadrupolar order parameter for a range of J/|J ab | values within region IV.Also shown is J/|J ab | = −1, which is the less constrained PHAFM and does not exhibit the nematic transition.Both observables show that T 0 increases as J/|J ab | increases, reaching a maximum around J/|J ab | = −8/7 and then decreases indicating non-monotonic behavior very close to the PHAFM boundary.The magnitude of |Q| in the sub-T 0 temperature regime also increases, indicating that a growing fraction of the sample nematically aligns.Defining Q l as the quadrupolar order parameter of the lth loop, this trend is reflected in the loop-loop quadrupolar correlation function shown in Fig. 3c.For J/|J ab | = −2 the two giant loops dominate the transition and only their Néel vectors become well aligned.However, as J/|J ab | approaches -1, smaller and smaller loops take part in the alignment and contribute to the bulk quadrupolar order parameter.This trend saturates around J/|J ab | = −8/7, where the smallest loops still remain only weakly aligned.
FIG. SM1.Breaking up the giant loops.a) Heat capacity and b) loop-loop quadrupolar correlation of an L = 8 charge-ice system with J/|J ab | = −2 in which the giant loops are decomposed into smaller winding loops.c) Visualization of its eight largest loops.

FIG. SM4 .
FIG. SM4.Density of normal mode frequencies of the transverse harmonic spin Hamiltonian.Data is shown for both aligned and randomly aligned (RLA) ground state configurations, respectively.Each loop contributes two zeromodes (not shown).In the low frequency regime, a power-law 1/ √ λ is seen, suggesting the dominance of 1D long-rangeordered AFM spin-chain behaviour.
is the quadrupolar mean-field felt by loops of length l of type a/b.Fig. SM5a displays the selfconsistent values of Q a/b z obtained upon iteration of the above, as a function of J/|J ab |, using the discrete loop distributions for the L = 8 ordered charge-ice (systems I and II) and the charge-ice realization used in Fig. 3 of the main text.The charge-ice configuration Q a/b z rapidly saturates to a maximum value, whereas for the chargeordered structure, Q a/b z indeed remains small.For the former, only a small difference is seen in the average value experienced by sites of type a and b reflecting the similar sizes of the two giant loops, whereas for the latter ordered structures the Q a/b z are identical in value.For the ordered charge-ice system II, Q a/b z grows most weakly reflecting the large number of hexagonal loops in the system.The Q a/b z and Eqn.SM18 may be used to calculate ⟨Tr FIG. SM5.Mean field description of nematic alignment.a) Self-consistent mean-field quadrupolar component for T → 0 as a function of J/|J ab | using the discrete loop distributions of the L = 8 charge-ice realization of Fig. 3 in the main text and the charge ordered systems of Fig. SM3.It is noted that values of J/|J ab | nearing -1 are outside the perturbative regime presently assumed.b) Mean-field prediction of ⟨Tr[Q l 1 Q l 2 ]⟩/(⟨|Q l 1 | 2 ⟩⟨|Q l 2 | 2 ⟩) 1/2for small loops (l1) embedded in the quadrupolar field of larger loops.Here we take l2 to be a giant loop.This is to be compared to Monte-Carlo simulations of the L = 8 charge-ice realization of Fig.3in the main text for the case of J/|J ab | = −2.
where −ϕ[ρ] (with ρ ≡ DJ/T 2 ) is the smallest eigenvalue of the quantum-mechanical hindered rotor Hamiltonian:H = − 1 2 L2 + ρ(cos 2 θ − 1/3), (SM30)L being the angular momentum operator in spherical coordinates.The ρ-independent term f 0 [T ] is immaterial for the discussion of the phase transition.It remains to minimize the mean field free energy per site,