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

Alec Cao  Roshan Sajjad  Hector Mas  Ethan Q. Simmons  Jeremy L. Tanlimco
Eber Nolasco-Martinez  Toshihiko Shimasaki  H. Esat Kondakci  Department of
Physics, University of California, Santa Barbara, California 93106, USA
Victor Galitski  Joint Quantum Institute and Department of Physics, University
of Maryland, College Park, MD 20742, USA  David M. Weld  Department of
Physics, University of California, Santa Barbara, California 93106, USA

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 incompletely-understood 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.

Fig. 1: Experimentally realizing an interacting quantum kicked rotor. (a)
Schematic of BEC in single pulsed optical lattice. (b) Experimental sequence.
After setting the scattering length the trap is removed and kicking is applied
with period $T$, pulse width $\tau$, and amplitude $V_{0}$ 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.  Fig.
2: Observing the interaction-induced emergence of quantum chaos. (a) Energy
versus kick number for varying $a$. Blue horizontal shaded region indicates
the measured single-rotor localization energy of
$E_{\mathrm{loc}}=2.5(4)E_{\mathrm{R}}$. Here $V_{0}=64E_{\mathrm{R}}$,
$T=1.2\,\upmu$s and $\tau=300$ ns ($K\approx 2.3$ and $\mathchar
22\relax\mkern-9.0muk\approx 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 $\left[-T/4,T/4\right]$.
The solid curve is noninteracting quantum theory and the dotted line is a
diffusion curve $4Dn/\mathchar 22\relax\mkern-9.0muk^{2}$ with $D\approx 0.19$
extracted from the classical standard map [14]. The red dot-dashed line is a
subdiffusive $\sqrt{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_{\mathrm{loc}}=\sqrt{E_{\mathrm{loc}}}\approx
1.6(1)k_{\mathrm{L}}$ [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_{\mathrm{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_{\mathrm{loc}}$ [13].

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 space [18,
19]. Although the QKR is a natural context for experimental probes of the
interplay between many-body 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 7Li 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 $\tau$ (see Fig. 1). We report
momentum and energy in units of $k_{\mathrm{L}}=\pi/d$ and
$E_{\mathrm{R}}=\hbar^{2}k_{\mathrm{L}}^{2}/2m$ with $m$ the mass of 7Li. The
single-particle QKR is defined by the 1-cycle Floquet map $U=e^{-i\mathchar
22\relax\mkern-9.0mukk^{2}/2}e^{-iK\cos z/\mathchar 22\relax\mkern-9.0muk}$
describing a sharp cosine potential impulse followed by free evolution. Here
$k$ and $z$ are momentum and position respectively, $K=\mathchar
22\relax\mkern-9.0mukV_{0}\tau/2\hbar$ is the stochasticity parameter
characterizing kicking strength and $\mathchar
22\relax\mkern-9.0muk=8E_{\mathrm{R}}T/\hbar$ 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
760$a_{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 240$a_{0}$ data
exhibit a clear steady-state behavior for 100 kicks, the 760$a_{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 $\sqrt{n}$
energy growth, and fitting the late-time data to $n^{\alpha}$ yields anomalous
diffusion exponents $\alpha$ in the range $\left[0.4,0.6\right]$. For
reference, 1D Gross-Pitaevskii simulations on a ring [22] predict
$\alpha\in\left[0.5,0.8\right]$, 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
$\alpha\in\left[0.3,0.4\right]$ [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?

Fig. 3: Effect of interactions on reversibility in Loschmidt echo experiments.
(a) Measured Loschmidt echo fidelity $F$ for a range of scattering lengths
$a=\left[0,1500\right]a_{0}$ for $N=6$ (blue circles) and $N=10$ (purple
triangles), where $N$ indicates the total number of kicks; a first set of
$N/2$ kicks propagates the system forward in time and a second time-reversal
set of $N/2$ kicks propagates it backwards. (b) Measured fidelity $F$ at three
different interaction strengths as a function of total number of kicks $N$ in
a Loschmidt echo experiment. (c-e) Averaged absorption images of a BEC after
the first $n$ kicks of an $N=10$ Loschmidt echo protocol, for three different
$a$.

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 for $\mathchar 22\relax\mkern-9.0muk=2\pi q$ with $q$
rational; in particular for $q=2$ ($T\approx 9.95$ $\upmu$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 [33, 34,
13]. This procedure would create exact time reversal for a single zero-
quasimomentum state in the absence of interactions. Due to finite
quasimomentum spread and non-reversed interactions, the reversal is imperfect,
yielding a Loschmidt fidelity
$F=\absolutevalue{\bra{\psi}U_{2}^{\dagger}U_{1}\ket{\psi}}^{2}$ where $U_{1}$
and $U_{2}$ are time-evolution 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 many-body 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. Funding:
D.M.W. acknowledges support from the Air Force Office of Scientific Research
(AFOSR FA9550-20-1-0240), the Army Research Office (PECASE W911NF1410154), and
the National Science Foundation (CAREER 1555313 and QLCI OMA-2016245). D.M.W.,
R.S., and E.N-M. acknowledge support from the UCSB NSF Quantum Foundry through
the Q-AMASE-i program (Grant No. DMR-1906325). V.G. was supported by US-ARO
Contract No.W911NF1310172, NSF DMR-2037158, and the Simons Foundation.
Authors’ contributions: A.C., R.S., H.M., E.Q.S., J.L.T., E.N-M., T.S. and
H.E.K. contributed to operating the experiment and performing the
measurements. A.C., R.S. and H.M. analyzed the data. A.C. conceptualized and
performed the theoretical simulations of rotor and spin models. V.G. and
D.M.W. developed the idea for the experiment. D.M.W. supervised the work.
A.C., R.S., H.M., V.G. and D.M.W. wrote the manuscript. All authors
contributed to discussion and interpretation of the results. Competing
interests: The authors declare no competing interests. Data and materials
availability: All data needed to evaluate the conclusions in this study are
presented in the paper and in the supplementary materials.

## References

  * Nandkishore and Huse [2015] R. Nandkishore and D. A. Huse, Annu. Rev. Condens. Matter Phys. 6, 15 (2015).
  * Abanin  _et al._ [2019] D. A. Abanin, E. Altman, I. Bloch, and M. Serbyn, Rev. Mod. Phys. 91, 021001 (2019).
  * Bernien  _et al._ [2017] H. Bernien, S. Schwartz, A. Keesling, H. Levine, A. Omran, H. Pichler, S. Choi, A. S. Zibrov, M. Endres, M. Greiner, V. Vuletić, and M. D. Lukin, Nature 551, 579 (2017).
  * Serbyn  _et al._ [2021] M. Serbyn, D. A. Abanin, and Z. Papić, Nature Physics 17, 675 (2021).
  * Singh  _et al._ [2019] K. Singh, C. J. Fujiwara, Z. A. Geiger, E. Q. Simmons, M. Lipatov, A. Cao, P. Dotti, S. V. Rajagopal, R. Senaratne, T. Shimasaki, M. Heyl, A. Eckardt, and D. M. Weld, Phys. Rev. X 9, 041021 (2019).
  * Rubio-Abadal  _et al._ [2020] A. Rubio-Abadal, M. Ippoliti, S. Hollerith, D. Wei, J. Rui, S. L. Sondhi, V. Khemani, C. Gross, and I. Bloch, Phys. Rev. X 10, 021044 (2020).
  * Peng  _et al._ [2021] P. Peng, C. Yin, X. Huang, C. Ramanathan, and P. Cappellaro, Nat. Phys. (2021).
  * Langen  _et al._ [2015] T. Langen, S. Erne, R. Geiger, B. Rauer, T. Schweigler, M. Kuhnert, W. Rohringer, I. E. Mazets, T. Gasenzer, and J. Schmiedmayer, Science 348, 207 (2015).
  * Ueda [2020] M. Ueda, Nature Reviews Physics 2, 669 (2020).
  * Agarwal  _et al._ [2015] K. Agarwal, S. Gopalakrishnan, M. Knap, M. Müller, and E. Demler, Phys. Rev. Lett. 114, 160401 (2015).
  * Weidinger and Knap [2017] S. A. Weidinger and M. Knap, Sci. Rep. 7, 1 (2017).
  * Menu and Roscilde [2020] R. Menu and T. Roscilde, Phys. Rev. Lett. 124, 130604 (2020).
  * [13] Materials and methods are available as supplementary materials.
  * Chirikov [1979] B. V. Chirikov, Phys. Rep. 52, 263 (1979).
  * Moore  _et al._ [1995] F. L. Moore, J. C. Robinson, C. F. Bharucha, B. Sundaram, and M. G. Raizen, Phys. Rev. Lett. 75, 4598 (1995).
  * Chabé  _et al._ [2008] J. Chabé, G. Lemarié, B. Grémaud, D. Delande, P. Szriftgiser, and J. C. Garreau, Phys. Rev. Lett. 101, 255702 (2008).
  * Gadway  _et al._ [2013] B. Gadway, J. Reeves, L. Krinner, and D. Schneble, Phys. Rev. Lett. 110, 190401 (2013).
  * Fishman  _et al._ [1982] S. Fishman, D. R. Grempel, and R. E. Prange, Phys. Rev. Lett. 49, 509 (1982).
  * Grempel  _et al._ [1984] D. R. Grempel, R. E. Prange, and S. Fishman, Phys. Rev. A 29, 1639 (1984).
  * Shepelyansky [1993] D. L. Shepelyansky, Phys. Rev. Lett. 70, 1787 (1993).
  * Gligorić  _et al._ [2011] G. Gligorić, J. D. Bodyfelt, and S. Flach, Europhys. Lett. 96, 30004 (2011).
  * Lellouch  _et al._ [2020] S. Lellouch, A. Rançon, S. De Bièvre, D. Delande, and J. C. Garreau, Phys. Rev. A 101, 043624 (2020).
  * Russomanno  _et al._ [2021] A. Russomanno, M. Fava, and R. Fazio, Phys. Rev. B 103, 224301 (2021).
  * Rajak  _et al._ [2019] A. Rajak, I. Dana, and E. G. Dalla Torre, Phys. Rev. B 100, 100302 (2019).
  * Rylands  _et al._ [2020] C. Rylands, E. B. Rozenbaum, V. Galitski, and R. Konik, Phys. Rev. Lett. 124, 155302 (2020).
  * Keser  _et al._ [2016] A. C. Keser, S. Ganeshan, G. Refael, and V. Galitski, Phys. Rev. B 94, 085120 (2016).
  * Oskay  _et al._ [2003] W. H. Oskay, D. A. Steck, and M. G. Raizen, Chaos, Solitons and Fractals 16, 409 (2003).
  * Pikovsky and Shepelyansky [2008] A. S. Pikovsky and D. L. Shepelyansky, Phys. Rev. Lett. 100, 094101 (2008).
  * Flach  _et al._ [2009] S. Flach, D. O. Krimer, and C. Skokos, Phys. Rev. Lett. 102, 024101 (2009).
  * Vuatelet and Rançon [2021] V. Vuatelet and A. Rançon, Effective thermalization of a many-body dynamically localized bose gas (2021), arXiv:2103.14388 [cond-mat.quant-gas] .
  * Rozenbaum  _et al._ [2017] E. B. Rozenbaum, S. Ganeshan, and V. Galitski, Phys. Rev. Lett. 118, 086801 (2017).
  * Gorin  _et al._ [2006] T. Gorin, T. Prosen, T. H. Seligman, and M. Žnidarič, Phys. Rep. 435, 33 (2006).
  * Ullah and Hoogerland [2011] A. Ullah and M. D. Hoogerland, Phys. Rev. E 83, 046218 (2011).
  * Martin  _et al._ [2008] J. Martin, B. Georgeot, and D. L. Shepelyansky, Phys. Rev. Lett. 100, 044106 (2008).
  * Lepers  _et al._ [2008] M. Lepers, V. Zehnlé, and J. C. Garreau, Phys. Rev. A 77, 043628 (2008).
  * Swingle  _et al._ [2016] B. Swingle, G. Bentsen, M. Schleier-Smith, and P. Hayden, Phys. Rev. A 94, 040302 (2016).
  * Toh  _et al._ [2021] J. H. S. Toh, K. C. McCormick, X. Tang, Y. Su, X.-W. Luo, C. Zhang, and S. Gupta, Observation of many-body dynamical delocalization in a kicked ultracold gas (2021), arXiv:2106.13773 [cond-mat.quant-gas] .
  * Pollack  _et al._ [2009] S. E. Pollack, D. Dries, M. Junker, Y. P. Chen, T. A. Corcovilos, and R. G. Hulet, Phys. Rev. Lett. 102, 090402 (2009).
  * Deuar and Drummond [2007] P. Deuar and P. D. Drummond, Phys. Rev. Lett. 98, 120402 (2007).
  * Tenart  _et al._ [2020] A. Tenart, C. Carcy, H. Cayla, T. Bourdel, M. Mancini, and D. Clément, Phys. Rev. Res. 2, 013017 (2020).
  * Gadway  _et al._ [2009] B. Gadway, D. Pertot, R. Reimann, M. G. Cohen, and D. Schneble, Opt. Express 17, 19173 (2009).
  * Ovchinnikov  _et al._ [1999] Y. B. Ovchinnikov, J. H. Müller, M. R. Doery, E. J. D. Vredenbregt, K. Helmerson, S. L. Rolston, and W. D. Phillips, Phys. Rev. Lett. 83, 284 (1999).
  * Klappauf  _et al._ [1999] B. Klappauf, W. Oskay, D. Steck, and M. Raizen, Physica D 131, 78 (1999).
  * Fishman  _et al._ [2003] S. Fishman, I. Guarneri, and L. Rebuzzini, J. Stat. Phys. 110, 911 (2003).
  * Dana and Dorofeev [2006] I. Dana and D. L. Dorofeev, Phys. Rev. E 73, 026206 (2006).
  * D’Alessio and Rigol [2014] L. D’Alessio and M. Rigol, Phys. Rev. X 4, 041048 (2014).
  * Boness  _et al._ [2006] T. Boness, S. Bose, and T. S. Monteiro, Phys. Rev. Lett. 96, 187201 (2006).
  * Boness  _et al._ [2010] T. Boness, K. Kudo, and T. S. Monteiro, Phys. Rev. E 81, 046201 (2010).
  * Rigol  _et al._ [2008] M. Rigol, V. Dunjko, and M. Olshanii, Nature 452, 854 (2008).
  * Monroe  _et al._ [2021] C. Monroe, W. C. Campbell, L.-M. Duan, Z.-X. Gong, A. V. Gorshkov, P. W. Hess, R. Islam, K. Kim, N. M. Linke, G. Pagano, P. Richerme, C. Senko, and N. Y. Yao, Rev. Mod. Phys. 93, 025001 (2021).
  * Browaeys and Lahaye [2020] A. Browaeys and T. Lahaye, Nat. Phys. 16, 132 (2020).
  * Salathé  _et al._ [2015] Y. Salathé, M. Mondal, M. Oppliger, J. Heinsoo, P. Kurpiers, A. Potočnik, A. Mezzacapo, U. Las Heras, L. Lamata, E. Solano, S. Filipp, and A. Wallraff, Phys. Rev. X 5, 021027 (2015).
  * Sieberer  _et al._ [2019] L. M. Sieberer, T. Olsacher, A. Elben, M. Heyl, P. Hauke, F. Haake, and P. Zoller, npj Quantum Inf. 5, 1 (2019).
  * Goldman and Dalibard [2014] N. Goldman and J. Dalibard, Phys. Rev. X 4, 031027 (2014).

Methods
Experimental platform and sequence. The experiments begin with a Bose-Einstein
condensate (BEC) of around $10^{5}$ 7Li atoms in a far-detuned optical dipole
trap with trapping frequencies $\omega_{x,z}/2\pi\approx 40$ Hz and
$\omega_{y}/2\pi\approx 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_{\mathrm{L}}=\pi/d$, and recoil energy
$E_{\mathrm{R}}=\hbar^{2}k_{\mathrm{L}}^{2}/2m$ with $m$ the mass of 7Li. The
full dynamics are then well described by the second-quantized Hamiltonian

$\begin{split}\mathcal{H}=\int\mathrm{d}^{3}r\,\hat{\Psi}^{\dagger}(\mathbf{r},t)H(\mathbf{r},t)\hat{\Psi}(\mathbf{r},t)+\\\
\frac{g}{2}\int\mathrm{d}^{3}r\,\hat{\Psi}^{\dagger}(\mathbf{r},t)\hat{\Psi}^{\dagger}(\mathbf{r},t)\hat{\Psi}(\mathbf{r},t)\hat{\Psi}(\mathbf{r},t)\end{split}$
(1) $\displaystyle
H(\mathbf{r},t)=\frac{p^{2}}{2m}+\frac{V_{0}}{2}\cos(2k_{L}z)I(x,y)\sum_{n}f_{\tau}(t-nT).$
(2)

The key kick parameters are the lattice depth $V_{0}$, effective pulse width
$\tau$, and kick period $T$. $V_{0}$ is calibrated through a standard Kapitza-
Dirac diffraction technique. $f_{\tau}(t)$ denotes a unit amplitude pulse
function beginning at $t=0$ of width $\tau$. 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 $\tau=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\pi\hbar^{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)\approx e^{-2(x^{2}+y^{2})/\sigma^{2}}$
with a measured 1/$e^{2}$ beam radius of $\sigma\approx 65$ $\upmu$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 7Li 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 $\pm 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 post-selection
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_{\mathrm{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\left(\langle z^{2}\rangle+2\langle
x^{2}\rangle\right)/2t_{\mathrm{TOF}}^{2}$ with $t_{\mathrm{TOF}}\approx 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_{\mathrm{loc}}$ is estimated by
averaging the noninteracting trace for $n\geq 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(-\absolutevalue{k}/k_{\mathrm{loc}})\sum_{j}\exp(-(k-2k_{\mathrm{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_{\mathrm{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_{\mathrm{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_{\mathrm{loc}})/(1-\exp(-k_{0}/k_{\mathrm{loc}}))^{2}$;
here $k_{0}\approx 9.85k_{L}$ is the half-width of our images which yields a
negligible correction factor of $\approx 1.006$. The width of the region is
computed through linearized error propagation.

In Fig. 2i, the plotted localization metric is
$\int_{-k_{0}}^{k_{0}}\max[r(k)-1,0]\differential k/2k_{0}$. Here,
$r(k)=\absolutevalue{\psi(k)}^{2}/\exp(-k/k_{\mathrm{loc}})$ is the ratio of
the measured axial density denoted $\absolutevalue{\psi(k)}^{2}$ and an
exponential localization envelope. Here the maximum of
$\absolutevalue{\psi(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_{\mathrm{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}\approx 50E_{\mathrm{R}}$, $\tau=300$ ns, and $T=9.93$ $\upmu$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^{\prime}$ 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 $\sim 0.1k_{\mathrm{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).

Supplementary Materials

###### Contents

  1. 1 Additional Delocalization Data
  2. 2 Monte Carlo distributions for exponential localization deviation
  3. 3 Systematics: Finite pulse width
  4. 4 Systematics: Position space dynamics and TOF conversion
  5. 5 Systematics: Beam-induced transverse dynamics
  6. 6 Systematics: Effects of scattering on measured Loschmidt fidelity
  7. 7 Mapping the QKR to a kicked spin chain model

##  1 Additional Delocalization Data

Fig. S1: Interaction-induced delocalization for a different set of kicking
parameters. The kick parameters are $V_{0}=70\,E_{\mathrm{R}}$, $T=2.2$
$\upmu$s and $\tau=300$ ns ($K\approx 4.6$ and $\mathchar
22\relax\mkern-9.0muk\approx 2.8$). The (a) energy, (b) 1D momentum-space IPR
and (c) deviation from exponential localization over time for varying
scattering lengths. In a, the shaded regions indicate the extracted initial
localization lengths for the three interaction strengths which we use for
computing c.

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$ $\upmu$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.

##  2 Monte Carlo distributions for exponential localization deviation

Fig. S2: Characterizing Monte Carlo distributions for deviation from
exponential localization. (a) Distribution for the 760$a_{0}$ data in Fig. 2
at $n=500$. Orange line indicates the fit to a log-normal distribution.
Vertical dashed line indicates the mean, and the vertical dotted lines
surrounding it indicate the interquartile range reported as the errorbars. (b)
Evolution of the distribution over time for noninteracting and interacting
samples. Note the difference in y-scale. The $n=42$ trace in the lower panel
is cut-off vertically for visual clarity on the larger $n$ distributions. (c)
Alternative characterizations of the exponential localization deviation in
terms of the mode and median of the simulated distributions (indicated by the
markers, the errorbars are left as the interquartile range), as opposed to the
mean shown in the main text Fig. 2i.

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_{\mathrm{loc}}$ drawn from a Gaussian centered at
1.58$k_{\mathrm{L}}$ and with standard deviation 0.12$k_{\mathrm{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 760$a_{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.

##  3 Systematics: Finite pulse width

Fig. S3: Effects of finite pulse width on noninteracting QKR localization. (a)
Time evolution of energy for $K=2.3$. The $\tau=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$.

The delta-function kicking assumption in the theoretical QKR model is not
perfectly realized in experiment owing to the finite atomic mass of 7Li. The
assumption corresponds to the Raman-Nath diffraction regime which is
approximately expressed by the condition
$2\sqrt{V_{0}E_{\mathrm{R}}}\tau/\hbar\ll 1$ [41]. For the experiment with
$\tau=0.3$ $\upmu$s and $V_{0}=64E_{\mathrm{R}}$ as in Fig. 2, this parameter
is approximately 0.76 (there is an ambiguity of a factor $2\pi$ 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 $\tau$ to the delta-kick Floquet map solution. To make this
comparison, we keep the effective stochasticity parameter $K\sim V_{0}\tau$
characterizing the kicking strength constant as we let $\tau\to 0$. 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\approx 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 $\pm 33.2k_{\mathrm{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.

##  4 Systematics: Position space dynamics and TOF conversion

Fig. S4: Effects of position-space dynamics on extracting energy from TOF
images. Noninteracting numerical simulations comparing different methods for
measuring energy in TOF and the true energy of the kicked rotor. (a) The
kicking parameters are comparable to the ones used in the main text
($T=1.2\,\upmu$s and $\tau=300\,$ns), but with lattice depth reduced to
$V_{0}=38E_{\mathrm{R}}$ 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\,\upmu$s roughly
corresponding to the supplementary delocalization data in Fig. S1 (with
adjusted lattice depth $V_{0}=52E_{\mathrm{R}}$). 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.

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 7Li 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 and b, revealing that in fact the simplest approach of using the TOF
to convert position to momentum works well for the parameters of the
experiments reported in the main text. For experiments at other parameter
values, however, the appropriate 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_{\mathrm{L}}$ momentum, though the added complexity
of this scheme is not justified by the results.

##  5 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.

Fig. S5: Transverse dynamics in experimental quantum kicked rotors. Measured
transverse kinetic energy versus time for the main delocalization data in Fig.
2. Each interaction strength undergoes an oscillation due to harmonic motion
in the time-averaged lattice potential. The difference in evolution among
different interaction strengths is due to scattering effects.

##  6 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 $\pm 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.

Fig. S6: Comparison of atom counting methods for Loschmidt fidelity. (A) Zero
mode fraction measured at the halfway point (F’) and at the end (F) of an
8-kick Loschmidt echo protocol for both raw and Gaussian counting methods. (B)
The difference in zero mode fraction between F and F’ in an 8-kick Loschmidt
echo protocol for both methods. Fig. S7: Experimentally measured scattered
fraction during Loschmidt experiments. (A) TOF absorption images of a BEC
after an $N=8$ Loschmidt protocol at various scattering lengths. (B)
Corresponding scattered fraction as computed using the Gaussian fitting
described in the methods section 1.3.

##  7 Mapping the QKR to a kicked spin chain model

Fig. S8: Kicked-rotor dynamics and many-body dynamical localization in kicked
spin chains. (a-b) Calculated single-spin-flip evolution for static and kicked
spin-1/2 XXX chains of length $L=13$ for $K=2.4$, $\mathchar
22\relax\mkern-9.0muk=1.5$ and $\beta=0$. (c) Spin chain quantum resonance
spectroscopy (solid) demonstrating equivalence to QKR (dashed) in the single
spin-flip sector with $K=3\mathchar 22\relax\mkern-9.0muk$, $L=51$ and
averaging over a Gaussian ensemble of $\beta$. Energy is after 10 kicks.
Zoomed-in panel compares experimental measurements on the atomic quantum
kicked rotor to predictions of spin chain numerics. Dashed vertical lines
indicate $q=$ 1/3, 3/8, 2/5, and 1/2. (d-e) The same as a-b but with 2 spin
flips. (f) Time-averaged staggered magnetization versus kick number, starting
from an initial Néel state with $\mathchar 22\relax\mkern-9.0muk=1$, $L=12$,
$\beta=0.1$ and varying $K$. (g) Staggered magnetization in the infinite-time
limit versus $K$, for the same parameters as panel f. (h) Gap-ratio statistic
of the $M=0$ sector at $\mathchar 22\relax\mkern-9.0muk=1$ for varying $L$
averaged over 100 values of $\beta$. Dashed lines indicate predictions of the
Poisson ($\langle r\rangle\approx 0.386$) and circular-orthogonal ensemble
($\langle r\rangle\approx 0.527$) [46].

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_{\mathrm{kickedXXZ}}=e^{-i\mathchar
22\relax\mkern-9.0mukH_{\mathrm{quad}}/4}e^{-iKH_{\mathrm{XXZ}}/4\mathchar
22\relax\mkern-9.0muk}$, with
$H_{\mathrm{quad}}=\sum_{j}(j+\beta)^{2}\sigma_{j}^{z}$ a quadratic field and
$H_{\mathrm{XXZ}}=\sum_{j}(\sigma_{j}^{x}\sigma_{j+1}^{x}+\sigma_{j}^{y}\sigma_{j+1}^{y}+\Delta\sigma_{j}^{z}\sigma_{j+1}^{z})$
an XXZ Hamiltonian with $\sigma^{x,y,z}$ the Pauli matrices and $\Delta$ the
anisotropy parameter. $\beta$ represents a field-center offset which serves as
an analogue to quasimomentum in the standard QKR problem. Total magnetization
$M=\langle\sum_{j}\sigma_{j}^{z}\rangle$ 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_{\mathrm{kickedXXZ}}$ on systems of length up to $L=14$. We center the
chains about $j=0$ (i.e. $j\in[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 ($\Delta=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_{\mathrm{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 $\ket{\uparrow\downarrow\uparrow\downarrow\ldots}$ 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=\sum_{j}(-1)^{j}\sigma_{j}^{z}/L$ after $n$ cycles, given by $\langle
O\rangle_{n}=\left(\sum_{i=1}^{n}\bra{\psi}U_{\mathrm{kickedXXZ}}^{\dagger
i}OU_{\mathrm{kickedXXZ}}^{i}\ket{\psi}\right)/n$. The infinite-time limit is
then calculated via the Floquet diagonal ensemble as $\lim_{n\to\infty}\langle
O\rangle_{n}=\sum_{\alpha}\absolutevalue{c_{\alpha}}^{2}\bra{\psi_{\alpha}}O\ket{\psi_{\alpha}}$.
Here $c_{\alpha}=\bra{\psi_{\alpha}}\ket{\psi}$ are the coefficients of the
initial state $\ket{\psi}$ in the basis of the many-body Floquet states
$\ket{\psi_{\alpha}}$ [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 $\langle
r\rangle$ [46, 13] (Fig. S8h). The gap ratio is defined as
$r_{\alpha}=\mathrm{min}\left(\delta_{\alpha},\delta_{\alpha+1}\right)/\mathrm{max}\left(\delta_{\alpha},\delta_{\alpha+1}\right)$,
where $\delta_{\alpha}=\epsilon_{\alpha+1}-\epsilon_{\alpha}$ is the gap
between consecutive quasi-energies $\epsilon_{\alpha}$, which we order on the
interval $[-\pi,\pi]$. To compute $\langle r\rangle$, we average $r_{\alpha}$
over the $M=0$ sector as well as for 100 values of $\beta$ drawn from a normal
distribution of standard deviation 0.1. We interpret the transition in
$\langle r\rangle$ 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 $\Delta=1$ point, as shown in Fig. S9. In the limit $\Delta=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$ and $\mathchar
22\relax\mkern-9.0muk$ in $U$ to the spin model $U_{\mathrm{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.

Fig. S9: Transition from MBDL to ergodicity for kicked XXZ model. (a) Gap-
ratio statistic for $\mathchar 22\relax\mkern-9.0muk=1$, $\Delta=0.5$ and
averaging over a Gaussian ensemble of $\beta$ with standard deviation 0.1 for
varying $K$ and $L$. (b) The same parameters except $\Delta=2$. Dashed
horizontal lines indicate the same Poisson and circular-orthogonal ensemble
predictions of Fig. S8h.