Many-body dynamics of holes in a driven, dissipative spin chain of Rydberg superatoms

Strongly interacting many-body systems subject to external driving and coupled to (possibly tailored) reservoirs offer a new route to create and stabilize interesting states of matter. As a simple example, a quantum state can be made immune to particle losses if it is the stationary state of an open system coupled to a particle reservoir. Furthermore, the competition between coherent driving and dissipation can lead to exotic steady states [1, 2] and phase transitions in open many-body systems [3–8].

Rydberg atoms [9] are well suited to study the interplay between strong interaction and coupling to coherent laser fields and dissipative environments. They are thus prime candidates to investigate many-body physics of driven, dissipative spin models. A prominent and well studied consequence of the strong, long-range interaction between atoms in Rydberg states is the so-called blockade phenomenon, whereby a Rydberg excited atom suppresses further excitations within a certain blockade distance [10]. Rydberg blockade in a dilute gas or in a regular array of single atoms leads to short-range spatial ordering of excitations, as was studied theoretically [11–14] and demonstrated experimentally [15, 16]. In the so-called anti-blockade regime, successive excitation of atoms at a certain distance from each other can be resonantly enhanced [17–23], which led to the lively debate on the possibility of attaining bistable steady states [22–26].

When many atoms are confined within the blockade distance from each other, they form an effective two-level system—Rydberg superatom—that can accommodate at most one Rydberg excitation [10, 13, 27]. The coupling of a superatom to the laser radiation is collectively enhanced, while the steady-state probability of a single Rydberg excitation can approach unity. This permits the level of control of single collective spins represented by superatoms far exceeding that for individual atoms. Moreover, being composed of many atoms, superatoms are relatively insensitive to atom number fluctuations and losses. Regular arrays of spins represented by superatoms can then be prepared with less experimental effort, which should be contrasted with the sophisticated dynamical preparation techniques used to realize defect-free arrays of individual Rydberg atoms [28, 29]. This, together with the strong, long-range interactions between the Rydberg excitations, makes superatoms ideal building blocks for realizing dissipative many-body spin models and analyzing their dynamics.

Single Rydberg superatoms have been observed in several experiments [16, 30–34]. A two-dimensional square lattice of superatoms with nearest-neighbor excitation blockade was studied in [35], demonstrating the possibility of phase transition to an anti-ferromagnetic steady state with spontaneously broken lattice symmetry. In the complementary interaction regime of the Rydberg anti-blockade, little is known about the many-body dynamics and the steady state of a lattice of superatoms. Here we study a one-dimensional lattice of Rydberg superatoms (see figure 1(a)), in which an already excited superatom facilitates resonant excitation of its neighbor, but the presence of two excited neighbors suppresses the excitation. This systems exhibits interesting excitation dynamics and a highly non-trivial steady state characterized by an almost universal density ${\rho }_{{\rm{R}}}=2/3$ of Rydberg excitations with strongly suppressed number fluctuations. We show that this behavior can be explained in terms of mobile excitation holes on the background of Rydberg excited lattice. The holes behave as a nearly incompressible liquid of hard rods with characteristic two-particle correlations (see figure 1(b)). We derive and verify an effective many-body model for holes. Varying the parameters of the effective model, we find a crossover between a liquid of holes with density–density correlations peaked at the distance of two lattice periods, $2a$, and the onset of crystalline order with period $3a$. In both cases the density of holes is ${\rho }_{h}=1/3$ with highly suppressed number fluctuations.

Zoom In Zoom Out Reset image size

The paper is organized as follows. In section 2, we formulate the model for a regular array of superatoms and derive the formalism for the efficient treatment of the system. In section 3 we present the results of numerical simulations of the dynamics of the chain of superatoms (driven spin chain) and introduce the effective hole model that leads to an intuitive physical picture for the equilibrium phase of the system. Section 4 summarizes our results. Technical derivations are deferred to the appendices A–C.

In this section, and in appendices A and B, starting from the fully quantum many-body master equation, we derive rate equations for the chain of laser-driven and mutually interacting Rydberg superatoms. These rate equations will then be used in section 3 for numerical simulations of the dynamics and steady state of the many-body system.

2.1. The microscopic model

We consider an ensemble of cold atoms trapped in a regular array of microtraps [32, 33, 36–38] or a long-wavelength optical lattice with the period a of several microns. Each lattice site j contains on average N atoms, see figure 1(a). A laser field of carrier frequency ω drives the atoms on the transition from the ground state $| g\rangle$ to the excited Rydberg state $| e\rangle$ with the Rabi frequency Ω and detuning ${\rm{\Delta }}=\omega -{\omega }_{{eg}}$. In the frame rotating with frequency ω, the coherent excitation dynamics of the atoms is described by the Hamiltonian (${\hslash }=1$)

$$\begin{eqnarray}&&{ \mathcal H }=\displaystyle \sum _{k}[{\rm{\Omega }}({\hat{\sigma }}_{{eg}}^{k}+{\hat{\sigma }}_{{ge}}^{k})-{\rm{\Delta }}{\hat{\sigma }}_{{ee}}^{k}]+\displaystyle \sum _{k\ne k^{\prime} }V({\vec{r}}_{k},{\vec{r}}_{k^{\prime} })\,{\hat{\sigma }}_{{ee}}^{k}\otimes {\hat{\sigma }}_{{ee}}^{k^{\prime} },\quad \end{eqnarray} \tag{ 1 }$$

where ${\hat{\sigma }}_{\mu \nu }^{k}\equiv | \mu {\rangle }_{k}\langle \nu |$ is the transition ($\mu \ne \nu$) or projection ($\mu =\nu$) operator for the kth atom, and $V({\vec{r}}_{k},{\vec{r}}_{k^{\prime} })=\tfrac{{C}_{p}}{| {\vec{r}}_{k}-{\vec{r}}_{k^{\prime} }{| }^{p}}$ is the interaction potential between pairs of atoms at positions ${\vec{r}}_{k},{\vec{r}}_{k^{\prime} }$ excited to the Rydberg state $| e\rangle$. The usual van der Waals interaction corresponds to p = 6, while the static dipole–dipole interaction to p = 3 [9].

Atoms excited to the Rydberg state $| e\rangle$ spontaneously decay to the ground state with the rate ${{\rm{\Gamma }}}_{s}$, and are dephased with a typically much larger rate ${{\rm{\Gamma }}}_{d}$ due to atomic collisions and motion in the inhomogeneous trapping potential, collisions of the Rydberg electron with the ground state atoms [39–41], or intermediate state decay if $| g\rangle \to | e\rangle$ is a two-photon transition [15, 42]. The dissipative dynamics is described by the master equation for the density matrix $\hat{\rho }$ of the system,

$$\begin{eqnarray}&&{\partial }_{t}\hat{\rho }=-{\rm{i}}[{ \mathcal H },\hat{\rho }]+\displaystyle \sum _{k}\displaystyle \sum _{\alpha =s,d}\left({L}_{\alpha }^{k}\hat{\rho }{L}_{\alpha }^{k\dagger }-\tfrac{1}{2}\{{L}_{\alpha }^{k\dagger }{L}_{\alpha }^{k},\hat{\rho }\}\right),\end{eqnarray} \tag{ 2 }$$

where ${L}_{\alpha }^{k}$ are the Lindblad jump operators for the two relaxation processes assumed acting independently on each atom k as ${L}_{s}^{k}=\sqrt{{{\rm{\Gamma }}}_{s}}{\hat{\sigma }}_{{ge}}^{k}$ and ${L}_{d}^{k}=\sqrt{{{\rm{\Gamma }}}_{d}}{\hat{\sigma }}_{{ee}}^{k}$.

2.2. Rate equations for Rydberg superatoms

In equation (1) we can split the sum over all the atoms into two parts: the sum over the lattice sites j, and the sum over the atoms k_j in each lattice site. We assume that all N atoms within each lattice site are confined within a small spatial interval $\delta r\ll a$, such that the interatomic interaction energy ${C}_{p}/{(\delta r)}^{p}$ exceeds all the relevant energy scales pertaining to the atoms, namely, the laser Rabi frequency Ω and detuning Δ, as well as the atomic spontaneous decay ${{\rm{\Gamma }}}_{s}$ and dephasing ${{\rm{\Gamma }}}_{d}$ rates and the resulting Rydberg state excitation linewidth $w\simeq 2{\rm{\Omega }}\sqrt{\gamma /{{\rm{\Gamma }}}_{s}}$ with $\gamma \equiv \tfrac{1}{2}({{\rm{\Gamma }}}_{d}+{{\rm{\Gamma }}}_{s})$ [13]. This allows us to neglect all the multi-atom states containing more than one Rydberg excitation per lattice site [10, 13]. If on the spatial scale $\delta r\lesssim 1\ \mu {\rm{m}}$ of such a Rydberg superatom the laser field can be assumed uniform, it would couple the collective ground state $| G\rangle \equiv | {g}_{1}{g}_{2}\ldots {g}_{N}\rangle$ only to the symmetric single excitation state $| {E}_{s}\rangle \equiv \tfrac{1}{\sqrt{N}}{\sum }_{k}| {g}_{1}{g}_{2}\ldots {e}_{k}\ldots {g}_{N}\rangle$ with the collectively enhanced Rabi frequency $\sqrt{N}{\rm{\Omega }}$, see figure 2(a). There are, in addition, $(N-1)$ non-symmetric single excitation states $| {E}_{{ns}}{\rangle }_{m}$, labeled by index m, decoupled from the laser field. We can then recast the Hamiltonian as

$$\begin{eqnarray}{ \mathcal H } & \simeq & \displaystyle \sum _{j}[\sqrt{N}{\rm{\Omega }}(| {G}^{j}\rangle \langle {E}_{s}^{j}| +| {E}_{s}^{j}\rangle \langle {G}^{j}| )]-{\rm{\Delta }}{\hat{P}}_{{EE}}^{j}\\ & & +\displaystyle \sum _{j\gt i}\displaystyle \frac{{C}_{p}}{| {\vec{r}}_{j}-{\vec{r}}_{i}{| }^{p}}{\hat{P}}_{{EE}}^{j}\otimes {\hat{P}}_{{EE}}^{i},\end{eqnarray} \tag{ 3 }$$

where ${\hat{P}}_{{EE}}^{j}\equiv | {E}_{s}^{j}\rangle \langle {E}_{s}^{j}| +{\sum }_{m}| {E}_{{ns}}^{j}{\rangle }_{m}\langle {E}_{{ns}}^{j}|$ is the projector onto the manifold of N single Rydberg excitation states of superatom at site j, and ${\vec{r}}_{j}=\tfrac{1}{N}{\sum }_{{k}_{j}}{\vec{r}}_{{k}_{j}}$ is its center of mass coordinate.

Zoom In Zoom Out Reset image size

Within each Rydberg superatom, the dephasing couples incoherently the symmetric single excitation state $| {E}_{s}\rangle$ to all $(N-1)$ non-symmetric states $| {E}_{{ns}}\rangle$ with the rate ${{\rm{\Gamma }}}_{d}(N-1)/N$ which approaches ${{\rm{\Gamma }}}_{d}$ for $N\gg 1$. In turn, the reverse coupling of the non-symmetric states to the symmetric state has the rate ${{\rm{\Gamma }}}_{d}/N$. All single excitation states decay back to the ground state $| G\rangle$ with rate ${{\rm{\Gamma }}}_{s}$, see figure 2(a). In appendix A we derive a simple rate equation model describing the excitation dynamics of a superatom in the limit of strong dephasing ${{\rm{\Gamma }}}_{d}\gtrsim {\rm{\Omega }}$. Starting from the density matrix equations for a single Rydberg superatom, we adiabatically eliminate all coherences and the population of the symmetric excited state, which scales as $\sim 1/N$. The superatom then reduces to an effective two-level system, see figure 2(b), and its dynamics is governed by the rate equations for the populations of the ground ${\rho }_{{GG}}^{j}$ and excited ${\rho }_{{EE}}^{j}$ states,

$$\begin{eqnarray}\displaystyle \frac{\partial }{\partial t}{\rho }_{{EE}}^{j} & = & {{\rm{\Gamma }}}_{\mathrm{ex}}({{\rm{\Delta }}}_{j}){\rho }_{{GG}}^{j}-{{\rm{\Gamma }}}_{\mathrm{de}}({{\rm{\Delta }}}_{j}){\rho }_{{EE}}^{j},\\ & & {\rho }_{{GG}}^{j}+{\rho }_{{EE}}^{j}\simeq 1,\end{eqnarray} \tag{ 4 }$$

where the excitation and de-excitation rates are given by

$$\begin{eqnarray}&&{{\rm{\Gamma }}}_{\mathrm{ex}}({{\rm{\Delta }}}_{j})=\displaystyle \frac{(N-1)\chi ({{\rm{\Delta }}}_{j})}{N\,\chi ({{\rm{\Delta }}}_{j})+2\gamma }{{\rm{\Gamma }}}_{d},\end{eqnarray} \tag{ 5a }$$

$$\begin{eqnarray}&&{{\rm{\Gamma }}}_{\mathrm{de}}({{\rm{\Delta }}}_{j})=\displaystyle \frac{\chi ({{\rm{\Delta }}}_{j})}{N\,\chi ({{\rm{\Delta }}}_{j})+2\gamma }{{\rm{\Gamma }}}_{d}+{{\rm{\Gamma }}}_{s},\end{eqnarray} \tag{ 5b }$$

with $\chi ({{\rm{\Delta }}}_{j})\equiv \tfrac{2{{\rm{\Omega }}}^{2}\gamma }{{\gamma }^{2}+{({{\rm{\Delta }}}_{j})}^{2}}$. Here

$$\begin{eqnarray}&&{{\rm{\Delta }}}_{j}={\rm{\Delta }}-\displaystyle \sum _{i\ne j}\displaystyle \frac{{C}_{p}}{{| {\vec{r}}_{i}-{\vec{r}}_{j}| }^{p}}{\rho }_{{EE}}^{i}\end{eqnarray} \tag{ 6 }$$

is the effective detuning of superatom j which includes the Rydberg level shift due to the interaction with all the other superatoms in the Rydberg state. In figure 2(c), and in more detail in appendix A, we compare the dynamics of a single superatom as obtained from the exact solution of the complete set of the density matrix equations and the solution of the rate equations. We observe that the rate equation model approximates well the relaxation timescale of a superatom and the steady-state population of the excited state, ${\rho }_{{EE}}=\tfrac{{{\rm{\Gamma }}}_{\mathrm{ex}}}{{{\rm{\Gamma }}}_{\mathrm{ex}}+{{\rm{\Gamma }}}_{\mathrm{de}}}$. Remarkably, unlike for a single two-level atom, the excitation probability of the superatom under continuous (near-)resonant driving and in the presence of strong dephasing $\gamma \gtrsim {\rm{\Omega }}$ can approach unity, ${\rho }_{{EE}}\simeq \tfrac{N}{N+1}\to 1$, with increasing N.

2.3. Facilitated excitation of superatoms

The laser irradiates continuously and uniformly the 1D chain of superatoms. We set the detuning Δ of the laser to be equal to the interaction strength $V(a)={C}_{p}/{a}^{p}\equiv {V}_{{\rm{N}}}$ between neighboring superatoms, ${\rm{\Delta }}={V}_{{\rm{N}}}$. We neglect the interaction between the next to nearest neighbors. The conditions for validity of our treatment of the chain of superatoms and the parameters for suitable experimental systems are discussed in appendix B. In this so-called facilitation regime [18], an already excited superatom shifts the frequency of its nearest neighbor into resonance with the laser. The excitation ${{\rm{\Gamma }}}_{\uparrow }$ and de-excitation ${{\rm{\Gamma }}}_{\downarrow }$ rates for the facilitated superatom, having one and only one, excited neighbor, are

$$\begin{eqnarray}&&{{\rm{\Gamma }}}_{\uparrow }={{\rm{\Gamma }}}_{\mathrm{ex}}({{\rm{\Delta }}}_{\mathrm{eff}}=0),\end{eqnarray} \tag{ 7a }$$

$$\begin{eqnarray}&&{{\rm{\Gamma }}}_{\downarrow }={{\rm{\Gamma }}}_{\mathrm{de}}({{\rm{\Delta }}}_{\mathrm{eff}}=0).\end{eqnarray} \tag{ 7b }$$

On the other hand, a superatom surrounded by non-excited neighbors is non-resonant with the laser and has a much lower excitation rate

$$\begin{eqnarray}&&{{\rm{\Gamma }}}_{\mathrm{seed}}={{\rm{\Gamma }}}_{\mathrm{ex}}({{\rm{\Delta }}}_{\mathrm{eff}}={\rm{\Delta }}).\end{eqnarray} \tag{ 8 }$$

But once excited, it will play the role of a seed for a rapid growth of a cluster of excited superatoms. Naively, such clusters will grow until they collide and nearly all the superatoms in a lattice will be excited. However, spontaneous decay of superatoms with rate ${{\rm{\Gamma }}}_{s}$ will produce excitation holes—ground state superatoms surrounded by excited superatoms. A hole cannot be resonantly excited as its Rydberg state is shifted by the interaction with two excited neighbors by $2{V}_{{\rm{N}}}$, leading to the effective detuning ${{\rm{\Delta }}}_{\mathrm{eff}}={\rm{\Delta }}-2{V}_{{\rm{N}}}=-{\rm{\Delta }}$. Hence, the rate to refill the hole turns out to be the same highly suppressed seed rate ${{\rm{\Gamma }}}_{\mathrm{seed}}$.

In the facilitation regime that we consider, the typical hierarchy of the relevant rates is ${{\rm{\Gamma }}}_{\mathrm{seed}}\ll {{\rm{\Gamma }}}_{s},{{\rm{\Gamma }}}_{\downarrow }\ll {{\rm{\Gamma }}}_{\uparrow }$.

Upon turning on the excitation laser, the chain of superatoms under the facilitation conditions described above will exhibit transient dynamics on the timescale of $t\sim {{\rm{\Gamma }}}_{s}^{-1}$ and then settle in a steady state. The excitation transients will be described later in this section. First we discuss the steady state characterized by the average Rydberg excitation density of ${\rho }_{{\rm{R}}}\simeq 2/3$ and highly non-trivial equilibrium dynamics of the excitation holes.

3.1. Steady-state distribution of holes

By definition, a hole is a ground state superatom surrounded by two excited superatoms. Holes originate from collisions of growing clusters of Rydberg excitations during the transient and spontaneous decay of excited superatoms inside the cluster. In figure 3 we show the density of holes in a long lattice of superatoms. After switching on the excitation laser, within a few lifetimes ${{\rm{\Gamma }}}_{s}^{-1}$ of Rydberg excitations, the density of holes reaches an equilibrium. For large enough values of the interaction strength ${V}_{{\rm{N}}}\gg w$, and thereby the laser detuning $({\rm{\Delta }}={V}_{{\rm{N}}})$, the steady-state density of holes approaches the value of ${\rho }_{h}=1/3$.

Zoom In Zoom Out Reset image size

In figure 4 we show the spatial correlation function ${g}^{(2)}(d)$ for holes in the steady state. The two-particle correlation function is defined via

$$\begin{eqnarray*}&&{g}^{(2)}({d}_{k})\equiv \displaystyle \frac{\langle {\hat{n}}_{h}^{(j)}{\hat{n}}_{h}^{(j+k)}\rangle }{\langle {\hat{n}}_{h}^{(j)}{\rangle }^{2}},\end{eqnarray*}$$

where ${\hat{n}}_{h}^{(j)}$ is the hole number operator for site j of the lattice, ${d}_{k}={ak}$ with $k\in {\mathbb{N}}$ is the distance, and we assume that the average hole density $\langle {\hat{n}}_{h}^{(j)}\rangle ={\rho }_{h}$ is spatially uniform. We observe that, to a good approximation, holes behave as hard rods of length $2a$, with the average density of rods being 1/3.

Zoom In Zoom Out Reset image size

3.2. Effective model for holes

To understand the results of the numerical simulations for the hole density and correlation function, we have derived an effective model for the equilibrium dynamics of the holes. The derivation, details of which are given in appendix C, is based on adiabatic elimination of short-lived configurations involving two or more neighboring superatoms in the ground state. Such configurations appear when an excited superatom next to a hole decays to the ground state with rate ${{\rm{\Gamma }}}_{s(\downarrow )}$, but the lifetime of these configurations is very short, $t\sim {{\rm{\Gamma }}}_{\uparrow }^{-1}$, due to the excitation facilitation with the fast rate ${{\rm{\Gamma }}}_{\uparrow }\gg {{\rm{\Gamma }}}_{s(\downarrow )}$. Hence, in the effective model, holes cannot be located on the neighboring sites, which gives the physically intuitive picture as to why they behave as hard rods of length $2a$.

There are three fundamental processes affecting the many-body dynamics of the holes on the lattice, as illustrated in the inset of figure 3. We describe these processes in terms of the Lindblad jump operators acting in the corresponding subspace for the holes.

(i) Holes are created by spontaneous decay of excited superatoms with rate ${{\rm{\Gamma }}}_{c}={{\rm{\Gamma }}}_{s}$. The corresponding Lindblad operator is given by

$$\begin{eqnarray}&&{L}_{c}^{(j)}=\sqrt{{{\rm{\Gamma }}}_{c}}\,{\hat{\sigma }}_{+}^{(j)}[1-{\hat{n}}_{h}^{(j+1)}][1-{\hat{n}}_{h}^{(j-1)}],\end{eqnarray} \tag{ 9 }$$

where ${\hat{\sigma }}_{\pm }^{(j)}$ is the hole creation/annihilation operator for site j of the lattice, and ${\hat{n}}_{h}^{(j)}\equiv {\hat{\sigma }}_{+}^{(j)}{\hat{\sigma }}_{-}^{(j)}$ is the number operator. The last two terms on the right-hand side ensure that a hole cannot be created next to an existing one.

(ii) When two holes are separated by one excited superatom, one of the holes can be annihilated. This process is initiated by spontaneous decay of the middle superatom which triggers fast excitation of two superatoms from one or both sides of the three-site region. The hole annihilation is described by

$$\begin{eqnarray}&&{L}_{a}^{(j)}=\sqrt{{{\rm{\Gamma }}}_{a}/2}\,{\hat{\sigma }}_{+}^{(j)}{\hat{\sigma }}_{-}^{(j+1)}{\hat{\sigma }}_{-}^{(j-1)},\end{eqnarray} \tag{ 10a }$$

$$\begin{eqnarray}&&{L}_{{a}_{\pm }}^{(j)}=\sqrt{{{\rm{\Gamma }}}_{a}/4}\,{\hat{\sigma }}_{-}^{(j\pm 1)}{\hat{n}}_{h}^{(j)},\end{eqnarray} \tag{ 10b }$$

with the total annihilation probability given by ${{\rm{\Gamma }}}_{a}={{\rm{\Gamma }}}_{s}$. The remaining hole can then occupy either the middle site with half of the total probability, or one of the side sites, each with quarter of the total probability.

(iii) Holes can hop from site to site. This process is mediated by de-excitation of the superatom next to a hole with rate ${{\rm{\Gamma }}}_{\downarrow }$ followed by facilitated excitation of one of the ground state superatoms. The hole transport is described by

$$\begin{eqnarray}&&{L}_{{t}_{\pm }}^{(j)}=\sqrt{{{\rm{\Gamma }}}_{t}}\,{\hat{\sigma }}_{+}^{(j\pm 1)}{\hat{\sigma }}_{-}^{(j)}[1-{\hat{n}}_{h}^{(j\pm 2)}],\end{eqnarray} \tag{ 11 }$$

where ${{\rm{\Gamma }}}_{t}={{\rm{\Gamma }}}_{\downarrow }/2$ is the transport rate and the last term ensured that the hole cannot hop to a site next to an existing hole.

In figure 4 we compare the spatial correlation function ${g}^{(2)}(d)$ for ground state superatoms obtained from the numerical simulations of the full superatom model and the effective hole model. We observe very good agreement between the full and effective models, including ${g}^{(2)}(1)\simeq 0$.

3.3. Liquid-crystal crossover for lattice holes

Although in the present setup the hole hopping rate cannot be made arbitrary small, ${{\rm{\Gamma }}}_{t}\gtrsim {{\rm{\Gamma }}}_{c,a}$, it is instructive to analyze how varying ${{\rm{\Gamma }}}_{t}$ would affect the many-body steady state.

Consider first the hypothetical case of no hole transport, ${{\rm{\Gamma }}}_{t}\to 0$. We can also neglect for now the probability of refilling the hole, due to smallness of the corresponding rate ${{\rm{\Gamma }}}_{\mathrm{seed}}$. Then the only stable configuration corresponds to holes on every third site of the lattice, since neither the hole creation nor annihilation processes of equations (9) and (10) affect the system. Due to translational invariance of the lattice, this configuration is triple degenerate. The hole density–density correlation function ${g}^{(2)}(d)$ is peaked at $d=3a$, and the steady-state approaches a pure state with long-range crystalline order. Since the density of holes is exactly ${\rho }_{h}=1/3$, their number does not fluctuate. We can characterize the number fluctuations of the holes by the Mandel Q parameter

$$\begin{eqnarray*}&&Q\equiv \displaystyle \frac{\langle {\hat{n}}_{h}^{2}\rangle -\langle {\hat{n}}_{h}{\rangle }^{2}}{\langle {\hat{n}}_{h}\rangle }-1,\end{eqnarray*}$$

where ${\hat{n}}_{h}\equiv {\sum }_{j}{\hat{n}}_{h}^{(j)}$ is the total number of holes. $Q\lt 0$ signifies sub-Poissonian number distribution, with $Q=-1$ corresponding to a precise number of holes with no fluctuations.

As we now increase ${{\rm{\Gamma }}}_{t}$, the holes become mobile and the peak of the correlation function ${g}^{(2)}(d)$ at $d=3a$ is progressively reduced, see figure 5. The holes can now approach each other and annihilate, followed by hole creation on the allowed sites, which causes finite fluctuation of the hole number. Yet, their statistics remains highly sub-Poissonian even for large ${{\rm{\Gamma }}}_{t}\gg {{\rm{\Gamma }}}_{c,a}$. This is attested by the asymptotic value of $Q\simeq -0.81$ for ${{\rm{\Gamma }}}_{t}/{{\rm{\Gamma }}}_{c}\to \infty$ extracted from the algebraic fit in figure 5(b). Since the hole creation and annihilation rates are the same, ${{\rm{\Gamma }}}_{c}={{\rm{\Gamma }}}_{a}={{\rm{\Gamma }}}_{s}$, their mean density stays close to ${\rho }_{h}=1/3$. Interestingly, the crystal with periodicity $d=3a$ does not simply melt into a liquid with the same period. Rather, for a large hopping rate, the correlation function ${g}^{(2)}(d)$ exhibit period $d=2a$ with short correlation length $\xi \lesssim 1a$, see figure 4. Since holes cannot come closer than two lattice sites, they start to behave as mobile hard rods of length $2a$.

Zoom In Zoom Out Reset image size

We finally note that when the interaction strength ${V}_{{\rm{N}}}$, and thereby the laser detuning ${\rm{\Delta }}={V}_{{\rm{N}}}$, are not sufficiently larger than the excitation linewidth, the seed rate ${{\rm{\Gamma }}}_{\mathrm{seed}}$ is not negligible and holes can refill, leading to ${\rho }_{h}\lt 1/3$, as can be seen in figure 3.

3.4. Transient excitation dynamics of the system

As promised above, we now discuss the dynamics of the system of superatoms initially in the ground state upon switching on the excitation laser. In figure 6 we show the density ${\rho }_{{\rm{R}}}$ of Rydberg excitations as a function of time, for different values of the interaction strength ${V}_{{\rm{N}}}$ (and laser detuning ${\rm{\Delta }}={V}_{{\rm{N}}}$). In the long time limit $t\gg {{\rm{\Gamma }}}_{s}^{-1}$, for large enough ${V}_{{\rm{N}}}\gg w$, the system reaches the steady state with the Rydberg excitation density ${\rho }_{{\rm{R}}}\simeq 2/3$, consistent with the hole density ${\rho }_{h}\simeq 1/3$ analyzed above. For smaller values of ${V}_{{\rm{N}}}$, we have in the steady state ${\rho }_{{\rm{R}}}\gt 2/3$, since holes can be refilled with appreciable seed rate ${{\rm{\Gamma }}}_{\mathrm{seed}}$.

Zoom In Zoom Out Reset image size

Consider now the transient regime $0\lt t\lesssim {{\rm{\Gamma }}}_{s}^{-1}$. In the simulations of figure 6, we put initially one seed excitation in order to start the facilitated excitation dynamics at a well-defined time, rather than at some uncertain (within $1/{{\rm{\Gamma }}}_{\mathrm{seed}}$) time when the seed excitation appears with small rate ${{\rm{\Gamma }}}_{\mathrm{seed}}$. As seen in that figure, after switching on the laser, the density of Rydberg excitations can peak at a large value and then relax to the lower steady-state density. At first sight, this observation is surprising, since such a behavior is reminiscent to partially coherent dynamics of a quantum system, such as, e.g., damped Rabi oscillations, while here our system is completely governed by rate equations and no coherences are involved. We now outline a macroscopic model that will explain the nature of the peak and the associated peak time t₀.

We consider three basic states of the system: the ground state (g), the fully excited state (e) and the final steady state (s). The corresponding probabilities are denoted by ${p}_{{\rm{g}}}$, ${p}_{{\rm{e}}}$, and ${p}_{{\rm{s}}}$. Initially the system is in the ground state, ${p}_{{\rm{g}}}=1$, and to start the dynamics we need a seed excitation. The probability of the seed p_seed is governed by the equation ${\partial }_{t}{p}_{\mathrm{seed}}(t)={p}_{{\rm{g}}}{{\rm{\Gamma }}}_{\mathrm{seed}}$. For short times, we can assume ${p}_{{\rm{g}}}\simeq 1$ and obtain linear growth of the seed probability, ${p}_{\mathrm{seed}}(t)\simeq {p}_{\mathrm{seed}}(0)+{{\rm{\Gamma }}}_{\mathrm{seed}}t$. For longer times, this approximation breaks down, but once ${p}_{{\rm{g}}}$ is depleted, the role of the seed becomes unimportant, as will become clear shortly.

The equations of motion for the probabilities of the three basic states are

$$\begin{eqnarray}&&{\partial }_{t}{p}_{{\rm{g}}}=-{p}_{\mathrm{seed}}\,2{{\rm{\Gamma }}}_{\uparrow }\,{p}_{{\rm{g}}},\end{eqnarray} \tag{ 12a }$$

$$\begin{eqnarray}&&{\partial }_{t}{p}_{{\rm{e}}}=+{p}_{\mathrm{seed}}\,2{{\rm{\Gamma }}}_{\uparrow }\,{p}_{{\rm{g}}}-3{{\rm{\Gamma }}}_{s}\,{p}_{{\rm{s}}}+{{\rm{\Gamma }}}_{\mathrm{seed}}\,{p}_{{\rm{s}}},\end{eqnarray} \tag{ 12b }$$

$$\begin{eqnarray}&&{\partial }_{t}{p}_{{\rm{s}}}=+3{{\rm{\Gamma }}}_{s}\,{p}_{{\rm{s}}}-{{\rm{\Gamma }}}_{\mathrm{seed}}\,{p}_{{\rm{s}}}.\end{eqnarray} \tag{ 12c }$$

Each seed excitation in the lattice triggers fast growth of facilitated excitations in both lattice directions. The ground state is then being depleted with the rate $2{{\rm{\Gamma }}}_{\uparrow }$ and its population is transferred to the fully excited state. In turn, the Rydberg excitations in the fully excited state decay with rate ${{\rm{\Gamma }}}_{s}$. Since the steady state corresponds to configurations with, on average, every third site non-excited, the decay rate of ${p}_{{\rm{e}}}$ to ${p}_{{\rm{s}}}$ is $3{{\rm{\Gamma }}}_{s}$, which we have also verified via numerical simulations of the full system with all the superatoms initially excited. We also include in the above equations the process of refilling the holes with the corresponding rate ${{\rm{\Gamma }}}_{\mathrm{seed}}$. Then the solution of these equations approximates remarkably well the mean density of Rydberg excitations ${\rho }_{{\rm{R}}}(t)={p}_{{\rm{e}}}(t)+\tfrac{2}{3}{p}_{{\rm{s}}}(t)$, even in the regime of sizable ${{\rm{\Gamma }}}_{\mathrm{seed}}$, as can be seen in the main panel of figure 6.

The peak in the excitation density ${\rho }_{{\rm{R}}}$ is reached when the ground state probability ${p}_{{\rm{g}}}$ is depleted and the majority of superatoms are transferred to the excited state, ${p}_{{\rm{e}}}\simeq 1$. Integration of ${\partial }_{t}{p}_{{\rm{g}}}$ suggests the peak time scaling as ${t}_{0}\propto 1/\sqrt{{{\rm{\Gamma }}}_{\mathrm{seed}}{{\rm{\Gamma }}}_{\uparrow }}\propto {V}_{{\rm{N}}}$. The scaling of the peak time with the interaction strength ${V}_{{\rm{N}}}/w$ is verified in the inset of figure 6. Note finally that our macroscopic model neglects cluster collisions, also producing holes, and thereby slightly overestimates the excitation density ${\rho }_{{\rm{R}}}$ obtained from the simulations of the complete microscopic model.

To summarize, we have studied the excitation dynamics and steady state of a lattice of Rydberg superatoms driven by a laser in the facilitation regime. We have shown that the steady state of the system has nearly universal Rydberg excitation density of ${\rho }_{{\rm{R}}}=2/3$. More interestingly, it corresponds to an equilibrium dynamics of mobile quasi-particles—excitation holes. We have derived an effective hole model which involves hole creation and pair annihilation of holes separated by two lattice sites. We have found that the number fluctuations of the holes are characterized by the Mandel Q parameter $Q\simeq -0.81$, and their spatial correlations decay on a distance of $\xi \lesssim 1a$ comparable to the lattice constant a.

That negative values of Q, which signify pronounced sub-Poissonian number distribution, have their origin in the hard rod constraint, has been pointed out already in [11] considering a Gibbs ensemble of Rydberg excitations. In our system, however, the pair annihilation of holes leads to a state with much stronger suppressed fluctuations. In the Gibbs state of [11], we would have to choose a rod length of about d = 0.611 to obtain the average hole density ${\rho }_{h}=1/3$, but this would then lead to $Q\simeq -0.63$ which is significantly larger than what we obtain from our simulations.

Our model corresponds to an experimentally amenable regime with a large hole transport rate ${{\rm{\Gamma }}}_{t}\sim {{\rm{\Gamma }}}_{c,a}$ comparable to the hole creation and annihilation rates ${{\rm{\Gamma }}}_{c,a}={{\rm{\Gamma }}}_{s}$, which, in turn, are determined by the spontaneous decay rate of Rydberg excitations. The resulting spatial correlations of the holes have a period of $d=2a$ and short correlation length $\xi \simeq a$. If we freeze the hole motion, ${{\rm{\Gamma }}}_{t}\ll {{\rm{\Gamma }}}_{s}$, we would obtain long-range order, $\xi \gg a$, with $d=3a$ periodicity. Such a regime can in principle be achieved, but at the expense of more complicated experimental setup involving additional lasers coupling the ground state of superatoms to a different long lived Rydberg or metastable state which would make the holes immobile.

Another interesting direction of research is to consider different lattice geometries, e.g., a two-dimensional square or triangular lattice. The latter might simulate dissipative frustrated spin models. Finally, the rate equation approach, amenable to large scale numerical simulations, is applicable in the regime of strong dephasing. Coherence effects might lead to interesting dynamics and yield long-range correlations and entanglement. However, fully quantum many-body simulations are limited to small system sizes.

This work was supported by DFG through SFB/TR185 (MF and FL) and the H2020 FET Proactive project RySQ (DP). FL is supported by a fellowship through the Excellence Initiative MAINZ (DFG/GSC 266). DP is grateful to the University of Kaiserslautern for hospitality and to the Alexander von Humboldt Foundation for support during his stay in Germany.

Here we consider in some detail the dissipative dynamics of a single Rydberg superatom. The superatom consists of N two-level atoms within a volume of linear dimension $\delta r$ smaller than the Rydberg blockade distance ${a}_{{\rm{B}}}$. We define ${a}_{{\rm{B}}}$ as the distance below which the interatomic interaction strength $V(r\leqslant {a}_{{\rm{B}}})\geqslant w$ starts to exceed the steady-state excitation linewidth of an atom, $w=\sqrt{4{{\rm{\Omega }}}^{2}\gamma /{{\rm{\Gamma }}}_{s}+{\gamma }^{2}}\simeq 2{\rm{\Omega }}\sqrt{\gamma /{{\rm{\Gamma }}}_{s}}$ [13, 14]. Here $\gamma \equiv \tfrac{1}{2}({{\rm{\Gamma }}}_{d}+{{\rm{\Gamma }}}_{s})$, and we assume ${{\rm{\Omega }}}^{2}\gg {{\rm{\Gamma }}}_{s}{{\rm{\Gamma }}}_{d}$. Due to the strong Rydberg blockade, the superatom can accommodate at most one Rydberg excitation.

The collective ground state of a superatom $| G\rangle \equiv | {g}_{1}{g}_{2}\ldots {g}_{N}\rangle$ is coherently coupled by the laser to the symmetric single excitation state $| {E}_{s}\rangle \equiv \tfrac{1}{\sqrt{N}}{\sum }_{k}| {g}_{1}{g}_{2}\ldots {e}_{k}\ldots {g}_{N}\rangle$ with Rabi frequency $\sqrt{N}{\rm{\Omega }}$, see figure 2(a). In addition, there are $N-1$ non-symmetric single excitation states $\{| {E}_{{ns}}{\rangle }_{m}\}$ which are not directly coupled to the ground state by the laser. All the excited states $\{| E\rangle \}$ spontaneously decay to the ground state $| G\rangle$ with rate ${{\rm{\Gamma }}}_{s}$. The dephasing ${{\rm{\Gamma }}}_{d}$ of the atomic Rydberg state $| e\rangle$, with respect to the ground state $| g\rangle$, leads to incoherent coupling of any single excitation state $| E\rangle$ of the superatom to any other such state $| E^{\prime} \rangle$ with rate ${{\rm{\Gamma }}}_{d}/N$. We may replace the manifold $\{| {E}_{{ns}}{\rangle }_{m}\}$ with a single aggregate non-symmetric state $| {E}_{{ns}}\rangle$, obtaining for the superatom an effective three-level system, $\{| G\rangle ,| {E}_{s}\rangle ,| {E}_{{ns}}\rangle \}$. We describe this dissipative system with a 'vector' of density matrix elements $\vec{\rho }={({\rho }_{{GG}},{\rho }_{{E}_{s}{E}_{s}},{\rho }_{{{GE}}_{s}},{\rho }_{{E}_{s}G},{\rho }_{{E}_{{ns}}{E}_{{ns}}})}^{T}$ which obeys the equation of motion

$$\begin{eqnarray}&&{\partial }_{t}\vec{\rho }={\rm{\Lambda }}\vec{\rho },\end{eqnarray} \tag{ A1 }$$

where

$$\begin{eqnarray}{\rm{\Lambda }}=\left(\begin{array}{ccccc}0 & {{\rm{\Gamma }}}_{s} & {\rm{i}}\sqrt{N}{\rm{\Omega }} & -{\rm{i}}\sqrt{N}{\rm{\Omega }} & {{\rm{\Gamma }}}_{s}\\ 0 & -{{\rm{\Gamma }}}_{s}-\displaystyle \frac{N-1}{N}{{\rm{\Gamma }}}_{d} & -{\rm{i}}\sqrt{N}{\rm{\Omega }} & {\rm{i}}\sqrt{N}{\rm{\Omega }} & \displaystyle \frac{1}{N}{{\rm{\Gamma }}}_{d}\\ {\rm{i}}\sqrt{N}{\rm{\Omega }} & -{\rm{i}}\sqrt{N}{\rm{\Omega }} & -{\rm{i}}{\rm{\Delta }}-\gamma & 0 & 0\\ -{\rm{i}}\sqrt{N}{\rm{\Omega }} & {\rm{i}}\sqrt{N}{\rm{\Omega }} & 0 & {\rm{i}}{\rm{\Delta }}-\gamma & 0\\ 0 & \displaystyle \frac{N-1}{N}{{\rm{\Gamma }}}_{d} & 0 & 0 & -{{\rm{\Gamma }}}_{s}-\displaystyle \frac{1}{N}{{\rm{\Gamma }}}_{d}\end{array}\right).\end{eqnarray} \tag{ A2 }$$

For N = 1, we retrieve the well known optical Bloch equations for a single two-level atom, with the decay rate ${{\rm{\Gamma }}}_{s}$ of the excited state population ${\rho }_{{ee}}$ and the relaxation rate $\gamma =\tfrac{1}{2}({{\rm{\Gamma }}}_{s}+{{\rm{\Gamma }}}_{d})$ of the coherence ${\rho }_{{ge}}$. For a superatom with $N\gt 1$, dephasing induces population ${\rho }_{{E}_{{ns}}{E}_{{ns}}}\gt 0$ of non-symmetric state(s), but their coherences remain decoupled, ${\rho }_{{E}_{{ns}}G}={\rho }_{{E}_{{ns}}{E}_{s}}=0$. By definition, the superatom can contain at most one Rydberg excitation, ${\rho }_{{GG}}+{\rho }_{{E}_{s}{E}_{s}}+{\rho }_{{E}_{{ns}}{E}_{{ns}}}=1$.

Solving equation (A1) in the steady state, ${\partial }_{t}\vec{\rho }=0$, we obtain the following expression for the total excited state population of the superatom,

$$\begin{eqnarray}{\rho }_{{EE}} & \equiv & {\rho }_{{E}_{s}{E}_{s}}+{\rho }_{{E}_{{ns}}{E}_{{ns}}}\\ & = & \displaystyle \frac{2N{{\rm{\Omega }}}^{2}\gamma }{{{\rm{\Omega }}}^{2}(N+1)\left({{\rm{\Gamma }}}_{d}+\tfrac{2N}{N+1}{{\rm{\Gamma }}}_{s}\right)+{{\rm{\Gamma }}}_{s}({\gamma }^{2}+{{\rm{\Delta }}}^{2})}.\end{eqnarray} \tag{ A3 }$$

For N = 1, this reduces to the excited state population of a single two-level atom [13],

$$\begin{eqnarray}&&{\rho }_{{ee}}=\displaystyle \frac{{{\rm{\Omega }}}^{2}}{2{{\rm{\Omega }}}^{2}+\tfrac{{{\rm{\Gamma }}}_{s}}{2\gamma }({\gamma }^{2}+{{\rm{\Delta }}}^{2})},\end{eqnarray} \tag{ A4 }$$

which is bounded by ${\rho }_{{ee}}\leqslant 0.5$ even for resonant excitation $| {\rm{\Delta }}| \ll \gamma$. On the other hand, for large N and small decay rate ${{\rm{\Gamma }}}_{s}\ll {{\rm{\Gamma }}}_{d},{{\rm{\Omega }}}^{2}/{{\rm{\Gamma }}}_{d}$, the excited state population of the superatom is approximately given by

$$\begin{eqnarray}&&{\rho }_{{EE}}\simeq \displaystyle \frac{N}{N+1},\end{eqnarray} \tag{ A5 }$$

which quickly approaches ${\rho }_{{EE}}\simeq 1$ with increasing N. We emphasize that the steady-state population inversion of the superatom is brought about by strong driving and dephasing, which tend to equalize populations of all $N+1$ states, i.e., N single excitation states and the ground state. In contrast, without dephasing (${{\rm{\Gamma }}}_{d}\ll {{\rm{\Gamma }}}_{s}$), the superatom reduces to a two-level system, with the ground state $| G\rangle$ coherently coupled to the symmetric excited state $| {E}_{s}\rangle$ which decays back to the ground state with rate ${{\rm{\Gamma }}}_{s}$, and the resulting steady-state Rydberg excitation probability is ${\rho }_{{EE}}\simeq {\rho }_{{E}_{s}{E}_{s}}\leqslant 1/2$.

From equation (A3), with strong dephasing ${{\rm{\Gamma }}}_{d}\gg {{\rm{\Gamma }}}_{s}$ and $N\gg 1$, we obtain the excitation linewidth of the superatom

$$\begin{eqnarray}&&{w}_{\mathrm{SA}}\simeq \sqrt{2N{{\rm{\Omega }}}^{2}\gamma /{{\rm{\Gamma }}}_{s}+{\gamma }^{2}}.\end{eqnarray} \tag{ A6 }$$

Comparing it with the single-atom excitation linewidth w which follows from equation (A4), we notice the analogy upon replacement ${{\rm{\Omega }}}^{2}\to \tfrac{1}{2}N{{\rm{\Omega }}}^{2}$, rather than ${{\rm{\Omega }}}^{2}\to N{{\rm{\Omega }}}^{2}$ as one would naively expect from the collective enhancement of the N-atom Rabi frequency $\sqrt{N}{\rm{\Omega }}$. This factor of 2 difference stems from the fact that the Rydberg excitation probability of a saturated superatom approaches unity, rather than $\tfrac{1}{2}$.

For strong dephasing ${{\rm{\Gamma }}}_{d}\gtrsim {\rm{\Omega }}$, the dynamics of a superatom can be described, to a good approximation, by rate equations which we now derive. In the equations for the density matrix elements, we adiabatically eliminate the coherences ${\rho }_{{E}_{s}G}$, ${\rho }_{{{GE}}_{s}}$, obtaining rate equations ${\partial }_{t}{\vec{\rho }}_{3}={{\rm{\Lambda }}}_{3}{\vec{\rho }}_{3}$ for populations ${\vec{\rho }}_{3}={({\rho }_{{GG}},{\rho }_{{E}_{s}{E}_{s}},{\rho }_{{E}_{{ns}}{E}_{{ns}}})}^{T}$ with

$$\begin{eqnarray}{{\rm{\Lambda }}}_{3}=\left(\begin{array}{ccc}-N\chi & N\chi +{{\rm{\Gamma }}}_{s} & {{\rm{\Gamma }}}_{s}\\ N\chi & -N\chi -{{\rm{\Gamma }}}_{s}-\tfrac{N-1}{N}{{\rm{\Gamma }}}_{d} & \tfrac{1}{N}{{\rm{\Gamma }}}_{d}\\ 0 & \tfrac{N-1}{N}{{\rm{\Gamma }}}_{d} & -{{\rm{\Gamma }}}_{s}-\tfrac{1}{N}{{\rm{\Gamma }}}_{d}\end{array}\right),\end{eqnarray} \tag{ A7 }$$

where

$$\begin{eqnarray*}&&\chi \equiv \displaystyle \frac{2{{\rm{\Omega }}}^{2}\gamma }{{\gamma }^{2}+{{\rm{\Delta }}}^{2}}.\end{eqnarray*}$$

The laser tends to equalize populations ${\rho }_{{GG}}$ and ${\rho }_{{E}_{s}{E}_{s}}$ of the ground and symmetric excited states. In turn, strong dephasing quickly transfers the population of the symmetric state to the non-symmetric states. For large N, the reverse transfer from ${\rho }_{{E}_{{ns}}{E}_{{ns}}}$ to ${\rho }_{{E}_{s}{E}_{s}}$ is suppressed by a factor of $1/N$. The symmetric state then plays the role of an intermediate state having small population ${\rho }_{{E}_{s}{E}_{s}}\sim {\rho }_{{E}_{{ns}}{E}_{{ns}}}/N$. Upon adiabatic elimination of ${\rho }_{{E}_{s}{E}_{s}}$ we finally obtain rate equations ${\partial }_{t}{\vec{\rho }}_{2}={{\rm{\Lambda }}}_{2}{\vec{\rho }}_{2}$ for the populations ${\vec{\rho }}_{2}={({\rho }_{{GG}},{\rho }_{{EE}})}^{T}$ of the ground and excited states (${\rho }_{{GG}}+{\rho }_{{EE}}\simeq 1$),

$$\begin{eqnarray}{{\rm{\Lambda }}}_{2}=\left(\begin{array}{cc}-{{\rm{\Gamma }}}_{\mathrm{ex}} & {{\rm{\Gamma }}}_{\mathrm{de}}\\ {{\rm{\Gamma }}}_{\mathrm{ex}} & -{{\rm{\Gamma }}}_{\mathrm{de}}\end{array}\right),\end{eqnarray} \tag{ A8 }$$

where the excitation and de-excitation rates are given by

$$\begin{eqnarray}&&{{\rm{\Gamma }}}_{\mathrm{ex}}=\displaystyle \frac{(N-1)\chi }{N\chi +2\gamma }{{\rm{\Gamma }}}_{d},\end{eqnarray} \tag{ A9 }$$

$$\begin{eqnarray}&&{{\rm{\Gamma }}}_{\mathrm{de}}=\displaystyle \frac{\chi }{N\chi +2\gamma }{{\rm{\Gamma }}}_{d}+{{\rm{\Gamma }}}_{s}.\end{eqnarray} \tag{ A10 }$$

This is illustrated in figure 2(b).

In figure A1(a) we compare the excitation dynamics of a resonantly driven superatom as obtained from the solution of the rate equations and the exact solution of the master equation for the density operator using the Monte Carlo wavefunction approach in the truncated single excitation Hilbert space. We observe that the rate equations provide accurate description of the dynamics of the system when $\gamma \gtrsim \sqrt{N}{\rm{\Omega }}$, while they cannot account for the (damped) Rabi oscillations between the ground and the excited states when $\sqrt{N}{\rm{\Omega }}\gg \gamma$. Nevertheless, the rate equations model approximates well the relaxation timescale and the steady-state population of the superatom.

Zoom In Zoom Out Reset image size

In figure A1(b) we show the spectrum of eigenvalues ${\lambda }_{k}$ of the matrix Λ in equation (A2). The eigenvalue ${\lambda }_{0}=0$ corresponds to the steady state of the system, while the (negative) real parts of the other eigenvalues, $\mathrm{Re}[{\lambda }_{k}]$ for $k=1,2,3,4$, characterize the relaxation rates of the superatom towards the steady state. In that figure, we also show the total relaxation rate of a superatom towards the steady state as given by the rate equations model, ${{\rm{\Gamma }}}_{\mathrm{tot}}={{\rm{\Gamma }}}_{\mathrm{ex}}+{{\rm{\Gamma }}}_{\mathrm{de}}$, which compares favorably with the exact relaxation rate for a broad range of parameters.

With the interatomic potential $V(r)={C}_{p}/{r}^{p}$ (assuming repulsive interaction ${C}_{p}\gt 0$) and the single-atom excitation linewidth w, the Rydberg blockade distance is defined as ${a}_{{\rm{B}}}\equiv {({C}_{p}/w)}^{1/p}$. Our model assumes that each superatom can accommodate at most one Rydberg excitation. We therefore require the spatial extent of a superatom to be small compared to blockade distance, $\delta r\ll {a}_{{\rm{B}}}$.

Clearly, the superatom excitation probability is maximal at resonant driving and is suppressed when the laser is detuned by more than its excitation linewidth ${w}_{\mathrm{SA}}\simeq w\sqrt{N/2}$. We consider a one-dimensional chain of superatoms driven by spatially uniform laser with detuning Δ. We set the laser detuning to be equal to the interaction strength between neighboring superatoms separated by the lattice constant a, ${\rm{\Delta }}=V(a)$. Then an already excited superatom will facilitate the excitation of the nearest neighbors by shifting the Rydberg energy level into resonance with the laser field. On the other hand, we require that ${\rm{\Delta }}\gt {w}_{\mathrm{SA}}$ in order to suppress excitation of superatoms that have either two non-excited neighbors or two excited neighbors. Since ${w}_{\mathrm{SA}}\gt w$, the facilitation distance ${r}_{\mathrm{fac}}={({C}_{p}/{\rm{\Delta }})}^{1/p}$, and thereby the lattice constant $a={r}_{\mathrm{fac}}$, are smaller than the blockade distance a_B.

The interaction potential V(r) is a convex monotonic function. Our approach is valid when the interaction-induced Rydberg level shifts of all the atoms of the facilitated superatom are within the superatom linewidth w_SA. Linearizing the interaction potential around the facilitation distance r_fac, we obtain the following condition on the spatial extent $\delta r$ of the superatom

$$\begin{eqnarray}&&\delta r\lesssim \displaystyle \frac{1}{p}\displaystyle \frac{{r}_{\mathrm{fac}}^{p+1}}{{C}_{p}/{w}_{\mathrm{SA}}}=\displaystyle \frac{1}{p}\displaystyle \frac{{w}_{\mathrm{SA}}}{{\rm{\Delta }}}a.\end{eqnarray} \tag{ B1 }$$

Finally, in our treatment of the 1D lattice of superatoms, we neglect the interaction shifts of the next neighbors, $V(2a)\lesssim {w}_{\mathrm{SA}}$. This can be rewritten as a condition on the power-law scaling p of the interaction potential,

$$\begin{eqnarray}&&p\mathrm{log}(2)\geqslant \mathrm{log}({\rm{\Delta }}/{w}_{\mathrm{SA}}).\end{eqnarray} \tag{ B2 }$$

To relate the forgoing discussion to a realistic experiment, we consider a lattice of microtraps each containing $N\simeq 50$ cold Rb atoms. The atomic ground state $| g\rangle \equiv 5{S}_{1/2}| F=2,{M}_{F}=-2\rangle$ is coupled, by a two-photon transition via the intermediate state $5{P}_{3/2}| F=3,{M}_{F}=-3\rangle$ (or $6{P}_{3/2}$) [15, 42], to the highly excited Rydberg state $| e\rangle \equiv {{nS}}_{1/2}$ with the principal quantum number $n\simeq 90$. The decay rate of state $| e\rangle$ is ${{\rm{\Gamma }}}_{s}\simeq 5\$ kHz. With the corresponding Rabi frequency ${\rm{\Omega }}=25{{\rm{\Gamma }}}_{s}=2\pi \times 20\$ kHz and the dephasing rate ${{\rm{\Gamma }}}_{d}=250{{\rm{\Gamma }}}_{s}$, we have the single-atom and superatom excitation linewidths $w=570{{\rm{\Gamma }}}_{s}=2\pi \times 450\$ kHz and ${w}_{\mathrm{SA}}=5w$. We chose the detuning ${\rm{\Delta }}=2\pi \times 10\$ MHz. With the van der Waals coefficient ${C}_{6}=2\pi \times 16.8\$ THz $\ \mu {\rm{m}}$⁶ [43], we have the lattice constant $a\simeq 11\ \mu {\rm{m}}$, while the size of the superatom should be $\delta r\lesssim 0.4\ \mu {\rm{m}}$ as per equation (B1) with p = 6. This implies large atomic density ${10}^{15}\$ cm⁻³, which would lead to strong Rydberg state dephasing [39–41] (assumed included in ${{\rm{\Gamma }}}_{d}$), and to collisional dephasing and loss of the ground state atoms from the microtraps (see below). Alternatively, we can employ the Rydberg state $| e\rangle$ corresponding to the stretched Stark eigenstate in the static electric field in the perpendicular to the array direction [9]. Such states possess large permanent dipole moments ($\sim \tfrac{3}{2}{n}^{2}{a}_{0}e$) leading to strong dipole–dipole interactions. With n = 60 we obtain the dipole–dipole coefficient ${C}_{3}\simeq 2\pi \times 28.5\$GHz $\ \mu {\rm{m}}$³ resulting in the lattice constant $a\simeq 14\ \mu {\rm{m}}$ for the same Δ as above. The size of the superatom can now be larger, $\delta r\sim 1\ \mu {\rm{m}}$ as per equation (B1) with p = 3, requiring lower atomic density $5\times {10}^{13}\$ cm⁻³ at which the collisional dephasing and losses are small. In both cases of p = 6 and 3 the condition (B2) is satisfied.

We finally note that loss of several atoms from the superatoms with $N\gg 1$ will not significantly affect the system dynamics on the equilibration timescale of $t\gtrsim {{\rm{\Gamma }}}_{s}^{-1}$.

Consider a small chain of superatoms shown in figure C1. We focus on a subset of the system containing three superatoms that can be in the ground state or excited to the Rydberg state. We denote the probability of each configuration $\{a,b,c\}$ by p_abc, where $a,b,c=0,1$ for non-excited (0) and excited (1) superatom on the corresponding site. Starting with all the superatoms excited, $\{1,1,1\}$, we can create a hole on any one site, e.g., $\{0,1,1\}$, with the spontaneous decay rate ${{\rm{\Gamma }}}_{s}$. The reverse process of exciting a hole can be neglected due to smallness of the transition rate ${{\rm{\Gamma }}}_{\mathrm{seed}}\ll {{\rm{\Gamma }}}_{\uparrow ,\downarrow },{{\rm{\Gamma }}}_{s}$. Next, an excited superatom neighboring a non-excited one can be de-excited, e.g., $\{0,1,1\}\to \{0,0,1\}$, with rate ${{\rm{\Gamma }}}_{\downarrow }$. An excited superatom between two non-excited ones can be de-excited, e.g., $\{0,1,0\}\to \{0,0,0\}$, with the spontaneous decay rate ${{\rm{\Gamma }}}_{s}$. Finally, the probability of configuration with two or more neighboring non-excited superatoms will quickly decay, e.g., $\{0,0,1\}\to \{1,0,1\}$ or $\{0,1,1\}$, with the the facilitated excitation rate ${{\rm{\Gamma }}}_{\uparrow }$. We thus obtain the following set of equations for the probabilities of various configurations (see figure C1 upper panel),

$$\begin{eqnarray}&&{\partial }_{t}{p}_{111}=-3{{\rm{\Gamma }}}_{s}{p}_{111},\end{eqnarray} \tag{ C1 }$$

$$\begin{eqnarray}&&{\partial }_{t}{p}_{011}=-({{\rm{\Gamma }}}_{s}+{{\rm{\Gamma }}}_{\downarrow }){p}_{011}+{{\rm{\Gamma }}}_{s}{p}_{111}+{{\rm{\Gamma }}}_{\uparrow }{p}_{001},\end{eqnarray} \tag{ C2 }$$

$$\begin{eqnarray}&&{\partial }_{t}{p}_{101}=-2{{\rm{\Gamma }}}_{\downarrow }{p}_{101}+{{\rm{\Gamma }}}_{s}{p}_{111}+{{\rm{\Gamma }}}_{\uparrow }{p}_{001}+{{\rm{\Gamma }}}_{\uparrow }{p}_{100},\end{eqnarray} \tag{ C3 }$$

$$\begin{eqnarray}&&{\partial }_{t}{p}_{110}=-({{\rm{\Gamma }}}_{s}+{{\rm{\Gamma }}}_{\downarrow }){p}_{110}+{{\rm{\Gamma }}}_{s}{p}_{111}+{{\rm{\Gamma }}}_{\uparrow }{p}_{100},\end{eqnarray} \tag{ C4 }$$

$$\begin{eqnarray}&&{\partial }_{t}{p}_{010}=-{{\rm{\Gamma }}}_{s}({p}_{010}-{p}_{011}-{p}_{110}),\end{eqnarray} \tag{ C5 }$$

$$\begin{eqnarray}&&{\partial }_{t}{p}_{001}=-(2{{\rm{\Gamma }}}_{\uparrow }+{{\rm{\Gamma }}}_{\downarrow }){p}_{001}+{{\rm{\Gamma }}}_{\uparrow }{p}_{000}+{{\rm{\Gamma }}}_{\downarrow }({p}_{011}+{p}_{101}),\end{eqnarray} \tag{ C6 }$$

$$\begin{eqnarray}&&{\partial }_{t}{p}_{100}=-(2{{\rm{\Gamma }}}_{\uparrow }+{{\rm{\Gamma }}}_{\downarrow }){p}_{100}+{{\rm{\Gamma }}}_{\uparrow }{p}_{000}+{{\rm{\Gamma }}}_{\downarrow }({p}_{110}+{p}_{101}),\end{eqnarray} \tag{ C7 }$$

$$\begin{eqnarray}&&{\partial }_{t}{p}_{000}=-2{{\rm{\Gamma }}}_{\uparrow }{p}_{000}+{{\rm{\Gamma }}}_{\downarrow }({p}_{001}+{p}_{100})+{{\rm{\Gamma }}}_{s}{p}_{010}.\end{eqnarray} \tag{ C8 }$$

The probability of full excitation p₁₁₁ decays to zero on the timescale $t\gt {{\rm{\Gamma }}}_{s}^{-1}$. The probabilities in the last three equations relax with the very fast rates $\sim 2{{\rm{\Gamma }}}_{\uparrow }$ and can be adiabatically eliminated. Setting there ${\partial }_{t}{p}_{{abc}}=0$, and substituting the resulting solutions in the other equations, we obtain

$$\begin{eqnarray}&&{\partial }_{t}{p}_{011}=-\left({{\rm{\Gamma }}}_{s}+\tfrac{1}{2}{{\rm{\Gamma }}}_{\downarrow }\right){p}_{011}+\tfrac{1}{2}{{\rm{\Gamma }}}_{\downarrow }{p}_{101}+\tfrac{1}{4}{{\rm{\Gamma }}}_{s}{p}_{010},\end{eqnarray} \tag{ C9 }$$

$$\begin{eqnarray}&&{\partial }_{t}{p}_{101}=-{{\rm{\Gamma }}}_{\downarrow }{p}_{101}+\tfrac{1}{2}{{\rm{\Gamma }}}_{s}{p}_{010}+\tfrac{1}{2}{{\rm{\Gamma }}}_{\downarrow }({p}_{011}+{p}_{110}),\end{eqnarray} \tag{ C10 }$$

$$\begin{eqnarray}&&{\partial }_{t}{p}_{110}=-\left({{\rm{\Gamma }}}_{s}+\tfrac{1}{2}{{\rm{\Gamma }}}_{\downarrow }\right){p}_{110}+\tfrac{1}{2}{{\rm{\Gamma }}}_{\downarrow }{p}_{101}+\tfrac{1}{4}{{\rm{\Gamma }}}_{s}{p}_{010},\end{eqnarray} \tag{ C11 }$$

$$\begin{eqnarray}&&{\partial }_{t}{p}_{010}=-{{\rm{\Gamma }}}_{s}{p}_{010}+{{\rm{\Gamma }}}_{s}({p}_{011}+{p}_{110}),\end{eqnarray} \tag{ C12 }$$

where we assumed $\tfrac{{{\rm{\Gamma }}}_{\downarrow }}{{{\rm{\Gamma }}}_{s}}{p}_{001},\tfrac{{{\rm{\Gamma }}}_{\downarrow }}{{{\rm{\Gamma }}}_{s}}{p}_{100}\ll {p}_{010}$ and ${{\rm{\Gamma }}}_{\downarrow }\ll 2{{\rm{\Gamma }}}_{\uparrow }$.

Zoom In Zoom Out Reset image size

By adiabatic elimination of the states with neighboring ground state superatoms, we projected the system onto a reduced configuration space. From the above rate equations, we can deduce three fundamental processes affecting the holes (see figure C1 lower panel):

• (i)  
  creation of holes with the rate ${{\rm{\Gamma }}}_{c}={{\rm{\Gamma }}}_{s}$,
• (ii)  
  annihilation of one of the two holes separated by one excitation with the rate ${{\rm{\Gamma }}}_{a}={{\rm{\Gamma }}}_{s}$,
• (iii)  
  transport of holes between neighboring sites with rate ${{\rm{\Gamma }}}_{t}=\tfrac{1}{2}{{\rm{\Gamma }}}_{\downarrow }$.

The corresponding jump operators, involving the constraints that no two neighboring holes are allowed, are given in section 3.1.