# Signatures of the Many-body Localized Regime in Two Dimensions

###### Abstract

Lessons from Anderson localization highlight the importance of dimensionality of real space for localization due to disorder. More recently, studies of many-body localization have focussed on the phenomena in one dimension using techniques of exact diagonalization and tensor networks. On the other hand experiments in two dimensions have provided concrete results going beyond the previously numerically accessible limits while posing several challenging questions. We present the first large-scale numerical examination of a disordered Bose-Hubbard model in two dimensions realized in cold atoms, which shows entanglement based signatures of many-body localization. By generalizing a low-depth quantum circuit to two dimensions we approximate eigenstates in the experimental parameter regimes for large systems, which is beyond the scope of exact diagonalization. A careful analysis of the eigenstate entanglement structure provides an indication of the putative phase transition marked by the peak in the fluctuations of entanglement entropy in a parameter range potentially consistent with experiments.

Many-body localization (MBL) is a paradigm shift in out-of-equilibrium quantum matter. This novel phase of matter is characterized by the absence of thermalization Anderson (1958); Basko et al. (2006); Gornyi et al. (2005); Pal and Huse (2010); Oganesyan and Huse (2007); Nandkishore and Huse (2015). Thus, unitary quantum dynamics influences only a few degrees of freedom even in an interacting system while it can retain quantum coherence at high energies Bardarson et al. (2012). By localizing the excitations MBL can also protect certain forms of topological and symmetry-breaking orders in excited states and provides an opportunity to process quantum information in a system driven far from equilibrium Huse et al. (2013); Pekker et al. (2014); Bahri et al. (2015); Chandran et al. (2014). The quantum phase transition separating the MBL and thermal phases poses a major challenge for developing a theory of dynamical critical phenomena described by many-body entanglement in highly excited states Kjäll et al. (2014); Vosk et al. (2015); Potter et al. (2015); Grover (2014); Chandran et al. (2015a); Lim and Sheng (2016); Khemani et al. (2017); Dumitrescu et al. (2017).

It is well-known that dimensionality of real space affects single particle Anderson localization where in one and two dimensions (without spin-orbit coupling and broken time-reversal symmetry), the entire spectrum of single particle eigenstates is localized for arbitrary weak disorder Lee and Ramakrishnan (1985). Although the defining properties of MBL in one dimension are firmly established both theoretically and experimentally Imbrie (2016); Schreiber et al. (2015), the existence of the phenomenon in two and higher dimensions is much debated Chandran et al. (2016); De Roeck and Huveneers (2017); De Roeck and Imbrie (2017); Agarwal et al. (2017); Luitz et al. (2017); Ponte et al. (2017). Experiments in cold atoms measuring local and global equilibration have shown the persistence of quantum memory for long times providing indications of an MBL phase in two dimensions Choi et al. (2016); Bordia et al. (2017). On the other hand, theoretical criteria suggest that the lifetime of local memory is finite, albeit extremely long Chandran et al. (2016); De Roeck and Huveneers (2017).

In this article we evaluate the eigenstates of bosons hopping in a disordered lattice in two dimensions with onsite interactions. We generalize the tensor network method developed earlier for one dimensional systems to approximate the eigenstates in the localized regime Wahl et al. (2017). The system sizes accessible by our method are much larger than those that can be currently achieved by exact diagonalization or quantum Monte Carlo Geraedts et al. (2016); Inglis and Pollet (2016). We study the model in parameter regimes directly relevant to the experiments. The approximate eigenstates evaluated by our method exhibit signatures of a phase with low entanglement. By studying the distribution of entanglement entropy of a single site, we probe the eigenstates in the localized regime at large disorder. On lowering the disorder strength the distribution of entanglement becomes bimodal with increasing weight at larger entanglement, indicating the transition into the thermal phase. We find a critical disorder strength which is reasonably close to the experimentally measured value. This is the first theoretical demonstration in a large-scale system of a transition between eigenstates with distinct statistical properties of entanglement in two dimensions.

Quantum memory in MBL

In one dimension, fully many-body localized (FMBL) systems are characterized by an extensive number of quasi-local operators, known as l-bits () Serbyn et al. (2013); Huse et al. (2014); Ros et al. (2015); Imbrie et al. (2017). These operators commute exactly with the Hamiltonian and also with each other. Therefore, all energy eigenstates are represented as a string of eigenvalues of the l-bits. Locality of the l-bits forces the eigenstates to satisfy the area law of entanglement Bauer and Nayak (2013). This feature greatly simplifies the spectral properties of an FMBL Hamiltonian. As a result, the unitary transformation () which diagonalizes the Hamiltonian can be efficiently approximated by a spectral tensor network (TN) Pekker and Clark (2017); Chandran et al. (2015b); Pollmann et al. (2016) where all eigenstates can be represented as matrix product states Friesdorf et al. (2015); Yu et al. (2017); Khemani et al. (2016). Throughout this work, we consider shallow quantum circuits, which are TNs which can be contracted efficiently in any dimension. They are composed of layers of local unitaries with each unitary acting on a finite contiguous block of spins Wahl et al. (2017). For the optimal set of local unitaries, the TN transforms the Hamiltonian into a predominantly diagonal form. The magnitude of the remaining off-diagonal matrix elements provides an estimate of the accuracy of the TN approximation. This method has been successfully implemented to evaluate the eigenstates of large systems in the fully many-body localized phase and estimate the location of the phase transition in one dimension Wahl et al. (2017). This method can also be used to construct the l-bits close to the MBL transition Kulshreshtha et al. (2017). Although exact diagonalization gives the entire spectrum of eigenstates, the construction of l-bits is non-trivial and the low-depth quantum circuit provides a systematic and efficient construction of the full set of operators.

In dimensions greater than one, the l-bit phenomenology can break down due to resonant interactions which destroy the exact conservation of the set of quasi-local operators Chandran et al. (2016); De Roeck and Huveneers (2017); De Roeck and Imbrie (2017). As a result, for a finite-size system thermalization can be restored at late times. However, if local operators remain approximately conserved, the system displays features analogous to MBL on experimentally relevant time scales. Unlike in one dimension where the l-bits are exactly conserved due to the suppression of resonances even for a finite system size, in higher dimensions the proliferation of resonances is believed to be exponentially sensitive to perturbations. Therefore, MBL giving rise to a stable quantum phase may only be well-defined as the system size tends to infinity, although there are arguments in favor of instability of MBL even in this asymptotic limit De Roeck and Huveneers (2017); De Roeck and Imbrie (2017). We use the semantics of quantum phase transitions to describe the features in our finite-size calculations, even though our numerical investigation may not capture this asymptotic limit. However, the error in our approximation provides a lower bound for the time scale over which local observables will decay in this system. Thus, our method approximates well the short-time dynamics. Hence, while not answering the question of the existence of exact l-bits in two dimensions, we expect our approach to describe the experimentally observed localization phenomena.

Note that even though the experimentally accessible time scales and system sizes are finite, there are still no theoretical descriptions of the observed phase diagram. Indeed, the conditons required to confirm the breakdown of conservation of these quasi-local operators likely require measurements for astronomically long times for a physical choice of parameters Choi et al. (2016). Nonetheless, our work sheds light in this experimentally relevant regime using our novel numerical technique.

In this article we evaluate eigenstates with an accuracy which grows with the range of the local unitaries. By studying the structure of entanglement in these eigenstates, we examine the system’s ability to serve as long-lived quantum memory over a range of disorder strengths which is consistent with state of the art experiments on MBL in two dimensions Choi et al. (2016); Bloch (). The low-depth quantum circuit captures the entanglement in local ensembles, but does not give the long-range entanglement of thermal systems.

Tensor network ansatz

We propose a shallow TN capable of approximating the two-dimensional localized phase. It encodes a unitary which approximately diagonalizes the localized Hamiltonian . We consider an lattice with periodic boundary conditions. The tensor network is composed of two layers of smaller unitaries acting on sites ( even) as shown in Fig. 1a. It is a natural extension of the ansatz in Ref. Wahl et al., 2017.
The approximately conserved operators are given by ( is the spin- operator at site ). These operators commute mutually by construction. The approximates , which is itself expected to be only approximately conserved, but acts as a local memory on the experimentally relevant time scales. Hence, the tensor network can be optimized by minimizing the (squared) trace norm of the commutator with the Hamiltonian,

(1) |

corresponds to the tensor network contraction shown in Fig. 1b, which is non-trivial only in the region covered by the unitaries , . As a result, for a nearest-neighbor Hamiltonian, Eq. (1) decomposes into local parts,

(2) |

The explicit representation of by tensor network contractions can be found in Appendix A. The overall computational cost scales as . In the general case, the length of the region within which is non-trivial is of order . Therefore, in the MBL regime, the error in representing is expected to decrease exponentially with . The same will be true for the error of expectation values of local observables. Hence, for fixed system size, this error is an inverse polynomial of the computational cost. The localization of the eigenstates in the MBL regime Chandran et al. (2016) implies that this error does not significantly increase with system size, making our approximation efficient. Note that even for , our approach can incorporate effects which are far beyond the reach of exact diagonalization: In that case, the operators are non-trivial within a region of . That implies non-vanishing correlations between a given site and other sites within a regime. Due to the strict locality of the operators, our approximation is not affected by the boundary conditions, i.e., it describes the MBL regime directly in the bulk.

MBL phase in 2D

In the following, we use our TN approach to investigate the MBL phase found experimentally Choi et al. (2016) in a two-dimensional optical lattice. The phase transition was determined by measuring the loss of memory of the initial conditions of bosonic Rubidium-87 atoms in a random optical potential. The dynamics of the atoms can be modeled by the disordered Bose-Hubbard Hamiltonian

(3) |

where the nearest neighbor links are counted once for each pair , () are bosonic annihilation (creation) operators and is the particle number operator. is random variable chosen from a Gaussian distribution with full-width half-maximum , on-site interaction and we fix the energy scale through . Note that we neglect the trapping potential, as we are interested in the phase transition in the bulk of the system, i.e., the center of the trap. The local bosonic Hilbert space is a priori infinite dimensional. We truncate the on-site occupation number to , which is an approximation suitable to describe the experiments where the higher occupation number states are rarely observed. (Experimentally, only up to 7 % of sites are occupied with doublons, and no triplons are observed Bloch ().) We perform calculations for (corresponding to hard core bosons or spin- particles) and (corresponding to spin- particles). As the Hamiltonian is particle number conserving and real, the unitaries can be parameterized as the exponentials of antisymmetric matrices, which also conserve particle number Wahl et al. (2017) (for details, see Appendix B). We verified numerically that no significant loss of accuracy occurs compared to a full parameterization.

A hallmark of the localized phase is the locality of entanglement even in the presence of interactions. A single spin is only weakly coupled to its environment and its reduced density matrix in an eigenstate (where represents the local spin and is the rest of the system), remains close to a pure state. On the other hand in the thermal phase, the single spin is maximally entangled with its environment. For a two level system (spin-), this maximal entanglement is equal to and is for a three-level system (spin-). We will first focus on the results. These come closest to the experiments Choi et al. (2016); Bloch (), which were carried out at an average filling . Since the overall boson number is conserved, the experimentally probed transition refers only to half filling. In contrast, our calculations describe a transition occurring at all occupation numbers (given the on-site truncation). For , the majority of simulated eigenstates are close to half-filled, whereas for , most eigenstates have on average approximately one boson per site.

The approximate eigenstates represented by the TN have a local structure which is reliable as the approximate integrals of motion given by the TN have a small commutator with the Hamiltonian. The degree of the commutation increases with increasing disorder strength (see Appendix E, Fig. 8b).

Signatures of the MBL transition

We calculated the on-site entanglement entropies () for 30 disorder realizations on a lattice as a function of disorder . Those entropies can be computed efficiently as explained in Appendix F. Their distribution displayed in Fig. 2 shows a sharp peak at . The distribution has a tail to large values of entanglement which decays quickly to zero.
The bimodal entanglement entropy distributions with a sharp peak at along with a broad peak at lower values of for provide evidence for the phase transition from the MBL to the thermal regime. The bimodality indicates the co-existence of highly localized states with low entanglement and delocalized states with larger entanglement in the critical regime. This resembles the behavior around the critical point of one-dimensional models exhibiting MBL Yu et al. (2016). In contrast, for , the distribution is mostly concentrated in a single peak at , where delocalized states are absent. Note that the peak at for appears because the map from the reduced density matrix to the entropy has a maximum at ( is diagonal due to symmetry.). Therefore, transforming the probability distribution of to produces a singularity due to the derivative of tending to zero at the maximum. As expected there is no peak in the histogram at (not shown).

We also computed the reduced density matrices of all plaquettes (using the contraction scheme presented in Appendix F) and extracted their entanglement spectra. The entanglement spectra are given by the eigenvalues of the entanglement Hamiltonian defined via . The corresponding distributions are shown in Fig. 2c,d. They have a double-peak structure throughout with a narrow maximum at small entanglement energies and a broad maximum at higher entanglement energies . The former moves towards zero and the latter towards infinite entanglement energies with increasing . The behavior of the narrow maximum is consistent with findings using exact diagonalization in one dimension Geraedts et al. (2017). However, in that reference, the broad maximum is observed only in the MBL phase and close to the transition. It is also found to move to larger entanglement energies with increasing disorder strength. It is thought to be a remnant of the critical regime even deep in the localized phase. We obtain a broad maximum for small , too, presumably because the optimized tensor network does not fully capture the delocalized phase, and instead displays critical features there also.

In order to locate the critical point more precisely, we evaluated the variance of the on-site entanglement entropies with respect to different disorder realizations. Finally, we averaged over the site positions to improve smoothness Wahl et al. (2017). The results for the entanglement entropy fluctuation are shown in Fig. 3. For , it peaks at which is potentially consistent with the latest experimentally measured transition point Bloch (). Note that in Ref. Choi et al., 2016, the MBL-to-thermal transition was measured at for ; however there are several experimental conditions which can lead to a lower estimate of the measured critical point which continue to be under intense investigation.

As already mentioned above, for , we do not obtain the same transition point as in the experiments, which were not carried out at average filling 1. The entropy fluctuations using 30 disorder realizations and lattices are also shown in Fig. 3. They indicate MBL-to-thermal transitions at and , respectively. Future experiments initialized with a charge density wave of doublons might verify this prediction. The trend towards lower transition points with decreasing interaction strength is consistent with Anderson localization at for

Conclusion

We have presented the first large-scale numerical study of eigenstates in the MBL regime of the disordered Bose-Hubbard model in two dimensions using shallow tensor networks (quantum circuits). The novel variational ansatz approximates the full spectrum of eigenstates as area-law entangled in this regime.
By characterizing the statistics of entanglement as a function of disorder, the location of the transition into the thermal phase has been estimated. The distribution of the entanglement spectrum of these eigenstates exhibits the characteristics of MBL which further supports our observations.

Our work has important consequences for experiments with ultra-cold atoms studying quench dynamics in the presence of disorder. It provides an estimate of the MBL transition point in two dimensions as a function of interaction strength, which can be verified in these experiments. The current experiments can access the transition using local dynamical properties of local observables in a finite system, which we expect to be well-described by our approximate eigenstates. This also opens the door to study the phenomenology of MBL in two dimensions, where the time-tested method of exact diagonalization is computationally prohibitively expensive. Large-scale simulations with higher will further elucidate the nature of the MBL regime in two dimensions.

Acknowledgments

We are grateful to Anushya Chandran and Chris Laumann for stimulating discussions and a careful reading of the manuscript, and to Antonio Abadal, Immanuel Bloch, Jae-Yoon Choi, and Christian Gross for detailed discussions related to the experiments. We would also like to thank David Huse, John Imbrie, Vadim Oganesyan, Wojciech De Roeck, Antonello Scardicchio, Shivaji Sondhi, and Thomas Spencer for fruitful discussions. S.H.S. and T.B.W. are both supported by TOPNES, EPSRC grant number EP/I031014/1. S.H.S. is also supported by EPSRC grant EP/N01930X/1. T.B.W. acknowledges usage of the University of Oxford Advanced Research Computing (ARC) facility in carrying out this work, http://dx.doi.org/10.5281/zenodo.22558. A.P. is supported by the Glasstone Fellowship and would like to thank the Aspen Center for Physics, which is supported by National Science Foundation grant PHY-1607611, and Simons Center for Geometry and Physics, Stonybrook University for their hospitality, where part of this work was performed. T.B.W. is grateful for support by the European Commission under the Marie Curie Programme. The contents of this article reflect only the authors’ views and not the views of the European Commission. Statement of compliance with EPSRC policy framework on research data: This publication is theoretical work that does not require supporting research data.

## Appendix A Decomposing the commutator norm into local parts

Here, we describe the contraction scheme to evaluate the commutator norm as defined in Eq. (1). We assume that the Hamiltonian is a sum of nearest neighbor terms, . In that case, the two terms in the final expression of Eq. (1) are

(4) |

and

(5) |

respectively. Using the representation for shown in Fig. 1b (which is non-trivial only in the region covered by the unitaries , ), the right side of Eq. (5) can be evaluated as shown in Fig. 4. The right side of Eq. (4) can be written as a similar tensor network contraction (this is most obvious if one inserts between and leading to an expression which is formally equivalent to Eq. (5)). Therefore, can be expressed as in Eq. (2).

The tensor network contraction of Fig. 4 is carried out most efficiently as follows. First, one combines unitaries and their adjoints (shown in the same color in Fig. 4) to new tensors. If there is no other tensor between them, this obviously results in the identity. In contrast, will be combined with and to form a new tensor . The other tensors might have to be combined with the Hamiltonian terms and , respectively. The resulting tensor network is shown in Fig. 5. A special case is the one of with the two Hamiltonian terms being connected by a separate line if they do not lie entirely within the region covered by the unitaries . This squares the dimension of the index contraction corresponding to one of the lines in Fig. 5. The leading contribution to the computational cost stems from the multiplication of matrices, i.e., it scales as . As the number of tensor network contractions of the type shown in Fig. 5 is independent of the system size, each term contributes to the overall computational cost, which is thus of order .

The derivative of our figure of merit can be calculated easily by cutting out the unitary matrix the derivative is taken of in Figs. 4 and 5, respectively, and contracting the tensors in such an order that the removed unitary would come last, cf. Ref. Wahl et al., 2017. This comes at a subleading increase of computational cost.

## Appendix B Optimization method

We take advantage of the symmetries of the considered model in order to reduce the number of parameters: The Hamiltonian is real and possesses symmetry, i.e., particle number conservation in the experimental model Eq. (3) and spin- conservation in its spin representations and the anti-ferromagnetic Heisenberg model, . In general, such global symmetries can be imposed on the individual tensors which make up the tensor network without much loss of accuracy Sanz et al. (2009); Pérez-García et al. (2010). We do so by taking real unitary matrices , i.e., orthogonal matrices, which on top of that preserve the spin- component individually, . Both taken together imply that the unitaries must have a block diagonal form (in the spin- basis),

(6) |

where , a skew-symmetric matrix. In each minimization run, the unitaries are initialized as identities and optimized using a quasi-Newtonian routine supplied by the derivative with respect to the parameters. We optimize one unitary at a time and sweep alternately over all the columns and all the rows of the system (starting each sweep at site ).

## Appendix C Comparison with exact diagonalization

We benchmark our approach by comparing with exact diagonalization for the hard-core Bose-Hubbard model, Eq. (1), on a lattice with open boundary conditions for . To that end, we place additional unitaries in the top layer, which pairwise connect the ones of the bottom layer as shown in Fig. 6. The open boundary conditions are meant to prevent artificial delocalization (as compared to periodic boundary conditions) for such a small system. In Fig. 7, we show the energy errors and overlaps of the approximate compared to the exact eigenstates for and 30 disorder realizations. More than 99.9 % of all approximate eigenstates have at least 50 % overlap with the actual ones, which is a very high accuracy given that the dimension of the Hilbert space is . Because of the local structure of the tensor network, it follows that the true eigenstates are relatively well localized for this disorder strength. If they stay well localized for larger system sizes, the tensor network approximation remains accurate and is able to capture local observables with a constant accuracy as the system size is increased. The average optimized figure of merit was . We also calculated the variance averaged over eigenstates,

(7) |

Its mean over disorder realizations was .

## Appendix D Increasing the system size

As the figure of merit decomposes into a sum of local parts, the expectation is it increases linearly with the system size on average. We confirmed this by studying 30 disorder realizations at for as a function of system size for periodic boundary conditions, see Fig. 8a.

In the MBL regime, eigenstates can be described by effective spin degrees of freedom Chandran et al. (2016) , which have exponentially decaying non-trivial matrix elements as a function of distance and commute with the Hamiltonian up to exponentially small corrections in the system size. Therefore, for sufficiently large systems, we gain the following picture: All effective spin degrees of freedom which are located at a distance from a given region contribute to local properties of that region by an amount of order . Hence, their contribution vanishes exponentially with distance and the error made by approximating the eigenstates in terms of strictly short-range tensor network -bits is independent of the system size. In other words, local observables can be approximated with a constant accuracy using our TN if the system size is increased. Note that at the same time the computational cost grows linearly with system size.

## Appendix E Accuracy as a function of disorder strength

We also investigated the accuracy of our approximation as a function of for 30 disorder realizations. For different , only the overall strength of the random magnetic fields was adjusted to enhance comparability. The optimized figure of merit is shown in Fig. 8b. We find a strong, monotonic increase of the error with decreasing , especially below . Notwithstanding, a concrete picture of the phase transition can only be gained from the distribution of entanglement entropies.

## Appendix F Calculation of the entanglement entropies and spectra of blocks of sites

In the following, we explain how the reduced density matrices of plaquettes of approximate eigenstates given by our tensor network can be calculated efficiently. The computation of single-site reduced density matrices is a special case thereof. We will focus on the case of plaquettes connecting four unitaries of the upper layer, which is the hardest one to evaluate. Other positions of the blocks can be treated in a simplified way.

A reduced density matrix of a given approximate eigenstate can be graphically represented by taking the tensor network in Fig. 1a, fixing its lower indices according to the -bit configuration of that approximate eigenstate and contracting it with its adjoint from above, leaving only the legs in the chosen region open. After cancelling unitaries with their Hermitian conjugates, one obtains the tensor network contraction shown in Fig. 9a (for ). Even though this tensor network is non-trivial in a Hilbert space of dimension ( Hilbert space dimension of a single site), its tensors can be contracted in such an order that the computational time is only of order , cf. Fig. 9b.

Finally, in Fig. 10 we present the data for the distribution of entropies for (after tracing out all but a single site) along with the histogram over the entanglement spectra for plaquettes. The data is based on 30 disorder realizations on a lattice and samples over 10,000 eigenstates per site.

## References

- Anderson (1958) P. Anderson, Phys. Rev. 109, 1492 (1958).
- Basko et al. (2006) D. Basko, I. Aleiner, and B. Altshuler, Annals of physics 321, 1126 (2006).
- Gornyi et al. (2005) I. Gornyi, A. Mirlin, and D. Polyakov, Phys. Rev. Lett. 95, 206603 (2005).
- Pal and Huse (2010) A. Pal and D. A. Huse, Phys. Rev. B 82, 174411 (2010).
- Oganesyan and Huse (2007) V. Oganesyan and D. Huse, Phys. Rev. B 75, 155111 (2007).
- Nandkishore and Huse (2015) R. Nandkishore and D. A. Huse, Annual Review of Condensed Matter Physics 6, 15 (2015).
- Bardarson et al. (2012) J. H. Bardarson, F. Pollmann, and J. E. Moore, Phys. Rev. Lett. 109, 017202 (2012).
- Huse et al. (2013) D. A. Huse, R. Nandkishore, V. Oganesyan, A. Pal, and S. L. Sondhi, Phys. Rev. B 88, 014206 (2013).
- Pekker et al. (2014) D. Pekker, G. Refael, E. Altman, E. Demler, and V. Oganesyan, Phys. Rev. X 4, 011052 (2014).
- Bahri et al. (2015) Y. Bahri, R. Vosk, E. Altman, and A. Vishwanath, Nature communications 6 (2015).
- Chandran et al. (2014) A. Chandran, V. Khemani, C. R. Laumann, and S. L. Sondhi, Phys. Rev. B 89, 144201 (2014).
- Kjäll et al. (2014) J. A. Kjäll, J. H. Bardarson, and F. Pollmann, Phys. Rev. Lett. 113, 107204 (2014).
- Vosk et al. (2015) R. Vosk, D. A. Huse, and E. Altman, Phys. Rev. X 5, 031032 (2015).
- Potter et al. (2015) A. C. Potter, R. Vasseur, and S. A. Parameswaran, Phys. Rev. X 5, 031033 (2015).
- Grover (2014) T. Grover, arXiv preprint arXiv:1405.1471 (2014).
- Chandran et al. (2015a) A. Chandran, C. R. Laumann, and V. Oganesyan, arXiv:1509.04285 (2015a).
- Lim and Sheng (2016) S. P. Lim and D. N. Sheng, Phys. Rev. B 94, 045111 (2016).
- Khemani et al. (2017) V. Khemani, S. P. Lim, D. N. Sheng, and D. A. Huse, Phys. Rev. X 7, 021013 (2017).
- Dumitrescu et al. (2017) P. T. Dumitrescu, R. Vasseur, and A. C. Potter, Phys. Rev. Lett. 119, 110604 (2017).
- Lee and Ramakrishnan (1985) P. A. Lee and T. V. Ramakrishnan, Rev. Mod. Phys. 57, 287 (1985).
- Imbrie (2016) J. Z. Imbrie, Journal of Statistical Physics 163, 998 (2016).
- Schreiber et al. (2015) M. Schreiber, S. S. Hodgman, P. Bordia, H. P. Lüschen, M. H. Fischer, R. Vosk, E. Altman, U. Schneider, and I. Bloch, Science 349, 842 (2015).
- Chandran et al. (2016) A. Chandran, A. Pal, C. R. Laumann, and A. Scardicchio, Phys. Rev. B 94, 144203 (2016).
- De Roeck and Huveneers (2017) W. De Roeck and F. m. c. Huveneers, Phys. Rev. B 95, 155129 (2017).
- De Roeck and Imbrie (2017) W. De Roeck and J. Z. Imbrie, arXiv preprint arXiv:1705.00756 (2017).
- Agarwal et al. (2017) K. Agarwal, E. Altman, E. Demler, S. Gopalakrishnan, D. A. Huse, and M. Knap, Annalen der Physik 529, 1600326 (2017).
- Luitz et al. (2017) D. J. Luitz, F. Huveneers, and W. De Roeck, arXiv preprint arXiv:1705.10807 (2017).
- Ponte et al. (2017) P. Ponte, C. Laumann, D. A. Huse, and A. Chandran, arXiv preprint arXiv:1707.00004 (2017).
- Choi et al. (2016) J.-y. Choi, S. Hild, J. Zeiher, P. Schauß, A. Rubio-Abadal, T. Yefsah, V. Khemani, D. A. Huse, I. Bloch, and C. Gross, Science 352, 1547 (2016).
- Bordia et al. (2017) P. Bordia, H. Lüschen, S. Scherg, S. Gopalakrishnan, M. Knap, U. Schneider, and I. Bloch, arXiv preprint arXiv:1704.03063 (2017).
- Wahl et al. (2017) T. B. Wahl, A. Pal, and S. H. Simon, Phys. Rev. X 7, 021018 (2017).
- Geraedts et al. (2016) S. D. Geraedts, R. Nandkishore, and N. Regnault, Phys. Rev. B 93, 174202 (2016).
- Inglis and Pollet (2016) S. Inglis and L. Pollet, Phys. Rev. Lett. 117, 120402 (2016).
- Serbyn et al. (2013) M. Serbyn, Z. Papić, and D. A. Abanin, Phys. Rev. Lett. 111, 127201 (2013).
- Huse et al. (2014) D. A. Huse, R. Nandkishore, and V. Oganesyan, Phys. Rev. B 90, 174202 (2014).
- Ros et al. (2015) V. Ros, M. Mueller, and A. Scardicchio, Nuclear Physics B 891, 420 (2015).
- Imbrie et al. (2017) J. Z. Imbrie, V. Ros, and A. Scardicchio, Annalen der Physik 529, 1600278 (2017), 1600278.
- Bauer and Nayak (2013) B. Bauer and C. Nayak, Journal Of Statistical Mechanics-Theory And Experiment 2013, P09005 (2013).
- Pekker and Clark (2017) D. Pekker and B. K. Clark, Phys. Rev. B 95, 035116 (2017).
- Chandran et al. (2015b) A. Chandran, J. Carrasquilla, I. H. Kim, D. A. Abanin, and G. Vidal, Phys. Rev. B 92, 024201 (2015b).
- Pollmann et al. (2016) F. Pollmann, V. Khemani, J. I. Cirac, and S. L. Sondhi, Phys. Rev. B 94, 041116 (2016).
- Friesdorf et al. (2015) M. Friesdorf, A. H. Werner, W. Brown, V. B. Scholz, and J. Eisert, Phys. Rev. Lett. 114, 170505 (2015).
- Yu et al. (2017) X. Yu, D. Pekker, and B. K. Clark, Phys. Rev. Lett. 118, 017201 (2017).
- Khemani et al. (2016) V. Khemani, F. Pollmann, and S. L. Sondhi, Phys. Rev. Lett. 116, 247204 (2016).
- Kulshreshtha et al. (2017) A. Kulshreshtha, A. Pal, T. B. Wahl, and S. H. Simon, arXiv:1707.05362 (2017).
- (46) I. Bloch, Private Communication.
- Yu et al. (2016) X. Yu, D. J. Luitz, and B. K. Clark, Phys. Rev. B 94, 184202 (2016).
- Geraedts et al. (2017) S. D. Geraedts, N. Regnault, and R. N. Nandkishore, arXiv:1705.00631 (2017).
- Sanz et al. (2009) M. Sanz, M. M. Wolf, D. Pérez-García, and J. I. Cirac, Phys. Rev. A 79, 042308 (2009).
- Pérez-García et al. (2010) D. Pérez-García, M. Sanz, C. Gonzalez-Guillen, M. M. Wolf, and J. I. Cirac, New Journal of Physics 12, 025010 (2010).