Spin-orbital glass transition in a model of a frustrated pyrochlore magnet without quenched disorder

TL;DR

This paper theoretically demonstrates that in a quenched-disorder-free pyrochlore magnet, spin and orbital degrees of freedom can undergo a simultaneous glass transition driven by their mutual interplay and Jahn-Teller distortions.

Contribution

It introduces a model showing disorder-free spin glass behavior arising from spin-orbital interactions mediated by lattice distortions.

Findings

Power-law divergence of relaxation times observed

Negative divergence of magnetic and dielectric susceptibilities

Simultaneous spin and orbital glass transition detected

Abstract

We show theoretically that spin and orbital degrees of freedom in the pyrochlore oxide Y2Mo2O7, which is free of quenched disorder, can exhibit a simultaneous glass transition, working as dynamical randomness to each other. The interplay of spins and orbitals is mediated by the Jahn-Teller lattice distortion that selects the choice of orbitals, which then generates variant spin exchange interactions ranging from ferromagnetic to antiferromagnetic ones. Our Monte Carlo simulations detect the power-law divergence of the relaxation times and the negative divergence of both the magnetic and dielectric non-linear susceptibilities, resolving the long-standing puzzle on the origin of the disorder-free spin glass.

Figures

Equations

$$H=\sum _{<ij>}{J}_{{\mathbf{\sigma }}_i,{\mathbf{\sigma }}_j}{\mathbf{S}}_i\cdot {\mathbf{S}}_j-ϵ\sum _{<ij>}{\mathbf{\sigma }}_i\cdot {\mathbf{\sigma }}_j(ϵ>0)$$

$$J_{{\mathbf{\sigma }}_i,{\mathbf{\sigma }}_j}=J[1+\delta ({\hat{\mathbf{r}}}_{ij}\cdot {\mathbf{\sigma }}_i+(-{\hat{\mathbf{r}}}_{ij})\cdot {\mathbf{\sigma }}_j)](\delta >0),$$

$${\chi }_3=\text{evaluated}{\frac{1}{3}\sum _{\mu }\frac{{\partial }^3⟨m_{\mu }{⟩}_{\mathrm{eq}}}{\partial h_{\mu }^3}}_{h_{\mu }\to 0},$$

$${\chi }_3^{\sigma }=\text{evaluated}{\frac{1}{4}\sum _{\nu }\frac{{\partial }^3⟨p_{\nu }{⟩}_{\mathrm{eq}}}{\partial E_{\nu }^3}}_{E_{\nu }\to 0},$$

$$\mathcal{H}=$$

$$+\sum _{i=\mathrm{A},\mathrm{B}}\sum _{a,a^{\prime }}\sum _{\sigma ,{\sigma }^{\prime }}\frac{K(a,a^{\prime })-{\delta }_{\sigma ,{\sigma }^{\prime }}J(a,a^{\prime })}{2}{c}_{ia\sigma }^{†}{c}_{ia\sigma }{c}_{i{a}^{\prime }{\sigma }^{\prime }}^{†}{c}_{i{a}^{\prime }{\sigma }^{\prime }}+\frac{U_o}{2}\sum _b{c}_{ob\sigma }^{†}{c}_{ob\sigma }{c}_{o{b}^{\prime }{\sigma }^{\prime }}^{†}{c}_{o{b}^{\prime }{\sigma }^{\prime }},$$

$${\chi }_3$$

$$\displaystyle=\frac{\beta^{3}N^{3}}{3}\sum_{\mu}\Bigl{(}\expectationvalue{m_{\mu}^{4}}_{\rm eq}-4\expectationvalue{m_{\mu}}_{\rm eq}\expectationvalue{m_{\mu}^{3}}_{\rm eq}-3\expectationvalue{m_{\mu}^{2}}_{\rm eq}^{2}+12\expectationvalue{m_{\mu}^{2}}_{\rm eq}\expectationvalue{m_{\mu}}_{\rm eq}^{2}-6\expectationvalue{m_{\mu}}_{\rm eq}^{4}\Bigr{)}.$$

$$\displaystyle\chi_{3}=\frac{\beta^{3}N^{3}}{3}\sum_{\mu}\Bigl{(}\expectationvalue{m_{\mu}^{4}}_{\rm eq}-3\expectationvalue{m_{\mu}^{2}}_{\rm eq}^{2}\Bigr{)}.$$

$$\displaystyle\chi^{(\sigma)}_{3}=\frac{\beta^{3}N^{3}}{4}\sum_{\nu}\Bigl{(}\expectationvalue{p_{\nu}^{4}}_{\rm eq}-3\expectationvalue{p_{\nu}^{2}}_{\rm eq}^{2}\Bigr{)}.$$

Full text

Glass formation is a generic phenomenon that emerges in strongly frustrated systems tarjus2005frustration, ranging from super-cooled molecular liquids and densely packed soft-matters angell2000relaxation; sciortino2005glassy; cavagna2009supercooled; berthier2011theoretical to geometrically frustrated spin systems on kagome martinez1992magnetic; ladieu2004relative; hamp2018supercooling; cugliandolo2019mean; upadhyay2017magnetic and pyrhoclore lattices alba1982very; greedan1986spin; reimers1991short; gingras1996nonlinear; gingras1997static; gardner1999glassy; zhou2008unconventional; taniguchi2009spin; ladieu2004relative; thygesen2017orbital; booth2000local. Frustration suppresses ordinary long-range orders and stabilizes liquid/paramagnetic states down to unusually low temperatures where the effect of many-body interactions become dominant. In sharp contrast to normal liquids at higher temperatures, molecules in super-cooled liquids can flow only through collective motions because of the strong inter-molecular interactions. In classical antiferromagnets on triangular, kagomé and pyrochlore lattices, the frustration among interactions caused by the lattice geometry precludes the possibility of all bonds to be simultaneously satisfied in energy diep2013frustrated. An outstanding open question there is whether truly thermodynamic glass transitions, for example the putative Kauzmann transition of super-cooled liquids Ka48, are possible. At present such ideal glass transitions is attained only theoretically in unphysical limits of infinitely large dimensions for some super-cooled liquids parisi2010mean; kurchan2012exact; charbonneau2014fractal or spin systems with infinitely large number of spin components yoshino2018disorder.

In realistic three-dimension, a true thermodynamic glass transition is established only in lattice systems with quenched disorder added on top of the frustration mydosh2014spin; binder1986spin; kawamura2015spin. Theoretically, this spin glass(SG) is found in the Edwards-Anderson model edwards1975theory; mezard1987spin based on both Ising ogielski1985dynamics; bhatt1988numerical; kawashima1996phase and Heisenberg spins hukushima2000chiral; lee2007large; viet2009numerical; ogawa2019monte.

The next fundamental question is whether the SG transition can be sustained even when such quenched disorder is switched off. It is known experimentally that some of the disorder-free pyrochlore magnets exhibit a sharp SG transition alba1982very; greedan1986spin; reimers1991short; gingras1996nonlinear; gingras1997static; gardner1999glassy; zhou2008unconventional; taniguchi2009spin; ladieu2004relative, which is indistinguishable from the conventional one with quenched disorder mydosh2014spin; binder1986spin; kawamura2015spin. The divergence of nonlinear susceptibility gingras1996nonlinear establishes this phase as SG. In contrast,, any of the theories available at present does not allow for the 3D SG without quenched disorder.

This Letter aims to clarify the origin of the SG transition in the disorder-free pyrochlore magnets by constructing a realistic model consisting of two key degrees of freedom of the system, i. e. spins and orbitals which are coupled through a Jahn-Teller (JT) lattice distortion which modifies the super-exchange interactions between the spins dynamically. The two degrees of freedom, if decoupled, live in highly degenerate energy landscapes of their own, due to the strong geometrical frustration, so that they remain in their liquid/paramagnetic states at all temperatures. Through the coupling they correlate producing dynamical disorder for each other and create an emergent lugged free-energy landscape. Our extensive MC simulation uncovers an unprecedented thermodynamic glass transition where the two degrees of freedom simultaneously freeze cooperatively.

The prototypical materials for our theory are the pyrochlore oxides, $A_{2}$Mo2O7 ($A=$ Ho, Y, Dy, Tb), which are insulating and show a SG transition at around 20 K gaulin1992spin; dunsiger1996muon; gingras1997static; gardner1999glassy; kezsmarki2004charge; hanasaki2007nature. The magnetic ions Mo${}^{4+}(4d^{2},S=1)$ sit on the vertices of the corner shared tetrahedra (see Fig. 1(a)), and the interactions between them are antiferromagnetic, whose underlying microscopic mechanism is the superexchange mediated by the O2- ions solovyev2003effects. It has been established that the classical spin model with purely antiferromagneic nearest neighbor interactions on the pyrochlore lattice remains nonmagnetic down to zero temperature because of the strong frustration effect reimers1991mean; moessner1998properties, not providing any hint of the SG transitions observed in experiments.

Very recently, some salient features of the local lattice distortions in Y2Mo2O7 thygesen2017orbital were uncovered experimentally. The analysis based on the neutron pair-distribution function data suggests that the Mo4+ ions are locally displaced towards (in) or away from (out) the centers of Mo4 tetrahedra (see Fig. 1(a)), explaining well a large variance of the Mo-Mo distances observed in the x-ray absorption study booth2000local. A natural mechanism of distortion is the JT effect that lowers the electrostatic energy of the Mo ion lifting its orbital degeneracy (Fig.1(c)) (See SI).

The displacements of the four Mo4+ ions on each tetrahedron possibly follow the 2-in-2-out rule, i. e. the ice-rule (see Fig. 1(a,b)), generating huge numbers of ionic configurations over the system equivalent in their JT energies. Such distortions thus do not lead to any long-range structural ordering.

It is empirically known 111The evaluation of the direct exchange interactions between Mo spins can be done based on the Hartree-Fock approach by assuming the orbital ordering, which gives a rough estimate of the sign of the interactions. We give a more detailed analysis by considering the dominant superexchange interactions for the insulating Y2Mo2O7. that the exchange interaction may vary significantly and even changes its sign, depending on the degree of the Mo-O-Mo angle $\alpha$ solovyev2003effects. Indeed, the angle evaluated from the variance of the displacements of the Mo ions distributes in the range $\alpha=116^{\circ}$ and $139^{\circ}$, which is large enough to change the sign of the interaction (Fig.1(d)) (See SI). Then, depending on the configurations of the local lattice distortion in space, the magnetic exchange interaction may acquire a random distribution with finite variances. This can be thought of as follows; If we consider the Jahn-Teller Hamiltonian including only the lattice degrees of freedom, the energy landscape should be such that all the ice-rules give the energy minima of the same height(see Fig.1(d)) pauling1960nature; bramwell2001spin. In putting spins on the ice-type displaced vertices, its energy landscape is modified to an irregular one with random valley structures, since the variance of the spin interactions from ferro to antiferromagnetic ones add different energies to each ice configuration. We verify the above picture showing numerically that the spin and lattice distortions simultaneously undergo a glass transition and freeze into the aperiodic, irregular types of configuration. Since the randomly frozen lattice distortions indicate the randomly selected orbital configurations, we call it a spin-orbital glass.

As a natural description of the material faithful to the experimental observations, we introduce a disorder-free spin model including not only the classical Heisenberg spins $\bm{S}_{i}~{}(i=1,2,3,...,N)$ for the magnetic moments of the Mo${}^{4+}(4d^{2},S=1)$ ions but also the degree of lattice distortion of the Mo ions $\bm{\sigma}_{i}$ represented as vectors taking two discrete values, $\bm{\sigma}_{i}=\sigma_{i}\hat{\bm{e}}_{\nu}$. Here, $\sigma_{i}=\pm 1$ corresponds to either in or out depending on the direction of $\bm{e}_{\nu}$ which is the unit vector in the $[111],[1\bar{1}\bar{1}],[\bar{1}1\bar{1}]$ and $[\bar{1}\bar{1}1]$ directions, respectively for the sub-lattices $\nu=1,2,3,4$ which the $i$-th spin belongs to (See Fig. 1.(a)). The Hamiltonian is given by

[TABLE]

with

[TABLE]

The first term in Eq. (1) represents the exchange interaction whose coupling constant $J_{\bm{\sigma}_{i},\bm{\sigma}_{j}}$ depends on the angle of Mo-O-Mo bond. This form is constructed in a way to reproduce the values of $J_{\rm eff}$ derived microscopically from the perturbation process on the Kanamori-type of Hamiltonian (see SI). As shown in Fig. 1 (d), for realistic values of the on-site Coulomb interactions on Mo4+ and O2- ions solovyev2003effects; shinaoka2013spin, $J_{\rm eff}$ changes its sign within the experimentally observed range of the Mo-O-Mo angle $\alpha=116^{\circ}$-$139^{\circ}$. This is modeled by Eq. (2) as such that $J_{\bm{\sigma}_{i},\bm{\sigma}_{j}}=(1+2\tilde{\delta})J\mbox{ ({in}, {in})},J\mbox{ ({in}, {out})},(1-2\tilde{\delta})J\mbox{ ({out}, {out})},$ where $\tilde{\delta}=\sqrt{6}\delta/3$ (See Fig. 1.(b)). The second term of Eq. (1) represents the elastic energy of the Mo4+ displacement. The elastic energy is minimized if a Mo4+ tetrahedron satisfies the ice-rule. More precisely, the elastic energy of each bond takes $\epsilon\bm{\sigma}_{i}\cdot\bm{\sigma}_{j}=-\epsilon/3$ if both displacements are out or in, otherwise $+\epsilon/3$. 222 This term is the simplified description of the energies of the classical elastic model whose lattice sites are connected by springs, for which one can easily find that the two-in two-out configurations give the lowest energy, consistent with the experimental findings.

There are three parameters in this system; $\tilde{T}=k_{\rm B}T/J$ is the dimensionless temperature, $\tilde{\epsilon}=\epsilon/3J$ the ratio of the energy scales between the exchange interaction and the elastic energy of the displacement, and $\tilde{\delta}=\sqrt{6}\delta/3$ is the amplitude of the displacement (hereafter we call them simply as $T,~{}\epsilon,~{}\delta$). At $\epsilon\rightarrow\infty$, the lattice distortion becomes static. In the following analysis, we mainly focus on a representative system at $\delta=1.5,\epsilon=0.6$.

We consider the periodic systems of cubic geometry consisting of $L^{3}$ unit cells with totally $N=16L^{3}$ spins, and perform 120 statistically independent runs for the system size $L=4,5,6,8$, evaluating the averages and mean-squared errors of observable. To equilibrate the system we made a combined use of the exchange Monte Carlo method hukushima1996exchange; hukushima1999domain, conventional and loop melko2001long update for lattice variables, Metropolis-Reflection pawig1998monte and over-relaxation PhysRevB.53.2537 update for spin variables. However, to study the relaxational dynamics, we switched to single-spin-flip method after the equilibration. (see SI) We took $3\times 10^{7}$ Monte Carlo steps for both equilibration and for taking thermal averages $\langle...\rangle_{\rm eq}$.

To explore spin and orbital freezing we examined the auto-correlation functions (ACFs) for both degrees of freedom $C_{s}(t)=N^{-1}\sum_{i=1}^{N}\bm{S}_{i}(t)\cdot\bm{S}_{i}(0)$ and $C_{\sigma}(t)=N^{-1}\sum_{i=1}^{N}\bm{\sigma}_{i}(t)\cdot\bm{\sigma}_{i}(0)$ measured in equilibrium. The behavior of the spin and orbital ACFs are shown in Figs. 2(a) and 2(c). Both ACFs indicate the slowing down of the dynamics. As a quantitative measure of dynamics, we define the relaxation times $\tau_{s}(T)$ and $\tau_{\sigma}(T)$ for each degree of freedom as the time-scale such that the ACFs decrease down to 0.5, which are evaluated as shown in Figs. 2(b) and (d). Both relaxation times diverge at the same temperature $T_{c}\approx 0.07$ with power laws although their exponents differ. This suggests that the spin and orbital are frozen simultaneously. We checked that there is no finite size effect within the time scale which we analyzed. We also analyzed data at $\epsilon=0.65$ and obtained consistent results. (See SI). In the insets of Figs. 2(b) and 2(d) we show the scaling plots of ACFs assuming a scaling form, $C(t,T)=\tilde{C}(t/\tau(T))$, which is clearly satisfied. The scaling functions turn out to be well fitted by simple exponentials. We note that the relaxation function in the Edwards-Anderson models in the paramagnetic phase is more complicated ogielski1985dynamics; yoshino1993monte presumably due to Griffith singularity bray1988dynamics.

To get further insights we analyzed the static-structure factors of spins and lattice distortions defined respectively as $S_{s}(\bm{k})=N^{-1}|\sum_{i=1}^{N}e^{-i\bm{k}\cdot\bm{r}_{i}}\bm{S}_{1}\cdot\bm{S}_{i}|$ and $S_{\sigma}(\bm{k})=(4N)^{-1}|\sum_{i=1}^{4}\sum_{j=1}^{N}e^{-i\bm{k}\cdot(\bm{r}_{j}-\bm{r}_{i})}\bm{\sigma}_{i}\cdot\bm{\sigma}_{j}|$, where $\bm{k}$ is a wave vector. In figures 3(a) and 3(b) we show $S_{s}(h,h,l)$ measured above and below $T_{c}$, which show no notable differences. No Bragg peaks are observed bellow $T_{c}$, indicating that there is no magnetic long range order. In Figs. 3(c) and 3(d), we also find no sign of long range order in the lattice distortions. The pinch-points similar to those observed in spin-ice systems bramwell2001spin333 Note however that here we don’t multiply any scattering factor $|f(\bm{k})|$ to the structure factor. can be seen. The pinch-points reflect the fraction of tetrahedra that satisfy the ice-rule which increases smoothly with decreasing temperature without any anomaly at the transition temperature $T_{\rm c}$. (see SI).

The present SG transition can be detected by the two thermodynamic quantities. One is the nonlinear magnetic susceptibility,

[TABLE]

where $m_{\mu}=N^{-1}\sum_{i=1}^{N}S_{i\mu}$ and $h_{\mu}$ are the magnetization and the magnetic field, respectively, in the $\mu$-direction ($\mu=x,y,z$). The other is the nonlinear dielectric susceptibility,

[TABLE]

where $p_{\nu}=N^{-1}\sum_{i=1}^{N}\bm{\sigma}_{i}\cdot\hat{\bm{e}}_{\nu}$ is the dielectric polarization in the $\hat{\bm{e}}_{\nu}$-direction and $E_{\nu}$ is the electric field. Fortunately, the since Mo4+ ion has an electric charge, the fluctuation of the orbital glass ordering can be detected electrically. Thus, the two quantities are measurable in laboratories and allow a direct comparison between theories and experiments. In conventional theories mezard1987spin; fisher1988equilibrium, one often analyzes the SG susceptibility, $\chi_{\rm SG}$, which is proportional to $\chi_{3}$. Although $\chi_{\rm SG}$ is much easier to calculate with better accuracy in the presence of quenched disorder by using the overlap between two independent replicas, it is not convenient for our model, since the Hamiltonian is translationally invariant so that the replicas usually do not overlap.

Figure 4(a) shows $\chi_{3}$ as a function of temperature for several system size $L$. In overall, it takes the negative value, and from slightly above $T_{c}$, its amplitude starts to grow rapidly in lowering the temperature. This growth is enhanced for larger $L$, strongly suggesting that the divergence of $\chi_{3}$ remains in the thermodynamic limit, which is consistent with the experimental observation in Y2Mo2O7. Similar behavior is observed for $\chi_{3}^{(\sigma)}$, indicating the freezing of the orbital degrees of freedom into a disordered state.

Let us discuss the relevance of our results with experiments. The existence of the SG transition itself is already established in Y2Mo2O7 and its families, and although its origin has not been understood, a recent analysis showed that the lattice distortions possibly plays a certain role in the glass transition. Our model is constructed based on this finding by including the direct coupling between the lattice distortion and the spin. The two degrees of freedom simultaneously undergo a glass transition, indicating that they generate an effect of dynamical randomness to each other. Experiments can test our scenario by examining whether the nonlinear dielectric susceptibility diverges in the vicinity of $T_{\rm c}$, where already the nonlinear magnetic susceptibility shows the divergence. Previously, the orbital-glass transition was observed in the FeCr2S4, where the dielectric spectroscopy measurement played a major role, detecting the slowing down of orbital dynamics fichtl2005orbital. Similar techniques shall be applied to the pyrochlore materials and in addition to that the ZFC/FC anomaly gingras1997static shall be tested in linear dielectric response of the pyrhochlore materials.

Let us finally comment on the scenario often posed experimentally, which expects the freezing of the orbital degrees of freedom at a higher temperature than the SG transition booth2000local; thygesen2017orbital. Once the orbital configuration is frozen, our model is reduced to the standard EA model of Heisenberg spins with quenched randomness on a pyrochlore lattice, which is known to exhibit a SG transition saunders2007spin; shinaoka2011spin. Actually, if we increase $\epsilon$ in Eq.(1), the ice-rule configuration of orbitals is selected at the higher temperature, which is detected by the broad peak in the heat capacity(see SI). However, this peak does not mean the transition but a crossover; Below that temperature, the orbitals are still dynamically fluctuating among huge numbers of quasi-degenerate ice-rule configurations. In lowering the temperature, the orbital dynamics naturally slows down but without any thermodynamic anomaly much as usual spin ices jaubert2009signature. Interestingly, at a lower temperature the heat capacity shows a second peak (see SI) which should be attributed to coupling between the spin and lattice distortions. In our scenario, we anticipate the thermodynamic spin-orbital glass transition at around that temperature.

Some other effect such as the spin-orbit interaction shinaoka2013spin, the spin-phonon coupling and longer range interactions, which are not included in our model, may be needed to reproduce more precisely the experimental measurements, whereas our findings proved that the interplay of nearest-neighbor interactions and the dynamical JT distortions can solely drive the system to the glass transition. The present model can be further applied to other SG frustrated magnets such as Tb2Mo2O7 whose magnetic ion is JT active.

To summarize, we introduced the realistic model without quenched randomness consisting of spin and orbital degrees of freedom, that shows the thermodynamic glass transition for the first time in the finite dimensional periodic lattice. We performed the dynamical simulations in the equilibrium and found that the relaxation times of both spins and orbitals show a power-law divergence at the same temperature, $T_{\rm c}$. The nonlinear magnetic susceptibility and the nonlinear dielectric susceptibility together show negative divergence at around $T_{\rm c}$, signaling the simultaneous glass transition of these two degrees of freedom. The former is consistent with the experiments already performed in the pyrochlore oxide, Y2Mo2O7, whereas the latter can be tested in the near future to fix the long standing problem on the origin of the SG transition in this disorder free crystalline solid.

Appendix A Microscopic justification of the effective model

In a pyrochlore oxide Y2Mo2O7, magnetic ion Mo${}^{4+}(4d^{2})$ is affected by trigonal crystal field, and its triply $t_{2g}$ orbitals split into a single $a_{1g}$ orbital and doubly $e_{g}^{\prime}$ orbitals, where $a_{1g}$ has the lower energy than $e_{g}^{\prime}$. The spins of two electrons on Mo4+ form an effective $S=1$ following Hund’s rule. This is because the magnitude of energy splitting between $a_{1g}$ and $e_{g}^{\prime}$ is smaller than the strength of Hund’s coupling. solovyev2003effects We further take account of the Jahn-Teller lattice distortion, which is observed experimentally booth2000local; thygesen2017orbital. Ref.thygesen2017orbital shows that each Mo4+ moves away from the center of the tetrahedron, and accordingly, the Mo-O-Mo angle varies from its equilibrium position $\alpha=$127∘ to 116∘ for the in-in bond and 139∘ for the out-out bond. We checked that the change of the angle is large enough to give the different sign of the exchange interaction as in the following.

We consider an extended Kanamori Hamiltonian consisting of $d$-orbitals $a=\ket{a_{1g}},\ket{e_{g+}^{\prime}},\ket{e_{g-}^{\prime}}$ of Mo4+ ions and $p$-orbitals $b=\ket{p(X)},\ket{p(Y)},\ket{p(Z)}$ on the O2- ion. As shown in Fig. S6, we consider two adjacent Mo4+ ions, A and B, and the relative position of the degenerate three $p$-orbitals whose quantization axis is defined as such that it is line symmetric with respect to A and B. The Hamiltonian is given as

[TABLE]

where $c_{la\sigma}/c_{la\sigma}^{\dagger}$ are the annihilation/creation operators of the $a$-orbital of either the Mo-site ($l={\rm A},{\rm B}$) or O-site ($l=o$) with spin $\sigma$. The transfer integral $t_{a,b}$ in the first term are those between $a$-orbital of Mo$(4d)$ and $b$-orbital of O$(2p)$, which depends on both the distances between ions and the Mo-O-Mo angle $\alpha$ harrison2012electronic. We calculated $t_{a,b}$ based on the Slater-Koster parameter, as presented in Table. 6. The hopping between $p(Y)$ and any Mo$(4d)$ orbitals are relatively small and does not depend much on $\alpha$, thus are neglected.

The second term represents the excitation energy from O$(2p)$ to Mo$(4d)$ and we adopt $\Delta e_{a_{1g}}=2.0$ and $\Delta e_{e_{g\pm}^{\prime}}=2.5$ taken from the first principles band calculation in Ref.solovyev2003effects. The third and fourth terms are the on-site Coulomb interactions on Mo4+ and O2-, respectively. $K(\alpha,\alpha^{\prime})$ and $J(\alpha,\alpha^{\prime})$ are the Coulomb and exchange integrations which are described by the Slater-Condon parameters $F^{0},F^{2}$ and $F^{4}$. We adopt $F^{2}=4.52$ and $F^{4}/F^{2}=0.63$ solovyev2003effects. We leave $F^{0}=U_{\rm Mo}$ and $U_{o}$, which represent the strength of the on-site Coulomb interactions on the Mo and O ions respectively, as parameters.

We start from the strong coupling limit, $t_{a,b}\ll U_{\rm Mo},U_{o},\Delta e_{a}$, where the Mott insulating ground state consists of fully occupied O$(2p)$ orbitals and the ${\rm A/B}$ ions each having two electrons that occupy $a_{1g}$ and either of $e_{g+}^{\prime}$/$e_{g-}^{\prime}$ (See Fig. S1 (a)). Since these two electrons on the same Mo$(4d)$ ions form the $S=1$ triplet due to Hund’s coupling, there are four different spin configurations $(S^{z}_{\rm A},S^{z}_{\rm B})=(\uparrow,\uparrow),(\uparrow,\downarrow),(\downarrow,\uparrow),(\downarrow,\downarrow)$ which we take as a low energy basis. The leading term which gives the effective exchange interaction $J_{\rm eff}\bm{S}_{\rm A}\cdot\bm{S}_{\rm B}$ is obtained at the fourth of perturbation in terms of $t_{a,b}$. By considering all the perturbation processes of either of the two electrons on the $p$-oribals each to hop to the empty orbitals on A and B ions and coming back, one can evaluate the energy gain $E_{\rm F/AF}^{(4)}$ for the parallel and anti-parallel spin orientations of A and B ions. The effective exchange interaction is given by $J_{\rm eff}=(E^{(4)}_{\rm F}-E^{(4)}_{\rm AF})/3$, where denominator 3 is from the effective spin $S=1$. The fact that both $|t_{p(X),a_{1g}}|^{2}$ and $|t_{p(Z),e_{g\pm}^{\prime}}|^{2}$ vary significantly with $\alpha$ affects the relative energy gain of $E^{(4)}_{\rm F}$ and $E^{(4)}_{\rm AF}$, the sign of $J_{\rm eff}$ also changes its sign at $\alpha\sim 127^{\circ}$, as we present in Fig. 1(d) in the main text.

Appendix B Details of the Monte Carlo simulation

Here we explain the details of our Monte Carlo simulation. For updating lattice displacements $\sigma_{i}$ $(i=1,2,3,...,N)$ we used the conventional single-spin-flip Metropolis method. However, it is hard to investigate the thermodynamic properties of the ice-type system at low temperature due to multiple energy minima with high energy barrier. Therefore, we also adopted a nonlocal update method called the loop algorithm. melko2001long The nonlocal update consists of two steps: (1) We look for the loop consisting of only in-out bonds. (2) We flip all the lattice displacements on the loop simultaneously with the acceptance rate determined by the Metropolis rule. Note that the energy contribution from the second term of the Hamiltonian remains unchanged in this update. For updating Heisenberg spins we used the Metropolis Reflection method pawig1998monte and the over-relaxation method. PhysRevB.53.2537 In the Metropolis Reflection method, we prepare a plane dividing spin space into two equal parts randomly. To update a spin $\bm{S}_{i}$, a new candidate spin configuration $\bm{S}_{i}^{\prime}$ is created by reflecting $\bm{S}_{i}$ by the plane. More precisely the sign of the perpendicular component to plane is changed. The reflection is performed with probability determined by Metropolis rule. In a unit Monte Carlo step (MCS), we try the update for all spins $\bm{S}_{i}$ $(i=1,2,3,...,N)$ and apply the same plane to reflect all spins. In the over-relaxation method, we first calculate the effective magnetic field $\bm{h}_{\rm eff}$ around a spin $\bm{S}_{i}$, and then rotate the spin around $\bm{h}_{\rm eff}$ by an angle of $\pi$. A unit Monte Carlo step is consisted of the following steps.

Sequential single-spin flip for all lattice displacements with the conventional Metropolis method 2. 2.

$L$ times loop update for lattice displacements 3. 3.

Calculating all the coupling constants $J_{\bm{\sigma}_{i},\bm{\sigma}_{j}}$ 4. 4.

Sequential single-spin flip for all Heisenberg spins with the Metropolis Reflection method 5. 5.

$L$ times sequential single-spin flip for all Heisenberg spins with the over-relaxation method

Here $L$ is the system size defined in the main text.

The equilibrium states were realized by the replica exchange method. hukushima1996exchange; hukushima1999domain We prepared 48 replicas for $L=4,5,6$ and 72 replicas for $L=8$. We determined the maximum temperature and the minimum temperature as $T_{\rm max}=0.400$ and $T_{\rm min}=0.0625$ respectively. We optimized the temperature set so that the exchange rate is independent of $T$. hukushima1999domain We performed the trial to exchange the replicas with the interval 15 (MCS). Initial spin and lattice configuration for each replica were prepared a fully random configuration. MCSs for thermalization was determined as $3.0\times 10^{7}$ (MCS) for all system sizes. Observations of physical quantities are made during additional $3.0\times 10^{7}$ (MCS) to evaluate their thermal averages $\langle...\rangle_{\rm eq}$ by their MCS averages.

In dynamical simulations in equilibrium, the initial configurations ($t=0$ (MCS)) are created by above equilibration. In order to observe the natural time evolution we apply only the sequential single-spin flip for both degrees of freedom. A unit Monte Carlo step is consisted of the following steps.

Sequential single-spin flip for all lattice displacements with the conventional Metropolis method 2. 2.

Calculating all the coupling constants $J_{\bm{\sigma}_{i},\bm{\sigma}_{j}}$ 3. 3.

Sequential single-spin flip for all Heisenberg spins with the Metropolis Reflection method

Appendix C Definitions of the non-linear magnetic and dielectric

susceptibilities

Here we present the detailed definition of the non-linear susceptibilities we measured in our simulations. From fluctuation formulae we obtain the non-linear magnetic susceptibility as

[TABLE]

Within our Monte Carlo simulations of finite sized systems, odd order moments of the magnetization vanish due to the rotational symmetry, e.g. $\expectationvalue{m_{\mu}}_{\rm eq}=\expectationvalue{m_{\mu}^{3}}_{\rm eq}=0$. Hence we can evaluate it simply as,

[TABLE]

Similarly we evaluate the the non-linear dielectric susceptibility in the same manner as,

[TABLE]

Appendix D Spin and orbital time auto-correlation functions

In main text, we present the results of the time auto-correlation functions of the spin and orbital at $\epsilon=0.60$. The data suggest that spin and orbital degrees of freedom simultaneously freeze at $T_{\rm c}\approx 0.07$. Here we show an additional data set obtained at $\epsilon=0.65$. The scaling analysis of the data suggests a simultaneous glass transitions at $T_{\rm c}\approx 0.086$. Notably, the exponents $z\nu$ agree with those obtained at $\epsilon=0.6$ within the numerical accuracy, suggesting a universality.

Appendix E Heat capacity

We show below the temperature dependence of heat capacity for various $\epsilon$. We can see a broad peak in the vicinity of spin glass transition temperature $T_{\rm c}$, similarly to the observations in usual spinglass systems (for example see Chap 3.1.2 of mydosh2014spin) and the pyrochlore SG raju1992magnetic. One see that the dependence on the system size is very weak (see Fig. S3 (a)).

For larger $\epsilon$, another broad peak emerges at higher temperatures (see Fig. S3 (c)), which can be interpreted as a reflection of the crossover from purely random phase to ice like phase. (see Fig. S4.)

Appendix F Fraction of ice like lattice distortions

Fig. S4 shows the temperature-dependence of the fraction of ice like 2-in 2-out structure of the lattice distortions for various $\epsilon$. It can be seen that the onset of ice like structure is gradual. The 2ndary peak of the heat-capacity observed at large $\epsilon$ (See Fig. S3(c)) can be interpreted as due to the onset of the ice like structures.

Appendix G Non-linear magnetic and dielectric susceptibilities

In the main text, we present the temperature dependence of magnetic and dielectric non-linear susceptibilities for $\epsilon=0.60$. Here, we show them for various $\epsilon$.