# Fluctuation control of non-thermal orbital order

Francesco Grandi  Department of Physics, University of Erlangen-Nürnberg,
91058 Erlangen, Germany  Martin Eckstein  Department of Physics, University of
Erlangen-Nürnberg, 91058 Erlangen, Germany

(August 29, 2024)

###### Abstract

Orbitally ordered states exhibit unique features which make them a promising
platform for exploring the ultrafast dynamics of long-range order in solids:
Their free energy typically has multiple discrete minima, and electric laser
fields or selectively excited phonons can exert effective forces that may be
used to steer the order parameter through these free energy landscapes.
Moreover, their free energy strongly depends on fluctuations, and in some
cases restoring forces close to a minimum are exclusively of entropic origin
(order-by-disorder mechanisms). This can open pathways to control the dynamics
of the order parameter via non-thermal fluctuations. In this work, we study
the laser-induced non-equilibrium dynamics in a $120^{\circ}$ compass model,
using time-dependent Ginzburg-Landau theory. We analyze protocols to switch
the order parameter between equivalent configurations, with a focus on the
interplay between the external force due to the driving field, and the non-
thermal entropic forces. In particular, we find that remanent non-thermal
fluctuations after some excitation can stabilize the high-symmetry phase even
when the homogeneous potential has retrieved its low-temperature form, which
facilitates laser-induced switching.

###### pacs:

##  I Introduction

The manipulation of matter with light Basov  _et al._ (2017); Giannetti  _et
al._ (2016) gives new perspectives on both fundamental theoretical questions
and future technological applications, such as magnetic memory devices that
operate at picosecond speed Kirilyuk  _et al._ (2010). Particularly appealing
is the ultrafast manipulation of multi-minima free-energy landscapes to
control competing phases or to reach hidden states along non-adiabatic
pathways Fausti  _et al._ (2011); Sun and Millis (2020); Stojchevska  _et al._
(2014); Ichikawa  _et al._ (2011); Först  _et al._ (2011); Li  _et al._
(2018). Such non-adiabatic transitions can be driven by generalized forces on
various order parameters, as obtained by the nonlinear effect of the electric
laser-field or selectively excited phonons on lattice distortions Först  _et
al._ (2011); Subedi  _et al._ (2014), superconducting pairing Murakami  _et
al._ (2017); Babadi  _et al._ (2017); Kennes  _et al._ (2017); Buzzi  _et al._
(2020), artificial magnetic fields Nova  _et al._ (2017), the Coulomb
interaction Singla  _et al._ (2015); Kaiser  _et al._ (2014); Grandi  _et al._
(2021), magnetic exchange interactions Wall  _et al._ (2009); Mentink  _et
al._ (2015), and more. A more indirect pathway to control the dynamics of an
order parameter relies on the excitation of fast degrees of freedom, such as
electronic variables, which then modify the phenomenological free energy
density $f(\phi_{0},T,...)$ of the order parameter $\phi_{0}$. Free energies
which depend on the fast degrees of freedom via few coarse-grained variables
only, such as an effective electronic temperature $T$ or a photo-doping
density, can indeed explain many observations regarding light-induced phase
transitions Yusupov  _et al._ (2010); Schäfer  _et al._ (2010); Huber  _et
al._ (2014); Beaud  _et al._ (2014).

On the other hand, assuming that $f(\phi_{0},T)$ is determined by a single
effective temperature even out of equilibrium would require all degrees of
freedom apart from $\phi_{0}$ to be in a quasi-equilibrium state, which is
clearly approximate. For example, recent observations show that the charge-
density wave (CDW) order in photo-excited TbTe3 apparently evolves in the low-
temperature potential, even if the electrons are still found to be above the
equilibrium $T_{c}$ Maklar  _et al._ (2020). In general, degrees of freedom
which are not equilibrated should be included explicitly in a phenomenological
description Zong  _et al._ (2019a, b); Kogar  _et al._ (2020). These non-
thermal fluctuations can renormalize the potential for $\phi_{0}$ in a
different way as compared to a temperature change. In particular, this
includes the long-wavelength spatial fluctuations of the order parameter
Dolgirev  _et al._ (2020), which, e.g., are found to delay the recovery of the
symmetry-broken state after an excitation for an $O(N)$-symmetric order
parameter.

A particularly important effect of fluctuations on the non-equilibrium
dynamics can be anticipated when fluctuations already have a qualitative
effect on the free energy landscape in equilibrium. This is the case for
systems which follow the so-called order-by-disorder mechanism. In this
scenario, the classical ground state potential has a continuous symmetry,
which is reduced to the lattice point group by fluctuations, i.e., by purely
entropic forces. The concept was introduced initially to describe the order
induced by quantum and thermal fluctuations in otherwise frustrated spin
models Henley (1989), and it is realized in particular in a class of spatially
anisotropic spin models known as compass models Nussinov and van den Brink
(2015). Such models are relevant for the description of orbital degrees of
freedom in transition metal compounds Kugel and Khomskii (1973) and for
applications in orbitronics Bernevig  _et al._ (2005) or topological quantum
computing Kitaev (2003). The time-dependent control of the orbital state in
solids might be achieved via an indirect coupling of the orbitals with the
electromagnetic field Eckstein  _et al._ (2017), or via their coupling to the
lattice Polli  _et al._ (2007); Miller  _et al._ (2015); Werner  _et al._
(2017); Beaulieu  _et al._ (2020). A measurement of the momentum-dependent
order parameter fluctuations could then be done using time-resolved inelastic
X-ray scattering (trRIXS) Chen  _et al._ (2019); Mitrano  _et al._ (2019);
Mitrano and Wang (2020); Parchenko  _et al._ (2020).

The focus of the present study is to investigate the effect of fluctuations on
the non-equilibrium order parameter dynamics in a paradigmatic model for
orbitally-ordered systems. We consider the $120^{\circ}$ compass model, an
orbital-only model representing Mott insulators with one electron (or one
hole) in an $\text{e}_{\text{g}}$ orbital doublet. In spite of the strong
directionality of the $\text{e}_{\text{g}}$ orbitals, the classical order
parameter potential at zero temperature has a continuous symmetry, which is
reduced to the cubic point group symmetry via the order-by-disorder mechanism.
We investigate the role of fluctuations on the ultrafast quench and recovery
after various excitation protocols. For a simple temperature quench of the
fast degrees of freedom, non-thermal fluctuations tend to stabilize the high
symmetry phase even when the temperature of the fast degrees of freedom goes
back to the original value, similar as for the $O(N)$-symmetric order
parameter Dolgirev  _et al._ (2020). As a second excitation protocol we
consider a time-dependent electric field, which can provide an effective
coherent force on the order parameter, and thus allows the system to switch
from one orbital state to another. We show how the spectrum of the
fluctuations evolves along the switching trajectory, emphasizing the interplay
between the external force due to the driving field, and the entropic forces.
Finally, combining the coherent driving with a temperature quench permits to
switch the orbital state much faster than at fixed temperature, similar to
heat-assisted switching in magnetic memory devices. In this protocol, the
transient trapping of the system in the high symmetry state due to non-thermal
fluctuations is essential to facilitate an efficient switching process.

The manuscript is structured as follows: In Sec. II, we introduce the model
and a time-dependent Ginzburg-Landau formalism implemented to take into
account the fluctuations of the order parameters. Sec. III presents the main
results of this work, particularly the effect of a temperature quench on the
dynamics of the order parameter (Sec. III.1), a coherent control of the order
parameter through the manipulation of the orbital exchange coupling (Sec.
III.2) and the combination of the two time-dependent perturbations (Sec.
III.3). In Sec. IV, we discuss possible experimental probes and establish a
connection between the correlation functions obtained within our formalism and
the structure factors measurable with trRIXS. Finally, Sec. V concludes the
work.

##  II Model and formalism

###  II.1 Ginzburg-Landau formalism for the $120^{\circ}$ compass model

In the following, we will analyze the $120^{\circ}$ compass model on a cubic
lattice as a specific representative for the general class of orbital-only
compass models Nussinov and van den Brink (2015). This model is obtained from
an underlying electronic Hamiltonian for a single electron or hole in the
$\text{e}_{\text{g}}$ orbital doublet, after freezing both charge degrees of
freedom (which is justified in the Mott regime) and the electron spin.
Freezing the spin can be achieved by strong magnetic fields, or alternatively,
the dynamics of the spin may be neglected to first order at temperatures where
magnetic order is absent while orbital order remains intact Nussinov  _et al._
(2004).

We define the orbital model in terms of a $\tfrac{1}{2}$-pseudospin
$\vec{\tau}=(\tau^{1},\tau^{2},\tau^{3})$ acting in the orbital
e${}_{\text{g}}$-doublet. Using a basis where the eigenstates of the
$\tau^{3}$ operators are the $|z^{2}-y^{2}\rangle$ and $|3x^{2}-r^{2}\rangle$
orbitals, $P^{x}=\frac{1}{2}(1-\tau^{3})$ is the projector on the
$|3x^{2}-r^{2}\rangle$ orbital. The projectors $P^{\gamma}$ on the orbital
$|3\gamma^{2}-r^{2}\rangle$ for $\gamma=x,y,z$ are then rotated with respect
to each other by $120^{\circ}$ around the $\tau^{2}$-direction in pseudospin
space, i.e., $P^{\gamma}=\tfrac{1}{2}(1-n^{\gamma}_{i}\tau^{i})$, where
summation over double indices is implied, and $n^{\gamma}_{i}$ are the
components of the unit vectors $\hat{n}^{x}=(0,0,1)$,
$\hat{n}^{y}=(\sqrt{3},0,-1)/2$, $\hat{n}^{z}=(-\sqrt{3},0,-1)/2$. The orbital
exchange Hamiltonian for a bond $(\mathbf{R},\mathbf{R}+\mathbf{e}_{\gamma})$
along a given direction $\gamma=x,y,z$ is
$-\bar{J}_{\gamma}P^{\gamma}_{\mathbf{R}}P^{\gamma}_{\mathbf{R}+\mathbf{e}_{\gamma}}$,
favoring both electrons in the orbital $|3\gamma^{2}-r^{2}\rangle$ with the
biggest overlap along the bond. The total Hamiltonian is therefore written as

$H=-\sum_{\mathbf{R}}\Big{[}\frac{1}{4}\sum_{\gamma}\bar{J}_{\gamma}n^{\gamma}_{i}n^{\gamma}_{j}\tau_{\mathbf{R}}^{i}\tau_{\mathbf{R}+\mathbf{e}_{\gamma}}^{j}-\sum_{\gamma}\bar{J}_{\gamma}n^{\gamma}_{i}\tau_{\mathbf{R}}^{i}\Big{]},$
(1)

up to irrelevant constants. In equilibrium, the exchange constants
$\bar{J}_{\gamma}$ for $\gamma=x,y,z$ are equal, but they can become different
when external fields are applied (see below). In particular, this implies that
the second term in the Hamiltonian, which corresponds to a pseudo-spin
magnetic field, is absent in equilibrium, because $\sum_{\gamma}n_{\gamma}=0$.

The classical ground state of model Eq. (1) is ferro-orbitally ordered, and
fully degenerate as function of the angle of the ordered moment in the
$\left(\tau^{3},\tau^{1}\right)$ plane. However, as soon as thermal or quantum
fluctuations are taken into account, the degeneracy is broken, leading to the
stabilization of six degenerate minima which correspond to equivalent ferro-
orbital order polarizations that are mapped onto each other by the cubic point
group operations van den Brink  _et al._ (1999); van den Brink (2004); Biskup
_et al._ (2005). This order-by-disorder mechanism can be understood already at
the level of spin-wave theory, because the excitation spectrum, and thus the
entropic contribution to the free energy at nonzero temperature, depends on
the direction of the ordered orbital moment.

To describe the physics of model (1) within the time-dependent Ginzburg-Landau
formalism, we write a free-energy functional which is consistent with the
continuum limit of the Hamiltonian (1),

$\displaystyle\mathcal{F}[\phi]$ $\displaystyle=\int
d^{3}\mathbf{x}\Big{[}\,\frac{r(T)}{2}\sum_{\gamma}J_{\gamma}n^{\gamma}_{i}n^{\gamma}_{j}\phi_{j}\phi_{i}-\sum_{\gamma}J_{\gamma}n^{\gamma}_{i}\phi_{i}$
$\displaystyle+u(\phi_{i}\phi_{i})^{2}+\frac{K}{2}\sum_{\gamma}J_{\gamma}n^{\gamma}_{i}n^{\gamma}_{j}(\partial_{\gamma}\phi_{j})(\partial_{\gamma}\phi_{i})\Big{]}.$
(2)

Here $(\phi_{1}\left(\mathbf{x},t\right),\phi_{2}\left(\mathbf{x},t\right))$
is a two-dimensional order parameter that corresponds to the coarse-grained
pseudospin components ($\tau^{3},\tau^{1})$. From now on, we restrict the
orbital indexes to the values $1,2$, imply summation over double orbital
indices, and redefine the unitary vectors $\hat{n}^{x}=(1,0)$,
$\hat{n}^{y}=(-1,\sqrt{3})/2$ and $\hat{n}^{z}=(-1,-\sqrt{3})/2$ by omitting
the irrelevant $\tau^{2}$ component. The coefficient of the quadratic
contribution in Eq. (II.1) is proportional to a temperature-dependent factor
$r(T)=-2(1-\frac{T}{T_{c}})$ which changes sign at the critical temperature
$T_{c}$. The quartic term $u$ implements a soft-spin approximation and
prevents the order parameter to grow arbitrarily large. The term proportional
to $K$, which has the form of a gradient approximation to the exchange
interaction in Eq. (1), is the free energy contribution from spatial order
parameter fluctuations. The couplings $J_{\gamma}$ in the free energy are
related to $\bar{J}_{\gamma}$ in Eq. (1) through some coarse-graining, and
must thus be understood as phenomenological parameters. However, one can
assume that they can be modified by external fields in the same way as the
bare parameters $\bar{J}_{\gamma}$. In particular, $J_{x}=J_{y}=J_{z}\equiv J$
in equilibrium, which implies that the free energy density arising from the
first three terms in Eq. (II.1) becomes
$\frac{3}{4}r(T)J\phi_{i}\phi_{i}+u(\phi_{i}\phi_{i})^{2}$, corresponding to
the standard $O(2)$-symmetric $\phi^{4}$ theory Dmitrašinović  _et al._
(1995). In contrast, the term proportional to $K$ is only invariant under the
cubic point group, implementing again the order-by-disorder scenario.

###  II.2 Time-evolution using model A dynamics

We compute the time evolution assuming an over-damped dissipative dynamics
(model A dynamics), which is applicable for anisotropic magnets and systems
without conserved fields Hohenberg and Halperin (1977); Stahl and Eckstein
(2021). The dissipative equations of motion read

$\dot{\phi}_{i}\left(\mathbf{x},t\right)=-\Gamma\frac{\delta\mathcal{F}}{\delta\phi_{i}\left(\mathbf{x},t\right)}+\eta_{i}\left(\mathbf{x},t\right)\;,$
(3)

where $\Gamma>0$ is the relaxation rate, and
$\eta_{i}\left(\mathbf{x},t\right)$ is a Gaussian white noise term that takes
into account the effect of the fast degrees of freedom, characterized by the
moments

$\begin{split}&\langle\eta_{i}\left(\mathbf{x},t\right)\rangle=0\;,\\\
&\langle\eta_{i}\left(\mathbf{x},t\right)\eta_{j}\left(\mathbf{x^{\prime}},t^{\prime}\right)\rangle=2T\Gamma\delta_{i,j}\delta\left(\mathbf{x}-\mathbf{x^{\prime}}\right)\delta\left(t-t^{\prime}\right)\;.\end{split}$
(4)

The parameter $T$ in Eq. (4) is the temperature of the fast degrees of
freedom, which are integrated out in the coarse-graining under the assumption
that they instantly thermalize on the timescale of the order parameter
dynamics. The noise temperature $T$ will therefore be assumed to be equal to
the temperature in the coefficient $r(T)$, even in an excitation protocol
where $T$ becomes time-dependent (see Sec. III).

With the Fourier transformation of the order parameter we define the
correlation functions
$D_{ij}\left(\mathbf{k},t\right)=\langle\phi_{i}\left(\mathbf{k},t\right)\phi_{j}\left(-\mathbf{k},t\right)\rangle$
in reciprocal space, and the total number of excitations in a particular
channel,

$n_{ij}\left(t\right)=\frac{1}{8\pi^{3}}\int_{-\Lambda}^{\Lambda}dk_{1}\int_{-\Lambda}^{\Lambda}dk_{2}\int_{-\Lambda}^{\Lambda}dk_{3}\
D_{ij}\left(\mathbf{k},t\right).$ (5)

Here $\Lambda$ is an ultraviolet cutoff. Note that the correlation matrix is
symmetric, $n_{ij}=n_{ji}$. The dynamical model given by Eqs. (3) and (4) can
then be recast as a Fokker-Planck-Smoluchowski equation for the probability
distribution of the fields $\phi_{i}\left(\mathbf{x},t\right)$. Under the
assumption that the probability distribution is Gaussian, we obtain a closed
set of equations for the correlation function and the statistically averaged
fields
$\bar{\phi}_{i}\left(t\right)=\langle\phi_{i}\left(\mathbf{x},t\right)\rangle$,

$\displaystyle\dot{\bar{\phi}}_{i}(t)$
$\displaystyle=-\Gamma\Big{[}A^{\text{eff}}_{ij}\bar{\phi}_{j}-\sum_{\gamma}J_{\gamma}n^{\gamma}_{i}\Big{]},$
(6) $\displaystyle\dot{D}_{ij}(\mathbf{k},t)$
$\displaystyle=\Gamma\Big{[}2\delta_{ij}T-M_{il}D_{lj}(\mathbf{k},t)-D_{il}(\mathbf{k},t)M_{lj}\Big{]},$
(7)

where

$\displaystyle A^{\text{eff}}_{ij}$
$\displaystyle=A^{0}_{ij}+4u\big{[}\delta_{ij}\sum_{l}(\bar{\phi}_{l}^{2}+n_{ll})+2n_{ij}\big{]},$
(8) $\displaystyle M_{ij}$
$\displaystyle=A^{\text{eff}}_{ij}+\kappa_{ij}+8u\bar{\phi}_{i}\bar{\phi}_{j}\;,$
(9)

with the bare coupling matrix

$A^{0}_{ij}(T)=r(T)\sum_{\gamma}J_{\gamma}n^{\gamma}_{i}n^{\gamma}_{j}\;,$
(10)

and the distortion contribution

$\kappa_{ij}=K\sum_{\gamma}J_{\gamma}n^{\gamma}_{i}n^{\gamma}_{j}k_{\gamma}^{2}\;.$
(11)

The second term in the first equation acts like a bare external force on the
order parameter,

$\displaystyle f^{0}_{i}=\sum_{\gamma}J_{\gamma}n^{\gamma}_{i}.$ (12)

The equations are solved numerically on an equidistant grid of k-points;
$N_{k}\sim\left(300\right)^{3}$ to $\left(400\right)^{3}$ is sufficient to
obtain converged results.

To analyze the results below, we will parametrrize the order parameter
($\bar{\phi}_{1},\bar{\phi}_{2})$ in polar coordinates, with amplitude and
phase defined as

$\begin{pmatrix}\bar{\phi}_{1}\\\
\bar{\phi}_{2}\end{pmatrix}=R\begin{pmatrix}\sin(\varphi)\\\
-\cos(\varphi)\end{pmatrix}\;.$ (13)

The zero-temperature stationary solution of Eqs. (6) and (7) implies
$n_{ij}=0$, and reflects the $\text{O}(2)$ symmetry of the problem. For
nonzero temperatures below the critical value, there are six equivalent
stationary solutions characterized by the same amplitude $R$ but different
phase $\varphi=\pi/6,\pi/2,5\pi/6,7\pi/6,3\pi/2$ and $11\pi/6$, each separated
by an angle $\pi/3$. More details about the equilibrium solution can be found
in App. A.

In the following, we take the parameters $J_{x}=J_{y}=J_{z}=5$, $u=1$, $K=1$,
$T_{c}=0.5$, $\Lambda=\pi$ and $\Gamma=0.5$. At zero temperature, the free
energy landscape (as well as $\Gamma$ and $\Lambda$) is the same as for the
totally $O(N)$-symmetric model defined in Ref. Dolgirev  _et al._ (2020) when
$N=2$. Unless otherwise specified, the equilibrium temperature of the system
is set to $T=0.125$ and, without any loss of generality, we initially place
the state of the system at $\varphi=\pi/6$.

###  II.3 Effective potential for the homogeneous order parameter

Looking at Eq. (6), one can formally write the dynamics for the homogeneous
order parameter $\bar{\phi}$ in terms of an effective potential
$\bar{F}(\bar{\phi})$, such that

$\displaystyle\dot{\bar{\phi}}_{i}$
$\displaystyle=-\Gamma\frac{\partial\bar{F}(\bar{\phi})}{\partial\bar{\phi}_{i}}.$
(14)

The corresponding free energy $\bar{F}$,

$\displaystyle\bar{F}(\bar{\phi})$
$\displaystyle=\frac{1}{2}A^{0}_{ij}(T)\bar{\phi}_{i}\bar{\phi}_{j}+u(\bar{\phi}_{i}\bar{\phi}_{i})^{2}-f_{i}^{0}\bar{\phi}_{i}+F_{\text{fl}}(\bar{\phi}),$
(15)

is given by the homogenous part of Eq. (II.1) and a quadratic fluctuation
potential
$F_{\text{fl}}(\bar{\phi})=\tfrac{1}{2}A_{ij}^{\text{fl}}\bar{\phi}_{i}\bar{\phi}_{j}$,
with the fluctuation matrix

$\displaystyle
A_{ij}^{\text{fl}}=4u\big{[}\delta_{ij}\sum_{l}n_{ll}+2n_{ij}\big{]}.$ (16)

During the dynamics, the fluctuations change in a nontrivial way, and the
matrix $A^{\text{fl}}$ will provide a convenient way to quantify their effect
on the dynamics of the order parameter.

##  III Results

Figure 1:  (a) Time-dependent profile of the temperature quench (17), for
$\tau_{QP}=0.3$, $T_{\text{s}}=0.125$ and $T_{\text{f}}=1.428$. (b) Time
dependence of the exchange parameter $J_{x}$ used in the switching protocol
[Eq. (20) with $\delta J=0.3$, $J_{x}^{0}=5$, $\bar{\Omega}=10$ and $t_{p}\sim
30$] together with its moving average $\langle J_{x}\rangle$.

In the following, we analyze three different time-dependent excitation
protocols: (i) A short temperature quench with subsequent recovery,
implemented by a time-dependent temperature $T(t)$ in Eqs. (4), (7), and (10)
(see Fig. 1(a)), (ii), a time-dependent modulation of the the renormalized
exchange interaction $J_{x}\left(t\right)$, as shown in Fig. 1(b), and (iii),
a combination of (i) and (ii). The excitation protocol (i) is used frequently
in the literature to describe the laser-induced heating of the fast degrees of
freedom, such as electrons, which couple directly to the light pulse and are
assumed to thermalize instantly on the timescale of the order-parameter
dynamics. The exponential recovery of the initial temperature after the pulse
reflects the energy dissipation to degrees of freedom which do not strongly
affect the order parameter (such as acoustic phonons). As mentioned, $T(t)$
enters both the free energy potential (II.1) and the noise (4).

The excitation protocol (ii) finds its rationale in the possibility of
controlling the time-dependent profiles of
$\bar{J}_{x}\left(t\right),\bar{J}_{y}\left(t\right)$ and
$\bar{J}_{z}\left(t\right)$ separately with suitable electric fields directed
along given lattice directions. Orbital exchange interactions can, e.g., be
changed directly by manipulating the exchange path with the electric field of
a laser. Explicit expressions for the exchange couplings in the field-driven
two-orbital Hubbard model have been obtained from a time-dependent Schrieffer-
Wolff transformation Eckstein  _et al._ (2017), indicating that few-percent
changes of $J$ may be possible in experiment. One can also envision other
symmetry-allowed processes, such as driving infrared-active phonons which are
nonlinearly coupled to the orbital pseudospin. In this study, we assume that
the effective parameters $J_{x}\left(t\right),J_{y}\left(t\right)$ and
$J_{z}\left(t\right)$ in the coarse-grained free energy are controlled in the
same way as the microscopic parameters. This enables a coherent control of the
order parameter, with a switching between different free energy minima.

Finally, protocol (iii) is motivated by the aim to explore how the
fluctuation-induced change of the free energy cooperates or competes with the
explicit control by the modulation of the exchange couplings. In practice, any
implementation of the coherent control protocol (ii) will go together with a
transient heating of the fast degrees of freedom.

###  III.1 Temperature quench

Figure 2: Solid lines show the time-dependent change of the amplitude of the
order parameter (normalized by the equilibrium value $R_{0}$) during
temperature quenches to different maximum temperatures $T_{\text{f}}$. For all
simulations, the initial temperature is $T_{\text{s}}=0.125$, and
$\tau_{QP}=0.3$. Dashed lines correspond to the dynamics of the order
parameter computed with the free energy (18) without including the
fluctuations (see text). The shaded region indicates times for which $T>T_{c}$
for the quench at $T_{\text{f}}=2.500$.

In this section, we focus on the effects of a temperature quench on the
dynamics of the system. We assume the time-dependent temperature

$T\left(t\right)=T_{\text{s}}+\left(T_{\text{f}}-T_{\text{s}}\right)\theta\left(t\right)e^{-t/\tau_{QP}}\;,$
(17)

where $\theta\left(t\right)$ is the Heaviside step function, $T_{\text{s}}$
the pre-pulse temperature, $T_{\text{f}}$ the maximum value of the temperature
reached at $t=0^{+}$, and $\tau_{QP}=0.3$ the relaxation rate to
$T_{\text{s}}$ (see Fig. 1(a)). The dynamics following such quenches is shown
in Fig. 2. The amplitude $R(t)$ of the order parameter is suppressed by the
temperature increase and subsequently recovers to the initial value. Upon
increasing $T_{\text{f}}$, the order parameter is suppressed to a value closer
to zero, and the recovery time increases.

The phase $\varphi=\pi/6$ remains constant after the temperature quench. To
characterize the dynamics, we can therefore analyze the effective potential
(15) at fixed angle $\varphi$ as a function of $R$. With
$J_{x}=J_{y}=J_{z}\equiv J$,

$\displaystyle\bar{F}(R)$
$\displaystyle=\tfrac{1}{2}[A^{0}_{R}(t)+A^{\text{fl}}_{R}(t)]R^{2}+uR^{4},$
(18)

where $A^{0}_{R}(t)=\frac{3}{2}r(T(t))J$, and the fluctuation component
$A^{\text{fl}}_{R}(t)$ is given by the projection of the matrix
$A_{ij}^{\text{fl}}$ in Eq. (16) along the unit vector
$\hat{e}^{R}=(\sin(\pi/6),-\cos(\pi/6))$ in radial direction in orbital space.
If one would consider only the dynamics generated by the bare potential
($\dot{R}=-\Gamma\partial_{R}\bar{F}(R)$, with $A^{\text{fl}}=0$), one obtains
a somewhat similar result as for the full solution (see dashed lines in Fig.
2). In this case, the delayed recovery after strong quenches has a trivial
interpretation: The quench brings the order parameter to a small value close
to $R=0$. After some time $t_{\text{min}}$ when the initial temperature is
recovered (which is very early on the scale of the plot, see shaded region in
Fig. 2), the order parameter undergoes an exponential growth phase $R(t)\sim
R_{\text{min}}e^{A^{0}(t-t_{\text{min}})}$ dominated by the quadratic term in
$\bar{F}$, starting from the small value $R_{\text{min}}$ at time
$t=t_{\text{min}}$. The stronger the quench, the smaller $R_{\text{min}}$, and
the longer it takes for the order parameter to leave the unstable region.
However, it is clear from the comparison of the dashed and full lines in Fig.
2, that this argument would largely underestimate the delay in the recovery of
the order parameter. The latter is therefore described only if the
fluctuations are taken into account. This delayed recovery is qualitatively
the same as given for the $O(N)$-symmetric model Dolgirev  _et al._ (2020). To
further analyze the behavior, we plot in Fig. 3(a) both the bare component
$A^{0}_{R}$ and the fluctuation-renormalized coefficient
$A^{0}_{R}+A_{R}^{\text{fl}}$ [c.f. Eq. (18)] of the quadratic term in
$\bar{F}(R)$. The comparison shows that while the bare potential would be
already unstable for most of the time ($A^{0}_{R}<0$), the fluctuations tend
to stabilize the high-symmetry state.

Figure 3:  (a) Analysis of the free energy (18) close to $R=0$: The bare
coefficient $A^{0}_{R}$ of the quadratic term and the fluctuation-renormalized
coefficient $A^{0}_{R}+A^{\text{fl}}_{R}$ are plotted as function of time for
the quench with $T_{\text{f}}=2.500$ shown in Fig. 2. The dashed line
indicates the zero of the $A_{R}$ parameter. (b) Excitation density in polar
coordinates for the same quench. The mixed component $n_{\varphi R}$ remains
equal to zero during the whole dynamics.

Finally, it is interesting to analyze the fluctuations themselves. For this
purpose, we chose a representation using polar coordinates in orbital space,

$\begin{split}&n_{\varphi\varphi}=\sin^{2}\left(\varphi\right)n_{22}+\cos^{2}\left(\varphi\right)n_{11}+\sin\left(2\varphi\right)n_{12}\;,\\\
&n_{RR}=\sin^{2}\left(\varphi\right)n_{11}+\cos^{2}\left(\varphi\right)n_{22}-\sin\left(2\varphi\right)n_{12}\;,\\\
&n_{\varphi
R}=\sin\left(\varphi\right)\cos\left(\varphi\right)\left(n_{11}-n_{22}\right)-\cos\left(2\varphi\right)n_{12}\;,\end{split}$
(19)

with $\varphi=\pi/6$. The components are obtained by projecting the
fluctuation matrix $n_{ij}$ along the given directions in orbital space, e.g.,
$n_{RR}=n_{ij}e^{R}_{i}e^{R}_{j}$, with the unit vector
$e^{R}=(\sin(\varphi),-\cos(\varphi))$ along the radial direction. In Fig.
3(b), a separation of the time-scales in the recovery of the radial and
angular fluctuations ($n_{RR}$ and $n_{\varphi\varphi}$, respectively) becomes
apparent. This observation simply reflects the fact that the free energy
landscape is more shallow in the angular than in radial direction. A similar
effect has been observed for the $O(N)$-symmetric model Dolgirev  _et al._
(2020). In the present case, the model does not have a continuous symmetry at
nonzero temperature, and the angular fluctuations are therefore stronger
confined. A measurement of the fluctuation dynamics in the two directions (see
also Sec. IV) could allow to map out this nontrivial potential landscape
arising from the order-by-disorder mechanism.

###  III.2 Time-dependent control of the exchange parameters

Figure 4:  (a) Time-dependence of the angle $\varphi$ when the system evolves
with the time-dependent exchange interactions shown in Fig. 1(b), for several
values of the variation $\delta J$. The pulse durations are $t_{p}\sim
65,40,30$ for $\delta J=0.1,0.2,0.3$, respectively. (b) Same as panel (a), but
for fixed pulse parameters $\delta J=0.1$ and $t_{p}\sim 65$, and different
initial temperatures.

In this section, we describe the time evolution of the system during an
orbital switching process, where the state moves from the free energy minimum
at $\varphi=\pi/6$ to the one at $\varphi=\pi/2$. To achieve the switching, we
modulate the value of the exchange coupling $J_{x}$ as shown in Fig. 1(b). The
precise time profile is

$J_{x}(t)=J_{x}^{0}+\delta J\cos(\bar{\Omega}t)^{2}s(t),$ (20)

where $s(t)$ is a smooth envelope of duration $t_{p}$. The oscillatory time-
dependence, with a large frequency $\bar{\Omega}$ is motivated by the shape
one would obtain when the exchange couplings are controlled with oscillatory
electric fields Eckstein  _et al._ (2017). Clearly, the fast oscillations are
almost irrelevant for the slow dynamics of the order parameter, and one could
equally well replace $\cos(\bar{\Omega}t)^{2}\to\tfrac{1}{2}$. The pulse is
therefore determined by the equilibrium coupling $J_{x}^{0}$, the duration
$t_{p}$, and the amplitude $\delta J$, which will generally taken to be of the
order of few percent of $J$.

The first effect of this protocol is to make the pseudo-spin magnetic fields
$f_{i}^{0}$ nonzero [c.f. Eq. (12)]. In addition, one has a direct
modification of the quadratic terms $A_{ij}^{0}$ in the potential, and the
indirect effect through the fluctuations. As shown in Fig. 4(a), the given
protocol can lead to a switch from one minimum to the other. As expected, with
larger $\delta J$, the final minimum is reached faster. We now fix the pulse
parameters $\delta J=0.1$ and $t_{p}\sim 65$ and compare the dynamics for
different temperatures (Fig. 4(b)). The resulting switching time is non-
monotonous as a function of $T$. As the temperature is increased, the
switching first becomes slower (compare the curves for $T=0.125$ and
$T=0.400$), and then again faster, as the critical temperature is approached.
As we will explain in the following, this non-monotonous behavior again
reflects the order-by-disorder mechanism.

The trajectories in the order parameter plane, corresponding to the switching
protocols described in Fig. 4(b), are shown in Fig. 5(a). The three curves, at
different amplitude of the order parameter, correspond to the three different
temperatures. The non-monotonous behavior of the switching speed can now be
understood by analyzing the different contributions to the effective force on
the order parameter. For this purpose, we project the external force
$f_{i}^{0}$ (Eq. (12)) at each point $\varphi$ on the trajectory along the
angular direction $e^{\varphi}=(\cos(\varphi),\sin(\varphi))$ in orbital
pseudo-spin space,

$f^{0}_{\varphi}=f^{0}\cdot
e^{\varphi}=f_{1}^{0}\cos(\varphi)+f_{2}^{0}\sin(\varphi)\;,$ (21)

and compare it with the force contribution from the fluctuating part of the
potential (15),

$f_{i}^{\text{fl}}=-\partial
F^{\text{fl}}/\partial\bar{\phi}_{i}=-A^{\text{fl}}_{ij}\bar{\phi}_{j}\;,$
(22)

also projected along the angular direction
$f^{\text{fl}}_{\varphi}=f^{\text{fl}}\cdot e^{\varphi}$. The result is shown
in Fig. 5(b). Because $f^{0}$ depends on the rapidly oscillating parameter
$J_{x}$, see Eq. (21), we show the moving average $\bar{f}^{0}_{\varphi}$ over
an oscillation period. In this representation, a positive force contributes to
the drag in the angular direction from $\pi/6$ to $\pi/2$, while a negative
force opposes the motion. The bare $\bar{f}^{0}_{\varphi}$, which directly
depends on the exchange constants, is almost unchanged with temperature. For
phases close to $\pi/6$, the fluctuating force is negative, indicating that
the fluctuations would pull the system back to the original minimum, while for
$\varphi$ larger than $\pi/3\sim 1.0472$, the fluctuation force pulls the
system towards the new minimum, and the system would evolve towards
$\varphi=\pi/2$ even without the external force $f^{0}_{\varphi}$. Apparently,
it is this fluctuation force which reflects the non-monotonous dependence of
the switching on $T$. This would be consistent with the order-by-disorder
mechanism: For small $T$, the increase of the temperature increases the
fluctuation-induced forces, while closer to the phase transition, where the
amplitude $R$ of the order parameter is small, the momentum fluctuations
$n_{ij}$ become isotropic. However, this argument holds formally only for
equilibrium, while for the non-equilibrium protocol one should evaluate the
fluctuations along the whole trajectory.

Figure 5:  (a) Trajectories of the order parameter for the manipulation of
$J_{x}$ at three different initial temperatures [fixed pulse parameters
$\delta J=0.1$ and $t_{p}\sim 65$ as in Fig. 4(b)]. The arrows indicate the
direction of motion. The shaded grey areas along the curve at $T=0.125$ show a
polar representation $n(\varphi^{\prime},t)$ of the non-equilibrium
fluctuation matrix around different points $\varphi$ along the trajectory [see
Eq. (23)]. The dashed green lines show the corresponding equilibrium
fluctuations distribution for the same points in the order parameter plane.
(b) Projections of the moving average of the external
($\bar{f}^{0}_{\varphi}$, solid lines) and fluctuating
($f^{\text{fl}}_{\varphi}$, dashed lines) contributions of the forces along
the angular direction as function of the angle $\varphi$ for the three
trajectories considered in panel (a).

To characterize the dynamics of the fluctuations during this time-dependent
process, we represent the fluctuation matrix $n_{ij}(t)$ in a polar plot,
projecting it along the direction
$e^{R}(\varphi^{\prime})=(\sin(\varphi^{\prime}),-\cos(\varphi^{\prime}))$ in
orbital pseudospin space

$n\left(t,\varphi^{\prime}\right)=e_{i}^{R}(\varphi^{\prime})e_{j}^{R}(\varphi^{\prime})n_{ij}\left(t\right)\;.$
(23)

The function $n(t,\varphi^{\prime})$ is shown for $T=0.125$ as polar plot
around several points on the trajectory (grey shaded areas). The plot shows a
characteristic double lobe structure, indicating a minimum of the fluctuations
in the longitudinal direction and a maximum in the transverse direction. In
agreement with the entropic origin of the free energy minima, the strongest
confinement of the fluctuations is observed around the thermodynamically most
unstable points (middle of the trajectory in Fig. 5(a)).

Figure 6:  Trajectory at the initial temperature $T=0.125$ for the driving of
the effective $J_{x}$ for pulse parameters $\delta J=1.0$ and $t_{p}\sim 13$.
As described in Fig. 5(a), the grey shaded areas and the green dashed lines
compare the angular distribution of the fluctuations in nonequilibrium and the
corresponding equilibrium distribution at different points of the trajectory.

By comparison one can see that the nonequilibrium angular distribution of the
fluctuations at different stages of the dynamics (shaded grey areas in Fig.
5(a)) and the corresponding equilibrium distribution for fixed values of the
angle and the radius of the order parameter (dashed green lines in Fig. 5(a))
are almost identical. This shows that the switching process is quasi-
adiabatic, as far as the integrated fluctuations are concerned. The quasi-
adiabatic switching can be contrasted with a process where the external
electric field is substantially larger. The trajectory for a switching with
$\delta J=1.0$ (corresponding to a $20$ percent change of the exchange
couplings) is shown in Fig. 6. The correspondingly faster switching indeed
leads to a mismatch between the time-dependent fluctuations and the
equilibrated fluctuations at each angle $\varphi$. In general, the time-
dependent fluctuations lag behind the corresponding equilibrium fluctuations.

Figure 7:  Spherical representation of the integral of the trace of the
correlation matrix Eq. (24) at several points of the trajectory shown in Fig.
5a for $T=0.125$ and $\delta J=0.1$. The five plots correspond to the points
for which the momentum integrated fluctuations are shown in Fig. 5a, at (a)
$\varphi=\pi/6$, (b) $\varphi=\pi/4$, (c) $\varphi=\pi/3$, (d)
$\varphi=5\pi/12$ and (e) $\varphi\sim\pi/2$.

As already the momentum averaged fluctuations show a non-trivial dynamics, it
would be interesting to map out the momentum-resolved fluctuations. These are
in principle accessible using scattering techniques (Sec. IV), and their
measurement would provide the most detailed characterization of the free
energy landscape along the switching path. In our model, the momentum-resolved
fluctuations are given by the momentum distribution of the correlation matrix
$D(\mathbf{k},t)$, whose entries are orbitally resolved. To simplify the
representation of this quantity, we analyze its trace
$D_{\text{Tr}}(\mathbf{k},t)=\sum_{i}D_{ii}(\mathbf{k},t)$. Moreover, since we
are mainly interested in the directionality of the correlation function in $k$
space, we show its integral over the modulus of the wave vector,

$D_{\text{Tr}}\left(\zeta,\chi,t\right)=\int_{0}^{\Lambda}dk\
k^{2}D_{\text{Tr}}\left(\mathbf{k},t\right)\;,$ (24)

having defined the vector $\mathbf{k}=\left(k_{x},k_{y},k_{z}\right)$ in
spherical coordinates $k_{x}=k\sin\left(\zeta\right)\cos\left(\chi\right)$,
$k_{y}=k\sin\left(\zeta\right)\sin\left(\chi\right)$ and
$k_{z}=k\cos\left(\zeta\right)$. At the initial point of the dynamics, Eq.
(24) shows peaks structure along the $\pm k_{x}$, see Fig. 7(a). During the
drag of the system from one minimum to the other, the k-resolved fluctuations
Eq. (24) continuously transfer weight from the original maxima along $\pm
k_{x}$ direction to the $\pm k_{y}$ direction, where the maxima of the final
state are observed, see Fig. 7(e). At intermediate steps, the distribution of
the fluctuations lies between these two limits, showing four well resolved
maxima in the k distribution when the system is at the maximum of the free
energy landscape, see Fig. 7(c). In Sec. IV, we describe a possible
experimental probe to measure the $k$-resolved dynamics of the fluctuations.

###  III.3 Heat-assisted switching

Figure 8:  (a) Time-dependent change of the amplitude of the order parameter
(normalized by the equilibrium value $R_{0}$) during a simultaneous modulation
of $J_{x}$ [pulse parameters $\delta J=0.1$ and $t_{p}\sim 13$] and a
temperature quench, for several maximum temperatures $T_{\text{f}}$. (b) Phase
of the order parameter $\varphi$ during the same time-dependent protocols. (c)
Corresponding trajectories in orbital space. Arrows in panel (c) show the
direction in which the state of the system moves along the trajectory.

In this section, we analyze the simultaneous action of the two time-dependent
perturbations considered in Sec. III.1 and Sec. III.2. Because we have seen
that in equilibrium a sufficiently high initial temperature simplifies the
switching protocol, one can ask whether the switching can be facilitated by
quenching the temperature in addition to the modification of the exchange
couplings. Our data confirm we can answer this question affirmatively, as we
are going to show.

In Fig. 8, we analyze the response of the system to temperature quenches of
increasing strength when we simultaneously apply a time-dependent $J_{x}$ with
fixed parameters $\delta J=0.1$ and $t_{p}\sim 13$. As apparent from Fig.
4(a), such a short pulse without the temperature quench would not be able to
induce the switch to the new orbital state. We generally notice that the
presence of the time-dependent $J_{x}$ makes the high-symmetry state at $R\sim
0$ less stable, compare the time scale for the recovery in Fig. 8(a) and Fig.
2. For small quenches in temperature, the dynamics of the phase of the order
parameter is not substantially perturbed to the case without the quench, see
Fig. 8(b) for $T_{\text{f}}=2.000$. In this case, $\varphi$ first increases
due to the action of $J_{x}(t)$, but it slowly recovers the original value
$\varphi=\pi/6$ once the pulse is over. If the quench is sufficiently strong,
the stability of the $R\sim 0$ state lasts enough for the time dependent
$J_{x}$ to make the switch, as observed in Fig. 8(b) for $T_{\text{f}}=2.857$.
Trajectories for different quenches are shown in Fig. 8(c).

This result is akin to the mechanism observed in heat-assisted recording
Kryder  _et al._ (2008), where the heating of magnetic memory devices reduces
the size of the magnetic field which is needed to switch the magnetization
state of the system against some anisotropy. However, there are some important
differences (apart from the fact that orbital order responds to different
forces). Firstly, the order-by-disorder mechanism implies that there is a
threshold for the temperature quench below which it would no longer facilitate
the switching, because for low temperatures the restoring forces would
increase with increasing fluctuations instead of decreasing. Moreover, we
emphasize that the transient stabilization of the high-symmetry state ($R\sim
0$) by the fluctuations also simplifies the switching because this prolongs
the time in which the order parameter can be manipulated with relatively weak
fields.

##  IV Measurement of the fluctuations

In this section, we briefly discuss the measurement of the $k$-resolved
fluctuation spectrum based on X-ray scattering. RIXS can give insights into
orbital correlations in the ground state of correlated materials Forte  _et
al._ (2008); Marra  _et al._ (2012); Ament  _et al._ (2011). Time-resolved
RIXS, based on free-electron laser sources, is therefore a potentially very
versatile probe for complex solids out of equilibrium Mitrano and Wang (2020).
In RIXS, an incoming X-ray photon with energy $\omega$ and momentum
$\mathbf{k}$ excites a core hole that subsequently decays through the emission
of an outgoing photon with energy $\omega^{\prime}$ and momentum
$\mathbf{k}^{\prime}$. The excess energy $\Omega=\omega-\omega^{\prime}$ and
momentum $\mathbf{q}=\mathbf{k}-\mathbf{k}^{\prime}$ are transferred to the
system, which allows to explore the energy momentum dispersion of the low
energy excitations. The polarization $\bm{\epsilon}$ and
$\bm{\epsilon}^{\prime}$ of the ingoing and outgoing photons can be used to
gain orbital selectivity. In general, under the assumption of a localized core
hole with a short lifetime $1/\Gamma$, the scattering can be expressed in
terms of a structure factor for transition operators $A$ and $B$ acting in the
low-energy manifold of the solid (valence excitations without core-hole),

$S_{A,B}\left(\mathbf{q},\Omega\right)\\!=\\!\sum_{f}\langle
g|A_{\mathbf{q}}^{\dagger}|f\rangle\langle
f|B_{\mathbf{q}}|g\rangle\delta(\Omega+E_{g}-E_{f}).$ (25)

Here $|g\rangle$, with energy $E_{g}$, is the initial ground state of the
system, and $|f\rangle$ an excited state in the low-energy valence manifold.
(The ground state average is replaced by a statistical average for finite
temperature.) Within this approximation, the RIXS intensity for scattering
between different orbital states in the $e_{g}$ manifold has been expressed in
terms of the orbital structure factor Marra (2016),

$\displaystyle
I_{\text{RIXS}}\propto\left|W^{A}_{\bm{\epsilon},\bm{\epsilon}^{\prime}}\right|^{2}S_{\tau^{3},\tau^{3}}\left(\mathbf{q},\Omega\right)+\left|W^{O}_{\bm{\epsilon},\bm{\epsilon}^{\prime}}\right|^{2}S_{\tau^{1},\tau^{1}}\left(\mathbf{q},\Omega\right)$
$\displaystyle+\left[(W^{A}_{\bm{\epsilon},\bm{\epsilon}^{\prime}})^{*}W^{O}_{\bm{\epsilon},\bm{\epsilon}^{\prime}}+W^{A}_{\bm{\epsilon},\bm{\epsilon}^{\prime}}(W^{O}_{\bm{\epsilon},\bm{\epsilon}^{\prime}})^{*}\right]S_{\tau^{1},\tau^{3}}\left(\mathbf{q},\Omega\right),$
(26)

with the polarization-dependent matrix elements,

$\begin{split}&W^{A}_{\bm{\epsilon},\bm{\epsilon}^{\prime}}=\frac{2}{i\Gamma}\left(\frac{1}{3}\bm{\epsilon}\cdot\bm{\epsilon}^{\prime}-\epsilon_{x}\epsilon_{x}^{\prime*}\right)\;,\\\
&W^{O}_{\bm{\epsilon},\bm{\epsilon}^{\prime}}=-\frac{1}{\sqrt{3}i\Gamma}\left(\epsilon_{y}\epsilon_{y}^{\prime*}-\epsilon_{z}\epsilon_{z}^{\prime*}\right)\;.\end{split}$
(27)

By properly choosing the incoming and the outgoing polarization of the
photons, one can therefore select specific components of the total RIXS
intensity. For instance, by selecting $\bm{\epsilon}$ and
$\bm{\epsilon}^{\prime}$ such that
$\frac{1}{3}\bm{\epsilon}\cdot\bm{\epsilon}^{\prime}=\epsilon_{x}\epsilon_{x}^{\prime*}$
and
$\epsilon_{y}\epsilon_{y}^{\prime*}\neq\epsilon_{z}\epsilon_{z}^{\prime*}$, we
obtain $I_{\text{RIXS}}\propto
S_{\tau^{1},\tau^{1}}\left(\mathbf{q},\Omega\right)$. Similarly, we can obtain
the other diagonal component
$S_{\tau^{3},\tau^{3}}\left(\mathbf{q},\Omega\right)$ and the off-diagonal
component $S_{\tau^{1},\tau^{3}}\left(\mathbf{q},\omega\right)$.

The static correlation functions discussed in Sec. III.2 are then obtained by
an energy integration in the energy window relevant for the orbital
excitations,

$D_{a,b}\left(\mathbf{k}\right)\propto\int d\Omega\
S_{A,B}\left(\mathbf{k},\Omega\right)\;,$ (28)

with $a,b=1,2$, $A,B=\tau^{3},\tau^{1}$ and the correspondences
$1\leftrightarrow\tau^{3}$, $2\leftrightarrow\tau^{1}$. Thus, through diffuse
X-ray scattering, one can probe the momentum distribution of the correlation
functions in a highly anisotropic system such as a compass model. Because the
signal is energy integrated, it should generalize directly to time-resolved
probes, with the same selection rules for the different components of the
orbital fluctuations.

##  V Conclusions

In this work, we designed a time-dependent Ginzburg-Landau formalism suitable
to study the non-equilibrium dynamics of the $120^{\circ}$ compass model, as a
paradigmatic representative for a system in which a multi-minima free energy
landscape for orbital order arises from the order by disorder scenario. We
find that the entropic forces that determine the free energy minima at
equilibrium play an important role also under nonequilibrium conditions: There
is a relatively long-time stabilization of the otherwise unstable unbroken
symmetry state once non-thermal fluctuations of the order parameters have been
activated. Moreover, we find it is generally possible to switch from one
minimum to another considering a time dependent exchange coupling that
emulates the action of a laser-pulse. If, together with this time dependent
exchange interaction, we consider a temperature quench, the orbital switching
is speeded up by one order of magnitude even for small applied external
fields. This mechanism can be viewed as a generalization of the heat-assisted
switching in magnetic memory devices to the orbitally-ordered states.

The present analysis is certainly highly idealized. For generic materials, one
cannot neglect the dynamics of the spin, and further parameters beyond the
effective temperature $T$ may be needed to describe the non-equilibrium state
of the electronic degrees of freedom, in particular for insulators Beaud  _et
al._ (2014). Moreover, microscopic theory (based on spin models, or the full
electronic model) would be needed to get an estimate of phenomenological
constants, such as the damping, and also to understand the ultra-fast regime
in which even the “fast” degrees of freedom are not thermalized Stahl and
Eckstein (2021). Nevertheless, general aspects discussed in the present study
should carry over to more realistic systems:

(1) In particular, the present work shows that for systems which follow the
order by disorder mechanism in equilibrium, forces arising from fluctuations
should also be taken into account in the interpretation of their dynamics.
Corrections to the free energy potential from non-thermal fluctuations become
increasingly more important for faster processes, as in the investigation of
the non-adiabatic switching in Fig. 6. (2) While control parameters such as
the electron temperature $T$ and order parameter fluctuations are rigidly
linked in equilibrium, they independently control the dynamics out of
equilibrium. In the present model the manifestation of this is relatively
simple (non-thermal fluctuations stabilize the symmetry unbroken state even
when $T$ is below the critical temperature).

This fact motivates a search for situations in which other, nontrivial hidden
states can be transiently stabilized by nonthermal fluctuations. Such
phenomena should be investigated, ideally by using scattering techniques to
analyze the time-dependent orbital fluctuations, in transition metal compounds
such as KCuF3. Estimates for a two-band Hubbard model suggest that for Mott-
based orbitally ordered systems electric field-induced changes of the exchange
interaction up to few percent should be possible Eckstein  _et al._ (2017),
potentially allowing to drive coherent dynamics of the order parameter.

## Acknowledgment

We gratefully acknowledge the work done by Aaron Müller at the early stages of
the project. We were supported by the ERC Starting Grant No. 716648.

##  Appendix A Equilibrium solution of the problem

Under equilibrium conditions, the solution of Eq. (7) for the correlation
functions reads:

$D_{ij}\left(\mathbf{k},t\right)=\left(-1\right)^{i+j}\frac{TM_{\bar{i}\bar{j}}\left(\mathbf{k},t\right)}{\text{det}\left[M\left(\mathbf{k},t\right)\right]}\;,$
(29)

where the $M\left(\mathbf{k},t\right)$ matrix has elements as defined in Eq.
(9), so:

$\text{det}\left[M\left(\mathbf{k},t\right)\right]=M_{11}\left(\mathbf{k},t\right)M_{22}\left(\mathbf{k},t\right)-M_{12}^{2}\left(\mathbf{k},t\right)\;,$
(30)

and $\bar{i}=3-i$. Here, we analyze the condition where
$J_{x}=J_{y}=J_{z}\equiv J$, leading to $A_{11}^{0}=A_{22}^{0}\equiv
A^{0}=\frac{3}{2}r\left(T\right)J$ and $A_{12}^{0}=0$. In passing, we notice
that, along specific directions in $k$ space, the correlation function does
not decay to zero when $|\mathbf{k}|\rightarrow\infty$, making the
fluctuations Eq. (5) cutoff dependent. However, the physics we describe is
qualitatively independent of the choice of $\Lambda$. We distinguish two cases
in our analysis, i.e. the zero-temperature condition $T=0$ and the limit of
high temperatures.

At zero temperature, all fluctuations vanish, so that the effective quadratic
couplings Eq. (8) read
$A^{\text{eff}}_{11}=A^{\text{eff}}_{22}=A^{\text{eff}}=A^{0}+4u\left(\bar{\phi}_{1}^{2}+\bar{\phi}_{2}^{2}\right)$
and $A^{\text{eff}}_{12}=0$. The stationarity condition for Eq. (6) implies
that $A^{\text{eff}}=0$, thus
$\bar{\phi}_{1}^{2}+\bar{\phi}_{2}^{2}=-\frac{A^{0}}{4u}$, reflecting the
$\text{O}\left(2\right)$ symmetry of the zero-temperature problem. At high
temperatures $T>T_{c}$, the order parameter is equal to zero, so
$\bar{\phi}_{1}=\bar{\phi}_{2}=0$. In this case, the momentum-integrated
fluctuations become isotropic, $n_{12}=0$ and $n_{11}=n_{22}=n$, so that
$A^{\text{eff}}_{11}=A^{\text{eff}}_{22}=A^{0}+16un$.

## References

  * Basov  _et al._ (2017) D. N. Basov, R. D. Averitt, and D. Hsieh, “Towards properties on demand in quantum materials,” Nature Materials 16, 1077–1088 (2017).
  * Giannetti  _et al._ (2016) Claudio Giannetti, Massimo Capone, Daniele Fausti, Michele Fabrizio, Fulvio Parmigiani, and Dragan Mihailovic, “Ultrafast optical spectroscopy of strongly correlated materials and high-temperature superconductors: a non-equilibrium approach,” Advances in Physics 65, 58–238 (2016).
  * Kirilyuk  _et al._ (2010) Andrei Kirilyuk, Alexey V. Kimel, and Theo Rasing, “Ultrafast optical manipulation of magnetic order,” Rev. Mod. Phys. 82, 2731–2784 (2010).
  * Fausti  _et al._ (2011) D. Fausti, R. I. Tobey, N. Dean, S. Kaiser, A. Dienst, M. C. Hoffmann, S. Pyon, T. Takayama, H. Takagi, and A. Cavalleri, “Light-induced superconductivity in a stripe-ordered cuprate,” Science 331, 189–191 (2011).
  * Sun and Millis (2020) Zhiyuan Sun and Andrew J. Millis, “Transient trapping into metastable states in systems with competing orders,” Phys. Rev. X 10, 021028 (2020).
  * Stojchevska  _et al._ (2014) L. Stojchevska, I. Vaskivskyi, T. Mertelj, P. Kusar, D. Svetin, S. Brazovskii, and D. Mihailovic, “Ultrafast switching to a stable hidden quantum state in an electronic crystal,” Science 344, 177–180 (2014).
  * Ichikawa  _et al._ (2011) Hirohiko Ichikawa, Shunsuke Nozawa, Tokushi Sato, Ayana Tomita, Kouhei Ichiyanagi, Matthieu Chollet, Laurent Guerin, Nicky Dean, Andrea Cavalleri, Shin-ichi Adachi, Taka-hisa Arima, Hiroshi Sawa, Yasushi Ogimoto, Masao Nakamura, Ryo Tamaki, Kenjiro Miyano, and Shin-ya Koshihara, “Transient photoinduced ‘hidden’ phase in a manganite,” Nature Materials 10, 101–105 (2011).
  * Först  _et al._ (2011) M. Först, R. I. Tobey, S. Wall, H. Bromberger, V. Khanna, A. L. Cavalieri, Y.-D. Chuang, W. S. Lee, R. Moore, W. F. Schlotter, J. J. Turner, O. Krupin, M. Trigo, H. Zheng, J. F. Mitchell, S. S. Dhesi, J. P. Hill, and A. Cavalleri, “Driving magnetic order in a manganite by ultrafast lattice excitation,” Phys. Rev. B 84, 241104 (2011).
  * Li  _et al._ (2018) Jiajun Li, Hugo U. R. Strand, Philipp Werner, and Martin Eckstein, “Theory of photoinduced ultrafast switching to a spin-orbital ordered hidden phase,” Nature Communications 9, 4581 (2018).
  * Först  _et al._ (2011) M. Först, C. Manzoni, S. Kaiser, Y. Tomioka, Y. Tokura, R. Merlin, and A. Cavalleri, “Nonlinear phononics as an ultrafast route to lattice control,” Nature Physics 7, 854–856 (2011).
  * Subedi  _et al._ (2014) Alaska Subedi, Andrea Cavalleri, and Antoine Georges, “Theory of nonlinear phononics for coherent light control of solids,” Phys. Rev. B 89, 220301 (2014).
  * Murakami  _et al._ (2017) Yuta Murakami, Naoto Tsuji, Martin Eckstein, and Philipp Werner, “Nonequilibrium steady states and transient dynamics of conventional superconductors under phonon driving,” Phys. Rev. B 96, 045125 (2017).
  * Babadi  _et al._ (2017) Mehrtash Babadi, Michael Knap, Ivar Martin, Gil Refael, and Eugene Demler, “Theory of parametrically amplified electron-phonon superconductivity,” Phys. Rev. B 96, 014512 (2017).
  * Kennes  _et al._ (2017) Dante M. Kennes, Eli Y. Wilner, David R. Reichman, and Andrew J. Millis, “Transient superconductivity from electronic squeezing of optically pumped phonons,” Nature Physics 13, 479–483 (2017).
  * Buzzi  _et al._ (2020) M. Buzzi, D. Nicoletti, M. Fechner, N. Tancogne-Dejean, M. A. Sentef, A. Georges, T. Biesner, E. Uykur, M. Dressel, A. Henderson, T. Siegrist, J. A. Schlueter, K. Miyagawa, K. Kanoda, M.-S. Nam, A. Ardavan, J. Coulthard, J. Tindall, F. Schlawin, D. Jaksch, and A. Cavalleri, “Photomolecular high-temperature superconductivity,” Phys. Rev. X 10, 031028 (2020).
  * Nova  _et al._ (2017) T. F. Nova, A. Cartella, A. Cantaluppi, M. Först, D. Bossini, R. V. Mikhaylovskiy, A. V. Kimel, R. Merlin, and A. Cavalleri, “An effective magnetic field from optically driven phonons,” Nature Physics 2, 132–136 (2017).
  * Singla  _et al._ (2015) R. Singla, G. Cotugno, S. Kaiser, M. Först, M. Mitrano, H. Y. Liu, A. Cartella, C. Manzoni, H. Okamoto, T. Hasegawa, S. R. Clark, D. Jaksch, and A. Cavalleri, “Thz-frequency modulation of the hubbard $u$ in an organic mott insulator,” Phys. Rev. Lett. 115, 187401 (2015).
  * Kaiser  _et al._ (2014) S. Kaiser, S. R. Clark, D. Nicoletti, G. Cotugno, R. I. Tobey, N. Dean, S. Lupi, H. Okamoto, T. Hasegawa, D. Jaksch, and A. Cavalleri, “Optical properties of a vibrationally modulated solid state mott insulator,” Scientific Reports 4, 5556 (2014).
  * Grandi  _et al._ (2021) Francesco Grandi, Jiajun Li, and Martin Eckstein, “Ultrafast mott transition driven by nonlinear electron-phonon interaction,” Phys. Rev. B 103, L041110 (2021).
  * Wall  _et al._ (2009) S. Wall, D. Prabhakaran, A. T. Boothroyd, and A. Cavalleri, “Ultrafast coupling between light, coherent lattice vibrations, and the magnetic structure of semicovalent ${\mathrm{lamno}}_{3}$,” Phys. Rev. Lett. 103, 097402 (2009).
  * Mentink  _et al._ (2015) J. H. Mentink, K. Balzer, and M. Eckstein, “Ultrafast and reversible control of the exchange interaction in mott insulators,” Nature Communications 6, 6708 (2015).
  * Yusupov  _et al._ (2010) Roman Yusupov, Tomaz Mertelj, Viktor V. Kabanov, Serguei Brazovskii, Primoz Kusar, Jiun-Haw Chu, Ian R. Fisher, and Dragan Mihailovic, “Coherent dynamics of macroscopic electronic order through a symmetry breaking transition,” Nature Physics 6, 681–684 (2010).
  * Schäfer  _et al._ (2010) H. Schäfer, V. V. Kabanov, M. Beyer, K. Biljakovic, and J. Demsar, “Disentanglement of the electronic and lattice parts of the order parameter in a 1d charge density wave system probed by femtosecond spectroscopy,” Phys. Rev. Lett. 105, 066402 (2010).
  * Huber  _et al._ (2014) T. Huber, S. O. Mariager, A. Ferrer, H. Schäfer, J. A. Johnson, S. Grübel, A. Lübcke, L. Huber, T. Kubacka, C. Dornes, C. Laulhe, S. Ravy, G. Ingold, P. Beaud, J. Demsar, and S. L. Johnson, “Coherent structural dynamics of a prototypical charge-density-wave-to-metal transition,” Phys. Rev. Lett. 113, 026401 (2014).
  * Beaud  _et al._ (2014) P. Beaud, A. Caviezel, S. O. Mariager, L. Rettig, G. Ingold, C. Dornes, S-W. Huang, J. A. Johnson, M. Radovic, T. Huber, T. Kubacka, A. Ferrer, H. T. Lemke, M. Chollet, D. Zhu, J. M. Glownia, M. Sikorski, A. Robert, H. Wadati, M. Nakamura, M. Kawasaki, Y. Tokura, S. L. Johnson, and U. Staub, “A time-dependent order parameter for ultrafast photoinduced phase transitions,” Nature Materials 13, 923–927 (2014).
  * Maklar  _et al._ (2020) J. Maklar, Y. W. Windsor, C. W. Nicholson, M. Puppin, P. Walmsley, V. Esposito, M. Porer, J. Rittmann, D. Leuenberger, M. Kubli, M. Savoini, E. Abreu, S. L. Johnson, P. Beaud, G. Ingold, U. Staub, I. R. Fisher, R. Ernstorfer, M. Wolf, and L. Rettig, “Nonequilibrium charge-density-wave order beyond the thermal limit,” (2020), arXiv:2011.03230 [cond-mat.mtrl-sci] .
  * Zong  _et al._ (2019a) Alfred Zong, Anshul Kogar, Ya-Qing Bie, Timm Rohwer, Changmin Lee, Edoardo Baldini, Emre Ergecen, Mehmet B. Yilmaz, Byron Freelon, Edbert J. Sie, Hengyun Zhou, Joshua Straquadine, Philip Walmsley, Pavel E. Dolgirev, Alexander V. Rozhkov, Ian R. Fisher, Pablo Jarillo-Herrero, Boris V. Fine, and Nuh Gedik, “Evidence for topological defects in a photoinduced phase transition,” Nature Physics 15, 27–31 (2019a).
  * Zong  _et al._ (2019b) Alfred Zong, Pavel E. Dolgirev, Anshul Kogar, Emre Ergecen, Mehmet B. Yilmaz, Ya-Qing Bie, Timm Rohwer, I-Cheng Tung, Joshua Straquadine, Xirui Wang, Yafang Yang, Xiaozhe Shen, Renkai Li, Jie Yang, Suji Park, Matthias C. Hoffmann, Benjamin K. Ofori-Okai, Michael E. Kozina, Haidan Wen, Xijie Wang, Ian R. Fisher, Pablo Jarillo-Herrero, and Nuh Gedik, “Dynamical slowing-down in an ultrafast photoinduced phase transition,” Phys. Rev. Lett. 123, 097601 (2019b).
  * Kogar  _et al._ (2020) Anshul Kogar, Alfred Zong, Pavel E. Dolgirev, Xiaozhe Shen, Joshua Straquadine, Ya-Qing Bie, Xirui Wang, Timm Rohwer, I-Cheng Tung, Yafang Yang, Renkai Li, Jie Yang, Stephen Weathersby, Suji Park, Michael E. Kozina, Edbert J. Sie, Haidan Wen, Pablo Jarillo-Herrero, Ian R. Fisher, Xijie Wang, and Nuh Gedik, “Light-induced charge density wave in late3,” Nature Physics 16, 159–163 (2020).
  * Dolgirev  _et al._ (2020) Pavel E. Dolgirev, Marios H. Michael, Alfred Zong, Nuh Gedik, and Eugene Demler, “Self-similar dynamics of order parameter fluctuations in pump-probe experiments,” Phys. Rev. B 101, 174306 (2020).
  * Henley (1989) Christopher L. Henley, “Ordering due to disorder in a frustrated vector antiferromagnet,” Phys. Rev. Lett. 62, 2056–2059 (1989).
  * Nussinov and van den Brink (2015) Zohar Nussinov and Jeroen van den Brink, “Compass models: Theory and physical motivations,” Rev. Mod. Phys. 87, 1–59 (2015).
  * Kugel and Khomskii (1973) K. I. Kugel and D. I. Khomskii, “Crystal-structure and magnetic properties of substances with orbital degeneracy,” Zh. Eksp. Teor. Fiz 64, 1429–1439 (1973).
  * Bernevig  _et al._ (2005) B. Andrei Bernevig, Taylor L. Hughes, and Shou-Cheng Zhang, “Orbitronics: The intrinsic orbital current in $p$-doped silicon,” Phys. Rev. Lett. 95, 066601 (2005).
  * Kitaev (2003) A.Yu. Kitaev, “Fault-tolerant quantum computation by anyons,” Annals of Physics 303, 2 – 30 (2003).
  * Eckstein  _et al._ (2017) Martin Eckstein, Johan H. Mentink, and Philipp Werner, “Designing spin and orbital exchange hamiltonians with ultrashort electric field transients,” (2017), arXiv:1703.03269 [cond-mat.str-el] .
  * Polli  _et al._ (2007) D. Polli, M. Rini, S. Wall, R. W. Schoenlein, Y. Tomioka, Y. Tokura, G. Cerullo, and A. Cavalleri, “Coherent orbital waves in the photo-induced insulator-metal dynamics of a magnetoresistive manganite,” Nature Materials 6, 643–647 (2007).
  * Miller  _et al._ (2015) Timothy A. Miller, Ravindra W. Chhajlany, Luca Tagliacozzo, Bertram Green, Sergey Kovalev, Dharmalingam Prabhakaran, Maciej Lewenstein, Michael Gensch, and Simon Wall, “Terahertz field control of in-plane orbital order in la0.5sr1.5mno4,” Nature Communications 6, 8175 (2015).
  * Werner  _et al._ (2017) Philipp Werner, Hugo U. R. Strand, Shintaro Hoshino, and Martin Eckstein, “Ultrafast switching of composite order in ${A}_{3}{\mathrm{c}}_{60}$,” Phys. Rev. B 95, 195405 (2017).
  * Beaulieu  _et al._ (2020) S. Beaulieu, J. Schusser, S. Dong, M. Schüler, T. Pincelli, M. Dendzik, J. Maklar, A. Neef, H. Ebert, K. Hricovini, M. Wolf, J. Braun, L. Rettig, J. Minár, and R. Ernstorfer, “Revealing hidden orbital pseudospin texture with time-reversal dichroism in photoelectron angular distributions,” Phys. Rev. Lett. 125, 216404 (2020).
  * Chen  _et al._ (2019) Yuan Chen, Yao Wang, Chunjing Jia, Brian Moritz, Andrij M. Shvaika, James K. Freericks, and Thomas P. Devereaux, “Theory for time-resolved resonant inelastic x-ray scattering,” Phys. Rev. B 99, 104306 (2019).
  * Mitrano  _et al._ (2019) Matteo Mitrano, Sangjun Lee, Ali A. Husain, Luca Delacretaz, Minhui Zhu, Gilberto de la Peña Munoz, Stella X.-L. Sun, Young Il Joe, Alexander H. Reid, Scott F. Wandel, Giacomo Coslovich, William Schlotter, Tim van Driel, John Schneeloch, G. D. Gu, Sean Hartnoll, Nigel Goldenfeld, and Peter Abbamonte, “Ultrafast time-resolved x-ray scattering reveals diffusive charge order dynamics in la2-xbaxcuo4,” Science Advances 5 (2019), 10.1126/sciadv.aax3346.
  * Mitrano and Wang (2020) Matteo Mitrano and Yao Wang, “Probing light-driven quantum materials with ultrafast resonant inelastic x-ray scattering,” Communications Physics 3, 184 (2020).
  * Parchenko  _et al._ (2020) Sergii Parchenko, Eugenio Paris, Daniel McNally, Elsa Abreu, Markus Dantz, Elisabeth M. Bothschafter, Alexander H. Reid, William F. Schlotter, Ming-Fu Lin, Scott F. Wandel, Giacomo Coslovich, Sioan Zohar, Georgi L. Dakovski, J. J. Turner, S. Moeller, Yi Tseng, Milan Radovic, Conny Saathe, Marcus Agaaker, Joseph E. Nordgren, Steven L. Johnson, Thorsten Schmitt, and Urs Staub, “Orbital dynamics during an ultrafast insulator to metal transition,” Phys. Rev. Research 2, 023110 (2020).
  * Nussinov  _et al._ (2004) Z Nussinov, M Biskup, L Chayes, and J. van den Brink, “Orbital order in classical models of transition-metal compounds,” Europhysics Letters (EPL) 67, 990–996 (2004).
  * van den Brink  _et al._ (1999) Jeroen van den Brink, Peter Horsch, Frank Mack, and Andrzej M. Oleś, “Orbital dynamics in ferromagnetic transition-metal oxides,” Phys. Rev. B 59, 6795–6805 (1999).
  * van den Brink (2004) Jeroen van den Brink, “Orbital-only models: ordering and excitations,” New Journal of Physics 6, 201–201 (2004).
  * Biskup  _et al._ (2005) M. Biskup, L. Chayes, and Z. Nussinov, “Orbital ordering in transition-metal compounds: I. the 120-degree model,” Communications in Mathematical Physics 255, 253–292 (2005).
  * Dmitrašinović  _et al._ (1995) V. Dmitrašinović, J. A. McNeil, and J. R. Shepard, “Goldstone theorem in the gaussian functional approximation to the scalar $\phi^{4}$ theory,” Zeitschrift für Physik C Particles and Fields 69, 359–363 (1995).
  * Hohenberg and Halperin (1977) P. C. Hohenberg and B. I. Halperin, “Theory of dynamic critical phenomena,” Rev. Mod. Phys. 49, 435–479 (1977).
  * Stahl and Eckstein (2021) Christopher Stahl and Martin Eckstein, “Electronic and fluctuation dynamics following a quench to the superconducting phase,” Phys. Rev. B 103, 035116 (2021).
  * Kryder  _et al._ (2008) M. H. Kryder, E. C. Gage, T. W. McDaniel, W. A. Challener, R. E. Rottmayer, G. Ju, Y. Hsia, and M. F. Erden, “Heat assisted magnetic recording,” Proceedings of the IEEE 96, 1810–1835 (2008).
  * Forte  _et al._ (2008) Filomena Forte, Luuk J. P. Ament, and Jeroen van den Brink, “Single and double orbital excitations probed by resonant inelastic x-ray scattering,” Phys. Rev. Lett. 101, 106406 (2008).
  * Marra  _et al._ (2012) Pasquale Marra, Krzysztof Wohlfeld, and Jeroen van den Brink, “Unraveling orbital correlations with magnetic resonant inelastic x-ray scattering,” Phys. Rev. Lett. 109, 117401 (2012).
  * Ament  _et al._ (2011) Luuk J. P. Ament, Michel van Veenendaal, Thomas P. Devereaux, John P. Hill, and Jeroen van den Brink, “Resonant inelastic x-ray scattering studies of elementary excitations,” Rev. Mod. Phys. 83, 705–767 (2011).
  * Marra (2016) Pasquale Marra, “Theoretical approach to direct resonant inelastic x-ray scattering on magnets and superconductors,” (2016), arXiv:1605.03189 [cond-mat.supr-con] .