Stability of many-body localization in two dimensions

TL;DR

This study investigates the stability of many-body localization in two-dimensional disordered quantum systems using ultracold atoms, revealing different behaviors for random and quasiperiodic disorder as system size increases.

Contribution

First experimental analysis of MBL stability in two dimensions with large system sizes, comparing quasiperiodic and random disorder effects.

Findings

MBL crossover shifts to higher disorder with system size for random disorder

No size dependence observed for quasiperiodic disorder, indicating potential stability

Results support the avalanche scenario for random disorder in 2D systems

Abstract

Disordered quantum many-body systems pose one of the central challenges in condensed matter physics and quantum information science, as their dynamics are generally intractable for classical computation. Many-body localization (MBL), hypothesized to evade thermalization indefinitely under strong disorder, exemplifies this difficulty. Here, we study the stability of MBL in two dimensions using ultracold atoms in optical lattices with variable system sizes up to $24\times 24$ sites, well beyond the classically simulable regime. Using the imbalance as a probe, we trace the long-time dynamics under two distinctive disorder potentials: quasiperiodic and random disorder. For random disorder, the MBL crossover point shifts to higher disorder strength with increasing system size, consistent with the avalanche scenario. In contrast, with quasiperiodic disorder, we observe no clear system size dependence, suggesting possible stability of MBL in two dimensions.

Equations

$$H=-J\sum _{⟨\mathbf{i},\mathbf{j}⟩}\left({\hat{b}}_{\mathbf{i}}^{†}{\hat{b}}_{\mathbf{j}}+\text{h.c.}\right)+\frac{U}{2}\sum _{\mathbf{i}}{\hat{n}}_{\mathbf{i}}({\hat{n}}_{\mathbf{i}}-1)+\sum _{\mathbf{i}}{W}_{\mathbf{i}}{\hat{n}}_{\mathbf{i}}$$

$$W_{\mathbf{i}=(x,y)}=\frac{W_q}{4}\left[\cos (2\pi {\beta }_{x}x+{\varphi }_x)+\cos (2\pi {\beta }_{y}y+{\varphi }_y)\right].$$

$$C(x,y)=\sum _{x^{\prime },y^{\prime }}\left[⟨W_{\mathrm{r},\mathrm{q}}(x,y)W_{\mathrm{r},\mathrm{q}}(x+x^{\prime },y+y^{\prime })⟩-⟨W_{\mathrm{r},\mathrm{q}}(x,y)⟩⟨W_{\mathrm{r},\mathrm{q}}(x+x^{\prime },y+y^{\prime })⟩\right]$$

Full text

Introduction

Statistical physics successfully links the microscopic motion of numerous particles to macroscopic observables. This connection relies on the ergodicity hypothesis, which assumes that all microstates are equally probable. In quantum systems, the concept of classical thermalization is reformulated by the eigenstate thermalization hypothesis (ETH) [1, 2, 3], according to which eigenstates are locally in Gibbs states. For quantum systems with many-body interactions, a naive expectation is that the system will obey ETH, as many-body quantum interference induces ergodicity [4, 5]. Many-body localization (MBL), which is expected to occur in strongly disordered systems, became known as a striking counter-example of this expectation [6, 7, 8], defying ETH even in the presence of many-body interactions. This emergent, non-ergodic phase is of fundamental interest not only for the breakdown of quantum thermalization but also as a potential application for protecting quantum information. Experimental studies of MBL, both in one and two dimensions, have been reported in a variety of experimental platforms, including ultracold atoms in optical lattices [9, 10, 11, 12, 13, 14, 15, 16], trapped ions [17], NV-centers [18], and superconducting qubit arrays [19, 20, 21, 22]. These experiments confirmed key features of the MBL phase, like logarithmic growth of entanglement entropy [14] and level statistics change across the MBL phase transition [19].

Despite these findings, the existence of the MBL phase under random disorder remains controversial [23, 24, 25, 26, 27]. For example, MBL can be unstable in a random potential because of the thermal avalanche scenario: rare and locally disorder-free regions can act as self-contained baths, gradually melting neighboring localized regions, which ultimately leads to full thermalization of the system in the thermodynamic limit [23, 24]. On the other hand, such thermal regions may occur so sparsely at strong disorder that they fail to connect into a system-spanning thermal cluster under certain conditions. In this view, the absence of a percolating network of inclusions prevents avalanches from destabilizing the entire system, allowing MBL to remain stable [28].

The stability of the MBL phase becomes even more unclear in higher dimensions. This is the case even in two dimensions (2D), the subject of this work. While experimental work on optical lattices [10, 29] and supporting numerical simulations [30, 31] suggest the existence of MBL in 2D, other numerical studies insist on the non-existence of the MBL phase [32]. It has been argued that the non-ergodic behavior in the experiments is a transient phenomenon in a finite time and will disappear in the thermodynamic limit. This ongoing debate is partly driven by the severe challenges in numerically simulating MBL, where the exponentially increasing complexity with system size forces state-of-the-art methods to rely on approximations that can fail in strongly interacting, disordered regimes. Consequently, numerical studies of 2D MBL are currently limited to systems no larger than $10\times 10$ sites [30, 33, 32, 34]. Nevertheless, recent theories suggest that MBL in 2D can be stable in quasiperiodic potentials [35, 36, 37], where thermal avalanches can be suppressed because of the absence of rare regions. Since the rare region effect of a random disorder potential is more prominent in higher dimensions [23, 38], a direct comparison of how these two disorder types behave in 2D can address the long-standing question about the stability of MBL.

Here we experimentally study the stability of the two-dimensional MBL phase under two different disorder potentials with ultracold atoms in optical lattices. Developing a potential engineering technique in a single-site resolution, we systematically change the system size from $8\times 8$ to $24\times 24$ lattice sites and create the random disorder and quasiperiodic potential independently. This allows for a finite-size analysis of the thermal-to-MBL crossover as well as a direct comparison of disorder types in a single platform. For random uniform disorder, the crossover shifts to larger disorder strengths as the system size increases [24], which is inconsistent with a stable two-dimensional MBL phase. In contrast, for quasiperiodic disorder, the crossover is insensitive to system size over the accessible range. These measurements suggest that quasiperiodic disorder is robust to the avalanche instability [35, 36, 37]. These observations are made possible by reaching large system sizes up to $24\times 24=576$ sites and very long measurement times up to $1500$ hoppings, well beyond the reach of current numerical methods. For example, the numerical simulations in this work (see below and Supplementary Material) required a computation time of approximately five days for each disorder realization for 100 hopping times on an $8\times 8$ lattice. For $24\times 24$, the half-system entanglement entropy near the transition would be nine times as large, increasing the required bond dimension by at least the same factor and computational time by $9^{3}$ [39] to $\sim 10$ years.

Experimental Setup

In our experiment, we investigate 2D MBL using ultracold atoms in optical lattices. This platform is highly scalable in system size and well isolated from the environment, enabling long-time dynamics and making it ideal for probing localization behavior [9, 10, 29, 15, 14, 40, 41, 42, 43]. Specifically, our system is described by the two-dimensional disordered Bose-Hubbard Hamiltonian

[TABLE]

where $\hat{b}_{\bf i}$ ($\hat{b}_{\bf i}^{\dagger}$) is the bosonic annihilation (creation) operator at site ${\bf i}=(x,y)$, $J$ is the tunneling amplitude, $U$ is the on-site interaction energy, and $W_{\bf i}$ is the random or quasiperiodic disorder potential. The tunneling is restricted to nearest-neighbor sites, as indicated by $\langle{\bf i},{\bf j}\rangle$.

The disorder is generated from a digital micromirror device (DMD), where we project an image generated from the DMD to an atomic plane by using a high-resolution imaging system [44, 45]. The optical projection convolutes the disorder potentials, where we characterize the disorder potentials from atomic responses from the local DMD potential (Fig. 1a and see Supplemental Material). The disorder strength is represented by the full width half maximum (FWHM) of the disorder distribution, $W_{\text{r}}$ and $W_{\text{q}}$, for random and quasiperiodic potential, respectively (Fig. 1b and c). The disorder potential on the DMD is set as follows: in the random case, $W_{\bf i}\in[0,W_{\text{r}}]$ is sampled from a uniform distribution, and in the quasiperiodic case, $W_{\bf i}=\frac{W_{\text{q}}}{4}\left[\cos(2\pi\beta_{x}x+\phi_{x})+\cos(2\pi\beta_{y}y+\phi_{y})\right],$ where $\phi_{x}$ and $\phi_{y}$ are realization-dependent random phase offsets, and $\beta_{x}=0.721,\beta_{y}=0.693$.

The experiment begins with preparing ultracold Bosonic ${}^{7}\text{Li}$ atoms in a 2D square optical lattice, where a quantum gas microscope is adopted to have a single site-resolved imaging [46]. We take a stripe pattern as an initial state for studying the MBL, where the total system size is controlled by imposing an optical barrier that confines the atoms within $L\times{}L$ sites during the dynamics (Supplementary Information). For a given system size $L$, we perform the two independent experiments under the random and quasi-periodic disorder to investigate the dynamical dependence of the localized phase on the disorder types (Fig. 1d and e). Even if limited to finite systems, studying the system size dependence of the localized phase is essential for understanding 2D MBL in the thermodynamic limit (Fig. 1g and h). For each type of disorder, we initiate the dynamics by ramping down the optical lattice, turning on the hopping process with interactions $U/J=5.7(2)$, and allowing the system to evolve under the disorder potential. Here, the atoms in the optical lattice evolve with the tunneling energy $J/h=84(2)~\text{Hz}$ with hopping timescale $\tau=\hbar/J=1.89~\text{ms}$. The strength of the disorder is controlled by changing the incident beam intensity shining onto the DMD.

Breakdown of thermalization

Following the time evolution, we measured the atomic distribution and used the imbalance $\mathcal{I}=(n_{\rm odd}-n_{\rm even})/(n_{\rm odd}+n_{\rm even})$ as an indicator for localization. Here, $n_{\rm odd}$ ($n_{\rm even}$) represents the parity-projected occupation number of atoms on odd (even) lattice columns. This observable can quantify the amount of initial state memory so that the imbalance shows a finite value in the localized regime [9, 10, 47].

Figure 2 displays the time evolution of the imbalance $\mathcal{I}(t)$ at system size $L\times L=24\times 24$ under random and quasiperiodic potential with various disorder strengths. For both cases, the imbalance shows a clear dependence on the disorder strength. In the absence of disorder, the imbalance rapidly drops to zero within a few tunneling times, while under the disorder potential, it decays slowly and eventually saturates to a finite value after a long time evolution. The dynamics of the imbalance under random disorder potential at $W_{\text{r}}=147{J}$ are similar to those in the previous 2D MBL experiment [29]. Right after the quench, the imbalance quickly drops due to the tunneling process, and then it gradually decays because of doublon formation. This can be seen in the dynamics of the atomic density, which decays on the same time scale as the imbalance. An exponential fit of the initial dynamics yields 3.6(7)$\tau$ and 3.1(5)$\tau$ for the decay time of the imbalance and density, respectively.

Dynamics under a disorder potential can be slow, and a complete study of the relaxation dynamics might not be accessible in real experiments [48]. In our experiment, the mean density shows a long decay time, which depends on the disorder strength. For example, for the quasiperiodic disorder (Fig. 2b inset), we find the mean density gradually decays over the entire hold time at strong disorder potential, $W_{\text{q}}=61J$. Note that the density decay cannot be attributed to atom loss during the time evolution, for example, from the off-resonant scattering of the DMD potential, since its absolute value is higher under a stronger disorder potential. The slow dynamics of the doublons are also reflected in the imbalance (Fig. 2b), where an exponential fit of the imbalance shows the relaxation time to be $t\approx 4000\tau$. When the doublon dynamics reach a quasi steady state (Fig. 2b at $W_{\text{q}}=30J$), the imbalance also stays at a nearly constant value. Therefore, we call this regime the “MBL regime”, where the relaxation time scale is so slow that we might not be able to discern the true localization within our experimental time window.

System size scaling

Within the experimentally accessible time window, investigating the system size dependence of the MBL transition becomes an effective strategy to assess the stability of MBL in two dimensions. By tracking how the transition point shifts with increasing system size, we can test whether the observed MBL is consistent with a stable phase that could persist in the thermodynamic limit, or instead is a prethermal phenomenon that thermalizes at longer times (Fig. 1g and h). To probe the transition point for a given system size, we set the evolution time to its maximum value of $t_{s}=1500\tau$ and measured the imbalance at $t_{s}$, $\mathcal{I}_{s}=\mathcal{I}(t_{s})$, under various disorder strengths. As the applied disorder strength increases, the onset of imbalance is observed, signaling the emergence of the “MBL regime” for both types of disorder. The transition curves for various system sizes ($L\times L=8\times 8,12\times 12,16\times 16,20\times 20,24\times 24$) are plotted in Fig. 3a (random disorder) and in Fig. 3b (quasiperiodic disorder), respectively.

To investigate the system size dependence of the finite time transition, we extract the critical disorder through an empirical fit to the transition curve, $\mathcal{I}_{s}(W_{\text{r(q)}})=a(W_{\text{r(q)}}-W_{c})^{2/3}+b$ for $W_{\rm r(q)}>W_{c}$ and $b$ otherwise, where $W_{c}$ is the critical disorder. Fig. 3c displays the extracted critical disorder for both disorder potentials as a function of the system size. For random disorder, the extracted transition points ($W_{c,r}$) exhibit a clear shift towards higher values as the system size increases. In contrast, for quasiperiodic disorder, the transition point ($W_{c,q}$) shows no significant dependence on system size across all measured system sizes up to $24\times{}24=576$ sites. These results are consistent with numerical simulations (for small system sizes, $L=6,8$) using the time-dependent variational principle (TDVP) [39] (see Supplemental Material for details).

The distinct system size scaling observed for these two disorder types suggests that the mechanism of thermalization in disordered systems can vary depending on the specific nature of the disorder. In the avalanche scenario uniform random disorder displays an asymptotic shift of the transition point as $W^{\mathrm{th}}_{c,\mathrm{r}}(L)\sim\text{exp}(c_{1}\ln^{1/3}L^{2})$ [24] with $c_{1}=2(\ln 2)^{2/3}\approx 1.57$ [32]. This asymptotic curve, depicted as a dotted line in Fig. 3, assumes a large system size and is consistent with our experimental data at larger system sizes. By following the asymptotic curve, the avalanche instability ultimately breaks the MBL phase in larger dimensions, leading to thermalization.

In contrast, the data for the quasiperiodic disorder distribution does not follow this avalanche asymptotic curve. Instead, it aligns with a constant value. This difference indicates that these two types of disorder possess intrinsically distinct natures of thermalization with increasing system size. Our results support the hypothesis that systems with quasiperiodic disorder in our quantum simulator remain stable against avalanche instability up to sizes of $L\times L=24\times 24$.

Quantum avalanche with time

In the avalanche mechanism for a random disorder potential, the thermalization time can depend on the system size. This is because the rare region with weak disorder can generate local thermal bubbles, which act as thermal baths. After a very long time, these thermal bubbles grow and separate bubbles can merge, then the entire system can eventually thermalize [23, 24]. Specifically, avalanche theory predicts that for fixed system size $L\times L$, an ergodic grain takes time $t_{\mathrm{deloc}}(L)\sim\exp(\mathrm{const}\cdot L)$ to delocalize the entire system [24]. This implies that if the observation time $t<t_{\mathrm{deloc}}(L)$ for all considered $L$, the experimentally measured transition point will only depend on the evolution time, $W_{c}(t)\sim\exp[\mathrm{const}\cdot(\ln\ln t)^{1/3}]$, and not on the system size.

Following this prediction, we perform the experiments that investigate the time- and size-dependence of the ergodic-to-MBL crossover in random potential and plot the results in Fig. 4. The characteristic transition point at short time ($t=150\tau$) is slightly increased with larger system size (Fig. 4a), and the critical disorder with $12\times 12$ systems is similar to the numerically obtained critical point. In contrast, at a later time ($t=1500\tau$), noticeable system size dependence is observed in the transition point (Fig. 4b and c). Our observation is consistent with the avalanche theory, where at short times the size of thermal bubbles is not sufficiently large to fully thermalize the system (Fig. 4c inset), such that the critical disorder dependence on the system size becomes weak.

Conclusions and outlook

We performed experiments on the stability of two-dimensional many-body localization (MBL), reaching system sizes well beyond the numerically accessible. Our results provide experimental evidence that random-disorder MBL in 2D is fragile to avalanche instability, by showing the upward drift of the transition point as a function of system size. However, we find that MBL induced by quasiperiodic disorder remains stable. Taken together, these observations support a disorder-dependent scenario: true random disorder allows for the appearance of rare local thermal regions that create avalanches that restore ergodicity, while quasiperiodicity suppresses such rare regions. Further work might include a more detailed experimental probe of the avalanche phenomenon–for instance, one might engineer a controlled thermal bubble, i.e., a locally weakened disorder patch, and track its growth through space and time [49]. These studies could also test the stability of the 2D MBL phase in quasiperiodic potential [35, 36, 37]. In particular, an intriguing coexistence of localization and transport can be studied in the Aubry–André model, where particle transport occurs along deterministic weak potential lines [36]. Moreover, our platform enables systematic finite-size scaling in two dimensions, which allows quantitative studies of critical phenomena at the ergodic-to-MBL phase transition in quasiperiodic models [50, 51].

Note added: While completing this manuscript, we learned of a recent complementary quantum simulation on a 70-qubit superconducting processor for uncorrelated random disorder [52].

Acknowledgements.

We acknowledge discussions with Soonwon Choi, Joonhee Choi, David. A. Huse, Johannes Zeiher, Seung-Jung Huh, Kyungtae Kim, and Elmer Doggen. J.-y.C. is supported by the Samsung Science and Technology Foundation (Grant No. BA1702-06) and the National Research Foundation of Korea (NRF) Grant under Project No. RS-2023-00207974, RS-2023-00218998, and 2023M3K5A1094812. T.B.W. is grateful for support from an EPSRC ERC underwrite Grant No. EP/X025829/1. J.L. acknowledges support from ERC Starting grant QARA (Grant No. 101041435) and the Austrian Science Fund (FWF): COE 1 and quantA. A.C. acknowledges support from the EPSRC Open Fellowship EP/X042812/1. The computational results presented here have been achieved (in part) using the LEO HPC infrastructure of the University of Innsbruck.

I Experimental systems

I.1 Experimental Sequence

As a first step to create the initial state, we prepared ultracold Bosonic ${}^{7}\text{Li}$ atoms in a 2D square optical lattice with a lattice constant of $a_{\rm lat}=752$ nm [46]. The vertical confinement is provided by an accordion lattice along $z-$axis, and the radial curvature is provided by an optical trap. Before ramping up the lattice depth to the Mott insulating regime, we apply a repulsive stripe-pattern potential, which is created by a digital micromirror device (DMD) [44, 45]. The stripe pattern has two-site periodicity ($\Lambda_{\text{Stripe}}=2a_{\text{lat}}$), and a square wall boundary was imposed to confine the atoms in a finite system size $L\times{}L$ [53]. Depending on the system size, we control the atom number by ramping down the radial trap, and then we transfer the atoms to the trap formed by the repulsive DMD potential. Finally, we ramp the lattice depth up to $V_{\text{lat}}=28.7E_{r}$ to form a Mott insulator state in the presence of the patterning potential. The $E_{r}/h=12.6$ kHz is the recoil energy of the lattice, and $h$ is Planck’s constant. Consequently, the atoms are arranged in an $n=1$ shell Mott insulator when the repulsive DMD potential is absent, whereas an $n=0$ state occurs when the DMD potential is present. The Mott insulator state with stripe pattern at $L\times L=24\times 24$ was optimized to achieve unity filling $\bar{n}=0.91(4)$ by tuning the Feshbach resonance and DMD patterning beam power. This optimization process was performed for each system size to achieve the maximum initial state imbalance. Furthermore, while taking data, we monitor the initial state every 5 minutes to obtain reliable measurements. We post-selected realizations with initial states characterized by sufficient imbalance $\mathcal{I}\geq 65-0.85$, which depends on the system size.

Then the magnetic field was adjusted to $B=676.9\,\text{G}$, setting the atomic scattering length to $a=31a_{B}$ with $a_{B}$ being the Bohr radius. The dynamics is started by ramping down the two-dimensional optical lattice to $18E_{r}$, initiating atomic tunneling with a hopping energy of $J/h=84(2)$ Hz. This configuration resulted in an on-site interaction energy of $U=478(6)$ Hz, yielding a $U/J=5.7(2)$. Details about the parameter calibration are described in Section I.3. The ramp-down time of the lattice was chosen to be 1 ms, which prevents unwanted inter-band heating during the quench process. During the evolution, the atoms interacted with a disorder potential generated by the DMD, where the potential realization was randomly chosen from a pre-generated disorder pattern realization. To ensure proper averaging over disorder realizations, we prepared 500 distinct disorder realizations, from which one was randomly selected for each experimental run. After the system evolved for various durations under these conditions, the lattice was rapidly ramped up again. Finally, the parity-projected occupation number of each site was measured using a quantum gas microscope [46].

I.2 Programmable potential

In our experiments, the key element is a programmable optical potential with single-site resolution, realized by projecting patterns from a digital micromirror device onto a two-dimensional square lattice through a high-magnification (150×) imaging system [53]. Each lattice site corresponds to 8.25×8.25 DMD pixels, enabling grayscale potentials via binary dithering. We use a blue-detuned 650 nm superluminescent diode as a light source, and amplified to have 35 mW after the fiber coupling. Residual phase drift of the optical lattice ($\sim$0.5 site/hour) is stabilized by monitoring lattice phase from fluorescence images and applying feedback to the DMD pattern at the end of each experimental run. In this subsection, we provide a detailed description about the disorder for random and quasi-periodic potential, and the optical barrier.

I.2.1 Disorder Potential

After preparing the initial state, we adjusted the system parameters for the main experiment while maintaining the lattice depth strong enough to prevent atomic tunneling. The magnetic field was tuned to modulate the on-site interaction potential so that the interaction-to-tunneling ratio is $U/J=5.7$. Simultaneously, the patterning potential was switched from the initial stripe pattern to a disorder pattern using an externally triggered transition. In the experiment, we use two different types of disorders. First, uniform random disorder was generated from a quasirandom generator, with potential values $W_{{\bf i}=(x,y)}$ uniformly distributed within the range, $W_{\bf i}\in[0,W_{\text{max}}]$.

Quasiperiodic Disorder was created by imposing another periodic potential with incommensurate period. The potential is given by

[TABLE]

Here, $\phi_{x}$ and $\phi_{y}$ are random phase offsets chosen independently for each experimental realization, and $\beta_{x}=0.721,\beta_{y}=0.693$. The strength of the disorder potential was controlled by either the intensity of the laser beam or the reflectivity of the DMD’s mirrors. In addition to the disorder potential, we added the harmonic curvature compensating potential to the disorder pattern, which can deconfine the harmonic curvature caused by the optical lattice beam. The specific disorder configuration for each experiment was generated by external software and changed for each experimental run.

I.2.2 Potential barrier

The overall pattern in the DMD is composed of a boundary wall and an interior disorder potential. The boundary walls height was set to maximum intensity, having the potential height of $V_{\rm wall}/h\geq 15$ kHz, which is much larger than the tunneling energy $V_{\rm wall}/J\geq 180$. Because of the finite imaging resolution, however, we have one site broadening of the potential wall with its height $V^{\prime}_{\rm wall}/h\geq 6$ kHz. The maximum intensity of this interior disorder potential was intentionally set lower than that of the boundary walls, creating a potential gradient between the interior and exterior regions. Typically, the maximum intensity for the interior potential was $20-40\%$ of the DMD’s maximum intensity. To avoid a homogeneous potential landscape at one site off from the wall, we also introduce the disorder pattern in this region. Overall, the disorder potential pattern is covered in the $(L+2)\times(L+2)$ region.

I.2.3 Calibration of the disorder potential

To characterize the disorder potential, we study the response of the atoms to the incident DMD beam power during the superfluid to Mott insulator transition and extract the point spread function of the optical projection system. We first prepared a Mott insulator without the DMD potential with an average density of $\bar{n}\approx 1$. At unity filling Mott insulator in low temperature, the local density can approximate to $\bar{n}(x,y)=1/[1+e^{-\beta\mu(x,y)}]$, where $\beta=1/k_{B}T$ with $k_{B}$ is the Boltzmann constant, $\mu(x,y)=\mu_{0}-V_{\rm trap}(x,y)$ is the local on-site potential with the harmonic potential $V_{\rm trap}(x,y)$. Fitting the radial profile of our Mott insulator, we determined the temperature $\beta$ and the local chemical potential $\mu(x,y)$ [54]. Then, we repeat the experiment with increasing DMD beam power, where the DMD light addresses a target lattice site (Fig. S2a). The density at a given site decreased with increasing DMD potential because the potential repels the atoms (Fig. S2b). The potential from the DMD light $V^{\rm DMD}(x,y)$ effectively reduces the local chemical potential to $\mu^{\prime}(x,y)=\mu(x,y)-V_{\rm DMD}(x,y)$. From the density profile at different beam intensity, we find the $V_{\rm DMD}(x,y)$ from the condition $\bar{n}=0.5$, which means $\mu^{\prime}(x,y)=0$ or $V_{\rm DMD}(x,y)=\mu(x,y)$. Away from the target site, the density response on the addressing beam is hardly visible, and we set $V_{\rm DMD}(x,y)=0$ in such a case.

After determining the potential $V_{\rm DMD}(x,y)$ from the single-site addressing, we calibrated the single-site potential, as shown in Fig. S2c. By fitting the calibration data to a two-dimensional Gaussian function $V_{\rm DMD}(x,y)=V_{0}e^{-[(x-x_{0})^{2}+(y-y_{0})^{2}]/2\sigma^{2}}$, we obtained the broadening parameters from the finite optical resolution as $\sigma=0.7a_{\text{lat}}$, where $a_{\text{lat}}=752~\text{nm}$ is the lattice spacing. Also, the central potential in terms of the applied beam power was determined to be $V_{0}=2.8~\text{kHz}/\text{mW}$. From this result, we can estimate disorder distribution in the atomic plane by convolving disorder images with $V_{\rm DMD}(x,y)$ (Fig. 1 b,c). Here, the disorder strength is characterized by the full width at half maximum (FWHM) of the disorder distribution. With the reflectivity of the DMD being 0.4, the calibrations for the disorders are $W_{\rm r}=5.3~J/\text{mW}$ and $W_{\rm q}=3.3~J/\text{mW}$ for the uniform random disorder and the quasiperiodic disorder, respectively.

Furthermore, we perform a one-dimensional quantum walk under a random disorder potential and complete the disorder calibration. In this experiment, we prepared a single atom in a 1D optical lattice [53], where it shows a quantum random walk in a clean system (Fig. S3a). We find the density distribution is well described by the equation, $n_{l}(t)=|\mathcal{J}_{l}(2J_{\rm{1D}}t/\hbar)|^{2}$, where $\mathcal{J}_{l}(x)$ is the Bessel function of the first kind, $l$ is the distance from the initial position, $J_{\text{1D}}$ is the 1D hopping energy, and $t$ is the evolution time. By fitting this function to the distribution at an evolution time of $t=1~\text{ms}$, we extracted the hopping parameter of the 1D lattice system as $J_{\rm 1D}/h=570~\text{Hz}$ (Fig. S3b). Then, we apply a random uniform disorder to the 1D system, where Anderson localization is expected with increasing disorder strength [55]. The localization is reflected in the density at the central position $n_{0}$, where the initial state information does not spread over the entire system and is localized (Fig. S3 c, d). The experimentally measured average density $\bar{n}_{0}$ shows a good agreement with the numerical curve, in which the convoluted disorder distributions have been used.

I.2.4 Correlation of disorders

Characterizing our potential projection system, we are able to characterize the disorder potential used in the experiments. The random disorder potential, for example, is convoluted with the experimentally determined point spread function $V_{\rm DMD}$, and the homogeneous disorder distribution is changed to have a flat-top shape distribution (Fig. 1b). Moreover, the finite resolution introduces a finite correlation length for the random disorder potential. The correlation length can be obtained by studying the two-point correlation function of the convoluted disorder potential $W_{\rm r,q}(x,y)$ with the lattice index $(x,y)$ in two dimensions,

[TABLE]

For a uniform random disorder (Fig. S4a), the spatial correlation function $C(x,y)$ was found to decay exponentially with distance. A fit to the exponential function $C(x,y)=C_{0}\exp(-\sqrt{x^{2}+y^{2}}/\xi)$ (Fig. S4b) yielded a correlation length of $\xi=0.89(2)a_{\text{lat}}$. In contrast, the spatial correlation for the quasiperiodic disorder (Fig. S4c,d) doesn’t follow a simple exponential decay. Instead, it shows a clear long-range periodic correlation from its quasiperiodic nature. The finite optical resolution of the system also affects this correlation, smoothing out the distribution as seen in Fig. S4d.

I.3 Hubbard parameters

The on-site interaction energy $U$ was measured using modulation spectroscopy [56]. While this technique is robust, its validity is contingent upon the tunneling parameter $J$ being sufficiently small compared to the on-site interaction $U$. Unfortunately, our experimental parameter $U/J=5.5$, did not meet this prerequisite for a direct measurement via modulation spectroscopy. Therefore, to accurately determine the on-site interaction energy $U$ for our experimental parameters, we first measure the interaction energy in the large $U/J$ limit at large scattering length or deep lattice and then estimate the interaction energy at the experimental condition. We first measured $U$ at a large lattice depth of $V_{\text{lat}}=26E_{r}$, where the tunneling is $J/h\simeq 23$ Hz. From the modulation spectroscopy, we find the onsite interaction energy $U/h=7.34(2)$ kHz and $U/h=571(7)$ Hz at two different magnetic fields $B_{1}=732$ G and $B_{2}=676.9$ G, respectively. Subsequently, using the same method, we measured the on-site energy at a shallower lattice depth $V_{\text{lat}}=18E_{r}$ to be $U/h=6.15(2)$ kHz. From these results, we can infer the interaction energy at $B_{2}$ at $18E_{r}$ as $U/h=478(6)$ Hz.

The hopping parameter $J$ is obtained from the dynamics of hard-core bosons in two dimensions, where the strong on-site interaction is obtained by ramping up the Feshbach magnetic field to have $a_{\text{lat}}=1100a_{B}$ [57]. After preparing the stripe pattern as described in Section I.1, we ramp down the lattice depth to $V_{lat}=18E_{r}$, and study the time evolution of the atoms without a disorder potential. Fig. S5b displays the dynamics of the imbalance, and we observe a Bessel-like oscillation even with exponential decay, which is consistent with the recent experiment in a 2D homogeneous optical lattice system [58]. From the fitting of the oscillation, we obtained the tunneling energy to be $J/h=84(2)\,\text{Hz}$.

I.4 Temperature after time evolution

During the experimental dynamics, atoms can acquire heat from various external sources, including pointing and phase fluctuations of the optical lattice, as well as collisions with residual vacuum atoms and blackbody radiation. To assess this heating effect, we measured the system’s temperature near the critical point ($V_{\text{lat}}=15E_{R}$), where the system is most susceptible to external thermal perturbations. In this measurement, the DMD potential is not applied, and the temperature is obtained from the radial density distribution as discussed in Sec. I.2.3. As shown in Fig. S6, a temperature increases at a rate of ${dT}/{dt}=0.8U/s$.

II Numerical simulations

Complementing the experimental work, we also performed numerical simulations using the time-dependent variational principle (TDVP) on matrix product states (MPS) [39]. We followed the methods of Ref. [32] and used the TeNPy library [59]. We considered $L\times L$ square lattices with open boundary conditions, focusing on the system sizes $L=6,8$ and 20 disorder realizations in each case. Following a similar protocol as in the experiment, we initialized the system in a charge-density-wave state with alternating full and empty columns and simulated real time dynamics to calculate the imbalance as a function of time. Transition points were extracted by performing a power law fit $\mathcal{I}(t)=A\,t^{-b}$ to the long-time imbalance $\mathcal{I}(t)=(n_{\mathrm{odd}}-n_{\mathrm{even}})/(n_{\mathrm{odd}}+n_{\mathrm{even}})$. We note that this fit is only valid in the thermal regime (resulting in a positive $b$) as opposed to the MBL regime, where the imbalance tends to oscillate around high values (giving rise to a negative $b$).

The time evolution was generated by the Hamiltonian in Eq. (1) for $J=1$, with the local Hilbert space restricted to the three lowest on-site occupation numbers (0, 1, and 2). As in the experiment, we considered both random and quasiperiodic disorder, with the disorder potential $W_{\bf i}$ as defined below Eq. (1). At each time step, the imbalance was evaluated without applying the parity projection. We found that for an MPS bond dimension of $\chi=340$ and time step size $\Delta t=0.05$, the results for the imbalance for $L=8$ were sufficiently well converged, even for long evolution times of $t=100J^{-1}$ hoppings. Figure S7 shows the time dependence of the imbalance $\mathcal{I}(t)$ for random and quasiperiodic disorder. The data are fit to $\mathcal{I}(t)=A\,t^{-b}$ for evolution times $t\geq 70$. Figure S8 plots the extracted decay exponent $b$ as a function of disorder strength for the two disorder types. The transition points are obtained from the disorder strengths where $b$ becomes zero within error bars. We found that for the accessible system sizes of $L=6$ and $L=8$, the obtained critical disorder strengths of $W_{c,\mathrm{r}}(L=6)=20.2\pm 4.5J$, $W_{c,\mathrm{r}}(L=8)=24.1\pm 3.6J$, $W_{c,\mathrm{q}}(L=6)=22.5\pm 0.7J$, and $W_{c,\mathrm{q}}(L=8)=23.4\pm 1.2J$ are consistent with the experimental data.

In the quasiperiodic case, we observe that $b$ is significantly smaller than zero for $W>W_{c,\mathrm{q}}(L)$. As discussed above, $\mathcal{I}(t)=At^{-b}$ is not a reliable fit in the MBL regime, and we therefore identify the transition point as the value of $W$ where $b$ first reaches zero with increasing $W$. The negative values of $b$ also indicate that the feasible simulation times remain insufficient to reach the true long-time regime. To ensure numerical reliability, we also performed convergence checks and confirmed that the imbalance dynamics and associated short-term oscillations are stable with respect to increasing bond dimension. Although full access to the asymptotic regime is computationally unfeasible, the observed persistence of a finite imbalance at late times is consistent with the phenomenology of many-body localization.