Interaction-driven breakdown of dynamical localization in a kicked quantum gas

Quantum interference can terminate energy growth in a continually kicked system, via a single-particle ergodicity-breaking mechanism known as dynamical localization. The effect of many-body interactions on dynamically localized states, while important to a fundamental understanding of quantum decoherence, has remained unexplored despite a quarter-century of experimental studies. We report the experimental realization of a tunably-interacting kicked quantum rotor ensemble using a Bose-Einstein condensate in a pulsed optical lattice. We observe signatures of a prethermal localized plateau, followed for interacting samples by interaction-induced anomalous diffusion with an exponent near one half. Echo-type time reversal experiments establish the role of interactions in destroying reversibility. These results quantitatively elucidate the dynamical transition to many-body quantum chaos, advance our understanding of quantum anomalous diffusion, and delimit some possibilities for protecting quantum information in interacting driven systems.

Quantum interference can terminate energy growth in a continually kicked system, via a single-particle ergodicity-breaking mechanism known as dynamical localization. The effect of many-body interactions on dynamically localized states, while important to a fundamental understanding of quantum decoherence, has remained unexplored despite a quarter-century of experimental studies. We report the experimental realization of a tunably-interacting kicked quantum rotor ensemble using a Bose-Einstein condensate in a pulsed optical lattice. We observe signatures of a prethermal localized plateau, followed for interacting samples by interaction-induced anomalous diffusion with an exponent near one half. Echo-type time reversal experiments establish the role of interactions in destroying reversibility. These results quantitatively elucidate the dynamical transition to many-body quantum chaos, advance our understanding of quantum anomalous diffusion, and delimit some possibilities for protecting quantum information in interacting driven systems.
Ergodicity breaking in quantum matter and relaxation dynamics of thermalizing phases are two aspects of a central question of non-equilibrium many-body physics: how and when do isolated quantum systems thermalize? A growing body of theoretical and experimental work suggests that ergodicity can be avoided or hindered by a variety of mechanisms, including many-body localization [1,2], quantum many-body scarring [3,4], and prethermalization [5][6][7][8][9]. Even without ergodicity breaking, the expected emergence of quantum chaos upon the addition of interactions to driven systems is not well understood. For example, one ubiquitous but incompletelyunderstood feature of the interface between localized and ergodic regimes is anomalous diffusion [10,11], which can potentially serve as an indicator of entanglement spreading [12]. A general predictive understanding of such phenomena remains an open challenge to theory and experiment.
The quantum kicked rotor (QKR) [15][16][17] is a paradigmatic model of dynamical ergodicity breaking. While strong, repeated kicking drives a classical rotor into chaotic diffusion, the corresponding quantum rotor stops absorbing energy after a finite time, signaling the onset of dynamical localization. Despite the complete absence of disorder, this phenomenon can be understood as a manifestation of Anderson localization in momentum After setting the scattering length the trap is removed and kicking is applied with period T , pulse width τ , and amplitude V0 for n cycles. The atoms are imaged after a time-of-flight expansion [13]. (c-d) Measured axial momentum distribution versus kick number n for noninteracting (c) and interacting (d) samples, revealing collisional momentum redistribution.
space [18,19]. Although the QKR is a natural context for experimental probes of the interplay between manybody interactions and dynamical localization, the quantum thermodynamics of interacting kicked rotors have not previously been experimentally explored. Depending on how interactions are introduced into the model, theoretical studies have predicted a variety of novel dynamical phenomena ranging from anomalous diffusion [20][21][22][23] to classical prethermalization [24] to many-body dynamical localization [25,26].
Here we report the first experimental study of dynamical localization in the presence of tunable contact interactions, which are nonlocal in momentum space. These experiments investigate a 7 Li Bose-Einstein condensate (BEC) kicked n times at period T by a far-detuned optical lattice of spacing d = 532 nm and depth V 0 for duration τ (see Fig. 1). We report momentum and energy in units of k L = π/d and E R = 2 k 2 L /2m with m the Here V0 = 64ER, T = 1.2 µs and τ = 300 ns (K ≈ 2.3 andk ≈ 1.5). The inset contrasts interaction-induced delocalization and anomalous diffusion with classical diffusion caused by random kicking, the latter achieved by adding random offsets to the average kick spacing T drawn uniformly from the interval [−T /4, T /4]. The solid curve is noninteracting quantum theory and the dotted line is a diffusion curve 4Dn/k 2 with D ≈ 0.19 extracted from the classical standard map [14]. The red dot-dashed line is a subdiffusive √ n law serving as a guide to the eye. (b) Momentum-space IPR with transverse dimensions integrated out. The shaded regions are predictions for two exponentially localized distributions with 1/e localization length k loc = √ E loc ≈ 1.6(1)kL [13]. (c-e) Normalized smoothed momentum space densities at various n. (f-h) The same densities on a logarithmic scale. The orange dotted and purple dashed lines are exponentially localized curves exp(−k/k loc ) with amplitudes normalized to match the peak of the measured distributions at the given n. (i) Deviation from exponential localization over time based on integrated ratio between measured and exponential distributions with error bars computed from uncertainty in k loc [13]. mass of 7 Li. The single-particle QKR is defined by the 1cycle Floquet map U = e −ikk 2 /2 e −iK cos z/k describing a sharp cosine potential impulse followed by free evolution. Here k and z are momentum and position respectively, K =kV 0 τ /2 is the stochasticity parameter characterizing kicking strength andk = 8E R T / is an effective Planck's constant determined by the kick period. Absorption imaging after free expansion is used to measure the momentum distribution; see the supplementary text for a full analysis of systematic effects in this procedure, such as dynamics transverse to the lattice direction. Interatomic interactions are varied by tuning the s-wave scattering length a (reported in units of the Bohr radius a 0 ) using a magnetic Feshbach resonance. While the kicking primarily couples discrete momentum states along a single dimension, the atoms are entirely unconfined between kicks; scattering between momentum modes thus couples the system to a bath of transverse free-particle states.
The main result of this work is presented in Fig. 2a.
While a noninteracting sample exhibits dynamical localization, saturating to a finite energy for over 800 kicks, interacting samples clearly demonstrate the destruction of the dynamically localized plateau with increasing scattering length. At intermediate interaction strength (a = 240a 0 ), we observe saturation to the same energy as non-interacting samples for approximately 300 kicks, suggesting the existence of a reasonably long-lived prethermal state. In contrast, the 760a 0 trace exceeds this localized energy after around 100 kicks; whether a quasiequilibrium dynamical state is truly established in this stronger-interacting sample is less clear. Fig. 2b shows another aspect of the same evolution, plotting the momentum space inverse participation ratio (IPR) versus kick number. The IPR characterizes the number of states over which the system is distributed, thereby also probing how collisional momentum redistribution washes out the originally discrete momentum modes, a process less easily inferred from energy measurements. While the 240a 0 data exhibit a clear steady-state behavior for 100 kicks, the 760a 0 IPR decreases monotonically for almost the entire experiment. A second key result of these measurements is that the observed delocalizing dynamics clearly exhibit anomalous diffusion: it appears that even interacting quantum kicked rotors absorb energy much more slowly than classical rotors. The inset of Fig. 2a compares the nature of the observed interaction-induced subdiffusive delocalization with linear energy growth in the classically chaotic model. We experimentally simulate classical dynamics by adding stochastic fluctuations to the kicking period T , making use of the known sensitivity of dynamical localization to timing noise [27]. These experimental data agree both with single-particle quantum numerics and with the linear energy growth predicted by the classical standard map [14], and stand in clear contrast to the measured interaction-induced anomalous diffusion away from the dynamically localized state. The dot-dashed red line indicates a √ n energy growth, and fitting the latetime data to n α yields anomalous diffusion exponents α in the range [0.4, 0.6]. For reference, 1D Gross-Pitaevskii simulations on a ring [22] predict α ∈ [0.5, 0.8], though a direct quantitative comparison to theory is challenging due to the significant condensate depletion and the three-dimensional nature of the experiment. Theoretical studies of the effect of local nonlinearity on real-space Anderson localization instead suggest α ∈ [0.3, 0.4] [28,29], but the long-range nature of contact interactions in momentum space similarly complicates comparison. This clear observation of anomalous diffusion in the interacting quantum kicked rotor raises a variety of fascinating questions for future exploration. What correlations are responsible for the anomalous diffusion dynamics? What feature of the interacting system prevents the interacting QKR from heating classically? What theoretical framework is appropriate for quantitatively predicting wavefunction evolution in this regime? What are the universality properties of the subdiffusive exponent?
For further insight into the dynamics of kicked interacting quantum systems we examine the evolution of the momentum distribution, shown in Figs. 2c-e. We observe a clear distinction between the noninteracting samples, which settle at a sharply-peaked dynamically-localized momentum distribution, and the interacting samples, which gradually smear out in momentum space due to scattering. Plotting these same densities on a logarithmic scale in Fig. 2f-h illuminates the destruction of dynamical Anderson localization by assessing the departure from exponentially-localized Floquet states. The smeared-out lower-energy modes actually appear to maintain the expected localization length, and thus do not trivially indicate a departure from exponential localization. This observation is also reflected in the fact that two predictions based on exponentially localized distributions bound the measured IPR in Fig. 2b. Instead, the departure from exponential localization manifests in the emergence of increased relative population in the tails of the distribution. It is interesting to note that recent theory suggests that even many-body dynamically localized phases are expected to exhibit universal power-law decaying tails [30]. In Fig. 2i we quantitatively characterize the overall deviation from exponential localization [13], revealing a break time near 200 kicks for both interaction strengths. These findings provide a second experimental signature of the destruction of the dynamically localized state by interactions, now both at the level of macroscopic observables and squared wavefunctions.
The onset of energy delocalization due to interactions indicates a transition to the regime of many-body quantum chaos, which can be probed directly by studying time-reversal dynamics [31,32]. In Fig. 3 we probe the onset of chaotic dynamics by measuring the effect of interactions on a Loschmidt echo time-reversal protocol [33,34]. The echo is realized using quantum resonances [35] which occur fork = 2πq with q rational; in particular for q = 2 (T ≈ 9.95 µs), the free evolution in U largely vanishes and effective time reversal can be achieved by setting q to 1 for a single kick halfway through the sequence [13,33,34]. This procedure would create exact time reversal for a single zeroquasimomentum state in the absence of interactions. Due to finite quasimomentum spread and non-reversed interactions, the reversal is imperfect, yielding a Loschmidt where U 1 and U 2 are timeevolution operators differing by some perturbation. Perhaps surprisingly, F initially increases as the scattering length a is turned up from zero. In this regime U 1 and U 2 are primarily distinguished by the failure to reverse kinetic energy, and thus the increase can be explained by Thomas-Fermi reduction of the initial state momentum spread. Eventually, for large enough a, the interaction becomes the primary perturbation and F begins to decrease with a, marking the transition to predominantly interaction-induced irreversibility. The decay of fidelity with total number of kicks in a Loschmidt echo experiment is shown in Fig. 3b. The use of Loschmidt echo techniques as a probe of many-body quantum chaos not only illuminates the origins of the delocalizing dynamics we observe, but opens up the possibility of extending these protocols to probe scrambling in many-body quantum chaotic systems [36].
In conclusion, we have experimentally realized a manybody ensemble of quantum kicked rotors. Following the evolution of interacting samples over hundreds of kicks, we observe signatures of an initial prethermal state, followed by an interaction-induced breakdown of dynamical localization via anomalous diffusion, signaling the onset of many-body quantum chaos. Characterization of the departure from the dynamically localized state indicates subdiffusive energy growth with an exponent near 0.5, easily distinguishable from classical Joule heating in a randomly kicked system, and reveals momentum space distributions which are not exponentially localized. Measuring Loschmidt echo time-reversal dynamics with a quantum resonance enabled us to directly probe the role of interaction-induced irreversibility in driving a transition to many-body quantum chaos. Together, these results demonstrate and quantitatively illuminate the emergence of interaction-driven quantum chaos in a paradigmatic localized system.
We acknowledge helpful conversations with Adam Rançon, Norman Yao, and Thomas Schuster. During the performance of the research described here we became aware of related efforts underway in another experimental group [37] which have reached similar conclusions using a complementary approach.

Methods
Experimental platform and sequence. The experiments begin with a Bose-Einstein condensate (BEC) of around 10 5 7 Li atoms in a far-detuned optical dipole trap with trapping frequencies ω x,z /2π ≈ 40 Hz and ω y /2π ≈ 56 Hz, where z is the axis of the optical lattice, y is the direction of gravity, and x is the remaining orthogonal axis. The condensate is produced by optical evaporation at an s-wave scattering length of a = 240a 0 , set by an applied magnetic field in the vicinity of the broad Feshbach resonance at 737 Gauss [38]. Immediately after evaporation, the fields are ramped to their desired value over 60-90 ms and maintained for the remainder of the experiment. The dipole trap is then extinguished and the BEC repeatedly subjected to a pulsed 1D optical lattice with lattice constant d = 532 nm, laser wave vector k L = π/d, and recoil energy E R = 2 k 2 L /2m with m the mass of 7 Li. The full dynamics are then well described by the second-quantized Hamiltonian The key kick parameters are the lattice depth V 0 , effective pulse width τ , and kick period T . V 0 is calibrated through a standard Kapitza-Dirac diffraction technique. f τ (t) denotes a unit amplitude pulse function beginning at t = 0 of width τ . The experimental pulse is approximated by a piecewise function with a linear rise and fall of 200 ns duration before and after a plateau of variable hold time. For the experimental data in the main text with τ = 300 ns between the half-maximum points, this hold duration is 100 ns. The scattering length a determines the two-body coupling coefficient g = 4π 2 a/m. Here I(x, y) denotes the transverse intensity profile of the lattice beams normalized to unity maximum; this is approximately Gaussian I(x, y) ≈ e −2(x 2 +y 2 )/σ 2 with a measured 1/e 2 beam radius of σ ≈ 65 µm. The total duration of kicking is at most 1 ms for our longest experiments, significantly shorter than the 4 ms it takes the BEC to fall under the influence of gravity through the lattice beam waist.
To measure the momentum distribution, we perform absorption imaging of the atoms after free expansion. The time-of-flight (TOF) duration is 3.5 ms for the delocalization data and 2 ms for the Loschmidt data. Due to the low mass of 7 Li and the breadth of the Feshbach resonance, coil inductance prevents sweeping the magnetic fields to the noninteracting regime for this expansion period. This means additional scattering occurs during expansion, which may lead to systematic errors in the measured quantities (see supplementary section 6). For the energy, we are able to account for this scattering in our analysis due to the energy-conserving nature of the collisions. For metrics such as the IPR, this systematic is challenging to avoid. However by tracking the evolution of these observables as a function of kick number at a fixed TOF duration, we can largely attribute the qualitative observed dynamics to the evolution under the Hamiltonian (1) as opposed to the expansion. At large n, the majority of scattering happens during the kicking duration so expansion effects become negligible.
Delocalization data analysis. This section discusses the analysis behind Fig. 2. Because the momentum distributions of the interacting samples change significantly over the course of the delocalization experiments, the quantities shown in Fig. 2 are computed directly from raw or averaged images as opposed to fitting procedures. However, this can make measurable quantities such as energy sensitive to noise, especially near the edge of the camera sensor due to the quadratic weighting. To maximize the signal-to-noise ratio in our measurement, we analyze raw images using an adaptive region-of-interest (ROI). First, a single base ROI capturing all detectable atoms at all times is created for each interaction strength. The integrated density in this ROI is used to post-select images with total atom numbers falling within a ±10% window of the mean, in order to reduce variations in the interaction energy, which depends directly on atom density. For these data we take 10 images at each kick number, of which typically 4-7 are discarded by this postselection procedure. The ROI boundaries at each kick number are then determined by the points at which the cumulative summed distributions of the averaged image outward from a center point reach a threshold value. The thresholds are set empirically and the boundaries obtained by the following procedure. First we compute the transverse bound by integrating out the entire axial direction to get the overall transverse distribution, find the point it crosses an 85% threshold and then expand the resulting boundary by a factor of 1.5 (1.2 in the delocalization data of supplementary section 1) to ensure all atoms are captured. We then compute an axial boundary going point by point along the transverse direction; at each transverse point we integrate over 10 neighboring transverse pixel rows to get a "local" axial distribution, find the point it crosses a 99.8% threshold and expand by a factor 1.15. Finally we smooth each ROI boundary and take a moving average across different kick numbers (4 on each side). Crucially, we have confirmed that the qualitative observation of delocalization is not significantly altered from the simple case where we use just the initial single base ROI across all shots. However, the details of the trends should be more accurately captured by the adaptive procedure because the signal-to-noise ratio over the ROI is optimized at each kick number. All measurable quantities are then calculated from the imaged densities within this region.
Since we do not observe any substantial atom loss during the kicking duration, we treat the imaged atomic densities as normalized distributions. For Figs. 2a-b, we compute the measured quantities from individual experimental runs and then average the results, with the reported error bar as the standard error of the mean. For Fig. 2i, we instead compute the averaged distributions first before computing the deviation from exponential localization; the errorbars are computed from a Monte-Carlo simulation of the uncertainty in k loc discussed later in this section. A smoothing filter is applied to the displayed densities in Figs. 2c-h for visual clarity, but not in the subsequent calculation of the localization deviation in Fig. 2i.
To measure the energy, we compute the post-expansion spatial variance of the distribution in both the kicking z and transverse x directions of the image. Assuming cylindrical symmetry, the kinetic energy is then calculated as m z 2 + 2 x 2 /2t 2 TOF with t TOF ≈ 3.5 ms (see section 2.4 for a discussion of possible corrections to the conversion of position to energy). For an accurate measurement of the interacting samples, inclusion of the transverse energy is necessary to account for energy-conserving scattering processes that occur both during the kicking and TOF. In addition, the inhomogeneous intensity profile of the beam I(x, y) leads to a transverse energy oscillation in all samples including the noninteracting ones (see supplementary section 5). Since we are not interested in this effect, we remove it to leading order by subtracting off the noninteracting transverse energy from each trace, so that the noninteracting energy is purely the kinetic energy along the kicking direction. To compute the error bars on the interacting data, we add the error of the total interacting energy and noninteracting transverse energy in quadrature. The single-particle localization energy E loc is estimated by averaging the noninteracting trace for n ≥ 100, and the reported uncertainty is based on the standard deviation of those points. We note that this uncertainty is not only due to experimental imperfections, but also due to natural dynamical fluctuations, as evidenced by the results of noninteracting simulations like those shown in Fig. S3.
We compute an effective 1D momentum-space IPR by first integrating out the transverse dimension and then summing the squares of the subsequent normalized axial density. We confirmed that this qualitatively matches the result of directly integrating the squared 2D distribution while largely eliminating the beam-induced transverse oscillation. Specifically for computing this metric, we apply a smoothing filter to the normalized densities consistently across all 3 interaction strengths. This suppresses high-frequency background noise which sets a lower bound on the measurable IPR due to the squaring procedure. The measured values are compared to two predictions based on an exponentially localized distribution. The blue shaded region is obtained by numerically computing the IPR for the momentum space distribution exp(−|k|/k loc ) j exp −(k − 2k L j) 2 /w 2 , which models a Gaussian comb with an exponential envelope. This is a reasonable expectation for a finite-size, localized noninteracting condensate occupying only discrete momentum modes. The width parameter w is measured from fitting the n = 0 noninteracting condensate and takes into account the momentum-space resolution of the TOF given the finite condensate spatial extent. The width of the region is based on Monte Carlo simulation of uncertainty in k loc , where the resulting distribution is fit to a Gaussian to extract the mean and standard deviation. The green shaded region is calculated analytically for a pure exponential distribution of infinite extent and is given by 1/4k loc . Taking into account the finite width of the imaging region changes the distribution normalization and leads to the following correction factor (1 − exp(−2k 0 /k loc )/(1 − exp(−k 0 /k loc )) 2 ; here k 0 ≈ 9.85k L is the half-width of our images which yields a negligible correction factor of ≈ 1.006. The width of the region is computed through linearized error propagation.
In Fig. 2i, the plotted localization metric is Here, r(k) = |ψ(k)| 2 / exp(−k/k loc ) is the ratio of the measured axial density denoted |ψ(k)| 2 and an exponential localization envelope. Here the maximum of |ψ(k)| 2 is set to unity. Taking the maximum of r(k) − 1 and 0 ensures that the result is only sensitive to regions of the distribution which decay more slowly than exponentially. That is, it interprets 0 as "at least exponentially" localized with respect to a given localization length, and thus characterizes departures from a given dynamically localized state in the traditional sense of exponentially localized wavefunctions. We note however that the system remaining exponentially localized but with a larger localization length would result in a non-zero value for this metric, which motivates the direct inspection of the distributions in Figs. 2f-h. The reported values and errorbars are extracted by propagating a Gaussian uncertainty in the measured k loc through a Monte-Carlo simulation. We find that the resulting distributions interpolate between sharply peaked at 0 with a rapid fall-off when well-localized, to positively skewed with non-zero peak in the delocalized regime. We empirically find that a log-normal distribution fits the Monte Carlo result well, and we use this fit to extract the data reported in Fig. 2i. In particular, the markers indicate the mean of the distribution and the errorbars represent the interquartile range containing the central 50% of the distribution. Because of the skewness, we investigate the Monte Carlo simulated distributions in more detail in supplementary section 2.
Loschmidt experimental sequence and data analysis. Here we discuss the methods and analysis used to produce Fig. 3. The Loschmidt experiments begin similarly to the previously described sequence; for an N kick Loschmidt sequence, the BEC is first kicked N/2 times near quantum resonance at the parameters V 0 ≈ 50E R , τ = 300 ns, and T = 9.93 µs. For this data, we adjusted the lattice depth V 0 for different interaction strengths to achieve the same amount of absorbed energy after the first N/2 kicks. This compensates for a decrease in energy absorption at the same lattice depth for higher interaction strengths, which we attribute to the increase of the Thomas-Fermi radius of the BEC relative to the lattice beam size. Neglecting this effect would artificially enhance the fidelity at very large interaction strengths due to a reduction in the effective stochasticity parameter K. We plot the zero mode fraction after the first N/2 kicks without time-reversal (denoted F Gauss) in Fig. S6 to benchmark this kicking amplitude normalization procedure.
After the first N/2 kicks, we wait a half period T /2 to shift the wavefunction spatially by half a lattice spacing, causing the sign of subsequent kicks to be reversed. We then apply another sequence of N/2 kicks using the same lattice parameters to complete the echo sequence. The time series in Fig. 3c-e show absorption images averaged over 5 shots for each kick number n in a N = 10 experiment. Since we begin with a zero-momentum condensate mode, to measure the Loschmidt fidelity we simply need to count the fraction of atoms remaining in this mode. While atoms in other momentum modes coupled by the lattice are easily distinguished, atoms that have undergone scattering events into a smeared-out background distribution are not always well-separated. Thus, to extract the return fraction we fit the axial atomic distribution around the zero-momentum mode with a pair of Gaussians of varying width. The narrower Gaussian accounts for atoms remaining in the zero-momentum condensate after expansion, while the broader Gaussian measures the atoms that have been collisionally ejected from the condensate [39,40]. In Fig. 3b, we show the fraction of atoms remaining in the narrow Gaussian and use this quantity as an estimate of the Loschmidt fidelity. Scattering during the expansion means that this necessarily underestimates the true fidelity, a possibility further addressed in supplementary section 6.
Noninteracting QKR numerics. One-dimensional simulations of the noninteracting kicked rotor problem for comparison with experimental data are executed in two ways. We either perform a split-step Fourier integration of the QKR Hamiltonian (2) (ignoring the transverse distribution I(x, y)) to model the finite-width pulse shapes, or iterate the QKR Floquet map described in the main text. The simulations are typically performed with periodic boundary conditions over a single lattice site (except when modeling the TOF readout; see supplementary section 4). We perform a Gaussian sampling of quasimomenta with standard deviation ∼ 0.1k L , in rough accordance with the measured BEC temperature of around 10-15 nK. For simulation of the stochastic kicking protocol, we use the same techniques and additionally average over 100 different realizations of the fluctuations (note this is slightly different than in the experiment where a single kick period disorder realization is used). To supplement the dynamical delocalization signals shown in Fig. 2 and demonstrate that this is not a particularly fine-tuned phenomenon in the kicking parameter space, in Fig. S1 we show the same metrics for a larger kicking period T = 2.2 µs. The overall picture is unchanged, as the interacting samples show starkly different behavior from the noninteracting traces, departing from the localized value of each metric after a variable break time. Here the energy delocalization is obscured slightly as the different interaction strengths seem to initially localize to different energies. We attribute this partly to Thomas-Fermi expansion which reduces both the effective lattice depth experienced by the condensate and the initial kinetic energy of the sample, though we do not entirely rule out the possibility of different early-time prethermal behavior across interaction strengths. The correlation between localization length and quasimomentum spread is observed in noninteracting numerics. The different-colored shaded regions indicate our best estimates for the different localization energies at the 3 interaction strengths by computing the mean energy (and standard deviation) over windows of n where the data are minimally changing. These values are used to compute the exponential localization deviation in Fig. S1c. We do note a small trend visible at the end of the noninteracting traces; numerics suggest that this is consistent with variations in the localization length that occur over time for certain kicking parameters.

MONTE CARLO DISTRIBUTIONS FOR EXPONENTIAL LOCALIZATION DEVIATION
In Fig. S2, we show further details on the Monte Carlo simulated distributions for quantitatively characterizing the deviation from exponential localization in Figs. 2i and S1c (in particular this data corresponds to Fig. 2i). The distributions are generated by computing the defined deviation parameter for 10 4 values of k loc drawn from a Gaussian centered at 1.58k L and with standard deviation 0.12k L . An example distribution for a sample which has delocalized is shown in Fig. S2a, clearly showing the skewed probability densities we obtain from this procedure. The solid orange line indicates the log-normal distribution fit we use to extract parameters such as the mean and interquartile range of the distribution. We note that the use of a log-normal distribution here is only motivated empirically as a systematic method to determine such quantities.
In Fig. S2b, we contrast how these simulated distributions evolve in time for localized noninteracting samples and delocalizing interacting ones. In the noninteracting case, the distributions are extremely sharply peaked at 0 and are relatively unchanging in-time, agreeing with the expectation of dynamical localization. In the latter, however, the distribution is only peaked at 0 for short times, a behavior indicative of the finite duration prethermal plateau we report, and gradually shifts away to non-zero values as the sample heats up. Importantly, at the later times the 760a 0 distribution has essentially vanishing probability density at 0 deviation, allowing us to confidently claim observation of departure from exponential localization. In Fig. S2c, we confirm that the reported behavior of deviation over time in Fig. 2I would not qualitatively change if we instead used the median or mode of the distribution instead of the mean.

SYSTEMATICS: FINITE PULSE WIDTH
The delta-function kicking assumption in the theoretical QKR model is not perfectly realized in experiment owing to the finite atomic mass of 7 Li. The assumption corresponds to the Raman-Nath diffraction regime which is approximately expressed by the condition 2 √ V 0 E R τ / 1 [41]. For the experiment with τ = 0.3 µs and V 0 = 64E R as in Fig. 2, this parameter is approximately 0.76 (there is an ambiguity of a factor 2π in defining the condition [42], which would reduce the parameter to 0.12). Either way, this suggests that the system is in between the Raman-Nath and Bragg diffraction regimes and thus finite-pulse-width effects require careful investigation.
We numerically explore the effects of realistic pulse duration on single-particle QKR localization by comparing square pulse simulations of varying pulse width τ to the delta-kick Floquet map solution. To make this comparison, we keep the effective stochasticity parameter K ∼ V 0 τ characterizing the kicking strength constant as we let τ → 0. The τ = 300 ns trace is comparable to the data in the main text. (b) Equivalent simulation for K = 5. The relative difference in localization energy between the achievable finite pulse durations in our experiment and the delta-kick limit becomes much more substantial at larger K.
Results for two values of K are shown in Fig. S3. In general, we find that larger pulse duration tends to decrease the localization energy, which from a classical perspective corresponds to a particle traversing a significant part of the cosine potential during the kick and thus feeling a smaller effective impulse. This effect depends on the value of K which determines the extent to which higher momentum modes are excited in the localized state. The many-body delocalization data in the main text were taken around K ≈ 2.3. These simulations indicate that for this data there is a roughly 20% decrease in the measured noninteracting localization energy with respect to the delta-kick limit.
The finite pulse duration leads to an effective kicking strength which decays with increasing momentum, causing even the classical phase space to localize above a certain momentum [43]. This is an important consideration in probing the destruction of the quantum dynamical localization which occurs in the classically chaotic regime. For our parameters, the estimate given in [43] for the momentum boundary between classically chaotic and integrable regions due to pulse width is roughly ±33.2k L , which is much larger than any excitation we observe in the data for Fig. 2 and S1. Thus we do not expect that the finite pulse duration qualitatively affects the observed delocalization dynamics.
We note that the reduction of absorbed energy by finite pulse width does play a practical role in determining which sets of system/kicking parameters are amenable to observation of interaction-induced delocalization. Because collisional processes are proportional to real-space density-density overlaps between different momentum modes, for a poor choice of kicking parameters a strong excitation of higher momentum modes in conjunction with real-space expansion (discussed in the next section) may rapidly dilute the system, and yield an effectively non-interacting sample before the interaction-induced delocalization break time. 2 µs and τ = 300 ns), but with lattice depth reduced to V0 = 38ER to compensate for modeling only delta function kicking. We expect that the strong oscillations shown here are largely washed out in the experimental data by effective averaging over kicking strengths due to transverse extent of the condensate and fluctuations in the beam intensity from run-to-run. The traces here are sampled every 30 kicks for visual clarity, and the energy actually fluctuates at a much higher frequency. This simulation reveals that for the main data, the correct energy is straightforwardly extracted by simply using the TOF time to convert position to momentum. The vertical line indicates the maximum kick number reached in the data. (b) A similar simulation but with T = 2.2 µs roughly corresponding to the supplementary delocalization data in Fig. S1 (with adjusted lattice depth V0 = 52ER). To illustrate the potential pitfalls of the readout method, here we model a TOF duration of only 2 ms, as opposed to the 3.5 ms used in the experiment. In this case, using only the TOF time in the velocity conversion leads to a false delocalization signal at late times. Instead, we show than an additional method using the TOF duration plus half the kicking time leads to the most faithful representation of the QKR energy for the longest time. We chose to use this conversion in the analysis of S1, though our simulations do not indicate a large difference between these two methods when modeling the full 3.5 ms TOF and restricting only to the max kick number indicated by the vertical line. Here the sampling is every 15 kicks for visual clarity.

SYSTEMATICS: POSITION SPACE DYNAMICS AND TOF CONVERSION
From a theoretical perspective, the difference between open and periodic boundary conditions in the single-particle QKR is resolved by Bloch's theorem [44,45]. Different quasimomenta evolve independently, manifesting different realizations of pseudo-randomness in the Anderson model mapping. The connection between the theoretical QKR and kicked quantum gas experiments is made by considering ensemble averages over quasimomenta. However, in practice, experimental readout of the kinetic energy even in the noninteracting case can be further complicated by spatial motion in a non-compact position variable during the course of the kicking. This effect was largely negligible for many previous QKR realizations using heavy atomic species and/or short kicking durations, but in these experiments using light 7 Li atoms and large kick numbers, careful consideration of the effect is required for an accurate energy measurement. Ideally, one would simply extend the TOF duration to suppress such effects, but technical limitations associated with the imaging procedure mean that this cannot be done indefinitely.
In interpreting the TOF absorption images as momentum space distributions, one must convert pixel position to velocity by dividing by an appropriate time. Without spatial motion, the correct time is trivially just the TOF duration. With spatial motion, a strict lower bound on the velocity conversion is set by the combined duration of the kicking and TOF which is equivalent to the assumption that each momentum mode propagates ballistically for the entire course of the experiment, ignoring the reshuffling of momentum modes by repeated kicking. To determine which conversion scheme leads to the most accurate energy measurement, we simulate the delta-function QKR model with an extended position space variable to model the TOF expansion explicitly. We are then able to compare the exact energy with various position-to-velocity conversions to determine the best metric.
Simulations for different kicking parameters and TOF durations are shown in Fig. S4a  conversion can change. In Fig. S4a corresponding to the main data, we show that simply using the TOF as a conversion factor matches the true energy. For the simulations in S4b, however, adding in half of the kicking duration gives a substantially more accurate approximation of the true energy than the simple TOF conversion, which produces a false delocalization signal at longer times. We also examine "box"-counting schemes where the image is instead binned into discrete modes which are multiples of 2k L momentum, though the added complexity of this scheme is not justified by the results.

SYSTEMATICS: BEAM-INDUCED TRANSVERSE DYNAMICS
While the main concern of QKR experiments is with momenta along the lattice direction, our experiments are three-dimensional and degrees of freedom transverse to the lattice beam cannot in general be ignored especially in the presence of scattering. In Fig. S5, we explicitly show the measured transverse kinetic energy for the main delocalization data in the text. Here we can see all three interaction strengths undergoing an oscillation in their transverse energy, which can be interpreted simply as harmonic motion in the time-averaged intensity distribution of the pulsed Gaussian lattice beam. The clear difference in the evolution between the different interaction strengths indicates the effects of 3D scattering for a system with uniform I(x, y). As discussed in the Methods section 1.2, this motivates inclusion of the difference between the noninteracting and interacting transverse energy traces in the plotted energy of Fig. 2. We have separately confirmed that ignoring the transverse dynamics altogether does not eliminate the observed delocalization signal.

SYSTEMATICS: EFFECTS OF SCATTERING ON MEASURED LOSCHMIDT FIDELITY
Readout of the momentum distribution of interacting samples can be complicated by scattering during the TOF, and this particularly impacts the measurements of metrics such as IPR and return probability F , where the signatures of scattering events occurring at different stages of the experiment (i.e. during the lattice pulse trains versus during the time-of-flight) are not easily extracted from the resulting distribution. Here we examine how scattering affects the reported Loschmidt echo return probability shown in Fig. 3. In Fig. S6 we present a comparison of two different methods for measuring the fidelity, which we argue should bound the true value and indicate the effect of this systematic. The Gaussian fitting method was described in the methods section 1.3 and presented in Fig. 3. The raw counting method computes the fidelity by integrating the raw distribution in a ±k L width around the central mode. If all the scattering occurs prior to the TOF, then the Gaussian fitting method is the appropriate counting procedure, as it discards all scattered atoms and only counts the remaining zero-order atoms. If, however, the majority of scattering  occurs during the TOF, then this population should be included in the return probability, and so the raw counting method would more accurately reflect the true fidelity.
In Fig. S6A, we compare these two methods before (red) and after (blue) the application of the time-reversal kicks, as a function of scattering length. In both cases we find that the raw counting method measures a higher fidelity than the Gaussian fitting method due to accounting for the scattered population. As expected, the two converge in the noninteracting limit where the overall scattered population vanishes but diverge as the scattering length and consequently the scattered fraction increase (the dependence of scattered fraction on scattering length is plotted in Fig. S7). This behavior of the Gaussian fitting and raw counting methods is consistent with the limits of validity expected for each, and supports the claim that the two methods bound the systematic measurement error in counting the zero-order population that results from scattering during the TOF.
Having established approximate bounds for the true fidelity as a function of scattering length, we further remark that both methods produce a time-reversed fidelity which exhibits a crossover in behavior as a function of scattering length. The fidelity as a function of a grows with weak interactions, and it saturates, or even decreases, as a function of a with stronger interactions. We note also that imperfect calibration of the effective kicking strength could give rise to errors in the measured return fidelity. In Fig. S6B, we attempt to account for this effect across the two methods by considering the difference in fidelity before and after the set of time-reversal kicks. This produces two curves of similar functional form, further supporting our observation of a crossover into interaction-induced irreversibility.

MAPPING THE QKR TO A KICKED SPIN CHAIN MODEL
The interplay between dynamical localization and quantum chaos is a topic of broad current interest, relevant in contexts well beyond the experimental model of the quantum kicked rotor which we explore here. To highlight this breadth, we describe and quantitatively explore a mapping from the QKR to a kicked Heisenberg spin chain, which can be used as a basis for generalization and further exploration. Specifically, we investigate the Floquet map U kickedXXZ = e −ikH quad /4 e −iKHXXZ/4k , with H quad = j (j + β) 2 σ z j a quadratic field and H XXZ = j (σ x j σ x j+1 + σ y j σ y j+1 + ∆σ z j σ z j+1 ) an XXZ Hamiltonian with σ x,y,z the Pauli matrices and ∆ the anisotropy parameter. β represents a field-center offset which serves as an analogue to quasimomentum in the standard QKR problem. Total magnetization M = j σ z j is a conserved quantity. This model is known to correspond to the QKR in single and few-body regimes [47,48]; for a single particle or spin flip the correspondence is essentially exact up to finite size effects. We perform exact diagonalization of U kickedXXZ on systems of length up to L = 14. We center the chains about j = 0 (i.e. j ∈ [1 − L, L − 1]/2, integer for L odd and half-integer for L even) and consider open boundary conditions. To first illustrate the mapping, single-particle dynamical localization and quantum resonance in the spin model are demonstrated in Fig. S8a-c by considering the single spin-flip sector M = −L + 2. In the zoom-in of Fig. S8c, we verify the QKR correspondence experimentally by comparing the kicked XXX (∆ = 1) numerics to our experimental observation of fractional q quantum resonances in the QKR, and find excellent agreement.
Given that the Heisenberg chain in a random magnetic field is a prototype model for traditional many-body localization (MBL), U kickedXXZ appears well suited to address the question of whether the emergent Floquet pseudo- randomness in a disorder-free kicked system is sufficient to reproduce MBL phenomenology. To numerically probe many-body dynamical localization in the kicked spin model, we study the evolution of an initial Néel state ordering |↑↓↑↓ . . . in the M = 0 sector with multiple spin-flips, and observe long-lived persistence of ordering at finite K values (Fig. S8f). As shown in Fig. S8g, this persistence lasts in the infinite time limit. These conclusions are obtained by computing the the time average of the staggered magnetization O = j (−1) j σ z j /L after n cycles, given by O n = n i=1 ψ| U †i kickedXXZ OU i kickedXXZ |ψ /n. The infinite-time limit is then calculated via the Floquet diagonal ensemble as lim n→∞ O n = α |c α | 2 ψ α | O |ψ α . Here c α = ψ α |ψ are the coefficients of the initial state |ψ in the basis of the many-body Floquet states |ψ α [49].
A transition from many-body dynamical localization (MBDL) to ergodicity at larger interaction strengths is indicated both by the Néel state persistence (Fig. S8f-g) and by the Floquet level-spacing gap ratio parameter r [13,46] (Fig. S8h). The gap ratio is defined as r α = min (δ α , δ α+1 ) /max (δ α , δ α+1 ), where δ α = α+1 − α is the gap between consecutive quasi-energies α , which we order on the interval [−π, π]. To compute r , we average r α over the M = 0 sector as well as for 100 values of β drawn from a normal distribution of standard deviation 0.1. We interpret the transition in r from the Poisson prediction to that of the circular-orthogonal ensemble as indication of distinct parameter regimes of MBDL and chaos in the many-body Floquet system. Similar results are found away from the isotropic ∆ = 1 point, as shown in Fig. S9. In the limit ∆ = 0, we checked that the statistics are Poisson for all values of K shown here. While extending to larger system sizes to extrapolate toward the thermodynamic limit remains an important task, these numerical results signal a true many-body dynamically localized state in kicked spin-chains with no disorder. This indicates a promising path towards experimentally observing MBDL in current-day quantum simulator platforms [50][51][52] and highlights connections between paradigmatic kicked spin models and the quantum kicked rotor. For the single-particle QKR, it is standard to consider periodic pulsing of the spatial potential separated by intervals of free kinetic energy evolution. In the spin-chain mapping, this corresponds to pulsed spin-exchange interactions and free evolution in a quadratic magnetic field. While here we used the exact mapping of the QKR parameters K andk in U to the spin model U kickedXXZ , physically we took the interpretation of an interacting XXZ Hamiltonian with a pulsed quadratic magnetic field, which we expect will be the most natural implementation in analog quantum platforms (the inherent Trotterization of a kicking Hamiltonian may also naturally lend itself to digital quantum simulation [53]). We remark that such a distinction is only manifest at the level of micromotion [54]. For the stroboscopic dynamics and Floquet level statistics analyzed in Fig. S8, the problems are identical and simply require a relabeling of parameters.