Engineering a flux-dependent mobility edge in disordered zigzag chains

There has been great interest in realizing quantum simulators of charged particles in artificial gauge fields. Here, we perform the first quantum simulation explorations of the combination of artificial gauge fields and disorder. Using synthetic lattice techniques based on parametrically-coupled atomic momentum states, we engineer zigzag chains with a tunable homogeneous flux. The breaking of time-reversal symmetry by the applied flux leads to analogs of spin-orbit coupling and spin-momentum locking, which we observe directly through the chiral dynamics of atoms initialized to single lattice sites. We additionally introduce precisely controlled disorder in the site energy landscape, allowing us to explore the interplay of disorder and large effective magnetic fields. The combination of correlated disorder and controlled intra- and inter-row tunneling in this system naturally supports energy-dependent localization, relating to a single-particle mobility edge. We measure the localization properties of the extremal eigenstates of this system, the ground state and the most-excited state, and demonstrate clear evidence for a flux-dependent mobility edge. These measurements constitute the first direct evidence for energy-dependent localization in a lower-dimensional system, as well as the first explorations of the combined influence of artificial gauge fields and engineered disorder. Moreover, we provide direct evidence for interaction shifts of the localization transitions for both low- and high-energy eigenstates in correlated disorder, relating to the presence of a many-body mobility edge. The unique combination of strong interactions, controlled disorder, and tunable artificial gauge fields present in this synthetic lattice system should enable myriad explorations into intriguing correlated transport phenomena.


INTRODUCTION
The idea that the transport of quantum particles in a random environment can be completely arrested due to the interference of multiple transport pathways was first pointed out by Anderson six decades ago [1]. While Anderson considered the localization of electrons in disordered solids, the presence of electron-phonon coupling and electron-electron interactions prohibit direct observation of most single-particle localization phenomena in such systems, even at low carrier density. In contrast, quantum simulation experiments using light [2] or atoms [3] have become an important testbed for disorder physics, since in these systems the issues of lattice phonons and interparticle interactions are either naturally unimportant or can be precisely controlled. For cold atoms, the abilities to tune system dimensionality, applied disorder, atomic interactions, artificial gauge fields, and lattice geometry open up myriad possibilities for exploring novel localization phenomena.
In the absence of interactions, Anderson localization is the generic fate of quantum states in lower-dimensional (d ≤ 2) systems featuring static, random potential energy landscapes and short-ranged tunneling [1,4]. In higher dimensions, the increasing density of states with increasing energy ensures the possibility of delocalization. The exploration of an energy-dependent localization transition, i.e., a mobility edge, has even been undertaken in atomic gases [5,6] in three dimensions through the precision control over disorder and atomic state energies. Cold atom techniques in principle also allow for the exploration of such physics in lower-dimensional systems, where mobility edges can be introduced by correlations in the applied disorder or modified lattice connectivities (e.g., through long-range tunneling).
Despite the exquisite control over cold atom systems and the observations of localization in one dimension (1D) over a decade ago, for both nearly random disorder [7] and correlated pseudodisorder [8], single-particle mobility edges (SPMEs) in lower dimensions have gone unexplored. The reasons for this are somewhat technical -it is quite difficult to modify lattice connectivities, and the varieties of engineered disorder that have been explored in experiment have either been practically random (speckle disorder [5][6][7]9], with short-range correlations due to diffraction) or of a particular form of correlated disorder which, due to a peculiar fine-tuning, does not admit a SPME. In the latter case, the pseudo-disorder that arises in a lattice system due to shifts of the site energies by an added, weaker incommensurate lattice is welldescribed by the Aubry-André model [8,[10][11][12]. While this form of correlated pseudodisorder allows for a localization transition in 1D, the fine-tuning of the cosinedistributed site energies and the cosine nearest-neighbor band dispersion results in an energy-independent metalinsulator transition, and thus the absence of a SPME. By deviating from this fine-tuned condition, either by modifying the band dispersion [13] or by modifying the form of the pseudodisorder [14], one can, in principle, controllably introduce a SPME in such a system.
In this work, we add multi-ranged tunneling pathways to a one-dimensional lattice that features site energy pseudodisorder described by the Aubry-André (AA) model. Specifically, we use our synthetic lattice sys-  The two-dimensional zigzag lattice representation, formed by a rearrangement of the one-dimensional picture of (a). A uniform clockwise flux φ through each triangular plaquette is generated via NNN tunneling phases φ with alternating sign. (c) Atomic dispersion indicating first-(black arrows) and second-order (dashed red arrows) Bragg transitions used to couple NN and NNN lattice sites, respectively, in the momentum-space lattice. The recoil energy is given by tem based on parametrically coupled atomic momentum states to engineer independently controllable nearestneighbor (NN) and next-nearest-neighbor (NNN) tunneling terms ( Fig. 1(a)). The combination of NN and NNN tunneling pathways results in closed tunneling loops that can support a nontrivial flux ( Fig. 1(b)), which we control directly through the complex phase of the various tunneling terms. This system realizes an effective zigzag chain with a tunable magnetic flux. With the combination of controlled pseudodisorder and tunable flux, we perform the first explorations of the interplay of disorder and artificial gauge fields.
We observe direct evidence for a flux-dependent SPME in this system, through measurement of the localization properties of the extremal energy eigenstates. In addition to the SPME that results from multi-ranged hopping, we observe asymmetric (with applied flux) localization behavior of the systems lowest-energy and highest-energy eigenstates, caused by the presence of effectively attractive interparticle interactions in the lattice of momentum states [15]. The influence of interactions is even more strongly evident in the case of the 1D AA with only NN tunneling, where a drastic shift in the localization transition is observed between low-and high-energy eigenstates, corresponding to a mobility edge driven purely by inter-particle interactions.

EXPERIMENTAL METHODS
To experimentally engineer effective zigzag chains, which are equivalent to a lattice model with NN and NNN tunneling terms, we coherently couple an array of discrete atomic momentum states with both first-and second-order Bragg transitions, as depicted in Fig. 1(c). Starting with atoms from a stationary Bose-Einstein condensate (BEC) of ∼10 5 87 Rb atoms, we apply a set of counter-propagating lattice laser beams with wavelength λ = 1064 nm, wavenumber k = 2π/λ, and frequency ω + = c/2πλ, allowing for quantized momentum transfer to the atoms in units of ±2 k. The parametric coupling of states separated in momentum by 2 k, which mimics NN tunneling, is realized by using a pair of acoustooptic modulators to write a controlled spectrum of frequency components onto one of the lattice beams. Starting with atoms at rest, the counter-propagating beams are able to couple the momentum states p n = 2n k as synthetic lattice sites. For example, to create a NN tunneling link between adjacent momentum states p = 0 and p = 2 k, a first-order Bragg resonance (solid black arrows in Fig. 1(c)) is fulfilled by matching the photon energy difference of the two laser fields to the added kinetic energy of an atom moving with momentum p = 2 k. More generally, there exists a unique energy difference between any pair of adjacent states with momenta p n and p n+1 , owing to the quadratic free-particle dispersion. In this way, the multiple frequency tones imprinted onto the one Bragg laser field enable the simultaneous addressing of many Bragg resonances.
In this study, we introduce the novel capability to engineer multi-range tunneling through the simultaneous addressing of first-and second-order Bragg transitions, shown in Fig. 1(c) as solid black and dashed red arrows, respectively. Because each of the spectral tones associated with a given NN or NNN tunneling term is unique, we are able to individually control each of the tunneling links in our synthetic lattice. Specifically, all of the site energies, tunneling amplitudes, and tunneling phases in our synthetic zigzag chains are individually controlled by the strength, phase, and frequency of a corresponding frequency component of the multi-frequency beam. For all of the studies described herein, a total of 21 synthetic lattice sites (momentum states) are coupled through firstand second-order Bragg transitions.
In addition to local parameter control, this system supports site-resolved detection by a simple time-of-flight expansion period where the momentum states separate in space according to their momenta, after which absorption imaging is used to determine the population at each site. A more detailed description of this momentum-space lattice scheme can be found in Refs. [16][17][18][19].

HOMOGENEOUS GAUGE FIELD STUDIES
We first demonstrate our control of a homogeneous synthetic gauge field in the zigzag lattice. We directly impose a synthetic magnetic flux φ on every three-site Chiral dynamics in the zigzag lattice. (a) Band structure for φ/π = ±0.5 considering a two-site unit cell (yellow boxes in lattice cartoon), for tunneling ratio t /t = 0.62. Color represents spin polarization σ , or the overlap of the quasimomentum eigenstate with the top (red, spin up) or bottom (blue, spin down) row of the lattice. Dashed black curves represent the folded band structure for t /t = 0. q should be considered "quasiposition" in our momentum-space lattice and is given in terms of the unit cell lattice spacing d = 4 k. (b) Population imbalance between sites 2 and −2 of the 21-site lattice, measured after 180 µs of dynamics (∼ 1.05 /t) with optical density (OD) images of atomic populations at φ/π = 0, ±0.5 above. Dashed and solid curves represent an ideal simulation of the experiment using Eq. (1) and a full simulation of experimental parameters, respectively. (c,d) Site population dynamics for applied flux (c) φ/π = 0.5 and (d) φ/π = −0.5. Left to right: data, full simulation, and ideal simulation of experiment. Arrows indicate direction of chiral motion. Data for (b-d) were taken with averaged NN tunneling time /t = 176(2) µs and tunneling ratio t /t = 0.622(3). All error bars denote one standard error of the mean. OD images in (b) and extracted site populations in (c,d) are plotted with the color scale in (b). plaquette using engineered tunneling phases. Because the plaquettes alternate pointing up and down, to generate a homogeneous positive flux φ we impose an alternating sign on the NNN tunneling phases, as shown in Fig. 1(a,b). The effective tight-binding Hamiltonian describing the 21-site zigzag lattice is then given bŷ where t (t ) is the NN (NNN) tunneling energy andĉ † n (ĉ n ) is the creation (annihilation) operator at site n. The synthetic gauge field, which can lead to the breaking of time-reversal symmetry, allows us to engineer an analog of spin-momentum locking in the zigzag lattice [20][21][22][23][24][25][26]. We consider the upper and lower rows of the lattice as an effective spin degree of freedom with (pseudo)spins σ = 1 and −1, respectively ( Fig. 2(a)). Under conditions of broken time-reversal symmetry (φ = 0, ±π) we expect to observe chiral trajectories for atoms "polarized" on one row of the lattice. The band structure (shown for the tunneling ratio t /t = 0.62 used in experiment) of the lattice shows this correlation between the sign of the group velocity and the (colored) spin/row degree of freedom [19]. The two bands here reflect the twosite unit cell of the lattice, highlighted in yellow boxes.
To explore this spin-momentum locking in experiment, we initialize atoms on the lower row at site 0 and quench on the tunnel couplings according to Eq. (1). With zero applied flux, the population delocalizes across the lattice symmetrically, as shown in the top middle optical density (OD) image of Fig. 2(b). For positive flux φ/π = +0.5 (right panel), population initially in site 0 moves towards lattice site 2, corresponding to counter-clockwise chiral motion. Under a negative flux φ/π = −0.5 (left panel), population moves in a clockwise fashion to lattice site −2. These observed chiral flows for φ = ±π/2 are clear signatures of spin-momentum locking.
By tuning the applied flux, we map out the entire range of chiral behavior, as shown in Fig. 2(b), bottom. Here we plot the population imbalance P 2 − P −2 between lattice sites 2 and −2, such that a positive (negative) value of imbalance indicates counter-clockwise (clockwise) motion. The data agree qualitatively with an ideal simulation of the experiment using only Eq. (1) (dashed curve), but agree more closely with a full simulation of the system parameters (solid curve) which considers the exact form of atomic coupling to the many laser frequency components, accounting for off-resonant Bragg couplings [19].
We are also able to directly observe the fully siteresolved chiral dynamics of initially localized atomic wave packets, as shown in Fig. 2(c,d). For positive flux, we see that atomic population moves counter-clockwise from site 0 to site 2, and further on to sites 4 and 6, remaining confined to the bottom row. Because the initial state (site 0) does not project entirely onto states with positive group velocity, a portion of the population stays near the center plaquette and oscillates between site 0 and sites ±1. Off-resonant Bragg coupling causes deviations from the ideal simulation (right), but these major qualitative features remain present in both the data (left) and full simulation (middle). For the case of negative applied flux, we observe the opposite chiral behavior, demonstrating that the nature of the spin-momentum locking can be controlled by the applied synthetic flux.

LOCALIZATION STUDIES
Localization phenomena in disordered quantum systems depend intimately on the properties of applied disorder and on the connectivity between regions of similar energy. For random potential disorder in three dimensions, a localization-delocalization transition is assured for states with energies beyond a critical value due to an increasing density of states. For a given disorder strength, a mobility edge, or energy-dependent localization transition, is found in such a system [5,6]. In lower dimensions, for truly random potential disorder, all energy states remain localized in the thermodynamic limit even for arbitrarily small strengths of disorder [4].
Considering instead the influence of correlated pseudodisorder, one finds that the localization physics is strongly modified, with delocalization and mobility edges permitted even in lower dimensions. One form of quasiperiodic pseudodisorder that has been of interest to quantum simulation studies with both light [27] and atoms [8] is that described by the diagonal AA model. Interest in this model has stemmed in part from its intriguing localization phenomenology and connections to the Hofstadter lattice model [10,28,29]. Experimental interest in this form of disorder has also been driven by the relative ease of its realization through the overlap of two incommensurate optical lattices [8].
The AA model of pseudodisorder has interesting properties in the context of SPMEs. The highly correlated disorder allows for the possibility of a metallic, delocalized states in lower dimensions. However, a subtlety arises due to a correspondence between the distribution of pseudodisorder -characterized by quasiperiodic, cosine-distributed site energies -and the cosine dispersion in a NN-coupled 1D lattice. This fine tuning results in a metal-insulator transition that occurs at the same critical disorder value (in units of the tunneling energy) for all energy eigenstates, and thus the absence of a mobility edge. By moving away from this fine-tuned scenario in any number of ways -by introducing longer-range hopping [13], by modifying the pseudodisorder correlations [14], or by adding nonlinear interactions [11,12,[30][31][32]] -a SPME can be introduced into the AA model.
The addition of longer-range tunneling, as in our zigzag lattice, allows for the band dispersion to be modified from its simple cosinusoidal form. For a flux of φ = ±π/2, as shown in Fig. 2(a), increasing the tunneling ratio t /t from zero leads to a deformation of the low-energy band structure from quadratic, to quartic, to forming a doublewell structure [33][34][35], with a symmetric modification of the band energies at high energy. The high ground state degeneracy of the quartic band in this system and of flat bands in similar multi-range hopping models has attracted great interest [36][37][38]. Such systems promise interesting localization properties under disorder [13], and the inherent high single-particle degeneracy allows for the study of emergent physics driven by interactions [36][37][38][39]. For all other flux values (φ = ±π/2) the dispersion of the bands at low and high energies is asymmetric, and this system permits the localization properties of the extremal energy eigenstates to be tuned through modification of the effective mass at low and high energies.
Here, we study the localization properties under the AA model on a 1D lattice and on the multi-range hopping zigzag lattice, observing evidence for an interactioninduced mobility edge as well as the emergence of a fluxdependent SPME.

1D Aubry-André localization transition
We first examine the localization properties of the onedimensional AA model, or the t /t = 0 limit of the zigzag lattice. Figure 3(a) shows this model's pseudodisordered distribution of site energies ε n = ∆ cos (2πβn + ϕ), for an irrational periodicity β = ( √ 5 − 1)/2 and a given value of the phase degree of freedom ϕ. Under this model, all energy eigenstates experience a transition from delocalized metallic states to localized insulating states at the same critical disorder, (∆/t) c = 2, for an infinite system size. To probe the crossover in our finite 21-site system, we initialize various energy eigenstates and explore their localization properties as a function of ∆/t. The experiment begins with population at site 0 (the BEC at rest) with all tunnelings turned off. In this initial limit of infinite disorder (∆/t) i = ∞, all eigenstates are trivially localized to individual sites of the lattice, with a vanishing localization length. We can initialize our atoms in a particular energy eigenstate of the system through choice of ϕ, as the eigenstates and eigenstate energies are solely determined by the site energies in this t = 0 limit. We then slowly ramp the magnitude of the tunneling energy to a final value, and probe the localization properties of the prepared eigenstate as a function of ∆/t. The ramp of t to its final strength t/ = 2π × 1013(9) Hz (corresponding to a tunneling time of /t = 157(1) µs, determined through two-site Rabi oscillations) is linear and performed over 1 ms, slow enough to largely remain within the prepared eigenstate. In each experiment, the disorder strength is fixed to a given value ∆, such that the tunneling ramp (always to the same t value) can be seen as traversing in parameter space from ∆/t = ∞ to the chosen final value (shown as an arrow in Fig. 3(b)). We expect that for final values with ∆/t > (∆/t) c , the population should largely remain localized to the initial site, whereas for ∆/t < (∆/t) c we should see population begin to delocalize across the lattice.
In Fig. 3(b), we plot the measured population outside the central three sites P out , averaged over four realizations of the AA phase ϕ/π = {0.96, 0.64, 1.35, 1.88} corresponding to energy eigenstates {|ψ 0 , |ψ 7 , |ψ 7 , |ψ 18 }, where |ψ 0 is the ground state and |ψ 20 is the highest excited state. As expected, the measured delocalized fraction is almost entirely absent for large disorder, and grows steeply for ∆/t < (∆/t) c . We find excellent agreement between our ϕ-averaged measurements and numerical simulation results based on our experimental ramp (dashed curve, idealized simulations ignoring off-resonant Bragg couplings) in Fig. 3(b), suggesting the observation of a localization crossover that is broadened due to finitesize effects as well as the finite ramp duration. This same behavior can also be seen in the integrated optical density data, shown in the inset, which directly shows the averaged site populations for each final disorder value ∆/t. For large disorder, population remains localized to the initial site, while the metallic regime shows population spreading out to sites n = ±7.
The data for individual energy eigenstates is also shown, both as integrated optical density images in Fig. 3(c) and the P out observable in Fig. 3(d). While all four data runs show localization crossovers, their positions in terms of a critical disorder-to-tunneling ratio (∆/t) c differ according to the state energies. Visually, the ground state |ψ 0 appears to localize for smaller disorders than the intermediate energy eigenstates, with the highly excited state |ψ 18 requiring the largest critical disorder strength for localization. While some of the broadening of the transition observed in Fig. 3(b) can be attributed to effects of finite size and finite ramp durations, to a large degree it is explained by this averaging over unique localization transitions of different energy eigenstates.
The difference in localization properties for different energy eigenstates runs counter to our expectations of an energy-independent transition for the NN-coupled AA model, but can be explained by the presence of nonlinear atomic interactions in our momentum-space lattice [15,19]. In particular, the interactions between indistinguishable bosons in momentum space are effectively attractive and site-local, in the sense that direct interactions are present for collisions between two atoms occupying any pair of momentum modes, while exchange interactions are present only when two identical bosons occupy distinguishable modes [40,41]. Thus, while the momentum-space interactions are physically long-ranged and repulsive, they give rise to an effective local attraction. For atoms initially prepared at the site with lowest energy, attractive interactions can be seen to bring atoms further away from tunneling resonance with other sites (Fig. 3(d), inset). Thus, such a state should remain localized even when the disorder drops below the singleparticle critical value. In contrast, for atoms prepared at the highest energy site, attractive interactions effectively lower the total site energy and bring the atoms closer to tunneling resonance with the unoccupied lower-energy sites of the lattice (Fig. 3(d), inset). Then, by filling the high-energy sites with attractively-interacting bosons, the disorder potential can be effectively smoothed out at high energies by atomic interactions [30].
This behavior for our effectively attractive momentumspace interactions is exactly the opposite of that found for real-space repulsive interactions, the influence of which has previously been studied on ground state localization properties of the AA model [30]. The simulation curves in Fig. 3(d) take into account the effective attractive interactions present in our system at an approximate, mean-field level (also ignoring the inhomogeneous atomic density and neglecting off-site contributions of the effective attraction, which arise due to partial indistinguishability of atoms in different momentum states resulting from superfluid screening [19]). The simulations assume a mean-field interaction based on our condensate's central mean-field energy U 0 / ≈ 2π × 860 Hz (as measured through Bragg spectroscopy), which is of the order of the single-particle tunneling energy t/ = 2π × 1013(9) Hz. To account for the inhomogeneous density distribution, we take a weighted average over homogeneous mean-field energies ranging from 0 to the peak mean-field energy U 0 to get an average mean-field energy of U/ ≈ 2π×500 Hz. We then use this average value as a homogeneous mean field energy in our simulations. These simplified simulation curves already reproduce well the observed shifts of the localization transitions for the low-(|ψ 0 ) and high-energy (|ψ 18 ) states. These direct observations of interaction-induced localization and delocalization for low and high-energy states, respectively, are indicative of a many-body mobility edge. Such measurements are enabled by our unique ability to stably prepare any particular eigenstate in our synthetic lattice.

Localization studies in zigzag chains
With the addition of longer-range tunneling, the energy-independent transition of the simple 1D AA model begins to depend critically on the eigenstate energy even at the single-particle level. By tuning the NNN tunneling strength and the artificial flux in our effective zigzag chains, we can introduce a tunable SPME through band structure engineering. While in the demonstration of control over flux and the observation of spinmomentum-locking in Fig. 2 we employed a tunneling ratio of t /t ≈ 0.6, here we work at a smaller value of t /t ≈ 1/4. Under this condition, a maximal difference in the band dispersion at low and high energies appears for flux values of 0 and π, where a quartic dispersion appears at high and low energies, respectively.
To probe the mobility edge, we prepare the two extremal energy eigenstates of the system, the ground state (GS) and the highest excited state (ES), and compare their localization properties. As in the 1D study, our experiment begins with all atomic population prepared at site 0 with all tunnelings turned off, i.e., in the infinitedisorder limit of the system (∆/t = ∆/t = ∞) where  all energy eigenstates are localized to individual sites of the lattice. To initialize the atoms in a particular energy eigenstate of the system, we simply vary the AA phase: ϕ = 0 for the GS and ϕ = π for the ES. In short, we track how the prepared eigenstate evolves as the parameters of the Hamiltonian, given bŷ are smoothly and slowly varied to some final desired conditions of ∆/t for fixed tunneling ratio t /t and fixed flux φ. To help ensure adiabaticity over a large part of the parameter ramp, an extra potential offset of strength V is added at the initial site n = 0, such that the modified site energies are given by ε n (V ) = ∆ cos (2πβn + ϕ) − V δ n,0 . By setting V > 0 (V < 0) for the GS (ES), we further separate the initial eigenstate from the rest of the spectrum by a potential well (hill). Starting from the initial limit of V /t = ∞ and ∆/t = ∞, we adiabatically load our desired eigenstate by linearly ramping up both tunneling terms (t and t ) over 2 ms while also smoothly removing the potential well by ramping V to zero [19].
We perform this procedure over parameter ranges 1 ≤ ∆/t ≤ 4.25 and 0 ≤ φ/π ≤ 1, mapping out the localization behavior of the GS and the ES in Fig. 4(a,d). We plot the standard deviation of the population distribution in the lattice, σ n (i.e., the momentum standard deviation σ p normalized to the spacing between sites of 2 k), where the values are resampled from the actual (∆/t, φ/π) points where data were taken (small black dots). The ∆/t values of the data have variations and uncertainties stemming from variations and measured uncertainties in calibrated tunneling rates for the experimental runs, with an overall averaged NN tunneling rate t/ = 493(2) Hz and tunneling ratio t /t = 0.247 (4).
For the ground state in Fig. 4(a), we see that the region of metallic, delocalized states (red region, corresponding to states with large σ n ) extends out to larger ∆/t values when the applied flux is near zero than for the case of an applied π flux. This can also be seen in the integrated optical density images at bottom: sites as far as n = ±2 remain populated even at large disorder ∆/t ∼ 3.5 at small flux φ/π = 0.05 (left), while for large flux φ/π = 0.95 (right) population fully localizes for ∆/t > 3. The top panel of Fig. 4(b) highlights that for a fixed disorderto-tunneling ratio of ∆/t ∼ 2.9, the GS can be driven from metallic to insulating by changing only the flux.
In the absence of flux, the shift of the GS localization transition to larger disorder values as compared to the t = 0 case is intuitive: simply adding longer-range tunneling increases the connectivity of the lattice, increasing the dispersion at low energy, and enhancing delocalization. As non-zero flux is added, however, the GS localization transition shifts towards smaller critical disorder values. This effect is perhaps surprising when considering effects such as the suppression of weak localization by broken time-reversal symmetry, as observed recently in measurements of coherent back scattering [42]. However, in the context of our zigzag flux chains, this fluxenhanced localization of the GS is easy to interpret. The shift of the GS localization transition towards smaller (∆/t) c is driven by a flattening of the low-energy band dispersion, owing to kinetic frustration of the different tunneling pathways. The system is maximally frustrated at φ = π for t /t = 1/4, corresponding to a nearly flat, quartic low-energy dispersion (Fig. 4(c), right). Under these conditions, the states at low energy become heavy (large effective mass) and easier to localize in the presence of disorder.
In considering the flux-dependent localization properties of the highest energy eigenstate, a similar line of argumentation holds, but with the opposite trend with applied flux. The high energy states of the band structure are maximally dispersive for φ = π, becoming flatter for decreasing flux, with a quartic band appearing for zero flux. The consequence of this modified band structure on the localization properties of the ES is reflected in the measured dependence of the ES localization proper-  Fig. 4(a,d). Non-interacting and interacting simulations (U/ = 2π × 500 Hz used in the latter) are shown as dashed and dotted lines, respectively. For flux values where no critical disorder is plotted, atomic population was determined to be delocalized (based on the set threshold value of the standard deviation) over the full range of disorder strengths. Vertical gray line at φ/π = 0.5 denotes flux value at which the GS and ES curves should cross in the absence of interactions and any off-resonant coupling terms. Error bars denote one standard error of the mean.
ties following the parameter ramp to final ∆/t values for different flux values (Fig. 4(d)). The flux-dependence of the localization transition is also seen in striking fashion in the integrated OD images at the bottom of Fig. 4 For both states, we empirically estimate the approximate "critical" disorder strength (normalized to t) relating to the metal-insulator transition by finding the ∆/t value at which σ n equals 0.68 lattice sites. This estimate is determined for each flux value of the data, and the extracted critical disorder strengths are shown as white circles in Fig. 4(a,d). We can compare these experimentallyextracted points to the predicted threshold values of disorder, based on numerical simulations of our experimental ramp protocol. These single-particle predictions are shown as dashed lines in Fig. 4(a,d), and show the same qualitative trend as the experimental points for both the GS and ES.
To better contrast the localization behavior of the GS and ES, we additionally plot both the experimentally determined transition points and the theory predictions for both extremal eigenstates together in Fig. 5. With the two datasets overlaid, one can more clearly see the direct evidence for a flux-dependent SPME. While this sampling of the two extremal eigenstates does not deter-mine the critical energy at which delocalization occurs for given values of ∆/t and φ, it does provide the first direct experimental evidence for a SPME in lower dimensions.
The behavior of the transition ∆/t values for the GS and ES are nearly opposite to one another. For flux values near zero, the disorder strength needed to localize the GS exceeds that of the ES by nearly t, due to kinetic frustration of the high energy states. The situation reverses for flux values near π: the GS becomes localized at lower disorder strengths ∆/t ∼ 2.3, and the ES remains delocalized even up to the highest disorder value used in experiment ∆/t ∼ 4.25. This apparent asymmetry, i.e., that a larger magnitude of shift between the GS and ES transition points is found for flux values near π than for flux values near 0, is in disagreement with the single-particle prediction. Moreover, at the singleparticle level the flux-dependence of the GS and ES localization properties should essentially be mirror images of one another (dashed lines, with a slight asymmetry resulting from effects due to off-resonant driving), such that their transitions points should cross very near to φ/π = 0.5 (vertical gray line in Fig. 5). However, the apparent crossing point is offset to lower flux values by nearly 0.1π. As in the previously discussed case of the 1D AA model with only NN interactions (Fig. 3(c,d)), the nonlinear interactions present in our atomic system are largely responsible for this asymmetry observed between the localization properties of our low and high energy eigenstates.
As described earlier in the context of the NN-coupled AA model, we can approximately capture the influence of the momentum-space interactions in this system by including a site-local mean-field attraction in a multisite nonlinear Schrodinger equation [19], with an interaction energy that is determined independently by calibration via Bragg spectroscopy. Including these interactions (dotted lines, also shown in Fig. 4(a,d)), the transition lines get shifted to lower (GS) and higher (ES) disorder values, so that they cross at lower flux values. The interacting simulation results better capture the localization properties of the ES, which was shifted to significantly higher disorder strengths than was predicted at the single-particle level. It also qualitatively captures the shift of the crossing of the critical disorder curves in Fig. 5 to lower flux values, although it predicts a slightly larger shift than seen in experiment. In the future, by studying fluctuations of the atomic number distribution and inter-site correlations in our synthetic lattice, or by more closely studying fine features of the localization properties, this simulation platform may enable unique explorations into the physics of interacting disordered systems, in particular related to the physics of many-body localization. It also offers a unique platform to study the interplay of disorder, artificial gauge fields, and interactions.

CONCLUSIONS
This work represents the first direct observation of a single-particle mobility edge in lower dimensions, which is enabled by the unique ability to stably prepare atoms in any energy eigenstate and explore their localization properties in a system with precisely controlled disorder and tunable artificial gauge fields. We also present the first direct quantum simulation evidence for a many-body mobility edge, studied through a shift of the localization properties of low-and high-energy eigenstates in the 1D AA model that arise due to many-body interactions. These interaction shifts are also observed in the localization transitions of a multi-range hopping AA model that admits a flux-dependent SPME, leading to the interplay of single-particle and many-body shifts of the localization transition for states at different energies.
This work also constitutes the first quantum simulation study combining synthetic gauge fields and disorder, and its extension to fully two-dimensional lattices beyond coupled chains promises to pave the way towards studies of disordered quantum Hall systems. In particular, by moving to a larger system containing bulk lattice sites, a robustness of the observed chiral-propagating modes to disorder (similar to the robustness to disorder observed recently for the bulk winding of chiral symmetric wires [18]) should be readily observable.

EXPERIMENTAL SETUP
Our experiments begin with a Bose-Einstein condensate (BEC) containing ∼10 5 87 Rb atoms, held in an optical dipole trap primarily formed by a single focused laser beam (wavelength λ = 1064 nm, wavenumber k = 2π/λ) with weak additional trapping provided by one other beam. To create the lattice, we allow this primary beam to pass through two acousto-optic modulators (AOMs) which, together, write onto the beam a spectrum of radiofrequency tones in the sub-MHz range. This multifrequency beam is sent back towards the atoms along the same path as the incoming single-frequency beam. The interference of these two beams creates a time-dependent optical potential comprised of multiple superimposed optical lattices moving at different velocities. Each of these velocities is determined by the frequency difference between one frequency component of the multi-frequency beam and the incoming beam, specifically tuned to address a particular Bragg transition between two discrete momentum states. The simultaneous driving of many Bragg transitions mimics tunneling between sites (momentum states) in our synthetic lattice [1,2].
By varying the amplitude and phase of each frequency component, along with the frequency detuning from Bragg resonance, we can control the amplitude and phase of each effective tunneling link, as well as each site-energy term. Because the single-particle dispersion relation is quadratic (E = p 2 /2M Rb ), all of the first-and secondorder Bragg transitions have unique frequencies and can be individually addressed.
Next-nearest-neighbor (NNN) tunnelings are realized via a four-photon, second-order Bragg process. Shown in Fig. S1(a) as dashed red lines, this involves virtually absorbing two photons from the single frequency beam (ω + ) and emitting two photons into the multi-frequency beam (ω 0,2 ), or vice versa. Control of the effective lattice parameters for these four-photon processes requires slightly different considerations compared to the nearestneighbor (NN) terms. For example, to enable a NNN tunneling phase of π, we apply half this phase (π/2, relative to incoming field) to the corresponding frequency component ω 0,2 . A similar consideration holds for the relationship between site energy and frequency detuning from resonance. More generally, these differences can be summarized in the relationships between the tunneling terms for NN (t nn e iϕnn ) and NNN (t nnn e iϕnnn ) processes. Taking into account the field strengths (assumed to be real)  (6) and applied flux φ/π = −0.5 (same data as in Fig. 2(d)). Dashed and solid curves in (c,d) represent results from an ideal simulation of the experiment and a full simulation accounting for off-resonant coupling, respectively. of the incoming beam (Ω I ) and a particular frequency component of the multi-frequency beam (Ω R ), the phases of these same fields (φ I and φ R ), the large single-photon detuning from atomic resonance ∆ (relating to roughly 100 THz for our laser wavelength λ = 1064 nm), and the recoil energy E R = h 2 /2M Rb λ 2 , these terms are given at resonance as: (S1) Here we present the exact NN tunneling times /t and NNN to NN tunneling ratios t /t for each individual data set shown in the main text. For the variable flux data of Fig. 2(b), /t = 172(1) µs and t /t = 0.633 (3). For the dynamics under φ/π = 0.5 of Fig. 2(c), /t = 182(1) µs and t /t = 0.605 (5). For the dynamics under φ/π = −0.5 of Fig. 2(d), /t = 174(2) µs and t /t = 0.628 (6). The tunneling ratio used to make the band structures of Fig. 2(a) was an average of these three values: t /t = 0.622 (3). For the 1D Aubry-André data of Fig. 3, /t = 157(1) µs. For the multi-range hopping Aubry-André data of Fig. 4 and Fig. 5, /t = 323(1) µs and t /t = 0.247(4), averaged over all data points taken.

OFF-RESONANT EXCITATIONS
While we seek to address individual Bragg transitions with single frequency components, in practice we apply the full spectrum of frequencies to the condensate, as shown in Fig. S1(b). Thus each transition is not only addressed by one resonant frequency, but also feels the effects of all of the other non-resonant frequency components. For a lattice with only NN tunnelings, the frequency components are equally spaced at 8E R / = 8 × k 2 /2M Rb ≈ 2π × 16.2 kHz. Adding NNN links halves the spacing of applied frequency components to 4E R / ≈ 2π × 8.1 kHz.
On a lattice with only NN tunnelings, the off-resonant couplings result in step-like intervals in the dynamics (Fig. S1(c)). For this data (from Ref. [3]), atomic population undergoes a continuous-time quantum walk on a 21-site lattice engineered with 20 equally-spaced frequency teeth. Due to the equal spacing of the frequency teeth, the off-resonant effects add up constructively. This is evident in the period T between these steps, which corresponds exactly to the spacing between adjacent frequency teeth, T = h/8E R . We have shown in a previous study [3] that adding random tunneling phases onto equally-spaced frequency teeth suppresses these steps, resulting in smoother dynamics.
We note that the magnitude of the tunneling rate plays a significant role in the magnitude of off-resonant effects. If the tunneling rate is comparable to or greater than the frequency spacing between the first-order Bragg resonances (t/ 8E R / ), then each frequency tooth may address multiple transitions. The extreme limit of non-resonant addressing is often encountered in cold atom experiments, e.g., when a deep stationary potential (two interfering fields with equal frequency) that is suddenly turned on results in Kapitza-Dirac diffraction [4] of atomic matter waves. With respect to synthetic lattices, this can be viewed as a system with constant NN tunneling terms in the presence of a quadratic poten-tial, set simply by the single-particle dispersion relation. This leads to expansion dynamics in momentum space at short times up until population reaches outer regions where the site (momentum state) energy roughly equals the effective tunneling bandwidth [5].
On the other hand, in the limit where the tunneling energy is much smaller than the energy spacing between relevant Bragg resonances (i.e., in the limit that the rotating wave approximation is valid), each transition will be ideally addressed by only a single spectral component, and the step-like behavior in Fig. S1(c) (data, solid curve) should approach the ideal smooth behavior (dashed curve). Intuitively, this occurs as the number of steps per tunneling period gets very large, or in other words as the tunneling time /t gets much larger than the step period T , such that the dynamics are spread out over many, many steps. For the data in Fig. S1(c), we are in the intermediate regime (with a tunneling time /t = 111.6(7) µs, corresponding to a tunneling rate of t/ = 2π × 1425(9) Hz), where off-resonant coupling primarily results in the observed dynamics with small steplike behavior.
For the zigzag lattice, we introduce NNN frequency teeth that halve the spacing to 4E R / ≈ 2π × 8.1 kHz (addition of dashed red peaks in Fig. S1(b)). This results in longer steps of exactly twice the duration, i.e., T = 2T (Fig. S1(d)) due to off-resonant first-order Bragg processes. The smaller structure on top of these long steps, with a spacing of T /2 relating to a frequency spacing of 16E R / , is due to off-resonant second-order Bragg processes. For this data, we used a tunneling ratio t /t = 0.628(6) and NN tunneling time /t = 174(2) µs, corresponding to a tunneling rate t/ = 2π × 917(9) Hz.
The "full" simulations in Fig. S1(c,d) (solid curves) and in the main text account for both resonant and offresonant driving on every Bragg transition, and thus retain these step-like features. The "ideal" simulations in Fig. S1(c,d) (dashed curves) and in the main text ignore off-resonant effects and consider only the smooth behavior of the idealized tight-binding Hamiltonians.

BAND STRUCTURE CALCULATIONS
The band diagrams shown in Fig. 2(a) and Fig. 4(c) were calculated using the same method as described in Sec. II of Ref. [6]. While the studies there focused on a semisynthetic lattice with one synthetic dimension and one real-space dimension, the same physics hold for our fully synthetic momentum-space lattice.
We consider the zigzag lattice in the absence of any applied disorder. We can take the rows of the lattice to be an effective spin degree of freedom, with a twosite unit cell comprised of one spin up (σ = +1) site from the top row and one spin down (σ = −1) site from the bottom row. We generate the spinful dispersion by calculating the 2×2 Hamiltonian introduced in Ref. [6] at each value of quasimomentum. 1 In terms of the creation and annihilation operators at spin σ and quasimomentum q (ĉ † σ, q andĉ σ, q ), the Hamiltonian is given bŷ where h jj = −2t cos qd + φ(−1) j and h 12 = h 21 = 2t cos [ qd/2] for lattice spacing d = 4 k (not 2 k due to the two-site unit cell). By diagonalizing this Hamiltonian for every value of q ∈ [−π/d, π/d], we generate the double-band dispersions. We note that for t /t = 0 (dashed black curves in Fig. 2(a) and Fig. 4(c)), the dispersion relation is cosinusoidal, but folded back at the edges of the Brillouin zone due to the two-site unit cell. The spin magnetization σ is simply the projection of the quasimomentum eigenvectors derived from Eq. (S3) onto the rows of the lattice. We take the difference between the projections onto the upper row and the lower row such that a positive (negative) σ corresponds to population on the upper (lower) row.

INFLUENCE OF INTERACTIONS
As mentioned in the main text, atomic interactions show effects on the localization properties of both the 1D Aubry-André data in Fig. 3 and the longer-range Aubry-André data in Fig. 4 and Fig. 5. Interactions in momentum-space lattices are described in detail in Ref. [7], where we show that effects like self-trapping can be observed when the mean-field energy becomes large compared to the tunneling.
Atoms in a particular momentum state experience an added positive self energy due to repulsive cold collisions, i.e., mean-field interactions. In addition, atoms overlapped in space but occupying distinct spatial eigenstates (i.e., distinguishable plane-wave momentum states) experience both a direct interaction as well as an added exchange energy, resulting in twice as large a repulsive energy [8,9]. For a fixed total density, this situation where atoms occupying the same momentum state have a weaker repulsive interaction energy may be recast as an effectively site-local attraction with a scale set by the mean-field interaction energy U . In reality, for an interacting degenerate Bose gas, superfluid screening can make distinct plane wave states partially indistinguishable, resulting in some off-site contribution to the effective attraction (although this vanishes as U becomes much less than 2E R ).
As discussed in Ref. [7], these interactions shift the Bragg resonance frequencies away from the single particle resonances. Under typical experimental conditions, we measured this shift to be ∼300 Hz, relating to a peak mean-field energy U 0 = gn 0 ≈ 2π × 860 Hz at the center of the harmonic trap, and an homogeneous mean-field energy U/ ≈ 2π × 500 Hz averaged over the entire trap (the measured shift is distinct from the average U value due to a combination of the aforementioned screening effects and the long duration of the Bragg pulses used in this determination). The central atomic density is n 0 ≈ 10 14 cm −3 and g = 4π 2 a/M Rb , where a is the scattering length [10].
We incorporate this mean-field energy U by considering an attractive interaction that depends on the population of atoms at each site [7], resulting in the curves shown in Fig. 3(d) and the dotted theory curves in Fig. 4. Specifically, the evolved state at the end of the tunnelingramps described in the main text are found by solving a time-dependent multi-site nonlinear Schrödinger equation that includes a local attractive self-nonlinearity −U |ψ n | 2 , where the ψ n are c-numbers (with normalization Σ n |ψ n | 2 = 1) relating to the atomic field terms at each site. This approach is approximate in a number of ways: it ignores quantum fluctuations, ignores off-site contributions to the effective atomic interaction, ignores energy-dependent corrections to the collisional scattering cross-section, ignores spatial variations in the atomic density in our trapped sample, ignores effects such as the loss of spatial overlap of the momentum wavepackets, and explicitly restricts the collisions to be mode-preserving (ignoring both s-wave collisions that may scatter atoms out from the considered set of 21 modes into many "halos" of additional states, as well as mode-changing collisions within our defined set of states, that would be energetically suppressed in the absence of our drive fields, but may be effectively enabled through higher-order, Braggmediated processes). We point out that a fuller quantum treatment of the problem (still restricted to being modepreserving, still ignoring spatial variations of the density, still ignoring loss of spatial overlap of momentum states, and still ignoring energy-dependent corrections to the scattering cross-section) would instead include an interaction term (U/N )Σ n,n c † n c † n c n c n in the effective tight-binding Hamiltonians of the main text. Even for our modest 21-mode system, the time-dependence of this problem would become intractable for particle numbers well below the ∼ 10 5 used in experiment.

RAMP PROCEDURE
For the localization studies of the zigzag lattice of Fig. 4 and 5, we slowly load into the extremal eigenstates (ground state or highest excited state) of the Hamiltonian we wish to explore (Eq. (2)). We begin, as described in the main text, by preparing with high fidelity the ground state of the system in the zero tunneling (t = t = 0) limit by shifting the Aubry-André site energy distribution such that the initial site has the lowest energy (ϕ = π for ∆ > 0). To ensure that there is a relatively large energy gap from this initially populated ground state to all other eigenstates even after finite tunneling is introduced, we add an effective potential well of depth V at the central site. 2 Then, over the course of 2 ms, we smoothly vary the system parameters until we reach the desired Hamiltonian. If these ramps are quasistatic adiabatic, such that the energy associated with the ramp rate is much smaller than the smallest energy gap encountered, then this procedure should prepare the desired ground state with high fidelity. First, we describe the different ramps used in experiment. The depth of the potential well at site n = 0 is ramped from V to zero over 2 ms (for comparison, the NN tunneling time for this experiment was /t = 304(4) µs). Over the same 2 ms duration, we also ramp both the NN and NNN tunneling amplitudes linearly from zero to their final magnitudes t and t , respectively. Over the course of this ramp we preserve the flux distribution, imposed by fixed tunneling phases. We additionally preserve the ratio of t /t by ramping the field strengths of the first-and second-order Bragg spectral components (∝ τ and ∝ √ τ ) according to their distinct scalings with the applied field strengths (Eqs. (S1) and (S2)). One complication arises from the small spacing between applied spectral components (frequency teeth), as mentioned in the "off-resonant excitations" section above