Bistability vs. Metastability in Driven Dissipative Rydberg Gases

We investigate the possibility of a bistable phase in an open many-body system. To this end we discuss the microscopic dynamics of a continuously off-resonantly driven Rydberg lattice gas in the regime of strong decoherence. Our experimental results reveal a prolongation of the temporal correlations with respect to the lifetime of a single Rydberg excitation and show strong evidence for the formation of finite-sized Rydberg excitation clusters in the steady state. We simulate our data using a simplified and a full many-body rate-equation model. The results are compatible with the formation of metastable states associated with a bimodal counting distribution as well as dynamic hysteresis. A scaling analysis reveals however, that the correlation times remain finite for all relevant system parameters. This suggest that the Rydberg aggregate is composed of many small clusters and all correlation lengths remain finite. This is a strong indication for the absence of a global bistable phase, previously suggested to exist in this system.


I. INTRODUCTION
Recent progress in experimental control offers unique possibilities to study the interplay between pure quantum mechanical systems and their coupling to reservoirs. In particular the competition between external drive and dissipation can result in interesting phases in the steady state [1][2][3][4]. One fascinating avenue is to study the properties of these phases and the possibility of phase transitions in the context of an open quantum system [1, [5][6][7][8]. In addition, including strong and nonlocal interactions may allow to prepare and study strongly correlated many-body systems [9][10][11].
Rydberg atoms are particularly interesting as they provide strong and quasi-long range interactions in an inherently dissipative environment. By means of an optical drive, ground state atoms can be excited to states with high principal quantum numbers. These Rydberg states have a large polarizability resulting in a strongly interacting quantum gas. Therefore, these systems are an ideal platform to study non-equilibrium many-body physics with strong dissipation [2,[12][13][14][15][16][17].
Depending on the driving scheme to the excited Rydberg state and the geometry of the atomic ensemble, Rydberg gases can exhibit very different phases. In the case of the so-called blockade regime, i.e. for resonant driving, the strong interaction shift of Rydberg states induced by an already present Rydberg excitation suppresses further excitations in a mesoscopic ensemble of surrounding atoms. This leads to the concept of a Rydberg superatom [18][19][20][21][22] and was initially proposed to realize quantum gates [23]. Here, interesting phases with antiferromagnetic, crystalline long-range order have been predicted [2] and phases with finite-length spatial correlations have been observed [24]. In the case of the so called anti-blockade regime, i.e. for off-resonant driving, the interaction shift may be compensated by detuning the frequency of the excitation laser from resonance. Now, instead of suppressing excitations within a certain volume, further Rydberg excitations are facilitated at a specific distance from an already excited atom. This leads to the formation of Rydberg excitation clusters [15,25,26].
One interesting aspect observed in recent anti-blockade experiments using a full counting analysis is a bimodal excitation number distribution [27,28] and a hysteresis in the excitation dynamics [13,17]. It has been under debate whether this indicates a transition to a bistable steady state of the open many body system [8,9,13,16,29,30]. We will show that both, hysteresis and bimodal counting statistics, are features of metastable states.
A metastable state manifests itself in a separation of timescales, where the relaxation time T is much larger than all other timescales in the system [6,31]. This may lead to a hysteresis behavior upon temporal changes of system parameters on a time-scale τ (sweep time) small compared to T . When the sweep time becomes much larger than T the hysteresis disappears. Associated with a finite value of T is an upper limit for all spatial correlation lengths, ξ ∼ T . We will argue that metastable systems with correlation length similar to the system size allow for a bimodal distribution. If T increases with the system size L and eventually diverges in the thermodynamic limit, the metastable states become degenerate, truly stationary states of the system. In this case the steady state is called bistable (multistable). As will be shown later bistability in an infinite, translationally invariant system implies long-range spatial order. Thus, different from metastability, the emergence of bistability can be associated with a phase transition to an ordered state.
In this paper, we experimentally investigate and theoretically model the excitation dynamics of a Rydberg lattice gas with off-resonant continuous driving and strong dissipation (see Fig. 1a). We continuously monitor the laser-induced excitation and subsequent ionization of Rydberg atoms and analyse its temporal correlations ( Rabi frequency Ω and detuning ∆ > 0. The Rydberg atoms are ionized at the rate Γion, which is much smaller than the natural decay rate. Therefore, the ions serve as a probe of the excitation dynamics. The color scheme of the atoms is chosen according to the rates in (c): Rydberg excitations are red, ground state atoms are gray and atoms, which can get facilitated are blue. (b) Pair correlation function g (2) (τ ) deduced from the ion signal (Ω = 2π × 77 kHz, ∆ = 2π × 13 MHz). The correlation function shows super Poissonian fluctuations with large bunching amplitude g (2) (0) 1 and relaxation time τ (2) τsp , where τsp = 20 µs is the natural lifetime of the excited Rydberg state 25P 1/2 . Inset: measured ion signal for a single (blue) and averaged over 40 (gray) experimental realizations. The shaded area around the gray line denotes Poissonian fluctuations. (c) Effective cluster model: A single seed excitation with rate NatΓ seed initializes the growth of a cluster. While a facilitation rate Γ ↑ Γ seed leads to an increase of cluster size m, the spontaneous decay rate Γsp limits the size of a cluster.
indicate the formation of small clusters with a long but finite lifetime, which originate from a correlated excitation cascade. We compare the experimental results to numerical simulations based on a simplified single-cluster model (Fig. 1c) and a more advanced rate equation model of extended systems. To check for the existence of bistability we perform a finite-size extrapolation in the numerical simulations. Our results indicate that the correlation time and length remain finite in the experimentally relevant regime amenable to a rate equation description. The bimodal character of the full counting distribution [28] disappears when the observed system size exceeds the correlation length, while a dynamic hysteresis can still be observed provided the sweep time is sufficiently small [32]. This suggests that the Rydberg aggregate in the anti-blockade regime consists of many small independent clusters and is incompatible with a global bistable regime.

II. BISTABLE AND METASTABLE STEADY STATES
Before discussing Rydberg lattices gases let us start with a general introduction to bistability and metastability in open quantum systems. Quantum optical systems coupled to external reservoirs and subject to external drives can often be described in terms of a Lindblad master equation for the many-body density matrix ρ [33][34][35] The Lindbladian superoperator L is the generator of the dynamics, determined by the Hamiltonian H describing the unitary evolution and the jump operators L µ describing the Markovian reservoir couplings L has complex eigenvalues λ k and corresponding right eigenvectors ρ k , i.e. Lρ k = −λ k ρ k . The non-negative real parts of the eigenvalues, Re[λ k ] ≥ 0, determine the characteristic relaxation rates towards a stationary state ρ 0 or a manifold of stationary states {ρ (n) 0 } with λ 0 = 0. It should be noted that all "excited" eigenvectors ρ k with λ k = 0 are no density matrices since Tr{ρ k } = 0. If the steady state is unique, it is called stable, if two orthogonal stationary solutions, ρ (1) 0 and ρ (2) 0 , exist, the steady state is called bistable. In the first case any initial state will eventually evolve into the unique steady state, in the second case the asymptotic state for large times is a mixture of the two solutions with coefficients determined by the initial state ρ(t = 0): is the left eigenvector of L corresponding to ρ (n) 0 and Tr{ρ In the cases we are interested in here, both the Hamiltonian H, as well as the jump operators L µ are local, i.e. they can be written as a sum of terms which act only locally. In such a case the steady state is generically unique if the system is finite unless the reservoir interactions are fine-tuned. General conditions for the uniqueness of steady states can be found e.g. in [36][37][38].
Since experiments can be performed only in finite time it is interesting to ask if states exist that, while not being true stationary states, appear to be stable over long but finite time scales. These states are called metastable and they occur in systems with a separation of time scales, i.e. where there is initial relaxation into long-lived states with the subsequent decay into the true stationary state on a much longer time scale [6,31]. Metastability is reflected in the spectrum of decay rates, Re[λ k ], as eigenvalues with small real part, clearly separated from all other eigenvalues, see Fig. 2a.
Bistability typically emerges only in the thermodynamic limit of infinite system size L → ∞, see Fig. 2b. An example in the context of driven Rydberg lattice gases is the formation of anti-ferromagnetically ordered steady states under conditions of Rydberg blockade [2]. Here a phase transition from an unordered steady state into one with two different checkerboard-like distribution of Rydberg excitations occurs above a certain critical value of the optical driving associated with a spontaneous breaking of translational lattice symmetry. For any finite system with linear dimension L < ∞ there is a unique steady state which is an equally weighted mixture of the two checkerboard configurations. Each of these two configurations corresponds however to a metastable state, which, as shown in [6] are admixtures of the true steady state ρ 0 and the first excited eigenvector ρ 1 . The metastable states eventually relax to the true steady state on a time scale Approaching the thermodynamic limit there is a critical slow down and the system becomes truly bistable.
FIG. 2. Generic spectrum of relaxation (damping) rates as a function of system size L for a) a metastable and b) a truly bistable system.
The transition from a unique to a bistable steady state is a true phase transition, which can be characterized by an order parameter. Bistability and metastability cannot be distinguished, however, in finite systems and one has to investigate the scaling behavior of the characteristic timescales for L → ∞. We will show now that bistability is always associated with long range order. To this end we consider a translationally invariant system and assume that two orthogonal steady states, ρ (1) 0 and ρ (2) 0 , exist. Then the asymptotic state at large times is given by eq.
(3). In order to detect bistability there must be a local observablex j that distinguishes the two states. Here the index j labels positions (or some compact part of space) wherex j acts on. Without loss of generality we assume Making use of translational invariance we find in the asymptotic state ρ ss , eq. (3): Obviously Ŷ ≥ 0 in any state. In both steady states, ρ (µ) 0 one finds j,k x jxk µ ≥ j x j µ k x k µ = L 2 . This yields for the general asymptotic state ρ ss j,k Now, considering correlations of the local observable Unless p = 0 or 1, this inequality can only be fulfilled in the thermodynamic limit L → ∞ if Thus bistability is only possible if there is long-range order. This result is very intuitive. If only two (or some finite number) of stationary states exists, which differ in the expectation value of some local observable, there must be long-range correlations in the system since one end of the system has to "know" in which state the other end is. For any system which is bistable in the thermodynamic limit we thus expect a scaling of the correlation length ξ ∼ L β>0 .
Since the speed at which correlations can spread in the system is finite, the transition into a bistable state is associated with a critical slow down, i.e. with characteristic time scales that diverge with the system size. If the correlation length ξ is finite, the number of stationary configurations is 2 N with N ∼ ξ/L. For L ξ the central limit theorem holds and global probability distributions become single-peaked Gaussians. We will argue in the following that under realistic experimental conditions correlation length and characteristic time scales in off-resonantly driven Rydberg lattice gases remain finite and there is no bistability in this system.

III. EXPERIMENT
Experimentally, we continuously probe the Rydberg excitation dynamics in the anti-blockade regime. To this end, we prepare N at = 20 000 87 Rb atoms in a 3D optical lattice in the Mott insulating phase at unit filling. The lattice constants in x-and y-direction are a x,y = 374 nm and in z-direction a z = 529 nm. We couple ground state atoms in the 5S 1/2 state continuously via a one-photon transition at a wavelength of λ = 297 nm to the excited Rydberg state 25P 1/2 with Rabi frequency Ω and blue detuning ∆ > 0. See Fig. 1a for a sketch of the experiment. The spontaneous decay rate from the Rydberg state is given by Γ sp = τ −1 sp = 50 kHz and we estimate a bare decoherence rate of γ 0 300 kHz which primarily originates from the laser linewidth.
The dipole trap lasers provide a weak photoionization rate Γ ion 2 kHz, which we use to observe the excitation dynamics in the system. The ions are guided towards a detector by a small electric field of 90 mV cm −1 . The detector efficiency is 40%. The retrieved ion signal with rate R ion (t) = Γ ion N Ry (t), corrected by the detector efficiency, serves as a weak probe for the full excitation number N Ry (t). We experimentally determine the Rabi frequency by analyzing the temporal statistics of the first detected ion. A more detailed discussion of the experimental setup and preparation can be found in [39] and in appendix A.
In the inset of Fig. 1b the measured ion rate R ion (t) for a Rabi frequency of Ω = 2π × 77 kHz and a detuning of ∆ = 2π × 13 MHz is shown. The recorded signal extends over several tens of ms, which is three orders of magnitude larger than the effective natural lifetime τ sp = 20 µs (reduced by black-body transitions into neighbouring states) of the 25P 1/2 state. Therefore, in contrast to previous studies [27,28], our measurements allow us to study the steady-state of the system [40].
The experimentally accessible quantities, which we use to characterize the formation and the dynamics of Rydberg clusters, are the average ion rate and the timedependent two particle correlation function R ion (t) and g (2) respectively. Already in a single experimental realization (see inset Fig. 1b) strongly bunched excitations are visible, indicating strong excitation number fluctuations. To quantify them, we extract g (2) (τ ) of the measured ion signals, see Fig. 1b. For typical parameters we find a value of g (2) (0) 1. Moreover, we find an exponential decay of the correlation signal with time scale τ (2) much larger than the lifetime of the 25P 1/2 Rydberg state. Crucially, the strong excitation bunching reflects two different timescales arising in the anti-blockade regime: While a first seed excitation is strongly suppressed (indications thereof can bee seen in the inset of Fig. 1b as large intervals of zero signal in a single experimental run), the excitation rate for a second atom can be enhanced in the presence of the first [25,26,41]. On top of that, further facilitated excitations lead to the correlated growth of Rydberg excitations which form a cluster. This leads to a large bunching amplitude g (2) (0) and a long corre- Inset: bunching amplitude g (2) (0) − 1 shows peaks at detunings which correspond to a cluster size of approximately two. Error bars correspond to the statistical uncertainty from fitting an exponential decay to the experimental retrieved correlation function. For simplicity they are not shown in the inset.
lation time τ (2) , which measures the lifetime of a single cluster τ cl , i.e. τ cl = τ (2) . The ion rate R ion and the cluster lifetime τ cl form the basis of our further analysis.

IV. SINGLE CLUSTER DYNAMICS
To gain insight into the cluster dynamics, we introduce a simplified cluster model and compare it to our experimental results. The model is capable to describe the growth of the excitation number in a single cluster and qualitatively captures all relevant features observed in the experiment. It characterizes the excitation cascade of independent, non-interacting excitation clusters, each of them being described by an effective rate equation model.
In the case of strong decoherence it has been shown that rate equation models are a good approximation [42][43][44][45] and were already successfully used to describe current experiments [15,20,27]. The effective cluster model is illustrated in Fig. 1c. It only considers the total number of excitations m and does not give insight into the internal microscopic structure of a cluster. The time evolution of the probability p m to be in a cluster of size m is governed by the equations Within the single-cluster model we can identify three main contributions to the dynamics: (i) a slow seed with rate Γ seed for generating new clusters, (ii) a fast rate Γ ↑ for increasing and decreasing the size of a cluster m and (iii) the spontaneous decay rate Γ sp = τ −1 sp limiting the size of a cluster.
The seed excitation with rate is strongly suppressed due to the large detuning ∆ in the anti-blockade regime. For a system with N at atoms, we expect independent seed events to occur on a timescale (N at Γ seed ) −1 . However, the presence of already excited Rydberg atoms strongly alters the summed excitation rate of all neighboring atoms where ∆ int is the interaction shift due to surrounding excited atoms. The competition between detuning ∆ and interaction shift ∆ int may lead to facilitated excitations [26]. The ratio between seed rate Γ seed and the rate for cluster growth Γ ↑ , already allows us to identify different regimes of correlated excitation dynamics. For Γ ↑ Γ seed , we expect a strong excitation cascade creating a cluster. In this regime each excitation triggers further excitations and we expect a substantial bunching amplitude g (2) (0). For Γ ↑ Γ seed , we expect primarily uncorrelated excitations leading to a bunching amplitude g (2) (0) 1 and clusters of size one. Finally we include spontaneous decay with rate mΓ sp increasing linearly with the cluster size m. Now, starting from an initial seed excitation, the cluster grows and shrinks with rate Γ ↑ similar to a random walk along the cluster size m ≥ 0. Although, the single-atom spontaneous decay rate Γ sp can be small, this is the limiting factor for large Rydberg excitation clusters. Microscopically, a decay event creates a defect in a consecutive chain of excited atoms and splits the cluster into smaller parts. The strong interactions prevent that these parts can merge again. Besides the geometry we account for the splittings introducing an effective coordination number z eff .
In appendix B we discuss the microscopic model in more detail.
The effective model allows us to compute quantitatively the typical sizem of a cluster. Moreover, it gives us an intuitive understanding of g (2) (0) and permits us to estimate the number of independent clusters N cl . Note, our further experimental analysis does not rely on the explicit evaluation of eqn. (13) and is thus independent of the exact interaction potential between two Rydberg atoms.

A. Cluster size
To begin with, we discuss the typical size of a cluster m. It is given by the ratio of the overall lifetime of a single cluster τ cl and the spontaneous lifetime of a single Rydberg excitation τ sp , The above expression can be understood in the following way: When a cluster containsm excitations it requires m consecutive decay events of timescale τ sp before the whole cluster vanishes. The facilitation mechanism leads to an correlated excitation growth. This is in contrast to a non interacting system, where all excitations are independent and decay on the same single atom timescale τ sp . We verified this intuitive picture by comparing full simulations to the simple cluster model and an analytic approach to the corresponding cluster size, see appendix B. We use equation (14) to extract the cluster size from the experimental data. Within the single-cluster model, the cluster lifetime τ cl and thus the typical cluster sizem only depend on the ratio Γ ↑ /Γ sp . In Fig. 3 the measured cluster sizem (corresponding to a certain cluster lifetime τ cl ) is shown for the full spectrum in the anti-blockade regime (∆ > 0) and different driving strength Ω. The data reveal a strong increase in the cluster size by increasing the driving strength Ω as one would naively expect. In the limit of large detuning ∆ ∆ int , all excitation rates are strongly suppressed and we obtain a cluster size of one corresponding to uncorrelated Rydberg excitations. In this case, the lifetime τ sp is the only relevant timescale. Close to the resonance, cascaded excitations lead to the formation of finite clusters of sizem 10. Note that for strong driving and/or small detuning, the lifetime of the entire sample approaches the cluster lifetime, and no reliable values for τ cl can be extracted. The preferred generation of clusters withm > 1 in this regime is associated with an excitation bunching. In the inset of Fig. 3 we plot g (2) (0) − 1. Experimentally, we find substantial bunching amplitudes up to g (2) (0) ∼ 13. The measurements show a peak in g (2) (0) which coincides with a cluster size ofm 2. The peak shifts towards larger detuning with increasing Rabi frequency Ω. The presence of a peak and its shift can be understood in the following way: For larger detuning, where the typical cluster size is smaller thanm = 2,  (15) for the different detunings and Rabi frequencies. For decreasing detuning the number of clusters increases steadily, exceeding the size of a cluster by roughly one order of magnitude. The color code corresponds to the same Rabi frequencies as in Fig. 3. The errorbars correspond to the statistical uncertainties from fitting the cluster lifetime τ cl as well as the multi particle seed rate NatΓ seed .
the influence of an uncorrelated background noise signal becomes important and diminishes the amplitude. For smaller detunings the cluster size as well as the rate of creating new independent clusters increases drastically. Both reduce the value of g (2) (0) due to an increase of uncorrelated ionization events.

B. Number of clusters
Next, we experimentally determine the seed rate Γ seed by an analysis of the first ionization event in each single measurement (details can be found in appendix A). Comparing the timescale with which new clusters are produced (N at Γ seed ) −1 to the lifetime of a single cluster τ cl , we can deduce approximately the number of clusters N cl The results are shown in Fig 4. Similar to the cluster sizem, the number of clusters N cl increases with increasing driving strength Ω and decreasing detuning ∆. This shows that for stronger driving and/or smaller detuning the Rydberg aggregate evolves into a steady state, which is characterized by the presence of a large fluctuating number of independent clusters with rather small size. This already suggests that in the experimentally observed parameter regime there is no global bistable phase with a large correlation length. For such a phase, we expect the system to be in a single cluster state with extent over the full system.

C. Validity of single-cluster model
The extraction of the cluster size and number from the experimental data performed in the previous section partially relies on the validity of the cluster model. In the following we make a consistency check to test the cluster model. To this end, we compare the fraction of Rydberg excitations retrieved from the cluster model ρ cl = N clm /N at = τ 2 cl Γ sp Γ seed to the excitation fraction ρ e directly extracted from the measured ionization rate R ion (t): ρ e = R max ion /(N at Γ ion ). While the first expression incorporates the equations (12) and (13) derived from the cluster model, the second solely relies on experimental data. Both quantities should be identical, independent on the Rabi frequency of the driving field and the detuning. The comparison is shown in Fig. 5. The dashed line indicates the ideal agreement ρ cl = ρ e . We see that all experimental curves obtained for different driving strength fall on top of each other. For small excitation fractions, we find furthermore good agreement between the results from the effective cluster model and the experimental ionization signal. In the strong driving regime, the experimentally retrieved ion signal results in a smaller excitation fraction than predicted by the cluster model. This may have two reasons: Firstly, we monitor the ionization signal in a single experimental run with an exponential atom loss and take the maximal ionization signal to reflect the steady state excitation fraction without atom loss. Especially in the strong driving regime, atom loss may lead to smaller maximum ionization rates reducing the extracted excitation fraction ρ e < ρ cl . Secondly, our simple single-cluster model relies on the assumption that clusters are independent. However, in appendix B we show by using a full rate equation simulation that cluster collisions become important for strong driving. Cluster collisions limit in particular the independent growth of clusters. This would also lead to a reduced excitation fraction ρ e < ρ cl . In summary we conclude that the cluster model is an appropriate description for the cluster dynamics, especially for small excitation fractions.

V. BEYOND THE SINGLE CLUSTER MODEL
The experimental protocol and the single cluster model introduced above only allow to measure and describe volume integrated quantities. To also study finite size effects, we perform in this section numerical simulations using a many-body rate equation approach, which goes beyond the simple cluster model and can be extended to large system sizes. In particular we numerically perform a finite size extrapolation of the correlation time. In contrast to the cluster model, which describes "point-like" multi level systems, here we take into account the spatial distribution of the atoms in the experiment and calculate the level shift for each atom individually. In this way, we have access to the full spatial distribution of excitations and are able to track the growth and decay of every single cluster. The microscopic description and numerical details of the system are discussed in the following section VI, while we here present the main results.

A. Finite-size scaling of the correlation time
We use the rate equation model to study the finite size scaling of the cluster lifetime. The results are shown in Fig. 6. We increase the linear system size L of a 3D Rydberg lattice gas and extract the cluster lifetime from the g (2) (τ ) correlation function for three different detunings and fixed Rabi frequency Ω = 2π × 160 kHz, using the same evaluation as for the experimental data. The data points strongly suggest a saturation in the limit L → ∞. The dashed lines in Fig. 6 correspond to an exponential saturation fit with a characteristic length scale matching the cluster sizem. This is already a first evidence that clusters, in general, do not extend over the whole system.

B. Full counting distribution and hysteresis
In a system where all timescales are finite, all correlation lengths must be finite unless the interactions are infinitely long-ranged. As a consequence the bimodality of the full counting distribution resulting from the presence of two distinguishable metastable states should only persist if the system size L is smaller than the correlation length ξ. When L exceeds ξ, averages are taken over independent regions and the bimodal distribution starts to wash out. In the limit L ξ the system behaves as N = (L/ξ) 3 independent sub-systems and the central limit theorem applies, eventually leading to an overall single-peaked Gaussian distribution. Fig. 7a shows the distribution function of the excitation number for different values of L after 1 ms excitation time. In a small system (L = 3), seed events are rare and the cluster size is limited by L. Therefore, the excitation number distribution is peaked at zero excitations. When the system size is comparable to the cluster extent (L = 6), a bimodal excitation number distribution occurs. Here, a dynamical switching between low and high total excitation number is observable. However, in the limit of large systems (L = 10) the bimodal structure disappears. Instead, a Gaussian distribution emerges with a peak at a large excitation number. Therefore, the bimodality seen in [28] may be a finite size effect. Our experimental results and theoretical investigations suggest that the correlation lengths remain short, on the order of the cluster extent and there is no phase transition to a global bistable regime. Now, let us discuss dynamic hysteresis features present in metastable systems [32]. While finite size effects vanish whenever the system size exceeds the correlation length L ξ, see Fig. 7a, a hysteresis behavior can still be seen depending on the ratio between relaxation time T and sweep time τ . Here, the sweep time τ is the duration of a continuous parameter change in the detuning ∆ from an initial to a final value. In our system, the relaxation towards the stationary state is determined by the growth of finite clusters. Therefore, we expect and numerically verified that T is on the order of the cluster lifetime τ cl . Exemplary, we show in Fig. 7b for two different sweep times τ the dynamic hysteresis. In the case of a small sweep time τ = 0.88 ms a large dynamic hysteresis area can be identified, while for τ = 88 ms the hysteresis area vanishes. The hysteresis behavior persists on a rather large timescale compared to the relaxation time, which is consistent with a simple model that describes the sweep as a sequence of successive small parameter quenches. Since the steady state excitation numberN Ry (∆) is different for different values of ∆, the coarse grained relaxation of the excitation number is no longer exponential, but attains algebraic corrections. In the limit of many small parameter quenches with approximately constant relaxation time T starting from an initial detuning ∆ 0 to ∆ 1 in time τ , the Rydberg excitation number is given by with x = (∆ 1 − ∆)/(∆ 1 − ∆ 0 ). ForN Ry = const we obtain the known exponential relaxation with timescale T . However, as seen in Fig. 7b, the mean Rydberg excitation number is typically not constant. Linearization of N Ry directly results in algebraic correction scaling with (τ /T ) −1 . This may have a tremendous effect on the hysterese relaxation. To give an example, in our system, the relaxation time is on the order of T ∼ 100 µs comparable to the cluster lifetime. However, the hysterese area vanishes only on a timescale of ∼ 100 ms, which is three orders of magnitude larger. The algebraic relaxation of the hysteresis area is consistent with other simulations [16,32] and agrees well with the findings in our system. For more details on the hysterese relaxation see appendix E.

VI. MICROSCOPIC DESCRIPTION & SIMULATION
Before we conclude, we will give some more details about the microscopic description of our experiment. For the relevant length scales we use an approximate interaction potential V (r) C 9 /r 9 with C 9 = 2π × 2.1 kHzµm 9 , which is motivated by exact diagonalization of a subspace of the two body interaction Hamiltonian including dipoledipole, dipole-quadrupole and quadrupole-quadrupole interactions [46]. For details regarding the calculation see appendix C. We do not take into account possible black-body induced transitions into neighbouring Rydberg states, which have recently been suggested to explain anomalous line broadening in dense samples [47] and which would introduce different types of interactions. While this effect probably is present on the timescales our experiments are performed on, we believe that because we continuously pump atoms into the 25P-state, the above given potential is still dominating the cluster growth dynamics. We did also check that the qualitative behavior of the system discussed here stays the same for a different type of interaction (e.g. for van der Waals-type potentials).
Besides the decoherence stemming mainly from laser noise γ 0 , we also include an inhomogeneous broadening γ inh motivated by the coupling to motional states in a lattice of finite width [48,49]. Specifically, we use an additional approximate decoherence rate where σ 60 nm is the width of the localized wavepacket in the optical lattice and the sum runs over all excited atoms j with distance r j . Therefore, the total decoherence rate is γ = γ 0 + γ inh . For our interaction potential this is a good approximation for ∆ 2π × 18 MHz.
(For a detailed discussion see appendix D.) Additionally, we take into account the ionization process by a rate Γ ion , which corresponds to a loss of atoms. Each ionization event alters the geometry of the sample of Rydberg atoms. To properly account for the detection background, we superimpose the extracted Rydberg excitation dynamics with an uncorrelated noise signal with rate Γ n . It originates from atoms which are trapped on the outer region of the optical lattice. Since the atomic density drops exponentially, these atoms do not contribute to the cascaded excitation dynamics. Although, the background noise rate Γ n may depend on the detuning ∆ and driving strength Ω, it is sufficient to approximate it with a constant value Γ n 1 kHz here. Note that all parameters used in our simulations are estimated from or determined by our experimental measurements. For a typical simulation using N at = 1000 atoms we already see strongly reduced finite size effects.
The numerical results obtained for the temporal correlation function g (2) (τ ) are shown in Fig. 8. Although, the simulations do not reproduce the experimental results on a quantitative level, we find good qualitative agreement. Firstly, the typical cluster sizes are consistent with the values obtained in the experiment and have approximately the same behavior when varying ∆ and Ω. Secondly, the structure of the numerically evaluated bunching amplitude g (2) (0) is similar to the one extracted from experiment. In particular, the g (2) (0) peak coincides with a cluster size ofm = 2 as seen experimentally. However, the actual size of g (2) (0) depends on the system size and therefore strongly deviates from the experiment. We believe this originates from e.g. the motion of the excited atoms due to the interaction or the decay to other Rydberg states due to black-body radiation. While we are aware that these effects could play a role it is numerically very challenging to incorporate all of them in a full many-body simulation. From the results obtained we do infer that the rate equation model used is nevertheless a sufficiently good approximation of the Rydberg excitation dynamics seen in the experiment.
After checking the validity of the numerical simulation, we can employ it to theoretically study the spatial cluster dynamics. In the case of a van-der-Waals interaction V (r) = C 6 /r 6 we observe directional excitation dynamics as indicated in Fig. 1a along a certain lattice direction. The directionality stems from the competition between the quasi long range interaction and the finite linewidth γ. To understand this, let us discuss the case of two neighboring Rydberg excitations on a lattice where a single excitation facilitates its nearest neighbors, i.e. ∆ = V (a). Now, for a third particle along the direction of the two excitations, the residual interaction shift in antiblockade configuration is V (a)/64, while in the orthogonal direction it is V (a)/8. Comparing to the linewidth γ, this may suppress the growth of the cluster in orthogonal directions while promoting the growth along the lattice direction. Crucially, the directionality depends on the range of the interaction. In the case of short range interaction, e.g. the interaction potential V (r) = C 9 /r 9 relevant for our experiment, the impact on the next nearest neighbors are typically negligible. Therefore, no directional Rydberg excitation dynamics is seen in the full numerical simulations using this potential. However, we believe that directional growth should be experimentally observable in different systems.

VII. CONCLUSION & OUTLOOK
To conclude, we discussed the excitation dynamics of a large dissipative Rydberg lattice gas probing the full anti-blockade regime. We showed that the increase in the relaxation time corresponds to the formation of small clusters within the Rydberg aggregate. Using a simplified cluster model, we could estimate the size of individual clustersm 10 as well as the number of independent clusters N cl 500. Furthermore, we experimentally tested the validity of the effective cluster model. Using a many-body rate equation model, we find qualitative agreement between our theoretical model and the experimental results. We performed an extrapolation of the cluster size with the linear system size indicating that the clusters remain small. The absence of long-range correlations indicates the absence of a true phase transition to a bistable regime. Nevertheless, a bimodal excitation number distribution can be observed for small system sizes, as well as a dynamic hysteresis for small sweep times.
Our results show that the experimental and theoretical identification of phase transitions in open quantum systems require a careful analysis of the system size scaling and of the resulting correlation functions. To fully understand the thermodynamic limit of such systems, numerical simulations should be benchmarked with experimental data. Given the huge numerical effort to calculate and the experimental challenge to measure higher order correlation functions with good spatial and temporal resolution, this is in general a non-trivial task. We believe that temporal correlations are very helpful in this respect as they directly reflect the intrinsic dynamics and fluctuations of the steady state. Spatial correlations instead are typically limited to single destructive snapshots of the system and are insensitive to the dynamics.
In our system a transition to a truly bistable stationary phase is expected to be of first order [9]. An experimental realization would allow to study these types of phase transitions in stationary states of open systems. We have shown that anti-blockaded Rydberg gases with van-der Waals type interactions do not lead to diverging correlation times and long-range order in the experimentally relevant regime of large decoherence. Thus there is no phase transition to a bistable phase. Whether or not long range order can occur in a coherent regime and taking into account atomic motion remains an open question. Numerical simulations in this regime can be performed only for very small system sizes as the coherent many-body dynamics is no longer accessible by classical Monte-Carlo simulations. Here experiments are needed, which however require a substantially reduced decoherence level.

ACKNOWLEDGMENTS
We acknowledge financial support by the DFG within the SFB/TR185. We thank Carsten Lippe for discussing and revising the manuscript. We thank Dominik Linzner and David Petrosyan for stimulating discussions. F.L. and O.T. are recipients of a fellowship through the Excellence Initiative MAINZ (DFG/GSC 266). This work was supported by the Allianz für Hochleistungsrechnen Rheinland-Pfalz (AHRP). An experimental cycle consists of loading a 3D magneto optical trap of 87 Rb for 2.5 s from a 2D magneto optical trap. The precooled atoms are loaded into a crossed optical dipole trap, generated from a Nd:YAG fiber amplifier at 1064 nm, and undergo forced evaporation within 4 s, resulting in a Bose-Einstein condensate of 20 000 atoms in an isotropic trap with a trap frequency of 2π × 64 Hz. After evaporation we perform a 2 ms long linear scan with a scanning electron microscope across the cloud, yielding an ion signal proportional to the atom number while only marginally influencing the atomic sample. This allows us to check for atom number fluctuations during preparation and to correct for long term instabilities. Because of a finite magnetic field gradient present during evaporation, the sample is spin polarized in the F = 1, m F = +1 state from which we transfer it to the F = 2, m F = +2 state by a microwave Landau-Zener sweep. We load the prepared condensate into a 3D optical lattice, created from a Ti:sapphire laser at 748 nm, with lattice constants a x,y = 374 nm and a z = 529 nm in an exponential ramp with time constant τ = 20 ms to a final lattice depth of S = 20E rec . This ensures that we are in the Mott insulating phase regime and have a maximum of one atom per lattice site [50].
To measure the cluster dynamics of Rydberg aggregates we couple the atomic cloud with a one-photon transition at 297 nm to the Rydberg state 25P 1/2 with fixed detuning ∆ and intensity I for 100 ms. The laser light is produced from a frequency doubled dye laser at 594 nm. During the light-matter interaction, Rydberg excitations generated in the cloud are photoionized by the trapping light and black-body radiation with a rate of Γ ion 2 kHz and guided by a small electric field of 90 mV cm −1 to a discrete dynode detector. The underlying lattice structure ensures hereby that effects from associative ionization [51] or molecular Rydberg states [50] are negligible. Nevertheless, we only probe with blue detuning where molecular states are only addressable through a spin-flip in the second atom [52], which is strongly suppressed compared to a normal molecular excitation. We analyse the detected ion signal by calculating the g (2) -correlation function binned with 1 µs. We correct for the overall decay of the sample by normalizing with the averaged ion rate of multiple runs. Because of detector ringing and ion repulsion during the time of flight to the detector we can't analyse the first two microseconds of the correlation function.
The seed rate Γ seed and Rabi frequency Ω for each experimental parameter are determined by an analysis of the arrival time of the first detected ion (see Fig. 9). Note that the first ion analysis is independent from the cluster dynamics analysis. By doing a statistical analysis for each laser power and detuning over multiple experimental runs we can extract the multi-particle seed rate, taking into account our finite detection and ionization efficiency. Additionally, we extract the underlying Rabi frequency by fitting the seed rates for different detuning with Γ seed = 2Ω 2 γ ∆ 2 (for ∆ γ). We checked that the extracted Rabi frequencies approximately follow the expected √ I scaling.

Appendix B: Cluster Model
Here, we study a simple 1D lattice of Rydberg atoms in anti-blockade configuration. Specifically, we set the lattice constant a equal to the facilitation radius r fac = (C 6 /∆) 1/6 assuming that Rydberg atoms interact via a typical van der Waals potential V (r) = C 6 /r 6 . We compare a full numerical simulation based on a manybody rate equation model [42][43][44][45] to our simplified cluster model introduced in Sec. IV.
As pointed out before, the first seed excitation is produced with a rate N at Γ seed . Now, this seed excitation triggers further excitations with a rate proportional to  9. We retrieve the seed rate Γ seed by fitting an exponential decay exp(−Γ seed t) to the probability of detecting a first ion after the time of flight t to our ion detector (inset). From a ∼ ∆ 2 fit to the extracted seed rates we can calculate the Rabi frequency Ω, taking into account the finite detection and ionization efficiency, as well as the bare decoherence rate.
In the case of a 1D system, the first excitation enhances the excitation rate of its two neighbors, see Fig. 10a.
We introduce a geometric coordination number z counting the number of atoms with enhanced excitation rate (B1). In the 1D case we identify a geometric coordination number z = 2 for the transition between cluster size m = 1 to m = 2, see Fig. 10b. First, let us neglect the effect of long-range interactions beyond the lattice constant a. Then, a cluster grows and shrinks with rate zΓ fac . This is in analogy to a random walk along the cluster size axis m with spreading ∝ √ Γ fac t. However, due to spontaneous decay with rate Γ sp , the size of a cluster m is limited. The competition between drive and decay leads to a finite size. A decay event may split the cluster into two parts, see Fig. 10b, leading to a doubling of the coordination number z. However, this requires that the individual parts of the cluster are spatially well separated and do not influence each other. For instance, the last configuration in Fig. 10a shows that the two parts of the cluster can not merge into a single one. Without the effect of two atoms blocking each other, the coordination number z 2n + 1 increases along the axis of splittings n. In Fig. 10b, the transition rates between all possible configurations labeled by cluster size m and number of splittings n are shown.
In the following, we discuss assumptions, which allow us to estimate the size and the lifetime of a cluster. Since we are mostly interested in the cluster size m, we firstly neglect all splitting processes fixing the dynamics to the case n = 0. This is correct in the regime of weak driving Γ fac Γ sp with coordination number z = 2 due to the geometry. For strong driving Γ fac ≥ Γ sp , we have to account for an increase in the coordination number z.
To do so, we determine an effective coordination number z eff self-consistently. This results in the simplified cluster model discussed in Sec. IV with Γ ↑ = z eff Γ fac , see Fig.  1c. Now, using a detailed balance ansatz and projecting to the case of having a cluster (m ≥ 1), we calculate the cluster size distribution p m . This yields the iterative formula with normalization condition m=1 p m = 1. Using the full width half maximum valuem of the cluster size distribution p m , we can estimate the cluster lifetime τ cl to be equal with equation (14): Note that the cluster lifetime agrees with the lifetime of a single Rydberg excitation in the limit Γ fac Γ sp , i.e. m = 1. Now, we compare the dynamical properties of a large lattice system using a many-body rate equation method to the simplified cluster model. As discussed in Sec. III, to obtain the lifetime of a single cluster, we calculate the second order temporal correlation function g (2) (τ ). We extract the relaxation time, which we identify with the cluster lifetime, and the bunching parameter g (2) (0). The results for different detunings ∆ and fixed facilitation radius r fac = a are shown in Fig. 11. As we would expect, with increasing driving strength Ω/Γ sp the lifetime of the cluster increases. This allows us to determine the effective coordination number z eff used in the simplified cluster model. Besides the geometric coordination number z = 2, we account for the number of splittings, which occur with increasing cluster lifetime, by setting where τ sp = Γ −1 sp is the lifetime of a single Rydberg excitation. Comparing the effective cluster model (dashed line) with the full rate equation simulations (bullet points), we find excellent agreement for small driving strength Ω/Γ sp . However, the full simulations show that the cluster lifetime decreases again with Ω/Γ sp depending on the detuning ∆. We can understand this as an effect of interacting clusters in the high excitation density regime, which is beyond the simple cluster model. Increasing the detuning ∆, we decrease the seed rate Γ seed with which new clusters are produced. Therefore, the Rydberg excitation density ρ e decreases, see inset Fig. 11a. In our numerical simulation, we identify the regime where cluster collisions are important near an excitation density ρ e 0.2. These collisions strongly reduce the independent growth of a cluster.
Moreover, using the effective cluster model and including the first seed excitation in the cluster size distribution, we can estimate the bunching parameter L ×  (g (2) (0) − 1) = m(m − 1) / m 2 . Comparing both, the full simulations with the simplified cluster model, we find good agreement. With increasing driving strength Ω/Γ sp , the bunching parameter, which measures deviations from Poisson counting statistics, decreases. Therefore, with growing Ω/Γ sp , the number of independent clusters N cl and the lifetime τ cl increases. Both contribute to a reduction of the bunching parameter.

Appendix C: Rydberg-Rydberg interaction potential
To calculate the interaction potential between two Rydberg atoms we perform an exact diagonalization of the interaction Hamiltonian as shown in [46]. We use a basis of pair states |n 1 l 1 j 1 m j1 n 2 l 2 j 2 m j2 in the energetic vicinity of the state of interest and include all interaction terms up to quadrupole-quadrupole interactions as well as a small finite magnetic field of 1 Gs. For the calculations shown here the basis consists of 1604 pair states with a maximum angular momentum quantum number of l = 3 in an energetic vicinity of 50 GHz to the pair state |25P 1/2 , 25P 1/2 . The results of the diagonalization are shown in Fig. 12. For the |25P 1/2 , 25P 1/2 state there are four different interaction channels depending on the m j state of the pair state. The symmetric superposition 1 √ 2 (|+1/2, −1/2 + |−1/2, +1/2 ) shows the strongest attraction, the two fully stretched states are nearly degenerate but still attractive, while the antisymmetric superposition gets blue shifted by an dipole-quadrupole interaction with the |22F 5/2 24D 3/2 m j -manifold. We fit the relevant repulsive interaction potential with a C 9 /r 9 -type potential and use it in the simulations presented here. Interaction potentials of the pair state |25P 1/2 , 25P 1/2 obtained from a diagonalization of the interaction Hamiltonian including up to quadrupole-quadrupole interactions. The different potentials correspond to different superpositions of the mj-states. The green dashed line is a fit to the only repulsive potential with V (r) = C9/r 9 and C9 = 2π × 2.1 kHzµm 9 . FIG. 13. Excitation rate Γex/Γsp for three different inter particle distances r0 using the parameters from the experiment and Ω/2π = 500 kHz. The excitation rate for r0 = 2ax already agrees well with the uncorrelated excitation rate r0 → ∞. The inset shows the same rates with linear scale.

Appendix D: Inhomogeneous Broadening
Here, we discuss an inhomogeneous broadening mechanism resulting from the finite width σ of a localized wave packet and the steep slope of an interaction potential V (r) = C α /r α . The excitation rate Γ ex for a ground state atom in the presence of an atom already excited to the Rydberg state is strongly suppressed and broadened. As already discussed in [48,49], this stems from the coupling to motional states.
Firstly, using the rate equation approximation [42], the bare excitation rate of an atom at distance r 0 to an ex- FIG. 14. Comparison of the bare excitation rate Γ 0 ex /Γsp without additional broadening mechanism, the full excitation rate Γex/Γsp eq. (D3) and the approximate excitation rate Γ L ex /Γsp. Here, we use the parameters from the experiment and the inter particle distance is r0 = ax and Rabi frequency Ω/2π = 500 kHz. The inset shows the same rates with linear scale.
cited Rydberg state is However, this is only true in the limit of vanishing width σ → 0. Now, consider an atom localized in an optical lattice. Using the harmonic oscillator approximation, we can approximate the extent of the localized Wannier wavefunction Ψ(r) = 1 πσ 2 with the harmonic oscillator width σ. In a semi-classical approximation, we interpret |Ψ(r)| 2 as the probability distribution of finding an atom at position r. Then, the total excitation rate Γ ex for an atom localized with width σ is where r 0 is the average distance to the already excited Rydberg atom. Note, in general we have to evaluate here a 3D integral over the positions of all excited atoms. In this case the effective detuning ∆ − j Cα | rj | α describes the interaction shift landscape generated from all excited Rydberg atoms.
In Fig. 13 we plot the excitation rate Γ ex for three different experimental relevant distances r 0 . For two neighboring atoms in a lattice (r 0 = a x ), the excitation rate is strongly suppressed and broadened compared to the uncorrelated excitation rate Γ ex already occurring at r 0 = 2a x . However, for excitation rates Γ ex /Γ sp 1 larger than the lifetime of a single Rydberg excitation, we  can still expect cascaded excitations. Here, we identify a large region up to ∆/2π 40 MHz where an excitation cascade can be possible.
The exact numerical evaluation of the excitation rate (D3) in a large scale 3D system using a stochastic manybody rate equation model is numerically challenging. However, to include the strong impact of the inhomogeneous broadening, we use an approximate decoherence rate γ = γ 0 + |∂ r V (r = r 0 )|σ/ √ π, which can be easily extended to many atoms, see Sec.
In Fig. 14, we compare the full excitation rate (D3) with the excitation rate Γ L ex using the decoherence rate (D4) for r 0 = a x . Moreover, we plot the bare excitation rate Γ 0 ex without additional broadening. The approximate rate Γ L ex nicely describes the suppression of excitation and the broadening compared to the bare rate Γ 0 ex . In particular, the regime ∆/2π 18 MHz agrees well with the full model (D3). However, for smaller detunings ∆/2π 18 MHz the simple approximate model is not valid anymore. Therefore, we believe that the approximate decoherence rate is an efficient and good approximation at least for detunings ∆/2π 18 MHz.

Appendix E: Hysterese Relaxation
Here, we discuss a simple model for the relaxation of the hysteresis area. In the case of a single quench from a parameter ∆ 0 to ∆ 1 , we typically have an exponential relaxation of some observable n with relaxation timescale T 1 : n(t) =n 1 + (n 0 −n 1 )e −t/T1 .
Here,n 1 is the equilibrium value of the observed parameter and n 0 its initial value, which is not necessarily at equilibrium.
In the case of a hysteresis, we perform N small quenches ∆ j to ∆ j+1 on a timescale τ . The total sweep time is τ sweep = N τ . We denoten j the equilibrium value at ∆ j and n j = n(jτ ) the non equilibrium value after time jτ . The relaxation time T j = T (∆ j ) might depend on the tuning parameter ∆ j . Then, after N consecutive quenches of the form eq. (E1), the non equilibrium value of the observable n is given by Now, let us perform a continuum limit, where we keep the sweep rate∆ = ∆ N −∆0 τsweep constant and increase the number of steps N → ∞. We assume equidistant steps in the quench parameter δ = ∆ j+1 − ∆ j . Since the sweep time τ sweep and the parameter regime ∆ N − ∆ 0 are constant, we have δ, τ → 0. This allows us to approximate the exponential 1 − e −τ /Tj τ /T j (E3) and replace the sum with an integral. We obtain the continuum counterpart of eq. (E2): For a constant relaxation time T = T (∆), we obtain eq. (16). In Fig. 15 we show the relaxation of the Rydberg excitation number with increasing sweep time τ sweep . The result shows the expected power law relaxation τ −1 sweep .