Non-equilibrium fluctuations and metastability arising from non-additive interactions in dissipative multi-component Rydberg gases

We study the out-of-equilibrium dynamics of dissipative gases of atoms excited to two or more high-lying Rydberg states. This situation bears interesting similarities to classical binary (in general $p$-ary) mixtures of particles. The effective forces between the components are determined by the inter-level and intra-level interactions of Rydberg atoms. These systems permit to explore new parameter regimes which are physically inaccessible in a classical setting, for example one in which the mixtures exhibit non-additive interactions. In this situation the out-of-equilibrium evolution is characterized by the formation of metastable domains that reach partial equilibration long before the attainment of stationarity. In experimental settings with mesoscopic sizes, this collective behavior may in fact take the appearance of dynamic symmetry breaking.

While single-component systems, where one Rydberg transition is driven, have been the focus of many efforts, the dynamics of multi-component Rydberg gases -i.e. systems with atoms excited to several Rydberg states-remains largely unexplored. As recent experiments are starting to probe multiple Rydberg states [29][30][31][32], it is important to achieve some understanding of the collective phenomena that can be expected to be found in such systems. One can anticipate that several competing length scales will arise from the interplay between intra-level and inter-level interactions (i.e. the interactions between atoms excited to the same or different levels, respectively). Indeed, a few theoretical studies have started exploring this competition [33,34].
The study of multi-component systems may help to further the strong analogies between the dynamics of dissipative Rydberg gases and soft-matter systems [1,12,15]. This connection ultimately originates from the Rydberg blockade effect [35,36], whereby an excitation of a given atom prevents that of neighboring atoms, an effect reminiscent of the excluded-volume interactions characteristic of soft-matter systems such as liquids and colloids [37]. These systems are often mixtures composed of more than one kind of particle, as such dispersity can give rise to interesting collective effects that are not present in the monodisperse case, see e.g. [38]. Furthermore, it is common when modelling soft matter computationally to consider "non-additive" mixtures, meaning mixtures where the cross interactions between different kinds of particles are not given by those between similar kinds: for example, if particles A and B interact among themselves with typical distances σ A and σ B , respectively, the distance for cross interaction is such that σ AB = (σ A + σ B )/2, as in e.g. Ref. [39]. An illustration is given in Fig. 1 (a). While non-additive interactions are unphysical in a classical setting where particles interact by excluded volume or similar effects, they are used to increase frustration in model liquids, thus precluding crystallization and promoting glass formation. In dissipative Rydberg gases the non-additivity of inter-atomic interactions is an experimentally realizable physical feature. Despite its quantum origin, such non-additivity survives in an effectively classical limit, giving a new handle on experimentally realizable binary (or generally p-ary) mixtures.
In this work, we elucidate the physics of multi-component dissipative Rydberg gases far from equilibrium. We first develop a general theory for the dynamics of systems with any number of components extending an approach that has been extensively validated in the one-component case [24,26]. We then perform an idealized numerical study, where different interaction strengths lead to a variety of length scales giving rise to strikingly different dynamical regimes. The phenomenology that emerges from non-additive interactions is characterized by the formation of domains, which are homogeneously populated by excitations of a given component when inter-level interactions dominate, and show an alternation of components in the opposite case. Homogeneous domains reach partial equilibration when detailed balance is achieved for the dominant atomic transition, leading to metastable behavior. In experimental settings with mesoscopic sizes, these domains will appear as non-equilibrium symmetry-broken states.

THEORY: EFFECTIVE DYNAMICS IN THE LIMIT OF STRONG DISSIPATION
We consider a system of N atoms, each of which can be in one of p + 1 levels, the ground state |0 , and p > 1 Rydberg states |1 , |2 , . . . |p , with energies E 0 < E 1 < E 2 < · · · < E p . See Fig. 1 (b) for an ilustration of the two-component case. Atoms in the Rydberg states |s and |s at positions r k and r m interact through a power-law potential V ss km = C ss α /|r k − r m | α with exponent α. For simplicity, we denote the intra-level interactions by V s km instead of V ss km . The value of the coefficients C ss α depends on the specific structure of the atomic spectrum and can be controlled through e.g. electric field induced Förster resonances [40] or microwave dressing [41,42]. Typically encountered exponents are α = 6 (van der Waals interaction) and α = 3 (dipole-dipole interaction) [43]. Each of the Rydberg states is resonantly coupled to the ground state by a laser field, and affected by dephasing noise [24][25][26]. The dynamics of the system is governed by a Master equation of Lindblad form ∂ t ρ = Lρ + D(ρ) [44]. The coherent part

ρ] includes an interaction Hamiltonian
and a driving term and Ω s is the Rabi frequency of the transition between |0 and |s . The dissipator is given by s , ρ , where γ s is the dephasing rate of |s w.r.t. |0 . Atomic decay is not considered, as we are especially interested in the short time dynamics that has been probed in experiments [24,26]. We deliberately focus on a situation where exchange interactions can be omitted, which can be achieved by a specific choice of Rydberg states [40].
In the limit of strong dissipation, Ω s γ s , which is relevant in a number of experimental settings [24,26], the time evolution is governed by an stochastic dynamics along the classical states represented in µ = diag(ρ) [1,45]. While the effective equations of motion of the multi-component Rydberg gas are crucially important for the rest of the paper, and simple enough as to provide physical insight into the phenomenology that is numerically observed (which would be very hard to infer from the quantum master equation), their derivation is relatively lengthy. We therefore include here only the main results, and give the technical details in Appendix A for the interested reader. The resulting rate equations are where I s + |0 k 0| projects on the subspace spanned by the ground state and the excited state |s of site k. For simplicity, we assume that the atoms sit in the sites of a chain with lattice constant a. A transition involving the excited level |s at site k, whether it is an excitation or a de-excitation, occurs with a rate wherer k = r k /a. The relevant length scales are given by the intra-level, R s = a −1 [2C s α /γ s ] 1/α , and the interlevel interaction parameters, R ss s = a −1 [2C ss α /γ s ] 1/α , which are the reduced distances at which the appearance of excitations of a given component correlate different sites. As in classical mixtures of particles (liquids, colloids, etc.) several components coexist and their interactions are characterized by different typical length scales depending on the components involved.
Experiments typically probe the dynamics starting from an initial state where all atoms are in the ground state, and this will be our choice as well. At the initial stages distant excitations to any level occur independently of each other with a rate that is O(1). This gives an "initial seed" for the correlated dynamics: as soon as the distance between excitations becomes comparable with R s and/or R ss s , the second term in Eq. (3) strongly correlates the atoms, and the transitions between the ground state and a particular level become less likely due to the presence of nearby excited particles [see Fig. 1 (c)].

PHENOMENOLOGY: NUMERICAL RESULTS
We next turn to a numerical exploration of the phenomenology that emerges in multi-component Rydberg gases.
As the dynamics is rich in collective effects, we start from the simplest possible case of a two-component system with symmetric interaction parameters, R 1 = R 2 ≡ R and R 12 1 = R 12 2 ≡ R c . The expression in brackets on the rhs of Eq. (3) then becomes m R α n (m) 1 /|r k −r m | α for transitions between |0 k and |1 k , and an equivalent expression for the transition between |0 k and |2 k is obtained by swapping n In keeping with the aim to simplify the parameter space as much as possible, we further assume Ω 2 1 /γ 1 = Ω 2 2 /γ 2 , and rescale the time variable by Ω 2 1 /γ 1 . We focus our study on three generic cases: (i ) R > R c , (ii ) R ∼ R c and (iii ) R < R c . Cases (i ) and (iii ) are examples of non-additive interactions, R = R c . Such interactions are of interest in the theoretical study of complex and glassy dynamics in classical mixtures, but their experimental realization remains challenging in those contexts, whereas they appear generically in Rydberg gases. We expect that in case (i ) the excitations of one component will be surrounded by excitations of the other component, in an anticorrelated pattern, as in the upper panel of Fig. 1 (c). By analogy, in case (iii ), one expects the clustering of excitations of a given component, as in the lower panel of Fig. 1

(c).
This phenomenology is indeed observed using kinetic Monte-Carlo simulations in a 1D chain of van der Waalsinteracting atoms (α = 6). We focus on a mesoscopic system of size N = 20, as such sizes are accessible by current experiments, and use periodic boundary conditions to prevent uncontrolled boundary effects. In Fig. 2 we show representative trajectories for cases (i ), (ii ) and (iii ) for fixed R = 2 and varying R c , where an atom appears in blue if it is in the excited state |1 , in red if it has been excited to |2 , and in white if it is in the ground state. Analogous results for four components are presented in Appendix B. For R c = 1/2 (R > R c ), the excitation pattern forms something that can be described as heterogeneous domains of alternating excitations of one and the other component, with some defects [ Fig. 2  that shown in Fig. 2 (d), domains of a given component dominate throughout the non-equilibrium regime. As for the situation in which R c = 2 [R = R c , Fig. 2 (b)], corresponding to additive interactions, excited atoms are as likely to be found close to excitations of either component throughout the evolution of the system. Indeed, the components act as labels that permit to distinguish different types of excitations, but they have no dynamical consequences. This is in stark contrast to situations in which the interactions are non-additive, where (as shown above) the configurations that emerge are highly dependent on the components of the excitations. As in mixtures of classical particles, non-additivity brings richness into the dynamics.
With increasing time, the lattice fills with more and more excitations. Eventually these highly structured configurations disappear and the system settles into the stationary state of the master equation, which is proportional to the identity, ρ st ≡ (p + 1) −N ⊗ k I k . Accordingly, in the effective dynamics the average number of atoms in each level becomes N/(p + 1), as can be seen from Eq. (2). In Fig. 3 (a) we show n 1 (t) ≡ (1/N ) k n (k) 1 (t) , i.e. the density of atoms in the excited state |1 , as a function of time. This observable gives us some important information of the generic aspects of the classes of dynamics illustrated in Fig. 2 for specific realizations. We again fix R = 2, and look at R c = 1/2, 2 and 8. The excitation density for the situations corresponding to the two extreme values of R c increases until it reaches a long plateau which has been highlighted with vertical arrows in the figure. The origin of these plateaus will be clarified below. Much later another increase leads the system towards the sationary state (see the black horizontal line). While the results for n (k) 2 (t) are identical, single trajectories fluctuate strongly, a situation reminiscent of dynamic symmetry breaking (see, e.g., [46,47]).
To gain insight into the relaxation behavior reported in Fig. 3 (a), and especially on the type of configurations that occur at different stages of the dynamics, it is useful to complement the study of the time evolution of the density of excitations with that of an observable that can help distinguish between different local patterns of excitations. For this purpose, we focus on the density of excitation blocks, i.e. excited atoms whose right neighbors are also excited. We consider separately homogeneous and heterogeneous blocks, which are made up of same-component or different-component excitations respectively. The former type is shown enclosed in a continuous-line box and the latter in a dashed-line box in the inset of Fig. 3 (b), where we show the density of homogeneous (continuous line) and heterogeneous blocks (dotted line). The block density in the stationary state is indicated by a black horizontal line. This is n /N = 2/(p + 1) 2 for homogeneous blocks, and the same value can be easily seen to apply to heterogenous blocks. For R c = 1/2 the density reaches the plateau in Fig. 3 (a) at the time the concentration of heterogeneous blocks gets close to the equilibrium value, and the final push into stationarity corresponds to an equivalent move on the part of the concentration of homogeneous blocks. This corresponds to a rapidly achieved alternating pattern of excitations, which persists for long times until it finally relaxes into the stationary state. For R c = 8, we see the opposite behavior: the plateau is reached first when the homogeneous blocks attain the equilibrium value, and stationarity is achieved after a long wait when the heterogeneous blocks reach that value as well. The interpretration is analogous to that of the R c = 1/2 case, but now the domains that appear at the time the plateau is reached are homogeneous. In which case it takes shorter or longer for the homogeneous or the heterogeneous blocks to reach the equilibrium value can of course be inferred from the rates in Eq. (3). The case where R c = R = 2 unsurprisingly shows a simultaneous equilibration of both types of blocks, and therefore stationarity is reached without an intermediate plateau.
These results suggest that the domain structure remains in place for very long times before reaching stationarity. To clarify this we consider the order parameter for the study of homogeneous domains. Bimodal distributions of this parameter indicate very strong dominance of one of the two components, while a narrow unimodal distribution that peaks at zero indicates the existence of configurations where both components are strongly mixed. We further define a similar parameter that assigns an alternating sign to consecutive excitations along the chain for the study of heterogeneous domains The probability distribution of P ± across time for R = 2, R c = 1/2, 2 and 8 is shown in Fig. 4 (a), (b) and (c), respectively. For R c = 1/2 [ Fig. 4 (a)], P − (t) has a relatively wide distribution that narrows down as the system approaches stationarity, indicating the loss of order. The presence of defects makes the distribution unimodal even for short times. For R c = 2 [ Fig. 4 (b)], P + (t) is narrowly distributed around zero, as the occurrence of both components is equally likely in all realizations. This case can be analyzed with P − (t) as well, yielding very similar results (not shown). A richer phenomenology occurs when R c > R [R c = 8, Fig. 4 (c)], with a clearly bimodal distribution throughout the non-equilibrium evolution of the system. Initially, two branches are formed symmetrically around zero, separated by a region of very low probability of occurrence. This corroborates the role of the initial seed in leading the system to domains of either component. The two peaks of P + (t) separate more and more until they saturate. Later on, the domain structure starts crumbling upon the appearance of excitations of the non-dominant component. This is illustrated in Fig. 4 (d), where the curves corresponding the distributions shown in (c) at times t = 10 2 , 10 4 and 10 8 are shown.
The saturation value for R c > R is |P + | 0.5, which corresponds to a metastable state, as we now explain. In the two-component case, the right-hand side of Eq. (2) contains two terms for each site. Within a homogeneous The corresponding term reaches a "partial equilibrium" when σ 1 µ is as likely to create excitations as de-excitations, i.e. when there are as many atoms in the ground state as in state |1 . Indeed, this state, in which detailed balance is satisfied for one of the transitions between the ground state and an excited state, would correspond to the stationary dynamics of a single-component Rydberg gas. In multi-component systems, however, excitations of the non-dominant component have to appear eventually in order for the system to reach the true stationary state, as shown in Fig. 4 (c) and (d). A similar behavior is observed in four-component systems (see Appendix B). The reader should note that such metastable states may not be achieved in experiments starting from an empty initial state if the times required to reach them exceed the lifetimes of the atoms. Starting from densely populated initial states can be helpful in probing this metastability.

CONCLUSIONS
We have derived an effective theory for a multi-component Rydberg gas in the presence of noise. For non-additive interactions, the emerging dynamics displays a domain structure that depends sensitively on the initial excitations. For large inter-species interactions this leads to a metastable dynamics when partial equilibration is reached for the dominant component, which corresponds to the stationary state of a single-component system where only that Rydberg transition is driven. To our knowledge, this could be the first system that is experimentally accessible in which non-additive interactions of the kind that are considered in classical mixtures of particles for the study of metastable dynamics can be naturally implemented, and are indeed expected to occur generically. Whether the phenomenology persists at the qualitative level when the dissipation is only moderately strong or even weak compared to the driving, as has been recently shown to occur in the case of single-component Rydberg gases [19], is an interesting question that remains to be studied, as is the general role of quantum fluctuations [48]. The possibility that the (to some extent) tunable exchange interaction of Rydberg gases [40] can open up new relaxation pathways in multi-component systems will be explored in the future.

APPENDIX A. DERIVATION OF THE EFFECTIVE EQUATIONS OF MOTION
We consider a gas of N atoms in a lattice. The ground state |0 of each atom is resonantly coupled by laser fields to the Rydberg states |1 , |2 , . . . , |p (with energies such that E 0 < E 1 < · · · < E p ). The Master equation is then given by ∂ t ρ = L 0 ρ + L 1 ρ, where L 0 contains the interaction Hamiltonian and the dissipator, and L 1 gives the time evolution due to the driving. For the derivation below, where L 1 will be treated as a perturbation, this is more convenient than the more physical decomposition into a coherent part and a dissipator that is used in the main text. The Liouvillian superoperator L 0 is defined as where n (k) s = |s k s| and γ s is the dephasing rate of |s w.r.t. |0 . Atoms in the Rydberg states |s and |s at positions r k and r m interact through a power-law potential V ss km = C ss α /|r k − r m | α with exponent α. For simplicity we denote the intra-level interactions by V s km instead of V ss km . As a result, the Hamiltonian H 0 can be written as The superoperator L 0 therefore consists of a Hamiltonian part and a dissipator whose individual terms commute. Additionally, we have the driving term, which in the rotating-wave approximation becomes where σ (k) sx = |s k 0| + |0 k s| and Ω s is the Rabi frequency of the coupling between |s and |0 . Our aim is to derive the effective dynamics in the limit of strong dissipation, Ω s γ s . We start by working out the effect of the dissipator on the dynamics. Using the notation, The first-order contribution of the action of e L k 0,d t = I + L k 0,d t + 1/2! (L k 0,d ) 2 t 2 + O(t 3 ) on the density operator is where the dissipative evolution of site k has been made explicit using the basis states |0 k , |1 k , |2 k , . . . , |p k , and ρ (k) mn are the p N −1 × p N −1 matrices defined by ρ (k) mn = k n|ρ |m k . By analogously deriving higher order terms, it can be shown that the action of the dissipator e L0t ρ is The off-diagonal entries of the density matrix are seen to decay exponentially, a fact that is not altered by the action of the coherenct dynamics given by H 0 , which is diagonal in the product basis formed by single particle states |0 k ,|1 k , |2 k , . . . , |p k . Therefore, the evolution under L 0 becomes, at time scales much larger than the inverse of the dephasing rates, a projector P on the diagonal of ρ in that same basis as happens in the case of just one Rydberg level [1]. The removal of all coherences leads to a diagonal density matrix, where each classically accessible configuration (e.g. |00100203 · · · 1 ) is given a certain probability of occurrence.
Using the projector operator P and its complement Q = 1 − P, we can formulate the effective evoluton equation for the diagonal density matrix µ = Pρ describing the slow evolution. To second order in L 1 , the general expression is given by In this case PL 1 µ = 0 and Qe L0t QL 1 P = e L0t L 1 P. We next calculate the integrand in Eq. (12), We have used the fact that e L0t does not shift matrix elements, and that the action of σ (k) sx = |s k 0| + |0 k s| followed by that of σ In the following, we explicitly work out the terms in Eq. (13). We focus on the contribution corresponding to level |1 for concreteness.
We will use V k s = m V s km n (m) s + s =s V ss km n (m) s as shorthand to refer to the increment in the interaction energy that one has to pay for the excitation of atom k to level |s . In the expressions above V k 1 appears in the oscillatory part of the diagonal elements of the matrix.
The term corresponding to s = 1 in Eq. (13) is therefore and the contributions due to the other levels take an analogous form. Thus, Eq. (12) can be rewritten as where the projection operator I s + |0 0| cancels all the elements in µ that do not correspond to the ground state or |s at site k. The effective dynamics is therefore given by with rates for a transition |0 → |s or |s → |0 To make explicit the power-law interactions, it is useful to refer to the atomic spatial arrangement in terms of reduced position vectorsr k = r k /a, where a is the lattice constant. We define an intra-level interaction parameter R s = a −1 [2C s α /γ s ] 1/α (for interactions between atoms in the same level, V s km = C s α n In some cases, the interaction exponent α could be different depending on the atomic levels involved. This more general case can be easily worked out from Eq. (21), but here we will assume that α is the same for all level pairs.

APPENDIX B. PHENOMENOLOGY OF A FOUR-COMPONENT DISSIPATIVE RYDBERG GAS
While the derivation of the effective equations of motion is valid for any number of species, in the numerical results reported in the main text we focus on the two-component case, p = 2. However, both the main observations on the phenomenology and the theoretical arguments given there can be extended without great difficulty to the p > 2 case. In this section we briefly report some results for p = 4. For the sake of simplicity, we again use a somewhat idealized parameter choice according to which all the intra-level interaction parameters, which we collectively denote as R, are equal to one another, while all the inter-level interaction parameters, denoted as R c , are also equal among themselves. We focus on the R c > R case, where homogeneous domains emerge, as it gives the richest phenomenology. More specifically, we consider R = 2 and R c = 8, which coincides with the parameter choice used in the main text.
In Fig. 5 (a), we see one representative trajectory of a system of N = 20 atoms with van der Waals interactions. The color coding is such that red corresponds to state |1 , green to |2 , cyan to |3 , magenta to |4 and white to ground state atoms. While there are some initial excitations to |1 , they finally de-excite and are replaced by excitations to |4 , which by then has become the dominant component. The large homogeneous |4 -domain that emerges is later replaced by a |3 -domain. At longer times, two domains, corresponding to |1 and |4 coexist. Eventually, when stationarity is approached, the system undergoes a strong mixing of all the components. To quantify the emerging dynamical order we focus on a complex order parameter that is an extension of the real order parameter P + that was proposed in the main text for p = 2 [see Eq. (4)]. It is defined as follows This order parameter, which has been inspired by the theory of the Potts model [49], can be easily extended to any number of components p. In Fig. 5 (b) we show P 4+ (t) at t = 10 2 , which corresponds to the time at which most of the trajectories inspected still show one domain that spans the whole chain. The existence of as many maxima as there are excited levels, all of them quite distant from the origin, indeed indicates that the formation of large domains of the kind seen in Fig. 2 (c) and (d) of the main text for p = 2 occurs generically in systems with a larger number of components as well. As in the two-component case [main text, Fig. 4 (d)], the four peaks reach the saturation value of |P 4+ | 0.5 at later times, and eventually subside into a unimodal distribution centered around the origin when the system approaches the stationary state.