Driven-dissipative many-body systems with mixed power-law interactions: Bistabilities and temperature-driven non-equilibrium phase transitions

We investigate the non-equilibrium dynamics of a driven-dissipative spin ensemble with competing power-law interactions. We demonstrate that dynamical phase transitions as well as bistabilities can emerge for asymptotic van der Waals interactions, but critically rely on the presence of a slower decaying potential-core. Upon introducing random particle motion, we show that a finite gas temperature can drive a phase transition with regards to the spin degree of freedom and eventually leads to mean-field behaviour in the high-temperature limit. Our work reconciles contrasting observations of recent experiments with Rydberg atoms in the cold-gas and hot-vapour domain, and introduces an efficient theoretical framework in the latter regime.

The idea that matter rapidly relaxes towards a thermal ensemble [1] has led to a profound understanding of many macroscopic phenomena within the powerful framework of equilibrium statistical physics. More recently, the experimental success in realizing synthetic many-body systems with controllable dissipation has motivated broad explorations of non-thermal steady states [2]. Examples include cold atoms in cavities [3], semiconductor excitonpolariton condensates [4], trapped ion crystals [5] and laser-driven Rydberg gases [6]. The interplay of coherent and dissipative dynamics in such driven-dissipative systems generates non-equilibrium phases and transitions that may have no equilibrium equivalent. An evident example of such distinct behaviour is the emergence of multiple steady states.
In this Letter we address this problem through numerical studies of the driven-dissipative dynamics in Isingspin ensembles with power-law interactions. We point  (color online) (a) An ensemble of interacting two-level systems is driven coherently with a coupling strength Ω and frequency detuning ∆ in the presence of decay and de-phasing with rates Γ and γ, respectively. (b) The potential eq.(1) (solid line) interpolates between different power-laws and accurately describes the actual interaction between excited Rubidium atoms (dots), shown exemplarily for Rb(70S 1/2 ) atoms. (d) Hysteresis with bistable steady states, whose typical spatial configurations and correlation functions are illustrated in panels (c) and (e). Blue spheres show excited particles, while the opacity of the red dots indicates the excitation rate γ ↑ of particles in state |g . out the importance of fluctuations for the topology of the non-equilibrium phase diagram and draw a direct connection to the form of the spin-spin interactions. In particular, bistability cannot occur under the common assumption of pure van der Waals (vdW) interactions, but instead requires a short-distance dipolar potentialcore. Such a short-distance behaviour is characteristic for Rydberg atoms [ Fig.1(b)], but typically neglected in theoretical models. Upon incorporating particle motion we reveal a temperature-driven phase transition to bistable non-equilibrium steady states. In the high-temperature limit we demonstrate a cross-over to MF-behaviour, offering an explanation of thermal-vapour experiments [8]. We consider a random ensemble of N two-level systems with states |g i and |e i (i = 1, ..., N ) coupled by a laser field with a Rabi frequency Ω and frequency detuning ∆ [ Fig.1(a)]. The associated unitary dynamics is governed by the HamiltonianĤ = Ω where V (r ij ) denotes the interaction potential of two particles at positions r i and r j and r ij = |r j − r i |. The N -body density matrix,ρ, of the system evolves according to the master describes one-body decoherence processes. We account for decay of the excited state (L i,0 = √ Γ|g i e i |) with a rate Γ and de-phasing of the laser-driven transition (L i,1 = √ γ|e i e i |) with a rate γ.
From now on we use dimensionless quantities, scaling time by Γ −1 and length by the critical radius r b , defined as the particle distance at which V (r b ) = Γ + γ. We consider a potential [22] V (r) which features a dipolar potential core (V ∼ 1/r 3 ) below a distance r vdW and vdW interactions (V ∼ 1/r 6 ) for r > r vdW . Herer = r/r b and ξ = r vdW /r b denotes the vdW distance relative to the blockade radius. Since r b defines the typical distance between Rydberg atoms limited by the excitation blockade, the value of ξ characterizes the importances of dipolar interactions. Eq.(1) reproduces the characteristics of Rydberg atom interactions [6], as illustrated in Fig.1(b) by a comparison to numerical results [23] for nS 1/2 states of Rubidium atoms. Provided that Ω/(Γ + γ) 1 one can adiabatically eliminate the dynamics of the off-diagonal densitymatrix elements and obtain a closed evolution equation for the diagonal ρ S,S [24][25][26][27]. The matrix elements ρ S,S describe the population of N -body configurations S ≡ (s 1 , . . . , s N ). Here, s i is an effective spin variable denoting the ground (s i = 0) and excited (s i = 1) state of the i th particle. Introducing the state vector S i ≡ (s 1 , . . . , 1 − s i , . . . , s N ), the resulting master equation can be written aṡ where the single-body (de)excitation rates are given by γ The rates are determined by two parameters: the scaled Rabi frequencyΩ = Ω/ Γ(Γ + γ) and the scaled frequency detuning∆ i (S) = ∆ i (S)/(Γ + γ). The latter consists of the laser detuning ∆ and the interaction-induced level shifts from nearby excited particles, [26]. Exact quantum-trajectory simulations for small systems [28] established the accuracy of this approach for Ω Γ + γ. Note that this condition does not restrict our parameters, and permitsΩ > 1 if γ > Γ, which is often the case in experiments [7].
The obtained effective master equation can be solved via kinetic Monte Carlo (MC) sampling [29]. To this end, we randomly sample N particle positions from a cubic volume with periodic boundary conditions and an edge length L, chosen to be much larger than r b and r vdW . The corresponding dimensionless density ρ = N r 3 b /L 3 defines the number of particles within a given blockade volume r 3 b . To calculate the excitation spectrum, we perform positive and negative scans of the detuning∆ with a corresponding chirp rate ±κ. Observables are calculated from an ensemble average over many realizations of particle disorder configurations.
We find two distinct steady states with a low and a high excitation density ρ e . The low-density phase corresponds to a dilute gas of excited pairs [ Fig.1(e)], formed by resonant sequential excitation of particles at a distance r ∆ for which the potential V (r ∆ ) = ∆ compensates the detuning. The correlations in high-density phase [ Fig.1(c)] do not feature strong ordering on the length scale r ∆ and resemble a liquid of repulsive excitations.
To investigate the stability of these two phases, we calculate the excitation spectrum for negative and positive scans of the detuning. For a proper choice of parameters, both phases are indeed found to coexist over a finite range of∆ where the excitation density shows hysteretic behaviour, showing qualitative resemblance to MF predictions [11]. However, in contrast to MF expectations [14], we find no evidence of bistability for pure vdW interactions (ξ → 0). This is illustrated in Fig.2(a) and (b) where we show typical excitation spectra for small and large values of the vdW radius. For small ξ, the excitation blockade prevents particles from exploring the dipolar region of the interaction potential and one finds a smooth resonance curve with a unique steady state. However, once the short-distance 1/r 3 -behaviour of V (r) starts to become significant the system develops bistable steady states beyond a critical driving strength [ Fig.2(b)].
This behaviour can be understood by considering the effect of the potential form on energy level fluctuations. Spontaneous decay inevitably causes |e → |g transitions and thereby temporal fluctuations of the corresponding interaction-induced level shifts ∆ i . For ξ 1, the total level shift, ∆ i , of an excited particle typically results from a small number of excitations in close proximity. Hence, a single decay event will cause a substantial change of ∆ i and disturb the excitation dynamics. The resulting large density fluctuations [ Fig.3(b)] prevent the formation of two distinct phases. For large ξ, a large number of excitations within a distance ξ collectively contribute to ∆ i , such that potential fluctuations are greatly reduced. To validate this picture, we have traced the microscopic steady state dynamics for two different values of ξ and otherwise identical parameters and average densities, ρ e [dots in Fig.2(a) and (b)]. By recording the maximum change, δ, of the level shift of excited particles due to a de-excitation, we construct the spectrum of potential fluctuations P (δ) from the long-time microscopic steady state dynamics. As seen in Fig.2(c), one indeed finds a broad distribution for ξ = 0.5 with extended tails well beyond the average potential shift V . On the contrary, for ξ = 2.5 [ Fig.2(d)] P (δ) is sharply peaked around small δ V and drops rapidly for larger values. It is this strong suppression of fluctuations [ Fig.3(d)] that facilitates the formation of bistable steady states in the limit of large ξ.
The microscopic steady state dynamics provides fur-ther insights about the transition between these two regimes. The non-equilibrium phase diagram shown in Fig.3(a) reveals a finite region of bistability at large ξ which ultimately closes upon decreasing ξ. In between these two limits, the low-and high-density phases, coexisting as long-lived metastable states, are connected by a first order phase transition over a finite range of ξ. This transition, generally obscured by the ensemble average over random particle configurations, is revealed by the counting statistics of a single N -body trajectory, as demonstrated in Fig.3(b)-(d) where we show the excitation-density distribution for a single particle configuration at a low chirp rate κ = 10 −8 . For 1.8 ξ 2 both phases dynamically coexist and yield a persistent bimodal counting statistics. The ensemble average of the corresponding transition point yields the first order transition line shown in Fig.3(a) which ends in a critical point around 1.7 < ξ < 1.8. For a broader characterization of the conditions leading to bistability we have calculated the asymptotic hysteresis area A 0 by extrapolating A(κ) to the limit κ → 0 [inset of Fig.4(b)]. Upon changing ξ, A 0 indicates a continuous transition with a critical exponent ∼ 1 [ Fig.4(b)]. The critical ξ expectedly decreases with particle density [ Fig.4(a)]. Yet, the apparent saturation of the transition line at large densities provides further indication for the absence of bistable behaviour for systems with dominant vdW interactions (ξ < 1).
We have also investigated the average switching time between the two steady states and found the expected power-law divergence upon approaching the transition points in Fig.3(a). Yet, the corresponding critical exponents strongly depend on ξ, which stands in contrast to thermal-vapour experiments [8] that observed a universal MF exponent of 0.5 [30]. To resolve this issue we now consider thermal particle motion which diminishes correlations, and thereby alters the spectrum of fluctuations. For simplicity we neglect inter-particle forces and adopt an ideal gas description with an equilibrium velocity distribution and dimensionless thermal velocity v th , measured in units of r b Γ. Tracking the evolving particle positions, now requires fixed-time-step MC simulations [29] of the spin dynamics.
As shown in Fig.5(a), thermal motion drives a continuous phase transition to bistability. At high temperatures the hysteresis area saturates to a finite value that can be understood within the following MF treatment. Assuming that rapid thermal motion completely randomizes any spatial excitation structures we can neglect correlations in eq.(2) to obtain a closed equation,ρ e =γ ↑ −(1+2γ ↑ )ρ e , for the average excitation density. Averaging the microscopic rates γ (i) ↑ over the uncorrelated ensemble yields a closed expression for the MF excitation ratē , our hightemperature results indeed approach this ensemble averaged mean-field (eaMF) limit. In contrast to corresponding lattice MF models [14,17,18], the functional ρ e -dependence ofγ ↑ depends strongly on the shape of the interaction potential. In particular, for ξ = 0 one finds Re(f ) = Im(f ) ∝ ρ e √ k, which implies that no phase transition can occur for pure vdW interactions. The numerically obtained eaMF transition line (Fig.4) demonstrates that this remains true for finite ξ < 1.
Finally, we put our findings into the context of recent experiments. Ref. [7] reports bimodal counting statis-tics of Rydberg excitations in a cold gas of Rb atoms excited to 70S 1/2 states. The quoted laser linewidth γ/2π ≈ 500 kHz and Rabi frequencies Ω < (Γ + γ) are within the regime of validity of the present theory. The vdW radius of ξ ≈ 0.3 implies that the conditions of [7] do not promote bistable steady states. Bimodality at short excitation times, τ , can, however, result from transient relaxation effects [28], while for larger τ dipolar state-mixing induced by black-body radiation [31,32] on a timescale τ bbr < τ [7,33] significantly affects the gas dynamics as observed in other experiments [34][35][36]. Note that the temperature corresponding to a thermal velocity of r b Γ, for typical values of r b ≈ 11µm and Γ −1 ≈ 200µs [7], can be as low as ∼ 30µK such that atom motion can be a factor even in cold-gas experiments. Consequently, thermal gases [8,37] are deep in the high-temperature limit, v th 1, and their measured excitation dynamics can be understood within the outlined eaMF approach. Importantly, the thermally activated transition to MF behaviour explains the emergence of the dynamical MF exponents [14,30] observed in [8].
In summary, we have investigated driven dissipative spin ensembles with competing power-law interactions. The steady-state of our Master equation (2) shares the same MF limit as that [9,14] obtained from the exact quantum evolution. Yet it accounts for classical correlations and fluctuations which turn out important for the non-equilibrium physics of such systems. As a striking consequence, the specific shape of the interaction potential was found to play a key role for the non-equilibrium phase diagram despite its general finite-range nature, dropping rapidly as ∼ 1/r 6 . The spatial extent of the inner dipolar potential can be tuned by external static [38] or microwave [39] fields, which should permit explorations of the predicted phase diagram (Figs. 3 and 4) in future experiments. We showed that thermal particle motion can drive a transition to bistablity and ultimately causes MF behaviour to emerge. The non-equilibrium phase transition takes place at a surprisingly low temperatures in the µK-to mK-domain, which should enable its observation in cold atom experiments. The established high-temperature MF limit (eaMF) provides a consistent explanation of recent thermal-vapour experiments [8] and, permits future analysis of multi-level atoms and more compex interactions that may occur in such systems [8,37].
We thank Rick van Bijnen, Simon Gardiner, Hendrik Weimer and Igor Lesanovsky for valuable discussions. This work was supported by the EU through the FET-Open Xtrack Project HAIRS, the FET-PROACT Project RySQ, and the Marie-Curie ITN COHERENCE and also by Durham University and the UK EPSRC grants EP/M014398/1 and EP/M013103/1.