Methods
Sample preparation
Our experiments are performed using a thermal gas of potassium-39
atoms, loaded directly from a magneto-optical trap into a crossed
optical dipole trap. The resulting cigar shaped atom cloud has a temper-
ature of 40 μK and e−1/2 radii of 10 μm × 100 μm. This should be compared
to the characteristic distance between facilitated Rydberg excitations
rfac = (C 6 /Δ)1/6, which for a detuning of Δ/2π = 30 MHz is about 1.7 μm.
The peak number of atoms in the |g⟩ state is 1.3 × 10^5 , and the density
determined by in situ imaging is 2.4 × 10^11 cm−3. To vary the density while
holding all other parameters fixed, we reduce the magneto-optical trap
loading time. The lifetime of the atoms in the trap without Rydberg
excitation is about 4 s—that is, much longer than the relevant timescales
for the SOC dynamics.
Excitation laser
To excite the atoms to the 39p3/2 Rydberg state we use a single photon
optical transition at a wavelength of 285 nm. This light is produced by a
frequency-doubled dye laser delivering up to 200 mW of single-mode
light and is frequency stabilized to a high-finesse cavity, resulting in an
independently measured linewidth of 400 kHz. The excitation beam is
aligned parallel to the long axis of the trap and weakly focused to a waist
much larger than the size of the atom cloud such that it is practically
uniform. We experimentally determine the Rabi frequency Ω for every
individual repetition of the experiment by logging the respective single-
shot-laser power on a photodiode and employing an independent Rabi
frequency calibration based on measuring the light shifts induced by
the laser via Ramsey interferometry^51.
Numerical simulation of the Langevin equation
Although the Langevin equation (equations ( 2 ) and (3)) is straightfor-
ward to solve in the mean-field approximation, in Fig. 1 we show exem-
plary numerical simulations that capture the effects of diffusion and
multiplicative noise terms in a one-dimensional setting. For these
simulations we make use of the XMDS2 (stochastic) differential-equa-
tion solver package^52 , assuming a transverse grid size of 128 points
and a timestep of 2.5 × 10−3. The noise term is implemented as a zero-
mean Wiener process with a standard deviation proportional to ρ.
However, to ensure numerical stability we found it necessary to impose
a noise cutoff by setting ξ = 0 when ρ < 0.0025n 0. For b = 0 the solutions
exhibit an absorbing-state phase transition at n 0 = 0.39 and power-law
scaling consistent with directed percolation universality (in one-dimen-
sion, βDP = 0.276). For b ≠ 0 we find that the individual timetraces
obtained from the full numerical solution are qualitatively very similar
to the corresponding mean-field solutions. By fitting the numerical
results in the same manner as performed for the experimental data,
we obtain slightly larger effective parameters κ and Γ.
Comparison of the power-law hypothesis to alternative
distributions
To test whether the avalanche data are indeed described by a power-law
distribution we employ the widely used Kolmogorov–Smirnov (KS) test
against several alternative distributions, including other heavy-tailed
distributions (following the definitions in ref.^6 ). The KS statistic is
defined as the maximum distance between the cumulative distribu-
tion of the empirical data and that of the hypothesized distribution,
with small values much less than 1 indicating good agreement. In all
cases we minimize the KS statistic as a function of the parameters of
the hypothesized distribution, restricting the data and the hypoth-
esized distributions to the range 20 ≤ s ≤ 400. For the data depicted
in Fig. 4 , the obtained KS-test statistics are: 0.015 (power law), 0.102
(exponential), 0.031 (log-normal), and 0.04 (gamma). This shows that
the power-law distribution provides a better fit to the data than the
alternative distributions. The power-law exponent α = −1.38 found via
KS minimization is in excellent agreement with the value obtained via
the maximum-likelihood estimation^53.Detuning dependence and further evidence for non-equilibrium
universality
In the following we present additional evidence for the universal nature
of the self-organized stationary state. For this we performed additional
measurements of the stationary density as a function of the driving
intensity but for different detunings of the excitation laser, as shown in
Extended Data Fig. 1. Each dataset shows qualitatively similar behaviour
to that presented in Fig. 3 , clearly showing the transition from an absorb-
ing phase to a self-organizing active phase. However, these data also show
that the location of the critical point depends on the laser detuning.
To further analyse these data we apply the scaling ansatz
nf0=(nFΩΔ^2 dn1/ 0 β′), where β′ = 0.869, and we have included as a new
parameter the detuning rescaling exponent, d. For d = −2.06(1) the data
again collapse onto a single universal curve. In this way we determine
the κ ∝ Ω^2 /Δ2.06 dependence of the spreading parameter, used elsewhere
in the paper to compare the data with mean-field theory.
Before analysing the scaling properties of the rescaled data, care-
ful inspection shows that it has a slightly different form to the scaling
function F(x) used to describe the data in Fig. 3b. This is evidenced by
the fit to F(x), shown as a blue dashed line in Extended Data Fig. 1b. The
deviation is most apparent in the normalized fit residuals (Extended
Data Fig. 1, inset) which, in contrast to Fig. 3b, exhibits some structure
(for example, the inverted U-shape of the black datapoints). Unless
properly accounted for, this deviation between the scaling form of
the data and the heuristic scaling function causes a systematic error in
the determination of the critical scaling exponent. To rectify this, we
model the detuning-dependent data by a generalized scaling function
F′(x) = [1 + (x/xa)vα + (x/xc)vβ]−1/v, where the newly introduced parameters
xa < xc and α < β empirically describe power-law scaling for intermediate
driving intensities. In the asymptotic regime x ≫ xc, the scaling function
once again reduces to a power law nf/n 0 ∝ x−β.Critical response
As additional evidence for the system reaching a critical state, we
have investigated the gapless response of the stationary state follow-
ing a parameter quench. Assuming the SOC state is indeed an attrac-
tor for the dynamics, on one hand we expect that small perturbations
(for example a sudden change of the spreading parameter κ),
should trigger avalanche-like processes that eventually bring the
system back to a new critical state corresponding to a lower station-
ary density. On the other hand, if the system evolves to a state that is
deep within the absorbing phase, then avalanches can only be trig-
gered by perturbations larger than a threshold value corresponding
to a non-zero dynamical gap. To measure this response we start from
the stationary state (reached after t = 10 ms) corresponding to differ-
ent driving intensities Ωi^2 (sketched in Extended Data Fig. 2a). We then
perturb the system by quenching the driving intensity to a new
value Ωf1^2 and then wait for a further 10 ms before measuring the
new stationary density. The whole procedure is then repeated
with a slightly larger final driving intensity Ω^2 f2≈+Ωπf1^22 (2 ×50kHz).
From these two measurements we estimate the susceptibility
χ=dnΩff/d^2 =[nΩff()^21 −(nΩff^22 )]/(ΩΩf1^2 −)f2^2.
Extended Data Fig. 2 shows the measured susceptibility as a function
of δ=(ΩΩf^2 −)i^2 /Ωc^2 for three different initial conditions corresponding
to Ωi < Ωc (absorbing), Ωi ≈ Ωc (critical) and Ωi > Ωc (active). For each of
these initial conditions we observe pronounced minima in χ corre-
sponding to the strongest system response. We interpret the leading
edge on the left side of each minimum as the point where the perturba-
tion is sufficient to bring the system back to the active phase, thereby
triggering avalanche-like dynamics and extra loss. When starting deep
in the absorbing phase (black circles) the onset occurs at a large value
of δ, which is a measure of the non-zero dynamical gap. By contrast,