Polynomially filtered exact diagonalization approach to many-body localization
PPolynomially filtered exact diagonalization approach to many-body localization
Piotr Sierant,
1, 2
Maciej Lewenstein,
2, 3 and Jakub Zakrzewski
1, 4 Institute of Theoretical Physics, Jagiellonian University in Krakow, Łojasiewicza 11, 30-348 Kraków, Poland ∗ ICFO- Institut de Sciences Fotoniques, The Barcelona Institute of Science and Technology,Av. Carl Friedrich Gauss 3, 08860 Castelldefels (Barcelona), Spain ICREA, Pg. Lluís Companys 23, 08010 Barcelona, Spain † Mark Kac Complex Systems Research Center, Jagiellonian University in Krakow, Łojasiewicza 11, 30-348 Kraków, Poland. ‡ (Dated: May 20, 2020)Polynomially filtered exact diagonalization method (POLFED) for large sparse matrices is intro-duced. The algorithm finds an optimal basis of a subspace spanned by eigenvectors with eigenvaluesclose to a specified energy target by a spectral transformation using a high order polynomial ofthe matrix. The memory requirements scale better with system size than in the state-of-the-artshift-invert approach. The potential of POLFED is demonstrated examining many-body localiza-tion transition in 1D interacting quantum spin-1/2 chains. We investigate the disorder strength andsystem size scaling of Thouless time. System size dependence of bipartite entanglement entropyand of the gap ratio highlights the importance of finite-size effects in the system. We discuss pos-sible scenarios regarding the many-body localization transition obtaining estimates for the criticaldisorder strength. Introduction.
Qantum many-body systems are gener-ically expected to approach equilibrium according toeigenstate thermalization hypothesis [1–3]. The phe-nomenon of many-body localization (MBL) [4–6] pro-vides a robust class of many-body systems which fail toreach thermal equilibrium [7–16]. Further examples ofnon-ergodic behavior include Stark localization [17, 18],persistent oscillations [19–23], the presence of confine-ment [24, 25], Hilbert space fragmentation [26–28] or lackof thermalization in lattice gauge theories [29–34].Classification of many-body systems according to theirergodic properties is a fascinating new direction of re-search, however, it poses serious technical challenges asexact methods are restricted either to small system sizes[35] or allow to trace time evolution only within a shorttime interval [36–38]. Hence, a fully consistent theory ofMBL transition is missing, with recent approaches point-ing towards Kosterlitz-Thouless scaling [39–43]. Thefinite-size effect strongly influence exact diagonalization(ED) results, leading to a recent debate [44–47] aboutdiscriminating between finite size effects and asymptoticfeatures of disordered many-body systems.The example of MBL transition shows that develop-ment of ED techniques allowing to study thermalizationproperties of possibly large many-body systems is in de-mand. In this letter, we introduce a polynomially filteredexact diagonalization (POLFED) as a tool to calculateeigenvectors with eigenvalues close to a specified energytarget of large sparse matrices. The polynomial spectraltransformation preserves the sparse structure of matricesavoiding the main bottleneck of shift-invert method of ex-act diagonalization (SIMED) [35]. We employ POLFEDin study of MBL transition in disordered quantum spinchains unveiling new aspects of system size scaling ofThouless time, entanglement entropy and level statistics.Our results provide novel qualitative and quantitative ar- − . − . . . . (cid:15) − R ( (cid:15) ) a) − . − . . . . (cid:15) . . . P ( (cid:15) ) b) Figure 1. Spectral transformation employed in a) SIMED;b) POLFED algorithm. The spectrum is transformed accord-ing to a) R ( (cid:15) ); b) P K =22 σ =0 ( (cid:15) ). Eigenvectors corresponding toeigenvalues at the edges of the transformed spectrum (shadedareas) are accessible for iterative methods. guments in favor of the existence of MBL transition inthe thermodynamic limit. Benchmark models.
We consider 1D disordered spinchains with Hamiltonian:ˆ H = X l =1 L X i =1 J l (cid:0) S xi S xi + l + S yi S yi + l + ∆ S zi S zi + l (cid:1) + L X i =1 h i S zi , (1)where ~S i are spin-1/2 matrices, L is the system size, J = 1 is fixed as the energy unit, periodic boundaryconditions are assumed and h i ∈ [ − W, W ] are indepen-dent, uniformly distributed random variables. The XXZmodel, widely studied in the MBL context [48–53], is ob-tained for J = 0 and ∆ = 1. The choice J = 1 and∆ = 0 .
55 leads to the J - J model studied in [44]. TheHamiltonian (1) becomes a real symmetric sparse ma-trix H ∈ R N ×N in basis of eigenstates of S iz operator;the matrix size, N , in the zero magnetization P i S zi = 0sector is given by N = (cid:0) LL/ (cid:1) ∝ e L ln 2 / √ L . a r X i v : . [ c ond - m a t . d i s - nn ] M a y Calculation of eigenpairs.
Hamiltonians of many-bodysystems are typically characterized by exponential scal-ing of matrix size, N , with the system size, L , and spar-sity in appropriately chosen basis. For a sparse matrixthe number of non-zero entries, N nz , is much smallerthan N implying that matrix vector multiplication re-quires much less operations than for a dense matrix. TheLanczos algorithm [54] utilizes this fact to find exterioreigenpairs (corresponding to highest/lowest eigenvalues).However, due to an increasing density of states and re-orthogonalization costs, Lanczos algorithm becomes in-efficient if many eigenpairs are requested. In contrast, afull ED procedure [55] allows one to determine all eigen-pairs of H but, with present day computers, it is limitedto N (cid:46) · corresponding to L = 18 in (1). Largermatrix sizes are tractable by SIMED [35]. The Hamil-tonian is transformed via H → R σ ( H ) = ( σ − H ) − sothat eigenvalues close to σ become exterior eigenvalues ofthe matrix R σ ( H ), see Fig. 1. Consequently, the Laczosalgorithm for the matrix R σ ( H ) converges to eigenpairsclose to the target σ . The Lanczos iteration with R σ ( H )is performed by calculating LU decomposition [56, 57]of the matrix H . That has a significant drawback: thesparsity pattern of H is lost resulting in a very severe forlarge N phenomenon of fill-in of the matrix. This wasidentified as the main bottleneck of SIMED when appliedto quantum many-body systems [35]. POLFED algorithm.
To avoid the fill-in phenomenonand utilize the sparsity of the H matrix in an efficientway, we use the polynomial spectral transformation H → P Kσ ( H ) = 1 D K X n =0 c σn T n ( H ) (2)where T n ( x ) denotes n -th Chebyshev polynomial, the co-efficients c σn = p − δ ,n cos( n arccos σ ) are obtainedfrom expanding a Dirac delta function centered at σ in Chebyshev polynomials and normalization D assuresthat P σ ( σ ) = 1. The eigenvalues close to the target en-ergy σ are the largest eigenvalues of the transformed ma-trix P σ ( H ) as shown in Fig. 1b). Hence, a block Lanczosmethod [58, 59] applied to matrix P Kσ ( H ) converges toeigenpairs close to the target σ . We note that eigen-solvers employing polynomial spectral transformationswere considered also in [60–63].The POLFED consists of the following steps. Lanc-zos algorithm is used to find the lowest (highest) eigen-value E ( E ) of matrix H which is then rescaled to˜ H = [2 H − ( E + E )] / ( E − E ). The order K of trans-formation (2) is specified by requiring that the numberof eigenvalues θ i of P Kσ ( ˜ H ) accessible to Lanczos algo-rithm (belonging to the shaded area in Fig. 1) is equal toa number of requested eigenvalues N ev – as the conditionwe take θ i (cid:62) p = 0 .
17. To find the value of K , an esti-mate of density of states ˜ ρ ( σ ) at energy σ of the matrix˜ H is needed. The ˜ ρ ( σ ) can be found efficiently for ar- bitrary sparse matrices using iterative methods [64, 65].For the benchmark models (1), the density of states isGaussian and is well approximated by an analytic ex-pression ˜ ρ (0) = ( E − E ) N / Γ at the center of spec-trum σ = 0 where Γ ∝ √ LW . Having found K , thePOLFED algorithm, starting with a matrix of orthonor-malized random vectors Q ∈ R N × s , performs the blockLanczos iteration U j = P Kσ ( ˜ H ) Q j − Q j − B Tj , A j = Q Tj U j (3) R j +1 = U j − Q j A j , Q j +1 B j +1 = R j +1 , (4)where Q = 0, B = 0 and the second operation in(4) is QR decomposition. The iteration is repeatedfor j = 1 , . . . , m resulting in Q j , U j , R j ∈ R N × s and A j , B j ∈ R s × s matrices. In exact arithmetic, columnsof Q j matrices form an orthonormal set of vectors. Thisproperty is gradually lost with increasing m during cal-culations with a finite precision. Hence, between (3) and(4), we perform a re-orthogonalization of columns of ma-trix U j against the columns of matrices { Q i } ji =1 . Theproduct of P Kσ ( ˜ H ) with each column of Q j in (3) is com-puted with the Clenshaw algorithm [66]. The orthogonalmatrix Q m = [ Q , . . . , Q m ] ∈ R N × ms defines a blocktridiagonal matrix T m = Q Tm P Kσ ( ˜ H ) Q m with A j matri-ces on the diagonal and B j ( B Tj ) below (above) the di-agonal. The eigenvectors t i ∈ R ms of T m are used tocalculate u i = Q m t i which converge, with increasing m ,to exterior eigenvectors of P Kσ ( ˜ H ) [67], that is to eigen-vectors of ˜ H with eigenvalues close to the target σ . Theconvergence is reached after m steps when the residualnorm || B j +1 ˜ t i || [59] (where ˜ t i are the last s componentsof the vector t i and || u || = √ u T u ) vanishes within thenumerical precision for each eigenvector t i correspondingto eigenvalue θ i (cid:62) p . The eigenvalues of the matrix ˜ H arefound as ε i = u Ti ˜ Hu i and the convergence is verified bya direct calculation of the residual norms || ˜ Hu i − ε i u i || .Each of our tests shows that eigenvalues ε i are, withinnumerical precision, equal to N ev eigenvalues of ˜ H clos-est to the target σ . For technical details of the algorithmsee [68].The POLFED is tailored for maximal efficiency in cal-culations for quantum many-body systems. The order K of the polynomial transformation (2) scales linearlywith the density of states ˜ ρ ( (cid:15) ) that increases exponen-tially with system size L . Thus, the product P Kσ ( ˜ H ) Q j in (3) is the most time consuming step of the calculation.POLFED offers high scalability as the product can beparallelized in two manners: i) it splits into independentmultiplications of subsequent columns of Q j by P Kσ ( ˜ H );ii) each of the matrix vector multiplications can be par-allelized. The re-orthogonalization step between (3) and(4) can be parallelized in a similar manner. The num-ber m of iterations after which the algorithm convergesis proportional to N ev . Hence, the memory consump-tion, dominated by Q m , scales as N ev N . The memory Figure 2. Thouless time t Th for system size L and disorderstrength W for J - J model. The dotted lines denote Heisen-berg time t H ; the dashed line denotes a scaling t Th ∝ L e W/ Ω broken by the L = 22 ,
24 data. requirements of SIMED are larger and scale as c ( L ) N where the factor c ( L ) is due to the fill-in of the ma-trix. For XXZ model c ( L ) ∝ L/ [35]. Moreover, c ( L )grows rapidly with number N nz of non-zero elements ofthe matrix significantly increasing the resources neededin calculations for J - J model. In contrast, computationtime of POLFED increases linearly with N nz – resourcesfor XXZ and J - J models are comparable, for detailedbenchmarks see [68]. POLFED allows to find larger num-ber of eigenpairs in a single run than the recently pro-posed eigensolver [69]. This reduces fluctuations of aver-ages over eigenstates and is essential in calculation of theThouless time. Thouless time.
The spectral form factor is defined as K ( τ ) = h| P N j =1 g ( E j )e − iE j τ | i /Z , where E j are eigen-values of H after an unfolding procedure [70], g ( (cid:15) ) is aGaussian function, the average is taken over disorder re-alizations. The spectral form factor of many-body system(with time reversal invariance) follows Gaussian Orthog-onal Ensemble (GOE) prediction K ( τ ) = K GOE ( τ ) onlyfor τ > τ T h defining the Thouless time t T h = τ T h t H ,where t H = 2 πρ (0) is the Heisenberg time.The Thouless time, t T h , calculated for J - J spinchains of length L (cid:54)
18 [44] scales as t T h ∝ L e W/ Ω where W is the disorder strength and Ω is constant. Ifthis scaling prevailed in L → ∞ limit, it would imply t T h /t H → J - J model with Thouless times obtained withPOLFED for L = 20 , ,
24 respectively for 800, 200,50 disorder realizations. Since we calculate N ev = 2500eigenvalues in the middle of spectrum ( σ = 0), the sum inthe definition of spectral form factor K ( τ ) is truncated.However, this does not influence the value of t T h as longas it is larger than a certain threshold value determinedby N ev [68]. The obtained Thouless times are shown inFig. 2. Data for L (cid:54)
20 follows the scaling t T h ∝ L e W/ Ω deviating from it at disorder strength ˜ W ( L ) which in- creases with the system size, for instance ˜ W (18) ≈ . W (20) ≈ .
6. This behavior changes qualitatively for L = 22 ,
24 data breaking the scaling t T h ∝ L e W/ Ω . Sim-ilar behavior heralds Anderson localization transition insingle particle disordered systems [45], hence, our datasuggest the presence of the transition to MBL phase in J - J model. Therefore, one has to reach a sufficientlylarge L to see the correct scaling of Thouless time, whichraises the question about the finite size effects at MBLtransition. Entanglement entropy and level statistics.
The en-tanglement entropy allows for insights in nature of MBLtransition [71–73]. The entanglement entropy of an eigen-state is defined as S E = − P i α i log( α i ), where α i areSchmidt basis coefficients (see e.g. [74]) associated withthe bipartition of the lattice into subsystems containingsites [ x, x + L/
2) and [ x + L/ , x + L ) (the sites are num-bered modulo L ). We average S E over the position of thecut x , over N ev (cid:54) min {N / , } eigenstates in themiddle of the spectrum ( σ = 0) of J - J model for systemsizes 12 (cid:54) L (cid:54)
24 (for L = 8 ,
10 we take N ev = 5) as wellas over more than 5000, 200, 50 disorder realizations re-spectively for L (cid:54) L = 22, L = 24. Finally, we obtainthe scaled entanglement entropy s E = S E /S RMT where S RMT ( L ) = ( L/
2) ln(2) + (1 / / / − / P i S zi = 0sector [75]. The resulting s E is shown in Fig. 3a). Foravailable system sizes, the scaled entanglement entropy s E : i) monotonically increases with L for W (cid:46) .
4; ii)monotonically decreases for W (cid:38)
11; iii) decreases forsmaller L and starts increasing for larger system sizes (asimilar reentrant behavior was observed e.g. in [47, 76]).The behavior i) clearly leads to an ergodic system at large L . In contrast, for large disorder strengths e.g. W = 15,an area law of entanglement entropy [77, 78] s E ∝ /L arises due to the emergent integrability of MBL phase[79–85]. Averaging r i = min { g i , g i +1 } / max { g i , g i +1 } (where g i = E i +1 − E i ) over eigenvalues correspondingto eigenstates from which s E was calculated, we obtain amean gap ratio r shown in Fig. 3b). The mean gap ratio r probes level statistics of the system, admitting valuescharacteristic for GOE and Poisson statistics for ergodicand localized systems [7, 86]. Similarly as for s E , themean gap ratio r follows the three types of behavior withsystem size depending on disorder strength W .To understand whether and at which disorder strengththe MBL transition takes place one has to study the in-terplay between the ii) and iii) trends. To this end, wefind the disorder strength W ∗ E ( L ) such that s E ( L −
1) = s E ( L + 1) for odd L and s E ( L −
2) = s E ( L + 2) for even L , for details see [68]. Smooth changes of s E with L and W assure that W ∗ E ( L ) is the largest disorder strength, fora given system size L , at which the volume-law S E ∝ L ,expected for an ergodic system, is still obeyed. Conse-quently, the disorder strength W ∗ E ( L ) is a lower boundfor the critical disorder strength W C of the transition to Figure 3. Finite size effects at MBL transition. a) The entanglement entropy s E of eigenstates of J - J model vs system size L for disorder strengths W = 1 . , ...,
15 (denoted on the color bar), dashed lines correspond to ergodic and MBL behavior; b)the same for the gap ratio r , dashed lines correspond to GOE and Poisson limits; c) W ∗ E,r and W T as function 1 /L for J - J model (see text); d) the same for XXZ model.
MBL phase. Fig. 3c) shows the relation between W ∗ E and 1 /L along with disorder strength W ∗ r obtained inanalogous manner for the average gap ratio r . Anotheraspect of finite size effects at MBL transition is revealedwhen, for given L , one finds a disorder strength W T ( L )for which the scaled entanglement entropy is close to theergodic limit, e.g. s E ( W T ) = 0 .
8. Such a criterion yields W T ∝ L . Equivalently, W T can be found as a disor-der strength for which the average gap ratio r departsfrom the GOE limit [44]. This allows us to identify thefollowing regimes: A) thermal, for W < W T , with en-tanglement entropy fulfilling the volume-law and closeto the value for chaotic spin chain S E ≈ S RMT ( L ) andlevel statistics well described by GOE; B) critical, for W ∗ E,r < W < W T , with S E < S RMT ( L ) but scalingsuper linearly with L and value of r increasing with L towards the GOE limit; C) MBL, for W < W ∗ E,r , withboth scaled entanglement entropy s E ( L ) and average gapratio r decreasing with system size L . Fig. 3d) shows thatbehavior of the XXZ model is similar (data for s E and r can be found in [68]). The three regimes resemble thequalitative picture of MBL transition proposed in [73].The asymptotic features of disordered spin chains de-pend on how W ∗ E,r and W T behave in thermodynamiclimit. For available system sizes, 8 (cid:54) L (cid:54)
24, the lin-ear scaling of W T with L as well as the linear scaling of W ∗ E,r with inverse of system size 1 /L , denoted by solidlines in Fig. 3 c), d), are accurately obeyed. Extrapo-lating the scalings (dashed lines in the same Fig.), leadsto the crossing W T = W ∗ E,r at L ≈
50 showing theincompatibility of the two scalings. Thus, it seems con-ceivable that studying eigenstates at at system size L would yield conclusive results about the L → ∞ limit,c.f. [47]. However, it is also possible that either of thescalings breaks down at smaller L achievable in the nearfuture with POLFED.The unveiled linear dependence of W ∗ E,r on 1 /L is con-sistently approached by data for all system sizes. Ex- trapolating to L → ∞ limit, we get estimates of criticaldisorder strength W J − J C ≈ . W XXZC ≈ . , (5)respectively for J - J and XXZ models. Our estimatefor W XXZC is larger than the value W C ≈ . W C ≈ . W C ≈ . s E ( W ) and r ( W ) curves, it does not rely on any finitesize scaling procedure. Our estimate for W XXZC is con-sistent with the lower bound W C > . W C > Conclusions.
The POLFED algorithm, thanks tothe employed polynomial spectral transformation, hasbetter scaling of computation time with matrix size N than the state-of-the-art SIMED algorithm. Avoiding thefill-in phenomenon, POLFED has a significantly lowermemory consumption than SIMED, moreover, its perfor-mance decreases only linearly with increasing the numberof non-zero off-diagonal matrix entries. For those rea-sons POLFED opens new pathways in studies of highlyexcited states of many-body systems, in particular withlong-range interactions. Understanding the relation ofPOLFED to alternative eigensolvers [69, 93, 94] is an in-teresting task for a further research.POLFED allowed us to study MBL transition in J - J model of size L (cid:54)
24. Such a system size is sufficient todemonstrate the breakdown of the scaling t T h ∝ L e W/ Ω of Thouless time [44]. Studying the system size scaling ofentanglement entropy S E of eigenstates we estimate thecritical disorder strength of transition to MBL phase. Acknowledgments.
We thank Fabien Alet andDominique Delande for insightful discussions. Thecomputations have been performed within PL-Grid In-frastructure, its support is acknowledged. M.L. ac-knowledges the Spanish Ministry MINECO (NationalPlan 15 Grant: FISICATEAMO No. FIS2016-79508-P, SEVERO OCHOA No. SEV-2015-0522,FPI), European Social Fund, Fundació Cellex, Fun-dació Mir-Puig, Generalitat de Catalunya (AGAURGrant No. 2017 SGR 1341, CERCA/Program), ERCAdG NOQIA, EU FEDER, MINECO-EU QUAN-TERA MAQS (funded by The State Research Agency(AEI) PCI2019-111828-2 / 10.13039/501100011033) ,and the National Science Centre, Poland-Symfonia GrantNo. 2016/20/W/ST4/00314. The support of Na-tional Science Centre, Poland under Unisono GrantNo. 2017/25/Z/ST2/03029 (Quantera: QTFLAG)is also acknowledged (J.Z.). P.S. acknowledges Na-tional Science Centre, Poland: ETIUDA grant No.2018/28/T/ST2/00401 . ∗ [email protected] † [email protected] ‡ [email protected][1] J. M. Deutsch, Phys. Rev. A , 2046 (1991).[2] M. Srednicki, Phys. Rev. E , 888 (1994).[3] L. D’Alessio, Y. Kafri, A. Polkovnikov, andM. Rigol, Advances in Physics , 239 (2016),https://doi.org/10.1080/00018732.2016.1198134.[4] R. Nandkishore and D. A. Huse, Ann. Rev. Cond. Mat.Phys. , 15 (2015).[5] F. Alet and N. Laflorencie, Comptes Rendus Physique , 498 (2018).[6] D. A. Abanin, E. Altman, I. Bloch, and M. Serbyn, Rev.Mod. Phys. , 021001 (2019).[7] V. Oganesyan and D. A. Huse, Phys. Rev. B , 155111(2007).[8] A. Pal and D. A. Huse, Phys. Rev. B , 174411 (2010).[9] J. A. Kjäll, J. H. Bardarson, and F. Pollmann, Phys.Rev. Lett. , 107204 (2014).[10] Y. Bar Lev, D. R. Reichman, and Y. Sagi, Phys. Rev. B , 201116 (2016).[11] R. Mondaini and M. Rigol, Phys. Rev. A , 041601(2015).[12] P. Prelovšek, O. S. Barišić, and M. Žnidarič, Phys. Rev.B , 241104 (2016).[13] P. Sierant, D. Delande, and J. Zakrzewski, Phys. Rev.A , 021601 (2017).[14] M. Kozarzewski, P. Prelovšek, and M. Mierzejewski,Phys. Rev. Lett. , 246602 (2018).[15] P. Sierant and J. Zakrzewski, New Journal of Physics ,043032 (2018).[16] N. Macé, N. Laflorencie, and F. Alet, SciPost Phys. ,50 (2019).[17] M. Schulz, C. A. Hooley, R. Moessner, and F. Pollmann,Phys. Rev. Lett. , 040606 (2019).[18] E. van Nieuwenburg, Y. Baum, and G. Refael, Pro-ceedings of the National Academy of Sciences , 9269(2019).[19] C. J. Turner, A. A. Michailidis, D. A. Abanin, M. Serbyn,and Z. Papić, Nature Physics , 745–749 (2018). [20] W. W. Ho, S. Choi, H. Pichler, and M. D. Lukin, Phys.Rev. Lett. , 040603 (2019).[21] V. Khemani, C. R. Laumann, and A. Chandran, PhysicalReview B , 161101 (2019).[22] T. Iadecola and M. Žnidarič, Phys. Rev. Lett. ,036403 (2019).[23] M. Schecter and T. Iadecola, Phys. Rev. Lett. ,147201 (2019).[24] A. J. A. James, R. M. Konik, and N. J. Robinson, Phys.Rev. Lett. , 130603 (2019).[25] T. Chanda, R. Yao, and J. Zakrzewski, “Coexistence oflocalized and extended phases: Many-body localizationin a harmonic trap,” (2020), arXiv:2004.00954.[26] P. Sala, T. Rakovszky, R. Verresen, M. Knap,and F. Pollmann, Physical Review X (2020),10.1103/physrevx.10.011047.[27] V. Khemani and R. Nandkishore, “Local constraints canglobally shatter hilbert space: a new route to quantuminformation protection,” (2019), arXiv:1904.04815.[28] T. Rakovszky, P. Sala, R. Verresen, M. Knap, andF. Pollmann, Physical Review B , 125126 (2020).[29] A. Smith, J. Knolle, R. Moessner, and D. L. Kovrizhin,Phys. Rev. Lett. , 176601 (2017).[30] M. Brenes, M. Dalmonte, M. Heyl, and A. Scardicchio,Phys. Rev. Lett. , 030601 (2018).[31] G. Magnifico, M. Dalmonte, P. Facchi, S. Pascazio, F. V.Pepe, and E. Ercolessi, “Real time dynamics and con-finement in the Z n Schwinger-Weyl lattice model for 1+1qed,” (2019), arXiv:1909.04821.[32] T. Chanda, J. Zakrzewski, M. Lewenstein, and L. Tagli-acozzo, Phys. Rev. Lett. , 180602 (2020).[33] G. Giudici, F. M. Surace, J. E. Ebot, A. Scardicchio, andM. Dalmonte, “Breakdown of ergodicity in disorderedu(1) lattice gauge theories,” (2019), arXiv:1912.09403.[34] F. M. Surace, P. P. Mazza, G. Giudici, A. Lerose,A. Gambassi, and M. Dalmonte, “Lattice gauge theoriesand string dynamics in Rydberg atom quantum simula-tors,” (2019), arXiv:1902.09551.[35] F. Pietracaprina, N. Macé, D. J. Luitz, and F. Alet,SciPost Phys. , 45 (2018).[36] T. Enss, F. Andraschko, and J. Sirker, Phys. Rev. B ,045121 (2017).[37] E. V. H. Doggen, F. Schindler, K. S. Tikhonov, A. D.Mirlin, T. Neupert, D. G. Polyakov, and I. V. Gornyi,Phys. Rev. B , 174202 (2018).[38] T. Chanda, P. Sierant, and J. Zakrzewski, Phys. Rev. B , 035148 (2020).[39] A. Goremykina, R. Vasseur, and M. Serbyn, Phys. Rev.Lett. , 040601 (2019).[40] A. Morningstar and D. A. Huse, Phys. Rev. B , 224205(2019).[41] P. T. Dumitrescu, A. Goremykina, S. A. Parameswaran,M. Serbyn, and R. Vasseur, Phys. Rev. B , 094205(2019).[42] N. Laflorencie, G. Lemarié, and N. Macé, “Chain break-ing and Kosterlitz-Thouless scaling at the many-body lo-calization transition,” (2020), arXiv:2004.02861.[43] J. Šuntajs, J. Bonča, T. Prosen, and L. Vidmar, “Ergod-icity breaking transition in finite disordered spin chains,”(2020), arXiv:2004.01719.[44] J. Šuntajs, J. Bonča, T. Prosen, and L. Vidmar, (2019),arXiv:1905.06345.[45] P. Sierant, D. Delande, and J. Zakrzewski, Phys. Rev.Lett. , 186601 (2020). [46] D. A. Abanin, J. H. Bardarson, G. D. Tomasi,S. Gopalakrishnan, V. Khemani, S. A. Parameswaran,F. Pollmann, A. C. Potter, M. Serbyn, and R. Vasseur,(2019), arXiv:1911.04501.[47] R. K. Panda, A. Scardicchio, M. Schulz, S. R. Taylor,and M. Žnidarič, EPL (Europhysics Letters) , 67003(2020).[48] K. Agarwal, S. Gopalakrishnan, M. Knap, M. Müller,and E. Demler, Phys. Rev. Lett. , 160401 (2015).[49] S. Bera, H. Schomerus, F. Heidrich-Meisner, and J. H.Bardarson, Phys. Rev. Lett. , 046603 (2015).[50] S. Bera, G. De Tomasi, F. Weiner, and F. Evers, Phys.Rev. Lett. , 196801 (2017).[51] L. Herviou, S. Bera, and J. H. Bardarson, Phys. Rev. B , 134205 (2019).[52] L. Colmenarez, P. A. McClarty, M. Haque, and D. J.Luitz, (2019), arXiv:1906.10701.[53] P. Sierant and J. Zakrzewski, Phys. Rev. B , 104201(2020).[54] C. Lanczos, Journal of Research of the National Bureauof Standards (1950).[55] G. H. Golub and C. F. Van Loan, Matrix computations ,Vol. 3 (JHU press, 2012).[56] P. R. Amestoy, I. S. Duff, J.-Y. L’Excellent, andJ. Koster, SIAM Journal on Matrix Analysis and Ap-plications , 15 (2001).[57] P. R. Amestoy, A. Guermouche, J.-Y. L’Excellent, andS. Pralet, Parallel Computing , 136 (2006), parallelMatrix Algorithms and Applications (PMAA’04).[58] J. Cullum and W. E. Donath, in (IEEE, 1974) pp. 505–509.[59] G. Golub and R. Underwood, in Mathematical Software ,edited by J. R. Rice (Academic Press, 1977) pp. 361 –377.[60] C. Bekas, E. Kokiopoulou, and Y. Saad, SIAM Journalon Matrix Analysis and Applications , 397 (2008).[61] H.-R. Fang and Y. Saad, SIAM Journal on ScientificComputing , A2220 (2012).[62] R. Li, Y. Xi, E. Vecharynski, C. Yang, and Y. Saad,SIAM Journal on Scientific Computing , A2512 (2016).[63] A. Pieper, M. Kreutzer, A. Alvermann, M. Galgon,H. Fehske, G. Hager, B. Lang, and G. Wellein, Jour-nal of Computational Physics , 226 (2016).[64] R. Silver and H. Röder, International Journal of ModernPhysics C , 735 (1994).[65] R. Silver, H. Roeder, A. Voter, and J. Kress, Journal ofComputational Physics , 115 (1996).[66] C. W. Clenshaw, Mathematics of Computation , 118(1955).[67] Y. Saad, SIAM Journal on Numerical Analysis , 687(1980).[68] See Supplemental Material at [URL will be inserted bypublisher].[69] R. Van Beeumen, G. D. Kahanamoku-Meyer, N. Y. Yao,and C. Yang, in Proceedings of the International Con-ference on High Performance Computing in Asia-PacificRegion (2020) pp. 179–187.[70] J. M. G. Gómez, R. A. Molina, A. Relaño, and J. Reta-mosa, Phys. Rev. E , 036209 (2002).[71] X. Yu, D. J. Luitz, and B. K. Clark, Phys. Rev. B ,184202 (2016).[72] V. Khemani, D. N. Sheng, and D. A. Huse, Phys. Rev.Lett. , 075702 (2017). [73] V. Khemani, S. P. Lim, D. N. Sheng, and D. A. Huse,Phys. Rev. X , 021013 (2017).[74] I. Bengtsson and K. Życzkowski, Geometry of Quan-tum States: An Introduction to Quantum Entanglement (Cambridge University Press, 2006).[75] L. Vidmar and M. Rigol, Phys. Rev. Lett. , 220603(2017).[76] M. Serbyn, Z. Papić, and D. A. Abanin, Phys. Rev. X , 041047 (2015).[77] B. Bauer and C. Nayak, Journal of Statistical Mechanics:Theory and Experiment , P09005 (2013).[78] M. Serbyn, Z. Papić, and D. A. Abanin, Phys. Rev. Lett. , 260601 (2013).[79] M. Serbyn, Z. Papić, and D. A. Abanin, Phys. Rev. Lett. , 127201 (2013).[80] D. A. Huse, R. Nandkishore, and V. Oganesyan, Phys.Rev. B , 174202 (2014).[81] V. Ros, M. Mueller, and A. Scardicchio, Nuclear PhysicsB , 420 (2015).[82] J. Z. Imbrie, Phys. Rev. Lett. , 027201 (2016).[83] T. B. Wahl, A. Pal, and S. H. Simon, Phys. Rev. X ,021018 (2017).[84] M. Mierzejewski, M. Kozarzewski, and P. Prelovšek,Phys. Rev. B , 064204 (2018).[85] S. J. Thomson and M. Schiró, Phys. Rev. B , 060201(2018).[86] Y. Y. Atas, E. Bogomolny, O. Giraud, and G. Roux,Phys. Rev. Lett. , 084101 (2013).[87] D. J. Luitz, N. Laflorencie, and F. Alet, Phys. Rev. B , 081103 (2015).[88] A. B. Harris, Journal of Physics C: Solid State Physics , 1671 (1974).[89] J. T. Chayes, L. Chayes, D. S. Fisher, and T. Spencer,Phys. Rev. Lett. , 2999 (1986).[90] A. Chandran, C. R. Laumann, and V. Oganesyan,(2015), arXiv:1509.04285.[91] N. Macé, F. Alet, and N. Laflorencie, Phys. Rev. Lett. , 180601 (2019).[92] T. Devakul and R. R. P. Singh, Phys. Rev. Lett. ,187201 (2015).[93] M. Bollhöfer and Y. Notay, Computer Physics Commu-nications , 951 (2007).[94] E. Polizzi, Phys. Rev. B , 115112 (2009).[95] N. Bell and M. Garland, Efficient sparse matrix-vectormultiplication on CUDA , Tech. Rep. (2008).[96] S. Acer, O. Selvitopi, and C. Aykanat, Parallel Comput-ing , 71 (2016).[97] S. Chen, J. Fang, D. Chen, C. Xu, and Z. Wang, arXivpreprint arXiv:1805.11938 (2018).[98] H. A. H. Baca and F. de Luz Palomino Valdivia, in (2019) pp. 1–4. SUPPLEMENTARYConvergence of POLFED
POLFED performs the iteration of the block Lanc-zos method for the transformed matrix P Kσ ( ˜ H ) untilall of the residual norms associated with eigenvalues θ i (cid:62) p = 0 .
17 vanish within numerical precision. Fig. 4shows the number n ev of converged (with the vanishingresidual norm) eigenpairs in few runs of POLFED. Re-
100 150 200 250 m n ev a)
300 350 400 450 m n ev b) m n ev c) m n ev d) Figure 4. Convergence of POLFED. Number of convergedeigenpairs n ev as a function of number of Lanczos steps m .Dotted lines correspond to 8 different disorder realizationswith disorder strength W = 5 . J - J model, system sizeis L = 20. The block size s = 1. Panels a), b), c) and d) cor-respond, respectively, to the number of requested eigenvalues N ev = 100 , , , gardless of the number N ev of requested eigenvalues, n ev increases rapidly, once the number m of Lanczos stepsexceeds a certain threshold value which is typically twicelarger than N ev . Therefore, the eigenpairs start to con-verge only when the Q m matrix contains the full basis ofsubspace of Hilbert space spanned by eigenvectors witheigenvalues close to the energy σ . It is beneficial to stopthe algorithm once n ev (cid:62) N ev . Further increase of m does not lead to an increase in n ev because of the largedensity of states of P Kσ ( ˜ H ) close to θ ≈ .
15 due to thethe secondary minima of P Kσ ( (cid:15) ) (see Fig. 1 of the maintext).POLFED follows the convergence pattern describedabove provided that the polynomial spectral transforma-tion P Kσ (in particular, its order K ) is chosen in such away that the number N ev of requested eigenvalues corre-sponds to the number n P of eigenvalues of P Kσ ( ˜ H ) that are accessible to the method (i.e. fulfill the θ > p = 0 . n P , POLFED uses the function P Kσ ( (cid:15) ) as well as the density of states of the ˜ H . As thedensity of states in the middle of the spectrum of thebenchmark models considered in this work we use an an-alytical expression ˜ ρ (0) = ( E − E ) N / Γ, whereΓ = q L [(1 + J ) / (1 + J ) /
16 + W / , (6)as obtained in [44]. The fluctuations of density of statesbetween disorder realizations lead to the fluctuations ofthe threshold value of m beyond which the convergenceoccurs – see Fig. 4. Those fluctuations are enhancedwhen disorder strength increases. However, our tests in-dicate that the convergence occurs for each of the consid-ered disorder values ( W (cid:54)
15) and disorder realizationsfor m < . N ev .Testing a variety of the polynomial spectral transfor-mation P Kσ ( ˜ H ) as well as different stopping criteria, wechecked that POLFED allows to minimize the time of cal-culation until the convergence is reached and, at the sametime, allows to keep a relatively large the total numberof eigenpairs obtained in a single run.When the block size s of the Lanczos method is in-creased, the total number of vectors generated in the it-eration, ms , required for the convergence of algorithm,is also increased. However, for the typical productionruns done in this work, i.e. with the block size s (cid:54) N ev (cid:62) ms < . N ev . Thus, duringits start, POLFED allocates 2 . N ev columns of the ma-trix Q m . The associated memory consumption is propor-tional to N ev N , which has the dominant contribution tototal memory occupation of POLFED. Technical details of POLFED
Calculation of the product of P Kσ ( ˜ H ) with subse-quent columns of Q j is the most time consumingstep of POLFED. The recurrence relation T n +2 ( x ) =2 xT n +1 ( x ) − T n ( x ) fulfilled by Chebyshev polynomialsreduces this product to multiplication of vectors by thesparse matrix ˜ H and basic linear algebra operations. TheClenshaw algorithm [66], allows us to reduce the numberoperations needed to calculate the product.The efficiency of computation of P Kσ ( ˜ H ) x where x ∈ R N is crucially dependent on efficiency of the singlesparse matrix vector multiplication ˜ Hx . In the currentversion of POLFED we store the ˜ H matrix in CSR for-mat. We do not store the off-diagonal Hamiltonian en-tries of H as they are all equal to H ij = 1 /
2. Onone hand this reduces the memory consumption asso-ciated with storing of the Hamiltonian matrix. On theother hand, POLFED does not access the values of H ij during the matrix-vector multiplication which increasesthe efficiency of the code. Throughout this work, weconsider block sizes s (cid:54)
24 calculating the products of P Kσ ( ˜ H ) with columns of matrix Q j independently. Eachof the products is calculated on a single core. Effec-tively, POLFED performs the computation in parallelon s cores. The re-orthogonalization of columns of ma-trix U j obtained in Lanczos step against the columns ofmatrices { Q i } ji =1 is parallelized similarly: each of the s cores orthogonalizes a single column of U j against thecolumns of { Q i } ji =1 .The matrix-vector multiplications ˜ Hx could be per-formed on multiple cores with use of external sparse ba-sic linear algebra libraries, resulting in higher degree ofparallelism in POLFED. Moreover, the promising way ofenhancing the performance of POLFED is to optimizethe sparse matrix-vector product, a subject that recentlyreceived attention both on CPUs as well as on GPUs[95–98]. Benchmark for disorder spin chains
In this section we compare performance of POLFEDwith state-of-the-art SIMED code for
XXZ and J - J models. Benchmark results are shown in Tab. I and inTab. II.The linear scaling of density of states ρ (0) with N im-plies that the order of the polynomial spectral transfor-mation K ∝ N . Therefore, up to a factor polynomialin L , the computation time of POLFED scales as N .The total CPU time t CP U for POLFED increases by afactor of ≈
16 both for
XXZ model (Tab. I) and for J - J model (Tab. II) when the system size L increasesby 2. Typically, the increase of t CP U is slightly largerwhen the number N cores of cores (equal to the blocksize s for POLFED) increases. The memory consump-tion of POLFED indeed scales as N ev N up to constanta additive factor due to the storing of the Hamiltonianmatrix as Tab. I and Tab. II show. The low memory con-sumption of POLFED allows for calculations on a singlenode ( N cores (cid:54)
24 on the supercomputer Prometheus,ACK Cyfronet AGH, Krakow) for both models as longas L (cid:54) LU decomposition of the Hamilto-nian, it scales as c ( L ) N (up to terms polynomial in thesystem size L ). The factor c ( L ) describes fill-in of thematrix. Tests performed in [35] indicate a phenomeno-logical scaling c ( L ) ∝ L/ for XXZ spin chain. Thisresults in total memory needed to store the LU factorsto be ≈ GB and 14000 GB respectively for L = 24and L = 26. The actual memory usage, due to peaksof allocated/de-allocated memory in the SIMED is sig-nificantly higher as shown in Tab. I. The rapid increaseof memory consumption with L forces one to use a verylarge number of nodes in calculations with SIMED, even- L t
CPU [ h ] N cores t W [ h ] RAM [ GB ] N ev P O L F E D
20 3.1 1 3.1 3.9 100022 62.2 4 15.5 21.2 140024 1503 24 62.6 114 200026 19870 24 828 488 2000 S I M E D
20 0.5 20 0.026 22 10022 20.2 120 0.17 244 10024 840 2880 0.23 12288 5026 36000 48000 0.75 204800 50Table I. POLFED vs SIMED for
XXZ spin chain: t CPU is total CPU time, N cores is the number of cores used in cal-culation, t W is total execution time, RAM is total memoryoccupation, N ev is the number of obtained eigenpairs in themiddle of the spectrum ( σ = 0). Tested on Intel Xeon E5-2680v3 (2.5GHz); SIMED data for L = 20 ,
22 obtained onIntel Ivybridge E5-2680 (2.8GHz), extracted from [35].
L t
CPU [ h ] N cores t W [ h ] RAM [ GB ] N ev P O L F E D
20 3.1 1 3.1 0.8 10020 3.6 1 3.6 3.9 100022 60.2 1 60.2 3.4 10022 63.2 2 31.6 4.5 20022 105 8 13.1 21.3 140024 3400 24 142 115 2000 S I M E D
20 2.4 36 0.067 70 10022 100 468 0.22 1840 10022 120 468 0.26 1840 200Table II. POLFED vs SIMED for J - J spin chain: t CPU is total CPU time, N cores is the number of cores used in cal-culation, t W is total execution time, RAM is total memoryoccupation, N ev is the number of obtained eigenpairs in themiddle of the spectrum ( σ = 0). POLFED tested on In-tel Xeon E5-2680v3 (2.5GHz); SIMED tested on Intel XeonGold 6140 CPU (2.3GHz), data provided by courtesy of F.Alet. tually making the calculations infeasible, even on largesupercomputers. Theoretically, the calculation time ofSIMED, dominated by the LU factorization, should beproportional to the number of elements in the factorsyielding the scaling of total CPU time t CP U ∝ c ( L ) N .However, as Tab. I shows, t CP U scales more rapidly withsystem size, increasing approximately 40 times when L increases by 2. Altogether, the system size scaling of t CP U is better for POLFED.Another aspect of the fill-in phenomenon of SIMEDis that it is quite unpredictable. For instance, the co-efficient c ( L ) may change after reordering of the basis.It is, however, clear that the fill-in becomes much moresevere as the number N nz of non-zero off-diagonal ele-ments increases. Tab. II shows that the total memoryconsumption of SIMED for J - J model is increased, incomparison to resources needed for XXZ model, by afactor of ≈ . ≈ . L = 20 and L = 22. The total CPU times t CP U for POLFED andSIMED for J - J model are very similar for L = 20 , L = 24 infeasible on presentday supercomputers. At the same time, POLFED allowsto obtain results for J - J model of size L = 24 with re-sources similar to the XXZ model – such a calculationfits in a single node of a supercomputer.Another advantage of POLFED is that it allows for asubstantial increase of the number of requested eigen-values, N ev , without a significant increase in the to-tal calculation time. This can be readily understood.When replacing N ev → αN ev where α >
1, the con-dition that N ev eigenvalues of P Kσ ( ˜ H ) are larger than p = 0 .
17 results in the order of the spectral transfor-mation K → K/α . Even though the total number ofLanczos iterations needed for the convergence of the al-gorithm increases by a factor of α , the cost of calcula-tion of a single polynomial spectral transformation de-creases α times. The re-orthogonalization performed byPOLFED is the only source of increase of total CPUtime when N ev → αN ev . This can be seen in Tab. IIas tests for L = 20 ,
22 were performed for few values of N ev . The change of total CPU time with number of re-quested eigenvalues is more significant for SIMED. Forinstance, Fig. 6. of [35] shows that increase of N ev from100 to 1000 results in approximately 6 times larger N ev for XXZ model of size L = 20.Typically, MBL calculations require averaging over dis-order realizations. The POLFED allows to find eigen-pairs of the disordered spin chains on relatively smallnumber of cores so that the averaging over disorder re-alizations can be done by performing calculations formany disorder realizations independently at the sametime. In contrast, SIMED requires much larger amountof resources. Ultimately, due to smaller total executiontimes, for XXZ model of sizes L (cid:54)
24 SIMED allows toget results for a comparable, but slightly larger number ofdisorder realizations using a fixed amount of CPU time(assuming than one is able to perform calculations for L = 24 simultaneously on 120 nodes of a cluster). For J - J model at system size L = 20 ,
22 POLFED has an ad-vantage. Moreover, while SIMED calculations for J - J model at L = 24 are infeasible, they can be readily doneby POLFED. The situation is similar for XXZ model at L = 26: single SIMED run requires 2000 nodes of a su-percomputer, and the calculation is performed in singleprecision [35] whereas POLFED calculation, in doubleprecision, requires only few nodes of a supercomputer. Extraction of Thouless time
To extract Thouless time t T h from spectral form factor K ( τ ) for system size L (cid:54)
18 we use data from full exactdiagonalization and follow the procedure outlined in [44].To this end we calculate the spectral form factor (SFF) according to its definition K ( τ ) = 1 Z *(cid:12)(cid:12)(cid:12)(cid:12)(cid:12)(cid:12) N X j =1 g ( E j )e − iE j τ (cid:12)(cid:12)(cid:12)(cid:12)(cid:12)(cid:12) + . (7)Subsequently, we perform the unfolding, during whichthe level staircase function σ ( (cid:15) ) = P i Θ( (cid:15) − ε i ) (obtainedfrom the set of eigenvalues of the system { ε i } ordered inascending manner) is separated into smooth and fluctu-ating parts σ ( E ) = σ ( E ) + δσ ( E ) and the eigenvaluesare mapped via ε j → E j = σ ( ε j ). As the smooth part σ ( E ) we take a polynomial of degree n p = 10 fitted to thelevel staircase function σ ( E ). To calculate SFF we use g ( (cid:15) ) ∝ exp( − ( (cid:15) − ¯ (cid:15) ) / ησ (cid:15) ), where ¯ (cid:15) denotes the averageof the unfolded eigenvalues for given disorder realization (cid:15) i , σ (cid:15) is the standard deviation of { (cid:15) i } and η = 0 .
3. Thischoice of parameters follows precisely [53]. Then, we cal-culate ∆ K ( t/t H ) = (cid:12)(cid:12)(cid:12)(cid:12) log (cid:18) K ( t/t H ) K GOE ( τ = t/t H ) (cid:19)(cid:12)(cid:12)(cid:12)(cid:12) , (8)where the spectral form factor for GOE is given by K GOE ( τ ) = ( τ − τ ln(1 + 2 τ ) for τ (cid:54) , − τ ln( τ +12 τ − ) for τ > . (9)The Thouless time t T h is the smallest positive time forwhich ∆ K ( t/t H ) < a . We choose the value of cut-off a = 0 . − − − τ − − K ( τ ) W=4.6W=5W=5.4W=5.8W=6.4W=7W=7.5
Figure 5. Spectral form factor K ( τ ) of J - J model, systemsize L = 24. Black dashed lines shows spectral form factorof GOE, K GOE ( τ ). The dashed lines show K F ( τ ) (9) fittedto K ( τ ); dots denote the obtained rescaled Thouless time: τ Th = t Th /t H , where t H is the Heisenberg time. For system sizes L = 20 , ,
24 we obtain N ev = 2500consecutive eigenvalues from the middle of spectrum.Firstly, we perform the unfolding procedure using fittingthe level staircase function with a polynomial of degree0 Figure 6. Determination of crossings s E ( L −
1) = s E ( L +1) for J - J model. The rescaled entropy s E , denoted bydashed line, is plotted as a function of disorder strength W for system sizes L = 10 , L = 14 , L = 18 , L =22 ,
24 respectively on panels a), b), c), d). s E ( W ) curves arefitted with polynomials of third degree (shown by solid lines)in vicinity of the crossing point. The crossing points of thepolynomials determine W ∗ E ( L ) for L = 11 , , , n p = 3. Then, we calculate the spectral form factor ac-cording to the definition (7) considering only the calcu-lated eigenvalues in the sum. Since the number of disor-der realizations we have for the largest system considered( L = 24) is only 50, we fit the spectral form factor withthe following formula K F ( τ ) = K GOE ( τ ) + c exp ( − c τ c ) , (10)where c , c , c are fit parameters. The results are shownin Fig. 5. The formula (10) provides very good fits of tothe spectral form factor for smaller system sizes ( L (cid:54) τ (cid:46) . τ T h . Thus, to extract τ T h for L = 22 ,
24 weuse K F ( τ ) in (8).Fig. 5 illustrates also an another aspect of calculationof the Thouless time when not all of the eigenvalues ofthe system are available. The number of N ev eigenvaluesdetermines the value of tau N ev below which the spectralform factor rapidly increases. In our case, as can be seenin Fig. 5, τ N ev ≈ · − . Once the extracted value of τ T h is significantly bigger than τ N ev , the value of τ T h isnot affected by the fact that N ev (cid:28) N . Extraction of W ∗ To extract the values of W ∗ E ( L ) we plot (for odd L ), s E ( L −
1) and s E ( L + 1) as functions of disorder strength Figure 7. The rescaled entanglement entropy s E of eigenstatesof XXZ model as a function of system size L for disorderstrengths W = 0 . , ..., s E = 1. W , and perform a fit with third order polynomial in thevicinity of the crossing point as shown in Fig. 6. Thecrossing point of the two polynomials is then the value of W ∗ E ( L ). A similar, analysis is performed for even L anddata for s E ( L − , s E ( L + 2).The values of W ∗ r ( L ) are found in analogous mannerfrom r ( L −
1) and r ( L + 1) ( r ( L −
2) and r ( L + 2)) asfunctions of disorder strength W for odd (even) L .We perform a similar analysis for XXZ model obtain-ing W ∗ E,r ( L ) for that system. The scaled entanglemententropy and average gap ratio for XXZ model are shownin Fig. 7 and Fig. 8. The scaled entanglement entropy s E is obtained when we average S E over the position ofthe cut x , over N ev (cid:54) min {N / , } eigenstates inthe middle of the spectrum ( σ = 0) of J − J modelfor system sizes 12 (cid:54) L (cid:54)
22 (whereas for L = 8 ,