Hydrodynamic stabilization of self-organized criticality in a driven Rydberg gas
K. Klocke, T. M. Wintermantel, G. Lochead, S. Whitlock, M. Buchhold
HHydrodynamic stabilization of self-organized criticality in a driven Rydberg gas
K. Klocke,
1, 2
T. M. Wintermantel,
3, 4
G. Lochead, S. Whitlock, and M. Buchhold
1, 5 Department of Physics and Institute for Quantum Information and Matter,California Institute of Technology, Pasadena, CA 91125, USA Department of Physics, University of California, Berkeley, California 94720, USA ISIS (UMR 7006), University of Strasbourg and CNRS, 67000 Strasbourg, France Physikalisches Institut, Universit¨at Heidelberg, 69120 Heidelberg, Germany Institut f¨ur Theoretische Physik, Universit¨at zu K¨oln, D-50937 Cologne, Germany (Dated: September 29, 2020)Signatures of self-organized criticality (SOC) have recently been observed in an ultracold atomic gas undercontinuous laser excitation to strongly-interacting Rydberg states [S. Helmrich et al. , Nature, , 481– 486(2020)]. This creates a unique possibility to study this intriguing dynamical phenomenon, e.g., to probe itsrobustness and universality, under controlled experimental conditions. Here we examine the self-organizingdynamics of a driven ultracold gas and identify an unanticipated feedback mechanism, which is especiallyimportant for systems coupled to thermal baths. It sustains an extended critical region in the trap center for anotably long time via hydrodynamic transport of particles from the flanks of the cloud toward the center. Thiscompensates the avalanche-induced atom loss and leads to a characteristic flat-top density profile, providing anadditional experimental signature for SOC and minimizing e ff ects of inhomogeneity on the SOC features. Introduction. – Many-body systems, may they be driven,open or excited by a sudden parameter quench, often evolvetoward steady or transient metastable states which can beclassified as far from thermal equilibrium [1–13]. Some-times these systems feature attractors for the non-equilibriumdynamics that give rise to emergent scale invariant proper-ties over a wide range of initial states or parameters [14–18]. One paradigmatic example is self-organized criticality(SOC), whereby a dissipative many-body system evolves to-ward a (non-equilibrium) critical state by an intrinsic feed-back mechanism. Since its first introduction by Bak, Tang,and Wiesenfeld in 1987 [12, 19], SOC has been intensivelystudied theoretically and associated with phenomena rangingfrom avalanches and earthquakes to solar flares and neuronalactivity to name a few [20–24].The range of phenomena found to exhibit SOC-like char-acteristics is at odds however with the relatively stringentconditions that are expected to lead to SOC in theory [25].For example, the typical requirements of a large separationof timescales between slow dissipation and fast, conserva-tive bulk dynamics will never be perfectly satisfied in prac-tice [26]. This has lead to the notion of self-organizedquasi-criticality (SOqC) where the system hovers around crit-icality with large excursions into the sub- and super-criticalphases [26, 27]. Nonetheless, key signatures of self-organizedcriticality including scale invariance of the stationary densityand power-law distributed excitation avalanches were recentlyobserved in the driven-dissipative dynamics of atomic Ryd-berg gases [28] (see also related experiments in driven thermalgases [29]). These experiments, however, lacked an obviousrefilling mechanism necessary to bring the system out of thesub-critical absorbing phase. This therefore raises importantquestions about how signatures of the SOC state persist forrelatively long times and whether possible universal character-istics of the SOC state [25] can be extracted from experimentsin a transient regime. excited state densityposition lossrearrangementwithoutrearrangement"ideal"SOCtrap (a) (b)(c) (ii-a) (ii-a)(ii-b) (ii-b)(iii) growth (i) A t o m i c den s i t y Figure 1. Mechanisms for self-organized criticality in an ultracoldatomic gas. (a) A trapped atomic gas with inhomogeneous densitydistribution is continuously driven by an o ff -resonant excitation laserto highly excited Rydberg states (blue disks). (b) Trajectory of theatom density n t and the excitation density ρ t driven by facilitatedexcitation, decay and hydrodynamic motion. Starting from the su-percritical phase n t = > n c the system undergoes: (i) rapid growthof Rydberg density; (ii-a) self-organization from the active phase to-ward the critical point due to gradual depletion of particles (causedby loss from the Rydberg state); (ii-b) refilling of the central densityfrom the sub-critical phase by atomic rearrangement (thermal mo-tion) from the wings to the trap center; (iii) stabilization close to thecritical point for an extended period of time. (c) upper panels: exper-imental absorption image ( n (cid:126) x , t integrated over z ) at di ff erent times.lower panels: Reconstructed three-dimensional atom density n (cid:126) x , t at y = z = Here we experimentally and theoretically demonstrate that a r X i v : . [ c ond - m a t . qu a n t - g a s ] S e p the mechanisms leading to the SOC state are remarkably ro-bust. We show that the slow motion of the particles providesan additional feedback mechanism which stabilizes the sys-tem close to the critical state over an extensive period of time.This is evidenced by the experimental observation of a sta-ble flat-top profile in the atomic gas, where the wings of thedistribution act as particle reservoirs that compensate parti-cle loss in the trap center (Fig. 1). To explain this result wedevelop a hydrodynamic Langevin equation which describesthe competition between thermalization of the gas (in the mo-tional degrees of freedom) and the driven-dissipative excita-tion dynamics leading to SOC. This allows the cloud to adaptby slowly refilling sub-critical regions back to a critical state,which plays a similar role to plasticity in biological neuralnetworks [30, 31]. Self-organization mechanism. – We consider a spatiallyinhomogeneous gas of ultracold atoms held in a harmonicoptical potential produced by a focused far-o ff -resonant laserbeam [28] (depicted in Fig. 1a). The atoms are continuouslydriven by a detuned laser field, which creates rare and isolatedRydberg excitations at random positions in the gas. Once anexcitation is present it will either spontaneously decay (whichis often accompanied by loss from the trap), or it can trig-ger secondary excitations through a process called Rydbergfacilitation [32–35]. This occurs at a characteristic distance r fac ≈ . µ m (for the present experiments) where the laserdetuning is compensated by the van der Waals interaction be-tween Rydberg pair states [36]. The self-organizing dynamicsare driven by the competition between facilitated excitation(with rate proportional to κ n (cid:126) x , t , where κ is the microscopicfacilitation rate and n (cid:126) x , t is the local density of ground-stateatoms) and density independent spontaneous decay or loss ofthe excited atoms with rate Γ . These two processes compete toproduce rich collective dynamics [28, 37–40] and become bal-anced at a critical atom density n c ≈ Γ /κ . For n (cid:126) x , t > n c (super-critical or active phase) individual excitations can grow intospatially extended clusters of excitations (avalanches) with ahigh degree of activity and particle loss. For n (cid:126) x , t < n c (absorb-ing phase), on the other hand, excitation avalanches are rareor vanishingly small.Figure 1 illustrates the mechanisms leading to SOC. Start-ing from the supercritical regions of the cloud the density ofexcitations ρ (cid:126) x , t undergoes a period of rapid growth [labeled (i)in Fig. 1b], followed by a slow decrease in both n (cid:126) x , t and ρ (cid:126) x , t owing to a gradual loss of excited atoms [labeled (ii-a)]. Inthe limit of vanishingly small loss rate (perfect separation oftimescales) the system will follow a characteristic trajectory(dashed blue curve in Fig. 1b) that terminates at the criticalpoint [orange cross at n (cid:126) x , t = n c and ρ (cid:126) x , t = t [ ms ] R y db e r ga t o m nu m b e r ( i )( ii ) ( iii ) → (a) −
500 500 x [ µm ] . . . n ~ x , t [ µ m − ] n c (b) t [ ms ] t [ ms ] R y db e r a t o m nu m b e r (c) −
100 100 x [ r fac ] . . . n ~ x , t [ µ m − ] n c (d) t [ ms ] t [ ms ] T o t a l a t o m nu m b e r −
100 100 x [ r fac ] n ~ x , t t = 10[ ms ] -85 85 x [ r fac ] -1010 y [ r f a c ] ExperimentTheory
Figure 2. Theory-experiment comparison showing the approach tothe SOC state (top: experiment, bottom: theory). (a) Instantaneousnumber of Rydberg excitations integrated over the cloud ( ∝ ρ t ). Eachdata point is obtained from a destructive measurement and corre-sponds to a distinct experimental realization. Inset: correspondingtotal atom number. (b) One-dimensional slices through the atomicdensity distribution n (cid:126) x , t at y = z =
0. (c) Simulated dynamics ofthe full time-evolution (red line: single trajectory, grey lines: over-lapped data of six di ff erent trajectories) showing temporally well sep-arated, extensive excitation avalanches that persist long after the ini-tial growth and self-organizing regimes (i) and (ii). (insets: snapshotsof the peak excitation density per avalanche at z = y = z = n (cid:126) x , t = n c analogous tothe experimental observations in (b) (see also inset). a refilling mechanism to escape the absorbing phase and ap-proach the critical point [red curve in Fig. 1a, labeled (ii-b)].This interplay of nonlinear excitation dynamics and atomicmotion explains how the system self-organizes close to thecritical point with a constant critical density across the cloudand sustains critical dynamics (e.g. avalanches) for long timescompared to the initial self-organization period [labeled (iii)]. Experimental approach. – Our experiments start with an ul-tracold gas of N = potassium-39 atoms trapped in a cigar-shaped optical potential with trap frequencies of ω x / π =
65 Hz and ω y , z / π =
950 Hz. The atomic cloud has a tem-perature of T = µ K and e − / radii of σ x = µ m, σ y , z = µ m with a peak density of n = . µ m − . Attime t =
0, we switch on an o ff -resonant ultraviolet (UV)laser coupling with Rabi frequency Ω ≈ ∆ / π =
30 MHz on the transition from the groundstate | g (cid:105) = | s / , F = (cid:105) to the Rydberg state | r (cid:105) = | p / (cid:105) .To strongly suppress single-particle excitations and to en-sure that many-body e ff ects dominate, we stay in the regime Γ (cid:28) Ω (cid:28) ∆ [37, 41–52]. Excitations decay with a calcu-lated rate of Γ / π = .
84 kHz, which either brings them backto the ground state | g (cid:105) or into states | (cid:105) which are not cou-pled or are lost from the trap. This loss of particles into inac-tive states | r (cid:105) → | (cid:105) provides the crucial feedback mechanismfor SOC [28]. After the laser exposure time t , we measureboth the number of Rydberg excitations in the cloud as wellas the spatial distribution of ground state atoms remaining inthe trap. For the former we first ionize the Rydberg atomsand detect them on a micro-channel plate detector. For thelatter we take an absorption image of the atom cloud, whichintegrates along the propagation of the light field [28].Example absorption images after di ff erent exposure times t are shown at the top of Fig. 1c, roughly coinciding with theones sketched in Fig. 1a. The line profiles shown in Fig. 1care reconstructed cross-sections of the three-dimensional den-sity distribution through the center of the cloud. These areobtained by an inverse Abel transformation [53], using radialsymmetry along the elongated axis of the three-dimensionalcigar-shaped cloud. Initially at t = t = t (cid:39)
15 ms, this dip has filled in and the cloud shows a flat-topcoinciding with the critical density n c . Theoretical description. – The collective dynamics of thedriven Rydberg ensemble is described by a nonequilibriumfield theory for the local density of particles n (cid:126) x , t and the den-sity of excitations ρ (cid:126) x , t [28]. Besides the facilitated spreadingof excitations and the dissipative decay, here we also accountfor the hydrodynamic motion of the atoms via two coupledstochastic evolution equations for ρ (cid:126) x , t , n (cid:126) x , t , including the in-ternal and external degrees of freedom.We label each atom with an index j , a set of operators σ αβ j = | α (cid:105) (cid:104) β | j where α, β label the states 0 , g , r , and a posi-tion (cid:126) x l (treated as classical variable). The equation of motion(EOM) for the internal degrees of freedom is given by the mi-croscopic Liouvillian ∂ t σ αβ l = i (cid:18)(cid:88) j (cid:44) l C σ rrj | (cid:126) x l − (cid:126) x j | − ∆ (cid:19) σ rrl + Ω σ rgl + σ grl , σ αβ l + δ αβ (cid:16) δ α r γ de + δ α g γ ↓ g + δ α γ ↓ (cid:17) σ αα l − Γ { σ rrl , σ αβ l } , with the anti-commutator {· , ·} , the commutator [ · , · ] and theKronecker symbol δ α,β . This includes coherent single-particleprocesses: laser driving with Rabi frequency Ω and detun-ing ∆ , and the van der Waals interaction between atoms l and j if both are in the Rydberg state. The dissipative single-particle processes are quantified by the dephasing rate γ de , thespontaneous decay rate γ ↓ g for the process | r (cid:105) → | g (cid:105) ( γ ↓ for | r (cid:105) → | (cid:105) ) and Γ = γ de + γ ↓ g + γ ↓ .In order to apply a coarse grained description, wedefine a unit cell as the sphere with radius r fac andvolume V fac . The densities per unit cell are [28] ρ (cid:126) x , t = (cid:80) (cid:48) j ,(cid:126) x (cid:104) σ rrj (cid:105) / V fac , n (cid:126) x , t = (cid:80) (cid:48) j ,(cid:126) x (cid:104) σ rrj + σ ggj (cid:105) / V fac where (cid:80) (cid:48) j ,(cid:126) x isrestricted to j with | (cid:126) x j − (cid:126) x | ≤ r fac . The EOM for the atomicdensity is evaluated by applying the chain rule ∂ t n (cid:126) x , t = (cid:48) (cid:88) j ,(cid:126) x ∂ t (cid:104) σ rrj + σ ggj (cid:105) V fac − ∇ (cid:48) (cid:88) j ,(cid:126) x (cid:104) σ rrj + σ ggj (cid:105) V fac ∂ t (cid:126) x l , (1)where ∇ = ( ∂ x , ∂ y , ∂ z ). It contains the EOM for the internaldegrees of freedom and for the position of the atoms. The sum over the velocities in Eq. (1) is by definition the coarsegrained current (cid:126) j .An equivalent computation for ∂ t (cid:126)ρ (cid:126) x , t yields the Langevinequation [40, 54] ∂ t ρ (cid:126) x , t = ( D ∇ − Γ ) ρ (cid:126) x , t + ( τ + κρ (cid:126) x , t ) (cid:0) n (cid:126) x , t − ρ (cid:126) x , t (cid:1) + ξ (cid:126) x , t . (2)Here the evolution within each unit cell is decomposed interms of facilitated (de-)excitation with rate κρ (cid:126) x , t ≈ Ω V fac ∆ ρ (cid:126) x , t and dissipative decay ∼ Γ . Excitations spread di ff usively be-tween unit cells with Dr fac ≈ κ . Rare, o ff -resonant single-particle excitations occur with rate τ r = κ Γ∆ ≈ − κ , act-ing as local seeds to prevent the system from getting stuck inan absorbing state. Local fluctuations in the excitation den-sity are described by a multiplicative Markovian noise ξ ( (cid:126) x , t )with auto-correlation function (cid:104) ξ ( (cid:126) x , t ) ξ ( (cid:126) y , t (cid:48) ) (cid:105) = δ ( (cid:126) x − (cid:126) y ) δ ( t − t (cid:48) ) (cid:0) Γ ρ (cid:126) x , t + τ (cid:1) [54].The EOM of the density n (cid:126) x , t yields ∂ t n (cid:126) x , t = −∇ (cid:126) j (cid:126) x , t − γ ↓ ρ (cid:126) x , t , (3)where the current (cid:126) j (cid:126) x , t = − ( D T ∇ + η ∇ V (cid:126) x ) n (cid:126) x , t fits the commonhydrodynamic form [55]. It includes di ff usion ( ∝ D T ) and anexternal force −∇ V (cid:126) x caused by the harmonic trapping poten-tial V (cid:126) x = M (cid:80) l = x , y , z ( ω l (cid:126) x l ) , for which we use the frequencies ω l and the atom mass M from the experiment. The mobility η is related to the di ff usion constant D T = η k B T via the Ein-stein relation. In the limit γ ↓ →
0, the steady state has zerocurrent (cid:126) j (cid:126) x , t = n (cid:126) x , t = n (eq) (cid:126) x = n exp( − V (cid:126) x k B T ). For the numerical simulationwe use the trap frequencies ω l and temperature T measured inthe experiment, and the initial spatial extension of the cloudis σ l = (cid:113) k B TM ω l , i.e., σ z = σ y = . r fac and σ x = r fac . For γ ↓ > η M ω x . In the EOMfor ρ (cid:126) x , t , particle motion is negligible compared to the facili-tated spreading of excitations, i.e., ∼ Dr − (cid:29) η M ω x .The Langevin equations are integrated numerically on a3 + ρ (cid:126) x , t , n (cid:126) x , t ) → (0 , n c ) as well as at large system sizes and longtimes. The parameters used in the simulations are chosento match the experimentally observed facilitation and decayrates [28, 58] as well as the real-space extension of the cloud.The mobility η is hard to quantify from microscopic param-eters alone and it was chosen such that the theoretical andexperimental thermalization times match. This respects theseparation of time scales between the fast spreading of exci-tations ∼ r / D = . ∼ / ( η M ω x ) = Dynamics. – Figure 2 shows comparable experiments andnumerical simulations for an initially Gaussian atomic cloudwith peak density n > n c . Looking at both compo-nents n (cid:126) x , t , ρ (cid:126) x , t provides insights into the di ff erent dynamicalregimes. This includes an initial growth regime (i), cover-ing the first few milliseconds of evolution and resulting ina macroscopic Rydberg population. The early time growthdynamics are interesting in their own right [59], but are notoverly important for the self-organizing behavior on longertimescales, apart from bringing the system into the initiallysupercritical phase. Subsequently the loss from the Rydbergstate begins to decrease the total atom number. This leads to aself-organizing regime (ii), evidenced by large bursts of Ryd-berg excitations (large activity seen in Fig. 2a,c) and a suddendrop in the central density of the atomic cloud. The densityapproaches a flat-top distribution with a central density givenby the critical value n (cid:126) x , t ≈ n c ( ≈ / n (cid:126) x , t = n c (red curves in Fig. 2b,d) and sporadicavalanche-like excitation events. This is reached after approx-imately 15 ms in the experiment and persists until at least150 ms. Simulations show that subsequent avalanches arewell separated in space and time, which implies that the ex-perimentally observed Rydberg excitation spikes correspondto individual avalanche events (Fig. 2c). In this regime eachavalanche event transiently imprints a slight extra depressionin the density profile such that n (cid:126) x , t < n c . However, particletransport from the flanks re-establishes n (cid:126) x , t ≈ n c between suc-cessive avalanches (Fig. 2) and pins n (cid:126) x , t = n c (witnessed byFig. 2d) thus sustaining a close to ideal critical SOC state overa large region of the system.In order to quantify the characteristic timescale associatedto this mechanism, we investigate the e ff ective refilling rate ofthe central region λ . A necessary condition for maintaininga SOC state is to satisfy a common separation of timescales γ ↓ (cid:29) λ (cid:29) τ [23, 57]. The refilling rate is determined bythe gradient of the particle current λ n (cid:126) x , t ≡ −∇ (cid:126) j (cid:126) x , t from thewings towards the center. In order to estimate λ , we apply amean-field approach based on our observation that the currentis dominated by particle flow along the elongated x -direction.Therefore, we assume a quasi-one-dimensional cloud with aflat-top of length L x and a constant density n x , t = ¯ n t ≥ n c in-side the center. Outside, the density is sub-critical and followsthe equilibrium profile n x , t = ¯ n t (cid:16) n (eq) x / n (eq) L x / (cid:17) , which minimizesthe current (cid:126) j x , t = ρ x , t = n t λ = − L x (cid:90) | x |≤ L x / dx ∂ x j x , t = η M ω x ¯ n t . (4)Using Eq. (4), we estimate γ ↓ /λ ≈
50 and λ/τ ≈ ∼ /γ ↓ ). It also ensures that the refilling of the central den-sity happens much faster than o ff -resonant excitations ( ∼ /τ ),leading to well-separated avalanches, fulfilling the necessaryconditions for SOC [23].Finally, we analyze the very late time dynamics, in whichthe cloud reaches thermal equilibrium. The thermal state isapproached when the particle reservoir represented by the −
500 500 x [ µm ] t [ m s ] −
100 100 x [ r fac ] t [ m s ] n ~x,t t [ms] − . − . . E x ce ss K u rt o s i s TheoryExperiment . . . a) b)c) Figure 3. Melting of the flat-top from the time-evolved densityprofile n (cid:126) x , t (projected onto the x -axis). (a) Experimental measure-ments at di ff erent times and (b) the simulated evolution show astable flat-top with a lifetime exceeding 200 ms and a boundary(white dashed curve), which slowly approaches the center. Bothplots extend over the same x -axis distance. (c) The equilibrationof the cloud profile is quantified by its time-dependent excess kur-tosis EK t = (cid:82) dx ( x /σ x ) n (cid:126) x , t −
3, where we integrate over a densityslice with y = z =
0. Here σ x is the width of the cloud in the x -direction. Starting from a Gaussian shape (EK = ≈ − flanks is continuously depleted, which in turn leads to themelting of the flat-top region (Fig. 3a,b). The approach tothermal equilibrium can be seen by the evolution of the excesskurtosis EK t , shown in Fig. 3c. A non-zero kurtosis serves asa measure for the deviation of the cloud shape from a thermalGaussian distribution, i.e., it measures relative flatness of thedistribution. Its relaxation monitors the melting of the flat-toptowards a robust, thermal equilibrium state without excitationoutbursts (corresponding to EK t = ≈
20 ms).
Conclusion. – We have identified an important additionalmechanism which explains how SOC can be sustained in adriven-dissipative ultracold atomic gas by nonequilibrium cur-rents. We show that this generates a flat-top density distri-bution at the SOC critical density, quantitatively confirmedby the hydrodynamic Langevin theory. This demonstratesan important additional signature for SOC that could helpidentify SOC-like behavior in other systems, such as room-temperature atomic vapours and cold molecular plasmas [29,60]. Similar mechanisms may also be at play in very di ff erentsystems including adaptive neural networks [30, 31]. The factthat the system naturally evolves to a stable, mostly homo-geneous shape combined with the e ff ectiveness of the hydro-dynamic Langevin theory will enable more stringent tests ofnon-equilibrium universality in SOC systems. Alternatively,the interplay between internal and external degrees of free-dom could lead to other rich dynamical regimes to test, suchas oscillatory behaviour associated with SOqC [26, 27, 61]. ACKNOWLEDGEMENTS
K.K. acknowledges support from the National ScienceFoundation through grant DMR-1723367. T.M.W. acknowl-edges the French National Research Agency (ANR) throughthe Programme d’Investissement d’Avenir under contractANR-17-EURE-0024. This project is part of and supportedby DFG SPP 1929 GiRyd through projects DI1745 / / [1] H. K. Janssen, Zeitschrift f¨ur Physik B Condensed Matter ,151 (1981).[2] P. Grassberger, Mathematical Biosciences , 157 (1983).[3] M. Babadi, E. Demler, and M. Knap, Phys. Rev. X , 041005(2015).[4] M. C. Tsatsos, P. E. Tavares, A. Cidrim, A. R. Fritsch, M. A.Caracanhas, F. E. A. dos Santos, C. F. Barenghi, and V. S.Bagnato, Physics Reports , 1 (2016).[5] H. Kadau, M. Schmitt, M. Wenzel, C. Wink, T. Maier, I. Ferrier-Barbut, and T. Pfau, Nature , 194 (2016).[6] E. Gillman, F. Carollo, and I. Lesanovsky, Phys. Rev. Lett. ,100403 (2020).[7] M. A. Mu˜noz, R. Dickman, A. Vespignani, and S. Zapperi,Phys. Rev. E , 6175 (1999).[8] K. A. Takeuchi, M. Kuroda, H. Chat´e, and M. Sano, Phys. Rev.Lett. , 234503 (2007).[9] G. Lemoult, L. Shi, K. Avila, S. V. Jalikop, M. Avila, andB. Hof, Nature Physics , 254 (2016).[10] M. Sano and K. Tamai, Nature Physics , 249 (2016).[11] B. Nowak, D. Sexty, and T. Gasenzer, Phys. Rev. B , 020506(2011).[12] P. Bak, C. Tang, and K. Wiesenfeld, Phys. Rev. Lett. , 381(1987).[13] M. Buchhold and S. Diehl, Phys. Rev. A , 013603 (2015).[14] J. Berges, A. Rothkopf, and J. Schmidt, Phys. Rev. Lett. ,041603 (2008).[15] O. Kinouchi and M. Copelli, Nature Physics , 348 (2006).[16] N. Bertschinger and T. Natschl¨ager, Neural Computation ,1413 (2004).[17] S. Bornholdt and T. Rohlf, Phys. Rev. Lett. , 6114 (2000).[18] E. Nicklas, M. Karl, M. H¨ofer, A. Johnson, W. Muessel,H. Strobel, J. Tomkoviˇcfi, T. Gasenzer, and M. K. Oberthaler,Phys. Rev. Lett. , 245301 (2015).[19] P. Bak, C. Tang, and K. Wiesenfeld, Phys. Rev. A , 364(1988).[20] E. Altshuler and T. H. Johansen, Rev. Mod. Phys. , 471(2004).[21] S. Field, J. Witt, F. Nori, and X. Ling, Phys. Rev. Lett. , 1206(1995). [22] S.-K. Jian, S. Yin, and B. Swingle, Phys. Rev. Lett. ,170606 (2019).[23] M. J. Aschwanden, N. B. Crosby, M. Dimitropoulou, M. K.Georgoulis, S. Hergarten, J. McAteer, A. V. Milovanov, S. Mi-neshige, L. Morales, N. Nishizuka, G. Pruessner, R. Sanchez,A. S. Sharma, A. Strugarek, and V. Uritsky, Space Science Re-views , 47 (2016).[24] D. L. Turcotte, Reports on Progress in Physics , 1377 (1999).[25] G. Pruessner, Self-organised criticality: theory, models andcharacterisation (Cambridge University Press, 2012).[26] J. A. Bonachela and M. A. Mu˜noz, Journal of Statistical Me-chanics: Theory and Experiment , P09009 (2009).[27] J. A. Bonachela, S. de Franciscis, J. J. Torres, and M. A.Mu˜noz, Journal of Statistical Mechanics: Theory and Experi-ment , P02015 (2010).[28] S. Helmrich, A. Arias, G. Lochead, T. M. Wintermantel,M. Buchhold, S. Diehl, and S. Whitlock, Nature , 481(2020).[29] D.-S. Ding, H. Busche, B.-S. Shi, G.-C. Guo, and C. S. Adams,Phys. Rev. X , 021023 (2020).[30] A. Levina, J. M. Herrmann, and T. Geisel, Nature Physics ,857 (2007).[31] J. Zierenberg, J. Wilting, and V. Priesemann, Phys. Rev. X ,031018 (2018).[32] H. Schempp, G. G¨unter, M. Robert-de Saint-Vincent, C. S. Hof-mann, D. Breyel, A. Komnik, D. W. Sch¨onleber, M. G¨arttner,J. Evers, S. Whitlock, and M. Weidem¨uller, Phys. Rev. Lett. , 13002 (2014).[33] N. Malossi, M. M. Valado, S. Scotto, P. Huillery, P. Pillet,D. Ciampini, E. Arimondo, and O. Morsch, Phys. Rev. Lett. , 023006 (2014).[34] E. A. Goldschmidt, T. Boulier, R. C. Brown, S. B. Koller, J. T.Young, A. V. Gorshkov, S. L. Rolston, and J. V. Porto, Phys.Rev. Lett. , 113001 (2016).[35] S. Helmrich, A. Arias, and S. Whitlock, Phys. Rev. A ,022109 (2018).[36] The facilitation radius can be approximated as r fac = ( C / ∆ ) / ,where ∆ is the detuning and C is the van der Waals coe ffi cient.[37] T. E. Lee, H. H¨a ff ner, and M. C. Cross, Phys. Rev. Lett. ,023602 (2012).[38] I. Lesanovsky and J. P. Garrahan, Phys. Rev. Lett. , 215305(2013).[39] M. Marcuzzi, E. Levi, W. Li, J. P. Garrahan, B. Olmos, andI. Lesanovsky, New Journal of Physics , 072003 (2015),1411.7984.[40] M. Marcuzzi, M. Buchhold, S. Diehl, and I. Lesanovsky, Phys.Rev. Lett. , 245701 (2016).[41] R. L¨ow, H. Weimer, J. Nipper, J. B. Balewski, B. Butscher, H. P.B¨uchler, and T. Pfau, Journal of Physics B: Atomic, Molecularand Optical Physics , 113001 (2012).[42] R. Heidemann, U. Raitzsch, V. Bendkowsky, B. Butscher,R. L¨ow, L. Santos, and T. Pfau, Phys. Rev. Lett. , 163601(2007).[43] A. Urvoy, F. Ripka, I. Lesanovsky, D. Booth, J. P. Sha ff er,T. Pfau, and R. L¨ow, Phys. Rev. Lett. , 203002 (2015).[44] H. Weimer, R. L¨ow, T. Pfau, and H. P. B¨uchler, Phys. Rev. Lett. , 250601 (2008).[45] C. Ates, T. Pohl, T. Pattard, and J. M. Rost, Phys. Rev. Lett. ,023002 (2007).[46] T. Amthor, C. Giese, C. S. Hofmann, and M. Weidem¨uller,Phys. Rev. Lett. , 013001 (2010).[47] M. G¨arttner, K. P. Heeg, T. Gasenzer, and J. Evers, Phys. Rev.A , 043410 (2013).[48] C. Simonelli, M. M. Valado, G. Masella, L. Asteria, E. Ari- mondo, D. Ciampini, and O. Morsch, Journal of Physics B:Atomic, Molecular and Optical Physics , 154002 (2016).[49] R. Faoro, C. Simonelli, M. Archimi, G. Masella, M. M. Valado,E. Arimondo, R. Mannella, D. Ciampini, and O. Morsch, Phys.Rev. A , 030701 (2016).[50] I. Lesanovsky and J. P. Garrahan, Phys. Rev. A , 011603(2014).[51] M. M. Valado, C. Simonelli, M. D. Hoogerland, I. Lesanovsky,J. P. Garrahan, E. Arimondo, D. Ciampini, and O. Morsch,Phys. Rev. A , 040701 (2016).[52] G. Pupillo, A. Micheli, M. Boninsegni, I. Lesanovsky, andP. Zoller, Phys. Rev. Lett. , 223002 (2010).[53] D. D. Hickstein, S. T. Gibson, R. Yurchak, D. D. Das, andM. Ryazanov, Review of Scientific Instruments , 065115(2019).[54] M. Buchhold, B. Everest, M. Marcuzzi, I. Lesanovsky, andS. Diehl, Phys. Rev. B , 014308 (2017).[55] See supplementary material appended to this manuscript.[56] I. Dornic, H. Chat´e, and M. A. Mu˜noz, Phys. Rev. Lett. ,100601 (2005).[57] K. Klocke and M. Buchhold, Phys. Rev. A , 053616 (2019).[58] For simulations we take a 30 × ×
600 spatial grid and a timediscretization ∆ t = .
002 with κ = . Γ = . γ ↓ = . D = τ = × − , D T = ffi ths e ff ects on an emergent network of excitedatoms,” (2020), arXiv:2007.07697.[60] R. Wang, M. Aghigh, K. L. Marroqun, K. M. Grant, J. Sous,F. B. V. Martins, J. S. Keller, and E. R. Grant, (2020),arXiv:2006.16412.[61] V. Buend´ıa, S. di Santo, J. A. Bonachela, and M. A. Mu˜noz,(2020), arXiv:2006.03020. Appendix A: Numerical integration scheme
We simulate the Langevin equations on a discrete spatio-temporal lattice by means of an operator-splitting update scheme.Consider starting from a general stochastic di ff erential equation (SDE) over the density field ρ (cid:126) x , t like ∂ t ρ (cid:126) x , t = D ∇ ρ (cid:126) x , t + a + b ρ (cid:126) x , t + c ρ (cid:126) x , t + σ (cid:112) ρ (cid:126) x , t + d η, (A1)where the Markovian noise kernel η has unit variance and zero mean. On a lattice the Laplacian is discretized so that it getsabsorbed into the coe ffi cients a and b . Then under appropriate change of variables ρ → ρ (cid:48) ≡ ρ + d we may eliminate the constanto ff set in the noise term. The temporal update is then decomposed into two steps: a stochastic step and a deterministic step. Forthe former we drop the quadratic term, yielding an SDE of the form ∂ t ρ (cid:126) x , t = α + βρ (cid:126) x , t + σ √ ρη. (A2)This class of SDEs with a multiplicative noise kernel admits an exact solution to the corresponding Fokker-Planck equation.Denoting the current value as ρ , the distribution of values ρ after a time step δ t is given by P ( ρ ) = λ e − λ ( ρ e βδ t + ρ ) (cid:32) ρρ e βδ t (cid:33) µ/ I µ (cid:16) λ (cid:112) ρ ρ e βδ t (cid:17) . (A3)Here we have denoted λ = βσ ( e βδ t − and µ = ασ −
1. This distribution may be e ffi ciently sampled by rewriting it as a mixedGamma distribution ρ = Γ [ µ + + Poisson[ λρ e βδ t ]] /λ. (A4)The stochastic evolution of the state at time t reduces to randomly sampling from the above distribution at each time step. For d (cid:44) ρ in terms of the transformed variables by taking all ρ (cid:48) < d to d . The remaining deterministic quadratic term may be treated by a standard Euler scheme. Similarly, we treat the evolution ofthe total density n (cid:126) x , t by an Euler scheme. Appendix B: Equations of motion