Synthetic lattices, flat bands and localization in Rydberg quantum simulators

The most recent manifestation of cold Rydberg atom quantum simulators that employs tailored optical tweezer arrays enables the study of many-body dynamics under so-called facilitation conditions. We show how the facilitation mechanism yields a Hilbert space structure in which the many-body states organize into synthetic lattices, which feature in general one or several flat bands and may support immobile localized states. We focus our discussion on the case of a ladder lattice geometry for which we analyze in particular the influence of disorder generated by the uncertainty of the atomic positions. The localization properties of this system are characterized through two localization lengths which are found to display anomalous scaling behavior at certain energies. Moreover, we discuss the experimental preparation of an immobile localized state, and analyze disorder-induced propagation effects.

The most recent manifestation of cold Rydberg atom quantum simulators that employs tailored optical tweezer arrays enables the study of many-body dynamics under so-called facilitation conditions. We show how the facilitation mechanism yields a Hilbert space structure in which the many-body states organize into synthetic lattices, which feature in general one or several flat bands and may support immobile localized states. We focus our discussion on the case of a ladder lattice geometry for which we analyze in particular the influence of disorder generated by the uncertainty of the atomic positions. The localization properties of this system are characterized through two localization lengths which are found to display anomalous scaling behavior at certain energies. Moreover, we discuss the experimental preparation of an immobile localized state, and analyze disorder-induced propagation effects.
Over the past few decades, advances in the manipulation of cold and ultra-cold atomic gases rendered them into a versatile quantum simulation platform [1,2]. Indeed, several paradigmatic many-body models have been studied experimentally, including the Luttinger liquid [3], the Tonks-Girardeau gas [4] as well as Bose-Hubbard [5,6] and Fermi-Hubbard Hamiltonians [7], permitting to directly observe several predicted phenomena, such as quantum revivals [8], Lieb-Robinson bounds [9], and topological phase transitions [10].
In quantum systems, it is well-established that transport can be heavily affected by the presence of quenched disorder, a phenomenon known as Anderson localization [36]. In the presence of randomly-distributed impurities in a metal, for example, different paths taken by an electron can interfere destructively, leading to localization. In one and two dimensions, this effect is so relevant that for arbitrarily small disorder all wavefunctions are localized and transport is effectively impossible [37,38]. Since their first prediction, these effects have been experimentally observed in a range of systems, spanning electron gases [39], cold atoms [40][41][42], thin films [43] or periodically-driven nitrogen molecules [44].
Apart from the case of quenched disorder, localized states can also arise in tight-binding models from particular lattice geometries. In these cases, destructive in-terference leads to the emergence of flat bands. Models with flat bands typically allow the construction of localized eigenstates, and have been experimentally realized with cold atoms [45], photonic lattices [46], and synthetic solid-state structures [47,48]. When disorder is introduced in such systems, these pre-existing localized states couple to the dispersive, system-spanning ones and start acting like scatterers, inducing a richer phenomenology, such as localization enhancement [6], Anderson transitions in lower-dimensional systems [50], and disorderinduced delocalization [51].
In this paper we demonstrate that Rydberg lattice quantum simulators [19,28,52] permit the exploration of disorder phenomena in the presence of flat bands. We show that under facilitation conditions -when the system parameters are set such that Rydberg states can only be excited next to an already existing excitation -the Hilbert space acquires a regular lattice structure featuring flat bands. In this picture, the uncertainty of atomic positions translates into a disordered potential on the newly created synthetic lattice. Scenarios similar to these were previously theoretically analyzed in [6,50]. Here we show that they emerge naturally in Rydberg quantum simulators employing optical tweezer arrays [28,52,53]. We illustrate our findings for the case of a so-called "Lieb ladder": we analyze the scaling of the localization length and discuss the spreading dynamics of a local flat-band eigenstate under the action of different disorder strengths.
Facilitation, Hilbert space structure and flat bands-We start by considering a regular [54] lattice of N optical tweezers, each loaded with a single Rydberg atom, and with nearest-neighbor distance R 0 . A laser is shone with a frequency detuned by ∆ with respect to an atomic transition between the electronic ground state |↓ and a Rydberg level |↑ . We work here in natural units = 1. Atoms in the Rydberg state |↑ interact, at distance d, via an algebraically-decaying potential V (d) = C α /d α , with α = 3(6) for dipole-dipole (van-der-Waals) interactions (without loss of generality we choose C α > 0). Within the rotating wave approximation the Hamiltonian of this system readŝ where Ω is the laser Rabi frequency, k and m are lattice indices, d km denotes the distance between atoms in sites k and m,σ The facilitation condition is obtained by setting ∆ = −V (R 0 ), so that an isolated excited atom makes the transitions of its neighbors resonant with the laser. In the following, we consider |∆| Ω, so that non-facilitated atoms are sufficiently off-resonant to neglect their excitation. Furthermore, we require V (2R 0 ) Ω which still ensures that an isolated excitation can facilitate the production of another on a neighboring site, but suppresses the creation of additional excitations in the neighborhood. For example, in one dimension ↑↑↑| e −iHt |↑↑↓ 2 ∼ O((Ω/V (R 1 )) 2 ). In the following, we neglect these strongly suppressed transitions, effectively splitting the Hilbert space into subspaces separated by energy scales Ω. Each subspace comprises a set of quasi-resonant states separated by scales ∼ O(Ω) (see Ref. [3] for more details on this structure). Intuitively, this means that a single excitation can at most produce one more in the neighborhood, after which either the former facilitates the de-excitation of the latter, or vice versa.
Amidst all various subspaces, the simplest non-trivial choice corresponds to the one consisting of all states with either a single excitation or a single pair of excitations on neighboring sites [3,56]. Hence, as sketched in Fig. 1 for a few planar examples, a lattice structure emerges in the Hilbert space which closely resembles the real-space geometry of the tweezer arrays. These synthetic lattices are constructed via the following rules: (i) in the original lattice structure, draw the links joining nearest neighbors; (ii) identify each site with the state having a single excitation on that site. This exhausts all "one-excitation" states in the subspace; (iii) each "pair" state can be straightforwardly associated to the link joining the two excited atoms; hence, place an additional site in the midpoint of each link and associate it with the corresponding "pair" state. The links in this new-found structure now effectively represent a pair of states connected by the Hamiltonian, which can be therefore seen as a tightbinding model on a generalized synthetic lattice. In the case of a square lattice, the new structure (see Fig. 1) is the Lieb lattice and is known to feature one flat and two dispersive bands which meet with a linear dispersion at the edges of the first Brillouin zone. However, this construction is general and can be extended to any kind of FIG. 1. Left column: geometry of a square, a triangular, and a honeycomb lattice in real space. The gray dots depict the position of the Rydberg atoms and the lines the interaction between neighboring atoms. R0 and R1 represent nearest and next-nearest neighbor distances, respectively. Middle column: respective "synthetic lattices" in the Hilbert space under facilitation conditions. The blue dots represent oneexcitation states while the red ones are pair states. Right column: Cut through the Brillouin zone for each lattice geometries at ky = 0. Each lattice features (at least) a flat band. The momentum scales for the three lattices (from top to bottom) are η = (1, 4 3 , 4 3 ).
regular [54] lattice. Most of these structures will support flat bands as well: It can be shown [57] that, calling n 1 (n 2 ) the number of one-excitation (pair) states in a unit cell, the number of flat bands n flat must be ≥ |n 1 − n 2 |.
For the examples of Fig. 1, the square, triangular and honeycomb lattices have (n 1 , n 2 , n flat ) = (1, 2, 1), (1,3,2) and (2, 3, 1) respectively. These flat bands constitute extensively-degenerate eigenspaces of the Hamiltonian; as such, it is often possible to recombine the usual (planewave-like) Bloch solutions to form a set of localized (or immobile) eigenstates. Disorder-Disorder enters the picture through the uncertainty in the atomic positions. Even small displacements from the center of the traps can significantly shift the atomic transitions off resonance from the laser frequency, thereby hindering the facilitation mechanism [3]. In fact, the interaction potential seen by an atom at a distance R = R 0 + δR from an excitation will be V (R) = V (R 0 + δR) ≡ V (R 0 ) + δV . At small disorder (δR R 0 and δV V (2R 0 )) the energy shifts can be approximated by δV ≈ −αC α /R α+1 0 δR [57]. These random variables only affect pair states, creating a disordered potential landscape over the pair (red) sites in Fig. 1.
In order to characterize the disorder, we denote by ω the optical tweezer trapping frequency (assumed hereafter to be isotropic in space), by m the atomic mass and by T the temperature. The probability distribution of a trapped atom can then be approximately described as a Gaussian of width σ around the trap center. We require now that (I) k B T ω: this implies that one can use the semiclassical estimate σ ≈ k B T /mω 2 and moreover that the thermal de Broglie wavelength of the atom is much smaller than the distribution width. In other words, the atom can be approximately considered localized somewhere within the trap according to a classical probability distribution. (II) ω∆t 1, with ∆t the duration of an experiment: this ensures that the atoms will not appreciably move from their positions in this time frame and thus the disorder is quenched. (III) Ω ω, or in other words the dynamics of the internal degrees of freedom is much faster than the one of the kinetic ones, so that within an experiment one can probe the action of the disordered Hamiltonian on the system while keeping the specific realization of the disorder fixed. The properties of the probability distribution of energy shifts are discussed in [57]; here we just mention that amplitudes of shifts over different pair sites are not independent, but correlated.
Disordered Lieb ladder-In the remainder of our discussion, we shall focus on a ladder configuration, i.e. a quasi-one-dimensional lattice formed by placing two linear chains parallel to each other at a lattice spacing R 0 . For this example, the synthetic lattice (a "1D Lieb lattice") in the Hilbert space is sketched in Fig. 2(a). The unit cell consists of five sites with n 1 = 2 and n 2 = 3 and the band structure features one zero-energy flat and four dispersive bands [ Fig. 2(d)].
This Lieb ladder constitutes one of the simplest examples where flat bands produce a non-trivial interplay with the on-site disorder [6]. In a Rydberg quantum simulator, however, the disorder only appears on pair states, i.e. all the blue sites of the synthetic lattice [ Fig. 1(a)] are unaffected by it. To investigate the effect of this unusual disorder scenario we study in the following the scaling behavior of the localization length ξ for small disorder strengths. This quantity encodes the localization properties of the energy eigenstates, whose amplitude is typically peaked in a specific area of the lattice and decays exponentially as e −r/ξ at large distances r.
For a ladder like the one under study, two different values of ξ can be extracted at any given energy, which we denote by ξ 1/2 and order according to ξ 1 < ξ 2 . To elucidate the reason, one can perform an appropriate change of basis ("detangling transformation" [6,58]) through which the Lieb ladder is mapped onto two uncoupled one-dimensional lattices [see Fig. 2(b)], a chain (in green, supporting the two innermost dispersive bands) and a stub lattice (in orange, supporting the flat and two outermost dispersive bands) [57]. At small disorder, one can thus associate either localization length to one of the two detangled chains.
The localization length ξ 1/2 are found numerically via a transfer matrix formalism and are displayed in Fig. 3(a) as a function of the disorder strength s ≡ σ/R 0 and the energy . In Fig. 3(b) we display log-log plots of the correlation lengths at selected energies as functions of s, which illustrate algebraic scaling ξ i ∼ s ν , for sufficiently small s. Where possible, we connect our findings to those presented in Ref. [6], where the same geometry is studied with independent disorder on all sites. The usual scaling for Anderson localization corresponds to ν = 0 at energies outside a band ("out"), ν = 2/3 at a band edge ("edge") and ν = 2 inside a band ("in"). The energies selected in Fig. 2 correspond to = 1 (out/in), √ 2 (edge/in), 1.8 (in/in), 2 (in/edge) and √ 6 (edge/out). Here the entries in the brackets refer to the two band structures depicted in Fig. 2(c,d): (orange/green). 6 In Ref. [6] an "anomalous" scaling ν = 4/3 was found at = √ 2 and 2. This was attributed to the fact that disorder, in the detangled picture, is not merely on-site but couples the two chains. This in turn may produce resonances between states in the middle of a band and states at the edge of the other when the latter displays vanishing group velocity. Comparing these values with the ones obtained for our situation, we observe reasonable agreement at = 1, = 1.8 and = √ 6, plus for the "edge" scaling at = √ 2. The anomalous "in" scaling at = √ 2 seems instead to be "cured" as we retrieve a result compatible with the usual Anderson one (ν ≈ 2). As we show in [57], this is likely to be due to the alternating structure of the disorder in the synthetic lattice, which in the detangled picture results in the absence of random couplings between Y ± n sites [see Fig. 2(b)], present instead in Ref. [6].
We find however discrepancies at = 2, where both localization lengths are close to 1.1 and do not seem to match with either of the expected values 2/3 (edge) or 4/3 (in, anomalous). An explanation for this behavior, which does not seem to be related simply to the alternating structure of the disorder [57], is currently lacking and requires further investigations.
Localized flat band state dynamics-Experimentally measuring the localization lengths studied above is challenging due to the required large systems size and small disorder amplitudes. However, one can probe the influence of disorder by initializing the system in a specific state and tracking the subsequent dynamics by measuring the on-site excitation probabilities [19,28,52]. A particularly interesting choice for an initial state is localized and an eigenstate of the flat band. Such state is not propagating in the absence of disorder. We show in [57] that it takes the form of this form can be prepared experimentally via single site addressing [57]. The time evolution of the excitation density is shown in Fig. 4(b). The effect of the disorder becomes apparent in the width ∆x [57] of the density packet which quickly reaches a stationary state. It is interesting to observe that, as shown in Fig. 4(c), the stationary value of ∆x shows a non-monotonic behavior as a function of s. This can be understood as follows: at very small (but finite) disorder strength s the initial state (energy ≈ 0) is almost a flat band eigenstate and it therefore only minimally spreads (see e.g. Refs. [51,59]). As s is increased, this picture breaks down and the state more and more strongly hybridizes with other states, allowing transport over larger distances. At the same time, however, the localization lengths at = 0 decrease. Hence, an interplay ensues: the spreading ∆x of the state increases with s as long as the localization length remains larger (∆x ξ i ). Once the decrease in the localization scale catches up with the increase of ∆x, the behavior is dominated by localization and, as expected, decreases with increasing disorder strength.
Conclusions and Outlook-We have shown that Ry-dberg quantum simulators allow to explore localization phenomena in synthetic lattices with flat bands and unconventional types of disorder (correlated, alternating). The current study focuses on the Lieb ladder and on a particular excitation sector. Higher-dimensional lattices hosting more excitations are straight-forwardly realizable in experiment. It is thus a future theoretical challenge to shed light on these intricate and unexplored scenarios.
Supplemental Material: Synthetic lattices, flat bands and localization in Rydberg quantum simulators

Approximate Gaussian distribution of the atoms
In order to show how the Gaussian distribution of the atomic positions arises, we consider here an atom of mass m sitting in a one-dimensional optical trap of frequency ω. The results will be straightforwardly generalizable to the three-dimensional case as the three Cartesian coordinates decouple and can be treated independently. We work in a regime of temperatures T much lower compared to the trap depth, but larger than the trap frequency, i.e., k B T ω. The first assumption allows us to treat the trap as an harmonic potential, yielding a Hamiltonian wherex andp are the quantum position and momentum operators, respectively. The thermal state of the system is described by the Gibbs form where β = 1/(k B T ) and Z is the partition function Employing the standard mappinĝ in terms of bosonic creation (b † ) and annihilation (b) operators b ,b † = 1 , one readily obtains Calling |x the position eigenvectorx |x = x |x , the probability density functions of the atomic position is defined as Its analytical form can be extracted from the Feynman propagator for the harmonic oscillator K(x, y, t) = x| e −itĤ/ |y , which, in the time interval t ∈ (0, π/ω), reads (see, e.g., [1,2] for detailed derivations) Substituting t → −iβ and y → x one finds Dividing by the partition function (4) one finally finds the Gaussian The variance σ can be read off directly and amounts to .
Since we assumed β ω, i.e., ωβ 1, we can expand this expression to lowest order, which yields as reported in the main text. The distribution (10) is straightforwardly generalized to three-dimensions and traps centered along a single linear chain at positions kR 0 = (0, 0, kR 0 ) with k an integer: For clarity, we remark here that the indices in the expression above distinguish between Cartesian components only, e.g, r 1 and r 2 are the components of the same atom along the x and y directions. In the following, when necessary the trap index will always appear before the component one, e.g., r k,i is the i-th component of the k−th atom's position. For a ladder, a second set of position distributions p Here we focus on a single one-dimensional chain as most of the properties which affect the results in the main text are due to the presence of an extended longitudinal direction. Still, the considerations made for the marginal distributions for pairs of atoms directly apply to any regular lattice configuration as well. The distribution of the differences d k = r k+1 − r k = (d k,1 , d k,2 , d k,3 ) can be found in the Supplemental Material of Ref. [3] and, for isotropic traps, reads where implying, Subsequent distances are therefore (anti-)correlated, and these correlations pass onto any (non-trivial) function of the distances, and in particular the energy displacements As a consistency check, we remark that C(L) is a (L − 1) × (L − 1) matrix, whose determinant satisfies the recursion relation with "seed" (or initial conditions) det C(2) = 2 and det C(3) = 3, which is solved by det C(L) = L. Consequently, the factor √ det A 3 produced by the Gaussian integration over all variables exactly cancels the L −3/2 appearing in the normalization factor, as expected.

Marginal distribution for a single pair of atoms
The d k s are identically distributed, so we can select any given one (and drop its index for brevity) and integrate over the other L − 2 variables from equation (14). This yields where d = |d| denotes the distance between a pair of neighboring atoms. The distribution for this new variable can be then obtained via a solid angle integration and reads The distribution of an energy shift δV is now just a change of variables (d → d(δV )) away, according to . For the sake of generality, we keep α generic in where Hence, the distribution of energy shifts for a pair is It is relatively simple to see that, if we define the dimensionless quantities δv = δV /V 0 and s = σ/R 0 , we can simplify this expression further: The probability distribution function in Eq. (23) is defined in the domain δ v ∈ [−1, +∞); for δv = −1 + ε, in the limit ε → 0 + it behaves as as the (vanishing) exponential factor dominates. In the opposite limit δv → ∞, instead, the distribution behaves asymptotically as This shows that this distribution is fat-tailed. In particular, all the distribution moments δv β with β ≥ 3/α are not defined and, for both α = 3 (dipole-dipole interactions) and α = 6 (van der Waals), this includes all integer moments (e.g., the mean and variance). These fat tails are the consequence of the approximation of an atom's position distribution as a Gaussian everywhere in space, i.e., including points much further away from the center of a trap than a few σs. In other words, it appears to be an artifact of the description, rather than something occurring in a real experiment. The result of this approximation is to allow for an extremely small (but not vanishing) probability that two atoms can be arbitrarily close, which, due to the algebraic scaling of the interactions, produces considerable energy shifts. Moments like the mean and variance are therefore dominated by the very rare events in which two atoms lie very close to each other. The rarity of such events is encoded in the exponential suppression e −(1/4s 2 ) in Eq. (25). In principle, these unphysical fat tails could affect our results, as it is known that, in the Anderson problem, the scaling of the localization length is modified when Cauchy-like distributions are chosen instead of more regular ones. However, as mentioned above, the fat tails in our case are strongly suppressed and one needs to assess how likely it is to actually probe them in a simulation or an experiment. For that purpose, let us first notice that the asymptotic behavior reported in (25) emerges when the argument of the sinh function in Eq. (23) is small, i.e., still assuming δv 1, for δv Let us calculate now the probability of generating an energy shift within the tails, i.e., Employing now the asymptotic expression (25) we obtain This result apparently does not depend on α, but we need to remember that the derivation assumes δv 1, and is therefore only consistent if (2s 2 ) −α 1. Considering α = 3 or 6, though, this is satisfied already for rather large disorder amplitudes, e.g., s = 0.3, which then yields P 0.3 ≈ 0.0013. Due to the exponential factor, these probabilities decrease very fast with s. For s = 0.1, for instance, we get P 0.1 ≈ 10 −14 and in the range spanned in the plots reported in the main text s ≤ 5 × 10 −4 this becomes P s 10 −400000 , which is clearly impossible to observe in any reasonable experiment or numerical procedure. Hence, we can safely assume the unphysical fat tails to be completely irrelevant in the determination of our numerical results in the regime considered.

Distribution bias towards positive energy shifts
As can be observed from Fig. 3 in the main text, there appears to be a bias of the distribution towards positive energy shifts for small disorder, which makes the features present in the localization length plots bend towards higher energies. Considering the marginal discussed in the previous section, this seems counter-intuitive, since the distribution (23) seems to shift, for increasing s, towards negative values instead, eventually becoming peaked very close to δv ≈ −1. It is also possible to provide an intuitive explanation for this, since, as two neighboring traps becomes wider and wider, it becomes much more likely for two atoms to lie at a larger distance than the one separating the two centers. The limiting case s R 0 is indicative, as one can imagine that the atoms' positions can be picked uniformly in space on length scales R 0 .
Although we do not hold at the moment a convincing explanation of why the shift seems to point in the opposite direction, we believe it is due to the correlated nature of the full distribution (which indeed would be consistent with not seeing the correct behavior in a marginal) and we provide a physical argument which partially supports this conjecture. For simplicity, we shall work here with a chain, rather than a ladder. We hypothesize that (A) the number L + 1 of atoms is large, i.e., L + 1 1 (this choice is such that the number of distances is L); (B) the disorder is extremely small s ≪ 1. A useful simplification from (B) is that we can effectively reduce the dimensionality of the problem and only consider the position displacement along the z direction: in fact, if we write d = (δx, δy, R 0 + δz), then Hence, we can extract the probability of the distances d k directly from Eq. (14): with the same A k,q (up to increasing L by 1). Clearly, the marginal for any given variable d k is also Gaussian and its mean and variance can be straightforwardly extracted from the expression above: These variables are therefore identically distributed and, due to the boundedness of the covariance (see matrix C in Eq. (15)), satisfy a generalized weak law of large numbers, as we demonstrate below for this very special case: let us define D = ( k d k )/L and consider the probability P(|D − R 0 | > ε) of a fluctuation ε around the mean value. By Chebyshev's inequality, where where C k,q are the elements of the matrix C = A −1 in Eq. (15). This sum is not difficult to calculate, since each row but the first and the last totals 0, whereas the first and last contribute 1 each, implying VarD = 2σ 2 /L 2 . As a consequence, by choosing a sufficiently large L the r.h.s. in Eq. (32) can be made arbitrarily small or, more precisely, ∀δ > 0 ∃L : ∀L >L and therefore D → R 0 in probability. Note that, had we been dealing with independent variables, we would have retrieved a result where VarD ∼ 1/L instead of 1/L 2 . On a less formal level, this can be understood as follows: consider that, in this effective one-dimensional picture, the sum of all distances corresponds to the distance between the first and last atoms. When generating the positions independently, this variable will not be affected by the random nature of the positions of all the intermediate ones; instead, if one were to generate the distances as independent variables, these would effectuate a random walk (with drift R 0 ) and thus the effective uncertainty in the position of the last atom, assuming knowledge of the first one, would be of order O( √ L) and increase with the length of the chain. Now, consider that, having chosen repulsive interactions (V (R) > 0), the interaction potential is a convex function. Hence, we can write down Jensen's inequality as By the weak law of large numbers, for large L 1 we can effectively replace the r.h.s. of the inequality above with V (R 0 ), which leaves us with the approximate statement i.e., i.e., k δV k 0.
Albeit not a rigorous proof, this argument provides a clear indication that, for sufficiently large system sizes, the correlations among the variables will make positive biases in the energy shifts preferable to negative ones, in agreement with the qualitative features observed in the main text. Accepting this claim, there must be at least a point s > 0 where this bias is strictly positive. It then follows that there exists a right neighborhood of s = 0 in which the bias increases with s.

Hilbert space reductions and restricted Hamiltonians
The Hamiltonian introduced in the main text readŝ where d km denotes the distance between the k-th and mth atoms. In order to exploit the large energy separations present in the system, we switch to the interaction picturê Recalling that σ (k) x ,n m = 0 for every k = m and that σ (k) x we can simplify the k-th addend in where in the first equality we singled out in the exponentials all the terms which depend uponn k ; all the remaining ones cancel out. The HamiltonianĤ I can then be written aŝ We apply a rotating-wave approximation to discard all terms which oscillate fast in time. This implies that the oscillation frequency ω should be Ω for a term to be neglected. Note that the frequency ω is however operator-valued: Since the prefactor −1 ≤ 2n k − 1 ≤ 1 is of order O(1), it is the second factor which is decisive for the selection.
We introduce now for every site k a projectorP k over all states where there is a single excitation among the neighbors of k and no additional one within a radius 2R 0 . Its specific structure depends clearly on the structure of the lattice, but if we define by F k the set of nearestneighboring sites of k and by S k the set of sites within a distance 2R 0 from k which are neither site k itself nor one of the sites in F k , then we can give an implicit definition according tô (1 −n q ). (44) Checking that the expression above satisfies P k 2 = P k is straightforward if one recalls thatn 2 q =n q and (1 −n q ) 2 = 1 −n q ∀ q. The relevance of the projectorP k is that it precisely identifies the constraints -identified in the main text -under which a spin (or atom) is able to flip (or being excited/de-excited). Slightly more formally, where we have neglected all contributions from excitations beyond a distance of 2R 0 . Furthermore, note that according to definition (44)P k acts trivially on site k and thus commutes with all local operators which instead exclusively act on that site; in particular, σ (k) x ,P k = 0. Defining for brevityQ k = 1 −P k the projector onto the orthogonal subspace (Q 2 k =Q k ,Q kPk = 0) we thus havê Hence, we can separate the interaction HamiltonianĤ I into two contributions: x .
The space of configurations onto whichQ k has support can be further split into three classes: (A) States where site k has two or more excited nearest neighbors; (B) States where site k has only one excited neighbor, but there is at least another excitation within a radius 2R 0 ; (C) States where no neighbors of k are excited.
In case (A) the interaction potential on site k is ≥ 2V (R 0 ); accounting for the facilitation condition ∆ = −V (R 0 ) we find ω V (R 0 ) Ω; these terms are thereby oscillating very fast and can be discarded. Terms of type (B) are facilitated by the single neighboring excitation, but the presence of an additional one within a distance 2R 0 implies that and therefore ω V (2R 0 ) Ω, which allows us to neglect all type-(B) contributions as well. Terms belonging to class (C) are instead more delicate, since an appropriate combination of the interactions with many excitations at different distances could approximately cancel out the detuning ∆. For instance, for dipole-dipole interactions (α = 3) the potential obeys V (γR 0 ) = V (R 0 )γ −3 ; considering a honeycomb lattice with 5 excited nextnearest neighbors at distance R 1 = √ 3R 0 and a single excited next-next-next-next-nearest (or fourth-nearest for brevity) neighbor at distance R 4 = 3R 0 one finds However, configurations such as the one described above always require a large local density of excitations, and hence can only affect Hilbert subspaces at higher energies than the ones considered in the main text, separated at least by some factors of V (R 1 ) Ω. As long as we consider the low-energy Hilbert subspaces, it is thus fine to neglect terms of type (C) as well. Overall, in the subspaces we are interested in we can approximatê Going back to the original Schrödinger representation is now straightforward and yieldŝ Note that in the specific subspace (let us call it H 1 ) considered in the main text, the one including all possible one-excitation states plus all possible pairs of neighboring ones, the diagonal partĤ 0 acts trivially as the null operator and can thus be discarded, implyinĝ We remark that the same derivation can be followed in the presence of weak disorder by changing the definition ofĤ 1 in Eq. (39) tô Since the second term is diagonal and commutes witĥ H 0 , the calculation of the interaction picture is straightforward: and one can follow the same steps outlined above.

Hilbert space lattice structure
Having derived the restricted Hamiltonian (52) we can now identify the geometric structure of the Hilbert space in the basis of eigenstates ofσ (k) z . To start with, we introduce the following definitions for the basis itself: we call |M k states with a single excitation present on site k, whereas we denote by |N kq states with a pair of excitations on sites k and q. Fixing the number N of tweezers, the Hilbert subspace we work in is therefore defined as where we recall that F k is the set of nearest neighbors of site k. Note that, since |N kq = |N qk the pair states are doubly counted; however, this clearly still leads to the generation of the same vector space. Alternatively, one can also define an equivalence relation |N kq ∼ |N ml ⇔ (k = m ∧ q = l) ∨ (k = l ∧ q = m) and take the quotient of the r.h.s. above. In the following, it is understood that the states |N kq are always taken from this space, i.e., we shall never consider states with two isolated excitations at distance d > R 0 .
By construction,Ĥ H1 H 1 ⊆ H 1 . Furthermore, we know that the action ofP kσ (k) x is to flip the spin in site k conditioned on the presence of a single excitation in F k and no additional one in S k . This implieŝ P kσ (k) Considering that l ∈ F k ⇔ k ∈ F l , one can see that Similarly, since by construction the only facilitated spins are in sites q and l. Hence, Collecting these considerations, we can find the Hamiltonian matrix elements: Now, there are as many states |M k as there are sites, so it is natural to make a connection: starting from the real-space geometry of the tweezer array, which defines the original lattice structure, we place for visual aid each state |M k on the corresponding site k. Crucially, each pair state |N kq is exclusively connected (via the Hamiltonian) to the two one-excitation states |M k and |M q , so it is placed as a mid-point between sites k and q, changing the structure to a generalized Lieb lattice. Now, by drawing a link between any pair of sites every time the corresponding states yield a non-zero Hamiltonian matrix element one precisely reconstructs the kind of lattices we displayed in Fig. 1 in the main text.

Bound on the number of flat bands
We provide an account of the lower bound of the number of flat bands n flat ≥ |n 1 − n 2 | mentioned in the main text. We recall that n flat denotes the number of flat bands in the model, n 1 the number of one-particle states per unit cell and n 2 the corresponding number of pair states per unit cell. Before doing that, however, we briefly comment on the fact that the spectrum of the hopping Hamiltonians (52) is always symmetric with respect to = 0. In fact, one can define the parity transformation which, in the subspace H 1 , acts according toÛ |M k = − |M k on all one-excitation states andÛ |N kq = |N kq on all pair states. Combined with Eqs. (60a)-(60c), this impliesÛ †Ĥ H1Û = −Ĥ H1 . Hence, if | is an eigenvector of the Hamiltonian at energy , thenÛ | is also an eigenvector, but at eigenvalue − , proving the symmetry of the spectrum under reflection → − .
We start directly from the synthetic lattice reconstructed in the Hilbert space according to the procedure described in the previous section. This structure is not in general a Bravais lattice and needs, as a first step, to be reduced to one by identifying an appropriate "basis". This is a standard procedure in crystallography and solid state physics and we refer the reader to any good introductory textbook (see e.g., [5]). For the reader's convenience, we however recall here just a few of the most basic concepts: a Bravais lattice is a lattice structure where the positions l of the lattice sites can be written as discrete combinations (62) of a set of d linearly-independent primitive lattice vectors a i (i = 1 . . . d), where d is the dimensionality of the system. If a site is located at the origin, all sites can be found this way and all points at positions l are lattice sites. Any lattice is, by definition, a periodically repeating pattern, and is therefore invariant under a certain set of translations by l for some specific choice of the primitive lattice vectors. However, in many cases an additional set of B vectors b 1 , . . . b B , called "basis", is required.
In such cases, and fixing conventionally b 1 = 0 which can be done without loss of generality, if one lattice point is located at the origin, every point at a position l is also a lattice site, but not all lattice sites are at positions l. All of them are instead found at positions l + b j with j = 1, . . . , B. We also remark that distances between sites in the synthetic lattice are not meaningful, being just a convenient way to visualize the structure of the Hilbert space. Hence, we are free to rescale the length of all (dimensionless) vectors a i , b j by a common factor.
In all the examples discussed below the primitive lattice vectors have the same length and we shall choose to normalize them to unit length (| a i | = 1). Also, for brevity in the following we refer to the R d space where these vectors live as the direct space.
We also introduce the reciprocal lattice vectors a * i , i = 1 . . . d which satisfy the defining relations The reciprocal Bravais lattice is then reconstructed by taking integer combinations of these vectors, i.e., We define a unit cell U * which contains only one reciprocal lattice point. All possible translations G of U * cover the whole space R d without any overlaps. It can be visualized as a tessellation with U * a tile. From a slightly different (but equivalent) perspective, one can define the equivalence relation between vectors k, q ∈ R d living in reciprocal space with G a reciprocal lattice vector. Hence, the unit cell may be defined as a quotient R d / ∼. By defining quasimomenta k as reciprocal space vectors belonging to a unit cell U * , one can define a Fourier series in the usual way for any generic quantity A l living on the direct-space Bravais lattice The corresponding inverse transform is also standard: as can be shown remembering that and using the Poisson-summation-derived distributional identity with α ∈ R and δ the Dirac delta. The choice of the unit cell is not unique; in the following we assume to be working in the first Brillouin zone B [5]. Clearly, the definitions above do not hinge upon working in a specific space and, indeed, one can analogously define a unit cell in direct space which contains a single Bravais lattice point. Hence, such a unit cell includes B synthetic lattice points. It is quite natural to subdivide them according to whether they are of the "oneexcitation" or "pair" kind. As done in the main text, we define n 1 the number of one-excitation states in a unit cell and n 2 = B − n 1 the number of pair ones. For example, • Synthetic square lattice (Lieb lattice): n 1 = 1, n 2 = 2, B = 3.
Since each synthetic lattice point can be uniquely associated to a given primitive lattice vector l and basis vector b i , we can unambiguously denote each state in the Hilbert subspace H 1 as a tensor product l ⊗ b i . For later convenience, we introduce now a new notation distinguishing between the basis vectors identifying oneexcitation states b i → |µ j , j = 1, . . . , n 1 and pair states b i → |ν j , j = 1, . . . , n 2 , so that the space of basis states is equivalently generated as Consequently, there is a bijective correspondence between states |M k and states l ⊗|µ i and between states |N kq and states l ⊗ |ν i .
We also define the lattice translation operator T j , where j is a Bravais lattice vector, which acts on the positional degrees of freedom according to where C j and D j are collections of connectivity matrices with elements 1 (if two states are linked) and 0 (if the two states are not). For instance, if the Hamiltonian can cause a hop from l to l + a 1 accompanied by a change |µ 1 → |ν 1 , then D a1,1,1 = 1. Note that these are, in general, rectangular matrices of size n 1 × n 2 . Furthermore, to ensure that H is hermitian they must satisfy where the last equality comes from the fact that they are defined to be real (their elements being either 0 or 1). Note that no terms ∝ |µ m µ n | or |ν m ν n | appear, as one-excitation states are exclusively connected to pair ones and vice versa (see Eqs. C j,m,n |µ m ν n | + C − j,m,n |ν n µ m | T j where again is, for every k ∈ B, a rectangular n 1 × n 2 matrix. Calling nowM we can represent it as a matrix in the basis {|µ 1 , . . . , |µ n1 , |ν 1 , . . . , |ν n2 }, which yields Due to this particular block structure, Furthermore, the rank of a rectangular matrix is never greater than its shortest side. In this case, Hence, if |n 1 − n 2 | ≥ 1 then for every k one can find a kernel vector v k in the basis such thatM k v k = 0.
Correspondingly,Ĥ H1 k ⊗ v k = 0 ∀ k and the set of all these states forms a zero-energy flat band. Clearly, if |n 1 − n 2 | > 1 then more than one choice of |v |k can be made per each quasi-momentum k, each identifying an independent flat band. Hence, calling the number of flat bands in the model n flat , consistently with the main text notation, which proves the bound. The general rules for filling the matrix elements of C k are the following: • Choose n-th column 1 ≤ n ≤ n 2 .
• Consider the two possible ways in which a particle can hop from the intermediate state |ν n within the basis to its neighbors |µ m and |µ p .
• Add e i k· jn to C − k,m,n and e i k· jp to C k,p,n , where j m/p are the lattice vectors pointing to the arrival lattice sites.
In the next sections we work out some examples among the ones displayed in the main text. For simplicity, we set Ω = 1.

Example: the triangular lattice
The triangular lattice is a two-dimensional Bravais lattice with primitive lattice vectors a (1, 0) and a 2 = a cos π 3 , sin π 3 , with a the real-space lattice spacing. In the Hilbert space, we have again a triangular structure where a new site is added on each link. It is not difficult to see that this reduces to a pure triangular lattice by choosing a basis of 4 sites, a single one-excitation one (n 1 = 1) and 3 pair ones (n 2 = 3). The primitive lattice vectors will be the same as above, where we fix for simplicity a = 1. The basis states can be chosen according to:    The matrix C − k is now a 1 × 3 matrix whose elements can be computed via the procedure outlined above: C − k,1,1 : from basis state |ν 1 one can reach state |µ 1 within the same Bravais lattice site (⇒ +1) or state |µ 1 at the neighboring site j = a 1 (⇒ +e i k· a1 ).
Collecting all terms, the matrix C − k reads C − k = 1 + e i k· a1 , 1 + e i k· a2 , 1 + e i k·( a1− a2) ≡ w † k (85) and is equivalent to a three-dimensional vector w k . Thus, the total matrix M k can be expressed as |µ 1 : a one-excitation site at b = 0. We thus see that the C − k are 2 × 3 matrices and that there is at least one flat zero-energy band. The matrix elements can be identified column by column as follows: |ν 1 : From |ν 1 one can jump to |µ 1 or to |µ 2 remaining in the same Bravais lattice site.
|ν 3 : From |ν 3 one can jump to |µ 1 in the same site or to |µ 2 changing site by j = − a 2 (⇒ e −i k· a2 ). Hence, with w k,1/2 three-dimensional vectors. The matrix M k is thus The kernel state is a five-dimensional vector (0, 0, v k ) which satisfies w † k,1/2 · v k = 0.
To identify the remaining non-zero bands, we again take the square of the total matrix M k : where the first block is 2 × 2 and the second one 3 × 3. We can now diagonalize the first block to find (see Fig. 1 in the main text) with λ k,± ≥ 0. The four non-trivial bands will thus correspond to ± λ k,+ and ± λ k,− . Working out the scalar products w k,1 and w † k,2 · w k,1 = 1 + e i k·( a1− a2) + e −i k· a2 we obtain by substitution λ k,± = 3 ± 1 + e i k·( a1− a2) + e −i k· a2 = 3 ± 3 + 2 cos k · ( a 1 − a 2 ) + 2 cos k · a 2 + 2 cos k · a 1 .
From the first equality we see that the second addend is always ≤ 3. It is 3 only when k = 0 (up to reciprocal lattice translations G, see (64)). Hence, λ − ( k = 0) = 0 is a minimum and λ + ( k = 0) = 6 is a maximum. The bands ± λ k,− touch at k = 0 with linear dispersion.
Second, the argument of the absolute value will vanish when k · ( a 1 − a 2 ) = ± 2π 3 + 2πn , − k · a 2 = ± 4π 3 + 2πm (100) where the signs must be chosen consistently. Up to reciprocal lattice translations, one can choose identifying the points at the vertices of the hexagonal first Brillouin zone (one can verify this point lies at the boundary of two of the conditions in (90)). Therefore, the two upper bands λ k,+ and λ k,− touch at the vertices of the first Brilluoin zone with linear dispersion and similarly do the lower bands − λ k,+ and − λ k,− .

Localization and scaling exponents
In Table I we list the scaling exponents ν i extracted from the localization lengths ξ i ∼ s νi , i = 1, 2 at a given disorder strength s as described in Fig. 3 of the main text. For comparison, and further to the discussion in the main text, we list in the second and third column scaling exponents obtained with flat disorder distribution, where the disorder energies δV affecting the sites of the (synthetic) Lieb ladder are drawn from a uniform interval [−W/2, W/2]. The second (third) column corresponds to a situation, where only the sites corresponding to the pair states (all sites of the Lieb ladder) are affected. Finally, we list in the last column the values presented in [6]. It is apparent from the Table I that the scaling exponents at energies = 1, 1.8 and √ 6 show a reasonable agreement corresponding to the "out" and "in" Anderson scalings 0 and 2. For = √ 2, we observe for the "in" scaling a different behaviour between the cases listed in Table  1. The anomalous value of 4/3 appears only when the disorder acts on all sites. However, when the disorder affects only the pair sites it corresponds to the Anderson value of 2 (both when the disorder is flat and drawn from the distribution (13)). In contrast, we remark that for = √ 2, for the flat distribution the result is independent on whether it acts on all or only on pair state sites. On the other hand, the values we obtain for the disorder distribution drawn from (13), i.e. ν( = √ 2) ≈ {1.1, 1.1} doesn't seem to be close to either the anomalous or the edge scaling exponents and, based solely on the present analysis, cannot be simply attributed to the disorder acting on only the pair state sites.

Initial state preparation and evolution
In this section we consider the preparation of the state ) localized at rungs i, i + 1 of the ladder. We assume that each atom in the ladder can be addressed individually with a laser pulse of Rabi frequency Ω R and duration τ so that the atomic spin evolves according to if the laser detuning ∆ is set such that it is resonant with the transition of the addressed atom. If ∆ Ω R , instead, it acts trivially like an identity operator. Specifically, we will distinguish two special cases, namely ∆ = 0 in addition to ∆ = −V (R 0 ) corresponding to the blockade and facilitation condition respectively. The state |ψ loc can be obtained by application of six pulses on initially all atoms in the spin-down state as |ψ loc = F 2 (2π)F 4 (2π)F 3 (π)F 2 ( π 2 )B 4 (π)B 1 ( π 2 ) |ψ ↓..↓ , where B j (θ), F j (θ) stand for the laser pulse of area θ = Ω R τ in the blockaded (B) and facilitated (F) regime applied at site j = 1, .., 4 labeling the effective plaquette formed by the four sites corresponding to the i-th and (i+1)-th rung of the ladder, see Eq. (102). Here, the first pulse creates an excitation at site 1, the second pulse then exploits the blockade mechanism to create a superposition of spin-up states at sites 1 and 4. Next, the pulse in the facilitated regime applied at site 2 creates a superposition of the form −i |↑ +|↓ if and only if a single nearestneighbor is already excited, and so forth. We have omitted the global −i factors in the second, and fourth lines of (102). In practice the choice of Ω R is a trade-off between the need to keep the state-preparation time to a minimum (implying higher values of Ω R ) and the upper bounds imposed for keeping the blockade and facilitation conditions preserved, see [7] for details of these issues.
Once the state |ψ loc has been prepared, it evolves according to H eff , Eq. (106). We define the probability of excitation at rung i is obtained as p α i = n α i / L i=1 n α i . Here n α i = ψ(t)|n α i |ψ(t) , α = u, l for the upper and lower leg of the ladder respectively. We then define the average positionx and the standard deviation ∆x of the excitations as  I. Scaling exponents ν for different energies obtained from fitting the behaviour of the localization lengths ξ. ξi ∼ s ν for the second and ξi ∼ W ν for the third and fourth columns, see text for details. The range of s and W in the first row denote the interval of the disorder parameter over which the fit was performed. Values in the second column obtained for α = 3 and N = 10 6 . N = 10 6 and 10 5 has been used in the third and fourth column.

Numerical simulation of the spin dynamics
We consider the full Hamiltonian Eq. ., E j and d Ξj is a shorthand for the spin separation in the given configuration Ξ j . We note that since configurations C, D correspond to single spin excitation, the associated disorder is vanishing by definition, δ Cj =δ Dj = 0, ∀j. The disorder energiesδ Ξj are generated from first drawing a specific realization of atomic positions at each site of the lattice in all three spatial directions with isotropic Gaussian distribution of width s. We then exactly evolve an initial state as |ψ(t) = exp −itĤ eff |ψ 0 , where b j are the elements of the basis (107) [strictly speaking there are only 5L − 2 non-trivial elements due to the OBC]. We note that the result of the evolution depends on two independent parameters, s and the ratio V (R 0 )/Ω, where the Rabi frequency should further satisfy Ω V (2R 0 ) for the effective Hamiltonian (106) to be valid. In Fig. 5a we present the results of the simulation analogous to that performed in Fig. 4, showing ∆x in the s − Ω/V (R 0 ) plane. We observe that the maximum of ∆x as a function of the disorder gets shifted towards higher disorder strength as Ω is increased.
The dependence of ∆x in Fig. 5a can be intuitively understood as follows. Smaller values of Ω/V (R 0 ) correspond to larger diagonal disorder elementsδ. Since it is the disorder which couples the flat and dispersive bands, the smaller the s, the smaller the Ω that is sufficient to cause the excitation hopping and thus the increase in ∆x. As s is increased, Anderson localization becomes more and more relevant and, correspondingly, the localization length at = 0 shrinks. Eventually, the state becomes capable of propagating over distances comparable to the localization length. Further increasing s then reduces this scale, corresponding to the decrease in ∆x. Clearly, by increasing Ω/V (R 0 ) the hopping amplitude becomes more relevant with respect to the typical energy shifts and the localization length is thus increased. Higher values of s are then required to localize the state again. In Fig. 5b,c we show a comparison between the exact evolution according to the full Hamiltonian (39), dashed line, and H eff , solid line. As expected, the predictions of the two models show an agreement in the regime where Ω V (2R 0 ), Fig.5c (V (R 0 )/Ω = 200). On the other hand for larger Ω, the two models start to differ as shown in Fig. 5b (V (R 0 )/Ω = 20).