We now present a semiclassical description of the particles’ center of mass motion and the cavity field. The coupled equations for the motion of the particles with position x j and momentum p j , and the real and imaginary part of the cavity field E r and E i evolve according to the following stochastic differential equations (Domokos et al., 2001) d x j = p j m d t , (1a) d p j = 2 ℏ k S E r ⁡ sin k x j d t , (1b) d E r = − Δ c E i − κ E r d t + d ξ r , (1c) d E i = Δ c E r − κ E i − N S Θ d t + d ξ i , (1d)

and j = 1,2 , … , N . The noise terms d ξ i , d ξ r have vanishing first moments, ⟨ d ξ i ⟩ = 0 = ⟨ d ξ r ⟩ , while the second moments fulfill ⟨ d ξ i d ξ i ⟩ = κ d t / 2 , ⟨ d ξ r d ξ r ⟩ = κ d t / 2 , and ⟨ d ξ r d ξ i ⟩ = 0 . Furthermore, the order parameter or magnetization Θ is defined by Θ = 1 N ∑ j = 1 N cos k x j . (2)

Equations 1a, b, c, and d describe the driven-dissipative dynamics of the particles that couple to a dissipative cavity mode.

To better understand the forces that are mediated by the cavity, it is useful to eliminate the cavity degrees of freedom from the dynamics. Here, we work in the limit where κ , | Δ c | ≫ k Δ p / m and κ 2 ≫ ω R N S , with ω R = ℏ k 2 / ( 2 m ) the recoil frequency. This implies that the cavity degrees evolve much faster and can be adiabatically eliminated (Schütz et al., 2013). Here, k Δ p / m is the Doppler-width and Δ p is the single-particle momentum width. Working in this regime allows us to simplify Equations 1a, b, c, and d by calculating the adiabatic stationary state of the cavity field. This is done by formally integrating the differential equations for E r and E i . The adiabatic solution is given by E r = Δ c N S Θ Δ c 2 + κ 2 , (3) E i = − κ N S Θ Δ c 2 + κ 2 . (4)

Using this result in Equations 1a, b, c, and d results in d x j = p j m d t , (5a) d p j = − 2 k V ⁡ sin k x j Θ d t , (5b)

with V = − ℏ Δ c N S 2 Δ c 2 + κ 2 .

We emphasize that V is positive since we assumed Δ c < 0 , which allows for self-organization and cavity cooling (Asbóth et al., 2005). The dynamics given by Equation 5a, b can be rewritten using an effective Hamiltonian H eff = ∑ j p j 2 2 m − N V Θ 2 ,

with d x j / d t = ∂ H eff / ∂ p j , d p j / d t = − ∂ H eff / ∂ x j . The term ∝ V is a long-range interaction potential which tries to maximize the value of Θ . The latter is the order parameter or magnetization and is used to distinguish between the self-organized and the spatially homogeneous phase. In this context, the values of cos ( k x j ) can be seen as a continuous magnetization for each atom which takes values between − 1 and + 1 . In the spatially homogeneous or paramagnetic phase, cos ( k x j ) takes random values between − 1 and + 1 such that Θ ≈ 0 . In the self-organized or ferromagnetic phase, the particles form a pattern with a periodicity that is determined by the wavelength λ = 2 π / k such that | Θ | > 0 , meaning that the individual spins fulfil either all cos ( k x j ) ≈ 1 or all cos ( k x j ) ≈ − 1 . In an experiment, the magnetization can be detected from the cavity output. This can be seen by finding the stationary state of Equations 1c, d that can be used to calculate the intra-cavity photon number I = E r 2 + E i 2 ≈ N 2 S 2 Δ c 2 + κ 2 〈 Θ 2 〉 , (6)

where we used Equations 3, 4 and the average runs over different initializations and trajectories.

After adiabatic elimination of the cavity degrees of freedom, we derived a dynamical description from a classical Hamiltonian. This implies that Equations 5a, b conserve the mean energy ⟨ H eff ⟩ for a time-independent interaction strength V . The description by means of Hamiltonian dynamics is, however, only true on a timescale where dissipative effects can be discarded (Schütz et al., 2015; Schütz et al., 2013; Schütz and Morigi, 2014; Jäger et al., 2016; Schütz et al., 2016).

In the following, we are interested in changing V very slowly such that the particles evolve mainly adiabatically but do so sufficiently fast such that dissipative effects are negligible.

We assume that the distribution function of the particles is given by a thermal state which can be seen as the stationary state of the system reached after sufficiently long times. This state is given by f t x , p = Z − 1 β t e − β t H eff

with single-particle kinetic energy E kin t = ⟨ p 2 ⟩ t 2 m = 1 2 β t , (7)

and partition function Z ( β t ) = ∫ d x ∫ d p e − β t H eff . Note that f t and β t are explicitly time-dependent. The expectation value is here defined by ⟨ h ( x , p ) ⟩ ( t ) = ∫ d x ∫ d p h ( x , p ) f t ( x , p ) with integrals ∫ d x = ∫ 0 λ d x 1 … ∫ 0 λ d x N and ∫ d p = ∫ − ∞ ∞ d p 1 … ∫ − ∞ ∞ d p N and for an arbitrary function h ( x , p ) of the atomic positions and momenta.

We assume a time dependent V and, in particular, that the temporal change of V is sufficiently slow such that the particles remain in a thermal state. With this assumption, a dynamical equation for the kinetic energy can be derived in the following.

Using Equations 5a, b and Equation 7, we obtain the dynamical evolution of the single-particle kinetic energy d E kin d t = V d ⟨ Θ 2 ⟩ d t . (8)

Furthermore, we may write ⟨ Θ 2 ⟩ = d F d y y = N V 2 E kin

with F y = ln ∫ d x e y Θ 2 .

Then, using V = 2 E kin y / N , we can rewrite Equation 8 as N d E kin d y = 2 E kin y d 2 F y d y 2 ,

and the integration of this equation leads to ∫ E 0 kin E 1 kin d E kin E kin = ∫ y 0 y 1 d y 2 y N d 2 F d y 2 , where y n = N V n / ( 2 E n kin ) with n = 0,1 . The latter can be solved using integration by parts to obtain ln E 1 kin E 0 kin = 2 N y d F y d y − F y y 0 y 1 , using the notation [ f ( y ) ] y 0 y 1 = f ( y 1 ) − f ( y 0 ) . Defining α n = y n / N = V n / ( 2 E n kin ) and performing the limit N → ∞ with α n = c o n s t we obtain the result E 1 kin E 0 kin = I 0 2 α 0 θ α 0 e − 2 α 0 θ 2 α 0 I 0 2 α 1 θ α 1 e − 2 α 1 θ 2 α 1 2 , (9) where I n is the n th modified Bessel function and θ ( α ) describes the stable solution of the equation θ = I 1 2 α θ I 0 2 α θ . (10) See Supplementary Appendix A for a detailed derivation. The value for θ calculated from Equation 10 is the mean magnetization of the particles for the given value of α . Thus, Equation 9 connects the magnetization before and after the ramp with the kinetic energy before and after the ramp. We now discuss how this result can be used to lower the kinetic energy of the particles.

In Figure 2a we plotted θ ( α ) as a function of α . It can be seen that θ ( α ) is 0 for α < 1 (paramagnetic phase) and increases for α > 1 (ferromagnetic phase) while it tends to 1 in α → ∞ . Note that only the positive solution θ > 0 has been shown, but there is also the solution − θ . This transition from spatially homogeneous (paramagnetic) to self-organized (ferromagnetic) has been described as a “phase transition.”

In the spatially homogeneous phase, for α 1 , α 0 ≤ 1 , the quotient of the kinetic energies in Equation 9 is always 1. This implies that any adiabatic change within the spatially homogeneous region will, to good approximation, not affect the kinetic energy.

However, when we assume that the coupling strength is initialized such that the particles are in the self-organized phase and ramped to a value where the particles are distributed spatially homogeneously—that is, α 0 > 1 and α 1 ≤ 1 —we obtain E 1 kin E 0 kin = I 0 2 α 0 θ α 0 e − 2 α 0 θ 2 α 0 2 . (11)

This result of the right-hand side of Equation 11 in Figure 2b is shown as black solid line. It is a monotonous decreasing function with α 0 . Therefore, we conclude that a potentially very low kinetic energy can be reached by starting the ramp from a high coupling strength V 0 . In this regime, for α 0 ≫ 1 , we obtain the asymptotic result E 1 kin E 0 kin = e 4 π α 0 , (12) where e is the Euler number. This shows that the ratio of the kinetic energies is proportional to 1 / α 0 and is plotted as a dashed gray line in Figure 2b.

We now discuss how this principle might be applicable to the driven-dissipative dynamics of particles in a cavity.

We first initialize the particles in a strongly self-organized thermal state with kinetic energy E kin Then, we ramp V exponentially as V t = V 0 ⋅ 1 0 − 5 t t ramp , (13)

for different ramping times t ramp . While the choice of an exponential ramp is a technical detail, it allows for a rather fast change of the interaction strength for large values of V ≈ V 0 and slow changes for V ≳ 0 . This seems a good compromise between being fast and remaining approximately adiabatic. We study the dynamics of the full system, including the cavity degrees of freedom (Equations 1a, b, c, d) and the dynamics where the cavity degrees of freedom are eliminated (Equations 5a, b). Figures 3a, b show the dynamics of the kinetic energy following the ramp for a ramping time of t ramp = 10 ω R − 1 and t ramp = 100 ω R − 1 , respectively. To put these values into actual numbers, we provide an explicit example and use the value ω R = 2 π × 2  kHz from Hosseini et al. (2017) for ¹³³Cs, which results in ramping times t ramp ≈ 1  ms and t ramp ≈ 10  ms , respectively. The dashed line shows the result using the conservative dynamics (Equation 5a, b), while the solid line represents the full dissipative dynamics (Equation 1a, b, c, d). Both curves for both ramping times show decreased kinetic energy. There is good agreement of both dynamics on short timescales while we observe discrepancies for the longer ramping time. Therefore, we expect that for sufficiently short times, dissipative effects are still negligible while they affect the dynamics on longer timescales. This relies on a timescale separation of dissipative and conservative forces that relies on (i) the number of particles and (ii) the typical timescale separation of motion and cavity relaxation—that is, k Δ p / m ≪ κ . This was studied by Jäger et al., (2016) and Schütz et al., (2016) and observed in Wu et al. (2023). In conclusion, this preliminary analysis demonstrates that there must be an optimal ramping time for which the lowest possible temperature can be achieved.

While our original assumption was that the ramp is close to adiabatic, we expect this assumption to fail, especially because the system parameters are ramped across a phase transition. An observable to test this is the kurtosis K t = 〈 p 4 〉 t 〈 p 2 〉 t 2 . (14)

The kurtosis is K = 3 for a Gaussian state and deviates from 3 for non-Gaussian states. Figures 3c, d show the kurtosis with the same labeling for the two different ramping times. In Figure 3c, we observe that the kurtosis remains close to 3 for times t ≲ 7 ω R − 1 , while it deviates for longer times as soon as the value of V crosses the phase transition line. In Figure 3d we observe the same for the simulations of the conservative dynamics (Equations 5a, b) and times t ≲ 70 ω R − 1 , while the simulation of the full dissipative dynamics (Equations 1a, b, c, d) shows values of K ≠ 3 on much shorter timescales. This finding supports our claim that the dynamics do not remain adiabatic across the phase transition. In addition, the discrepancies between the conservative and dissipative dynamics predict a dissipation-assisted creation of non-Gaussian states that is closely related to the dissipation-assisted stabilization of non-Gaussian states predicted in Schütz et al. (2016).

We now analyze the dependence of the minimum achievable temperature on the ramping time t ramp . Figure 4 compares the values of the final kinetic energies E 1 kin = E kin ( t ramp ) for different ramping times t ramp , different particle numbers, and different ratios of κ / ω R . The black line with symbols are calculated using simulations of Equations 1a, b, c, d with N = 50 (circles), N = 100 (crosses), and N = 200 (pluses), where we show κ = 400 ω R in Figure 4a and κ = 40 ω R in Figure 4b. These are realistic values, and the lower value of κ = 40 ω R is close to that realized in Hosseini et al. (2017). For both simulations, the system has been initialized with a pumping strength of α 0 = V 0 / ( 2 E 0 kin ) = 50 and E 1 kin = ℏ κ / 4 . The results of the simulations predict a local minimum of the kinetic energy in the range ω R t ramp = 10 − 100 , thus showing that there is an optimal ramping time. In general, this optimal ramping time is shorter for smaller particle numbers. In addition, the minimum achievable kinetic energy is larger for smaller particle numbers. This should be due to the timescale separation between the conservative and dissipative forces that becomes larger for increasing particle numbers.

We also find that the result of the kinetic energy for κ = 40 ω R (b) appears to be slightly displaced to larger ramping times with respect to the simulations for κ = 400 ω R (a). We expect that this is due to a violation of the adiabaticity criteria, κ t ramp ≫ 1 and k Δ p t ramp / m ≫ 1 , that are not fulfilled for smaller values of k Δ p / m and κ and short ramping times in Figure 4b. For completeness, we have also included a simulation of Equations 5a and b that does not include dissipation and noise. The results are visible as gray dashed lines in Figure 4. Those curves are monotonically decreasing, thus showing that noise and dissipation are the origins for the local minima in the kinetic energy in the full simulations. The theoretical minimum of the achievable kinetic energy is shown as a gray dashed–dotted line. It is calculated using Equation 12 and α 0 = 50 , resulting in E 1 kin / E 0 kin ≈ 0.04 . We observe that the simulation without dissipation and noise (gray dashed line) converges to this theoretical minimum in the limit ω R t ramp → ∞ .

We now use this gained insight to minimize the kinetic energy of particles that are initially in a spatially homogeneous configuration. We thus assume that the initial state is a thermal state with temperature k B T in = ℏ κ / 2 that can be reached by cavity cooling for Δ c = − κ . The actual choice of this state is rather arbitrary but should be sufficiently cold such that the ensemble can be cavity (laser) cooled.

In a first stage, we perform a quench in the driving-laser intensity determined by S such that V has a value V fer for the system to reach a state well inside the self-organized phase. Over a very long time, the system again reaches a stationary state which is thermal. To be consistent, for large laser intensities, we need to take corrections of the final temperature into account which come from the laser driving power. This final kinetic energy is given in the well-organized regime (Niedenzu et al., 2011; Grießer et al., 2012) by E fer kin = k B T fer 2 = ℏ Δ c 2 + κ 2 + 4 ω 0 2 − 8 Δ c , (15) where ω 0 2 = 4 ω R V fer ℏ is the effective trapping frequency.

In a second stage, we consider a ramp from V fer back to a value close to 0 such that both magnetization and kinetic energy are adiabatically reduced. If we assume that the system remains adiabatic in a thermal state, the optimum final kinetic energy can be approximated using Equation 12 by E par kin = e 2 π E fer kin 2 V fer . (16)

By minimizing Equation 16 with respect to V fer , we find E min kin = e 2 π ℏ ω R Δ c 2 + κ 2 Δ c 2 (17)

at an optimum value of V fer opt = ℏ Δ c 2 + κ 2 16 ω R . (18)

This minimum kinetic energy is of the order of the recoil energy ℏ ω R .

We now consider simulations that test this prediction. Following the procedure of the first stage, Figure 5a shows the dynamics of the kinetic energy after a quench from V ≈ 0 to V = V fer opt (Equation 18). Initially, a rapid increase is observed in the magnetization determined by Θ 2 . This can be seen in Figure 5c, which plots the cavity field determined by Equation 6. Since energy is conserved on short timescales, the kinetic energy is also exponentially increasing, reaching a maximum of E kin ≈ 6.5 × 1 0 3   ℏ ω R . On longer timescales, dissipation guides the system toward a thermal state with a temperature given by Equation 15 (horizontal gray dashed line in Figure 5a) with corresponding magnetization (horizontal gray dashed line in Figure 5c). The steady-state magnetization has been calculated using ⟨ Θ 2 ⟩ ≈ θ 2 where θ is the solution of Equation 10. Figure 5e plots the kurtosis (Equation 14) that starts and ends at a value close to K ≈ 3 , suggesting that both the final and initial state are thermal.

Following the second stage, we ramp the coupling strength according to Equation 13 with a ramping time of t ramp = 10 ω R − 1 that we have found to be close to optimal for the choice of parameters in the previous subsection. Figure 5b shows the decrease of the kinetic energy that eventually reaches a value that is of the order of the recoil energy. We observe a final kinetic energy of E final kin = 1.3 ℏ ω R , whereby Equation 17 predicts a similar value of E min kin ≈ 0.9 ℏ ω R . During this process, the field intensity and mean magnetization of the system are decreased (Figure 5d), reminiscent of adiabatic demagnetization. We emphasize that the adiabatic process is much faster than cavity cooling in the first stage. This is visible by comparing the time axes in Figures 5a and b. This difference in the timescales comes from the fact that the adiabatic process is collective. The kurtosis, visible in Figure 5f, remains 3 during the ramp until the phase transition line is crossed and the process is no longer adiabatic.