Universal non-equilibrium properties of dissipative Rydberg gases

We investigate the out-of-equilibrium behavior of a dissipative gas of Rydberg atoms that features a dynamical transition between two stationary states characterized by different excitation densities. We determine the structure and properties of the phase diagram and identify the universality class of the transition, both for the statics and the dynamics. We show that the proper dynamical order parameter is in fact not the excitation density and find evidence that the dynamical transition is in the"model A"universality class, i.e. it features a non-trivial $\mathbb{Z}_2$ symmetry and a dynamics with non-conserved order parameter. This sheds light on some relevant and observable aspects of dynamical transitions in Rydberg gases. In particular, it permits a quantitative understanding of a recent experiment [C. Carr et al., Phys. Rev. Lett. 111, 113901 (2013)] which observed bistable behaviour as well as power-law scaling of the relaxation time. The latter emerges not due to critical slowing down in the vicinity of a second order transition, but from the non-equilibrium dynamics near a so-called spinodal line.

In this context, a class of systems that offers a rich and intricate physics is represented by so-called Rydberg gases [26][27][28][29][30][31], i.e., atomic clouds in which atoms are laser-excited to high-lying energy levels. The main consequence of the population of such orbitals is a considerable increase [27,28] in the interaction strength. This is at the heart of several non-trivial dynamical phenomena, both for closed systems undergoing coherent evolution and showing enhanced spatial (anti-)correlations [32][33][34][35], and for open ones, in which the interplay between driving and dissipation leads instead to intermittency [36], glassy behavior [37] and bistable behavior [38].
The dissipative case has been recently studied via a mean-field approach [39][40][41], numerical calculations in one dimension [42][43][44] and an approximate rate equation description in higher dimensions [31,[45][46][47]. These investigations highlighted the presence of various stationary regimes and the existence of first and second order phase transitions. In addition, experiments have started to probe the static and dynamic features of these systems revealing a bimodal behavior of the excitation density [29] and the presence of an optical bistability [30].
The aim of this work is to shed light on the bistable transition in a dissipative Rydberg gas with particular fo-cus on its dynamics and to connect the findings to recent experimental studies. For the stationary state, the transition is related to the spontaneous breaking of a Z 2 symmetry and falls into the Ising universality class. The effective static order parameter is an appropriately shifted Rydberg excitation density. The dynamics is found to be of "model A" type according to the standard classification of Ref. [4], i.e., it is akin to a classical Ising model subject to a noise-induced spin-flipping process which does not preserve the total magnetization. However, within the dynamical framework it becomes clear that the dynamical order parameter is not formally identical to the Rydberg excitation density and the Z 2 symmetry identified in the static case must be non-trivially generalized. Linking to recent experimental studies [30], we note that the dynamic transitions observed there take in fact place near the so-called spinodal lines of the mean-field phase diagram. The connection established to "model A" physics allows us, moreover, to extract a universal scaling law for relaxation times for which quantitative agreement with experiment is found. We believe that this perspective will be useful for analyzing and understanding the dynamical phenomena observed in other related experiments, such as the one presented in Ref. [29].
The model.-We employ the standard description of a Rydberg gas in terms of (fictitious) interacting spin-1/2 particles, where the states |↓ and |↑ correspond to the atomic ground and Rydberg states respectively. The many-body dynamics of the system's density matrix ρ is governed by the quantum master equation (QME) Here Ω is the Rabi frequency and ∆ the detuning of the excitation laser with respect to the ground state -Ryd-arXiv:1406.1015v3 [cond-mat.stat-mech] 10 Jun 2014 berg state transition. The interaction between two atoms positioned at r k and r p is, e.g., of van der Waals type V kp = C 6 | r k − r p | −6 . Moreover, we have defined the excitation density n k = (1 k + σ z k )/2, with { σ x k , σ y k , σ z k } being the usual quantum spin operators acting on the kth site. Dissipation within our model is described by the dissipator of Lindblad form Relating to previous experimental observations [27,48], we account for two dissipation mechanisms: One is independent atomic decay (at rate Γ) from the Rydberg state to the ground state, with the corresponding jump opera- The second one is dephasing of the Rydberg state relative to the ground state, occurring at rate K with L 2k = √ K n k . Mean-field equations of motion.-A mean-field treatment of the Rydberg gas has been already conducted to some extent in other works, see, e.g., [39]; here we just briefly summarize the derivation of the equations of motion. We consider the complete set of one-atom observables {1 k , σ x k , σ y k , n k } and calculate their respective averages 1, S ≡ {1, S x , S y , n} according to (·) = tr { ρ(·)}. Applying the QME, assuming spatial uniformity and factorising all quadratic expectations yields the closed set of dynamical equations with V = 2 p V kp the mean-field interaction energy. Stationary regime.-Introducing the effective parameters allows us to formulate the problem in a concise way. We can eliminate S x and S y from the stationary solutions of (3), thus obtaining an algebraic equation for the stationary average number of excitations n, This expression is a cubic real polynomial in n and therefore the number of real roots may vary from 1 to 3 depending on the specific values taken by (a, b, c) within the physically allowed space {a ≥ 2, b ≥ 0}. In Fig. 1 we report the corresponding phase diagram in the a − b plane for different choices of c. The stable phase of the system corresponds to the parameter domain displaying only one acceptable solution. Complementary to this domain is the bistable regime [36,38,39], with Eq. (5) featuring three solutions, only two of which are stable. The boundaries between stable and bimodal regimes are the spinodal lines, where at least two solutions coincide. For any value of c, the spinodal lines coalesce into a critical point identified by a c = −9/(8c) and b c = −27/(8c 3 ), which corresponds to having three coincident real solutions for Eq. (5). This point moves along the curve b = (4a/3) 3 shown in Fig. 1 and lies within the aforementioned physical parameter space only when −9/16 ≤ c ≤ 0. This results in a constraint on the physical parameters: a minimal value b min = 512/27. From them, one can work out the threshold value V min = 4(Γ + K) below which the transition cannot be found by just varying the laser parameters Ω, ∆.
We investigate now the universal features near the critical point. To this end we expand Eq. (5) to leading order in a perturbation of the parameters around their critical values: a = a c + δa and b = b c + δb. We then study the corresponding variation of the stable solutions n st = n c + δn = −2c/3 + δn and identify a special direction δb = (−9/c 2 )δa [in the following referred to as symmetry line, see Fig. 1(a)] along which the solution is invariant under the transformation δn → −δn (with a more complicated one holding for S x and S y ). Thus, a Z 2 symmetry for the stationary value of the excitation density n emerges, which is spontaneously broken FIG. 2. (Color online) Stylized sketch of the mean-field potential V, which gives rise to the dynamics governed by Eqs.
(3), in the S x − n plane. The transformation from the stationary basis of observables {S x , S y , n} to the dynamical one ( S x , S y , n ), where the Z2 symmetry becomes manifest also in the dynamic structure, is in general nonlinear. On the right we show stylized sections of this potential obtained by moving along the directions n , S x and S y respectively. Crucially, the double-well structure (responsible for the "model A" physics) is only felt by the critical n , whereas along the other "massive" directions the system only probes single quadratic wells which play no role in the transition.
in the bistable phase [see Fig. 1(b.1)]. When approaching the critical point along the symmetry line we find δn ∼ (−δa) 1/2 . For any other direction [e.g., the red dashed line in Fig. 1(a)] the system does not switch phases when crossing the critical point. The corresponding behavior, portrayed in Fig. 1(b.2), is described by |δn| ∼ |δa| 1/3 . We can thus conclude that this transition belongs to the (static) Ising universality class with order parameter δn = n − n c : In fact, the magnetization m of an Ising model, as a function of the temperature T , the critical temperature T c and the magnetic field h is known to obey m(T, h = 0) ∼ |T − T c | β and m(T c , h) ∼ h 1/δ , with mean-field exponents β = 1/2 and δ = 3 [3,49]. In analogy, we associate the symmetry line (b.1) to the thermal direction and any deviation from it to the presence of a magnetic field which explicitly breaks the Z 2 symmetry. Finally, a generic choice of the parameters will lead to probing the spinodal behavior shown in Fig. 1(b.3), which has indeed been highlighted in previous theoretical and experimental studies [30,40].
Dynamical vs. static order parameter.-The discussion so far is incomplete as it does not include the dynamical aspects of the system. As a first step, we perform an analysis of the stability of the stationary points. In their neighborhood, we expand the r.h.s. of Eqs. (3) to linear order in the deviations (e.g., δn = n − n st ), which obey the differential equation δ˙ S = M δ S; the eigenvalues of the stability matrix M constitute the rates of approach or escape from the stationary point. Whenever the solution is unique, it is stable as well. When three solutions are present, the two extremal ones are stable, while the one in the middle is unstable, cf. [39].
Here, however, we focus on the universal properties that emerge near the critical point and the spinodal lines, where null eigenvalues appear. The latter are related to (leading) algebraic decays δn ∼ t −1/ζ towards stationarity. The corresponding exponent is ζ = 1 on the spinodal lines and ζ = 2 at the critical point. In the latter case, an algebraic law of the form t −β/(νz) is expected on the basis of scaling arguments, with z being the dynamical critical exponent. The determination of the static universality class (Ising) provides us with the mean-field exponents β = ν = 1/2. Thus, we conclude that z = 2, describing the dynamics of a diffusive system. In addition to the null direction, the stability matrix displays two non-vanishing ("massive") eigenvalues, which identify directions that are not involved in the critical physics (see qualitative sketch in Fig. 2). Thus, the effective order parameter for the low frequency (i.e., long time) dynamics has only one component, and is described by the real variable δn . A discrete Z 2 invariance of the equations of motion under the transformation δn → −δn emerges in a specific direction (i.e., the symmetry line in Fig. 1) emanating from the critical point. This can be verified beyond the linear analysis: Applying a quadratic transformation δΣ i = R ij δS j + Q ijk δS j δS k , we find that along the previously found null direction, the quadratic term vanishes as well. The absence of any apparent conservation law strongly suggests that the dynamics of the system at hand belongs to the (one component) "model A" universality class. We remark that this is similar to the critical point in the driven open Dicke model [10][11][12][13], which however constitutes a zero dimensional model where the mean-field exponents are exact. Consistently with this picture, along the symmetry line the equation of motion reads δṅ ∝ (δn ) 3 at leading order. In contrast, along the spinodal lines where this symmetry is not present, the leading-order equation is of the form δṅ ∝ (δn ) 2 and, consequently, gives rise to an exponent ζ = 1. We remark that the emergent symmetry introduced above lies not among those identified in Ref. [39] (i.e., {∆, V, S x } → {−∆, −V, −S x } and {Ω, S x , S y } → {−Ω, −S x , −S y }), which are unbroken in both phases. Metastable dynamics and connection to experiments.-We now connect these findings to recent experiments [29,30] that have investigated the dynamics of dissipative Rydberg gases. The work presented in Ref. [29] has explored the phase diagram in the Ω − ∆ plane, shown in Fig. 3(a). Another experiment [30] has demonstrated a bistable behavior similar to the one presented in Fig. 1(b.3). Moreover, a power-law behavior of the relaxation time close to a "critical value" of the excitation laser strength was reported.
In order to gain some intuitive insight on the origin of these phenomena, we exploit our knowledge of the universality class and introduce a phenomenological mean-field potential V(n ) = αn − β(n ) 2 + γ(n ) 4 which reflects the profile reported in the topmost panel on the r.h.s. of  1(a)]. We furthermore display the qualitative structure of the mean-field potential V(n ) for parameters corresponding to inside ( ), on ( ) and outside ( ) the spinodal lines (see text). In panel (b) we show an example for the relaxation of the excitation density n (from an initial value n0) towards the stationary value nst, for parameters near one of these boundaries. Here we observe a metastable plateau whose lifetime τ is determined by the first crossing time of the midpoint n = (nst + n0)/2. Panel (c, d) display a power-law divergence of τ as a function of the reduced Rabi frequency ω = (Ω − Ωc)/Ωc varied along the red and purple curves in panel (a) [see for comparison the experimental data shown in Fig. 4 of Ref. [30]]. The power changes from θ = 1/2 in the bistable region [diamonds in panels (a,c,d)] to θ = 2/3 when intersecting the critical point [disks in panels (a,c,d)]. Fig. 2. The corresponding mean-field dynamics is given byṅ = −∂ n V(n ). For fixed β > 0 and γ > 0, this equation portrays a stable (one minimum) to bistable (two minima) transition at a threshold value α c = (2β/3) 3/2 γ. The experiments mentioned above are performed such that initially no excited atoms are present and subsequently the excitation laser is switched on at given values of ∆ and Ω. Within our qualitative framework, in the bistable phase (α < α c ) this may lead to a fast relaxation towards the nearest local minimum of V(n ), which is not necessarily the global one [see ( ) in Fig. 3(a)]. Accounting for fluctuations (beyond mean-field) may in general introduce an additional time scale beyond which this picture is no longer valid and a different physics emerges. However, the agreement between our predictions and experimental observations highlighted below suggests that these features are quite robust in three dimensions and that current experiments indeed probe this "short time physics".
When α = α c , i.e., on a spinodal line, one of the minima becomes an inflection point [see ( ) in Fig. 3(a)]. For α α c [see ( ) in Fig. 3(a)], in the proximity of the disappeared minimum one can identify a region of vanishing slope of the potential, which leads to a characteristic slow dynamics in a flat landscape. This is reflected in the evolution of observables by the appearance of long-lived plateaus (see, e.g., Fig. 3(b) and Fig. (4) in Ref. [30]) whose lifetime τ diverges when approaching α c . By analytically solving the phenomenological equations of motion one obtains τ ∼ (α − α c ) −θ with θ = 1/2, which agrees with the experimental estimate θ = 0.53 ± 0.10 of Ref. [30]. If, instead of a spinodal line, one crosses the critical point (i.e., β = 0 in V (n )), a different exponent θ = 2/3 is found. It should be possible to test this prediction in the experimental setting of Ref. [30].
Note, that instead of n the standard experimental ob-servable is the excitation density n. Nonetheless, critical scaling behaviors characterize the latter as well: in fact, although n is a (non-linear) combination of n and the two "massive" variables, the latter decay exponentially fast (cfr. Fig. 2). Hence, on long time scales the dynamics is determined by the n component. This and the results of the phenomenological model are confirmed by the numerical solution of the full dynamical equations (3): Mimicking the experimental procedure [50], i.e., computing τ for different values of Ω in the proximity of the spinodal lines while keeping V, ∆, Γ, K fixed [see Fig. 3 Conclusions and outlook.-We have found strong evidence for the non-equilibrium dynamics of the dissipative Rydberg gas being governed by the "model A" universality class. This is an extensively studied class [4,6,[51][52][53][54] and in particular it is known that the lower critical dimension is two. This would exclude the presence of a phase transition in dimension one -a question that was raised by the authors of Ref. [40] who have found numerical evidence for long-ranged correlations.
We moreover observed the emergence of a metastable regime in close proximity to the spinodal lines, whose lifetime τ diverges algebraically, in agreement with recent experimental results. Surprisingly enough, our meanfield approach is quantitatively accurate in determining the exponent of this power-law. The accord found between experimental, mean-field and phenomenological approaches suggests that this feature, even far from the critical point, is universal, i.e., it is dictated by the largescale features of the system, rather than by its specific microscopic description. This hints at the presence of a more subtle transition taking place at the spinodal lines. The investigation of this phenomenon -in particular in different dimensions -constitutes a matter of future experimental and theoretical investigation.