Effective dynamics of strongly dissipative Rydberg gases

We investigate the evolution of interacting Rydberg gases in the limit of strong noise and dissipation. Starting from a description in terms of a Markovian quantum master equation we derive effective equations of motion that govern the dynamics on a"coarse-grained"timescale where fast dissipative degrees of freedom have been adiabatically eliminated. Specifically, we consider two scenarios which are of relevance for current theoretical and experimental studies --- Rydberg atoms in a two-level (spin) approximation subject to strong dephasing noise as well as Rydberg atoms under so-called electromagnetically induced transparency (EIT) conditions and fast radiative decay. In the former case we find that the effective dynamics is described by classical rate equations up to second order in an appropriate perturbative expansion. This drastically reduces the computational complexity of numerical simulations in comparison to the full quantum master equation. When accounting for the fourth order correction in this expansion, however, we find that the resulting equation breaks the preservation of positivity and thus cannot be interpreted as a proper classical master rate equation. In the EIT system we find that the expansion up to second order retains information not only on the"classical"observables, but also on some quantum coherences. Nevertheless, this perturbative treatment still achieves a non-trivial reduction of complexity with respect to the original problem.


Introduction
Gases of interacting highly excited Rydberg atoms are becoming an increasingly popular theoretical and experimental platform for the investigation of the physics of strongly interacting many-body systems [1,2]. The main distinction between these systems and "traditional" ones restricted to low-lying excited states lies in the huge enhancement in the interaction between the atoms, which can be several orders of magnitude stronger for typical experimental parameters. Indeed, two atoms in a Rydberg state usually experience extremely strong dipole-dipole or van der Waals forces (see e.g., [3]). These in turn considerably affect both the static [4,5,6,7,8] and dynamic [9,10,11,12,13,14] properties of the system.
The rather long lifetime of these Rydberg states allows for probing of the coherent evolution of these many-body systems up to relatively long time-scales. Here, however, we will focus on situations (currently studied with great interest) where accounting for dissipative processes leads to interesting changes in both the dynamics and the stationary properties. These processes emerge due to the coupling of the system to external degrees of freedom, which produce e.g. decay via spontaneous emission (fluorescence) and noise-induced loss of coherence (dephasing). The evolution is typically well described in terms of a quantum master equation with Markovian noise (an overview on methods for treating the non-Markovian case as well can be found in [15]). While such a modelling is certainly among the simplest descriptions of an open quantum system it still poses severe challenges when trying to conduct a numerical treatment for large system sizes N. This is due to the very fast increase in the dimension of the many-body Liouville space (b 2N with b being the dimension of the single-particle Hilbert space, e.g., b = 2 for an Ising spin). Some procedures to address this issue have been developed which divide the system into two parts, a subsystem of interest and an environment which is traced away [16,17].
A somewhat different framework arises when it is possible to identify degrees of freedom that evolve on vastly different timescales. Adiabatically eliminating fastevolving ones might then allow the derivation of an effective equation of motion for the remaining slow degrees of freedom that portrays a reduced complexity and might be amenable to numerical treatment (see Appendix A or [16] for a general description of the method, based on the Nakajima-Zwanzig projection formalism [18,19]). One of the first works to apply an idea along those lines to Rydberg systems was that of Ates et al. [20,21], where it was shown that -in an appropriate limit -properties of the stationary state of an interacting Rydberg gas can be extracted via a classical Monte-Carlo method. This made the simulation of stationary properties of large scale Rydberg systems feasible and several recent works employ variations of the same method [14,22,23,24,25,26,27].
Besides studying the properties of the stationary state there is great current interest in the understanding of the out-of-equilibrium dynamics of these interacting Rydberg gases. Here, effective equations of motion that describe the systems' dynamics in terms of a classical rate equation have been put forward and used by several authors both in purely theoretical works [14,28,29] and to model actual experimental data [30,31].
The purpose of this paper is to give a detailed account of the derivation of effective equations of motion that describe the many-body dynamics of interacting Rydberg gases in the limit of strong dissipation. Specifically, we will discuss the two scenarios depicted in figure 1 that are directly relevant to current experiments: The first one is that of a Rydberg gas in which atoms are modelled by coherently driven and interacting two-level systems. Here dissipation is present in the form of dephasing noise that quickly destroys coherent superpositions between the two states. This corresponds approximately to the experimental situations discussed in Refs. [30,31]. We show that in the limit of strong dephasing the dynamics of the interacting two-level systems is described -up to second order in the relevant perturbative expansion -by a classical rate equation; the corresponding stochastic process is described by single spin flips subject to kinetic constraints [32]. We then proceed further and calculate the fourth-order corrections, showing that they result in new processes, such as simultaneous two-spin flips. However, it turns out that, unlike for the second order case, there are domains in the space of physical parameters for which some of the "rates" become negative, thereby breaking the conservation of positivity. Hence, a standard treatment in terms of a classical stochastic dynamics is not always possible. Yet, the perturbative expansion is formally correct and our numerical analysis shows that said breakdown only affects the initial stage of the dynamics, whereas for long times a good agreement with numerically exact data is still found. We conclude the discussion by including (radiative) decay from the upper to the lower atomic level in the rate equation treatment.
The second scenario we are considering is that of Rydberg gases under so-called electromagnetically induced transparency (EIT) conditions. In this regime, which has been studied experimentally in Refs. [33,34,35,36,37,38,39], atoms are modelled by coherently driven and interacting three-level systems. Dissipation enters through a fast (radiative) decay of the middle level to the lower one. In this case, the fast dissipative dynamics is of a different nature and does not necessarily project onto a classical spin configuration. Therefore, information on some quantum coherences must be retained. Despite the fact that one does not gain as simple a description as in the case above, the resulting reduced equation is still in Lindblad form and offers therefore a simplified alternative to the study of the one acting on the whole Hilbert space. In particular, the interatomic interaction needs to be taken into account coherently, while the elimination of the fast decay leads to an unusual form of effective dissipation that drives the system into coherent superposition states. We conclude the discussion of the three-level system by considering the limit of infinitely strong nearest-neighbour interaction, for which we find an effective purely-dissipative quantum dynamics that is reminiscent of that of quantum analogues of kinetically constrained spin models [40].
The paper is organized such that the central results and conclusions are presented in the main text, whilst a more detailed formal derivation is provided in the final appendices V (b) V (a) Figure 1. Atomic level schemes considered in this paper. Two atoms effectively only interact when they are both excited to the Rydberg state |↑ ; here we denote with V the strength of such an interaction. Most of the current experiments can be modelled by the following two descriptions: (a): Two-level atoms driven by a laser with Rabi frequency Ω and detuning ∆. The main dissipation mechanism we consider here is dephasing (at rate γ) of superpositions between the states |↑ and |↓ . We will also take into account the decay from |↑ to |↓ with rate Γ ryd , which is usually small compared to γ and thus will be treated perturbatively. (b): Three-level atoms in an EIT configuration, i.e. where the excitation of the Rydberg states is performed via a transition between the ground state |↓ and the intermediate state |← . Within the time-scales of interest, spontaneous decay processes are assumed to be relevant only for the intermediate state, whose inverse lifetime is Γ.
to which we refer in the appropriate sections.

Two-level Rydberg atoms in the presence of strong dephasing noise
We consider here a gas of N atoms with two relevant internal levels as shown in figure 1(a). We assume that the timescale of the external motion of the atoms is much larger than the one in which the electronic dynamics takes place. This "frozen gas" picture has been shown to be adequate in a vast number of theoretical and experimental works. The ground state |↓ and the Rydberg state |↑ are coupled by a laser with Rabi frequency Ω and detuning ∆. In order to avoid having an explicit angular dependence of the interaction, the excited state is typically chosen to have spherical symmetry, i.e., to be an S orbital. Sharing the same parity of the ground state, however, it is not possible to reach it via a single dipole transition; in practice, this excitation must be achieved by means of a two-photon process via a far off-resonant excitation of an intermediate state. In this case the two-level approximation is adequate. Later we will, however, also account for the case of near-resonant excitation of such an intermediate level, which makes it necessary to include it as well in the description.
When two atoms (k and m) are simultaneously excited to the Rydberg state they interact due to the electrostatic coupling of the respective (permanent or induced) dipole moments. The strength of this interaction V km is therefore of the form [3] V km = C p |r k − r m | p , where r k denotes the position of the k-th atom, C p is the dispersion coefficient and p characterizes the interaction type: p = 3 stands for dipole-dipole and p = 6 for van der Waals forces. The dynamics of the density matrix ρ of the system is described by a quantum master equatioṅ The first term describes the coherent evolution of the system which (within the rotating wave approximation) is governed by the many-body Hamiltonian The generator of the dissipative dynamics is modeled in terms of a dissipator in Lindblad form which in case of the dephasing noise considered here is given by with γ being the dephasing rate. Note, that this relies on the assumption that the noise can be considered white and spatially uncorrelated, i.e., acting independently on each atom. In practice, this is not always the case: for instance, dephasing noise can be introduced by fluctuations in the laser fields with a finite correlation length; if the typical interatomic distance is smaller than this correlation length, the noise experienced by nearby atoms will be spatially-correlated. Nevertheless, to consider independent fluctuations represents a reasonable approach for dilute Rydberg ensembles and, moreover, recent experimental work [30,31,41] suggests that this approximation captures the essential physics of the problem.

Second order effective evolution
We are now interested in deriving an effective equation of motion for the system in the regime where the dephasing rate is large or more precisely γ ≫ Ω. In this limit coherent superpositions of the local atomic states |↓ and |↑ will dephase exponentially fast on timescales of the order of γ −1 . On timescales longer than this dephasing time the density matrix will therefore no longer show coherences and we can thus describe the system's state by a reduced density matrix µ which includes only the diagonal elements of ρ, i.e., only the probabilities of the classical spin configurations (direct products of the form |· · · ↑↑↓↑ · · · ). Note that, because of this, the only observables which can be calculated within this scheme are the diagonal ones, i.e., those which can be written as combinations of n k -s and the identity. A more formal version of the arguments above is given in Appendix B.
The coherent flipping induced by the laser provides a much slower dynamics which, due to the fast action of the dephasing, can be effectively projected onto the stationary subspace of the dephasing and accounted for in a perturbative expansion in powers of Ω. Accordingly, the effective equation of motion can be cast in the forṁ where L (α) is the evolution operator of order Ω α . All odd terms of this series identically vanish, hence the first non-vanishing term is of second order, i.e. L (2) . The corresponding truncated evolution of µ at this level readṡ as was derived in [28]. Since µ is the diagonal matrix of probabilities associated to classical spin configurations, (3) describes a continuous-time stochastic process which flips the k-th spin with operator-valued rate Γ k . More precisely, the rate depends via the interaction term on the configurations of all sites but the k-th one. The expression (3) is therefore equivalent to a kinetically-constrained rate equation [32,42,43], i.e., it can be regarded in terms of a trivial evolution (flipping one spin at a time) subject to a non-trivial constraint (here the number of excitations present in the neighbourhood). This can be more clearly understood by introducing the probability vector v = diag(µ) which evolves according tȯ Here each term describes the incoherent state change of an atom from |↓ to |↑ , and viceversa, with rate Γ k . This representation has the distinct advantage of allowing for numerical investigations of large scale systems by virtue of classical Monte-Carlo simulations. This fact has been exploited in a number of recent works [14,28,29].
In order to assess the validity of this effective description, we have computed the evolution of small systems (up to N = 9 atoms) accounting for both the full quantum master equation (numerically simulated via Quantum Jump Monte-Carlo [44,45]) and the resulting classical second order equation [numerically exact integration of (3)]. For our simulations we have considered the atoms trapped in a one dimensional chain with lattice constant a, with one atom per site, periodic boundary conditions and van der Waals interaction. We have chosen a value for the dephasing compatible with the perturbative requirement, namely γ = 10Ω, and fixed the value of the nearest neighbour interaction to V = C 6 /a 6 = 10Ω. In figure 2 we show for a system of N = 9 atoms the  Figure 2. Time evolution of the density of excitations n , its associated fluctuations n 2 − n 2 and the density-density correlations g 2 (d) for nearest (d = 1) and nextnearest neighbours (d = 2). In all cases the initial state is the one without excitations ⊗ k |↓ k . We compare the results obtained from the Quantum Jump Monte-Carlo simulation of the full quantum system and the numerically exact integration of the effective Master equation obtained up to second [given by (3)] and fourth order [adding the contribution given by (5)] for N = 9 atoms. The parameters used in the simulations shown are ∆ = −10Ω, 0, 10Ω, γ = 10Ω and V = 10Ω.
short-time evolution of the mean density of excitations n = k n k /N, its fluctuations n 2 − n 2 and the density-density correlations between nearest (d = 1) and next-nearest (d = 2) neighbours, with d being the distance in units of a. We choose to focus on the short time behaviour as here the difference between the exact and effective dynamics becomes most visible. For longer times, both dynamics reach the same stationary state, which is completely mixed, i.e., µ and ρ become proportional to the identity. This is a consequence of the fact that the dissipator (1) is constructed solely upon Hermitian jump operators n k . In general there is good agreement between the approximate and exact dynamics for the chosen parameter sets and observables. The initial increase in the density of excitations for very short times (smaller than Ω −1 ) is not expected to be well-captured by the approximation introduced above. Considering that the initial state is the completely polarised one k |↓ k , which belongs to the stationary space of the dephasing, in the full quantum problem the early stage of the dynamics is approximately driven by the coherent part and thus must be reversible. Therefore the density can only start from 0 with vanishing slope, so that the initial increase of the excitation density is always proportional to (at least) t 2 . This behavior cannot be captured by the classical rate equation (3). This can be understood by determining the equation of motion of the density of excitations n k at site k, which reads Due to the factorised nature of the initial state, Γ k n k = Γ k n k . Hence, starting from n k = 0, we observe that the gradient at t = 0 is indeed different from 0 and the density's initial increase is proportional to t (see e.g. [31]). This difference is most obvious in the resonant case ∆ = 0, where a magnified view of the very short-time regime is shown (see second row in figure 2). For times beyond Ω −1 one can still observe small amplitude oscillations of the numerically exact solution; these are due to the dampening effects of the coherences and are thus not captured by the approximate dynamics. Compatibly with our assumptions, they become less and less pronounced as γ increases and off-diagonal terms are quickly damped out. Within the parameter regime analysed here, it also becomes apparent that the agreement is enhanced for positive values of the detuning. At an intuitive level, this can be related to the form of the rates (4), which can be interpreted as the effective perturbative parameters. It is clear that, for negative ∆, the effects of the detuning and the interaction are competing and one can generally obtain rates which are smaller than in the case of positive detuning.

Fourth order corrections
The next non-vanishing order in the perturbative expansion (2) of the effective master equation is the fourth one (∝ Ω 4 ). Its structure is considerably more involved than the second order contribution. It can be divided into five terms: where each superoperator "G" can be represented in general as with R i and R ′ i being hermitian operator-valued coefficients. Analogously to the second order rates above, the structure of these coefficients is diagonal, in the sense that they can be written as non-linear combinations of the local density operators n q . Their specific form is, however, much more complicated in this case: Note that the preservation of the trace is here ensured by the fact that every R ′ km i except R ′ km 1 commutes with σ x k , whilst R km 1,4 commute with σ x m and R km 2,3 with σ x k . The last term of (5) is instead of the form where β k is again an operator-valued rate (i.e., a kinetic constraint) commuting with σ x k , which reads A detailed derivation of the fourth order contribution to the perturbative expansion (5) and the specific forms of the rates are given in Appendix B. We now wish to give a stochastic interpretation to the terms resulting from this fourth order expansion. The action of the "F " superoperator (7) displays the same structure as the second order terms (3). Thus, it represents simply a perturbative correction of order Ω 4 to these processes. The action of the "G" superoperators (6) is more involved. Here, we start by separating terms that lead to a single spin-flip from the ones that lead to two correlated spin flips. The former constitute an additional fourth-order correction to (3), whereas the latter introduce novel dynamical processes. Collecting all terms up to fourth order, the effective equation for the probability vector v reads nowv Here, the single-flip Γ s k and double-flip rates Γ d km are given by Conservation of probability is ensured by the preservation of the trace discussed above. The remaining requirement for obtaining a proper classical rate equation is that it must preserve positivity as well, i.e., probabilities cannot become negative. This is equivalent to requiring that every stochastic rate must be positive. Within the perturbative regime, this is automatically satisfied for all single spin-flip processes, since the second-order rates Γ k constitute the leading terms and are strictly positive for all k. On the other hand, the two spin-flip ones can become negative for some choices of the parameters, signalling the breakdown of this simplified stochastic interpretation.
In figure 3 we show these two regimes emerging from a numerical analysis in the V -∆ plane for γ = 10Ω and two system sizes: N = 4 and N = 9. The white area corresponds to parameter choices for which all rates are positive, whereas within the black one at least one is negative. We note that typically both V and ∆ must be quite large compared to Ω in order for the stochastic interpretation to formally hold. We also observe that the boundaries of the "negative (black) region" shift towards larger values of both V and ∆ as the system size is increased. In passing, we remark that, as it should be expected, in the non-interacting case V = 0 the rates for two spin-flip processes vanish, leaving finite only the corrections to the second-order term.
It is worth mentioning that, even when positivity is not ensured for all times, this approach yields a non-negligible reduction of the degrees of freedom. We have numerically integrated the fourth-order equation (8) and compared it with the full quantum evolution. The corresponding results for the evolution of the average density, its fluctuations and the density-density correlations are shown in figure 2. We can see here that, as expected within our perturbative scheme, the fourth-order terms give rise to very small corrections in the dynamics.

Perturbative treatment of the radiative decay
Additional processes can in principle be included in this treatment as long as they do not violate the separation of time-scales. For instance, spontaneous radiative decay from the excited (|↑ ) to the ground state (|↓ ) with a rate Γ ryd can be modelled via a Markovian dissipator which commutes with the one in (1), i.e., Therefore, the typical time-scale 1/γ due to dephasing is unaffected and one can still analyse the projected dynamics in the corresponding stationary (diagonal) ensemble. However, decay does not commute with the interaction term and thus introduces an evolution which is not easy to account for analytically. This issue can be overcome by also considering the decay as a slow process or a perturbation, i.e. Γ ryd ≪ γ. Up to the first non-trivial order in both the decay rate and the coherent flipping amplitude Ω the effective rate equation readṡ with Γ k as in (4). Note that Γ ryd in (9) must be positive and therefore (10) constitutes a proper classical master equation.

Three-level Rydberg atoms in a EIT configuration
The second scenario we consider is a frozen gas of N atoms with three internal levels subject to EIT conditions, as shown in figure 1(b). Here the ground state |↓ is resonantly coupled to an intermediate state |← via a laser field with Rabi frequency Ω p . A second laser couples |← to |↑ with Rabi frequency Ω c and detuning ∆. Once again, we only account for interactions V km between pairs of atoms in the Rydberg state.
Overall, the Hamiltonian that governs the coherent evolution of this many body system can be expressed as H = H 0 + H 1 , with where n k denotes the occupation number of the k-th Rydberg level, and Atoms excited to a Rydberg level are typically quite stable and display mesoscopic lifetimes of the order of tens of µs [2]. Thus, on microscopic time-scales the main process causing loss of energy is spontaneous radiative decay of the intermediate state |← to the ground state |↓ , which occurs with rate Γ. We model such a source of dissipation as where we have again neglected spatial and temporal correlations, i.e., each atom decays independently of the state of the others.

Second order effective evolution
We assume now that Γ is much larger than both Rabi frequencies Ω c and Ω p . In this limit the population of the intermediate state |← will decay on a fast timescale Γ −1 and thus one can adiabatically eliminate it. We can then describe the system's state by means of a reduced density matrix µ which includes only the two internal states |↓ and |↑ . Note that in this case coherences between the Rydberg and ground states are preserved. Despite the fact that a classical interpretation is no longer possible, this approach yields a considerable reduction in the growth of the Hilbert space dimension with the system size (from 3 N to 2 N ). This can prove useful for numerical approaches focussing on the aforementioned subspace. In this case the observables one is effectively restricted to are those which can be written as combinations of A detailed discussion on how to implement this approximation can be found in Appendix C.
First, let us note that H 0 is entirely written in terms of the operators n k = |↑ k ↑ k |, which are not directly affected by the dissipation. Hence, defining However, in contrast with the previous case, the stationary subspace of D is not entirely included in the one of H 0 , which implies that H 0 generates a nontrivial dynamics within it. Because of this fact, it becomes more difficult to account for its presence (although it is still possible to do it analytically, as we mention at the end of Appendix C). In the following, we shall treat both H 0 and H 1 perturbatively, which yields, up to second order, the effective equation for the reduced density matrix μ This is a Lindblad quantum master equation with coherent part governed by H 0 and dissipation provided by the jump operators with p k = |↓ k ↓ k | being the projector onto the ground state in the reduced space. Note that here, whilst the coherent part H 0 simply acts onto the "classical" configurations by associating a given energy to each of them, it is the jump operators that tend to bring the system into coherent superposition states. The engineering of such a type of dissipation, namely one that leads to dynamics and stationary states featuring quantum coherence and many-body superpositions, has been attracting an increasing amount of attention recently [46,47,48,49,50].
In order to assess the validity of this effective equation of motion, we resort again to the numerical simulation of a one-dimensional system on a periodic lattice of spacing a and van der Waals interaction, with the nearest-neighbour interaction denoted by V = C 6 /a 6 . We have numerically solved both the full quantum many-body and the second order effective dynamics for small systems of up to N = 6 atoms. In both cases we have employed a numerically exact direct integration of the corresponding master equations, which quickly becomes very demanding in terms of computational resources as the system size is increased. In figure 4 we display the time evolution of the expectation value of the Rydberg density and its fluctuations for the resonant (∆ = 0) case, Γ = V = 100Ω c and three different values of Ω p /Ω c = 0.1, 1, and 10. As discussed above, in contrast with the previous two-level system [figure 1(a)] here we keep also track of the dynamics of some off-diagonal observables, such as the operator σ x = k σ x k /N (see third column in figure 4). Moreover, we display not only the shorttime behaviour but also the long-time one, as in this case the stationary state is not trivial and generically depends on the parameters of the system (in particular on the Rabi frequencies Ω c and Ω p ). As a consequence, in this case the steady state will in principle only be reproduced in a perturbative fashion. We observe in general good agreement between the approximate and the exact results. In particular, the curves are hardly distinguishable for Ω p /Ω c = 0.1 and 1. Small deviations show up instead in the case Ω p /Ω c = 10, which are related to the fact that Γ is only 10 times larger than Ω p and therefore we are approaching the limits of applicability of our perturbative scheme. As long as one chooses parameters in a range compatible with the latter, however, we can conclude that (11) provides a good approximation for the description of the dynamics of an interacting Rydberg gas under EIT conditions.

Nearest-neighbour exclusion
A further simplification can be obtained in a particularly simple case: a one-dimensional chain of atoms with resonant excitation (∆ = 0) where we approximate the interaction as a hard-wall repulsion between neighbouring excitations, i.e., This effectively yields a projection of the dynamics onto the set of ground states of H 0 , i.e., the portion of the Hilbert space spanned by classical states without neighbouring pairs of excitations. This nearest neighbour exclusion approximation has been often used for gaining insight on the underlying physics of strongly-interacting Rydberg gases [51,52,25,23]. One can show that, for any initial condition with overlap only on allowed (zeroenergy) states, the effective equation for the dynamics of the reduced density matrix µ has a purely dissipative form [40], i.e. with Here we have introduced the operator P k = p k−1 p k+1 , which indeed ensures that an excitation |↑ is never created next to an already existing one. In figure 5 we assess the validity of this approximation by numerical methods. This time we use as a measure the trace distance of two density matrices ρ and µ In particular, we calculate the trace distance between the stationary state of the full three-level many-body system with nearest neighbour interaction V , ρ ss , and the one obtained from the approximate equation (12), µ ss . In figure 5 we plot T (ρ ss , µ ss ) as a function of the system size N for Ω p /Ω c = 10 and Γ/Ω c = 100 [panel (a)] and 1000 [panel (b)], as well as for different values of the interaction V . First, we observe that the validity of the approximation appears to get slightly worse as the system size is increased. Secondly, while for low values of the interaction V the approximation is poor (as expected), the trace distance rapidly reaches a saturation value when increasing V . Finally, by comparing panels (a) and (b) we infer that this saturation value tends to vanish when the expansion parameters Ω p/c /Γ of the perturbative series are made smaller, which is indeed compatible with the fact that we expect the stationary state to be only perturbatively reproduced to second order in our treatment.

Conclusions
In this paper, we have applied the Nakajima-Zwanzig projection technique to the study of strongly-interacting many-body dynamics, particularly in the context of Rydberg gases. By relying on a time-scale separation between a fast and a slow dynamics, effectively integrating out the fast degrees of freedom and using perturbation theory, we have obtained effective equations of motion that approximately describe the dynamics of a selected set of observables of the system within a reduced subspace. This yields a reduction of the complexity of the corresponding problem, which allows for a numerical treatment of larger systems. Via numerical simulations of small systems we have verified that the obtained effective dynamics yields indeed an excellent approximation to both the stationary state and the relaxation towards it. We have focussed in particular on two models which describe many-body systems which are currently intensely investigated and experimentally realised with stronglyinteracting Rydberg gases. The first one considers each atom as a subsystem with only two physically relevant levels coupled by a laser and subject to strong dephasing. Here, we have found that a classical effective description of the problem is possible, thereby making very large systems amenable to numerical treatment. In the second case we have considered a Rydberg gas under EIT conditions where the fast evolving timescale is provided by the rapid decay of the intermediate level. In this latter case we have also obtained a reduction of the complexity of the problem, whilst not as significant as in the previous case, since part of the Hilbert space structure is retained.
The main aim of this work is to provide a formal framework for the effective description of these strongly interacting systems in the limit of strong dissipation. We hope this effort to contribute towards unifying different results already obtained in the literature (e.g. [20,21,28,40]) and, moreover, provide some degree of guidance and a reference for future efforts employing these techniques in the context of interacting Rydberg gases.

Appendix A. Adiabatic elimination of fast degrees of freedom
In this appendix we establish the general formalism and notation we have employed to derive the reduced dynamical equations in the main text. It is based on the Nakajima-Zwanzig projection formalism [18,19], which relies on finding a criterion to divide the degrees of freedom in relevant and irrelevant, and effectively keeping track only on the former. In a quantum setting, the most intuitive application of this frame of thought would be to focus the attention on a subsystem, while treating the rest as an effective "external" bath [15,16]. Here, however, we shall take a slightly different perspective and hinge instead upon a clear time-scale separation for different dynamical processes.
Let us consider a general system whose Markovian dynamics is described by the von Neumann equatioṅ where L 0 and L 1 are two time-independent Liouville operators acting on the density matrix ρ. Our aim is to project the evolution onto a reduced subspace by adiabatically eliminating the fast degrees of freedom, which we assume to be entirely described by L 0 . We also assume for simplicity that L 0 includes a dissipative part. In other words, we are trying to obtain a coarse-grained equation of motion which effectively captures the dynamics of the system on time scales longer than the typical ones of L 0 . Within this frame of thought, we introduce the projector onto the stationary subspace of L 0 (i.e., its null eigenspace, or kernel). The existence of the limit is ensured by the fact that positivity must be preserved by L 0 and therefore all its eigenvalues must have non-positive real part. The reduced density matrix whose dynamics we want to describe is thus µ = P ρ. We correspondingly define the complementary projector Q = 1−P , where 1 is the identity superoperator, and χ = Qρ. Applying P and Q to (A.1), we can rewrite it as where we have used the fact that, by construction, P L 0 = L 0 P = 0. An implicit solution of the first equation is given by This allows us to write down an integro-differential equation for µ which does not depend on χ(t), but only on its initial value. If we further assume that the initial condition entirely lies within the kernel of L 0 , i.e., χ(0) = 0, we can rewrite the equation for µ aṡ The equation above is still exact, but its integro-differential nature makes it difficult to approach analytically. We thereby proceed by applying a Laplace transform L[•] = ∞ 0 dte −st (•), as it readily yields a perturbative expansion. This leads to where we have employed some of the general properties of Laplace transforms, i.e., where f * g denotes the convolution t 0 dτ f (τ )g(t − τ ). By assuming that the amplitude of L 1 is much smaller than the other energy scales (which makes the corresponding dynamics much slower), we can expand and subsequently truncate the fraction in (A.4) as a power series in L 1 : Note that, since both L 0 and the sum L 0 + L 1 are assumed to be proper evolution superoperators, their spectra (and thus the singularities of the Laplace transform above) lie to the left of the imaginary axis. It is therefore possible to easily choose a contour on which to define the Laplace anti-transform. By exploiting the same general properties seen above, one can rewrite (A.3) in powers of L 1 aṡ where τ 0 ≡ t.
In order to obtain a differential equation (i.e., an expression forμ(t) which only depends on µ(t)), we perform the trivial substitution µ(τ α−1 ) = µ(t) and express the second addend as Note that the difference above is of the same order of the derivative, i.e., at least O(L 1 ), and that it can be made time-local to any finite perturbative order by iteration: for instance, up to second order it can be cast as This constitutes substantially an adaptation to the present case of the time convolutionless technique (see e.g., [15]). In fact, we can iteratively apply (A.6) to obtain a time-local representation of the corrections to higher orders in the expansion. This leads to a redefinition of the evolution operators L (α) ; the first four orders read where we have applied the set of changes of variables τ k+1 → τ k − τ k+1 and τ 1 → t − τ 1 . Note that the projectors Q in front of each exponential e τ k L 0 ensure that the latter can only act on states not belonging to its kernel. We furthermore assume that there are no eigenvalues of L 0 which are purely imaginary. This implies that Q projects onto eigenspaces corresponding to eigenvalues with strictly negative real parts, such that the action of e τ k L 0 introduces an exponential dampening typically dictated by the eigenvalue λ with the largest non-trivial real part. Our assumption of a clear time-scale separation implies that −ℜ(λ) must be large compared with the typical energy scales of L 1 or, more precisely, larger than those associated to operators which couple the stationary subspace of L 0 to its complement. In the light of this, extending the integration domains to the whole real axis should introduce only a small correction. After this approximation, the final form of the terms up to fourth order is Before applying these expressions to the specific systems mentioned in the main text, let us briefly discuss which observables can be calculated within this scheme. Exploiting the fact that the trace defines a scalar product in the superoperatorial space (ρ, σ) ≡ Tr ρ † σ we can rewrite the relation above as which should be valid for every possible choice of the density matrix ρ. Hence, we conclude that only observables that satisfy O = P † O can be calculated within the reduced-space formalism discussed here.

Appendix B. Two-level Rydberg atoms in the presence of strong dephasing
In this Appendix we will give detailed account of the derivation of an effective equation of motion of a system of N driven two-level atoms strongly interacting and in the presence of strong dephasing, as described in Section 2.
The dynamics of the system is modelled via the master equation V km n k n m and the dissipator Dρ = γ k n k ρn k − 1 2 {n k , ρ} . The stationary space of L 0 is formed here by all matrices µ which are diagonal in the basis of all possible classical spin configurations. In the following, we will always use the terms "diagonal" and "off-diagonal" referring to this basis. Strictly speaking, the fast dynamics is provided by the dephasing and thus one would have to define L 0 = D in order to directly connect to the results of Appendix A. However, in this case the commutator −i [H 0 , •] not only commutes with the dephasing dissipator, but actually its stationary subspace includes the one of D. Therefore, this term can be conveniently included in L 0 without altering the structure of said subspace. On a different note, the perturbation L 1 does not connect states belonging to the kernel of L 0 , i.e., it cannot map any state which is stationary under the action of L 0 into another one. This implies that P L 1 P = 0, and hence the first four terms of the expansion are notably simplified to Let us first calculate the second order contribution to the dynamics. To do so, and in order to be able to compute any term in the perturbative expansion, it is fundamental to understand how the operator O = ∞ 0 dτ e τ L 0 QL 1 acts on a diagonal matrix. Note that, according to (B.1), one can actually reduce this problem to studying the action of the single-site components We start by decomposing the diagonal matrix µ as matrices acting on all sites but the k-th one. This representation allows us to more easily calculate the action of L k 1 on µ, which reads As this matrix is entirely off-diagonal, the operator Q effectively acts as the identity when applied to it. We then have now to compute the action of the superoperator e τ L 0 on a generic matrix of the form µ k 1 ,...km ι 1 ,...ιm ⊗ σ y k 1 ⊗ . . . ⊗ σ y km , with ι n =↑, ↓ and k n = 1 . . . N. To this end, we first notice that the action of the Hamiltonian H 0 and the dissipator D commute, which allows us to factorize the exponential as which amounts simply to the multiplication by a damping factor. The action of the coherent part is more involved. In order to give an expression for it as well, we divide the Hamiltonian H 0 into two operators: one which does not depend on the indices k 1 . . . k m and therefore inconsequentially commutes with the density matrix above, and a part which instead depends on them, i.e., The action of the Hamiltonian on a matrix of the form (B.3) can be obtained by the realisation that h k 1 ,...km (n k 1 , . . . , n km ) σ + k j = h k 1 ,...,km n k 1 , . . . , n k j = 1, . . . n km σ + k j and σ + k j h k 1 ,...,km (n k 1 , . . . , n km ) = σ + k j h k 1 ,...,km n k 1 , . . . , n k j = 0, . . . , n km , We emphasize here that h k 1 ,...,km (n k j = 1) stands for n k j = 1 h k 1 ,...,km n k j = 1 and consists of a reduced matrix which does not act on the k j -th subspace. Note that this implies that it commutes with every local operator acting only on it, e.g., σ x k j . Introducing the function N p with p = ± such that N + = 1 and N − = 0 we find ) With these expressions we can already calculate the contribution to the second order. To do so, first we realize that the action of O k (τ ) on the diagonal matrix µ yields an off-diagonal form. As a consequence, due to the presence of a projector P as a last step, one needs the subsequent action of L k 1 σ ± k = ±iΩσ z k -i.e., specifically of the k-th component of L 1 -in order to recover a diagonal matrix and get a non-vanishing outcome. Thus, the second order contribution yields which after the integration over time reads as shown in (3), where we have used that h k (0) = 0 and h k (1) = ∆ + q =k V kq n q .
Let us now look into the calculation of the next orders. Here, we again note that the action of the operator O on the diagonal matrix µ yields an off-diagonal one. The action of L 0 does not modify the matrix structure, and thus only the subsequent action of L 1 can recover a diagonal matrix that is not annihilated by the final application of the projector P . This in turn means that any odd number of L 1 s applied to µ will always render an off-diagonal matrix, yielding and, in particular, L (3) µ = 0. Moreover, as seen above, for every occurrence of e.g., L k 1 , a second k-th component must be present, since no other L q 1 with q = k acts on the k-th subspace and is able to recover the diagonal structure. Hence, the different components of L 1 always appear in pairs.
We look now into the calculation of the fourth order contribution to the perturbative expansion. We first split it into two terms as Let us focus on the first term, A 4 µ: from the discussion above, we know that the only non-zero terms are the ones where the L 1 operators come in pairs with equal indices, i.e., On the other hand, the off-diagonal projectors Q included in each O (note that O = QO) prevent the matrix structure from being diagonal at any intermediate step before the last one, which means that any subsequence of Os which appears on the right (i.e., directly acts on µ), is strictly shorter than the full sequence, and in which all indices of the L 1 components can be paired up with each other identically vanishes. Thus, (B.7) can be simplified to We calculate now step by step the action of O m (τ 2 )O k (τ 1 )µ, common to both terms in (B.8). The first step, O k (τ 1 )µ, we already calculated for the second order in (B.6).
We have now to apply O m (τ 2 ) = e τ 2 L 0 QL m 1 . The action of L m 1 here does not involve solely the difference µ k − ≡ µ k ↓ − µ k ↑ , but also each operator-valued prefactor. In order to calculate its action we employ the same decomposition used in (B.2) and rewrite µ k − as We now make use of the general identity which, recalling that we are defining Hence, the overall effect of O m (τ 2 ) reads where we have applied (B.5) to each addend.
Let us now calculate the first term of the sum in (B.8), which means that the next step involves the action of a superoperator O k (τ 3 ). Following the same procedure outlined above, applying the identities L k 1 (M ⊗ σ ± k ) = ±iΩM ⊗ σ z k and introducing the shorthand notation one obtains after some algebraic manipulation As a final step for the calculation of this first term, we have to apply L m 1 , which renders In order to calculate the second contribution to A 4 in (B.8), we need to go back to (B.10) and apply O m (τ 3 ) to it, thus obtaining Finally, we apply the operator L k 1 and arrive at where we have used the definition (B.11). Note that the remaining time integrations involve only diagonal matrices (in the sense defined at the beginning, i.e., in the classical basis) and can be therefore straightforwardly evaluated. We thus obtain where the first summand comes from the first term in (B.8) and the other two from the second one. In this expression the operator-valued rates R km i and R ′ km as already shown in (5) in Section 2 of the paper. Note that, as stated above, the indices of the operators "h" in these expressions denote the subspaces on which they do not act and, as a consequence, [h k (1), σ x k ] = 0 and h k,m (n k , n m ), σ x k/m = 0. This property then trivially transmits to the corresponding Γ-s. Therefore, as reported in the main text, R km 2/3 and R ′ km 2/3 commute with σ x k . The calculation of the second term of (B.8) follows analogous steps; the final expression reads Note that, from the second order calculation, we already know the action of P L k 1 O k (τ 1 ) on µ, which is indeed diagonal. For the next step (the action of L m 1 O m (τ 2 )), we distinguish the cases m = k and m = k. Let us now first consider the case m = k, which yields On the other hand, when m = k we obtain Finally, integrating over time the last two expressions we obtain the contribution for B 4 µ to fourth order , which indeed coincide with the expressions shown in (5) in Section 2 of the paper, as one can check via the definitions (B.13).

Appendix B.1. Radiative decay
We now aim to include the effect of radiative decay from the Rydberg state with rate Γ ryd , described by the dissipator D dec ρ = Γ ryd k σ − k ρσ + k − 1 2 {n k , ρ} . The dynamics of the system are thus described by the von Neumann equatioṅ where we set the decay rate to be much smaller than the dephasing rate (Γ ryd ≪ γ), so that we can consider D dec as a perturbation together with the coherent driving represented by L 1 .
By following the same procedure as in Appendix A we can now writė For simplicity, we now restrict ourselves to the lowest non-trivial order in both processes, i.e., the decay (expansion in powers of Γ ryd ) and the coherent spin-flipping (expansion in powers of Ω). We also consider the order Γ ryd Ω 2 to be negligible, which allows us to disregard the corrections coming from the substitution µ(τ ) → µ(t) in the expression above. Since the dephasing and decay dissipators commute, our calculation can still hinge on the fact that the long-time behaviour outside of the classical subspace will portray an exponential decay ≈ e −tγ/2 (or a faster one). Thus, after the change of variables τ → t − τ , we can also bring the upper bound of the integral to infinity, which yieldsμ Note that the second addend is simply L (2) . Thus, we can straightforwardly obtain the expressionμ up to first order in Γ ryd and second order in Ω. This equation is equivalent to (10) in the main paper.

Appendix C. Three-level Rydberg atoms in a EIT configuration
In this Appendix we will give detailed account of the derivation of an effective equation of motion of a system of N driven three-level atoms which display strong interactions between excited states and in the presence of fast decay processes from the intermediate to the ground state, as described in Section 3. The dynamics of the system is described by the master equation (A.1), where Note that H 0 and D act on different subspaces and therefore their actions on the state of the system commute. Therefore, once again, we can include the commutator with H 0 in the "fast" term L 0 while considering only D for the determination of the stationary subspace.
While in the previous case the action of the projector P was equivalent to a projection onto the diagonal of the density matrix, here its action is slightly more involved. In order to obtain it, we use the fact that e tD = e t k D k = k e tD k .
Each D k non-trivially acts on the k-th (three-dimensional) subspace and its action on a generic matrix A can be represented as and then define for this case µ = P ρ.
We are now in position to calculate the effective dynamics for µ up to second order of perturbation in L 1 using the results in Appendix A L (1) µ = P L 1 µ(t) L (2) µ = P L 1 ∞ 0 dτ e τ L 0 QL 1 µ(t).
Let us start by calculating the first order correction. Here, we separate the action of L 1 into the one associated to H 0 and to H 1 . The latter can be written as a sum H 1 = k H k 1 and we can therefore restrict our analysis here to a generic (k-th) subspace. In particular, H k 1 acts on an element of the kernel of D as which constitutes a matrix orthogonal to the kernel of D and, hence, implies that P H k 1 , µ = 0 and Q H k 1 , µ = H k 1 , µ . Thus, the first order contribution to the effective equation of motion reads We now use this knowledge as well to calculate the second order contribution L (2) µ. The first thing we realize is that the presence of a projector Q after the application of L 1 to µ leaves only the contribution of H 1 , which we know already from (C.2), as [H 0 , P ] = 0. The next step is the application of e τ L 0 = k e τ D k to a matrix of the form (C.2). From (C.1) we can extract that this action amounts simply to a multiplication of the matrix by e −τ Γ/2 . The last step is thus the application of L 1 to the matrix (C.2). First we realize that, as the interaction Hamiltonian H 0 keeps the matrix within the subspace orthogonal to the kernel, its contribution vanishes as P is subsequently applied. Hence, we only need to understand the action of H 1 , which yields Note that this matrix does not generally belong to the kernel of D, and it is only after applying P k that one gets Thus, one can obtain now the final form of the second order contribution, which yields with ǫ k = µ k ↑↓ + µ k ↓↑ and δ k = µ k ↓↓ + µ k ↑↑ , and where we have eliminated the intermediate level and hence used a 2 × 2 matrix for the description of the k-th atom.
Note that this contribution can be also written out as a purely dissipative Lindblad equation of the form where the jump operators have a non-classical form where p k = |↓ k ↓ k | and σ k − = |↓ k ↑ k | are spin-1/2 operators. It is worth mentioning that in this simple case one can actually treat the "classical" part of the Hamiltonian H 0 in a non-perturbative fashion. This can be done employing an interaction representation for L 0 and L 1 : L 0 → L 0 (t) = e −tH 0 L 0 e tH 0 and L 1 → L 1 (t) = e −tH 0 L 1 e tH 0 .
One then finds L 0 (t) = L 0 and L k where h k (1) = ∆ + q =k V kq n q is the same object defined in (B.4). The procedure outlined in these appendices can be then carried on in a similar manner; the main difference being that the exponentials e t(L 0 +L 1 ) must be replaced by their time-ordered counterparts T e t 0 dτ [L0+ L 1 (τ )] . At second order in L 1 one eventually findṡ .
We have verified numerically that this expression generally yields negligible corrections to the dynamics with respect to (C.3) in the perturbative regime Γ ≫ Ω c/p .