Control of effective free energy landscape in a frustrated magnet by a field pulse

TL;DR

This paper demonstrates how a magnetic field pulse can manipulate the free energy landscape of a frustrated magnet, specifically a pyrochlore XY antiferromagnet, by exciting fast modes to control state stability.

Contribution

It introduces a method to control the non-equilibrium free energy landscape in frustrated magnets using field pulses, linking fast mode excitation to slow ground state dynamics.

Findings

Field pulses can create tunable free energy minima at unstable states.

Fast mode excitation influences slow ground state motion.

Control of non-equilibrium landscapes in frustrated magnets.

Abstract

Thermal fluctuations can lift the degeneracy of a ground state manifold, producing a free energy landscape without accidentally degenerate minima. In a process known as order by disorder, a subset of states incorporating symmetry-breaking may be selected. Here, we show that such a free energy landscape can be controlled in a non-equilibrium setting as the slow motion within the ground state manifold is governed by the fast modes out of it. For the paradigmatic case of the classical pyrochlore XY antiferromagnet, we show that a uniform magnetic field pulse can excite these fast modes to generate a tunable effective free energy landscape with minima at thermodynamically unstable portions of the ground state manifold.

Figures

Equations

$$H=$$

$$-J^{zz}{S}_i^z{S}_j^z.$$

$${\ddot{\varphi }}_l=-{\omega }_0^2{f}_l(\varphi ){\varphi }_l^2;\ddot{\varphi }=-{\tilde{\omega }}_0^2\sum _{l=a,b,c}{\varphi }_l^2\sin (2\varphi +{\Theta }_l),$$

$$V(\varphi )=\frac{3(1-j_{zz})}{\sqrt{3+j_{zz}}}\sum _{l=a,b,c}{I}_l\sqrt{f_l(\varphi )},$$

$$\frac{M_{a,b,c}}{\mathcal{N}{\mu }_{\mathrm{B}}S}=-\frac{2g_{\parallel }}{\sqrt{3}}{z}_{a,b,c}-\frac{4g_{\perp }}{\sqrt{6}}\sin (\varphi -{\Theta }_{a,b,c}){\varphi }_{a,b,c},$$

$${\hat{x}}_0$$

$${\hat{y}}_0$$

$${\hat{z}}_0$$

$${\psi }_{01}={\psi }_{23}=0;{\psi }_{02}={\psi }_{13}=\frac{2\pi }{3};{\psi }_{03}={\psi }_{12}=-\frac{2\pi }{3}.$$

$$L=\frac{\mathcal{N}}{4}\sum _i{z}_i\dot{{\varphi }_i}-E,$$

$$E=-\frac{\mathcal{N}}{2}\sum _{i>j}\sqrt{(1-z_i^2)(1-z_j^2)}[\cos ({\varphi }_i-{\varphi }_j)-j_{\pm \pm }\cos ({\varphi }_i+{\varphi }_j+{\psi }_{ij})]+j_{zz}{z}_i{z}_j.$$

$$L=L_{\varphi }+L_a+L_b+L_c,$$

$$L_{\varphi }$$

$$L_a$$

$$L_b$$

$$L_c$$

$${\ddot{\varphi }}_a$$

$${\ddot{\varphi }}_b$$

$${\ddot{\varphi }}_c$$

$$\ddot{\varphi }$$

$$\overline{{\varphi }_a^2}=\frac{I_a}{{\omega }_0\sqrt{1-j_{\pm \pm }/2\cos 2\varphi }};\overline{{\varphi }_b^2}=\frac{I_b}{{\omega }_0\sqrt{1-j_{\pm \pm }/2\cos (2\varphi +2\pi /3)}};\overline{{\varphi }_c^2}=\frac{I_c}{{\omega }_0\sqrt{1-j_{\pm \pm }/2\cos (2\varphi -2\pi /3)}}.$$

$$\ddot{\varphi }=-\frac{\partial V(\varphi )}{\partial \varphi },$$

$$V(\varphi )=\frac{3(1-j_{zz})}{\sqrt{3+j_{zz}}}\left(I_a\sqrt{1-\frac{j_{\pm \pm }}{2}\cos 2\varphi }+I_b\sqrt{1-\frac{j_{\pm \pm }}{2}\cos (2\varphi +\frac{2\pi }{3})}+I_c\sqrt{1-\frac{j_{\pm \pm }}{2}\cos (2\varphi -\frac{2\pi }{3})}\right).$$

$$\frac{E}{\mathcal{N}{J}_{\pm }{S}^2}=3(1-j_{zz})z^2+\frac{V(\varphi )}{6(1-j_{zz})}=\mathrm{constant}.$$

$$\frac{2\pi \hslash }{4S\sqrt{J_{\pm }(3J_{\pm }+J_{zz})}}\approx \frac{2\pi \hslash }{0.5\mathrm{meV}}\approx 8\mathrm{ps}.$$

$${\dot{\varphi }}_a$$

$${\dot{\varphi }}_b$$

$${\dot{\varphi }}_c$$

$$\dot{\varphi }$$

$${\varphi }_a(0^{+})=2n_a{a}_{\parallel };{\varphi }_b(0^{+})=2n_b{a}_{\parallel };{\varphi }_c(0^{+})=2n_c{a}_{\parallel };\varphi (0^{+})=0;$$

$$z_a(0^{+})=0;z_b(0^{+})=\sqrt{3}{n}_b{a}_{\perp };z_c(0^{+})=-\sqrt{3}{n}_c{a}_{\perp };z(0^{+})=0.$$

$${\dot{\varphi }}_a(0^{+})=0;{\dot{\varphi }}_b(0^{+})=2\sqrt{3}(3+j_{zz})n_b{a}_{\perp };{\dot{\varphi }}_c(0^{+})=-2\sqrt{3}(3+j_{zz})n_c{a}_{\perp };\dot{\varphi }(0^{+})=0.$$

$$I_{a,b,c}=\frac{1}{2{\omega }_{a,b,c}}{\dot{\varphi }}_{a,b,c}^2+\frac{{\omega }_{a,b,c}}{2}{\varphi }_{a,b,c}^2;{\omega }_a={\omega }_0\sqrt{1-\frac{j_{\pm \pm }}{2}};{\omega }_b={\omega }_0\sqrt{1+\frac{j_{\pm \pm }}{4}};{\omega }_c={\omega }_0\sqrt{1+\frac{j_{\pm \pm }}{4}}.$$

$$\displaystyle\left(\begin{array}[]{c}S^{x}_{i}(0^{+})\\ S^{y}_{i}(0^{+})\\ S^{z}_{i}(0^{+})\end{array}\right)=R(\frac{\mathbf{h}}{h},h\tau)\left(\begin{array}[]{c}S^{x}_{i}(0^{-})\\ S^{y}_{i}(0^{-})\\ S^{z}_{i}(0^{-})\end{array}\right).$$

Full text

The emergence of a thermodynamic landscape from microscopic interactions, and its consequences for macroscopic behavior, is a central theme of condensed matter physics. The essential role of fluctuation contributions to the free energy landscape is underlined by entropic interactions in soft matter Kirkwood (1939); Guth and James (1941); Onsager (1949), infinite temperature phase transitions Alder and Wainwright (1957); Frank H. Stillinger Jr. and Zevi W. Salsburg (1967), and, at low temperatures, the phenomenon of order by disorder (ObD) Villain et al. (1980); Shender (1982); Kawamura (1984); Henley (1989); Moessner and Chalker (1998).

This raises the question of how one can control a free energy landscape in a condensed matter system by manipulating fluctuations instead of changing its Hamiltonian. In equilibrium (EQ), the Boltzmann distribution fully dictates the free energy landscape and hence allows little room for control. In this work, we show that such control is feasible if one drives the system out of equilibrium, where fluctuations dynamically produce a free energy landscape which may be tuned by changing the non-equilibrium conditions.

To illustrate this, we sketch a simple but generic scenario. It is based on a crisp distinction between energetic and entropic contributions to the free energy as exists at low energies in a broad class of geometrically frustrated magnets. This arises because geometrical frustration often results in an accidentally degenerate ground state (GS) manifold, i.e. a large, continuous family of degenerate GSs not related by any symmetry operation. In EQ, the entropy due to the thermal fluctuations near each GS can vary along the manifold. Parametrizing the ground state manifold by coordinates $\mathcal{Q}$, one thus obtains an effective free energy $V(\mathcal{Q})$ from integrating out the fluctuations, which can thus lift the accidental degeneracy and tend to stabilize the GS(s) with maximal entropy.

We focus on systems with Hamiltonian dynamics such as frustrated vector spin models. We separate variables into the pseudo-Goldstone modes describing the drift motion within the GS manifold, and the other normal modes involving deviations out of it (Fig. 1) Moessner (2001). In the low temperature limit, the drift motion within GS manifold is vanishingly slow as they experience no linear restoring force. By contrast, the other modes are fast owing to their finite stiffness. Integrating these out (a controlled procedure provided their frequencies are bounded above zero) yields a non-equilibrium contribution to the effective free energy of the form $V(\mathcal{Q})\propto\sum_{i}I_{i}\sqrt{K_{i}(\mathcal{Q})}$, with action variable $I_{i}>0$ of the $i$-th excited normal mode, and $K_{i}$ its $\mathcal{Q}$-dependent stiffness Douçot and Simon (1998). Importantly, $V(\mathcal{Q})$ can be tuned by adjusting the weights $I_{i}$. This we show can be achieved straightforwardly by a magnetic field pulse.

We flesh out this scenario with a concrete model, for which we establish properties of the effective landscape, its tunability and resulting dynamics, as well as its interplay with the thermal fluctuations. We also develop an analytic toy model that transparently explains the central phenomenon.

Consider the classical pyrochlore XY antiferromagnet, which has gained prominence thanks to its likely realization of ObD in the rare earth magnet Er2Ti2O7 Champion et al. (2003); Champion and Holdsworth (2004); Zhitomirsky et al. (2012); Savary et al. (2012); Wong et al. (2013); Oitmaa et al. (2013); McClarty et al. (2014); Zhitomirsky et al. (2014):

[TABLE]

Here, spin $\mathbf{S}_{i}$ of length $S$ resides on pyrochlore lattice site $i$. The summation runs over nearest-neighbor bonds. $S^{x,y,z}_{i}$ are the Cartesian components of the spin in the local frame $\{\hat{x}_{i},\hat{y}_{i},\hat{z}_{i}\}$, and $S^{\pm}_{i}\equiv S^{x}\pm iS^{y}_{i}$. $\psi_{ij}$ are link-dependent phase angles due to local spin frames (See SM for details). $J_{\pm},J_{zz}$ parametrize anisotropic Heisenberg exchange interaction, whereas $J_{\pm\pm}$ originates from a Dzyaloshinskii-Moriya or pseudo-dipolar interaction. We assume $J_{\pm}>J_{\pm\pm},J_{zz}>0$. For simplicity, we have omitted a symmetry-allowed term known as $J_{z\pm}$ as it will not change the physics discussed in this work.

We briefly review the ObD following from Eq. (1) Zhitomirsky et al. (2012); Savary et al. (2012); Wong et al. (2013); Oitmaa et al. (2013); McClarty et al. (2014); Zhitomirsky et al. (2014). The GSs, which preserve lattice translation symmetry, show Néel order with complex order parameter $O\equiv\sum_{i}S^{+}_{i}/(\mathcal{N}S)$ ($\mathcal{N}$ is the number of sites). GS spin configurations are $\mathbf{S}_{i}=S(\cos\phi\hat{x}_{i}+\sin\phi\hat{y}_{i})$ with $\phi=\arg O$ (Fig. 2a). Crucially, the GS energy is independent of $\phi$: Eq. (1) possesses an accidental $U(1)$ degeneracy although it only exhibits discrete symmetries.

At small but finite temperature, the entropy due to the spin wave fluctuations lifts the accidental degeneracy, yielding six symmetry-equivalent maxima known as $\psi_{2}$ states, located at $\phi=m\pi/3$, $m\in\mathbb{Z}$: thermal ObD corresponds to an emergent entropic six-state clock anisotropy in the $U(1)$ manifold.

To study the non-EQ dynamics of Eq. (1), we endow the spins with precessional Landau-Lifshitz (LL) dynamics. We begin with a toy model that captures the essential features: taking the spins on each of the four sublattices $i=0\sim 3$ with the same orientation, $\mathbf{S}_{i}/S=\sqrt{1-z^{2}_{i}}(\cos\phi_{i}\hat{x}+\sin\phi_{i}\hat{y})+z_{i}\hat{z}_{i}$ reduces the degrees of freedom to the four spins in a unit cell. In other words, we freeze spin wave modes with wave vector $\mathbf{k}\neq 0$.

We define $\phi\equiv(\phi_{0}+\phi_{1}+\phi_{2}+\phi_{3})/4$, $\phi_{a}\equiv(\phi_{0}+\phi_{1}-\phi_{2}-\phi_{3})/2$, $\phi_{b}\equiv(\phi_{0}+\phi_{2}-\phi_{1}-\phi_{3})/2$, and $\phi_{c}\equiv(\phi_{0}+\phi_{3}-\phi_{1}-\phi_{2})/2$, along with corresponding $z$ and $z_{a,b,c}$. $\phi$ corresponds to the pseudo-Goldstone mode, whereas $\phi_{a,b,c}$ correspond to finite-frequency optical magnons. Linearizing Eq. (1) in $\phi_{a,b,c}$ yields SM :

[TABLE]

with rescaled variables such that $J_{\pm}S\to 1$, $j_{\pm\pm}\equiv J_{\pm\pm}/J_{\pm}$, $j_{zz}\equiv J_{zz}/J_{\pm}$. $l$ runs over labels $a,b,c$. $f_{l}(\phi)=1-j_{\pm\pm}/2\,\cos(2\phi+\Theta_{l})$. $\Theta_{a}\equiv 0$, $\Theta_{b}\equiv 2\pi/3$, and $\Theta_{c}\equiv-2\pi/3$. $\omega_{0}\equiv 4\sqrt{3+j_{zz}}$ and $\widetilde{\omega}_{0}\equiv\sqrt{6j_{\pm\pm}(1-j_{zz})}$ are constant frequencies.

Since $\phi_{l}\ll 1$, the optical magnons are fast harmonic oscillators parametrically driven by the slow motion of the pseudo-Goldstone mode $\phi$. We thus proceed to integrate out optical magnons by using the method of averaging Arnol’d (1989); Douçot and Simon (1998). We obtain SM : $\ddot{\phi}=-\partial{}V(\phi)/\partial\phi$, where

[TABLE]

is the effective potential in the $U(1)$ degenerate manifold. $I_{l}$, the action variable of mode $l$, is an adiabatic invariant Landau and Lifshitz (1976), so that $V(\phi)$ is approximately time independent.

For $I_{a}=I_{b}=I_{c}$, $V(\phi)\propto-j^{3}_{\pm\pm}\cos(6\phi)+o(j^{3}_{\pm\pm})$ with minima at the states selected by ObD in EQ. For generic values of $I_{a,b,c}$, however, the minima are two-fold degenerate, $V(\phi)=V(\phi+\pi)$ and, crucially, located elsewhere.

Selectively exciting optical magnons thus permits control of the individual values of $I_{a,b,c}$ and thereby $V(\phi)$. This can in fact be achieved simply via the polarization of an applied magnetic field pulse. To see this, note that the magnetization $M_{a,b,c}$ along the crystallographic $a,b,c$ axes, to leading order in $\phi_{a,b,c}$ and $z_{a,b,c}$ is Maryasin et al. (2016); SM ,

[TABLE]

with $g_{\parallel}(g_{\perp})$ being the Landé g-factor along the local $\hat{z}$ axis ($\hat{x},\hat{y}$ axes): the $a,b,c$ magnons respectively carry magnetic dipole moments in the three crystallographic axes.

We next confirm numerically such generation and control of an effective free energy landscape. We first equilibrate the model at temperature $T$ and then apply a short magnetic pulse whose temporal profile is a Dirac-$\delta$ function: $\mathbf{B}(t)=B_{\mathrm{max}}\hat{n}\delta(t/\tau)$, where $B_{\mathrm{max}}$ is the peak strength, $\tau$ is the duration, and $\hat{n}$ is the polarization. After the pulse, we remove the bath and let the system evolve according to the LL equation at $t>0$.

The simulation uses a periodic lattice of $8\times 8\times 8$ unit cells. We do not observe significant system size dependence in dynamics SM . We integrate the LL equation via the 4th order Runge-Kutta method (RK4). For nonzero $T$, we generate $2^{10}$ initial states for the LL equation from canonical Monte Carlo (MC) and average over microcanonical trajectories. The RK4 step width is chosen such that the relative error of energy $\epsilon<10^{-5}$. The integration stops at $10^{3}\hbar/(J_{\pm}S)$, corresponding to $>10^{3}$ oscillation cycles of optical magnons. Model parameters $j_{\pm\pm}=0.646$, $j_{zz}=0.192$, $g_{\parallel}=2.45$, and $g_{\perp}=6.0$ are similar to those for Er2Ti2O7 Savary et al. (2012). In EQ, the $\psi_{2}$ states occur with equal probability. We assume the model is initially in a single domain with $\phi=0$. Only the area of the $\delta$ peak $B_{\mathrm{max}}\tau$ enters the equation of motion, which we set to 0.1 ps$\cdot$T SM . The pulse polarization lies in the $ab$ plane, i.e. $\hat{n}=\cos\theta_{\mathrm{B}}\hat{a}+\sin\theta_{\mathrm{B}}\hat{b}$ (Fig. 2a, inset). In this setup, the pulse only excites $a$ and $b$ magnons. The energy deposited by the pulse is $\sim 10^{-3}\,J_{\pm}S^{2}$ per spin.

First consider the $T=0$ limit. $\phi$ initially rests at [math] for $t<0$. The pulse generates a $V(\phi)$ whose minima are located elsewhere. $\phi$ thus oscillates around a nearby minimum of $V(\mathrm{\phi})$ (see Fig. 2b for $\theta_{\mathrm{B}}=0.15\pi$). To extract $V(\phi)$ from simulation, we use the conservation of energy SM , $V(\phi)/(6(1-j_{zz}))+K=\mathrm{const.}$, where the first and second terms are respectively the potential and kinetic energy density of the pseudo-Goldstone mode. $K$ can be evaluated from data. The result is in good agreement with Eq. (3), with $I_{a,b,c}$ extracted from the initial condition (Fig. 2c).

We next show that the accidental degeneracy underpinning ObD is constitutive to the controllability of the free energy landscape in our present scheme. To do this, we contrast Eq. (1) with a clock model where ObD is mimicked by microscopic interactions yielding the same degeneracy lifting Rau et al. (2016): we add a six-state clock anisotropy term, $-\Delta/S^{4}[(S^{+}_{i})^{6}+h.c.]$, to Eq. (1) and set $j_{\pm\pm}=0$. For the clock model, $V(\phi)\propto-\Delta\cos(6\phi)$, and the pulse only produces a small renormalization of $\Delta$: for $\Delta/J_{\pm}=10^{-4}$, $\phi$ oscillates around [math] (Fig. 2b). The oscillation amplitude is small as the Néel order couples to the magnetic field only nonlinearly.

We further demonstrate the control on $V(\phi)$ by scanning the field polarization $\hat{n}$. As $\hat{n}$ rotates from the $a$ to the $b$ axis, $I_{b}$ increases while $I_{a}$ decreases. The $I_{a}$ term of Eq. (3) favors $\phi=0,\pi$, whereas the $I_{b}$ term favors $\phi=-\pi/3,2\pi/3$. Thus, the minimum positions $\phi_{\mathrm{min}}$ continuously shift from $0,\pi$ to $-\pi/3,2\pi/3$. In numerical simulation, the shift of $\phi_{\mathrm{min}}$ is manifest as the change of the center position of $\phi$ oscillation (Fig. 2d). The small discrepancy between the toy model and the simulation is likely due to nonlinear effects neglected in Eq. (2).

At nonzero temperature, the effective potential produced by optical magnons is also subject to the thermal fluctuations in $\mathbf{k}\neq 0$ modes. In EQ, these modes produce an entropic effective potential with the magnitude of $\sim 10^{-3}\,k_{B}T$ per spin. We now show that the above non-EQ phenomena are nonetheless robust at low $T$. We start with $k_{B}T/(J_{\pm}S^{2})=10^{-3}$. Fig. 3a,b show two representative cases. For $\theta_{\mathrm{B}}=0.1\pi$, $\mathrm{arg}\langle{}O\rangle$, the argument of the average Néel order parameter, exhibits a persistent oscillation with small damping. On the other hand, for $\theta_{\mathrm{B}}=0.5\pi$, in addition to oscillation in the argument, the modulus $|\langle{}O\rangle|$ decreases quickly. However, the Néel order persists as the average modulus squared $\langle|O|^{2}\rangle\approx 1$ within the simulation time window (Fig. 3b). Thus, the decrease in $|\langle O\rangle|$ is due to the loss of anisotropy, i.e. the argument of $O$ explores the $U(1)$ degenerate manifold. This is made clear by the histogram of $\mathrm{arg}\,O$ for $\theta_{\mathrm{B}}=0.5\pi$, which spreads out in the $U(1)$ manifold at late time (Fig. 3d).

The difference between these two types of behavior originates dynamically as follows. Consider the $a$ and $b$ magnons underpinning $V(\phi)$. Their ensemble-averaged action variables $\langle I_{a,b}\rangle$ gradually decay in time, indicating that the energy is being slowly transferred to other modes (Fig. 3e). Note $I_{a,b,c}$ must equalize and $\phi$ must return to $\psi_{2}$ positions after the thermal equilibration time $\tau_{\mathrm{eq}}$. Yet, the persistent difference in $I_{a,b}$ indicates $\tau_{\mathrm{eq}}>10^{3}\hbar/{J_{\pm}S}$. Comparing $\langle{}I_{a,b}\rangle$ for $\theta_{\mathrm{B}}=0.1\pi$ with $\theta_{\mathrm{B}}=0.5\pi$, the latter exhibits more pronounced decrease. This slow decay implies that of $V(\phi)$ as well (Eq. (3)). Now, the motion of the pseudo-Goldstone mode $\phi$ is oscillatory as long as its total energy $E$ is less than the potential maxima $V_{\mathrm{max}}$. As $V_{\mathrm{max}}$ decreases, $\phi$ may therefore overcome the potential barrier and enter an open orbit. As a qualitative test of this picture, we extract $E$ and $V_{\mathrm{max}}$ from data and count the number of trajectories for which $E>V_{\mathrm{max}}$ (Fig. 3f). For $\theta_{\mathrm{B}}=0.1\pi$, the number of such trajectories is negligible, whereas the count steadily grows in time for $\theta_{\mathrm{B}}=0.5\pi$.

Having gained a qualitative understanding the two representative cases, we consider the systematic $\theta_{\mathrm{B}}$-dependence of the order parameter dynamics (Fig. 4a). Throughout, $\mathrm{arg}\langle{}O\rangle$ exhibits damped oscillation, while $|\langle{}O\rangle|$ decreases at an larger rate as $\theta_{\mathrm{B}}$ approaches $\pi/2$. Similar to $T=0$, the center of $\mathrm{arg}\langle{}O\rangle$ oscillation gradually shifts from 0 to $-\pi/3$. The center position at early times agrees with the potential minima of Eq. (3) with $I_{a,b,c}$ extracted from initial conditions. The small drift in oscillation center at late time is due to a small change in the relative weight of $I_{a,b,c}$ and the increasing importance of thermal fluctuations. The higher oscillation frequency at larger $\theta_{\mathrm{B}}$ results from the larger curvature of $V(\phi)$. The beating of the $\mathrm{arg}\langle{}O\rangle$ observed at $\theta_{\mathrm{B}}=0.5\pi$ is likely due to a finite-width distribution of the oscillation frequencies.

Finally, we study the temperature dependence of this phenomenon. Fig. 4b shows the order parameter dynamics for fixed $\theta_{\mathrm{B}}=0.15\pi$. Damping of the oscillation in $\mathrm{arg}\langle{}O\rangle$ increases with temperature $T$, while the oscillation center moves toward 0. Both signify the growing influence of thermal fluctuations. At $k_{B}T/(J_{\pm}S^{2})=0.2$, $\mathrm{arg}\langle{}N\rangle$ is nearly stationary and very close to [math], which suggests the thermal fluctuations dominate over the non-EQ potential from optical magnons. This is clearly seen in the histogram of $\mathrm{arg}\,O$ (Fig. 4c): initially at a single $\psi_{2}$ state, the system escapes to other $\psi_{2}$ states at late time. Since the model remains XY-ordered at this temperature and hence no symmetry changes occur, we infer a crossover at temperature $T^{\ast}$, above which the landscape essentially reverts to being thermal. Simple dimensional analysis indicates $T^{\ast}$ scales with the energy density deposited by the pulse.

Looking ahead, our identification of a tunable effective free energy landscape leads to a number of interesting directions for future research. Firstly, a more detailed analysis of the finite temperature dynamics is needed for a complete picture of the aforementioned crossover. Secondly, while we have focused on classical models, it is intriguing to see to what extent the physics discussed here has a natural analogue in the quantum realm. This would be directly applicable to quantum magnets such as Er2Ti2O7 and NaCaCo2F7 Krizan and Cava (2015); Ross et al. (2016); Sarkar et al. (2016), which hold the promise of allowing a detailed and quantitative study of thermal and/or quantum ObD effects.

Appendix A Local spin frames and phase factors

We define the following local spin frames for the four sublattices, labeled 0, 1, 2, and 3. All vectors are given in terms of Cartesian components in the crystallographic axes $\hat{a},\hat{b},\hat{c}$:

[TABLE]

In particular, $\hat{z}_{i}$ is parallel with the 3-fold axis that passes through site $i$.

The phase factors $\psi_{ij}$ arise due to the above local spin frames. $\psi_{ij}=\psi_{ji}$. Thanks to translation symmetry, $\psi_{ij}$ depends only on the sublattice labels of site $i$ and $j$. There are six independent $\psi_{ij}$ values:

[TABLE]

Appendix B Equations of motion of the toy model

It is convenient to use the Lagrangian formalism. The Lagrangian of the toy model is given by:

[TABLE]

Here, $i$ and $j$ are sublattice indices running from 0 to 3. We have rescaled time such that $J_{\pm}S\to 1$. $j_{\pm\pm}\equiv J_{\pm\pm}/J_{\pm}$, and $j_{zz}\equiv J_{zz}/J_{\pm}$. $\mathcal{N}$ is the number of spins.

Using the new canonical variables $\phi$, $\phi_{a,b,c}$, $z$, and $z_{a,b,c}$ defined in the main text, and expanding $E$ to the quadratic order in all variables except for $\phi$, we find,

[TABLE]

Using Euler-Lagrangian equations, we obtain,

[TABLE]

Here, $\omega_{0}\equiv 4\sqrt{3+j_{zz}}$ is a constant frequency.

To derive an effective equation of motion for $\phi$, we approximate $\phi^{2}_{a,b,c}$ in Eq. (9d) by their respective time average $\overline{\phi^{2}}_{a,b,c}$ in an oscillation cycle and, then, relate them to the action variables $I_{a,b,c}$:

[TABLE]

Substituting the above into Eq. (9d), we find an effective equation of motion for $\phi$:

[TABLE]

This is the result given in the main text.

The above equation of motion for $\phi$ implies an effective energy conservation law:

[TABLE]

We recognize the first term as the kinetic energy of the pseudo-Goldstone mode and the second term as the effective potential energy due to the excited optical magnons.

Appendix C Modeling the magnetic pulse

In this work, we model the magnetic pulse as a Dirac-$\delta$ function for the sake of simplicity: $\mathbf{B}(t)=B_{\mathrm{max}}\hat{n}\delta(t/\tau)=B_{\mathrm{max}}\tau\hat{n}\delta(t)$, where $\tau$ is the duration of the pulse, $B_{\mathrm{max}}$ is the peak field strength, and $\hat{n}$ is the unit polarization vector. In our simulation, we chose the pulse area $B_{\mathrm{max}}\tau=0.1\mathrm{T}\cdot\mathrm{ps}$.

Using a Dirac-$\delta$ function is a sound approximation provided that the duration of the pulse $\tau$ is much shorter than oscillation period of the optical magnons. From Eq. (9), omitting a small dependence on $\phi$, the oscillation frequency of the optical magnon is $\omega_{0}=4\sqrt{3+j_{zz}}$ . Therefore, the oscillation period is approximately,

[TABLE]

Here, we have restored the dimension of time. In the second equality, we have used the model parameters similar to Er2Ti2O7. Setting a short duration $\tau=1\mathrm{ps}$ and the peak strength typical for THz laser $B_{\mathrm{max}}=0.1\mathrm{T}$ yields the peak area $B_{\mathrm{max}}\tau=0.1\mathrm{T}\cdot\mathrm{ps}$.

We now analyze the impact of the Dirac $\delta$ pulse on spin dynamics. We begin with the toy model. The equations of motion are modified in the presence of the $\delta$ pulse:

[TABLE]

Note that we have omitted the regular terms that do not involve $\delta(t)$. The dimensionless parameters $a_{\perp}\equiv 2g_{\perp}\mu_{0}B_{\mathrm{max}}\tau/\sqrt{6}$ and $a_{\parallel}\equiv g_{\parallel}\mu_{0}B_{\mathrm{max}}\tau/\sqrt{3}$. $n_{a,b,c}$ are the Cartesian components of the polarization vector $\hat{n}$ in the crystallographic (cubic) axis. Assuming $\phi_{a,b,c}(0^{-})=\phi(0^{-})=z_{a,b,c}(0^{-})=z(0^{-})=0$ before the pulse, immediately after the pulse, we have

[TABLE]

The above serve as the initial conditions for the ensuing free evolution govern by Eq. (9). In particular, one may evaluate the action variables $I_{a,b,c}$ by using the following:

[TABLE]

Evidently, $I_{a,b,c}$ are respectively related to the values of $\phi_{a,b,c}$ and $\dot{\phi}_{a,b,c}$, which in turn depend on the polarization vector $\hat{n}$.

In simulating the full lattice model, the Dirac-$\delta$ magnetic pulse is treated in the same spirit. Let $\mathbf{S}_{i}=S^{x}_{i}\hat{x}_{i}+S^{y}_{i}\hat{y}_{i}+S^{z}_{i}\hat{z}_{i}$ be the spin vector on a pyrochlore lattice site $i$. $S^{x,y,z}_{i}$ are the Cartesian components in the local spin frame $\{\hat{x}_{i},\hat{y}_{i},\hat{z}_{i}\}$. The pulse induces a pure rotation:

[TABLE]

Here, $S^{x,y,z}_{i}(0^{\pm})$ are the Cartesian components immediately after and before the pulse. $R(\hat{u},\theta)$ is the $3\times 3$ Euclidean rotation matrix with respect to the rotational axis $\hat{u}$ by an angle $\theta$. The vector $\mathbf{h}=(g_{\perp}\mu_{0}B\hat{n}\cdot\hat{x}_{i},g_{\perp}\mu_{0}B\hat{n}\cdot\hat{y}_{i},g_{\parallel}\mu_{0}B\hat{n}\cdot\hat{z}_{i})$ is the Zeeman field in the local spin frame, and $h$ is its length.

Appendix D Monte Carlo simulation

In our finite temperature molecular dynamics simulation, we use canonical Monte Carlo (MC) method to generate initial conditions for the Landau-Lifshitz equation. In the MC simulation, we equilibrate the system with $2^{15}$ MC steps. We skip $2^{6}$ MC steps before we record a new configuration as a sample in order to reduce the sample correlation. In thermal equilibrium, the six $\psi_{2}$ states occur with equal probability, giving rise to six magnetic domains. We assume the system is initially in a single domain of the $\psi_{2}$ states with $\phi=0$. This is done by discarding samples whose order parameter argument $\mathrm{arg}\,O$ falls out of the interval $\mathrm{arg}\,O\in[-\pi/6,\pi/6)$. We generate $2^{10}$ samples and feed them to the Runge-Kutta integrator as the initial conditions.

The details of the MC update scheme are given as follows. Each MC step consists of 1 spin flip update and 8 over-relaxation updates:

• •

To flip a spin, we draw with uniform probability a new orientation from a spherical cap that centers on the previous orientation of the spin. This procedure is performed on each and every spin in a spin flip update. The size of the cap is a parameter of the update, which we adjust to reduce the autocorrelation time. The updated configuration is accepted according to the Metropolis rule.

• •

An over-relaxation update consists of four steps of rotations. In each step, we rotate each and every spin belonging to a given sublattice $i$ with respect to its local molecular field by angle $\pi$. The rotations within the same sublattice are independent from each other. As a result, the precise order of the rotations within a sublattice does not affect the outcome. However, rotations for two different sublattices do not commute. We first rotate the spins belonging to sublattice 0, and then proceed to sublattice 1, 2, and 3. The updated configuration is always accepted. By construction, the over-relaxation update conserves energy and satisfies detailed balance.

We benchmark our MC update scheme by measuring the low temperature thermodynamic properties of the model. To do this, we perform a separate MC simulation. The results are summarized in Fig. 5, where the data are obtained by averaging over $2^{20}$ successive MC steps after equilibrating the system for $2^{18}$ MC steps. Note we no longer discard MC samples that fall outside of the $\psi_{2}$ domain with $\phi=0$.

In agreement with the equipartition theorem, the low temperature specific heat clearly approaches $k_{B}$ (Fig. 5a). To characterize thermal order by disorder, we use two order parameters: $\rho\equiv\sqrt{\langle|O|^{2}\rangle}$, and $\zeta\equiv\langle\cos 6\phi\rangle$, where $|O|$ and $\phi$ are respectively the modulus and argument of the complex Néel order parameter $O$, and $\langle\cdots\rangle$ denotes the ensemble average. $\rho$ characterizes the Néel order, whereas $\zeta$ characterizes the six-state clock order. In particular, $\zeta=1$ if the system falls precisely on a $\psi_{2}$ state. Fig. 5b shows the Néel order is well established in the temperature window $k_{B}T<J_{\pm}S^{2}$. As $T$ decreases, the six-state clock order starts to develop and quickly approaches the value predicted by linear spin wave theory.

While $\rho$ exhibits little dependence on system size (hence the $\rho$ plots of different $L$ lie on top of each other), $\zeta$ shows significant finite size effect. The latter is due to the entropic nature of the six-state clock anisotropy. In the $T\to 0$ limit, the probability distribution of the order parameter argument $\phi$ is given by $\rho(\phi)\propto\exp[\mathcal{N}s(\phi)]$, where $s(\phi)$ is the thermal entropy density due to spin wave fluctuations and $\mathcal{N}$ is the number of spins. The peaks of $\rho(\phi)$ are located at $\psi_{2}$ states (Fig. 5c). Crucially, for finite but large $\mathcal{N}$, the width of these peaks is inversely proportional to $\sqrt{\mathcal{N}}$. As $\mathcal{N}$ increases, more and more spin configurations fall on $\psi_{2}$ states, and hence $\zeta$ increases towards 1.

Appendix E Finite size effect of the order parameter dynamics

In our simulation, we use a system of $L\times{}L\times{}L$ unit cells with linear dimension $L=8$. Here, we assess the finite size effect of the order parameter dynamics by monitoring its dependence on $L$. Fig. 6 shows the argument and the modulus of the average Néel order parameter $\langle O\rangle$ as a function of time for $L=6,7,8$. The temperature $k_{B}T=10^{-3}J_{\pm}S^{2}$. The pulse polarization angle $\theta_{\mathrm{B}}=0.15\pi$.

For all three linear dimensions, $\mathrm{arg}\langle{}O\rangle$, the argument of the average Néel order parameter $\langle{}O\rangle$, exhibit qualitatively similar behavior within the time window of simulation. At early times, the oscillation amplitude of $\mathrm{arg}\langle{}O\rangle$ increases rather than decreases with $L$. Meanwhile, the center position of the oscillation, which indicates the position of a non-equilibrium potential minimum, changes with $L$ by a small amount. Likewise, the modulus of the average order parameter, $|\langle O\rangle|$, decays in time for all three $L$s. We also note that the decay rate of the modulus $|\langle O\rangle|$ shows non-monotonic dependence on $L$: the slope of the $L=7$ plot is greater than those of $L=6$ and $L=8$.

Without a finite-size scaling theory, we are unable to perform controlled extrapolation toward $L\to\infty$. Nevertheless, given that the order parameter dynamics shows only small dependence on $L$ within the simulated time window, the qualitative features of the non-equilibrium dynamics should be robust in the thermodynamic limit.