From classical to quantum non-equilibrium dynamics of Rydberg excitations in optical lattices

The glass phase and its quantum analog are prominent challenges of current non-equilibrium statistical mechanics and condensed matter physics. As a model system to study the transition from classical to quantum glassy dynamics, we propose a setup of laser driven three-level atoms trapped in an optical lattice. Tuning the strength of the laser driving to the intermediate level allows one to study the transition from a classical Kinetically Constrained Model to the coherent regime. For strong driving, Rydberg excitations evolve analogously to defects in the One-Spin Facilitated Model, a minimal model known to exhibit glassy dynamics. In our setup, the constraints result from the interplay between Rydberg interactions and the laser detuning from the Rydberg state. The emerging heterogeneous relaxation timescales are tuneable over several orders of magnitudes. In the opposite limit of weak driving of the intermediate level, we find an effective cluster model which describes the dynamics in a reduced subspace of the allowed number and positions of Rydberg excitations. This subspace is uniquely determined by the initial state and is characterized by a fixed number of clusters of Rydberg excitations. In addition, we investigate the influence of random fields on the classical relaxation. We find that the glassy dynamics can relax faster in the presence of weak random fields.


Introduction
When disordered materials are quenched to low temperatures, the dynamics can become non-ergodic and the system may fall out of equilibrium [1,2,3,4,5]. Understanding this phenomenon, known as the glass transition, is one of the open questions in current statistical mechanics. While a complete theory which satisfactorily explains all features of glass forming materials is missing, spin models have been proposed to understand particular aspects of the glass problem. For example, static properties of glasses are captured by Ising models with random interactions, the so-called spin glasses [6,7,8]. These models have recently gained renewed interest in the context of both quantum adiabatic optimization [9,10] and, in presence of external random fields, of the many-body localization phase transition [11,12,13,14]. Spin models which study dynamical properties of glasses are so-called Kinetically Constrained Models [15] (KCMs). In these KCMs, separation of relaxation timescales as well as space-time dynamical heterogeneities [5,16] result from transitions (spin flips) which are allowed or forbidden depending on the state of neighboring spins. Quantum effects in KCMs have been recently studied in Refs. [17,18,19,20,21] in ultracold Rydberg atoms. Thanks to their strong and long-range interactions, as well as to their long lifetimes [22,23,24], Rydberg atoms offer a promising framework to study quantum spin models in-and out-of-equilibrium [25,26,27,28,29,30,31,32,33,34,35,36]. With the considerable recent experimental development in atomic, molecular and optical physics, Rydberg atoms can in fact be trapped, laser-excited and detected with single particle resolution in state-of-the-art experiments [37,38,39,40,41,42,43,44].
In this work we propose a model system based on V-shaped three-level atoms, that allows one to study the transition between a classical dissipation-induced glass model and a fully coherent quantum system. Atoms are resonantly driven from the electronic ground-state |g to a short-lived low-lying excited state |e and far-off-resonantly driven to a metastable Rydberg state |r (see Fig. 1). Tuning the driving strength Ω e of the |g ↔ |e transition, allows one to study the crossover from classical to quantum dynamics of Rydberg excitations. For large Ω e , the relaxation of Rydberg excitations resembles the one of defects in a classical KCM, the One-Spin Facilitated Model [15,18] (1-SFM). Here, kinetic constraints are a consequence of the interplay between strong interactions of Rydberg-excited atoms [45] and the laser detuning from |r . In contrast, for vanishing Ω e , coherences of the |g ↔ |r transition dominate the dynamics within the lifetime of the Rydberg states [37]. In the case of large laser detuning from the Rydberg state the system resembles a transverse Ising model for a particular choice of the detuning [46,47]. We present an effective simple model to describe the evolution in this regime which is based on a fixed number of Rydberg clusters with variable number of Rydberg excitations.
In addition to the study of the transition from a glass model to a fully coherent integrable model, the proposed framework opens up the possibility to investigate the interplay between internal dynamical constraints and external disorder arising from random fields. In our system, static random fields can be implemented with sitedependent random detunings from the Rydberg state. Our results indicate that the glassy relaxation of Rydberg excitations toward the steady state can be accelerated by the presence of weak random fields. Figure 1. The setup we have in mind consists of three-level atoms, each one schematized with a ground-state |g , a Rydberg state |r and an intermediate excited state |e . Atoms are coherently laser-driven with Rabi frequencies Ω e (Ω r ) between |g -|e (|g -|r ) and γ e is the decay rate due to spontaneous emission from |e . The antiblockade condition V = ∆ r , i.e. when the laser detuning ∆ r with respect to the |g -|r atomic transition is equal to the interactions between neighbouring Rydberg-excited atoms V (see Subsec. 2.2), ensures that whenever an atom is in |r (defect, blue sphere), an atom next to it, initially in |g (facilitated atom, red sphere), will also be favoured to be excited to |r . An atom in |g far from a defect (non-facilitated atom, black sphere) is blocked, because its |g -|r transition is far-off-resonant (see Subsec. 4.1). The case with additional site-dependent random detunings ∆ ri is studied in Subsec. 4

.3.
The remainder of the paper is organized as follows: In Sec. 2 the level-scheme and the model are introduced; in Sec. 3 we give a brief overview of SFMs. Numerical results from quantum trajectory simulations are presented in Sec. 4. In Subsec. 4.2 we study the coherent regime and the transition to the rate equation limit and in Subsec. 4.3 we analyze the case of finite random detunings from the Rydberg state. The validity of the assumptions and approximations of Sec. 4 are discussed in Sec. 5. Finally, in Sec. 6 we summarize our conclusions.

Model
We consider N atoms individually trapped in a one-dimensional (1d) large spacing optical lattice [37,48,49,50,51] or in a magnetic trap array [52,53]. Atomic motion due to tunnelling between adjacent sites is negligible within experimental timescales. Each atom is approximated as a three-level system in the V-configuration, as schematically depicted in Fig. 1: the electronic ground-state |g is laser driven to (i) a Rydberg state |r with a small Rabi frequency Ω r , and (ii) to an excited state |e with a large Rabi frequency Ω e . The intermediate low-lying excited state |e is rapidly decaying to |g and its short lifetime is given by the inverse of the decay rate γ e . Two atoms in their Rydberg states interact via dipole-dipole V (x) = V /x 3 or van der Waals V (x) = V /x 6 potentials, depending on the presence or absence of an applied external field [23], respectively. We approximate interactions between Rydberg-excited atoms to be nearest-neighbor, i.e. V (1) = V and V (x 2) = 0. Results under this assumption are in qualitative good agreement with simulations of the full model including long-range interactions. Furthermore, we neglect processes of spontaneous emission from the Rydberg state, setting the associated decay rate γ r to zero. Restrictions on the validity of these approximation are discussed in detail in Sec. 5.
The Hamiltonian in the rotating frame (with = 1) is Here, the sum runs over all N lattice sites, while ∆ r and ∆ ri are the homogeneous and site-dependent random laser detunings from |r , respectively. The random detuning can, for example, be induced by a spatially inhomogeneous ac-Stark shift acting on the Rydberg-or ground-state using either a random light pattern or differential lightshifts between atoms in different traps, in close analogy to experiments performed in the context of Anderson localization [54,55].

A single-atom in the V-configuration
Even though a single three-level systems in the V-configuration has been extensively studied in literature, both theoretically [56,57,58,59,60,61,62] and experimentally [63,64], let us briefly discuss its properties. For a more complete and detailed discussion, we refer the reader to Appendix A. The single-atom Hamiltonian of the three-level scheme in Eq. (1), where the Rydberg-state is far-off-resonantly laser driven (∆ r = 0) and no random detuning is From classical to quantum non-equilibrium dynamics of Rydberg excitations in optical lattices5 In the Born-Markov approximation, the master equation describing the dissipative evolution of the reduced system density matrix ρ(t) is in the Lindblad form [65] where c = |g e| is the quantum jump operator. Assuming Ω r Ω 2 e /γ e and ∆ r Ω e , Ω r , γ e (the reason behind this choice is explained in Appendix A), we imagine to measure photons emitted by our three-level atom. A typical photon counting sequence registered by a photodetector would then look like an alternation of time intervals with many photons detected (bright periods, which we label with 0) and no photons detected (dark periods, which we label with 1). The transition probability (rate) in going from 0 to 1 is given by whereas the transition probability from 1 to 0 is The analytical derivation of these formulas and the physical mechanism behind these transitions, known as quantum jumps, are given in Appendix A.

Many-body model
In this Subsection we study the many-body model (∆ r , V = 0) with no random detunings. In analogy to the single-atom case [Eq. (3)], the many-body master equatioṅ with c i = |g e| i the quantum jump operator at site i, describes the driven-dissipative dynamics of the many-body system. Furthermore, let us assume: (i) Ω r Ω 2 e /γ e , as in Subsec. 2.1, for well defined quantum jumps to occur and (ii) ∆ r Ω e , Ω r , γ e and (iii) V = ∆ r , which we call the anti-blockade condition. The many-body transition rates, corresponding to the single-particle rates in Eqs. (4) and (5), can be determined as follows. In a 1d system, due to the assumption V (x 2) = 0, every process at a certain site i can be influenced at most by its two nearest-neighbors, i ± 1. The many-body rates can be calculated from Eqs.  (5), ∆ r appears only in even powers, both off-resonant rates = 0, 2 are equally suppressed with respect to the resonant rate = 1, as highlighted in Tab. 1. ∆ r = 0 ∆ r = ±3γ e ∆ r = ±10γ e Γ 0→1 (∆ r )/γ e 3.00 · 10 −4 8.80 · 10 −6 7.54 · 10 −7 Γ 1→0 (∆ r )/γ e 9.00 · 10 −4 7.14 · 10 −7 5.64 · 10 −9  (5), for ∆ r = 0 (resonant rates) and ∆ r = ±3γ e , ±10γ e (far-off-resonant rates). Other parameters are: Ω e = γ e and Ω r = 0.03γ e . The difference between resonant and far-off-resonant rates is of several order of magnitudes and increases with |∆ r |.
In absence of conditions (i), (ii) and (iii), the strong transition |g ↔ |e might significantly increase the probability Γ 0(1)→1(0) (±∆ r ) of (originally) suppressed processes of excitation (deexcitation) to (from) |r to take place. Since rates differ by several order of magnitudes (see Tab. 1), the dynamics shows a large separation of timescales (cfr. Sec. 4).

Spin Facilitated Models
Spin Facilitated Models (SFMs) are a class of KCMs that exhibit dynamical heterogeneities and slowing down of relaxation timescales [15]. Even though they can be implemented in a variety of lattice geometries and dimensions, from now on we will limit ourselves to 1d. Let us now consider the energy function of non-interacting spins in an external field, introduced by Fredrickson and Andersen in Refs. [67,68], where n i = 0, 1 is a classical spin variable. In our setup, n i = 0 (1) corresponds to the i-th atom being in a bright (dark) period. The energy function in Eq. (7) does not feature a finite temperature phase transition for a field strength K = 0. The equilibrium concentration of defects, which we define as spins in the n i = 1 state, reads where β = 1/(k B T ), k B is the Boltzmann's constant and T is the equilibrium temperature of the system. Furthermore, for every K > 0, lim T →0 d eq = 0, in line with the intuition that at low T the concentration of defects should be small. From now on, without loss of generality, we set k B = 1 and we fix the temperature scale setting K = 1. Despite the simple ground-state of Eq. (7), the non-equilibrium dynamics towards it can be non-trivial. In general, the evolution of a classical system of spins is described by the following rate equation for the probabilities p(n, t) of being in the many-body spin state specified by the vector n = (n 1 , ..., n i , ..., n N ) at time t: Here, Γ(m → n) is the rate of the transition from m to n = m. The equilibrium probabilities π(m, t) satisfy the detailed balance Γ(m → n) π(m, t) = Γ(n → m) π(n, t) with respect to Eq. (7). Under the assumption of local spin-flips and Glauber dynamics [69], the transition rate from (n 1 , ..., n i , ..., n N ) to (n 1 , ..., 1 − n i , ..., n N ) at site i, Γ n i →1−n i , is equal to 1/(1 + e β∆E ), with ∆E = 1 − 2n i the energy difference between the final and initial spin configurations. Since in our case ∆E = ±1, from Eq. (8) single-particle rates can be written as Thus, the asymmetry between the ratio of the rates of creation over destruction of a defect fixes the temperature T at which the system of non-interacting spins thermalizes. (Note that Γ 0→1 → Γ 1→0 for T → ∞.) Up to now, all transitions between different configurations are allowed. The key idea behind SFMs is to introduce restrictions in these transitions. The general rule can be stated as: a spin flip is allowed only if there are enough defects in the neighbourhood that can facilitate the process. In the 1d 1-SFM, these so-called kinetic constraints require that a spin can flip provided at least one of its two neighbors is a defect. This facilitation condition transforms Eq. (10) into The association between rates in Eq. (12) and the many-body rates calculated in Subsec. 2.2 is straightforward: n i−1 = n i+1 = 0 corresponds to = 0 for which no spin flip at site i is allowed, while either {n i−1 = 1, n i+1 = 0} or {n i−1 = 0, n i+1 = 1} corresponds to = 1 and i is facilitated to flip its spin. The only difference arises in the case of {n i−1 = n i+1 = 1}: while rates in the 1-SFM are doubled with respect to the case {n i−1 = 1, n i+1 = 0}, for = 2, rates are suppressed as much as in the case of = 0. This will lead to different steady states which we will discuss in Sec. 4.

Semi-classical regime
In the limit of strong driving to the intermediate state Ω e ∼ γ e , the dynamics of the system is dominated by dissipation. As a method to study the system in the presence of both coherent driving and dissipation, we use quantum trajectory simulations  [15,70]. This dynamics can be intuitively understood from three different species of atoms, which we label, in analogy to the language of SFMs, as: • defects, that is atoms in |r ; • facilitated atoms, that is atoms in either |g or |e with one defect as nearestneighbor; • non-facilitated atoms, that is atoms in either |g or |e with no defect as nearestneighbor.
The dynamics of defects in the 1d 1-SFM is known to be mainly diffusive [70]. Let us illustrate the diffusion process of defects in more detail. A typical trajectory is shown in Fig. 2(b) and marked by a blue ellipsis. A defect can facilitate the excitation of a neighboring atom. In this example, the atom to the right is excited. For a short time, both are in the excited state but soon the original atom is de-excited: therefore, the defect has effectively moved to the right by one site. Another possible move is highlighted in Fig. 2(b) with a red ellipsis. In this second example, the left neighbor is excited but later on immediately de-excited. Here, the defect did effectively not move. The diffusion process can also be understood quantitatively in some limits. In the regime T 1 (or, equivalently, d eq 1), Γ 0→1 (0)/Γ 1→0 (0) ∼ d eq ∼ e −β and the diffusion constant is d eq /2 [15]. In the intermediate regime we studied in this work, the diffusion constant does not have such a simple analytical expression, since Γ 0→1 (0)/Γ 1→0 (0) = 0.33.
In the following, we study the relaxation of the three different species individually. Numerical simulations of N traj stochastically different trajectories are averaged in order to approximate the dynamics of the system described by the master equation (6). Results plotting the trajectory averaged Rydberg population R i (t) (see Appendix A) for each atom i are depicted in Fig. 3. Facilitated atoms (red lines) show a dynamical behaviour analogous to the one of resonantly driven single atoms, due to the antiblockade condition which acts as a facilitating term in the 1-SFM [cfr. Eq. (12)]. In fact, the inset of Fig. 3 shows that, in proximity of t = 0, the profiles of facilitated atoms are very well approximated by the resonant (∆ r = 0) single-atom Rydberg population ρ rr (t) (see Appendix A), the associated transition rate being Γ 0→1 (0).
Contrarily, the dynamics of defects (blue lines) and non-facilitated atoms (black lines) is slowed down with respect to the one of facilitated atoms. At t = 0, the weak transition |g ↔ |r of a defect is far-off-resonant and Γ 1→0 (∆ r ) is orders of magnitude suppressed with respect to Γ 1→0 (0) [cfr. Tab. 1]. A defect can in fact relax only after one of its neighboring facilitated atoms is excited to the Rydberg state. Similarly, at t = 0, the weak transition |g ↔ |r of a non-facilitated atom is far-off-resonant and Γ 0→1 (∆ r ) is orders of magnitude suppressed with respect to Γ 0→1 (0). The inset of Fig. 3 shows that profiles of non-facilitated atoms agree with the far-off resonant (∆ r = 0) single atom Rydberg population ρ rr (t): non-facilitated atoms, in order to be excited to the Rydberg state, have to wait for one of their two neighbors to be Rydberg-excited too. Note that for larger N we predict the existence of non-facilitated atoms with profiles of R i (t) even slower than the ones in Fig. 3. This can be regarded as a form of hierarchical relaxation [20], in that the more distant a non-facilitated atom is with respect to a defect at t = 0, the more slowed-down the evolution of R i (t) will be.
We now consider the concentration of Rydberg excitations The dynamics of r(t) for different N and V , starting with the same initial state is depicted in Fig. 4. After an initial transient relaxation, r(t) reaches a steady state constant value r ≡ lim t→∞ r(t), except in the case of the largest system size N = 10 and smaller interaction strength V = 3γ e , where it increases in time with a small but clear trend. We determine the origin of this anomalous slow growth of r(t) in the excess number of defects created from non-facilitated atoms within (originally suppressed) = 0 processes, which give a negative temperature contribution in the SFM language [compare Eq. (8) with Eq. (11)]. The reason why all reverse = 0 processes of destruction of isolated Rydberg excitations are not relevant, is in the asymmetry between off-resonant rates Γ 0→1 (∆ r ) > Γ 1→0 (∆ r ), which increases with |∆ r |, as depicted in the inset of Fig. 9. More quantitatively, the difference δr(t) between the excess in the concentration of Rydberg excitations at time t and the steady state value r can be  estimated as where N nf is the number of non-facilitated atoms at t = 0, which we approximate constant in time and equal to its value in the initial state. In our case, at the final time of the trajectory t max , δr(t max ) ∼ 0.06, which well agrees with numerical results plotted in Fig. 4, within statistical errors due to finite N traj simulated. Contrarily, the difference −δr (t) between the lack in the concentration of Rydberg excitations at time t and the steady state value r is where N d is the number of isolated defects at t = 0, which we also assume constant in time. Within our parameters, δr (t max ) ∼ 0.003, which is negligible with respect to both finite N traj errors and δr(t). From Eqs. (14), and (15) it is clear that, for ∆ r sufficiently large, l = 0 off-resonant processes are inhibited due to the suppression of Γ 0→1 (∆ r ) and Γ 1→0 (∆ r ) and consequently r(t) levels-off, reaching a steady state (see ∆ r 10γ e ). All far-off-resonant = 2 processes in Fig. 4 of merging and splitting of Rydberg excitations can be limited by a 'clever' choice of the initial concentration r(0) and configuration of Rydberg excitations. This is the case for, e.g., the initial states of one defect or two far defects in Fig. 4, in which the initial number N b of blocked atoms is zero. Blocked atoms are atoms with both left and right nearest-neighbour in the Rydberg state: the choice of such initial state ensures that, most probably, N b = 0 also at later times and we can thus limit the conditions that lead to the occurrence of merging processes. Plateaus in r(t) at long times in Fig. 4 show that a steady state is always characterized by the general result r > d eq , where and not r = d eq , as expected from the theory of 1-SFM. (It can be verified that d eq , in the limit Ω r Ω e , γ e , is well approximated by the single-atom steady state Rydberg population ρ rr (t → ∞) calculated in Appendix B. This confirms that the equilibrium ground-state of the 1-SFM is trivial, despite its spatially correlated and slowed down transient dynamics.) Our system thus does not thermalize in the long-time limit, but it converges to a different steady state. Nevertheless, the trend in Fig. 4 suggests that for larger N , r decreases approaching d eq .
Furthermore let us note that, due to finite size effects, r crucially depends on the initial configuration of Rydberg excitations, as depicted in Fig. 5(a): while the best approximation of the 1d 1-SFM result is given, as expected, by the lowest initial concentration of defects (blue line), for higher initial concentrations (black line), avoided merging processes become relevant and other corrections to r(t), similar to those in Eqs. (14) and (15), come into play. Contrarily, in the case in which both = 2 merging and splitting processes are present and in the limit of large N , our system would thermalize for every initial state in the long-time limit.
We finally mention that, ideally, one could think of simulating the full dynamics of the 1d 1-SFM introducing a third laser which weakly couples |g ↔ |r with a small Rabi frequency Ω r and detuning 2∆ r from |r . In this configuration, processes of merging of Rydberg excitations would be allowed also for the finite N we simulated numerically.

Quantum Regime
In the limit of vanishing driving of the |g ↔ |e transition, Ω e γ e , Ω r , the dynamics, up to times t γ e /Ω 2 e , is dominated by |g ↔ |r coherences. Starting from the fully coherent limit Ω e = 0, one can study, increasing Ω e , the crossover to the classical rate equation limit of Subsec. 4.1, i.e. Ω e ∼ γ e with Ω r Ω 2 e /γ e . In the case of Ω e = 0, the system Hamiltonian is that of N interacting two-level atoms The aim here is to study the non-equilibrium dynamics of the average concentration of Rydberg excitations r(t). The Hamiltonian in Eq. (17), assuming V = ∆ r and PBCs, can be mapped exactly to the so-called transverse field Ising model [71]. We note that the time evolution of the magnetization, which in our formalism can be easily related to r(t), following a quantum quench within the ordered (ferromagnetic) phase and across the phase transition to the disordered (paramagnetic) phase, has been extensively studied in Refs. [46] and [47]. There, the authors provide explicit analytical results in the thermodynamic limit N → ∞, corroborated by numerical computations showing that finite-size effects are negligible.
In the following, we show that for V = ∆ r and in the limit Ω r ∆ r , one can approximate the dynamics of excitations of the transverse field Ising model [or of Eq. (17)], whose Hilbert space scales exponentially (2 N ), with the dynamics of excitations in a simplified model with only polynomial scaling (N 2 ). The Hilbert space of such a simplified model (Fig. 6) is just a subspace of the transverse field Ising model total Hilbert space, and is univocally determined by the initial state of the dynamics.
If the initial state has energy −V , the dynamics of the system is limited, for times t ∆ r /Ω 2 r , to the subspace |rrr , |ggr , |grg , |grr , |rgg , |rrg of the full Hilbert space. The matrixH N =3 tl , enclosed within the dotted lines in the down-right corner of H N =3 tl in Eq. (18), governs the dynamics of the system in the reduced subspace. In our example, the subspace is characterized by the presence of one cluster of Rydberg excitations with varying size. Here, a cluster is a group of contiguous Rydberg excitations without atoms in |g in between them. The total number of Rydberg excitations N r in the cluster changes during the dynamics (see Fig. 6), but the number of clusters does not. At this point we remark that, at times t ∆ r /Ω 2 r , second (and higher) order perturbative processes in the small parameter Ω r /∆ r will drive the dynamics outside the subspace. This discussion can be extended also for a general N and for every initial condition. The number of clusters will be a constant of motion and it will always be determined by the initial state. For example, in Fig. 7, we compare the evolution of r(t) governed by H N tl (green line) andH N tl (red line) for N = 9. We restrict the dynamics of the system in the subspace with energy −V starting with one Rydberg excitation in the center of the open chain and all the remaining atoms in |g . Numerical simulations of the full (H N =9 tl ) and effective (H N =9 tl ) dynamics are in excellent agreement for t ∆ r /Ω 2 r . In this time window, the peaks of oscillations of r(t) have a periodicity 2π/(Ω r /2).
Remarkably, the dimension ofH N tl scales as N (N +1)/2 contrarily to the exponential growth 2 N of H N tl . Such quadratic scaling for large N is understood from the the fact that the effective dynamics of our system, in that limit, can be interpreted as hopping processes of a single particle in a 2d square lattice, which is the well-known 2d tightbinding model. We will comment on this in a few lines.
It is worthwhile first to study properties ofH N tl more in detail. Setting ∆ r = V = 0 without loss of generality (in the reduced subspace ofH N tl , V is just a constant energy offset) and considering again the example N = 3, the eigenvalues E N =3 λ ofH N =3 tl , written in ascending order, are In Fig. 8, we plot the spectrum of eigenvalues E N λ for different N . Two features emerge: for Ω e = 0. The red curve reproduces with high accuracy the exact dynamics evolving the reduced matrixH N =9 tl for Ω e = 0 (see discussion in the main text). For all Ω e = 0, N traj = 500 trajectories have been simulated.
(i) Independently of N , there is always at least one dark state, that is an eigenstate ofH N tl with E N λ = 0. We verified that the total number of degenerate dark states increases linearly in N . More specifically, for N 3 [N 2] odd [even], it scales with (N + 1)/2 [N/2]. In the case of N = 3, the two (non-normalized) dark states can be expressed as In particular, |E N exists for every N , and it is obtained from the linear combination of all states consisting of only one Rydberg excitation with alternating plus and minus signs, in analogy to |E N =3 . Due to the linear scaling in N of the degeneracy of dark states, and since bulk properties of the system scale quadratically, we can think of dark states also as edge states in the Hilbert space ofH N tl : for example, |E N =3 is localized at the lower boundary of the scheme depicted in Fig. 6.
(ii) The eigenvalue energy spectrum, shown in Fig. 8 for different N , is bounded by −2Ω r and 2Ω r for large N . This is a consequence of the fact that, as already mentioned, our model, in the limit N → ∞, can be mapped exactly to the 2d tight binding model. Its dispersion relation looks like where ε( q) is the energy of the particle at quasi-momentum q = q 1 b 1 +q 2 b 2 and J is the hopping amplitude between adjacent lattice sites with spacing s. b 1,2 = (2π/s 2 ) a 1,2 are the vectors which span the reciprocal lattice generated by a 1 = s(1, 0) and a 2 = s(0, 1) [(1, 0) and (0, 1) is the standard basis in the Euclidean plane]. Furthermore, due to PBCs, q 1,2 s = 2πm 1,2 /N and m 1,2 ∈ [−N/2, N/2]: from Eq. (21), ε( q) has thus a maximum width of 8J. In our case J ≡ Ω r /2, and the spectrum width of ε( q), 4Ω r , coincides with the one observed in Fig. 8 for large N . We stress that this mapping, which simplifies a d−dimensional many-body problem with nearest-neighbour interactions to a (d + 1)−dimensional single particle problem, is exact only in the limit N → ∞, where boundary effects in our model are absent, i.e. PBCs, assumed in Eq. (21), hold. Let us finally comment on the transition from the fully coherent limit to the classical limit switching on and increasing Ω e , while keeping fixed all other parameters. Even for small Ω e = (0.03, 0.1)γ e (blue and orange lines, respectively, in Fig. 7), the system dynamics deviates from the coherent case, since Ω e Ω r , and oscillations of r(t) are damped by dissipative processes due to the coupling to the short-lived intermediate level |e . For t γ e /Ω 2 e , r(t) is driven to its steady state value r, which is well approximated by ρ rr (t → ∞) [see Appendix B]. Upon further increase of Ω e , the rate equation limit is recovered, in which r(t) increases slowly with time, oscillations are completely suppressed and, as already discussed in Subsec. 4.1, r does not coincide with ρ rr (t → ∞).

External disorder
The external disorder, arising from static spatially-random detunings from |r , ∆ ri = 0 in Eq. (1), considerably changes the dynamics of the relaxation. In the quantum regime, they are responsible of the many-body localization transition [11,12,13,14]. Here, we focus on the effects of random detunings on the relaxation of our system in the rate equation limit of the model. We randomly choose the values of ∆ ri within the interval [−A, A] according to a uniform distribution. Moreover, we assume 0 < A ∆ r , in order to avoid accidental balancing of the homogeneous and random detunings, ∆ ri ∼ ∆ r . Quantum trajectories are simulated according to the following procedure: (i) construct N rnd different sets of random detunings (∆ r1 , ..., ∆ ri , ..., ∆ rN ) to be added to the homogeneous detuning ∆ r ; (ii) for each random field configuration (∆ r1 +∆ r , ..., ∆ ri +∆ r , ..., ∆ rN +∆ r ), simulate N traj stochastically independent quantum trajectories; (iii) average the desired observable over the total number of trajectories N traj N rnd . Note that, for sufficiently large N rnd , no extra energy is added to the system, because [−A, A] is centered around 0. Fig. 9 shows that increasing N rnd (for a sufficiently large N traj ), all atoms relax faster toward the steady state with respect to the case of zero random detuning. This tendency can be understood looking at Tab. 2, which summarizes the effect of a finite random detuning on each species.
As depicted in the inset of Fig 9, for A ∆ r , rates with a finite random detuning Γ n i →1−n i (∆ ri ) of facilitated atoms are larger than the corresponding rates with no random detuning Γ n i →1−n i (0) (see the two local maxima of both rates within the beige area). As a consequence, at each time, the relaxation of facilitated atoms tends to be sped up by any finite random detuning, independently from its sign. Contrarily, the relaxation of defects/non-facilitated atoms and blocked atoms is either slowed down or sped up, depending on the sign and the magnitude of ∆ ri (see shaded green and blue areas, respectively). As a general rule, blocked atoms and defects/non-facilitated atoms are affected in the opposite way by the same random detuning: e.g., for ∆ ri > 0 defects/non-facilitated atoms are slowed down, whereas blocked atoms are sped up.
The result of an averaged (over trajectories and random field realizations) faster relaxation to the steady state can be understood more deeply as a consequence of the breaking of integrability of the transverse field Ising model the Hamiltonian in Eq. (17) 0 − − − Figure 9. Evolution of the trajectory averaged Rydberg population R i (t) for different finite random detunings from the Rydberg state, ∆ ri = 0. The initial configuration and parameters are the same as in Fig. 3 and A = 0.5γ e (see main text). Different colors refer to different numbers of random detuning configurations N rnd , for a fixed number of trajectories simulated in each random configuration, N traj = 100. Inset: Semi-log plot of single particle rates Γ ni→1−ni (∆ r ) [with n i = 0, 1] from Eqs. (5) and (4) as a function of ∆ r . The shaded areas delimit the intervals in which the total detuning ∆ r ± |∆ ri | ranges for facilitated atoms (beige), defects/non-facilitated atoms (green) and blocked atoms (blue).
can be mapped to. We would like to point out that, in the different limit of large random detunings, i.e. A ∼ ∆ r , the dynamics is contrarily slowed down, in that facilitated atoms have much higher probabilities of turning non-facilitated.

Experimental considerations
In the following, we discuss and validate the various assumptions made in Secs. 2 and 4. These are: (i) neglecting the decay from the Rydberg state during the relevant experimental timescales, (ii) truncating long-range van der Waals or dipolar interactions to nearest-neighbor interactions, (iii) assuming atoms in ground-and Rydberg states to be trapped simultaneously in the same optical lattice, and (iv) neglecting heating due to photons emitted from the |g ↔ |e transition.
(i) The lifetime of Rydberg states is determined by both radiative decay (spontaneous emission) to lower-lying electronic energy levels and blackbody radiation which induces transitions to higher-and lower-lying neighbouring Rydberg states [23,24]. While for low-lying electronic levels spontaneous decay to the ground-state dominates, for Rydberg states the interaction with thermally occupied infrared photons, which is the main contribution in the overall lifetime, leads to strong mechanical effects [72]. In general, (alkali-metal) Rydberg states have, compared to low-lying excited states, exceedingly long lifetimes ∝ n 3 (∝ n 5 ) for low (high) angular momentum states. For n 80, these can be of the order of hundreds of µs [73]. In our setup, the rate of decay from the Rydberg state γ r can be neglected when holds [62], implying that stimulated emission processes from |r have to be much larger than spontaneous ones, which is typically the case for Rydberg states with large enough principal quantum numbers n and ∆ r = 0. For ∆ r = 0, nonetheless, care has to be taken since, contrarily to the case of the general requirement for well-defined quantum jumps to occur [62] the condition in Eq. (22) becomes stricter with respect to the resonant case ∆ r = 0.
(ii) In Fig. 5(b) we check, for a certain initial state and no random detunings, the reliability of the approximation of truncating Rydberg interactions to nearestneighbour. While for van der Waals potentials there is no significative quantitative disagreement, dipolar interactions result to be detrimental, despite the qualitative slow growth of the number of Rydberg excitations is recovered in the long time limit. Nevertheless, Rydberg atoms at distances of a few µm, typical of large spacing optical lattices [48,49,50,51,52,53], are known to interact with van der Waals potentials [23], and our approximation is thus not only qualitatively but also quantitatively justified.
(iii) It is important, in our setup, that atoms are trapped in the same optical lattice during both bright and dark periods. While ground-state atoms can be trapped using different state-of-the-art techniques, it has been suggested to optically trap Rydberg atoms via ponderomotive potentials from spatially modulated light fields [74] or in dcelectric or magnetic fields [75]. In general, electronic ground-states and Rydberg states will experience different optical trapping potentials due to their different polarizabilities. A promising route to trap ground-state atoms and Rydberg-atoms in the same lattice is also to use alkaline-earth atoms as suggested in Ref. [76].
(iv) Heating due to photons emitted from the strong transition can be considerably reduced adding a finite detuning from |e , ∆ e = 0, in order to Doppler-cool atoms with the same fluorescence laser. Results of quantum trajectory simulations (not shown) indicate that the typical slowed relaxation of atoms is qualitatively unaltered, although larger trajectory times with respect to the case ∆ e = 0 are needed in order to observe the convergence of r(t) to its steady state value r.

Conclusions
We have shown that the non-equilibrium dynamics of Rydberg excitations in three-level atoms can be tuned from the classical to the quantum regime with the laser driving to a short-lived intermediate excited state. The anti-blockade condition, where the nearestneighbor interaction among Rydberg levels equals the single atom laser detuning from the Rydberg state, implements a form of kinetic constraint. This phenomenon resembles the assisted processes of creation and destruction of defects in the One-Spin Facilitated Model. The absence of merging processes in our system prevents the thermalization to the ground-state of the One-Spin Facilitated Model, despite the concentration of Rydberg excitations can be, in the long time limit, stabilized to some steady state value for sufficiently large interactions and detunings from the Rydberg state. In the intermediate regime of interactions, we observed instead a (tuneable) slow increase of the concentration of Rydberg excitations, caused by the appearance of processes which are not present in the One-Spin Facilitated Model. In this regime, an additional slow relaxation timescale emerges.
In the coherent limit, the dynamics of the resulting interacting two-level system is limited to a subspace of the Hilbert space. This subspace is determined by the choice of the initial state, whose number of clusters of Rydberg excitations is conserved. However, the size of these clusters changes during the evolution. We derived an effective reduced Hamiltonian and verified excellent agreement with simulations of the exact dynamics of the full Hamiltonian.
Finally, we have shown that, in the rate equation limit, the competition of the internal disorder of glasses with an externally imposed disorder can result in the slowing down but also in the speeding up of the overall relaxation towards the steady state. The setup we presented allows for the investigation of the quantum dynamics of Spin Facilitated Models in the presence of random fields. As an outlook, we plan to study the dynamics of the coherent cluster model in the spirit of many-body localization.

Appendix A. Review of the properties of V-shaped three-level atom
In this Appendix we review properties of a single atom in the V-configuration and describe both quantitatively and qualitatively the physical mechanism of quantum jumps in such a setup.
Let us imagine to measure photons emitted by the three-level atom described by the Hamiltonian in Eq. (2) of the main text with a photodetector, which we assume to be 100% efficient. The theory of continuos measurement [62] ensures that, under continuous monitoring of the strong transition |g ↔ |e with our photodetector and during time intervals where no photons are detected, the system wave function |ψ(t) evolves according to the Schrödinger equation i d dt |ψ(t) = H eff |ψ(t) , with an effective non-Hermitian Hamiltonian and H given by Eq. (2) in the main text. The general solution for the system wave function is where λ n and |u n are the eigenvalues and associated eigenvectors of H eff , respectively [66]. The coefficients c n can be determined from the initial condition, e.g. |ψ(0) = |g . In the limit Ω r Ω e , γ e , ∆ r , λ n and |u n can be computed in perturbation theory.
Appendix A.1. Brillouin-Wigner perturbation theory calculation of eigenvalues and eigenvectors of H eff In the following, we summarize the steps for the calculation of λ n and |u n of H eff in Eq. (A.1) within a Brillouin-Wigner perturbation theory approach. In particular, we focus on the calculation of λ e and |u 3 (the other eigenvalues and eigenvectors can be straightforardly computed following the same procedure). In the {|g , |e , |r } basis, H eff has the following matrix form where U † is the Hermitian conjugate of U , E m = i 4 −γ e + γ 2 e − 4Ω 2 e and E p = − i 4 γ e + γ 2 e − 4Ω 2 e . Applying the same transformation on H 1 ,H 1 = U H 1 U † , we get where Ω c = Ω r cosh { 1 2 arctanh 2Ωe γe } and Ω s = Ω r sinh { 1 2 arctanh 2Ωe γe }. The unitary transformation U defines a change of basis from {|g , |e , |r } to, We now define two operators P and Q which project onto |r and its complementary subspace, respectively. Since From classical to quantum non-equilibrium dynamics of Rydberg excitations in optical lattices23 Here, [A] i,j indicates the i-th row and j-th column of the matrix A and R 0 = (E 0 1 − QH 0 Q) −1 is the so-called resolvent. Up to second order in Ω r , λ 3 can be approximated with We can now calculate the first-order correction to |r originated by the other two eigenstates ofH 0 , |− and |+ , both mixed to |r by the transformed perturbation HamiltonianH 1 . For example, the eigenvector |u 3 associated to λ 3 , is: As mentioned above, the coefficients c n can be calculated from the initial condition. For example, setting |ψ(0) = |g , we get (A.13) The eigenvalues λ n are complex numbers all with negative imaginary parts Im{λ n } [see, e.g., Eq. (A.11)]. As a consequence, |ψ(t) does not conserve the norm during the evolution. For Ω r Ω 2 e /γ e , one of the eigenvalues of H eff , λ 3 in Eq. (A.11), has a much less negative imaginary part with respect to the other two, |Im{λ 3 }| |Im{λ 1 }|, |Im{λ 2 }| [66]. This determines long tails in the delay function D 0 (t) = 3 n=1 |c n | 2 e 2tIm{λn} , (A.14) defined as the conditional probability that no photon has been detected until time t, given a photon count has occurred at t = 0. Note that D 0 (t) is exactly the norm of |ψ(t) , and its monotonous decay is dominated, in the long time limit, by Im{λ 3 }. Thus Im{λ 3 } is responsible of the existence of a small but finite probability of a long time window with no photon counts. When such a long interval of no photon counts occurs, the atom is prepared in the Rydberg state and |ψ(t) is dominated by |u 3 , which in turn can be expressed as |r plus a perturbatively small admixture of |g and |e (see Appendix A.1).
The wavefunction of trajectories, or system realizations, |ψ α (t) (where α labels each of the N traj stochastically independent trajectories), can be simulated following the procedure below [58]: (i) Prepare the system in a certain (normalized) state |ψ(t 0 ) α at the initial time t 0 .
(ii) Compute the probability of the occurrence of a quantum jump at t 0 from p α (t 0 ) = γ e ψ (t 0 )|c † c|ψ(t 0 ) α dt 1, where dt is an ideally infinitesimal time step. The corresponding probability of no quantum jump is 1 − p α (t 0 ) ≈ 1.
(iii) According to these probabilities, randomly choose a no jump/jump event.
(v) Repeat the procedure from (i) to (iv), until the desired final trajectory time t max is reached.
A typical experimental sequence of photon detections looks like an alternation of bright and dark periods. During a bright period, labelled by 0, the atomic population is rapidly cycled between |g and |e and the time between two successive photon detections is of the order of the lifetime of |e , γ −1 e . In contrast, during a dark period, labelled with 1, the atomic population is shelved in |r and no photons are detected. An approximation of the exact time evolution of ρ(t) governed by Eq. In the following, the observable we are interested in is the Rydberg population, that is, the expectation value of the Rydberg projector |r r| in |ψ(t) α , Averaging r α (t) over N traj leads to In analogy to Eq. (A.15), lim N traj →∞ R(t) = ρ rr (t), where ρ rr (t) ≡ r|ρ(t)|r is the Rydberg population whose evolution is governed by the exact master equation in Eq. (3). Fig. A1 shows that R(t) converges to ρ rr (t) by increasing N traj . In the following we illustrate three independent methods for the estimation of the average duration of bright and dark periods T 0 and T 1 , respectively, or equivalently, of the transition probabilities (rates) from 0 to 1, Γ 0→1 ≡ T −1 0 , and from 1 to 0, Γ 1→0 ≡ T −1 1 . In Ref. [60], these rates are found to be together with ρ rr (0) = 1. Rates can thus be evaluated by linearization of the short-time evolution of ρ rr (t) with the proper initial condition (e.g., from Fig. A1, we can estimate Γ 1→0 ). For comparison, we show the analytical derivation of the rates in Eqs. (4) and (5) of the main text. In the long time limit D 0 (t) ∼ p e −t/T 1 , where p is the (small) probability to be in a dark period [61,62]. Comparing it with Eq. (A.14) one gets where both results are obtained substituting λ 3 and c 3 calculated in Appendix A.1. In order to compute Γ 0→1 (∆ r ), we note that, during a bright period, |r is almost never populated and the effective emission rate will be that of the (resonantly driven) twolevel system |g ↔ |e , i.e. γ e Ω 2 e /(γ 2 e + 2Ω 2 e ). Then, Γ 0→1 (∆ r ) is simply the two-level emission rate multiplied by p [66], . (A.22) From classical to quantum non-equilibrium dynamics of Rydberg excitations in optical lattices26 The third independent method relies on the direct numerical calculation of the average over trajectories of the time-averaged length of dark and bright periods in each trajectory. Results for Γ 1→0 , obtained with the three different methods, are compared in Tab. A1: