Quenches across the self-organization transition in multimode cavities

# Quenches across the self-organization transition in multimode cavities

## 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 but below the cavity resonance. When sufficiently illuminated, moreover, the atoms’ steady state undergoes a phase transition from homogeneous density to crystalline order. 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, their intensity is either suddenly switched or slowly ramped across the self-ordering transition. Numerical simulations for different ramp protocols predict that the system exhibits long-lived metastable states, whose occurrence strongly depends on initial temperature, ramp speed, and number of atoms.

## I 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 Horak et al. (1997); Domokos and Ritsch (2003); Black et al. (2003); Baumann et al. (2010). 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 Asbóth et al. (2005); Baumann et al. (2010). 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 Asbóth et al. (2005); Schütz and Morigi (2014); Schütz et al. (2015); O’Dell et al. (2003); Münstermann et al. (2000). This force has also 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 Vuletić and Chu (2000); Black et al. (2003). Theoretical studies with single-mode resonators predicted that this dissipation can establish long-range correlations and support the onset of metastable ordered structures Schütz et al. (2016); Jäger et al. (2016).

In a multimode cavity and for several illumination frequencies, competing ordering processes are present and lead to a richer phase dynamics. In a two-mode cavity, like the one depicted in Fig. 1(a), the transition to self-organization can be a phase-transition of first or second order depending on the laser intensities and on their relative strength Keller et al. (2017). The corresponding self-ordered phases can exhibit superradiant scattering either in one or in both cavity modes, as illustrated in Fig. 1(b), while the asymptotic distribution of the atoms can be thermal provided that the lasers’ frequencies are suitably chosen Keller et al. (2017). In our example the particles can order in a lattice at a given length scale and/or on a lattice with half the period . 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 Ref. Schütz et al. (2013) for a similar system. We find that even at very long times the atoms’ spatial distribution strongly depends on the initial temperature, ramp speeds, and on the 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 when compared to higher initial temperatures.

Our work is organized as follows: In Section II we introduce the system and the semi-classical equations describing the dynamics. The atoms’ stationary properties are then summarized in a phase diagram, which was derived in Ref. Keller et al. (2017). In Sec. III we numerically study the real time dynamics when the parameters are varied within the phase diagram according to different quench protocols. In Sec. IV we analyze the crystallization dynamics starting from spatially homogeneous distributions with different momentum widths. In Sec. V 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 Refs. Domokos et al. (2001); Torggler and Ritsch (2014). The conclusions are drawn and future perspectives are discussed in Sec. VI.

## Ii Semiclassical dynamics

The system we consider consists of a gas of cold atoms with mass , 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 axis) by a tight external dipole trap Brennecke et al. (2007); Landig et al. (2016) and is here described in the semi-classical limit.

The geometry of the setup is illustrated in Fig. 1. Lasers with (rescaled) intensities propagate in a direction orthogonal to the cavity axis and are quasi resonant with the standing wave cavity modes with frequency and wave number () foo (); Léonard et al. (2017). The lasers have frequency and linear polarization which is parallel to the one of the corresponding cavity mode. Each pair of laser and cavity mode couples to an atomic dipolar transition at frequency , where and are the laser and the vacuum Rabi frequency, respectively. Spontaneous scattering processes are suppressed when the absolute value of the detuning exceeds the coupling strengths and the detuning between laser and cavity mode by orders of magnitude: . The relevant dissipative process are given by cavity decay, and we denote by the loss rate of cavity mode .

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 a 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 which, at all instants of time, is orders of magnitude larger than the photon recoil : Dalibard and Cohen-Tannoudji (1985); Schütz et al. (2013); Keller et al. (2017). The corresponding stochastic differential equations read

 dxj =pjmdt, (1a) dpj =(Fj,ad+Fj,ret)dt+dW(1)j+dW(2)j, (1b)

where

 Fj,ad= −∑n=1,22ℏnkαnℏβnsin(nkxj)Θn, (2) Fj,ret= −∑n=1,2ℏ(nk)2mαnκn−Δnsin(nkxj)1NN∑l=1plsin(nkxl), (3)

and

 αn= 4NS2nΔ2n(Δ2n+κ2n)2, (4) βn= −4Δnℏ(Δ2n+κ2n). (5)

Here, is the amplitude of coherent scattering by a single atom and has the dimension of a frequency, while and in Eq. (1b) describe Wiener processes, which fulfill and ( and ). Here,

 Dnij=(ℏnk)2αnℏβnκn−Δnsin(nkxi)sin(nkxj). (6)

Finally, the parameter

 Θn=1NN∑i=1cos(nkxi) (7)

quantifies Bragg ordering of the atoms in the cavity mode with wave number . In particular, when the atoms are localized either at the maxima or at the minima of , which is the configuration which maximizes the intracavity field intensity. We identify with the order parameter for self-organization in the corresponding cavity mode Keller et al. (2017). Below, we will denote by ”long-wavelength order” a configuration with non-vanishing value of , corresponding to a Bragg grating with period . Similarly, ”short-wavelength order” refers to a configuration with , corresponding to a Bragg grating with . 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 . For details we refer to Ref. Keller et al. (2017).

### ii.1 Stationary states

An analysis of the Fokker-Planck equation at the basis of Eq. (1) allows to identify the conditions for the existence of a stationary state. The latter exists provided that and , see Eq. (5). In this case the atoms’ distribution at steady state reads Keller et al. (2017)

 fst(x1,p1,…,xN,pN)=exp(−βHeff)Z(β) (8)

where is the effective Hamiltonian derived after eliminating the cavity field variables,

 Heff=N∑j=1p2j2m−∑n=1,2NαnβnΘ2n, (9)

while denotes the partition function:

 Z(β)=1ΔN∫∞−∞dp1...∫∞−∞dpN∫λ0dx1...∫λ0dxNe−βHeff, (10)

with and with the single particle unit phase space volume. In the following we will assume that the cavity decay rates are equal,

 κ1=κ2=:κ, (11)

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

 Δ1=Δ2=:Δc<0. (12)

The phase diagram of the system can be determined by using that the steady state, Eq. (8), has the form of a thermal state. On the basis of this observation we introduce the temperature of the stationary state, which is defined as

 kBT=β−1=ℏ(Δ2c+κ2)−4Δc, (13)

with the Boltzman constant . The steady-state temperature has the same functional dependence on and as for a single-mode cavity Schütz et al. (2013, 2015). We can further define the free energy per particle using the formal equivalence with the canonical ensemble of equilibrium statistical mechanics Schütz et al. (2015):

 F=−1Nβln(Z(β)). (14)

Following the procedure detailed in Refs. Schütz et al. (2015); Keller et al. (2017); Jäger et al. (2016) we determine the global minima of in an appropriately defined thermodynamic limit, which consists in keeping constant for . The global minima are the resulting stationary phases. The corresponding order parameters and , in particular, are determined by and . When the fields are sufficiently weak, then the density is homogeneous and there is no structural order. We denote this phase by paramagnetic, borrowing the notation of the generalized Hamiltonian mean-field model (GHMF) Campa et al. (2009); Teles et al. (2012); Pikovsky et al. (2014) to which this model can be mapped. The possible ordered phases at steady state are illustrated in Fig. 1(b) and take one of four set of values. In particular, the ferromagnetic phase is characterized by (i) , and (ii) , , exhibiting Bragg order in both cavity modes. The nematic phases (iii) and (iv), instead, are characterized by no order in the long-wavelength mode, while can be either negative or positive.

The resulting phase diagram in the plane is shown in Fig. LABEL:fig:phase_diagram and reproduces the one of Ref. Keller et al. (2017). The phases are separated by either first- or second-order transitions, which have been determined using Ehrenfest’s criterion Pikovsky et al. (2014). 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) phase. Examples of the free-energy landscape in the plane are shown in subplots (b) and (c). Subplot (b) corresponds to the parameters of the red bullet labeled by (b) in subplot (a): Here, the free energy exhibits two symmetric global minima which correspond to the ferromagnetic phase. In subplot (c), corresponding to the parameters of the red bullet labeled by (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 . We denote this region by bistable referring to the existence of a second, metastable state in which the system can be dynamically trapped.

## Iii Dynamics of self-organization

We now examine the dynamics of the system when the values of and are varied as a function of time. Experimentally, this corresponds to vary the pump laser intensities or their detuning with respect to the cavity mode frequencies. At time we assume that the system is prepared in the stationary state of a paramagnetic phase, described by the distribution in Eq. (8) by setting in Eq. (9) (). The values appearing in the equations of motion (1) are then varied in time, by performing either (i) a sudden quench, i.e. suddenly switching the two values of and , or (ii) a slow quench, consisting in varying monotonously and continuously in time towards the final values and . We choose the final values in the ferromagnetic phase. The quench protocols we consider are illustrated by the green lines in Fig. LABEL:fig:phase_diagram(a): for sudden quenches, the initial and final values are two points connected by the 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.

In what follows we perform numerical simulations of Eq. (1) using the parameters of a gas of Rb atoms. In particular, we take with the wavelength of the line. The corresponding recoil frequency is . The cavity linewidth is taken to be , so that . A possible realization of the two-mode setup here considered has been discussed in Refs. Keller et al. (2017); Léonard et al. (2017).

### iii.1 Sudden quench into the ferromagnetic phase

We first consider sudden quenches from in the paramagnetic phase to in the ferromagnetic phase, keeping . 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, Eq. (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 displays the distribution for the order parameters and as a function of time for . It is defined as a time sequence of normalized histograms

 Pt(Θ)=\# trajectories with Θ(t)∈[Θ−ΔΘ/2,Θ+ΔΘ/2]\# trajectories×ΔΘ, (15)

where is calculated on each trajectory of the simulations with the stochastic different equations and its value is determined according to the precision of the grid in . We observe that at a given time scale of the order of , 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 Asbóth et al. (2005). The order parameter of the short-wavelength mode , which is weakly pumped, substantially grows to a positive value long after the symmetry breaking. The fact that vanishes for negative values comes from the ordering of the atoms close to the anti-nodes of the dominant long-wavelength mode field (see Fig. 1).

The distributions at the asymptotics are reported in the right panels of Fig. 3. They are obtained by averaging over times , where a stable configuration has been reached. Formally

 P(Θ)=Nt∑i=1Pti(Θ)/Nt, (16)

where is the number of instants of times at which the distribution is sampled in the interval , with and . Comparing the widths of the distributions in the right panel of Fig. 3 one observes that after sufficiently long times the long-wavelength order parameter fluctuates less than the short-wavelength order parameter.

Figures LABEL:fig:ferro_quench:Theta(a)-(b) display the dynamics of the mean absolute value of the order parameters for different values of . Figure LABEL:fig:ferro_quench:Theta (c) shows the time evolution of the fluctuations of the order parameters for particles. The order parameters asymptotically tend to the values predicted by the free energy, indicated by the horizontal dashed line, for a time scale of the order of . Meanwhile the fluctuation relaxes to a much smaller value than the asymptotic value of reproducing the widths of the distributions in the right panel in Fig. 3. The time evolution of , in particular, is reminiscent of the one observed for quenches into the ferromagnetic phase in a single-mode resonator Schütz et al. (2016). It can be separated into three stages which we denote by (in order of their temporal appearance) (i) violent relaxation, corresponding to an exponential increase of the absolute value of the order parameter ; (ii) transient dynamics, corresponding to power-law scaling with time, and (iii) relaxation phase, where the mean values tend exponentially towards the asymptotic value. The violent relaxation can be described by a mean-field model Jäger et al. (2016), in the transient stage the coherent dynamics is prevailing, while the relaxation stage is dominated by dissipation Schütz et al. (2016). The transient and relaxation stages are characterized by time scales which increase with but with different functional dependence Jäger et al. (2016). The time scale can here be identified as the one at which the asymptotic state is reached for , 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 timescale. The metastable phase of the transient dynamics can be therefore denoted by ”nematic”, its lifetime increases with and for it is of the order of . 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 the quasi-stationary state due to the long-range coherent dynamics analogous to the ones reported in Ref. Teles et al. (2012). This conjecture is also supported by the behaviour of the single-particle kinetic energy and of the kurtosis , which are shown in Fig. LABEL:fig:ferro_quench:p. The latter quantifies the deviation of the momentum distribution from a Gaussian, for which it takes the value . 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 starting from the same initial values of but with different endpoints . We take different endpoints ranging from the paramagnetic to the ferromagnetic phase, under the constrain . The circles in Fig. 6 correspond to the numerical results for 100 particles at time , where we expect that the system has reached the steady state. These are in good agreement with the analytical results (dashed lines) based on evaluating the corresponding observables at the global minimum of the free energy. The interval where grows monotonically from to the value of the ferromagnetic phase is expected to shrink as 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 that is negative at :

 Pt(Θ2<0)=∫0−1dΘ2Pt(Θ2). (17)

We note that in the paramagnetic phase (homogeneous spatial distribution) we expect . In contrast, due to the given mode structure we expect that for long-wavelength ordering in the ferromagnetic phase. Indeed, as increases across the critical value, , quickly drops down to zero.

### iii.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 Fig. LABEL:fig:phase_diagram(a), which consists in equal effective pumping . In this parameter region (bistable phase) the free energy exhibits a local minimum, which is nematic. As in the previous case, the initial values are vanishingly small and the atoms are at the corresponding stationary distribution, whose temperature is determined by the detuning and whose spatial density is homogeneous. The quench is performed by switching the laser power keeping the detuning constant, thus the atoms should reach a thermal distribution with the same temperature as the initial state.

Figure 7 displays the time evolution of the trajectories’ -distribution for and . As opposed to the previous section, here a finite fraction of trajectories gets trapped in the nematic phase with vanishing value of and finite probability that takes negative values. This is visible in the small extra peaks in and (right panels). The trapping occurs at the time scale of the violent relaxation, and it seems stable over times of the order of . We conjecture that it persists also at asymptotic times. In Fig. LABEL:fig:mixed_quench:Theta the time evolution of the mean absolute value of the order parameters is shown for different numbers of particles. While reaches the same stationary value (in reality its value decreases slightly with ), instead the asymptotic value of decreases as 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 in subplot (c) reflects the contribution of these trajectories.

Mean single-particle kinetic energy and kurtosis are shown in Fig. LABEL:fig:mixed_quench:p. From their behaviour we infer that the metastable nematic state does not significantly deviate from a thermal distribution with the expected asymptotic temperature (Eq. (13)).

Peculiar features of these dynamics become visible when inspecting the probability at the asymptotics and as a function of in Fig. 10. As in Fig. 6, it vanishes identically upon leaving the paramagnetic phase, but increases again as are chosen deeper in the bistable phase of Fig. LABEL:fig:phase_diagram(a). Correspondingly, starts to decrease as increases, which suggests that from this point on the depth of the local minimum grows. The value of the order parameter at which starts to grow again identifies a threshold, above which the local minimum is sufficiently deep to stably trap particles.

### iii.3 Slow ramp into the bistable phase

We now consider linear ramps of across the transition region separating the paramagnetic from the bistable region. The ramps protocols have duration and sweep between the values , with . In particular, if , while for then is constant and equal to . Note that a sudden quench is the limit of a linear quench. We choose to vary the values of along the rightmost green line in Fig. LABEL:fig:phase_diagram(a), so that at all instants of time, with in the bistable phase. We further keep constant, and solely vary the pump intensity. This means that the asymptotic temperatures at each value of are equal.

Figure LABEL:fig:ramp:Theta displays the dynamics of the mean absolute value of the order parameters for for linear ramps with different durations . The dynamics following the sudden quench (cif. Fig. LABEL:fig:mixed_quench:Theta (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 and . This behaviour seems to be initiated at the instant of time when the parameters cross the critical point of the phase diagram. Moreover, for sufficiently slow ramps approaches the asymptotic value of the free energy’s global minimum, signaling stationary long-wavelength order.

We further note that for the orders parameters undergo a three-stage dynamics, as for the sudden quench (we attribute the fluctuations to the statistics of trajectories). For slower ramps, instead, the mean value of the order parameter tends exponentially towards the steady state, which approaches the free energy’s global minimum of Eq. (14) for . We believe that this behaviour is determined by the ramp duration with respect to the time scale of the transient dynamics, and thus by the time the parameters 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 Fig. 12, corresponding to the curves in Fig. LABEL:fig:ramp:Theta. For faster ramps it is similar to the sudden quench, exhibiting first a violent relaxation followed by a time interval where the dynamics is prevailingly coherent, and finally an exponential decay to the steady state value due to cavity cooling. In contrast, upon increasing the ramp duration towards slower ramps this transient regime disappears. Particularly, for the slowest ramp considered here, dissipation leads to quasi-adiabatic dynamics. Figure 13 shows the order parameters and at , where the curves of Fig. LABEL:fig:ramp:Theta have reached an asymptotic behaviour. Self-organization in the long-wavelength grating depends on the ramp duration and is found for . Note that short-wavelength order quantified by , in contrast, 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 (), 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 Sec. III.2 into two subsequent quenches: one at from a paramagnetic to a ferromagnetic bistable phase, but close to the transition line: . This quench shows a vanishing value of for sufficiently long times as in Fig. 10. The second sudden quench occurs after an elapsed time and goes from this intermediate point into . The detuning is kept constant during the evolution.

Figure LABEL:fig:jump_ramp:Theta displays the time evolution of the mean absolute values of the order parameters for different time intervals elapsed between the two quenches. The order parameters undergo a first violent relaxation at , 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 larger the time elapsed between the two quenches, the closer the asymptotic value is to the one of the global minimum. Inspecting the dynamics of the kinetic energy in Fig. 15 we observe that for large the atoms are cooled into the stationary state at . 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 , and remains trapped there after the second quench.

Figure 16 displays the mean absolute value of the order parameters, as extracted from the numerical data at , as a function of the time elapsed between the two quenches. Their value is compared to the predictions of the global minimum of the free energy at and . The behaviour is quite similar to the one observed when performing a linear ramp of corresponding duration, Fig. 13. Dynamical ordering in the long-wavelength mode seems thus to require that the atoms are initially cooled close to the global minima. This is realised by means of the sufficiently large time spent close to the transition point.

## Iv Cooling into crystalline order

We now analyse sudden quenches of the parameter starting with different initial single-particle momentum widths. A possible realization is a quench in the detuning since controls the steady state temperature, see Eq. (13). By these means we consider quenches which could either lead to heating or cooling of the system to the stationary temperature ,

 kBT0=ℏκ2, (18)

namely, the minimal temperature achieved by cavity cooling, corresponding to setting . Thereby we also consider initial thermal distributions which are spatially uniform and with temperature . The initial momentum distribution we consider are Gaussian and their width is .

Figure LABEL:fig:mpemba:Theta shows the time evolution of the mean absolute values of the order parameters for different values of ranging from up to . The asymptotic value of increases with the initial temperature: The hotter is initially the system, the smaller is the fraction of trajectories which remain trapped in the metastable, nematic state. The corresponding time evolution of mean kinetic energy per particle is displayed in Fig. 18 and shows that for (and even more for ) the system stays relatively hot over time scales of the order of . For lower initial temperatures, instead, the system is heated by the energy released by the sudden quench before relaxation cools the atoms.

As visible in Fig. LABEL:fig:mpemba:Theta, for samples initially cold 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 the sudden quenches of the laser intensity, when a transient long-range order is established for times and . For dissipation becomes important and increases to a stationary value. This relaxation stage is also present for samples with initial temperatures larger than , however, in this hotter case it is significant faster. Taking a threshold value , we observe that buildup of long-wavelength order can take up to a hundred times shorter than for a cold initial state. This is reminiscent of the Mpemba effect in supercooled water Auerbach (1995); Brownridge (2011); Tao et al. (2017); Jin and Goddard III (2015); Zhang et al. (2013). Its origin could be traced to a suppression of long-wavelength order if short-wavelength order is already established on a much faster time scale, visible in Fig. LABEL:fig:mpemba:Theta(b).

In Fig. LABEL:fig:mpemba:Theta (a) we observe that the final value of does not coincide with its predicted stationary value even after very long cooling times. This can also be seen in Fig. 19, which shows the mean absolute value of the order parameters at as a function of the initial temperature for . One would expect that should have reached a constant value corresponding to the stationary state. Apparently this is not the case and even for finite a significant fraction of trajectories converges to and remains in the local minimum. This behavior gets much less pronounced, if the initial temperature lies above a certain threshold set by the energy released by the quench itself.

## V Comparison between different numerical approaches

The discussion of this paper has been based on results obtained by numerical integration of the 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 fastest time scale, and in particular on treating retardation as a small parameter in the dynamics. This regime allows one to systematically describe the quantum fluctuations of the cavity degrees of freedom by eliminating the cavity variables from the equations of motion of the external degrees of freedom. We now compare these predictions with the ones of stochastic differential equations derived in Ref. Domokos et al. (2001) where the cavity degrees of freedom are treated in the semi-classical limit but included at all orders of the retardation expansion. These stochastic differential equations are here extended to our setup composed of two cavity modes Torggler and Ritsch (2014) and read:

 dxj=pjmdt, (19a) dpj=∑n=1,22ℏnkSnEn,rsin(nkxj)dt, (19b) dEn,r=(−ΔnEn,i−κnEn,r)dt+dξn,r, (19c) dEn,i=(ΔnEn,r−κnEn,i−NSnΘn)dt+dξn,i, (19d)

where and are the real and imaginary part of the positive-frequency component of the cavity field mode . The Wiener processes have vanishing first moment, , while the second moments fulfill , , and .

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 Eq. (19) predict for certain parameters samples whose temperature is 10% hotter than the one obtained with Eq. (1). Small differences are found also for the order parameters after the quenches into the bistable phase.

Figure LABEL:Fig:comparison shows a representative result of the discrepancies found after the quench protocol discussed in Sec. III.2.

The two simulations predict different stationary values for both kinetic energy and the order parameters. We believe that this discrepancy is due to retardation effects, which are neglected in the approach of Eq. (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 Refs. Niedenzu et al. (2011); Grießer et al. (2012), where the temperature of the stationary thermal distribution was corrected by the contribution due to the atoms’ localization at the minima of the self-organized lattice,

 kB~T=ℏΔ2c+κ24|Δc|+ℏω20|Δc|. (20)

Here, is the frequency of oscillation about the lattice minima in the harmonic approximation. It can be estimated using Eq. (19b) and imposing the equality

 dpj≈∑n=1,22ℏ(nk)2SnEn,rxjdt≡−mω20xjdt.

This delivers an analytic estimate of the frequency

 ω20=ωrΔ2c+κ2−Δc(α1Θ1+4α2Θ2),

where we used Eq. (4). For the parameters of the quench in Fig. LABEL:Fig:comparison, with and , we obtain , where is the temperature given in Eq. (18). Indeed, this corrected value of the final temperature is in good agreement with the discrepancy observed in Fig. LABEL:Fig:comparison (a).

This hypothesis is also consistent with the discrepancy observed in the asymptotic values of the order parameters. In fact, stationary temperature and the final values of the order parameters are related: the stationary values of the order parameters are determined by the parameters (cif. Keller et al. (2017)) and thus depend on both field intensities as well as detunings, see Eq. (4). According to this hypothesis, the asymptotic values of the order parameters for the simulation using Eq. (19) should be the ones corresponding to the system’s parameters with the corrected temperature , hence we shall minimize the free energy of Eq. (14) using , Eq. (20), instead of . This is equivalent to rescale the phase diagram in Fig. LABEL:fig:phase_diagram(a) using the prescription , and results in a smaller stationary value of the order parameter which is consistent with the discrepancies visible in Fig. LABEL:Fig:comparison (b) and (c).

## Vi 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 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 Ref. Keller et al. (2017). 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 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 . They share analogies with metastable configurations found in the GHMF when performing quenches in the microcanonical ensemble Teles et al. (2012). Since the phase diagrams of GHMF and the model here considered can be formally mapped into one another Keller et al. (2017), 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 as the one performed in Ref. Jäger et al. (2016) 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 systems remains trapped in metastable configurations is further increased. This reminds of the behavior of supercooled water Auerbach (1995); Brownridge (2011); Tao et al. (2017); Jin and Goddard III (2015); Zhang et al. (2013).

Here we have considered 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 wavelength of the cavity modes are incommensurate Habibian et al. (2013), so that the ordering mechanisms are much more strongly competing and a multitude of meta-stable states can form. A further interesting direction is operation with much colder temperatures or in the side-band resolved cooling regime Krämer and Ritsch (2014). Here it is intriguing to consider in which form meta-stable 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 Maschler et al. (2007), a process which should also be relevant in closely related schemes of simulated quantum annealing Torggler et al. (2017).

###### Acknowledgements.
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 network ”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.

### References

1. P. Horak, G. Hechenblaikner, K. M. Gheri, H. Stecher,  and H. Ritsch, Phys. Rev. Lett. 79, 4974 (1997).
2. P. Domokos and H. Ritsch, J. Opt. Soc. Am. B 20, 1098 (2003).
3. A. T. Black, H. W. Chan,  and V. Vuletić, Phys. Rev. Lett. 91, 203001 (2003).
4. K. Baumann, C. Guerlin, F. Brennecke,  and T. Esslinger, Nature 464, 1301 (2010).
5. J. K. Asbóth, P. Domokos, H. Ritsch,  and A. Vukics, Phys. Rev. A 72, 053417 (2005).
6. S. Schütz and G. Morigi, Phys. Rev. Lett. 113, 203002 (2014).
7. S. Schütz, S. B. Jäger,  and G. Morigi, Phys. Rev. A 92, 063808 (2015).
8. D. H. J. O’Dell, S. Giovanazzi,  and G. Kurizki, Phys. Rev. Lett. 90, 110402 (2003).
9. P. Münstermann, T. Fischer, P. Maunz, P. W. H. Pinkse,  and G. Rempe, Phys. Rev. Lett. 84, 4068 (2000).
10. V. Vuletić and S. Chu, Phys. Rev. Lett. 84, 3787 (2000).
11. S. Schütz, S. B. Jäger,  and G. Morigi, Phys. Rev. Lett. 117, 083001 (2016).
12. S. B. Jäger, S. Schütz,  and G. Morigi, Phys. Rev. A 94, 023807 (2016).
13. T. Keller, S. B. Jäger,  and G. Morigi, J. Stat. Mech. , 064002 (2017).
14. S. Schütz, H. Habibian,  and G. Morigi, Phys. Rev. A 88, 033427 (2013).
15. P. Domokos, P. Horak,  and H. Ritsch, J. Phys. B 34, 187 (2001).
16. V. Torggler and H. Ritsch, Optica 1, 336 (2014).
17. F. Brennecke, T. Donner, S. Ritter, T. Bourdel, M. Köhl,  and T. Esslinger, Nature 450, 268 (2007).
18. R. Landig, L. Hruby, N. Dogra, M. Landini, R. Mottl, T. Donner,  and T. Esslinger, Nature 532, 476 (2016).
19. This can be realised by assuming , giving . Another possible realisation, where , has been discussed in Ref. Keller et al. (2017), and uses two optical single-mode cavities crossing at an angle of . For a similar experimental setup see also Ref. Léonard et al. (2017).
20. J. Léonard, A. Morales, P. Zupancic, T. Esslinger,  and T. Donner, Nature 543, 87 (2017).
21. J. Dalibard and C. Cohen-Tannoudji, J. Phys. B 18, 1661 (1985).
22. A. Campa, T. Dauxois,  and S. Ruffo, Phys. Rep. 480, 57 (2009).
23. T. N. Teles, F. P. d. C. Benetti, R. Pakter,  and Y. Levin, Phys. Rev. Lett. 109, 230601 (2012).
24. A. Pikovsky, S. Gupta, T. N. Teles, F. P. d. C. Benetti, R. Pakter, Y. Levin,  and S. Ruffo, Phys. Rev. E 90, 062141 (2014).
25. D. Auerbach, Am. J. Phys. 63, 882 (1995).
26. J. D. Brownridge, Am. J. Phys. 79, 78 (2011).
27. Y. Tao, W. Zou, J. Jia, W. Li,  and D. Cremer, J. Chem. Theory Comput. 13, 55 (2017).
28. J. Jin and W. A. Goddard III, J. Phys. Chem. C 119, 2622 (2015).
29. X. Zhang, Y. Huang, Z. Ma,  and C. Q. Sun, preprint arXiv:1310.6514  (2013).
30. W. Niedenzu, T. Grießer,  and H. Ritsch, Europhys. Lett. 96, 43001 (2011).
31. T. Grießer, W. Niedenzu,  and H. Ritsch, New J. Phys. 14, 053031 (2012).
32. H. Habibian, A. Winter, S. Paganelli, H. Rieger,  and G. Morigi, Phys. Rev. Lett. 110, 075304 (2013).
33. S. Krämer and H. Ritsch, Phys. Rev. A 90, 033833 (2014).
34. C. Maschler, H. Ritsch, A. Vukics,  and P. Domokos, Opt. Commun. 273, 446 (2007).
35. V. Torggler, S. Krämer,  and H. Ritsch, Phys. Rev. A 95, 032310 (2017).
You are adding the first comment!
How to quickly get a good reply:
• Give credit where it’s due by listing out the positive aspects of a paper before getting into which changes should be made.
• Be specific in your critique, and provide supporting evidence with appropriate references to substantiate general statements.
• Your comment should inspire ideas to flow and help the author improves the paper.

The better we are at sharing our knowledge with each other, the faster we move forward.
The feedback must be of minumum 40 characters