Dynamics of the antiferromagnetic spin ice phase in pyrochlore spinels

,


I. INTRODUCTION
In frustrated magnets, no single magnetic ground state is able to satisfy all competing interactions.Such materials can show a wide range of exotic physical properties, including extensive ground-state degeneracy, fractional excitations, topological order, and other hallmarks of spin-liquid physics [1][2][3].Such behaviour is, however, usually limited to fine-tuned models: In most real magnetic materials, additional (e.g., furtherneighbour or spin-lattice) interactions and disorder tend to suppress spin liquids in favour of an ordered state or a spin glass [4].
Spinel pyrochlores of formula  2 O 4 , where the magnetic  ions form a lattice of corner-sharing tetrahedra, exhibit a variety of such mechanisms.For example, the magneto-structural order of ZnCr 2 O 4 at low temperature involves both a lowering of symmetry from cubic to tetragonal and a doubling of the unit cell, as well as a non-collinear spin arrangement [5].This complexity arises from the interplay between furtherneighbour interactions and spin-lattice coupling [6][7][8].In the Zn 2−  Cd  Cr 2 O 4 solid solution, on the other hand, the introduction of bond disorder by chemical substitution produces an apparent spin glass ground state for even small values of  [9].
Perhaps the most intriguing mechanism for relieving frustration in spinels is due to spin-lattice coupling.The simplest approach to model this interaction is to consider a bondlength-dependent Heisenberg model coupled to phonons that modulate the bond lengths independently: Integrating out the phonons then results in the bilinear-biquadratic (BLBQ) effective Hamiltonian [6,10,11] In the large- limit (a justified approximation for the  = 3/2 Cr 3+ ions), the pure Heisenberg model  = 0 is a classical spin liquid [12].Spin-lattice coupling introduces a finite  > 0, which causes the model to develop collinear nematic order at  ≃  [13].This order, however, retains residual frustration: Once all spins align collinearly along a nematic director to optimise the  term, the  term is equivalent to nearestneighbour spin ice, which is optimised by exponentially many two-up-two-down states [14].While the BLBQ model has been very successful in describing the magnetisation plateaux of several pyrochlore spinels [11,[15][16][17][18], longer-range interactions not captured by (1) cause most of them to exhibit full magnetic ordering rather than the predicted nematic spin-ice state.In particular, the BLBQ model decouples the length modulation of different bonds around the same site: in the more realistic site-phonon model [19,20], phonons mediate further-neighbour multi-spin couplings that explain such orders as the plateau phases of CdCr 2 O 4 and HgCr 2 O 4 [18].
In the past decade, much of the attention on pyrochlore spinels has shifted to the breathing pyrochlores  ′ Cr 4 O 8 , where the ordering of the  and  ′ cations cause translationinequivalent (up-and down-pointing) tetrahedra in the pyrochlore lattice to have different sizes and hence exchange couplings.High-temperature susceptibility measurements on the two best-known materials in the family, LiGaCr 4 O 8 and LiInCr 4 O 8 , indicate a ratio of Heisenberg couplings  ′ / of around 0.6 and 0.1 between the inequivalent tetrahedra, respectively [21].At low temperatures, both materials show magneto-structural ordering driven by spin-lattice coupling [20,22,23], alongside phase separation due to site disorder.
Mixing Ga and In on the  ′ site, however, quickly suppresses this ordering, with a possible gapped spin liquid on the Inrich side of the phase diagram, and a spin glass on the Garich side [24].For LiGa 0.95 In 0.05 Cr 4 O 8 (with  ′ / close to that of LiGaCr 4 O 8 ), however, neutron-diffraction and specificheat measurements indicate that the nematic state predicted by the BLBQ model ( 1) is stabilised [25].This indicates that the bond disorder introduced by the Ga-In mixing is weaker than the leading biquadratic interaction generated by spinlattice coupling, allowing nematic ordering, but it is strong enough to disrupt full ordering, stabilising instead a glassy spin-ice state [13], similar to dipolar spin ice [26].Inelasticneutron-scattering experiments on LiGa 0.95 In 0.05 Cr 4 O 8 [27] found a broad peak in the dynamical structure factor S(, ) at ℏ ≈ 5.5 meV and wave vector  ≈ 1.6 Å −1 [corresponding approximately to the (200) reciprocal lattice vector for the cubic lattice constant  0 = 8.25 Å], tentatively ascribed to excitations localised on antiferromagnetic hexagon loops.
In this paper, we explore the dynamics of LiGa 0.95 In 0.05 Cr 4 O 8 and the nematic state of breathing pyrochlores in general.Our main focus will be the classical BLBQ model, generalised to the breathing lattice: (2) where ↑, ↓ stand for up and down tetrahedra, respectively.For brevity, we will use units in which the spin magnitude |ì   | is 1 [28] and introduce  = ( +  ′ )/2,  = ( +  ′ )/2.In Section II, we discuss simulations of the semiclassical Landau-Lifshitz dynamics under (2).The numerically obtained dynamical structure factor is dominated by a broad peak at ℏ = 16, a sharp linearly dispersing mode at small , , and a weak continuum extending up to about 4.We discuss the fate of these features in the presence of Landau-Lifshitz damping and structural disorder; in Section III, we explain them in terms of linear "spin-wave" theory around the disordered spin-ice ground states of (2).We find that the 16 peak is caused by out-of-phase precession around ferromagnetic loops of spins, while the linearly dispersing feature originates in long-wave fluctuations of the nematic director.We perform the same analysis for the more complex site-phonon model as well: we find that all qualitative features of the BLBQ dynamics survive, albeit with strong effective disorder and an approximately halved effective .We discuss our findings in the context of experimental results in Section IV.

II. DYNAMICAL SIMULATIONS
We drew initial configurations from the thermal ensemble  − H using single-spin-flip Metropolis Monte Carlo on 16 × 16 × 16 pyrochlore cubic unit cells (65 536 spins).Similar to Ising spin ice [26], we expect Monte Carlo dynamics to slow down substantially in the nematic phase.Therefore, to avoid getting stuck in local minima that do not satisfy the ice rules, we used simulated annealing from a temperature well above the nematic transition [at least 2 max(,  ′ )] down to  = 0.01, well below the transition in every case, where we performed all dynamical simulations.
We then computed the time evolution of the initial spin configurations under the stochastic Landau-Lifshitz dynamical equation where ì   = −H /ì   is the effective field acting on spin ì   and ì   is a stochastic field satisfying the fluctuation-dissipation relation We note that the sign of the precession term in (3) is flipped compared to its usual presentation [29], as we take the negative gyromagnetic ratio of the electron into account through the definition of ì   .Details of the numerical method are described in Appendix A; we benchmarked the simulations by ensuring that thermodynamic properties, such as the average energy, match the Monte Carlo results within statistical error.
We first focused on the case  =  ′ ,  =  ′ = 0.1,  = 0.01.The powder-averaged dynamical correlation function (, ) is plotted in Fig. 1.Similar to the experimental powder neutron-scattering pattern of Ref. [27], we see a prominent scattering maximum at energy transfer ℏ ≈ 1.6, with the highest intensity at  ≈ 4/ 0 .In addition, we observe a linearly dispersing branch, extending from the origin to about the frequency of the intensity maximum.
We also plotted the static structure factor (), as well as (, ) integrated over two frequency windows, in Fig. 2. The static structure factor shows sharp pinch points in the pattern seen in Ref. [30] for spin ice; this is expected as the effective Ising spins in the nematic order obey the same ice rules.The same pinch points, albeit much broadened, are seen at finite  as well; below the intensity maximum at 1.6, we also see sharp circular features corresponding to the linearly dispersing mode in Fig. 1.
As shown in Fig. 3, these features of the scattering pattern are quite robust in parameter space.While the model with  ′ =  ′ = 0 is qualitatively different from the symmetric case (it is made up of disconnected tetrahedra), most features of the latter are already recovered for / ′ = 0.2 and the powder patterns at  ′ / ≥ 0.4 only differ in quantitative details.Likewise, the scattering pattern of the pure Heisenberg model  =  ′ = 0  is qualitatively different due to the lack of nematic ordering, but any finite  is enough to bring about both the linearly dispersing mode and a broad scattering maximum at 16.
The intensity maximum at 16 in our simulations is much sharper than the analogous experimental feature [27].This may either be because interactions with other dynamical degrees of freedom (e.g., magnon-magnon interactions) induce stronger damping than the Landau-Lifshitz damping term  = 0.01 used above, or because the exchange couplings ,  are disordered due to structural disorder or magnetoelastic distortions [13].To distinguish these possibilities, we performed dynamical simulations at  =  ′ ,  =  ′ = 0.1 for two additional values of  = 0.1, 0.001, as well as for a disordered model in which ,  for each bond is independently drawn from a Gaussian distribution of 10% standard deviation relative to the mean.The results are summarised in Fig. 4. The 16 feature is broadened by a similar amount both for  = 0.1 [Fig.4(a)] and on introducing disorder [Fig.4(b)]: the line shapes of the -integrated structure factor [Fig. 4(c)] appear Lorentzian and Gaussian, respectively, although the actual line shape is difficult to distinguish due to the background intensity.The linearly dispersing mode, however, behaves qualitatively differently: it does not blur noticeably on introducing disorder but becomes broad and faint beyond the point of clear detection on increasing  [see also Fig. 4(d)].Lowering  causes the low-frequency structure factor to decompose into discrete normal modes of the finite simulation box, resulting in an array of sharp peaks in Fig. 4(c).Finally, we considered dynamics under the site-phonon Hamiltonian where the inner sum runs over all pairs ,  of nearest neighbours of site  ( ,  and ,  are both counted),    =  ( ′ ) if the bond   is part of an up (down) tetrahedron, and ê  is the unit vector pointing from site  to .Since this Hamiltonian has a fully ordered ground state [19,20], we emulated the glassy nematic order of LiGa 0.95 In 0.05 Cr 4 O 8 by first preparing low-temperature states of the BLBQ Hamiltonian and annealing them under (5) before measuring dynamical correlation functions.At  =  ′ ,  =  ′ = 0.1, we obtained the powderaveraged structure factors shown in Fig. 5; see Appendix A 4 for a wider range of parameters.The general structure of the powder pattern remains unchanged and, in particular, the nematic state appears to be metastable even without quenched disorder.However, the finite-frequency peak becomes much broader, and the peak frequency is reduced substantially, from 16 to about 7.The linearly dispersing mode remains sharp, but its velocity is reduced, too.
Remarkably, we see a strong modulation of the intensity with frequency that appears to split the peak into a number of fringes.Similar, albeit weaker, fringes have already appeared in the BLBQ model [cf.Fig. 4(b)]: those are caused by the finite gap between normal modes in the linearly dispersing mode in a finite sample.To rule out this origin for the modulation seen in the site-phonon model, we performed dynamical simulations on clusters of  3 cubic unit cells for every 12 ≤  ≤ 17.
After averaging the powder pattern for the different clusters, several fringes at low frequencies (where finite-size effects are the most pronounced) indeed disappear; however, the peak remains split into three.Nevertheless, we expect that these fringes would be washed out either by stronger damping or quenched disorder, similar to the peak broadening seen in the BLBQ model in Fig. 4.

III. LINEAR SPIN-WAVE THEORY (LSWT)
Spin dynamics in the nematically ordered phase consists of small thermal fluctuations around the equilibrium configuration, a spin-ice configuration with an arbitrary Ising axis.Without loss of generality, we choose this axis to be ±  , so we can write where the   = ±1 satisfy the ice rules, and  +  ∼ √︃ / ≪ 1. Substituting this into (2) and expanding to quadratic order in  ± gives where the nonzero matrix elements    are in an ice-like arrangement of   .Likewise, substituting (6) into the dynamical equation (3) (without the stochastic fields ì   ) and expanding to linear order gives (see Appendix B) The dynamical modes of ( 9) and their frequencies are given by the eigenvalue equation where we introduce bra-ket notation for the vectors comprised of all  +  and define for convenience the diagonal matrix  with the Ising configuration   along the diagonal.For  = 0, all eigenfrequencies of (10) are real (see Appendix C), as expected for energy-conserving dynamics near a ground state.Likewise, for  > 0, all modes decay exponentially.
We diagonalised (10) for a cluster of 12 × 12 × 12 cubic unit cells (27 648 spins), for both  = 0 and 0.01.As explained in Appendix D, these eigenvalues and eigenvectors can be used to compute the dynamical structure factor (, ) within the linear-spin-wave approximation.We find excellent quantitative agreement in the -integrated structure factor (Fig. 6) as well as the powder pattern (not shown).The two curves in Fig. 6 differ in two ways.First, the low-frequency oscillations show a different pattern due to the different system sizes (and thus different low-frequency modes).Second, higher-frequency features of the LSWT spectrum are consistently shifted to slightly higher frequencies.This is due to spin-wave interactions, most of which can be accounted for in a simple mean-field picture: As the length |ì   | of spins is fixed to 1, transverse fluctuations cause ⟨   ⟩ =:  0 to shorten, which renormalises the coefficients of (9) as  ↦ →  0 ,  ↦ →  3 0 [27,28].From (D6), we estimate  0 = √︁ 1 − ⟨ +   −  ⟩ ≈ 0.9916; scaling LSWT frequencies by a factor of  3 0 indeed causes the 16 peaks of the two curves to overlap perfectly.
In summary, spin waves give a full, quantitative account of the inelastic spin dynamics, and nonlinear effects affect the spin-wave spectrum very weakly at low temperatures.In the following sections therefore, we will explain the salient features of the dynamical structure factor in terms of particular eigenmodes of the dynamical equation (9).

A. Exact eigenmodes at 16𝑄
In the simulations, we see an accumulation of eigenfrequencies near 16.More interestingly, a number (140 for the pattern of   we used) of them is equal to ±16 within numerical accuracy for  = 0. We also observe that these exact eigenmodes live exclusively on up spins (for  > 0) or down spins (for  < 0).
Since the LSWT matrix  only has on-site and nearestneighbour matrix elements, we can decompose it into terms acting on a single tetrahedron only.These terms have the form for up tetrahedra; for down tetrahedra,  →  ′ ,  →  ′ .The first two and last two rows and columns of the matrix correspond to up and down spins, respectively. ↑ has two eigenvectors with eigenvalue 8: they are orthogonal to both (1, 1, 1, 1) (i.e., they respect the ice rules) and the spin configuration of the tetrahedron.By enforcing these constraints on every tetrahedron, we can construct a number of eigenmodes of ; the corresponding eigenvalue is 8( +  ′ ) = 16 as every spin belongs to one up and one down tetrahedron.
To obtain an eigenvector of the dynamical matrix  from this construction, we need them to be eigenvectors of  as well, so they must be constrained to up ( = +1) or down ( = −1) spins in the nematic Ising configuration.On each tetrahedron, there are two configurations that obey all of these requirements: out-of-phase fluctuations of either the two up or the two down spins.We can build joint eigenstates of  and  from these by following closed loops of up or down spins and giving nearest neighbours out-of-phase fluctuations.In periodic boundary conditions, the loops always close, so the resulting fluctuation vectors are exact eigenvectors of both  and the dynamical matrix .Therefore, they are also eigenmodes of the damped dynamics (10) with complex frequency ℏ = 16(±1 − ), as we also found in exact diagonalisation.
Since each spin has precisely two neighbours of the same spin in a spin-ice configuration (one on each tetrahedron), the loops that give rise to these states are uniquely defined.By constructing them on simulated ice configurations, we found that their lengths have a very broad distribution: a few loops cover almost all spins, while the remaining ones form very small loops, often as small as a single hexagon.The exact eigenmodes on the latter resemble the "weathervane modes" proposed in Ref. [27] (inset of Fig. 7).
The exact eigenmodes described above, however, do not account for the full intensity of the 16 peak in the dynamical structure factor, or the high density of LSWT eigenmodes near this frequency.On a long loop, however, we can consider "excited loop states," in which the exact eigenmode with alternating phases is modulated with a standing wave along the loop.Locally, this pattern is very similar to the exact eigenmode, thus we expect them to be eigenmodes to a good approximation, with a frequency very close to ±16.A few numerically obtained eigenmodes of  follow this pattern closely (Fig. 7) and most of those near  = ±16 show similar features, albeit obscured somewhat by local interference between different loops.Furthermore, the fact that these modes live on loops explains the singular cusp, characteristic of one-dimensional van Hove singularities, in the structure factor.

B. Low-frequency dispersive modes
For  = 0, we found that the lowest-magnitude eigenvalues of both  and the dynamical matrix  organise themselves in approximate multiplets (Fig. 8).Their multiplicities match those of the reciprocal lattice vectors {100}, {110}, {111}, {200}, . . . of the cubic simulation box, while the eigenvalues of  and  scale with the same wave vectors as  ∝ ±| |,  ∝  2 , respectively.This indicates a linearly dispersing dynamical mode, consistent with the dynamical simulations.The corresponding eigenvectors do not show a particularly high overlap with the plane waves |⟩, but rather with the vectors At  = 0, this mode corresponds to rotating the Ising axis of the nematic order; for small , it captures long-wavelength fluctuations of the director, which we anticipate to cost little energy.
To explain these findings, we first note that both |⟩ and |⟩ are approximate eigenmodes of  for small .The first closely resembles the all-in-all-out configuration (1, 1, 1, 1) on each tetrahedron, which is an eigenvector of the singletetrahedron Hamiltonian (11) with eigenvalue 4 + 8.Since each spin belongs to one up and one down tetrahedron, these contributions add up to give In the limit  → 0, |⟩ is an exact eigenvector of : the matrix elements of  have no factors of   , so it is translation symmetric and its eigenvectors are plane waves [31] with eigenvalue () =  2 /4 +  ( 4 ).A finite  adds disordered terms to  that penalise fluctuations proportional to the spin configuration   on each tetrahedron.However, the weight of such fluctuations is only  ( 2 ) as the plane wave |⟩ is proportional to (1, 1, 1, 1) +  () on each tetrahedron.That is, even if  ≫  causes these fluctuations to be completely projected out, the resulting eigenmode of  is still |⟩ up to  ( 2 ) corrections.That is, where () ∝  2 .Now, up to  ( 2 ) corrections, the dynamical equation ( 10) can be written for a mode |⟩ = |⟩ + |⟩ as For  = 0, we obtain the eigenmodes The numerically obtained  and  plotted in Fig. 8 match (16a) closely.Since () ≪ , the modes (16b) are dominated by |⟩, but the small admixture of |⟩ is enough to yield a visible dispersing mode, as it is not diffuse in -space.
For  ≠ 0, the complex eigenfrequencies of ( 15) are Even for small , the decay rate Im  is independent of , thus coherent oscillations at the longest wavelengths are always disrupted (the square root in (17) becomes imaginary, purely decaying modes).For  = 0.1, this becomes the case for most values of , that is, the linearly dispersing modes blur completely, as seen in the dynamical simulations.By contrast, we expect that they remain robust against disorder: since they are dominated by long-wave modulations of the nematic director, they are only sensitive to coarse-grained averages of the exchange couplings , , which are affected far less by disorder.

C. Site-phonon model
Finally, we consider the site-phonon Hamiltonian (5).Expanding the Hamiltonian to quadratic order is somewhat more complicated than in the BLBQ case, and is described in detail in Appendix E. We find that, despite the further-range quartic terms in (5), the matrix  in (7) only contains nearestneighbour terms.In a spin-ice configuration, the coefficient of if the bond   is on an up tetrahedron, where  ′ , , ,  ′ are consecutive sites along a ⟨110⟩ chain (cf.Fig. 11); for a down tetrahedron,  →  ′ ,  ↔  ′ ; the diagonal terms   are determined by the constraint |⟩ = 0 imposed by spin-rotation symmetry.The first three terms of ( 18) only renormalise ,  ′ : unless  ′ ≪  or  is large compared to , this does not affect the stability of the dynamics.The fourth term in analogous to the  term of ( 8), but is halved in magnitude: this accounts for most of the reduction in the inelastic-peak frequency.
The last term depends on spins outside of the bond, so it acts as disorder on top of this renormalised BLBQ quadratic Hamiltonian.For  ≈  ′ , its magnitude is comparable to the renormalised  term, so it is expected to strongly broaden the finite-frequency peak, as indeed seen in Fig. 5.To account for the remaining discrepancy in the peak frequency, we consider a simple mean-field picture, where we replace the last two terms with In nearest-neighbour spin ice at zero temperature, the correlator of two spins in this position is ⟨    ′ ⟩ = ⟨    ′ ⟩ ≈ 0.0883; at  =  ′ , this predicts a further renormalisation of  that brings the peak down from 8 to ≈ 7.3, in good quantitative agreement with Fig. 5.
A detailed account of the splitting of the renormalised 16 peak is beyond the scope of this work.However, we speculate that due to the strong but discrete [the last term of (18) can only be ± √  ′ or 0] disorder, the eigenmodes living on long ferromagnetic chains discussed in Sec.III A break up into short segments with equal disorder terms, leading to three peaks.

IV. DISCUSSION
To summarise, dynamical simulations of the bilinearbiquadratic Hamiltonian (2) on the breathing pyrochlore lattice with  ′ / ≳ 0.3 show a dynamical structure factor made up of three components: (1) a broad inelastic peak at 16, (2) a sharp linearly dispersing mode, and (3) a broad, weakly dispersing continuum extending to about 4.We accounted for this spectrum quantitatively in terms of small fluctuations around the spin-ice-like ground states of the nematically ordered model, similar to linear spin-wave theory around conventionally ordered magnets.In particular, feature (1) is due to collective spin precession around long ferromagnetic loops, while feature (2) originates in long-wavelength fluctuations of the nematic director.We developed a similar small-fluctuation theory for the more accurate site-phonon model [19,20], which shows the same qualitative features, albeit with a renormalised dispersion relation: in particular, the position of the finite-frequency peak is renormalised down to about 7.As shown in Fig. 9, the theoretically predicted spectrum matches inelastic-neutron-scattering experiments on the nematically ordered spinel LiGa 0.95 In 0.05 Cr 4 O 8 [27].Comparing the positions of the inelastic peaks, we estimate  ≈ 0.35 meV assuming the BLBQ model and  ≈ 0.75 meV assuming the site-phonon model.Taking the length of the  = 3/2 Cr 3+ moments into account [28], the latter corresponds to  ≈ 0.22 meV.Estimates of  and  for LiGa 0.95 In 0.05 Cr 4 O 8 in the literature vary widely between  ≈ 40 K [21] to  ≈ 80 K [32,33].For the former, our estimate yields  ≈ 0.07, comparable to typical figures in the literature.On the other hand, the parameters proposed in Ref. [32] yield  ≈ 1.0 meV, almost an order of magnitude too high; furthermore, their estimate of  ′ / ≈ 0.04 leads to low-temperature dynamics desribed by isolated tetrahedra, qualitatively different to the experiment.Earlier estimates of  ′ / ≈ 0.6 [21,33] appear more consistent with the neutronscattering results.
The idealised BLBQ spectrum shown in, e.g., Fig. 1 differs from the experimental results in two key ways: (1) the experimental inelastic peak is far broader than the simulated 16 peak; (2) the linearly dispersing mode is absent in the experiments.We found that the former can be explained either by finite excitation lifetime due to dynamical processes (modelled using strong Landau-Lifshitz damping ) or by static disorder in the Hamiltonian, cf.Figs. 4 and 9.The experimentally observed Gaussian shape of the peak [27], however, agrees better with the latter scenario.The additional couplings of the site-phonon model can also be regarded as strong disorder on top of the BLBQ dynamics, which indeed broaden the peak to a similar extent as seen in the experiment.We believe that the additional structure of this peak would also be washed out by either static disorder or dynamical damping.In future work, it will be interesting to extend inelastic neutron-scattering measurements to smaller values of (, ), where the fate of the linearly dispersing mode could be studied directly.
Our findings are also potentially relevant to spinel materials beyond LiGa 0.95 In 0.05 Cr 4 O 8 .For example, the solid solution Zn 2− Cd  Cr 2 O 4 with  = 0.05 shows a similar phenomenology to LiGa 0.95 In 0.05 Cr 4 O 8 in both magnetic susceptibility and inelastic neutron scattering [9], raising the possibility that the nematic state is generic to chromium spinels with light disorder.The inelastic excitation energy 4 meV measured for this system [9] implies  BLBQ ≈ 0.25 meV and  sp ≈ 0.55 meV in the BLBQ and site-phonon models, respectively.Given the value of  = 3.5 meV for the end-member ZnCr 2 O 4 [34], the latter implies  ≈ 0.04, consistent with the value  ≈ 0.02 suggested by high-field magnetisation measurements on ZnCr 2 O 4 [16].
We finally note that the order-by-disorder-induced nematic phase of the classical kagome Heisenberg model [35] also exhibits sharp linearly dispersing dynamical modes on top of a partially ordered nematic background.The similarity of the ice-like ground states, as well as the dynamics, of these two systems raises the tantalising possibility of a deeper analogy between them.In future work, therefore, it will be interesting to study the long-time relaxation dynamics of the nematic order in our models: In the kagome case, this dynamics is governed by qualitatively different processes from the LSWTlike precession dynamics studied in this work.The relaxation dynamics may also be affected in exotic ways by kinematic constraints, possibly analogous to the fractal dynamics recently uncovered in pyrochlore spin ice [36].colour maps based on Ref. [38].Computing resources were provided by STFC Scientific Computing Department's SCARF cluster. A. Sz.gratefully acknowledges the ISIS Neutron and Muon Source and the Oxford-ShanghaiTech collaboration for support of the Keeley-Rutherford fellowship at Wadham College, Oxford.For the purpose of open access, the authors have applied a Creative Commons Attribution (CC-BY) licence to any author accepted manuscript version arising.
Appendix A: Details of the dynamical simulations

Monte Carlo sampling
We used Metropolis Monte Carlo with single-spin updates to draw spin configurations from the thermal ensemble of the Hamiltonian (2).In each step, the proposed spin update was constructed as where each component of ì   is drawn independently from a Gaussian distribution of variance  2 = /( + 6), chosen to match the thermal fluctuations under the on-site part of the LSWT Hamiltonian (8).We note that, since there are no interactions between spins on the same sublattice, the proposalacceptance cycle can be performed in parallel for all spins of the same sublattice, allowing for efficient vectorisation.
This protocol results in a temperature-independent acceptance rate at low temperatures, indicating that the thermal fluctuations around the ordered state are captured well.However, the nematically ordered moments become frozen sufficiently below the ordering temperature, resulting in rather noisy static structure factors from a single run.Therefore, to obtain the static structure factor shown in Fig. 2(a), we initialised the Monte Carlo with 256 independent spin-ice configurations obtained from a variant of the codes in Refs.[39,40].

Stochastic dynamics
To solve the Landau-Lifshitz dynamical equation (3), we implemented the semi-implicit integrator SIB proposed in Ref. [41].This algorithm achieves comparable accuracy to fully symplectic solvers (such as the implicit midpoint method) at a fraction of the computational cost by exploiting the sparseness of the dynamical equation.We performed 65 536 time steps of size Δ = ℏ/(16) for a total simulation time  = 4096ℏ/, resulting in a frequency resolution Δ = 2/ ≈ 1.53 × 10 −3 /ℏ.We only saved the spin configuration after every fourth step, as this still allowed us to resolve the full dynamical spectrum.
In stochastic differential equation solvers, the noise term ì   of (3) is implemented using a noise vector ì   , whose components should be unit Gaussian random numbers.Due to the implicitness of the solver, however, using unbounded ì   can lead to instabilities.Ref. [41] proposes to simply apply a cutoff to Gaussian ì   components: we found that this results in strong numerical damping and equilibrium energies well below that obtained for the same temperature from Monte Carlo at any temperature, time step size, or value of .By contrast, drawing ì   uniformly from the surface of a sphere of radius √ 3 (such that the standard deviation of each component is 1) resulted in energies that match the Monte Carlo results within statistical error.We believe that matching the (co)variances of the ideal Gaussian noise in any projection (not only along the Cartesian axes) is crucial for this.
For most parameter values, we ran a single dynamical simulation, as the dynamical fluctuations appear to remain selfaveraging even in a frozen spin-ice background.For the parameters  =  ′ ,  =  ′ = 0.1,  = 0.01 used in Figs. 1, 2 and 6, we averaged four independent runs to improve statistics.

Powder averaging
To compute powder averages of -dependent quantities, we broadened every -point obtained from FFT with a Gaussian of standard deviation   = √ 2/ 0 , where  is the number of cubic unit cells along each Cartesian direction ( 0 is the linear size of the simulation box) and integrated the result over bins of width Δ: We can think of this as spreading out the discrete -points into three-dimensional Gaussians in reciprocal space and averaging the result over spherical shells; the denominator 4 2 Δ is the volume of such a shell.The width   of the Gaussian was chosen such that the effective volume ( √ 2  ) 3 taken up by them in reciprocal space match the volume around each allowed -point, (2/ 0 ) 3 .Both this choice and integrating over bins of  reduce spurious fluctuations due to the discrete -points available on the finite-size system, allowing us to use the relatively narrow bin width Δ = / 0 .

Dynamical structure factor of the site-phonon model
We performed dynamical simulations of the site-phonon model (5) for all the parameter sets used in Fig. 3.The resulting dynamical structure factors are shown in Fig. 10.Similar to the BLBQ model, we find that the finite-frequency peak and the linearly dispersing mode form clearly for all  ′ / ≳ 0.3, and all quantitative features are essentially the same for all  ′ ≳ 0.5.The position of the main inelastic peak also appears to remain proportional to . spin-wave dynamics from having any exponentially decaying or exploding modes, thus all eigenvalues of the dynamical matrix  must be real.To prove that this is the case, we multiply (10) with  from the left: Now, both  and  are Hermitian matrices, and  is positive definite [43]: this implies that the eigenvalues   are all real [44].As both matrices are also real, the eigenvectors are real, too.We can extend these arguments to show that the eigenfrequencies of the  > 0 dynamics have Im  ≤ 0, that is, they all decay.We multiply (10) from the left by  1/2 to get Multiplying on the left by ⟨ r |, we find that the eigenfrequency is given by the Rayleigh quotient Every matrix in the second form of (C3) is Hermitian, so Im  is entirely due to the second term.Furthermore, both  and  2 are positive definite [43], so ⟨  | 2 |  ⟩/⟨  ||  ⟩ > 0 for any |  ⟩, thus Im   < 0 for all but the zero mode.
The matrix  1/2 (−) 1/2 in (C2) is symmetric.Complex symmetric matrices are generally not normal, so their left and right eigenvectors are different.However, taking the transpose (not conjugate transpose) of (C2) gives that is, ⟨ r *  | is a left eigenvector corresponding to the same eigenvalue as | r ⟩.The resolvent of  1/2 ( − ) 1/2 can therefore be written as To incorporate the stochastic fluctuation terms of (3) into linear spin-wave theory, we note that the leading-order components of ì   × ì   and ì   × ì   × ì   are those involving the -component of ì   .Therefore, where | + ⟩ is the vector of    +    .The dynamical structure factor in real space is given by the thermal average of | + ()⟩⟨ + ()|.To perform the average, we note that In the second line, we use that  is a diagonal matrix with ±1 as entries.The fourth line uses the identity In the last two lines, we substitute the fluctuation-dissipation relation (4) and the spectral decomposition (C7), making the normalisation (C6) explicit.The two terms in the last two lines are manifestly complex conjugate symmetric matrices, so S() is real and symmetric, as expected from its definition.Eq. (D3) assumes that the eigenvectors |  ⟩ satisfy the orthogonality condition (C6).For degenerate modes (such as the exact 16 modes), the eigenvectors returned by nonhermitian eigensolvers do not satisfy any such relation, so blindly applying (D3) leads to incorrect results.For the results shown in Fig. 6, we added very weak (Δ/ = 10 −6 ) bond disorder to lift all degeneracies without perceptibly changing S().
The -integrated and -resolved dynamical structure factors can be expressed from (D3) as

FIG. 1 .
FIG. 1. Simulated powder-averaged dynamic structure factor for  =  ′ ,  =  ′ = 0.1,  = 0.01.The pattern is dominated by a sharp maximum at ℏ ≈ 1.6 and a linearly dispersing mode that connects this maximum to the origin.

FIG. 2 ./ = 0. 20 FIG. 3 .
FIG.2.Static and frequency-window-integrated dynamical structure factors in the (ℎℎℓ) plane for the same parameters as in Fig.1.The frequency windows used in (b, c) are highlighted with lighter background in Fig.1.The "arbitrary units" are consistent across the three plots: the static structure factor, dominated by the  = 0 component, is much larger than dynamical correlations, consistent with (nematic) ordering.

FIG. 4 .
FIG. 4. (a, b) Simulated powder-averaged dynamic structure factor for the same parameters as Fig. 1 but with  = 0.1 (a) or 10% Gaussian disorder in  and  (b).(c) -integrated structure factors for the disordered model and three values of .(d) Cuts of the powder pattern for the disordered model and two values of  at  = /(2 0 ) [dashed green lines in panels (a, b)].The "arbitrary units" are not consistent between panels (c, d).

FIG. 5 .
FIG.5.Simulated powder-averaged (left) and -integrated (right) dynamic structure factor for the same parameters as Fig.1in the site-phonon model(5).The powder pattern is averaged from simulations of clusters between 12 3 and 17 3 cubic unit cells; the integrated structure factor is shown for each cluster as well as the average.

FIG. 7 .
FIG. 7.An eigenvector |⟩ of the  = 0 dynamical matrix with frequency ℏ/ = −1.6 + 3.01 × 10 −8 , restricted on the longest closed loop of down spins, where 93% of the statistical weight falls.Inset: illustration of the exact 16 eigenmode on the shortest possible ferromagnetic loop (blue atoms).Magnetic moments around the loop (blue arrows) precess around their equilibrium Ising direction (cones); the fluctuations of nearest neighbours are out of phase.

FIG. 9 .
FIG. 9. (a) Simulated inelastic neutron-scattering pattern for 10% Gaussian disorder in  and  in the BLBQ model [same parameters as Fig. 4(b)].(b) Inelastic neutron-scattering intensity of LiGa 0.95 In 0.05 Cr 4 O 8 at  = 5.2 K, measured with 16 meV incident neutron energy[27].(c) Simulated inelastic-neutron scattering pattern for the site-phonon model (same parameters as Fig.5).The data in panels (a,c) are multiplied with the Cr 3+ magnetic form factor, and the wave vector scaled by 2/ 0 (with lattice parameter  0 = 8.25 Å), to aid comparison.Frequency ranges are chosen to approximately match the peak positions.