Many-Body Localization: Transitions in Spin Models
John Schliemann, Joao Vitor I. Costa, Paul Wenk, J. Carlos Egues
MMany-Body Localization: Transitions in Spin Models
John Schliemann , Jo˜ao Vitor I. Costa , Paul Wenk , and J. Carlos Egues Institute for Theoretical Physics, University of Regensburg, Regensburg, Germany Instituto de Fısica de Sao Carlos, Universidade de Sao Paulo, Sao Carlos, Brazil (Dated: November 24, 2020)We study the transitions between ergodic and many-body localized phases in spin systems, subjectto quenched disorder, including the Heisenberg chain and the central spin model. In both casessystems with common spin lengths 1 / s r ) of the averaged consecutive-gap ratio (cid:104) r (cid:105) for different disorder realizations. For both types of systemsand spin lengths we find a maximum in ∆ s r as a function of disorder strength, accompanied byan inflection point of (cid:104) r (cid:105) , signaling the transition from ergodicity to many-body localization. Thecritical disorder strength is found to be somewhat smaller than the values reported in the recentliterature. Further information about the transitions can be gained from the probability distributionof expectation values within a given disorder realization. I. INTRODUCTION
Many-body localization has become in the last yearsone of the most intensively growing areas of research incondensed matter physics and beyond . It denotes theabsence of thermalization in an isolated interacting quan-tum system in the presence of typically strong disorder.In the opposite ergodic phase the eigenstate thermaliza-tion hypothesis is fulfilled stating that any appropriatesubsystem of the isolated total system (being in a purestate) is accurately described by equilibrium statisticalmechanics . On the other hand, the presence of inter-actions distinguishes many-body localization from tradi-tional Anderson localization .In this work we revisit the transition from the ergodicto the many-body localized phase in disordered Heisen-berg spin chains and compare it with the behavior ofcentral spin models also subject to quenched disorder.For both types of systems we consider the spin lengths1 / (cid:104) r (cid:105) , suggesting aclose analogy to classic phase transitions with (cid:104) r (cid:105) beingan order parameter. These central observations are sum-marized in Figs. 2, 7 for the Heisenberg chain and thecentral spin model, respectively.An additional tool in the analysis of the transition fromthe ergodic to the many-body localized phase is the prob-ability distribution of expectation values within a givendisorder realization. The corresponding data is containedin Figs 5 and 9.This paper is further organized as follows: In section IIwe introduce the spin models to be studied and summa-rize the underlying theoretical techniques. Our numericalresults are presented in section III, and we close with a summary and an outlook in section IV. II. MODEL AND APPROACHA. Spin Hamiltonian
We study a periodic Heisenberg spin chain interactingwith an additional central spin and being subject to anuniaxial quenched disorder field on each site of the chain, H = J K (cid:88) i =1 (cid:126)I i · (cid:126)I i +1 + AK (cid:126)S · K (cid:88) i =1 (cid:126)I i + 2 S K (cid:88) i =1 h i I zi , (1)where the parameter J describes the coupling of the K chain (or bath) spins (cid:126)I i = (cid:126)I i + K to their nearest neigh-bors. The coupling to the central spin (cid:126)S is parametrizedby A , and the random magnetic field h i is chosen froma uniform distribution within the interval [ − h, h ]. Inwhat follows all spins will have length S = I = 1 / S = I = 1.For S = I = 1 /
2, the Hamiltonian (1) with A = 0,i.e. the Heisenberg spin-1 / . The factor 2 S in front of the dis-order term makes contact to the usual parametrizationfor S = 1 / S (cid:55)→ qS , I (cid:55)→ qI leads to H (cid:55)→ q H . Thus,disorder-induced effects on the dynamics should occur atthe same disorder strength h . As we shall see in sec-tion III, this remains approximately true when switchingbetween S = I = 1 / S = I = 1.Recent work by Hetterich et al. studied the full model(1) for S = I = 1 / J = 1. Asthese authors argue, dividing the interaction parameter A by the number of central spins K , as done in Eq. (1), en-sures that the spectral bandwidth of that coupling term a r X i v : . [ c ond - m a t . d i s - nn ] N ov is approximately independent of the system size. We willsee in the following that this stipulation has indeed ad-vantages when comparing data for different numbers ofspins. B. Random Matrix Theory
An important method to distinguish a many-body lo-calized phase from an ergodic phase is random matrixtheory, which is a theory for statistical fluctuations ofenergy levels of a given quantum system . A moderntool to analyze the energy level statistics is the consecu-tive gap ratio defined as r n = min { s n , s n − } max { s n , s n − } = min (cid:26) ¯ r n , r n (cid:27) , ¯ r n = s n s n − , (2)where s n = e n +1 − e n is the difference between two neigh-boring energy levels e n +1 , e n . In the strictly many-bodylocalized (or integrable) phase, characterized by an ex-tensive number of independent conserved quantities, thedifferences s n obey Poisson statistics, and the probabilitydistribution for the random variable r = r n can easily bedetermined to be p ( r ) = 2(1 + r ) . (3)On the other hand, in the fully ergodic phase, a sys-tem of the type (1) is generally assumed to be describedby the Gaussian orthogonal ensemble (GOE) of randommatrices . Here the analysis of small random matricesmimicking the Hamiltonian predicts the pertaining prob-ability distribution to be p ( r ) = 274 r + r (1 + r + r ) / , (4)which can be seen as an analog of the classic Wignersurmises for probability distributions governing the tra-ditional random variable s = s n .As a consequence, the lowest moments of the probabil-ity distribution (3) in the integrable case are given by (cid:104) r (cid:105) p = 2 ln 2 − ≈ . , (5) (cid:104) r (cid:105) p = 3 − ≈ . , (6)∆ p r = (cid:113) (cid:104) r (cid:105) p − (cid:104) r (cid:105) p ≈ . , (7)with (cid:104)·(cid:105) p = (cid:90) dr ( p ( r ) · ) , (8)whereas in the ergodic situation (4) we have (cid:104) r (cid:105) p = 4 − √ ≈ . , (9) (cid:104) r (cid:105) p = 274 ln (cid:18) √ (cid:19) − − √ ≈ . , (10)∆ p r ≈ . , (11) C. Statistical Data Analysis
1. Disorder ensemble
Consider an ensemble of Q realizations of the local dis-order field h i , i ∈ { , . . . , K } . Each disorder realization,or sample, labeled by α ∈ { . . . , Q } leads to a prob-ability distribution p α ( r ) for the consecutive gap ratio r ∈ [0 , f ( r ), these dis-tributions determine the realization-dependent averages(or expectation values) (cid:104) f (cid:105) α = (cid:90) drp α ( r ) f ( r ) (12)with variances (∆ α f ) = (cid:104) f (cid:105) α − (cid:104) f (cid:105) α . (13)The disorder-averaged probability distribution p ( r ) atgiven disorder strength (and other system parameters)is p ( r ) = lim Q →∞ Q Q (cid:88) α =1 p α ( r ) , (14)and a good numerical estimate for this quantity is p ( r ) ≈ Q Q (cid:88) α =1 p α ( r ) , Q (cid:29) , (15)for sufficiently large Q . The disorder-averaged expecta-tion values of f ( r ) reads (cid:104) f (cid:105) p = (cid:90) drp ( r ) f ( r ) = lim Q →∞ Q Q (cid:88) α =1 (cid:104) f (cid:105) α . (16)
2. Probability distribution for average within sample
On the other hand, we can view the numbers x = (cid:104) f (cid:105) α as random variables according to the distribution s ( x ) = 1(2 h ) K (cid:90) h − h dh · · · (cid:90) h − h dh K · δ (cid:18) x − (cid:90) drp ( r ; h , . . . , h K ) f ( r ) (cid:19) (17)where p ( r ; h , . . . , h K ) = p α is the probability distribu-tion within a system with local disorder fields h , . . . , h K forming the disorder realization α . Thus, the disorder-averaged expectation value of f ( r ) can be formulated as (cid:104) f (cid:105) s = (cid:90) dxs ( x ) x = lim Q →∞ Q Q (cid:88) α =1 (cid:104) f (cid:105) α = (cid:104) f (cid:105) p , (18)where the integration goes over all values of x = f ( r )for r ∈ [0 , s ( x ) does in general notcoincide with p ( r ) even for x = f ( r ) = r . Moreover s ( x )will of course have a dependence on the function f , whichwe, however, shall suppress in the notation.The finite average¯ f = 1 Q Q (cid:88) α =1 (cid:104) f (cid:105) α ≈ (cid:104) f (cid:105) s , (19)is, again for appropriately large Q , an approximation tothe expression (18). On the other hand, it is a sum ofstochastically independent (and therefore uncorrelated)random variables with identical probability distribution,so that the resulting joint distribution is π ( x , . . . , x Q ) = Q (cid:89) α =1 s ( x α ) . (20)Thus, the expectation value of ¯ f is (cid:104) ¯ f (cid:105) π = 1 Q Q (cid:88) α =1 (cid:104) f (cid:105) s = (cid:104) f (cid:105) s . (21)
3. Sample-to-sample variance
For the variance pertaining to the expectation value(21) one finds (cid:68)(cid:0) ¯ f − (cid:104) f (cid:105) s (cid:1) (cid:69) π = Q (cid:88) α,β =1 (cid:104) ( (cid:104) f (cid:105) α − (cid:104) f (cid:105) s ) ( (cid:104) f (cid:105) β − (cid:104) f (cid:105) s ) (cid:105) π Q = 1 Q Q (cid:88) α =1 (cid:68) ( (cid:104) f (cid:105) α − (cid:104) f (cid:105) s ) (cid:69) s = (∆ s f ) Q (22)with (∆ s f ) = (cid:104) ( f − (cid:104) f (cid:105) s ) (cid:105) s = lim Q →∞ Q Q (cid:88) α =1 ( (cid:104) f (cid:105) α − (cid:104) f (cid:105) s ) . (23)Therefore, the fluctuations of the finite average (18)around its expectation value (21) are characterized bythe standard deviation∆ π ¯ f = (cid:114)(cid:68)(cid:0) ¯ f − (cid:104) f (cid:105) s (cid:1) (cid:69) π = ∆ s f √ Q (24)and shows the familiar decay ∝ / √ Q . This result fol-lows of course also from the Lindeberg-Levy central limittheorem applied to the sum (19) of random variables. Anapproximate expression for the variance (23) is(∆ s f ) ≈ Q Q (cid:88) α =1 ( (cid:104) f (cid:105) α − ¯ f ) , Q (cid:29) . (25) Note that the variance (23) also occurs as a contribu-tion to the variance of f ( r ) calculated from the disorder-averaged probability distribution (14),(∆ p f ) = (cid:90) drp ( r )( f ( r ) − (cid:104) f (cid:105) p ) (26)= lim Q →∞ Q Q (cid:88) α =1 (cid:90) drp α ( r )( f ( r ) − (cid:104) f (cid:105) p ) = lim Q →∞ Q Q (cid:88) α =1 ( (cid:104) f (cid:105) α − (cid:104) f (cid:105) α )+ lim Q →∞ Q Q (cid:88) α =1 ( (cid:104) f (cid:105) α − (cid:104) f (cid:105) p )= lim Q →∞ Q Q (cid:88) α =1 (∆ α f ) + (∆ s f ) . (27)Hence, the variance (26) with respect to the averagedprobability distribution (14) is the average of all vari-ances within the disorder realizations around their in-dividual expectation value of f , plus the variance (23)describing the fluctuations of these expectation valuesaround their mean. An approximate expression for theabove result is(∆ p f ) ≈ Q Q (cid:88) α =1 (cid:90) drp α ( r )( f ( r ) − ¯ f ) (28)= 1 Q Q (cid:88) α =1 (∆ α f ) + 1 Q Q (cid:88) α =1 ( (cid:104) f (cid:105) α − ¯ f ) (29)where again Q (cid:29)
4. Variance of the variance
Let us now analyze the statistical fluctuations of ther.h.s. of the approximate quantity (25) for finite Q anddefine¯ g := 1 Q Q (cid:88) α =1 ( (cid:104) f (cid:105) α − ¯ f ) = 1 Q Q (cid:88) α =1 (cid:0) (cid:104) f (cid:105) α − ¯ f (cid:1) . (30)An important difference between the above expressionand the quantity (19) is that above the summands de-pend, again via Eq. (19), on all expectation values (cid:104) f (cid:105) α , α ∈ { , . . . , Q } , and are therefore distributed accordingto the joint probability distribution (20). The expecta-tion value of the random variable (30) with respect tothe latter distribution is (cid:104) ¯ g (cid:105) π = 1 Q Q (cid:88) α =1 (cid:10) (cid:104) f (cid:105) α (cid:11) π − Q Q (cid:88) α,β =1 (cid:104)(cid:104) f (cid:105) α (cid:104) f (cid:105) β (cid:105) π = Q − Q (cid:0) (cid:104) f (cid:105) s − (cid:104) f (cid:105) s (cid:1) = Q − Q (∆ s f ) (31)which approaches the variance (23) for large Q . Notealso that, in contrast to Eq. (18), it holds (cid:104) f (cid:105) s = (cid:90) dxs ( x ) x = lim Q →∞ Q Q (cid:88) α =1 (cid:104) f (cid:105) α (cid:54) = (cid:104) f (cid:105) p . (32)It is now straightforward to establish that (cid:10) ( (cid:104) f (cid:105) α − ¯ f − (cid:104) ¯ g (cid:105) π )( (cid:104) f (cid:105) β − ¯ f − (cid:104) ¯ g (cid:105) π ) (cid:11) π = δ αβ (cid:0) (cid:104) f (cid:105) s − (cid:104) f (cid:105) s (cid:1) + O (cid:18) Q (cid:19) , (33)and therefore, analogously as in Eq. (22),(∆ π ¯ g ) = Q (cid:88) α,β =1 (cid:68) ( (cid:104) f (cid:105) α − ¯ f − (cid:104) ¯ g (cid:105) π )( (cid:104) f (cid:105) β − ¯ f − (cid:104) ¯ g (cid:105) π ) (cid:69) π Q = (∆ s f ) Q + O (cid:18) Q (cid:19) (34)where(∆ s f ) = (cid:104) f (cid:105) s − (cid:104) f (cid:105) s = (cid:68)(cid:0) f − (cid:104) f (cid:105) s (cid:1) (cid:69) s = (cid:0) ∆ s (cid:0) (∆ s f ) (cid:1)(cid:1) . (35)
5. Consecutive gap ratio
In what follows we will be mainly concerned with thefunction f ( r ) = r where we have(∆ s r ) = lim Q →∞ Q Q (cid:88) α =1 ( (cid:104) r (cid:105) α − (cid:104) r (cid:105) s ) = lim Q →∞ Q Q (cid:88) α =1 (cid:0) (cid:104) r (cid:105) α − (cid:104) r (cid:105) s (cid:1) (36)with (cid:104) r (cid:105) s = (cid:90) dxs ( x ) x = lim Q →∞ Q Q (cid:88) α =1 (cid:104) r (cid:105) α = (cid:104) r (cid:105) p , (37)where s ( x ) is, in accordance with Eq. (17), the prob-ability distribution for the random variable x = (cid:104) r (cid:105) α .Finally, as seen in Eq. (35), the variance of the variance(36) is determined by(∆ s r ) = (cid:104) r (cid:105) s − (cid:104) r (cid:105) s = (cid:68)(cid:0) r − (cid:104) r (cid:105) s (cid:1) (cid:69) s = (cid:0) ∆ s (cid:0) (∆ s r ) (cid:1)(cid:1) . (38) III. NUMERICAL RESULTS
The Hamiltonian (1) obviously conserves the z -component of the total spin, (cid:126)J = (cid:126)S + K (cid:88) i =1 (cid:126)I , [ H, J z ] = 0 . (39) Thus, in order to apply random matrix theory, the spec-tra of each invariant subspace of J z have to be analyzedseparately . In this section we present accumulatedexact-diagonalization data from a separate evaluation ofall subspaces of J z except for the four subspaces of small-est dimension where | J z | is maximal or differs from itsmaximal value by 1. The number of disorder realizationsvaries, depending on system size, between several hun-dreds and 2 · . A. Heisenberg Chain
For A = 0 the central spin (cid:126)S becomes obsolete,and for I = 1 / h = 0 theresulting Heisenberg chain is integrable via the Betheansatz . However, this is a rather isolated point in thephase diagram as seen in Fig. 1 showing the disorder-averaged probability distribution (15) obtained fromexact-diagonalization data of Heisenberg chain with spinlengths I = 1 / I = 1 at different disorder strengths.For even small disorder such as h = 0 . J the systemshows ergodic statistics (4) while upon increasing h itchanges to the Poisson-type distribution (3). This tran-sition occurs for both spin lengths at about the samedisorder strength, which is a consequence of the scalingfactor 2 S = 2 I in the disorder term of the Hamiltonian(1). The fact that both spin lengths show such a simi-lar behavior is good news for semiclassical approaches tomany-body localization in spin chains .In Fig. 2 we show the expectation value (cid:104) r (cid:105) p = (cid:104) r (cid:105) s as a function of disorder strength for Heisenberg chainsof different sizes along with the standard deviation ∆ p r (top panels). The data shows a transition between theergodic phase at small h characterized by Eqs. (9), (11)to the values (5), (7) of the many-body localized phase.The bottom panels display the sample-to-sample stan-dard deviation ∆ s r according to Eqs. (25), (36). As seenfrom the figures, ∆ s r amounts only to about ten percentof ∆ p r , which demonstrates via Eq. (27) that (∆ s r ) isonly a tiny contribution to the variance (∆ p r ) . On theother hand, ∆ s r shows a pronounced maximum (∆ s r ) max which grows rapidly with system size, as displayed in theinsets of the lower panels. Moreover, in close vicinityto the corresponding position h = h max , (cid:104) r (cid:105) s . has aninflection point at h = h inf . In Fig. 3 we have plottedboth disorder strengths for I ∈ { / , } as functions ofsystems size K , which shows that both quantities seemto converge to a common value for large K . Thus, theexpectation value (cid:104) r (cid:105) s and the standard deviation ∆ s r show as a function of disorder strength typical featuresof a phase transition with the former quantity playingthe role of an order parameter.The crossing points of the data shown in the top panelsof Fig. 2 are also often considered as indications for aphase transition. Therefore, following Refs. , we alsoplot in Fig. 3 the positions h = h cr ( K ) where two curvesof (cid:104) r (cid:105) s with consecutive system sizes K and K + 1 cross. Figure 1. The probability distribution (15) for the consecutive gap ration r of a Heisenberg chain ( A = 0) of K = 18 spinsof length I = 1 / K = 11 spins of length I = 1 (right) for different disorder strength h obtained from exact-diagonalization data. At small disorder the system is ergodic and well described by the distribution (4) (red), while withincreasing h a transition to the Poisson-typed distribution (3) (green) sets in. This data set clearly deviates from h max , h inf and growsto larger disorder strengths, an observation known as the“drifting of the critical disorder strength” with systemsite . It is an interesting speculation whether h max ≈ h inf and h cr correspond, for large systems, to two distincttransitions occurring in the same systems.For the finite-size data depicted in Fig. 3 we estimatethe transition point to h max ≈ h inf ≈ . J . . . . J forspin length I = 1 /
2, and h max ≈ h inf ≈ . J . . . . J for I = 1. These values for the critical disorder strength forthe transition from the ergodic to the many-body lpcal-ized phase are somewhat smaller than the ones reportedin other works , which favor, for I = 1 / h/J ≈ concentrate on chains with an even numberof spins and the subspace with total spin J z = 0, whereashere we also take into account odd numbers of spins andall subspaces except for those with | J z | ∈ { IK, IK − } .Moreover, we introduce a new and different criterion tolocate the transition given by the position of (∆ s r ) max .Fig. 4 shows the standard deviation ∆ s r along withthe square root of “variance of the variance” (38) asa function of disorder strength for both spin lengths I = 1 / I = 1. Remarkably, both quantities arevery close to each other, especially at small disorderstrength h . This observation should be taken as an indi-cation that the underlying probability distribution s ( x )is rather narrow since both quantities become strictly equal, ∆ s r = ∆ s r = 0, for a δ -type distribution. Thisconjecture is confirmed by the data of Fig. 5 which dis-plays the probability distribution (17) for the realization-specific average (cid:104) r (cid:105) α with the disorder strengths beingthe same as in Fig. 1. The probability distribution ismuch narrower than the distribution p ( r ) and broadenssignificantly in the transition region. The latter result issimilar to an obsevation by Pal and Huse who foundnear the transition a maximum in the width of the prob-ability distribution of a long-ranged spin correlator.We also note that a closer analysis of the data ofthe lower panels of Fig. 4 suggests that ∆ s r extrapo-lates to zero for K → ∞ deep in the ergodic phase( h/J ≈
1) as well as deep in the many-body localizedphase ( h/J (cid:38) δ -function in thislimit and the above range of disorder strengths, which isconsistent with the data of Fig. 5. Note that the state-ment lim K →∞ ∆ s r = 0 implies, according to Eq. (27),that the variance (∆ p r ) ) is entirely given by the aver-aged variances within the individual disorder realizations. B. Central Spin Model
For J = 0 the central spin model resulting from theHamiltonian (1) is also integrable via an appropriateBethe ansatz and known as the Gaudin model . Com- Figure 2. Top panels: The expectation value (cid:104) r (cid:105) p as a function of disorder strength in Heisenberg chains of spin length I = 1 / I = 1 (right) for various system sizes and pertaining numbers of disorder realizations. The error bars are determinedby Eqs. (24), (25). Inset: The standard deviation ∆ p r as a function of disorder strength. The horizontal lines indicate theexpected values for GOE and Poissonian statistics as given in Eqs. (5)-(11).Bottom panels: The standard deviation ∆ s r according to Eq. (25) for the same parameters as in the top panels. For a moredetailed view on the data at large system sizes see also Fig. 4. The insets show the maximum of the standard deviation as afunction of system size where the error bars follow Eqs. (35), (38). pared to the Heisenberg chain, the coupling to the centralspin provides an alternative mechanism of introducing in-teraction among the bath spins, which are subject to arandom magnetic field.Fig. 6 shows data analogous to Fig. 1 now for centralspin models of spin length S = I = 1 / S = I = 1.For small disorder, the system clearly deviates from thePoisson-type distribution (3) and shows level repulsion.However, differently from the case of the Heisenbergchain, the level statistics do not fully reach the Gaus-sian orthogonal ensemble but change back before to anintegrable or many-body localized phase.Fig. 7 displays data for the central spin model analo-gous to Fig. 2 for the Heisenberg chain. Comparing thetop panels of both figures suggests that a transition fromthe (approximately) ergodic to the many-body localizedphase at the inflection point h = h inf (cid:46)
1, which is con-sistent with the findings of Ref. for A, h (cid:29) J . Here thefact that the transition occurs at about the same disorder strength for different system sizes depends on the scalingfactor 1 /K in front of the second term in the Hamiltonian(1). Also the sample-to-sample standard deviation ∆ s r plotted in the bottom panels of Fig. 7 behaves similarlyas for the Heisenberg chain: For large enough systemssizes this quantity developes a maximum (∆ s r ) max near h = h inf whose value increases monotonously with systemsize, as shown in the insets. Thus, we have qualitativelythe same situation as for the Heisenberg chain.Moreover, as seen in Fig. 8, the square root of the“variance of the variance” (38) follows, similarly as forthe Heisenberg chain, closely the standard deviation ∆ s r as a function of disorder strength for both spin lengths I = 1 / I = 1. This is consistent with the prob-ability distribution (17) for the realization-specific aver-age (cid:104) r (cid:105) α shown in Fig. 9. As already seen in the caseof the Heisenberg chain, the probability distribution ismuch narrower than the distribution p ( r ) and becomessignificantly broader in the transition region. Heisenberg chain, I = / h max / J h inf / J h cr / J12 14 16 18 20 222.02.53.03.54.04.5 K h ( J ) Heisenberg chain, I = h max / J h inf / J h cr / J6 8 10 12 14 16 18 20K
Figure 3. Finite-size transition data for the Heisenberg chainof spin length I = 1 / I = 1 (right). The panelsshow as a function of system size K the position h = h max of(∆ s r ) max , the position h = h inf of the inflection point of (cid:104) r (cid:105) s ,and the crossing point h = h cr ( K ) of two curves of the latterquantity with consecutive systems sizes K and K + 1. Thedashed lines are linear fits to h max and h inf , and the shadedregions estimate the transition where h max ≈ h inf Figure 4. The sample-to-sample standard deviation ∆ s r alongwith the square root of “variance of the variance” (38) as afunction of disorder strength h for Heisenberg chains of spinlength I = 1 / I = 1 (bottom panels). Thedata sets are remarkably close to each other, in particular atsmall disorder strength h IV. SUMMARY AND OUTLOOK
We have compared the transitions between ergodicand many-body localized phases in disodered Heisenbergchains as well as central spin models composed of spinsof length 1 / s r of the ex-pectation value (cid:104) r (cid:105) α of the consecutive-gap ratio in anindividual disorder realization (sample) α . This quan-tity assumes, for both types of systems and spin lengths,a maximum as a function of disorder strength, accom-panied by an inflection point of (cid:104) r (cid:105) . These are typicalfeatures of a phase transition where the latter quantityplay the role of an order parameter. The critical disorderstrength deduced from these observations turn out to besmaller than those reported in the recent literature.Further information about the transitions is containedin the probability distribution of the expectation valueswithin a given disorder realization. We expect the studyof this probability distribution and its moments to be auseful tool in the investgation of phenomena related tomany-body localization also in other systems. ACKNOWLEDGMENTS
We thank F. Evers and F. G¨ohmann for useful dis-cussions, and M. Trivelato for collaboration on an ear-lier stage of this project. J.S. acknowledges supportby FAPESP and the hospitaliy of the University of SaoPaulo at Sao Carlos and of IIP Natal. J.C.E. acknowl-edges support from the Sao Paulo Research Foundation(FAPESP) Grants No. 2016/08468-0, No. 2018/19017-4,No. 2020/00841-9, and and from Conselho Nacional dePesquisas (CNPq), Grant No. 306122/2018-9.
Figure 5. The probability distribution (17) for the realization-specific average (cid:104) r (cid:105) α for disordered Heisenberg chains of spinlength I = 1 / I = 1 (right). The disorder strengths h in both panels are the same as in Fig. 1.Figure 6. The probability distribution (15) for the consecutive gap ration r of a central spin model ( J = 0) of N = K + 1 = 11spins of length S = I = 1 / N = K + 1 = 10 spin of length S = I = 1 (right) for different disorder strength h .Similarly as for the Heisenberg chain, the data follows for small but finite disorder approximately the GOE distribution (4)(red), while with increasing h a transition to the Poisson-typed distribution (3) (green) sets in. Figure 7. Top panels: The expectation value (cid:104) r (cid:105) p as a function of disorder strength in central spin modelss of spin length S = I = 1 / S = I = 1 (right) for various system sizes and pertaining numbers of disorder realizations. The errorbars are determined by Eqs. (24), (25). Inset: The standard deviation ∆ p r as a function of disorder strength. The horizontallines indicate the expected values for GOE and Poissonian statistics as given in Eqs. (5)-(11).Bottom panels: The standard deviation ∆ s r according to Eq. (25) for the same parameters as in the top panels. For a moredetailed view on the data at large system sizes see also Fig. 8. The insets show the maximum of the standard deviation as afunction of system size.Figure 8. The sample-to-sample standard deviation ∆ s r alongwith the square root of “variance of the variance” (38) as afunction of disorder strength h for central spin models of spinlength S = I = 1 / S = I = 1 (bottompanels). Figure 9. The probability distribution (17) for the realization-specific average (cid:104) r (cid:105) α for disordered cenral spin systems of spinlength S = I = 1 / S = I = 1 (right). The disorder strengths h in both panels are the same as in Fig. 6. R. Nandkishore and D. A. Huse, Ann. Rev. Cond. Mat.Phys. , 15 (2015). E. Altman and R. Vosk, Ann. Rev. Cond. Mat. Phys. ,383 (2015). J. Z. Imbrie, V. Ros, and A. Scardicchio, Ann. Phys.(Berlin) , 1600278 (2017). K. Agarwal, E. Altma, E. Demler, S. Gopalakrishnan,D. A. Huse, and M. Knap, Ann. Phys. (Berlin) ,1600326 (2017). D. J. Luitz and Y. Bar Lev, Ann. Phys. (Berlin) ,1600350 (2017). A. Haldar and A. Das, Ann. Phys. (Berlin) , 1600333(2017). D. A. Abanin and Z. Papic, Ann. Phys. (Berlin) ,1700169 (2017). D. A. Abanin, E. Altman, I. Bloch, and M. Serbyn,arXiv:1804.11065. J. M. Deutsch, Phys. Rev. A , 2046 (1991). M. Srednicki, Phys. Rev. E
888 (1994). P. W. Anderson, Phys. Rev. , 1492 (1958). F. Evers and A. D. Mirlin, Rev. Mod. Phys. , 1355(2008). L F. Santos, G. Rigolin, and C. O. Escobar, Phys. Rev. A , 042304 (2004), A. Pal and D. A. Huse, Phys. Rev. B. , 174411 (2010). J. H. Bardarson, F. Pollmann, and J. E. Moore, Phys. Rev.Lett. , 017202 (2012). D. J. Luitz, N. Laflorencie, and F. Alet, Phys. Rev. B ,081103 (2015). A. Chandran, I. H. Kim, G. Vidal, and D. A. Abanin,Phys. Rev. B , 085425 (2015). K. Agarwal, S. Gopalakrishnan, M. Knap, M. Muller, andE. Demler, Phys. Rev. Lett. , 160401 (2015). E. Baygan, S. P. Lim, and D. N. Sheng, Phys. Rev. B ,195153 (2015). T. Devakul and R. R. P. Singh, Phys. Rev. Lett. ,187201 (2015). D. J. Luitz, N. Laflorencie, and F. Alet, Phys. Rev. B ,060201 (2016). D. J. Luitz, Phys. Rev. B , 134201 (2016). S. D. Geraedts, R. M. Nandkishore, and N. Regnault,Phys. Rev. B V. Khemani, A. Lazarides, R. Moessner, and S. L. Sondhi,Phys. Rev. Lett. , 250401 (2016). S. P. Lim and D. N. Sheng, Phys. Rev. B , 045111 (2016). M. Serbyn, A. A. Michailidis, D. A. Abanin, and Z. Papic,Phys. Rev. Lett. , 160601 (2016). T. Enss, F. Andratschko, and J. Sirker, Phys. Rev. B ,045121 (2017). O. L. Acevedo, A. Safavi-Naini, J. Schachenmayer, M. L.Wall, R. Nandkishore, and A. M. Rey, Phys. Rev. A , 033604 (2017). V. Khemani, D. N. Sheng, and D. A. Huse, Phys. Rev.Lett. , 075702 (2017). J. L. C. da C. Filho, A Saguia, L. F. Santos, and M. S.Sarandy, Phys. Rev. B , 014204 (2017). S. D. Geraedts, N. Regnault, and R. M. Nandkishore, NewJ. Phys. , 113021 (2017). K. Xu, J.-J. Chen, Y. Zeng, Y.-R. Zhang, C. Song, W.Liu, Q. Guo, P- Zhang, D. Xu, H. Deng, K. Huang, H.Wang, X. Zhu, D. Zheng, and H. Fan, Phys. Rev. Lett. , 050507 (2018). E. V. H. Doggen, F. Schindler, K. S. Tikhonov, A. D. Mir-lin, T. Neupert, D. G. Polyakov, and I. V. Gornyi, Phys.Rev. B , 174202 (2018). J. Suntajs, J. Bonca, T. Prosen, and L. Vidmar,arXiv:1905.06345. R. K. Panda, A. Scardicchio, M. Schulz, S. R. Taylor, andM. Znidaric, EPL , 67003 (2020). T. Chanda, P. Sierant, and J. Zakrzewski, Phys. Rev. B , 035148 (2020). J. Suntajs, J. Bonca, T. Prosen, and L. Vidmar,arXiv:2004.01719. N. Laflorencie, G. Lemarie, and N. Mace,arXiv:2004.02861. P. Sierant, M. Lewenstein, and J. Zakrzewski,arXiv:2005.09534. S. Dhara, A. Hamma, and E. R. Mucciolo, Phys. Rev. B , 045140 (2020). T. Chanda, P. Sierant, and J. Zakrzewski,arXiv:2006:02860. R. E. Throckmorton and S. Das Sarma, arXiv:2009:04457. D. Hetterich, N. Y. Yao, M. Serbyn, F. Pollmann, and B.Trauzettel, Phys. Rev. B , 161122 (2018). T. Guhr, A. Muller-Groeling, and H. A. Weidenmuller,Phys. Rep. , 189 (1998). V. Oganesyan and D. A. Huse, Phys. Rev. B , 155111(2007). Y. Y. Atas, E. Bogomolny, O. Giraud, and G. Roux, Phys.Rev. Lett. , 084101 (2013). H. Bethe, Z. Phys. A , 205 (1931). B. Craps, M. De Clerck, D. Janssens, V. Luyten, and C.Rabideau Phys. Rev. B , 174313 (2020). K. S. Tikhonov, A. D. Mirlin, and M. A. Skvortsov, Phys.Rev. B , 220203 (2016). M. Gaudin, J. Phys. (Paris) , 1087 (1976). B. Erbe and J. Schliemann, Phys. Rev. Lett.105