# Quenches across the self-organization transition in multimode cavities

Tim Keller, Valentin Torggler, Simon B Jäger, Stefan Schütz, Helmut Ritsch, Giovanna Morigi

► To cite this version:

Tim Keller, Valentin Torggler, Simon B Jäger, Stefan Schütz, Helmut Ritsch, et al. Quenches across the self-organization transition in multimode cavities. New Journal of Physics, 2018, 20 (2), pp.025004. 10.1088/1367-2630/aaa161. hal-04876972

HAL Id: hal-04876972

https://hal.science/hal-04876972v1

Submitted on 9 Jan 2025

HAL is a multi-disciplinary open access
archive for the deposit and dissemination of sci-
entific research documents, whether they are pub-
lished or not. The documents may come from
teaching and research institutions in France or
abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire HAL, est
destinée au dépôt et à la diffusion de documents
scientifiques de niveau recherche, publiés ou non,
émanant des établissements d'enseignement et de
recherche français ou étrangers, des laboratoires
publics ou privés.

<figure><img src="image_1.png" /></figure>

Distributed under a Creative Commons Attribution 4.0 International License
PAPER • OPEN ACCESS

Quenches across the self-organization transition in multimode cavities

To cite this article: Tim Keller et al 2018 *New J. Phys.* 20 025004

View the [article online](#) for updates and enhancements.

You may also like

- *Self-ordering and cavity cooling using a train of ultrashort pulses*
  Valentin Torggler, Ivor Kreši, Ticijana Ban et al.

- *Self-Ordering of Anodic Porous Alumina Induced by Local Current Concentration: Burning*
  Sachiko Ono, Makiko Saito and Hidetaka Asoh

- *Extended self-ordering regime in hard anodization and its application to make asymmetric AAO membranes for large pitch-distance nanostructures*
  Minwoo Kim, Yoon-Cheol Ha, Truong Nhat Nguyen et al.
New Journal of Physics

The open access journal at the forefront of physics

Deutsche Physikalische Gesellschaft DPG

IOP Institute of Physics

Published in partnership
with: Deutsche Physikalische
Gesellschaft and the Institute
of Physics

<figure><img src="image_2.png" /></figure>

CrossMark

OPEN ACCESS

RECEIVED

22 September 2017

REVISED

30 November 2017

ACCEPTED FOR PUBLICATION

13 December 2017

PUBLISHED

7 February 2018

Original content from this work may be used under the terms of the [Creative Commons Attribution 3.0 licence](https://creativecommons.org/licenses/by/3.0/).

Any further distribution of this work must maintain attribution to the author(s) and the title of the work, journal citation and DOI.

<figure><img src="image_3.png" /></figure>

PAPER

# Quenches across the self-organization transition in multimode cavities

Tim Keller<fnref n="1" />,<fnref n="4" />, Valentin Torggler<fnref n="2" />,<fnref n="4" />, Simon B Jäger<fnref n="1" />, Stefan Schütz<fnref n="1" />,<fnref n="3" />, Helmut Ritsch<fnref n="2" /> and Giovanna Morigi<fnref n="1" />

E-mail: simonjaeger1@aol.com

**Keywords:** cavity QED, long-range interactions, out-of-equilibrium dynamics, cavity cooling, stochastic differential equations, metastability, self-organization

## Abstract

A cold dilute atomic gas in an optical resonator can be radiatively cooled by coherent scattering processes when the driving laser frequency is tuned close to but below the cavity resonance. When the atoms are sufficiently illuminated, their steady state undergoes a phase transition from a homogeneous distribution to a spatially organized Bragg grating. We characterize the dynamics of this self-ordering process in the semi-classical regime when distinct cavity modes with commensurate wavelengths are quasi-resonantly driven by laser fields via scattering by the atoms. The lasers are simultaneously applied and uniformly illuminate the atoms; their frequencies are chosen so that the atoms are cooled by the radiative processes, and their intensities are either suddenly switched or slowly ramped across the self-ordering transition. Numerical simulations for different ramp protocols predict that the system will exhibit long-lived metastable states, whose occurrence strongly depends on the initial temperature, ramp speed, and the number of atoms.

## 1. Introduction

Laser light creates an attractive optical potential for cold atoms when far detuned below an optical transition. Such potential can be significantly enhanced if the light is confined by an optical resonator [1–4]. In addition, if the laser illuminates the atoms, trapping is induced by a dynamical optical potential emerging from the interference between the scattered light and the laser, which tends to order the particles at the maxima of the intensity [4, 5]. The interference contrast and thus the trapping depends on the relative positions of the scattering atoms. Therefore, this phenomenon can be also understood in terms of an effective long-range force, which is mediated by the collectively scattered photons [5–9]. This force also has a dissipative component, which is due to the dissipative nature of the resonator and which cools the atoms when the pump is tuned below the cavity resonance [3, 10]. Theoretical studies with single-mode resonators have predicted that this dissipation can establish long-range correlations and support the onset of metastable ordered structures [11, 12].

In a multimode cavity and for several illumination frequencies, competing ordering processes are present and lead to richer phase dynamics. In a two-mode cavity, like the one depicted in figure 1(a), the transition to self-organization can be a phase transition of the first or second order depending on the laser intensities and on their relative strength [13]. The corresponding self-ordered phases can exhibit superradiant scattering either in one or in both cavity modes, as illustrated in figure 1(b), while the asymptotic distribution of the atoms can be thermal provided that the lasers' frequencies are suitably chosen [13]. In our example the particles can order in a lattice at a given length scale $\lambda$ and/or on a lattice with half the period $\lambda/2$. For these settings we numerically analyze the semi-classical dynamics following sudden quenches or slow ramps of the laser intensities across the thresholds separating the homogeneous from one of the self-organized phases. We describe the evolution by stochastic differential equations, which correspond to the Fokker-Planck equation derived in [14] for a similar system. We find that even at very long times the atoms' spatial distribution strongly depends on the initial

<footnote n="1">¹ Theoretische Physik, Universität des Saarlandes, D-66123 Saarbrücken, Germany</footnote>

<footnote n="2">² Institut für Theoretische Physik, Universität Innsbruck, A-6020 Innsbruck, Austria</footnote>

<footnote n="3">³ icFRC, IPCMS (UMR 7504) and ISIS (UMR 7006), University of Strasbourg and CNRS, F-67000 Strasbourg, France</footnote>

<footnote n="4">⁴ First authors</footnote>
<figure><img src="image_4.png" /><figcaption>**Figure 1.** (a) Cold atoms are confined within an optical cavity and move along the cavity axis (x axis). They coherently scatter photons from transverse lasers with rescaled amplitudes $\alpha_1$ and $\alpha_2$ into the (correspondingly) resonant cavity modes with spatial mode functions $\cos(kx)$ (red) and $\cos(2kx)$ (green) and loss rates $\kappa_1$ and $\kappa_2$, respectively. (b) Sketch of the atomic density distribution $n(x)$ along the cavity axis (in units of $1/k$) for the four possible stationary self-organized orders. On the right we report the corresponding values of the quantities $\Theta_1$ and $\Theta_2$, which signal the Bragg order in modes 1 and 2, respectively. See text and [13] for details.</figcaption></figure>

temperature, ramp speed, and quench protocol, such that the system gets trapped in long-lived metastable states. In particular, for quenches starting with ensembles at low temperatures, the buildup of long-range order requires longer times than that for higher initial temperatures does.

Our work is organized as follows: in section 2 we introduce the system and the semi-classical equations describing the dynamics. The atoms' stationary properties are then summarized in a phase diagram, which is derived from [13]. In section 3 we numerically study the real-time dynamics when the parameters are varied within the phase diagram according to different quench protocols. In section 4 we analyze the dynamics of the distributions from the spatially homogeneous to the organized ones with different momentum widths. In section 5 we compare the predictions of the stochastic differential equations we employ with an extended approach including the dynamical evolution of the field modes introduced in [15, 16]. Conclusions are drawn and future perspectives are discussed in section 6.

## 2. Semi-classical dynamics

The system we consider consists of a gas of $N$ cold atoms with mass $m$, which are trapped inside a high-finesse optical resonator and coherently scatter laser light into the cavity modes. The atomic motion is confined along the cavity axis (here the x axis) by a tight external dipole trap [17, 18] and is here described in the semi-classical limit.

The geometry of the setup is illustrated in figure 1. Lasers with (rescaled) intensity $\alpha_n$ propagate in a direction orthogonal to the cavity axis and are quasi-resonant with the standing wave cavity modes $\cos(nkx)$ with frequency $\omega_{c,n}$ and wave number $nk$ ($n = 1, 2$)$^5$. The lasers have frequency $\omega_{p,n}$ and linear polarization, which is parallel to that of the corresponding cavity mode. Each pair of laser and cavity mode couples to an atomic dipolar transition at frequency $\omega_{a,n}$, where $\Omega_{p,n}$ and $g_n$ are the laser and vacuum Rabi frequency, respectively. Spontaneous scattering processes are suppressed when the absolute value of the detuning $|\Delta_{a,n} = \omega_{p,n} - \omega_{a,n}|$ exceeds the coupling strength and the detuning $|\Delta_n = \omega_{p,n} - \omega_{c,n}|$ between the laser and cavity mode by orders of magnitude: $|\Delta_{a,n}| \gg \Omega_{p,n}$, $g_n$, $|\Delta_n|$. The relevant dissipative process is given by cavity decay, and we denote by $\kappa_n$ the loss rate of cavity mode $n = 1, 2$.

In the so-called bad cavity limit, assuming that the cavity field loss rates are faster than the rate of the dynamics of the atomic motion, one can eliminate the cavity field variables from the equations of motion of the atoms by means of coarse graining in time. This gives rise to an effective model, where the atoms experience a long-range interaction mediated by the cavity photons, while retardation effects and fluctuations of the cavity field are responsible for friction forces and diffusion. In the semi-classical limit one can derive a Fokker-Planck equation for the atoms' position and momentum distribution, assuming that the single-atom momentum distribution has a width $\Delta p$ which, at all instants of time, is orders of magnitude greater than the photon recoil $\hbar k: \Delta p \gg \hbar k$ [13, 14, 20]. The corresponding stochastic differential equations read as

<footnote>$^5$This can be realised by assuming $\omega_{c,2} = 2\omega_{c,1} \equiv 2\omega_c$, which yields $k_2 = 2k_1 \equiv 2k$. Another possible realisation, where $\omega_{c,1} \approx \omega_{c,2}$, has been discussed in [13]; it uses two optical single-mode cavities crossing at an angle of 60°. For a similar experimental setup see also [19].</footnote>
$$dx_j = \frac{p_j}{m}dt, \qquad (1a)$$

$$dp_j = (F_{j,ad} + F_{j,ret})dt + dW_j^{(1)} + dW_j^{(2)}, \qquad (1b)$$

where

$$F_{j,ad} = - \sum_{n=1,2} 2\hbar nk \frac{\alpha_n}{\hbar \beta_n} \sin(nkx_j) \Theta_n, \qquad (2)$$

$$F_{j,ret} = - \sum_{n=1,2} \frac{\hbar(nk)^2}{m} \alpha_n \frac{\kappa_n}{-\Delta_n} \sin(nkx_j) \frac{1}{N} \sum_{l=1}^{N} p_l \sin(nkx_l), \qquad (3)$$

and

$$\alpha_n = \frac{4NS_n^2 \Delta_n^2}{(\Delta_n^2 + \kappa_n^2)^2}, \qquad (4)$$

$$\beta_n = \frac{-4\Delta_n}{\hbar(\Delta_n^2 + \kappa_n^2)}. \qquad (5)$$

Here, $S_n = g_n\Omega_n/\Delta_{a,n}$ is the amplitude of coherent scattering by a single atom and has the dimension of a frequency, while $dW_j^{(1)}$ and $dW_j^{(2)}$ in equation (1b) describe Wiener processes, which fulfill $\langle dW_i^{(n)} \rangle = 0$ and $\langle dW_i^{(n)} dW_j^{(m)} \rangle = 2D_{ij}^n \delta_{nm} dt$ ($n, m = 1, 2$ and $i, j = 1, \dots, N$). Here,

$$D_{ij}^n = (\hbar nk)^2 \frac{\alpha_n}{\hbar N \beta_n} \left( \frac{\kappa_n}{-\Delta_n} \right) \sin(nkx_i) \sin(nkx_j). \qquad (6)$$

Finally, the parameter

$$\Theta_n = \frac{1}{N} \sum_{i=1}^{N} \cos(nkx_i) \qquad (7)$$

quantifies the Bragg ordering of the atoms in the cavity mode with wave number $nk$. In particular, $|\Theta_n| = 1$ when the atoms are localized either at the maxima or at the minima of $\cos(nkx)$, which is the configuration which maximizes the intracavity field intensity. We identify $\Theta_n$ with the order parameter for self-organization in the corresponding cavity mode [13]. Below, we denote by 'long-wavelength order' any configuration with a non-vanishing value of $\Theta_1$, corresponding to a Bragg grating with period $\lambda = 2\pi/k$. Similarly, 'short-wavelength order' refers to a configuration with $\Theta_2 \neq 0$, corresponding to a Bragg grating with $\lambda/2$. Note that here and in the rest of the paper we discard the dynamical Stark shift of the cavity frequency assuming that this is much smaller than the cavity mode linewidth $Ng_n^2/|\Delta_{a,n}| \ll \kappa_n$. For details refer to [13].

## 2.1. Stationary states

Analysis of the Fokker-Planck equation on the basis of equation (1) allows one to identify the conditions for the existence of a stationary state. This state exists provided that $\Delta_n < 0$ and $\beta_1 = \beta_2 \equiv \beta$ (see equation (5)). In this case the atoms' distribution in the steady state reads as [13]

$$f_{\text{st}}(x_1, p_1, \dots, x_N, p_N) = \frac{\exp(-\beta H_{\text{eff}})}{Z(\beta)} \qquad (8)$$

where $H_{\text{eff}}$ is the effective Hamiltonian derived after eliminating the cavity field variables,

$$H_{\text{eff}} = \sum_{j=1}^{N} \frac{p_j^2}{2m} - \sum_{n=1,2} N \frac{\alpha_n}{\beta_n} \Theta_n^2, \qquad (9)$$

while $Z(\beta)$ denotes the partition function:

$$Z(\beta) = \frac{1}{\Delta^N} \int_{-\infty}^{\infty} dp_1 \dots \int_{-\infty}^{\infty} dp_N \int_0^{\lambda} dx_1 \dots \int_0^{\lambda} dx_N e^{-\beta H_{\text{eff}}}, \qquad (10)$$

with $\lambda = 2\pi/k$ and $\Delta = 2\pi\hbar$ as the single particle unit phase space volume. In the following we assume that the cavity decay rates are equal,

$$\kappa_1 = \kappa_2 := \kappa, \qquad (11)$$

so that the condition for the existence of the stationary state in equation (8) becomes

$$\Delta_1 = \Delta_2 := \Delta_c < 0. \qquad (12)$$

The phase diagram of the system can be determined by using the steady state, equation (8), in the form of a thermal state. On the basis of this observation we introduce the temperature $T$ of the stationary state, which is defined as
$$k_B T = \beta^{-1} = \frac{\hbar(\Delta_c^2 + \kappa^2)}{-4\Delta_c}, \qquad (13)$$

with $k_B$ as the Boltzmann constant. The steady-state temperature $T$ has the same functional dependence on $\Delta_c$ and $\kappa$ as for a single-mode cavity [7, 14]. We can further define the free energy per particle $\mathcal{F}$ using formal equivalence with the canonical ensemble of equilibrium statistical mechanics [7]:

$$\mathcal{F} = -\frac{1}{N\beta} \ln(\mathcal{Z}(\beta)). \qquad (14)$$

Following the procedure detailed in [7, 12, 13] we can calculate $\mathcal{F}$ as

$$\mathcal{F} = \inf_{\theta_1, \theta_2} F(\theta_1, \theta_2) \qquad (15)$$

with

$$\beta F(\theta_1, \theta_2) = \frac{\ln(\pi \hbar \omega_r \beta)}{2} + \alpha_1 \theta_1^2 + \alpha_2 \theta_2^2 - \ln(I(\theta_1, \theta_2)) \qquad (16)$$

and

$$I(\theta_1, \theta_2) = \int_0^{2\pi} dq \ e^{2\alpha_1 \theta_1 \cos(q) + 2\alpha_2 \theta_2 \cos(2q)}. \qquad (17)$$

Here $\omega_r = \hbar k^2/(2m)$ is the recoil frequency. We determine the value of $\mathcal{F}$ in an appropriately defined thermodynamic limit, which consists in keeping $\alpha_n$ constant for $N \to \infty$. The global minima of $F$ are the resulting stationary phases. The corresponding points $(\theta_{1,\min}, \theta_{2,\min})$ where $F$ achieves its minimum are the stationary values for the order parameters $\Theta_1 = \theta_{1,\min}$ and $\Theta_2 = \theta_{2,\min}$; they are determined by $\alpha_1$ and $\alpha_2$. When the fields are sufficiently weak, then $\Theta_1 = \Theta_2 = 0$, the density is homogeneous and there is no structural order. We call this phase paramagnetic, borrowing the notation of the generalized Hamiltonian mean-field model (GHMF) [21–23] to which this model can be mapped. The possible ordered phases in the steady state are illustrated in figure 1(b) and take one of four sets of values. In particular, the ferromagnetic phase is characterized by (i) $\Theta_1 > 0$, $\Theta_2 > 0$, and (ii) $\Theta_1 < 0$, $\Theta_2 > 0$, exhibiting Bragg order in both cavity modes. In contrast, the nematic phases, (iii) and (iv), are characterized by no order in the long-wavelength mode, $\Theta_1 = 0$, while $\Theta_2$ can be either negative or positive.

We note that the spatial distributions depicted in figure 1(b) corresponding to phases (i), (ii), (iii), and (iv) are only possible configurations out of many. For example configuration (i) can also correspond to all atoms sitting at one site $kx = 0$. Order here refers to Bragg gratings corresponding to the long- and short-wavelength modes. No long-wavelength order is found in the case of $\Theta_1 = 0$, where photons scattered into the long-wavelength mode destructively interfere. However, $\Theta_1$ cannot give more detailed information about the positions of the atoms in the long-wavelength mode. The same is true for short-wavelength order and $\Theta_2$.

The resulting phase diagram in the $\alpha_1 - \alpha_2$ plane, shown in figure 2, reproduces that in [13]. The phases are separated by either first- or second-order transitions, which have been determined using Ehrenfest’s criterion [23]. The shaded areas show stability regions in which the free energy has a local minimum that corresponds to the paramagnetic (dark gray region) and nematic (light gray region) phases. Examples of the free-energy landscape in the $\Theta_1 - \Theta_2$ plane are shown in subplots (b) and (c). Subplot (b) corresponds to the parameters of the red bullet labeled (b) in subplot (a): Here, the free energy exhibits two symmetric global minima which correspond to the ferromagnetic phase. In subplot (c), which corresponds to the parameters of the red bullet labeled (c), there is an additional local minimum corresponding to a nematic phase. In the latter there is only ordering in the short-wavelength lattice, while $\Theta_1 = 0$. We call this region *bistable*, which refers to the existence of a second, metastable state in which the system can be dynamically trapped.

### 3. Dynamics of self-organization

We now examine the dynamics of the system when the values of $\alpha_1$ and $\alpha_2$ are varied as a function of time. Experimentally, this corresponds to varying the pump laser intensities or their detuning with respect to the cavity mode frequencies. At time $t = 0$ we assume that the system is prepared in the stationary state of a paramagnetic phase, described by the distribution in equation (8), by setting $\alpha_n = \alpha_{ni} \ll 1$ in equation (9) ($n = 1, 2$). The values $\alpha_n$ appearing in the equations of motion (1) are then varied in time, by performing either (i) a sudden quench, i.e. a sudden switching of the values to $\alpha_{1f}$ and $\alpha_{2f}$, or (ii) a slow quench, which consists in varying $\alpha_n(t)$ monotonically and continuously in time towards the final values $\alpha_{1f}$ and $\alpha_{2f}$. We choose the final values $\alpha_{nf}$ in the ferromagnetic phase. The quench protocols we consider are illustrated by the green lines in figure 2(a): for sudden quenches, the initial and final values are two points connected by a green line. A slow quench sweeps across the intermediate points along the line. We are interested in determining the dynamics leading to the steady state.
<figure><img src="image_5.png" /><figcaption>Figure 2. (a) Phase diagram of the stationary phases, corresponding to the global minima of equation (16), in plane $α_1 - α_2$. Blue (yellow) lines indicate second (first) order phase transitions. The light gray area (dark gray area) within the ferromagnetic phase indicates the parameter region where the nematic (paramagnetic) phase is the local minimum of the free energy. The red circles labeled (b) and (c) indicate the parameters for which the contour plots of free energy in subplots (b) and (c) are shown in the $Θ_1 - Θ_2$ landscape. The free energy in subplot (b) exhibits two global minima at $Θ_1 = ±0.92$ and $Θ_2 = 0.73$; in (c) the two global minima of the ferromagnetic phase are at $Θ_1 = ±0.97$ and $Θ_2 = 0.88$, and the local minimum in the nematic phase is at $Θ_1 = 0$ and $Θ_2 = -0.83$ (the contour of $F$ is below a convenient threshold). The green dash-dot lines in subplot (a) illustrate the paths of the quench protocols discussed in section 3. Circle (d) indicates the parameters of the quench discussed in section 3.3.</figcaption></figure>

In what follows we perform numerical simulations of equation (1) using the parameters of a gas of $^{85}$Rb atoms. In particular, we take $k = 2\pi/\lambda$ with $\lambda = 780$ nm as the wavelength of the $D_2$ line. The corresponding recoil frequency is $\omega_r = 2\pi \times 3.86$ kHz. The cavity linewidth is taken to be $\kappa = 2\pi \times 1.5$ MHz, so that $\kappa \approx 388.6\omega_r$. A possible realization of the two-mode setup here considered has been discussed in [13, 19].

### 3.1. Sudden quench into the ferromagnetic phase

We first consider sudden quenches from $α_{1i}$, $α_{2i}$ in the paramagnetic phase to $α_{1f}$, $α_{2f}$ in the ferromagnetic phase, keeping $α_{1i}/α_{2i} = α_{1f}/α_{2f} = 5$. The initial values are vanishingly small and the atoms are at the corresponding stationary distribution, which is a thermal distribution at the temperature determined by the corresponding detuning, equation (13), with homogeneous density. The detuning before and after the quench is taken to be equal; thus it is expected that the atoms reach a thermal distribution with the same temperature as the initial state.

Figure 3 shows the distribution $P_t(Θ)$ for order parameters $Θ_1$ and $Θ_2$ as a function of time for $(α_{1f}, α_{2f}) = (2.5, 0.5)$. It is defined as a time sequence of normalized histograms:

$$P_t(\Theta) = \frac{\# \text{ trajectories with } \Theta(t) \in [\Theta - \Delta_{\Theta}/2, \Theta + \Delta_{\Theta}/2]}{\# \text{ trajectories} \times \Delta_{\Theta}}, \quad (18)$$

where $Θ$ is calculated on each trajectory of the simulations with the stochastic differential equations and its value is determined according to the precision $\Delta_{\Theta}$ of the grid in $Θ$. We observe that at a given time scale of the order of $10^2/\kappa$, $P_t(Θ_1)$ splits into two branches corresponding to two possible orders in the long-wavelength lattice. This symmetry breaking is well known from the single-mode case [5]. The order parameter of the short-wavelength mode $Θ_2$, which is weakly pumped, substantially grows to a positive value long after the symmetry breaking. The fact that $P_t(Θ_2)$ vanishes for negative $Θ_2$ values comes from the ordering of the atoms close to the anti-nodes of the dominant long-wavelength mode field $\cos(kx)$ (see figure 1).

The distributions $P(Θ_n)$ at the asymptotics are reported in the right panels of figure 3. They are obtained by averaging $P_t(Θ_n)$ over times $t \ge 10^6/\kappa$, where a stable configuration is reached. Formally
<figure><img src="image_6.png" /><figcaption>Figure 3. Dynamics following a sudden quench from $\alpha_{1i}$, $\alpha_{2i} \ll 1$ to $\alpha_{1f} = 2.5$, $\alpha_{2f} = 0.5$ with the detuning kept constant $\Delta_c = -\kappa$. The left panels display the contour plots of distribution $P(\Theta)$ (equation (18)) for $\Theta = \Theta_1$ and $\Theta = \Theta_2$ as a function of time (in units of $1/\kappa$). The distribution is extracted from numerical simulations using equation (1) for $N = 100$ atoms and 1000 trajectories. The grid in $\Theta$ for the left panels has a minimum step $\Delta\Theta = 2/111$; the grayscale gives the relative weight. The right panels display the distributions $P(\Theta_1)$ and $P(\Theta_2)$ as a function of $\Theta_1$ and $\Theta_2$, respectively (see equation (19)). Here, the time average is obtained for $N_t = 113$ instants of time chosen between $t_1 = 10^6/\kappa$ and $t_f = 3.77 \times 10^6/\kappa$. The vertical dashed lines in the left panels mark the instant of time $t_1 = 10^6/\kappa$.</figcaption></figure>

$$P(\Theta) = \sum_{i=1}^{N_t} P_{t_i}(\Theta)/N_t, \qquad (19)$$

where $N_t$ is the number of instants of times at which the distribution is sampled in the interval $[t_1, t_f]$, with $t_1 = 10^6/\kappa$ and $t_f = t_{N_t} > t_1$. Comparing the widths of the distributions in the right panels of figure 3 one observes that after sufficiently long times the long-wavelength order parameter fluctuates less than the short-wavelength order parameter.

Figures 4(a)–(b) display the dynamics of the mean absolute value of the order parameters for different values of $N$. Figure 4(c) shows the time evolution of the fluctuations of the order parameters $\delta\Theta_n = \sqrt{\langle\Theta_n^2\rangle} - \langle|\Theta_n|\rangle^2$ for $N = 100$ particles. The order parameters asymptotically tend to the values predicted by the free energy, indicated by the horizontal dashed lines, for a time scale of the order of $10^6/\kappa$. Meanwhile the fluctuation $\delta\Theta_1$ relaxes to a much smaller value than the asymptotic value of $\delta\Theta_2$ reproducing the widths of the distributions in the right panels in figure 3. The time evolution of $\langle|\Theta_1|\rangle$, in particular, is reminiscent of the one observed for quenches into the ferromagnetic phase in a single-mode resonator [11]. It can be separated into three stages which we call (in order of their temporal appearance) (i) violent relaxation, which corresponds to an exponential increase of the absolute value of the order parameter $\langle|\Theta_1|\rangle$; (ii) transient dynamics, which corresponds to power-law scaling with time; and (iii) the relaxation phase, where the mean values tend exponentially towards the asymptotic value. Violent relaxation can be described by a mean-field model [12]; in the transient stage coherent dynamics prevails, while the relaxation stage is dominated by dissipation [11]. The transient and relaxation stages are characterized by time scales which increase with $N$ but have different functional dependence [12]. The time scale $10^6/\kappa$ can here be identified as the one at which the asymptotic state is reached for $N \lesssim 200$, while for larger numbers of particles longer time scales shall be considered.

Interestingly, in the transient phase there is ordering only in the long-wavelength mode of the cavity, while ferromagnetic order is finally established by dissipation on a longer time scale. The metastable phase of the transient dynamics can therefore be denoted by ‘nematic’: its lifetime increases with $N$ and for $N \sim 200$ it is of the order of $t \sim 10^4/\kappa$. However, this metastable nematic state cannot be understood in terms of the landscape of the free energy, but rather seems to exhibit the features of a quasi-stationary state due to long-range coherent dynamics analogous to that reported in [22]. This conjecture is also supported by the behavior of the single-particle kinetic energy and of the kurtosis $K = \langle p^4 \rangle / \langle p^2 \rangle^2$, which are shown in figure 5. The latter quantifies the deviation of the momentum distribution from a Gaussian one, for which it takes a value of $K_{Gauss} = 3$. For these quantities we observe that in the metastable nematic phase the kinetic energy grows, while the distribution is non-thermal. Ordering in the second, short-wavelength lattice is accompanied by cooling into a thermal distribution.

We now compare the numerical results with the analytic theory for different quenches with the same initial values of $\alpha_{1i}$, $\alpha_{2i} \ll 1$ but with different endpoints $\alpha_{1f}$, $\alpha_{2f}$. We take different endpoints from the paramagnetic to the ferromagnetic phase, under the constraint $\alpha_{1f}/\alpha_{2f} = 5$. The circles in figure 6 correspond to the numerical results for 100 particles at time $t_f > 10^6/\kappa$, where we expect the system to reach the steady state.
<figure><img src="image_7.png" /><figcaption>Figure 4. Dynamics of (a) $\langle|\Theta_1|\rangle$, (b) $\langle|\Theta_2|\rangle$, and (c) their fluctuations $\delta\Theta_n = \sqrt{\langle\Theta_n^2\rangle - \langle|\Theta_n|\rangle^2}$ as a function of time (in units of $1/\kappa$). The parameters and quench protocol are the same as in figure 3; the curves are however evaluated for different numbers of atoms and of trajectories. In (a) and (b) the data correspond to $N = 25, 50, 100, 200$ particles (see legends for color code) and respectively 1000, 500, 250, 125 trajectories. The horizontal dashed lines indicate the values predicted by the global minimum of the free energy in equation (16). The finite values of the order parameters at $t = 0$ are due to finite size effects, $\langle|\Theta_n(0)|\rangle = 1/\sqrt{\pi N}$. The curves in (c) are calculated for $N = 100$ and 250 trajectories.</figcaption></figure>

These are in good agreement with the analytical results (dashed lines) based on an evaluation of the corresponding observables at the global minimum of the free energy. The interval where $\langle|\Theta_n|\rangle$ grows monotonically from $\sim 1/\sqrt{N}$ to the value of the ferromagnetic phase is expected to shrink as $N$ is increased, in agreement with a second-order phase transition at the thermodynamic limit. Further information on the onset of this ferromagnetic order can be gained by the probability $P_t(\Theta_2 < 0)$ that $\Theta_2$ is negative at $t$:

$$P_t(\Theta_2 < 0) = \int_{-1}^{0} d\Theta_2 P_t(\Theta_2). \qquad (20)$$

We note that in the paramagnetic phase (homogeneous spatial distribution) we expect $P_t(\Theta_2 < 0) \simeq 0.5$. In contrast, due to the given mode structure we expect that $P_t(\Theta_2 < 0) \simeq 0$ for long-wavelength ordering in the ferromagnetic phase. Indeed, as $\alpha_{1f}$ increases across the critical value, $P_t(\Theta_2 < 0)$ quickly drops down to zero.

## 3.2. Sudden quenches into the bistable phase

We now turn to the dynamics following sudden quenches from the paramagnetic to the ferromagnetic phase but following the right path of figure 2(a), which consists in equal effective pumping $\alpha_{1i}/\alpha_{2i} = \alpha_{1f}/\alpha_{2f} = 1$. In this parameter region (the bistable phase) the free energy exhibits a local minimum, which is nematic. As in the
<figure><img src="image_8.png" /></figure>

<figure><img src="image_9.png" /><figcaption>**Figure 5.** Dynamics of (a) the single-particle kinetic energy $\langle p^2/2m \rangle$ (in units of $\hbar\omega_r$) and (b) kurtosis $K = \langle p^4 \rangle / \langle p^2 \rangle^2$ for $N = 25, 50, 100, 200$ particles (see legends) and respectively 1000, 500, 250, 125 trajectories. The other parameters and quench protocol are the same as in figure 4. The horizontal dashed line in (a) indicates the asymptotic value predicted by equation (13).</figcaption></figure>

<figure><img src="image_10.png" /><figcaption>**Figure 6.** Asymptotic values of $\langle |\Theta_1(t_f)| \rangle$, $\langle |\Theta_2(t_f)| \rangle$, and $P_{t_f}(\Theta_2 < 0)$, equation (20), as a function of $\alpha_{1f}$. The quenches start from the same initial values in the paramagnetic phase ($\alpha_{1i}$, $\alpha_{2i} \ll 1$ and $\Delta_c = -\kappa$) and end up with different values $\alpha_{1f}$, $\alpha_{2f}$ with $\alpha_{1f} = 5\alpha_{2f}$ (lying along the left green line in figure 2(a)) and $\Delta_c = -\kappa$. The circles correspond to the results of the numerical simulations at $t_f = 3.77 \times 10^6/\kappa$ with $N = 100$ particles and 250 trajectories. The dashed lines indicate the predictions of the global minima of equation (16).</figcaption></figure>

previous case, the initial values $\alpha_{1i}$, $\alpha_{2i}$ are vanishingly small and the atoms are at the corresponding stationary distribution, whose temperature is determined by the detuning $\Delta_c$ and whose spatial density is homogeneous. The quench is performed by switching the laser power while keeping the detuning constant; thus the atoms should reach a thermal distribution with the same temperature as the initial state.

Figure 7 shows the time evolution of the trajectories' $\Theta$-distribution for $\alpha_{1f} = \alpha_{2f} = 2$ and $\Delta_c = -\kappa$. In contrast to the previous section, here a finite fraction of trajectories gets trapped in the nematic phase with a vanishing value of $\Theta_1$ and a finite probability that $\Theta_2$ takes negative values. This is visible in the small extra peaks in $P(\Theta_1)$ and $P(\Theta_2)$ (right panels). The trapping occurs at the time scale of the violent relaxation, and it seems stable over times of the order of $10^6/\kappa$. We conjecture that it also persists at asymptotic times. In figure 8 the time evolution of the mean absolute value of the order parameters is shown for different numbers of particles. While $\langle |\Theta_2| \rangle$ reaches the same stationary value (in reality its value decreases slightly with $N$), the asymptotic value
<figure><img src="image_11.png" /><figcaption>**Figure 7.** Dynamics following a sudden quench from $\alpha_{1f}$, $\alpha_{2f} \ll 1$ to $\alpha_{1f} = \alpha_{2f} = 2$ with the detuning kept constant $\Delta_c = -\kappa$ (red circle (c) in figure 2(a)). The left panels display the contour plot of the distribution $P_t(\Theta)$ (equation (18)) for $\Theta = \Theta_1$ and $\Theta = \Theta_2$ as a function of time (in units of $1/\kappa$). The distribution is extracted from the numerical simulations using equation (1) for $N = 100$ atoms and 1000 trajectories. The right panels display distributions $P(\Theta_1)$ and $P(\Theta_2)$ as a function of $\Theta_1$ and $\Theta_2$, respectively (see equation (19)). See figure 3 for further details.</figcaption></figure>

of $\langle|\Theta_1|\rangle$ decreases as $N$ grows. This suggests that the probability that the dynamics gets trapped in the local minimum increases with the number of particles. The asymptotic value of $\delta\Theta_1 = \sqrt{\langle\Theta_1^2\rangle - \langle|\Theta_1|\rangle^2}$ in subplot (c) reflects the contribution of these trajectories.

The mean single-particle kinetic energy and kurtosis are shown in figure 9. From their behavior we infer that the metastable nematic state does not significantly deviate from a thermal distribution with the expected asymptotic temperature (equation (13)).

Peculiar features of these dynamics become visible when inspecting the probability $P_t(\Theta_2 < 0)$ at the asymptotics and as a function of $\alpha_{1f}$ in figure 10. As in figure 6, it vanishes upon leaving the paramagnetic phase, but increases again as the $\alpha_{1f}$, $\alpha_{2f}$ chosen are deeper into the bistable phase of figure 2(a). Correspondingly, $\langle|\Theta_1|\rangle$ starts to decrease as $\alpha_{1f}$ increases, which suggests that from this point on the depth of the local minimum grows. The value of the order parameter $\langle|\Theta_2|\rangle$ at which $P_t(\Theta_2 < 0)$ starts to grow again signifies a threshold, above which the local minimum is sufficiently deep to stably trap particles.

### 3.3. Slow ramp into the bistable phase

We now consider linear ramps of $\alpha_n(t)$ across the transition region separating the paramagnetic from the bistable region. The ramp protocols have duration $\tau$ and sweep between the values $[\epsilon, \alpha_{nf}]$, with $\epsilon \ll 1$. In particular, $\alpha_n(t) = \epsilon + \alpha_{nf} \frac{t}{\tau}$ if $t \in [0, \tau]$, while for $t > \tau$, $\alpha_n(t)$ is constant and equal to $\alpha_{nf}$. Note that a sudden quench is the limit $\tau \to 0$ of a linear quench. We choose to vary the values of $\alpha_n(t)$ along the rightmost green line in figure 2(a), so that $\alpha_1(t) = \alpha_2(t)$ at all instants of time, with $\alpha_{nf}$ in the bistable phase. We further keep $\Delta_c$ constant, and vary only the pump intensity. This means that the asymptotic temperatures at each value of $\alpha_n$ are equal.

Figure 11 shows the dynamics of the mean absolute value of the order parameters for $\alpha_{1f} = \alpha_{2f} = 2$ for linear ramps with different durations $\tau$. The dynamics following the sudden quench (figures 8(a) and (b)) is shown for comparison (blue curve). We observe that the dynamics of the order parameters exhibits an exponential increase which occurs almost simultaneously for both $\langle|\Theta_1|\rangle$ and $\langle|\Theta_2|\rangle$. This behavior seems to be initiated at the instant of time when the parameters $\alpha_n(t)$ cross the critical point of the phase diagram. Moreover, for sufficiently slow ramps $\langle|\Theta_1|\rangle$ approaches the asymptotic value of the free energy's global minimum, signaling stationary long-wavelength order.

We further note that for $\tau \lesssim 10^3/\kappa$ the order parameters undergo a three-stage dynamics, as for the sudden quench (we attribute the fluctuations to the statistics of the trajectories). For slower ramps, the mean value of the order parameters tends exponentially towards the steady state, which approaches the free energy's global minimum in equation (16) for $\tau > 10^4/\kappa$. We believe that this behavior is determined by the ramp duration $\tau$ with respect to the time scale of the transient dynamics, and thus by the time the parameters $\alpha_n(t)$ spend close to the transition point. This conjecture is supported by the analysis of the time evolution of the single-particle kinetic energy shown in figure 12, which corresponds to the curves in figure 11. For faster ramps it is similar to the sudden quench, exhibiting first a violent relaxation followed by a time interval when the dynamics is predominantly coherent, and then an exponential decay to the steady-state value due to cavity cooling. In
<figure><img src="image_12.png" /></figure>

<figure><img src="image_13.png" /></figure>

<figure><img src="image_14.png" /><figcaption>**Figure 8.** Dynamics of (a) $\langle|\Theta_1|\rangle$, (b) $\langle|\Theta_2|\rangle$, and (c) their fluctuations $\delta\Theta_n = \sqrt{\langle\Theta_n^2\rangle - \langle|\Theta_n|\rangle^2}$ as a function of time (in units of $1/\kappa$). The parameters and quench protocol are the same as in figure 7; the curves are evaluated for different numbers of atoms and of trajectories. In (a) and (b) the data correspond to $N = 25, 50, 100, 200$ particles (see legends for color code) and respectively 1000, 500, 250, 125 trajectories. The horizontal dashed lines indicate the values predicted by the global minimum of the free energy in equation (16). The finite values of the order parameters at $t = 0$ are due to finite size effects, $\langle|\Theta_n(0)|\rangle = 1/\sqrt{\pi N}$. The curves in (c) are calculated for $N = 100$ and 250 trajectories.</figcaption></figure>

contrast, upon increasing the ramp duration towards slower ramps this transient regime disappears. In particular, for the slowest ramp considered here, dissipation leads to quasi-adiabatic dynamics. Figure 13 shows order parameters $\langle|\Theta_1(t)|\rangle$ and $\langle|\Theta_2(t)|\rangle$ at $t = 3.77 \times 10^6/\kappa$, where the curves of figure 11 reach an asymptotic behavior. Self-organization in the long-wavelength grating depends on the ramp duration $\tau$ and is found for $\tau > 10^4/\kappa$. Note that short-wavelength order quantified by $\langle|\Theta_2(t)|\rangle$ only slightly depends on the ramp duration.

On a microscopic scale, it seems that the reason for better long-wavelength ordering after slower ramps is that more time is spent close to the transition line ($\alpha_1 = \alpha_2 \sim 1$), where the local minimum of the free energy is not deep enough to stably trap the system. In order to test this conjecture, we consider a two-step quench protocol which splits the sudden quench of section 3.2 into two subsequent quenches. One occurs at $t = 0$ from a paramagnetic to a ferromagnetic bistable phase, but close to the transition line: $\alpha_{1int} = \alpha_{2int} = 1.1$. This quench shows a vanishing value of $P_i(\Theta_2 < 0)$ for sufficiently long times as in figure 10. The second sudden quench occurs after an elapsed time $\tau$ and goes from this intermediate point to $\alpha_{1f} = \alpha_{2f} = 2$. The detuning $\Delta_c$ is kept constant during the evolution.
<figure><img src="image_15.png" /></figure>

<figure><img src="image_16.png" /><figcaption>**Figure 9.** Dynamics of (a) the single-particle kinetic energy $\langle p^2/2m \rangle$ (in units of $\hbar\omega_r$) and (b) the kurtosis $K = \langle p^4 \rangle / \langle p^2 \rangle^2$ for $N = 25, 50, 100, 200$ particles (see legends) and respectively 1000, 500, 250, 125 trajectories. The other parameters and initial conditions are the same as in figure 8. The horizontal dashed line in (a) indicates the asymptotic value predicted by equation (13).</figcaption></figure>

<figure><img src="image_17.png" /><figcaption>**Figure 10.** Asymptotic values of $\langle |\Theta_1(t_f)| \rangle$, $\langle |\Theta_2(t_f)| \rangle$, and $P_{t_f}(\Theta_2 < 0)$, equation (20), as a function of $\alpha_{1f}$. The quenches start from the same initial values in the paramagnetic phase ($\alpha_{1i}, \alpha_{2i} \ll 1$ and $\Delta_c = -\kappa$) and end up with different values $\alpha_{1f}, \alpha_{2f}$ with $\alpha_{1f} = \alpha_{2f}$ (lying along the right green line in figure 2(a)) and $\Delta_c = -\kappa$. The circles correspond to the results of the numerical simulations at $t_f = 3.77 \times 10^6/\kappa$ with $N = 100$ particles and 250 trajectories. The dashed lines indicate the predictions of the global minima of equation (16).</figcaption></figure>

Figure 14 shows the time evolution of the mean absolute values of the order parameters for different time intervals $\tau$ between the two quenches. The order parameters undergo an initial violent relaxation at $t = 0$, when the first sudden quench occurs, and a second one immediately after the second quench (which looks like a jump in logarithmic scale). As expected, the longer the time between the two quenches, the closer the asymptotic value to that of the global minimum. Inspecting the dynamics of the kinetic energy in figure 15 we observe that for large $\tau$ the atoms are cooled into the stationary state at $\alpha_n \sim 1$. At this point of the phase diagram the free energy has two ferromagnetic global minima, while the nematic local minimum is very shallow. The system thus gets cooled close to the global minima of the free energy at $\alpha_n = 2$, and remains trapped there after the second quench.

Figure 16 shows the mean absolute values of the order parameters, as extracted from the numerical data at $t = 3.77 \times 10^6/\kappa$, as a function of the time between the two quenches. These values are compared to the
<figure><img src="image_18.png" /></figure>

<figure><img src="image_19.png" /><figcaption>**Figure 11.** Mean value of order parameters (a) $\langle|\Theta_1|\rangle$ and (b) $\langle|\Theta_2|\rangle$ as a function of time (in units of $1/\kappa$) for $N = 100$ and $\Delta_c = -\kappa$, evaluated numerically with 250 trajectories. The curves are the time evolution during and after linear ramps of duration $\tau = 0$ (blue), $\tau = 5.5 \times 10^2/\kappa$ (red), $\tau = 6.8 \times 10^3/\kappa$ (yellow), $\tau = 8.5 \times 10^4/\kappa$ (purple), and $\tau = 2 \times 10^6/\kappa$ (green). The ramps occur from the paramagnetic to the bistable phase, specifically from $\alpha_{1i} = \alpha_{2i} = \varepsilon \ll 1$ to $\alpha_{1f} = \alpha_{2f} = 2$. As before, at $t = 0$ the initial state of the atoms is the steady state, equation (8), for $\alpha_n = \alpha_m$ and $\Delta_c = -\kappa$. The dashed horizontal lines show the steady-state values of the global minima of the free energy, equation (16).</figcaption></figure>

<figure><img src="image_20.png" /><figcaption>**Figure 12.** Mean value of the single-particle kinetic energy $\langle p^2/(2m) \rangle$ (in units of $E_{\text{kin},0} = \hbar\kappa/4$) as a function of time (in units of $1/\kappa$) for the same parameters and color codes as in figure 11.</figcaption></figure>

predictions of the global minimum of the free energy at $\alpha_1 = \alpha_2 = 2$ and $\Delta_c = -\kappa$. The behavior is quite similar to that observed when performing a linear ramp of corresponding duration (figure 13). Dynamical ordering in the long-wavelength mode thus seems to require that the atoms are initially cooled close to the global minima. This is realized by means of sufficiently long time $\tau$ spent close to the transition point.

4. Cooling into organized structures

We now analyze sudden quenches of the parameter $\alpha_n$ starting with different initial single-particle momentum widths. A possible realization is a quench in the detuning since $\Delta_c$ controls the steady-state temperature (see equation (13)). Using these we consider quenches which could lead to either heating or cooling of the system to the stationary temperature $T_0$,
<figure><img src="image_21.png" /><figcaption>Figure 13. Values of $\langle |\Theta_1(t_f)| \rangle$ (blue) and $\langle |\Theta_2(t_f)| \rangle$ (red) at $t_f = 3.77 \times 10^6/\kappa$ and as a function of the ramp duration $\tau$ (in units of $1/\kappa$) for the same parameters as in figure 11. The dashed horizontal lines show the steady-state value predicted by the free energy, equation (16).</figcaption></figure>

<figure><img src="image_22.png" /></figure>

<figure><img src="image_23.png" /><figcaption>Figure 14. Dynamics of $\langle |\Theta_1| \rangle$ and $\langle |\Theta_2| \rangle$ as a function of time (in units of $1/\kappa$) for the two-step quench protocol. Here, the parameters $\alpha_n$ suddenly ramp at $t = 0$ from the initial values $\alpha_{1i}, \alpha_{2i} \ll 1$ to $\alpha_{1int} = \alpha_{2int} = 1.1$; after a time interval $\tau$, there is a second quench from $\alpha_{1int} = \alpha_{2int} = 1.1$ to $\alpha_{1f} = \alpha_{2f} = 2$. The parameters are $\Delta_c = -\kappa$, $N = 100$ with 250 trajectories, and $\tau = 1/\kappa$ (blue), $10/\kappa$ (red), $10^2/\kappa$ (yellow), $10^3/\kappa$ (purple), $10^4/\kappa$ (green), $10^5/\kappa$ (light blue) and $10^6/\kappa$ (dark red). The dashed horizontal lines show the values predicted by the global minimum of the free energy, equation (16), at $\alpha_1 = \alpha_2 = 2$ and $\Delta_c = -\kappa$.</figcaption></figure>

$$k_B T_0 = \frac{\hbar \kappa}{2}, \qquad (21)$$

namely, the minimal temperature achieved by cavity cooling, which is equivalent to setting $\Delta_c = -\kappa$. Therefore, we also consider initial thermal distributions which are spatially uniform and have a temperature $T_{ini} < T_0$. The initial momentum distributions we consider are Gaussian and their width is $\Delta p^2 = mk_B T_{ini}$.

Figure 17 shows the time evolution of the mean absolute values of the order parameters for different values of $T_{ini}$ ranging from $0.1 T_0$ up to $5 T_0$. The asymptotic value of $\langle |\Theta_1| \rangle$ increases with the initial temperature: The hotter the system is initially, the smaller the fraction of trajectories which remain trapped in the metastable, nematic state is. The corresponding time evolution of the mean kinetic energy per particle is displayed in figure 18 and it shows that for $T_{ini} = 2 T_0$ (and even more for $T_{ini} = 5 T_0$) the system stays relatively hot over time
<figure><img src="image_24.png" /><figcaption>**Figure 15.** Mean kinetic energy per particle ($\langle p^2/2m \rangle$) in units of the asymptotic value $E_{\text{kin,0}} = \hbar\kappa/4$) as a function of time (in units of $1/\kappa$) for the same parameters and color codes as in figure 14.</figcaption></figure>

<figure><img src="image_25.png" /><figcaption>**Figure 16.** The symbols correspond to $\langle|\Theta_1(t_f)|\rangle$ (blue) and $\langle|\Theta_2(t_f)|\rangle$ (red) at $t_f = 3.77 \times 10^6/\kappa$ as a function of the time between the two quenches, $\tau$ (in units of $1/\kappa$), for the same parameters as in figure 14. The dashed horizontal lines show the steady-state values predicted by the global minimum of the free energy, equation (16), at $\alpha_1 = \alpha_2 = 2$ and $\Delta_c = -\kappa$.</figcaption></figure>

scales of the order of $10^4/\kappa$. For lower initial temperatures, the system is instead heated by the energy released by the sudden quench before relaxation cools the atoms.

As shown in figure 17, for initially cold samples a long-wavelength Bragg grating is formed faster than for hotter samples. In this case we recognize a three-stage dynamics like the one observed for sudden quenches of the laser intensity, when a transient long-range order is established for times $t > 10/\kappa$ and $t < 10^3/\kappa$. For $t > 10^3/\kappa$ dissipation becomes significant and $\langle|\Theta_1|\rangle$ increases to a stationary value. This relaxation stage is also present for samples with initial temperatures higher than $T_0$; however, in this hotter case it is significantly faster. Taking a threshold value $\langle|\Theta_1|\rangle|_{\text{thres}} = 0.5$, we observe that buildup of long-wavelength order can be up to a hundred times shorter than for a cold initial state. This is reminiscent of the Mpemba effect in supercooled water [24–28]. Its origin could be traced to a suppression of long-wavelength order if short-wavelength order is already established on a much shorter time scale, as shown in figure 17(b).

In figure 17(a) we observe that the final value of $\langle|\Theta_1|\rangle$ does not coincide with its predicted stationary value even after very long cooling times. This can also be seen in figure 19, which shows the mean absolute value of the order parameters at $t = 3.77 \times 10^6/\kappa$ as a function of the initial temperature for $N = 100, 200$. One would expect that $\langle|\Theta_1|\rangle$ should have reached a constant value corresponding to the stationary state. Apparently this is not the case and even for finite $N$ a significant fraction of trajectories converge to and remain in the local minimum. This behavior gets much less pronounced when the initial temperature lies above a certain threshold set by the energy released by the quench itself.

## 5. Comparison of numerical approaches

The discussion in this paper is based on results obtained by numerical integration of stochastic differential equations (1) and on their comparison with the corresponding analytical model. Both rely on the validity of the so-called bad cavity limit, where cavity damping is the shortest time scale, and particularly on treating retardation as a small parameter in the dynamics. This regime allows one to systematically describe quantum fluctuations of the cavity degrees of freedom by eliminating the cavity variables from the equations of motion of
<figure><img src="image_26.png" /><figcaption>Figure 17. Dynamics of (a) $\langle|\Theta_1|\rangle$ and (b) $\langle|\Theta_2|\rangle$ as a function of time (in units of $1/\kappa$) after a sudden quench at $t = 0$ from $\alpha_{1i} = \alpha_{2i} = 0$ and $T_{ini}$ to $\alpha_{1f} = \alpha_{2f} = 2$ and $\Delta_c = -\kappa$ (corresponding to the asymptotic temperature $T_0$) for $N = 100$ and 250 trajectories. The different curves correspond to $T_{ini} = 5T_0$ (blue), $2T_0$ (red), $T_0$ (yellow), $0.5T_0$ (purple) and $0.1T_0$ (green). The dashed horizontal lines show the steady-state values predicted by the global minimum of the free energy, equation (16).</figcaption></figure>

<figure><img src="image_27.png" /><figcaption>Figure 18. Mean kinetic energy per particle $\langle p^2/2m \rangle$ (in units of $E_{kin,0} = \hbar\kappa/4$) as a function of time (in units of $1/\kappa$) for the same parameters and color codes as in figure 17. The horizontal dashed line corresponds to the asymptotic value $\langle p^2/2m \rangle = E_{kin,0}$.</figcaption></figure>

the external degrees of freedom. We now compare these predictions with those of the stochastic differential equations derived in [15], where the cavity degrees of freedom are treated in the semi-classical limit but included at all orders of retardation expansion. These stochastic differential equations are here extended to our setup composed of two cavity modes [16]:

$$dx_j = \frac{p_j}{m}dt, \qquad (22a)$$

$$dp_j = \sum_{n=1,2} 2\hbar nkS_n \mathcal{E}_{n,r} \sin(nkx_j) dt, \qquad (22b)$$

$$d\mathcal{E}_{n,r} = (-\Delta_n \mathcal{E}_{n,i} - \kappa_n \mathcal{E}_{n,r}) dt + d\xi_{n,r}, \qquad (22c)$$

$$d\mathcal{E}_{n,i} = (\Delta_n \mathcal{E}_{n,r} - \kappa_n \mathcal{E}_{n,i} - NS_n \Theta_n) dt + d\xi_{n,i}, \qquad (22d)$$

where $\mathcal{E}_{n,r} = \text{Re}\{\mathcal{E}_n\}$ and $\mathcal{E}_{n,i} = \text{Im}\{\mathcal{E}_n\}$ are the real and imaginary parts of the positive-frequency component of the cavity field mode $n = 1, 2$. The Wiener processes $d\xi_{n,i}$, $d\xi_{n,r}$ have a vanishing first moment,
<figure><img src="image_28.png" /><figcaption>Figure 19. The symbols correspond to the values of $\langle|\Theta_1(t_f)|\rangle$ (blue, yellow) and $\langle|\Theta_2(t_f)|\rangle$ (red, purple) at $t_f = 3.77 \times 10^6/\kappa$ as a function of $T_{ini}$ (in units of $T_0$) for the same parameters as in figure 14 but for $N = 200$ particles (125 trajectories). See inset for the color code. The dashed horizontal lines show the steady-state values predicted by the global minimum of the free energy, equation (16).</figcaption></figure>

$\langle d\xi_{n,i} \rangle = 0 = \langle d\xi_{n,r} \rangle$, while the second moments fulfill $\langle d\xi_{n,i}d\xi_{m,i} \rangle = \delta_{nm}\kappa/2dt$, $\langle d\xi_{n,r}d\xi_{m,r} \rangle = \delta_{nm}\kappa/2dt$, and $\langle d\xi_{n,r}d\xi_{m,i} \rangle = 0$.

The results of the simulations based on the two approaches for a single-mode cavity show good agreement. For the two-mode cavity we generally find qualitative agreement. Quantitative discrepancies are found in general for the momentum distribution: The simulations based on equation (22) predict for certain parameters samples whose temperatures are 10% hotter than the ones obtained with equation (1). Small differences are also found for the order parameters after the quenches into the bistable phase.

Figure 20 shows a representative result of the discrepancies found after the quench protocol discussed in section 3.2. The two simulations predict different stationary values for both the kinetic energy and the order parameters. We believe that this discrepancy is due to retardation effects, which are neglected in the approach of equation (1) and become relevant when the atoms are trapped at tight minima.

In order to test our conjecture we use the prediction of the kinetic theory of [29, 30], where the temperature of the stationary thermal distribution was corrected by the contribution from the atoms’ localization at the minima of the self-organized lattice,

$$k_B \tilde{T} = \hbar \frac{\Delta_c^2 + \kappa^2}{4|\Delta_c|} + \hbar \frac{\omega_0^2}{|\Delta_c|} \quad (23)$$

Here, $\omega_0$ is the frequency of oscillation about the lattice minima in the harmonic approximation. It can be estimated by using equation (22b) and imposing the equality

$$d p_j \approx \sum_{n=1,2} 2\hbar(nk)^2 S_n \varepsilon_{n,r} x_j dt \equiv -m\omega_0^2 x_j dt.$$

This delivers an analytic estimate of the frequency

$$\omega_0^2 = \omega_r \frac{\Delta_c^2 + \kappa^2}{-\Delta_c}(\alpha_1 \Theta_1 + 4\alpha_2 \Theta_2),$$

where we use equation (4). For the parameters of the quench in figure 20, with $\Delta_c = -\kappa$ and $\alpha_1 = \alpha_2 = 2$, we obtain $k_B \tilde{T} \approx 1.1k_B T_0$, where $T_0$ is the temperature given in equation (21). Indeed, this corrected value of the final temperature is in good agreement with the discrepancy observed in figure 20(a).

This hypothesis is also consistent with the discrepancy observed in the asymptotic values of the order parameters. In fact, the stationary temperature and the final values of the order parameters are related: the stationary values of the order parameters are determined by the parameters $\alpha_1$, $\alpha_2$ [13] and thus depend on both field intensities and detunings (see equation (4)). According to this hypothesis, the asymptotic values of the order parameters for the simulation using equation (22) should be the ones corresponding to the system’s parameters with the corrected temperature $\tilde{T}$; hence we shall minimize the free energy of equation (16) using $\tilde{\beta} = 1/(k_B \tilde{T})$, equation (23), instead of $1/(k_B T_0)$. This is equivalent to rescaling the phase diagram in figure 2(a) using the prescription $\tilde{\alpha}_n = \alpha_n T_0 / \tilde{T} < \alpha_n$, and results in a smaller stationary value of the order parameter, which is consistent with the discrepancies visible in figures 20(b) and (c).
<figure><img src="image_29.png" /></figure>

<figure><img src="image_30.png" /></figure>

<figure><img src="image_31.png" /></figure>

**Figure 20.** Dynamics of (a) single-particle kinetic energy (in units of $E_{\text{kin},0}$), (b) $\langle|\Theta_1|\rangle$, and (c) $\langle|\Theta_2|\rangle$ as a function of time (in units of $1/\kappa$) following a quench at $t=0$ from $\alpha_{1i} = \alpha_{2i} \ll 1$ to $\alpha_1 = \alpha_2 = 2$ for $\Delta_c = -\kappa$ and $N=100$. The blue (red) lines correspond to simulations using equation (1) (equation (22)). The black dashed lines denote the values of the order parameters obtained by the free energy, equation (16). In (a) the blue (red) line corresponds to 250 (500) trajectories. In (b) and (c) the blue and red lines correspond to 1000 (500) trajectories. Note that a quench from $\alpha_1, \alpha_2$ as performed in section 3.2 for fixed $\Delta_c = -\kappa$ corresponds for the simulation of equation (22) to a quench in the pumping strength $S_n$ such that $NS_n^2 = \alpha_n \kappa^2$. At $t_f = 3.77 \times 10^6/\kappa$ we observe that 16.9% (15.4%) of the trajectories are in the nematic phase using equation (1) (equation (22)).

## 6. Conclusions

In this work we have studied the semi-classical dynamics of atoms interacting with two cavity modes after quenches of the intensity and/or frequency of the pumping lasers. In the quench protocols the laser parameters were varied across transition lines separating a disordered phase from an ordered self-organized phase. We could verify numerically that the states reached at the asymptotics of the dynamics correspond to the minima of the free energy of a corresponding thermodynamic description developed in [13]. This picture is further confirmed by the comparison with numerical simulations based on different initial assumptions. This analysis shows, in particular, that trapping of the system in the local minima of the free energy crucially depends on the initial temperature and on the cooling rate.

We observe, in addition, that the system can be trapped in metastable configurations for transient times which cannot be understood in terms of the effective thermodynamic description. For hundreds of particles the lifetime of these states is about four orders of magnitude longer than the cavity lifetime, and is expected to increase with $N$. They share analogies with metastable configurations found in the GHMF when performing quenches in the microcanonical ensemble [22]. Since the phase diagrams of the GHMF and the model here considered can be formally mapped onto each other [13], we conjecture that these metastable configurations
could be due to the coherent dynamics. This conjecture can be tested by means of a mean-field analysis such as
the one performed in [12] for a single-mode cavity.

Interestingly, when the initial temperature of the atomic ensemble is different from the stationary temperature of cavity cooling, we observe that the final magnitude of asymptotic order changes. In particular when the initial temperature is even lower than the predicted cavity cooling temperature, the probability that the system remains trapped in metastable configurations is further increased. This is similar to the behavior of supercooled water [24–28].

Here we consider the very special case of two commensurate modes. While this already highlights many generic properties of the dynamics, future considerations certainly should include the case in which the wavelengths of the cavity modes are incommensurate [31], so that the ordering mechanisms are much more intensely competing and a multitude of metastable states can form. A further interesting direction is operation with much colder temperatures or in the side-band resolved cooling regime [32]. Here it is intriguing to consider in which form metastable states survive deep in the quantum regime. Besides diffusion they could be depleted via tunneling and atom–field entanglement plays an important role in this dynamics [33], a process which should also be relevant in closely related schemes of simulated quantum annealing [34].

## Acknowledgments

The authors thank Tobias Donner, Sebastian Krämer, and Francesco Rosati for stimulating and helpful discussions. This work was supported by the German Research Foundation (DFG, DACH project ‘Quantum crystals of matter and light’) and by the European Commission (ITN ‘ColOpt’). V T and H R are supported by Austrian Science Fund Project No. I1697-N27. T K and V T contributed equally to this work.

## ORCID iDs

Helmut Ritsch https://orcid.org/0000-0001-7013-5208

## References

[1] Horak P, Hechenblaikner G, Gheri K M, Stecher H and Ritsch H 1997 *Phys. Rev. Lett.* **79** 4974

[2] Domokos P and Ritsch H 2003 *J. Opt. Soc. Am. B* **20** 1098

[3] Black A T, Chan H W and Vuletić V 2003 *Phys. Rev. Lett.* **91** 203001

[4] Baumann K, Guerlin C, Brennecke F and Esslinger T 2010 *Nature* **464** 1301

[5] Asbóth J K, Domokos P, Ritsch H and Vukics A 2005 *Phys. Rev. A* **72** 053417

[6] Schütz S and Morigi G 2014 *Phys. Rev. Lett.* **113** 203002

[7] Schütz S, Jäger S B and Morigi G 2015 *Phys. Rev. A* **92** 063808

[8] O'Dell D H J, Giovanazzi S and Kurizki G 2003 *Phys. Rev. Lett.* **90** 110402

[9] Münstermann P, Fischer T, Maunz P, Pinkse P W H and Rempe G 2000 *Phys. Rev. Lett.* **84** 4068

[10] Vuletić V and Chu S 2000 *Phys. Rev. Lett.* **84** 3787

[11] Schütz S, Jäger S B and Morigi G 2016 *Phys. Rev. Lett.* **117** 083001

[12] Jäger S B, Schütz S and Morigi G 2016 *Phys. Rev. A* **94** 023807

[13] Keller T, Jäger S B and Morigi G 2017 *J. Stat. Mech.* **6** 064002

[14] Schütz S, Habibian H and Morigi G 2013 *Phys. Rev. A* **88** 033427

[15] Domokos P, Horak P and Ritsch H 2001 *J. Phys. B* **34** 187

[16] Torggler V and Ritsch H 2014 *Optica* **1** 336

[17] Brennecke F, Donner T, Ritter S, Bourdel T, Köhl M and Esslinger T 2007 *Nature* **450** 268

[18] Landig R, Hruby L, Dogra N, Landini M, Mottl R, Donner T and Esslinger T 2016 *Nature* **532** 476

[19] Léonard J, Morales A, Zupancic P, Esslinger T and Donner T 2017 *Nature* **543** 87

[20] Dalibard J and Cohen-Tannoudji C 1985 *J. Phys. B* **18** 1661

[21] Campa A, Dauxois T and Ruffo S 2009 *Phys. Rep.* **480** 57

[22] Teles T N, Benetti F P d C, Pakter R and Levin Y 2012 *Phys. Rev. Lett.* **109** 230601

[23] Pikovsky A, Gupta S, Teles T N, Benetti F P d C, Pakter R, Levin Y and Ruffo S 2014 *Phys. Rev. E* **90** 062141

[24] Auerbach D 1995 *Am. J. Phys.* **63** 882

[25] Brownridge J D 2011 *Am. J. Phys.* **79** 78

[26] Tao Y, Zou W, Jia J, Li W and Cremer D 2017 *J. Chem. Theory Comput.* **13** 55

[27] Jin J and Goddard W A III 2015 *J. Phys. Chem. C* **119** 2622

[28] Zhang X *et al* 2014 *Phys. Chem. Chem. Phys.* **16** 22995–3002

[29] Niedenzu W, Grießer T and Ritsch H 2011 *Europhys. Lett.* **96** 43001

[30] Grießer T, Niedenzu W and Ritsch H 2012 *New J. Phys.* **14** 053031

[31] Habibian H, Winter A, Paganelli S, Rieger H and Morigi G 2013 *Phys. Rev. Lett.* **110** 075304

[32] Krämer S and Ritsch H 2014 *Phys. Rev. A* **90** 033833

[33] Maschler C, Ritsch H, Vukics A and Domokos P 2007 *Opt. Commun.* **273** 446

[34] Torggler V, Krämer S and Ritsch H 2017 *Phys. Rev. A* **95** 032310