Griffiths Phase in a Facilitated Rydberg Gas at Low Temperature

TL;DR

This paper investigates the dynamics of facilitated Rydberg gases at different temperatures, revealing a Griffiths phase at low temperatures and phase transition behaviors at high temperatures, with implications for understanding many-body quantum systems.

Contribution

It introduces a detailed analysis of Rydberg facilitation dynamics across temperature regimes, highlighting the emergence of a Griffiths phase and extending Langevin models for better accuracy.

Findings

At high temperatures, a homogeneous mean-field behavior is observed.

In frozen gases, a Griffiths phase replaces the traditional phase transition.

The extended Langevin equation improves modeling of Rydberg atom density.

Abstract

The spread of excitations by Rydberg facilitation bears many similarities to epidemics. Such systems can be modeled with Monte-Carlo simulations of classical rate equations to great accuracy as a result of high dephasing. In this paper, we analyze the dynamics of a Rydberg many-body system in the facilitation regime in the limits of high and low temperatures. While in the high-temperature limit a homogeneous mean-field behaviour is recovered, characteristic effects of heterogeneity can be seen in a frozen gas. At large temperatures the system displays an absorbing-state phase transition and, in the presence of an additional loss channel, self-organized criticality. In a frozen or low-temperature gas, excitations are constrained to a network resembling an Erd\"os-Renyi graph. We show that the absorbing-state phase transition is replaced with an extended Griffiths phase, which we accurately describe by a susceptible-infected-susceptible model on the Erd\"os-Renyi network taking into account Rydberg blockade. Furthermore, we expand upon an existing macroscopic Langevin equation to more accurately describe the density of Rydberg atoms in the frozen and finite temperature regimes.

Equations

$$\frac{\text{d}}{\text{d}t}\hat{\rho }=i[\hat{\rho },\hat{\mathcal{H}}]+\sum _l{\hat{L}}_l\hat{\rho }{\hat{L}}_l^{†}-\frac{1}{2}\{{\hat{L}}_l^{†}{\hat{L}}_l,\hat{\rho }\}.$$

$$\hat{\mathcal{H}}$$

$$\displaystyle\Gamma_{\mathrm{f}}(\Sigma)=\frac{2\Omega^{2}\gamma_{\perp}}{\gamma_{\perp}^{2}+\Delta^{2}\big{(}\sum_{\begin{subarray}{c}j\neq i\\ j\in\Sigma\end{subarray}}\frac{r_{\text{f}}^{6}}{r_{ij}^{6}}-1\big{)}^{2}},$$

$$\frac{\text{d}}{\text{d}t}{P}_{\mathrm{r}}^{(i)}$$

$$\frac{\text{d}}{\text{d}t}{P}_{\mathrm{g}}^{(i)}$$

$$\tau =\frac{2{\Omega }^2{\gamma }_{\perp }}{{\gamma }_{\perp }^2+{\Delta }^2}.$$

$$\delta r_{\text{f}}=\frac{{\gamma }_{\perp }}{2\Delta }{r}_{\text{f}}.$$

$$\kappa ={\Gamma }_{\mathrm{f}}{V}_s.$$

$$\rho (\vec{r},t)=\frac{1}{\Delta V}\sum _{i:{\vec{r}}_i\in \Delta V}⟨{\hat{\sigma }}_i^{rr}⟩,$$

$$\displaystyle n(\vec{r},t)=\frac{1}{\Delta V}\sum_{i:\vec{r}_{i}\in\Delta V}\Bigl{(}\langle\hat{\sigma}^{rr}_{i}\rangle+\langle\hat{\sigma}^{gg}_{i}\rangle\Bigr{)}.$$

$$\frac{\text{d}}{\text{d}t}\rho$$

$$n$$

$$n_{\text{crit}}=\frac{\gamma }{\kappa },$$

$$P(k)=\frac{(n{V}_s{)}^k}{k!}\exp (-n{V}_s).$$

$$n_{\text{perc}}=\frac{1}{V_s}.$$

$$P(s)$$

$$={\text{e}}^{-n{V}_s}(1-{\text{e}}^{-n{V}_s}{)}^{s-1}$$

$$\propto {\text{e}}^{-cs},$$

$$\overline{v}\gg \gamma n^{-1/3},$$

$$\overline{v}\ll {\Gamma }_{\mathrm{f}}\delta r_{\text{f}},$$

$$\lambda$$

$$\mu$$

$${\lambda }_c^{(1)}=\mu$$

$${\lambda }_c^{(2)}=\mu \frac{⟨k⟩}{⟨k⟩-1},$$

$$n_{\text{crit}}^{(2)}$$

$$\equiv n_{\text{perc}}+n_{\text{crit}}^{(1)}.$$

$$\displaystyle\tau(s)\propto\sqrt{\frac{2\pi}{s}}\frac{\lambda}{(\lambda-1)^{2}}\text{exp}\Big{\{}s\Big{(}\text{ln}(\lambda)-1+\frac{1}{\lambda}\Big{)}\Big{\}}.$$

$$\rho (t)=\int \text{d}ssP(s){\text{e}}^{-t/\tau (s)}.$$

$$\rho \sim t^{-c/a},$$

$$\nu \equiv -\frac{c}{a}=-\frac{\ln \left(1-e^{-n{V}_s}\right)}{\ln (\lambda )-1+{\lambda }^{-1}},$$

$$p_{\text{exc}}$$

Full text

I Introduction

Rydberg atoms have gained a lot of attention in recent years due to their strong interactions over large distances [1]. This, paired with their long lifetimes in the order of milliseconds, creates a platform to explore quantum many-body physics of strongly interacting spin systems [2, 3, 4, 5, 6, 7, 8, 9] and to implement key elements for quantum information processing [10, 11, 12, 13, 14]. Moreover, optically driven Rydberg atoms (see Fig. 1a) can be used to investigate many-body dynamics of spin systems in inherently dissipative environments [15, 16, 17, 18, 19], as the laser excitation into high lying Rydberg states is often accompanied by strong dephasing. The latter includes important dynamical phenomena such as an absorbing-state phase transition (see Fig. 1b), one of the simplest classical non-equilibrium phase transitions displaying critical behavior and universality [20, 21].

Absorbing-state phase transitions are of general interest as they occur in many phenomena outside of physics such as population dynamics, epidemics, or the spreading of information in social media [22, 23, 24, 25]. Systems with this phase transition are believed fall into the universality class of directed percolation (DP) [20]. The unambiguous experimental observation of DP universal behavior is however challenging and has only been achieved in few systems in recent years [26, 27, 28, 29, 30, 31]. More recently, experimental signatures of such a transition have been reported in optically driven Rydberg gases [32].

In Rydberg systems, dissipation can give rise to another important dynamical phenomenon: self-organized criticality (SOC) [33, 34], which is believed to be a cause for the abundance of scale-invariance in nature [35, 36, 37, 38]. An SOC system dynamically evolves to the critical point of a phase transition by itself due to dissipation and without the need for parameter fine tuning (see Fig. 1c). Since the dissipation is strongly reduced once the critical point is reached, further evolution into the absorbing phase happens on much longer time scales. Recent experiments on Rydberg facilitation have shown some evidence of SOC through the use of ionization, or a decay into an auxiliary inert/”dead” state as a loss mechanism (see Fig. 1a) [39] (see [40] for related experiments).

However, the directed percolation transition is known to be susceptible to disorder [41] and more recent experiments on Rydberg facilitation in a trapped ultra-cold gas of atoms gave indications for an emergent heterogeneity in the system [42]. In such a heterogeneous system, the critical point of the absorbing-state phase transition is replaced by an intermediate extended Griffiths phase. Griffiths phases are characterized by generic scale invariance and the lack of universal behavior. This is in contrast to an absorbing state phase transition where scale invariance is only expected at the critical point. As a result (e.g. in the Rydberg gas), one expects a power law decay in active density over time with continuously varying exponents depending on the driving strength [43].

In [42], it was experimentally shown that a Rydberg system in the facilitation regime produces signatures of such a Griffiths phase for short times compared to the lifetime of the Rydberg state. A power-law decay in Rydberg density over time was observed with the decay exponents varying with driving strength and a phenomenological susceptible-infected-susceptible (SIS) network model was put forward to describe the observations. The model included a fitting function for the node weights of the network depending on the excitation rate $\kappa$. The interpretation being, that in the network model heterogeneity originates from a velocity selective excitation mechanism, where only atoms with relative velocities smaller than the Landau-Zener velocity $v_{\textrm{LZ}}(\kappa)$ could participate in facilitation dynamics. Above this velocity all further excitations are exponentially suppressed.

In the present paper, we present experimental indications for generic scale-invariance and strong theoretical indications for a Griffiths phase in a Rydberg facilitation gas by Monte-Carlo simulations.

In the experiment we continuously monitor the number of Rydberg excitations in a trapped ultra-cold gas of 87Rb atoms. We show that the size distribution of the Rydberg excitation number follows a power-law distribution, i.e. shows a scale free behavior, over an extended parameter regime, which is a key characteristic of a Griffiths phase.

In order to understand and to quantitatively describe the emergence of the Griffiths phase, we theoretically analyze two limiting cases: (i) a frozen gas and (ii) a gas with high temperature. While we recover a direct absorbing-state phase transition in the high-temperature limit with no signs of a velocity induced heterogeneity, we can identify a Griffiths phase in the frozen gas limit as a result of the finite paths along which facilitated excitations can spread. We give a quantitative analysis of the factors contributing to the emergence of a Griffiths phase and provide an estimate for the characteristic exponents of the power-law decay of Rydberg activity in this phase.

The facilitation of Rydberg excitations in a gas of optically driven atoms can be microscopically described by a Lindblad master equation [44] for the density matrix $\hat{\rho}$, which takes the form

[TABLE]

Here, the atom-light interaction Hamiltonian $\hat{\mathcal{H}}$ is given by

[TABLE]

where $\hat{\sigma}_{j}^{\mu\nu}=|\mu\rangle_{jj}\langle\nu|$ is the transition operator between states $|\nu\rangle$ and $|\mu\rangle$ of the $j$th atom. The strength of the laser driving shifted from the ground-Rydberg resonance frequency by the detuning $\Delta$ is described by the Rabi-frequency $\Omega$, and there is a van der Waals interaction proportional to $c_{6}/r_{ij}^{6}$, with $r_{ij}=|\vec{r}_{i}-\vec{r}_{j}|$ being the distance between atoms $i$ and $j$. Dissipative processes are taken into account by the Lindblad jump operators ${\hat{L}_{1}^{(i)}=\sqrt{(1-b)\gamma}\hat{\sigma}^{gr}_{i}}$, ${\hat{L}_{2}^{(i)}=\sqrt{b\gamma}\hat{\sigma}_{i}^{r0}}$ describing spontaneous decay of the Rydberg state into the ground state $\ket{g}$ and the inert state $\ket{0}$, with the branching parameter $b$. Finally, dephasing, attributed to laser phase noise and Doppler broadening [39] as well as the spread of the atomic wave packet over the van-der-Waals potential [45], is described by ${\hat{L}_{\perp}^{(i)}=\sqrt{\gamma_{\perp}}\hat{\sigma}^{rr}_{i}}$.

The strong van der Waals interaction of a Rydberg atom shifts energy levels of surrounding atoms significantly up to distances of multiple $\text{\,}\mathrm{\SIUnitSymbolMicro m}$. When the atoms are resonantly coupled to a laser field, this will block further excitations into Rydberg states from occurring for all atoms within a finite distance, a phenomenon known as Rydberg blockade [46]. On the other hand, if the laser excitation is strongly detuned, the excitation of isolated atoms is suppressed while atoms close to the facilitation distance $r_{\text{f}}\equiv\sqrt[6]{\frac{c_{6}}{\Delta}}$ are shifted into resonance (Fig. 1d) and are excited with a greatly increased rate. This process, termed Rydberg facilitation, leads to a cascade of excitations quickly spreading through the system following a single (off-resonant) excitation [47, 48]. It is important to note that Rydberg blockade still occurs in this regime. The excitation of atoms with distances ${r<r_{\text{f}}}$ is greatly suppressed (red zone in Fig. 1d).

II Experimental observation of scale-free behavior in a driven Rydberg gas

To experimentally test scale-invariance, we investigate the excitation density in a trapped gas of 87Rb atoms. To this end, we prepare a sample containing $150\text{\times}{10}^{3}$ atoms at a temperature of $1\text{\,}\mathrm{\SIUnitSymbolMicro K}$ in a crossed optical dipole trap. The sample has a density on the order of 10^{12}/$${\mathrm{cm}}^{3}. From the $5\mathrm{S}_{1/2}$ ground state, a UV laser at $297\text{\,}\mathrm{nm}$ continuously couples to the $40\mathrm{P}_{3/2}$ Rydberg state with a detuning of +$40\text{\,}\mathrm{MHz}$ and a resonant Rabi frequency of $2\pi\times$100\text{,}\mathrm{kHz}$$. The temperature of the gas corresponds to a most probable speed $\hat{v}=0.7\,r_{\text{f}}\gamma$ with the facilitation radius $r_{\text{f}}$ and decay rate $\gamma$.

Atoms in the $40\mathrm{P}_{3/2}$ state are ionized because of multiple intrinsic processes [49, 50], which we use to continuously monitor the excitation number. To this end we guide the resulting ions to a detector using a small electric field. This yields a time-resolved signal proportional to the number of Rydberg excitations in the sample (Fig. 2a).

At the beginning of the continuous laser exposure, which lasts $100\text{\,}\mathrm{ms}$, there are no excitations in the sample. As soon as the first off-resonant excitation is created, activity spreads through the system via facilitation, setting it up in the active phase. Due to the continuous atom loss caused by the ionization of excited atoms, the sample density decreases, reducing the effective driving strength. The sample thus approaches the phase transition.

We divide the ion signal in segments of $5\text{\,}\mathrm{ms}$ to account for the temporally varying effective driving. For each of these segments, we analyze the ion count distribution in $10\text{\,}\mathrm{\SIUnitSymbolMicro s}$ bins and average over $1000$ experimental runs. After about $10\text{\,}\mathrm{ms}$ the average activity has dropped more than an order of magnitude compared to its maximum value, while in individual runs it is dominated by avalanches. Therefore, we assume that at this time the sample is leaving the active phase.

Our measurement data shows persistent power-law behavior in the distribution of avalanche sizes over a wide range of densities (Fig. 2b). Power laws are a clear signature of scale invariance, which is expected only at the critical point of an absorbing-state phase transition or in a Griffiths phase characterizing a heterogeneous system. The extracted exponent of the power-law distribution is not fixed but varies with the density, strongly indicating non-universal behavior. While these observations are not an experimental proof of heterogeneity we use them as motivation to theoretically investigate possible origins of heterogeneity and a related Griffiths phase in the system.

III Microscopic Model of Rydberg facilitation

After having shown indications of scale-invariant behavior in the Rydberg facilitation gas, however with varying exponents in the experiments, we now turn to a theoretical modelling of the microscopic dynamics.

In the limit of large dephasing, the dynamics of a many-body Rydberg gas are effectively governed by classical rate equations [51]. As such, we will simulate a gas of atoms governed by (1) using classical Monte-Carlo simulations of a set of rate equations derived from (1) in the limit of large dephasing. After adiabatic elimination of coherences, eq. (1) reduce to classical rate equations between ground, Rydberg, and inert states (see Fig. 1a), with the stimulated rate $\Gamma_{\mathrm{f}}(\Sigma)$ given as

[TABLE]

where $\Sigma$ is the set of indices of Rydberg-excited atoms. In order to ensure numerical stability in the simulation, the singularity of the potential in eq. (3) is truncated at a cutoff value.

Using eq. (3) we can formulate a set of classical rate equations for the probability of the $i$th atom being in the Rydberg state $P_{\mathrm{r}}^{(i)}$ or the ground state $P_{\mathrm{g}}^{(i)}$ as

[TABLE]

If no other Rydberg atom exists in the gas or their distance is much larger than $r_{\text{f}}$, $\Gamma_{\mathrm{f}}(\Sigma)$ reduces to the off-resonant excitation rate of an isolated atom

[TABLE]

As a result of the broadening of the ground-Rydberg transition, given by the dephasing rate $\gamma_{\perp}$, facilitation can occur in a smeared out region around the facilitation distance $r_{\text{f}}$, given by

[TABLE]

Therefore, each Rydberg atom spans a facilitation shell around it at the radius $r_{\text{f}}$ and with the width ${\delta r_{\text{f}}}$ (white disks in Fig. 1d). Inside this shell, the stimulated rate takes its maximal value ${\Gamma_{\mathrm{f}}=\frac{2\Omega^{2}}{\gamma_{\perp}}}$ referred to as the facilitation rate. Relevant for later mappings to epidemic models is this rate integrated over the volume $V_{s}$ of the facilitation shell given by

[TABLE]

The relevant quantities of interest here are the coarse grained Rydberg density (in a small volume $\Delta V$)

[TABLE]

and the total active density of ground-state and Rydberg atoms

[TABLE]

In the following, $n$ will be referred to as the total density of the gas, for simplicity. As atoms in the state $\ket{0}$ do not participate in the dynamics of the system (see Fig. 1a), a decay into this state corresponds to a reduction of the total density, i.e. atom loss.

With this, $n\kappa$ corresponds to the rate with which excitations spread through the cloud.

The gas is simulated in a cube with size $L^{3}$ and periodic boundary conditions, typically ${L=7\,r_{\text{f}}}$. Atom positions are chosen randomly and velocities are sampled from the Maxwell-Boltzmann distribution with the temperature parameter $\hat{v}$, corresponding to the most probable atom velocity in the gas. After choosing a fixed time-step (${dt=1/400\;\gamma}$), the time evolution of the system is given by a fixed time step Monte-Carlo (ftsMC) algorithm [52]. We choose an ftsMC algorithm as opposed to a kinetic Monte-Carlo algorithm [53] as atomic movement, paired with long-range interactions leads to quickly changing transitional rates in the system.

In Ref. [39] Langevin equations have been derived to macroscopically describe the density of Rydberg atoms $\rho$ and the total density $n$ in the system. As is shown in [39], the homogeneous mean-field solution, in which diffusion terms are neglected, is sufficient to model the system. These equations then take the form

[TABLE]

with the off-resonant excitation rate $\tau$, and a noise term $\xi$. The parameter $b$ characterizes the percentage of Rydberg atoms which spontaneously decay into the dead state $\ket{0}$ (see Fig. 1a). As mentioned above, atoms that decay into this state are effectively removed from the system.

Assuming a gas with a heterogeneous density, diffusion results in a stabilization of the critical point over long times. For details pertaining to this see [54].

In the absence of decay into $|0\rangle$, i.e. for ${b=0}$, and in the absence of an off-resonant excitation, i.e. $\tau=0$, the dynamics described by eqs. (10), feature an absorbing-state phase transition at the critical atom density

[TABLE]

when the facilitation rate is fixed, or alternatively, at the critical facilitation rate $\kappa_{\textrm{crit}}=\gamma/n_{0}$ for fixed density. Below the critical point, any initially existing excitations in the system will eventually decay and the steady state of the system is one where all atoms are in the ground state (absorbing phase). Above the critical point, any arbitrarily small number of excitations initially present in the system will facilitate further excitations cascading through the system until a steady state with finite excitation density ${\rho(t\to\infty)>0}$ is reached (active phase).

Off-resonant excitations, with the rate $\tau$, will seed an excitation cascade in the active phase; whereas, in the absorbing phase, they cause fluctuations in the excitation number. As a result, the true absorbing state ${\rho=0}$ can only be approximately reached experimentally through a large separation of the off-resonant and facilitation time-scales, suppressing off-resonant excitations on the experimentally relevant facilitation time-scales.

Finally, the (slow) decay into a dead state $|0\rangle$ with rate $b\gamma$ is responsible for the self-organized approach to the critical point when starting in the active phase as indicated in Fig. 1c. Starting at an initial density $n_{0}$ above the critical value $n_{\textrm{crit}}$, i.e. in the active phase, the large number of atoms in the Rydberg state causes a fast loss of atoms into the dead state. As a consequence, the total density of atoms $n$ effectively participating in the facilitation process, i.e. atoms in states $|g\rangle$ and $|r\rangle$ decreases quickly and approaches the critical value. This loss continues at the critical density and drives the system further into the absorbing state. However, this happens on a much slower time-scale, as fewer Rydberg excitations are present at the critical point.

In Fig. 3 we have plotted the time evolution of the total density $n$, initially ten times higher than the critical density $n_{\textrm{crit}}$, and the Rydberg density $\rho$ for a frozen gas as well as a high-temperature gas with otherwise identical conditions, obtained from Monte-Carlo simulations. Here, all atoms in the system are initially in the ground state until one atom is off-resonantly excited to the Rydberg state. For comparison we also show the solution of the mean-field Langevin equations (10), which capture the long-time SOC dynamics of the high-temperature gas, but fail to describe the frozen gas outside of very short times (see Figs. 3b and c). The discrepancy in the peak values of $\rho$ can be attributed to Rydberg blockade, which truncates the maximum number of Rydberg excitations simultaneously present in the gas.

Qualitatively, the Rydberg density in the frozen gas displays a similar time dynamic to that of the high temperature gas, albeit with substantial quantitative differences in the long-time limit. We will show that in the low temperature regime of the Rydberg gas the absorbing state phase transition is replaced with an extended Griffiths phase, whose characteristic features become visible when off-resonant excitations and decay into state $\ket{0}$ are negligible.

It is important to note that for ${b>0}$ the decay into $\ket{0}$ dominates the dynamics at times ${t>1/b\gamma}$. Thus in order to experimentally observe a Griffiths phase by monitoring the long-time dynamics with this system, ionization and loss of Rydberg atoms (manifested in the parameter $b$) must be reduced as much as possible.

In Ref. [42] it was argued that a Rydberg atom moving at an average velocity larger than the Landau-Zener velocity ${v_{\textrm{LZ}}=2\pi^{2}\Omega^{2}r_{\text{f}}/(3\Delta)}$, effectively decouples from the excitation cascade. As a result, it was argued that this system features an emerging heterogeneity at high temperatures. Considering the two limiting cases of a frozen gas and a high temperature gas we argue that the Griffiths phase, which originates from spatial inhomogeneity, disappears when the atoms average velocity is increased above a certain limit, resulting in a direct absorbing-state phase transition.

A quantitative discussion of the crossover between a frozen system with an extended Griffiths phase and a high-temperature gas with a direct absorbing-state phase transition is beyond the scope of the present paper and is subject to future work. Instead we will focus on the quantitative description of the facilitation dynamics in a low-temperature or frozen gas.

IV Network Structure of Facilitation Paths in a Frozen Gas

The emergence of a Griffiths phase results from facilitation events being constrained to a network structure. In the limit of a frozen gas, atoms have random but fixed positions. If we regard the system at the time scale of facilitated excitations, off-resonant excitations can be neglected. Therefore, the dynamics are described by the facilitated spreading of Rydberg excitations, which is only possible if atomic distances are approximately $r_{\text{f}}$. As a result, we can regard the structure of atom positions and the paths along which excitations can spread as a random graph with edges where atoms have the distance ${r\in[r_{\text{f}}-\frac{\delta r_{\text{f}}}{2},r_{\text{f}}+\frac{\delta r_{\text{f}}}{2}]}$.

Assuming a uniform distribution of atom positions in the gas, the probability that a randomly selected atom has $k$ atoms in its facilitation shell (see Fig. 1d), meaning the atom is of degree $k$, is given by the Poissonian distribution

[TABLE]

As the degree distribution is Poissonian we can map this problem to a random Erdős–Rényi (ER) network [55]. In contrast, the network structure of atoms trapped by an optical lattice or tweezer array would be given by a regular lattice network.

Of particular interest in random graph theory is the question if a system percolates. In a percolating system, the probability $p$ that a bond between two randomly selected atoms exists is high enough, such that a path exists which runs through the entire system, i.e. there almost surely exists a single cluster (i.e. a single connected set of vertices) with its size in the order of the system size. If, however, the connectivity is below a critical threshold for bond connectivity $p<p_{c}$, the system is composed of many small, disconnected clusters [55, 56]. For $p=p_{c}$ the percolation transition occurs. A 2D network with ${p=p_{c}}$ from Monte-Carlo sampling is illustrated in Fig. 4.

If $N$ is the number of atoms and $s_{1}(N)$ is the size of the largest connected cluster (LCC), then the system percolates if ${\lim_{N\to\infty}{s_{1}(N)}/{N}>0}$. For an ER network the percolation transition occurs when the average network degree is $\langle k\rangle=1$ [56, 57]. Using eq. (12), the density at which the percolation transition occurs is therefore

[TABLE]

This density is a factor $\Gamma_{\mathrm{f}}/\gamma$ larger than the critical density $n_{\text{crit}}$ of the absorbing state phase transition, given by eq. (11). We can verify that eq. (13) corresponds to the correct percolation density by calculating the size of the LCC $s_{1}$ depending on the density of the gas (Fig. 4). In the thermodynamic limit, ${s_{1}/N=0}$ for all densities ${n<n_{\textrm{perc}}}$. As numeric simulations are restricted to a finite system size however, we instead consider the percolation transition to occur when $s_{1}$ grows faster than linear with the density $n$ (the black dashed line in Fig. 4 corresponds to linear growth).

Of relevance for the Griffiths phase is the size distribution of clusters in the network. Using Monte-Carlo simulations we can verify that the lengths of clusters follow a geometric distribution ${P(s)\sim\text{e}^{-cs}}$ under the assumption that clusters are made of linear chains of $s$ atoms. This assumption holds true for small cluster sizes and an average network degree ${\langle k\rangle\ll 1}$.

We can then approximate the decay constant $c$, with $p_{0}$ being the probability of an atom having the degree ${k=0}$, as

[TABLE]

with ${c=-\text{ln}(1-\text{e}^{-nV_{s}})}$. In Fig. 5 a comparison between cluster sizes in Monte-Carlo simulations and eqs. (14) is shown. The agreement is very good for small densities as almost all clusters in the gas are composed of linear chains. As the density of the gas increases, the probability that at least one atom in the cluster has more than two connections, i.e. $k\geq 3$, increases. While the distribution remains exponential, the probability for larger clusters to exist in the system greatly increases compared to the prediction by eqs. (14).

V Epidemic Dynamics on the Network

It is known that Rydberg systems in the facilitation regime bear close similarities to epidemics [58, 42]. In the following, we will systematically analyze the Rydberg facilitation dynamics on the random network formed by atoms within their respective facilitation shells. For this, we will map the dynamics to the Susceptible-Infected-Susceptible (SIS) epidemic model [59, 60, 61]. We will (i) disregard the decay of Rydberg atoms into the dead state $\ket{0}$ (parameter ${b=0}$, see Fig. 1a), reducing the dynamics of each atom to a two-level system. Additionally, we will (ii) neglect off-resonant excitations by setting ${\tau=0}$, meaning excitations can only be created by means of facilitation. We will refer to simulations carried out with these two constraints as the SIS approximation.

One major difference between our Rydberg system in the SIS approximation and a classical SIS system remains with Rydberg blockade. Atoms excited to the Rydberg state do not only facilitate the spread of excitations, they can also block the spreading in adjacent clusters. This will be analyzed more systematically later.

The network structure of cluster of atoms in facilitation distance of each other strongly depends on the temperature of the gas. If the RMS average relative velocity $\overline{v}$ is large, such that each excited Rydberg atom meets many ground state atoms during its lifetime $\gamma^{-1}$, i.e. if in a 3D gas

[TABLE]

any network structure is effectively washed out and the system is homogeneous. Close to the critical point of the absorbing-state phase transition, the above condition is equivalent to $\overline{v}\gg\gamma\,r_{f}$. If, on the other hand, the average velocity of atoms is very small, such that during a facilitation time $\Gamma_{\mathrm{f}}^{-1}$ they do not move out of the facilitation shell, i.e. if

[TABLE]

the atoms form a finitely connected network. We will now discuss these two limits.

V.1 High temperature limit

In a high-temperature gas with RMS average relative velocity $\overline{v}\gg\gamma n^{1/3}$, we can map the system to the Susceptible-Infected-Susceptible (SIS) epidemic model [59, 60, 61]. The SIS model is characterized by the infection and recovery rates, $\lambda$ and $\mu$ respectively, which for the Rydberg gas read

[TABLE]

The SIS model predicts an active/absorbing phase transition when

[TABLE]

where excitation spread equals spontaneous decay. This corresponds to the critical density (11) of the absorbing state phase transition discussed before.

In Fig. 6a, Monte-Carlo simulations of the Rydberg system with the SIS approximation and ${\rho(t=0)=n}$ are shown for the high temperature gas for different values of $n$ and a fixed facilitation rate $\Gamma_{\mathrm{f}}$. We note that as shown in Appendix A the excitation probabilities following from Monte-Carlo simulations of rate equations and those from full coherent density matrix simulations agree, showing that the rate equation approach remains valid also in the high-temperature limit. One recognizes that an active/absorbing-state phase transition occurs for $\lambda=\lambda_{c}^{(1)}$, with the Rydberg density either exponentially decaying at the time-scale $\mu$ (for $\lambda<\lambda_{c}^{(1)}$), or decaying to a steady-state active density (for $\lambda>\lambda_{c}^{(1)}$). At the critical density (green curve in Fig. 6a) the system should decay with $\rho\sim t^{-1}$ [62], however, this decay is truncated by an exponential decay due to finite system size.

V.2 Frozen gas limit

In the limit of an effectively frozen gas the atoms that can participate at the facilitation process form a network. The dynamics of an SIS epidemic strongly depend on the structure of this underlying network. For example, in the case of a heterogeneous but scale-free network, i.e. ${P(k)\sim k^{-\nu}}$, the absorbing phase can disappear altogether, leaving the system in an endemic phase regardless of the infection rate [63, 64, 65].

For the case of a heterogeneous ER network, which describes the frozen gas of atoms, an active phase can only occur if the network is above the percolation threshold (i.e. ${\langle k\rangle>1}$). However, for a finitely connected (but percolating) ER network, the threshold for the active phase is modified since activity occurs in localized regions and thus the effective infection rate is reduced. One finds [62]

[TABLE]

with $\langle k\rangle$ being the average degree of the network. For ${\langle k\rangle\to\infty}$ the threshold given by eq. (18) is recovered. For a fixed facilitation rate and facilitation volume this threshold can be expressed in terms of a critical density of atoms using $\langle k\rangle=nV_{s}$

[TABLE]

If the network is below the percolation threshold, the finite size of clusters truncates the spread of activity through the system. Therefore, the network cannot support an active phase, and instead, a Griffiths phase emerges above the critical infection rate $\lambda_{c}^{(1)}$ [62]. One of the most distinguishing characteristics of a Griffiths phase is the presence of rare regions with above average activity which lead to a slow, algebraic decay of excitations [43].

In the non-percolating network (i.e. ${\langle k\rangle<1}$), for ${\lambda\leq\mu}$ decay dominates, leading to very short times until all activity disappears in clusters as excitations cannot sustain themselves. If however ${\lambda>\lambda_{c}^{(1)}=\mu}$, the time until activity disappears in clusters increases exponentially with cluster size $s$ and is given by [66]

[TABLE]

In the following we will refer to $\tau(s)$ as the extinction time of activity in a cluster. As a result of the convolution of exponentially rare cluster sizes ${P(s)\sim\text{e}^{-cs}}$ and a cluster lifetime increasing exponentially with cluster size ${\tau(s)=\text{e}^{as}}$, the activity in the Griffiths phase decays with a power law:

[TABLE]

Using eqs. (14) and (21), the integral in (22) can be approximated with Laplace’s method and results in an algebraic decay

[TABLE]

with the coefficient $a$ given by (21) as ${a=\text{ln}(\lambda)-1+\frac{1}{\lambda}}$. If the network is above the percolation threshold, i.e. if $\langle k\rangle\geq 1$, but the driving strength is below the critical value for the active phase $\lambda_{c}^{(2)}$ the decay of activity is expected to follow a stretched exponential. A qualitative phase diagram of the facilitation dynamics in the frozen Rydberg gas is shown in Fig. 7.

Fig. 6b shows the results of Monte Carlo simulations for a frozen gas in the SIS approximation for the same parameters and color code as in the high-temperature case of Fig. 6a. For ${n<n_{\textrm{crit}}}$ all initial excitations decay exponentially (curves 1 and 2 from left to right), corresponding to the absorbing phase. The behavior changes at and above the critical point but below the percolation threshold ${n_{\textrm{crit}}\leq n<n_{\textrm{perc}}}$ (curves 3-7). Here, the system is in an extended Griffiths phase with a power-law decay with varying exponents. Above the percolation threshold but below the threshold of the active phase $n_{\textrm{perc}}\leq n<n_{\textrm{crit}}^{(2)}$ the decay is expected to become a streched exponential [62], which we cannot resolve however in our simulations due to the very long time scale of this decay. Finally, for ${n\geq n_{\textrm{crit}}^{(2)}}$ the system enters the active phase where excitations simply decay to a steady state.

In the following, we want to give a quantitative estimate for the power law decay coefficient in the Griffiths phase based on the SIS model and compare them with those from the Monte-Carlo simulations. In contrast to the standard SIS model, a Rydberg system features Rydberg blockade and facilitated de-excitation, making it unclear if analytic predictions from an SIS model would be accurate. To check this, we compare the extinction time of activity in clusters in a linear excitation chain, given by eq. (21), using the spreading rate ${\lambda=\Gamma_{\mathrm{f}}V_{s}\cdot 1r_{\text{f}}^{-3}}$, with Monte-Carlo simulations of the SIS approximation in Fig. 8. For this, we simulate a 1D cluster of length $s$ where each atom is initially in the Rydberg state and measure the average time until all atoms are decayed. Here, we assume the above mentioned SIS approximation (no decay to $\ket{0}$ and no off-resonant excitations). One recognizes that eq. (21) gives a good approximation of the extinction time.

Using (14) and (21) we can approximate the power-law exponent $\nu$ in the Griffiths phase dependent on the density and internal rates. We receive

[TABLE]

with $\lambda=4\pi\Gamma_{\mathrm{f}}\frac{\delta r_{\text{f}}}{r_{\text{f}}}$. The comparison with exponents fitted from the power law decay of Rydberg density in Monte-Carlo simulations of the frozen gas under the SIS approximation (seen in Fig. 6) can be seen in Fig. 9b. Our very rough approximation of the Griffth-phase decay exponents qualitatively fits with Monte-Carlo data.

A fundamental difference between Rydberg facilitation and classical SIS activity spreading is Rydberg blockade. Considering the frozen gas limit, two effects arise from Rydberg blockade: first, if an atom is surrounded by two Rydberg atoms in the facilitation distance, i.e. the atom is in the middle of a cluster, then it cannot be facilitated, as it receives twice the dipole shift and is pushed out of resonance again. If this atom decays or is in the ground state at the beginning, it cannot be excited resulting in a hole splitting the cluster [19]. Additionally, Rydberg atoms can block excitations from spreading through adjacent clusters. However, neither of these effects change the actual structure of the network, instead they effectively retard the timescale at which excitations spread. For a quantitative comparison we simulate the Rydberg gas and compare the decay of excitations in the SIS approximation with and without Rydberg blockade (Fig. 9).

As blockade allows fewer Rydberg atoms to be present in the system, the steady state Rydberg density of the active phase, and therefore the density at which the power law decay of the Griffiths phase begins is much lower. However, as seen in Fig. 9, the exponents of the power-law decay in the Griffiths phase show no qualitative change depending on the presence of Rydberg blockade in the system.

VI Conclusion

We studied the facilitation dynamics of Rydberg excitations in an ultra-cold gas of atoms. While in the homogeneous limit, the system is expected to show a phase transition between an absorbing phase and an active phase, and – in the presence of an additional loss channel from the Rydberg state – self-organized criticality (SOC). However, experiments with a gas of trapped 87Rb atoms at low temperatures show signs of scale invariant dynamics in an extended parameter regime, which is a signature of a Griffiths phase replacing the critical point of the absorbing-state phase transition.

To understand the emergence of scale invariance in the experiment in an extended parameter regime, we numerically simulated the many-body Rydberg gas in the facilitation regime through the use of Monte-Carlo simulations in the classical rate-equation approximation. We showed that the latter is well justified for the large dephasing characteristic for the experiment even for a high-temperature gas. Since a Griffiths phase originates from heterogeneity in the system, we numerically and theoretically analyzed two limiting cases: (i) a high-temperature gas and (ii) a frozen gas. While in the high-temperature limit, a homogeneous mean-field behavior is recovered, with a clear absorbing-state phase transition and SOC, the facilitation dynamics in a low-temperature or frozen gas is governed by the presence of a network structure of atoms that can participate in the excitation spread. Numerical simulations show characteristic power-law decay of Rydberg excitations in time if off-resonant excitations and atom losses are neglected.

We have shown that in the frozen gas the spread of excitations is constrained to a network resembling a random Erdős–Rényi (ER) graph. Increasing the density of atoms the ER network has a percolation transition from a fragmented phase, in which the maximum cluster size of connected atoms remains finite, to a phase where the size of the largest cluster scales with the size of the system. A theoretical explanation of the Rydberg facilitation dynamics observed in Monte-Carlo simulations can then be given by mapping to a susceptible-infected-susceptible (SIS) epidemic model on such an ER graph taking into account the effects of Rydberg blockade, which truncates the maximum Rydberg excitation density. An active phase of self-sustained Rydberg excitations is only possible above the percolation threshold. Below this threshold, an extended Griffiths phase emerges in the place of the (for homogeneous systems) expected absorbing-state phase transition. We showed that the modified SIS model quantitatively explains the observed power-law decay exponents as well as the overall dynamics of the Rydberg density.

While the limits of a high-temperature and a frozen gas are well captured with our model, it does not yet allow the study of the crossover between the two regimes. To this end, the Rydberg facilitation process needs to be mapped to a dynamical network, which is beyond the scope of the present work and will be the subject of future work. Furthermore, in order to quantitatively understand the power-law exponents in the number distribution of Rydberg atoms in a given time interval observed in the experiment, it is necessary to extend the microscopic simulations to much larger system sizes matching those used in the experiments. To this end different approaches, e.g. using machine-learning algorithms might be useful [67]. Finally, the interplay between coherent quantum dynamics and dynamical network structures in Rydberg facilitation under conditions where dephasing is much less dominant could give rise to very different dynamics [68, 69]. The latter requires, however, the development of new microscopic simulation techniques capable of incorporating quantum coherences in 3D Rydberg gases, at least in an approximate way [70].

Acknowledgement

The authors thank Johannes Otterbach for fruitful discussions. Financial support from the DFG through SFB TR 185, project number 277625399, is gratefully acknowledged.

Authors contributions

DB and MF developed the theoretical models, DB performed all numerical simulations and developed the mapping to the ER network with support by SO. MF, TN, and HO conceived the project. JB, PM, and DB performed the experiments guided by TN and HO. DB and MF wrote the initial version of the manuscript with support by JB. All authors discussed the results and contributed to the writing of the manuscript.

Appendix A Effects of relative motion between atoms

The rate equation approximation used for the Monte-Carlo simulations (e.g. eq. 3) is valid as long as the population dynamics are slow compared to the dephasing rate. In a frozen gas, the relevant time scales are solely determined by the internal dynamics of an individual atom for a given (fixed) configuration of Rydberg atoms in its vicinity. If however, the gas of atoms has a finite temperature, a ground state atom can fly in and out of the facilitation volume $V_{s}$ of a Rydberg atom, which can amount to a fast sweep of the detuning of the ground state atom. Thus there is an additional time scale given by the crossing time $\sim\delta r_{f}/v$.

To analyze the effects of atomic motion onto the facilitation process we consider the two-body problem of a ground state atom moving with velocity $v$ and with the impact parameter $d$ relative to a Rydberg atom (see inset of Fig. 10). For ${d>r_{\text{f}}}$ the ground state atom is not shifted into resonance and no facilitation occurs. For ${d\leq r_{\text{f}}}$ one has to distinguish two cases depending on the impact parameter: (i) ${d<r_{\text{f}}}$ the ground state atom flies through the facilitation shell twice (blue case in Fig. 10), and (ii) ${d\approx r_{\text{f}}}$ the ground state atom grazes the facilitation shell and is only briefly shifted into resonance (orange case in Fig. 10).

In case (i), using the excitation rate from eq. (3) as $\Gamma_{\uparrow}(t)$, we find the excitation probability after a single pass of the ground state atom through the facilitation shell as

[TABLE]

Note that this expression assumes a short passage time through the facilitation shell, so that the facilitated de-excitation process can be ignored. For longer passage times the excitation probability approachs the steady state value of ${1/2}$, as can be seen Fig. 10.

Linearizing the time dependent detuning $\Delta(t)$ for times $t_{i}<t<t_{f}$, while passing through the facilitation shell, we receive $\Delta(t)\approx\dot{\Delta}\times(t-t_{0})$ yielding

[TABLE]

where we have assumed that $|\Delta_{i,f}|=|\Delta(t_{i,f})|\gg\gamma_{\perp}$, which is exactly the same expression as given by the Landau-Zener formula.

If $p_{\text{exc}}$ is small, the asymptotic excitation probability after two passages is just ${p_{\text{exc}}\approx 1-\exp\{-4\pi\Omega^{2}/\dot{\Delta}\}}$. From this discussion we expect the rate equations to accurately describe the facilitation process even for large atom velocities as long as the impact parameter $d$ is different from $r_{f}\pm\delta r_{f}$.

In case (ii) however, i.e. for ”grazing incidence”, the Landau-Zener formula is no longer valid and there could be a difference between the rate-equation approximation and the solution of the full two-particle density matrix equations. This is indeed the case, as can be seen from Fig. 10, where we have plotted the asymptotic excitation probability of the ground state atom as function of relative velocity and impact parameter both from a simulation of the full density-matrix equations (dashed lines), the analytic Landau Zener formula (solid line), and by a Monte-Carlo simulation of the rate equation in the large-dephasing limit with time step ${dt=1/400\,\gamma}$ (dots). One recognizes perfect agreement except for large relative velocities and impact parameters close to the facilitation radius ${d\approx r_{\text{f}}}$, where the rate equations predict up to an order of magnitude higher excitation probabilities than the full simulation. Since ${\delta r_{\text{f}}\ll r_{\text{f}}}$ the contribution of these ”grazing-incidence” cases is negligibly small, allowing us to accurately describe high gas temperatures with a fixed time step Monte-Carlo algorithm.

Furthermore, at high temperatures (as can be seen from Fig. 10) the excitation probability above the Landau-Zener velocity ${v_{\textrm{LZ}}=2\pi^{2}\Omega^{2}r_{\text{f}}/(3\Delta)}$ indeed quickly drops and scales as $1/v$, the number of ground state atoms seen by a moving Rydberg atom in a given time increases linearly with its velocity $v$, too. This compensates the former effect, and thus does not lead to an emerging heterogeneity in phase space as argued in [42], as long as the Rydberg atom does not move out of the gas sample.