Quenches across the self-organization transition in multimode cavities

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.


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][2][3][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][6][7][8][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 λ and/or on a lattice with half the period l 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

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 a n propagate in a direction orthogonal to the cavity axis and are quasi-resonant with the standing wave cavity modes ( ) nkx cos with frequency w n c, and wave number nk ( = n 1, 2) 5 . The lasers have frequency w n p, 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 w n a, , where W n p, and g n are the laser and vacuum Rabi frequency, respectively. Spontaneous scattering processes are suppressed when the absolute value of the detuning w w D = - . The relevant dissipative process is given by cavity decay, and we denote by k 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 Dp which, at all instants of time, is orders of magnitude greater than the photon recoil k:  D  p k [13,14,20]. The corresponding stochastic differential equations read as 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 a 1 and a 2 into the (correspondingly) resonant cavity modes with spatial mode functions ( ) kx cos (red) and ( ) kx cos 2 (green) and loss rates k 1 and k 2 , respectively. (b) Sketch of the atomic density distribution ( ) n x along the cavity axis (in units of k 1 ) for the four possible stationary self-organized orders. On the right we report the corresponding values of the quantities Q 1 and Q 2 , which signal the Bragg order in modes 1 and 2, respectively. See text and [13] for details. 5 This can be realised by assuming . Another possible realisation, where w w » c,1 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]. Here, = W D S g n n n n a, is the amplitude of coherent scattering by a single atom and has the dimension of a frequency, while Finally, the parameter quantifies the Bragg ordering of the atoms in the cavity mode with wave number nk. In particular, Q = | | 1 n when the atoms are localized either at the maxima or at the minima of ( ) nkx cos , which is the configuration which maximizes the intracavity field intensity. We identify Q 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 nonvanishing value of Q 1 , corresponding to a Bragg grating with period l p = k 2 . Similarly, 'short-wavelength order' refers to a configuration with Q ¹ 0 2 , corresponding to a Bragg grating with l 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 . For details refer to [13]. (see equation (5)). In this case the atoms' distribution in the steady state reads as [13]  while  b ( )denotes the partition function: and  p D = 2 as the single particle unit phase space volume. In the following we assume that the cavity decay rates are equal, 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 with k B as the Boltzmann constant. The steady-state temperature T has the same functional dependence on D c and κ as for a single-mode cavity [7,14]. We can further define the free energy per particle  using formal equivalence with the canonical ensemble of equilibrium statistical mechanics [7]: Following the procedure detailed in [7,12,13] we can calculate  as is the recoil frequency. We determine the value of  in an appropriately defined thermodynamic limit, which consists in keeping a n constant for  ¥ N . The global minima of F are the resulting stationary phases. The corresponding points q q ( ) , 1,min 2,min where F achieves its minimum are the stationary values for the order parameters q Q = 1 1 , m i n and q Q = 2 2 , m i n ; they are determined by a 1 and a 2 . When the fields are sufficiently weak, then Q = Q = 0 1 2 , 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][22][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 , 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, Q = 0 1 , while Q 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 Q = 0 1 , where photons scattered into the long-wavelength mode destructively interfere. However, Q 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 Q 2 .
The resulting phase diagram in the a a -1 2 plane, shown in figure 2, reproduces that in [13]. The phases are separated by either firstor 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 Q -Q 1 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 Q = 0 1 . We call this region bistable, which refers to the existence of a second, metastable state in which the system can be dynamically trapped.

Dynamics of self-organization
We now examine the dynamics of the system when the values of a 1 and a 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 In what follows we perform numerical simulations of equation (1) using the parameters of a gas of 85 Rb atoms. In particular, we take p l = k 2 with l = 780 nm as the wavelength of the D 2 line. The corresponding recoil frequency is w p =2 3.86kHz r . The cavity linewidth is taken to be k p =2 1.5MHz, so that k w » 388.6 r . A possible realization of the two-mode setup here considered has been discussed in [13,19].

Sudden quench into the ferromagnetic phase
We first consider sudden quenches from a a , 1i 2i in the paramagnetic phase to a a , 1f 2f in the ferromagnetic phase, keeping a a a a = =5 1i 2i 1f 2f . 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  Q ( ) t for order parameters Q 1 and Q 2 as a function of time for . It is defined as a time sequence of normalized histograms: where Θ is calculated on each trajectory of the simulations with the stochastic differential equations and its value is determined according to the precision D Q of the grid in Θ. We observe that at a given time scale of the order of k 10 2 ,  Q ( ) 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 Q 2 , which is weakly pumped, substantially grows to a positive value long after the symmetry breaking. The fact that  Q ( ) t 2 vanishes for negative Q 2 values comes from the ordering of the atoms close to the anti-nodes of the dominant long-wavelength mode field ( ) kx cos (see figure 1). The distributions  Q ( ) n at the asymptotics are reported in the right panels of figure 3. They are obtained by averaging  Q ( ) t n over times  k t 10 6 , where a stable configuration is reached. Formally 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 shortwavelength order parameter.  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 k 10 6 . Meanwhile the fluctuation dQ 1 relaxes to a much smaller value than the asymptotic value of dQ 2 reproducing the widths of the distributions in the right panels in figure 3. The time evolution of á Q ñ | | 1 , 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 á Q ñ | | ; 1 (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 k 10 6 can here be identified as the one at which the asymptotic state is reached for  N 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Ñ 200 it is of the order of k t 10 4 . 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 singleparticle kinetic energy and of the kurtosis  = á ñ á ñ p p 4 2 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  = 3 Gauss . 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 a a  , 1 1i 2i but with different endpoints a a , 1f 2f . We take different endpoints from the paramagnetic to the ferromagnetic phase, under the constraint a a = 5 1f 2f . The circles in figure 6 correspond to the numerical results for 100 particles at time , where we expect the system to reach the steady state. 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 Q < ( ) P 0 t 2 that Q 2 is negative at t: We note that in the paramagnetic phase (homogeneous spatial distribution) we expect Q <  ( . In contrast, due to the given mode structure we expect that Q <  ( ) P 0 0 t 2 for long-wavelength ordering in the ferromagnetic phase. Indeed, as a 1f increases across the critical value, Q < ( ) P 0 t 2 quickly drops down to zero.

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 a a a a = =1 1i 2i 1f 2f . In this parameter region (the bistable phase) the free energy exhibits a local minimum, which is nematic. As in the 2i are vanishingly small and the atoms are at the corresponding stationary distribution, whose temperature is determined by the detuning D 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' Θ-distribution for a a = = 2 1f 2f and k D =c . In contrast to the previous section, here a finite fraction of trajectories gets trapped in the nematic phase with a vanishing value of Q 1 and a finite probability that Q 2 takes negative values. This is visible in the small extra peaks in  Q ( ) 1 and  Q ( ) 2 (right panels). The trapping occurs at the time scale of the violent relaxation, and it seems stable over times of the order of k 10 6 . 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 á Q ñ | | 2 reaches the same stationary value (in reality its value decreases slightly with N), the asymptotic value  of á Q ñ | | 1 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 dQ = áQ ñ -á Q ñ | | 1 1 2 1 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 Q < ( ) P 0 t 2 at the asymptotics and as a function of a 1f in figure 10. As in figure 6, it vanishes upon leaving the paramagnetic phase, but increases again as the a a , 1f 2f chosen are deeper into the bistable phase of figure 2(a). Correspondingly, á Q ñ | | 1 starts to decrease as a 1f increases, which suggests that from this point on the depth of the local minimum grows. The value of the order parameter á Q ñ | | 2 at which Q < ( ) P 0 t 2 starts to grow again signifies a threshold, above which the local minimum is sufficiently deep to stably trap particles.

Slow ramp into the bistable phase
We now consider linear ramps of a ( ) t n across the transition region separating the paramagnetic from the bistable region. The ramp protocols have duration τ and sweep between the values e a [ ] , nf , with e  1. In particular, a e a = + , a ( ) t n is constant and equal to a nf . Note that a sudden quench is the limit t  0 of a linear quench. We choose to vary the values of a ( ) t n along the rightmost green line in figure 2(a), so that a a = ( ) ( ) t t 1 2 at all instants of time, with a nf in the bistable phase. We further keep D c constant, and vary only the pump intensity. This means that the asymptotic temperatures at each value of a n are equal. Figure 11 shows the dynamics of the mean absolute value of the order parameters for a a = = 2 1f 2f for linear ramps with different durations τ. 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 á Q ñ | | 1 and á Q ñ | | 2 . This behavior seems to be initiated at the instant of time when the parameters a ( ) t n cross the critical point of the phase diagram. Moreover, for sufficiently slow ramps á Q ñ | | 1 approaches the asymptotic value of the free energyʼs global minimum, signaling stationary long-wavelength order.
We further note that for  t k 10 3 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 t k > 10 4 . We believe that this behavior is determined by the ramp duration τ with respect to the time scale of the transient dynamics, and thus by the time the parameters a ( ) t n 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 with the detuning kept constant k D =c (red circle (c) in figure 2(a)). The left panels display the contour plot of the distribution  Q ( ) t (equation (18)) for Q = Q 1 and Q = Q 2 as a function of time (in units of k 1 ). The distribution is extracted from the numerical simulations using equation (1) for N=100 atoms and 1000 trajectories. The right panels display distributions  Q ( ) 1 and  Q ( ) 2 as a function of Q 1 and Q 2 , respectively (see equation (19)). See figure 3 for further details.
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 á Q ñ | ( )| t 1 and á Q ñ | ( )| t 2 at k =t 3.77 10 6 , where the curves of figure 11 reach an asymptotic behavior. Self-organization in the long-wavelength grating depends on the ramp duration τ and is found for t k > 10 4 . Note that short-wavelength order quantified by á Q ñ | ( )| t 2 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 (a a =~1 1 2 ), 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: a a = =1.1 . The detuning D c is kept constant during the evolution.  Figure 14 shows the time evolution of the mean absolute values of the order parameters for different time intervals τ 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 τ the atoms are cooled into the stationary state at a~1 n . 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 a = 2 n , 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 k =t 3.77 10 6 , as a function of the time between the two quenches. These values are compared to the  and k D =c . 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 τ spent close to the transition point.

Cooling into organized structures
We now analyze sudden quenches of the parameter a n starting with different initial single-particle momentum widths. A possible realization is a quench in the detuning since D 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 , to a a = =2 1f 2f . As before, at t=0 the initial state of the atoms is the steady state, equation (8), for a a =    Figure 17 shows the time evolution of the mean absolute values of the order parameters for different values of T ini ranging from T 0.1 0 up to T 5 0 . The asymptotic value of á Q ñ | | 1 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 T 2 ini 0 (and even more for = T T 5 ini 0 ) the system stays relatively hot over time and as a function of the ramp duration τ (in units of k 1 ) for the same parameters as in figure 11. The dashed horizontal lines show the steady-state value predicted by the free energy, equation (16). scales of the order of k 10 4 . 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 k > t 10 and k < t 10 3 . For k > t 10 3 dissipation becomes significant and á Q ñ | | 1 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 á Q ñ = | | | 0.5 1 thres , 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][25][26][27][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 á Q ñ | | 1 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 k =t 3.77 10 6 as a function of the initial temperature for = N 100, 200. One would expect that á Q ñ | | 1 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.

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

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 for k D =c 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 a a , we observe that 16.9% (15.4%) of the trajectories are in the nematic phase using equation (1)(equation (22)).