Non-equilibrium universality in the dynamics of dissipative cold atomic gases
The theory of continuous phase transitions predicts the universal collective properties of a physical system near a critical point, which for instance manifest in characteristic power-law behaviours of physical observables. The well-established concept at or near equilibrium, universality, can also characterize the physics of systems out of equilibrium. The most fundamental instance of a genuine non-equilibrium phase transition is the directed percolation universality class, where a system switches from an absorbing inactive to a fluctuating active phase. Despite being known for several decades it has been challenging to find experimental systems that manifest this transition. Here we show theoretically that signatures of the directed percolation universality class can be observed in an atomic system with long range interactions. Moreover, we demonstrate that even mesoscopic ensembles — which are currently studied experimentally — are sufficient to observe traces of this non-equilibrium phase transition in one, two and three dimensions.
In strongly interacting systems, long range correlations can lead to the emergence of collective behaviours which can be distinctively different from the single-body physics governing the constituents. A remarkable example can be found within the theory of continuous phase transitions [1, 2, 3], which indeed predicts characteristic singular behaviours of macroscopic observables, in contrast to the usually smooth microscopic description. At a so-called critical point, the correlation length of the system diverges and the overall behaviour is fundamentally determined by certain properties — such as dimensionality, range of interactions and symmetries — that do not depend on the specific scale. All systems sharing these few coarse-grained features display the same qualitative macroscopic physics  and form a universality class, which is in turn characterized by the corresponding critical exponents and ratios [1, 2, 3, 5].
Phase transitions have been extensively studied in equilibrium — both in the classical regime [6, 2, 3, 7] and in the quantum one [8, 9] — within the general framework of statistical mechanics, based on the properties of a few macroscopic thermodynamic potentials. Out of equilibrium, the lack of an equivalent unified picture has prevented reaching a similar degree of insight. Despite this, critical phenomena are known to occur out of equilibrium as well, e.g., in the statistics of the work [10, 11, 12], the temporal evolution [13, 14, 15, 16] and the spectral properties  of quenched quantum systems. In particular, phase transitions can take place in the properties of the steady state; this typically leads to an enrichment of the stationary phase diagram which depends upon the coarse-grained aspects of the dynamics, such as symmetries and conservation laws (see, e.g., [18, 19]). Equilibrium conditions, for example, are specifically related to the microreversibility symmetry [20, 21, 22, 23]. Systems in which the latter is not recovered at large length and time scales undergo genuine non-equilibrium phase transitions [24, 25, 26].
One of the most fundamental and well-studied universality classes of this kind is directed percolation (DP). A DP-like transition is defined as a continuous transition between a unique absorbing state and a fluctuating one observed for a positive, one-component order parameter. An intuitive description can be provided in terms of the “contact process”, i.e., a classical stochastic process on a chain of Ising spins in which an up spin (or excitation) can always flip down (self-destruction), whereas a down spin can only flip up if another up spin is present in its neighbourhood (offspring production), see Fig. 2. Due to the fact that the latter of these processes require a nearby up spin to take place, the “all-down” state constitutes not only a stationary state of the dynamics, but an absorbing one as well (i.e. once the system is in this state it cannot abandon it), being completely free of stochastic fluctuations. When the rate at which self-destruction occurs is much larger than the rate of offspring production, the system invariably ends up in the absorbing state, no matter what the initial conditions were. At a critical ratio between the rates, the system switches continuously from the “all-down” steady state to a fluctuating active one with a finite density of excitations. Despite its simplicity and robustness , it has been very difficult to identify clear signs of DP universality in physical systems . Only recently an experiment focussing on two distinct topological phases of nematic liquid crystals provided the first clean realization of DP in a two-dimensional physical system [28, 29].
In this work we show that strongly interacting ensembles of highly excited (so-called Rydberg) atoms [30, 31, 32] feature a dynamics which — in a certain limit — is governed by the elementary rules of a DP process . Beyond revealing insights into the out-of-equilibrium behaviour of this currently much studied system, our work highlights an alternative approach for the experimental exploration of DP universal features not only in two but in one and three dimensions as well. Remarkably, our results show that even for relatively small (mesoscopic) system sizes which compare to those currently studied in experiment [33, 34, 35, 36], clear signatures of DP are observable.
This paper is structured as follows. In Section 2 we briefly introduce the open quantum many-body system under consideration and focus on a parameter regime where its dynamics can be described via a classical rate equation. Section 3 shows how the parameters of the system can be tuned appropriately so that this dynamics is very close to that of a DP process, which is then further tested at the mean-field level in Section 4. Numerical simulations of mesoscopic systems in one and two dimensions indicate that indeed strong signatures of the DP universality class can be observed in this long range interacting system (Section 5). Our concluding remarks can be found in Section 6.
2 The setup
The specific setup we consider (figure 1) is an ensemble of atoms trapped in a lattice with spacing . Note that while we focus here on the one- and two-dimensional cases, the arguments below can be analogously applied to three dimensions as well. We describe the internal level structure of the atoms with two relevant levels: the ground state and a highly excited (Rydberg) state [31, 32], coupled by a laser with Rabi frequency and detuning . When two atoms (at positions and ) are in the Rydberg state they experience a long range interaction of strength which is parameterized as . We will focus here on Rydberg -states that interact via van der Waals (vdW) forces (). The corresponding Hamiltonian of this system, expressed in a rotating frame and after discarding fast-oscillating terms (i.e. in the Rotating Wave Approximation) yields
where , and . Moreover, we consider that the system is subject to two sources of Markovian noise (or dissipation), one that leads to spontaneous radiative decay from the Rydberg state to the ground state at a rate and the other causing dephasing of atomic superposition states at a rate (due to uncorrelated noise) [37, 38]. These two dissipation processes are described via the Lindblad superoperators
In the limit of strong dephasing, i.e. , a perturbative treatment of the full dissipative quantum dynamics allows for the description of the dynamics of this system by means of a classical rate equation for the probability vector whose components are the statistical weights of the classical spin configurations (e.g., ). The derivation of said equation is discussed in a more general framework in Refs. [37, 38, 39, 40], to which we refer the interested reader. Here we report the one corresponding to the setup outlined above:
This depicts a classical stochastic process in which the -th spin flips up with rate and down with rate (see figure 1). Here, the rate
which depends on the state of all spins but the -th one, is analogous to those that appear in dynamically-facilitated models of glasses [41, 42]. In A we compare the full quantum dissipative dynamics and the one described by (4) for small systems, thus providing further evidence of the validity of this approach in the current case.
3 Emergent DP process
It has been conjectured that the defining conditions for the emergence of DP universality are: (i) a local dynamics with a unique absorbing state, (ii) a continuous phase transition with a positive, one-component order parameter and (iii) absence of additional symmetries [43, 25]. In a spin chain, the simplest setup that meets these conditions consists of two fundamental processes: flipping down spins (self-destruction) and flipping them up provided there is an up spin in the neighbourhood (offspring production) (see two topmost rows in figure 2). On the other hand, processes that can create isolated excitations (self-activation) destroy the absorbing property of the “all-down” state. Hence, they constitute a relevant perturbation away from the DP critical behaviour (fourth row in figure 2). The presence of other processes, such as the flipping down of a spin provided there is another up spin in the neighbourhood (coagulation, third row in figure 2) do not modify the critical properties of the system.
Let us now study the dynamical processes in our physical system and the corresponding underlying rates in more detail. While the self-destruction of excitations is simply provided by the radiative decay (rate ), the facilitated excitation and de-excitation of the -th spin (offspring production and coagulation) occur at rates and , respectively (see equation (5)). These processes can be favoured by tuning the laser parameters adequately. In particular, given an existent excitation, fixing the detuning such that it cancels exactly the nearest neighbour interaction , i.e. , effectively brings on resonance the excitation and de-excitation of the atoms sitting next to it [44, 45], i.e., it maximizes the value of .
Clearly, due to the long-distance tails of the potential , the presence of an excitation not only affects its neighbours, but slightly modifies the rates of spins lying far from it as well. As we argue in B, these changes are akin to introducing the possibility of Lévy flights in the offspring production. However, the van der Waals potential decays fast enough for the resulting corrections not to affect the critical properties of the DP dynamics [46, 47]. We also refer to B for the derivation of the main dynamical processes of the Rydberg system. Here we report them in figure 2 along with the corresponding rates, re-expressed for simplicity in terms of the dimensionless parameters and and measured in units of the decay rate . Focussing on the fourth row, we note that the parameter can be increased (limit of large detuning ) to reduce the impact of self-activation on the system. By comparing the first two rows, instead, we note that within the same regime measures the relative strength of the offspring production and self-destruction processes and hence drives the DP phase transition. However, here the presence of long range interactions imposes an additional constraint on the actual emergence of DP universality, as it affects the growth of clusters of excitations. For illustration purposes, let us consider only the presence of next-to-nearest neighbour interaction of strength with : While the offspring production rate (i.e. ) is given by , processes of the type occur at a rate . This can be seen calculating the corresponding energy difference appearing at the denominator of (5) for such a process. The latter rate, always smaller than , becomes extremely small in the limit of large , eventually hindering the growth of clusters. To prevent this one needs to impose the constraint (). In 1D and for a vdW potential we have which allows one to make reasonably large. In contrast, in higher dimensions we find the more restrictive , as next-to-nearest neighbours lie at a distance from one another. This restriction can be softened, as we discuss in C, by modulating the interaction between the atoms via external microwave fields.
4 Mean-field analysis
As a simple consistency check, our first step towards understanding the critical properties of this non-equilibrium system is a mean-field analysis. The time evolution of the average local density of excitations , after factorizing all spatial correlations (), is governed by the equation . Employing the representation (8) and assuming the density to be uniform, i.e., , one finds the mean-field equation
where for convenience we have defined , and is the coordination number (i.e., number of nearest neighbours) of the lattice. Here, one can recognize in the first two terms spontaneous decay (self-destruction) and processes occurring where there is a single excitation around (offspring production and coagulation) (three first rows in figure 2). The term, , represents the self-activation process, which is the one driving the system away from criticality. This can be understood by noticing that it is the only surviving term when , hence allowing the system to leave the “all-down” state.
The stationary mean-field solution as a function of and is shown in figure 3(a) for a 1D chain (). For finite we observe a smooth crossover from a low-density region to a high-density one. As we increase [see figure 3(b)], criticality emerges and one observes a sharp transition at . In this regime, the last terms of (4) affect the evolution only on long timescales — dictated for a single spin by . The dynamics at times shorter than those is thus effectively captured by
which is similar to the mean-field equations of other known DP processes . In fact, it features a DP-like phase transition: while for the only stationary state is the absorbing one with density , for the latter becomes unstable and a new, stable one with a non-vanishing density appears, identifying as a critical point.
The mean-field critical exponents can also be extracted: the approach to stationarity can be studied at the critical point by introducing a small perturbation from the all-down state and linearising (7), which yields and therefore agrees with the prediction for DP with . Focusing on the properties of the steady state as a function of in the proximity of the critical point () one finds, at leading order, , which is again compatible with the DP prediction with .
5 Numerical results
We have performed Monte Carlo dynamic simulations of small 1D and 2D systems with open boundary conditions for (4) (see figure 4). While (4) can be simulated for very large lattices, we focus here in system sizes and boundary conditions relevant to current Rydberg experiments. We have used and the initial state is a random classical spin configuration with fixed global density .
We first focus on the 1D case, where we have accounted for both vdW () and nearest neighbour (NN) interactions () in a system of atoms. Signatures of the transition are most effectively highlighted in the dynamics. By increasing (e.g., by acting on the Rabi frequency ), one encounters a qualitative change of the curves from an exponential decay towards (lower curves) to an exponential relaxation to a finite value (upper curves) with a characteristic bending down and up, respectively [see figures 4(a) and (d)]. The intermediate, critical regime () displays instead an algebraic decay in reasonable agreement with the DP one  ( with is shown for comparison). We then analyze the behaviour of the stationary density in figures 4(b) and (e). This yields again an algebraic law close to which resembles the DP one [the expected value is again shown for comparison in panels (c) and (f)]. As a further check, we have assumed scaling with the DP critical exponents  and plotted as a function of for different values of close to the critical point as an inset in panels (a) and (d). This shows that the curves tend to collapse, as expected in the presence of a continuous phase transition. This collapse is not present at short times, since we have not accounted for initial slip effects, and at long times , where the relevant processes perturbing DP become non-negligible. As discussed earlier, the main difference between the dynamics resulting from the two interactions (NN and vdW) is the rate of growth from a two- to a three-excitation cluster () which, with our choice of , is times smaller in the vdW case. Correspondingly, we expect a shift in the critical point by approximately the same factor, as can be indeed observed by comparing figures 4(b) and 4(e). Note that the profiles of these curves are analogously stretched by the same factor, making the nearest-neighbour crossover look sharper than the vdW one. We have verified, however, that the two curves plotted as functions of the rescaled parameter look reasonably similar. Moreover, as the value of is increased in the vdW case, the time at which the self-activation mechanism becomes relevant (proportional to ) is shorter. Despite this, one can still observe a collapse of the curves in the intermediate regime as discussed above.
The dynamics in 2D with vdW interactions is actually more complicated: at the rather small sizes accessible in experiments and numerically investigated here, one observes the emergence of a spurious one-dimensional, rather than two-dimensional critical behaviour. This can be intuitively understood by analyzing once again the growth from a two- to a three-excitation cluster: The two neighbouring excitations identify a direction (say, along a row or a column of the square lattice). Further growth along this direction occurs at a rate . In contrast, the rate for growing in the orthogonal direction is , since in this case the next-to-nearest excitation lies at distance from the first excitation (i.e., along the diagonal of a square cell). Hence, it is much more likely to grow linear structures than it is to percolate in both directions. At the same time, this anisotropy is not a geometrical property of the system, but it is stochastically determined for a cluster by the first excitation it forms. In small systems, clusters growing vertically and horizontally hinder the respective growths, effectively reducing the available amount of space and thus, in a sense, amplifying finite-size effects. In order to recover proper 2D directed percolation behaviour in this setting one would have to increase the system size beyond the current experimental possibilities. We remark that, although lowering the detuning reduces the effective anisotropy, it also increases the self-activation rate and is therefore not a valid strategy. A seemingly more viable approach is to shape the atomic interaction potential. In C we discuss a way to modulate the potential and induce a faster decay of the tail at large distances. This allows to employ the nearest-neighbour case as an approximate description in higher dimensions. The 2D data shown in figures 4(g), (h) and (i) are therefore obtained for a system of atoms considering NN interactions only. They show that the exponents and fit again fairly well the data from the numerical simulations. Accordingly, the collapse also displays the emergence of two distinct master curves.
We have shown that interacting gases of Rydberg atoms feature a microscopic dynamics that leads to an emergent out-of-equilibrium behaviour displaying specific features characteristic of DP. Surprisingly this is the case already for rather small mesoscopic system sizes which are currently studied experimentally. This renders these cold atomic systems into a viable candidate for the experimental observation of the simple — yet intriguingly hard to observe — DP non-equilibrium universality class. Such experimental studies could furthermore become a platform for a systematic exploration of the role of quantum effects on the dynamics (a very challenging task from the theoretical point of view), and shed light on the role of the assumed ordered configuration (lattice) by considering continuous gases.
Appendix A Effectively classical dynamics
Here, we provide evidence of the validity of the effectively classical equation of motion (4) in the specific parameter regime considered in this work. A numerically exact solution of the full quantum Master equation underlying this problem is only possible for small systems sizes. Here we use atoms in a 1D chain with NN interactions. Initially there is a single excitation in the center of the lattice and the corresponding dissipative quantum dynamics is simulated with Quantum Jump Monte Carlo. A comparison to the dynamics obtained from (4) is displayed in figure 5(a) for , and and . Indeed an excellent agreement between the full quantum and the approximate classical dynamics is found .
Appendix B Adjusting the rates of dynamical processes in the Rydberg lattice
We introduce the projectors (, with being the coordination number of the lattice) on the subspaces where exactly nearest neighbours of the -th atom are excited. As these projectors commute with , we can decompose with , which yields
where denotes all sites included in the neighbourhood of site . Consequently, restricts the sum on atoms lying further than nearest neighbours from site . Since the interaction is quickly decaying with the distance we have which allows us to approximate . This represents the rate at which a spin flips up if it has a single excited neighbour (second row in figure 2). It is important to remark that this approximation, made at the microscopic scale, is expected not to influence the critical behaviour emerging on mesoscopic scales: in fact, the leading correction can be expressed as
Exploiting the fact that the operator inside the brackets is diagonal in the basis, we see that for every possible classical state it can be bounded from above according to
which is reminiscent of a term causing offspring production at arbitrary distances from the parent excitation with an algebraically-decreasing rate. Such corrections have been investigated in the past [46, 47] and have been found to modify the scaling behaviour only if the exponent of the aforementioned algebraic decay lies below a certain threshold. The latter was estimated to be in 1D and in 2D (see, e.g., equation (42) in ). For van der Waals interactions, the exponent appearing in equation (10) is and thus definitely higher than these thresholds and hence the dynamics is deep within the “ordinary” DP regime.
Correspondingly, provides an estimate of the rate at which spins flip up when far away from excitations. These processes can create isolated excitations, thereby destroying the absorbing property of the “all-down” state. They therefore constitute a relevant perturbation away from the DP critical region (fourth row in figure 2). Using the fact that the Rydberg interaction is quickly decaying, i.e. , the rate of these processes can be estimated as
Hence, the magnitude of the “DP-breaking” processes is strongly suppressed at large laser detuning: when one finds that where denotes the largest eigenvalue of in modulus and, in the simple case of (11), it can be straightforwardly obtained by replacing the projectors with .
Appendix C Potential shaping
The Rydberg states we are considering here display long range interactions of the vdW type. As we have discussed in the main text, in one dimension they allow for the observation of DP in 1D. In 2D and 3D, however, the presence of long range interactions induces critical properties different from the 2D and 3D DP ones due to the preference of the clusters to grow linearly. In order to observe DP in Rydberg gases in dimensions two and three, it is desirable to have a potential closer to a NN interaction.
This can be achieved through potential shaping via microwave dressing of Rydberg states [48, 49, 50]. To illustrate this we consider a Rydberg -state (orbital quantum number ) with principal quantum number coupled to a state via a linearly polarized microwave field with Rabi frequency and detuning [see figure 5(b)]. This results in a “hybridization” of the respective interactions, which consequently modifies the bare purely algebraic decay of vdW potential. For example, in figure 5(b) we show the result of coupling the and states of Rubidium with Rabi frequency MHz and detuning MHz. Considering the lattice constant to be , the ratio between nearest (distance ) and next-to-nearest neighbour ( in 2D and 3D) interaction increases from (in the vdW case) to .
Note that for the coherent superposition between the two states ( and ) to be maintained throughout the experiment the Rabi frequency of the microwave field must be the largest energy scale in the problem, i.e. .
-  Ma S K 1976 Modern theory of critical phenomena Frontiers in Physics (Reading, MA: W. A. Benjamin)
-  Goldenfeld N 1992 Lectures on Phase Transitions and the Renormalization Group Frontiers in Physics (Westview Press)
-  Zinn-Justin J 2002 Quantum Field Theory and Critical Phenomena International Series of Monographs on Physics (Clarendon Press)
-  Back C H, Würsch C, Vaterlaus A, Ramsperger U, Maier U and Pescia D 1995 Nature 378 597
-  Pelissetto A and Vicari E 2002 Phys. Rep. 368(6) 549
-  Huang K 2009 Introduction to Statistical Physics 2nd ed (Chapman and Hall/CRC)
-  Bellac M L 1991 Quantum and Statistical Field Theory (Oxford Science Publications)
-  Sachdev S 1999 Quantum Phase Transitions (Cambridge University Press)
-  Mussardo G 2009 Statistical Field Theory Oxford Graduate Texts (Oxford University Press)
-  Gambassi A and Silva A 2012 Phys. Rev. Lett. 109(25) 250602
-  Sotiriadis S, Gambassi A and Silva A 2013 Phys. Rev. E 87(5) 052129
-  Gambassi A and Silva A 2011 preprint arXiv:1106.2671
-  Heyl M, Polkovnikov A and Kehrein S 2013 Phys. Rev. Lett. 110(13) 135704
-  Lesanovsky I, van Horssen M, Guţă M and Garrahan J P 2013 Phys. Rev. Lett. 110(15) 150401
-  Karrasch C and Schuricht D 2013 Phys. Rev. B 87(19) 195104
-  Buchhold M and Diehl S 2014 preprint arXiv:1404.3740
-  Torre E G D, Diehl S, Lukin M D, Sachdev S and Strack P 2013 Phys. Rev. A 87(2) 023831
-  Hohenberg P C and Halperin B I 1977 Rev. Mod. Phys. 49(3) 435
-  Sieberer L M, Huber S D, Altman E and Diehl S 2013 Phys. Rev. Lett. 110(19) 195301
-  Agarwal G 1973 Z. Phys. 258 409–422 ISSN 0044-3328
-  Chetrite R and Mallick K 2012 J. Stat. Phys. 148 480–501 ISSN 0022-4715
-  Aron C, Biroli G and Cugliandolo L F 2010 J. Stat. Mech. 2010 P11018
-  Sieberer L M, Chiocchetta A, Gambassi A, Täuber U C and Diehl S 2015 preprint arXiv:1505.00912
-  Henkel M, Hinrichsen H and Lübeck S 2009 Non-Equilibrium Phase Transitions (Theoretical and Mathematical Physics vol 1) (Springer)
-  Hinrichsen H 2000 Adv. Phys. 49 815
-  Ódor G 2004 Rev. Mod. Phys. 76(3) 663–724
-  Hinrichsen H 2000 Brazilian Journal of Physics 30 69 – 82 ISSN 0103-9733
-  Takeuchi K A, Kuroda M, Chaté H and Sano M 2007 Phys. Rev. Lett. 99(23) 234503
-  Takeuchi K A, Kuroda M, Chaté H and Sano M 2009 Phys. Rev. E 80(5) 051116
-  Bloch I, Dalibard J and Zwerger W 2008 Rev. Mod. Phys. 80(3) 885
-  Löw R, Weimer H, Nipper J, Balewski J B, Butscher B, Büchler H P and Pfau T 2012 J. Phys. B: At. Mol. Opt. Phys. 45 113001
-  Saffman M, Walker T G and Mølmer K 2010 Rev. Mod. Phys. 82(3) 2313–2363
-  Viteau M, Bason M G, Radogostowicz J, Malossi N, Ciampini D, Morsch O and Arimondo E 2011 Phys. Rev. Lett. 107(6) 060402
-  Schauß P, Cheneau M, Endres M, Fukuhara T, Hild S, Omran A, Pohl T, Gross C, Kuhr S and Bloch I 2012 Nature 491 87
-  Schauß P, Zeiher J, Fukuhara T, Hild S, Cheneau M, Macrì T, Pohl T, Bloch I and Gross C 2014 preprint arXiv:1404.0980
-  Barredo D, Labuhn H, Ravets S, Lahaye T, Browaeys A and Adams C S 2014 preprint arXiv:1408.1055
-  Lesanovsky I and Garrahan J P 2013 Phys. Rev. Lett. 111(21) 215305
-  Marcuzzi M, Schick J, Olmos B and Lesanovsky I 2014 J. Phys. A: Math. Theor. 47 482001
-  Ates C, Pohl T, Pattard T and Rost J M 2006 J. Phys. B: At. Mol. Opt. Phys. 39 L233–L239
-  Ates C, Pohl T, Pattard T and Rost J M 2007 Phys. Rev. A 76(1) 013413
-  Ritort F and Sollich P 2003 Adv. Phys. 52 219
-  Chandler D and Garrahan J P 2010 Annual Review of Physical Chemistry 61 191–217
-  Janssen H K 1981 Z. Phys. B 42 151–154 ISSN 0722-3277
-  Lesanovsky I and Garrahan J P 2014 Phys. Rev. A 90(1) 011603
-  Urvoy A, Ripka F, Lesanovsky I, Booth D, Shaffer J, Pfau T and Löw R 2015 Phys. Rev. Lett. 114(20) 203002
-  Janssen H K, Oerding K, van Wijland F and Hilhorst H J 1999 Eur. Phys. J. B 7(1) 137–145
-  Hinrichsen H 2007 J. Stat. Mech. 2007(7) P07006
-  Petrosyan D and Mølmer K 2014 Phys. Rev. Lett. 113(12) 123003
-  Kiffner M, Li W and Jaksch D 2013 Phys. Rev. Lett. 111(23) 233003
-  Kiffner M, Li W and Jaksch D 2013 Phys. Rev. Lett. 110(17) 170402