Coherence and decoherence in the Harper-Hofstadter model

We quantum-simulated the 2D Harper-Hofstadter (HH) lattice model in a highly elongated tube geometry -- three sites in circumference -- using an atomic Bose-Einstein condensate. In addition to the usual transverse (out-of-plane) magnetic flux, piercing the surface of the tube, we threaded a longitudinal flux $\Phi_{\rm L}$ down the axis of the tube This geometry evokes an Aharonov-Bohm interferometer, where noise in $\Phi_{\rm L}$ would readily decohere the interference present in trajectories encircling the tube. We observe this behavior only when transverse flux is a rational fraction of the flux-quantum, and remarkably find that for irrational fractions the decoherence is absent. Furthermore, at rational values of transverse flux, we show that the time evolution averaged over the noisy longitudinal flux matches the time evolution at nearby irrational fluxes. Thus, the appealing intuitive picture of an Aharonov-Bohm interferometer is insufficient. Instead, we quantitatively explain our observations by transforming the HH model into a collection of momentum-space Aubry-Andr\'{e} models.

Understanding how and when closed quantum systems lose or retain coherence is a central intellectual and practical question for quantum technologies. For example, modern optical atomic clocks operate in highly optimized decoherence-free subspaces created by using "clock" states [1] that are insensitive to the environment, as well as using lasers at magic wavelengths and polarizations that give only common mode energy shifts. In rare cases, such as collisional narrowing [2] or environment assisted tunneling [3], random processes can enhance coherence. Here we add to this list the quasi-periodic lattice described by the Harper-Hofstadter (HH) model [4,5] in a highly-elongated tube geometry-a 1D quasi-crystalby showing that the dynamics can be made immune to environmental noise.
Ultracold atomic gases in optical lattices can mimic Aharonov-Bohm (AB) phase factors using the optical phase of interfering laser beams [6,7]. Even in units of the magnetic flux quantum Φ 0 = h/q (for Planck's constant h and charged q), these systems realize large tunable magnetic fluxes Φ = a 2 B/Φ 0 per lattice plaquette (lattice constant a). Planar geometries [8][9][10], narrow Hall ribbons [11,12] and even tubes [13,14] have been realized in experiment. In the tube geometry, the longitudinal flux Φ L threading the tube has significant physical consequences. For example, adiabatically ramping Φ L by one flux quantum would drive one cycle of Laughlin's topological charge pump [15] that can probe both non-interacting and many-body topological systems [16].
The HH model, initially formulated to describe electrons moving in a 2D crystalline lattice subject to a transverse magnetic field; in terms of the AB phase Φ the HH model takes the form For isotropic tunneling J s = J x , the resulting "Hofstadter-butterfly" energy spectrum was one of the first quantum fractals ever predicted [5]. In planar geometries the additional uniform Peierls phase φ has no physical consequence, however, for a tube M sites in circumference [M = 3 depicted in Fig. 1(a)], the uniform Peierls phase contributes M φ/(2π) to the longitudinal flux Φ L [17]. Using the synthetic dimension approach [11,12,18], we assembled our 2D lattice by combining the sites of a 1D optical lattice with three internal atomic states to respectively define the longitudinal (e x ) and azimuthal (e s ) directions of our tube. We perform AB-like interference experiments [ Fig. 1(a)] in which particles prepared in a wavepacket at site m = 1 along e s , but extended along e x , are released and potentially interfere as they rapidly encircle the tube.
In this Letter, we report three key observations summarized in Fig. 1(b), where · denotes the average over φ uniformly sampled from [0, 2π), and · t marks the timeaverage. (i) For rational transverse flux Φ = P/Q (expressed in reduced form), the time-evolving population in each m-site depends strongly on φ and therefore exhibits large uncertainties, as one would expect for an AB interferometer. (ii) This dependence decreases with increasing Q and vanishes for irrational Φ. (iii) The φ-averaged dynamics at rational Φ is equal to the that at nearby irrational Φ. In all cases, our numerical simulations are in excellent agreement with both mean evolution and φsensitivity. While our experimentally probed Φ near 2/3, these observations are generalized by the numerical simulation shown in Fig. 1(c) that plots the time-averaged sensitivity to φ; this curve is approximated by an everywhere discontinuous Thomae-like function [19]. The spatial extent w limits the degree to which Φ can be distinguished to ≈ a/w, broadening the otherwise singular peaks.
Implementation We performed these experiments using 87 Rb BECs in the |5S 1/2 , F = 1 electronic ground state manifold in a crossed optical dipole trap [20]. The longitudinal Thomas-Fermi radius R TF = 11.5(5) µm was averaged over φ at Φ = 2/3 (i) and Φ = 2/3 − 1/141 (ii), along with their difference (iii). The pink uncertainty marks denote the associated standard deviation from the mean using 14-22 experimental runs and result from the sensitivity to φ along with technical noise from: variations in magnetic field, Raman coupling strength, and rf field strength, along with statistical detection uncertainties. The dark blue curves and error bands plot our numerical simulations. (c) Numerically simulated phase noise sensitivity quantified by the normalized variance, and time-averages were performed over the range [0, 1.8] ms. (d) Schematic of Raman coupling scheme and top view of the experimental setup. The optical lattice beam was retro-reflected to form a standing wave. The Raman beams were orthogonally polarized and colinear with the counterpropagating lattice beams.
obtained both by direct in-situ imaging and mean-field driven expansion [21]. We used continuous dynamical decoupling (CDD) to eliminate the magnetic field sensitivity [22,23] of this manifold's three |m F states, giving three dressed |m states that served as sites along e s . Our experiments took place in a B 0 ≈ 31.47 G bias field, along with a resonant 22.1 MHz radio-frequency (rf) magnetic field B rf with 150 kHz Rabi frequency for CDD. The resulting energy differences δ m ≡ m − m+1 in the CDD basis were (δ 0 , δ 1 , δ 2 ) = h×(−308.3, 118.6, 189.7) kHz. We adiabatically loaded BECs into the ground state of a 1D optical lattice generated by a retro-reflected λ L = 532.008(5) nm laser beam along e x . This lattice defined our tube's longitudinal sites |n with spac-ing a = λ L /2, and the V = 5.0(1)E L lattice depth set the tunneling matrix element J x = 0.066(2)E L . Here E L = 2 k 2 L /(2M a ) and k L = 2π/λ L are the single-photon recoil energy and wavevector respectively, for particles with atomic mass M a , A pair of "Raman" laser beams counter-propagating along e x [ Fig. 1(d)] (wavelength λ R , and recoil wavevector k R = 2π/λ R ) resonantly Raman coupled the |m states with strength Ω R = 0.296(6)E L , providing hopping J s = 0.111(3)E L along e s . Because |m = 2 was Raman-coupled to |m = 0 , we adopt periodic labels, i.e., |m ≡ |mod(m, 3) . Since Φ = k R /k L [24], we tuned Φ = 2/3 + ∆Φ by varying the Raman wavelength from 770.94(1) nm (∆Φ = 2/87) to 806.46(1) nm (∆Φ = −1/141), with range limited by the increasing power requirement as the detuning from the excited states increased. While φ could be tuned by changing the phase of any of the Raman lasers or displacing the optical lattice, we instead phase-shifted the CDD field B rf , giving the same effect in the CDD basis [25]. Each experimental run randomly sampled a Peierls phase φ uniformly distributed from 0 to 2π.
In the data presented here, we began with BECs in |m = 1 and initiated dynamics by abruptly introducing J s for a time t up to 1.8 ms, at which time the lattice, Raman, and dipole trap lasers were simultaneously extinguished. During the subsequent 21 ms time-of-flight (TOF) we applied a magnetic field gradient to spatially separate atoms in the three |m states. We then used absorption imaging to detect the resulting density distribution, yielding the longitudinal momentum distributions of each |m state. Figure 2(b) shows representative TOF data after time evolution for ∆Φ = 1/84 and initial state |m = 1 (middle row). Each crystal momentum state of the longitudinal lattice consists of momentum states that are imaged as horizontally spaced diffraction orders with spacing 2 k L . The synthetic lattice sites are resolved vertically; the diffraction orders in these sites are shifted by 2 k L × Φ for each synthetic hopping event as indicated by Fig. 2(a). The elongated diffraction orders in m = 0, 2 are present only for ∆Φ = 0, and result from the incommensurate lattice and Raman recoils.
Physical picture We present an intuitive picture explaining when the dynamics depends on φ, focusing on the case relevant to our experiment, with M = 3 and initial state m = 1. For each m state, Fig 2(a) plots the sinusoidal band structure of the longitudinal lattice potential as a function of crystal momentum q. Given the BEC's RMS spatial extent w ≈ 0.46R TF = 20(1)a (see SM), the corresponding momentum-space density distribution had RMS width w q ≈ /(2w) = 4.0(2) × 10 −3 × 2 k L , schematically indicated by the grey Gaussian in the left panel. The azimuthal tunneling from Raman coupling induces transitions (tan arrows) that change the crystal momentum by 2 k L Φ, as well as imparting the phase φ. After evolving under this coupling a period of time, some amplitude will remain in the ini- tial state, while some will undergo three transitions, returning to the initial |m state (right panel), with crystal momentum shifted 3∆Φ × 2 k L . This closes a φsensitive momentum-space interferometer, provided the final wavepacket has non-negligible overlap with the initial state, i.e., 3|∆Φ| × 2 k L w q .
Model We make the intuitive picture quantitive by Fourier transforming the 2D HH Hamiltonian along e x givingĤ = q0Ĥ AA (q 0 ) over Hamiltonianŝ each labeled by crystal momentum q 0 . EachĤ AA (q 0 ) is a realization of the 1D Aubry-Andrey (AA) lattice [26] with nearist-neighbor hopping strength J s and sinusoidal potential with depth 4J x , and phase set by q 0 . The sites of this AA lattice |j ≡ |m 0 + j, q 0 + jΦ × 2 k L are labeled by azimuthal site-index m along with longitudinal crystal momentum q. As shown in Fig. 2(a), the sinusoidal potential originates from Raman transitions changing the crystal momentum by 2 k L Φ as m is incremented, in effect sampling the lowest band of the longitudinal lattice. For rational Φ, eacĥ H AA (q 0 ) describes a ring [ Fig. 2(c)] of size N AA = LCM(M, Q) (LCM denotes the least common multiple) since |m 0 + N AA , q 0 + N AA Φ × 2 k L coincides with the initial state |m 0 , q 0 . For irrational Φ (incommensurate potential) the AA model is 1D quasicrystal with a prototypical metal (J s > J x ) to insulator (J s < J x ) transition [27,28], and at criticality J s = J x it is a quantum fractal showing features of quantum chaos [29].
Together with Φ, the momentum space width w q specifies the initial occupation of AA lattice sites: the mapping in Fig. 2(a) shows that every third AA site potentially samples the initial wavepacket until |3j∆Φ| × 2 k L w q when the site's crystal momentum falls outside the initial wavepacket. A representative AA lattice is shown in Fig. 2(c) where the grey bars result from sampling the filled grey curve in Fig. 2(a). In the AA model's metallic phase (J s > J x , where our experiments took place) amplitude ballistically expands from each initially occupied site, while in the insulating phase amplitude remains exponentially localized near each initial sites [28]. When 3|∆Φ| × 2 k L w q , the AA ring contains multiple initially occupied sites; as the system evolves in the metallic phase, amplitude originating in different sites can overlap and interfere. In contrast, for 3|∆Φ| × 2 k L w q only isolated sites are initially occupied and the interference is absent. Any interference depends strongly on φ, which appears in the tunneling term of Eq. (2).
We quantify the degree of interference by the normalized variance var[P m (t)]/ P m (t) 2 in each HH m-site. The analytic treatment in the SM shows that for 3|∆Φ|× 2 k L w q this normalized variance is proportional to |f (3∆Φ × 2 k L )| 2 , where f (δq) ≡ dqψ * (q + δq)ψ(q) is the form factor expressing the overlap integral anticipated by our intuitive description, and ψ(q) is the momentum-space wavefunction.
Discussion Figure 3 summarizes our data for three ∆Φ and compares the experimental data with our numerical modeling. Figure 3(a) depicts a control case where we converted our tube geometry into a ribbon by removing the Raman coupling between |m = 0 and |m = 2 . The top row depicts the AA lattice along with the initial state (grey bars), showing the transition from multiple sources for ∆Φ ≈ 0, to a single isolated source for ∆Φ = 1/84. In the ribbon case, multiple sources are present, but the disconnected links in the lattice prevent any potential interference, and the φ-sensitivity vanishes for all transverse flux Φ. The accompanying numerical simulations used parameters obtained from fits to the data in Fig. 1(b), lie within the calibration uncertainties (see SM). To obtain quantitative agreement with our data, we include a phenomenological dephasing parameter obtained from fitting the decaying sinusoid of Fig. 3(a).
The second row of Fig. 3 plots the three diffraction peaks in the range −1/2 < q/(2 k L ) < 1/2 observed after an evolution time t = 1.506 ms. When ∆Φ was small enough to allow interference between atoms originating FIG. 3. Representative data Each column depicts a different value of ∆Φ, from 0 (showing significant φ-dependence) to 1/84 (showing negligible φ-dependence), along with a ribbon geometry as a control case. The top row depicts representative AA lattices |j = 0 = |m0 = 1, q0 = 0 for each ∆Φ, along with the initial wavepackets calculated from the experimentally measured Thomas-Fermi radius (grey bars). The AA lattice schematic for ∆Φ = 0 was calculated with ∆Φ = 1/10000 since the AA ring only has 3 sites at ∆Φ = 0. The middle row shows TOF data for each case measured at t = 1.506 ms. Since diffraction orders with q/(2 kL) / ∈ [−1/2, 1/2) are replicas of those with q/(2 kL) ∈ [−1/2, 1/2), with amplitude governed by the Wannier orbitals of the longitudinal lattice [30], we focus on this regime. The bottom row plots the time evolution (top: data, bottom: numerics), and each pink data point results from an experimental run. The shaded regions marks t < 0.185 ms after which atoms have tunneled in excess of one site.
from different AA lattice sites [ Fig. 3(b)], the corresponding crystal momentum difference between different paths to the same order is too small to be resolved. When 3|∆Φ| × 2 k L w q , the momentum difference associated with the different AA became resolvable, causing each order in the TOF image to fragment into multiple overlapping sub-orders. In the ribbon case, only three orders are present irrespective of ∆Φ.
The third row of Fig. 3 compares the observed time dependence (top) with the prediction of our model (bottom). Each symbol marks the observed fractional population in the initial state P m=1 for a single measurement (no averaging), and the bottom panel is a histogram of the predicted trajectories for all φ. The experimental data is highly variable only for small ∆Φ with a manifestly non-Gaussian distribution. The noisy dynamics begin at t ≈ π /(3J s ) = 0.185 ms when population was significantly transferred out of the initial state, and interference became possible [31]. The numerical model fully captures the observed spread of data, its timedependence and even the non-Gaussian distribution, including singular features, i.e., caustics [32]. The ribbon case lacks interference and exhibits only technical noise. Figure 4(a) plots the time-averaged normalized vari-ance for fully connected rings (black data) as well as ribbons (grey data). A sharp peak is present at Φ = 2/3 for the tube geometry, while the signal is featureless for the ribbon geometry. The measured time dependence of the normalized variance, shown in Fig. 4(b), is peaked at ∆Φ = 0 for all times, but with variable amplitude. In both cases, the numerical calculations [ Fig. 4(c), and solid curve in Fig. 4(a)] are nearly indistinguishable from the data. Lastly, Fig. 4(d,e) plots the φ-averaged time evolution both for experiment (d) and numerics (e); neither show any feature at ∆Φ = 0 where the noise feature is maximal in Fig. 4(b,c). This can be directly understood in the Hall-tube picture as follows: the HH describing the system in the neighborhood of site n can be transformed to that at n = 0 by a suitable change of φ. For irrational transverse fluxes, the system uniformly samples φ, leading to spatially self-averaged time evolution. Although individual systems at rational flux lack the spatial selfaveraging effect, averaging over φ recovers the uniform sampling in the irrational case, resulting in similar mean time evolution as nearby irrational fluxes.
Outlook The realization of the HH model in the highlyasymmetric geometry proved an ideal testbed of the phe- Our AA model makes additional predictions that go beyond the experimental observations presented here. At longer evolution times than in our experiments, AA rings in the metallic phase with just one initially occupied site can exhibit long-time interference as amplitude fully encircles the ring. While this work did not address the question of giving φ explicit time-dependence, a simple argument suggests that sufficiently slowly varying φ(t) will still leave the dynamics unchanged for irrational ∆Φ near Φ = 2/3 and in the AA insulating phase. In the momentum space AA picture, φ(t) is equivalent to the potential jdφ/dt (i.e., a uniform electric field). We might expect that the dynamics will be unchanged, provided that both J s and J x are large compared to the energy shift from the electric field ξdφ/dt over the extent of the the AA localization length ξ. If true, this would be a new type of quasidisordered-stabilised decoherence-free subspace.

ACKNOWLEDGMENTS
We thank R. P. Anderson for assistance in the early stages of the experiment. We are grateful to Q. Zhou for providing Ref. [28]  Each experiment began with a Bose-Einstein condensate (BEC) of 87 Rb in a far-detuned crossed dipole trap, whose trapping frequencies were (f x , f y , f z ) ≈ (44, 53, 159) Hz in |F = 1, m F = −1 sublevel, where F and m F are the total atomic and magnetic angular momentum, respectively. Owing to the central role of the system size in our discussion, we independently determined the longitudinal Thomas-Fermi radius R TF of the BEC using both TOF and in-situ images. We used the Castin-Dum equations [S1] to obtain R TF = 11.6(5) µm from TOF images. We used digitally refocused in-situ images to give R TF = 11.4(6) µm. Because these two measurements share the same magnification, their uncertainties are correlated, and take the average R TF = 11.5(5) µm.
The rf field for the dynamical decoupling B rf [ Fig. 1(d) in the main text] was assigned a randomized phase right before the atoms were adiabatically loaded into the |m = 0 pseudo-spin state. This rapid adiabatic passage maped the magnetic sublevels |m F = −1, +1, 0 to |m = 0, 1, 2 pseudo-spin states [S2], which are insensitive to the magnetic field and facilitate cyclic Raman coupling. The magnetic field B 0 and rf strength Ω rf were stabilized using the procedure detailed in Sec. S1 C. After that, we applied an rf π-pulse, along the direction orthogonal to both B rf and B 0 , to prepare the atoms in the |m = 1 state. We removed any remaining atoms in |m = 0 by transferring them to the F =2 hyperfine manifold and blowing them away with a resonant light pulse. Meanwhile, the atoms were adiabatically loaded into the 1D optical lattice with a ramp duration of 200 ms. We then switched on the Raman coupling for various time t, followed by simultaneous sudden extinction of the Raman, optical lattice and dipole trap beams, to study the dynamics.
To achieve the momentum and pseudo-spin resolved imaging, we reversed the process of loading the atoms into the rf dressed states. Namely, during a time-of-flight of 21.3 ms, we adiabatically ramped down the magnetic field, followed by ramping down Ω rf to zero. As a result, states |m = 0, 1, 2 mapped back to |m F = −1, +1, 0 , respectively, which were then separated using the Stern-Gerlach effect with a bias field perpendicular to the lattice beam propagation direction. The absorption image in Fig.2(b) in the main text does not show the |m states in the same order as observed in the experiments (|m = 1 and |m = 2 states were switched). Additionally, due to the high magnetic field gradient applied, a harmonic potential perpendicular to the bias field stretches |m F = +1 and compresses |m F = −1 states by 3 ∼ 4%. The optical density profiles of |m F = ±1 sublevels were rescaled with linear interpolation such that the distances between the neighboring orders, which should all be the two-lattice-photon recoil momenta, is the same for all the sublevels.

B. Raman setup
The Raman beams were almost colinear with the counter-propagating optical lattice beams, with the lattice beam bisecting the angle β = 0.34 (6) • between the two Raman beams. In this configuration, the transverse flux per lattice plaquette is Φ = k R cos[β/2]/k L ≈ k R /k L . The tiny angle β was introduced to avoid retro-reflected beams. The Raman beams were carefully aligned such that there was no momentum transfer along the transverse direction so the dynamics of our system are essentially one-dimensional.
For phase coherence between Raman beams, optical lattice beams and B rf , we detected the beatnote of the two Raman beams in the vicinity of the retro mirror of the optical lattice. There were two or three frequency components in the beatnote corresponding to the two or three Raman transitions. The |m = 0 ↔ |m = 1 transition has the largest intensity and is relatively far away from the other peaks, and was therefore chosen to be the locking frequency. The local oscillator was taken from the same Direct Digital Synthesized (DDS) signal generator as B rf , such that the rf field for the dynamical decoupling and the Raman beatnote can be in phase. The feedback was applied to the the Raman beam that only contained one frequency component [Raman 1 in Fig. 1(d) in the main text]. In the rotating frame, the frequency differences between the two Raman beams are δω m,m+1 , where m = 0, 1, 2. In the lab frame, they are δω m,m+1 − ω rf , so the frequencies of the beatnotes are ∼ 22 MHz [S3].

C. Magnetic field and rf locks
We then applied two microwave pulses of duration T µw = 100 µs, separated by 1/60 s to partially transfer the atoms to |F = 2, m F = −1 sublevel. At the locking point, the two pulses are blue and red detuned by 1/(2T µw ) from the resonance, and each pulse transfers ∼ 5% of the atoms. The lock point B lock ≈ 31.37 G was chosen such that the rf frequency ω rf /(2π) = 22.1 MHz is resonant with the |F = 1, m F = −1 ↔ |F = 1, m F = 0 transition. The error signal was the imbalance of the two transfers, OD int , i = 1, 2, is the integrated optical density of the two images, and was fed to the next experimental run, setting its magnetic fields. Our procedure is extends earlier work [S4], and for further specifics see Ref. [S5].
After the atoms were transferred to |m = 0 state, we used a similar technique to stabilize the strength of the rf field Ω rf . The microwave pulses drove the transition that has the largest sensitivity to Ω rf and almost least sensitive to the magnetic field, namely, from |m = 0 to the state that asymptotically goes to |F = 2, m F = −2 with decreasing magnetic field.

S2. SIMULATIONS AND REDUCTION TO HH MODEL
We start with the full Hamiltonian describing the light-matter interaction of our three-state atoms of mass M a including a state independent confining potential V ext (x). In our experiment, we measured the V = 5.0(1)E L lattice depth by suddenly applying the lattice potential and fitting the resulting Kapitza-Dirac time evolution [S6]. We obtained the Raman coupling strength Ω R,m = Ω R = 0.296(6)E L by separately measuring the Rabi frequencies of each transition and adjusting them to be equal within our uncertainties. This process was repeated each time we changed the wavelength of the titanium sapphire laser producing the Raman laser beams. The parameters δ m = 0.00(2)E L describe detuning from Raman resonance. The stated uncertainties include the variation across the measured transverse fluxes. The ground-band behavior ofĤ full can be approximated by the tight binding HH Hamiltonian where in analogy to the relation between the Planck constant h and ≡ h/(2π), we defineφ ≡ φ/(2π). We connected these two models by numerically solving the first term in Eq. (S1) to obtain the ground-band Wannier orbitals w 0 (x − na) for each lattice site n, from which we obtained the longitudinal tunneling strength We then projected second term in Eq. (S1) to the lowest band subspace giving the synthetic dimension tunneling strength where the integral defines a Lamb-Dicke suppression factor, equal to 0.75 for our lattice depth.

A. Displacement property of HH Hamiltonian
The Hamiltonian Eq. (S2) at longitudinal site n is identical to the Hamiltonian at site n = 0 with the Peierls phase factor changed toφ = Φn +φ. In terms of the displacement operator that shifts the longitudinal sites as |n + δn =D(δn) |n , the displacement property of the Hamiltonian iŝ This property was used in the main text to explain the similarity between the time evolution at rational transverse flux Φ and nearby irrational fluxes, as well as the suppression of uncertainty in the time evolution at irrational fluxes.

B. Simulations
All of the numerical simulations presented in the main manuscript were complete real space 1D discrete variable representation (DVR) simulations of the full light-matter Hamiltonian [S7]. Here we compare full-system simulations to those of the reduced Hofstadter model, and show: (1) their qualitative time-evolution is the same, but after long times their behavior differs quantitatively; and (2) the predicted peaks in normalized variance as a function of Φ are indistinguishable. We used a real space DVR method rather than momentum space band structure approach for two reasons: (1) band structure simulations are not possible for irrational flux, since there is no periodic potential; and (2) our observed noise signature depended critically on the spatial extent of our system.
Our simulations of non-interacting atoms included a harmonic potential V ext that served to define the spatial extent of the ground state wavefunction with harmonic oscillator length HO . Because the squared form factor |f (q)| 2 governs the coherence peak width, we selected the frequency of the harmonic potential frequency so the RMS width of |f (q)| 2 from the numerical simulation was equal to that resulting from n(x) ∝ [1 − (x/R TF ) 2 ] 2 , the 1D profile resulting from integrating a 3D Thomas-Fermi profile along both transverse directions. This yields the harmonic oscillator length HO = 20 √ 2πR TF /77 ≈ 0.651R TF giving a spatial density profile with RMS width w = HO / √ 2 ≈ 0.460R TF . The associated widths in momentum space are κ HO = 1/ HO with a momentum-density RMS width w q = 1/(2w). With the definition of the potential, every term in both the light matter Hamiltonian and the HH Hamiltonian are fully defined.
We validated our calibrations by performing least squares fits to the data in the middle panel of Fig. 1(b), including P 0 (t) and P 2 (t), for a ring-coupling geometry at ∆Φ = −1/141, as well as the data in Fig. 3(d) for the ribbon geometry. The resulting best fit coefficients (Ω R,0 , Ω R,1 , Ω R,2 ) = (1.04, 0.98, 0.98)Ω R and (δ 0 , δ 1 , δ 2 ) = (0, 0, 0.02)E L are consistent with the uncertainties of our independent calibrations. All of our simulations include a single phenomenological fitting parameter τ = 1.5(2) ms, obtained only from the fit to the ribbon geometry data, to capture the slow decay of the observed coherent evolution.
We conclude with a side-by-side comparison of the dynamics of these two descriptions, as depicted in Fig. S1, with the full light matter simulation in (a) and the HH description in (b). Firstly, the time-averaged noise variance, plotted in the first row as a function of Φ, is indistinguishable confirming that our key observation is a property both of our true physical system as well as the reduced HH model. The second and third rows display the time evolution of the initial state probability P m=1 (t) and the noise variance, respectively. These data show that the exact time evolution of these simulations share the same dynamical time scale and qualitative features, but differ markedly in the quantitative evolution. Still in both cases, the probability P m=1 (t) shows no feature associated with the peak in the noise variance.
To quantitatively understand the variance of the time evolution in our experiments, we reduce the 2D HH model to 1D AA model. We begin by Fourier transforming the real dimension to realize a series of AA lattices in momentum space: In the main text, we approximated irrational transverse fluxes by nearby rational numbers within the finite resolution of the system. Here, we instead approximate rational numbers by neary irrational numbers. In this case, each AA HamiltonianĤ AA (q i ) describes an infinite chain, and the total HamiltonianĤ sums over an infinite set of such AA chains, where the crystal momentum q i is written in units of two-photon recoil momenta 2 k L . Here, both the pseudospin and crystal momentum are labeled periodically, namely |m + M = |m and |q + 1 = |q . Each chain contains a countably infinite set of sites labeled by |j; q i ≡ |m, q = |m i + j, q i + jΦ , where j ∈ Z. In what follows, we take m i = 0 without any loss of generality. Figure S2 shows an example of an AA chain with Φ = P r /Q r + ∆Φ, and ∆Φ = 1/84 + ( √ 5 − 1)/15000, close to the simple rational fraction P r /Q r = 2/3. This chain is nearly indistinguishable from the corresponding AA ring, with

C. Mean time evolution
The initial state is a wavepacket of width w q in momentum space with initial pseudo-spin state |m i = 0 : which is unraveled in the AA lattices as whereψ AA 0,j (q i ) =ψ 0 (q i + jΦ) denotes the wavefunction in the chain containing site |j = 0; q i = |0, q i . The sum over j includes all sites within a chain. D(q i ) is an everywhere discontinuous function and takes the value either 1 or 0, reminiscent of the Dirichlet function, such that the integral over q i (summing over infinite chains) includes all the momenta within [0, 1) without repeating the momenta already connected within a chain.
At time t, an initial site |j; q i evolves to where α j,j (t; q i ) = j ; q i |Û AA (t; q i ) |j; q i describes the time evolution from site j to j with φ = 0. We can therefore interpret α j,j (t; q i ) as the wavefunction at time t for a particle starting in |j ; q i expressed in the |j; q i basis, with φ = 0. When φ = 0, the time evolution only contributes a phase factor φ(j − j), as can be seen from the expansion of the evolution operator:Û The initial state has pseudo-spin M j 1 , j 1 ∈ Z. Consequently, the total wavefunction of chain q i at time t evolves to with the probability to arrive in the final site |j; q i equal tõ Inserting exp[iφM (j 1 − j 2 )] φ = δ j1,j2 into the above equation, leads to the averaged probabilities: This is an incoherent sum over the probabilities of the atoms starting from each initial site weighted by the probability of traveling from that site to the final site.
Assuming the initial wavepacket is very narrow, the Hamiltonian in the vicinity of an initially occupied site |j 1 ; q i is close to that at |j 1 = 0; q i and independent of q i . This allows us to replace the evolution from site j 1 to site j with that from site 0 to j − j 1 , i.e., α j1,j (t; q i ) → α 0,j−j1 (t), leading to probability depending only on the distance traveled.
Our experiments measured the probability to arrive at the final state with pseudo-spin m f + M , ∈ Z: The final state with pseudo-spin m f + M has AA lattice site indices m f + M .
This equation depends only on the evolution of a particle starting at AA lattice site |j = 0 . We consider a small range of transverse fluxes in the vicinity of a rational fraction Φ r of the flux quantum. Within the range of propagation, the small deviation of the transverse flux is not resolved, i.e. α Φ 0,j−j1 (t) ≈ α Φr 0,j−j1 (t). We hereby conclude that the mean time evolution in the vicinity of Φ r is a smooth function of Φ.

D. Variance of the time evolution
The variance of the time evolution for uniformly sampled Peierls phase φ is defined as We calculate var m f (t) beginning We will insert Eq. (S13).
We will change the variables as ji − j i → ji and i − j i → i, i = 1, 2 We will use e iM φ(j 1 +j 2 ) The expression in the first square bracket of the last step in the above equation is

[·]
As noted below Eq. (S9), the wavefunction was unraveled to the AA chain according to the momentum distribution.
We will reverse the unraveling, i.e. dqiD(qi) (S18) where f is the atomic form factor, and S is the static structure factor [S8]. From the normalization of the wavefunction, we have f (0) = 1 and S(0) = 1.
Combining Eqs. (S17, S18) we arrive at (S19) Importantly, if we also perform a time average we de-phase the cross-terms in the double sum The Kronecker delta functions locate those cases that match conjugated terms with non-conjugated terms, and the last term avoids a double counting error. We now shift the first term in the last line to the LHS to get the time-averaged variance having used f (0) = 1. In the second line, we defined the total variance by summing over final states. Due to the narrow wavepacket assumption, the sites that could possibly contribute to the form factor f are separated by LCM(M, Q r ), rather than M (Fig. S3), where the nearby rational transverse flux Φ r = P r /Q r with P r and Q r being co-prime. In the main text and the experiments, P r /Q r = 2/3, while here we derive a general case. Therefore, Only the non-integer part of the variable matters because the crystal momentum q is periodically labeled. If LCM(M, Q r )|∆Φ| w q , then the form factor and hence the variance approach zero for j = 0. When LCM(M, Q r )|∆Φ| w q , the closest distance between initially occupied sites on the chain is LCM(M, Q r ) (Fig. S3). This leads to the final form of the total variance (3) We extend (2) by assuming from the beginning that enough time has passed so the wavepackets have spread over a range wide compared to the initial distribution, giving where we made the far-field replacement |α 0,i−LCM(M,Qr)j | 2 ≈ |α 0,i | 2 . Note that i |α 0,i (t)| 4 t is the inverse participation ratio in the AA lattice.
(4) for ∆Φ → 0, we get in the second line we assumed that the wavefunction is uniform on the scale of LCM(M, Q r ) sites (a poor assumption in the insulating phase, and even in our data there is a ≈ 2× difference in typical averaged populations) and broke S9 the first term into LCM(M, Q r ) individual terms |α 0,i−LCM(M,Qr)j−m (t)| 2 |α 0,i (t)| 2 with m from 0 to LCM(M, Q r ) − 1.
With that we had a new sum over m, that we used to complete the j sum to be on every lattice site. At which point we used the normalization condition to do this double sum, and the → is assuming a spreading wavepacket for which the i will fall to zero.

S3. AA LATTICES IN REAL SPACE
In accordance with our TOF data, we explained our observations in a momentum space picture. In this section, we will derive an analogous description in real space that gives the dual AA model, where the noise-suppression for irrational transverse flux results from spatial self-averaging. Inserting the Fourier expansion

A. Mean time evolution
The initial state, with pseudo-spin m i = 0, is a superposition of q m state and a wavepacket of width w in real space: The atoms evolve along all M chains, and the wavefunction at time t is whereÛ AA (φ) is the evolution operator governed by the HamiltonianĤ AA (φ). The probability of the final state to have pseudo-spin m f at time t is We will insert Eqs. (S32) and (S35). Also noteÛAA does not change |qm and n |n n| = I The spatial AA model and therefore its evolution operator have the same displacement property [Eq. (S5)] as the HH model. We will shift site n to 0.
The mean time evolution is We will change the variable asφ − q m + Φn →φ, whose implicit dependence on n and q m will be averaged out overφ. We will also change the variable q m − qm → δqm and complete the summation over qm.
We will make the approximation that the spatial extent of our atomic cloud is so large such that the post-propagation wavefunction overlaps with the original wavefunction.
The quantity in the bracket describes the interference of atoms initially located in sites n = 0 and δn. The small wavefunction spread approximation echoes the narrow momentum distribution approximation made in the momentumspace picture. Next, we will quantitatively analyze how the deviation of the transverse flux ∆Φ = Φ − Φ r affects the mean time evolution. If we denote the farthest site from n = 0 that has nonzero contribution to Eq. (S37) as δn max , then we require the atoms not to travel too far to resolve the difference in the transverse flux yet ∆Φδn max 1. The cosine term that modulates the on-site energy in the AA HamiltonianĤ AA (φ) [Eq. (S33)] can then be approximated as cos 2π (Φ r + ∆Φ) n +φ = cos(2π∆Φn) cos [2π(Φ r n +φ)] − sin(2π∆Φn) sin[2π(Φ r n +φ)] Plugging this into Eq. (S37), we have The approximation in Eq. (S39) requires 2J s t/ 1 in order for the higher order terms to be small. Together with the earlier approximation ∆Φδn max 1, which validates Eq. (S38), the phase factor within the bracket, linear with both t and ∆Φ in Eq. (S40), remains small and smooth.
To summarize, we obtained the mean time evolution by shifting the sites to the vicinity of n = 0. This is possible because the dynamics on different sites of the AA chains is equivalent to that on the same site under differentφ. The global transverse phase dependence (∆Φn) is averaged out by uniformly samplingφ. Consequently, the system can only resolve the transverse phase within the range of the propagation (∆Φδn max ), leading to similar mean time evolution observed across the full range of experimentally probed transverse fluxes. S11 B. Variance of the time evolution Without averaging overφ, when shifting from site n to 0,φ is changed to Φn +φ. At irrational fluxes Φ, for a large system, summing over n randomly samples the phases, equivalent to averaging overφ. Therefore, each individual time evolution is expected to evolve as the averaged evolution-a spatial self-averaging effect. On the other hand, at rational fluxes, averaging over n only allows sampling Q r different phases in each chain and a total of LCM(M, Q r ) phases including all chains. For example, at P r /Q r = 2/3 andφ = 0, only three phases 0, 2π/3 and 4π/3 are sampled, inequivalent to averaging overφ. Thus, each time evolution is expected to have its own unique trajectory.
To calculate the variance of the time evolution, we calculate the two terms P m f (t) 2 φ and P m f (t) 2 φ separately. The former follows straightforwardly from Eq. (S37): For the latter, we have × δn1,δn2 ∆qm δn 1 |Û † AA (φ)Û AA (φ + δq m,1 ) |0 δn|Û † AA (φ + Φ∆n + ∆q m )Û AA (φ + Φ∆n + ∆q m + δq m,2 ) |0 φ (S42) from Eq. (S36), where we performed an analogous procedure as in Eq. (S37). Under the assumption of small spread of the wavepacket, we made the following approximation: n ψ * 0 (n + δn 1 )ψ 0 (n)ψ * 0 (n + ∆n + δn 2 )ψ 0 (n + ∆n) ≈ n |ψ 0 (n)| 2 |ψ 0 (n + ∆n)| 2 ≡S(∆n) (S43) whereS(∆n) is the autocorrelation function, which is the Fourier transform of the static structure factor introduced in Eq. (S18). This can be explicitly seen as follows: As phases wrap around 2π, only the non-integer part of the variable inÛ AA matters. To take out the integer part of the variable and highlight the role of ∆Φ, We will change the variable ∆n → aQ r + b, where a, b ∈ Z and b = 0, 1, · · · , Q r − 1. Therefore, Φ∆n = P r a + bP r /Q r + ∆ΦaQ r + ∆Φb ≈ bP r /Q r + ∆ΦaQ, where we ignored the ∆Φb S12 term because it is small under our assumption of small ∆Φ. Ignoring the contribution of b toS, as Q r w, leads to  w, whereS is almost independent of Θ, and thereby Θ only showing up in the bracket. The relative phases in the two brackets average out by the integral over Θ, allowing averaging over the Peierls phase to act on each bracket individually rather than their product. As a result, we have where ∆Φ M,Qr ≡ ∆ΦLCM(M, Q r ). The function σ(x) ≡ nS (n/x) /x − 1 is defined so that var m f (t)/ P m f (t) 2 φ = σ(∆Φ M,Qr ). When 1/x → 0, σ(x) → dx S (x ) − 1 = 0, which can be seen from the definition ofS [Eq. (S43)] and the normalization of the wavefunction. As a result, when 1/∆Φ M,Qr w, the uncertainty of the time evolution becomes vanishingly small. This is the same condition as what we obtained in the momentum-space picture. S13 FIG. S5. Time evolution uncertainty plotted in (a) linear and (b) log scale for a Gaussian wavefunction. In both cases, the thick red curve comes from exactly evaluating Eq. (S51), while the black curve results from the small w∆ΦM,Q r expansion and the black curve results from the large w∆ΦM,Q r expansion. The vertical grey line marks 2πw∆ΦM,Q r = π/2 were we expect the validity of these expansions to crossover. The pink lines mark the maximum possible noise for quantum states drawn completely at random: the top line at 2 is for M → ∞, and the bottom line at 0.8 is for M = 3.