Glassy dynamics in a disordered Heisenberg quantum spin system

A. Signoles*,1, 2 T. Franz*,1 R. Ferracini Alves,1 M. Gärttner,3 S. Whitlock,1, 4 G. Zürn,1 and M. Weidemüller1, 5 Physikalisches Institut, Universität Heidelberg, Im Neuenheimer Feld 226, 69120 Heidelberg, Germany. Laboratoire Charles Fabry, Institut d’Optique Graduate School, CNRS, Université Paris-Saclay, 91127 Palaiseau cedex, France Kirchhoff-Institut für Physik, Universität Heidelberg, Im Neuenheimer Feld 227, 69120 Heidelberg, Germany IPCMS (UMR 7504) and ISIS (UMR 7006), University of Strasbourg and CNRS, 67000 Strasbourg, France Hefei National Laboratory for Physical Sciences at the Microscale and Department of Modern Physics, and CAS Center for Excellence and Synergetic Innovation Center in Quantum Information and Quantum Physics, University of Science and Technology of China, Hefei, Anhui 230026, China. (Dated: January 10, 2020)


I. INTRODUCTION
The far-from-equilibrium behavior of isolated quantum systems and in particular their relaxation towards equilibrium still evades a unifying description. It has been conjectured that these systems generically relax to a state of local thermal equilibrium according to the eigenstate thermalization hypothesis (ETH) [1]. However, the ETH does not explain how the equilibrium state will be reached, or even if it will be reached in experimentally accessible timescales. This appears to be especially important in disordered quantum systems where a strong interplay between interactions and many-body entanglement can give rise to new and intrinsically non-equilibrium effects such as pre-thermalization [2][3][4], many-body localization [5][6][7], Floquet time crystals [8,9], and quantum scars [10,11].
In contrast, most natural systems (e.g., in condensed matter) are not fully isolated from their environment and hence always relax to thermal equilibrium imposed by the external bath [12]. But it is known that disorder and frustration effects can lead to a dramatic slowdown of thermalization, associated with the onset of glassy behavior [13]. A key signature of this behavior is that macroscopic observables relax in characteristically nonexponential way, as encountered for example, in doped semiconductors [14] and organic superconductors [15], quasi crystals [16], atoms in optical lattices [17] or diamond color centers [18,19]. This opens the question whether slow relaxation, which appears to be ubiquitous in open disordered systems, also emerges in iso-lated quantum systems, and if so, what are the underlying mechanisms responsible for quantum glassiness.
A prototypical system for studying far-fromequilibrium quantum dynamics is the Heisenberg XXZ Hamiltonian for spin-1/2 particles.
Compared to the Ising Hamiltonian, this system has fewer conserved quantities and shows complex, chaotic far-fromequilibrium dynamics which are difficult to describe theoretically. Here, we experimentally realize a disordered quantum Heisenberg-XXZ spin-1/2 model in an ultracold atomic gas by encoding the spin degree of freedom in two electronically excited (Rydberg) states of each atom. Spin-spin interactions arise naturally through the state-dependent dipolar interactions between Rydberg states, while disorder originates from the random positions of each atom in the gas which gives rise to distance dependent couplings [20]. Using a strong microwave field pulse that couples the two Rydberg states, we initialize the spins in a far-from-equilibrium state and probe their time evolution, thus employing our system as a quantum simulator for unitary spin dynamics out of equilibrium [10,[21][22][23][24][25].
Exploring the dynamics in a large ensemble of Rydberg spins, we observe that the magnetisation follows a sub-exponential dependence characterised by a universal exponent that is independent of the strength of disorder up to a critical value. Our experiments and supporting numerical simulations that go beyond the mean-field approximation suggest that slow dynamics described by stretched exponential decay is a generic feature of disordered quantum spin systems, hinting to-wards a unifying effective theory description.
The paper is organized as follows: In Sec. II, we give a qualitative physical picture by solving the timedependent Schrödinger equation for a few spins exactly. We then describe in Sec. III how to implement the Heisenberg-XXZ spin model in a gas of ultracold atoms that are excited to Rydberg states. In Sec. IV we experimentally characterize the relaxation dynamics, and we discuss its universal character in Sec. V. Finally, we theoretically investigate the dependence on disorder strength and character in Sec. VI.

II. QUALITATIVE PICTURE OF THE QUANTUM DYNAMICS
We consider an ensemble of spin-1/2 particles randomly positioned in space and all initialized in the |→ y, z}) refers to the spin-1/2 operator of the ith spin. The experimental protocol is illustrated in Figure 1 (a). The unitary dynamical evolution of the system is governed by the Heisenberg XXZ-Hamiltonian in the absence of magnetic fields (in units where = 1), where δ is the anisotropy parameter and J ij are the interaction couplings between the spins i and j. To remain consistent with the experimental implementation (see Section III), we focus on an anisotropy parameter δ = −0.73 and spin-spin interactions that decay as a power law J ij = C 6 /r 6 ij with the inter-particle distance r ij .
To obtain a qualitative understanding of the quantum dynamics in this system, we perform a full quantum mechanical simulation on a small ensemble of 12 spins i. In Fig. 1 (b) we show the time evolution of the magnetization x k for a single disorder realization k (grey curve). Due to the spatial disorder, spin-spin interactions give rise to complex many-body dynamics on strongly varying energy scales. This is in stark contrast to an effectively classical, mean-field prediction for this Hamiltonian [dotted line in Fig. 1(c)], which assumes each spin to evolve in the average field generated by all other spins thus neglecting quantum correlations. The initial fully magnetized state is an eigenstate of the mean-field Hamiltonian which explains the total absence of relaxation. Therefore, in the manybody case the loss of magnetization is not caused by classical dephasing, but by the build-up of entanglement between spins, witnessed by the decrease of the local purity Tr(ρ 2 i ) for each spin. Since the dynamics are unitary and therefore the full system remains pure (Tr(ρ 2 ) = 1), we can quantify entanglement by the second order Rényi entropy which increases to S (2) i = 1 on a similar timescale as the relaxation of the magnetization (see Fig. 2).
After ensemble and disorder averaging, the magnetization approaches a fully randomized state with S x,y,z = 0 [see black curve in Fig. 1(b)], consistent with the ETH prediction. However, this relaxation occurs very slowly compared to the timescales associated with spin-spin interactions, and appears to follow an almost logarithmic dependence (dashed line in Figure 2. Build-up of entanglement: Time evolution of the second order Rényi entropy of the few-particle simulation from Fig. 1(b). The plot shows the ensemble averaged entropy for a single realization (gray curve) and the disorder average over 1000 realizations (black curve). Since the full system remains pure, the Rényi entropy is a measure of entanglement, that increases on similar timescales to the maximal value of S (2) i = 1 as the magnetization relaxes to zero. Fig. 1(b)). From a few-body perspective, one may wonder whether such glassy dynamics can actually be observed in fully isolated quantum many-body systems and to what extent it shares common features with classical spin glasses.

III. REALIZING A HEISENBERG-XXZ SPIN SYSTEM WITH RYDBERG ATOMS
We address this question experimentally using a gas of ultracold rubidium atoms prepared in a superposition of two different Rydberg states: |↓ = |48s and |↑ = |49s . The electric dipole-dipole coupling between these two states realizes the Hamiltonian (1) [26], with C 6 /2π = 59 GHz · µm 6 characterizing the strength of the beyond nearest-neighbor interactions (for the derivation of the Hamiltonian, see Appendix A).
The experimental procedure [ Fig. 1 (a)] starts with a gas of 87 Rb atoms prepared in their electronic ground state |g = 5S 1/2 , F = 2, m F = 2 in an optical dipole trap and with a temperature T ∼ 50 µK, low-enough to freeze the motional degrees of freedom over the time scale of the experiment. A laser pulse of variable duration brings a controllable number of atoms N ≤ 1200 to the |↓ = 48S 1/2 , m j = +1/2 Rydberg state. For this we use a two-photon laser excitation at 780 nm and 480 nm, with a detuning ∆/2π = −100 MHz from the intermediate state |e = 5P 3/2 , F = 3, m F = 3 and an effective Rabi frequency of Ω/2π = 150 kHz. To individualy address two specific Rydberg states, including Zeeman substructure, we continuously apply a magnetic field of 6 G. To characterize the resulting threedimensional Rydberg density distribution we perform depletion imaging, where the Rydberg density is deduced by absorption imaging of the ground-state atoms before and after the Rydberg excitation laser pulse [27].
A two-photon microwave field is then used to couple the |↓ state to the |↑ = 49S 1/2 , m j = +1/2 Rydberg state. The single photon frequency ν = 35.2 GHz is detuned from the intermediate state 48P 3/2 by 170 MHz, far enough to guarantee the population in this state due to off-resonant coupling is smaller than 4%. The atoms can therefore be considered as two-level systems described by a pseudo-spin degree of freedom. In this description the microwave field acts as an external field described by the Hamiltonian with Ω the Rabi frequency, ∆ the detuning and φ the phase of the field. To vary these three parameters at the 10 ns time scale, we use frequency up-conversion with a radio-frequency field of frequency 400 MHz. The Rabi frequency is calibrated from the period of Rabi oscillations between the two spin states. To perform the Ramsey sequence shown in Fig. 1 (a), a resonant microwave π/2 pulse at an effective Rabi frequency Ω/2π = 3.00(1) MHz rotates all spins to the fully-magnetized state |→ ⊗N x with all spins pointing along the x-direction on the Bloch sphere. Uncertainties in the duration and amplitude of the pulses as well as interaction effects lead to imperfect initial spin state preparation. Based on simulations, we estimate the fidelity to be higher than 96%.
After a free evolution time t in the absence of the microwave field (Ω = 0), a second microwave π/2 pulse with adjustable phase φ is applied to rotate the equatorial magnetization components to the detection basis ( Fig. 1 (c)): In this way we effectively read out the S x and S y magnetizations from population measurements of the two spin states |↓ and |↑ using electric field ionization (see appendix B).
To ensure unitary Hamiltonian dynamics, we restrict the experimental time scales to a maximum of 10 µs which is short compared to the spontaneous decay time and redistribution by black-body radiation (113 µs and 121 µs respectively for the chosen Rydberg states [28]). To verify that the single-spin phase coherence is preserved during experimental time, we perform a Ramsey measurement with finite detuning ∆ at low spin densities where interactions can be neglected [29,30]. The full contrast oscillation shows that the single-spin phase coherence is preserved over the duration of the experiment, as shown in Fig. 3).  Figure 3.
Ramsey oscillations showing a high degree of phase coherence for a detuning of ∆/2π = 0.47 MHz for a peak spin density of ρ 0 S = 6.0(15) × 10 7 cm −3 . The solid line shows dTWA simulations for that density.

IV. EXPERIMENTAL OBSERVATION OF RELAXATION DYNAMICS
We now study the relaxation dynamics due to spinspin interactions for increasing spin densities. Fig. 4 shows the experimentally observed relaxation of the magnetization using tomographic spin-resolved readout of S x and S y . Starting from the almost fully magnetized state S x = 1/2 we observe that the magnetization decays towards the unmagnetized state within ≈ 10 µs. This is much shorter than the singlespin phase coherence time measured in Fig. 3 but still slower than the characteristic timescale governing microscopic physics, C 6 /(2π)/a 6 0 −1 = 0.7(3) µs. Here, a 0 = 4πρ 0 S /3 −1/3 is the Wigner-Seitz radius of the Gaussian spin distribution with peak density ρ 0 S [49]. The time evolution seen in the experiment is qualitatively similar to the dynamics obtained by exact diagonalisation in Fig. 1 (b). Even when accounting for imperfect preparation of the initial state S x 1/2 the mean field prediction shows essentially no relaxation on the experimentally relevant timescales (see Appendix C). This absence of dynamics at the mean-field level implies that the relaxation seen in the experiment is closely related to build-up of entanglement also apparent in the exact diagonalisation calculations.
We empirically observe that the relaxation is well described by a stretched exponential, in apparent similarity with glassy dynamics in classical disordered media [18,[31][32][33]: where β is the stretching exponent and γ J defines an effective relaxation rate. The exponent β characterizes the deviation from a simple exponential (β = 1) towards a purely logarithmic decay ( The experimental data are well described by this phenomenological function (dashed line in Fig. 4) yielding an exponent β = 0.32 (2). This value clearly rules out a pure exponential decay, i.e. β = 1, that could be expected on the basis of single-particle dephasing. It is also inconsistent with an exponent of 1/2 that is expected for an Ising model [33,34] or a fluctuator model for open quantum systems [18,19]. Our observation therefore indicates a close similarity between classical glassy dynamics and quantum glassy dynamics in the form of subexponential relaxation of the macroscopic magnetization [35].

V. EVIDENCE FOR UNIVERSAL RELAXATION DYNAMICS
To further investigate how slow relaxation and the characteristic exponent depends on microscopic details, we control the degree of spatial disorder by taking advantage of the Rydberg blockade effect in the state preparation stage [36,37]. During laser excitation the strong van der Waals interactions between Rydberg states prevents two spins to be prepared at distances smaller than the Rydberg blockade radius R bl . The degree of disorder is thus controlled by the ratio between blockade radius R bl and Wigner-Seitz radius a 0 (Fig. 5a). For a 0 R bl the blockade effect has little influence and the spins are randomly distributed, whereas the limit a 0 ≈ R bl corresponds to a strongly ordered configuration. In between the short distance cutoff imposed by the Rydberg blockade effect effectively reduces the strength of the disorder compared to fully uncorrelated random spin positions.
In the experiment we can tune the disorder strength by changing the peak spin density ρ 0 S and thus the mean number of spins per blockade sphere (a 0 /R bl ) −3 from 0.20(5) to 0.7(2) (two-dimensional representations of corresponding distributions are depicted in Fig. 5 (a)). Remarkably, we find the stretching exponent β to be almost constant over this range (inset in Figure 5b). Furthermore, after rescaling the time axis by the characteristic energy scale C 6 /a 6 0 , the time-dependent data collapse onto a single line (Fig. 5 (b)). From this we conclude that the dynamics are universal, i.e. independent of the microscopic details of the system, an observation that will be explored further in numeric simulations.
To simulate the three-dimensional Rydberg gas on a classical computer, we model the system by an ensemble of spin-1/2 particles with interactions described by the Heisenberg XXZ-Hamiltonian. Due to the large number of spins in the experiment, the full simulation of the unitary dynamics of the Hamiltonian Eq. (1) is not possible. Instead, quantum effects can be partially taken into account by applying the semiclassical discrete Truncated Wigner Approximation (dTWA, see Appendix C) [38,39] which has recently been shown to describe the dynamics of Rydberg interacting spin systems very well [20]. To model the present experiments, all relevant parameters for the simulation are determined through independent measurements, such as the spatial density distribution, total number of spins and the microwave coupling strength Ω used in the preparation and readout stages. The initial spin distribution is generated from a random excitation model of the Rydberg atoms, including a cut-off distance to account for the blockade effect (see appendix E).
The numerical simulations describe the glassy dynamics and its universal character very well (solid line in Fig. 4), further confirming the validity of the dTWA approximation for treating the dynamics of disordered quantum systems. The ability to reproduce the experimental observations using dTWA simulation also suggests a possible interpretation in analogy to classical spin glasses: quantum, instead of thermal fluctuations, allow the system to explore a rugged energy landscape leading to a hierarchical separation of timescales [35,40], resulting in stretched exponential relaxation. 0 the data collapses on a single curve described by a stretched exponential function with stretching exponent β = 0.32 (2). The inset shows β as a function of the corresponding ratio (a0/R bl ) −3 for the experimental data and discrete Truncated Wigner Approximation (dTWA) calculations (solid gray line).

VI. ROLE OF DISORDER AND BREAKDOWN OF UNIVERSALITY
Theoretical modeling using the dTWA allows to further test the role of disorder on the universal behavior, while excluding possible effects of the inhomogeneous spin density resulting from the optical trap. The simulated spin dynamics for a uniform density distribution ρ S are shown by the solid lines in Fig. 6(a). For early times where t ≤ 2πR 6 bl /C 6 , the leading order quadratic Hamiltonian evolution is clearly visible. For times beyond the perturbative regime, the relaxation is well described by a stretched exponential [dashed lines in Fig. 6(a), see Appendix D], proving that the glassy dynamics is an intrinsic many-body effect and not a result of the Gaussian spin distribution.
The fitted relaxation rate γ J does not scale with C 6 /a 6 [a = (4πρ S ) −1/3 ] but with the median of the mean field interaction strengths J mf [see Fig 6(b)], which we conclude is the characteristic energy scale in the system (see Appendix E). This scale J mf = C 6 /ã 6 defines an effective Wigner-Seitz radiusã, that coincides with the usual Wigner-Seitz radius a for small densities but deviates at larger densities when spatial correlations induced by the Rydberg blockade effect become important. Rescaling the time with C 6 /ã 6 , the simulated data collapse on a single stretched exponential curve [see Fig. 6(c)], similarly to the experimental observations in Fig. 5(b), hence confirming the universal relaxation dynamics. This is further confirmed by the fitted stretching exponent β that is shown in the inset of Fig. 6(c). We find it to be approximately constant for large disorder strengths where spatial correlations are weak (i.e. (ã/R bl ) −3 0.7), with a value of 0.36 close to the experimental one of 0.32 (2).
The simulations allow to access even higher densities than those accessible in the experiment, corresponding to more correlated spatial configurations of the atoms. The dTWA simulations shows significantly different dynamics [dotted line in Fig. 6(a)] that translates into a stretching exponent β that is sensitive to the strength of the disorder above a threshold (see inset in Fig. 6 (c)). This is linked to an increasing importance of the cutoff at large interaction strengths induced by the blockade effect which modifies the spatial distribution function. It can be seen in the probability distribution of nearest-neighbor interaction strengths, which becomes more peaked compared to a random spin distribution without blockade effect (see Fig. 7).
On the other hand, in the regime of large disorder where we observe universal, glassy dynamics, the cutoff is much beyond the median interaction strength, and the functional form of the spin distribution function is significantly modified only at larger interactions strengths. Only short time dynamics is influenced by the existence of this cutoff, and one might be tempted to conjecture that universal glassy dynamics emerges as a consequence of the interaction strengths following distribution functions of similar shape. However, this hypothesis is in stark contrast to the results of a rate model analogue to the fluctuator model [18], which assumes that each spin decays with a local rate γ i sampled from the probability distributions from Fig. 7: In this oversimplified model, one finds non-universal relaxation dynamics with a stretched exponential varying as a function of disorder strength, as shown by the dashed line the inset of Fig. 6 (c). In this model the difference in distribution functions matters across the whole range of disorder. From this we conclude that the universal stretching exponent β we observe in ex- , the numerical data collapses on one curve after rescaling time with C6/ã 6 , wherẽ a plays the role of an effective distance that takes into account the roles of disorder and power-law interactions. For densities as large as ρS = 2.11 × 10 9 cm −3 (dotted line) the spatial order introduced by the blockade is so large, that the dynamics does not follow the universal behavior observed for smaller densities. Inset: Stretching exponent β as a function of the order in the system expressed by (ã/R bl ) −3 . Below a critical value of (ã/R bl ) −3 0.7 the exponent becomes constant. The dots denote the exponent β derived from the curves depicted in Fig. 6. periment and simulation emerges from many-body effects.
For the experimental conditions the inhomogeneous density distribution leads to an average over a range of disorder strengths. However, the densities are sufficiently low that the experiments fall well within the window where a universal stretching exponent is observed. This is confirmed by additional dTWA simulations performed for the actual inhomogeneous spin distribution used in the experiments, which results in a density-independent stretched exponent β which is in excellent agreement with the measurements, as shown by the inset in Fig. 5 (b).

VII. CONCLUSION AND OUTLOOK
In this work we observed glassy dynamics in close analogy to the subexponential relaxation known from classical spin glasses. While the latter is driven by thermal fluctuations, the dynamics of the disordered isolated quantum system is governed by quantum fluctuations and spreading of entanglement going beyond mean-field approximations. Remarkably, the stretching exponent β takes on a universal value below a certain disorder strength, as confirmed by both experiment and semi-classical simulations. In the experiment, disorder is changed by controlling of the density of spins, which shifts the upper cutoff scale in the distribution of interaction strengths due the excitation blockade effect. The long time evolution is therefore insensitive to the details of the system parameters on high energy scales which hints at a low-energy effective theory valid at long times.
In addition, the observation that the dynamics of the magnetization is well described by semi-classical truncated Wigner simulations suggests that quantum interference effects become less important as the system approaches its equilibrium state. This is in line with previous findings that the long-time dynamics of generic thermalizing quantum many-body systems simplifies also in the sense that states can be represented efficiently due to limited entanglement [41].
We interpret the observed independence of the cut-off at large interaction strengths and the validity a of semiclassical description as a strong hint that the dynamics of many-body quantum systems might be amenable to a simplified description of the late-time dynamics in terms of effective low energy degrees of freedom. Concretely, this could be approached is the spirit of the strong disorder renormalization group, iteratively integrating out the highest energy degrees of freedom resulting from most strongly interacting spins or clusters of spins [42]. Furthermore, spin glasses in the aging regime have been found to show certain quasi-thermal properties. For examples a fluctuation dissipation theorem has been found to hold [43] and the spin-glass transition shows similarities to a thermal phase transition [44]. Thus, the similarities to classical glassy dynamics observed in this work are encouraging to extend such an effective thermal-like description to quantum spin glass dynamics.
In order to describe the interaction between two Rydberg excitations, the Hamiltonian is expanded in multipoles. This is well justified, as the minimal distance between the Rydberg atoms that is determined by the blockade radius R bl is much larger than the LeRoy radius R LR that describes the typical spread of the electron wave function. The leading order term of this expansion is the dipole-dipole interaction Hamiltonian that couples Rydberg atoms with different angular moment quantum number l. For dipolar forbidden transitions, the second order term in perturbation theory needs to be calculated giving rise to the van der Waals Hamiltonian where the Foerster defect ∆ F = E m − E i is the energy difference between the intermediate and initial state. Aiming for a simpler notation, the two different Rydberg states can be identified as spin states |↑ and |↓ . In the pair state basis (|↑↑ , |↑↓ , |↓↑ , |↓↓ ), the total Hamiltonian describing the interaction between two atoms i and j can be written in matrix form aŝ The additional single-spin detuning ∆ vdW is an order of magnitude smaller than the interaction strength J ij and thus negligible. The matrix elements E ↑↑ , E ↓↑ , E ↓↓ and J ex were calculated using the python module ARC [28]. For the Rydberg states 48S differential equations is solved by Tsitouras 5/4 Runge-Kutta method [46] using the Julia Differential Equations package [47]. For a perfect initial state where all spins are aligned in x-direction, mean-field theory does not predict any dynamics. However, the interactions present during the first π/2-pulse of the Ramsey protocol induce small fluctuations in the initial state. We take these imperfections into account by including the preparation and read-out pulses into the simulations which leads to the dynamics shown by the dotted line in Fig. 4. For the relevant time-scale of the experiment, these dynamics are negligible. For the considered dynamics, the initial state is an eigenstate of the mean-field Hamiltonian. The relaxation is thus triggered by the initial quantum fluctuations, meaning that mean-field approaches fails in this case. Instead, we use a discrete Truncated Wigner Approximation (dTWA), which still performs classical evolution of the spins following eq. (C2) but includes the initial quantum fluctuations into statistical ensembles of the initial state by sampling Monte-Carlo trajectories on discrete phase-space [38,39]. For the simulations in this paper, we sample over 100 initial conditions which is sufficient for the magnetization to be converged. Simulating far-from equilibrium dynamics of disordered spin systems has been successfully applied to spin systems in recent work [20]. Imperfections of the preparation and readout are taken into account by simulating the whole Ramsey sequence including those two microwave pulses. We also compared dTWA with an approximate quantum mechanical model, the so-called Moving Average Cluster Expansion [48]  tively gives similar results (see Fig. 9).

Appendix D: Quantification of slow dynamics by a stretched exponential
A phenomenological approach to describe slow dynamics in disordered systems is a fit of the magnetization with a stretched exponential with relaxation rate γ J and stretched exponent β. This was already proposed by Kohlrausch in 1847 [51], a review on the stretched exponent in numerical simulations and in experimental data of various materials can be found in [31].
For β = 1 the stretched exponential describes an exponential decay. In the limit β → 0, the stretched exponential approaches the logarithmic decay which can be seen by performing a Taylor expansion at small β  is an eigenstate of S x , the expectation value of the commutator [H, S x ] vanishes for this state and we expect the initial dynamics to be quadratic in time. However, this doesn't hold for the stretched exponential function that is a power low with exponent β for short times t Therefore, we exclude the very early dynamics from the fit where t < 1/J max (see Fig. 6).

Appendix E: Spatial spin distribution and disorder strength
To model the experimental 3D spin distribution, we simulate the Rydberg excitation dynamics in a cloud of ground state atoms with Gaussian density distribution (measured radii at 1/e 2 : σ x = 203(3) µm, σ y = σ z = 35(1) µm). We excite each atom randomly one after the other with a certain probability if it is not within the blockade sphere of an already excited atom. The excitation probability takes into account the √ N b enhancement due to Rydberg interactions, with N b the number of atoms within a blockade sphere. It also consider the profile of the laser excitation, characterized by a Gaussian distribution of the two-photon Rabi frequency with measured radius σ = 70.6(3) µm (1/e 2 ). Its amplitude is adjusted such that the total number of excited atoms equals the one measured by field ionization.
For simulations of a homogeneous system, the spins are randomly distributed in a uniform box taking into account a blockade radius of 5 µm until the desired densityρ S is reached. In the limit of no blockade effect, the nearest-neighbor distribution for a given Wigner-Seitz radius a = (4πρ/3) −1/3 would be given by [49] h(r) = 3 a (r/a) yielding the distribution of coupling strengths h(J) = h(r) ∂r ∂J . Instead, the blockade effect modifies the nearest-neighbour distribution, resulting in a different coupling distribution g(J) (see Fig. 7 (a)). We quantify the disorder strength of the spin system with the Kullback-Leibler divergence [50] D KL (g h) = g(J) log g(J) h(J) dJ, i.e. the amount of information that is gained by updating from the distribution h(J). Indeed, the Kulback-Leibler divergence increases almost linearly with density (see Fig. 10. This confirms that (a/R bl ) −3 is a relevant scale to describe the disorder strength.
When investigating the universal character of the spin dynamics for homogeneous system, we have concluded that the typical energy scale should be determined by C 6 /ã 6 (see Fig. 6). The effective distanceã is defined as the median of the non-Gaussian distribution {r i } with r −6 i = j r −6 ij . The energy scale C 6 /ã 6 is thus the median of the mean-field energies C 6 /r 6 i .
Appendix F: Summary of parameters Table II summarizes the parameters of the individual measurements shown in figure 2 and 3 of the main text.