Probing many-body localization in a disordered quantum dimer model on the honeycomb lattice
SSciPost Physics Submission
Probing many-body localization in a disordered quantumdimer model on the honeycomb lattice
Francesca Pietracaprina , Fabien Alet † Laboratoire de Physique Th´eorique, IRSAMC, Universit´e de Toulouse, CNRS, UPS,France School of Physics, Trinity College Dublin, College Green, Dublin 2, Ireland* [email protected] † [email protected] Abstract
We numerically study the possibility of many-body localization transition in adisordered quantum dimer model on the honeycomb lattice. By using the pecu-liar constraints of this model and state-of-the-art exact diagonalization and timeevolution methods, we probe both eigenstates and dynamical properties and con-clude on the existence of a localization transition, on the available time and lengthscales (system sizes of up to N = 108 sites). We critically discuss these resultsand their implications. Contents a r X i v : . [ c ond - m a t . d i s - nn ] J un ciPost Physics Submission Localization in disordered, interacting quantum systems [1, 2] is a topic that has recentlyreceived wide attention due to the very peculiar phenomenology [3–6], the foundational issuesabout quantum integrability and ergodicity involved [7, 8], and the increased precision andcontrol 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 thermal-ization hypothesis (ETH) and one at high disorder which exhibits no transport, no thermal-ization [11–14] and emergent integrability due to an extensive number of quasi-local integralsof motion [15–19]. Furthermore, localized states have low entanglement at any energy andobey an area law, a property usually valid for ground states only [20, 21]. Finally, localizationin interacting systems is characterized by the very slow spreading of information, namely theentanglement [22–24], and the total absence of transport for local observables [1, 2]. All thesefeatures have contributed to make MBL a compelling physical phenomenon, including withrespect to quantum information processing protocols [20, 25–27].In the context of the study of MBL transitions, a wide range of results, outlining thephenomenology described above, have been produced for one-dimensional (1D) systems [3–6].Remarkably, a proof of the existence of the MBL transition has been obtained for a 1Dquantum Ising model with a transverse field [28, 29]. In higher dimensions, however, no suchproof exists. One generally expects that in higher dimensions delocalization is favoured dueto the increase in channels for the delocalizing terms, similarly to the phenomenology ofAnderson localization in higher dimensions. More specifically, general arguments based onthe existence and size-scaling of thermalizing bubbles support the absence of localization forlarge enough times [30, 31], even though no rigorous proof was obtained either.A number of results on 2D systems have notably been presented. Experimental resultsobtained in cold atoms setups interestingly show absence of dynamics and localization at highdisorder [10,32]. At present, this experimental evidence is arguably of higher quality than theanalytical and numerical modeling of MBL in 2D. Numerically, a number of approaches havebeen explored in 2D lattice models, using both unbiased and biased methods, and showingindications of a localized phase [33–38]. Other simulations conclude in favor of absence ofMBL [39]. However, the main limit of numerical approaches is the small system sizes and/ortime scales that are reachable in the computations. The size of the Hilbert space and thusof the quantum problem grows exponentially with the number of particles N in the systemwhile the physical length scale of the sample grows as a square root of N . For unbiasedmethods this is an especially strong constraint, effectively limiting the analysis to systems upto around 20 spins 1/2. While in one dimension several different lattice sizes can fulfill thisrequirement, thus allowing in principle finite size scaling to be performed, this is no longerthe case in two dimensions where the number of system sizes are greatly limited. While largersystem sizes can be reached using methods geared towards capturing properties of an MBLphase [33–35,40–42], these methods are not unbiased and by construction will miss the ergodicphase or the phase transition.Here, we aim to investigate an MBL transition in a specific system up to a real-space sizeas large as possible and with unbiased methods. We do this by considering a highly con-strained model and state-of-the-art numerically exact methods [43]. Specifically we considera disordered quantum dimer model (QDM) on a honeycomb lattice, where each lattice link iseither free or occupied by a dimer with the constraint that each lattice site is touched by one2 ciPost Physics Submission and only one dimer [44–46]. An immediate consequence for this is that the dynamics of such amodel is very constrained: single dimer moves are not allowed and the simplest move involvesan hexagonal plaquette. Moreover, this constraint also automatically encodes strong interac-tions which for the honeycomb lattice already imply long-range correlations in the statisticalensemble of dimer coverings. The interplay between a constrained dynamics, which favorsslow dynamics and localization [47], and the strong interactions, which favor delocalization,creates an ideal situation for an MBL transition to exist. Finally, we note that such modelsare based on Hilbert spaces that, due to the constraints, have considerably lower dimensioncompared to spin systems: for N N , while it scales only as ≈ . N [44] for a dimer system on a N -sites honeycomb lattice, giving an obvious numericaladvantage for large system sizes. A previous work has analyzed a similar disordered QDM ona square lattice [48]. Here, we substantially push forward this analysis, almost doubling themaximum system size reached, by turning to the honeycomb lattice instead.The article is structured as follows. In Section 2 we detail the model Hamiltonian, thesymmetry sectors and the lattices used as well as the procedures used to obtain the numericalresults. Such results are outlined in Section 3, first considering observables within exact mid-spectrum eigenstates, and, secondly, the dynamical properties obtained with Krylov timeevolution. Finally, we provide conclusions in Section 4. In the appendix we discuss in detailthe lattice clusters used in the numerical analysis (Appendix A), further energy-resolvedquantities (Appendix B) and comparisons with the entanglement properties of specific states(Appendix C). We consider the following quantum dimer model on the honeycomb lattice [45, 46, 49] with arandom potential: H QDM = − τ (cid:88) (cid:55) p (cid:16) | (cid:105)(cid:104) | + | (cid:105)(cid:104) | (cid:17) + (cid:88) (cid:55) p v p (cid:16) | (cid:105)(cid:104) | + | (cid:105)(cid:104) | (cid:17) The first term, an hexagon “flip”, is a kinetic term. The second term is a disorderedpotential 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 showthe N = 72 lattice and we refer the reader to the Appendix for more details on the otherclusters. On the honeycomb lattice with periodic boundary conditions, the constraints dueto the dimers and to the allowed plaquette moves are such that two conserved quantities, thewinding numbers, exist. The winding numbers are defined as the sum along a line parallel tothe x or y axis; having labeled the honeycomb lattice sites with binary alternating symbols A (yellow in Fig. 1) and B (green) and orienting all links from A → B , we add a +1 valueto the sum if the line crosses a dimer with an arrow in the positive direction, − w x = w y = 0, which is thelargest one. We remark that, for finite lattices, not all lattice shapes allow the existence ofthis zero winding sector; we discard lattice shapes that do not satisfy this requirement [50].Table 1 displays the number of allowed coverings in the zero winding sector, which cor-respond to the size N H of the Hilbert space. The number of nonzero elements in the matrix3 ciPost Physics Submission A B BB AA ba c
Figure 1: Dimers on the honeycomb lattice. (a) Untilted N = 72 cluster. We consider periodicboundary conditions. The lattice is split in two parts (here shown in red and blue) for thecomputation of the bipartite entanglement. (b) Example of dimer covering, namely the starconfiguration that is used as a reference for the computation of the imbalance (Sec. 3.2) andas initial state for the quenched dynamics (Sec. 3.4). The definition of the phases φ p used todefine the imbalance in Sec. 3.2 are shown for each plaquette. (c) Computation of the twoindependent winding numbers on a single plaquette.is also noted, which, in addition to matrix size, contributes to limiting the feasibility of thenumerical calculations.Cluster size Coverings w x = w y = 0 Nonzero elements N N H nnz
42 1 032 8 04654 7 311 69 51972 131 727 1 596 92778 349 326 4 536 28896 6 460 809 100 676 169108 45 649 431 791 275 167Table 1: Matrix size N H and number of nonzero elements nnz for the clusters that have beenconsidered.We perform exact diagonalization on some of these lattices (up to size 78). We use eitherfull diagonalization or shift-invert methods [43] to obtain around 100 eigenstates at the centerof the spectrum. We also study the dynamics of nonequilibrium initial states though Krylovsubspace time evolution methods for all lattice sizes [51]. In all cases, we average over disorderrealizations of the random potential (at least 1000 for most system sizes and around 100 forthe dynamics on the largest one). 4 ciPost Physics Submission 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 thesystem. N = = = = r V = = = =
20 V = = = P ( s ) Figure 2:
Right panel.
Gap ratio r as a function of the disorder strength V for differenthoneycomb samples. The two limiting values at approximately 0 .
53 and 0 .
39 (dashed lines)correspond to the ones obtained for a Wigner-Dyson and Poisson distribution respectively.
Left panel.
Distribution of the energy gaps s i of the unfolded spectrum. The Wigner-Dysonand Poisson reference distributions are shown in dashed lines. Spectral gap ratio
We start by analyzing the spectral properties of the two phases. Specif-ically, we consider the energy level gap ratio [14]: (cid:104) r (cid:105) = (cid:28) min( s i , s i +1 )max( s i , s i +1 ) (cid:29) i , (1)where s i = E i +1 − E i is the gap between two adjacent eigenvalues. We average in a smallwindow of about 100 eigenstates around the center of the spectrum as well as over disorderrealizations. Depending on the level gap statistics, (cid:104) r (cid:105) ≈ .
39 for a Poisson distribution inthe localized phase and (cid:104) r (cid:105) ≈ .
53 [52] for a Wigner-Dyson distribution corresponding to theETH phase.In Fig. 2, top panel, we show the value of (cid:104) r (cid:105) as a function of the disorder for varioussystem sizes. It appears that both localized and ETH phases are captured with the availablecluster sizes. The transition value can typically be inferred by where the curves for increasingsize cross, as it denotes opposite flows in the system size scaling in the two phases. We notethat 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 theunfolded spectrum for various values of the disorder V , showing excellent agreement with aPoissonian 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 afunction of the energy density (not just for the middle of the spectrum), see Appendix B.5 ciPost Physics Submission Kullback-Leibler divergence for energy-adjacent eigenstates
We now consider quan-tities characterizing eigenstate properties which have been shown to be good indicators oflocalization. In the localized phase, eigenstates and local observables close in energy are verydifferent in structure, as opposed to the ETH phase. Thus, we consider the Kullback-Leiblerdivergence for two consecutive eigenstates | ψ (cid:105) and | ψ (cid:48) (cid:105) in the spectrum, defined asKL = (cid:88) i |(cid:104) ψ | b i (cid:105)| ln |(cid:104) ψ | b i (cid:105)| |(cid:104) ψ (cid:48) | b i (cid:105)| , (2)where the sum runs over the N H elements | b i (cid:105) of the Hilbert space basis. We expect KL toapproach to KL GOE = 2 (the value obtained for the Gaussian orthogonal ensemble of randommatrices) in an ETH phase and to diverge with system size in a localized phase [12]. N = = = = K L Figure 3: Kullback-Leibler divergence KL of eigenstates adjacent in the many-body spectrumas a function of the disorder strength V and for different samples.We show the results for KL in Fig. 3 as a function of the disorder strength V . The limitvalue 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 tothe small size), with some drift due to finite-size scaling, suggesting a localization transitionaround V ≈ − Eigenstate participation entropy
In a similar manner, we consider the participationentropy of the eigenstates, which gives information about localization in the Hilbert space [12,53]. It is defined as S p = − (cid:88) i |(cid:104) ψ | b i (cid:105)| ln |(cid:104) ψ | b i (cid:105)| . (3)For a state which is localized in the Hilbert space, S p is of O (1). For many-body localizedstates, a multifractal behavior is expected in this computational basis [53], with a participation6 ciPost Physics Submission 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 thecoefficient a up to higher order corrections), as a function of the disorder V . At low disorderwe see that a has a high value which is likely to scale to 1 with increasing size. A differentbehavior onsets at around V ≈ −
25: the curves for different system sizes join and collapse,suggesting a finite a < N = = = = S p / Log H Figure 4: Eigenstates participation entropy rescaled by the logarithm of the Hilbert spacesize, S p / ln N H , as a function of disorder strength V and for different samples. V = = = = | ℐ P ( | ℐ ) V = = = | ℐ P ( | ℐ ) Figure 5: Probability distribution of modulus-squared imbalance for system size N = 78.The imbalance of the configuration basis states are shown as dashed lines. Left panel.
Smalldisorders ( V = 5 to 20) show a distribution peaked in 0 for small disorder ( V = 5, 10)and broadening and development of peaks corresponding to the basis vectors imbalance (for V = 20). Right panel.
Large disorders display a clear structure reproducing the dimerconfiguration basis states, while maintaining a continuous distribution.7 ciPost Physics SubmissionEigenstate imbalance
We next consider the imbalance of the eigenstates with respect toa specific configuration where the state used as a reference is chosen as the basis element ofthe so-called star configuration displayed in Fig. 1b. We define the (complex) imbalance as I = 1 N p (cid:88) (cid:55) p e iφ p (cid:16) | (cid:105)(cid:104) | + | (cid:105)(cid:104) | (cid:17) . (4)The phases φ p assume three possible values: 0, π and − π , depending on the dimer config-uration on the plaquette p (see Fig. 1b). With this definition, the imbalance of the referencebasis state in Fig. 1b has (cid:60) I = 1 and maximum amplitude | I | = 1. Delocalized eigen-states will have a probability distribution for the modulus squared imbalance which is sharplypeaked in 0. On the other hand, if states are localized and close to basis states, the imbalancewill 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 twopanels of Fig. 5 respectively for low (top panel) and high (bottom panel) values of disorder forthe largest system size ( N = 78). As expected, the imbalance distribution is sharply peakedat 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| > V ≥ Eigenstate dimer bond occupation
We finally consider a local observable, the dimerbond occupation, and specifically the probability distribution of O k = (cid:104) ψ | n k | ψ (cid:105) (5)where the operator n k acts on the basis vectors | b i (cid:105) as n k | b i (cid:105) = 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 thesame probability to be occupied, i.e. 1 /
3. As shown in the main panel of Fig. 6, for adelocalized eigenstate, this translates into a probability distribution of O sharply peaked at1 / 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. Inthe limit of infinite disorder, where the eigenstates coincide with the configuration states, thedistribution is 2 / δ (0) + 1 / δ (1), given that one bond per lattice site is occupied. In the insetof Fig. 6 the expected behavior is further evidenced by the computation of the integral of thepeaks in small intervals near 0 (solid lines) and 1 (dashed lines) respectively; for increasingsystem size and disorder strength, the peaks approach 2 / / Next, we consider the entanglement properties of eigenstates through their von Neumannentanglement entropy S = − Tr ρ A ln ρ A (6)8 ciPost Physics Submission ● ● ● ● ● ● ●■ ■ ■ ■ ■ ■ ■ ◆ ◆ ◆ ◆ ◆ ◆ ◆ ▲ ▲ ▲ ▲ ▲ ▲ ▲ ● ● ● ● ● ● ●■ ■ ■ ■ ■ ■ ■ ◆ ◆ ◆ ◆ ◆ ◆ ◆ ▲ ▲ ▲ ▲ ▲ ▲ ▲ ● N = ■ N = ◆ N = ▲ N = V = = = = = = = P ( O ) Figure 6: Probability distribution of bond occupation for the system size N = 78. Inset.
Integral of P ( O ) over a small interval ( ≤ O = 0 (solid lines) and O = 1 (dashedlines), compared to the respective limit values 2 / / N = = = = S N = = = = S / Figure 7:
Left panel.
Entanglement entropy S as a function of the disorder V and for differentsample sizes. The dashed lines are the bipartite entanglement entropy for random states,averaged over 10 realizations, and represent the limiting value of the entanglement entropyof eigenstates as V → Right panel.
Entanglement entropy rescaled with the size of theboundary A of the bipartition, showing a collapse, and thus an area-law scaling, at highdisorder values. 9 ciPost Physics Submission (cid:1) (cid:1) (cid:1) (cid:1) (cid:1) (cid:1) (cid:1) (cid:1) (cid:1) (cid:1) (cid:1)(cid:1) (cid:1) (cid:1) (cid:1) (cid:1) (cid:1) (cid:1) (cid:1) (cid:1) (cid:1) (cid:1)(cid:1) (cid:1)(cid:1) (cid:1) (cid:1) (cid:1) (cid:1) (cid:1)(cid:2) (cid:2) (cid:2) (cid:2) (cid:2) (cid:2) (cid:2) (cid:2) (cid:2) (cid:2) (cid:2)(cid:2) (cid:2) (cid:2) (cid:2) (cid:2) (cid:2) (cid:2) (cid:2) (cid:2) (cid:2) (cid:2)(cid:2) (cid:2)(cid:2) (cid:2) (cid:2) (cid:2) (cid:2) (cid:2) (cid:3) (cid:3) (cid:3) (cid:3) (cid:3) (cid:3) (cid:3) (cid:3) (cid:3) (cid:3) (cid:3)(cid:3) (cid:3) (cid:3) (cid:3) (cid:3) (cid:3) (cid:3) (cid:3) (cid:3) (cid:3) (cid:3)(cid:3) (cid:3)(cid:3) (cid:3) (cid:3) (cid:3) (cid:3) (cid:3) (cid:4) (cid:4) (cid:4) (cid:4) (cid:4) (cid:4) (cid:4) (cid:4) (cid:4) (cid:4) (cid:4)(cid:4) (cid:4) (cid:4) (cid:4) (cid:4) (cid:4) (cid:4) (cid:4) (cid:4) (cid:4) (cid:4)(cid:4) (cid:4)(cid:4) (cid:4) (cid:4) (cid:4) (cid:4) (cid:4) (cid:1) N = (cid:2) N = (cid:3) N = (cid:4) N = (cid:5) S (cid:1) (cid:1) (cid:1) (cid:1)
45 50 55 60 65 70 750.150.170.19 N V m a x / N Figure 8:
Top panel.
Standard deviation of the entanglement entropy as a function of disorderfor different sample sizes.
Bottom panel.
Position of the maximum entanglement entropystandard deviation V max scaled by the system size, as a function of the system size. For therange of accessible sizes, V max has an approximately linear increase with system size.where A is a region comprising half of the sample and ρ A = Tr B ρ is the reduced densitymatrix obtained from an eigenstate by tracing out the complementary region B . The analysisof the entanglement entropy has been especially useful in the study of MBL transitions giventhe low, area law entanglement of all localized states, to be compared with a volume lawscaling in the extended phase [12, 20, 54, 55]. In the clusters taken into consideration, thereis some freedom in the choice of the two regions A and B ; here, where possible, we considera cut that runs parallel to the lattice vectors. The two regions are shown in red and bluerespectively 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 thedisorder strength V for different sample sizes. For low- V values, we see that S approachesthe value obtained for random states (shown in the figure as a dashed line) as V →
0, thusmaking evident a volume law entanglement. At high disorder, on the other hand, we observean area law growth; specifically, by considering S/ A where A is the length of the boundarybetween 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 bothsets of curves collapsing only for larger V .Given the relatively arbitrary choice of the boundary of the bipartition, as an additionalcomparison and justification for adequateness of the use of volume and A area laws, weconsidered the entanglement entropy of some special states. One class is the already mentionedrandom states with volume law entanglement growth. We additionally considered the uniform10 ciPost Physics Submission (‘Rokshar-Kivelson’ [45]) state ψ RK = 1 / √N H | . . . (cid:105) and the ground state of the model(1) with no disorder and a small constant field V c , which both have an area law entanglementscaling. The entanglement entropy computed for the clusters and the cuts under considerationindeed 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 varianceof the entanglement entropy distribution as a function of disorder. The variance is expectedto 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 inthe 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 peakrescaled with cluster size shows an approximately linear increase with respect to the systemsize (see bottom panel). Thus, for the entanglement entropy, system sizes up to N = 78 donot show convergence to a finite transition value. This might be an indication that the systemsizes that we considered are still within the non-universal scaling regime or that the transitiondoes not hold asymptotically in the thermodynamic limit. 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 disorderedmodel: | ψ ( t ) (cid:105) = exp( − iHt ) | ψ (0) (cid:105) . (7)The chosen initial state | ψ (0) (cid:105) is the same as the reference state for the imbalance calculationin Sec. 3.2. In the 1D MBL phase, transport of local quantities is absent and entanglementhas a well-understood slow logarithmic growth [22, 56, 57]. We look for these markers oflocalization in the present model at high disoder. We consider the same clusters that havebeen used in the exact diagonalization analysis, that is N = 42, 54, 72 and 78, with theaddition of the N = 96 and 108 clusters. The time evolution is performed through full exactdiagonalization for the clusters N = 42 and 54, and with the Krylov method for the largerones. We average over 10 ÷ disorder realizations for clusters up to N = 96 and around100 realizations for the largest cluster N = 108. V = = =
10 V = = =
25 V = | ℐ N = V = = =
10 V = = =
25 V = | ℐ N = Figure 9: Modulus-squared imbalance | I | as a function of time for the two clusters N = 78(left) and N = 108 (right). 11 ciPost Physics Submission ● ● ● ● ● ● ●■ ■ ■ ■ ■ ■ ■ ◆ ◆ ◆ ◆ ◆ ◆ ◆ ▲ ▲ ▲ ▲ ▲ ▲ ▲▼ ▼ ▼ ▼ ▼ ▼ ▼○ ○ ○ ○ ○ ○ ○ ● ● ● ● ● ● ● - I ● N = ■ N = ◆ N = ▲ N = ▼ N = ○ N = | ℐ ( t m a x ) Figure 10: Asymptotic value of modulus-squared imbalance at long times, as a function ofdisorder.
Inset.
Coefficient I of the scaling function |I| = I + a/N as a function ofdisorder. Imbalance
We start by considering the imbalance of the time-evolved state with respectto the initial state, as defined in Sec. 3.2 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 systemsizes N = 78 and N = 108. An imbalance value | I | > N = 78 cluster for which longertimes are available). For this reason, we look at the asymptotic value, estimated from the lastavailable point of the time evolution, and look at its scaling with the system size. We remarkthat, while for the smallest clusters N = 42 and N = 54 we are able to obtain the evolvedstates at very large times, for the larger sizes and with the Krylov time evolution we are onlyable 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 attributedto 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 theirdependence on the system size. For finite size systems it is expected that | I | > I ( V ) from a scaling function of the form |I| = I + a/N , and we observe (see inset) that itis 0, or reasonably close to it, for V (cid:46)
20, while it increases to non-zero values for
V >
Entanglement entropy
A known remarkable feature of the localized phase in one dimen-sion is a slow growth of the entanglement, which spreads logarithmically in time as opposedto a ballistic (linear in time) spread in the extended phase. In finite systems the growth iseventually limited by the corresponding volume law in the two phases [22–24, 56–58]. We re-mark that given the geometry imposed by the entanglement cut of the 2d system (see Fig. 1b12 ciPost Physics Submission V = = = =
15 V = = = S N = V = = = =
15 V = = = S N = Figure 11: Bipartite entanglement entropy as a function of time for two system sizes, N = 78(top panel) and N = 108 (bottom panel), showing a logarithmic growth at high disorder( V ≥ N = =
54 N = =
78 N = =
100 200 300 400 5000.00.20.40.60.81.0 t - Log R / Log H N = =
54 N = =
78 N = =
100 200 300 400 5000.00.20.40.60.8 t - Log R / Log H N = =
54 N = =
78 N = =
100 200 300 400 5000.00.10.20.30.4 t - Log R / Log H N = =
54 N = =
78 N = = - Log R / Log H N = =
54 N = =
78 N = = - Log R / Log H N = =
54 N = =
78 N = = - Log R / Log H Figure 12: Time dynamics of the return probability of the initial product state, rescaled bythe system size − ln R/ ln N H . Top row.
Values for low disorder: from left to right, V = 1, V = 10 and V = 15. Bottom row.
Values for high disorder, showing a logarithmic decrease:from left to right, V = 20, V = 25 and V = 30.and 14), entanglement can spread only in the direction perpendicular to the cut, and we thusexpect 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 functionof time, for the cluster sizes N = 78 and N = 108. For low disorder, a fast saturation to thevolume law value can be readily observed. As disorder increases, the entanglement entropycontinues to quickly reach a size-dependent limiting value. For disorders V (cid:38)
20, a logarithmicgrowth appears to be present, consistent with the existence of a MBL phase. We note thatthis feature is only visible in the largest clusters, N = 96 and N = 108, highlighting the needof analysing very large system sizes in order to obtain evidence of a localized phase.13 ciPost Physics SubmissionReturn probability We then consider the return probability R = |(cid:104) ψ ( t ) | ψ (0) (cid:105)| . Being anoverlap 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 localizedphases, but its time dependence may reveal non-trivial differences. To account for the systemsize scaling, we consider (minus) the logarithm of the return probability rescaled with the (logof the) Hilbert space size − ln R/ ln N H which is displayed as a function of time in Fig. 12for six values of the disorder V . For low disorder ( V = 1, V = 10 and V = 15, top row inFig. 12) the rescaled return probability quickly reaches a limiting value, which is smaller inabsolute value as the disorder increases. At larger disorder ( V (cid:38) N = 42 and N = 54). N = =
54 N = =
78 N = =
100 200 300 400 5000.00.20.40.60.81.0 t S p / Log H N = =
54 N = =
78 N = =
100 200 300 400 5000.00.20.40.60.81.0 t S p / Log H N = =
54 N = =
78 N = =
100 200 300 400 5000.00.20.40.60.81.0 t S p / Log H N = =
54 N = =
78 N = = S p / Log H N = =
54 N = =
78 N = = S p / Log H N = =
54 N = =
78 N = = S p / Log H Figure 13: Participation entropy rescaled by the size of the Hilbert space S p / ln N H as afunction of time. Top row.
Values for low disorder: from left to right, V = 1, V = 10 and V = 15. Bottom row.
Values for high disorder, showing a logarithmic growth: from left toright, V = 20, V = 25 and V = 30. 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 thelogarithm 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 saturationto system-size dependent values can be readily observed, with notably a saturation to a valuevery close to 1 for very small disorder; for higher, V ≥
20 disorders, shown in logarithmicscale in the bottom row of Fig. 4, a slow, logarithmic growth suggesting localization becomesapparent for the two largest system sizes. The behavior of the participation entropy thusclosely resembles the one of the bipartite entanglement entropy.14 ciPost Physics Submission
The analysis of the eigenstates and of the dynamics after a quench suggest the same con-clusions than the ones reached in a similar disordered QDM on the square lattice [48]: thepresence of an extended and a many-body localized phase at low and high disorder respec-tively. This conclusion is justified by the study of a large system size, as large as possiblewith unbiased numerical methods, up to N = 78 for exact diagonalization and N = 108 fordynamics. The analysis of large sizes is essential in finding some characteristic features ofMBL, such as the slow logarithmic growth of entanglement entropy.From our analysis on finite systems at finite times, however, it cannot be excluded that, inthe 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 entanglemententropy variance (even though this quantity does not accurately locate the transition, alreadyin the standard model of MBL in 1D [12]). We remark that in that case the time scales forthermalization at high disorder would likely still be so long for the system to be effectivelylocalized for practical purposes, in particular potential experimental platforms.We also attempted a scaling analysis (done through the bayesian method [59]) on someof the quantities presented in Sec. 3. With the available system sizes, it was not possible toobtain a collapse. This tends to suggest that, if a true transition exists, the finite systemsconsidered 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 otherwords, we use a very constrained model in order to numerically study the 2D MBL problemin the largest physical system attainable with the current numerical capabilities. Consideringthe current lack of theoretical arguments for MBL in 2D, alternative opportunities come frompossible experimental realizations in specifically arranged experimental setups. There hasbeen a lot of recent effort devoted to perform analog quantum simulations of lattice gaugetheories (see Ref. [60] for a recent review), in order to implement experimentally e.g. theGauss law equivalent to the dimer constraint. Let us for instance highlight explicit proposalsfor implementing QDMs with different possible setups using Rydberg atoms [61–63]. Finally,it would be interesting to see whether the constraints and the non-tensor product structure inQDM could allow the existence of quantum scar states [64], similar e.g. to what happens inthe 1D constrained PXP model. These scar states have been argued to realize intermediatescenarios between the extended and localized paradigms. Acknowledgments
This work benefited from the support of the project THERMOLOC ANR-16-CE30-0023-02 ofthe French National Research Agency (ANR) and by the French Programme Investissementsd’Avenir under the program ANR-11-IDEX-0002-02, reference ANR-10-LABX-0037-NEXT.This project has received funding from the European Union’s Horizon 2020 research andinnovation programme under the Marie Sklodowska-Curie grant agreement No 838773. Weacknowledge PRACE for awarding access to HLRS’s Hazel Hen computer based in Stuttgart,Germany under grant number 2016153659, as well as the use of HPC resources from CALMIP(grants 2017-P0677 and 2018-P0677) and GENCI (grant x2018050225). The computer codesthat allowed to obtain the results presented in Sec. 3 make use of the libraries PETSc [65,66],15 ciPost Physics Submission
SLEPc [67, 68] and Strumpack [69, 70].
References [1] D. M. Basko, I. L. Aleiner and B. L. Altshuler,
Metal–insulator transition in a weaklyinteracting many-electron system with localized single-particle states , Ann. Phys. ,1126 (2006), doi:10.1016/j.aop.2005.11.014.[2] I. V. Gornyi, A. D. Mirlin and D. G. Polyakov,
Interacting Electrons in Disordered Wires:Anderson Localization and Low- $ T $ Transport , Phys. Rev. Lett. (20), 206603 (2005),doi:10.1103/PhysRevLett.95.206603.[3] R. Nandkishore and D. A. Huse, Many-Body Localization and Thermalization in Quan-tum Statistical Mechanics , Annual Review of Condensed Matter Physics (1), 15 (2015),doi:10.1146/annurev-conmatphys-031214-014726.[4] F. Alet and N. Laflorencie, Many-body localization: An introduction and selected topics ,Comptes Rendus Physique (6), 498 (2018), doi:10.1016/j.crhy.2018.03.003.[5] D. A. Abanin and Z. Papi´c, Recent progress in many-body localization , Annalen derPhysik (7), 1700169 (2017), doi:10.1002/andp.201700169.[6] D. A. Abanin, E. Altman, I. Bloch and M. Serbyn,
Many-body localiza-tion, thermalization, and entanglement , Rev. Mod. Phys. (2), 021001 (2019),doi:10.1103/RevModPhys.91.021001.[7] C. Gogolin and J. Eisert, Equilibration, thermalisation, and the emergence of statisticalmechanics in closed quantum systems , Reports on Progress in Physics (5), 056001(2016), doi:10.1088/0034-4885/79/5/056001.[8] F. Pietracaprina, C. Gogolin and J. Goold, Total correlations of the diagonal ensemble asa generic indicator for ergodicity breaking in quantum systems , Phys. Rev. B , 125118(2017), doi:10.1103/PhysRevB.95.125118.[9] M. Schreiber, S. S. Hodgman, P. Bordia, H. P. L¨uschen, M. H. Fischer, R. Vosk,E. Altman, U. Schneider and I. Bloch, Observation of many-body localization of in-teracting fermions in a quasirandom optical lattice , Science (6250), 842 (2015),doi:10.1126/science.aaa7432.[10] J.-y. Choi, S. Hild, J. Zeiher, P. Schauß, A. Rubio-Abadal, T. Yefsah, V. Khemani, D. A.Huse, I. Bloch and C. Gross,
Exploring the many-body localization transition in twodimensions , Science (6293), 1547 (2016), doi:10.1126/science.aaf8834.[11] A. Pal and D. A. Huse,
Many-body localization phase transition , Phys. Rev. B (17),174411 (2010), doi:10.1103/PhysRevB.82.174411.[12] D. J. Luitz, N. Laflorencie and F. Alet, Many-body localization edge in the random-fieldHeisenberg chain , Phys. Rev. B , 081103 (2015), doi:10.1103/PhysRevB.91.081103.16 ciPost Physics Submission [13] A. De Luca and A. Scardicchio, Ergodicity breaking in a model showing many-bodylocalization , Europhys. Lett. , 37003 (2013), doi:10.1209/0295-5075/101/37003.[14] V. Oganesyan and D. A. Huse,
Localization of interacting fermions at high temperature ,Phys. Rev. B , 155111 (2007), doi:10.1103/PhysRevB.75.155111.[15] D. A. Huse, R. Nandkishore and V. Oganesyan, Phenomenology of certain many-body-localized systems , Phys. Rev. B (17), 174202 (2014), doi:10.1103/PhysRevB.90.174202.[16] J. Z. Imbrie, V. Ros and A. Scardicchio, Local integrals of motion in many-body localizedsystems , Annalen der Physik (2017), doi:10.1002/andp.201600278.[17] V. Ros, M. M¨uller and A. Scardicchio,
Integrals of motion in the many-body localizedphase , Nucl. Phys. B , 420 (2015), doi:10.1016/j.nuclphysb.2014.12.014.[18] M. Serbyn, Z. Papi´c and D. A. Abanin,
Local conservation laws and the struc-ture of the many-body localized states , Phys. Rev. Lett. (12), 127201 (2013),doi:10.1103/PhysRevLett.111.127201.[19] A. Chandran, I. H. Kim, G. Vidal and D. A. Abanin,
Constructing local integralsof motion in the many-body localized phase , Phys. Rev. B (8), 085425 (2015),doi:10.1103/PhysRevB.91.085425.[20] B. Bauer and C. Nayak, Area laws in a many-body localized state and its implica-tions for topological order , J. Stat. Mech. , P09005 (2013), doi:10.1088/1742-5468/2013/09/P09005.[21] M. Friesdorf, A. H. Werner, W. Brown, V. B. Scholz and J. Eisert,
Many-body localizationimplies that eigenvectors are matrix-product states , Phys. Rev. Lett. , 170505 (2015),doi:10.1103/PhysRevLett.114.170505.[22] M. ˇZnidariˇc, T. Prosen and P. Prelovˇsek,
Many-body localization in the Heisen-berg XXZ magnet in a random field , Phys. Rev. B , 064426 (2008),doi:10.1103/PhysRevB.77.064426.[23] M. Serbyn, Z. Papi´c and D. A. Abanin, Universal slow growth of entanglementin interacting strongly disordered systems , Phys. Rev. Lett. , 260601 (2013),doi:10.1103/PhysRevLett.110.260601.[24] A. Nanduri, H. Kim and D. A. Huse,
Entanglement spreading in a many-body localizedsystem , Phys. Rev. B , 064201 (2014), doi:10.1103/PhysRevB.90.064201.[25] S. Choi, N. Y. Yao, S. Gopalakrishnan and M. D. Lukin, Quantum control of many-bodylocalized states (2015), .[26] N. Y. Yao, C. R. Laumann and A. Vishwanath,
Many-body localization protected quantumstate transfer (2015), .[27] E. A. A. V. Yasaman Bahri, Ronen Vosk,
Localization and topology protectedquantum coherence at the edge of hot matter , Nat. Commun. , 7341 (2015),doi:10.1038/ncomms8341. 17 ciPost Physics Submission [28] J. Imbrie, On many-body localization for quantum spin chains , J Stat Phys , 998(2016), doi:10.1007/s10955-016-1508-x.[29] J. Z. Imbrie,
Diagonalization and many-body localization for a disordered quantum spinchain , Phys. Rev. Lett. , 027201 (2016), doi:10.1103/PhysRevLett.117.027201.[30] W. De Roeck and F. Huveneers,
Stability and instability towards delocal-ization in many-body localization systems , Phys. Rev. B , 155129 (2017),doi:10.1103/PhysRevB.95.155129.[31] I.-D. Potirniche, S. Banerjee and E. Altman, Exploration of the stability of many-bodylocalization in d >
1, Phys. Rev. B , 205149 (2019), doi:10.1103/PhysRevB.99.205149.[32] P. Bordia, H. L¨uschen, S. Scherg, S. Gopalakrishnan, M. Knap, U. Schneider and I. Bloch, Probing slow relaxation and many-body localization in two-dimensional quasiperiodic sys-tems , Phys. Rev. X , 041047 (2017), doi:10.1103/PhysRevX.7.041047.[33] P. A. Wahl, T.B. and S. Simon, Signatures of the many-body localized regime in twodimensions , Nature Phys , 164 (2019), doi:10.1038/s41567-018-0339-x.[34] S. J. Thomson and M. Schir´o, Time evolution of many-body localized systems with the flowequation approach , Phys. Rev. B , 060201 (2018), doi:10.1103/PhysRevB.97.060201.[35] G. De Tomasi, F. Pollmann and M. Heyl, Efficiently solving the dynamics ofmany-body localized systems at strong disorder , Phys. Rev. B , 241114 (2019),doi:10.1103/PhysRevB.99.241114.[36] S. D. Geraedts, R. Nandkishore and N. Regnault, Many-body localization and thermal-ization: Insights from the entanglement spectrum , Phys. Rev. B , 174202 (2016),doi:10.1103/PhysRevB.93.174202.[37] E. van Nieuwenburg, Y. Baum and G. Refael, From bloch oscillations to many-body localization in clean interacting systems , PNAS (19), 9269 (2019),doi:10.1073/pnas.1819316116.[38] S. Inglis and L. Pollet,
Accessing many-body localized states through the generalized gibbsensemble , Phys. Rev. Lett. , 120402 (2016), doi:10.1103/PhysRevLett.117.120402.[39] E. V. H. Doggen, I. V. Gornyi, A. D. Mirlin and D. G. Polyakov,
Slow many-bodydelocalization beyond one dimension (2020), .[40] A. Kshetrimayum, M. Goihl and J. Eisert,
Time evolution of many-body localized systemsin two spatial dimensions (2019), .[41] D. M. Kennes,
Many-body localization in two dimensions from projected entangled-pairstates (2018), .[42] Y. Bar Lev and D. R. Reichman,
Slow dynamics in a two-dimensional anderson-hubbardmodel , Europhys. Lett. (4), 46001 (2016), doi:10.1209/0295-5075/113/46001.[43] F. Pietracaprina, N. Mac´e, D. J. Luitz and F. Alet,
Shift-invert diagonaliza-tion of large many-body localizing spin chains , SciPost Physics (5), 045 (2018),doi:10.21468/SciPostPhys.5.5.045. 18 ciPost Physics Submission [44] P. W. Kasteleyn, Dimer Statistics and Phase Transitions , Journal of MathematicalPhysics (2), 287 (1963), doi:10.1063/1.1703953.[45] D. S. Rokhsar and S. A. Kivelson, Superconductivity and the quantum hard-core dimergas , Phys. Rev. Lett. , 2376 (1988), doi:10.1103/PhysRevLett.61.2376.[46] R. Moessner, S. L. Sondhi and P. Chandra, Phase diagram of the hexagonal lattice quan-tum dimer model , Phys. Rev. B , 144416 (2001), doi:10.1103/PhysRevB.64.144416.[47] J. Feldmeier, F. Pollmann and M. Knap, Emergent glassy dynamics in a quantum dimermodel , Phys. Rev. Lett. , 040601 (2019), doi:10.1103/PhysRevLett.123.040601.[48] H. Th´eveniaut, Z. Lan and F. Alet,
Many-body localization transition in a two-dimensional disordered quantum dimer model (2019), .[49] T. Schlittler, T. Barthel, G. Misguich, J. Vidal and R. Mosseri,
Phase diagram of anextended quantum dimer model on the hexagonal lattice , Phys. Rev. Lett. , 217202(2015), doi:10.1103/PhysRevLett.115.217202.[50] O. C´epas,
Colorings of odd or even chirality on hexagonal lattices , Phys. Rev. B ,064405 (2017), doi:10.1103/PhysRevB.95.064405.[51] A. Nauts and R. E. Wyatt, New approach to many-state quantum dynam-ics: The recursive-residue-generation method , Phys. Rev. Lett. , 2238 (1983),doi:10.1103/PhysRevLett.51.2238.[52] Y. Y. Atas, E. Bogomolny, O. Giraud and G. Roux, Distribution of the ratio of consec-utive level spacings in random matrix ensembles , Phys. Rev. Lett. , 084101 (2013),doi:10.1103/PhysRevLett.110.084101.[53] N. Mac´e, F. Alet and N. Laflorencie,
Multifractal scalings across themany-body localization transition , Phys. Rev. Lett. , 180601 (2019),doi:10.1103/PhysRevLett.123.180601.[54] J. A. Kj¨all, J. H. Bardarson and F. Pollmann,
Many-body localization ina disordered quantum ising chain , Phys. Rev. Lett. , 107204 (2014),doi:10.1103/PhysRevLett.113.107204.[55] V. Khemani, S. P. Lim, D. N. Sheng and D. A. Huse,
Critical properties of the many-bodylocalization transition , Phys. Rev. X , 021013 (2017), doi:10.1103/PhysRevX.7.021013.[56] J. H. Bardarson, F. Pollmann and J. E. Moore, Unbounded Growth of Entangle-ment in Models of Many-Body Localization , Phys. Rev. Lett. , 017202 (2012),doi:10.1103/PhysRevLett.109.017202.[57] D. J. Luitz, N. Laflorencie and F. Alet,
Extended slow dynamical regimeclose to the many-body localization transition , Phys. Rev. B , 060201 (2016),doi:10.1103/PhysRevB.93.060201.[58] H. Kim and D. A. Huse, Ballistic spreading of entanglement in a diffusive nonintegrablesystem , Phys. Rev. Lett. , 127205 (2013), doi:10.1103/PhysRevLett.111.127205.19 ciPost Physics Submission [59] K. Harada,
Bayesian inference in the scaling analysis of critical phenomena , Phys. Rev.E , 056704 (2011), doi:10.1103/PhysRevE.84.056704.[60] M. C. Ba˜nuls, R. Blatt, J. Catani, A. Celi, J. I. Cirac, M. Dalmonte, L. Fallani, K. Jansen,M. Lewenstein, S. Montangero, C. A. Muschik, B. Reznik et al. , Simulating lattice gaugetheories within quantum technologies (2019), .[61] A. W. Glaetzle, M. Dalmonte, R. Nath, I. Rousochatzakis, R. Moessner and P. Zoller,
Quantum spin-ice and dimer models with rydberg atoms , Phys. Rev. X , 041037 (2014),doi:10.1103/PhysRevX.4.041037.[62] A. W. Glaetzle, M. Dalmonte, R. Nath, C. Gross, I. Bloch and P. Zoller, Designingfrustrated quantum magnets with laser-dressed rydberg atoms , Phys. Rev. Lett. ,173002 (2015), doi:10.1103/PhysRevLett.114.173002.[63] A. Celi, B. Vermersch, O. Viyuela, H. Pichler, M. D. Lukin and P. Zoller,
Emerging 2dgauge theories in rydberg configurable arrays (2019), .[64] C. J. Turner, A. A. Michailidis, D. A. Abanin, M. Serbyn and Z. Papi´c,
Weak er-godicity breaking from quantum many-body scars , Nature Physics (7), 745 (2018),doi:10.1038/s41567-018-0137-5.[65] S. Balay, W. D. Gropp, L. C. McInnes and B. F. Smith, Efficient management ofparallelism in object oriented numerical software libraries , In E. Arge, A. M. Bruasetand H. P. Langtangen, eds.,
Modern Software Tools in Scientific Computing , pp. 163–202. Birkh¨auser Press, doi:10.1007/978-1-4612-1986-6 8 (1997).[66] S. Balay, S. Abhyankar, M. F. Adams, J. Brown, P. Brune, K. Buschelman, L. Dalcin,V. Eijkhout, W. D. Gropp, D. Kaushik, M. G. Knepley, L. C. McInnes et al. , PETSc usersmanual , Tech. Rep. ANL-95/11 - Revision 3.8, Argonne National Laboratory (2017).[67] V. Hernandez, J. E. Roman and V. Vidal,
SLEPc: A scalable and flexible toolkit forthe solution of eigenvalue problems , ACM Trans. Math. Software (3), 351 (2005),doi:10.1145/1089014.1089019.[68] J. E. Roman, C. Campos, E. Romero and A. Tomas, SLEPc users manual , Tech.Rep. DSIC-II/24/02 - Revision 3.8, D. Sistemes Inform`atics i Computaci´o, UniversitatPolit`ecnica de Val`encia (2017).[69] P. Ghysels, X. Li, F. Rouet, S. Williams and A. Napov,
An efficient multicore implemen-tation of a novel HSS-structured multifrontal solver using randomized sampling , SIAMJ. Sci. Comput. (5), S358 (2016), doi:10.1137/15M1010117.[70] P. Ghysels, X. S. Li, C. Gorman and F. H. Rouet, A robust parallel preconditionerfor indefinite systems using hierarchical matrices and randomized sampling , In , pp. 897–906, doi:10.1109/IPDPS.2017.21 (2017).[71] M. D. Schulz, S. Dusuel, G. Misguich, K. P. Schmidt and J. Vidal,
Ising anyons with astring tension , Phys. Rev. B , 201103 (2014), doi:10.1103/PhysRevB.89.201103.20 ciPost Physics Submission AppendixA Overview of the finite-size lattices used
Figure 14: Overview of the honeycomb lattice clusters, with periodic boundary conditions,that have been used for exact diagonalization and dynamics. The red and blue subsystemscorrespond to the ones used for entanglement entropy computations.In this work we have used the honeycomb lattices with N = 42, 54, 72, 78, 96 and 108sites shown in Fig. 14. These were all considered with periodic boundary conditions and areconstructed with the following basis vectors [71], written in the basis { u , u } where u = (1 , u = (1 / , √ / N v v
42 (1 4) (5 − − − − −
8) (4 − − { u , u } where u = (1 ,
0) and u = (1 / , √ / N = 54 and N = 78)this was not exactly possible but was chosen as close as possible to the parallel boundary line. B Mobility edge
We present here an additional analysis of the gap ratio defined in Sec. 3.1, this time resolvedin energy. The purpose is to identify a possible dependence of the localization transition valuefrom the energy, i.e. the presence of a so-called mobility edge [12].21 ciPost Physics Submission ϵ V N = ϵ V N = Figure 15: Color map plot of the gap ratio r as a function of the disorder strength V and therescaled energy (cid:15) for the N = 42 (left panel) and N = 54 (right panel) cluster sizes, showingan energy-dependent mobility edge.For the smallest system sizes, we consider full exact diagonalization. As customary, weintroduce the parameter (cid:15) = E − E min E max − E min ∈ [0 , . (8)Thus, from the whole spectrum, we compute the gap ratio for ten (cid:15) windows of fixed widthand 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 cannotdefinitively conclude the existence of a mobility edge in the model (1), although Fig. 15 doesshow an indication of an enhanced localization at the spectrum extrema, which appears moremarked for N = 54 than N = 42. C Area law entanglement of selected states
Given the different symmetries and aspect ratio of the clusters, dividing them in two sub-systems for the purpose of the computation of the bipartite entanglement entropy should bedone respecting the vectors of each cluster, as outlined in Appendix A. In order to checkthat the chosen cut is sufficiently general, we computed the entanglement entropy of somereference states which are known to have an area law as the system size increases. The en-tanglement entropy, rescaled by the area of the cut, is shown in Fig. 16. The reference statesare: the ground state | ψ GS (cid:105) of the nondisordered model with constant potential V e = 0 . | ψ RK (cid:105) = 1 / √N H | . . . (cid:105) ; two localized states at highdisorder, 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 lawscaling for the selected cut in all the clusters shown in Fig. 14.22 ciPost Physics Submission ● ● ● ●■ ■ ■ ■ ◆ ◆ ◆ ◆ ◆ ▲ ▲ ▲ ▲ ● V = ■ V = ◆ | ψ rk > ▲ | ψ GS >
50 60 70 80 900.00.10.20.30.4 N S / A Figure 16: Bipartite entanglement entropy of some reference states in the honeycomb latticerescaled by the size of the boundary between the two subsystems, S/ A , showing an area lawscaling. Shown are entanglement scalings for the ground state | ψ GS (cid:105) of the nondisorderedmodel with constant field V e = 0 .