Transition to a many-body localized regime in a two-dimensional disordered quantum dimer model

Many-body localization is a unique physical phenomenon driven by interactions and disorder for which a quantum system can evade thermalization. Although the existence of a many-body localized phase is now well established in one-dimensional systems, its fate in higher dimensions is an open question. We present evidence for the occurrence of a transition to a many-body localized regime in a two-dimensional quantum dimer model with interactions and disorder. Our analysis is based on the results of large-scale simulations for static and dynamical properties of a consequent number of observables. Our results pave the way for a generic understanding of occurrence of a many-body localization transition in a dimension larger than one and highlight the unusual quantum dynamics that can be present in constrained systems.


I. INTRODUCTION
The emergence of thermalization in isolated quantum systems, which are not in contact with a bath and solely follow unitary Hamiltonian dynamics, has long been a central issue in statistical mechanics [1][2][3].This question has been revived thanks to an extraordinary increase in the degree of control and isolation reached in cold-atom experiments [4,5].The central concept of the eigenstate thermalization hypothesis (ETH) [6,7], a conjecture on expectation values of local operators in the eigenstates of the isolated system, allows to provide a physical understanding and justification of this thermalization.The ETH is found to be valid for the generic majority of many-body systems [8,9].Recent years have witnessed a vast body of research concentrating on a situation where the ETH fails, namely, the phenomenon of many-body localization (MBL) [10,11], which occurs most noticeably in the presence of strong disorder.
In a many-body localized system, eigenstates at the same energy density do not form a statistical ensemble, and local expectation values strongly differ from one eigenstate to the other.MBL states sustain a low level of quantum entanglement in the form of an area law [12] (in contrast with ETH states which obey a volume law for entanglement entropy) and have been found to be able to sustain long-range order in cases where equilibrium states are forbidden to [13][14][15][16][17]. Dynamics in MBL systems markedly differ from those of thermalizing systems, with, e.g., long-time memory of the initial state, absence of transport for conserved quantities (localization), and a slow propagation of quantum information with logarithmic growth of entanglement [18,19].The latter example allows to distinguish MBL from Anderson localization [20] (where entanglement is bounded) which occurs in noninteracting systems.
The existence of many-body localization in one dimension is now established and reasonably well understood in one-dimensional (1D) quantum lattice systems with shortrange interactions, thanks to a concerted effort of approaches, including numerical simulations [21][22][23], phenomenological renormalization [24][25][26][27][28] approaches, rigorous results [29,30], and cold-atom experiments [31][32][33].Crucial to this understanding is the provable [29] existence of local integral of motions [34,35], denoted as l-bits (for localized bits): These emergent localized operators, which diagonalize the Hamiltonian of the system at strong disorder, explain most of the salient features of the MBL phase [36].In the typical setting of lattice models, the MBL phase, indeed, emerges in the presence of strong enough disorder and is separated from a ETH phase at low disorder by a many-body localization transition involving an exponential number of eigenstates [21].The reviews [37][38][39][40][41][42] provide an introduction and recent insights on this alive field of research.

A. Many-body localization in higher dimensions
An outstanding question which remains open is the fate of many-body localization in dimensions higher than one.A basic observation in this case is that localization is more difficult to achieve as there are more transport channels than in one dimension.The Anderson localization survives in two dimensions (2D) for all disorder strengths but with a much higher localization length and needs an increasingly stronger critical disorder in higher dimensions [43].
In the many-body case, analytical arguments were presented [44][45][46] that suggest the intrinsic instability of a MBL phase in dimensions higher than one.They are based on the inevitable existence of thermal "bubbles" (created by, e.g., local regions of space with low disorder) which grow in time and eventually thermalize any MBL sample.The thermalization conveyed in this scenario is asymptotic, meaning that it may require an asymptotically long time before signs of localization (for instance, memory of the initial state) disappear.In these works, the MBL region is represented in terms of the l-bits, the exponentially localized integral of motions which form a complete set of operators spanning the Hilbert space.
Despite these arguments, recent remarkable numerical simulations [47][48][49][50] analyze possible signatures of a MBL phase in two-dimensional lattice models of spins or bosons in either static [47] or dynamical quantities [48][49][50].Whereas results of Refs.[47][48][49] are in favor of the existence of a true thermodynamic MBL phase at finite disorder strength in 2D, Ref. [50] concluded oppositely (albeit evidencing slow dynamics even at moderate disorder strength).These results are very intriguing, but the methods used in these works are not unbiased: In Refs.[47][48][49], an approximate unitary diagonalization transformation based on the lowly entangled nature or localized nature of MBL eigenstates is used, which misses by construction the ETH phase (we note, nevertheless, that Ref. [47] presents evidence for the transition in 2D but without a finite-size or finite-entanglement scaling analysis).The finite-time dynamical results of Ref. [50] also rely on a method which is entanglement limited.Other exact diagonalization studies in 2D [51][52][53] are limited to small systems (lattices with, at most, 16 lattice sites) such that finite-size scaling is not possible.We also note a recent interesting investigation using exact numerics on a 2D model in the continuum compatible with MBL, albeit with a transition point shifting with the number of particles [54].A quantum Monte Carlo treatment of a Hamiltonian constructed from l-bits also argues for the existence of MBL in a 2D spin system [55].On the experimental side, recent advances in cold-atom experiments present clear dynamical signs of localization (in terms of memory of an initial imbalance of atoms) on the longest accessible timescale in two-dimensional systems either for bosonic atoms in random potential [56] or fermionic gases in quasiperiodic potentials [57].
In this paper, we present a comprehensive study of a two-dimensional strongly interacting disordered model, based on extensive numerical simulations on fairly large systems (square lattices with up to 64 sites).These simulations are exact, unbiased, and allow to study both ETH and MBL regimes and the transition between them.We study an important range of clusters of different sizes and shapes, allowing to follow the evolution of many ETH or MBL indicators with size.This is made possible thanks to two main specificities of our paper: We push large-scale parallel numerical methods to their limits (see Sec. II B), and we work on a constrained model (described in Sec.II A) that presents specific features appropriate for localization and numerical studies (Sec.I B).

B. Searching for MBL in 2D constrained models
Constraints exist or emerge in the description of a wide variety of physical systems, ranging from gauge theories, surface physics (adsorption of molecules [58]), glass physics [59][60][61], quantum or classical magnets [62,63], or atomic physics (such as with the Rydberg blockade [64]).Theoretical descriptions are written in terms of models with adapted degrees of freedom, such as dimer, loop, or ice models, which exhibit physical features not present in unconstrained counterparts.Consider, for instance, the dimer model (Fig. 1) that will be used in this paper where the degrees of freedom are dimers located on bonds of a given lattice (or graph).The constraint is given by the following rule: Each site of a lattice must be touched by one and only one dimer.Already, at the classical static level, this local rule enforces long-range algebraic correlations between dimers on bipartite lattices [65,66].Taking again the dimer example, it is not possible to obtain another dimer configuration by simply moving one dimer: Only constrained moves (the smallest of which is represented in Fig. 1) are possible.This can result in very complex dynamics (see Refs. [67][68][69] for classical examples).
Recent work highlighted the also unusual quantum dynamics exhibited by quantum constrained models either in theory [70][71][72][73][74] or in experiments [64] as well as for related lattice gauge theories [75][76][77][78][79][80][81][82].On one hand, we would, thus, expect that constrained models are good natural candidates to exhibit localization due to their enhanced slow dynamics.On the other hand, the local constraint is such that these models are already in a strongly interacting limit (even without interactions encoded in a Hamiltonian), which would rather favor delocalization and thermalization.This competition provides strong motivation to investigate the possibility of MBL in quantum constrained models.The interplay among constraints, interactions, and disorder has been explored in recent works on 1D models [83,84] (for Ref. [83], equivalent to a dimer model on a ladder) where constraints are relevant for the physics of cold-atomic chains in the regime of the Rydberg blockade [64], resulting in rich phase diagrams, with, e.g., different mechanisms to produce MBL [83].
Another important aspect (crucial for the numerical approach taken in this paper) worth mentioning is the fact that the constraints overall reduce the size of the phase space for allowed configurations.For instance, the size of the Hilbert space scales as ∼1.34 N for dimer coverings of square lattices with N sites, instead, of 2 N for spin-1/2 models (or 4 N if no constraints were present).This greatly helps numerics based on exact approaches as this enlarges the set of available finite samples in 2D.Given that exact diagonalization simulations have been instrumental in the study and understanding of MBL, this is also an important positive feature in the search for MBL in quantum constrained models.
The plan of the paper is as follows.We first describe the model and the 2D lattices used in this paper in Sec.II A along with a brief description of the exact numerical methods employed (Sec.II B).We then describe in detail physical properties of the eigenstates of this model in Sec.III, including spectral and eigenstate statistics, expectation values of local observables, as well as half-system entanglement entropy.Section IV is devoted to a machine-learning analysis

Potential term
Quantum dimer model with random potential  1).Plaquettes p of the lattice where a pair of parallel dimers is present (highlighted in green) contribute a potential V p (different for each plaquette) to the potential energy of this configuration.The kinetic term flips dimers on such plaquettes such that the orientation of the two dimers on this plaquette is changed.In the bottom right panel, the flip of the plaquette p = 20 (surrounded by red lines) is represented, resulting in a new configuration with a different potential energy.(c) On the square (bipartite) lattice with periodic boundary conditions, all dimer configurations can be classified by a pair of winding numbers (W x , W y ) defined as follows.Consider a horizontal (vertical) line-such as represented in the figure-perpendicular to lattice bonds.W x (W y ) is then equal to the algebraic number (i.e., ±1 where the sign is defined by the sublattice on a fixed side of the line) of dimers in this configuration crossing this line.The sign ±1 is represented by arrows on lattice bonds.Note that the winding numbers do not depend on the position of the line, and they cannot be changed by a local reconfiguration of dimers.The represented configurations all belong to the (W x , W y ) = (0, 0) sector in which all calculations are performed in this paper.Periodic boundary conditions are represented with the help of duplicated gray sites on the boundaries. of entanglement spectra.Section V then presents an analysis of dynamical properties after a quench.Finally, we summarize our results and present perspectives in Sec.VI.

A. Model
This paper focuses on a quantum dimer model (QDM) on a square lattice with random potential terms.Quantum dimer models have a long history in condensed matter as low-energy effective models.They have been first introduced by Rokhsar and Kivelson [85] in the context of the resonating valence bond theory of high-temperature superconductivity, an interest which has been revived with the use of QDMs to describe experiments on underdoped cuprates [86].QDMs are also prominent representations of highly frustrated magnets where dimers represent nearest-neighbor spin singlets [87].Dimer models are, indeed, often more tractable than frustrated spin models from which they can be systematically derived [88], while capturing the essential physics.For instance, the existence of 2D quantum spin liquids exhibiting topological order was first demonstrated in a QDM on a nonbipartite lattice [89].
Quantum dimer models are defined in a configuration space where only configurations which fulfill the dimer constraint are allowed [dimer coverings, see Fig. 1(a)].Their Hamiltonian contains a potential term that attributes an energy to every dimer covering and a kinetic term that flips dimers along a loop, allowing to pass from one dimer configuration to another.The specific QDM Hamiltonian that we consider in this paper is based on the square lattice and contains a potential term with a random component, taking the following form: (1) The sums run over all plaquettes p of the square lattice, V p is a random potential different for each plaquette which is drawn uniformly from a box distribution V p ∈ [−V, V ].The kinetic-energy scale is set to t = 1.The various terms are illustrated in Fig. 1(b).At V = 0 (no disorder), Ref. [90] showed that eigenstates in the middle of the spectrum obey the eigenstate thermalization hypothesis.Quite importantly, there is no limit where Eq. ( 1) can be mapped to a model of free particles moving in a random potential (no Anderson localization limit).
Quantum dimer models based on coverings of bipartite lattices (such as the square lattice) have an important property: They display conserved topological numbers, known as winding numbers, when considering systems with periodic boundary conditions.The definition of the two conserved winding numbers (denoted W x and W y ) is given in Fig. 1(c).It is easy to see that they are conserved by the Hamiltonian Eq. ( 1), and by any local moves (i.e., that do not go through periodic boundary conditions), hence, their topological nature.W x and W y both range from −L/2 to L/2 on a square lattice.The Hamiltonian is block diagonal with sectors identified by (W x , W y ).We note that the resulting Abelian U (1) symmetry of the model Eq. ( 1) is compatible with many-body localization [91].Without loss of generality, we will restrict to the largest sector (W x = 0, W y = 0) in what follows.

B. Lattices and numerical methods
This paper uses two-dimensional lattices (represented in Fig. 14 in the Appendix) with N sites and 2N bonds of the following form: regular square lattices with N = L 2 sites (with L = 4, 6, 8), rectangular lattices with N = L x × L y sites (with L x,y = 4, 6, 8, 10), as well as tilted square lattices with N = 32, 40, and 52 sites (see the Appendix for details).Periodic boundary conditions are taken.All the lattices used in this paper are presented in Table I, together with the corresponding Hilbert space sizes in the (W x = 0, W y = 0) sector.
We use different numerical methods to treat exactly the Hamiltonian Eq. ( 1).In Sec.III, we first use methods which allow to obtain exact eigenstates of H.As usually performed in the MBL context, we concentrate on eigenstates with eigenvalue E located right in the middle of the many-body spectrum, that is, at = 0.5 with = (E − E min )/(E max − E min ) where E min/max are the extremal energies of H. Full numerical diagonalization is possible up to N = 36 sites.We also use massively parallel shift-invert diagonalization [92] that allows to target eigenstates at = 0.5 on larger samples up to N = 52.We typically obtain about 100 eigenstates near = 0.5.Results are averaged over more than 500 (up to 10 000) realizations of disorder.To show the interest of working with constrained models, the Hilbert space of the largest cluster on which we can obtain exact eigenstates (N = 52) is almost twice as large as the one for N = 24 spins 1/2 (in the zero magnetization sector), which constitutes the state of the art for the analysis of MBL in spin chains [92][93][94].Obtaining ten eigenstates in the middle of the spectrum for one disorder sample at N = 52 requires about 20 min on 10 000 CPU cores.
In Sec.V, we also consider dynamical properties of various observables after a quench.The knowledge of all eigenstates allows to consider properties at arbitrary time up to N = 36.We also use Krylov expansion time-evolution techniques [95], which allow simulations up to time of the order of thousand plaquette flips on large clusters up to N = 64.Note that simulations on the largest sample N = 64 are extremely demanding (see sizes of Hilbert space in Table I).In the study of dynamics, we average over between 30 and 1000 realizations of disorder for each simulation set.

A. Spectral statistics
We first consider the simplest property of a many-body system, namely, its energy spectrum.The nature of the distribution of energy levels can already discriminate between different regimes: Ergodic systems display level repulsion whereas localized systems display a random distribution of level spacing characterized by Poisson statistics.In the manybody context, and in order to counter the fact that the density of states is not uniform, it is useful to consider the distribution of consecutive ratios as introduced in Ref. [96].We define gaps in the many-body spectrum as s n = E n − E n−1 and consider the consecutive reduced gap ratio r n = min (s n ,s n+1 ) max (s n ,s n+1 ) for which 0 r 1.The distribution P(r) and average value of r = 1 0 rP(r)dr, averaged over eigenstates and disorder realizations, have been computed for random matrix ensembles (corresponding to ergodic systems) and for independent energy levels (Poisson statistics) [97].For real Hamiltonians, such as Eq. ( 1), the random matrix ensemble to be considered is the Gaussian orthogonal ensemble (GOE).Figure 2 displays both the distribution P(r) and its average r for various values of disorder and system sizes.We compare them to predictions for the GOE ensemble [97]   , we find that the agreement with P GOE (r) for low values of disorder (V = 3) and P Poisson (r) for high values (V = 25, 30) is excellent with no adjustable parameter.Intermediate values of disorder V = 12, 17 present distributions which display a crossover behavior between the two limit distributions for the system size considered N = 40.The average gap ratio r also displays an interesting crossover between the two limiting cases with different system sizes showing an apparent crossing point (see the inset in Fig. 2(a) for square samples N = 32, 36, 40) around V 15-20.We note that rectangular samples also show the same crossover between the GOE and the Poisson limits but do not cross exactly at the same value which can be attributed to their different aspect ratios [98].The critical value of the gap ratio r * 0.392 is smaller than for the 1D standard MBL model [22], indicating that the putative transition point looks even closer to the Poissonian localized limit.Note that, for the largest disorder V = 30 considered here, we also checked on the largest system available to full diagonalization (N = 36) that no bandlike structure is present in the energy spectrum, a feature which was found [99] to mimic MBL on too small system sizes.

B. Eigenstate statistics
We next consider the similarity between eigenstates at the same energy density.It is well known that eigenstates in an ergodic phase look "similar" (e.g., they present similar expectation values for local observables), whereas eigenstates in a MBL phase are very different (e.g., local observables can differ strongly from one eigenstate to the following one in the many-body spectrum).To quantify the degree of similarity, it is useful to consider [22] the Kullback-Leibler divergence between two eigenstates |n and |n located nearby in the spectrum (here, at the same energy density = 0.5).It is defined as KL = i p i ln(p i /p i ) where p i and p i are the participation coefficients, i.e., the square amplitudes eigenstates in a given basis (here, the computational basis where dimer occupations are diagonal).The sum runs over the full Hilbert space, and the eigenstates are normalized such that i p i = i p i = 1.Eigenstates at "infinite temperature" ( = 0.5) in an ergodic phase should be very close to eigenstates of GOE random matrices for which one has KL GOE = 2. On the other hand, for eigenstates which are many-body localized, we expect KL to diverge with system size as the eigenstates are increasingly different.Figure 3 precisely shows this behavior with expectation values KL converging to their limiting GOE value at low disorder and diverging for large disorder.We also observe a crossing point at V 15, albeit slightly drifting, for different system sizes.This drifting behavior is also observed for 1D MBL in the random-field Heisenberg model [22].The value of KL at the transition point is larger than for the 1D spin model [22], pointing again to a putative critical point closer to the localized limit.around a plaquette in some eigenstates (e.g., left-middle plaquette for the leftmost eigenstate at V = 30).At intermediate disorder (V = 17), whereas, for the two first left eigenstates, dimer occupations look similar (at a superficial level), the right eigenstate configuration is markedly different.
Columnar imbalance.We now consider expectation values for another local observable, which will be useful later when discussing dynamics in Sec.V. We define the imbalance with respect to some specific state, the columnar configuration represented in the inset of Fig. 5 as (2) Phase ϕ p carried out by each plaquette is simply chosen to be 0 or π on a columnar arrangement (see the inset of Fig. 5) such that the imbalance of this perfect columnar state is maximal and equal to I = 1.At thermal equilibrium and infinite temperature, the imbalance vanishes on average.
We display, in Fig. 5, the probability distribution P(I ) of n|I|n taken over eigenstates |n located at = 0.5 for different values of disorder and computed on the tilted square lattice with N = 40 sites.For small values of disorder (left column), the probability distribution is Gaussian around the mean value 0 as expected for an ETH phase.The typical width of the distribution (which is very small for the smallest disorder V = 1) broadens with disorder.Close to V 15, peaks for specific values of the imbalance start to appear in the distribution which become clearly visible for all strong values of disorder V 20 (right panels).The distribution of imbalance for pure dimer configurations P conf (I ) (obtained by  computing I for all dimer configurations in the Hilbert space) is represented in the V = 30 panel for comparison (there are 1/N values possible for the imbalance I).The location of the peaks of P(I ) matches precisely those of P conf (I ) with, however, different amplitudes.This indicates that, at large V , eigenstates are close (but not exactly equal) to pure dimer configurations, in agreement with the qualitative view provided by the snapshots of Fig. 4.

D. Participation entropies
To understand the localization properties of the eigenfunctions in the Hilbert (configuration) space, we next consider the participation entropy, where p i = | n|i | 2 and the sum runs over all basis states |i of the computational basis.At low disorder, we expect eigenfunctions to be fully delocalized with a participation entropy scaling as S p = a ln(H) − • • • (with a = 1 for a fully delocalized state and • • • denoting finite-size corrections).In the large disorder regime, on the other hand, based on data for the MBL phase in the 1D random-field Heisenberg [94], we do not expect a true localization in configuration space but, instead, a volumic multifractal behavior where S p = a ln(H) with a strictly less than 1 (and increasingly small as disorder is increased).Figure 6 presents the numerical data for the average participation entropy divided by the logarithm of the size of the Hilbert space S p / ln(H).Two regimes can be clearly distinguished: one at low V where data tend to their maximal value (a = 1) with system size and another at larger disorder where curves tend to group together towards sensibly smaller values of the volume law coefficient a.This is in accordance with the results for the MBL phase in one dimension [94].We find that data for S p / ln H for different sizes tend to regroup themselves around V 15-20, consistent with previous estimates of the possible transition point.

E. Entanglement entropy
Many-body eigenstates in ergodic or localized phases have strikingly different entanglement properties.Consider the von Neumann entanglement entropy of an eigenstate |n for a system divided in two parts A and B,

S(|n ) = −Tr
with ρ A as the reduced density matrix obtained by tracing out bond degrees of freedom in B ρ A = Tr B |n n|.For our samples, we choose equal parts A and B containing identical number of lattice bonds and sites (see the Appendix for a graphical representation of the entanglement cut).The entanglement entropy S, averaged over eigenstates and realizations of disorder, is represented in Fig. 7.In the ETH phase, a volume law of the entanglement entropy is expected with S scaling with the number of degrees of freedom in A. This is what is found in the low-V regime of Fig. 7 with the following caveat: Due to the constrained nature of the Hilbert space, not all configurations of A are compatible with those of B (see a discussion in Ref. [100] and in the Appendix), and the number of these "constraint sectors" in ρ A depends, in particular, on the area between A and B. This results in different samples with the same volume having different constraints on ρ A when bipartitioned.This is clearly seen when computing entanglement entropy on random states in the different samples denoted as S random and represented by the star symbol in Fig. 7: For instance, the 8 × 4 rectangular lattice has a much higher entropy than the 4 × 8 lattice, even though both have the same volume (number of bonds and sites) for A. Figure 7 shows that the low-V entanglement entropies of all samples converge to the random-state values, confirming clearly that the low-V phase fulfills the ETH.For larger values of V , the entanglement entropy drops down for all samples.The inset of Fig. 7 shows a reasonable scaling of S with area A at large enough V , in agreement with the expected area law scaling expected in a localized phase.Finite-size effects combined with the constraint effect on S discussed above prevent us from pinpointing a transition/crossover point between the volume and the area law scaling with great precision.
In previous studies of the behavior of entanglement entropy close to the MBL transition, the variance of entanglement entropy was shown [16,101] to peak when approaching the transition, which can be rationalized as the result of the crossover between the volume (area) law of entanglement entropy in the ETH (MBL) phase.This is also seen in Fig. 8 where the standard deviation std(S) of the entanglement entropy (computed over all eigenstates and realizations of disorder), scaled by the entanglement entropy of a random state S random , is represented.Focusing on the square samples, a right shift of the peak with system size towards larger values of V is apparent: A similar shift was always found to be present in numerical studies of the MBL transition in 1D [16,101].For the largest square (rectangular) sample, the peak is located at V = 10 (V = 12), providing a lower bound for the possible transition.

IV. MACHINE LEARNING ANALYSIS OF ENTANGLEMENT SPECTRA
As a complementary approach, we study the quantum dimer model Eq.(1) using machine learning techniques and follow the supervised scheme introduced in Ref. [102].As inputs representative of the two phases, we provide entanglement spectra (i.e., the eigenvalues λ i of the reduced density matrix ρ A ) obtained deep in the ETH and the MBL phases, and we train a neural network to classify them accordingly.The supervised learning has been used in multiple occasions in the context of delocalization-localization transitions [103][104][105][106][107][108] producing results in good qualitative agreement with more conventional analyses.In particular, FIG. 9. Entanglement spectra (eigenvalues λ i of the reduced density matrix ordered by amplitude) used in the neural network approach to label the ETH (V = 1, blue colors) and MBL phase (V = 30, orange colors) for the N = 36 square sample.For both cases, 100 spectra obtained from different realizations of disorder are represented.The dashed line represents the average of all spectra for V = 30.The inset: Higher part of the average entanglement spectra for different sample sizes for V = 30.The log-log scale highlights a power-law behavior for the larger eigenvalues λ i .
Refs. [103,105,106] showed the interest of using labeled entanglement spectra.
For this scheme, we consider the N = 32, 36, 40 square samples and the largest rectangular sample N = 6 × 8.For each of them, we provide 10 000 entanglement spectra (2000 for N = 40, 6 × 8) including between 100 and 200 disorder realizations per disorder strength at V = 1 for the ETHlabeled phase and V = 30 for the MBL-labeled phase.The spectra being rather large (1972 values for N = 32, up to 21 286 for N = 6 × 8), we found that sorting them allowed for both perfect training and test accuracies for all system sizes.Figure 9 represents the form of the entanglement spectra used in these two limits for the sample N = 36 for 100 different realizations of disorder.The very similar form for various disorder samples in the ETH phase clearly contrasts with the stronger dispersion observed in the MBL phase.In the latter, superposing the entanglement spectra for various samples sizes (inset of Fig. 9) highlights that the larger values of the spectrum decay as a power law, similar to what is found in the 1D MBL [109].
We used a fully connected neural network consisting of one hidden layer of 32 neurons followed by two softmax output neurons.We follow a cross-validation procedure where we randomly selected half of the dataset to form the training dataset, the rest being assigned to the test set.This process is repeated multiple times, generating new training and test partitions each time.This allowed us to track whether the neural networks were overfitting depending on the training conditions.Namely, we checked that data from V = 1 and V = 30 give perfect training and test accuracies for each system size.
Figure 10 displays features that are consistent with the previous analysis.The left panel displays the average output of the neural network for the different samples as a function of disorder strength.At low V , the machine learning analysis validates a fully ETH phase (i.e., where all samples are classified ETH) that extends up to a value V = V 1 , and at large V , a fully MBL phase (with more than 99% accuracy) at large for V V 2 (V 1 6 and V 2 20 for the largest N = 6 × 8 sample).Note how these bounding values (in particular, V 1 ) shift to higher values of disorder with system size.This reflects that ETH is easier to disrupt on a too small sample, in perfect agreement with the trend in all other observables discussed in Sec.III.We find no crossing point with system size in the current data, different from what is observed with the 2D Ising model [102].
The right panel of Fig. 10 shows the standard deviation of the neural network output as a function of disorder.The standard deviation is low in both limits where the phases are well distinguished (at low and large V 's) and peaks at an intermediate value of disorder.The location of the peak (which shifts with system size) is the point where the neural network has most difficulties to classify phases.It can be interpreted as a finite-size estimate of a possible transition point.Note the similarity between the standard deviation of the neural network output and the standard deviation of the entanglement entropy (Fig. 8), in particular, the positions of the peaks are almost the same for both quantities for the different sample sizes.
In conclusion of this section, we find that a neural network only fed with entanglement spectra is able to learn how to correctly distinguish the ETH and MBL phases for the quantum dimer model with random potential Eq. ( 1) as well as to provide finite-size estimates of the transition point between the two.This automated method gives results in very good agreement with the analysis based on more standard featureengineered estimators of the phases presented in the previous Sec.III.Our results also indicate that the machine learning analysis is also subject to finite-size effects, indicating that the use of numerical data coming from a single size may not be able to provide quantitative results for the study of phase transitions (see a discussion in Ref. [110]).One noticeable interest of the neural network analysis (already pinpointed earlier [103,105]) is that the required amount of data and overall computational cost is considerably lower than with more traditional observables to obtain approximately similar quality of prediction: For instance, good statistics on the gap ratio (Fig. 2) requires 100 times more realizations of disorder than with the machine learning analysis (Fig. 10).

V. DYNAMICAL PROPERTIES AFTER A QUENCH
We consider here dynamical properties after a quench of the model Eq. ( 1) in order to probe for the arrested transport dynamics or slow propagation of information expected in a MBL phase.We perform quenches from an initial product state (with no entanglement) and take, here, the columnar configuration depicted in the inset of Fig. 5 (a product state).This columnar configuration is the most flippable state and is well connected to other dimer configurations by the action of the Hamiltonian.We, thus, intuitively expect this state to rapidly thermalize if it does.The choice of a columnar configuration echoes the experimental protocols for which the system is often initialized in a charge-density state [31,57] or numerical simulations which choose a Néel state in spin chains, both states sharing the same properties as the columnar configuration.

A. Dynamics of imbalance from the initial state
We start with the dynamics of the columnar imbalance defined in Eq. (2).By definition, I (t = 0) = 1, whereas I = 0 when averaged over all eigenstates (see Fig. 5).In the long-time limit, I (t ) is, thus, expected to vanish in a thermal phase, whereas a nonvanishing long-time value would signal a localized phase.Note that in finite systems, I (t → ∞) never strictly vanishes.Figure 11 displays the columnar imbalance as a function of time with panels corresponding to different disorder strengths V .In the low-disorder limit, the imbalance quickly decays to a vanishing value on short timescales (this is most readily seen at V = 1).When disorder is increased, decays become slower, and for intermediate disorder values V = 10, 15, 20, the saturation value is not reached for the largest clusters within the timescale studied (for N = 32, 36 data not shown at longer time show that the saturation is reached within less than a percent for t = 1000 for all disorders considered).We can, nevertheless, observe that the imbalance clearly decreases with system size for the largest time available, in a manner compatible with a vanishing value in the thermodynamic limit for V 15.For larger values of disorder, one clearly observes that the saturation value is much larger and, within available time and cluster sizes, the imbalance does not tend to vanish: This is most speaking for V = 30 where more than 80% of the initial imbalance is retained even on the largest clusters.Pinpointing exactly a transition point with dynamics is difficult due to the timescale limitations for the largest clusters, but the present data appear to be compatible with a transition/crossover for values of V

B. Entanglement and participation entropies dynamics
The nature of the growth of entanglement entropy S(t ) after a quench is a hallmark feature to distinguish MBL in one dimension where it takes a logarithmic form with time [18,19] from an ETH regime (with ballistic growth S(t ) [111]) or an Anderson localized regime [where it saturates to a constant O (1)].In both MBL and ETH regimes, the saturation value at long times is expected to fulfill a volume law with a prefactor equal to (smaller than) the corresponding thermal entropy for the ETH (MBL) case [112,113].
We represent, in Fig. 12, the dynamics after a quench from the columnar configuration of the average entanglement entropy S(t ) normalized by the entanglement entropy of a random state S random [which is the thermal entropy up to O( 1) corrections] for different disorder strengths.The ETH regime with a rapid growth of entanglement is readily identified at low disorder V = 1, where S(t ) reaches its saturation value (close to the maximum value reached by a random state) in less than ten plaquette flips.In the top panels where the same scale is used for the vertical axis, the entanglement entropy is observed to grow slower towards its saturation value with visible finite-size effects.Again, the time and size limitations do not strictly allow to draw definitive conclusions, but it appears very likely that for large enough clusters, the entanglement entropy will reach its thermal value at a long time for V 15.In the bottom panels (disorders V = 20, 25, 30), on the other hand, the entanglement entropy develops at an even much slower pace.The logarithmic horizontal scale chosen for the bottom panels highlights the striking feature that the growth of entanglement entropy appears to behave as ln(t ) for the largest clusters which are not hampered by finite-size saturation effects.Note that for the largest data at V = 15, a logarithmic growth is also compatible with the data (not shown).Naively exporting the argument (based on the exponentially decaying interactions between local integral of motions) that explains the logarithmic growth of entanglement in 1D MBL, we would also expect a logarithmic growth in 2D for the entanglement cut geometry chosen here as entanglement can grow only in one direction.
We finally present the evolution after a quench of the participation entropy S p (t ) in Fig. 13.The dynamics of the participation entropy is almost identical to the one of the entanglement entropy discussed above: At small disorder, rapid saturation to the maximum attainable value ln(H), slower growth at intermediate values and logarithmic growth for the largest samples for strong disorder (note again the logarithmic scale on the bottom panels of Fig. 13).The only noticeable difference in the dynamics of both entropies comes from different finite-size effects at short times with participation data showing an improved nesting of curves as system size is increased which is not observed for entanglement.From participation entropy, we conclude again on slow dynamics starting from V 15-20 consistent with the existence of a localized phase.

VI. DISCUSSION AND PERSPECTIVES
We presented a comprehensive and extensive large-scale exact numerical study of the localization properties of a twodimensional constrained quantum many-body system with disorder.All data both for eigenstate or dynamical properties can be interpreted as consistent with the existence of two distinct phases: An ETH phase at low disorder and a manybody localized phase at strong disorder, which are separated by a transition located around disorder strength V 15-20.Our evidence for a MBL transition in the 2D quantum dimer model comes from numerical simulations on finite lattices.Of course, finite-size simulations can always be argued to artificially detect a MBL phase even when only a ETH phase occurs in the thermodynamic limit.We would like to emphasize that the level of numerical evidence for a MBL transition in the 2D quantum dimer model is, in our opinion and experience, similar to the one obtained for the standard model of MBL in one dimension [21,22] with similar or larger Hilbert space sizes and timescales probed in the current paper.This is, of course, not a definite proof that a MBL transition occurs in the thermodynamic limit.
What is perhaps more important with respect to potential experiments is the fact that even if the behavior at large values of disorder (say V 25) is ultimately ergodic, the timescales and/or system sizes needed to probe ergodicity would be extremely large.An experimental realization of Eq. ( 1) or of a similar constrained model in 2D would see localization for all practical purposes on the timescales available in the laboratory.We deliberately avoided to attempt to perform a data collapse analysis to extract, e.g., critical exponents due to the two following limiting points: The moderate linear system sizes that we can reach as well as the nonuniform aspect ratio between available samples would likely provide biased estimates for critical exponents.For one-dimensional MBL, the larger linear length scales that can be reached have been argued (see, e.g., Ref. [114]) not to be large enough to provide correct estimates of asymptotic critical behavior.The situation is likely to be the same here.With these numerical limitations in mind, we can, nevertheless, observe that critical values of the gap ratio or of the Kullback-Leibler divergence of eigenstates are closer to their Poisson than their ETH limits, indicating that the putative transition point is even less ergodic than for the one-dimensional MBL transition in the random-field Heisenberg spin chain [22].
Going back to the l-bits arguments which predict the asymptotic unstability of MBL in two dimensions [44][45][46], we remark that we do not know how to simply construct l-bits for the quantum dimer model.The crucial point is that there is no tensor product structure for the Hilbert space of quantum dimer models, thus, constructing a basis of commuting operators with local support appears ill defined.As the bubble argument is explicitly based on l-bits, it cannot strictly speaking apply here.
Regarding experiments, the relative strong disorder needed to induce the MBL phase does not provide immediate relevance of our results in frustrated magnets where QDMs naturally emerge.However, the artificial realizations of quantum dimer models (and other related 2D constrained models) have been extensively studied in cold-atomic and mesoscopic systems.Precise experimental protocols have been proposed using Ryberg atoms [115,116], trapped ions [117], large-spin ultracold atoms [118], Josephson-junction arrays [119], or superconducting quantum circuits [120].Further explorations are needed to understand which experimental platform is the most suited to implement random potential or interactions.Constrained models in 1D have already been implemented (with no randomness) exploiting the Rydberg blockade in atomic chains [64], already exhibiting rich unusual quantum dynamics.
There are several perspectives opened up by our paper.First, the road map to two-dimensional MBL can be exploited using the QDM on other two-dimensional lattices, allowing to test for universality and to search for other features.A recent investigation [121] considered the honeycomb lattice, which has an effective smaller local Hilbert space allowing to reach larger samples and obtained results similar to ours.The QDM on the kagome lattice is also an interesting case worth pursuing at it possesses conserved Z 2 quantum numbers, allowing the exciting possibility of 2D topological order in MBL states [89].The possibility of MBL in other quantum constrained models, such as quantum ice or loop models also provide an interesting follow-up.The large-scale numerical techniques that we use in the present paper would be useful to investigate other aspects of the model Eq. ( 1) not addressed here, in particular, to make connection with 1D models where MBL is best understood.The current investigation focused on the "infinite temperature limit," that is, excited states in the middle of the spectrum ( = 0.5).It would be interesting to construct a full energy-disorder phase diagram to see if a many-body mobility edge exists as in 1D spin chain models [22].Also the detailed dynamics occurring in the ETH phase remains to be investigated.The relatively slow dynamics exhibited around V ∼ 10, 15 in Figs.11-13 could be signs of a precursor subdiffusive regime, similar to what is observed in 1D systems [122][123][124].Finally, it would be intriguing to see if the machine learning techniques used in Sec.IV would be able to distinguish 1D from 2D MBL, and if not, to use neural networks trained on 1D spin chain models to probe 2D MBL (transfer learning).
On the analytical side, the absence of an obvious scheme for constructing l-bits points towards the interest of building an alternative effective description for the many-body localized phase in such 2D constrained models.Recall as well that there is no free-particle limit to start with in the model Eq. ( 1).The height description [125] of dimer coverings (and of constrained models in general) appears as a good candidate to develop such a framework.Finally, from a purely classical perspective, we note that the strong disorder limit V → ∞ constitutes an interesting classical statistical mechanics problem (with, for instance, nontrivial degeneracies), unsolved to the best of our knowledge and distinct from other interacting dimer models [126].We use regular and tilted square lattices as well as rectangular lattices with N sites and 2N bonds.Regular square lattices have N = L 2 sites, and we used samples with L = 4, 6, 8. Rectangular lattices, which have N = L x × L y sites (with L x,y = 4, 6, 8, 10), have a nonunity aspect ratio, which can affect the analysis of finite-size effects.Finally, tilted square lattices host N = p 2 + q 2 sites, and we use the combination p = 4, q = 4 (with N = 32 sites), p = 6, q = 2 (N = 40), and p = 6, q = 4 (N = 52).Tilted square lattices have a unity aspect ratio, but they do not have all symmetries of the full regular square: This is, however, irrelevant as the Hamiltonian Eq. (1) does not allow any lattice symmetries due to disorder.All these lattices, which are represented in Fig. 14, have in common that they host a (W x , W y ) = (0, 0) winding sector: This is true when L, L x , L y , p, and q are even.One can also easily define the columnar state on these samples (see the example for N = 52 in Fig. 14).

Entanglement cut
When computing bipartite entanglement properties, there is some freedom in choosing where is located the entanglement cut, i.e., which degrees of freedom belong to A or B in Eq. (4).We choose the cuts displayed in Fig. 14 where all lattices are divided in two parts with equal number of bonds.
With this cut, the area A between A and B is expected to scale with L for regular square lattices, L y for rectangular lattices, and p 2 + q 2 for tilted square lattices.
This choice of cuts for the regular square, rectangular, and N = 32 tilted square lattice also induces the following property: One of the two winding numbers (the one associated with a line parallel to the cut, e.g., W y for the rectangular lattices) is conserved by the cut but not the other one.The reduced density matrix on A is block diagonal with blocks labeled by winding numbers in A (W A y which is forced to be 0 and W A x ), which we take advantage of in numerical computations.The existence of blocks in the reduced density matrices and its effect on entanglement has been discussed in other models with conservation laws [100,[135][136][137].
We note also a supplementary specificity of constrained models for entanglement properties (see the discussion in Ref. [100]).Due to the constraint and even within the conserved-number blocks (e.g., configurations of A with W A x and configurations of B with W B x such that W A x + W B x = W x = 0), not all configurations of A are compatible with those of B. This property has an important consequence, for instance, when computing volume-law entanglement as different samples with the same volume can have different constraints on configurations in A and B when bipartitioned.This is illustrated in the main text with the example of the 8 × 4 and 4 × 8 rectangular samples.This property can also cause difficulties in comparing scaling of entanglement with system size for different samples.

FIG. 1 .
FIG.1.Quantum dimer models.(a) A configuration of dimers is valid when each lattice site is touched by one and only one dimer.(b) Illustration of the quantum dimer model on the square lattice with a disordered potential Eq.(1).Plaquettes p of the lattice where a pair of parallel dimers is present (highlighted in green) contribute a potential V p (different for each plaquette) to the potential energy of this configuration.The kinetic term flips dimers on such plaquettes such that the orientation of the two dimers on this plaquette is changed.In the bottom right panel, the flip of the plaquette p = 20 (surrounded by red lines) is represented, resulting in a new configuration with a different potential energy.(c) On the square (bipartite) lattice with periodic boundary conditions, all dimer configurations can be classified by a pair of winding numbers (W x , W y ) defined as follows.Consider a horizontal (vertical) line-such as represented in the figure-perpendicular to lattice bonds.W x (W y ) is then equal to the algebraic number (i.e., ±1 where the sign is defined by the sublattice on a fixed side of the line) of dimers in this configuration crossing this line.The sign ±1 is represented by arrows on lattice bonds.Note that the winding numbers do not depend on the position of the line, and they cannot be changed by a local reconfiguration of dimers.The represented configurations all belong to the (W x , W y ) = (0, 0) sector in which all calculations are performed in this paper.Periodic boundary conditions are represented with the help of duplicated gray sites on the boundaries.

FIG. 2 .
FIG. 2. Statistics of the reduced gap ratio r for eigenstates located in the middle of the spectrum = 0.5.(a) Average gap ratio r for different 2D samples as a function of disorder.Square samples are color highlighted.(b) Probability distribution P(r) for N = 40 sites tilted square lattice for different values of disorder.

FIG. 3 .
FIG.3.Average Kullback-Leibler divergence KL as a function of disorder for different system sizes.Square samples are color highlighted.

C.
FIG. 4. Dimer occupation n b for three consecutive eigenstates (from left to right) in the middle of the spectrum, for a single disorder realization for different values of disorder strength (from top to bottom).The width of orange dimers located on each bond is proportional to n b .Data for the N = 52 tilted square lattice.

5 .
FIG. 5. Probability distribution of columnar imbalance I for different values of disorder.On the bottom right panel is also represented (circles) the distribution for pure dimer configurations P conf (I ) (normalized such that P conf (I = 0) = P V =30 (I = 0).Data for the tilted square lattice with N = 40 sites with more than 10 000 eigenstates for each value of disorder.The inset (top right): definition of the phase φ p used in the definition of the columnar imbalance, on top of the columnar configuration for which I = 1.

FIG. 6 .
FIG.6.Average participation entropy divided by the logarithm of the size of Hilbert space as a function of disorder for different system sizes.Square samples are color highlighted.

FIG. 7 .
FIG.7.Average entanglement entropy S of subsystem A as a function of disorder for different system sizes.The star symbol on the V = 0 axis is the entanglement entropy of a random state S random .The inset: Area law scaling of the average entanglement entropy divided by the size of the surface A between A and B as a function of disorder (see the Appendix for definition of A for tilted clusters).Square samples are color highlighted.

FIG. 8 .
FIG.8.Standard deviation of entanglement entropy std(S) of subsystem A as a function of disorder for different system sizes, scaled by the entanglement entropy of a random state S random .Square samples are color highlighted.

FIG. 10 .
FIG.10.As a function of disorder strength V and for different system sizes: (a) mean and (b) standard deviation (right panel) of the neural network output p, defined such that the labels are p = 0 (p = 1) in the ETH (MBL) samples at V = 1 (V = 30).The neural network predictions involved 5000 entanglement spectra per disorder (including 100 disorder realizations per disorder).Error bars show standard deviation over ten instances of neural networks from the cross-validation procedure, stopped after 100 epochs.Square samples are color highlighted.

30 FIG. 11 .
FIG. 11.Decay with time of columnar imbalance I(t ) after a quench from the columnar configuration for different values of disorder V and different system sizes N of square clusters.between 15 and 20, similar to what is found with the analysis of eigenstate data in Sec.III.

30 FIG. 12 .
FIG.12.Average entanglement entropy S(t ) after a quench from the columnar configuration, scaled by the (extensive) entanglement entropy of a random state S random for different values of disorder V and different system sizes N of square clusters.The top four panels are represented using the same vertical scale and a linear scale for the time axis.The bottom three panels use a logarithmic scale for the time axis, highlighting the slow growth of entanglement for the largest values of disorder.

30 FIG. 13 .
FIG.13.Growth of the average participation entropy S p (t ), normalized by its maximal value ln(H) after a quench from the columnar configuration for different values of disorder V and different system sizes N. As in Fig.12, the top four panels use the same vertical scale, and a linear scale for the time axis whereas the bottom three panels use a logarithm scale for the time axis.

FIG. 14 .
FIG. 14. Lattices used in this paper: regular square lattices with N = L 2 sites (top left panel with L = 6), rectangular lattices with N = L x × L y (top middle panel with L x = 6 and L y = 8), and tilted square lattices (bottom row with N = 32, 40, 52 sites).The red/blue bond decomposition is used for the entanglement partition.Top right: columnar state for the tilted N = 52 square lattice.

TABLE I .
Two-dimensional lattices used in this paper with their shape [regular or tilted square (sq.) lattices, or rectangular] and the size of their Hilbert space in the (W x = 0, W y = 0) winding sector.