Many Body Localization Transition in the strong disorder limit : entanglement entropy from the statistics of rare extensive resonances

The space of one-dimensional disordered interacting quantum models displaying a Many-Body-Localization Transition seems sufficiently rich to produce critical points with level statistics interpolating continuously between the Poisson statistics of the Localized phase and the Wigner-Dyson statistics of the Delocalized Phase. In this paper, we consider the strong disorder limit of the MBL transition, where the critical level statistics is close to the Poisson statistics. We analyse a one-dimensional quantum spin model, in order to determine the statistical properties of the rare extensive resonances that are needed to destabilize the MBL phase. At criticality, we find that the entanglement entropy can grow with an exponent $0<\alpha<1$ anywhere between the area law $\alpha=0$ and the volume law $\alpha=1$, as a function of the resonances properties, while the entanglement spectrum follows the strong multifractality statistics. In the MBL phase near criticality, we obtain the simple value $\nu=1$ for the correlation length exponent. Independently of the strong disorder limit, we explain why for the Many-Body-Localization transition concerning individual eigenstates, the correlation length exponent $\nu$ is not constrained by the usual Harris inequality $\nu \geq 2/d$, so that there is no theoretical inconsistency with the best numerical measure $\nu = 0.8 (3)$ obtained by D. J. Luitz, N. Laflorencie and F. Alet, Phys. Rev. B 91, 081103 (2015).


INTRODUCTION
The thermalization of isolated many-body quantum systems is nowadays discussed in terms of the Eigenstate Thermalization Hypothesis (E.T.H.) [1][2][3][4] (see the review [5] and references therein) : the idea is that each manybody excited eigenstate |ψ > is 'thermal', i.e. the reduced density matrix ρ A of a sub-region A corresponds the thermal density matrix where the inverse temperature β selects the correct average energy corresponding to the initial quantum state |ψ >.
Then the Von Neumann entanglement entropy of the region A with the complementary region coincides with the thermal entropy, which is extensive with respect to the volume L d A of the region A However the presence of strong disorder can prevent this thermalization and lead to the phenomenon of Many-Body-Localization (MBL) (see the recent reviews [6,7] and references therein), where the disorder-averaged entanglement entropy of Eq. 2 follows instead the area-law [8] So in the MBL phase, excited states are somewhat similar to ground states, with efficient representation via Density-Matrix-RG or Matrix Product States [9][10][11][12] and Tensor Networks [13]. The Fisher Strong Disorder Real Space RG to construct the ground states of random quantum spin models [14][15][16] has been extended into the Strong Disorder RG procedure for the unitary dynamics [17,18], and into the RSRG-X procedure in order to construct the whole set of excited eigenstates [19][20][21][22][23][24]. This construction is actually possible only because the MBL phase is characterized by an extensive number of emergent localized conserved operators [25][26][27][28][29][30][31][32], but breaks down as the MBL transition towards delocalization is approached. As a consequence, the current RG descriptions of the MBL transition are based on different types of RG rules concerning either the entanglement [33] or the resonances [34].
From the point of view of the entanglement entropy of Eq. 2, a natural question is whether at criticality, it is possible to obtain an entanglement power-law growth intermediate between the area-law and the volume law d − 1 ≤ α ≤ d, or with possibly some logarithmic corrections. Let us now concentrate on the one-dimensional case d = 1. If one assumes some standard finite-size-scaling form in the critical region in terms of the diverging correlation length ξ one obtains that the matching with the area-law of the MBL phase (Eq. 4 in d = 1) requires the divergence as whereas the matching with the volume-law of the delocalized phase in d = 1 requires a vanishing coefficient of the volume-law As discussed in detail in [35,36], one should then distinguish two possibilities : (i) if the transition is directly towards the thermal ergodic phase satisfying E.T.H., the continuous vanishing of the coefficient 1 ξ 1−α of the volume-law of Eq. 8 is actually forbidden by the strong-subadditivity property [35], and the only possibility is that the critical point is itself thermal, i.e. it should satisfy the volume law α = 1 with the finite thermal coefficient given by Eq. 3. This scenario was found numerically [37,38] and via the RG based on entanglement [33] or resonances [34]. Then the difference between the critical point and the delocalized phase is not visible in the disorder-averaged entanglement entropy, but in its variance [33,37,38] and in the dynamical properties [33,34].
(ii) if the transition is towards a delocalized non-ergodic phase, i.e. a phase satisfying the volume law, but not the E.T.H. with the thermal coefficient fixed by Eq. 3, then the continuous vanishing of the coefficient of the volume-law of Eq. 8 is possible, and the exponent α is not a priori fixed by the strong-subadditivity property (see [35] for more details). This scenario is suggested by the point of view that the MBL transition is somewhat similar to an Anderson Localization transition in the Hilbert space of 'infinite dimensionality' as a consequence of the exponential growth of the size of the Hilbert space with the volume [39][40][41][42][43][44][45].
In this paper, we consider the strong disorder limit of a one-dimensional quantum spin model in order to analyze the statistical properties of resonances that are needed to destabilize an MBL eigenstate satisfying the area-law S A = O(1) and to determine the entanglement properties of the obtained critical state. The paper is organized as follows. In section II, we introduce the one-dimensional quantum spin model and derive the entanglement spectrum in terms of the couplings. In section III, the statistical properties of the entanglement spectrum are studied in terms of Lévy variables. The behavior of the entanglement entropy is analyzed in Section IV. Finally in section V, the multifractal statistics of the entanglement spectrum is obtained. Our conclusions are summarized in section VI. Two appendices contain complementary computations.

II. SINGULAR PERTURBATION FROM THE MANY-BODY LOCALIZED PHASE
Before we introduce the MBL model, we need to explain the motivations coming from strong disorder limit of the Anderson Localization Transition. A.
Motivation : strong disorder limit of the Anderson Localization Transition At Anderson localization transitions, the critical states can be more or less multifractal (see the review [46] and references therein) : (i) for the short-ranged tight-binding model in dimension d, there is a continuous interpolation between the 'weak multifractality' regime in d = 2 + ǫ and the 'strong multifractality' in high dimension d.
(ii) for the tight-binding model with long-ranged hopping criticality is known to occur exactly at a c = d, as can be understood from the scaling of the difference between energy levels Here there is also a continuous family of critical points as a function of the amplitude V , that interpolates between the 'weak multifractality' regime for large V → +∞ and the 'strong multifractality' regime for small V → 0 [46]. From the point of view of the statistics of the energy levels, the 'weak multifractality' of wavefunctions corresponds to a level statistics close to the Wigner Dyson statistics of the delocalized phase, with a strong level repulsion and a vanishing level compressibility χ → 0, whereas the 'strong multifractality' of wavefunctions corresponds to a level statistics close to the Poisson statistics of the localized phase, with a very weak level repulsion and a level compressibility close to unity χ ≃ 1 (see more details in [46]).
In a previous work [61], we have described how these strong multifractality results for small aamplitude V in the Long-Ranged-Hopping model can be reproduced via the first-order perturbation theory with respect to the completely localized basis, provided one takes into account the broad Lévy statistics of the terms involved in the perturbative expansion.

B. Model and notations
In the present paper, we wish to adapt the above strategy for the MBL transition as follows : instead of the standard nearest-neighbor MBL models, we wish to consider a toy model where the Hamiltonian contains direct couplings between any pair of configurations in the Hilbert space, with a vanishing small amplitude and with an appropriate decay as a function of the distance in Hilbert space in order to reach an MBL critical point. This critical point will be then as close as possible to the Many-Body-Localized Phase, and can be studied via singular perturbation theory around the completely localized limit.
So we consider a one-dimensional model involving L quantum spins with the following properties. The Hamiltonian H 0 is chosen as the simplest possible Many-Body-Localized Hamiltonian where the fields h i are random variables drawn with the Gaussian distribution of variance W 2 The 2 L eigenstates are simply given by the tensor products with the random energies For instance, the ground state corresponds to the choice S i = sgn(h i ) and has the extensive energy In the following, we focus on the middle of the spectrum, where the density of states follows the Gaussian of zero mean and of variance (LW 2 ) Since there are 2 L levels, the level spacing near zero energy scales as In analogy with Eq. 9, we wish to introduce a small perturbation V that produces a direct coupling to all other (2 L − 1) states of the Hilbert space The couplings J i1,..,i k are assumed to be of small amplitude, but they should be able to produce resonances at all scales. Their scaling should thus be directly related to the level spacing of Eq. 18, and in particular they should decay exponentially in space. However the precise conditions will be discussed later, and it is clearer to write the perturbation theory in the arbitrary small couplings J i1,..,i k , first for the eigenstates, then for the corresponding reduced density matrices, and finally for the entanglement spectrum.
and the eigenstates read To simplify the notations , let us now focus on the particular state whose energy is arbitrary in the spectrum of the unperturbed Hamiltonian H 0 , as a consequence of the random fields of Eq. 13. Let us label the other (2 L − 1) states by the 1 ≤ k ≤ L positions (i 1 , .., i k ) of flipped spins with respect to this reference configuration Then the perturbed eigenstate of Eq. 21 reads , it is convenient to introduce the following basis in each region with the same notations above and Then the eigenstate of Eq. 25 can be decomposed into It is thus convenient to use instead the following perturbed basis for the region A and the similar basis for the region B Then Eq. 28 can be rewritten at first order in the perturbation as The corresponding density matrix leads to the reduced density matrix for the region A In particular, the diagonal elements involve the positive coefficients and their sum corresponds to Eq. 33

E. Entanglement spectrum
As we will see in the next section, the diagonal coefficients D i1,..,i k A of Eq. 36 and their sum Σ 1 of Eq. 37 are very broadly distributed with a Lévy law of index µ = 1/2. The physical meaning is that they are dominated by the few biggest terms corresponding to the smallest denominators 1/ 2 kA q=1 h iq + 2 kB p=1 h jp 2 that can be interpreted as very rare resonances involving spins of both regions A and B. The technical consequence is that the diagonal coefficients D i1,..,i k A of Eq. 36 and their sum Σ 1 give contributions of first order O(|J|) in the couplings. On the contrary, the off-diagonal coefficients of Eq. 35 involving two different denominators have a finite averaged value of second order O(J 2 ) (see Appendix A).
In conclusion, at first order in the couplings, the off-diagonal coefficients of Eq. 35 do not contribute, so that the diagonal elements directly represent the entanglement spectrum, i.e. the 2 LA eigenvalues of the reduced density matrix read To characterize the statistical properties of these weights, it is convenient to introduce where the sum generalizes Eq. 37. Then the Rényi entanglement entropy of index q allows to recover the usual entanglement entropy in the limit q → 1 In particular, to obtain the disorder-averaged value of the Renyi entropy of Eq. 41 one only needs to study separately the probability distributions of Σ q and Σ 1 , as done in the next section. Other statistical properties involving correlations between Σ q and Σ 1 can also be computed (see for instance Appendix B for the disorder-averaged values Y q and for the variance of the entanglement entropy S 1 ).

III. STATISTICAL PROPERTIES OF THE ENTANGLEMENT SPECTRUM
In this section, we analyze the statistical properties of the entanglement spectrum obtained in Eq. 38.
A. Probability distribution of the variable Di 1 ,..,i k A Let us first consider the probability distribution P i1,..,i k A (D i1,..,i k A ) of the positive random variable D i1,..,i k A of Eq. 36 by evaluating its Laplace transform at lowest order in the couplings The variable is a sum of (k A + k B ) Gaussian variables (Eq. 13) and is thus distributed with the Gaussian of zero mean and of variance ( So using the change of variable x = t i1,..,i k A so that Eq 44 becomes at first order in the couplings J i1,.
It is thus convenient to introduce the notation and the Lévy positive stable-law of index µ = 1/2 and of parameter Ω with its Laplace transformL The parameter Ω directly measures the weight of the singularity in t 1 2 of the Laplace transform, and thus also the weight of the power-law behavior at large D of the distribution of Eq. 50. Lévy variables [62] actually appear very often in the studies of disordered systems (see for instance [63][64][65][66] for more details on their properties).
Then Eq. 48 means that at lowest order in the couplings J, the variable D i1,..,i k A is distributed with the Lévy law L 1 B. Probability distribution of Σ1 The probability distribution P 1 (Σ 1 ) of the positive random variable Σ1 of Eq. 37 can be evaluated similarly via its Laplace transform Eq. 47 then yields that with Eq. 54 means that at lowest order in the couplings J, the variable Σ 1 is distributed with the Lévy law L 1 2 ;Ω1 (Eqs 50 and 51) of index µ = 1/2 and of parameter Ω 1 In particular, using the identity one obtains C. Probability distribution of the weight p0 The probability distribution P(p 0 ) of the weight 0 ≤ p 0 ≤ 1 of the unperturbed state (Eq 38) can be obtained directly from Eq 56 by the change of variable In particular, its averaged value 1 0 dp 0 p 0 P(p 0 ) = 1 − Ω 1 e Ω 2 1 4 +∞ Ω 1 2 is close to unity for small Ω 1 D.

Probability distribution of Σq
Let us now consider the variable Σ q of Eq. 40.
Here one has to distinguish two regions for the index q : (i) For 0 < q < 1 2 , the moment of order q of D i1,..,i k A distributed with the Lévy law of Eq. 50 is finite and reads As a consequence, the average of Σ q is also finite and given by so that at lowest order in the couplings one has (ii) For q > 1 2 , the average of Σ q is infinite, and one needs to estimate the Laplace transform of its probability distribution P q (Σ q ) as above Using Eq. 50 and the change of variable x = tD q , one obtains so that Eq. 66 becomes It is thus convenient to introduce the notation that generalizes Eq. 55, and the Lévy positive stable-law L µ;Ω of index 0 < µ < 1 and of parameter Ω defined by its Laplace transformL Eq. 68 means that for q > 1 2 , the variable Σ q is distributed with the Lévy stable law L µq;Ωq of index µ q = 1 2q and parameter Ω q that generalizes Eq 56. In particular, using the identity of Eq. 58 for Σ q , one obtains E.

Disorder-averaged values of the Renyi entropies
The disorder-averaged entanglement entropy of Eq. 43 can be obtained from the previous results : (i) in the region 0 < q < 1/2, the contribution of Eq. 65 of order |J| 2q dominates over the contribution of Eq. 58, so that the leading order reads (ii) in the region q > 1 2 one obtains using Eq. 58 and Eq. 72 and in particular in the limit q → 1 In conclusion, the parameter Ω 1 of Eq. 55 directly represents the scale of the disorder-averaged entanglement entropy of Eq. 75. As shown in Eq. B15 of Appendix B, the scale Ω 1 also governs the variance of the entanglement entropy S 1 .

IV. SCALING OF THE ENTANGLEMENT ENTROPY WITH THE LENGTH LA
In this section, we study how the scale Ω 1 (Eq. 55) of the disorder-averaged entanglement entropy of Eq. 75 depends on the length L A of the region A. Here we need to make more precise assumptions on the couplings J i1,..,i k .

A.
Statistical properties of the couplings Ji 1 ,..,i k As explained in Section II B, within the toy model of Eq. 11 that we consider for the MBL case, the couplings J i1,..,i k involved in the perturbation of Eq. 19 are the analog of the long-ranged hoppings of Eq. 9 for the Anderson model. In the Anderson case, the hopping between two points depends only on the distance, but in the MBL case, we have to choose the properties for all the couplings J i1,..,i k . To simplify the discussion, let us make the simplest possible choice and assume that the couplings J i1,..,i k are independent Gaussian random variables of zero mean J i1,..,i k = 0 (77) and of variance depending only on the spatial range r ≡ i k − i 1 , i.e. on the distance between the two extremal spins In addition, we consider the following size dependence with three parameters where b governs the exponential decay, a governs the power-law prefactor, and v is the small amplitude of the perturbation theory described in the previous sections. From the exponential decay of the level spacing of Eq. 18, one may already guess that the critical value for the exponential decay will be b c = 1, as found indeed below.
Since the perturbation V of Eq. 11 is a part of the full Hamiltonian of Eq. 19, we wish to impose that it remains extensive with respect to the number of spins. This means that the local field in the x direction on a given spin, for instance the first one should remain a finite random variable as L → +∞. Its variance may be evaluated as When i k = 1 + r is fixed, the number of ways to place the remaining q = (k − 2) points (i 2 , .., i k−1 ) into the (i k − 2) = r − 1 possible positions (i 1 + 1, .., i k − 1) is given by the binomial coefficient r−1 q , so that Eq. 80 yields So the extensivity of the perturbation requires the convergence of the sum +∞ r=1 ∆ 2 (r)2 r < +∞ (83) With the form of Eq. 79, the convergence of the sum requires either (i) b > 1/2 , or (ii) b = 1/2 with a > 1/2.

B. Study of the scale Ω1
The average value of the absolute value of the coupling distributed with Eq. 76 is also governed by the scale ∆(r) as a function of the spatial range r.
As a consequence, the scale Ω 1 of Eq. 55 can be evaluated as above, i.e. once the extreme points i 1 and j kB have been chosen, one has binomial coefficients to take into account the number of ways to place the (k A − 1) points (i 2 , .., i kA ) in the (L A − i 1 ) remaining possible positions (i 1 + 1, .., L A ), and to place the (k B − 1) points (j 1 , .., j kB −1 ) in the (j kB − L A − 1) remaining possible positions (L A + 1, .., j kB − 1) where n = L A − i 1 and m = j kB − L A − 1 represent the distances with respect to the frontier between the regions A and B.
Of course the scale Ω 1 has always a finite contribution O(1) coming from the finite distances (n, m) of order O(1) with respect to the frontier. Now we wish to evaluate the contribution of large distances (n, m). Using the identity As could have been anticipated, this leading behavior means that the sum is dominated by the regions k A ∼ n 2 and k B ∼ m 2 that maximize the binomial coefficients, i.e. the dominant resonances on the large distance (n + m) are extensive resonances involving a number of spins of order k A + k B ≃ n+m 2 . As a consequence, after the multiplication by ∆(1 + n + m) with Eq. 79, one obtains the asymptotic behavior that we need to integrate over n and m to obtain the scale Ω 1 of Eq. 86. So we arrive at the following conclusions : (i) for b > 1, there is an exponential convergence 2 −(1−b)(n+m) at large distance, so the Many-Body-Localized phase is stable and displays a finite entanglement entropy.
(ii) for b c = 1, there is no exponential factor anymore in Eq. 89, but only the power-law . We wish that the integral over m converges in order to have a well-defined thermodynamic limit L B → +∞ for the region B : this is the case for a > 1/2. Then for 1/2 < a < 3 2 , the scale Ω 1 is dominated by the contribution of the long distance

C. Entanglement growth at criticality
In summary, for the value b c = 1 in Eq. 79 for the coupling that matches exactly the exponential decay of the level spacing of Eq. 18, we have obtained that, as the parameter a of the power-law of Eq. 79 varies in the interval 1/2 < a < 3 2 , it is possible to construct a critical state with an entanglement growth governed by the exponent anywhere between the area law α = 0 and the volume law α = 1.
Of course within the present approach, we have to stop at the critical point and we have no access to the delocalized phase, but as recalled in the Introduction near Eq. 8, the values 0 < α < 1 are only possible if the transition is towards a delocalized non-ergodic phase [35].

D.
Correlation length exponent ν of the MBL phase In the MBL phase b < b c = 1, the exponential convergence 2 −(1−b)(n+m) in Eq. 89 corresponds to the correlation length that diverges near the transition b → b c = 1 with the simple correlation length exponent This value simply reflects the crossing of the exponential decay of the level spacing of Eq. 18 and of exponential decay of the couplings of Eq. 79. It coincides with the exact correlation length exponent ν loc = 1 for the Anderson localization on the Bethe lattice [69][70][71], where the Hilbert space also grows exponentially with the distance. Since this value ν = 1 does not satisfy the usual Harris [67] or Chayes-Chayes-Fisher-Spencer [68] inequality that has been rediscussed recently for the specific case of the MBL transition [36], it is important to understand why. The Harris inequality of Eq. 94 is based on the fact that in a volume L d of a disordered sample, there are of order L d random variables defining the disorder, so that there will be fluctuations of order L d 2 as a consequence of the Central Limit Theorem. In the present model for instance, the system of L spins contains L random fields h i (Eq. 13). However, the MBL transition is very different from other types of phase transitions, because it is a transition concerning an individual arbitrary eigenstate in the middle of the spectrum, in our case the state |0 > of Eq. 22. Then the full enumeration of possible resonances involve the other (2 L − 1) states of energies So in some sense, the eigenstate |0 > does not see only the L random variables h i , but it sees effectively the (2 L − 1) other random energies that are build from the L variables h i by the various choices of the spin values S i , as can be seen also in the exponentially small level spacing of Eq. 18. As a consequence, independently of the strong disorder limit considered in the present paper, we believe that the correlation length exponent ν has no reason to satisfy the inequality of Eq. 94, so that there is no theoretical inconsistency in the numerical works where the correlation length exponent ν is found to violate the Harris inequality [37,38].

V. MULTIFRACTAL STATISTICS OF THE ENTANGLEMENT SPECTRUM
As in Anderson transitions where the critical point is characterized by multifractal eigenfunctions [46], one expects that the MBL transition is related to some multifractal properties [45,[74][75][76][77][78][79]. Besides the entanglement entropy described above, it is thus interesting to characterize the statistics of the whole entanglement spectrum of Eq. 38, both in the Many-Body-Localized phase b > 1 and at the critical point b c = 1.

A.
Multifractal exponents τ (q) In the region 0 < q < 1/2, the disorder-averaged value of the Renyi entropy of Eq. 73 can be evaluated by generalizing the previous calculations of Section IV B to obtain For the bipartite case L A = L B = L, this yields the following power-laws in terms of the size M = 2 2L of the Hilbert space with the multifractal exponents In the region q > 1 2 , Eq. 74 yields in continuity with the result of Eq. 98.
More recently, the strong multifractality spectrum of Eq. 101 has been found to describe the statistics of matrix elements of local operators and the statistics of hybridization Ratios in nearest-neighbors Many-Body-Localized models at criticality [45]. Our conclusion is thus that the present toy model is not so far from more realistic short-ranged models, since it is able to reproduce the same critical multifractal spectrum.

C.
Multifractal spectrum in the Many-Body-Localized Phase b > bc = 1 In the many-Body-Localized phase b > b c = 1, the multifractal exponents of Eqs 98 and 99 is the Legendre transform of the multifractal spectrum that is a very simple deformation of the critical result of Eq. 101. The fact that multifractality occurs also in the Many-Body-Localized Phase is in agreement with the analysis of matrix elements of local operators and the statistics of hybridization Ratios in nearest-neighbors Many-Body-Localized models [45].

VI. CONCLUSIONS
In this paper, in analogy with the strong disorder limit of the Anderson transition for a single particle (recalled in II A), we have considered the strong disorder limit of the MBL transition, defined as the limit where the level statistics at the MBL critical point is close to the Poisson statistics of the MBL phase. For a quantum one-dimensional toy model, we have analyzed the statistical properties of the rare extensive resonances that are needed to destabilize the Many-Body Localized phase. At criticality, we have found that the entanglement entropy can grow with an exponent 0 < α = 3/2 − a < 1 anywhere between the area law α = 0 and the volume law α = 1, as a function of the power-law exponent a of the couplings (Eq. 79), while the entanglement spectrum follows the strong multifractality statistics of Eq. 101, well-known as the 'strong-multifractality' spectrum in the context of Anderson Localization Transition [46], and found recently for nearest-neighbor MBL models at criticality [45].
For an initial short-ranged model, we thus expect that the important extended rare resonances are described by some effective renormalized couplings J i1,...,i k and it would be interesting to compute their properties in terms of the initial parameters, either via the forward approximation [32,39,43,72,73] or via some RG procedure.
The main difference between the present approach and the existing RG on resonances [34] is as follows : (i) here, in the strong disorder limit, we have considered the resonances concerning a single eigenstate in the middle of the spectrum.
(ii) on the contrary, Ref. [34] is based on the notion of resonant blocks, so that when two blocks are declared to be resonant, all eigenstates of the two blocks are strongly mixed and exhibit level repulsion. So we feel that this assumption should be valid in the opposite weak-disorder limit of the transition, defined as the limit where the level statistics at the critical point is close to the Wigner-Dyson statistics of the delocalized phase.
So we believe that (i) and (ii) are actually the two extreme theories of a continuous family of MBL critical points, where the level statistics interpolates between the Poisson statistics and the Wigner-Dyson statistics.
Finally, in the MBL phase near criticality, we have obtained the simple value ν = 1 for the correlation length exponent, that simply reflects the crossing of the exponential decay of the level spacing and of the exponential decay of the couplings. More generally, and independently of the strong disorder limit, we have explained why the correlation length exponent ν of the MBL transition has no reason to satisfy the usual Harris inequality ν ≥ 2/d, so that there is actually no theoretical inconsistency in the numerical works where the correlation length exponent ν is found to violate the Harris inequality [37,38].
Let us consider the probability distribution P i1,..,i k A ;i ′ 1 ,..,i ′ k ′ A (R(i 1 , .., i kA ; i ′ 1 , .., i ′ The variables are three independent Gaussian variables (Eq. 13) of zero mean and of variances 4W 2 k A , 4W 2 k ′ A and 4W 2 k B respectively. As a consequence, the probability distribution Q(z) of the variable z = (E A + E B )(E ′ A + E B ) involved in Eq. A1 reads using the change of variables y = E A + E B and in terms of the Bessel function K 0 (u) displaying the logarithmic singularity near the origin and an exponential decay for large u.
As a consequence, the average of the inverse of the variable z remains finite ) has the finite average value which is of second order O(J 2 ) with respect to the couplings. As a consequence, at the first order in the couplings considered in the main text, the off-diagonal elements are negligible.

Appendix B: Disorder-averaged values Yq
Using the identity the disorder-averaged value of Eq. 39 can be decomposed into the two contributions Eq. 54 yields that the first contribution reads To evaluate the second contribution, one needs to consider separately the two cases 0 < q < 1 2 and q > 1 2 .

1.
Region 0 < q < 1 2 For 0 < q < 1 2 , the average of Σ q is finite (Eq. 64) so that the second contribution at leading order Y q<1/2 | second = Σ q = Γ 1 2 − q 4 q √ π L kA=1 1≤i1<i2..<i k A ≤L Ω 2q i1,..,i k A is bigger than the correction of order Ω 1 = O(J) of the first contribution of Eq. B3, so that the sum of the two contributions reads For q > 1 2 , the average of Σ q is infinite, so one needs to evaluate the divergence for small t of The sum with the first contribution of Eq. B3 yields the disorder-average of Y q in the region q > 1/2 This disorder-averaged value is thus different from the typical value computed in Eq. 74 In the vicinity of q = 1, the difference found above between the averaged value of Eq. B11 Y q> 1 2 = e (q−1)Sq = 1 + (q − 1)S q + (q − 1) 2 2 S 2 q + O(q − 1) 3 (B13) and the typical value of Eq. B12 Y typ q>1/2 = e (1−q)Sq = 1 + (q − 1)S q + (q − 1) 2 2 (S q ) 2 + O(q − 1) 3 (B14) yields the variance of the entanglement entropy S 1 So the scale Ω 1 that governs the scale of the disorder-averaged entanglement entropy (Eq. 75) also determines its variance (Eq. B15).