# A first-principles study on the lattice thermal conductivity of irradiated
glassy states of the Ge2Sb2Te5 phase-change memory material

Felix C. Mocanu<EMAIL_ADDRESS>Department of Chemistry, University of
Cambridge, Lensfield Road, CB2 1EW, Cambridge, United Kingdom  Konstantinos
Konstantinou  Department of Chemistry, University of Cambridge, Lensfield
Road, CB2 1EW, Cambridge, United Kingdom  Stephen R. Elliott  Department of
Chemistry, University of Cambridge, Lensfield Road, CB2 1EW, Cambridge, United
Kingdom

###### Abstract

An analysis of thermal transients from non-equilibrium _ab initio_ molecular-
dynamics simulations can be used to calculate the thermal conductivity of
materials with a short phonon mean-free path. We adapt the approach-to-
equilibrium methodology to the three-dimensional case of a simulation that
consists of a cubic core region at higher temperature approaching thermal
equilibrium with a thermostatted boundary. This leads to estimates of the
lattice thermal conductivity for the glassy state of the phase-change memory
material, Ge2Sb2Te5, which are close to previously reported experimental
measurements. Self-atom irradiation of the material, modelled using thermal
spikes and stochastic-boundary conditions, results in glassy models with a
significant reduction of the lattice thermal conductivity compared to the
pristine glassy structure. This approach may prove to be useful in
technological applications, e.g. for the suppression of thermal cross-talk in
phase-change memory and data-storage devices.

The cost of computation is significantly determined by the energy needed to
keep key components cool, whether it is central processing units,
accelerators, memory or data storage. Even when components are kept at
reasonable temperatures during normal operation, there may still be thermal-
throttling under peak-load which affects their performance and lifetime.
Hence, the ability to predict the thermal properties of solids is critical for
the rational design of materials and the management of devices that are
sensitive to thermal fluctuations. Modestino and Hausssener (2015); Zebarjadi
(2016)

Phase-change memory devices, which are widely believed to be a part of at
least some commercial non-volatile solid-state drives, rely on the fast and
reversible switching between a conductive crystalline phase (the “1” state of
a bit), and a resistive glassy phase (the “0” state of a bit), of the same
chalcogenide material, induced via Joule heating from the application of
appropriate voltage pulses. Elliott (2015) Consequently, thermal cross-talk
between adjacent memory cells can limit the size down-scaling in phase-change
memory devices. Fong, Neumann, and Wong (2017) An understanding of the
thermal-transport properties of these materials at the atomic level is
therefore essential for their design and optimization. Lencer, Salinga, and
Wuttig (2011)

The lattice thermal conductivity of the phase-change memory material, GeTe,
has been predicted successfully with quasi-static calculations, Sosso  _et
al._ (2012) as well as with classical equilibrium, Campi  _et al._ (2015) and
non-equilibrium, Campi  _et al._ (2015); Sosso  _et al._ (2018) molecular-
dynamics methods by using a linear-scaling neural-network interatomic
potential. Sosso  _et al._ (2012) First-principles calculations were also
employed to estimate the thermal conductivity for the crystalline phases of
different chalcogenide phase-change materials, including Ge2Sb2Te5, Campi  _et
al._ (2017) in good agreement with experimental observations. However, there
are no previously reported modelling studies of the lattice thermal
conductivity for the glassy phase of Ge2Sb2Te5 which is of equal importance
with the crystalline counterpart for technological applications.

Methods for simulating thermal transport at the atomic level have seen a rapid
evolution and are getting closer to quantitative agreement with experimental
measurements. Baroni  _et al._ (2018); Sosso  _et al._ (2018) _Ab initio_
methods, based on the Boltzmann transport equation, Broido  _et al._ (2007);
Carg  _et al._ (2011) or on Green-Kubo dynamical formulations, Marcolongo,
Umari, and Baroni (2016); Carbogno, Ramprasad, and Scheffler (2017) have been
employed recently in the literature to provide a first-principles description
of thermal transport, and they represent significant advances in the field.
Nevertheless, these approaches are computationally very demanding within a
density-functional-theory framework, while some of them are specifically
designed for harmonic solids near equilibrium.

Therefore, there is the necessity for an efficient first-principles molecular-
dyanamics approach to model the thermal properties of glasses and in
particular the glassy phase of Ge2Sb2Te5. An efficient and quantitative
assessment of lattice thermal conductivity can be obtained from the approach-
to-equilibrium molecular-dynamics (AEMD) method, Lampin  _et al._ (2012, 2013)
which has been successfully deployed in _ab initio_ (as well as with empirical
force fields) molecular-dynamics simulations of several different materials.
Puligheddu, Gygi, and Galli (2017)

In this Letter, the AEMD methodology, which belongs to a larger class of non-
equilibrium molecular dynamics (NEMD) methods, has been adapted to the case of
a cubic core region in contact with a thermostatted boundary shell in order to
calculate the thermal conductivity of irradiated glassy Ge2Sb2Te5. Below, we
briefly describe the computational details of the NEMD simulation protocol
using stochastic-boundary conditions that is used to simulate the energetic
thermal spikes. Numerical results are presented, starting with the fitting of
thermal transients, the predicted thermal conductivity and the effect of
irradiation on the thermal properties. We also compare our findings to
experimental data and previous simulation studies, and we discuss some of the
limitations of the approach.

The simulated system is a 315-atom melt-quenched model of glassy Ge2Sb2Te5.
Lee and Elliott (2017) The simulation box, which has a length of 21.65 Å, was
divided into a core cubic region and an outer boundary shell of thickness 1 Å
on each side. This type of separation has been dubbed stochastic-boundary
conditions, and was originally used to investigate thermal transport at
interfaces. Tenenbaum, Ciccotti, and Gallico (1982) The core region samples a
micro-canonical ensemble (NVE) while the boundary undergoes Langevin dynamics
in the canonical ensemble (NVT) and dissipates the heat generated in the core
region, during an ionic cascade. Toton  _et al._ (2010) The thermostat
parameters were chosen such that the damping period was 100 fs. _Ab initio_
molecular-dynamics simulations were carried out using the CP2K code,
Vandevondele  _et al._ (2005) in which the stochastic-boundary-conditions
approach is implemented, based on the Generalized Langevin Equation
formulation. Kantorovich and Rompotis (2008) We modelled radiation-induced
non-equilibrium cascades by performing thermal-spike simulations with initial
kinetic energies in the range of $15-200$ eV. Further details related to the
computational set-up, the radiation-damage cascades and the _ab initio_
molecular-dynamics simulations can be found in our previous work, Konstantinou
_et al._ (2018a, b) where some of the initial trajectories used in this work
have been taken from.

In the approach-to-equilibrium methodology, heat conduction is usually
modelled by Fourier’s law. Lampin  _et al._ (2013) Based on our set-up, the
simulation box is cubic and is assumed to be approximately homogeneous and
isotropic. Additionally, we have considered that the boundary region acts as a
“thermal wall” at a temperature of $T_{0}$ ($300$ K) and the core region will
rapidly come into equilibrium with it after the cascade. In practice, the
temperature of the boundary region will oscillate significantly and there will
be an artificial thermal boundary resistance at the interface with the core
region due to the thermostat that is employed in the outer shell. Kapitza
(1941); Tenenbaum, Ciccotti, and Gallico (1982); Singh and Tadmor (2015) The
Cartesian coordinates of atoms in the core region, $x$, $y$ and $z$, reside in
the real interval $[0,L]$ where $L=a-2r$, $a$ is the size of the periodic
cubic simulation box and $r$ is the thickness of the boundary region. Under
these assumptions, the resulting heat equation can be written as:

$\dfrac{\partial T}{\partial t}=-\alpha\left(\dfrac{\partial^{2}T}{\partial
x^{2}}+\dfrac{\partial^{2}T}{\partial y^{2}}+\dfrac{\partial^{2}T}{\partial
z^{2}}\right)$ (1)

where $T$ is the temperature and $t$ is the time.

The thermal diffusivity $\alpha$ is defined as:

$\alpha=\dfrac{\kappa}{C_{v}\rho}$ (2)

where $\kappa$ is the thermal conductivity, $C_{v}$ the constant-volume heat
capacity and $\rho$ the density of the system. The general solution, assuming
a separable form, is then given by:

$\begin{split}T(x,y,z,t)=T_{0}+\sum_{n=1}^{\infty}\sum_{m=1}^{\infty}\sum_{l=1}^{\infty}a_{nml}e^{-\alpha\lambda_{nml}t}\\\
\sin\left(\dfrac{n\pi x}{L}\right)\sin\left(\dfrac{m\pi
y}{L}\right)\sin\left(\dfrac{l\pi z}{L}\right)\end{split}$ (3)

The Fourier-series coefficients in the general solution, $\lambda_{nml}$ and
$a_{nml}$, can be inferred from the boundary conditions and are given in the
two equations below:

$\lambda_{nml}=\left(\dfrac{\pi}{L}\right)^{2}\left(n^{2}+m^{2}+l^{2}\right)\\\
$ (4)

$\begin{split}a_{nml}=-\Delta
T_{0}\left(\dfrac{2}{\pi}\right)^{3}\dfrac{\left(1-(-1)^{n}\right)\left(1-(-1)^{m}\right)\left(1-(-1)^{l}\right)}{nml}\end{split}$
(5)

The difference between the spatially averaged temperature of the core region
(undergoing NVE dynamics) and the target temperature of the thermostatted
boundary region (undergoing Langevin NVT dynamics) corresponds to: $\Delta
T(t)=\bar{T}(t)-T_{0}$. It has an initial value, $\Delta T_{0}$, at the start
of the thermal quench and it decays to zero during the quench.

$\Delta T\sim e^{-\dfrac{t}{\tau}}$ (6)

The exponential temporal decay of this temperature difference in equation 6,
has a dominant contribution from the leading term $n=m=l=1$ of the Fourier
series in equation 3. Once the relaxation time $\tau$ is obtained from
simulations, it can be inserted into the time-dependent part of this dominant
term. Hence, the thermal conductivity, $\kappa$, can be calculated from the
expression:

$\kappa=\dfrac{L^{2}}{3\pi^{2}}\dfrac{C_{v}\rho}{\tau}$ (7)

The thermal transients of the non-equilibrium ion-irradiation simulations for
glassy Ge2Sb2Te5 can be directly fitted from the thermal quench of the core
region as it reaches thermal equilibrium with the boundary layer. The
trajectory of the radiation-induced cascade can be split into three intervals,
based on the time evolution of the kinetic temperature, shown in figure 1 for
different initial thermal-spike kinetic energies: (a) The high-energy cascade
generated by the thermal spike; (b) An approach-to-equilbrium transient that
is reasonably well described by an exponential temporal decay of the
temperature difference between the core and the boundary; and (c) An
equilibrium region where the system as a whole fluctuates around the target
temperature of the thermostatted boundary layer (300 K).

Figure 1: Temperature (K) versus time (fs), on a linear-log graph, for a model
of glassy Ge2Sb2Te5, for different thermal-spike energies, shown as shades of
red (darker means higher energy). The characteristic dips in temperature
correspond to collision events during the non-equilibrium cascade. This is
followed by an exponential decay of the temperature as the system approaches
equilibrium, while at the end, the kinetic temperature of the system
fluctuates around 300 K.

By examining the approach-to-equilibrium of the system after the high-energy
cascade, thermal-conductivity estimates are obtained “on the fly” without
having to run a separate simulation for this purpose. It is important to only
fit the relaxation time $\tau$ using data from a restricted time interval in
which the temperature decays exponentially. In practice this means including
data only after the shock of the thermal spike has been absorbed by the
boundary and just before reaching equilibrium. In order to avoid including
data from the high-energy cascade at the beginning of the simulation, a time
period of $1-2$ ps from the start of the simulation needs to be removed from
the fitting interval, based on the initial energy of the thermal spike, as
indicated by our kinetic analysis of the approximate cascade duration.
Konstantinou  _et al._ (2018b) Examples of exponential fits of the thermal
relaxation time from the asymptotic regime of the ion-irradiation simulations
are shown in figure 2 for 50 eV, 100 eV and 200 eV initial thermal-spike
energies.

Figure 2: The temperature difference $\Delta T$ between the core region and
the target temperature of the boundary thermostat, viz. 300 K, as a function
of time. Thermal transients in the approach-to-equilibrium are shown in
different shades of red as a function of the energy of the thermal spike
(darker is higher). The corresponding exponential fits are shown as black
lines.

After obtaining the transients for all the thermal-spike simulations, we
examined the thermal conductivity of the initial pristine glassy Ge2Sb2Te5
structure. The computational procedure used in this case comprised the
following steps: (1) Fix the atoms in the boundary region; (2) Initialize and
equilibrate the velocities in the core region at a higher temperature; (3)
Release the constraints for the boundary region; and (4) Remove the thermostat
from the core region. In this case, there is no radiation-induced cascade and
the simulation consists simply of a thermal quench and equilibration with the
boundary. The core-region initial maximum kinetic temperature was chosen to be
700 K in order to avoid any intermixing between the core and the boundary.

The thermal relaxation times corresponding to the different ion-irradiation
simulations and to the pristine glassy Ge2Sb2Te5 structure were calculated,
and are shown in figure 3. It can be observed that the relaxation time $\tau$
scales as a power-law function of the initial energy of the thermal spike. The
temperature evolution is regularly analyzed after simulations of thermal
spikes and a mechanical model predicts that the thermal relaxation time will
scale as a power-law function of the thermal-spike energy, with an exponent
around $2/3$. Marks (1997) However, values for the power-law exponent below
$2/3$ have been reported from computer simulations, suggesting that the
exponent depends on the structure of the material. Buchan  _et al._ (2015)
From our simulations, a best-fit power-law exponent of $0.354$ was obtained
for glassy Ge2Sb2Te5, indicating a gentle increase of the thermal relaxation
time with the energy of the thermal spike. We expect that this trend will not
be significantly changed if the calculation is repeated for several
independent amorphous models or indeed a larger model to accommodate higher
thermal-spike energies.

Figure 3:  Thermal relaxation times versus the initial kinetic energy of the
thermal-spike atom. Data points are shown as black circles. A power-law fit to
the ion-irradiation data is shown as the solid black line and the shaded
region is the 95% confidence interval of the fit.

An additional ingredient for the estimation of the lattice thermal
conductivity is the heat capacity of the material. In order to obtain the heat
capacity, energy fluctuations, $\delta E=E-\langle E\rangle$, were sampled
from equilibrium _ab initio_ molecular-dynamics simulations. The starting
point for each trajectory corresponds to the structural models at the end of
each non-equilibrium thermal-spike simulation. In these subsequent molecular-
dynamics runs, there is no longer a separation between core and boundary, and
the entire system undergoes NVT dynamics with a single Langevin equation
(GLE), or colored-noise, thermostat. Ceriotti, Bussi, and Parrinello (2008,
2010) Trajectories of 40 ps were generated with a time-step of 1 fs for each
glassy sample. The first 8 ps were discarded and the heat capacity was
estimated from the remaining 32000 time steps.

The calculated heat capacity at 300 K for the pristine glassy phase of
Ge2Sb2Te5 was 272 $\pm$ 60.19 J kg-1 K-1, which is above the Dulong-Petit
limit value of 219 J kg-1K-1 ( $3Nk_{B}$). This is in good agreement with
experimental data from differential scanning calorimetry, viz. $220-230$ J
kg-1 K-1, Kuwahara  _et al._ (2007); Zalden  _et al._ (2014) and with the fact
that the Debye temperature of glassy (as-deposited) Ge2Sb2Te5 was found to be
below 111 K from inelastic neutron-scattering experiments. Zalden  _et al._
(2014) In the Supplemental Material, details are provided for the calculation
of the heat capacity and thermal conductivity in all of the simulated systems.

The lattice thermal conductivity of glassy Ge2Sb2Te5 was calculated from
equation 7, using a relaxation time corresponding to the intercept from our
power-law extrapolation. This estimated thermal conductivity was found to be
0.16$\pm$0.04 W K-1 m-1 in good agreement with the experimentally reported
values, which are in the range $0.19-0.3$ W K-1 m-1. Peng, Cheng, and
Mansuripur (1997); Giraud  _et al._ (2005); Lyeo  _et al._ (2006) From figure
3, it can be seen that the irradiated glassy Ge2Sb2Te5 models exhibit an
increased thermal-relaxation time which ultimately leads to a decreased
thermal conductivity. The estimated lattice thermal conductivity of the self-
irradiated configurations was found to be in the range of 0.067 to 0.113 W K-1
m-1, depending on the thermal-spike energy, further revealing that the thermal
conductivity could be significantly lowered by controlled irradiation. It is
noted that controlled irradiation (with He ions) has been used successfully
before to lower the lattice thermal conductivity of Si nanowires. Zhao  _et
al._ (2017)

In a recent AEMD study for a related glassy chalcogenide material, namely
GeTe4, the authors reported a value of 0.013$\pm$0.003 W K-1m-1 for the
thermal conductivity, Bouzid  _et al._ (2017) which is an order of magnitude
lower than the experimental measurement for this material (0.1 W K-1 m-1).
Zhang  _et al._ (2009) This discrepancy is likely due to the small model
system size (185 atoms), which illustrates the limitations of tractable _ab
initio_ molecular-dynamics simulations. The same authors, in a different
study, Martin  _et al._ (2018) doubled the simulated system size (up to 370
atoms) and obtained a value of 0.044$\pm$0.001 W K-1 m-1 for the thermal
conductivity of glassy GeTe4, closer to, but still smaller than, the
experimental value, revealing the influence of size effects.

The relatively low estimated value for the thermal conductivity of glassy
Ge2Sb2Te5 in this work could suggest that the system size is still somewhat
too small to fully account for the contribution of any long-wavelength
vibrational modes to the lattice thermal conductivity. Allen and Feldman
(1993); Sosso  _et al._ (2018) While not a substitute for analysing the
finite-size effects, the influence of the boundary region Langevin-thermostat
damping time was included in our estimates of the lattice thermal conductivity
and is discussed in the Supplemental Material.

In conclusion, a non-equilibrium molecular-dynamics methodology is proposed
for calculating the lattice thermal conductivity of a cubic-core region as it
approaches equilibrium with a thermostatted-boundary layer. This approach has
been applied to estimate the change in the lattice thermal conductivity with
thermal-spike energy in self-irradiated glassy Ge2Sb2Te5 models. As the
irradiated core approached equilibrium with the thermostatted boundary, the
thermal relaxation time is fitted using an appropriate time interval from the
asymptotic regime of the simulation. The good agreement obtained between the
calculated value for the lattice thermal conductivity of pristine glassy
Ge2Sb2Te5 and the results of experimental measurements provides validation for
this approach.

The thermal relaxation time is found to increase as a sublinear power-law
function of the thermal-spike energy. This results in an overall decrease of
the lattice thermal conductivity compared to that of the pristine glassy
structure. Our simulations suggest that radiation-induced cascades can
therefore reduce the lattice thermal conductivity of glassy Ge2Sb2Te5 by as
much as 60%. Given the remarkable recovery of the electronic structure of the
glass after irradiation, as already demonstrated for this material in our
previous work, Konstantinou  _et al._ (2018a) ion irradiation can therefore be
a potential strategy for improving the performance of phase-change memory and
data-storage devices by reducing thermal cross-talk between memory cells.
While the qualitative trends should stay the same, in the future, these
effects should be explored in larger models which allow for better statistics
and a wider range of higher thermal-spike energies, as well as in models of
the crystalline phases of Ge2Sb2Te5.

See Supplemental Material for details about the thermal profiles established
in the asymptotic regime in the irradiated glass, the calculation of the heat
capacity in the simulated systems and the effect of the Langevin-thermostat
damping time during irradiation on the calculation of the thermal
conductivity.

F.C.M. acknowledges financial support from the UK Engineering and Physical
Sciences Research Council (EPSRC) Centre for Doctoral Training in
Computational Methods for Materials Science under grant EP/L015552/1, and
resources provided by the “Cambridge Service for Data Driven Discovery” (CSD3,
http://csd3.cam.ac.uk) system operated by the University of Cambridge Research
Computing Service (http://www.hpc.cam.ac.uk ) and funded by EPSRC Tier-2
capital grant EP/P020259/1. Via our membership of the UK’s HEC Materials
Chemistry Consortium, which is funded by EPSRC (EP/L000202, EP/R029431), this
work used the ARCHER UK National Supercomputing Service
(http://www.archer.ac.uk). K.K. acknowledges financial support from the EPSRC
grant EP/N022009/1 (“Development and Application of Non-Equilibrium Doping in
Amorphous Chalcogenides”), and also acknowledges the use of the High
Performance Computing Facility (Grace@UCL), and associated support services,
in the completion of this work.

## References

## References

  * Modestino and Hausssener (2015) M. A. Modestino and S. Hausssener, Annu. Rev. Chem. Biomed. Eng. 6, 13 (2015).
  * Zebarjadi (2016) M. Zebarjadi, Sci. Rep. 6, 20951 (2016).
  * Elliott (2015) S. R. Elliott, Int. J. Appl. Glass Sci. 6, 15 (2015).
  * Fong, Neumann, and Wong (2017) S. W. Fong, C. M. Neumann, and H. S. Wong, IEEE Transactions on Electron Devices 64, 4374 (2017).
  * Lencer, Salinga, and Wuttig (2011) D. Lencer, M. Salinga, and M. Wuttig, Adv. Mater. 23, 2030 (2011).
  * Sosso  _et al._ (2012) G. C. Sosso, D. Donadio, S. Caravati, J. Behler, and M. Bernasconi, Phys. Rev. B 86, 104301 (2012).
  * Campi  _et al._ (2015) D. Campi, D. Donadio, G. C. Sosso, J. Behler, and M. Bernasconi, J. Appl. Phys. 117, 015304 (2015).
  * Sosso  _et al._ (2018) G. C. Sosso, V. L. Deringer, S. R. Elliott, and G. Csányi, Mol. Sim. 44, 866 (2018).
  * Campi  _et al._ (2017) D. Campi, L. Paulatto, G. Fugallo, F. Mauri, and M. Bernasconi, Phys. Rev. B 95, 024311 (2017).
  * Baroni  _et al._ (2018) S. Baroni, R. Bertossa, L. Ercole, F. Grasselli, and A. Marcolongo, in  _Handbook of Materials Modeling_ (Springer International Publishing, Cham, 2018) pp. 1–36, arXiv:1802.08006 .
  * Broido  _et al._ (2007) D. A. Broido, M. Malorny, G. Birner, N. Mingo, and D. A. Stewart, Appl. Phys. Lett. 91, 231922 (2007).
  * Carg  _et al._ (2011) J. Carg, N. Bonini, B. Kozinsky, and N. Marzari, Phys. Rev. Lett. 106, 045901 (2011).
  * Marcolongo, Umari, and Baroni (2016) A. Marcolongo, P. Umari, and S. Baroni, Nat. Phys. 12, 80 (2016).
  * Carbogno, Ramprasad, and Scheffler (2017) C. Carbogno, R. Ramprasad, and M. Scheffler, Phys. Rev. Lett. 175901, 118 (2017).
  * Lampin  _et al._ (2012) E. Lampin, Q. H. Nguyen, P. A. Francioso, and F. Cleri, Appl. Phys. Lett. 100, 131906 (2012).
  * Lampin  _et al._ (2013) E. Lampin, P. L. Palla, P. A. Francioso, and F. Cleri, J. Appl. Phys. 114, 033525 (2013).
  * Puligheddu, Gygi, and Galli (2017) M. Puligheddu, F. Gygi, and G. Galli, Phys. Rev. Mater. 1, 060802 (2017).
  * Lee and Elliott (2017) T. H. Lee and S. R. Elliott, Adv. Mater. 29, 1700814 (2017).
  * Tenenbaum, Ciccotti, and Gallico (1982) A. Tenenbaum, G. Ciccotti, and R. Gallico, Phys. Rev. A 25, 2778 (1982).
  * Toton  _et al._ (2010) D. Toton, C. D. Lorenz, N. Rompotis, N. Martsinovich, and L. Kantorovich, J. Phys.: Condens. Matter 22, 074205 (2010).
  * Vandevondele  _et al._ (2005) J. Vandevondele, M. Krack, F. Mohamed, M. Parrinello, T. Chassaing, and J. Hutter, Comput. Phys. Commun. 167, 103 (2005).
  * Kantorovich and Rompotis (2008) L. Kantorovich and N. Rompotis, Phys. Rev. B 78, 094305 (2008).
  * Konstantinou  _et al._ (2018a) K. Konstantinou, T. H. Lee, F. C. Mocanu, and S. R. Elliott, Proc. Natl. Acad. Sci. USA 115, 5353 (2018a).
  * Konstantinou  _et al._ (2018b) K. Konstantinou, F. C. Mocanu, T. H. Lee, and S. R. Elliott, J. Phys.: Condens. Matter 30, 455401 (2018b).
  * Kapitza (1941) P. L. Kapitza, Physical Review 60, 354 (1941).
  * Singh and Tadmor (2015) A. Singh and E. B. Tadmor, J. Appl. Phys. 117, 185101 (2015).
  * Marks (1997) N. A. Marks, Phys. Rev. B 56, 2441 (1997).
  * Buchan  _et al._ (2015) J. T. Buchan, M. Robinson, H. J. Christie, D. L. Roach, D. K. Ross, and N. A. Marks, J. Appl. Phys. 117, 245901 (2015).
  * Ceriotti, Bussi, and Parrinello (2008) M. Ceriotti, G. Bussi, and M. Parrinello, Phys. Rev. Lett. 102, 020601 (2008).
  * Ceriotti, Bussi, and Parrinello (2010) M. Ceriotti, G. Bussi, and M. Parrinello, J. Chem. Theory Comput. 6, 1170 (2010).
  * Kuwahara  _et al._ (2007) M. Kuwahara, O. Suzuki, Y. Yamakawa, N. Taketoshi, T. Yagi, P. Fons, T. Fukaya, J. Tominaga, and T. Baba, Jap. J. Appl. Phys. 46, 3909 (2007).
  * Zalden  _et al._ (2014) P. Zalden, K. S. Siegert, S. Rols, H. E. Fischer, F. Schlich, T. Hu, and M. Wuttig, Chem. Mater. 26, 2307 (2014).
  * Peng, Cheng, and Mansuripur (1997) C. Peng, L. Cheng, and M. Mansuripur, J. Appl. Phys. 82, 4183 (1997).
  * Giraud  _et al._ (2005) V. Giraud, J. Cluzel, V. Sousa, A. Jacquot, A. Dauscher, B. Lenoir, H. Scherrer, and S. Romer, J. Appl. Phys. 98, 13520 (2005).
  * Lyeo  _et al._ (2006) H. K. Lyeo, D. G. Cahill, B. S. Lee, J. R. Abelson, M. H. Kwon, K. B. Kim, S. G. Bishop, and B. K. Cheong, Appl. Phys. Lett. 89, 151904 (2006).
  * Zhao  _et al._ (2017) Y. Zhao, D. Liu, J. Chen, L. Zhu, A. Belianinov, O. S. Ovchinnikova, R. R. Unocic, M. J. Burch, S. Kim, H. Hao, D. S. Pickard, B. Li, and J. T. Thong, Nat. Comm. 8, 15919 (2017).
  * Bouzid  _et al._ (2017) A. Bouzid, H. Zaoui, P. Luca Palla, G. Ori, M. Boero, C. Massobrio, F. Cleri, and E. Lampin, Phys. Chem. Chem. Phys. 19, 9729 (2017).
  * Zhang  _et al._ (2009) S. N. Zhang, J. He, T. J. Zhu, X. B. Zhao, and T. M. Tritt, J. Non-Cryst. Solids 355, 79 (2009).
  * Martin  _et al._ (2018) E. Martin, P. L. Palla, F. Cleri, A. Bouzid, G. Ori, S. Le Roux, M. Boero, and C. Massobrio, J. Non-Cryst. Solids 498, 190 (2018).
  * Allen and Feldman (1993) P. B. Allen and J. L. Feldman, Phys. Rev. B 48, 12581 (1993).


††thanks: Equal contribution.††thanks: Equal contribution.

# Measuring a localization phase diagram
controlled by the interplay of disorder and driving

Peter Dotti  Yifei Bai  Toshihiko Shimasaki  Anna R. Dardia  David Weld
<EMAIL_ADDRESS>Department of Physics, University of California, Santa Barbara,
California 93106, USA

###### Abstract

The interplay of localizing and delocalizing mechanisms such as disorder,
interactions, and driving is a central topic of modern condensed matter
physics. In this work we experimentally explore the interplay between
quasiperiodic disorder and periodic driving, each of which in isolation is
capable of driving a metal-insulator phase transition. Using a 1D
quasiperiodic cold-atom chain we measure transport across the full phase
diagram varying both drive strength and quasi-disorder strength. We observe
lobes of metallic phases bounded by quantum phase transitions which depend on
both drive and disorder. While these observations are broadly consistent with
expectations from a high-drive-frequency theoretical model, we also observe
clear departures from the predictions of this model, including anomalous
changes in localization behavior at lower drive frequency. We demonstrate
experimentally and theoretically that understanding the full measured phase
diagram requires an extension to commonly-used approximate theories of Floquet
matter.

††preprint: APS/123-QED

Due to the extended character of Bloch wavefunctions, quantum lattice systems
tend to be metallic unless they are subjected to some localizing influence.
Such influences include interactions (Mott localization), disorder (Anderson
localization), quantum statistics (band insulators), and periodic driving
(dynamic localization). The interplay of various localizing and delocalizing
phenomena is a topic which remains at the forefront of the study of condensed
matter, as evidenced for example by the recent intense investigation of the
possible many-body localized state [1, 2, 3, 4].

In this work, we experimentally and theoretically investigate the interplay of
two independent localizing influences: periodic driving and quasiperiodic
pseudo-disorder. A time-periodic driving force applied at specific amplitudes
can cause dynamic localization by resonantly flattening out a Floquet band
[5]. Separately, quasiperiodic pseudo-disorder can drive a localization-
delocalization phase transition via Anderson-insulator-like interference
effects, as expressed in the 1D Aubry-André-Harper (AAH) model [6, 7, 8]. This
work aims to explore open questions concerning the application of _both_ of
these localizing influences to a quantum system: does the driving destroy the
Anderson-like insulator? Can the quasi-disorder disrupt dynamical
localization? What phase diagram arises from the interplay of the two effects?
How is the phase diagram modified as the drive frequency approaches the
bandwidth?

Figure 1: Observing two kinds of localization. (a) Primary (blue) and
secondary (orange) lattices combine to create a quasiperiodic potential
(green) which can be translated via frequency modulation. (b) Atoms imaged
after 500 ms hold time for different secondary lattice strengths demonstrate
AAH localization at $\Delta/J>2$. (c) Atoms imaged after 500 ms hold time for
different drive strengths $K$ demonstrate dynamic localization. Top panel
indicates theoretically predicted $|J_{\mathrm{eff}}|$. The lattice wavevector
is in the vertical direction in (b) and (c).

Previous work has theoretically explored transport in a driven disordered
lattice [9, 10, 11, 12] and experimentally investigated the effects of
disorder in the vicinity of dynamic localization in waveguides [13]. These
studies were generally restricted to the high-frequency limit.

The experiments we describe begin by loading a Bose-Einstein condensate of
$8\times 10^{4}$ 84Sr atoms into a pair of superimposed 1D optical lattices
with different spatial periods aligned along an axis perpendicular to gravity.
The primary (secondary) lattice is formed by two counter-propagating beams of
wavelength $\lambda_{p}=1064$ nm ($\lambda_{s}=874.61$ nm
$\equiv\lambda_{p}/\beta$). All four beams pass through individual acousto-
optic modulators. This allows for control of the lattice depths, denoted
$V_{p}$ and $V_{s}$ for the primary and secondary lattices, and for
synchronized translation of the lattices by setting frequency differences
$\Delta\nu_{p}$ ($\Delta\nu_{s}$) between the two beams forming the primary
(secondary) lattice (Fig. 1(a)). The resulting ability to rapidly tune the
depth and spatial phase of both lattice potentials enables the simultaneous
experimental realization of dynamic and quasidisorder-induced localization.
More experimental details appear in [14].

Synchronized sinusoidal lattice phase modulation is conveniently treated by
transforming to the co-moving frame, in which the modulation appears as a
time-varying electric field. We focus on the case where the secondary lattice
is at rest in the co-moving frame of the primary lattice, implemented by
fixing $\Delta\nu_{s}=\beta\Delta\nu_{p}$. This admits a description by the
tight-binding Hamiltonian [11]

$\begin{split}H=&\sum_{l}-J\left(\ket{l+1}\\!\bra{l}+\mathrm{h.c.}\right)+\Delta\cos(2\pi\beta
l-\delta)\ket{l}\\!\bra{l}\\\ &+\hbar\omega K_{0}\cos(\omega
t)l\ket{l}\\!\bra{l},\end{split}$ (1)

where $J$ is the tunneling energy, $\Delta$ is the quasidisorder strength,
$K_{0}$ is the dimensionless drive strength,
$\beta=\lambda_{p}/\lambda_{s}=1064/874.61$, and the phason parameter $\delta$
is constant. The tight-binding approximation is justified provided that
$V_{p}$ and $V_{p}/V_{s}$ are sufficiently large, as is the case in our
experiments [15, 16]. $V_{p}$ is chosen to be 9$E_{R,p}$, leading to
$J=0.0242E_{R,p}$, where $E_{R,p}=2\pi^{2}\hbar^{2}/m\lambda_{p}^{2}$ is the
recoil energy of the primary lattice. We explore a range of $V_{s}$ values,
from 0 up to 0.592$E_{R,p}$.

Two special cases of Hamiltonian (1) correspond to the two localization
mechanisms previously mentioned. When $K_{0}=0$, the Hamiltonian becomes the
static AAH model. For irrational $\beta$ [17], all eigenstates of this
Hamiltonian exhibit a localization-delocalization quantum phase transition at
$\Delta/J=2$. Fig. 1(b) experimentally demonstrates the phase transition of
the AAH model, showing _in situ_ absorption images taken after 500 ms of
expansion in the bichromatic lattice potential for $K_{0}=0$ with $\Delta/J$
values ranging from 0 to 6.7. As in Ref. [7], the phase transition is observed
in the suppression of expansion for values of $\Delta/J>2$. When $\Delta=0$ in
Hamiltonian (1), dynamic localization can occur [18, 19, 20, 21]. Applying
Floquet theory [22, 23] results in a modified tight-binding Hamiltonian with
effective tunneling energy $J_{\mathrm{eff}}=J\mathcal{J}_{0}(K_{0})$ [18].
Crucially, hopping is entirely suppressed for dimensionless drive amplitudes
corresponding to zeros of the Bessel function $\mathcal{J}_{0}$, for example
at $K_{0}\approx 2.405$ and $5.520$. Fig. 1(c) shows _in situ_ absorption
images taken after 500 ms of expansion for varying drive amplitudes at
$\Delta=0$ and $\omega=2\pi\times(1\mathrm{\ kHz})=19.7J/\hbar$. Expansion is
initiated by sudden removal of a confining optical dipole trap, simultaneous
with the beginning of sinusoidal phase modulation. These data clearly
demonstrate the suppression of transport due to dynamic localization for
$K_{0}$ close to the zeroes of $\mathcal{J}_{0}(K_{0})$.

Figure 2:  Interplay of dynamic localization and quasidisorder-induced
localization in two different drive frequency regimes. (a) Measured width
$\sigma$ (standard deviation of a Gaussian fit) of the atomic density
distribution after a 500 ms hold, for a high drive frequency
$\omega=2\pi\times(1\mathrm{\ kHz})=19.7J/\hbar$. (b) Same quantity for a
lower drive frequency $\omega=2\pi\times(200\mathrm{\ Hz})=3.93J/\hbar$, near
the width of the ground band. The orange dot marks the experimental condition
in Fig. 4. (c) and (d) show numerically calculated average Floquet eigenstate
IPR for the 1 kHz drive and the 200 Hz drive respectively. In all panels
dashed lines indicate phase boundaries predicted by the approximate model
described in the text. In the region $\Delta/J>2$, we observe clear
differences between (a) and (b) and between (c) and (d), indicating differing
Floquet eigenstate properties in the two frequency regimes. $V_{p}=9E_{R,p}$
for all panels. The white arrow in (d) indicates the observed region of
enhanced localization, and corresponds to the arrow in Fig. 3(a).

Having demonstrated that the experimental system can exhibit both
quasidisorder-induced localization and drive-induced localization, we proceed
to investigate the main topic of this work, which is the interplay of these
two phenomena. We measured the width of the atomic density distribution after
500 ms of expansion, with quasidisorder strength $\Delta$ and drive amplitude
$K_{0}$ independently varied to map out a 2D parameter space. Figs. 2(a) and
(b) present experimental measurements of the localization phase diagram for a
drive frequency of 1 kHz ($\hbar\omega=19.7J$) and 200 Hz
($\hbar\omega=3.9J$), respectively. For comparison, Figs. 2(c) and 2(d)
present numerical calculations of the average inverse participation ratio
$\langle\mathrm{IPR}\rangle$ of the Floquet eigenstates.

The basic structure of Fig 2(a) and (c) can be best understood by the time-
averaged Hamitonian:

$\begin{split}H_{\text{eff}}\approx&\sum_{l}-J_{\mathrm{eff}}\left(\ket{l+1}\\!\bra{l}+\mathrm{h.c.}\right)+\Delta\cos(2\pi\beta
l-\delta)\ket{l}\\!\bra{l}\end{split}$ (2)

which holds in the high-frequency approximation $\hbar\omega\gg 4J,2\Delta$.
Eq (2) is the AAH model with a Bessel-renormalized hopping amplitude
$J_{\mathrm{eff}}$. Thus, the localization phase transition occurs at
$\Delta_{c}=2|J_{\mathrm{eff}}|=2|J\mathcal{J}_{0}(K_{0})|$. This
theoretically predicted phase boundary is indicated by dashed curves in Fig.
2.

The first main result of this work is the experimental observation of the
predicted delocalized lobes in the region $|\Delta/J_{\mathrm{eff}}|<2$. The
data clearly exhibit a drive-amplitude-dependent localization-delocalization
phase transition, featuring metallic lobes separated by fingers of an
insulating phase touching the y axis at the zeroes of the Bessel function
$\mathcal{J}_{0}(K_{0})$ where $J_{\mathrm{eff}}=0$. The analytic
$\Delta/|J_{\mathrm{eff}}|=2$ curves correspond very well to the transition
observed both in experimental expansion data (Fig. 2(a) and Fig. 2(b)) and in
numerically calculated $\langle\mathrm{IPR}\rangle$ (Fig. 2(c) and Fig. 2(d)).
The transition between phases is intrinsically somewhat less sharp when
measuring width $\sigma$ rather than $\langle\mathrm{IPR}\rangle$, because the
measured $\sigma$ scales strongly with the effective hopping time
$\hbar/|J_{\mathrm{eff}}|$ even in delocalized regions. The form of the
measured high-frequency phase diagram of Fig. 2(a) agrees with both numerical
theory (Fig. 2(b)) and analytic approximation, indicating that these data
constitute observation of the phase boundaries arising from the interplay
between quasi-disorder and driving.

Figure 3:  Measured ratio between high and low frequency phase diagrams. (a)
Ratio between the expansion with drive frequency
$\omega=2\pi\times(200\mathrm{\ Hz})=3.93J/\hbar$ and
$\omega=2\pi\times(1\mathrm{\ kHz})=19.7J/\hbar$. The data used to generate
(a) is the same as that presented in Fig. 2. (b) Ratio between data taken at
drive frequency $\omega=2\pi\times(300\mathrm{\ Hz})=5.90J/\hbar$ and the same
1 kHz data. The white arrow in (a) is the same as that in Fig. 2(d) and
indicates the observed region of enhanced localization.

For lower drive frequencies, comparable to the bandwidth, new features appear
in both the experimental and numerical phase diagrams which are not captured
by the model of Eq. (2). The observation and characterization of these new
features constitutes the second main result of this work. These differences
are apparent in Figs. 2(b) and 2(d). To visualize the differences, in Fig.
3(a) we show the ratio between the data presented in Fig. 2(b) and Fig. 2(a),
denoted $\sigma_{200\mathrm{Hz}}/\sigma_{1\mathrm{kHz}}$. We observe the
dominant deviation from the high-frequency behavior to occur when
$\Delta/J\gtrsim 2$. We note in particular two features, which appear in both
the experimental data of Fig. 2(b) and 3(a) and the numerical results of Fig.
2(d): a narrow band of anomalously low expansion near the line $K_{0}\approx
3.8$ for larger values of $\Delta/J$, and a broad area of anomalously high
expansion in the lower-right of the parameter space. We interpret these
features as a signature of nontrivial phenomena not captured by the simple
high-frequency approximation.

To further investigate the drive frequency dependence of these features of the
phase diagram, we took similar measurements of the expansion with a drive
frequency of 300 Hz ($\hbar\omega=5.90J$), but otherwise identical
experimental conditions. The ratio between the measured widths at 300 Hz and
the 1 kHz data set, denoted $\sigma_{300\mathrm{Hz}}/\sigma_{1\mathrm{kHz}}$,
appears in Fig. 3(b). We observe the same general features at this slightly
higher drive frequency, but the extent and magnitude of the deviation is
reduced, indicating the approach towards the high-frequency limit.

The first anomalous feature of the low-drive-frequency phase diagram, the
observed enhanced localization near $K_{0}\approx 3.8$, corresponds closely to
a zero of the _first_ Bessel function $\mathcal{J}_{1}(K_{0})$. To better
understand this observation, we first note that a gauge transformation takes
the oscillating linear force into an oscillating phase in the tunneling
strength [18]. The transformed Hamiltonian reads

$\begin{split}H\rightarrow H^{\prime}=&\sum_{l}-J\left(e^{iK_{0}\sin\omega
t}\ket{l+1}\\!\bra{l}+\mathrm{h.c.}\right)\\\ &+\Delta\cos(2\pi\beta
l-\delta)\ket{l}\\!\bra{l}.\end{split}$ (3)

The effect of the oscillating phase can be captured by the Jacobi-Anger
expansion

$e^{iK_{0}\sin\omega
t}=\sum_{m=-\infty}^{\infty}\mathcal{J}_{m}(K_{0})e^{im\omega t}.$

In the high frequency limit, the nonzero-frequency Bessel terms average away
over many cycles. However, when $\hbar\omega\approx 4J$, higher-order terms
contribute, with the $\mathcal{J}_{1}(K_{0})$ term the most relevant term past
leading order. We thus interpret the observed localized stripe near
$K_{0}\approx 3.8$, indistinguishable within experimental error from the first
zero of the $\mathcal{J}_{1}$ Bessel function near 3.8317, as the effect of a
beyond-leading-order localization process not captured in the commonly-used
high-frequency approximate description of dynamical localization.

Figure 4:  Measured expansion versus time in various regimes. (a) Expansion in a purely periodic lattice (diamond markers) is much more rapid than that in the quasiperiodic lattice in the localized phase (circular markers). (b) Same data as in (a) with y axis expanded to show variation among the expansions in quasiperiodic potentials. Different driving frequencies are represented by the color of the marker. Each data point is averaged from 5 experimental realizations. Here $\Delta/J=6.35$ and $K_{0}\approx 2.4046$. Error bars show standard error of the mean. $\omega/2\pi\,(\mathrm{kHz})$ | 0 | 0.2 | 0.3 | 0.4 | 1
---|---|---|---|---|---
$\gamma$ | 0.24(5) | 0.26(5) | 0.26(7) | 0.12(6) | 0.03(21)
Table 1: Transport exponents $\gamma$ extracted by fitting the expansion of
the condensate in a quasiperiodic lattice under various driving frequencies.
Error bars represent 95% confidence interval.

The second anomalous feature, the apparent decrease in localization in the
lower-right quadrant of the low-drive-frequency phase diagram, requires a
separate explanation. This feature could in principle arise from true
delocalization of the Floquet eigenstates or even an extended critical phase.
However, as discussed below, additional measurements of transport dynamics and
numerical calculations of the IPR of the Floquet eigenstates indicate that
this region in fact remains a localized phase in the presence of driving, but
with an increased localization length.

To explore the nature of this region of the phase diagram, we experimentally
measure the evolution of the density distribution of the condensate at various
driving frequencies, for hold times up to 2 seconds ($\approx 639.4\hbar/J$).
In particular, we focus on the case of a strong quasiperiodic potential
($\Delta/J\approx 6.35$) and a driving amplitude of $K_{0}\approx 2.4046$.
Here, in the absence of periodic driving, the non-interacting system is in the
localized phase; without the quasiperiodic potential, the system is near
dynamic localization ($\mathcal{J}_{0}(K_{0})=1.3\times 10^{-4}$). The results
of transport measurements shown in Fig. 4 characterize the interplay between
quasidisorder-induced and dynamic localization in this regime. We extract the
transport exponent $\gamma$ by fitting the measured evolution to the form
$\sigma(t)=\sigma_{0}(1+t/t_{0})^{\gamma}$ for the quasiperiodic cases [24].
Here $\sigma_{0}$ is the initial width, and $t_{0}$ acts as an activation
time. Localization, diffusion, and ballistic expansion correspond to
$\gamma=0$, $0.5$, and $1$, respectively. We do not fit the periodic case
($\Delta=0$) to this form, since the residual axial trapping from the lattice
beam complicates the dynamics at long evolution time. We benchmark the
expansion of the condensate in two cases: a static periodic lattice and a
static quasiperiodic lattice. The evolution in the periodic lattice (diamond
markers in Fig. 4(a)) shows the typical fast expansion. In the quasiperiodic,
localized case (Fig. 4(b)), the expansion exhibits slow sub-diffusion with a
transport exponent $\gamma=0.24(5)$. This sub-diffusion in the localized
regime is caused by the finite mean-field interaction inherent to our
condensate, as found in previous studies [24, 25].

Fig. 4(b) shows the expansion dynamics in the quasiperiodic lattice with
various driving frequencies. When the driving frequency is low, as in the
$\omega=2\pi\times 200\,\mathrm{Hz}$ case, the measured evolution follows that
of the static quasiperiodic system. As the driving frequency increases, the
width evolution remains subdiffusive and becomes slower. Notably, the
expansion is totally arrested within the experimental time in the high-
frequency case ($\omega=2\pi\times 1\,\mathrm{kHz}$). In all the cases
reported, the measured expansion is slower than the undriven quasiperiodic
case, which is in the localized phase.

These experimental results imply that the anomalous region in the lower-right
of the phase diagram is still a localized phase. We corroborate this with
numerical studies of the single-band model (Eq.1). For these numerical
studies, we fix $K_{0}\approx 2.4046$, the same as the experimental condition
in Fig.4, but vary the quasi-disorder strength $\Delta$. We numerically
diagonalize the time-dependent Hamiltonian (Eq.1) via the Floquet formalism,
and calculate the IPR of the eigenstates in the central Floquet-Brillouin
zone. The scaling of the IPR versus system size distinguishes the phases. Our
calculation of the IPR scaling in the Supplementary Information [14] shows
only a localized phase without any mobility edges.

What then distinguishes the anomalous region observed in experiment? We argue
that this region features localization lengths which are both quantitatively
larger and much more quasi-energy dependent than in the high-frequency phase
diagram. Numerically calculated values of IPR in the low-frequency case are
smaller than those in the high-frequency case (Fig. 5b), meaning that the
localization lengths in the low-frequency regime are larger. Past numerical
studies showed that localization lengths increase under low-frequency driving
for a tight-binding lattice with a truly random disorder [26]. Our
calculations thus generalize this result to a quasiperiodic system in the
localized phase. The IPR of the low-frequency Floquet eigenstates also varies
with quasienergy, as shown in Fig. 5(a) and discussed from a different
perspective in [14]. This is in stark contrast to the static AAH model or the
high-frequency approximation of the driven case in which the localization
length of the eigenstates, and thus their IPR, are independent of the
(quasi-)energy [6]. Therefore, although the expansion of the condensate in the
low-frequency regime resembles that of the static quasiperiodic lattice, their
mechanisms are completely different.

Figure 5:  IPR of Floquet eigenstates in the central Floquet-Brillouin zone in
(a) low and (b) high-frequency regimes. Values of IPR of each eigenstate
fluctuate heavily for the low-frequency case in (a), but they are independent
of each eigenstate in the high-frequency case in (b) up to numerical precision
and finite-size effects. The red dashed lines correspond to the quasiperiodic
strength $\Delta/J$ used in the experiments in Fig. 4.

In summary, we have experimentally and theoretically investigated the
interplay of dynamic and quasidisorder-induced localization. In the high-
frequency regime, the measured phase diagram is consistent with both numerical
and time-averaged analytic theory, featuring metallic Bessel lobes separated
from insulating regions by a phase boundary which depends on both drive
strength and quasidisorder strength. For lower frequency driving, comparable
to the natural energy scale of the system, we observe localization behavior
not captured by the simple time-averaged analytic model, including areas of
both enhanced and suppressed transport, and show that it can be understood via
extensions to that simple model. This work touches on some of the rich
structures that can arise when quantum systems are driven at frequencies near
their natural bandwidth, which may prove relevant for future Floquet
engineering applications where optimal system parameters deviate from the
simple high frequency model. A natural direction for future work is the
addition of tunable interactions, enabling experimental probes of the
interplay among three distinct fundamental localizing mechanisms of many-body
quantum matter: interactions, driving, and disorder [27].

###### Acknowledgements.

We acknowledge helpful conversations with Martin Holthaus, and research
support from the Air Force Office of Scientific Research (FA9550-20-1-0240)
and the NSF QLCI program (OMA-2016245). A.R.D. acknowledges support from the
UC Santa Barbara NSF Quantum Foundry funded via the Q-AMASE-i program under
Grant DMR1906325. The optical lattices used herein were developed in work
supported by the U.S. Department of Energy, Office of Science, National
Quantum Information Science Research Centers, Quantum Science Center.

## References

  * Alet and Laflorencie [2018] F. Alet and N. Laflorencie, Comptes Rendus Physique 19, 498 (2018).
  * Moessner and Sondhi [2017] R. Moessner and S. L. Sondhi, Nature Physics 13, 424 (2017).
  * Abanin  _et al._ [2019] D. A. Abanin, E. Altman, I. Bloch, and M. Serbyn, Rev. Mod. Phys. 91, 021001 (2019).
  * Sierant  _et al._ [2024] P. Sierant, M. Lewenstein, A. Scardicchio, L. Vidmar, and J. Zakrzewski, Many-body localization in the age of classical computing (2024), arXiv:2403.07111 [cond-mat.dis-nn] .
  * Holthaus [1992] M. Holthaus, Phys. Rev. Lett. 69, 351 (1992).
  * Aubry and André [1980] S. Aubry and G. André, Ann. Israel Phys. Soc. 3, 18 (1980).
  * Roati  _et al._ [2008] G. Roati, C. D’Errico, L. Fallani, M. Fattori, C. Fort, M. Zaccanti, G. Modugno, M. Modugno, and M. Inguscio, Nature (London) 453, 895 (2008).
  * Lahini  _et al._ [2009] Y. Lahini, R. Pugatch, F. Pozzi, M. Sorel, R. Morandotti, N. Davidson, and Y. Silberberg, Phys. Rev. Lett. 103, 013901 (2009).
  * Holthaus  _et al._ [1995] M. Holthaus, G. H. Ristow, and D. W. Hone, Phys. Rev. Lett. 75, 3914 (1995).
  * Holthaus and Hone [1996] M. Holthaus and D. W. Hone, Philos. Mag. B 74, 105 (1996).
  * Drese and Holthaus [1997] K. Drese and M. Holthaus, Phys. Rev. Lett. 78, 2932 (1997).
  * Arlinghaus  _et al._ [2011] S. Arlinghaus, M. Langemeyer, and M. Holthaus, in  _Dynamical Tunneling: Theory and Experiment_ , edited by S. Keshavamurthy and P. Schlagheck (CRC Press, 2011) 1st ed., Chap. 12.
  * Guzman-Silva  _et al._ [2020] D. Guzman-Silva, M. Heinrich, T. Biesenthal, Y. V. Kartashov, and A. Szameit, Opt. Lett. 45, 415 (2020).
  * [14] See Supplementary Material for experimental and theoretical details, which includes Ref. [28, 29, 30, 27].
  * Li  _et al._ [2017] X. Li, X. Li, and S. Das Sarma, Phys. Rev. B 96, 085119 (2017).
  * Modugno [2009] M. Modugno, New Journal of Physics 11, 033023 (2009).
  * Jitomirskaya [1999] S. Y. Jitomirskaya, Annals of Mathematics 150, 1159 (1999).
  * Eckardt  _et al._ [2005] A. Eckardt, C. Weiss, and M. Holthaus, Phys. Rev. Lett. 95, 260404 (2005).
  * Graham  _et al._ [1992] R. Graham, M. Schlautmann, and P. Zoller, Phys. Rev. A 45, R19 (1992).
  * Madison  _et al._ [1998] K. W. Madison, M. C. Fischer, R. B. Diener, Q. Niu, and M. G. Raizen, Phys. Rev. Lett. 81, 5093 (1998).
  * Lignier  _et al._ [2007] H. Lignier, C. Sias, D. Ciampini, Y. Singh, A. Zenesini, O. Morsch, and E. Arimondo, Phys. Rev. Lett. 99, 220403 (2007).
  * Holthaus [2016] M. Holthaus, J. Phys. B 49, 013001 (2016).
  * Le  _et al._ [2022] C. M. Le, R. Akashi, and S. Tsuneyuki, Phys. Rev. A 105, 052213 (2022).
  * Lucioni  _et al._ [2011] E. Lucioni, B. Deissler, L. Tanzi, G. Roati, M. Zaccanti, M. Modugno, M. Larcher, F. Dalfovo, M. Inguscio, and G. Modugno, Phys. Rev. Lett. 106, 230403 (2011).
  * Shimasaki  _et al._ [2024] T. Shimasaki, M. Prichard, H. E. Kondakci, J. E. Pagett, Y. Bai, P. Dotti, A. Cao, A. R. Dardia, T.-C. Lu, T. Grover, and D. M. Weld, Nature Physics 10.1038/s41567-023-02329-4 (2024).
  * Martinez and Molina [2006] D. F. Martinez and R. A. Molina, Phys. Rev. B 73, 073104 (2006).
  * Bairey  _et al._ [2017] E. Bairey, G. Refael, and N. H. Lindner, Phys. Rev. B 96, 020201 (2017).
  * Walters  _et al._ [2013] R. Walters, G. Cotugno, T. H. Johnson, S. R. Clark, and D. Jaksch, Phys. Rev. A 87, 043613 (2013).
  * Thouless [1972] D. J. Thouless, Journal of Physics C: Solid State Physics 5, 77 (1972).
  * Kramer and MacKinnon [1993] B. Kramer and A. MacKinnon, Reports on Progress in Physics 56, 1469 (1993).