Observation of Stark many-body localization without disorder
W. Morong, F. Liu, P. Becker, K. S. Collins, L. Feng, A. Kyprianidis, G. Pagano, T. You, A. V. Gorshkov, C. Monroe
OObservation of Stark many-body localization without disorder
W. Morong, ∗ F. Liu, P. Becker, K. S. Collins, L. Feng, A.Kyprianidis, G. Pagano, T. You, A. V. Gorshkov, and C. Monroe Joint Quantum Institute and Joint Center for Quantum Information and ComputerScience, University of Maryland and NIST, College Park, MD 20742 USA Department of Physics and Astronomy, Rice University, Houston, TX 77005 USA
Thermalization is a ubiquitous process of statistical physics, in which details of few-body observ-ables are washed out in favor of a featureless steady state. Even in isolated quantum many-bodysystems, limited to reversible dynamics, thermalization typically prevails [1]. However, in thesesystems, there is another possibility: many-body localization (MBL) can result in preservation ofa non-thermal state [2, 3]. While disorder has long been considered an essential ingredient for thisphenomenon, recent theoretical work has suggested that a quantum many-body system with a uni-formly increasing field—but no disorder—can also exhibit MBL [4], resulting in ‘Stark MBL’ [5].Here we realize Stark MBL in a trapped-ion quantum simulator and demonstrate its key properties:halting of thermalization and slow propagation of correlations. Tailoring the interactions betweenionic spins in an effective field gradient, we directly observe their microscopic equilibration for a vari-ety of initial states, and we apply single-site control to measure correlations between separate regionsof the spin chain. Further, by engineering a varying gradient, we create a disorder-free system withcoexisting long-lived thermalized and nonthermal regions. The results demonstrate the unexpectedgenerality of MBL, with implications about the fundamental requirements for thermalization andwith potential uses in engineering long-lived non-equilibrium quantum matter.
Many-body localization was first formulated as a gen-eralization of the non-interacting Anderson transition [6–9]. With disorder, quantum particles can experience de-structive interference through multiple scattering, caus-ing a transition to exponentially localized wavepackets.Over time, a cohesive picture of MBL in interacting sys-tems has also developed [10, 11]. In this description,the MBL regime has extensive local conserved quantitiesthat generalize the particle occupancies in Anderson lo-calization. However, interactions result in additional slowspreading of correlations via entanglement. Strikingly,MBL creates a phase of matter that is non-ergodic: for arange of parameters, local features of the initial state arepreserved for all times, preventing thermalization [2].In considering MBL, it is natural to ask whether ran-dom disorder is a requirement. A partial answer has longbeen known: MBL is possible with incommensurate pe-riodic potentials [12]. However, the question of whetheran MBL phase might exist which preserves translationalsymmetry, for instance in a system with gauge invariance[13] or multiple particle species [14, 15], has continued togenerate extensive discussion [16]. Recently, this prob-lem has been approached from a different starting point:the Bloch oscillations and Wannier-Stark localization ofnon-interacting particles in a uniformly tilted lattice [17].From this, it has been predicted that interacting systemswith a strong linear tilt can also display MBL-like behav-ior [4, 5]. This effect, sometimes called Stark MBL, hasattracted considerable theoretical and experimental in-terest [18–27]. However, clear experimental realization ofStark MBL has been complicated by exact degeneraciesbetween states with the same center of mass that occur ∗ [email protected] in the limit of short-range interactions [4, 5, 26]. The set-ting of a trapped-ion quantum simulator with long-rangeinteractions naturally overcomes this complication. EXPERIMENTAL SETUP
Investigation of many-body localization has beendriven in part by the development of isolated quantumsimulator platforms with site-resolved control and detec-tion [28–31]. Our experimental apparatus (Fig. 1a) exem-plifies these capabilities. It consists of a chain ( N = 15–25) of Yb + ions, with pseudospin states |↑ z (cid:105) and |↓ z (cid:105) encoded in hyperfine ground-state levels. The Hamilto-nian has two ingredients. The first is an overall spin-spininteraction, mediated by global laser beams coupling spinand motion using the Mølmer-Sørensen scheme [32]. Thesecond, a tightly-focused beam creating a programmableeffective B z magnetic field at each ion using the AC Starkeffect [33]. A key feature of this platform is its high de-gree of controllability. In addition to turning on or off ei-ther Hamiltonian term, we use the tightly-focused beamto initialize spins in any desired product state, and wemeasure arbitrary local observables with state-dependentfluorescence collected onto a CCD camera.Combining the global spin-spin interactions with a pro-grammable local field set to a linear gradient results inthe tilted long-range Ising Hamiltonian ( (cid:126) = 1): H = (cid:88) j
0. In practice, we apply the terms inthis Hamiltonian sequentially in time, using a Trotter-ization scheme that reduces decoherence while still re-sulting in evolution closely following the Hamiltonian inEq. 1 (see Methods). The bias field B z is set to be large( B z /J > (cid:80) j (cid:104) σ zj (cid:105) is approximately conserved. With this constraint, andneglecting edge effects, J jj (cid:48) = J | j − j (cid:48) | and this Hamilto-nian is translationally invariant: the operation j → j + n for integer n is equivalent to a shift in B z , which hasno effect in the bulk. For an initial state of definite to-tal magnetization, this model can then be mapped to achain of hard-core bosons with long-range hopping in atilted lattice (see Methods), indicating that it has similaringredients to models previously shown to realize StarkMBL [4, 5]. This system has also been used previouslyto study MBL in a disordered field [28].A useful numeric diagnostic of whether a model ex-hibits an MBL regime can be found in the level statistics,which feature similar behavior in regular MBL [34] andStark MBL [4, 5]. A generic thermalizing ergodic systemhas energy levels following the Wigner-Dyson distribu-tion that characterizes random matrices, while a genericmany-body localized system has a Poissonian level dis-tribution [34]. This difference can be quantified by theaverage ratio of adjacent energy level gaps, defined as (cid:104) r (cid:105) = 1 n (cid:88) n min( E n +1 − E n , E n − E n − )max( E n +1 − E n , E n − E n − ) . (2)The quantity (cid:104) r (cid:105) varies from 0.53 for the Wigner-Dysoncase to 0.39 for the Poissonian case [4, 5, 34]. Diagonaliz-ing the Hamiltonian (Eq. 1) for N = 15, we find that (cid:104) r (cid:105) moves from 0.50 to 0.39 as the gradient g/J is increased,suggesting localization (Fig. 1b). While Fig. 1b showsthe exact experimental Hamiltonian, including deviations from uniform interactions near the edges of the chain,this behavior persists for the ideal power-law Hamilto-nian (see Methods). Unlike the first studies of StarkMBL, where a small amount of disorder or curvature wasrequired to create a state with generic Poissonian levelstatistics [4, 5], Eq. 1 exhibits these properties withoutany terms perturbing the translational symmetry.We probe the degree of localization using a quench pro-cedure, shown schematically in Fig. 1c. The initial state,such as a N´eel state of staggered up and down spins, ishighly excited and far-from-equilibrium. If it thermal-izes, it will result in a high-temperature equilibrium inwhich each spin has nearly equal probabilities of being upor down. Many-body localization will instead result inpersisting memory of the initial configuration, breakingergodicity. ERGODICITY BREAKING IN STARK MBL
Performing the quench experiment, we see the ex-pected signature of localization: a low gradient resultsin quick equilibration of the spins (Fig. 2a), while in astrong gradient they are nearly unchanged from their ini-tial values (Fig. 2b). The experimental data correspondclosely to exact numerics for the system evolution.To quantify the amount of initial state memory as thegradient is increased, it is useful to define a measure thatcan serve as an effective order parameter. We choose ageneralized imbalance, I ( t ), which reflects the preserva-tion of the local magnetizations of the initial state. Thisobservable is similar to other previously used measuresof initial state memory, such as the imbalance [35] or theHamming distance [28], but is advantageous for compar-ing different initial states (see Methods). For an initialstate with M spins that are up, and N − M down, I isequal to the subsequent difference between the average FIG. 2. Ergodicity breaking in Stark MBL. a , Ion-resolved dynamics for an initial N´eel state ( N = 15) at g/J = 0 .
24, and b ,at g/J = 2 .
4, corresponding to the red and blue points on Fig. 1b. While the state quickly relaxes to a uniform magnetizationin the low gradient, the large gradient results in a persisting memory of the initial state. The top row are experimental data,and the bottom row are exact numerics. c , Memory of the initial state, here a N´eel state ( N = 15), can be quantified by thegeneralized imbalance I . For a state of frozen up and down spins, I = 2, and for complete relaxation to a uniform state, I = 0.As the gradient is increased (light to dark), the imbalance crosses from quick relaxation towards zero to a persistent finite value.Points are experimental data at g/J = { } , with statistical error bars smaller than the symbol size, and lines areexact numerics for the lowest and highest gradient. d , For various initial states, shown at top, we see a similar value of thelate-time imbalance at large gradient, suggesting complete localization. From top to bottom, the three initial states correspondto the { triangle, square, round } points. e , Dependence of the late-time imbalance on system size is shown, using an initial N´eelstate with N = 15 (a subset of the data in panel b) and N = 25. The overall increase of late-time imbalance with gradientis robust to the system size increase. The pronounced dip in I near g/J = 1 . magnetizations of the two groups: I ( t ) = (cid:80) Mj (cid:104) σ zj ( t ) (cid:105) M − (cid:80) N − Mj (cid:48) (cid:104) σ zj (cid:48) ( t ) (cid:105) N − M (3)where the sums are respectively over the spins initially upand initially down. In general, |I ( t ) | reaches a maximumvalue of 2 for perfect memory of an initial state with upand down spins, and is zero for a uniform state as atthermal equilibrium.The imbalance shows a clear trend as we increase thegradient (Fig. 2c). At lower gradients, it quickly relaxesto a decaying oscillation about zero, indicating quickthermalization. However, as the gradient is increased, theimbalance instead settles to a progressively higher value.Compared to exact numerics, we observe a damping of the imbalance oscillations, resulting in a lower steady-state imbalance. Furthermore, at the highest gradientvalues, decoherence causes a slow decay of I over time.These are both attributed primarily to intensity fluctua-tions in the tightly-focused beam, as well as to residualcoupling to ion-chain motion from the Mølmer-Sørensenbeams. However, the separation between this decoher-ence time and the fast relaxation dynamics allows us tocharacterize the late-time imbalance.To study initial-state memory for different gradients,we average I ( t ) over a time window tJ from 5 to 7. Thiswindow is chosen to be late enough that transient oscilla-tions have largely decayed, while early enough that deco-herence is limited. This late-time imbalance, ¯ I , capturesthe amount of initial-state memory after fast relaxationhas subsided, and thus the approximate degree of local-ization (Fig. 2d). ¯ I is consistent with zero at the low-est gradient: averaging over the initial states shown inFig. 2d we have ¯ I = 0 . ± . I becomes clearly distinct from zero and progressivelyincreases, reflecting an increasing memory of the initialstate. Crucially, this memory does not show strong de-pendence on the specific initial state chosen: for stateswith different numbers of initial spin flips and differentsymmetry properties, similar behavior is observed. Theinitial state insensitivity observed here is consistent withmany-body localization, which can have some energy de-pendence [20] but is a robust mechanism for breakingergodicity that can span the entire spectrum. This insen-sitivity distinguishes our observations from other effectsthat cause thermalization to have a strong dependenceon the initial state, such as quantum many-body scars[36] and domain wall confinement [37].A key further test of the stability of Stark MBL is tocharacterize the dependence of the observed behavior onincreasing system size. This is especially relevant to lo-calization in systems with long-range interactions, wherefinite-size effects may be particularly important [28, 38].Increasing the spin chain length to N = 25, we see arise in the imbalance at low g/J that is similar to the N = 15 case (Fig. 2e). This length reaches a regime thatis challenging for numerical simulation, and beyond ourability to compute exact dynamics. While we are unableto reach the deeply localized regime for N = 25, due tothe scaling of the experimentally achievable maximumgradient with N (see Methods), the small nonzero valueof I that we observe indicates the persistence of a StarkMBL regime. REVEALING THE CORRELATED STARK MBLSTATE
Probes of the local magnetization, as in Fig. 2, can es-tablish non-ergodicity, but they do not reveal the corre-lations that characterize a localized phase. The structureof the regular MBL phase is understood as being definedby emergent local conserved quantities [10, 11]. Theseconservation laws result in localization, but the localizedregions still have interactions with one another, result-ing in slow spreading of correlations via entanglement af-ter a quench from a product state (typically logarithmicspreading in time, but potentially faster for long-rangeinteractions [39]). While the existence of similar con-served quantities in Stark MBL is debated [23, 24], thereare indications that it can display similar entanglementdynamics [5, 18].Some observables have been established to directlyprobe this correlation spreading, such as quantum Fisherinformation [27, 28] (see Methods and Extended DataFig. 10) or techniques to measure subsystem entangle-ment entropy [30, 31]. Here we instead adopt a localinterferometric scheme, the double electron-electron res-
FIG. 3. DEER Protocol. a , In the spin-echo procedure (darkgreen), a single probe spin undergoes a spin-echo sequence,while the rest of the spins experience normal evolution under H for total time t . In the DEER procedure (dark and lightgreen), there are additional perturbing π/ R spins away. Thedifference in the probe magnetization following these proce-dures reflects the ability of the DEER region to influence thedynamics at the probe spin. We study this protocol using aninitial N´eel state ( N = 15). b , At intermediate times, beforethe spin-echo signal approaches zero due to decoherence, adifference develops between the spin-echo (dark green) andDEER (light green) signals. We quantify this by taking theaverage difference (DEER-spin echo) between tJ = 2 and4 (shaded region) after imbalance dynamics have stabilized.These data are for R = 1 and g/J = 0 . c , As R is in-creased (at g/J = 0 . d , As g is increased (at R = 2), thedifference signal also decreases with increasing gradient, con-sistent with the expectation that within the Stark MBL phase,increasing localization leads to progressively slower develop-ment of correlations. Points in c. and d. are the experimentaldata, and solid lines are exact numerics. onance (DEER) protocol, to reveal the spread of corre-lations controlled by the structure of the localized state[18, 29, 40]. This protocol, shown in Fig. 3a, comparestwo experimental sequences: one that is a standard spin-echo sequence on a probe spin within a system of in-terest, and one that combines this with a set of π/ tJ =2–4, we character-ize the structure of these spreading correlations by tak-ing the average difference between the signals over thistime, ∆ (cid:104) σ z (cid:105) . This time window is slightly earlier thanthe window used for the steady-state imbalance, as theDEER signal is more sensitive to experimental decoher-ence. Varying the DEER spin distance, R , we see thatthis difference signal drops as the DEER spins move pro-gressively farther from the probe, reflecting the local na-ture of correlation propagation (Fig. 3c). Similarly, bysitting at a fixed separation and increasing the gradient,we observe the reduction of the difference signal at agiven time, confirming that the correlation spread is con-trolled by the degree of localization (Fig. 3d). The depen-dences of the difference signal on both R and g/J trackexact numerics, with an overall scaling difference due todecoherence reducing the experimental signal. Taken to-gether, these probes identify the Stark MBL regime asone in which correlations spread slowly through the sys-tem despite persisting memory of the initial state. Thesecorrelations capture the role that interactions play inStark many-body localization, distinguishing it from non-interacting localization. DISORDER-FREE MBL BEYOND A LINEARFIELD
If many-body localized effects are possible in the sim-ple setting of a linearly increasing field, might theyalso appear in a more general class of smoothly vary-ing fields? Utilizing the high degree of tunability ofthis simulator, we investigate a natural generalization:a quadratic, rather than linear, potential. We parame-terize the Hamiltonian as: H = (cid:88) j We have seen the signatures of many-body localiza-tion in a system without disorder, suggesting that theconcept of MBL may be relevant in settings well beyondthe original considerations [8, 9]. Further work could an-alyze these observations in terms of Hilbert space frag-mentation (or shattering) [23, 24, 26, 44], clarifying theconnection between Stark and disordered MBL. Our re-alization of Stark MBL would not appear to naturallyextend to the thermodynamic limit, as this results in in-finite energy differences between different parts of thesystem. However, the Stark MBL Hamiltonian (Eq. 1) isequivalent via a gauge transformation to a Hamiltonianwithout a linear potential that is time-dependent, whichhas a well-defined thermodynamic limit [4].Beyond these conceptual questions, from the perspec-tive of near-term quantum devices, our results suggestthat Stark MBL retains key aspects of the disorderedMBL phase while offering certain advantages, such as notrequiring a fine-grained field and being free of rare-regioneffects or the need for disorder averaging of observables.We summarize some aspects of the comparison in Ta-ble I. Stark MBL may be a useful resource for such de-vices, serving as a tool to stabilize driven non-equilibriumphases [19, 45], or as a means of making a quantum mem-ory [3] with each site spectroscopically resolved. FIG. 4. Relaxation in a quadratic field. a , We reconfigure the site-resolved field from a linear gradient to a quadratic,characterized by the maximum slope γ . For clarity, we show N = 7. b , Dynamics are split into a thermalizing region nearthe center of the system and localized regions near the edges, with the approximate boundaries indicated by the dashed lines.As the maximum gradient is increased, the fraction of the system in the thermalizing regime shrinks. c , Ion-resolved traces ofthe dynamics for max g/J = 1 . 8, showing separation of the spins into localizing regions (bright hues with round points) andthermalizing regions (faded hues with square points). Colors reflect the local field strength at each ion.Disordered MBL Stark MBLErgodicity breaking Yes [2] Yes [4, 5]Slow entanglement growth Yes [2] Yes [5]Max. potential O ( J ) O ( NJ )Requires site-resolved field Yes NoRare-region effects Yes [41, 46] No [4]TABLE I. Comparison of disordered MBL and Stark MBLrequirements, focusing on applications with near-term quan-tum devices. Quasi-periodic MBL occupies an intermediateposition from this perspective, with some of the advantages ofboth disordered and disorder-free localization. For all types ofMBL, questions about the conditions for asymptotic stabilityof localization remain, particularly in long-range interactionsor more than one dimension [4, 41, 46]. ACKNOWLEDGEMENTS This work is supported by the DARPA Driven andNon-equilibrium Quantum Systems (DRINQS) Program(D18AC00033), NSF Practical Fully-Connected Quan-tum Computer Program (PHY-1818914), the DOE Ba-sic Energy Sciences: Materials and Chemical Sci-ences for Quantum Information Science program (DE-SC0019449), the DOE High Energy Physics: Quan-tum Information Science Enabled Discovery Program (DE-0001893), the DoE Quantum System Accelera-tor, the DOE ASCR Quantum Testbed Pathfinder pro-gram (DE-SC0019040), DOE Award DE-SC0019449, theDoE ASCR Accelerated Research in Quantum Comput-ing program (DE-SC0020312), and the AFOSR MURIon Dissipation Engineering in Open Quantum Systems(FA9550-19-1-0399). DATA AVAILABILITY The data that support the findings of this study areavailable from the corresponding author upon request. CODE AVAILABILITY The code used for analyses is available from the corre-sponding author upon request. AUTHOR CONTRIBUTIONS F.L., L.F., and W.M. proposed the experiment. W.M.,P.B., K.S.C., A.K., G.P., T.Y., and C.M. contributed toexperimental design, data collection, and analysis. F.L.and A.V.G. contributed supporting theory and numerics.All authors contributed to the manuscript. COMPETING INTERESTS The authors declare competing financial interests:C.M. is Co-Founder and Chief Scientist at IonQ, Inc. [1] M. Rigol, V. Dunjko, and M. Olshanii, Thermalizationand its mechanism for generic isolated quantum systems,Nature , 854 (2008).[2] D. A. Abanin, E. Altman, I. Bloch, and M. Serbyn, Col-loquium: Many-body localization, thermalization, andentanglement, Reviews of Modern Physics , 021001(2019).[3] R. Nandkishore and D. A. Huse, Many-Body Localiza-tion and Thermalization in Quantum Statistical Mechan-ics, Annual Review of Condensed Matter Physics , 15(2015).[4] E. Van Nieuwenburg, Y. Baum, and G. Refael, FromBloch oscillations to many-body localization in clean in-teracting systems, Proceedings of the National Academyof Sciences of the United States of America , 9269(2019).[5] M. Schulz, C. A. Hooley, R. Moessner, and F. Pollmann,Stark Many-Body Localization, Physical Review Letters , 040606 (2019).[6] P. W. Anderson, Absence of Diffusion in Certain RandomLattices, Physical Review , 1492 (1958).[7] P. A. Lee, Disordered electronic systems, Reviews ofModern Physics , 287 (1985).[8] I. V. Gornyi, A. D. Mirlin, and D. G. Polyakov, Interact-ing Electrons in Disordered Wires: Anderson Localiza-tion and Low-T Transport, Physical Review Letters ,206603 (2005).[9] D. Basko, I. Aleiner, and B. Altshuler, Metal–insulatortransition in a weakly interacting many-electron systemwith localized single-particle states, Annals of Physics , 1126 (2006).[10] M. Serbyn, Z. Papi´c, and D. A. Abanin, Local Conserva-tion Laws and the Structure of the Many-Body LocalizedStates, Physical Review Letters , 127201 (2013).[11] D. A. Huse, R. Nandkishore, and V. Oganesyan, Phe-nomenology of fully many-body-localized systems, Phys-ical Review B , 174202 (2014).[12] S. Iyer, V. Oganesyan, G. Refael, and D. A. Huse, Many-body localization in a quasiperiodic system, Physical Re-view B , 134202 (2013).[13] M. Brenes, M. Dalmonte, M. Heyl, and A. Scardicchio,Many-Body Localization Dynamics from Gauge Invari-ance, Physical Review Letters , 030601 (2018).[14] T. Grover and M. P. A. Fisher, Quantum disentangledliquids, Journal of Statistical Mechanics: Theory and Ex-periment , P10010 (2014).[15] N. Y. Yao, C. R. Laumann, J. I. Cirac, M. D. Lukin,and J. E. Moore, Quasi-Many-Body Localization inTranslation-Invariant Systems, Physical Review Letters , 240601 (2016). [16] F. Alet and N. Laflorencie, Many-body localization:An introduction and selected topics, Comptes RendusPhysique , 498 (2018).[17] G. H. Wannier, Dynamics of Band Electrons in Electricand Magnetic Fields, Reviews of Modern Physics , 645(1962).[18] S. R. Taylor, M. Schulz, F. Pollmann, and R. Moess-ner, Experimental probes of Stark many-body localiza-tion, Physical Review B , 054206 (2020).[19] A. Kshetrimayum, J. Eisert, and D. M. Kennes, Starktime crystals: Symmetry breaking in space and time,Physical Review B , 195116 (2020).[20] L. Zhang, Y. Ke, W. Liu, and C. Lee, Mobility edge ofStark many-body localization, arXiv:2009.08357 (2020).[21] T. Chanda, R. Yao, and J. Zakrzewski, Coexistence oflocalized and extended phases: Many-body localizationin a harmonic trap, Physical Review Research , 032039(2020).[22] D. S. Bhakuni and A. Sharma, Stability of electric fielddriven many-body localization in an interacting long-range hopping model, Physical Review B , 085133(2020).[23] E. V. H. Doggen, I. V. Gornyi, and D. G. Polyakov, Starkmany-body localization: Evidence for Hilbert-space shat-tering, arXiv:2012.13722 (2020).[24] V. Khemani, M. Hermele, and R. Nandkishore, Localiza-tion from Hilbert space shattering: From theory to phys-ical realizations, Physical Review B , 174204 (2020).[25] E. Guardado-Sanchez, A. Morningstar, B. M. Spar, P. T.Brown, D. A. Huse, and W. S. Bakr, Subdiffusion andHeat Transport in a Tilted Two-Dimensional Fermi-Hubbard System, Physical Review X , 011042 (2020).[26] S. Scherg, T. Kohlert, P. Sala, F. Pollmann, B. H. M.,I. Bloch, and M. Aidelsburger, Observing non-ergodicitydue to kinetic constraints in tilted Fermi-Hubbard chains,arXiv:2010.12965 (2020).[27] Q. Guo, C. Cheng, H. Li, S. Xu, P. Zhang, Z. Wang,C. Song, W. Liu, W. Ren, H. Dong, R. Mondaini, andH. Wang, Stark many-body localization on a supercon-ducting quantum processor, arXiv:2011.13895 (2020).[28] J. Smith, A. Lee, P. Richerme, B. Neyenhuis, P. W.Hess, P. Hauke, M. Heyl, D. A. Huse, and C. Mon-roe, Many-body localization in a quantum simulator withprogrammable random disorder, Nature Physics , 907(2016).[29] B. Chiaro, C. Neill, A. Bohrdt, M. Filippone, F. Arute,K. Arya, R. Babbush, D. Bacon, J. Bardin, R. Barends,S. Boixo, D. Buell, B. Burkett, Y. Chen, Z. Chen,R. Collins, A. Dunsworth, E. Farhi, A. Fowler, B. Foxen,C. Gidney, M. Giustina, M. Harrigan, T. Huang, S. Isakov, E. Jeffrey, Z. Jiang, D. Kafri, K. Kechedzhi,J. Kelly, P. Klimov, A. Korotkov, F. Kostritsa, D. Land-huis, E. Lucero, J. McClean, X. Mi, A. Megrant,M. Mohseni, J. Mutus, M. McEwen, O. Naaman, M. Nee-ley, M. Niu, A. Petukhov, C. Quintana, N. Rubin,D. Sank, K. Satzinger, A. Vainsencher, T. White,Z. Yao, P. Yeh, A. Zalcman, V. Smelyanskiy, H. Neven,S. Gopalakrishnan, D. Abanin, M. Knap, J. Martinis, andP. Roushan, Direct measurement of non-local interac-tions in the many-body localized phase, arXiv:1910.06024(2019).[30] A. Lukin, M. Rispoli, R. Schittko, M. E. Tai, A. M. Kauf-man, S. Choi, V. Khemani, J. L´eonard, and M. Greiner,Probing entanglement in a many-body-localized system,Science , 256 (2019).[31] T. Brydges, A. Elben, P. Jurcevic, B. Vermersch,C. Maier, B. P. Lanyon, P. Zoller, R. Blatt, and C. F.Roos, Probing R´enyi entanglement entropy via random-ized measurements, Science , 260 (2019).[32] K. Mølmer and A. Sørensen, Multiparticle Entanglementof Hot Trapped Ions, Physical Review Letters , 1835(1999).[33] A. C. Lee, J. Smith, P. Richerme, B. Neyenhuis, P. W.Hess, J. Zhang, and C. Monroe, Engineering large Starkshifts for control of individual clock state qubits, PhysicalReview A , 042308 (2016).[34] V. Oganesyan and D. A. Huse, Localization of interact-ing fermions at high temperature, Physical Review B ,155111 (2007).[35] M. Schreiber, S. S. Hodgman, P. Bordia, H. P. Luschen,M. H. Fischer, R. Vosk, E. Altman, U. Schneider, andI. Bloch, Observation of many-body localization of inter-acting fermions in a quasi-random optical lattice, Science , 842 (2015).[36] H. Bernien, S. Schwartz, A. Keesling, H. Levine, A. Om-ran, H. Pichler, S. Choi, A. S. Zibrov, M. Endres,M. Greiner, V. Vuleti´c, and M. D. Lukin, Probing many-body dynamics on a 51-atom quantum simulator, Nature , 579 (2017).[37] W. L. Tan, P. Becker, F. Liu, G. Pagano, K. S.Collins, A. De, L. Feng, H. B. Kaplan, A. Kypriani-dis, R. Lundgren, W. Morong, S. Whitsitt, A. V. Gor-shkov, and C. Monroe, Observation of Domain WallConfinement and Dynamics in a Quantum Simulator,arXiv:1912.11117 (2019).[38] Y.-L. Wu and S. Das Sarma, Understanding analog quan-tum simulation dynamics in coupled ion-trap qubits,Physical Review A , 022332 (2016).[39] M. Pino, Entanglement growth in many-body localizedsystems with long-range interactions, Physical Review B , 174204 (2014).[40] M. Serbyn, M. Knap, S. Gopalakrishnan, Z. Papi´c, N. Y.Yao, C. R. Laumann, D. A. Abanin, M. D. Lukin, andE. A. Demler, Interferometric Probes of Many-Body Lo-calization, Physical Review Letters , 147204 (2014).[41] W. De Roeck and F. Huveneers, Stability and instabil-ity towards delocalization in many-body localization sys-tems, Physical Review B , 155129 (2017).[42] J. L´eonard, M. Rispoli, A. Lukin, R. Schittko, S. Kim,J. Kwan, D. Sels, E. Demler, and M. Greiner, Signaturesof bath-induced quantum avalanches in a many-body–localized system, arXiv:2012.15270 (2020).[43] S. S. Kondov, W. R. McGehee, W. Xu, and B. De-Marco, Disorder-Induced Localization in a Strongly Cor- related Atomic Hubbard Gas, Physical Review Letters , 083002 (2015).[44] P. Sala, T. Rakovszky, R. Verresen, M. Knap, andF. Pollmann, Ergodicity Breaking Arising from HilbertSpace Fragmentation in Dipole-Conserving Hamiltoni-ans, Physical Review X , 011047 (2020).[45] D. V. Else, C. Monroe, C. Nayak, and N. Y. Yao, Dis-crete Time Crystals, Annual Review of Condensed Mat-ter Physics , 467 (2020).[46] K. Agarwal, E. Altman, E. Demler, S. Gopalakrishnan,D. A. Huse, and M. Knap, Rare-region effects and dy-namics near the many-body localization transition, An-nalen der Physik , 1600326 (2017).[47] R. Islam, E. E. Edwards, K. Kim, S. Korenblit, C. Noh,H. Carmichael, G. D. Lin, L. M. Duan, C. C. JosephWang, J. K. Freericks, and C. Monroe, Onset of a quan-tum phase transition with a trapped ion quantum simu-lator, Nature Communications , 1 (2011).[48] J. Zhang, P. W. Hess, A. Kyprianidis, P. Becker, A. Lee,J. Smith, G. Pagano, I.-D. Potirniche, A. C. Potter,A. Vishwanath, N. Y. Yao, and C. Monroe, Observationof a discrete time crystal, Nature , 217 (2017).[49] J. Zhang, G. Pagano, P. W. Hess, A. Kyprianidis,P. Becker, H. Kaplan, A. V. Gorshkov, Z. X. Gong, andC. Monroe, Observation of a many-body dynamical phasetransition with a 53-qubit quantum simulator, Nature , 601 (2017).[50] C. Monroe, W. C. Campbell, L. M. M. Duan, Z. X. X.Gong, A. V. Gorshkov, P. Hess, R. Islam, K. Kim,G. Pagano, P. Richerme, C. Senko, and N. Y. Yao, Pro-grammable quantum simulations of spin systems withtrapped ions, arXiv:1912.07845 (2019).[51] D. J. Luitz and Y. B. Lev, The ergodic side of the many-body localization transition, Annalen der Physik ,1600350 (2017).[52] A. Nauts and R. E. Wyatt, New approach to many-statequantum dynamics: The recursive-residue- generationmethod, Physical Review Letters , 2238 (1983).[53] B. P. Lanyon, C. Hempel, D. Nigg, M. Muller, R. Ger-ritsma, F. Zahringer, P. Schindler, J. T. Barreiro,M. Rambach, G. Kirchmair, M. Hennrich, P. Zoller,R. Blatt, and C. F. Roos, Universal Digital QuantumSimulation with Trapped Ions, Science , 57 (2011).[54] P. Ponte, Z. Papi´c, F. Huveneers, and D. A. Abanin,Many-Body Localization in Periodically Driven Systems,Physical Review Letters , 140401 (2015).[55] P. Richerme, Z. X. Gong, A. Lee, C. Senko, J. Smith,M. Foss-Feig, S. Michalakis, A. V. Gorshkov, andC. Monroe, Non-local propagation of correlations inquantum systems with long-range interactions, Nature , 198 (2014).[56] B. Neyenhuis, J. Zhang, P. W. Hess, J. Smith, A. C.Lee, P. Richerme, Z.-X. X. Gong, A. V. Gorshkov,and C. Monroe, Observation of prethermalization inlong-range interacting spin chains, Science Advances ,e1700672 (2017).[57] Y. Y. Atas, E. Bogomolny, O. Giraud, and G. Roux,Distribution of the Ratio of Consecutive Level Spacingsin Random Matrix Ensembles, Physical Review Letters , 084101 (2013).[58] P. Hyllus, W. Laskowski, R. Krischek, C. Schwem-mer, W. Wieczorek, H. Weinfurter, L. Pezz´e, andA. Smerzi, Fisher information and multiparticle entan-glement, Physical Review A , 022321 (2012). METHODSEXPERIMENTAL APPARATUSState preparation and readout Our apparatus has been previously described in [47–50]. We employ a three-layer Paul trap to confine Yb + ions in a harmonic pseudopotential with trapping fre-quencies f x,y = 4.64 MHz and either f z = 0.51 MHz( N = 15) or 0.35 MHz ( N = 25). There is a 1-2% day-to-day variation in these frequencies. Pseudospins areencoded in the two clock ground hyperfine states, with | F = 0 , m F = 0 (cid:105) = |↓ z (cid:105) and | F = 1 , m F = 0 (cid:105) = |↑ z (cid:105) . Wedrive coherent global rotations between these spin statesusing stimulated Raman transitions. Long-range spin-spin interactions are generated via a bichromatic beat-note that couples these states via motional modes alongthe ˆ x direction. This is generated by three Raman beamsfrom a pulsed 355 nm laser driving a symmetric pair oftransitions, with average detunings of µ/ π = 200 kHzfrom the red and blue sideband transitions of the high-est frequency (center-of-mass) transverse motional modealong ˆ x . The resulting distribution of J jj (cid:48) couplings hasa best-fit power law of α = 1 . 28 for N = 15 and α = 1 . N = 25, and a best-fit J / π between 0.25 and 0.33kHz, depending on day-to-day variations in laser power.This value of J , calibrated for a given day, is used toscale energies and times in the main text.Each experimental cycle begins with state initializationvia optical pumping and Doppler and resolved-sidebandcooling, which prepares the spin state |↓ z (cid:105) with fidelity > . 99 and the ground motional state with fidelity > . > . | ↑ z (cid:105) → P / transition collected on a CCDcamera, with typical detection errors of 3%. All mea-surements presented in the main text, except for theDEER measurements, are repeated at each setting 200times for statistics. For the DEER measurements, weinstead average over 2000 repetitions, which are takenalternating between DEER and spin-echo sequences ev-ery 100 measurements so that to a very good approxi-mation both sample any noise profile equally. The datapresented have not been corrected for state preparationand measurement (SPAM) errors. Calibration of Hamiltonian parameters The experimental J jj (cid:48) matrix is determined by mea-surements of motional sideband Rabi frequencies andtrap parameters. Past work has validated this model against direct measurements of the matrix elements [28].We directly measure and calibrate the linear field foreach spin individually. As this calibration process is im-perfect, each spin has a finite amount of deviation fromthe ideal linear gradient and thus there is a finite amountof effective site-by-site disorder in the experimental re-alization, with δ B zj gj ≈ . 02. While a small amount ofdisorder can be crucial in simulations of Stark MBL withshort-ranged interactions, because it breaks the exact de-generacies of that problem [4], in the context of long-range interactions the level statistics are already generic,and this disorder does not have a substantial effect onthe system in numerics. As such, we call our system‘disorder-free’ in the sense that we only have small, tech-nical and well-understood imperfections limiting our re-alization of the ideal disorder-free Hamiltonian. Anyreal quantum simulator can only hope to asymptoticallyapproach a perfectly uniform environment, just as anyquantum simulator can only hope to approximately re-alize MBL because there will always be some residualcoupling to the environment that restores ergodicity atsufficiently long times. NUMERICS Studies of Hamiltonian level statistics with (cid:104) r (cid:105) use ex-act diagonalization of the Hamiltonian. For simulationsof dynamics we solve the Schroedinger equation using theKrylov space technique [51, 52].For all numerics, except those shown in the subse-quent Methods sections ‘Numerical studies of the idealpower-law Hamiltonian’ and ‘Scaling of I with systemsize,’ we use the experimentally determined J jj (cid:48) ma-trix. These couplings show some inhomogeneity acrossthe chain, with the nearest-neighbor hopping varying 7%for N = 15. At large ion-ion separation they also showdeviations from power-law behavior, with the couplingsfalling off faster than the best-fit power law [50]. Thecomparison to power-law numerics shows that each ofthese effects does not strongly alter the dynamics. TROTTERIZED M-S HAMILTONIAN We generate two types of Hamiltonian terms in thiswork. The first is the Mølmer-Sørensen Hamiltonian inthe resolved sideband and Lamb-Dicke limits [50], cre-ated with a pair of detuned bichromatic beatnotes: H ( t ) = (cid:88) j,ν σ + j (cid:20) − i Ω η ν b νj a ν e − iω ν t + a † ν e iω ν t )( e − iδ B t − e − iδ R t ) (cid:3) + h.c. (5)Here j is the ion index and ν is the normal mode index, a ν is the destruction operator of a phonon of motion for agiven normal mode of the ion chain, Ω is the carrier Rabi0 Extended Data Figure 5. Trotterization scheme. Left top: Comparison of the imbalance dynamics for the averaged Hamiltonianof Eq. 12 (solid blue line) with the full Trotter evolution (dashed orange), for the case of an initial N´eel state ( N = 15) andparameters corresponding to the strongest experimental field gradient. Left bottom: difference (averaged - Trotter), showingthe the error over experimental timescales is on the order of one percent. Right: experimental examples (top row) of continuousand Trotterized evolution, both at g/J = 1 . 5, compared to simulations (bottom row) using the (slightly different) parameters ofthe individual experimental realizations. Although the Trotterized evolution lasts nearly twice as much time in absolute units,since the averaged J is roughly half as large, it nonetheless shows a substantial reduction in decoherence and improvement infidelity to the desired Hamiltonian. An initial state with one spin flip is chosen for this comparison, as it makes the effect ofdecoherence due to phonons more pronounced compared with a state near zero net magnetization. rate, η ν is the Lamb-Dicke parameter, b νj is the mode am-plitude for ion j , ω ν is the mode frequency, and δ B ( R ) isthe red(blue) detuning. This term generates spin-motionentanglement, and in the limit η ν Ω (cid:28) | δ R,B − ω ν | themotion can be adiabatically eliminated for an effectivespin-spin interaction.The second Hamiltonian term is the local field gener-ated by the individual addressing beam. This beam onlyaddresses one ion at a time, and is rastered across thechain to create an overall field landscape. A single cycleof this term can be written as: H ( t ) = N (cid:88) j B zj σ zj Θ( t − ( j − t pulse )Θ( jt pulse − t ) , (6)with Θ( t ) as the Heaviside theta and t pulse the time fora pulse of the beam on one ion, which we experimentallyfix at t pulse = 0 . µ s.When these terms are applied simultaneously, in thelimit | δ R,B − ω ν | (cid:29) η ν Ω (cid:29) B zj , the transverse Ising Hamiltonian is approximately realized: H T F IM = (cid:88) j,j (cid:48) J jj (cid:48) σ xj σ xj (cid:48) + (cid:88) j B zj N σ zj . (7)However, the validity of this Hamiltonian is limited tosmall B zj . Therefore, when realizing a linear field gra-dient, B zj = gj , this results in the constraint gN (cid:28) η ν Ω, which prevents the simultaneous attainment of longchains and large linear field gradients. For example, fortypical experimental parameters of N = 15, η Ω = 2 π · J = 2 π · 250 Hz, this would require that g/J (cid:28) . 5. When this is not satisfied, additionalphonon terms are present in the Hamiltonian that re-sult in undesired spin-motion entanglement, or effectivedecoherence of the dynamics when measuring only spin.We can reduce these constraints by applying a Trot-terized Hamiltonian [53]. The evolution under this time-varying Hamiltonian can be analyzed using the Magnusexpansion, to find the dominant contributions to time-averaged dynamics [50]. Within this framework, the un-desired effects arise from the commutator [ H ( t ) , H ( t )].1Intuitively, when these terms are no longer applied simul-taneously the effect of this commutator is reduced.Consider unitary evolution of a single Trotter cycle,using the lowest-order symmetrized sequence: U = e − i (cid:82) ∆ t / H ( t ) dt × e − i (cid:82) ∆ t t / t / H ( t ) dt e − i (cid:82) ∆ t t t t / H ( t ) dt (8)The Hamiltonians governing each part of the unitaryevolution may be approximately replaced by their time-averaged values, simplifying both. For H we have (cid:90) ∆ t / H ( t ) dt = (cid:90) ∆ t / (cid:88) j B zj σ zj Θ( t − ( j − t pulse )Θ( jt pulse − t ) dt = ∆ t N (cid:88) j B zj σ zj , (9)an exact identity since each of the terms in H ( t ) com-mute with one another. For H ( t ) we have (cid:90) ∆ t dt (cid:88) j,ν σ + j (cid:20) − i Ω η ν b νj a ν e − iω ν t + a † ν e iω ν t )( e − iδ B t − e − iδ R t ) (cid:3) + h.c. (10)However, this is just the usual M − S Hamiltonian, andin the limit that | δ R,B − ω ν | t (cid:29) δ R = − δ B this results in the pure σ x σ x interaction. Wheninstead a small rotating frame transformation is appliedwe generate the Ising Hamiltonian with a small overalltransverse field [50]: (cid:90) ∆ t dtH ( t ) ≈ ∆ t (cid:88) j,j (cid:48) J jj (cid:48) σ xj σ xj (cid:48) + B z (cid:88) j σ zj . (11)The combined evolution of the full Trotter cycle isthen, to lowest order, described by the Hamiltonian H = ∆ t ∆ t + ∆ t (cid:88) j,j (cid:48) J jj (cid:48) σ xj σ xj (cid:48) + (cid:88) j σ zj (cid:18) B z + ∆ t ∆ t + ∆ t B zj N (cid:19) + O (∆ t ) . (12)We program B zj to the desired functional form and ab-sorb the factors with ∆ t and ∆ t into re-definitions of J and B zj , leading to Eqs. 1 and 4 of the main text. Theconstant term B z does not depend on these times, be-cause it is created by moving into a rotating frame that isapplied to the entire time evolution. This approximationrequires that | δ R,B − ω ν | ∆ t (cid:29) | δ R,B − ω ν | min = µ = 2 π · 200 kHz and ∆ t ≥ µ s, whose product is 22.6. Addition-ally, ∆ t and ∆ t must not be so long that the Trotterapproximation (Eq. 12) breaks down. However, the lowenergy scale of J and the use of the symmetrized Trot-ter form make this limit less constraining than the limitfor continuous evolution, allowing us to reach g/J = 2 . t and ∆ t , because[ g/J ] avg = (∆ t / ∆ t ) g/J . This capability allows us toscan over a range of gradient values with a single calibra-tion, and it makes any errors on the gradient calibrationcommon to all these scans. In the data presented here,we fix the instantaneous values of g and J and vary ∆ t (see Trotterized Hamiltonian parameters). In addition,we ramp the interactions up and down over 9 µ s witha shaped Tukey profile to reduce adiabatic creation ofphonons [48].This implementation of Trotterized Stark MBL dy-namics would be difficult to extend to more than tens ofspins, as the maximum instantaneous shift required onthe edge ion scales as N , leading to the requirement ofan increasingly fast drive. However, given the unboundednature of a linear gradient, any large- scale simulation ofStark MBL is likely to be challenged by the required fielddifference between the two ends.Throughout this discussion, we have taken the per-spective of a Trotterized quantum simulation of a desiredHamiltonian. We could also understand this experimentin terms of Floquet theory. From this perspective, thisdriven system is described stroboscopically by a FloquetHamiltonian, which to lowest order is the Hamiltonian(12), and the steady-state equilibration that we see rep-resents prethermal evolution under this effective Hamil-tonian that is expected be altered at long times by Flo-quet heating arising from the higher-order terms. Whilethis picture offers a complementary way to understandthese results, and interesting connections to studies ofdriven localization [54], for simplicity we focus on theTrotterized perspective. Trotterized Hamiltonian parameters For imbalance measurements at N = 15, we calibrateto g/J of 2.5 for ∆ t = ∆ t . To scan the gradientstrength, ∆ t is fixed at 18 µ s and ∆ t is varied from18-180 µ s. In addition, there is an extra 9 µ s of effectivedead time per Trotter step associated with the Tukeypulse shaping. We fix B z at 2 π · γ = 2 . t = ∆ t , and vary∆ t from 10-180 µ s, with all other settings kept the sameas in the linear gradient.For N = 25, we instead set g/J to 1.25 for ∆ t = ∆ t .∆ t is fixed at 30 µ s, and ∆ t is varied between 25 and190 µ s, again with an extra 9 µ s of effective dead timeper cycle due to pulse shaping. B z is again fixed at 2 π · g/J of 2.0.∆ t is fixed at 18 µ s and ∆ t is varied from 18-180 µ s,plus an extra 9 µ s of dead time associated with Tukeypulse shaping. We fix B z at values varying for differentdatasets between 2 π · MAPPING TO BOSON MODEL Our experimental Hamiltonian, from Eq. 1 of the maintext, is: H = (cid:88) j 2, resultingin the Hamiltonian: H HC = (cid:88) j A typical ergodic system has a single-particle densitymatrix with support throughout the bulk, and thus hasa high degree of overlap between particles. This resultsin level repulsion in the many-body spectrum, leadingto a Wigner-Dyson energy level distribution characteris-tic of random matrices [34]. A typical localized system,on the other hand, has single particles that are spatiallyconfined, and thus have little overlap, resulting in a Pois-sonian distribution of the many-body spectrum. In Ex-tended Data Fig. 6 we show the full distribution of r , theratio of adjacent energy level spacings, for the experimen-tal Hamiltonian at selected values of g/J . We compareit to the probability density distributions resulting fromPoisson and Wigner-Dyson statistics [5]: P p ( r ) = 2(1 + r ) (Poisson), (16) P W D ( r ) = 27( r + r )4(1 + r + r ) / (Wigner-Dyson), (17)where Eq. 17 is an analytic approximation to the Gaus-sian Orthogonal Ensemble based on the Wigner Surmise[57].3While a small field gradient is needed to break the ap-proximate integrability of the Hamiltonian [56] in the lim-its of g = 0 and B z (cid:29) J , over the range of tilts studiedexperimentally the level statistics cross from being closeto the Wigner-Dyson limit, with an evident dip at low r due to the proliferation of avoided crossings, to veryclose to the Poisson limit at high gradients. This shouldbe contrasted with the case of short-range hopping, inwhich the level statistics may be highly non-generic dueto exact degeneracies associated with center-of-mass con-servation, making concepts of Hilbert space shatteringespecially relevant [4, 5, 24, 26, 44]. However, we cannotexclude the possibility of a ‘quasi-ergodic’ regime at lowgradient [23]. Although the level statistics shown here arefor an experimentally measured Hamiltonian, featuringsmall deviations from a perfectly linear gradient, thesedeviations do not substantially affect the level statistics,as the long-range interactions already lift the degenera-cies. In the next section we show this explicitly, usingthe ideal power-law Hamiltonian to study more generalfeatures of Stark MBL with long-range interactions suchas the scaling behavior. NUMERICAL STUDIES OF THE IDEALPOWER-LAW HAMILTONIAN The experimental system is approximately describedby a Hamiltonian with a power-law hopping. However, asthe exact experimental couplings feature inhomogeneityacross the chain and deviations from power-law scalingfor large ion separations, all numerics shown in the maintext (as well as the previous section) use the exact Hamil-tonian as determined by experimental measurements ofmode structure and detuning. Nonetheless, to study thegeneral behavior of the system it is useful to also look atthe power-law Hamiltonian, which captures the dominantbehavior while being translation-invariant and thereforehaving a more natural scaling with size. We study thisnumerically to characterize the behavior of (cid:104) r (cid:105) with re-spect to α and g/J , and to study the finite-size depen-dence. Dependence of (cid:104) r (cid:105) on α and g/J Extended Data Fig. 7 shows the dependence of thelevel statistics (cid:104) r (cid:105) on the Hamiltonian parameters α and g/J . The primary features of the experimental Hamil-tonian statistics are retained, such as non-generic statis-tics for very low gradient values and a crossover from (cid:104) r (cid:105) ≈ . g/J between 0.1 and 2.0. For α < 1, the ergodic regime progressively increases, as theinteraction energy is superextensive in this regime andthus delocalization is always expected for a sufficientlylarge system. For large α , (cid:104) r (cid:105) generally decreases, whichmay reflect an approach to the exact degeneracies thatare present in the short-range limit. The general features Extended Data Figure 7. Dependence of (cid:104) r (cid:105) on power-lawrange α and g/J (N=13, B z /J = 5). In the experimentspresented in the main text α ≈ . N = { } (light todark), for α = 1 . B z /J = 5. observed are consistent with a recent study of long-rangehopping in a tilt [22] that also found persistence of acrossover in (cid:104) r (cid:105) up to N = 18 and for α > Dependence of (cid:104) r (cid:105) on system size Using the power-law Hamiltonian, we can study thedependence of the level statistics on system size. Ex-tended Data Fig. 8 shows this for N ranging from 9 to15. In general, the curves do not exhibit a simple finite-size scaling. This may be due to the long-range interac-tions, which also cause a system size-dependent shift inthe transition in numerics for the disordered MBL case[38]. We see that the gradient-driven localization persistsup to the largest systems we can diagonalize, coincidingwith the size used for most of the data presented in themain text, with a full study of the scaling left as an in-4teresting subject for future work. GENERALIZED IMBALANCE The generalized imbalance used in the main text isdefined as: I ( t ) = (cid:80) j (cid:104) σ zj ( t ) (cid:105) (1 + (cid:104) σ zj (0) (cid:105) ) (cid:80) j (1 + (cid:104) σ zj (0) (cid:105) ) − (cid:80) j (cid:104) σ zj ( t ) (cid:105) (1 − (cid:104) σ zj (0) (cid:105) ) (cid:80) j (1 − (cid:104) σ zj (0) (cid:105) ) (18)For an initial state that is a product of up and downspins along z , this reduces to a simple form: the averagemagnetization of the spins initialized up minus the aver-age magnetization of the spins initialized down. For aninitial state that is fully polarized this imbalance is un-defined, which may be considered as a drawback to thismeasure, but such a state is already near equilibrium andthus is not useful for quantifying equilibration.This definition is similar to many other variations ofthe imbalance. For an initial N´eel state with an evennumber of spins it is identical up to scaling factors to boththe imbalance and the Hamming distance. However, ingeneral this definition offers a few advantages: • Unlike the imbalance, it is exactly zero for a ther-malized system with an odd number of spins. • It does not require any knowledge of the initial stateto be added in by hand, unlike alternative observ-ables in which the initially flipped spins are tracked. • Unlike the Hamming distance, this generalized im-balance is zero for a thermalized system, and hasunits of magnetization difference (therefore rangingfrom -2 to 2). • Finally, this generalized imbalance is less sensitiveto noise than the Hamming distance. An exampleis useful: consider an initial state of one flippedspin ( (cid:104) σ z (cid:105) = 1), with N = 10, and a backgroundof spin-down ( (cid:104) σ z (cid:105) = − Extended Data Figure 9. Scaling of I with system size. Topleft: As the system increases from N = 9 to N = 15, thelargest change is in a sharpening feature near g/J = 1 thatshifts downward and towards higher gradient. Top right:while we cannot solve for I for N = 25, experimentally wesee a similar dip (reproduced from Fig. 2e of the main text).Bottom left: expanded view of I for N = 9 and N = 15,showing similar localization beyond g/J = 1. While the Hamming distance is always 1 at time zero,this generalized imbalance only starts at 2 for an initialstate in which each spin is in a definite state of σ z . InFig. 2c the experimental imbalances do not start exactlyat 2, reflecting SPAM errors. SCALING OF I WITH SYSTEM SIZE Extended Data Fig. 9 shows a comparison of our datafor I varying system size (Fig. 2e) with numerics. Wecannot present an exact comparison, due to the compu-tational resources needed to solve the Schroedinger equa-tion for a 25-spin system. However, we instead presentdata for N = 9 and N = 15, corresponding to the samescaling factor and the lower of the two experimental sys-tem sizes. To facilitate system size comparison, we usethe ideal power-law Hamiltonian for these numerics.For the most part, I only shows a slight shift withincreasing N . However, there is a sharp feature near g/J = 1 . N = 25. Interpretation of this feature inexperimental data is complicated by decoherence that in-creases both with g/J and with N .The dip feature seen here is initial-state dependent,and may reflect a few-body delocalization process that5 Extended Data Figure 10. Quantum Fisher information. Nor-malized quantum Fisher information for a N´eel state ( N = 15)with g/J = 0 . 24 (white) and g/J = 2 . tJ ≈ 1, the QFI isconsistent with saturation for the low gradient, and with slowentanglement growth for the high gradient, with behavior verysimilar to that previously observed in disordered MBL [28]. is especially favorable for the N´eel state. This illustratesthe challenge of determining the onset of localization infinite size systems, and in quenches from a particular ini-tial state [23]. However, any such few-body process wouldonly occur for g/J < 1, the regime where such reso-nances are possible, and we expect (consistent with thebottom panel of Extended Data Fig. 9) that for g/J > QUANTUM FISHER INFORMATION Quantum Fisher information (QFI) has gained atten-tion as a scalable entanglement witness [28, 58]. For apure state, it is nothing more than the variance of the wit-ness operator O : f Q = 4( (cid:104)O (cid:105) − (cid:104)O(cid:105) ) /N . For f Q > f Q =1 N (cid:88) jj (cid:48) ( − j + j (cid:48) (cid:104) σ zj σ zj (cid:48) (cid:105) − ( (cid:88) j ( − j (cid:104) σ zj (cid:105) ) . (19) The results are shown in Extended Data Fig. 10. We seea significant difference between f Q with weak and strongfield gradients. In a weak gradient, entanglement buildsup rapidly before slowly tapering off. In a strong gradient f Q instead grows slowly, exhibiting similar behavior asexpected for entanglement in an MBL phase and in StarkMBL [5].A few shortcomings limit the value of the QFI. First,it is only easily calculated when assuming a pure state.Second, it can only be interpreted as an entanglementwitness when it exceeds one, challenging in a strongly lo-calized phase. Finally, unlike the DEER protocol it doesnot give spatially resolved information. Still, within itslimits the QFI behavior is consistent with the expecta-tions for an MBL phase. The QFI dynamics also closelyresemble previous observations for disordered MBL [28],consistent with expectations that disorder or strong gra-dients result in similar entanglement spreading. ADDITIONAL DEER DATA Additional data for the DEER protocol difference sig-nal (∆ (cid:104) σ z (cid:105) ) is shown in Extended Data Fig. 11. Lookingat the DEER difference signal, we see that correlationsdevelop more slowly as the DEER region R is moved pro-gressively away from the source. For R = 2, these cor-relations are only visible after the imbalance dynamicshave reached a steady state. This rules out attributionof the correlations to the transient population dynamics,and instead resembles the slow correlation dynamics thatoccur in a disordered MBL system after populations havereached a steady state [10, 11, 30]. CRITICAL SLOPE IN QUADRATIC FIELD Extended Data Fig. 12 presents the dependence of thecritical value of g/J for a quadratic field with differentvalues of the curvature γ. The critical value is determinedby the innermost pair of spins that are both separatedfrom the center spin by more than their mutual errorbars, judged by taking the mean and standard deviationof the average magnetizations for the last five time points.The data are largely consistent in suggesting a criti-cal gradient value on the order of g/J = 0 . 5. However,the strongest curvature is notably different, possibly re-flecting a breakdown of the local gradient approximationfor this case. For curvatures less than this, we concludethat the system seems roughly consistent with a pictureof localization that is determined by the local Stark MBLfield slope at any given spin.6 Extended Data Figure 11. DEER Difference signal for R = { } (light to dark), compared with the imbalance I ( t ) for thesame parameters. Data are offset for clarity but otherwise share the same axes. I is taken from the same dataset as the R = 1 spin-echo data, with the probe spin excluded from the imbalance calculation. After tJ ≈ 2, the imbalance is essentiallyconstant at the low but finite steady-state value corresponding to this gradient strength. However, correlation dynamics arestill progressing- in particular, correlations as measured by the difference signal only begin to develop for R = 2 after thispoint. This is similar to the disordered MBL state, in which slow entanglement dynamics continue after the locally conservedpopulations have reached a steady state [10, 11, 30]. Extended Data Figure 12. Dependence of the critical slope separating thermalizing and non-thermalized regions on the curvature γ . As the quadratic curvature is varied, the division between thermalizing and nonthermal regions is largely consistent with acritical slope near g/J = 0 . 5. However, the strongest curvature of γ = 3 . γ the system was completely delocalized, and thus only the lower bound is meaningful. Error bars (aside from the first twopoints) denote a variation of ±±