Giant Magnetoelastic-Coupling Driven Spin-Lattice Liquid State in Molybdate Pyrochlores

TL;DR

This paper introduces a novel spin-lattice liquid state in molybdate pyrochlores driven by strong magnetoelastic coupling, explaining experimental observations in Y$_2$Mo$_2$O$_7$ and expanding understanding of disordered magnetic systems.

Contribution

It proposes the existence of a spin-lattice liquid state in molybdate pyrochlores caused by giant magnetoelastic effects, supported by microscopic analysis.

Findings

Identification of a spin-lattice liquid state in molybdate pyrochlores.

Explanation of experimental features in Y$_2$Mo$_2$O$_7$.

Demonstration of strong spin-lattice coupling effects.

Abstract

We propose the idea of a spin-lattice liquid, in which spin and lattice degrees of freedom are strongly coupled and remain disordered and fluctuating down to low temperatures. We show that such a state arises naturally from a microscopic analysis of a class of molybdate pyrochlore compounds, and is driven by a giant magnetoelastic effect. Finally, we argue that this could explain some of the experimental features of Y$_2$Mo$_2$O$_7$.

Figures

Equations

$${\mathcal{H}}_{\mathsf{sl}}=\sum _{⟨ij⟩}\left[\left(J_{\mathsf{is}}({\alpha }_{\mathsf{av}})-g\delta {\alpha }_{ij}\right){\sigma }_i{\sigma }_j+K\delta {\alpha }_{ij}^2\right],$$

$${\mathcal{H}}_{\mathsf{loc}}={\Delta }_{\mathsf{trig}}\left[\frac{2}{3}-(L_{\mathsf{eff}}^{{\mathsf{z}}^{\prime }}{)}^2\right]-{\lambda }_{\mathsf{SO}}\mathbf{S}\cdot {\mathbf{L}}_{\mathsf{eff}},$$

$$R_0^T$$

$$R_2^T$$

$$\displaystyle P_{\sf x}^{\pm}=\pm V_{pd\pi}\left(\begin{array}[]{ccc}0&0&0\\ 0&0&1\\ 0&1&0\\ \end{array}\right),\quad P_{\sf y}^{\pm}=\pm V_{pd\pi}\left(\begin{array}[]{ccc}0&0&1\\ 0&0&0\\ 1&0&0\\ \end{array}\right),\quad P_{\sf z}^{\pm}=\pm V_{pd\pi}\left(\begin{array}[]{ccc}0&1&0\\ 1&0&0\\ 0&0&0\\ \end{array}\right),$$

$${\Lambda }_{ij}=t{P}_{ij}{R}_i{R}_j^T{P}_{ji},$$

$$\displaystyle\Lambda_{13}=\Lambda_{20}=t\left(\begin{array}[]{ccc}0&0&0\\ 0&-\frac{1}{9}&\frac{8}{9}\\ 0&\frac{8}{9}&-\frac{1}{9}\\ \end{array}\right),\quad\Lambda_{10}=\Lambda_{32}=t\left(\begin{array}[]{ccc}-\frac{1}{9}&0&\frac{8}{9}\\ 0&0&0\\ \frac{8}{9}&0&-\frac{1}{9}\\ \end{array}\right),\quad\Lambda_{03}=\Lambda_{12}=t\left(\begin{array}[]{ccc}-\frac{1}{9}&\frac{8}{9}&0\\ \frac{8}{9}&-\frac{1}{9}&0\\ 0&0&0\\ \end{array}\right),$$

$${\mathcal{H}}_{\text{hop}}=\sum _{⟨ij⟩}\sum _{\alpha ,\beta ,s}{\alpha }_{is}^{†}{\Lambda }_{ij}^{\alpha \beta }{\beta }_{js}+h.c.,$$

$${\mathcal{H}}_{\mathsf{site}}=$$

$$+J_H\left[\sum _{\alpha }{n}_{\alpha \uparrow }{n}_{\alpha \downarrow }+\sum _{\alpha <\beta }\left({\alpha }_{\uparrow }^{†}{\alpha }_{\downarrow }^{†}{\beta }_{\downarrow }{\beta }_{\uparrow }+{\beta }_{\uparrow }^{†}{\beta }_{\downarrow }^{†}{\alpha }_{\downarrow }{\alpha }_{\uparrow }\right)\right],$$

$${\mathcal{H}}_{\text{SO}}=\frac{t^2}{U}\sum _{\gamma }\sum _{⟨ij⟩\parallel \gamma }\left[\left({\text{S}}_i\cdot {\text{S}}_j+1\right){\hat{J}}_{ij}^{(\gamma )}+{\hat{K}}_{ij}^{(\gamma )}\right],$$

$${\hat{J}}_{ij}^{(\gamma )}(\eta )=(1+2\eta R){\hat{A}}_{ij}^{(\gamma )}-\eta r{\hat{B}}_{ij}^{(\gamma )}-\eta R{\hat{C}}_{ij}^{(\gamma )},{\hat{K}}_{ij}^{(\gamma )}(\eta )=2\eta R{\hat{A}}_{ij}^{(\gamma )}+2\eta r{\hat{B}}_{ij}^{(\gamma )}-(1+\eta R){\hat{C}}_{ij}^{(\gamma )},$$

$${\hat{A}}_{ij}^{(\gamma )}$$

$${-\frac{8}{81}\left(n_i({\tau }_j^{+}+{\tau }_j^{-})+({\tau }_i^{+}+{\tau }_i^{-})n_j\right)]}^{(\gamma )}$$

$${\hat{B}}_{ij}^{(\gamma )}$$

$${-\frac{8}{81}\left(n_i({\tau }_j^{+}+{\tau }_j^{-})+({\tau }_i^{+}+{\tau }_i^{-})n_j\right)]}^{(\gamma )}$$

$${\hat{C}}_{ij}^{(\gamma )}$$

$$\displaystyle P_{\sf z}^{+}=\left(\begin{array}[]{ccc}\sqrt{3}lmnV_{pd\sigma}-2lmnV_{pd\pi}&\sqrt{3}l^{2}nV_{pd\sigma}+n(1-2l^{2})V_{pd\pi}&\sqrt{3}l^{2}mV_{pd\sigma}+m(1-2l^{2})V_{pd\pi}\\ \sqrt{3}m^{2}nV_{pd\sigma}+n(1-2m^{2})V_{pd\pi}&\sqrt{3}lmnV_{pd\sigma}-2lmnV_{pd\pi}&\sqrt{3}m^{2}lV_{pd\sigma}+l(1-2m^{2})V_{pd\pi}\\ \sqrt{3}n^{2}mV_{pd\sigma}+m(1-2n^{2})V_{pd\pi}&\sqrt{3}n^{2}lV_{pd\sigma}+l(1-2n^{2})V_{pd\pi}&\sqrt{3}lmnV_{pd\sigma}-2lmnV_{pd\pi}\\ \end{array}\right),$$

$${\text{r}}_{\text{O}}=\left(\frac{1}{\sqrt{2}}\sin \left[\frac{{\alpha }_{\mathsf{cub}}-\alpha }{2}\right],\frac{1}{\sqrt{2}}\sin \left[\frac{{\alpha }_{\mathsf{cub}}-\alpha }{2}\right],-\cos \left[\frac{{\alpha }_{\mathsf{cub}}-\alpha }{2}\right]\right),$$

$$\displaystyle\Lambda_{03}(127^{\circ})=t\left(\begin{array}[]{ccc}-0.63&0.36&0.0025\\ 0.36&-0.63&0.0025\\ 0.0025&0.0025&0.0029\\ \end{array}\right).$$

$${\mathcal{H}}_{\mathsf{is}}=J_{\mathsf{is}}(\alpha ,\eta )\sum _{⟨ij⟩}{\sigma }_i{\sigma }_j,$$

$$J_{\mathsf{is}}(\alpha ,\eta )=$$

$$-\left(-\frac{2R(\eta )}{9}+\frac{r(\eta )}{3}+\frac{5}{9}\right)\frac{|t_{a_{1g}e}(\alpha ){|}^2}{U},$$

$$\displaystyle\Lambda^{\textsf{trig}}_{ij}(\alpha)=\left(\begin{array}[]{ccc}t_{a_{1g}}&t_{a_{1g}e}&t_{a_{1g}e}^{*}\\ t_{a_{1g}e}^{*}&t_{d}&t_{od}\\ t_{a_{1g}e}&t_{od}^{*}&t_{d}\\ \end{array}\right)=\textbf{U}\cdot\Lambda_{ij}\cdot\textbf{U}^{\dagger},\quad\textbf{U}=\frac{1}{\sqrt{3}}\left(\begin{array}[]{ccc}1&1&1\\ 1&\omega&\omega^{*}\\ 1&\omega^{*}&\omega\\ \end{array}\right),\quad\omega=e^{-\frac{2i\pi}{3}},$$

$${\mathcal{O}}_{\mathsf{loop}}=\sum _{\mathsf{loops}}\left|\sum _{i\in \mathsf{loop}}{\sigma }_i{\sigma }_i^{\mathsf{loop},\mathsf{gs}}\right|,$$

$${\mathcal{O}}_{\mathsf{mis}}=\frac{N_{4+4m}-N_{2+4m}}{N},$$

$${\mathcal{S}}_{\mathsf{sp}}(\mathbf{q})=\sum _{\alpha }⟨S^{\alpha }(\mathbf{q})S^{\alpha }(-\mathbf{q})⟩,$$

$${\mathcal{S}}_{\mathsf{lat}}(\mathbf{q})=⟨{\left|b_{\mathsf{Mo}}{e}^{i\mathbf{q}.{\mathbf{r}}_{\mathsf{Mo}}}+b_{\mathsf{O}}{e}^{i\mathbf{q}.{\mathbf{r}}_{\mathsf{O}}}\right|}^2⟩,$$

Full text

Frustration is at the heart of the search for unusual states of matter, and can promote spin-liquid behaviour, in which the spins remain disordered and fluctuating down to low temperature Lacroix et al. (2011); Balents (2010); Savary and Balents (2016). However, the interaction of spins with lattice degrees of freedom is typically expected to result in long-range order via an “order-by-disorder”-type mechanism Villain, J. et al. (1980); Shender (1982). For example, magnetoelastic coupling can drive a structural distortion that selects a small subset of the otherwise extensively degenerate spin configurations Yamashita and Ueda (2000); Tchernyshyov et al. (2002a, b); Penc et al. (2004); Di Matteo et al. (2004, 2005); Bergman et al. (2006); Shannon et al. (2010). Nevertheless, there is the intriguing possibility that spin and lattice degrees of freedom can become tightly coupled, but, rather than ordering, remain fluctuating at low temperature, and we dub the resultant state a spin-lattice liquid.

In order to identify candidate materials for such a spin-lattice liquid state, there are a number of characteristic features that are required. These include strong magnetoelastic coupling, an absence of spin ordering at temperatures well below the interaction energy and evidence for incoherently disordered lattice degrees of freedom at low temperature. All of these characteristics have been experimentally observed in the pyrochlore molybdates, R2Mo2O7 (*R$=$*Y,Tb,Dy,Ho,Er,Tm,Yb,Lu Clark et al. (2013); Shinaoka et al. (2018)), of which Y2Mo2O7 is the best studied example Greedan et al. (1986); Reimers et al. (1988); Raju et al. (1992); Dunsiger et al. (1996); Gingras et al. (1997); Gardner et al. (1999); Booth et al. (2000); Keren and Gardner (2001); Sagi et al. (2005); Greedan et al. (2009); Ofer et al. (2010); Silverstein et al. (2014); Thygesen et al. (2017). However, many aspects of their low-temperature behaviour remain mysterious, despite more than three decades of study.

In this Letter we analyse the Mo pyrochlores and show that at intermediate temperatures they may support a state in which spin and lattice degrees of freedom become strongly coupled and explore an extensive manifold of degenerate configurations. This state, which we refer to as a spin-lattice liquid, both serves as a first example of a class of states, as well as providing an explanation for the unusual behaviour of the R2Mo2O7 materials. We use microscopic analysis to show that spin-orbital and lattice degrees of freedom would individually be expected to show behaviour similar to the classical spin-liquid “spin ice” Bramwell and Gingras (2001). However, due to strong magnetoelastic coupling, a “spin-lattice” liquid phase emerges at low temperature. Since the magnetoelastic coupling plays a decisive and nonperturbative role in selecting the low-energy configurations, we refer to it as a “giant” magnetoelastic effect.

The R2Mo2O7 materials crystallise into the well-known Fd$\bar{3}$m space group, where Mo4+ ions form a pyrochlore lattice of corner-sharing tetrahedra and are surrounded by oxygen octahedra Reimers et al. (1988) (see Fig. 1a,b). At low temperature ($T_{\sf f}\approx 20K$) a spin-glass transition is observed Greedan et al. (1986); Gingras et al. (1997), whose origin remains an open question Shinaoka et al. (2018). However, this occurs well below the Curie-Weiss temperature ($\theta_{\sf cw}=-200K$ for Y2Mo2O7 Gardner et al. (2004)), which approximately measures the strength of magnetic interactions. As such there is a wide temperature range, $T_{\sf f}<T<\theta_{\sf cw}$ in which spin correlations are expected to be significant, but no spin ordering is observed.

At the same time, extensive studies of the crystal structure have shown that there are significant low-temperature local distortions away from the average Fd$\bar{3}$m structure that survive up to temperatures in excess of $\theta_{\sf cw}$, but no structural phase transition Booth et al. (2000); Keren and Gardner (2001); Sagi et al. (2005); Ofer et al. (2010); Greedan et al. (2009); Thygesen et al. (2017). Recent pair-distribution-function analysis has shown that for each tetrahedron two Mo4+ ions move towards and two away from the tetrahedral centre, forming a highly-degenerate “2-in-2-out” distortion pattern Thygesen et al. (2017) (see Fig. 1c). There is a concomitant distortion of the O octahedra, such that both the Mo-O bond length and the local crystal field remains almost invariant. However, the nearest-neighbour Mo-Mo distances and the Mo-O-Mo bond angles change significantly.

Magnetoelastic coupling provides a link between the spin and lattice degrees of freedom, and there is evidence from NMR and $\mu$sr measurements that it is significant in Y2Mo2O7 Keren and Gardner (2001); Sagi et al. (2005); Ofer et al. (2010). In consequence, the experimental prerequisites appear to be in place to think the R2Mo2O7 materials may realise a spin-lattice liquid state. In order to explore whether this is the case, we use microscopic arguments to understand the interplay between spin, orbital and lattice, and derive a simple effective model in an attempt to capture the essence of their low-temperature behaviour.

Local electronic states: The starting point is to consider the local physics of the Mo4+ ions. Mo4+ has two electrons in the $4d$-shell, and these are localised by an on-site Coulomb repulsion, $U$, that is large compared to the bandwidth Shinaoka et al. (2018). The local cubic crystal field due to the surrounding O2- octahedron (see Fig. 1a,b) splits the $4d$ energy levels into a low-energy $t_{2g}$ manifold of xy, yz and zx orbitals that point between the O ions, and a high-energy $e_{g}$ manifold of orbitals pointing towards O ions Fazekas (1999). The $t_{2g}$ manifold is split by Hund’s coupling ($J_{\sf H}\approx 0.5$eV Shinaoka et al. (2013)), which favours parallel spins, and therefore selects a 9-fold degenerate set of low-energy states, labelled by spin $S=1$ and effective orbital angular momentum $L_{\sf eff}=1$ Abragam and Bleaney (1970). This is further split by a non-cubic component of the crystal field with trigonal symmetry that results from a compression of O octahedra along the local $z^{\prime}$-axes pointing into/out of Mo tetrahedra (see Fig. 1b) and spin-orbit coupling. The characteristic energy scales are repectively $\Delta_{\sf trig}\approx 160$meV and $\lambda_{\sf SO}\approx 40$meV Shinaoka et al. (2013). The result is a low-energy doublet, labelled as $J^{z^{\prime}}_{\sf eff}=\pm 2$, where ${\bf J}_{\sf eff}={\bf S}+{\bf L}_{\sf eff}$, that is separated from the first excited state by an energy gap $\Delta E=\lambda_{\sf SO}[1-\lambda_{\sf SO}/\Delta_{\sf trig}+\dots]\approx 30$meV$\approx 300$K ($k_{B}=1$) sup .

Exchange interaction: An effective low-energy model can be derived by considering superexchange interactions between neighbouring $J^{z^{\prime}}_{\sf eff}=\pm 2$ doublets. This is different from the approach taken by a number of other theories, which assume an isotropic $S=1$ spin on every Mo site Saunders and Chalker (2007); Andreanov et al. (2010); Shinaoka et al. (2014). The geometry is simplest in the case of idealised undistorted oxygen octahedra, where local cubic axes can be aligned with the O octahedra, as shown in Fig. 1d. Superexchange occurs via a single intermediate oxygen with Mo-O-Mo bond angle $\alpha=2\arctan 2\sqrt{2}\approx 141^{\circ}$. Within the local cubic axes the dominant hopping channel is off-diagonal, bond-dependent and follows paths of the type $d_{1}^{\sf yz}$-$(p_{1}^{\sf y}$:$p_{2}^{\sf x})$-$d_{2}^{\sf zx}$, where $p_{1}^{\sf y}$ and $p_{2}^{\sf x}$ denote orbitals on the same oxygen ion but in the local axes of the two different Mo ions (1 and 2 in Fig. 1d). It can be seen that the $d$-$p$ bonds involve a lateral overlap of orbitals, and this is known as $\pi$-bonding Harrison (1989).

We work in the Mott insulator picture, and for two electrons in the $t_{2g}$ manifold, the derivation of a spin-orbital Hamiltonian from a multiorbital Hubbard model is a standard procedure sup ; Khaliullin et al. (2001); Horsch et al. (2003); Khaliullin (2005). In order to arrive at an effective low-temperature model, it is required to project the full spin-orbital Hamiltonian for an $S=1$, $L_{\sf eff}=1$ configuration onto the $J^{z^{\prime}}_{\sf eff}=\pm 2$ doublet. Since no matrix elements connect $J^{z^{\prime}}_{\sf eff}=2$ to $J^{z^{\prime}}_{\sf eff}=-2$, the effective Hamiltonian takes the form of a classical Ising model $\mathcal{H}_{\sf is}=J_{\sf is}\sum_{\langle ij\rangle}\sigma_{i}\sigma_{j}$ sup , where $\sigma=\pm 1$ labels the doublet states, and can be thought of as an Ising “spin” pointing into or out of tetrahedra (see Fig. 1b). $J_{\sf is}>0$ favours antiferromagnetic alignment in the local trigonal axes, which corresponds to ferromagnet exchange in crystallographic axes, and low-energy states consist of the extensively-degenerate spin-ice configurations with a 2-in-2-out arrangement. $J_{\sf is}<0$ favours ferromagnetic alignment in the local trigonal axes, corresponding to antiferromagnet exchange in the crystallographic axes and therefore unfrustrated all in/out arrangements on tetrahedra. In the case of undistorted cubic octahedra we find $J_{\sf is}<0$ for any value of $J_{\sf H}/U$.

In reality the oxygen octahedra are trigonally distorted, and in Y2Mo2O7 the Mo-O-Mo bond angle in the average Fd$\bar{3}$m structure is $\alpha=\alpha_{\sf av}=127^{\circ}$ Reimers et al. (1988); Greedan et al. (2009); Thygesen et al. (2017). The displacement of the O relative to the undistorted case opens up a $\sigma$-type hopping channel, where orbitals directly overlap one another Harrison (1989), and this is exemplified by the path $d_{1}^{\sf yz}$-$(p_{1}^{\sf z}$:$p_{2}^{\sf z})$-$d_{2}^{\sf yz}$ shown in Fig. 1e. The $\sigma$ channel becomes comparable to the original $\pi$-bonding channel for a relatively small trigonal distortion, and as a result $J_{\sf is}>0$ for $\alpha=127^{\circ}$, favouring a 2-in-2-out set of low-energy configurations sup .

Due to the competition between the $\pi$ and $\sigma$ channels, the value of $J_{\sf is}$ is extremely sensitive to the bond angle, as shown in Fig. 2a. Serendipitously, the bond angle at which the exchange interaction cancels is close to $\alpha=\alpha_{\sf av}=127^{\circ}$, both in our simple analytical calculations and also in more involved band structure calculations Shinaoka et al. (2013). Deviations of the Mo-O-Mo bond angle from its average value thus lead to very large relative changes in $J_{\sf is}$.

In addition to the trigonal distortion, it is also necessary to take into account the 2-in-2-out distortions of the Mo ions, and concomitant movement of the O octahedra, which leaves the Mo-O bond length and local crystal field almost invariant, but significantly changes the Mo-Mo distances and the Mo-O-Mo bond angles Thygesen et al. (2017). The result is three classes of Mo-Mo bonds with significantly different exchange couplings (see Fig. 1c). Combining the pair-distribution-function analysis with the calculations yields on each tetrahedron one short bond with $\alpha<\alpha_{\sf av}$ and $J_{\sf is}>0$, four medium bonds with $\alpha\approx\alpha_{\sf av}$ and $J_{\sf is}>0$ and one long bond with $\alpha>\alpha_{\sf av}$ and $J_{\sf is}<0$. This conclusion remains robust, even if the bond angle deviations are smaller than the $\sim 11^{\circ}$ reported in Thygesen et al. (2017) and shown in Fig. 2b. The magnitude of $J_{\sf is}$ is significantly larger on the short and long bonds than the medium bonds, resulting in a very strong magnetoelastic coupling effect.

Spin-lattice model: Next we consider the minimal effective model that captures the essence of the microscopic analysis. Mo lattice sites are allowed to displace in a 2-in-2-out pattern (Fig. 1c), corresponding to the triplet normal mode of a tetrahedron Yamashita and Ueda (2000); Tchernyshyov et al. (2002b); Di Matteo et al. (2005). The resultant bond angles are $\alpha=\alpha_{\sf av}+\delta\alpha$ on the long bond, $\alpha=\alpha_{\sf av}$ on the medium bonds and $\alpha=\alpha_{\sf av}-\delta\alpha$ on the short bond, with $\delta\alpha\geq 0$ (see Fig. 2b). For simplicity, $\delta\alpha$ is allowed to vary globally but not locally. Spins are constrained to point into or out of tetrahedra and described by $\sigma=\pm 1$. The resultant spin-lattice Hamiltonian is,

[TABLE]

where $J_{\sf is}(\alpha_{\sf av})$ is the superexchange interaction for undistorted bonds, $g=\partial J_{\sf is}/\partial\alpha|_{\alpha=\alpha_{\sf av}}$ describes at first order the angular dependence of the interaction, $\delta\alpha_{ij}\in\{-\delta\alpha,0,\delta\alpha\}$ and $K$ is the elastic energy cost of distorting the lattice.

The nature of the ground state can be understood starting from a single tetrahedron, and depends on the ratio $g^{2}/(4KJ_{\sf is}(\alpha_{\sf av}))$. For $g^{2}/(4K)<J_{\sf is}(\alpha_{\sf av})$ it is energetically favourable to have $\delta\alpha(T=0)=\delta\alpha_{\sf gs}=0$ combined with a 2-in-2-out spin configuration on each tetrahedra, and as a result a spin ice forms at low $T$. For $g^{2}/(4K)>J_{\sf is}(\alpha_{\sf av})$ the lowest energy is reached with $\delta\alpha_{\sf gs}=g/(2K)$, and writing $J_{\sf is}(\alpha)=J_{\sf is}(\alpha_{\sf av})\pm\delta J_{\sf is}$, one finds $\delta J_{\sf is}>2J_{\sf is}(\alpha_{\sf av})$ on long and short bonds. For a given lattice distortion a tetrahedron has four degenerate low-energy 3-1 spin configurations, where 3 spins point in and 1 out or vice versa, as shown in the inset of Fig. 3 111The mechanism by which a 3-1 ground state appears is quite different from previous work, where a partially ordered monopole crystal has been proposed Brooks-Bartlett et al. (2014); Jaubert (2015); Jaubert and Moessner (2015) and multiferroic behaviour discussed Khomskii (2012).. The combined manifold of spin and lattice configurations is extensive sup , and so a spin-lattice liquid may be realised at low $T$.

The full phase diagram can be explored using Monte Carlo simulation, and is shown in Fig. 3 (see supplemental material sup for details of the simulations). We set $g/(2K)=0.2$, since this corresponds to the experimental findings at low $T$ ($\delta\alpha_{\sf gs}=0.2$ $\approx 11^{\circ}$) Thygesen et al. (2017). As expected, for $g^{2}/(4K)<J_{\sf is}(\alpha_{\sf av})$ there is no lattice distortion at any temperature, and, below a specific heat peak at $T\approx 0.8J_{\sf is}$, the system shows the characteristic behaviour of spin ice 222Since analytical calculations show $J_{\sf is}(\alpha_{\sf av})$ is close to 0, we have checked that the phase diagram for $J_{\sf is}(\alpha_{\sf av})<0$ is qualitatively similar to that for $J_{\sf is}(\alpha_{\sf av})>0$. The main difference is that the value of $g^{2}/(4K)$ at which all-in-all-out order turns into a spin-lattice liquid is larger than the value at which a spin-ice becomes a spin-lattice liquid..

For $g^{2}/(4K)>J_{\sf is}(\alpha_{\sf av})$ there are two transitions as temperature is reduced, and this can be seen in simulations of the heat capacity shown in Fig. 4a. At $T=T_{1}$ there is a first-order liquid-gas-like transition in which $\delta\alpha(T)$ jumps discontinuously from $\delta\alpha(T>T_{1})\approx 0$ to $\delta\alpha(T=T_{1}^{-})\approx J_{\sf is}(\alpha_{\sf av})/g$. Below this there are three classes of bonds, long and short bonds with coupling $J_{\sf is}(\alpha_{\sf av})\pm g\delta\alpha(T)$ and medium bonds with coupling $J_{\sf is}(\alpha_{\sf av})$ (see inset to Fig. 3). For a fixed lattice configuration, the long and short bonds form non-intersecting loops that visit every lattice site and on which the sign of the coupling alternates. Strong spin correlations form on these loops at the transition temperature, favouring alternating pairs of $\sigma=\pm 1$ sup . The correlations increase on further decreasing $T$, and this is related to the increasing density of tetrahedra with a 3-1 spin configuration (see Fig. 4b). The loops themselves are not static, but constantly rearranging themselves subject to the 2-in-2-out constraint on the lattice displacements. The distribution of loop lengths follows the same power-law scaling as is found in spin-ice Jaubert et al. (2011), despite the interaction of lattice displacements via spin degrees of freedom (see supplemental material sup for more details about loops).

At lower $T$ there is a transition at $T=T_{2}$ that is likely second order, and in which the system excludes loops with length $l_{\sf loop}=4m+2$ ($m=1,2,...$) in favour of those with length $l_{\sf loop}=4m+4$, which we call loop symmetry breaking (see Fig. 4c and sup ). This occurs because loops with length $l_{\sf loop}=4m+2$ cannot simultaneously minimise the energy of every bond, as they are constrained to have an odd number of unsatisfied bonds, while loops with length $l_{\sf loop}=4m+4$ can satisfy all the bonds (or have an even number of unsatisfied bonds) sup . Below $T=T_{2}$ essentially all tetrahedra adopt a 3-1 spin configuration (Fig. 4b) and the energy saturates, but the number of allowed configurations remains extensive.

The essence of the resultant spin-lattice liquid state is that it is built by nonperturbatively coupling two classical spin liquids, which in the present case are both spin ice. As with classical spin liquids, spin-lattice liquids can be described by a coarse-grained gauge theory, and, while this inherits some of the structure of the underlying spin-liquid gauge theories Henley (2010); Castelnovo et al. (2012), the nonperturbative coupling provides additional constraints on the allowed configurations. A detailed exploration of such an effective field theory is left for future work.

Relation to experiment: In order to determine an approximate value of the control parameter $g^{2}/(4KJ_{\sf is}(\alpha_{\sf av}))$ for Y2Mo2O7, the calculated value of $g/J_{\sf is}(\alpha_{\sf av})\sim 20$ can be combined with the experimentally determined $\delta\alpha_{\sf gs}=g/(2K)\approx 0.2$ Thygesen et al. (2017), to give $g^{2}/(4KJ_{\sf is}(\alpha_{\sf av}))\sim 2$. This places the material on the right of the phase diagram in Fig. 3, firmly within the spin-lattice liquid region.

The predictions of the model are predominantly aimed at understanding the intermediate temperature regime of Y2Mo2O7, above the spin-glass freezing temperature, $T_{\sf f}$, but below the energy of the first crystal-field excitation, $\Delta E$. This should be compared with the model in the region $T_{2}<T<T_{1}$. Here, the discreteness of the Ising degrees of freedom, coupled with the discrete distortions of the Mo ions provide an explanation for the multiple peaks observed in NMR spectra of Y nuclei Keren and Gardner (2001). It is also consistent with the finding in $\mu$sr and NMR $1/T_{2}$ measurements that the lattice distortions are temperature dependent and respond strongly to magnetic field Sagi et al. (2005); Ofer et al. (2010).

At lower temperatures the degeneracy of the model means that other subdominant effects are likely to play a role. To understand the experimentally observed spin-glass region ($T\leq T_{\sf f}$) it is likely necessary to take into account some combination of disorder, further-neighbour interactions and additional exchange couplings arising from higher-energy crystal-field levels. In particular, small local changes in the Mo-O-Mo bond angle, (i.e. in $\delta\alpha$) will result in additional disorder in the exchange couplings. One promising feature of the model is that its highly-correlated low-temperature dynamics likely make it very susceptible to spin-glass freezing, particularly in the vicinity of $T_{2}$, where Monte-Carlo simulations become difficult to equilibrate 333Using Monte Carlo to simulate the effect of disorder, and test whether it drives a spin-glass transition at temperatures close to $T_{2}$, is challenging, since it is very difficult to equilibrate the simulations due to the non-local nature of the dynamics and the strong correlation between spin and lattice degrees of freedom. We leave such a detailed investigation to future work..

At high temperatures ($T\sim\Delta E\sim 300$K) the excited crystal-field levels become thermally populated, and therefore the approximation of being in the low-energy doublet breaks down. Since $\Delta E$ is smaller than the temperature at which NMR linewidth extrapolations determine that the lattice distortion dies away ($\sim 430$K) Keren and Gardner (2001), the first-order transition predicted by the model is likely smeared out due to the additional local degrees of freedom and the opening of new hopping pathways.

A particular mystery in Y2Mo2O7 is why susceptibility measurements show a negative (AFM) Curie-Weiss temperature of $\theta_{\sf cw}\approx-200$K in the high temperature range $500-800$K Gardner et al. (2004), a different Curie-Weiss temperature of $\theta_{\sf cw}\approx-41$K in the intermediate range $50-300$K Silverstein et al. (2014) but ring features surrounding the $\Gamma$ point in low-temperature neutron-scattering measurements that can only be explained by including FM coupling Silverstein et al. (2014). Our model provides a resolution by proposing a mechanism by which AFM and FM interactions can coexist. Simulations in the temperature range $T_{2}<T<T_{1}$ show an emergent Curie-Weiss behaviour with $\theta_{\sf cw}<0$ (see Fig. 4d), despite the presence of strong correlations. Above $T\sim 300K$ the local electronic state is no longer confined to the ground state doublet, and a change in the Curie-Weiss behaviour is unsurprising.

One way to test for the existence of a spin-lattice liquid state in Y2Mo2O7 would be to use magnetic field to control the lattice or uniaxial pressure to control the magnetism. Large magnetic fields will suppress the 2-in-2-out lattice displacements by selecting spin configurations (either 2-in-2-out or 3-1) that are incompatible with the displacement. Intermediate fields will have a more subtle effect, splitting the low-energy spin-lattice configurations, and could be an interesting direction for further study.

In conclusion, we have proposed the idea of a spin-lattice liquid state, arising from strong coupling between two “spin” liquids, one associated with the spin-orbital and the other with the lattice degrees of freedom. Furthermore, we have shown that it may be realised in Y2Mo2O7 and related compounds, driven by a giant magnetoelastic coupling effect. According to our theory, the strength of the magnetoelastic coupling is the main thing that distinguishes the molybdate pyrochlores from $f$-electron spin-ice compounds such as Dy2Ti2O7 Bramwell and Gingras (2001). Finally, we end with the hope that other materials may be found, in which the dance of spin and lattice extends down to even lower temperature.

Acknowledgments. We thank Reinhard Kremer, Ludovic Jaubert, Joe Paddison and Hidenori Takagi for useful comments on the manuscript.

I Local electronic state of $d^{2}$ ions

Here we consider in more detail the local electronic state of an ion with two electrons in the $d$-shell ($d^{2}$ configuration), residing in a trigonally distorted oxygen octahedra and subject to spin-orbit coupling. The largest energy scale is the cubic crystal-field splitting of the $4d$ electron states, and as a result we consider only the $t_{\sf 2g}$ levels in which the $d$ electron orbitals point in between the surrounding oxygen sites Fazekas (1999). Further splitting occurs due to Hund’s coupling, which favours parallel spin alignment, resulting in a low-energy manifold of states labelled by spin $S=1$ and effective angular momentum $L_{\sf eff}=1$. The 9-fold degeneracy of these states is split by the combination of a trigonal flattening of the surrounding oxygen octahedra and spin orbit coupling, and we consider the local Hamiltonian,

[TABLE]

where ${\bf S}$ and ${\bf L}_{\sf eff}$ are measured relative to the axis of the trigonal distortion, which points into/out-of Mo tetrahedra and defines the local $z^{\prime}$ axis (see Fig. 1 in the main text).

Diagonalisation of $\mathcal{H}_{\sf loc}$ gives the energy levels shown in Fig. S1. In general, the presence of a trigonal distortion means that the levels cannot be labelled by their total angular momentum, ${\bf J}_{\sf eff}={\bf S}+{\bf L}_{\sf eff}$. However, the lowest energy doublet is independently an eigenstate of both parts of $\mathcal{H}_{\sf loc}$ and can therefore be written in terms of the total angular momentum as $|J_{\sf eff}=2,J_{\sf eff}^{z^{\prime}}=\pm 2\rangle=|S^{z^{\prime}}=\pm 1,L_{\sf eff}^{z^{\prime}}=\pm 1\rangle$.

In the main text we consider this lowest-energy doublet to be the only accessible state at low temperature, and this is valid as long as $T\ll\Delta E=\lambda_{\sf SO}[1-\lambda_{\sf SO}/\Delta_{\sf trig}+\dots]$, where $k_{B}=1$ and the energy difference between the levels has been expanded in terms of $\lambda_{\sf SO}/\Delta_{\sf trig}$. In Y2Mo2O7 band structure calculations have estimated that $\Delta_{\sf trig}\approx 160$meV and $\lambda_{\sf SO}\approx 40$meV Shinaoka et al. (2013), in reasonable agreement with the $\lambda_{\sf SO}\approx 58$meV found experimentally for isolated Mo4+ ions Abragam and Bleaney (1970). In consequence one finds $\Delta E\approx 30$meV$\approx 300$K.

II Superexchange Hamiltonian on the pyrochlore lattice

Here we derive the superexchange Hamiltonian used in the main text. First we consider the idealised case of a pyrochlore lattice surrounded by regular oxygen octahedra. While this is never exactly the case in materials such as Y2Mo2O7, due to a tension between the oxygen octahedra surrounding the Mo ions and the oxygen cubes surrounding the Y ions, it serves as a useful starting point. We then go on to consider the effect of trigonally disorting the oxygen octahedra, and show that the distortion can have a significant effect on the superexchange Hamiltonian.

II.1 Superexchange in idealised case with regular octahedra

In the Y2Mo2O7 structure, oxygen octahedra surrounding each of the 4 basis sites of the Mo pyrochlore lattice have different orientations. As such it is useful to set-up a local cubic coordinate system for each site, as shown in Fig. S2. Rotation matrices from local cubic to global cubic (i.e. crytallographic) coordinates are given by,

[TABLE]

We consider oxygen mediated $d$-$p$-$d$ hopping between neighbouring sites. The $d$-$p$ hopping matrices taking $(d_{yz},d_{zx},d_{xy})$ to $(p_{x},p_{y},p_{z})$ are given by Slater and Koster (1954); Harrison (1989) (see also Ref. Pesin and Balents (2010)),

[TABLE]

where $V_{pd\pi}$ parametrises the $\pi$-bond hopping, the $\pm$ refers to a positive or negative oxygen coordinate and ${\sf x,y,z}$ labels the axis on which the oxygen lies. On “integrating out” the oxygen orbitals one finds the $d$-$d$ hopping matrices,

[TABLE]

where $t=V_{pd\pi}^{2}/\Delta_{pd}$, $\Delta_{pd}$ gives the energy difference between the $d$ and $p$ orbitals and $P_{ij}$ is chosen appropriately for the bond in question. For the regular octahedra shown in Fig. S2, this gives the bond-dependent hopping matrices,

[TABLE]

and it can be seen that the dominant hopping process is off-diagonal in the local coordinates. The nearest-neighbour hopping Hamiltonian is given by,

[TABLE]

with $\alpha,\beta\in\{d_{\sf yz},d_{\sf zx},d_{\sf xy}\}$ measured in the local cubic coordinates and $s\in{\uparrow,\downarrow}$.

It is also necessary to specify the on-site Coulomb interaction, and in the local octahedral coordinates this is given by Castellani et al. (1978),

[TABLE]

where $U$ is the on-site Coulomb repulsion between electrons in the same orbital, $J_{\sf H}$ is Hund’s coupling, $n_{\alpha s}=\alpha^{\dagger}_{s}\alpha^{\phantom{\dagger}}_{s}$, $n_{\alpha}=n_{\alpha\uparrow}+n_{\alpha\downarrow}$ and $\textbf{S}_{\alpha}=(\frac{1}{2}[\alpha^{\dagger}_{\uparrow}\alpha^{\phantom{\dagger}}_{\downarrow}+\alpha^{\dagger}_{\downarrow}\alpha^{\phantom{\dagger}}_{\uparrow}],\frac{1}{2i}[\alpha^{\dagger}_{\uparrow}\alpha^{\phantom{\dagger}}_{\downarrow}+\alpha^{\dagger}_{\downarrow}\alpha^{\phantom{\dagger}}_{\uparrow}]),\frac{1}{2}[\alpha^{\dagger}_{\uparrow}\alpha^{\phantom{\dagger}}_{\uparrow}-\alpha^{\dagger}_{\downarrow}\alpha^{\phantom{\dagger}}_{\downarrow}])$.

A spin-orbital Hamiltonian can be derived by the standard method of treating $t^{2}/U$ as a small parameter and considering virtual hopping processes of the type $d^{2}d^{2}\rightarrow d^{3}d^{1}\rightarrow d^{2}d^{2}$ Khaliullin et al. (2001); Horsch et al. (2003); Khaliullin (2005). This results in,

[TABLE]

where $\gamma$ refers to the three types of bonds (01 & 23 ; 02 & 13 ; 03 & 12) shown in Fig. S2, $\textbf{S}_{i}$ is a spin-1 operator on site $i$ measured in the global coordinate system,

[TABLE]

$\eta=J_{H}/U$, $R=1/(1-3\eta)$, $r=1/(1+2\eta)$ and,

[TABLE]

where the orbital operators are bond dependent pseudospins and given by $\tau^{z}=\frac{1}{2}[d_{\sf yz}^{\dagger}d_{\sf yz}-d_{\sf zx}^{\dagger}d_{\sf zx}]$, $\tau^{+}=d_{\sf yz}^{\dagger}d_{\sf zx}$, $\tau^{-}=d_{\sf zx}^{\dagger}d_{\sf yz}$ and $n=d_{\sf yz}^{\dagger}d_{\sf yz}+d_{\sf zx}^{\dagger}d_{\sf zx}$ on 03 and 12 bonds and by related expressions on other bonds.

II.2 The effect of trigonal distortion and spin-orbit coupling

The effect of the trigonal distortion on the spin-orbital Hamiltonian is twofold, since it is necessary to take into account a modified hopping matrix, as well the way in which it combines with spin-orbital coupling to split the otherwise 9-fold degenerate local energy levels.

When considering trigonally distorted oxygen octahedra, we continue to use a local cubic coordinate system that is aligned with the idealised undistorted octahedra (see Fig. 1 in the main text). As a consequence, oxygen ions no longer sit on the axes of the local coordinates, and thus the $d$-$p$ hopping matrices become more complicated. For example, using the tables given in Ref. Slater and Koster (1954); Harrison (1989), the $P_{\sf z}^{+}$ matrix given in Eq. S25 should be replaced by,

[TABLE]

where $\textbf{r}_{\textsf{O}}=(l,m,n)$ is the position of the oxygen ion in the local cubic coordinate system. Other hopping matrices follow in a similar manner. For Mo4+ ions it has been determied that $V_{pd\sigma}/V_{pd\pi}\approx-2.2$, and in fact this ratio is only weakly ion dependent Harrison (1989).

For trigonal distortions of the type found to exist in Y2Mo2O7, the position of the oxygen ion on a 03 bond can be written as,

[TABLE]

where $\alpha_{\sf cub}=2\arctan 2\sqrt{2}\approx 141^{\circ}$ is the bond angle in the case of regular octahedra and $\alpha$ is the bond angle of the trigonally distorted structure, which for the average $Fd\bar{3}m$ structure of Y2Mo2O7 is given by $\alpha\approx 127^{\circ}$. Inserting these values into Eq. S26 for the $d$-$d$ hopping matrix, one finds for the 03 bond and for $\alpha=127^{\circ}$,

[TABLE]

Comparison to Eq. S36 for the case of regular oxygen octahedra shows that the balance between diagonal and off-diagonal hopping of $d_{\sf yz}$ and $d_{\sf zx}$ orbitals has changed significantly, while hopping to/from $d_{\sf xy}$ orbitals remains unimportant. For regular octahedra the off-diagonal $d_{\sf yz}\leftrightarrow d_{\sf zx}$ hopping dominates, while for trigonally distorted octahedra the diagonal $d_{\sf yz}\leftrightarrow d_{\sf yz}$ and $d_{\sf zx}\leftrightarrow d_{\sf zx}$ hopping becomes more important. The reason for this is the appearance of a $\sigma$-bond hopping channel that competes with the original $\pi$-bond hopping channel as soon as the oxygen octahedra are not regular (see Fig. 1 in the main text). Due to the relative importance of the $\sigma$ channel ($V_{pd\sigma}^{2}\approx 4.8V_{pd\pi}^{2}$) even a relatively small change in the bond angle away from the case of regular octahedra can significantly alter the hopping matrix, and this is at the heart of the giant magnetoelastic effect proposed in the main text.

The combination of the trigonal distortion and the spin-orbit coupling split the 9-fold degenerate Hund’s-rule coupled $t_{2g}$ levels ($S=1$, $L_{\sf eff}=1$), resulting in a ground state doublet labelled by $J_{\sf eff}^{z^{\prime}}=\pm 2$, where the $z^{\prime}$ axis points into or out of the tetrahedra and defines a trigonal coordinate system (see Fig. S2). If the energy gap between the doublet and the higher energy levels is considerably larger than the exchange coupling strengh, $\sim t^{2}/U$, then the higher energy levels can be ignored. We assume this is the case and therefore project a modified version of $\mathcal{H}_{\textsf{SO}}$ [Eq. S39] that takes the bond distortion into account into the manifold of local $J_{\sf eff}^{z^{\prime}}=\pm 2$ states. Since none of the operators appearing in $\mathcal{H}_{\textsf{SO}}$ [Eq. S39] connect $J_{\sf eff}^{z^{\prime}}=2$ and $J_{\sf eff}^{z^{\prime}}=-2$, the Hamiltonian takes the simple Ising form,

[TABLE]

where $\sigma=\pm 1$ represents moments pointing into or out of tetrahedra along the local trigonal $z^{\prime}$ axis.

The value of $J_{\sf is}(\alpha,\eta)$ can be determined by comparing the different configurational energies of a pair of neighbouring moments, resulting in,

[TABLE]

where $t_{a_{1g}}(\alpha)$, $t_{d}(\alpha)$, $t_{od}(\alpha)$ and $t_{a_{1g}e}(\alpha)$ are hopping parameters in the trigonal coordinate system, defined for 03 bonds according to,

[TABLE]

where the trigonal basis is $(a_{1g},e_{g+}^{\prime},e_{g-}^{\prime})$. For other bonds the U must be modified, but this does not affect the value of $J_{\sf is}(\alpha,\eta)$. Fig. 2 in the main text is plotted using the expression for $J_{\sf is}(\alpha,\eta)$ given in Eq. S52.

III Physical properties of $\mathcal{H}_{\sf sl}$

Here we describe in more detail the properties of the spin-lattice Hamiltonian $\mathcal{H}_{\sf sl}$ [Eq. 1 in the main text].

III.1 Monte Carlo simulation details

We used Monte Carlo simulations to determine the properties of $\mathcal{H}_{\sf sl}$. Monte Carlo updates were performed consecutively for spin and lattice degrees of freedom. For the spin degrees of freedom a single spin-flip algorithm was employed. Updates of the lattice degrees of freedom (displacements into or out of a tetrahedra) were split into two steps. First a worm algorithm was used to move between lattice configurations in which each tetrahedra has a 2-in-2-out pattern Jaubert (2009). The worm was formed taking into account the interaction of the lattice with the (fixed) spin degrees of freedom, and could therefore be made rejection free. Second, a global change to the magnitude of the lattice distortion, $\delta\alpha$, was proposed, and accepted or rejected according to the Metropolis criteria. Parallel tempering was employed to further improve equilibration.

Typically a Monte Carlo step involved $N$ attempts to make a single spin flip, one worm step and 100 attempts to change $\delta\alpha$, with parallel tempering after each 10 Monte Carlo steps. The number of spins was $N=16L^{3}$, where each unit cell has 16 sites, $L$ is the linear size of a cubic cluster and periodic boundary conditions were used. Simulations were performed for $L=3,4,\dots,12$ (i.e. $N=432,1024,\dots,27648$). Equilibriation at the lower transition ($T=T_{2}$) proved difficult except for the smallest cluster sizes ($L=3,4$), due to the strong correlation between spin and lattice degrees of freedom. Overcoming this equilibration difficulty would almost certainly require a combined spin and lattice update. Typical acceptance ratios of the single spin-flip algorithm are shown in Fig. S3, and it can be seen that they become very small in the neighbourhood of $T=T_{2}$.

III.2 Thermodynamic properties

The Monte Carlo simulations were used to determine the thermodynamic properties of $\mathcal{H}_{\sf sl}$ and thus map out the phase diagram shown in Fig. 3 of the main text.

The first-order nature of the transition at $T=T_{1}$ can be ascertained from energy histogram analysis, as shown in Fig. S4. The double-peak structure is a sign of phase coexistence at the transition temperature. That the transition is of the liquid-gas type can be seen from the behaviour of the lattice displacement parameter, $\delta\alpha$, which plays the role of density, and jumps discontinuously at the transition (see Fig. S4). It can also be seen in Fig. S4 that the transition temperature is relatively insensitive to system size.

The nature of the lower transition is more difficult to determine, since the Monte Carlo simulations encounter equilibration problems in the vicinity of the transition for larger system sizes. Energy histogram analysis of small system sizes ($L=3,4,6$) appear to show a second-order transition. However, finite size scaling analysis suffers from uncertainty due to the equilibration difficulties.

III.3 Loops

Here we explain in more detail the microscopic structure in the spin-lattice liquid phase (see Fig. 3 in the main text for the phase diagram).

The building blocks of the pyrochlore lattice are tetrahedra, and in the spin-lattice liquid phase these have two ions displaced towards and two away from their centres. Each tetrahedron thus has one long bond, four medium-length bonds and one short bond, and the long and short bonds are opposite one another (see Fig. S5). The effective spin model on a tetrahedron with a fixed lattice displacement consists, in the local Ising basis, of one ferromagnetic bond (long), four weakly antiferromagnetic bonds (medium) and one more strongly antiferromagnetic bond (short). As a result, an isolated tetrahedron with fixed lattice displacement has a fourfold-degenerate spin ground state, with three spins pointing in and one out or vice versa (see Fig. 3 in the main text). In these configurations, the four medium-length bonds host two satisfied interactions and two unsatisfied interactions, and therefore do not contribute to the energy.

On the pyrochlore lattice, the requirement that lattice displacements follow a 2-in-2-out rule results in the formation of loops consisting of alternating long and short bonds. Each lattice site is visited by exactly one loop, and the loops are non-intersecting and are closed in the presence of periodic boundary conditions. The shortest possible loop has length $l_{\sf loop}=6$, and, due to the structure of the pyrochlore lattice, loops have to be even in length. Some example loops are shown in Fig. S6.

Configurations of the system are labelled by a combination of the loop structure (or equivalently the in/out lattice displacements) and the spin directions. There is an extensive set of degenerate ground-state configurations in which all tetrahedra have a 3-1 (three in and one out or vice versa) spin configuration of the type shown in Fig. S5 (see also Fig. 3 in the main text). These ground-state (gs) configurations consist of loops of length $l_{\sf loop}=4+4m$ with $m=1,2,3,\dots$ (i.e. loops of length $l_{\sf loop}=2+4m$ are excluded) and Ising spin configurations on the loops of $\boldsymbol{\sigma}^{\sf loop,gs}=(1,1,-1,-1,1,1,-1,-1,\dots)$. There are two types of local move that leave the system within the ground state manifold. One possibility is that the spins on a loop are all flipped simultaneously, taking $\boldsymbol{\sigma}^{\sf loop,gs}\to-\boldsymbol{\sigma}^{\sf loop,gs}$ (see Fig. S7). A second possibility is that the bond displacements are exchanged around a closed set of medium-length bonds thus rearranging the loop structure of the lattice (see Fig. S8 for an example of a non-ground-state loop update). Typically it is simultaneously necessary to flip some spins to ensure that the reconstructed loops retain the correct $\boldsymbol{\sigma}^{\sf loop,gs}$ structure.

One example of an excitation out of the ground state is the introduction of kinks into the otherwise periodic Ising configuration on a loop (see Fig. S7). For loops of length $l_{\sf loop}=4+4m$ the number of kinks has to be even, and each kink is associated with a 2-in-2-out (or potentially all-in-all-out) spin configuration on the relevant tetrahedron. Excitations can also involve the creation of loops with length $l_{\sf loop}=2+4m$, and the reason that these cost energy is that they have to contain at least one kink (see Fig. S6a). Further kinks can be created in pairs, and so the number of kinks is constrained to be odd.

Monte Carlo simulations show two phase transitions. At the high-temperature first-order phase transition ($T=T_{1}$) the lattice displacement parameter, $\delta\alpha$, acquires a finite value. For $T<T_{1}$ the bonds can therefore be divided into long, medium and short lengths with differing spin interactions, and it becomes meaningful to discuss the system in terms of loops of alternating long and short bonds. A non-local loop-spin order parameter can be constructed as the product of Ising spins around the loop projected onto its ground state value and given by,

[TABLE]

where $\mathcal{O}_{\sf loop}=0$ for loops with uncorrelated spins and $\mathcal{O}_{\sf loop}=1$ in the ground state. The behaviour of this order parameter in Monte Carlo simulations is shown in Fig. S9, and it can be seen that it jumps from $\mathcal{O}_{\sf loop}\approx 0$ for $T>T_{1}$ to a finite value for $T<T_{1}$. Further reducing the temperature builds up spin correlations on the loops and therefore increases the value of $\mathcal{O}_{\sf loop}$.

The distribution of loop lengths can be measured by Monte Carlo simulations. In the bulk of the intermediate phase ($T_{2}<T<T_{1}$) we find a power-law distribution of loop lengths, with no sign of a symmetry breaking between loops of length $l_{\sf loop}=4+4m$ and $l_{\sf loop}=2+4m$. This is shown in Fig. S10, where it can be seen that there is a crossover between different power laws at $l_{\sf loop}=L^{2}$. The crossover is the effect of performing simulations on finite size clusters, where loops can be divided into those that wind or don’t wind the system Jaubert et al. (2011). The power-law distribution matches well to that already determined for a 2-in-2-out manifold of states without further interactions Jaubert et al. (2011). As such, the fact that lattice displacements are effectively coupled via the spin degrees of freedom, doesn’t restrict their exploration of the full 2-in-2-out manifold in the bulk of the spin-lattice-liquid phase.

At lower temperature, Monte Carlo simulations show that there is a phase transition at $T=T_{2}$ at which the system begins to favour loops of length $l_{\sf loop}=4+4m$ over those with $l_{\sf loop}=2+4m$. This can be seen from studying the distribution of loop lengths, and a simple way of doing this is by considering,

[TABLE]

where $N_{i+4m}$ is the number of sites belonging to loops of lengths $l_{\sf loop}=i+4m$. $\mathcal{O}_{\sf mis}$ measures the mismatch in the number of sites belonging to the two loop-length classes, and simulation results for its behaviour can be seen in Fig. S11. It takes a relatively low value for $T>T_{2}$ before rapidly increasing at $T_{2}$ and saturating at $\mathcal{O}_{\sf mis}=1$ approaching $T=0$. It can be seen in Fig. S11 that increasing the system size makes the transition sharper, but the reliability of finite size scaling analysis is hampered by equilibration difficulties.

III.4 Structure factor

For comparison with neutron scattering studies it is useful to calculate the structure factor. Since neutrons scatter both off nuclei and from magnetic moments, we calculate both the spin and lattice structure factors.

The spin structure factor is defined by,

[TABLE]

where $S^{\alpha}({\bf r})$ denotes the $\alpha$th component of spin, measured in the crystallographic coordinates and at position ${\bf r}$, and $S^{\alpha}({\bf q})$ is its Fourier transform. It should be remembered that in the spin-lattice liquid region, the lattice positions of the spins are not fixed but dynamical, and it is necessary to take this into account. Fig. S12a shows a cut through the spin structure factor in the spin-lattice liquid state, and displays structure arising from spin correlations. The cut taken is the $hhl$ plane, where $hkl$ are used to label $(q_{\sf a}/2\pi,q_{\sf b}/2\pi,q_{\sf c}/2\pi)$.

Due to the disordered lattice displacements in the spin-lattice liquid state, nuclear scattering also contributes to the diffuse scattering Paddison (2019). The total nuclear scattering from Mo ions and their surrounding O octahedra is given by,

[TABLE]

where $b_{\sf Mo}=6.715$fm and $b_{\sf O}=5.803$fm are the neutron scattering lengths of Mo and O and ${\bf r}_{\sf Mo}$ and ${\bf r}_{\sf O}$ denote the positions of the Mo and O nuclei. An example structure factor in the spin-lattice liquid state is shown in Fig. S12b and includes both Bragg peaks arising from the average nuclear positions, as well as diffuse scattering arising from the correlated but disordered displacements.

The total neutron scattering signal will have contributions from both the spin and lattice structure factors, with the relative weight dependent on material parameters.