Controlling excitation avalanches in driven Rydberg gases

TL;DR

This paper proposes a control mechanism for avalanche dynamics in driven Rydberg gases, enabling exploration of self-organized criticality and related phenomena in cold atom systems.

Contribution

It introduces a time-dependent drive amplitude to control avalanche behavior, expanding the understanding of SOC in Rydberg ensembles.

Findings

Control over avalanche regimes achieved via drive modulation

Demonstrated dependence of dynamics on external scales and dimensionality

Potential for exploring extended SOC phenomenology in cold atom systems

Abstract

Recent experiments with strongly interacting, driven Rydberg ensembles have introduced a promising setup for the study of self-organized criticality (SOC) in cold atom systems. Based on this setup, we theoretically propose a control mechanism for the paradigmatic avalanche dynamics of SOC in terms of a time-dependent drive amplitude. This gives access to a variety of avalanche dominated, self-organization scenarios, prominently including self-organized criticality, as well as sub- and supercritical dynamics. We analyze the dependence of the dynamics on external scales and spatial dimensionality. It demonstrates the potential of driven Rydberg systems as a playground for the exploration of an extended SOC phenomenology and their relation to other common scenarios of SOC, such as e.g. in neural networks and on graphs.

Equations

$$\Omega \to {\Omega }_t={\Omega }_0(1+t\frac{\lambda }{2n_{c,0}}),$$

$${\partial }_t\hat{\rho }=i[\hat{\rho },H]+\sum _l{\mathcal{L}}_l\hat{\rho }$$

$$\displaystyle H=\sum_{l}\left\{\left(\sum_{l^{\prime}\neq l}\frac{C_{6}}{2|\vec{r}_{l}-\vec{r}_{l^{\prime}}|^{6}}\sigma^{rr}_{l^{\prime}}-\Delta\right)\sigma^{rr}_{l}+\frac{\Omega_{t}}{2}\left(\sigma^{rg}_{l}+\sigma^{gr}_{l}\right)\right\},\ \$$

$${\mathcal{L}}_l\hat{\rho }$$

$${\rho }_{\vec{x},t}$$

$$n_{\vec{x},t}$$

$$\displaystyle\partial_{t}\rho_{\vec{x},t}=D\nabla^{2}\rho_{\vec{x},t}+(\kappa_{t}\rho_{\vec{x},t}+\tau_{t})(n_{\vec{x},t}-2\rho_{\vec{x},t})-\Gamma\rho_{\vec{x,t}}+\xi_{\vec{x},t}.\ \ \$$

$$\overline{{\xi }_{\vec{x},t}{\xi }_{\vec{y},t^{\prime }}}=\delta (\vec{x}-\vec{y})\delta (t-t^{\prime })\left[(\tau +\kappa {\rho }_{\vec{x},t})n_{\vec{x},t}+2\Gamma {\rho }_{\vec{x},t}\right].$$

$${\partial }_t{n}_{\vec{x},t}=D_n{\nabla }^2{n}_{\vec{x},t}-{\gamma }_{\downarrow 0}{\rho }_{\vec{x},t}+{\eta }_{\vec{x},t}$$

$$0\stackrel{!}{=}{\partial }_t⟨{\sigma }_l^{r\alpha }⟩=\text{Tr}\left[{\sigma }_l^{r\alpha }\left(i[\hat{\rho },H]+\sum _l{\mathcal{L}}_l\hat{\rho }\right)\right].$$

$${\partial }_t{\sigma }_l^{gg}$$

$${\partial }_t{\sigma }_l^{rr}$$

$$\overline{({\xi }_l^r{)}^2}$$

$$\frac{{\Omega }_t^2\Gamma }{{\Gamma }^2+4(\Delta -{\sum }_{l^{\prime }\ne l}{V}_{l,l^{\prime }}^I{\sigma }_{l^{\prime }}^{rr}{)}^2}=\underset{={\tau }_t}{\underbrace{\frac{{\Omega }_t^2\Gamma }{{\Gamma }^2+4{\Delta }^2}}}+$$

$$\sum _{l^{\prime }\ne l}\left(\frac{{\Omega }_t^2\Gamma }{{\Gamma }^2+4(\Delta -V_{l,l^{\prime }}^I{)}^2}-{\tau }_t\right){\sigma }_{l^{\prime }}^{rr}+O({\sigma }_{l^{\prime }}^{rr}{\sigma }_{l^{\prime \prime }}^{rr}).$$

$${\partial }_t{\sigma }_l^{rr}$$

$${\partial }_t{\rho }_{\vec{x},t}=\sum _l\left(⟨{\partial }_t{\sigma }_l^{rr}⟩+⟨{\sigma }_l^{rr}⟩{\partial }_t{\vec{r}}_l\vec{\nabla }\right)\Theta (r_{\text{fac}}-|\vec{x}-\overrightarrow{r_l}|)$$

$${\partial }_t{n}_{\vec{x},t}$$

$${\partial }_t{\rho }_{\vec{x},t}$$

$${\mathcal{F}}_{\vec{x}}({\rho }_{\vec{z},t})={\mathcal{F}}_{\vec{x}}(1){\rho }_{\vec{x},t}+\frac{{\mathcal{F}}_{\vec{x}}({\vec{z}}^2)}{2}{\nabla }^2{\rho }_{\vec{x},t}+O({\nabla }^4{\rho }_{\vec{x},t}).$$

$${\partial }_t{n}_{\vec{x},t}$$

$${\partial }_t{\rho }_{\vec{x},t}$$

$$n_{c,t}=\frac{\Gamma }{{\kappa }_t}=n_{c,0}-{\frac{\partial {\kappa }_t}{\partial t}|}_{t=0}\frac{t\cdot n_{c,t}}{{\kappa }_t}=n_{c,0}-\lambda t,$$

$${\rho }_{\vec{x},t}=\frac{\lambda }{{\gamma }_{\downarrow 0}}\text{ and }{n}_{\vec{x},t}=n_{c,t}+\frac{2\lambda }{{\gamma }_{\downarrow 0}}+\frac{{\gamma }_{\downarrow 0}{\tau }_t}{\kappa \lambda n_{c,t}}.$$

$$\text{(i) }\frac{{\tau }_t}{\lambda }\to 0^{+}\text{ and (ii) }\frac{\lambda }{{\gamma }_{\downarrow 0}}\to 0^{+},$$

$${\tau }_t\ll \lambda \ll {\gamma }_{\downarrow 0}.$$

$$R_t\equiv n_0+\lambda t-{\int }_0^{t}d{t}^{\prime }{\gamma }_{\downarrow 0}⟨{\rho }_{\vec{x},t^{\prime }}{⟩}_V,$$

$${\partial }_t{\rho }_{\vec{x},t}$$

$$P(\rho ,\delta t)=\lambda e^{-\lambda \left({\rho }_0{e}^{\beta \delta t}+\rho \right)}{\left(\frac{\rho }{{\rho }_0{e}^{\beta \delta t}}\right)}^{\mu /2}{I}_{\mu }\left(2\lambda \sqrt{{\rho }_0\rho e^{\beta \delta t}}\right),$$

$$\rho =\Gamma [\mu +1+\text{Poisson}[\lambda {\rho }_0{e}^{\beta \delta t}]]/\lambda ,$$

$${\partial }_t{n}_{\vec{x},t}$$

$${\partial }_t\overrightarrow{p_l}=-\vec{\nabla }V({\vec{r}}_l)-\frac{1}{t_{\text{rel}}}{\vec{p}}_l.$$

$${\partial }_t{n}_{\vec{x},t}$$

Full text

I Introduction

Away from thermal equilibrium and in the absence of detailed balance, (quasi-) stationary states emerge from ordering principles different from the equipartition of energy. Outstanding amongst such out-of-equilibrium ordering mechanism is self-organized criticality. Introduced in the seminal paper of Bak, Tang and Wiesenfeld (BTW) Bak et al. (1987, 1988) to explain the emergence of flicker noise in electrical circuits, SOC has since then been observed in a variety of diverse, mainly large scale systems, ranging from earth quakes Sornette and Sornette (1989); Chen et al. (1991); Bak et al. (2002), forest fires  Schenk et al. (2000); Drossel and Schwabl (1992); Malamud et al. (1998); Turcotte (1999) and solar flares Lu and Hamilton (1991); Aschwanden et al. (2016) to vortex dynamics in superconductors Field et al. (1995); Altshuler and Johansen (2004) and turbulence Chapman and Nicol (2009). Only recently, SOC was recognized as a possible mechanism to establish optimal conditions for information spreading de Arcangelis et al. (2006); Shew et al. (2015); Hesse and Gross (2014); Marković and Gros (2014); Kinouchi and Copelli (2006); Levina et al. (2007).

The phenomenon of SOC can be described by simple means, by balancing dissipation and external drive, a many-body system is attracted, i.e., it self-organizes, towards a state with scale invariant correlations Lu (1995); Watkins et al. (2016); Dickman et al. (2000). In thermal equilibrium, scale invariance is associated with dynamics at a critical point signaling a continuous phase transition Zinn-Justin (1996). Compared to a fine tuned critical point, scale invariance due to SOC is believed to occur in an extended parameter regime, commonly enabled by a separation of time scales between drive and dissipation Vespignani and Zapperi (1997); Dickman et al. (2000); Watkins et al. (2016). While this makes SOC robust to changes in the external conditions, the interplay of interactions, drive and dissipation obscure its origin and only few microscopic models are found in the literature.

Apart from the sandpile model of BTW, manifestations of SOC in nature are mostly approached via phenomenological models Malamud et al. (1998); Rybarsch and Bornholdt (2014), either because the microscopic description is too complex or the elementary building blocks are unknown Rhodes and Anderson (1996); Shew et al. (2015); Aschwanden et al. (2018). Unfortunately, many realizations of SOC don’t match the energy conserving dynamics of BTW’s sandpile model. This makes both the microscopic understanding and, even more, the controllability of SOC extremely challenging Watkins et al. (2016); Dickman et al. (2000). This applies especially to neural network dynamics, which are entirely based on effective models Strogatz (2001); Barzel and Barabási (2013). Consequently, a setup exploring an extended SOC phenomenology on the one hand and featuring the knowledge and a large degree of controllability of its basic elements on the other hand represents a promising tool to study aspects of SOC in generic nonequilibrium settings.

Only recently a promising candidate has been introduced in an experiment with a gas of driven Rydberg atoms Helmrich et al. (2018); above a certain driving threshold, the atomic pseudo-spins self-organize towards a transient, scale invariant state, featuring common signatures of SOC Dickman et al. (2000); Aschwanden et al. (2018). Our work builds up on this basic setting for SOC in cold Rydberg gases.

We propose the implementation of a control mechanism for excitation avalanches in driven Rydberg ensembles and explore the corresponding many-body dynamics. We show how this gives access to an extended SOC phenomenology, including subcritical and supercritical avalanche dynamics. By adjusting the proposed mechanism to common control parameters such as the laser intensity and the detuning, one can access the paradigms of SOC: a scale invariant avalanche distribution Turcotte (1999) and a $\frac{1}{\omega}$-noise pattern Bak et al. (1987, 1988).

II Facilitated Rydberg dynamics

We consider the many-body dynamics in a gas of interacting Rydberg atoms Schauß et al. (2012); Günter et al. (2013); Gorshkov et al. (2013); Helmrich et al. (2018); Letscher et al. (2017); Thomas et al. (2018), which move freely inside a trap. Each Rydberg atom is modeled as an effectively three level system, consisting of a non-interacting ground state $|g\rangle$, a highly polarizable Rydberg state $|r\rangle$ with large principle quantum number $n\gg 1$ Gallagher (1984); Saffman et al. (2010); Löw et al. (2012) and an auxiliary, removed state $|0\rangle$. The latter is a container state representing a set of internal states that can be reached via dissipative decay but are otherwise decoupled from the $|g\rangle-|r\rangle$ sector Helmrich et al. (2018, 2018). Each atom obtains a label $l$ and a set of operators $\sigma^{ab}_{l}\equiv|a\rangle\langle b|_{l}$ acting on its internal states.

The ensemble is subject to a laser, coherently driving the $|g\rangle-|r\rangle$ transition with a Rabi frequency $\Omega$ and detuning from resonance $\Delta$. The highly excited Rydberg state is subject to dissipation originating from dephasing as well as spontaneous decay into both the ground state $|g\rangle$ and the removed state manifold $|0\rangle$ with effective rates labeled by $\gamma_{\text{de}},\gamma_{\downarrow g},\gamma_{\downarrow 0}$ Helmrich et al. (2018). Due to their polarizability, two atoms, labeled $l,l^{\prime}$, in the Rydberg state experience a mutual van-der-Waals repulsion. Its potential form is $V^{I}_{l,l^{\prime}}=C_{6}|\vec{r}_{l}-\vec{r}_{l^{\prime}}|^{-6}$, where $C_{6}$ is the van-der-Waals coefficient and $\vec{r}_{l},\vec{r}_{l^{\prime}}$ are the atomic positions Baluktsian et al. (2013)111The interaction might as well acquire a dipole-dipole form, $V\sim|\vec{r}_{l}-\vec{r}_{l^{\prime}}|^{-3}$, e.g., due to Foerster resonances Li et al. (2005). This does, however, not modify the structure of Eq. (7)..

As a simple but crucial innovation we consider here a time-dependent Rabi frequency

[TABLE]

with an initial frequency $\Omega_{0}$, a dimensionless density $n_{c,0}$, which we define later, and a ramp parameter $\lambda\ll n_{c,0}\Omega_{0}$. This corresponds to a slow, linear increase of the pump laser intensity $I_{t}\sim\Omega_{0}^{2}(1+\frac{\lambda t}{n_{c,0}}+O(\lambda^{2}t^{2})$. It gives rise to a continuously increasing excitation probability for the $|g\rangle\leftrightarrow|r\rangle$ transition, counteracting the decay into the removed state and balancing the system at a fixed, non-zero density of excited states for transient times $t<\frac{n_{c,0}}{\lambda}$.

The microscopic dynamics of the $d$-dimensional gas are given by the master equation ($\hbar=1$)

[TABLE]

for the ensemble density matrix $\hat{\rho}$. The coherent atom-light and atom-atom interaction is captured by the Hamiltonian

[TABLE]

while dissipative processes are described by the Liouvillian

[TABLE]

where $\Gamma=\gamma_{\text{de}}+\gamma_{\downarrow g}+\gamma_{\downarrow 0}$ is the sum of all dissipative rates. In typical experiments Helmrich et al. (2018, 2018); Valado et al. (2016), the product of the atomic mass $M$ and temperature $T$ is ’large’ compared to the density $n_{0}$, causing a thermal de-Broglie wavelength $\lambda_{\text{th}}=\frac{h}{\sqrt{2\pi Mk_{\text{B}}T}}$ much smaller than the mean free path $d_{a}\sim n_{0}^{-1/d}$. The motional degrees of freedom $\vec{r}_{l,l^{\prime}}$ thus cannot maintain coherence between two subsequent scattering events and are treated as classical variables undergoing thermal motion, see below and Ref. Helmrich et al. (2018).

We focus on a very large detuning $\Delta/\Gamma\sim O(10^{2}-10^{3})$ Urvoy et al. (2015); Gärttner et al. (2013); Lee et al. (2012), leading to strongly suppressed, off-resonant single particle transitions $|g\rangle\leftrightarrow|r\rangle$ at a rate $\tau_{t}\equiv\frac{\Gamma\Omega_{t}^{2}}{\Gamma^{2}+4\Delta^{2}}$. Due to interactions, an atom in the Rydberg state, however, creates a facilitation shell of radius $r_{\text{fac}}=(C_{6}/\Delta)^{1/6}$ and width $\delta r_{\text{fac}}\sim r_{\text{fac}}\Gamma/\Delta$. Inside the shell, the Rydberg repulsion compensates the detuning in Eq. (3), yielding an effective resonant excitation rate $\kappa_{t}\approx\Omega_{t}^{2}/\Gamma$ with $\kappa_{t}\gg\tau_{t}$ Ates et al. (2007); Amthor et al. (2010); Lesanovsky and Garrahan (2014); Faoro et al. (2016); Marcuzzi et al. (2016).

In the limit of strong dephasing, the atom coherences decay rapidly in time and the relevant dynamical degrees of freedom are the Rydberg state density and the density of ’active’ states, i.e., of atoms in the Rydberg and in the ground state. Their coarse grained values, averaged over a ’facilitation cluster’ of volume $V_{\text{fac}}=\frac{\pi^{\frac{d}{2}}r_{\text{fac}}^{d}}{\Gamma_{\text{Euler}}(\frac{d}{2}-1)}$ are

[TABLE]

The evolution equations for $\rho_{\vec{x},t}$ and $n_{\vec{x},t}$ are obtained by adiabatically eliminating the atom coherences from the Heisenberg equations of motion Marcuzzi et al. (2015, 2016); Valado et al. (2016); Buchhold et al. (2017); Pérez-Espigares et al. (2017); Gutiérrez et al. (2017). This yields the Langevin equation (see  Buchhold et al. (2017); Helmrich et al. (2018))

[TABLE]

Equation (7), describes four different processes on a coarse grained time scale $t\sim\gamma_{\text{de}}^{-1}$. It covers the average over Rabi oscillations inside each cluster, which occur with rate $\kappa_{t}\rho_{\vec{x},t}+\tau_{t}$ and prefer (averaged over time) a semi-excited state $\rho_{\vec{x},t}=\frac{n_{\vec{x},t}}{2}$. The rate combines the off-resonant oscillation rate $\tau_{t}$ and the resonant, facilitated rate $\kappa_{t}\rho_{\vec{x},t}$, which is proportional to the number of facilitating atoms $\rho_{\vec{x},t}$. This process competes with the linear decay channel $\sim\Gamma$, which prefers the ground state $\rho_{\vec{x},t}=0$.

The spreading of excitations from cluster to cluster is described by the diffusion term $\sim D\nabla^{2}\rho_{\vec{x},t}$, with $D=\kappa_{t}S$ being proportional to the facilitation rate and the surface $S$ of the clusters Helmrich et al. (2018). In a dissipative environment, each cluster experiences fluctuations of $\rho_{\vec{x},t}$, which are proportional to the oscillation rate Marcuzzi et al. (2015, 2016); Valado et al. (2016); Buchhold et al. (2017) and covered by the Markovian noise kernel (overline indicating noise average)

[TABLE]

Before turning to the evolution of $n_{\vec{x},t}$, we discuss the mean-field solution of Eq. (7) in the limit where $\tau_{t}\ll\Gamma,\kappa_{t}n_{\vec{x},t}$ by setting $D=\xi_{\vec{x},t}=0$. Defining a critical density $n_{c,t}\equiv\frac{\Gamma}{\kappa_{t}}$, one distinguishes two different regimes: an inactive regime for $n_{\vec{x},t}<n_{c}$, where the Rydberg density is suppressed and evolves towards $\rho_{\vec{x},t}\rightarrow\frac{\tau}{\Gamma}n_{\vec{x},t}$, and an active regime for $n_{\vec{x},t}>n_{c}$, where it evolves towards $\rho_{\vec{x},t}\rightarrow\frac{1}{2}(n_{\vec{x},t}-n_{c})$. The crossover between the two regimes at $n_{\vec{x},t}=n_{c,t}$ features a maximal correlation length of $\xi_{||}=\sqrt{\frac{D}{\sqrt{8\Gamma\tau_{t}}}}$. It turns into a sharp, second order phase transition in the limit $\tau_{t}\rightarrow 0$ Janssen (1981); Hinrichsen (2000); Marcuzzi et al. (2016, 2015); Buchhold et al. (2017).

The above discussed mean-field solution illustrates the dynamics in the regimes $\Gamma\ll\kappa_{t}n_{\vec{x},t}$ and $\Gamma\gg\kappa_{t}n_{\vec{x}t}$. In the presence of spatial fluctuations, i.e., for $D>0$, the asymptotic values for $\rho_{\vec{x},t}$ and $n_{\vec{x},t}$ above remain good approximations far away from the critical point $|\kappa_{t}n_{\vec{x},t}-\Gamma|\gg 1$. For $n_{\vec{x},t}\rightarrow n_{c,t}$, however, spatial fluctuations, manifesting via propagating avalanches with strongly fluctuating density, become increasingly strong and lead to deviations of the uniform behavior. In addition, the critical density is generally shifted towards larger values $n_{c}>\Gamma/\kappa$. In order to determine $n_{c}$ in $d=1,2$, we compute the location of the critical point in Eq. (7) numerically, e.g., we find $n_{c}=3.86$ in $d=1$.

The evolution of the density $n_{\vec{x},t}$ is governed by thermal motion of the atoms, the decay into the removed state and density fluctuations. It is summarized in the Langevin equation  Helmrich et al. (2018)

[TABLE]

with a Markovian noise kernel $\langle\eta_{\vec{x},t}\eta_{\vec{y},t^{\prime}}\rangle=\delta(\vec{x}-\vec{y})\delta(t-t^{\prime})\gamma_{\downarrow 0}\rho_{\vec{x},t}$ and a thermal diffusion constant $D_{n}$. It has minor impact on the dynamics but reduces geometrical constraints due to rare, inhomogeneous configurations of $n_{\vec{x},t}$ 222Any rare configuration with $n_{\vec{x},t}=0$ would otherwise block the spreading of excitations forever..

III Derivation of the Langevin equations

In this section, we present the detailed derivation of the Langevin equations (7) and (9) from the master equation (2). Readers interested in the effective dynamics may continue with its discussion in the following section.

Due to the exponential growth of the Hilbert space, the master equation Eq. (2) becomes too complex to solve for realistic, macroscopic system sizes. In order to reduce the complexity, the dynamics are projected onto the relevant long-wavelength degrees of freedom, i.e., the Rydberg density $\rho$ and the active density $n$ as defined in Eqs. (5) and (6). This procedure has been discussed for the case of $\lambda=\gamma_{\downarrow 0}=0$ in Refs. Marcuzzi et al. (2016); Buchhold et al. (2017) and for the case $\lambda=0,\gamma_{\downarrow 0}\neq 0$ in Ref. Helmrich et al. (2018).

For strong dephasing $\gamma_{\text{de}}\gg\Omega_{t}$ the decay of the atomic coherences $\sigma_{l}^{rg},\sigma_{l}^{r0}$ towards their steady state value is the fastest process in the quantum master equation. They can be adiabatically eliminated by formally solving the steady state equation for the average ($\alpha=g,0$)

[TABLE]

Inserting the solution of Eq. (10) and the completeness relation $\sigma_{l}^{rr}+\sigma_{l}^{gg}+\sigma_{l}^{00}=\mathds{1}$ into the full Heisenberg-Langevin equations for $\sigma^{rr}_{l},\sigma^{gg}_{l}$ yields

[TABLE]

The Markovian noise operators $\xi_{l}^{r,g}$ are added in order to enforce the fluctuation-dissipation relation of the driven dissipative master equation. They are local in space and time and fulfill the generalized Einstein relation (overline indicating noise average)

[TABLE]

This noise average leads to the $\delta$-correlated Markovian noise in Eq. (8) for the Langevin equation after the coarse graining procedure. Its crucial property for the realization of SOC is the scaling of the noise $\sim\rho_{\vec{x},t}$ (except for the tiny fluctuations $\sim\tau_{t}n_{\vec{x},t}$), which is responsible for a well defined, fluctuationless inactive phase.

Since the operators $\sigma^{rr}_{l},\sigma^{gg}_{l}$ are projection operators with eigenvalues $0,1$, any function $f$ of, say $\sigma^{rr}_{l}$, can be expressed as $f(\sigma^{rr}_{l})=f(0)+(f(1)-f(0))\sigma^{rr}_{l}$. Extending this to the whole set of $\{\sigma^{rr}_{l},\sigma^{gg}_{l}\}$, one rewrites

[TABLE]

This expression is exact up to second order powers in the projection operators. It separates off-resonant single particle transitions with rate $\tau$ and facilitated, two-particle transitions. For $2|\Delta-V^{I}_{l,l^{\prime}}|<\Gamma$, the facilitation rate deviates significantly from zero. Depending on the interaction potential, this defines the facilitation radius $r_{\text{fac}}$, i.e., for a typical van der Waals potential $V_{l,l^{\prime}}=\frac{C_{6}}{r^{6}}$ one finds $r_{\text{fac}}\equiv(\rm{C}_{6}/\Delta)^{1/6}$ and the facilitation shell $|\vec{r}_{l}-\vec{r}_{l^{\prime}}|\in[r_{\text{fac}}-\Delta r_{\text{fac}},r_{\text{fac}}+\Delta r_{\text{fac}}]$ with $\Delta r_{\text{fac}}=r_{\text{fac}}\frac{\Gamma}{12\Delta}$. We introduce a real space projector $\Pi_{ll^{\prime}}$ with $\Pi_{ll^{\prime}}=1$ if $|\vec{r}_{l}-\vec{r}_{l^{\prime}}|$ is inside the facilitation shell and zero otherwise. This yields

[TABLE]

This provides a good approximation for the facilitation rate when the density of excitations is small. For a number of $m\geq 1$ excited states inside a single shell, however, the exact solution shows a growth of the shell radius as $r_{\text{fac}}^{(m)}=m^{1/6}r_{\text{fac}}$ (in $d=3$ dimensions). This scaling behavior could be either taken into account by expanding Eq. (14) up to higher orders in the $\sigma^{rr}$ operators, which would account for a larger number $m>1$ of excitations per cluster, or by including the scaling of the facilitation volume for $m>1$ particles compared to the case of $m=1$. In both cases, the facilitation rate for $m>1$ then grows $\propto\sqrt{m}$, compared to the $\propto m$ prediction of Eqs. (15) and (14).If one bears in mind, however, the weak off-resonant excitation rate, configurations of $m\geq 1$ are suppressed by a factor $o(10^{-4})$. Our simulations show that $\rho_{\vec{x},t}\leq 1$ in most cases, which validates the restriction to $m=0,1$ in Eqs. (15) and (14).

The equation of motion for $\rho_{\vec{x},t}=\sum_{l}\Theta(r_{\text{fac}}-|\vec{x}-\vec{r_{l}}|)\langle\sigma^{rr}_{l}\rangle$ and $n_{\vec{x},t}=\sum_{l}\Theta(r_{\text{fac}}-|\vec{x}-\vec{r_{l}}|)\langle\sigma^{rr}_{l}+\sigma^{gg}_{l}\rangle$ yields

[TABLE]

and similar for $n_{\vec{x},t}$. For a homogeneous density, the drift term $\sim\partial_{t}\vec{r}_{l}$ can be approximated to be zero (see below for an inhomogeneous setting). This yields

[TABLE]

where $\mathcal{F}_{\vec{x}}(\rho_{\vec{z}},t)$ is some linear, quasi-local functional of $\rho_{\vec{x},t}$.

$\mathcal{F}_{\vec{x}}(\rho_{\vec{z},t})$ has support only around $|\vec{x}-\vec{z}|=r_{\text{fac}}$, enabling a Taylor expansion of the density. Since the Rydberg facilitation mechanism is isotropic in space, the expansion contains only even powers of derivatives. It reads as [cf. Eq. (14) in Ref. Helmrich et al. (2018)]

[TABLE]

The noise $\langle\xi_{\vec{x},t}\xi_{\vec{y},s}\rangle=\sum_{l,m}\Theta(r_{\text{fac}}-|\vec{x}-\vec{r_{l}}|)\Theta(r_{\text{fac}}-|\vec{y}-\vec{r_{m}}|)\langle\xi_{l,t}\xi_{m,s}\rangle=\delta(s-t)\delta(|\vec{x}-\vec{y}|)\left[\kappa_{t}\rho_{\vec{x},t}+\tau_{t}\right]$ remains Markovian and $\delta$-correlated on length scales of the facilitation radius.

Making a conservative estimate for the temperature of the motional degrees of freedom $T=O(10\mu K)$ and the atomic mass $M=O(20$u) Helmrich et al. (2018), one finds a thermal de Broglie wavelength $\lambda_{T}=\frac{h}{\sqrt{2\pi Mk_{\text{B}}T}}\approx 200$nm. For an atomic density of $n_{0}\approx 10^{11}\text{cm}^{-3}$ the mean free path in three dimensions amounts to d_{a}=\Big{(}\frac{6}{\pi n_{0}}\Big{)}^{1/3}\sim 2\mum, which is at least one order of magnitude larger than $\lambda_{T}$. Consequently, coherence in the motional degrees of freedom is lost between two subsequent scattering events and they can be treated classically. In the absence of an external trapping potential, the particles perform Brownian motion, i.e., thermal diffusion in a dilute van der Waals gas. This allows us to treat the atomic positions as slowly diffusing and uniformly distributed in space.

Including Brownian motion with diffusion constant $D_{n}$ the final form of the Langevin equations is

[TABLE]

Here $\kappa_{t}=\mathcal{F}_{\vec{x}}(1)\frac{\Omega_{t}^{2}}{\Gamma}$ is the facilitation rate. The diffusion constant $D=\mathcal{F}_{\vec{x}}(\vec{z}^{2})\frac{\Omega_{t}^{2}}{2\Gamma}(n_{\vec{x},t}-2\rho_{\vec{x},t})+D_{n}\approx\mathcal{F}_{\vec{x}}(\vec{z}^{2})\frac{\Omega_{t}^{2}}{2\kappa_{t}}$ is dominated by the facilitated spreading, which is proportional to the average density, i.e., $n_{\vec{x},t}-2\rho_{\vec{x},t}\approx\frac{\Gamma}{\kappa_{t}}$. This makes $D$, apart from local density fluctuations, time independent.

IV Self-organized criticality and avalanche dynamics

In order to observe self-organization towards a long-range correlated state, the dynamics should push any initial density $n_{\vec{x},0}$ close towards $n_{\vec{x},t}\rightarrow n_{c,t}$ and thereby maximize the correlation length $\xi_{||}$. This is achieved by the combination of loss into the auxiliary state $\sim\gamma_{\downarrow 0}$ and the continuously growing pump strength $\sim\lambda$.

Their interplay is best understood by expanding the critical density $n_{c,t}$ up to first order in $\lambda t$, yielding

[TABLE]

which is valid for $\lambda t<n_{c,0}$. For active densities $n_{\vec{x},t}\approx n_{c,t}$, the Rydberg state density $\rho_{\vec{x},t}$ experiences a large correlation length, leading to long-lived and and far spreading excitations, i.e., the formation of avalanches. Once an avalanche has formed, parts of it decay into the removed state, leading to a decrease of $n_{\vec{x},t}$. It reaches a stationary point when the decay of both $n_{\vec{x},t}$ and $n_{c,t}$ compensate each other, i.e., for $\lambda=\gamma_{\downarrow 0}\rho_{\vec{x},t}$.

On times $t<\frac{\lambda}{n_{c,0}}$, this is the only homogeneous solution of Eqs. (7) and (9) with

[TABLE]

It is reached after a time $t\approx\max\{\kappa_{t}^{-1},\gamma_{\downarrow 0}^{-1}\}$ and it survives up to times of order $t\approx\frac{n_{c,0}}{\lambda}$. On larger times, effects of order $\lambda^{2}t^{2}$ set in and the active density depletes to zero, i.e., $\rho_{\vec{x},t},n_{\vec{x},t}\rightarrow 0$.

Imposing a double separation of time scales on the dynamics via

[TABLE]

Eq. (21) predicts the self-organization towards a long-lived and long-range correlated state with $\rho_{\vec{x},t}=0^{+}$, $n_{\vec{x},t}=n_{c,t}+0^{+}$ and $\xi_{||}\rightarrow\infty$. We thus call (i) + (ii) the conditions for SOC in our driven Rydberg setup. The degree up to which both conditions are met, i.e., SOC is realized, can be adjusted experimentally via the Rabi frequency $\Omega_{t}$, the detuning $\Delta$ or the decay $\gamma_{\downarrow 0}$.

Such double separation of scales is a common requirement for realizations of SOC without energy conservation Bonachela and Muñoz (2009); Bonachela et al. (2010) 333This is contrasted with SOC in energy conserving systems, e.g., the sandpile model, which requires only a single pair of separated scales Bonachela and Muñoz (2009).. Since both our Hamiltonian and the Lindblad dynamics do not conserve the energy, the conditions (i)+(ii) can be seen as the present manifestations of this phenomenon. One may now argue that such strict requirements do not really differ from parameter fine tuning in conventional criticality. We, however, show that the dynamics of Eqs. (7) and (9) display SOC even for very weak realizations of (i) and (ii), e.g., for $\frac{\tau_{t}}{\lambda}\sim 10^{-4}$ and $\frac{\lambda}{\gamma_{\downarrow 0}}\sim 0.1$, making it accessible to experiments.

We emphasize that for $t<\frac{n_{c,0}}{\lambda}$ the increase of $\Omega_{t}$ with $\lambda$ is identical to loading ground state atoms with rate $\lambda$ into the system. The excitation avalanches of $\rho_{\vec{x},t}$ depend only on the difference $n_{\vec{x},t}-n_{c,t}=n_{\vec{x},t}+\lambda t-n_{c,0}$ and cannot distinguish between $n_{c,t}$ being decreased and $n_{\vec{x},t}$ being increased with rate $\lambda$. Experimentally, however, a controlled repopulation with rate $\lambda$ is often less feasible than adjusting the drive strength.

In order to confirm the prediction of emergent SOC from the homogeneous treatment above and to observe its paradigmatic avalanche dynamics, we simulate the full time evolution of the Rydberg density via Eqs. (7) and (9) in spatial dimensions $1\leq d\leq 3$. The equations are integrated on a $d$-dimensional grid of linear lattice spacing $\Delta x$ and we use dimensionless rates, expressed in units of $\Delta x^{2}/D$. The integration scheme is a derivative of the splitting scheme for stochastic differential equations with multiplicative noise Dornic et al. (2005), adapted to the noise kernel of Eq. (7), see Appendix A.

For the simulations we set $\kappa_{0}=\Gamma=\frac{\Delta x^{2}}{2D}$, $\tau_{0}=10^{-7}\Gamma$ and $\gamma_{\downarrow 0}=10^{-2}\Gamma$, which is consistent with recent experiments Helmrich et al. (2018, 2018); Gutiérrez et al. (2017). Different degrees of scale separation are realized by varying $\lambda$ within the interval $\lambda\in[0,0.2\Gamma]$. We point out that, as for our choice of parameters, any realistic experiment will realize the conditions (i) and (ii) only on an approximate level.

Our simulations reveal an extended dynamical regime, which is governed by the formation, propagation and decay of avalanches containing a significant number of excitations, $\rho_{\vec{x},t}\gg\frac{\lambda}{\gamma_{\downarrow 0}}$, (see Fig. 1c). Parametrically it coincides well with the criterion $\tau_{t}<\lambda<\gamma_{\downarrow 0}$, matching (i) and (ii). In general, the distribution $P_{\text{ava}}(s)$ of avalanche sizes $s$ varies with $\lambda$. In the vicinity of a critical value $\lambda\approx\lambda_{soc}$ it, however, approaches a scale invariant form $P_{\text{ava}}(s)\sim s^{-\alpha}$ with an exponent $\alpha>0$.

In $d=1$, we obtain $\alpha=1.44\pm 0.1$, which is consistent with results obtained from other SOC models, e.g., the forest fire model Schenk et al. (2002) or activity patterns in the cortex V. Stewart and Plenz (2006), and is associated with the underlying directed percolation universality class Hesse and Gross (2014). Its statistical error results from our sampling procedure, which dynamically counts avalanches from a finite number of patches of $10^{4}\times 10^{4}$ sites (time and space). For $d>1$, we predict $\alpha\approx 1.5$, however, with larger errors due to our avalanche counting scheme.

The scale invariant avalanche distribution is the hallmark of SOC Lu (1995); Watkins et al. (2016); Dickman et al. (2000). It is accompanied by fractal spatio-temporal Rydberg excitation patterns (see Fig. 1c) and paradigmatic $\frac{1}{\omega}$-fluctuations Bak et al. (1987, 1988) in the Rydberg density $\rho_{\vec{x},\omega}\equiv\int\rho_{\vec{x},t}e^{i\omega t}dt\sim\omega^{-\beta}$, with $\beta\lessapprox 1$ (see Fig. 2b). This clearly demonstrates a dynamical regime with SOC in the driven Rydberg gas. Its location at $\lambda\approx\lambda_{soc}$ can be understood as a trade-off in optimizing (i) and (ii) simultaneously for fixed values of $\tau_{t},\gamma_{\downarrow 0}$. For dimensions $d>1$ it approaches the estimate $\lambda_{soc}\sim\sqrt{\tau_{t}\gamma_{\downarrow 0}}$.

Moving $\lambda$ away from $\lambda_{soc}$, $P_{\text{ava}}(s)$ remains scale invariant in a finite range $|\lambda-\lambda_{soc}|<\eta$. We found $\eta\approx 0.2\lambda_{soc}$ for system sizes of $N=10^{6}$ lattice sites and our set of parameters. For larger deviations $|\lambda-\lambda_{soc}|>\eta$, the algebraic form of $P_{\text{ava}}(s)$ persists only for avalanche sizes $s<s_{||}(\lambda)$, i.e., below a $\lambda$-dependent cutoff scale $s_{||}(\lambda)$. Estimating the cutoff scale from the mean-field correlation length, i.e., $s_{||}(\lambda)=\xi_{||}$, which is justified far away from the SOC regime, one finds $s_{||}(\lambda)\sim\sqrt{\frac{D\gamma_{\downarrow 0}}{2\kappa_{t}\lambda}}$ for $\lambda\gg\tau_{t}$ and $s_{||}(\lambda)\sim\sqrt{\frac{D\lambda}{\kappa_{t}\tau_{t}}}$ for $\lambda\ll\gamma_{\downarrow 0}$.

The behavior on distances above $s_{||}$ in the two regimes $\lambda\lessgtr\lambda_{soc}$ manifestly differs from each other. For supercritical values $\lambda\gg\lambda_{soc}$, the critical density $n_{c,t}$ decreases rapidly, leading to a large avalanche triggering rate and a high density of avalanches. On sizes $s>s_{||}(\lambda)$ different avalanches start to overlap, which makes them indistinguishable and generates a random excitation pattern (displayed in Fig. 1c), revealing the underlying avalanches only for $s<s_{||}(\lambda)$, (squares in Fig. 1d).

The slow decrease of $n_{c,t}$ in the subcritical regime, $\lambda\ll\lambda_{soc}$ makes two subsequently following avalanches unfavorable and enforces a relative delay. It destroys the scale invariance above $s_{||}(\lambda)$ in favor of periodically triggered avalanches with increasing length $s\gg s_{||}(\lambda)$. This transforms the fractal real space structure found in the SOC regime into a time-periodic pattern, which is dominated by thermodynamically large excitation avalanches, shown in Fig. 1c. The period between two subsequent avalanches appears to be the time by which $n_{c,t}$ decreases by an integer value, i.e., $\delta t\approx\lambda^{-1}$.

Our simulations reveal that the conditions (i) and (ii) do not have to be fulfilled exactly in order to realize avalanche dominated dynamics and self-organized criticality. We find SOC also for a broader parameter regime, which is approximately described by the condition

[TABLE]

This condition can serve as a rule of thumb for the realization of self-organized criticality in experiments on driven Rydberg ensembles.

V Experimental observability

While the real space evolution of excitation avalanches is hard to access in experiments, the statistics of excitations, i.e., $\rho_{\vec{x},t}$ and $n_{\vec{x},t}$, can be measured via the particle loss rate $\propto\gamma_{\downarrow 0}\rho_{\vec{x},t}$ Helmrich et al. (2018, 2018). A robust, time-translational invariant observable is the integrated density

[TABLE]

where $n_{0}$ is the total initial density and $\langle...\rangle_{V}=\frac{1}{V}\int_{V}d^{d}x$ denotes the spatial average over the system volume. Its meaning becomes clear when comparing it with the initial critical density $n_{c,0}$ at times $t\lambda\ll n_{c,0}$, yielding $R_{t}-n_{c,0}=\langle n_{\vec{x},t}\rangle_{V}-n_{c,t}$.

Both $\rho_{\vec{x},\omega}$ and $R_{t}$ display very characteristic features in the three different regimes. For subcritical $\lambda$, the real time evolution of $R_{t}$ shows large, periodic amplitude fluctuations, reflecting individual, periodically triggered, extended avalanches. Instead, both the SOC and the supercritical regime feature much smaller amplitude fluctuations around $R_{t}\approx n_{c,o}$ (SOC) or $R_{t}\gg n_{c,o}$ (supercritical) as shown in Fig. 2a. In the subcritical (supercritical) regime, $\rho_{\vec{x},\omega}$ departs from its scale invariant form at SOC and one finds instead suppressed (pronounced) density fluctuations at intermediate frequencies, see Fig. 2b.

Significant information is encoded in the statistics of $R_{t}$, especially its mean $\bar{R}\equiv\lambda\int_{0}^{\lambda^{-1}}R_{t}dt$ and fluctuations $\sigma_{R}^{2}\equiv\lambda\int_{0}^{\lambda^{-1}}R_{t}^{2}dt-\bar{R}^{2}$ as displayed in Fig. 2c. For subcritical $\lambda$ both $\bar{R}$ and $\sigma_{R}$ increase with $\lambda$ faster than the linear mean-field prediction. At the onset of SOC, however, both $\bar{R}$ and $\sigma_{R}$ experience a sharp drop, manifest in a non-analytic kink in their $\lambda$-dependence. While $\bar{R}\rightarrow n_{c,0}$ rapidly approaches the critical density, the fluctuations decrease by several orders of magnitude. Upon further increasing $\lambda$, $\bar{R}$ reaches a valley at $\approx n_{c,0}$ and subsequently increases again into the supercritical regime. $\sigma_{R}$ is featureless at the SOC-supercritical transition.

In order to reason the observability of SOC for realistic conditions, where the atomic cloud is confined inside a trap, we expose $n_{\vec{x},t}$ to a potential of the form $V_{\text{trap}}(\vec{x})=V_{0}\exp(-|\vec{x}|^{2}/\xi_{\text{trap}}^{2})$, e.g., resulting from a Gaussian trapping laser with beam waist $\xi_{\text{trap}}$ Helmrich et al. (2018). For a mean free path $d_{a}\ll\xi_{\text{trap}}$, the effect of $V_{\text{trap}}(\vec{x})$ can be treated within the relaxation time approximation, see Appendix B. This adds a drift $\sim-\vec{v}_{\vec{x}}\cdot\vec{\nabla}n_{\vec{x},t}$ to the right-hand side of Eq. (9). Here $\vec{v}_{\vec{x}}=\frac{d_{a}}{\sqrt{Mk_{\text{B}}T}}\vec{\nabla}V_{\text{trap}}$ is the relaxation velocity. The dynamics following this drift at low temperatures $T$ ($\frac{V_{0}d_{a}}{\sqrt{Mk_{\text{B}}T}}=0.7D$) is displayed in Fig. 3a. On distances $|\vec{x}|<\xi_{\text{trap}}$, avalanches remain well defined and both their fractal real space pattern and the scale invariant statistics are observable below the trap scale, see Fig. 3b.

VI Effect of the spatial dimension

Apart from Rydberg atoms, the continuum model in Eq. (7) may also serve as a coarse grained description for activity spreading in sparse networks Hesse and Gross (2014). In this picture, each Rydberg atom represents a node and the parameters $\kappa,\tau,\Gamma$ describe its reaction to external stimuli and the decay of information. The density $n_{\vec{x},t}$ represents a ’node energy’, which is consumed by active nodes with rate $\gamma_{\downarrow 0}\rho_{\vec{x},t}$ and recharged with rate $\lambda$.

Optimal networks are expected to operate close to SOC de Arcangelis et al. (2006); Shew et al. (2015); Hesse and Gross (2014); Marković and Gros (2014); Kinouchi and Copelli (2006); Levina et al. (2007). Their natural tuning parameter is the average connectivity $z$ of the nodes, which is adjusted to match external conditions Levina et al. (2009); Bornholdt and Rohlf (2000); Bertschinger and Natschläger (2004); Kinouchi and Copelli (2006); Levina et al. (2007). Figure 2c confirms that here the dimensionality $d$ acts as a second ’control parameter’. Changing $d$ from $d=1$ to $d=2$ shifts the scale invariant regime (shaded region) and increases its range. For a given set $\tau,\lambda,\gamma_{\downarrow 0}$, there may exist an ’optimal’ $d$, for the system to display SOC. In Rydberg experiments $d$ can be controlled by adjusting the trapping geometry. Combined with the tuneability of $\lambda$ and $\tau$, this offers many possibilities to study self-organized criticality in network-like setups.

VII Conclusion

We propose and study an experimentally feasible mechanism to control excitation avalanches in driven Rydberg setups Helmrich et al. (2018). On large, transient times, one can observe subcritical, supercritical and self-organized critical avalanche dynamics, depending on the control parameter. Each regime features unique signatures, including a scale invariant avalanche distribution and $\frac{1}{\omega}$-noise, both paradigmatic signals for SOC. This motivates driven Rydberg ensembles Helmrich et al. (2018) as viable platforms for the study of SOC and the conditions under which simple dynamical rules, as imposed by the facilitation condition, can establish and maintain self-ordering towards complex dynamics structures.

While the crossover from the SOC to the supercritical regime does not produce a pronounced feature in the integrated density, Fig. 2 reveals a developing non-analyticity in both the integrated density as well as its fluctuations as $\tau_{t}$ is decreased. It hints towards an underlying critical point, on the one hand such a critical point might describe the SOC universality class, including avalanche and correlation exponents. On the other hand, it could be a remnant of the directed percolation critical point, which would be reached for $\lambda,\tau\rightarrow 0$. In both cases, the investigation of this conjectured critical point and its relation to the SOC universality seems worthwhile for future work.

Based on the similarity of the corresponding master equations, we conjecture a relation between driven Rydberg gases and self-organizing neural networks. The analogy is strengthened by frequently observed periodic or random activity patterns in non-optimal operating networks Prinz (2008); Ong et al. (2012). Exploring this connection, especially for the role that is played by scale separation, appears a promising direction to connect driven Rydberg systems with neurosciences.

Acknowledgements.

We thank G. Refael, S. Diehl and S. Whitlock for valuable comments on the manuscript. K. K. was supported by the J. Weldon Green SURF fellowship and M. B. acknowledges support from the Alexander von Humboldt foundation.

Appendix A Numerical integration scheme

Numerical integration of Eqs. (7) and (9) is performed by an operator-splitting update scheme Dornic et al. (2005). At each time step, the evolution is decomposed into a stochastic evolution step and a deterministic step. The former is designed to solve a stochastic differential equation of the form:

[TABLE]

Here $\eta$ is a Markovian noise kernel with mean zero and unit variance. For small $\gamma_{\downarrow 0},\kappa,\tau$, we may approximate $\alpha$ and $\beta$ to be constant over each time step. The corresponding Fokker-Planck equation has the exact solution

[TABLE]

where we set $\rho\equiv\rho_{\vec{x},t+\delta}$ and $\rho_{0}\equiv\rho_{\vec{x},t}$ as well as $\lambda=\frac{2\beta}{\sigma^{2}(e^{\beta t}-1)}$ and $\mu=\frac{2\alpha}{\sigma^{2}}-1$ and $I_{\mu}(x)$ is the modified Bessel function of the first kind with index $\mu$ and argument $x$. This can be expressed via a mixed Gamma distribution which allows for efficient sampling:

[TABLE]

which is shorthand notation for a random variable which is drawn from a Gamma distribution with argument $\mu+1+x$, whereas $x$ was drawn from a Poisson distribution with argument $\lambda\rho_{0}e^{\beta\delta t}$.

Given the values of $\rho_{\vec{x},t}$ at time $t$, its stochastic evolution $\rho_{\vec{x},t+\delta t}$ after a step $\delta t$ can be drawn from the above distribution. The deterministic part of the equation of motion has a purely polynomial form and can also be solved exactly. The time discretization error is therefore only caused by the splitting of the evolution into a stochastic and a deterministic part.

A non-zero $\tau_{t}$ can be incorporated by using the same procedure with a simple change of variables: $u=\rho+\tau_{t}/\kappa_{t}$. The non-negativity of $\rho$ is enforced after sampling by resetting any value of $u<\tau_{t}$ to $\tau$. The well-behaving evolution of $n_{\vec{x},t}$ is performed via an Euler scheme.

Appendix B Relaxation time approximation in a trap

In the presence of an inhomogeneous background potential $V(\vec{r})$ for the particles, the drift term in Eq. (16) becomes significant. For the active density it yields

[TABLE]

where we applied the chain rule and inserted the momentum $\vec{p}_{l}=M\partial_{t}\vec{r}_{l}$. In the relaxation time approximation, the momentum $\vec{p}$ is reset after a characteristic scattering time $t_{\text{rel}}=d_{a}{\sqrt{\frac{M}{2\pi k_{\text{B}}T}}}$, where $d_{a}$ is the mean free path and $T$ is the temperature. This yields the equation of motion

[TABLE]

It is stationary for $\vec{p}_{l}=-t_{\text{rel}}\vec{\nabla}V(\vec{r}_{l})$ and induces an average drift for times $t>t_{\text{rel}}$. Inserting this result in Eq. (B) and neglecting the variation of $V$ on length scales $\sim r_{\text{fac}}$, i.e., $V(\vec{r}_{l})\approx V(\vec{x})$, one finds

[TABLE]

where $...$ describes the dynamics of the internal states of the atoms. This approximation works well if both the facilitation shell and the mean free path are much smaller than the typical length scale of the potential $V$.