# Many-body-localization transition in the strong disorder limit: entanglement entropy from the statistics of rare extensive resonances

Cécile Monthus

► To cite this version:

Cécile Monthus. Many-body-localization transition in the strong disorder limit: entanglement entropy from the statistics of rare extensive resonances. Entropy, 2016, 18, pp.122. 10.3390/e18040122. cea-01321700v2

HAL Id: cea-01321700

https://cea.hal.science/cea-01321700v2

Submitted on 3 Oct 2022

HAL is a multi-disciplinary open access
archive for the deposit and dissemination of sci-
entific research documents, whether they are pub-
lished or not. The documents may come from
teaching and research institutions in France or
abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire HAL, est
destinée au dépôt et à la diffusion de documents
scientifiques de niveau recherche, publiés ou non,
émanant des établissements d'enseignement et de
recherche français ou étrangers, des laboratoires
publics ou privés.

<figure><img src="image_1.png" /></figure>

Distributed under a Creative Commons Attribution 4.0 International License
Article

# Many-Body-Localization Transition in the Strong Disorder Limit: Entanglement Entropy from the Statistics of Rare Extensive Resonances

Cécile Monthus¹,²,³

¹ Institut de Physique Théorique, Commissariat à l'Energie Atomique (CEA), 91191 Gif-sur-Yvette, France; cecile.monthus@cea.fr; Tel.: +33-1-69-08-74-73

² Centre National de la Recherche Scientifique (CNRS), Institut de Physique (INP), Unité Mixte de Recherches (UMR) 3681, Avenue de la Terrasse, 91190 Gif-sur-Yvette, France

³ Université Paris Saclay, 91190 Saint-Aubin, France

Academic Editor: Jay Lawrence

Received: 25 January 2016; Accepted: 29 March 2016; Published: 1 April 2016

**Abstract:** The space of one-dimensional disordered interacting quantum models displaying a many-body localization (MBL) 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 level statistics at the MBL critical point is close to the Poisson statistics. We analyze 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 \ge 2/d$, so that there is no theoretical inconsistency with the best numerical measure $\nu = 0.8(3)$ obtained by Luitz et al. (2015).

**Keywords:** many-body-localization; entanglement; random quantum spin chains

## 1. Introduction

The thermalization of isolated many-body quantum systems is nowadays discussed in terms of the Eigenstate thermalization hypothesis (E.T.H.) [1-4] (see the review [5] and the references therein): the idea is that each many-body excited eigenstate $|\psi\rangle$ is "thermal", i.e., the reduced density matrix $\rho_A$ of a sub-region A corresponds the thermal density matrix:

$$ \rho_A^{(\beta)} = \frac{e^{-\beta H_A}}{\mathrm{Tr}_A(e^{-\beta H_A})} \qquad (1) $$

where the inverse temperature $\beta$ selects the correct average energy corresponding to the initial quantum state $|\psi\rangle$. Then, the von Neumann entanglement entropy of the region A with the complementary region:

$$ S_A = -\mathrm{Tr}_A(\rho_A \ln \rho_A) \qquad (2) $$
coincides with the thermal entropy, which is extensive with respect to the volume $L_A^d$ of the region A:

$$S_A \underset{ETH}{\approx} -\mathrm{Tr}_A(\rho_A^{th} \ln \rho_A^{th}) = s^{th}(\beta)L_A^d \quad (3)$$

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 the references therein), where the disorder-averaged entanglement entropy of Equation (2) follows instead the area law [8]:

$$\overline{S_A}_{MBL} \propto L_A^{d-1} \quad (4)$$

Therefore, in the MBL phase, excited states are somewhat similar to ground states, with efficient representation via the Density Matrix Renormalization Group (DMRG) or matrix product states [9–12] and tensor networks [13]. The Fisher strong disorder real space RG to construct the ground states of random quantum spin models [14–17] has been extended into the strong disorder RG procedure for the unitary dynamics [18,19] and into the Real-Space-RG-for-Excited-States (RSRG-X) in order to construct the whole set of excited eigenstates [20–25]. This construction is actually possible only because the MBL phase is characterized by an extensive number of emergent localized conserved operators [26–33], 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 [34] or the resonances [35].

From the point of view of the entanglement entropy of Equation (2), a natural question is whether, at criticality, it is possible to obtain an entanglement power law growth:

$$\overline{S_A}_{criti} \approx L_A^{\alpha} \quad (5)$$

intermediate between the area law and the volume law $d-1 \le \alpha \le 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 $\zeta$:

$$\overline{S_A}_{FSS} \approx L_A^{\alpha} \Phi\left(\frac{L_A}{\zeta}\right) \quad (6)$$

one obtains that the matching with the area law of the MBL phase (Equation (4) in $d=1$) requires the divergence as:

$$\overline{S_A}_{critical,loc} \propto \zeta^{\alpha} \quad (7)$$

whereas the matching with the volume law of the delocalized phase in $d=1$ requires a vanishing coefficient of the volume law:

$$\overline{S_A}_{critical,deloc} \propto L_A \left(\frac{1}{\zeta}\right)^{1-\alpha} \quad (8)$$

As discussed in detail in [36,37], 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 $(\frac{1}{\zeta})^{1-\alpha}$ of the volume law of Equation (8) is actually forbidden by the strong subadditivity property [36], and the only possibility is that the critical point is itself thermal, i.e., it should satisfy the volume law $\alpha = 1$ with the finite thermal coefficient given by Equation (3). This scenario was found numerically [38,39] and via the RG based on entanglement [34] or resonances [35]. Then, the difference between the critical point and the delocalized phase is
not visible in the disorder-averaged entanglement entropy, but in its variance [34,38,39] and in the dynamical properties [34,35].

(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 Equation (3), then the continuous vanishing of the coefficient of the volume law of Equation (8) is possible, and the exponent $\alpha$ is not a priori fixed by the strong-subadditivity property (see [36] 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 [40–46].

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 $\overline{S}_A = O(1)$ and to determine the entanglement properties of the obtained critical state. The paper is organized as follows. In Section 2, we introduce the one-dimensional quantum spin model and derive the entanglement spectrum in terms of the couplings. In Section 3, 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 4. Finally in Section 5, the multifractal statistics of the entanglement spectrum is obtained. Our conclusions are summarized in Section 6. Two appendices contain complementary computations.

## 2. Singular Perturbation from the Many-Body Localized Phase

Before we introduce the MBL model, we need to explain the motivations coming from the strong disorder limit of the Anderson localization transition.

### 2.1. 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 [47] and the references therein):

(i) For the short-range tight-binding model in dimension $d$, there is a continuous interpolation between the “weak multifractality” regime in $d = 2 + \epsilon$ and the “strong multifractality” in high dimension $d$.

(ii) For the tight-binding model with long-range hopping:

$$V(r) = \frac{V}{r^d} \qquad (9)$$

criticality is known to occur exactly at $a_c = d$, as can be understood from the scaling of the difference between energy levels:

$$\Delta E(L^d) \propto L^{-d} \qquad (10)$$

Here, there is also a continuous family of critical points as a function of the amplitude $V$, which interpolates between the “weak multifractality” regime for large $V \to +\infty$ and the “strong multifractality” regime for small $V \to 0$ [47].

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 $\chi \to 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 $\chi \simeq 1$ (see more details in [47]).

The “strong multifractality” regime $V \to 0$ has been analyzed via the Levitov renormalization method that takes into account the resonances occurring at various scales [48–55] or other methods [56–62].
In a previous work [63], we have described how these strong multifractality results for small amplitude V in the long-range 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.

## 2.2. 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.

Therefore, we consider a one-dimensional model involving L quantum spins

$$H = H_0 + V \qquad (11)$$

with the following properties. The Hamiltonian $H_0$ is chosen as the simplest possible many-body localized Hamiltonian:

$$H_0 = - \sum_{i=1}^{L} h_i \sigma_i^z \qquad (12)$$

where the fields $h_i$ are random variables drawn with the Gaussian distribution of variance $W^2$:

$$G_W(h) = \frac{1}{\sqrt{2\pi W^2}} e^{-\frac{h^2}{2W^2}} \qquad (13)$$

The $2^L$ eigenstates are simply given by the tensor products:

$$|\psi_{S_1, \dots, S_L}^{(0)}\rangle \equiv |S_1\rangle \otimes |S_2\rangle \dots \otimes |S_L\rangle \qquad (14)$$

with the random energies:

$$E_{S_1, \dots, S_L}^{(0)} = - \sum_{i=1}^{L} h_i S_i \qquad (15)$$

For instance, the ground state corresponds to the choice $S_i = \text{sgn}(h_i)$ and has the extensive energy:

$$E_{GS}^{(0)} = - \sum_{i=1}^{L} |h_i| \qquad (16)$$

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$):

$$\rho_0(E) = \frac{1}{\sqrt{2\pi LW^2}} e^{-\frac{E^2}{2LW^2}} \qquad (17)$$

Since there are $2^L$ levels, the level spacing near zero energy scales as:

$$\Delta E(L) \propto L^{\frac{1}{2}} 2^{-L} \qquad (18)$$
In analogy with Equation (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:

$$V = - \sum_{k=1}^{L} \sum_{1 \le i_1 < i_2 < \dots < i_k \le L} J_{i_1 \dots i_k} \sigma_{i_1}^x \sigma_{i_2}^x \dots \sigma_{i_k}^x \quad (19)$$

The couplings $J_{i_1 \dots 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 Equation (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_{i_1 \dots i_k}$, first for the eigenstates, then for the corresponding reduced density matrices and, finally, for the entanglement spectrum.

## 2.3. First Order Perturbation Theory for the Eigenstates

At first order in the perturbation V, the eigenvalues of Equation (15) are unchanged:

$$E_{S_1, \dots, S_L}^{(1)} = E_{S_1, \dots, S_L}^{(0)} + \langle \psi_{S_1, \dots, S_L}^{(0)} | V | \psi_{S_1, \dots, S_L}^{(0)} \rangle = E_{S_1, \dots, S_L}^{(0)} \quad (20)$$

and the eigenstates read:

$$|\psi_{S_1, \dots, S_L}^{(1)}\rangle = |\psi_{S_1, \dots, S_L}^{(0)}\rangle + \sum_{\{S'_i\}} |\psi_{S'_1, \dots, S'_L}^{(0)}\rangle > \frac{\langle \psi_{S'_1, \dots, S'_L}^{(0)} | V | \psi_{S_1, \dots, S_L}^{(0)} \rangle}{E_{S_1, \dots, S_L}^{(0)} - E_{S'_1, \dots, S'_L}^{(0)}} \quad (21)$$

To simplify the notations, let us now focus on the particular state:

$$|0\rangle \equiv |S_1 = +1, S_2 = +1, \dots, S_L = +1\rangle \quad (22)$$

whose energy:

$$E_0^{(0)} = - \sum_{i=1}^{L} h_i \quad (23)$$

is arbitrary in the spectrum of the unperturbed Hamiltonian $H_0$, as a consequence of the random fields of Equation (13). Let us label the other $(2^L - 1)$ states by the $1 \le k \le L$ positions $(i_1, \dots, i_k)$ of flipped spins with respect to this reference configuration:

$$|i_1, \dots, i_k\rangle \equiv \sigma_{i_1}^x \sigma_{i_2}^x \dots \sigma_{i_k}^x |0\rangle \quad (24)$$

Then, the perturbed eigenstate of Equation (21) reads:

$$|\psi^{(1)}\rangle = |0\rangle + \sum_{k=1}^{L} \sum_{1 \le i_1 < i_2 < \dots < i_k \le L} \left( \frac{J_{i_1, i_2, \dots, i_k}}{2 \sum_{q=1}^{k} h_{i_q}} \right) |i_1, \dots, i_k\rangle \quad (25)$$
## 2.4. Reduced Density Matrix of the Region A

To evaluate the entanglement between the two regions $A = [1, L_A]$ and $B = [L_A + 1, L = L_A + L_B]$, it is convenient to introduce the following basis in each region with the same notations above:

$$ egin{aligned} |0\rangle_A &\equiv |S_1 = +1, S_2 = +1, \dots, S_{L_A} = +1\rangle \\ |i_1, \dots, i_k\rangle_A &\equiv \sigma_{i_1}^x \sigma_{i_2}^x \dots \sigma_{i_k}^x |0\rangle_A \quad \text{for } 1 \le i_1 < i_2 < \dots < i_k \le L_A \end{aligned} \quad (26) $$

and:

$$ egin{aligned} |0\rangle_B &\equiv |S_{L_A+1} = +1, S_{L+2} = +1, \dots, S_L = +1\rangle \\ |j_1, \dots, j_k\rangle_B &\equiv \sigma_{j_1}^x \sigma_{j_2}^x \dots \sigma_{j_k}^x |0\rangle_B \quad \text{for } L_A + 1 \le j_1 < j_2 < \dots < j_k \le L_A + L_B \end{aligned} \quad (27) $$

Then, the eigenstate of Equation (25) can be decomposed into:

$$ egin{aligned} |\psi^{(1)}\rangle &= |0\rangle_A \otimes |0\rangle_B \\ &+ \sum_{k=1}^{L_A} \sum_{1 \le i_1 < i_2 < \dots < i_k \le L_A} \frac{J_{i_1, i_2, \dots, i_k}}{2 \sum_{q=1}^k h_{i_q}} |i_1, \dots, i_k\rangle_A \otimes |0\rangle_B \\ &+ \sum_{k=1}^{L_B} \sum_{1 \le j_1 < j_2 < \dots < j_k \le L_B} \frac{J_{i_1, j_1, \dots, j_k}}{2 \sum_{q=1}^k h_{j_q}} |0\rangle_A \otimes |i_1, \dots, i_k\rangle_B \\ &+ \sum_{k_A=1}^{L_A} \sum_{1 \le i_1 < i_2 < \dots < i_{k_A} \le L_A} \sum_{k_B=1}^{L_B} \sum_{1 \le j_1 < j_2 < \dots < j_{k_B} \le L_B} \frac{J_{i_1, \dots, i_{k_A}, j_1, \dots, j_{k_B}}}{2 \sum_{q=1}^{k_A} h_{i_q} + 2 \sum_{p=1}^{k_B} h_{j_p}} |i_1, \dots, i_{k_A} \rangle_A \otimes |j_1, \dots, j_{k_B} \rangle_B \end{aligned} \quad (28) $$

It is thus convenient to use instead the following perturbed basis for the region A:

$$ egin{aligned} |\phi_0^A\rangle &\equiv |0\rangle_A + \sum_{k=1}^{L_A} \sum_{1 \le i_1 < i_2 < \dots < i_k \le L_A} \frac{J_{i_1, i_2, \dots, i_k}}{2 \sum_{q=1}^k h_{i_q}} |i_1, \dots, i_k\rangle_A \\ |\phi_{i_1, \dots, i_k}^A\rangle &\equiv |i_1, \dots, i_k\rangle_A - \frac{J_{i_1, i_2, \dots, i_k}}{2 \sum_{q=1}^k h_{i_q}} |0\rangle_A \quad \text{for } 1 \le i_1 < i_2 < \dots < i_k \le L_A \end{aligned} \quad (29) $$

and the similar basis for the region B:

$$ egin{aligned} |\phi_0^B\rangle &\equiv |0\rangle_B + \sum_{k=1}^{L_B} \sum_{1 \le j_1 < j_2 < \dots < j_k \le L_B + L_A} \frac{J_{j_1, j_2, \dots, j_k}}{2 \sum_{q=1}^k h_{j_q}} |j_1, \dots, j_k\rangle_B \\ |\phi_{j_1, \dots, j_k}^B\rangle &\equiv |j_1, \dots, j_k\rangle_B - \frac{J_{j_1, j_2, \dots, j_k}}{2 \sum_{q=1}^k h_{j_q}} |0\rangle_B \quad \text{for } L_A + 1 \le j_1 < j_2 < \dots < j_k \le L_A + L_B \end{aligned} \quad (30) $$

Then, Equation (28) can be rewritten at first order in the perturbation as:

$$ egin{aligned} |\psi^{(1)}\rangle &= |\phi_0^A\rangle \otimes |\phi_0^B\rangle \\ &+ \sum_{k_A=1}^{L_A} \sum_{1 \le i_1 < i_2 < \dots < i_{k_A} \le L_A} \sum_{k_B=1}^{L_B} \sum_{L_A+1 \le j_1 < j_2 < \dots < j_{k_B} \le L_A + L_B} \frac{J_{i_1, \dots, i_{k_A}, j_1, \dots, j_{k_B}}}{2 \sum_{q=1}^{k_A} h_{i_q} + 2 \sum_{p=1}^{k_B} h_{j_p}} |\phi_{i_1, \dots, i_{k_A}}^A\rangle \otimes |\phi_{j_1, \dots, j_{k_B}}^B\rangle \end{aligned} \quad (31) $$
The corresponding density matrix:

$$
\begin{equation}
\begin{split}
\rho \equiv \frac{|\psi^{(1)}\rangle\langle\psi^{(1)}|}{\langle\psi^{(1)}|\psi^{(1)}\rangle} = \frac{1}{1+\Sigma_1} (|\phi_0^A\rangle\otimes|\phi_0^B\rangle \\
+ \sum_{k_A=1}^{L_A} \sum_{1\le i_1<i_2<\dots<i_{k_A}} \sum_{k_B=1}^{L_B} \sum_{L_A+1\le j_1<j_2<\dots<j_{k_B}\le L_A+L_B} \frac{J_{i_1\dots,i_{k_A},j_1\dots,j_{k_B}}}{2\sum_{q=1}^{k_A} h_{i_q} + 2\sum_{p=1}^{k_B} h_{j_p}} |\phi_{i_1\dots,i_{k_A}}^A\rangle\otimes|\phi_{j_1\dots,j_{k_B}}^B\rangle) \\
\times (\langle\phi_0^A|\otimes\langle\phi_0^B| \\
+ \sum_{k'_A=1}^{L_A} \sum_{1\le i'_1<i'_2<\dots<i'_{k'_A}} \sum_{k'_B=1}^{L_B} \sum_{L_A+1\le j'_1<j'_2<\dots<j'_{k'_B}\le L_A+L_B} \frac{J_{i'_1\dots,i'_{k'_A},j'_1\dots,j'_{k'_B}}}{2\sum_{q=1}^{k'_A} h'_{i'_q} + 2\sum_{p=1}^{k'_B} h'_{j'_p}} \langle\phi_{i'_1\dots,i'_{k'_A}}^A\otimes\langle\phi_{j'_1\dots,j'_{k'_B}}^B|)
\end{split}
\tag{32}
\end{equation}
$$

with:

$$
\Sigma_1 \equiv \sum_{k_A=1}^{L_A} \sum_{1 \le i_1 < i_2 < \dots < i_{k_A}} \sum_{k_B=1}^{L_B} \sum_{L_A+1 \le j_1 < j_2 < \dots < j_{k_B} \le L_A + L_B} \left( \frac{J_{i_1 \dots i_{k_A}, j_1 \dots j_{k_B}}}{2 \sum_{q=1}^{k_A} h_{i_q} + 2 \sum_{p=1}^{k_B} h_{j_p}} \right)^2 \quad (33)
$$

leads to the reduced density matrix for the region A:

$$
\begin{align}
\rho_A &\equiv \mathrm{Tr}_B(\rho) = \langle \phi_0^B | \rho | \phi_0^B \rangle + \sum_{k_B=1}^{L_B} \sum_{L_A+1 \le j_1 < j_2 < \dots < j_{k_B} \le L_A + L_B} \langle \phi_{j_1 \dots j_{k_B}}^B | \rho | \phi_{j_1 \dots j_{k_B}}^B \rangle \nonumber \\
&= \frac{1}{1+\Sigma_1} |\phi_0^A\rangle\langle\phi_0^A| \nonumber \\
&\quad + \sum_{k_A=1}^{L_A} \sum_{1 \le i_1 < i_2 < \dots < i_{k_A} \le L_A} \sum_{k'_A=1}^{L_A} \sum_{1 \le i'_1 < i'_2 < \dots < i'_{k'_A} \le L_A} \frac{R(i_1, \dots, i_{k_A}; i'_1, \dots, i'_{k'_A})}{1+\Sigma_1} |\phi_{i_1 \dots i_{k_A}}^A\rangle\langle\phi_{i'_1 \dots i'_{k'_A}}^A| \tag{34}
\end{align}
$$

with:

$$
R(i_1, \dots, i_{k_A}; i'_1, \dots, i'_{k'_A}) \equiv \sum_{k_B=1}^{L_B} \sum_{L_A+1 \le j_1 < j_2 < \dots < j_{k_B} \le L_A + L_B} \left( \frac{J_{i_1 \dots i_{k_A}, j_1 \dots j_{k_B}}}{2 \sum_{q=1}^{k_A} h_{i_q} + 2 \sum_{p=1}^{k_B} h_{j_p}} \right) \left( \frac{J_{i'_1 \dots i'_{k_A}, j'_1 \dots j'_{k_B}}}{2 \sum_{q=1}^{k'_A} h'_{i'_q} + 2 \sum_{p=1}^{k'_B} h'_{j'_p}} \right) \quad (35)
$$

In particular, the diagonal elements involve the positive coefficients:

$$
D_{i_1, \dots, i_{k_A}} \equiv R(i_1, \dots, i_{k_A}; i_1, \dots, i_{k_A}) \equiv \sum_{k_B=1}^{L_B} \sum_{L_A+1 \le j_1 < j_2 < \dots < j_{k_B} \le L_A + L_B} \left( \frac{J_{i_1, \dots, i_{k_A}, j_1, \dots, j_{k_B}}}{2 \sum_{q=1}^{k_A} h_{i_q} + 2 \sum_{p=1}^{k_B} h_{j_p}} \right)^2 \quad (36)
$$

and their sum corresponds to Equation (33):

$$
\Sigma_1 \equiv \sum_{k_A=1}^{L_A} \sum_{1 \le i_1 < i_2 < \dots < i_{k_A} \le L_A} D_{i_1, \dots, i_{k_A}} \quad (37)
$$

2.5. Entanglement Spectrum

As we will see in the next section, the diagonal coefficients $D_{i_1, \dots, i_{k_A}}$ of Equation (36) and their sum $\Sigma_1$ of Equation (37) are very broadly distributed with a Lévy law of index $\mu = 1/2$. The physical meaning is that they are dominated by the few biggest terms corresponding to the smallest denominators $1/(2^{\sum_{q=1}^{k_A} h_{i_q}} + 2^{\sum_{p=1}^{k_B} h_{j_p}})^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_{i_1, \dots, i_{k_A}}$ of Equation (36) and their sum $\Sigma_1$ give contributions of first order $O(|J|)$ in the couplings. On the contrary, the off-diagonal coefficients of Equation (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 Equation (35) do not contribute, so that the diagonal elements directly represent the entanglement spectrum, i.e., the 2$^{L_A}$ eigenvalues of the reduced density matrix read:

$$ 
\begin{aligned}
p_0 &= \frac{1}{1 + \Sigma_1} \\
p_{i_1, \dots, i_{k_A}} &= \frac{D_{i_1, \dots, i_{k_A}}}{1 + \Sigma_1} \quad \text{for } 1 \le i_1 < i_2 < \dots < i_{k_A} \le L_A
\end{aligned} 
\quad (38) $$

To characterize the statistical properties of these weights, it is convenient to introduce:

$$ Y_q(L_A) = \mathrm{Tr}_A(\rho_A^q) = p_0^q + \sum_{k_A=1}^{L_A} \sum_{1 \le i_1 < i_2 < \dots < i_{k_A} \le L_A} p_{i_1, \dots, i_{k_A}}^q = \frac{1 + \Sigma_q}{(1 + \Sigma_1)^q} \quad (39) $$

where the sum:

$$ \Sigma_q = \sum_{k_A=1}^{L_A} \sum_{1 \le i_1 < i_2 < \dots < i_{k_A} \le L_A} D_{i_1, \dots, i_{k_A}}^q \quad (40) $$

generalizes Equation (37).

Then, the Rényi entanglement entropy of index q:

$$ S_q(L_A) = \frac{\ln Y_q(L_A)}{1-q} = \frac{\ln(1+\Sigma_q) - q \ln(1+\Sigma_1)}{1-q} \quad (41) $$

allows one to recover the usual entanglement entropy in the limit $q \to 1$:

$$ S_1(L_A) = -\mathrm{Tr}_A(\rho_A \ln \rho_A) = -p_0 \ln p_0 - \sum_{k_A=1}^{L_A} \sum_{1 \le i_1 < i_2 < \dots < i_{k_A} \le L_A} p_{i_1, \dots, i_{k_A}} \ln p_{i_1, \dots, i_{k_A}} \quad (42) $$

In particular, to obtain the disorder-averaged value of the Rényi entropy of Equation (41):

$$ \bar{S}_q = \frac{\overline{\ln(1 + \Sigma_q)} - q \overline{\ln(1 + \Sigma_1)}}{1 - q} \quad (43) $$

one only needs to study separately the probability distributions of $\Sigma_q$ and $\Sigma_1$, as done in the next section. Other statistical properties involving correlations between $\Sigma_q$ and $\Sigma_1$ can also be computed (see, for instance, Appendix B for the disorder-averaged values $\bar{Y}_q$ and for the variance of the entanglement entropy $S_1$).

### 3. Statistical Properties of the Entanglement Spectrum

In this section, we analyze the statistical properties of the entanglement spectrum obtained in Equation (38).
### 3.1. Probability Distribution of the Variable $D_{i_1, \dots, i_{k_A}}$

Let us first consider the probability distribution $P_{i_1, \dots, i_{k_A}}(D_{i_1, \dots, i_{k_A}})$ of the positive random variable $D_{i_1, \dots, i_{k_A}}$ of Equation (36) by evaluating its Laplace transform at lowest order in the couplings:

$$ 
\begin{aligned}
\overline{e^{-tD_{i_1, \dots, i_{k_A}}}} &= \int_0^{+\infty} dD_{i_1, \dots, i_{k_A}} P_{i_1, \dots, i_{k_A}}(D_{i_1, \dots, i_{k_A}}) e^{-tD_{i_1, \dots, i_{k_A}}} \\
&= e^{-t \sum_{k_B=1}^{L_B} \left( \sum_{A=1}^{k_A} \sum_{q=1}^{k_B} \frac{J_{i_1, \dots, i_{k_A}, j_1, \dots, j_{k_B}}}{2 \sum_{q=1}^{k_A} h_{iq} + 2 \sum_{p=1}^{k_B} h_{jp}} \right)^2} \\
&= \prod_{k_B=1}^{L_B} \prod_{L_A+1 \le j_1 < j_2 < \dots < j_{k_B} \le L_A + L_B} \left[ 1 - \left( 1 - e^{-t \left( \frac{J_{i_1, \dots, i_{k_A}, j_1, \dots, j_{k_B}}}{2 \sum_{q=1}^{k_A} h_{iq} + 2 \sum_{p=1}^{k_B} h_{jp}} \right)^2} \right)^2 \right] \\
&= 1 - \sum_{k_B=1}^{L_B} \sum_{L_A+1 \le j_1 < j_2 < \dots < j_{k_B} \le L_A + L_B} \left( 1 - e^{-t \left( \frac{J_{i_1, \dots, i_{k_A}, j_1, \dots, j_{k_B}}}{2 \sum_{q=1}^{k_A} h_{iq} + 2 \sum_{p=1}^{k_B} h_{jp}} \right)^2} \right) + o(J) \quad (44)
\end{aligned}
$$

The variable:

$$ E = 2 \sum_{q=1}^{k_A} h_{iq} + 2 \sum_{p=1}^{k_B} h_{jp} \quad (45) $$

is a sum of $(k_A + k_B)$ Gaussian variables (Equation (13)) and is thus distributed with the Gaussian of zero mean and of variance $(4W^2(k_A + k_B)):

$$ G_{4W^2(k_A+k_B)}(E) = \frac{1}{2\sqrt{2(k_A+k_B)\pi W^2}} e^{-\frac{E^2}{8(k_A+k_B)W^2}} \quad (46) $$

Therefore, using the change of variable $x = t_{i_1, \dots, i_{k_A}} \frac{J_{i_1, \dots, i_{k_B}}^2}{E^2}$, one obtains:

$$ 
\begin{aligned}
& \frac{1 - e^{-t \left( \frac{J_{i_1, \dots, i_{k_A}, j_1, \dots, j_{k_B}}}{2 \sum_{q=1}^{k_A} h_{iq} + 2 \sum_{p=1}^{k_B} h_{jp}} \right)^2}}{1 - e^{-t \left( \frac{J_{i_1, \dots, i_{k_A}, j_1, \dots, j_{k_B}}}{2 \sum_{q=1}^{k_A} h_{iq} + 2 \sum_{p=1}^{k_B} h_{jp}} \right)^2}} \\
&= \int_0^{+\infty} \frac{dE}{\sqrt{2(k_A + k_B)\pi W^2}} e^{-\frac{E^2}{8(k_A + k_B)W^2}} \left( 1 - e^{-t \left( \frac{J_{i_1, \dots, i_{k_A}, j_1, \dots, j_{k_B}}}{E} \right)^2} \right) \\
&= |J_{i_1, \dots, i_{k_A}, j_1, \dots, j_{k_B}}| \ t^{\frac{1}{2}} \frac{1}{2\sqrt{2(k_A + k_B)\pi W^2}} \int_0^{+\infty} \frac{dx}{x^{\frac{3}{2}}} e^{-\frac{J_{i_1, j_2, \dots, i_{k_A}, j_1, \dots, j_{k_B}}^2}{8(k_A + k_B)W^2 x}} (1 - e^{-x}) \\
&= |J_{i_1, \dots, i_{k_A}, j_1, \dots, j_{k_B}}| \ t^{\frac{1}{2}} \frac{1}{2\sqrt{2(k_A + k_B)\pi W^2}} \int_0^{+\infty} \frac{dx}{x^{\frac{3}{2}}} (1 - e^{-x}) + o(J_{i_1, \dots, i_{k_A}, j_1, \dots, j_{k_B}}) \\
&= |J_{i_1, \dots, i_{k_A}, j_1, \dots, j_{k_B}}| \ t^{\frac{1}{2}} \frac{1}{\sqrt{2(k_A + k_B)W^2}} + o(J_{i_1, \dots, i_{k_A}, j_1, \dots, j_{k_B}}) \quad (47)
\end{aligned}
$$
Therefore, Equation (44) becomes at first order in the couplings $J_{i_1, \dots, i_{k_A}, j_1, \dots, j_{k_B}}$:

$$ \begin{aligned} & \overline{e^{-t_{i_1 \dots i_{k_A}} D_{i_1 \dots i_{k_A}}}} \\ &= 1 - t^{\frac{1}{2}} \sum_{k_B=1}^{L_B} \sum_{L_A+1 \le j_1 < j_2 < \dots < j_{k_B} \le L_A + L_B} \overline{|J_{i_1 \dots i_{k_A}, j_1 \dots j_{k_B}}|} \frac{1}{\sqrt{2(k_A + k_B)W^2}} + o(J) \\ &= e^{-t^{\frac{1}{2}} \sum_{k_B=1}^{L_B} \sum_{L_A+1 \le j_1 < j_2 < \dots < j_{k_B} \le L_A + L_B} \overline{|J_{i_1 \dots i_{k_A}, j_1 \dots j_{k_B}}|} \frac{1}{\sqrt{2(k_A + k_B)W^2}} + o(J)} \end{aligned} \quad (48) $$

It is thus convenient to introduce the notation:

$$ \Omega_{i_1, \dots, i_{k_A}} \equiv \sum_{k_B=1}^{L_B} \sum_{L_A+1 \le j_1 < j_2 < \dots < j_{k_B} \le L_A + L_B} \overline{|J_{i_1, \dots, i_{k_A}, j_1, \dots, j_{k_B}}|} \frac{1}{\sqrt{2(k_A + k_B)W^2}} \quad (49) $$

and the Lévy positive stable law of index $\mu = 1/2$ and of parameter $\Omega$:

$$ L_{\frac{1}{2}; \Omega}(D) \equiv \frac{\Omega}{2\sqrt{\pi D}^{\frac{3}{2}}} e^{-\frac{\Omega^2}{4D}} \quad (50) $$

with its Laplace transform:

$$ \hat{L}_{\frac{1}{2}, \Omega}(t) \equiv \int_{0}^{+\infty} dL \hat{L}_{\frac{1}{2}; \Omega}(D)e^{-tD} = e^{-\Omega t^{\frac{1}{2}}} \quad (51) $$

The parameter $\Omega$ directly measures the weight of the singularity in $t^{\frac{1}{2}}$ of the Laplace transform, and, thus, also the weight of the power-law behavior at large $D$ of the distribution of Equation (50). Lévy variables [64] actually appear very often in the studies of disordered systems (see, for instance, [65–68] for more details on their properties).

Then, Equation (48) means that at lowest order in the couplings $J$, the variable $D_{i_1, \dots, i_{k_A}}$ is distributed with the Lévy law $L_{\frac{1}{2}; \Omega_{i_1, \dots, i_{k_A}}}$:

$$ P_{i_1, \dots, i_{k_A}}(D_{i_1, \dots, i_{k_A}}) = L_{\frac{1}{2}; \Omega_{i_1, \dots, i_{k_A}}}(D_{i_1, \dots, i_{k_A}}) \quad (52) $$

## 3.2. Probability Distribution of $\Sigma_1$

The probability distribution $P_1(\Sigma_1)$ of the positive random variable $\Sigma_1$ of Equation (37) can be evaluated similarly via its Laplace transform:

$$ \begin{aligned} \overline{e^{-t\Sigma_1}} &\equiv \int_0^{+\infty} d\Sigma_1 P_1(\Sigma_1) e^{-t\Sigma_1} \\ &= e^{-t\sum_{k_A=1}^{L_A} \sum_{i_1 < i_2 < \dots < i_{k_A} \le L_A} \sum_{k_B=1}^{L_B} \sum_{L_A+1 \le j_1 < j_2 < \dots < j_{k_B} \le L_A + L_B} \left( \frac{l_{i_1 \dots i_{k_A}, j_1 \dots j_{k_B}}}{2^{\sum_{q=1}^k h_{iq}} + 2^{\sum_{p=1}^k h_{jp}}} \right)^2} \\ &= \prod_{k_A=1}^{L_A} \prod_{1 \le i_1 < i_2 < \dots < i_{k_A} \le L_A} \prod_{k_B=1}^{L_B} \prod_{L_A+1 \le j_1 < j_2 < \dots < j_{k_B} \le L_A + L_B} \left[ 1 - \left( 1 - e^{-t \left( \frac{l_{i_1 \dots i_{k_A}, j_1 \dots j_{k_B}}}{2^{\sum_{q=1}^k h_{iq}} + 2^{\sum_{p=1}^k h_{jp}}} \right)^2} \right) \right] \\ &= 1 - \sum_{k_A=1}^{L_A} \sum_{1 \le i_1 < i_2 < \dots < i_{k_A} \le L_A} \sum_{k_B=1}^{L_B} \sum_{L_A+1 \le j_1 < j_2 < \dots < j_{k_B} \le L_A + L_B} \left( 1 - e^{-t \left( \frac{l_{i_1 \dots i_{k_A}, j_1 \dots j_{k_B}}}{2^{\sum_{q=1}^k h_{iq}} + 2^{\sum_{p=1}^k h_{jp}}} \right)^2} \right) + o(J) \end{aligned} \quad (53) $$
Equation (47) then yields that:

$$ 
\begin{aligned}
\overline{e^{-t\Sigma_1}} &= 1 - t^{\frac{1}{2}} \sum_{k_A=1}^{L_A} \sum_{1 \le i_1 < i_2 < \dots < i_{k_A} \le L_A} \sum_{k_B=1}^{L_B} \sum_{L_A+1 \le j_1 < j_2 < \dots < j_{k_B} \le L_A+L_B} \overline{|J_{i_1, \dots, i_{k_A}, j_1, \dots, j_{k_B}}|} \frac{1}{\sqrt{2(k_A + k_B)W^2}} + o(J) \\
&= e^{-\Omega_1 t^{\frac{1}{2}} + o(J)}
\end{aligned} 
\quad (54) $$

with:

$$ 
\begin{aligned}
\Omega_1 &= \sum_{k_A=1}^{L_A} \sum_{1 \le i_1 < i_2 < \dots < i_{k_A} \le L_A} \Omega_{i_1, \dots, i_{k_A}} \\
&= \sum_{k_A=1}^{L_A} \sum_{1 \le i_1 < i_2 < \dots < i_{k_A} \le L_A} \sum_{k_B=1}^{L_B} \sum_{L_A+1 \le j_1 < j_2 < \dots < j_{k_B} \le L_A+L_B} \overline{|J_{i_1, \dots, i_{k_A}, j_1, \dots, j_{k_B}}|} \frac{1}{\sqrt{2(k_A + k_B)W^2}}
\end{aligned} 
\quad (55) $$

Equation (54) means that at lowest order in the couplings $J$, the variable $\Sigma_1$ is distributed with the Lévy law $L_{1/2;\Omega_1}$ (Equations (50) and (51)) of index $\mu = 1/2$ and of parameter $\Omega_1$:

$$ P_1(\Sigma_1) = L_{1/2;\Omega_1}(\Sigma_1) \quad (56) $$

In particular, using the identity:

$$ \ln(1 + \Sigma_1) = \int_0^{+\infty} \frac{dt}{t} e^{-t} (1 - e^{-t\Sigma_1}) \quad (57) $$

one obtains:

$$ 
\begin{aligned}
\overline{\ln(1 + \Sigma_1)} &= \int_0^{+\infty} d\Sigma_1 P_1(\Sigma_1) \ln(1 + \Sigma_1) \\
&= \int_0^{+\infty} \frac{dt}{t} e^{-t} \left(1 - \overline{e^{-t\Sigma_1}}\right) \\
&= \int_0^{+\infty} \frac{dt}{t} e^{-t} \Omega_1 t^{\frac{1}{2}} + o(J) \\
&= \Omega_1 \sqrt{\pi} + o(J)
\end{aligned} 
\quad (58) $$

### 3.3. Probability Distribution of the Weight $p_0$

The probability distribution $P(p_0)$ of the weight $0 \le p_0 \le 1$ of the unperturbed state (Equation (38)):

$$ p_0 = \frac{1}{1 + \Sigma_1} \quad (59) $$

can be obtained directly from Equation (56) by the change of variable:

$$ P(p_0) = \frac{1}{p_0^2} L_{\frac{1}{2};\Omega_1} \left(\frac{1}{p_0} - 1\right) = \frac{\Omega_1 p_0^{\frac{1}{2}}}{2\sqrt{\pi}(1-p_0)^{\frac{3}{2}}} e^{-\frac{\Omega_1^2 p_0}{4(1-p_0)}} \quad (60) $$

In particular, its averaged value:

$$ \int_0^1 dp_0 p_0 P(p_0) = 1 - \Omega_1 e^{\frac{\Omega_1^2}{4}} \int_{\frac{\Omega_1}{2}}^{+\infty} dt e^{-t^2} \quad (61) $$
is close to unity for small $\Omega_1$:

$$ \int_{0}^{1} dp_0 p_0 \mathcal{P}(p_0) = 1 - \frac{\sqrt{\pi}}{2} \Omega_1 + O(\Omega_1^2) \quad (62) $$

### 3.4. Probability Distribution of $\Sigma_q$

Let us now consider the variable $\Sigma_q$ of Equation (40).

Here, one has to distinguish two regions for the index $q$:

(i) For $0 < q < \frac{1}{2}$, the moment of order $q$ of $D_{i_1, \dots, i_{k_A}}$ distributed with the Lévy law of Equation (50) is finite and reads:

$$ \overline{D_{i_1, \dots, i_{k_A}}^q} = \int_0^{+\infty} dD \frac{\Omega_{i_1, \dots, i_{k_A}}}{2\sqrt{\pi}} D^{q-\frac{3}{2}} e^{-\frac{\Omega_{i_1, \dots, i_{k_A}}^2}{4D}} = \frac{\Omega_{i_1, \dots, i_{k_A}}^{2q} \Gamma(\frac{1}{2}-q)}{4^q \sqrt{\pi}} \quad (63) $$

As a consequence, the average of $\Sigma_q$ is also finite and given by:

$$ \overline{\Sigma_q} = \sum_{k_A=1}^{L} \sum_{1 \le i_1 < i_2 < \dots < i_{k_A} \le L} \overline{D_{i_1, \dots, i_{k_A}}^q} = \frac{\Gamma(\frac{1}{2}-q)}{4^q \sqrt{\pi}} \sum_{k_A=1}^{L} \sum_{1 \le i_1 < i_2 < \dots < i_{k_A} \le L} \Omega_{i_1, \dots, i_{k_A}}^{2q} \quad (64) $$

Therefore, at lowest order in the couplings, one has:

$$ \overline{\ln(1 + \Sigma_q)} \approx \overline{\Sigma_q} = \frac{\Gamma(\frac{1}{2} - q)}{4^q \sqrt{\pi}} \sum_{k_A=1}^{L} \sum_{1 \le i_1 < i_2 < \dots < i_{k_A} \le L} \Omega_{i_1, \dots, i_{k_A}}^{2q} + o(|J|^{2q}) \quad (65) $$

(ii) For $q > \frac{1}{2}$, the average of $\Sigma_q$ is infinite, and one needs to estimate the Laplace transform of its probability distribution $P_q(\Sigma_q)$ as above:

$$ \begin{aligned} \overline{e^{-t\Sigma_q}} &\equiv \int_0^{+\infty} d\Sigma_q P_q(\Sigma_q) e^{-t\Sigma_q} = \prod_{k_A=1}^{L} \prod_{1 \le i_1 < i_2 < \dots < i_{k_A} \le L} \left[ e^{-tD_{i_1, \dots, i_{k_A}}^q} \right] \\ &= \prod_{k_A=1}^{L} \prod_{1 \le i_1 < i_2 < \dots < i_{k_A} \le L} \frac{1 - \left(1 - e^{-tD_{i_1, \dots, i_{k_A}}^q}\right)}{\overline{1 - e^{-tD_{i_1, \dots, i_{k_A}}^q}}} \\ &= 1 - \sum_{k_A=1}^{L} \sum_{1 \le i_1 < i_2 < \dots < i_{k_A} \le L} \overline{\left(1 - e^{-tD_{i_1, \dots, i_{k_A}}^q}\right)} + o(J) \end{aligned} \quad (66) $$

Using Equation (50) and the change of variable $x = tD^q$, one obtains:

$$ \begin{aligned} \overline{(1 - e^{-tD^q})} &= \int_{0}^{+\infty} dD \frac{\Omega}{2\sqrt{\pi}D^{\frac{3}{2}}} e^{-\frac{\Omega^2}{4D}} (1 - e^{-tD^q}) \\ &= \frac{\Omega}{2q\sqrt{\pi}} t^{\frac{1}{2q}} \int_{0}^{+\infty} \frac{dx}{x^{1+\frac{1}{2q}}} (1 - e^{-x}) + o(J) \\ &= \frac{\Omega\sqrt{\pi}}{\sin\left(\frac{\pi}{2q}\right) \Gamma\left(\frac{1}{2q}\right)} t^{\frac{1}{2q}} + o(J) \end{aligned} \quad (67) $$
Therefore, Equation (66) becomes:

$$ 
\begin{aligned}
\overline{e^{-t\Sigma_q}} &= 1 - \frac{\sqrt{\pi}}{\sin\left(\frac{\pi}{2q}\right) \Gamma\left(\frac{1}{2q}\right)} t^{\frac{1}{2q}} \sum_{k_A=1}^{L} \sum_{1 \le i_1 < i_2 < \dots < i_{k_A} \le L} \Omega_{i_1, \dots, i_{k_A}} + o(J) \\
&= e^{-\frac{\sqrt{\pi}\Omega_1}{\sin\left(\frac{\pi}{2q}\right) \Gamma\left(\frac{1}{2q}\right)} t^{\frac{1}{2q}}} + o(J)
\end{aligned}
\quad (68) $$

It is thus convenient to introduce the notation:

$$ \Omega_q \equiv \frac{\sqrt{\pi\Omega_1}}{\sin\left(\frac{\pi}{2q}\right) \Gamma\left(\frac{1}{2q}\right)} \quad (69) $$

that generalizes Equation (55) and the Lévy positive stable law $L_{\mu;\Omega}$ of index $0 < \mu < 1$ and of parameter $\Omega$ defined by its Laplace transform:

$$ \hat{L}_{\mu, \Omega}(t) \equiv \int_{0}^{+\infty} dD L_{\mu; \Omega}(D) e^{-tD} = e^{-\Omega t^{\mu}} \quad (70) $$

Equation (68) means that for $q > \frac{1}{2}$, the variable $\Sigma_q$ is distributed with the Lévy stable law $L_{\mu_q;\Omega_q}$ of index $\mu_q = \frac{1}{2q}$ and parameter $\Omega_q$:

$$ P_q(\Sigma_q) = L_{\frac{1}{2q}; \Omega_q}(\Sigma_q) \quad (71) $$

That generalizes Equation (56). In particular, using the identity of Equation (58) for $\Sigma_q$, one obtains:

$$ 
\begin{aligned}
\overline{\ln(1 + \Sigma_q)} &= \int_{0}^{+\infty} d\Sigma_q P_q(\Sigma_q) \ln(1 + \Sigma_q) \\
&= \int_{0}^{+\infty} \frac{dt}{t} e^{-t} \left(1 - \overline{e^{-t\Sigma_q}}\right) \\
&= \int_{0}^{+\infty} \frac{dt}{t} e^{-t} \Omega_q t^{\frac{1}{2q}} + o(J) \\
&= \Omega_q \Gamma\left(\frac{1}{2q}\right) + o(J) \\
&= \Omega_1 \frac{\sqrt{\pi}}{\sin\left(\frac{\pi}{2q}\right)} + o(J)
\end{aligned}
\quad (72) $$

### 3.5. Disorder-Averaged Values of the Rényi Entropies

The disorder-averaged entanglement entropy of Equation (43) can be obtained from the previous results:

(i) In the region $0 < q < 1/2$, the contribution of Equation (65) of order $|J|^{2q}$ dominates over the contribution of Equation (58), so that the leading order reads:

$$ \overline{S_q} = \frac{\Gamma\left(\frac{1}{2} - q\right)}{(1-q)4^q \sqrt{\pi}} \sum_{k_A=1}^{L} \sum_{1 \le i_1 < i_2 < \dots < i_{k_A} \le L} \Omega_{i_1, \dots, i_{k_A}}^{2q} + o(|J|^{2q}) \quad (73) $$
(ii) In the region $q > \frac{1}{2}$, one obtains using Equations (58) and (72):

$$ \overline{S}_q = \frac{\Omega_1 \sqrt{\pi}}{1-q} \left[ \frac{1}{\sin\left(\frac{\pi}{2q}\right)} - q \right] + o(J) \quad (74) $$

and in particular, in the limit $q \to 1$:

$$ \overline{S}_1 = \Omega_1 \sqrt{\pi} + o(J) \quad (75) $$

In conclusion, the parameter $\Omega_1$ of Equation (55) directly represents the scale of the disorder-averaged entanglement entropy of Equation (75). As shown in Equation (B15) of Appendix B, the scale $\Omega_1$ also governs the variance of the entanglement entropy $S_1$.

## 4. Scaling of the Entanglement Entropy with the Length $L_A$

In this section, we study how the scale $\Omega_1$ (Equation (55)) of the disorder-averaged entanglement entropy of Equation (75) depends on the length $L_A$ of the region A. Here, we need to make more precise assumptions on the couplings $J_{i_1, \dots, i_k}$.

### 4.1. Statistical Properties of the Couplings $J_{i_1, \dots, i_k}$

As explained in Section 2.2, within the toy model of Equation (11) that we consider for the MBL case, the couplings $J_{i_1, \dots, i_k}$ involved in the perturbation of Equation (19) are the analog of the long-range hoppings of Equation (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 of the couplings $J_{i_1, \dots, i_k}$. To simplify the discussion, let us make the simplest possible choice and assume that the couplings $J_{i_1, \dots, i_k}$ are independent Gaussian random variables:

$$ P(J_{i_1, \dots, i_k}) = \frac{1}{\sqrt{2\pi\Delta^2(i_k - i_1)}} e^{-\frac{J_{i_1, \dots, i_k}^2}{2\Delta^2(i_k - i_1)}} \quad (76) $$

of zero mean:

$$ \overline{J_{i_1, \dots, i_k}} = 0 \quad (77) $$

and of variance depending only on the spatial range $r \equiv i_k - i_1$, i.e., on the distance between the two extremal spins:

$$ \overline{J_{i_1, \dots, i_k}^2} = \Delta^2(r \equiv i_k - i_1) \quad (78) $$

In addition, we consider the following size dependence with three parameters:

$$ \Delta(r) = \frac{v}{r^d} 2^{-br} \quad (79) $$

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 Equation (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 Equation (11) is a part of the full Hamiltonian of Equation (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:

$$B_1^x = - \sum_{k=2}^{L} \sum_{1<i_2<\dots<i_k \le L} J_{1,\dots,i_k} S_{i_2}^x \dots S_{i_k}^x \qquad (80)$$

should remain a finite random variable as $L \to +\infty$. Its variance may be evaluated as:

$$\overline{(B_1^x)^2} \approx \sum_{k=2}^{L} \sum_{1<i_2<\dots<i_k \le L} \overline{J_{1,\dots,i_k}^2} = \sum_{k=2}^{L} \sum_{1<i_2<\dots<i_k \le L} \Delta^2(i_k - 1) \qquad (81)$$

When $i_k = 1 + r$ is fixed, the number of ways to place the remaining $q = (k-2)$ points $(i_2, \dots, i_{k-1})$ into the $(i_k - 2) = r - 1$ possible positions $(i_1 + 1, \dots, i_k - 1)$ is given by the binomial coefficient $\binom{r-1}{q}$, so that Equation (80) yields:

$$\overline{(B_1^x)^2} \approx \sum_{r=1} \Delta^2(r) \sum_{q=0}^{r-1} \binom{r-1}{q} = \sum_{r=1} \Delta^2(r) 2^{r-1} \qquad (82)$$

Therefore, the extensivity of the perturbation requires the convergence of the sum:

$$\sum_{r=1}^{+\infty} \Delta^2(r) 2^r < +\infty \qquad (83)$$

With the form of Equation (79), the convergence of the sum:

$$\sum_{r=1}^{+\infty} \Delta^2(r) 2^r = v^2 \sum_{r=1}^{+\infty} \frac{1}{r^{2a}} 2^{(1-2b)r} \qquad (84)$$

requires either (i) $b > 1/2$ or (ii) $b = 1/2$ with $a > 1/2$.

## 4.2. Study of the Scale $\Omega_1$

The average value of the absolute value of the coupling distributed with Equation (76):

$$\overline{|J_{i_1, \dots, i_k}|} = \int_0^{+\infty} dJ J \frac{2}{\sqrt{2\pi\Delta^2(i_k - i_1)}} e^{-\frac{j^2}{2\Delta^2(i_k - i_1)}} = \sqrt{\frac{2}{\pi}} \Delta(i_k - i_1) \qquad (85)$$

is also governed by the scale $\Delta(r)$ as a function of the spatial range $r$.

As a consequence, the scale $\Omega_1$ of Equation (55) can be evaluated as above, i.e., once the extreme points $i_1$ and $j_{k_B}$ have been chosen, one has binomial coefficients to take into account the number of ways to place the $(k_A - 1)$ points $(i_2, \dots, i_{k_A})$ in the $(L_A - i_1)$ remaining possible positions $(i_1 + 1, \dots, L_A)$, and to place the $(k_B - 1)$ points $(j_1, \dots, j_{k_B - 1})$ in the $(j_{k_B} - L_A - 1)$ remaining possible positions $(L_A + 1, \dots, j_{k_B} - 1)$:

$$\begin{aligned} \Omega_1 &= \sum_{k_A=1}^{L_A} \sum_{1 \le i_1 < i_2 < \dots < i_{k_A} \le L_A} \sum_{k_B=1}^{L_B} \sum_{L_A+1 \le j_1 < j_2 < \dots < j_{k_B} \le L_A+L_B} \overline{|J_{i_1, \dots, i_{k_A}, j_1, \dots, j_{k_B}}|} \frac{1}{\sqrt{2(k_A + k_B)W^2}} \\ &= \frac{\sqrt{2}}{\pi W} \sum_{i_1=1}^{L_A} \sum_{j_{k_B}=L_A+1}^{L_A+L_B} \Delta(j_{k_B} - i_1) \sum_{k_A=1}^{L_A-i_1+1} \sum_{k_B=1}^{L_B-L_A} \frac{1}{\sqrt{(k_A + k_B)}} \binom{L_A-i_1}{k_A-1} \binom{j_{k_B}-L_A-1}{k_B-1} \\ &= \frac{\sqrt{2}}{\pi W} \sum_{n=0}^{L_A-1} \sum_{m=0}^{L_B-1} \Delta(1+n+m) \sum_{k_A=1}^{n+1} \sum_{k_B=1}^{m+1} \frac{1}{\sqrt{(k_A + k_B)}} \binom{n}{k_A-1} \binom{m}{k_B-1} \end{aligned} \qquad (86)$$
where $n = L_A - i_1$ and $m = j_{k_B} - L_A - 1$ represent the distances with respect to the frontier between the regions A and B.

Of course, the scale $\Omega_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:

$$ \frac{1}{\sqrt{k_A + k_B}} = \frac{1}{\sqrt{\pi}} \int_0^{+\infty} du \, u^{-\frac{1}{2}} e^{-(k_A + k_B)u} \qquad (87) $$

one obtains:

$$ \begin{aligned} & \sum_{k_A=1}^{n+1} \sum_{k_B=1}^{m+1} \frac{1}{\sqrt{(k_A + k_B)}} \binom{n}{k_A - 1} \binom{m}{k_B - 1} \\ &= \frac{1}{\sqrt{\pi}} \int_0^{+\infty} du \, u^{-\frac{1}{2}} \left[ \sum_{k_A=1}^{n+1} \binom{n}{k_A - 1} e^{-k_A u} \right] \left[ \sum_{k_B=1}^{m+1} \binom{m}{k_B - 1} e^{-k_B u} \right] \\ &= \frac{1}{\sqrt{\pi}} \int_0^{+\infty} du \, u^{-\frac{1}{2}} e^{-2u} [1 + e^{-u}]^{1+n+m} \\ &= \frac{1}{\sqrt{\pi}} \int_0^{+\infty} du \, u^{-\frac{1}{2}} e^{-2u} [2 - (1 - e^{-u})]^{1+n+m} \\ &= 2^{1+n+m} \frac{1}{\sqrt{\pi}} \int_0^{+\infty} du \, u^{-\frac{1}{2}} e^{-2u} \left[ 1 - \frac{(1 - e^{-u})}{2} \right]^{1+n+m} \\ &\approx \frac{2^{1+n+m}}{n+m \rightarrow +\infty} \frac{\sqrt{2}}{\sqrt{n+m}} \end{aligned} \qquad (88) $$

As could have been anticipated, this leading behavior means that the sum is dominated by the regions $k_A \sim \frac{n}{2}$ and $k_B \sim \frac{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 \approx \frac{n+m}{2}$.

As a consequence, after the multiplication by $\Delta(1 + n + m)$ with Equation (79), one obtains the asymptotic behavior:

$$ \Delta(1 + n + m) \sum_{k_A=1}^{n+1} \sum_{k_B=1}^{m+1} \frac{1}{\sqrt{(k_A + k_B)}} \binom{n}{k_A - 1} \binom{m}{k_B - 1} \stackrel{\alpha}{\propto} \frac{v}{(n+m)^{a+\frac{1}{2}}} 2^{-(1-b)(n+m)} \qquad (89) $$

that we need to integrate over $n$ and $m$ to obtain the scale $\Omega_1$ of Equation (86).

Therefore, we arrive at the following conclusions:

(i) for $b > 1$, there is an exponential convergence $2^{-(1-b)(n+m)}$ at a 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 Equation (89), but only the power law $\frac{1}{(n+m)^{a+\frac{1}{2}}}$.

We wish that the integral over $m$ converges in order to have a well-defined thermodynamic limit $L_B \to +\infty$ for the region B: this is the case for $a > 1/2$. Then, for $1/2 < a < \frac{3}{2}$, the scale $\Omega_1$ is dominated by the contribution of the long distance:

$$ \begin{aligned} \Omega_1 \left( b_c = 1; 1/2 < a < \frac{3}{2} \right) &\propto v \sum_{n=1}^{L_A} \frac{1}{\left(a - \frac{1}{2}\right) n^{a - \frac{1}{2}}} \\ &\propto \frac{v}{\left(\frac{3}{2} - a\right) \left(a - \frac{1}{2}\right)} L_A^{\frac{3}{2} - a} \end{aligned} \qquad (90) $$
### 4.3. Entanglement Growth at Criticality

In summary, for the value $b_c = 1$ in Equation (79) for the coupling that matches exactly the exponential decay of the level spacing of Equation (18), we have obtained that, as the parameter $a$ of the power law of Equation (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:

$$0 < \alpha = \frac{3}{2} - a < 1 \qquad (91)$$

anywhere between the area law $\alpha = 0$ and the volume law $\alpha = 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 Equation (8), the values $0 < \alpha < 1$ are only possible if the transition is towards a delocalized non-ergodic phase [36].

### 4.4. Correlation Length Exponent $v$ of the MBL Phase

In the MBL phase $b < b_c = 1$, the exponential convergence $2^{-(1-b)(n+m)}$ in Equation (89) corresponds to the correlation length:

$$\zeta = \frac{1}{b-1} \qquad (92)$$

that diverges near the transition $b \to b_c = 1$ with the simple correlation length exponent:

$$v = 1 \qquad (93)$$

This value simply reflects the crossing of the exponential decay of the level spacing of Equation (18) and of the exponential decay of the couplings of Equation (79). It coincides with the exact correlation length exponent $v_{loc} = 1$ for the Anderson localization on the Bethe lattice [69–71], where the Hilbert space also grows exponentially with the distance.

Since this value $v = 1$ does not satisfy the usual Harris [72] or Chayes-Chayes-Fisher-Spencer [73] inequality:

$$v \ge \frac{2}{d} \qquad (94)$$

that has been rediscussed recently for the specific case of the MBL transition [37], it is important to understand why. The Harris inequality of Equation (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$ (Equation (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\rangle$ of Equation (22). Then, the full enumeration of possible resonances involves the other $(2^L - 1)$ states of energies:

$$E_{S_1, \dots, S_L}^{(0)} = - \sum_{i=1}^{L} h_i S_i \qquad (95)$$

Therefore, in some sense, the eigenstate $|0\rangle$ 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 Equation (18).

As a consequence, independently of the strong disorder limit considered in the present paper, we believe that the correlation length exponent $v$ has no reason to satisfy the inequality of Equation (94), so
that there is no theoretical inconsistency in the numerical works where the correlation length exponent $\nu$ is found to violate the Harris inequality [38,39].

## 5. Multifractal Statistics of the Entanglement Spectrum

As in Anderson transitions where the critical point is characterized by multifractal eigenfunctions [47], one expects that the MBL transition is related to some multifractal properties [46,74–79]. Besides the entanglement entropy described above, it is thus interesting to characterize the statistics of the whole entanglement spectrum of Equation (38), both in the many-body localized phase $b > 1$ and at the critical point $b_c = 1$.

### 5.1. Multifractal Exponents $\tau(q)$

In the region $0 < q < 1/2$, the disorder-averaged value of the Rényi entropy of Equation (73) can be evaluated by generalizing the previous calculations of Section 4.2 to obtain:

$$ 
\begin{aligned}
\bar{S}_q &= \frac{\Gamma(\frac{1}{2}-q)}{(1-q)4^q \sqrt{\pi}} \sum_{k_A=1}^{L} \sum_{1 \le i_1 < i_2 < \dots < i_{k_A} \le L} \Omega_{i_1, \dots, i_{k_A}}^{2q} + o(|J|^{2q}) \\
&\propto \sum_{n=0}^{L_A-1} \sum_{m=0}^{L_B-1} \Delta^{2q} (1+n+m) \frac{2^{n+m}}{\sqrt{n+m}} \\
&\propto \sum_{n=0}^{L_A-1} \sum_{m=0}^{L_B-1} \left[ \frac{v}{(n+m)^a} 2^{-b(n+m)} \right]^{2q} \frac{2^{n+m}}{\sqrt{n+m}} \\
&\propto v^{2q} \sum_{n=0}^{L_A-1} \sum_{m=0}^{L_B-1} \frac{2^{(n+m)(1-2qb)}}{(n+m)^{\left(\frac{1}{2}+2qa\right)}}
\end{aligned} 
\quad (96) $$

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 for the generalized moments:

$$ Y_q^{\text{typ}}(L) = e^{(1-q)\bar{S}_q} \propto v^{2q} M^{-\tau_b(q)} \quad (97) $$

with the multifractal exponents:

$$ 
\begin{aligned}
\tau_b(q) &= 2qb - 1 \quad \text{for } q \le \frac{1}{2b} \\
\tau_b(q) &= 0 \quad \text{for } \frac{1}{2b} \le q \le \frac{1}{2}
\end{aligned} 
\quad (98) $$

In the region $q > 1/2$, Equation (74) yields:

$$ \tau_b(q) = 0 \quad \text{for } q \ge \frac{1}{2} \quad (99) $$

in continuity with the result of Equation (98).

### 5.2. Multifractal Spectrum $f^{\text{criti}}(\alpha)$ at the Critical Point $b_c = 1$

At the critical point $b_c = 1$, the multifractal exponents of Equations (98) and (99):

$$ 
\begin{aligned}
\tau^{\text{criti}}(q) &= 2q - 1 \quad \text{for } q \le \frac{1}{2} \\
\tau^{\text{criti}}(q) &= 0 \quad \text{for } \frac{1}{2} \le q
\end{aligned} 
\quad (100) $$
correspond to the Legendre transform of the multifractal spectrum:

$$f^{\text{criti}}(\alpha) = \frac{\alpha}{2} \quad \text{for } 0 \le \alpha \le 2 \qquad (101)$$

that is well known as the “strong-multifractality” spectrum in the context of Anderson localization transition in the limit of infinite dimension $d \to +\infty$ [47] or in the long-ranged power law Anderson model in one-dimension of Equation (9) [48–63].

More recently, the strong multifractality spectrum of Equation (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 [46]. Our conclusion is thus that the present toy model is not so far from more realistic short-range models, since it is able to reproduce the same critical multifractal spectrum.

### 5.3. Multifractal Spectrum in the Many-Body Localized Phase $b > b_c = 1$

In the many-body localized phase $b > b_c = 1$, the multifractal exponents of Equations (98) and (99):

$$\begin{aligned} \tau_b^{\text{loc}}(q) &= 2qb - 1 && \text{for } q \le \frac{1}{2b} \\ \tau_b^{\text{loc}}(q) &= 0 && \text{for } \frac{1}{2b} \le q \end{aligned} \qquad (102)$$

is the Legendre transform of the multifractal spectrum:

$$f_b^{\text{loc}}(\alpha) = \frac{\alpha}{2b} \quad \text{for } 0 \le \alpha \le 2b \qquad (103)$$

that is a very simple deformation of the critical result of Equation (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 [46].

## 6. Conclusions

In this paper, in analogy with the strong disorder limit of the Anderson transition for a single particle (recalled in Section 2.1), 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 < \alpha = 3/2 - a < 1$ anywhere between the area law $\alpha = 0$ and the volume law $\alpha = 1$, as a function of the power law exponent $a$ of the couplings (Equation (79)), while the entanglement spectrum follows the strong multifractality statistics of Equation (101), well-known as the “strong-multifractality” spectrum in the context of Anderson localization transition [47] and found recently for nearest-neighbor MBL models at criticality [46].

For an initial short-ranged model, we thus expect that the important extended rare resonances are described by some effective renormalized couplings $J_{i_1,...,i_k}$, and it would be interesting to compute their properties in terms of the initial parameters, either via the forward approximation [33,40,44,80–82] or via some RG procedure.

The main difference between the present approach and the existing RG on resonances [35] 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, [35] 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. Therefore, 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.

Therefore, 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 $v = 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 $v$ of the MBL transition has no reason to satisfy the usual Harris inequality $v \ge 2/d$, so that there is actually no theoretical inconsistency in the numerical works where the correlation length exponent $v$ is found to violate the Harris inequality [38,39].

**Acknowledgments:** It is a pleasure to thank Romain Vasseur for an illuminating discussion on various aspects of the MBL transition and, in particular, for his explanations concerning the RG procedure on resonances [35].

**Conflicts of Interest:** The author declares no conflict of interest.

# Appendix

## Appendix A. Off-Diagonal Elements of the Reduced Density Matrix

In this Appendix, we explain why the off-diagonal elements of Equation (35) of the reduced density matrix of Equation (34) can be neglected at first order in the couplings.

Let us consider the probability distribution $P_{i_1, \dots, i_{k_A} i'_1, \dots, i'_{k'_A}} (R(i_1, \dots, i_{k_A}; i'_1, \dots, i'_{k'_A}))$ of the variable $R(i_1, \dots, i_{k_A}; i'_1, \dots, i'_{k'_A})$ of Equation (35) when the indices $(i_1, \dots, i_{k_A}; i'_1, \dots, i'_{k'_A})$ are all distinct (when some indices are the same, the generalization is straightforward). The Fourier transform at lowest order in the couplings reads:

$$ 
\begin{aligned}
\overline{e^{itR(i_1, \dots, i_{k_A} i'_1, \dots, i'_{k'_A})}} &= \int_{-\infty}^{+\infty} dR (i_1, \dots, i_{k_A}; i'_1, \dots, i'_{k'_A}) P_{i_1, \dots, i_{k_A} i'_1, \dots, i'_{k'_A}} (R(i_1, \dots, i_{k_A}; i'_1, \dots, i'_{k'_A})) e^{itR(i_1, \dots, i_{k_A} i'_1, \dots, i'_{k'_A})} \\
&= e^{it \sum_{k_B=1}^{L_B} \left( \sum_{A+1 \le j_1 < j_2 < \dots < j_{k_B} \le L_A + L_B} \left( \frac{i_{1, \dots, k_A} j_1 \dots j_{k_B}}{2 \sum_{q=1}^{k_A} h_{iq} + 2 \sum_{p=1}^{k_B} h_{jp}} \right) \left( \frac{i'_{1, \dots, k_A} j_1 \dots j_{k_B}}{2 \sum_{q=1}^{k'_A} h'_{iq} + 2 \sum_{p=1}^{k'_B} h'_{jp}} \right) \right)} \\
&= \prod_{k_B=1}^{L_B} \prod_{L_A+1 \le j_1 < j_2 < \dots < j_{k_B} \le L_A + L_B} \left[ 1 - e^{it \left( \frac{i_{1, \dots, k_A} j_1 \dots j_{k_B}}{2 \sum_{q=1}^{k_A} h_{iq} + 2 \sum_{p=1}^{k_B} h_{jp}} \right) \left( \frac{i'_{1, \dots, k_A} j_1 \dots j_{k_B}}{2 \sum_{q=1}^{k'_A} h'_{iq} + 2 \sum_{p=1}^{k'_B} h'_{jp}} \right)} \right] \\
&= 1 - \sum_{k_B=1}^{L_B} \sum_{L_A+1 \le j_1 < j_2 < \dots < j_{k_B} \le L_A + L_B} \left( 1 - e^{it \left( \frac{i_{1, \dots, k_A} j_1 \dots j_{k_B}}{2 \sum_{q=1}^{k_A} h_{iq} + 2 \sum_{p=1}^{k_B} h_{jp}} \right) \left( \frac{i'_{1, \dots, k_A} j_1 \dots j_{k_B}}{2 \sum_{q=1}^{k'_A} h'_{iq} + 2 \sum_{p=1}^{k'_B} h'_{jp}} \right)} \right) + \dots \quad (A1)
\end{aligned}
$$

The variables:

$$ 
\begin{aligned}
E_A &= 2 \sum_{q=1}^{k_A} h_{iq} \\
E'_A &= 2 \sum_{q=1}^{k'_A} h'_{iq} \\
E_B &= 2 \sum_{p=1}^{k_B} h_{jp}
\end{aligned}
\qquad (A2) 
$$
are three independent Gaussian variables (Equation (13)) of zero mean and of variances $4W^2k_A$, $4W^2k'_A$ and $4W^2k_B$, respectively. As a consequence, the probability distribution $Q(z)$ of the variable $z = (E_A + E_B)(E'_A + E_B)$ involved in Equation (A1) reads using the change of variables $y = E_A + E_B$ and $y' = E_A + E_B$:

$$ 
\begin{aligned}
Q(z) &= \int_{-\infty}^{+\infty} \frac{dE_B e^{-\frac{E_B^2}{8k_B W^2}}}{2\sqrt{2k_B \pi W^2}} \int_{-\infty}^{+\infty} \frac{dE_A e^{-\frac{E_A^2}{8k_A W^2}}}{2\sqrt{2k_A \pi W^2}} \int_{-\infty}^{+\infty} \frac{dE'_A e^{-\frac{(E'_A)^2}{8k'_A W^2}}}{2\sqrt{2k'_A \pi W^2}} \delta(z - (E_A + E_B)(E'_A + E_B)) \\
&= \int_{-\infty}^{+\infty} \frac{dE_B e^{-\frac{E_B^2}{8k_B W^2}}}{2\sqrt{2k_B \pi W^2}} \int_{-\infty}^{+\infty} \frac{dy e^{-\frac{(y-E_B)^2}{8k_A W^2}}}{2\sqrt{2k_A \pi W^2}} \int_{-\infty}^{+\infty} \frac{dy' e^{-\frac{(y'-E_B)^2}{8k'_A W^2}}}{2\sqrt{2k'_A \pi W^2}} \delta(z - yy') \\
&= \int_{-\infty}^{+\infty} \frac{dy e^{-\frac{y^2}{8k_A W^2}}}{\sqrt{8k_A \pi W^2}} \int_{-\infty}^{+\infty} \frac{dy' e^{-\frac{(y')^2}{8k'_A W^2}}}{\sqrt{8k'_A \pi W^2}} \frac{\delta(y' - \frac{z}{y})}{|y|} \int_{-\infty}^{+\infty} \frac{dE_B}{\sqrt{8k_B \pi W^2}} e^{-\frac{E_B^2}{8W^2} \left( \frac{1}{k_B} + \frac{1}{k'_A} + \frac{1}{k'_A} \right)} e^{\frac{(y+y')E_B}{4k_A W^2}} \\
&= e^{\frac{z}{4W^2 \left( k_A + k'_A + \frac{k_A k'_A}{k_B} \right)}} \frac{1}{4\pi W^2 \sqrt{k_A k'_A + k_B (k_A + k'_A)}} K_0 \left( |z| \frac{\sqrt{(k_B + k_A)(k_B + k'_A)}}{4W^2 \left( k_A + k'_A + \frac{k_A k'_A}{k_B} \right)} \right) \quad (A3)
\end{aligned}
$$

in terms of the Bessel function $K_0(u)$ displaying the logarithmic singularity near the origin:

$$ K_0(u) = -\ln u + (\ln 2 - \gamma_{Euler}) + O(u^2) \quad (A4) $$

and an exponential decay for large $u$.

As a consequence, the average of the inverse of the variable $z$ remains finite:

$$ 
\begin{aligned}
\int_{-\infty}^{+\infty} \frac{dz}{z} Q(z) &= \int_{0}^{+\infty} \frac{dz}{z} [Q(z) - Q(-z)] \\
&= \int_{0}^{+\infty} \frac{dz}{z} \frac{\sinh\left(\frac{z}{4W^2\left(k_A + k'_A + \frac{k_A k'_A}{k_B}\right)}\right)}{2\pi W^2 \sqrt{k_A k'_A + k_B(k_A + k'_A)}} K_0\left(|z| \frac{\sqrt{(k_B + k_A)(k_B + k'_A)}}{4W^2\left(k_A + k'_A + \frac{k_A k'_A}{k_B}\right)}\right) \\
&= \frac{1}{4W^2 \sqrt{k_A k'_A + k_B(k_A + k'_A)}} \arcsin\left(\frac{1}{(k_B + k_A)(k_B + k'_A)}\right) \quad (A5)
\end{aligned}
$$

Therefore, Equation (A1) yields that the variable $R(i_1, ..., i_{k_A}; i'_1, ..., i'_{k'_A})$ has the finite average value:

$$ 
\begin{aligned}
\overline{R(i_1, ..., i_{k_A}; i'_1, ..., i'_{k'_A})} &= \sum_{k_B=1}^{L_B} \sum_{L_A+1 \le j_1 < \dots < j_{k_B} \le L_A+L_B} \frac{J_{i_1, \dots, i_{k_A}; j_1, \dots, j_{k_B}} J'_{i'_1, \dots, i'_{k_A}; j'_1, \dots, j'_{k_B}}}{4W^2 \sqrt{k_A k'_A + k_B(k_A + k'_A)}} \arcsin \left( \frac{1}{(k_B + k_A)(k_B + k'_A)} \right) + o(J^2) \quad (A6)
\end{aligned}
$$

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 $\overline{Y_q}$

Using the identity:

$$ \frac{1}{a^q} = \frac{1}{\Gamma(q)} \int_0^{+\infty} dt \, t^{q-1} e^{-at} \quad (B1) $$
the disorder-averaged value of Equation (39) can be decomposed into the two contributions:

$$
\begin{align*}
\overline{Y}_q &= \overline{\frac{1 + \Sigma_q}{(1 + \Sigma_1)^q}} = \overline{Y_q}_{first} + \overline{Y_q}_{second} \\
\overline{Y_q}_{first} &= \overline{\frac{1}{(1 + \Sigma_1)^q}} = \frac{1}{\Gamma(q)} \int_0^{+\infty} dt \, t^{q-1} e^{-t} \overline{e^{-t\Sigma_1}} \\
\overline{Y_q}_{second} &= \overline{\frac{\Sigma_q}{(1 + \Sigma_1)^q}} = \frac{1}{\Gamma(q)} \int_0^{+\infty} dt \, t^{q-1} e^{-t} \overline{\Sigma_q e^{-t\Sigma_1}} \tag{B2}
\end{align*}
$$

Equation (54) yields that the first contribution reads:

$$
\begin{align}
\overline{Y}_q|_{first} &= \frac{1}{\Gamma(q)} \int_0^{+\infty} dt \, t^{q-1} e^{-t} (1 - \Omega_1 t^{\frac{1}{2}} + o(J)) \nonumber \\
&= 1 - \Omega_1 \frac{\Gamma\left(q + \frac{1}{2}\right)}{\Gamma(q)} + o(J) \tag{B3}
\end{align}
$$

To evaluate the second contribution, one needs to consider separately the two cases $0 < q < \frac{1}{2}$ and $q > \frac{1}{2}$.

Appendix B.1. Region $0 < q < \frac{1}{2}$

For $0 < q < \frac{1}{2}$, the average of $\Sigma_q$ is finite (Equation (64)), so that the second contribution at leading order:

$$
\overline{Y_{q<1/2}|_{second}} = \overline{\Sigma_q} = \frac{\Gamma\left(\frac{1}{2} - q\right)}{4^q \sqrt{\pi}} \sum_{k_A=1}^{L} \sum_{1 \le i_1 < i_2 < \dots < i_{k_A} \le L} \Omega_{i_1, \dots, i_{k_A}}^{2q} + o(|J|^{2q}) \quad (B4)
$$

is bigger than the correction of order $\Omega_1 = O(J)$ of the first contribution of Equation (B3), so that the sum of the two contributions reads:

$$
\overline{Y_{q<1/2}} = 1 + \frac{\Gamma\left(\frac{1}{2} - q\right)}{4^q \sqrt{\pi}} \sum_{k_A=1}^{L} \sum_{1 \le i_1 < i_2 < \dots < i_{k_A} \le L} \Omega_{i_1, \dots, i_{k_A}}^{2q} + o(|J|^{2q}) \quad (B5)
$$

This disorder-averaged value thus coincides at leading order with the typical value computed in Equation (73):

$$
Y_{q<1/2}^{typ} \equiv e^{\overline{\ln Y_{q<1/2}}} = e^{(1-q)\overline{S_q}} = 1 + \frac{\Gamma\left(\frac{1}{2} - q\right)}{4^q \sqrt{\pi}} \sum_{k_A=1}^{L} \sum_{1 \le i_1 < i_2 < \dots < i_{k_A} \le L} \Omega_{i_1, \dots, i_{k_A}}^{2q} + o(|J|^{2q}) \quad (B6)
$$

Appendix B.2. Region $q > \frac{1}{2}$

For $q > \frac{1}{2}$, the average of $\Sigma_q$ is infinite, so one needs to evaluate the divergence for small $t$ of:

$$
\begin{align}
\overline{\Sigma_q e^{-t\Sigma_1}} &= \overline{\sum_{k_A=1}^{L_A} \sum_{1 \le i_1 < i_2 < \dots < i_{k_A} \le L_A} D_{i_1, \dots, i_{k_A}}^q e^{-t \sum_{k'_A=1}^{L_A} \left( \sum_{1 \le i'_1 < i'_2 < \dots < i'_{k'_A} \le L_A} D_{i'_1, \dots, i'_A}^{i'_1} \right)}} \nonumber \\
&= \overline{\sum_{k_A=1}^{L_A} \sum_{1 \le i_1 < i_2 < \dots < i_{k_A} \le L_A} D_{i_1, \dots, i_{k_A}}^q e^{-t D_{i_1, \dots, i_{k_A}}}} \left[ 1 + O(t^{\frac{1}{2}}) \right] \tag{B7}
\end{align}
$$
The Lévy distribution of Equations (50) and (51) for the variable $D_{i_1, \dots, i_{k_A}}$ yields the singularity:

$$ 
\begin{aligned}
\overline{D_{i_1, \dots, i_{k_A}}^q e^{-tD_{i_1, \dots, i_{k_A}}}} &= \frac{\Omega_{i_1, \dots, i_{k_A}}}{2\sqrt{\pi}} \int_0^{+\infty} dD D^{q-\frac{3}{2}} e^{-\frac{\Omega_{i_1, \dots, i_{k_A}}^2}{4D}} e^{-tD} \\
&= \frac{\Omega_{i_1, \dots, i_{k_A}}}{2\sqrt{\pi}} t^{\frac{1}{2}-q} \int_0^{+\infty} dx x^{q-\frac{3}{2}} e^{-\frac{\Omega_{i_1, \dots, i_{k_A}}^2}{4x}} t e^{-x} \\
&= \frac{\Omega_{i_1, \dots, i_{k_A}}}{2\sqrt{\pi}} t^{\frac{1}{2}-q} \Gamma\left(q-\frac{1}{2}\right) + o(|J|)
\end{aligned}
\quad (B8) $$

Therefore, Equation (B7) displays the same singularity in terms of the parameter of Equation (55):

$$ \overline{\Sigma_q e^{-t\Sigma_1}} = \frac{\Omega_1}{2\sqrt{\pi}} t^{\frac{1}{2}-q} \Gamma\left(q-\frac{1}{2}\right) + o(|J|) \quad (B9) $$

and the second contribution of Equation (B2) reads at leading order:

$$ 
\begin{aligned}
\overline{Y_{q>\frac{1}{2}}}_{\text{second}} &= \frac{1}{\Gamma(q)} \int_0^{+\infty} dt \, t^{q-1} e^{-t} \overline{\Sigma_q e^{-t\Sigma_1}} \\
&= \frac{\Omega_1}{2} \frac{\Gamma\left(q-\frac{1}{2}\right)}{\Gamma(q)} + o(|J|)
\end{aligned}
\quad (B10) $$

The sum with the first contribution of Equation (B3) yields the disorder-average of $Y_q$ in the region $q > 1/2$:

$$ 
\begin{aligned}
\overline{Y_{q>\frac{1}{2}}} &= 1 + \frac{\Omega_1}{2\Gamma(q)} \left[ \Gamma\left(q-\frac{1}{2}\right) - 2\Gamma\left(q+\frac{1}{2}\right) \right] + o(J) \\
&= 1 - \Omega_1 \frac{\Gamma\left(q-\frac{1}{2}\right)}{\Gamma(q-1)} + o(J)
\end{aligned}
\quad (B11) $$

This disorder-averaged value is thus different from the typical value computed in Equation (74):

$$ Y_{q>1/2}^{\text{typ}} \equiv e^{\overline{\ln Y_{q>1/2}}} = e^{(1-q)\overline{S}_q} = 1 + \Omega_1 \sqrt{\pi} \left[ \frac{1}{\sin\left(\frac{\pi}{2q}\right)} - q \right] + o(J) \quad (B12) $$

## Appendix B.3. Consequence for the Variance of the Entanglement Entropy

In the vicinity of $q = 1$, the difference found above between the averaged value of Equation (B11):

$$ \overline{Y_{q>\frac{1}{2}}} = e^{(q-1)\overline{S}_q} = 1 + (q-1)\overline{S}_q + \frac{(q-1)^2}{2}\overline{S}_q^2 + O(q-1)^3 \quad (B13) $$

and the typical value of Equation (B12):

$$ Y_{q>1/2}^{\text{typ}} = e^{(1-q)\overline{S}_q} = 1 + (q-1)\overline{S}_q + \frac{(q-1)^2}{2}(\overline{S}_q)^2 + O(q-1)^3 \quad (B14) $$
yields the variance of the entanglement entropy $S_1$:

$$ 
\begin{aligned}
\overline{S_1^2} - (\overline{S_1})^2 &= 2 \lim_{q \to 1} \left[ \frac{\overline{Y_{q>\frac{1}{2}}} - Y_{q>1/2}^{typ}}{(q-1)^2} \right] \\
&= 2\Omega_1 \lim_{q \to 1} \left[ \frac{\sqrt{\pi} \left[ q - \frac{1}{\sin\left(\frac{\pi}{2q}\right)} \right] - \frac{\Gamma(q-\frac{1}{2})}{\Gamma(q-1)}}{(q-1)^2} \right] + o(J) \\
&= \Omega_1 \sqrt{\pi} \left( 4 \ln 2 - \frac{\pi^2}{4} \right) + o(J) \qquad (B15)
\end{aligned}
$$

Therefore, the scale $\Omega_1$ that governs the scale of the disorder-averaged entanglement entropy (Equation (75)) also determines its variance (Equation (B15)).

## References

1. Deutsch, J.M. Quantum statistical mechanics in a closed system. Phys. Rev. A **1991**, *43*, 2046–2049.

2. Srednicki, M. Chaos and quantum thermalization. Phys. Rev. E **1994**, *50*, 888–901.

3. Rigol, M.; Dunjko, V.; Olshanii, M. Thermalization and its mechanism for generic isolated quantum systems. Nature **2008**, *452*, 854–858.

4. Goldstein, S.; Huse, D.A.; Lebowitz J.L.; Tumulka, R. Thermal equilibrium of a macroscopic quantum system in a pure state. Phys. Rev. Lett. **2015**, *115*, 100402.

5. D’Alessio, L.; Kafri, Y.; Polkovnikov A.; Rigol, M. From Quantum Chaos and Eigenstate Thermalization to Statistical Mechanics and Thermodynamics. **2015**, arxiv:1509.06411.

6. Nandkishore, R.; Huse, D.A. Many body localization and thermalization in quantum statistical mechanics. Annu. Rev. Condens. Matter Phys. **2015**, *6*, 15–38.

7. Altman, E.; Vosk, R. Universal dynamics and renormalization in many-body-localized Systems. Annu. Rev. Condens. Matter Phys. **2015**, *6*, 383–409.

8. Bauer, B.; Nayak, C. Area laws in a many-body localized state and its implications for topological order. J. Stat. Mech. **2013**, 2013, doi:10.1088/1742-5468/2013/09/P09005.

9. Pekker, D.; Clark, B.K. Encoding the structure of many-body localization with matrix product operators. **2014**, arxiv:1410.2224.

10. Yu, X.; Pekker, D.; Clark, B.K. Finding matrix product state representations of highly-excited eigenstates of many-body localized Hamiltonians. **2015**, arxiv:1509.01244.

11. Friesdorf, M.; Werner, A.H.; Brown, W.; Scholz, V.B.; Eisert, J. Many-body localization implies that eigenvectors are matrix-product states. Phys. Rev. Lett. **2015**, *114*, 170505.

12. Khemani, V.; Pollmann, F.; Sondhi, S.L. Obtaining highly-excited eigenstates of many-body localized Hamiltonians by the density matrix renormalization group. **2015**, arxiv:1509.00483.

13. Chandran, A.; Carrasquilla, J.; Kim, I.H.; Abanin D.A.; Vidal, G. Spectral tensor networks for many-body localization. Phys. Rev. B **2015**, *92*, 024201.

14. Fisher, D.S. Random antiferromagnetic quantum spin chains. Phys. Rev. B **1994**, *50*, 3799–3821.

15. Fisher, D.S. Random transverse field Ising spin chains. Phys. Rev. Lett. **1992**, *69*, 534–537.

16. Fisher, D.S. Critical behavior of random transverse-field Ising spin chains. Phys. Rev. B **1995**, *51*, 6411–6461.

17. Fisher, D.S. Phase transitions and singularities in random quantum systems. Physica A **1999**, *263*, 222–233.

18. Vosk, R.; Altman, E. Many-body localization in one dimension as a dynamical renormalization group fixed point. Phys. Rev. Lett. **2013**, *110*, 067204.

19. Vosk, R.; Altman, E. Dynamical quantum phase transitions in random spin chains. Phys. Rev. Lett. **2014**, *112*, 217204.

20. Pekker, D.; Refael, G.; Altman, E.; Demler, E.; Oganesyan, V. Hilbert-glass transition: New universality of temperature-tuned many-body dynamical quantum criticality. Phys. Rev. X **2014**, *4*, 011052.
21. Huang, Y.; Moore, J.E. Excited-state entanglement and thermal mutual information in random spin chains. *Phys. Rev. B* **2014**, *90*, 220202(R).

22. Vasseur, R.; Potter, A.C.; Parameswaran, S.A. Quantum criticality of hot random spin chains. *Phys. Rev. Lett.* **2015**, *114*, 217201.

23. Pouranvari, M.; Yang, K. Entanglement spectrum and entangled modes of highly excited states in random XX spin chains. *Phys. Rev. B* **2015**, *92*, 245134.

24. You, Y.Z.; Qi, X.L.; Xu, C. Entanglement Holographic Mapping of Many-Body Localized System by Spectrum Bifurcation Renormalization Group. **2016**, arxiv:1508.03635.

25. Monthus, C. Many-body localization: Construction of the emergent local conserved operators via block real-space renormalization. *J. Stat. Mech. 2016*, 2016, 033101.

26. Swingle, B. A simple model of many-body localization. **2013**, arxiv:1307.0507.

27. Serbyn, M.; Papić, Z.; Abanin, D.A. Local conservation laws and the structure of the many-body localized states. *Phys. Rev. Lett.* **2013**, *111*, 127201.

28. Huse, D.A.; Nandkishore, R.; Oganesyan, V. Phenomenology of fully many-body-localized systems. *Phys. Rev. B* **2014**, *90*, 174202.

29. Nanduri, A.; Kim, H.; Huse, D.A. Entanglement spreading in a many-body localized system. *Phys. Rev. B* **2014**, *90*, 064201.

30. Imbrie, J.Z. On Many-Body Localization for Quantum Spin Chains. **2014**, arxiv:1403.7837.

31. Serbyn, M.; Papić, Z.; Abanin, D.A. Quantum quenches in the many-body localized phase. *Phys. Rev. B* **2014**, *90*, 174302.

32. Chandran, A.; Kim, I.H.; Vidal, G.; Abanin, D.A. Constructing local integrals of motion in the many-body localized phase. *Phys. Rev. B* **2015**, *91*, 085425.

33. Ros, V.; Müller, M.; Scardicchio, A. Integrals of motion in the many-body localized phase. *Nucl. Phys. B* **2015**, *891*, 420-465.

34. Vosk, R.; Huse, D.A.; Altman, E. Theory of the many-body localization transition in one-dimensional systems. *Phys. Rev. X* **2015**, *5*, 031032.

35. Potter, A.C.; Vasseur, R. Universal properties of many-body delocalization transitions. *Phys. Rev. X* **2015**, *5*, 031033.

36. Grover, T. Certain general constraints on the many-body localization transition. **2014**, arxiv:1405.1471.

37. Chandran, A.; Laumann, C.R.; Oganesyan, V. Finite size scaling bounds on many-body localized phase transitions. **2015**, arxiv:1509.04285.

38. Kjäll, J.A.; Bardarson, J.H.; Pollmann, F. Many-body localization in a disordered quantum Ising chain. *Phys. Rev. Lett.* **2014**, *113*, 107204.

39. Luitz, D.J.; Laflorencie, N.; Alet, F. Many-body localization edge in the random-field Heisenberg chain. *Phys. Rev. B* **2015**, *91*, 081103.

40. Altshuler, B.L.; Gefen, Y.; Kamenev, A.; Levitov, L.S. Quasiparticle lifetime in a finite system: A nonperturbative approach. *Phys. Rev. Lett.* **1997**, *78*, 2803-2806.

41. Gornyi, I.V.; Mirlin, A.D.; Polyakov, D.G. Interacting electrons in disordered wires: Anderson localization and low-*T* transport. *Phys. Rev. Lett.* **2005**, *95*, 206603.

42. Oganesyan, V.; Huse, D.A. Localization of interacting fermions at high temperature. *Phys. Rev. B* **2007**, *75*, 155111.

43. Monthus, C.; Garel, T. Many-body localization transition in a lattice model of interacting fermions: Statistics of renormalized hoppings in configuration space. *Phys. Rev. B* **2010**, *81*, 134202.

44. De Luca, A.; Altshuler, B.L.; Kravtsov, V.E.; Scardicchio, A. Anderson localization on the Bethe lattice: Nonergodicity of extended states. *Phys. Rev. Lett.* **2014**, *113*, 046806.

45. Gornyi, I.V.; Mirlin, A.D.; Polyakov, D.G. Many-body delocalization transition and relaxation in a quantum dot. *Phys. Rev. B* **2016**, *93*, 125419.

46. Monthus, C. Many-body-localization transition: Strong multifractality spectrum at the critical point towards a non-ergodic delocalized phase. **2016**, arxiv:1603.04701.

47. Evers, F.; Mirlin, A.D. Anderson transitions. *Rev. Mod. Phys.* **2008**, *80*, 1355-1417.

48. Levitov, L.S. Absence of localization of vibrational modes due to dipole-dipole interaction. *Europhys. Lett.* **1989**, *9*, 83-86.

49. Levitov, L.S. Delocalization of vibrational modes caused by electric dipole interaction. *Phys. Rev. Lett.* **1990**, *64*, 547-550.
50. Altshuler, B.L.; Levitov, L.S. Weak chaos in a quantum Kepler problem. *Phys. Rep.* **1997**, *288*, 487-512.

51. Levitov, L.S. Critical Hamiltonians with long range hopping. *Annalen der Physik* **1999**, *8*, 697-706.

52. Evers, F.; Mirlin A.D. Fluctuations of the inverse participation ratio at the Anderson transition. *Phys. Rev. Lett.* **2000**, *84*, 3690-3693.

53. Mirlin, A.D.; Evers, F. Multifractality and critical fluctuations at the Anderson transition. *Phys. Rev. B* **2000**, *62*, 7920-7933.

54. Fyodorov, Y.V.; Ossipov, A.; Rodriguez, A. The Anderson localization transition and eigenfunction multifractality in an ensemble of ultrametric random matrices. *J. Stat. Mech.* **2009**, L12001.

55. Fyodorov, Y.V.; Kupiainen, A.; Webb, C. Towards rigorous analysis of the Levitov-Mirlin-Evers recursion. *J. Phys. A* **2015**, arXiv:1509.01366.

56. Yektushenko, O.; Kratsov, V.E. Virial expansion for almost diagonal random matrices. *J. Phys. A* **2003**, *36*, 8265-8289.

57. Yektushenko, O.; Ossipov, A. Supersymmetry approach to almost diagonal random matrices. *J. Phys. A* **2007**, *40*, 4691-4716.

58. Kronmüller, S.; Yektushenko, O.M.; Cuevas, E. Supersymmetric virial expansion for time-reversal invariant disordered systems. *J. Phys. A* **2010**, *43*, 075001.

59. Kratsov, V.E.; Ossipov, A.; Yektushenko, O.M.; Cuevas, E. Dynamical scaling for critical states: Is Chalker's ansatz valid for strong fractality? *Phys. Rev. B* **2010**, *82*, 161102(R).

60. Bogomolny, E.; Giraud, O. Perturbation approach to multifractal dimensions for certain critical random-matrix ensembles. *Phys. Rev. E* **2012**, *84*, 036212.

61. Bogomolny, E.; Giraud, O. Multifractal dimensions for all moments for certain critical random matrix ensembles in the strong multifractality regime. *Phys. Rev. E* **2012**, *85*, 046208.

62. Bogomolny, E.; Giraud, O. Eigenfunction entropy and spectral compressibility for critical random matrix ensembles. *Phys. Rev. Lett.* **2011**, *106*, 044101.

63. Monthus, C.; Garel, T. Anderson localization transition with long-ranged hoppings: Analysis of the strong multifractality regime in terms of weighted Levy sums. *J. Stat. Mech.* **2010**, P09015.

64. Lévy, P. *Théorie de l'addition des Variables Aléatoires*; Gauthier-Villars: Paris, France, 1937. (In French)

65. Bouchaud, J.P.; Georges, A. Anomalous diffusion in disordered media: Statistical mechanisms, models and physical applications. *Phys. Rep.* **1990**, *195*, 127-293.

66. Derrida, B. From random walks to spin glasses. *Physica D* **1997**, *107*, 186-198.

67. Derrida, B.; Flyvbjerg, H. Statistical properties of randomly broken objects and of multivalley structures in disordered systems. *J. Phys. A* **1987**, *20*, 5273-5288.

68. Monthus, C.; Garel, T. Critical weight statistics of the random energy model and of the directed polymer on the Cayley tree. *Phys. Rev. E* **2007**, *75*, 051119.

69. Kunz, H.; Souillard, B. The localization transition on the Bethe lattice. *J. Phys. Lett.* **1983**, *44*, 411-414.

70. Mirlin, A.D.; Fyodorov, Y.V. Localization transition in the Anderson model on the Bethe lattice: Spontaneous symmetry breaking and correlation functions. *Nucl. Phys. B* **1991**, *366*, 507-532.

71. Monthus, C.; Garel, T. Anderson transition on the Cayley tree as a traveling wave critical point for various probability distributions. *J. Phys. A* **2009**, *42*, 075002.

72. Harris, A.B. Effect of random defects on the critical behaviour of Ising models. *J. Phys. C* **1974**, *7*, 1671-1692.

73. Chayes, J.T.; Chayes, L.; Fisher, D.S.; Spencer, T. Finite-size scaling and correlation lengths for disordered systems. *Phys. Rev. Lett.* **1986**, *57*, 2999-3002.

74. De Luca, A.; Scardicchio, A. Ergodicity breaking in a model showing many-body localization. *Europhys. Lett.* **2013**, *101*, 37003.

75. Torres-Herrera, E.J.; Santos, L.S. Dynamics at the many-body localization transition. *Phys. Rev. B* **2015**, *92*, 014208.

76. Chen, X.; Yu, X.; Cho, G.Y.; Clark, B.K. Eduardo Fradkin Many-body localization transition in Rokhsar-Kivelson-type wave functions. *Phys. Rev. B* **2015**, *92*, 214204.

77. Serbyn, M.; Papic, Z.; Abanin, D.A. Criterion for many-body localization-delocalization phase transition. *Phys. Rev. X* **2015**, *041047*.

78. Serbyn, M.; Moore, J.E. Spectral statistics across the many-body localization transition. *Phys. Rev. B* **2016**, *93*, 041424.

79. Monthus, C. Level repulsion exponent $\beta$ for many-body localization transitions and for Anderson localization transitions via Dyson Brownian motion. *J. Stat. Mech.* **2016**, *3*, 033113.
80. Pietracaprina, F.; Ros, V.; Scardicchio, A. Forward approximation as a mean-field approximation for the Anderson and many-body localization transitions. *Phys. Rev. B* **2016**, *93*, 054201.

81. Laumann, C.R.; Pal, A.; Scardicchio, A. Many-body mobility edge in a mean-field quantum spin glass. *Phys. Rev. Lett.* **2014**, *113*, 200405.

82. Baldwin, C.L.; Laumann, C.R.; Pal, A.; Scardicchio, A. The many-body localized phase of the quantum random energy model. *Phys. Rev. B* **2016**, *93*, 024202.

<figure><img src="image_2.png" /></figure>

© 2016 by the author; licensee MDPI, Basel, Switzerland. This article is an open access article distributed under the terms and conditions of the Creative Commons by Attribution (CC-BY) license (http://creativecommons.org/licenses/by/4.0/).