Probing many-body localization in a disordered quantum dimer model on the honeycomb lattice

We numerically study the possibility of many-body localization transition in a disordered quantum dimer model on the honeycomb lattice. By using the peculiar constraints of this model and state-of-the-art exact diagonalization and time evolution methods, we probe both eigenstates and dynamical properties and conclude on the existence of a localization transition, on the available time and length scales (system sizes of up to N=108 sites). We critically discuss these results and their implications.


I. INTRODUCTION
Localization in disordered, interacting quantum systems 1,2 is a topic that has recently received wide attention due to the very peculiar phenomenology [3][4][5][6] , the foundational issues about quantum integrability and ergodicity involved 7,8 , and the increased precision and control on experimental realizations 9,10 . Systems with a many-body localization (MBL) transition typically exhibit two phases, one at low disorder which obeys the eigenstate thermalization hypothesis (ETH) and one at high disorder which exhibits no transport, no thermalization [11][12][13][14] and emergent integrability due to an extensive number of quasi-local integrals of motion [15][16][17][18][19] . Furthermore, localized states have low entanglement at any energy and obey an area law, a property usually valid for ground states only 20,21 . Finally, localization in interacting systems is characterized by the very slow spreading of information, namely the entanglement [22][23][24] , and the total absence of transport for local observables 1,2 . All these features have contributed to make MBL a compelling physical phenomenon, including with respect to quantum information processing protocols 20,25-27 . In the context of the study of MBL transitions, a wide range of results, outlining the phenomenology described above, have been produced for one-dimensional (1D) systems [3][4][5][6] . Remarkably, a proof of the existence of the MBL transition has been obtained for a 1D quantum Ising model with a transverse field 28,29 . In higher dimensions, however, no such proof exists. One generally expects that in higher dimensions delocalization is favoured due to the increase in channels for the delocalizing terms, similarly to the phenomenology of Anderson localization in higher dimensions. More specifically, general arguments based on the existence and size-scaling of thermalizing bubbles support the absence of localization for large enough times 30,31 , even though no rigorous proof was obtained either.
A number of results on 2D systems have notably been presented. Experimental results obtained in cold atoms setups interestingly show absence of dynamics and localization at high disorder 10,32 . At present, this exper-imental evidence is arguably of higher quality than the analytical and numerical modeling of MBL in 2D. Numerically, a number of approaches have been explored in 2D lattice models, using both unbiased and biased methods, and showing indications of a localized phase [33][34][35][36][37][38] . Other simulations conclude in favor of absence of MBL 39 . However, the main limit of numerical approaches is the small system sizes and/or time scales that are reachable in the computations. The size of the Hilbert space and thus of the quantum problem grows exponentially with the number of particles N in the system while the physical length scale of the sample grows as a square root of N . For unbiased methods this is an especially strong constraint, effectively limiting the analysis to systems up to around 20 spins 1/2. While in one dimension several different lattice sizes can fulfill this requirement, thus allowing in principle finite size scaling to be performed, this is no longer the case in two dimensions where the number of system sizes are greatly limited. While larger system sizes can be reached using methods geared towards capturing properties of an MBL phase [33][34][35][40][41][42] , these methods are not unbiased and by construction will miss the ergodic phase or the phase transition.
Here, we aim to investigate an MBL transition in a specific system up to a real-space size as large as possible and with unbiased methods. We do this by considering a highly constrained model and state-of-the-art numerically exact methods 43 . Specifically we consider a disordered quantum dimer model (QDM) on a honeycomb lattice, where each lattice link is either free or occupied by a dimer with the constraint that each lattice site is touched by one and only one dimer [44][45][46] . An immediate consequence for this is that the dynamics of such a model is very constrained: single dimer moves are not allowed and the simplest move involves an hexagonal plaquette. Moreover, this constraint also automatically encodes strong interactions which for the honeycomb lattice already imply long-range correlations in the statistical ensemble of dimer coverings. The interplay between a constrained dynamics, which favors slow dynamics and localization 47 , and the strong interactions, which favor delocalization, creates an ideal situation for arXiv:2005.10233v1 [cond-mat.dis-nn] 20 May 2020 an MBL transition to exist. Finally, we note that such models are based on Hilbert spaces that, due to the constraints, have considerably lower dimension compared to spin systems: for N 1/2-spins, the Hilbert space size is 2 N , while it scales only as ≈ 1.175 N 44 for a dimer system on a N -sites honeycomb lattice, giving an obvious numerical advantage for large system sizes. A previous work has analyzed a similar disordered QDM on a square lattice 48 . Here, we substantially push forward this analysis, almost doubling the maximum system size reached, by turning to the honeycomb lattice instead.
The article is structured as follows. In Section II we detail the model Hamiltonian, the symmetry sectors and the lattices used as well as the procedures used to obtain the numerical results. Such results are outlined in Sec-tion III, first considering observables within exact midspectrum eigenstates, and, secondly, the dynamical properties obtained with Krylov time evolution. Finally, we provide conclusions in Section IV. In the appendix we discuss in detail the lattice clusters used in the numerical analysis (Appendix A), further energy-resolved quantities (Appendix B) and comparisons with the entanglement properties of specific states (Appendix C).

II. MODEL
We consider the following quantum dimer model on the honeycomb lattice 45,46,49 with a random potential: The first term, an hexagon "flip", is a kinetic term. The second term is a disordered potential on each flippable hexagon; the v p are drawn from a uniform distribution in [−V, V ].
We construct lattices with N = 42, 54, 72, 78, 96 and 108 sites; in Fig. 1 (a) we show the N = 72 lattice and we refer the reader to the Appendix for more details on the other clusters. On the honeycomb lattice with periodic boundary conditions, the constraints due to the dimers and to the allowed plaquette moves are such that two conserved quantities, the winding numbers, exist. The winding numbers are defined as the sum along a line parallel to the x or y axis; having labeled the honeycomb lattice  Fig.  1 (c)). Among the sectors with conserved total winding number, we select the one for which w x = w y = 0, which is the largest one. We remark that, for finite lattices, not all lattice shapes allow the existence of this zero winding sector; we discard lattice shapes that do not satisfy this requirement 50 . Table I displays the number of allowed coverings in the zero winding sector, which correspond to the size N H of the Hilbert space. The number of nonzero elements in the matrix is also noted, which, in addition to matrix size, contributes to limiting the feasibility of the numerical calculations. We perform exact diagonalization on some of these lattices (up to size 78). We use either full diagonalization or shift-invert methods 43 to obtain around 100 eigenstates at the center of the spectrum. We also study the dynamics of nonequilibrium initial states though Krylov subspace time evolution methods for all lattice sizes 51 . In all cases, we average over disorder realizations of the random potential (at least 1000 for most system sizes and around 100 for the dynamics on the largest one).

III. RESULTS
We consider various quantities with known different behaviors in the MBL and ETH phases. We analyze spectral, eigenstate, and entanglement properties as well as the dynamics of the system.

A. Spectral properties
Spectral gap ratio We start by analyzing the spectral properties of the two phases. Specifically, we consider the energy level gap ratio 14 : where s i = E i+1 − E i is the gap between two adjacent eigenvalues. We average in a small window of about 100 eigenstates around the center of the spectrum as well as over disorder realizations. Depending on the level gap statistics, r ≈ 0.39 for a Poisson distribution in the localized phase and r ≈ 0.53 52 for a Wigner-Dyson distribution corresponding to the ETH phase. In Fig. 2, top panel, we show the value of r as a function of the disorder for various system sizes. It appears that both localized and ETH phases are captured with the available cluster sizes. The transition value can typically be inferred by where the curves for increasing size cross, as it denotes opposite flows in the system size scaling in the two phases. We note that here the crossing point has a noticeable drift towards higher V values.
In the bottom panel of Fig. 2 we show the probability distributions of the gaps s of the unfolded spectrum for various values of the disorder V , showing excellent agreement with a Poissonian or a Wigner-Dyson distribution (shown in black) for high and low V respectively.
For the smallest sizes N = 42 and N = 54, we additionally computed the gap ratio as a function of the energy density (not just for the middle of the spectrum), see Appendix B.

B. Eigenstates
Kullback-Leibler divergence for energy-adjacent eigenstates We now consider quantities characterizing eigenstate properties which have been shown to be good indicators of localization. In the localized phase, eigenstates and local observables close in energy are very different in structure, as opposed to the ETH phase. Thus, we consider the Kullback-Leibler divergence for two consecutive eigenstates |ψ and |ψ in the spectrum, defined as where the sum runs over the N H elements |b i of the Hilbert space basis. We expect KL to approach to KL GOE = 2 (the value obtained for the Gaussian orthogonal ensemble of random matrices) in an ETH phase and to diverge with system size in a localized phase 12 . We show the results for KL in Fig. 3 as a function of the disorder strength V . The limit value KL = 2 is well captured at small disorders V , as well as a crossing point between the N = 54, N = 72 and N = 78 clusters (N = 42 appears to show stronger deviations due to the small size), with some drift due to finite-size scaling, suggesting a localization transition around V ≈ 22 − 25. Eigenstate participation entropy In a similar manner, we consider the participation entropy of the eigenstates, which gives information about localization in the Hilbert space 12,53 . It is defined as For a state which is localized in the Hilbert space, S p is of O(1). For many-body localized states, a multifractal behavior is expected in this computational basis 53 , with a participation entropy behaving as S p ∝ a ln N H , with a < 1. For extended states in the ETH regime, S p will scale as ln N H , with a = 1.
In Fig. 4 we show the participation entropy, rescaled by ln N H (i.e. this ratio is the coefficient a up to higher order corrections), as a function of the disorder V . At low disorder we see that a has a high value which is likely to scale to 1 with increasing size. A different behavior onsets at around V ≈ 20 − 25: the curves for different system sizes join and collapse, suggesting a finite a < 1 asymptotically for disorders larger than this value.
Eigenstate imbalance We next consider the imbalance of the eigenstates with respect to a specific configuration where the state used as a reference is chosen as the basis element of the so-called star configuration displayed in Fig. 1b. We define the (complex) imbalance as The phases φ p assume three possible values: 0, 2 3 π and − 2 3 π, depending on the dimer configuration on the plaquette p (see Fig. 1b). With this definition, the imbalance of the reference basis state in Fig. 1b has I = 1 and maximum amplitude |I| 2 = 1. Delocalized eigenstates will have a probability distribution for the modulus squared imbalance which is sharply peaked in 0. On the other hand, if states are localized and close to basis states, the imbalance will be peaked around the values corresponding to dimer configurations of the basis states. The probability distribution of the modulus squared of the imbalance is shown in the two panels of Fig. 5 respectively for low (top panel) and high (bottom panel) values of disorder for the largest system size (N = 78). As expected, the imbalance distribution is sharply peaked at 0 for very small values of disorder with an increasing variance for higher disorders. Between V = 15 and V = 20 the distribution broadens and develops peaks at values |I| 2 > 0 which, at higher disorder (V ≥ 25), are shown to closely correspond to the imbalance of configuration basis states (shown in dashed lines).
As for 1D systems, the imbalance is a quantity especially useful for characterizing the dynamical properties of the system in different phases. We will further analyze dynamics of the imbalance after a quench in Sec. III D.
Eigenstate dimer bond occupation We finally consider a local observable, the dimer bond occupation, and specifically the probability distribution of where the operator n k acts on the basis vectors |b i as n k |b i = 1 if bond k is occupied in b i and 0 otherwise. In the limit of a uniformly extended state, all three bonds belonging to a site have the same probability to be occupied, i.e. 1/3. As shown in the main panel of Fig. 6, for a delocalized eigenstate, this translates into a probability distribution of O sharply peaked at 1/3 for low disorder values, with an increasing variance for higher disorders. At disorder between V = 15 to V = 20 the distribution becomes bimodal with two peaks at O = 0 and O = 1, meaning that the eigenstates start to resemble some given dimer configurations. In the limit of infinite disorder, where the eigenstates coincide with the configuration states, the distribution is 2/3 δ(0)+1/3 δ(1), given that one bond per lattice site is occupied. In the inset of Fig. 6 the expected behavior is further evidenced by the computation of the integral of the peaks in small intervals near 0 (solid lines) and 1 (dashed lines) respectively; for increasing system size and disorder strength, the peaks approach 2/3 and 1/3 respectively.

C. Half-system entanglement entropy
Next, we consider the entanglement properties of eigenstates through their von Neumann entanglement entropy where A is a region comprising half of the sample and ρ A = Tr B ρ is the reduced density matrix obtained from an eigenstate by tracing out the complementary region B. The analysis of the entanglement entropy has been especially useful in the study of MBL transitions given the low, area law entanglement of all localized states, to be compared with a volume law scaling in the extended phase 12,20,54,55 . In the clusters taken into consideration, there is some freedom in the choice of the two regions A and B; here, where possible, we consider a cut that runs parallel to the lattice vectors. The two regions are shown in red and blue respectively in Fig. 1 of the main text and in Fig. 14 in the appendix.
In the top panel of Fig. 7 we show the entanglement entropy Eq.(6) as a function of the disorder strength V for different sample sizes. For low-V values, we see that S approaches the value obtained for random states (shown in the figure as a dashed line) as V → 0, thus making evident a volume law entanglement. At high disorder, on the other hand, we observe an area law growth; specifically, by considering S/A where A is the length of the boundary between the two subsections, we observe a collapse (see Fig. 7 bottom panel). Interestingly, as seen for other quantities, the curves for different system sizes collapse in pairs, at around V = 18 for sizes N = 42 and N = 54 and at V = 20 for sizes N = 72 and N = 78, with both sets of curves collapsing only for larger V .
Given the relatively arbitrary choice of the boundary of the bipartition, as an additional comparison and justification for adequateness of the use of volume and A area laws, we considered the entanglement entropy of some special states. One class is the already mentioned random states with volume law entanglement growth. We additionally considered the uniform ('Rokshar-Kivelson' 45 ) state ψ RK = 1/ √ N H |1 1 . . . 1 and the ground state of the model (1) with no disorder and a small constant field V c , which both have an area law entanglement scaling. The entanglement entropy computed for the clusters and the cuts under consideration indeed scales with A as expected (see Appendix C and Fig. 16).
In order to better understand the position of a transition point, we consider the variance of the entanglement entropy distribution as a function of disorder. The vari- ance is expected to have a peak at the transition value (with possibly strong finite-size corrections) 54,55 . In the main panel of Fig. 8 we show the standard deviation σ S of S for the eigenstates in the energy window around E = 0 and for different disorder realizations. A peak is present, although with a substantial drift towards higher disorder values. The position of the peak rescaled with cluster size shows an approximately linear increase with respect to the system size (see bottom panel). Thus, for the entanglement entropy, system sizes up to N = 78 do not show convergence to a finite transition value. This might be an indication that the system sizes that we considered are still within the non-universal scaling regime or that the transition does not hold asymptotically in the thermodynamic limit.

D. Dynamics
We finally consider the dynamical properties of the system. Starting from a product state, which is taken as an element of the computational basis, we perform a quench to the disordered model: The chosen initial state |ψ(0) is the same as the reference state for the imbalance calculation in Sec. III B. In the 1D MBL phase, transport of local quantities is absent and entanglement has a well-understood slow logarithmic growth 22,56,57 . We look for these markers of localization in the present model at high disoder. We consider the same clusters that have been used in the exact diagonalization analysis, that is N = 42, 54, 72 and 78, with the addition of the N = 96 and 108 clusters. The time evolution is performed through full exact diagonalization for the clusters N = 42 and 54, and with the Krylov method for the larger ones. We average over 10 4 ÷ 10 3 disorder realizations for clusters up to N = 96 and around 100 realizations for the largest cluster N = 108. Imbalance We start by considering the imbalance of the time-evolved state with respect to the initial state, as defined in Sec. III B and Eq. (4). In Fig. 9 we show the modulus-squared imbalance as a function of time for various values of the disorder V for the system sizes N = 78 and N = 108. An imbalance value |I| 2 > 0 indicates that some memory of the initial state is kept after the time evolution. From Fig. 9, we see a decrease and, for most disorder values, a saturation of the imbalance (namely for the N = 78 cluster for which longer times are available). For this reason, we look at the asymptotic value, estimated from the last available point of the time evolution, and look at its scaling with the system size. We remark that, while for the smallest clusters N = 42 and N = 54 we are able to obtain the evolved states at very large times, for the larger sizes and with the Krylov time evolution we are only able to reach times of order t = 1000 (in units of the inverse of the plaquette flip energy scale τ ) according to the system size and the disorder strength. An appropriate error is attributed to the points that are not sufficiently close to the saturation value. In Fig. 10, we show the asymptotic values as a function of the disorder highlighting their dependence on the system size. For finite size systems it is expected that |I| 2 > 0 and one should therefore look at the thermodynamic limit. We extrapolate the infinite-size imbalance I 0 (V ) from a scaling function of the form |I| 2 = I 0 +a/N , and we observe (see inset) that it is 0, or reasonably close to it, for V 20, while it increases to non-zero values for V > 20, indicating a localized state where some memory of the initial state is kept at infinite time.
Entanglement entropy A known remarkable feature of the localized phase in one dimension is a slow growth of the entanglement, which spreads logarithmically in time as opposed to a ballistic (linear in time) spread in the extended phase. In finite systems the growth is eventu- ally limited by the corresponding volume law in the two phases [22][23][24][56][57][58] . We remark that given the geometry imposed by the entanglement cut of the 2d system (see Fig. 1b and 14), entanglement can spread only in the direction perpendicular to the cut, and we thus expect a spread similar to a 1D localized phase in this case. In Fig. 11, we show the bipartite entanglement entropy, as defined in Eq. (6), as a function of time, for the cluster sizes N = 78 and N = 108. For low disorder, a fast saturation to the volume law value can be readily observed. As disorder increases, the entanglement entropy continues to quickly reach a size-dependent limiting value. For disorders V 20, a logarithmic growth appears to be present, consistent with the existence of a MBL phase. We note that this feature is only visible in the largest clusters, N = 96 and N = 108, highlighting the need of analysing very large system sizes in order to obtain evidence of a localized phase.
Return probability We then consider the return probability R = | ψ(t)|ψ(0) |. Being an overlap of two vectors in the Hilbert space, one expects that it will be exponentially small (scaling as the inverse of the Hilbert space size) at long times in both ergodic and localized phases, but its time dependence may reveal non-trivial differences. To account for the system size scaling, we consider (minus) the logarithm of the return probability rescaled with the (log of the) Hilbert space size − ln R/ ln N H which is displayed as a function of time in Fig. 12 for six values of the disorder V . For low disorder (V = 1, V = 10 and V = 15, top row in Fig. 12) the rescaled return probability quickly reaches a limiting value, which is smaller in absolute value as the disorder increases. At larger disorder (V 20), a logarithmic increase appears for larger system sizes, indicating a slow spreading in a range consistent with the one obtained from entanglement entropy and the participation entropy results shown below. We finally note that there is reasonable collapse between different system sizes (except the smallest two N = 42 and N = 54).
Participation entropy Finally, we consider the participation entropy, as defined in Eq. (3), of the time evolved state. In Fig. 13 we show the participation entropy, rescaled by the logarithm of the Hilbert space size, as a function of time, for six values of the disorder strength. For the small disorders V = 1, V = 10 and V = 15, shown in the top row, a quick saturation to system-size dependent values can be readily observed, with notably a saturation to a value very close to 1 for very small disorder; for higher, V ≥ 20 disorders, shown in logarithmic scale in the bottom row of Fig. 4, a slow, logarithmic growth suggesting localization becomes apparent for the two largest system sizes. The behavior of the participation entropy thus closely resembles the one of the bipartite entanglement entropy.

IV. CONCLUSIONS
The analysis of the eigenstates and of the dynamics after a quench suggest the same conclusions than the ones reached in a similar disordered QDM on the square lat- tice 48 : the presence of an extended and a many-body localized phase at low and high disorder respectively. This conclusion is justified by the study of a large system size, as large as possible with unbiased numerical methods, up to N = 78 for exact diagonalization and N = 108 for dynamics. The analysis of large sizes is essential in finding some characteristic features of MBL, such as the slow logarithmic growth of entanglement entropy. From our analysis on finite systems at finite times, however, it cannot be excluded that, in the thermodynamic limit, there is no transition but a crossover to increasingly slow dynamics. This is hinted by e.g. the linear scaling with system size of the maximum of the entanglement entropy variance (even though this quantity does not accurately locate the transition, already in the standard model of MBL in 1D 12 ). We remark that in that case the time scales for thermalization at high disorder would likely still be so long for the system to be effectively localized for practical purposes, in particular potential experimental platforms.
We also attempted a scaling analysis (done through the bayesian method 59 ) on some of the quantities presented in Sec. III. With the available system sizes, it was not possible to obtain a collapse. This tends to suggest that, if a true transition exists, the finite systems considered are not large enough to be in the universal scaling regime.
The analysis presented in this work makes use of dimers on a peculiar lattice: in other words, we use a very constrained model in order to numerically study the 2D MBL problem in the largest physical system attainable with the current numerical capabilities. Considering the current lack of theoretical arguments for MBL in 2D, alternative opportunities come from possible experimental realizations in specifically arranged experimental setups. There has been a lot of recent effort devoted to perform analog quantum simulations of lattice gauge theories (see Ref. 60 for a recent review), in order to implement experimentally e.g. the Gauss law equivalent to the dimer constraint. Let us for instance highlight explicit proposals for implementing QDMs with different possible setups using Rydberg atoms [61][62][63] . Finally, it would be interesting to see whether the constraints and the non-tensor product structure in QDM could allow the existence of quantum scar states 64 , similar e.g. to what happens in the 1D constrained PXP model. These scar states have been argued to realize intermediate scenarios between the extended and localized paradigms.  In this work we have used the honeycomb lattices with N = 42, 54, 72, 78, 96 and 108 sites shown in Fig. 14. These were all considered with periodic boundary conditions and are constructed with the following basis vectors 71 , written in the basis {u 1 , u 2 } where u 1 = (1, 0) and The separation into two subsystems used for the calculation of the bipartite entanglement entropy is shown in different colors in each cluster. The boundary has been chosen parallel to one of the basis vectors. We note that in some cases (namely, clusters N = 54 and N = 78) this was not exactly possible but was chosen as close as possible to the parallel boundary line.

Appendix B: Mobility edge
We present here an additional analysis of the gap ratio defined in Sec. III A, this time resolved in energy. The purpose is to identify a possible dependence of the localization transition value from the energy, i.e. the presence of a so-called mobility edge 12 . For the smallest system sizes, we consider full exact diagonalization. As customary, we introduce the parameter Thus, from the whole spectrum, we compute the gap ratio for ten windows of fixed width and average on around 1000 disorder realizations. The result for cluster sizes N = 42 and N = 54 is shown in Fig. 15. Having only the two smallest system sizes available, we cannot definitively conclude the existence of a mobility edge in the model (1), although Fig. 15 does show an indication of an enhanced localization at the spectrum extrema, which appears more marked for N = 54 than N = 42. Given the different symmetries and aspect ratio of the clusters, dividing them in two subsystems for the purpose of the computation of the bipartite entanglement entropy should be done respecting the vectors of each cluster, as outlined in Appendix A. In order to check that the chosen cut is sufficiently general, we computed the entanglement entropy of some reference states which are known to have an area law as the system size increases. The entanglement entropy, rescaled by the area of the cut, is shown in Fig. 16. The reference states are: the ground state |ψ GS of the nondisordered model with constant potential V e = 0.1; the 'Rokshar-Kivelson' 45 state |ψ RK = 1/ √ N H |1 1 . . . 1 ; two localized states at high disorder, respectively obtained at disorder strength V = 30 and V = 50. For all states, S/A is approximately constant with respect to system size N , showing thus the correct area law scaling for the selected cut in all the clusters shown in Fig. 14.