{'fraction_non_alphanumeric': 0.03653772079458537, 'fraction_numerical': 0.013078009501274785, 'mean_word_length': 4.41679085941381, 'pattern_counts': {'":': 0, '<': 14, '<?xml version=': 0, '>': 0, 'https://': 0, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 5, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
Slow Thermalization of Few-Body Dipole-Dipole Interactions

We simulate the dynamics of Rydberg atoms resonantly exchanging energy via two-, three-, and four-body dipole-dipole interactions in a one-dimensional array. Using a simplified model of a realistic experimental system, we study the initial state survival probability, the level spacing statistics, the spread of entanglement, and the properties of the energy eigenstates. By exploring a range of disorders and interaction strengths, we find regions in parameter space where the three- and four-body dynamics exhibit nonergodic behavior and either fail to reach thermodynamic equilibrium or do so slowly. The interplay between the always-resonant and field-tuned interactions gives rise to quantum many-body scar states, which play a critical role in slowing the dynamics of the three- and four-body interactions.


I. INTRODUCTION
The thermodynamic evolution of strongly interacting, closed quantum systems offers insight into behavior that cannot be described by conventional statistical mechanics. Unlike classical systems, which explore all of their states with equal probability in the thermodynamic limit, the evolution of an isolated quantum system is unitary. The eigenstate thermalization hypothesis (ETH) proposes that if a quantum system thermalizes for any initial state, then it must likewise reach thermal equilibrium when initialized in any of the many-body energy eigenstates [1][2][3]. For a system that obeys the ETH, entanglement can spread rapidly while the memory of the initial state becomes inaccessible to local measurements. In the presence of disorder, an alternative fate is many-body localization (MBL), in which some signature of the initial state remains at long times [4][5][6]. Nonergodic dynamics are also possible if the initial state has large overlap with quantum many-body scar states, which form a small subset of energy eigenstates disconnected from the rest of the Hilbert space [7]. Understanding these dynamics will provide insight into the fundamental physics of energy transport, localization, collective behavior, and entanglement growth. Furthermore, control of closed quantum systems is necessary for engineering quantum devices and for the analog modeling of condensed matter.
The development of experimental platforms using superconducting qubits [31][32][33][34], trapped ions [35,36], nitrogen-vacancy centers [37], and ultracold atoms [38][39][40][41][42][43][44][45] has made it possible to explore these dynamics in the laboratory. The ability to precisely control their geometry and interactions makes ultracold Rydberg atoms an excellent candidate for studying the quantum dynamics of closed systems. Systems of ultracold Rydberg atoms have been used to measure the thermalization of a dipolar spin system [38], observe a nonergodic extended phase [39], explore the possibility of MBL in an ultracold plasma [40,41], measure slow relaxation in a Heisenberg XXZ spin system [42], measure the spread of entanglement in a many-body localized system [43], and probe the localization-delocalization transition in a three-dimensional system with random hopping interactions [44]. Nonergodic dynamics have been observed in an interacting Rydberg atom array [45] and quantum many-body scar states have been proposed as the mechanism [7].
Resonant Förster, or dipole-dipole, interactions among ultracold Rydberg atoms have been studied extensively [46][47][48][49][50][51][52][53][54][55][56]. In a two-body interaction, the initial and final atomic pair states are degenerate and coupled by the dipole-dipole operatorV dd ∝ 1/R 3 . If both final electronic states are different from both initial electronic states, then the interaction is typically tuned into resonance by using an electric field to adjust the Stark energy levels. In rubidium, for example, the 32d + 32d ↔ 34p 3/2 + 30g exchange is resonant at about 0.5 V/cm [48,52]. There are also the alwaysresonant interactions, in which the electronic states are swapped between initial and final states. For example, the 39p 3/2 + 39s 1/2 ↔ 39s 1/2 + 39p 3/2 exchange is always resonant in a spatially homogeneous field [57]. Alwaysresonant interactions are often referred to as hopping interactions, since one can imagine the upper energy state to be an excitation that hops from atom to atom. This type of interaction maps to the Heisenberg XY model and thus has been well-studied [38,44,58,59].
While some recent work has shown that localiza-tion can persist in the presence of three-body interactions [60,61], most studies of thermalization to date focus on two-body interactions. Three-and four-body interactions have been observed in both cesium and rubidium [62][63][64][65]. Three-body interactions have been considered for use in quantum gates [66,67] and a new class of three-body interactions in rubidium has been proposed [68]. The recent observation of few-body dipoledipole interactions among Rydberg atoms presents the opportunity for further study and is the subject of our numerical investigation.
Recently, Liu et al., including some of the authors of this work, measured the time evolution of two-, three-, and four-body interactions in rubidium for the maximally disordered case of a frozen Rydberg gas [65]. The experiment was a typical quantum quench in which the atoms were excited to the same initial Rydberg state. An electric field then tuned dipole-dipole interactions into resonance, which coupled the atoms to nearby final Rydberg states. The fraction of atoms in the upper final Rydberg state was measured and appeared to grow more slowly than expected for the three-and four-body interactions.
Motivated by these experimental results, we study the dynamics of a one-dimensional array of Rydberg atoms that exchange energy via these few-body interactions. We do so by solving the Schrödinger equation using exact diagonalization. We have developed a simplified model that allows us to simulate the dynamics to longer times, which is necessary to determine if the system eventually thermalizes. While our studies are prompted by ultracold Rydberg experiments, they also present an in-depth analysis of thermalization in a few-body interacting system that is more broadly applicable.
We start by examining the dynamics and eigenstates in an ordered system. We then add disorder to the system and study the effect of varying both disorder and interaction strength. We numerically determine the decay of the initial state survival probability, the growth of entanglement entropy, and the level spacing statistics. We graph the results for the survival probability and the level spacing statistics as a function of disorder and interaction strength, which shows regions of nonergodic evolution.
While the two-body interaction thermalizes rapidly in all of our simulations, our results reveal slow thermalization and nonergodic evolution for the three-and fourbody interactions. In those cases, we find numerical evidence for quantum many-body scar states. The scar states arise when the hopping interactions are significantly stronger than the field-tuned interactions, as can be the case for the three-and four-body resonances. Furthermore, increasing the disorder significantly slows the dynamics and suggests that these systems could be candidates for MBL studies.  Figure 1. The two-, three-, and four-body interactions for (a) a realistic four-level experimental system and (b) the simplified three-level model presented in the text. The horizontal axis is the applied electric field while the vertical axis is energy. The pairs of solid red, dashed blue, and dotted green arrows represent individual detuned steps in the three-and four-body interactions. Open circles indicate initial states and closed circles indicate final states. While these circles are horizontally offset to allow space for the transition arrows, the individual resonant interactions occur at only one electric field. A more in-depth comparison of the two models is presented in the appendix.

II. MODEL
Rubidium Rydberg atoms in np states with 32 < n < 38 can resonantly exchange energy through a set of few-body dipole-dipole interactions [64,65]. An applied electric field tunes these interactions into resonance via the Stark effect. For example, Liu et al. initially excited a sample of atoms in a magneto-optical trap to the 36p 3/2,|mj|=1/2 state, which we label as the p state. The energetically nearby 36p 3/2,|mj|=3/2 , or p ′ state, was initially unpopulated. As shown in Fig. 1(a) by the solid black arrows, two p atoms can resonantly exchange energy by the field-tuned interaction where s refers to the 37s state and s ′ refers to the 36s state. At a slightly larger electric field, the two-body interaction is tuned out of resonance by the energy E p −E p ′ as shown by the dotted red arrows in Fig. 1(a). A third p atom can be recruited to account for this defect via the s + p → p ′ + s interaction so that one atom ends up in the p ′ state, as shown by the dashed blue arrows in Fig. 1(a). By further increasing the electric field, the energy defect can be set to 2(E p − E p ′ ) so that a fourth p atom is needed to bring the interaction into resonance. These three-and four-body interactions are Borromean in nature and require all atoms to participate simultaneously [64]. The simulation results presented in Liu et al. used exact diagonalization to solve for the time evolution, based on a model that included the four energy levels s, p, p ′ , and s [65]. The number of states N for this Hamiltonian increases rapidly with the number of Rydberg atoms such that the simulations were limited to nine atoms, since at ten atoms N = 184 756. This presents a challenge to studying the dynamics at longer times. The timescale on which we expect the finite size of the system to become relevant is given by the Heisenberg time t H = 2π/∆E. We estimate ∆E from the average energy level spacing and find that t H is a few 10s of µs for nine atoms. While this is sufficient for comparison of the simulation results to the experimental data in [65], we cannot numerically study thermalization using this model.
To address this problem, we have developed a simplified model that involves only three energy levels: an initial state p, an upper final state s, and a lower final state s ′ . This model presents a similar set of field-tuned resonant two-, three-, and four-body interactions, as shown in Fig. 1(b), given by The three-and four-body interactions require s ↔ s ′ transitions which are forbidden by angular momentum selection rules in a real atomic system. In our model, we allow these transitions, resulting in three possible alwaysresonant, or hopping, interactions We find that our simplified model shows essentially the same dynamics as a model including all four energy levels. See the appendix for a more detailed comparison of the two models. In our three-level model, the three-and four-body matrix elements can be calculated perturbatively by considering a series of detuned two-body interactions, yielding similar results to [65]. For example, in the three-body energy exchange of Eq. (3) shown in Fig. 1(b), the initial two-body interaction p+p → s+s ′ is detuned by E p −E s ′ . The subsequent two-body interaction s ′ + s → s + s ′ accounts for this detuning and brings the three-body interaction into resonance. Similarly, the four-body interaction is detuned by 2(E p − E s ′ ) and can be made resonant with two additional two-body exchanges.
In our simulations, we solve the Schrödinger equation by exact diagonalization to obtain the time evolution of the initial state. Our Hamiltonian includes only the fewbody dipole-dipole interactions. We ignore the angular dependence and spontaneous decay. For example, the three-body Hamiltonian iŝ whereσ i ps is an operator that takes the i th atom from the p to the s state and the sums are performed over distinct triplets or pairs of atoms. The second term represents the two-body hopping interactions while the first term is the field-tuned resonant three-body interaction with where R ij is the distance between atoms i and j. The detuning from the two-body resonance is δ and µ, ν, and γ are the dipole matrix elements coupling p → s, p → s ′ , and s → s ′ , respectively. In Eq. (8) we have suppressed the diagonal matrix elements as they are all zero since we perform all simulations on resonance. The two-and fourbody Hamiltonian operators are shown in the appendix. The dependence on the distance between atoms in Eq. (9) results from summing over all possible paths from the initial state to the final state. While the two-body field-tuned resonant and hopping matrix elements are proportional to r −3 , the field-tuned resonant matrix elements of the three-and four-body interactions effectively scale as r −6 and r −9 , respectively.
At the Rydberg atom densities studied in [65], typical values for the two-body matrix elements were a few hundreds of kHz and for the three-and four-body matrix elements a few kHz. In this work, we simulate a range of values for the matrix elements. From the largest to the smallest linear array spacings, the three-and four-body matrix elements range from a few kHz to a few hundred kHz.
The Hamiltonian matrix for the new model grows less rapidly, enabling us to include more atoms in the simulation. We can simulate up to 12 atoms for two-body interactions (N = 73 789), 13 atoms for three-body interactions (N = 93 094), and 14 atoms for four-body interactions (N = 108 109). Since the Heisenberg time increases exponentially with the size of the system, we find that it is generally greater than a few hundred µs in our simulations. Where possible in the discussion of our results, we display t H directly on the plots.
Panda et al. have recently argued that the system sizes in current numerical studies are generally too small to model the transition from thermalization to localization [69]. The metrics we use to measure the dynamics of our system generally stabilize well before t H . However, we exercise caution when drawing definitive conclusions about the presence of MBL in our results.
In all of our simulations, we arrange the atoms in a one-dimensional array, with the spacing between atoms d chosen to be experimentally realistic. We vary this spacing to adjust the typical strength of the interactions. Results are presented using atomic units for energy and, in order to compare the different few-body interactions, we use a natural time unit. The time unit is the reciprocal of the matrix element for N nearest neighbors interacting via the N -body interaction that are separated by the nominal spacing for that simulation run. The values of µ, ν, γ, and δ are chosen to be similar to the real experimental values. We assume that the atoms are stationary during the time of the simulation.
All simulations are run on a supercomputer with four Nvidia A100 graphical processing units operating in parallel. Linear algebra operations are performed with the MAGMA software package [70][71][72].

III. RESULTS
In order to study the thermalization of few-body dipole-dipole interactions among Rydberg atoms, we examine some relevant dynamical variables in our simulations. We calculate the fraction of atoms excited to the s state as a function of time, which has the advantage of being straightforward to measure in experiment. We also determine the survival probability of the initial state, the energy level spacing statistics, and the entanglement entropy. All three of these are used widely in numerical studies to distinguish between thermalizing and nonergodic systems. For ordered linear arrays of atoms, we also study the properties of the energy eigenstates.

A. Ordered System
We start by studying the ordered system, where we can gain insight by examining the many-body eigenstates. When the initial state is composed entirely of p atoms, we find that the dynamics for three-and four-body interactions are significantly slower than for two-body interactions. This is true even when the field-tuned matrix elements are similar in magnitude. For example, the three-body matrix element of a triplet of nearestneighbor atoms with an array spacing of 11 µm is about the same as the two-body matrix element of a pair of nearest-neighbor atoms with a spacing of 50 µm.
The fraction of atoms that end up in the s state as a function of time is displayed in Fig. 2(a), with the twobody in solid black and the three-body in dotted blue. For comparison, the data are displayed as a fraction of the expected saturation level. The saturation level is calculated under the assumption that the ETH is obeyed and all many-body eigenstates are thermodynamically equivalent. The two-body saturation s fraction is 0.326, the three-body is 0.213, and the four-body is 0.153. The twobody s fraction reaches saturation quickly compared to the three-body s fraction, which eventually approaches saturation after t = 10.
The survival probability or fidelity, F (t), is defined as and it gives the probability of finding the system in its initial state |Ψ(t = 0) at some later time t. In Fig. 2(b), the two-body survival probability decays rapidly to zero and remains there. The three-body survival probability decays more slowly and exhibits periodic partial revivals, indicating that some coherence is preserved in this system. We also compare the entanglement entropy (EE) for each system. We divide the linear array into two subsystems, the left and right halves, and calculate the bipartite Von Neumann EE, S(t), which measures the degree of entanglement between the two subsystems. The EE is S(t) = − i λ i ln λ i where λ i are the eigenvalues of the reduced density matrix of one half of the system. The EE are shown in Fig. 2(c) as a fraction of the maximum value.
The EE has been found to grow only logarithmically in the MBL phase for short-range systems [12,13,73], while it grows algebraically for the MBL phase in long-range systems [43,58,74]. The interactions 1/r β are shortrange if β is greater than the dimension of the system, as in our case. We see in Fig. 2(c) that the EE for the twobody field-tuned interactions grows exponentially before beginning to plateau near its maximum value. The EE for the three-body field-tuned interactions, after an initial rapid increase, grows only logarithmically with time. The initial rapid increase is likely due to the fact that the hopping interactions are stronger in the more closelyspaced three-body system.
An essential difference between the two-and threebody cases shown in Fig. 2 is the relationship between the always-resonant hopping and the field-tuned interac-tions. In the two-body case, these matrix elements are approximately equal and both scale as r −3 . However, in the three-body case the field-tuned matrix element scales as r −6 and is thus significantly smaller than the hopping matrix element. We can numerically test the effect of the hopping interactions in the three-body case by artificially disabling them in our simulation. The results are shown by the dashed red line in Fig. 2.
In Fig. 2(a), we see that the artificially modified threebody field-tuned interaction now approaches saturation more quickly, similar to the two-body dynamics. In Fig. 2(b), the survival probability collapses rapidly like that of the two-body case. The EE, on the other hand, still shows logarithmic growth but no longer outpaces the two-body EE at early times. We can conclude that the always-resonant hopping interaction is an important factor in the relatively slow dynamics of the three-body case.
One might not expect the presence of hopping interactions to slow the dynamics. However, Celardo et al. have also studied a model in which hopping interactions lead to localized features [15]. They add long-range hopping interactions to a one-dimensional Anderson model [75] and find that an energy gap forms between the ground state and excited states. Beyond a critical value of the hopping strength, the size of the gap grows with the hopping strength. When the gap is present, this results in "cooperative shielding" in which certain eigensubspaces evolve as if there were no hopping interactions [17]. Additionally, localization has been found in one-dimensional systems when the hopping interactions are correlated [23,76].
In order to study the strong hopping interactions apart from any effects due to the three-or four-body nature of the field-tuned interactions, we construct a two-body model in which we artificially weaken the field-tuned interactions. The two-body field-tuned matrix element is then αµν/R 3 , where 0 < α ≤ 1. The alwaysresonant matrix elements are either µ 2 /R 3 , ν 2 /R 3 , or γ 2 /R 3 . Excluding the last, which represents the unphysical s + s ′ ↔ s ′ + s hopping, does not have a significant effect on the results.
We start by running our two-body model for only three atoms so that we can calculate the full angular dependence. To do so, we include the m j sublevels and modify the matrix elements to add the appropriate angular prefactors, depending on the total ∆m j for each transition [77]. Figure 3 shows the mean probability, P i , that the system remains in the initial state as a function of 1/α, where α is the numerical factor that modifies the strength of the field-tuned interaction. The mean probability approaches one for values of α that are typical for the three-body interaction at experimentally relevant spacings of around 10 µm. We have included as many as six atoms in this model and we find the angular dependence does not play a significant role in our onedimensional system. The effect of the anisotropy on the thermalization of two-and three-dimensional systems has 0.8 Figure 3. The mean probability Pi of remaining in the initial state as a function of 1/α plotted on a semi-log scale, where α is a numerical factor that modifies the strength of the fieldtuned interaction. The probability was averaged over many Rabi periods at each value of α for a three-atom model including angular dependence.
We can gain further insight by extending our two-body model to larger systems, where we again ignore the angular dependence, and examining the many-body eigenstates. The local density of states (LDOS) is the overlap of each eigenstate with the initial state as a function of energy, given by where E i are the energy eigenvalues, |ψ i are the eigenstates, and the delta function selects the energy. For the unmodified case where α = 1, Fig. 4(a) shows the LDOS, which is binned by energy. The LDOS fits well to a Gaussian, shown in dashed red, which implies rapid thermalization and collapse of the survival probability [19]. We also graph log 10 of the overlap as function of energy on a scatter plot as shown in Fig. 4(b). In this figure, darker colors indicate a lower density of points while lighter colors indicate a higher density. For the unmodified two-body case where α = 1, no eigenstate or set of eigenstates dominates and all of the overlaps are smaller than 10 −2 . The initial state is delocalized and ergodically spread across the energy eigenbasis.
Figures 4(c) and 4(d) show the results when α ≪ 1. Here, the field-tuned interactions are 60 times weaker than the always-resonant interactions, which is a typical ratio for the three-body interaction. The LDOS is extremely narrow compared to the unmodified case in Fig. 4(a) (note the zoomed in horizontal scale in Fig. 4(c)). The initial state has significant overlap with very few eigenstates. A sparse LDOS can indicate nonergodic behavior [19].
In the scatter plot of Fig. 4(d), there are a few eigenstates near zero energy with large overlap. There are also eigenstates, in bands of energies to the left and right of zero, that have atypically high overlap for that energy. Special eigenstates like these, dubbed quantum manybody scar states by Turner et al. [7], have recently been discovered in a number of model systems [81][82][83][84][85][86][87][88] and proposed as the mechanism behind observed revivals in  In (c) and (d), the field-tuned interaction matrix elements have been decreased by a factor of 60 relative to the hopping interactions, as is typical for the three-body case. In (e) and (f), the value of α is the same as in (c) and (d) but we start in the alternative initial state |ss ′ p · · · p .
a Rydberg atom lattice [7,45]. These scar states form a subset that has minimal overlap with the rest of the Hilbert space and thus fail to satisfy the ETH. If the initial state strongly overlaps one of these scar states, the subsequent time evolution features revivals as the system oscillates among this small subset. Indeed, we see such oscillations in the three-body simulation results of Fig. 2.
The appearance of a special class of eigenstates in our system is broadly similar to another interference phenomenon that arises in a strongly interacting Rydberg system. In the dipole blockade, the optical excitation of more than one Rydberg atom is suppressed by the addition of strong dipole-dipole interactions among the atoms [47,49,89,90]. The energies of states with multiple Rydberg atoms are shifted out of resonance with the excitation laser. These strong interactions effectively decouple the initial state from any multiple excitation eigenstates, leaving only couplings to singly excited states.
Interestingly, Iñarrea et al. have recently studied the analogous classical system of a chain of interacting dipoles [91]. They perturb their system by adding a small amount of energy ∆K to one of the dipoles and simulate the subsequent evolution in time. They find that for certain initial conditions and values of ∆K, the chain thermalizes very slowly due to "chaotic breather" modes. These modes persist for long times and keep most of the system's energy localized in a few dipoles.
The nonergodic behavior visible in Figs. 2 and 4 depends on the high overlap of the initial state with the scar states. Instead of the initial state |pp · · · p used in Fig. 4(a)-(d), consider the initial state |ss ′ p · · · p , in which one nearest-neighbor pair of atoms starts in ss ′ . The LDOS shown in Fig. 4(e), while not Gaussian, is only slightly narrower than in the unmodified case of Fig. 4(a). The scatter plot of Fig. 4(f) does not show evidence of scar states. The survival probability collapses to near zero and the EE grows more rapidly, confirming that we see faster thermalization for this initial state.
Leaving our two-body model, we can visualize the eigenstates for the three-and four-body interactions in the same way. A scatter plot of the overlap as a function of energy for the three-body example of Fig. 2 reveals similar quantum many-body scar states, as shown in Fig. 5(a). In Fig. 5(b), we show the results for a linear array with a 7 µm spacing undergoing four-body fieldtuned interactions, for which the scar states are even more apparent. The field-tuned matrix elements are similar in magnitude for these two cases, but the hopping interactions are significantly stronger for the more closely spaced atoms in the four-body simulation. The emergence of quantum many-body scar states in the threeand four-body interactions is the primary driver of their slow dynamics.

B. Disordered System
The presence of disorder in low-dimensional systems diminishes spin and energy transport, effectively driving the system towards the MBL phase [4,5,43,92]. Disorder can be added to ultracold atomic lattices, for example, by superimposing an incommensurate lattice [43]. The frozen Rydberg gas which motivates our work, on the other hand, is already disordered by the random po- sitions of the atoms. Disorder is added to our one-dimensional model by shifting the position of the i th atom by w i d along the array axis, where w i is a random number drawn from the uniform distribution [−w, w]. We simulate our system for five disorders ranging from w = 0.05 to w = 0.45 in increments of 0.10. For the three-body case, we examine seven linear array spacings ranging from 7 µm to 13 µm, and for the four-body case, we survey five spacings ranging from 4 µm to 8 µm. In both cases, the range of d was chosen to demonstrate the entire range of behaviors while also being experimentally realizable. For each case we average over at least 100 samples.
The archetypal model for studying thermalization is the disordered interacting spin chain, in which the spins are fixed on a one-dimensional lattice and the disorder is introduced via random on-site potentials. In such models, the disorder and the interaction strength can be adjusted independently. However, in this model, disorder and interaction strength are coupled, since perturbing the positions of the atoms necessarily changes the dipoledipole matrix elements. This is an important feature of the amorphous Rydberg systems we model.
We use four metrics to characterize the dynamics: algebraic fits of the fidelity decay profile at long times, s saturation, level spacing statistics, and the entanglement entropy. We observe slow thermalization for the threeand four-body interactions with high disorder and with large array spacings. As d increases, the ratio of fieldtuned to hopping matrix elements decreases rapidly. Like that of the zero disorder case, this results in scar states that slow the dynamics. The addition of disorder further slows the dynamics and, for systems with high disorder, we observe very slow thermalization or even failure to thermalize.
We first examine the long-time behavior of the fidelity. the onset of thermalization in isolated many-body systems based on fits of the decay profile of F (t) at long times [19]. They apply their method to an interacting spin-1 2 model in a one-dimensional lattice. After a rapid exponential collapse, the fidelity is best described by algebraic behavior F (t) ∝ t −γ . The value of the exponent γ, which is determined by the shape and filling of the LDOS, is shown to predict the onset of thermalization.
Távora et al. identify three regimes of behavior depending on the value of γ. Decay profiles consistent with γ ≥ 2 indicate a delocalized phase, while profiles matching 0 < γ < 1 reflect a lack of ergodicity. For the delocalized phase, thermalization is guaranteed as the LDOS is ergodically-filled, implying that the eigenstates are not correlated and the initial state is strongly delocalized. For the nonergodic phase, the eigenstates are correlated as a consequence of a sparse LDOS. They conclude that the intermediate regime, for which 1 ≤ γ < 2, requires further investigation.
Applying the criteria developed by Távora et al. to our disordered system, we examine a typical algebraic decay fit of the fidelity as a function of time in Fig. 6. In this example, the system is undergoing three-body interactions with w = 0.45 and d = 9 µm. The horizontal axes show time in both natural units and in µs. The initial rapid exponential collapse of the survival probability is shown in blue, while the long-time power-law decay regime is shown in red. In this example, the power-law decay fit yields γ = 0.2 which is consistent with slow thermalization. The power-law behavior begins well before the Heisenberg time t H , shown by the vertical dashed gray line near t = 10.
The simulated timescale in Fig. 6 extends approxi-mately to the Heisenberg time t H , near 10 3 µs. This is beyond the typical lifetime of a Rydberg state; for example, the 36p initial state studied in Liu et al. [65] has a lifetime of about 40 µs [93,94], though this could be extended using a cryostat [95]. Extending our simulated time in this way was necessary to obtain a good fit to the power-law decay region of F (t). We note, however, that the onset of the power-law decay behavior is earlier than 10 µs and therefore accessible to current experiments. Figures 7(a) and 7(b) show plots of the fidelity and the upper state saturation level as functions of natural time and approximate real time in µs. The conversions between real time in µs and natural time are within a factor of two for the two-, three-, and four-body cases shown. The two-body case with d = 32 µm is shown in black, the three-body case with d = 8 µm is shown in blue, and the four-body case with d = 6 µm is shown in red. The dashed and solid lines correspond to w = 0.05 and w = 0.45 respectively. In this example, the array spacings were chosen so the field-tuned matrix elements were of similar magnitude. The vertical dashed gray lines represent the Heisenberg times for both four-body cases and the high disorder three-body case. The remaining Heisenberg times are beyond the simulated time.
First, examining the two-body cases in Fig. 7, we observe F (t) in Fig. 7(a) to rapidly collapse to nearly zero. This rapid thermalization is also evident in the two-body s saturation levels shown in Fig. 7(b), both of which approach the maximum fraction of the expected saturation level. We find that the fidelity rapidly approaches zero for all two-body simulations. Since all of these were observed to thermalize, we do not perform the fits to F (t).
The three-body cases in Fig. 7 reveal slower dynamics. Figure 7(a) shows that while F (t) for w = 0.05 does approach zero before the corresponding t H , this decay occurs on a slower timescale than that of both two-body cases. Increasing the disorder further slows the dynamics; the three-body case with w = 0.45 yields γ = 0.3. These behaviors are reflected in the upper state saturation plotted in Fig. 7(b). While the low disorder three-body case approaches saturation prior to t H , this occurs on a much slower timescale than that of both two-body cases. The high disorder three-body case fails to reach saturation before t H .
We observe the slowest dynamics in the four-body cases of Fig. 7. The fidelity of both four-body cases fails to reach zero before their corresponding Heisenberg times. Fits to the long-time algebraic decay yield γ = 0.2 and γ = 0.05 for w = 0.05 and w = 0.45 respectively, reflecting slow thermalization. For both four-body cases, the s saturation levels remain significantly lower than the two-and three-body cases.
While the field-tuned matrix elements for each of the cases shown in Fig. 7 are of similar magnitude, the dynamics are considerably slower in the three-and fourbody results. When w = 0.05 for the three-and fourbody interactions, the behavior should be similar to systems with zero disorder. Thus, one expects scar states  The scale indicates the value of γ, with redder colors reflecting larger values of γ and bluer colors corresponding to smaller values of γ. All cases where γ ≥ 2 or F (t) rapidly collapses to zero are colored red. For the three-body case, we observe in Fig. 8(a) that only the smallest values of d with w ≤ 0.15 are predicted to thermalize. Conversely, all spacings with w ≥ 0.35 reflect slow thermalization. For larger array spacings where d ≥ 11 µm, we observe slow thermalization regardless of disorder. This is expected for large enough d as the hopping matrix elements are significantly greater than the field-tuned matrix elements in the three-body case, giving rise to scar states.
In the four-body case, the hopping matrix elements dominate even at small array spacings; we observe in Fig. 8(b) that thermalization is guaranteed only at d = 4 µm with w ≤ 0.15. Similar to the three-body case, we observe that for d ≥ 6 µm, the system is nonergodic regardless of disorder. Since the four-body field-tuned matrix elements effectively scale as r −9 , they are significantly smaller than the hopping matrix elements at almost all spacings. The combined effect of scar states and disorder leads to γ < 1 for most cases and γ < 0.1 for many.
We next consider the level spacing statistics in Fig. 9, a computational metric developed by Oganesyan and Huse in which the distribution of nearby energy levels indicates system behavior [8]. For adjacent many-body energy levels, the gap is δ n = E n+1 − E n ≥ 0 with eigenvalues E n listed in ascending order. The ratio of two consecutive gaps is defined as 0 ≤ r n = min{δ n , δ n−1 } max{δ n , δ n−1 } ≤ 1.
Oganesyan and Huse conclude that for the nonergodic phase, Poissonian statistics best describe the energy spectrum with the average value r n ≈ 0.386. A Poissonian distribution for this uncorrelated phase is expected as nearby energy levels show no level repulsion and thus are nearly randomly distributed. Conversely, they conclude that for the delocalized, thermalizing phase in which level repulsion is present, a Wigner-Dyson distribution in the Gaussian orthogonal ensemble (GOE) best describes the energy spectrum, with r n ≈ 0.5295. Oganesyan and Huse apply this metric to onedimensional lattices of spinless interacting fermions to map the localization transition. This metric has been widely used to draw numerical conclusions about MBL phase transitions. Pal and Huse estimate the location of the transition to MBL for a random-field spin-1 2 chain [9]. Burin analyzes the shape of the MBL phase transition in a one-dimensional chain of spins with long-range interactions and finds that the transitions becomes narrower as the system size is increased [96]. Bar Lev et al. study the phase diagram of a one-dimensional model and find that Wigner-Dyson statistics do not rule out sub-diffusive behavior [97].
We plot r n as a function of array spacing and disorder for the three-and four-body cases in Fig. 9(a) and Fig. 9(b) respectively. As shown by the scale, the color of the intensity plot indicates the value of r n . Redder colors are closer to Wigner-Dyson statistics consistent with thermalization and blue colors are closer to nonergodic dynamics described by Poissonian statistics.
For the three-body cases plotted in Fig. 9(a), the computed values of r n range from 0.531 to 0.450. Increasing disorder has the strongest effect on the level spacing statistics, as r n begins to deviate from a Wigner-Dyson distribution for all d with w ≥ 0.35. It is interesting to compare Fig. 9(a) with Fig. 8(a). The level statistics are very close to Wigner-Dyson for all w ≤ 0.25. However, particularly for w = 0.25, in many cases fits are consistent with γ < 1. This combination points to a delocalized phase that nonetheless thermalizes slowly. A similarly intermediate many-body critical phase was recently discovered and analyzed by Wang et al. [20,98].
The four-body cases shown in Fig. 9(b) are more consistent with the values of γ shown in Fig. 8(b). When d = 7 µm or 8 µm, we have both γ < 0.1 and r n approaching Poissonian statistics. Together, these indicate that an MBL phase is likely for the four-body interactions.
We also calculate the entanglement entropy, which has been widely used to distinguish between thermalizing and nonergodic dynamics. The EE has been found to grow only logarithmically for nonergodic cases [12,13,73], while it grows faster for thermalizing systems [99]. Figure 10 depicts a plot of the EE for the same cases shown in Fig. 7. The left side of Fig. 10 shows the early growth of the EE as a function of natural time up to t = 1.5, while the right side depicts the entanglement entropies plateauing near the end of the simulation time. The EE for each case has been scaled to its maximum value, with the approximate real time in µs depicted on the upper horizontal axis.
After the short, early period of more rapid growth, we perform numeric fits to the EE curves of all of our simulations to determine their growth rate. When γ ≥ 2, we find the EE grows faster than logarithmically, as in both two-body cases and the low disorder three-body case in Fig. 10. When γ < 1, as in the high disorder threebody and both four-body cases, we find the growth is logarithmic or slower. For the intermediate cases with 1 ≤ γ < 2, which are mostly orange-colored in Fig. 8, we find a mix of behaviors among all three metrics.
On the left portion of Fig. 10, the entanglement entropies of the three-and four-body cases initially outpace that of the two-body. This is because the spacing is much smaller and therefore the hopping matrix elements are significantly larger for those cases. While a single field-tuned interaction can entangle a small cluster of atoms, the resulting s and s ′ atoms are now able to hop and entangle more distant parts of the system. Similar behavior was observed in the ordered system of Fig. 2(c), in which the three-body case initially outpaces that of the two-body.
One might expect the long-time EE, as shown in the right half of Fig. 10, to saturate when the growth profile is The time scale is broken in two with early times on the left, and times near the end of the simulation on the right. Following a short period of rapid initial growth, numeric fits determine that both two-body cases and the three-body case with w = 0.05 grow faster-than-logarithmically. The remaining cases are consistent with slower-than-logarithmic growth. The real times for the two-and four-body cases are approximately two times larger than those of the three-body shown.
consistent with thermalization. Both the two-and threebody cases with w = 0.05 indeed approach the maximum value, while the remaining cases, including that of the thermalizing two-body case with w = 0.45, fail to do so. When disorder is added to the system by varying the positions of the atoms within the lattice, we create close clusters of atoms. The atoms in these clusters interact strongly and partially decouple from the rest of the system, inhibiting the spread of entanglement [13,43,100]. Both four-body cases in Fig. 10 show slower-thanlogarithmic growth of the EE and fail to saturate. Combining the three metrics we examined for our disordered system for the high disorder case, we find γ = 0.05, r n = 0.423, and slower-than-logarithmic growth of the EE. All three metrics point towards slow thermalization and suggest this case is a candidate for MBL.
We can use each of the three metrics to make independent conclusions about the thermodynamic behavior of our simulations. We find that in most cases, the conclusions we draw from these metrics agree with the criteria for the algebraic decay parameter described by Távora et al. We observe that for all simulations in which γ ≥ 2, the EE and r n are consistent with thermalization. We also identify intermediate cases, with 1 < γ ≤ 2. These cases show a variety of behaviors, and form the boundary between the thermalizing and the nonergodic regions in Fig. 8. We find that for the majority of simulations that fit γ < 1, the EE and the mean level spacing statistics point toward slow thermalization or failure to thermalize entirely. These slowly thermalizing systems are excellent candidates for future studies of MBL and nonergodic behavior.

IV. CONCLUSION
We have presented simulated results for onedimensional chains of Rydberg atoms exchanging energy via few-body dipole-dipole interactions. We analyzed the dynamics by studying the properties of the energy eigenstates, the decay of the initial state fidelity, the energy level spacing statistics, and the growth of entanglement entropy. While our model was motivated by ultracold Rydberg experiments, our results give broader insight into the rich dynamics of few-body interacting systems.
By studying the zero disorder case, we found numerical evidence for quantum many-body scar states when the hopping matrix elements were significantly larger than the field-tuned matrix elements. This condition is met for a wide range of experimentally typical interatomic spacings for the three-and four-body interactions; thus, we expect that scar states play an important role in slowing their dynamics.
When disorder is added to the linear array of atoms by randomly perturbing their positions, we find that thermalization is further inhibited. Power-law fits to the decay of the initial state fidelity, the mean energy level spacing ratio, and the slow growth of entanglement entropy all point to the possibility of many-body localization in our model. While we would need to include more atoms in simulation to draw a definitive conclusion about MBL, our model is accessible in current experiments.
Density-matrix renormalization group [101] and matrix product state simulation methods [102] can be effective for one-dimensional lattice systems if the EE spreads slowly. If the positions of the atoms are fixed in our model, such methods may allow for larger system sizes to be simulated in the future.
We plan to further explore the dynamics of these fewbody interactions in experiment. For lower values of the principal quantum number, the three-body resonance is well-resolved from the two-body resonance so that they can be compared in experiment. Using an initial state that does not map so closely to a scar state may result in significantly altered dynamics for the three-body case. Additionally, varying the dimensionality of the excited sample of Rydberg atoms should produce different dynamics for the two-and three-body interactions.

V. APPENDIX
In this appendix, we provide more details of our simplified three-level sps ′ model along with comparisons to a four-level spp ′ s ′ model. While the two models are not in exact numerical agreement, they share the same fundamental dynamics. In particular, both models have few-body dipole-dipole field-tuned resonant interactions that exist side-by-side with two-body hopping interactions. Furthermore, the hopping interactions have significantly larger matrix elements than the three-and fourbody field-tuned interactions in both models. This gives rise to quantum many-body scar states that greatly determine the long-time thermodynamic evolution of the system.
The most relevant comparison between the two models would be at long times when we can attempt to ascertain whether or not the system has thermalized. However, the Heisenberg times for the four-level model are too short for those long-time results to be trustworthy and we cannot therefore make such a comparison. That is, of course, the primary motivation for introducing the three-level model.

A. Two-and four-body Hamiltonians
The three-body dipole-dipole Hamiltonian operator is given in the main text by Eq. (8). We write the two-body dipole-dipole Hamiltonian operator aŝ using the same notation as in Eq. (8), where the first term is the field-tuned interaction and the second term represents the hopping interactions.
where the second term again represents the hopping interactions and the first term is the field-tuned interaction Just as for the three-body interaction in Eq. (9), V ijkℓ 4 results from summing over all possible paths from initial state to final state.
B. Comparison of the sps ′ and spp ′ s ′ models We compare the two models directly by examining the calculated probabilities as a function of time. With a few simplifying assumptions, we can find analytical solutions for the time evolution when including only three atoms. For example, consider the three-body interaction for three atoms located at the vertices of an equilateral triangle, in which we continue to ignore angular dependence.
If we continue to assume that all transition dipole moments are the same, then all field-tuned resonant Hamiltonian matrix elements will be identical. Similarly, the hopping matrix elements will also be identical to one another. These assumptions are reasonable, as all of the dipole moments for the experiments of [65] are in the neighborhood of 600 ea 0 . Beginning with the three-level sps ′ model, the four possible states are where the placement of the state labels within the ket corresponds to each atom. These states yield the Hamiltonian where v is the field-tuned resonant matrix element and u is the hopping matrix element. Diagonalizing this Hamiltonian and solving for the probability of remaining in the initial state |ppp , for the three-level model we find g sps ′ (t) = 2u 2 + 3v 2 + 3v 2 cos 2 √ u 2 + 3v 2 t 2 (u 2 + 3v 2 ) .
For the four-level spp ′ s ′ model, the seven possible states are These yield the Hamiltonian where we ignore any off-resonant hopping, like s + p → p ′ + s. Diagonalizing this Hamiltonian and solving for the probability of remaining in the initial state |ppp , we find for the four-level model. The behavior of these two models is evidently quite similar as there is little difference between g sps ′ (t) and g spp ′ s ′ (t). For the three-body energy exchange, it is typical that u ≫ v; that is, the hopping interaction is much stronger than that of the field-tuned. In this limit, the solutions for the two models approach each other asymptotically.
These analytical solutions suggest that the models should behave similarly when disorder is added. Returning to a linear array of three atoms, Fig. 11 shows the simulated p and s fractions as a function of time for a spacing of d = 5 µm. The spp ′ s ′ model is shown in dashed red and the sps ′ model in solid blue; and disorder increases in increments of 0.1 from w = 0.05 in Fig. 11(a) up to w = 0.45 in Fig. 11(e). The upper two lines initially decreasing from the maximum state fraction on each plot correspond to the p fraction, whereas the lower two lines increasing from the minimum represent the s fraction. While the numerical results for the two models are not identical, they do display the same general behavior, with similar timescales for the increase and subsequent plateau of the s fraction. Figure 12 shows the results for the spp ′ s ′ model for a linear array of nine atoms with a spacing of 8 µm. The simulated p and s fractions are plotted as a function of time, with the spp ′ s ′ and sps ′ models shown in dashed red and solid blue. Disorder increases from w = 0.05 in Fig. 12(a) up to w = 0.45 in Fig. 12 representing the p fraction and the lower two increasing from the minimum corresponding to the s fraction. Comparing these results to Fig. 11, it is clear that the state fractions of the two models become more similar as we include more atoms.
More importantly, we can also look at the longertime thermodynamic behavior of the spp ′ s ′ model and compare it to the results shown for the sps ′ model in Sec. III B. Figure 13 shows the fidelity, or initial state survival probability, and the entanglement entropy (EE) plotted as a function of time for nine atoms undergoing three-body interactions with the spp ′ s ′ model. The atoms are in a linear array with d = 8 µm and w = 0.45. In Fig. 13(a), the decay of the fidelity F (t) has the same form as that of the sps ′ model in Fig. 6 and a similar value of the fitting parameter γ. After a brief initial period of rapid growth, Fig. 13(b) shows that the EE grows logarithmically before saturating due to the finite size of the system. The long-time dynamics of the sps ′ model are significantly determined by the emergence of scar states. The scar states arise when the hopping matrix elements become much larger than the field-tuned matrix elements, as is the case for the three-and four-body interactions. Relatively large hopping matrix elements are also a property of the spp ′ s ′ model and thus also give rise to scar states, as we show here. Figure 14 shows a scatter plot of log 10 of the overlap of the energy eigenstates with the initial state as a function of energy for the spp ′ s ′ model, where the initial state is |p · · · p . In this case, we simulate nine atoms arranged in an ordered linear array with a spacing of 11 µm, just as for the sps ′ model in Fig. 5. Similar to the sps ′ model, there are a few states with high overlap near zero energy. These are the quantum many-body scar states. There is also a similar, though much more pronounced, banded structure to the overlaps.
In conclusion, the sps ′ model retains the essential features of the spps ′ model while also displaying the same dynamics.

C. Counting States
Our primary motivation for introducing the simplified model is the more favorable scaling of the Hamiltonian with the number of atoms. As the Heisenberg time grows exponentially, even small modifications that reduce the size of the Hamiltonian are beneficial.
We can count the states for the sps ′ model in the following way, using the three-body interaction as an example. Consider n total atoms and choose m of them to interact, where m is a multiple of three. Of the m interacting atoms, choose one third of them to be in the s state since the interaction takes p + p + p → s + s ′ + s ′ . Multiplying these two combinations, m n × m m/3 , yields the total number of states when m atoms interact. All that remains is to sum over all possible values of m so that where N is the total number of states. The calculations for two-and four-body interactions proceed similarly. The state counting for the spp ′ s ′ model is not too different. However, we must account for the fact that the non-interacting atoms can either be in the p or p ′ state if we allow off-resonant hopping interactions like s + p ↔ p ′ + s. In this case, the number of states for two-, three-, and four-body interactions is the same and given by We could also disregard the off-resonant hopping interactions to lower the number of states. Again using the three-body interaction as an example, we proceed just as in the sps ′ case. We choose m interacting atoms and m/3 of them to be in the s state. However, now we must choose m/3 of the remaining 2m/3 atoms to be in the p ′ state since p + p + p → s + p ′ + s ′ . This yields for the number of states. Figure 15(a) compares the number of states N for each of these three models as a function of the number of atoms included in the simulation. On our hardware, using exact diagonalization, we can simulate up to about 120 000 states. The two versions of the spp ′ s ′ model would only allow us to include nine or ten atoms. However, we can accommodate 13 atoms using the sps ′ model.
By increasing from nine or ten atoms to 13 atoms, we can dramatically lengthen the Heisenberg time t H . Figure 15(b) shows the estimated t H for a linear array with w = 0.45, d = 12 µm, and averaged over 50 simulation runs. For this example, t H ≈ 100 µs when including nine atoms. This time is insufficient to study the long-time thermodynamic behavior of this system. However, when including 13 atoms, t H is increased to a few thousand µs. We find that this is generally a sufficient amount of time; that is, the decay of the initial state fidelity has stabilized to the power law behavior predicted in [19].