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

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 ν 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\left(3\right)$ obtained by Luitz et al. (2015).

1. 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 the references therein): the idea is that each many-body excited eigenstate $|\psi >$ is “thermal”, i.e., the reduced density matrix ${\rho }_A$ of a sub-region A corresponds the thermal density matrix:

$$\begin{array}{c}\hfill {\rho }_A^{\left(\beta \right)}=\frac{e^{-\beta H_A}}{T{r}_A\left(e^{-\beta H_A}\right)}\end{array}$$

(1)

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

$$\begin{array}{c}\hfill S_A\equiv -T{r}_A({\rho }_{A}ln{\rho }_A)\end{array}$$

(2)

coincides with the thermal entropy, which is extensive with respect to the volume $L_A^d$ of the region A:

$$\begin{array}{c}\hfill S_A\underset{ETH}{\simeq }-T{r}_A({\rho }_A^{th}ln{\rho }_A^{th})=s^{th}\left(\beta \right)L_A^d\end{array}$$

(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]:

$$\begin{array}{c}\hfill \overline{S_A}\underset{MBL}{\propto }{L}_A^{d-1}\end{array}$$

(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,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,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,21,22,23,24,25]. This construction is actually possible only because the MBL phase is characterized by an extensive number of emergent localized conserved operators [26,27,28,29,30,31,32,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:

$$\begin{array}{c}\hfill \overline{S_A}\underset{criti}{\simeq }{L}_A^{\alpha }\end{array}$$

(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 ξ:

$$\begin{array}{c}\hfill \overline{S_A}\underset{FSS}{\simeq }{L}_A^{\alpha }\mathsf{\Phi }\left(\frac{L_A}{\xi }\right)\end{array}$$

(6)

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

$$\begin{array}{c}\hfill \overline{S_A}\underset{critical,loc}{\propto }{\xi }^{\alpha }\end{array}$$

(7)

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

$$\begin{array}{c}\hfill \overline{S_A}\underset{critical,deloc}{\propto }{L}_A{\left(\frac{1}{\xi }\right)}^{1-\alpha }\end{array}$$

(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 ${\left(\frac{1}{\xi }\right)}^{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 α 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,41,42,43,44,45,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\left(1\right)$ 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+ϵ$ and the “strong multifractality” in high dimension d.

(ii)

For the tight-binding model with long-range hopping:

$$\begin{array}{c}\hfill V\left(r\right)=\frac{V}{r^a}\end{array}$$

(9)

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

$$\begin{array}{c}\hfill \mathsf{\Delta }E\left(L^d\right)\propto L^{-d}\end{array}$$

(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,49,50,51,52,53,54,55] or other methods [56,57,58,59,60,61,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

$$\begin{array}{c}\hfill H=H_0+V\end{array}$$

(11)

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

$$\begin{array}{c}\hfill H_0=-\sum _{i=1}^L{h}_i{\sigma }_i^z\end{array}$$

(12)

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

$$\begin{array}{c}\hfill G_W\left(h\right)=\frac{1}{\sqrt{2\pi W^2}}{e}^{-\frac{h^2}{2W^2}}\end{array}$$

(13)

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

$$\begin{array}{ccc}\hfill |{\psi }_{S_1,..,S_L}^{\left(0\right)}>& \equiv & |S_1>\otimes |S_2>...\otimes |S_L>\hfill \end{array}$$

(14)

with the random energies:

$$\begin{array}{c}\hfill E_{S_1,..,S_L}^{\left(0\right)}=-\sum _{i=1}^L{h}_i{S}_i\end{array}$$

(15)

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

$$\begin{array}{c}\hfill E_{GS}^{\left(0\right)}=-\sum _{i=1}^L\left|h_i\right|\end{array}$$

(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 $\left(L{W}^2\right)$:

$$\begin{array}{c}\hfill {\rho }_0\left(E\right)=\frac{1}{\sqrt{2\pi L{W}^2}}{e}^{-\frac{E^2}{2L{W}^2}}\end{array}$$

(17)

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

$$\begin{array}{c}\hfill \mathsf{\Delta }E\left(L\right)\propto L^{\frac{1}{2}}{2}^{-L}\end{array}$$

(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:

$$\begin{array}{c}\hfill V=-\sum _{k=1}^L\sum _{1\le i_1<i_2..<i_k\le L}{J}_{i_1,..,i_k}{\sigma }_{i_1}^x{\sigma }_{i_2}^x...{\sigma }_{i_k}^x\end{array}$$

(19)

The couplings $J_{i_1,...,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,...,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:

$$\begin{array}{c}\hfill E_{S_1,...,S_L}^{\left(1\right)}=E_{S_1,...,S_L}^{\left(0\right)}+<{\psi }_{S_1,...,S_L}^{\left(0\right)}\left|V\right|{\psi }_{S_1,...,S_L}^{\left(0\right)}>=E_{S_1,...,S_L}^{\left(0\right)}\end{array}$$

(20)

and the eigenstates read:

$$\begin{array}{ccc}\hfill |{\psi }_{S_1,...,S_L}^{\left(1\right)}>& =& |{\psi }_{S_1,...,S_L}^{\left(0\right)}>+\sum _{\left\{S_i^{\prime }\right\}}|{\psi }_{S_1^{\prime },...,S_L^{\prime }}^{\left(0\right)}>\frac{<{\psi }_{S_1^{\prime },...,S_L^{\prime }}^{\left(0\right)}\left|V\right|{\psi }_{S_1,...,S_L}^{\left(0\right)}>}{E_{S_1,...,S_L}^{\left(0\right)}-E_{S_1^{\prime },...,S_L^{\prime }}^{\left(0\right)}}\hfill \end{array}$$

(21)

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

$$\begin{array}{ccc}\hfill |0>& \equiv & |S_1=+1,S_2=+1,...,S_L=+1>\hfill \end{array}$$

(22)

whose energy:

$$\begin{array}{c}\hfill E_0^{\left(0\right)}=-\sum _{i=1}^L{h}_i\end{array}$$

(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,...,i_k)$ of flipped spins with respect to this reference configuration:

$$\begin{array}{ccc}\hfill |i_1,...,i_k>& \equiv & {\sigma }_{i_1}^x{\sigma }_{i_2}^x...{\sigma }_{i_k}^x|0>\hfill \end{array}$$

(24)

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

$$\begin{array}{ccc}\hfill |{\psi }^{\left(1\right)}>& =& |0>+\sum _{k=1}^L\sum _{1\le i_1<i_2...<i_k\le L}\left(\frac{J_{i_1,i_2,...,i_k}}{{\displaystyle 2\sum _{q=1}^k{h}_{i_q}}}\right)|i_1,...,i_k>\hfill \end{array}$$

(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:

$$\begin{array}{ccc}\hfill |0{>}_A& \equiv & |S_1=+1,S_2=+1,...,S_{L_A}=+1>\hfill \\ \hfill |i_1,...,i_k{>}_A& \equiv & {\sigma }_{i_1}^x{\sigma }_{i_2}^x...{\sigma }_{i_k}^x|0{>}_A\mathrm{for}1\le i_1<i_2...<i_k\le L_A\hfill \end{array}$$

(26)

and:

$$\begin{array}{ccc}\hfill |0{>}_B& \equiv & |S_{L_A+1}=+1,S_{L+2}=+1,...,S_L=+1>\hfill \\ \hfill |j_1,...,j_k{>}_B& \equiv & {\sigma }_{j_1}^x{\sigma }_{j_2}^x...{\sigma }_{j_k}^x|0{>}_B\mathrm{for}{L}_A+1\le j_1<j_2...<j_k\le L_A+L_B\hfill \end{array}$$

(27)

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

$$\begin{array}{ccc}& & |{\psi }^{\left(1\right)}>=|0{>}_A\otimes |0{>}_B\hfill \\ & & +\sum _{k=1}^{L_A}\sum _{1\le i_1<i_2...<i_k\le L_A}\frac{J_{i_1,i_2,...,i_k}}{2{\sum }_{q=1}^k{h}_{i_q}}|i_1,...,i_k{>}_A\otimes |0{>}_B\hfill \\ & & +\sum _{k=1}^{L_B}\sum _{L_A+1\le j_1<j_2...<k_k\le L_A+L_B}\frac{J_{i_1,i_2,...,i_k}}{2{\sum }_{q=1}^k{h}_{i_q}}|0{>}_A\otimes |i_1,...,i_k{>}_B\hfill \\ & & +\sum _{k_A=1}^{L_A}\sum _{1\le i_1<i_2...<i_{k_A}\le L_A}\sum _{k_B=1}^{L_B}\sum _{L_A+1\le j_1<j_2...<j_{k_B}\le L_A+L_B}\frac{J_{i_1,...,i_{k_A},j_1,...,j_{k_B}}}{2{\sum }_{q=1}^{k_A}{h}_{i_q}+2{\sum }_{p=1}^{k_B}{h}_{j_p}}|i_1,...,i_{k_A}{>}_A\otimes |j_1,...,j_{k_B}{>}_B\hfill \end{array}$$

(28)

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

$$\begin{array}{ccc}\hfill |{\varphi }_0^A>& \equiv & |0{>}_A+\sum _{k=1}^{L_A}\sum _{1\le i_1<i_2...<i_k\le L_A}\frac{J_{i_1,i_2,...,i_k}}{2{\sum }_{q=1}^k{h}_{i_q}}|i_1,...,i_k{>}_A\hfill \\ \hfill |{\varphi }_{i_1,...,i_k}^A>& \equiv & |i_1,...,i_k{>}_A-\frac{J_{i_1,i_2,...,i_k}}{2{\sum }_{q=1}^k{h}_{i_q}}|0{>}_A\mathrm{for}1\le i_1<i_2...<i_k\le L_A\hfill \end{array}$$

(29)

and the similar basis for the region B:

$$\begin{array}{ccc}\hfill |{\varphi }_0^B>& \equiv & |0{>}_B+\sum _{k=1}^{L_B}\sum _{L_A+1\le j_1<j_2...<j_k\le L_A+L_B}\frac{J_{j_1,j_2,...,j_k}}{2{\sum }_{q=1}^k{h}_{j_q}}|j_1,...,j_k{>}_B\hfill \\ \hfill |{\varphi }_{j_1,...,j_k}^B>& \equiv & |j_1,...,j_k{>}_B-\frac{J_{j_1,j_2,...,j_k}}{2{\sum }_{q=1}^k{h}_{j_q}}|0{>}_B\mathrm{for}{L}_A+1\le j_1<j_2...<j_k\le L_A+L_B\hfill \end{array}$$

(30)

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

$$\begin{array}{ccc}& & |{\psi }^{\left(1\right)}>=|{\varphi }_0^A>\otimes |{\varphi }_0^B>\hfill \\ & & +\sum _{k_A=1}^{L_A}\sum _{1\le i_1<i_2...<i_{k_A}\le L_A}\sum _{k_B=1}^{L_B}\sum _{L_A+1\le j_1<j_2...<j_{k_B}\le L_A+L_B}\frac{J_{i_1,...,i_{k_A},j_1,...j_{k_B}}}{2{\sum }_{q=1}^{k_A}{h}_{i_q}+2{\sum }_{p=1}^{k_B}{h}_{j_p}}|{\varphi }_{i_1,...,i_{k_A}}^A>\otimes |{\varphi }_{j_1,...,j_{k_B}}^B>\hfill \end{array}$$

(31)

The corresponding density matrix:

$$\begin{array}{ccc}\hfill \rho & \equiv & \frac{|{\psi }^{\left(1\right)}><{\psi }^{\left(1\right)}|}{<{\psi }^{\left(1\right)}|{\psi }^{\left(1\right)}>}=\frac{1}{1+{\mathsf{\Sigma }}_1}\left(\right|{\varphi }_0^A>\otimes |{\varphi }_0^B>\hfill \\ & & +\sum _{k_A=1}^{L_A}\sum _{1\le i_1<i_2...<i_{k_A}\le L_A}\sum _{k_B=1}^{L_B}\sum _{L_A+1\le j_1<j_2...<j_{k_B}\le L_A+L_B}\frac{J_{i_1,...,i_{k_A},j_1,...,j_{k_B}}}{2{\sum }_{q=1}^{k_A}{h}_{i_q}+2{\sum }_{p=1}^{k_B}{h}_{j_p}}|{\varphi }_{i_1,...,i_{k_A}}^A>\otimes |{\varphi }_{j_1,...,j_{k_B}}^B>)\hfill \\ & & (<{\varphi }_0^A|\otimes <{\varphi }_0^B|\hfill \\ & & +\sum _{k_A^{\prime }=1}^{L_A}\sum _{1\le i_1^{\prime }<i_2^{\prime }..<i_{k_A^{\prime }}^{\prime }\le L_A}\sum _{k_B^{\prime }=1}^{L_B}\sum _{L_A+1\le j_1^{\prime }<j_2^{\prime }..<j_{k_B^{\prime }}^{\prime }\le L_A+L_B}\frac{J_{i_1^{\prime },..,i_{k_A}^{\prime },j_1^{\prime },...j_{k_B^{\prime }}^{\prime }}}{2{\sum }_{q=1}^{k_A^{\prime }}{h}_{i_q^{\prime }}+2{\sum }_{p=1}^{k_B^{\prime }}{h}_{j_p^{\prime }}}<{\varphi }_{i_1^{\prime },..,i_{k_A^{\prime }}^{\prime }}^A|\otimes <{\varphi }_{j_1^{\prime },..,j_{k_B^{\prime }}^{\prime }}^B\left|\right)\hfill \end{array}$$

(32)

with:

$$\begin{array}{ccc}\hfill {\mathsf{\Sigma }}_1& \equiv & \sum _{k_A=1}^{L_A}\sum _{1\le i_1<i_2...<i_{k_A}\le L_A}\sum _{k_B=1}^{L_B}\sum _{L_A+1\le j_1<j_2...<j_{k_B}\le L_A+L_B}{\left(\frac{J_{i_1,...,i_{k_A},j_1,...j_{k_B}}}{2{\sum }_{q=1}^{k_A}{h}_{i_q}+2{\sum }_{p=1}^{k_B}{h}_{j_p}}\right)}^2\hfill \end{array}$$

(33)

leads to the reduced density matrix for the region A:

$$\begin{array}{ccc}\hfill {\rho }_A& \equiv & T{r}_B\left(\rho \right)=<{\varphi }_0^B\left|\rho \right|{\varphi }_0^B>+\sum _{k_B=1}^{L_B}\sum _{L_A+1\le j_1<j_2...<j_{k_B}\le L_A+L_B}<{\varphi }_{j_1,...,j_{k_B}}^B\left|\rho \right|{\varphi }_{j_1,...,j_{k_B}}^B>\hfill \\ & =& \frac{1}{1+{\mathsf{\Sigma }}_1}|{\varphi }_0^A><{\varphi }_0^A|\hfill \\ & & +\sum _{k_A=1}^{L_A}\sum _{1\le i_1<i_2...<i_{k_A}\le L_A}\sum _{k_A^{\prime }=1}^{L_A}\sum _{1\le i_1^{\prime }<i_2^{\prime }...<i_{k_A}^{\prime }\le L_A}\frac{R(i_1,...,i_{k_A};i_1^{\prime },...,i_{k_A^{\prime }}^{\prime })}{1+{\mathsf{\Sigma }}_1}|{\varphi }_{i_1,...,i_{k_A}}^A><{\varphi }_{i_1^{\prime },...,i_{k_A^{\prime }}^{\prime }}^A|\hfill \end{array}$$

(34)

with:

$$\begin{array}{c}\hfill R(i_1,...,i_{k_A};i_1^{\prime },...,i_{k_A^{\prime }}^{\prime })\equiv \sum _{k_B=1}^{L_B}\sum _{L_A+1\le j_1<j_2...<j_{k_B}\le L_A+L_B}\left(\frac{J_{i_1,...,i_{k_A},j_1,...,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^{\prime },...,i_{k_A}^{\prime },j_1,...,j_{k_B}}}{2{\sum }_{q=1}^{k_A^{\prime }}{h}_{i_q^{\prime }}+2{\sum }_{p=1}^{k_B}{h}_{j_p}}\right)\end{array}$$

(35)

In particular, the diagonal elements involve the positive coefficients:

$$\begin{array}{c}\hfill D_{i_1,...,i_{k_A}}\equiv R(i_1,...,i_{k_A};i_1,...,i_{k_A})\equiv \sum _{k_B=1}^{L_B}\sum _{L_A+1\le j_1<j_2...<j_{k_B}\le L_A+L_B}{\left(\frac{J_{i_1,...,i_{k_A},j_1,...,j_{k_B}}}{2{\sum }_{q=1}^{k_A}{h}_{i_q}+2{\sum }_{p=1}^{k_B}{h}_{j_p}}\right)}^2\end{array}$$

(36)

and their sum corresponds to Equation (33):

$$\begin{array}{c}\hfill {\mathsf{\Sigma }}_1\equiv \sum _{k_A=1}^{L_A}\sum _{1\le i_1<i_2...<i_{k_A}\le L_A}{D}_{i_1,...,i_{k_A}}\end{array}$$

(37)

2.5. Entanglement Spectrum

As we will see in the next section, the diagonal coefficients $D_{i_1,...,i_{k_A}}$ of Equation (36) and their sum ${\mathsf{\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/{\left(2{\sum }_{q=1}^{k_A}{h}_{i_q}+2{\sum }_{p=1}^{k_B}{h}_{j_p}\right)}^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,...,i_{k_A}}$ of Equation (36) and their sum ${\mathsf{\Sigma }}_1$ give contributions of first order $O\left(\right|J\left|\right)$ 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\left(J^2\right)$ (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{array}{ccc}\hfill p_0& =& \frac{1}{1+{\mathsf{\Sigma }}_1}\hfill \\ \hfill p_{i_1,...,i_{k_A}}& =& \frac{D_{i_1,...,i_{k_A}}}{1+{\mathsf{\Sigma }}_1}\mathrm{for}1\le i_1<i_2...<i_{k_A}\le L_A\hfill \end{array}$$

(38)

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

$$\begin{array}{c}\hfill Y_q\left(L_A\right)=T{r}_A\left({\rho }_A^q\right)=p_0^q+\sum _{k_A=1}^{L_A}\sum _{1\le i_1<i_2...<i_{k_A}\le L_A}{p}_{i_1,...,i_{k_A}}^q=\frac{1+{\mathsf{\Sigma }}_q}{{(1+{\mathsf{\Sigma }}_1)}^q}\end{array}$$

(39)

where the sum:

$$\begin{array}{c}\hfill {\mathsf{\Sigma }}_q\equiv \sum _{k_A=1}^{L_A}\sum _{1\le i_1<i_2...<i_{k_A}\le L_A}{D}_{i_1,...,i_{k_A}}^q\end{array}$$

(40)

generalizes Equation (37).

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

$$\begin{array}{c}\hfill S_q\left(L_A\right)\equiv \frac{ln{Y}_q\left(L_A\right)}{1-q}=\frac{ln(1+{\mathsf{\Sigma }}_q)-qln(1+{\mathsf{\Sigma }}_1)}{1-q}\end{array}$$

(41)

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

$$\begin{array}{c}\hfill S_1\left(L_A\right)\equiv -T{r}_A({\rho }_{A}ln{\rho }_A)=-p_{0}ln{p}_0-\sum _{k_A=1}^L\sum _{1\le i_1<i_2...<i_{k_A}\le L}{p}_{i_1,...,i_{k_A}}ln{p}_{i_1,...,i_{k_A}}\end{array}$$

(42)

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

$$\begin{array}{c}\hfill \overline{S_q}=\frac{\overline{ln(1+{\mathsf{\Sigma }}_q)}-q\overline{ln(1+{\mathsf{\Sigma }}_1)}}{1-q}\end{array}$$

(43)

one only needs to study separately the probability distributions of ${\mathsf{\Sigma }}_q$ and ${\mathsf{\Sigma }}_1$, as done in the next section. Other statistical properties involving correlations between ${\mathsf{\Sigma }}_q$ and ${\mathsf{\Sigma }}_1$ can also be computed (see, for instance, Appendix B for the disorder-averaged values $\overline{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,...,i_{k_A}}$

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

$$\begin{array}{ccc}\hfill \overline{e^{-t{D}_{i_1,...,i_{k_A}}}}& =& {\int }_0^{+\infty }d{D}_{i_1,...,i_{k_A}}{P}_{i_1,...,i_{k_A}}\left(D_{i_1,...,i_{k_A}}\right)e^{-t{D}_{i_1,...,i_{k_A}}}\hfill \\ & =& \overline{e^{-t{\sum }_{k_B=1}^{L_B}{\sum }_{L_A+1\le j_1<j_2...<j_{k_B}\le L_A+L_B}{\left(\frac{J_{i_1,...,i_{k_A},j_1,...,j_{k_B}}}{2{\sum }_{q=1}^{k_A}{h}_{i_q}+2{\sum }_{p=1}^{k_B}{h}_{j_p}}\right)}^2}}\hfill \\ & =& \overline{\prod _{k_B=1}^{L_B}\prod _{L_A+1\le j_1<j_2...<j_{k_B}\le L_A+L_B}\left[1-\left(1-e^{-t{\left(\frac{J_{i_1,...,i_{k_A},j_1,...,j_{k_B}}}{2{\sum }_{q=1}^{k_A}{h}_{i_q}+2{\sum }_{p=1}^{k_B}{h}_{j_p}}\right)}^2}\right)\right]}\hfill \\ & =& 1-\sum _{k_B=1}^{L_B}\sum _{L_A+1\le j_1<j_2...<j_{k_B}\le L_A+L_B}\overline{\left(1-e^{-t{\left(\frac{J_{i_1,...,i_{k_A},j_1,...,j_{k_B}}}{2{\sum }_{q=1}^{k_A}{h}_{i_q}+2{\sum }_{p=1}^{k_B}{h}_{j_p}}\right)}^2}\right)}+o\left(J\right)\hfill \end{array}$$

(44)

The variable:

$$\begin{array}{c}\hfill E=2\sum _{q=1}^{k_A}{h}_{i_q}+2\sum _{p=1}^{k_B}{h}_{j_p}\end{array}$$

(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 $\left(4W^2(k_A+k_B)\right)$:

$$\begin{array}{c}\hfill G_{4W^2(k_A+k_B)}\left(E\right)=\frac{1}{2\sqrt{2(k_A+k_B)\pi W^2}}{e}^{-\frac{E^2}{8(k_A+k_B)W^2}}\end{array}$$

(46)

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

$$\begin{array}{ccc}& & \overline{1-e^{-t{\left(\frac{J_{i_1,...,i_{k_A},j_1,...,j_{k_B}}}{2{\sum }_{q=1}^{k_A}{h}_{i_q}+2{\sum }_{p=1}^{k_B}{h}_{j_p}}\right)}^2}}\hfill \\ & & =\overline{{\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,...,i_{k_A},j_1,...,j_{k_B}}}{E}\right)}^2}\right)}\hfill \\ & & =\overline{|J_{i_1,...,i_{k_A},j_1,...,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,i_2,...,i_k}^2{t}_A}{8(k_A+k_B)W^{2}x}}\left(1-e^{-x}\right)}\hfill \\ & & =\overline{|J_{i_1,...,i_{k_A},j_1,...,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}}}\left(1-e^{-x}\right)+o\left(J_{i_1,...,i_{k_A},j_1,...,j_{k_B}}\right)\hfill \\ & & =\overline{|J_{i_1,...,i_{k_A},j_1,...,j_{k_B}}|}{t}^{\frac{1}{2}}\frac{1}{\sqrt{2(k_A+k_B)W^2}}+o\left(J_{i_1,...,i_{k_A},j_1,...,j_{k_B}}\right)\hfill \end{array}$$

(47)

Therefore, Equation (44) becomes at first order in the couplings $J_{i_1,...,i_{k_A},j_1,...,j_{k_B}}$:

$$\begin{array}{ccc}& & \overline{e^{-t_{i_1,...,i_{k_A}}{D}_{i_1,...,i_{k_A}}}}\hfill \\ & & =1-t^{\frac{1}{2}}\sum _{k_B=1}^{L_B}\sum _{L_A+1\le j_1<j_2...<j_{k_B}\le L_A+L_B}\overline{|J_{i_1,...,i_{k_A},j_1,...,j_{k_B}}|}\frac{1}{\sqrt{2(k_A+k_B)W^2}}+o\left(J\right)\hfill \\ & & =e^{-t^{\frac{1}{2}}{\sum }_{k_B=1}^{L_B}{\sum }_{L_A+1\le j_1<j_2...<j_{k_B}\le L_A+L_B}\overline{|J_{i_1,...,i_{k_A},j_1,...,j_{k_B}}|}\frac{1}{\sqrt{2(k_A+k_B)W^2}}+o\left(J\right)}\hfill \end{array}$$

(48)

It is thus convenient to introduce the notation:

$$\begin{array}{c}\hfill {\mathsf{\Omega }}_{i_1,...,i_{k_A}}\equiv \sum _{k_B=1}^{L_B}\sum _{L_A+1\le j_1<j_2...<j_{k_B}\le L_A+L_B}\overline{|J_{i_1,...,i_{k_A},j_1,...,j_{k_B}}|}\frac{1}{\sqrt{2(k_A+k_B)W^2}}\end{array}$$

(49)

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

$$\begin{array}{c}\hfill L_{\frac{1}{2};\mathsf{\Omega }}\left(D\right)\equiv \frac{\mathsf{\Omega }}{2\sqrt{\pi }{D}^{\frac{3}{2}}}{e}^{-\frac{{\mathsf{\Omega }}^2}{4D}}\end{array}$$

(50)

with its Laplace transform:

$$\begin{array}{c}\hfill {\widehat{L}}_{\frac{1}{2},\mathsf{\Omega }}\left(t\right)\equiv {\int }_0^{+\infty }dD{L}_{\frac{1}{2};\mathsf{\Omega }}\left(D\right)e^{-tD}=e^{-\mathsf{\Omega }{t}^{\frac{1}{2}}}\end{array}$$

(51)

The parameter Ω 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,66,67,68] for more details on their properties).

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

$$\begin{array}{c}\hfill P_{i_1,...,i_{k_A}}\left(D_{i_1,...,i_{k_A}}\right)=L_{\frac{1}{2};{\mathsf{\Omega }}_{i_1,...,i_{k_A}}}\left(D_{i_1,...,i_{k_A}}\right)\end{array}$$

(52)

3.2. Probability Distribution of ${\mathsf{\Sigma }}_1$

The probability distribution $P_1\left({\mathsf{\Sigma }}_1\right)$ of the positive random variable $\mathsf{\Sigma }1$ of Equation (37) can be evaluated similarly via its Laplace transform:

$$\begin{array}{ccc}\hfill \overline{e^{-t{\mathsf{\Sigma }}_1}}& \equiv & {\int }_0^{+\infty }d{\mathsf{\Sigma }}_1{P}_1\left({\mathsf{\Sigma }}_1\right)e^{-t{\mathsf{\Sigma }}_1}\hfill \\ & =& \overline{e^{-t{\sum }_{k_A=1}^{L_A}{\sum }_{1\le i_1<i_2...<i_{k_A}\le L_A}{\sum }_{k_B=1}^{L_B}{\sum }_{L_A+1\le j_1<j_2...<j_{k_B}\le L_A+L_B}{\left(\frac{J_{i_1,...,i_{k_A},j_1,...,j_{k_B}}}{2{\sum }_{q=1}^{k_A}{h}_{i_q}+2{\sum }_{p=1}^{k_B}{h}_{j_p}}\right)}^2}}\hfill \\ & =& \overline{\prod _{k_A=1}^{L_A}\prod _{1\le i_1<i_2...<i_{k_A}\le L_A}\prod _{k_B=1}^{L_B}\prod _{L_A+1\le j_1<j_2...<j_{k_B}\le L_A+L_B}\left[1-\left(1-e^{-t{\left(\frac{J_{i_1,...,i_{k_A},j_1,...,j_{k_B}}}{2{\sum }_{q=1}^{k_A}{h}_{i_q}+2{\sum }_{p=1}^{k_B}{h}_{j_p}}\right)}^2}\right)\right]}\hfill \\ & =& 1-\sum _{k_A=1}^{L_A}\sum _{1\le i_1<i_2...<i_{k_A}\le L_A}\sum _{k_B=1}^{L_B}\sum _{L_A+1\le j_1<j_2...<j_{k_B}\le L_A+L_B}\overline{\left(1-e^{-t{\left(\frac{J_{i_1,...,i_{k_A},j_1,...,j_{k_B}}}{2{\sum }_{q=1}^{k_A}{h}_{i_q}+2{\sum }_{p=1}^{k_B}{h}_{j_p}}\right)}^2}\right)}+o\left(J\right)\hfill \end{array}$$

(53)

Equation (47) then yields that:

$$\begin{array}{ccc}\hfill \overline{e^{-t{\mathsf{\Sigma }}_1}}& =& 1-t^{\frac{1}{2}}\sum _{k_A=1}^{L_A}\sum _{1\le i_1<i_2...<i_{k_A}\le L_A}\sum _{k_B=1}^{L_B}\sum _{L_A+1\le j_1<j_2...<j_{k_B}\le L_A+L_B}\overline{|J_{i_1,...,i_{k_A},j_1,...,j_{k_B}}|}\frac{1}{\sqrt{2(k_A+k_B)W^2}}+o\left(J\right)\hfill \\ & =& e^{-{\mathsf{\Omega }}_1{t}^{\frac{1}{2}}+o\left(J\right)}\hfill \end{array}$$

(54)

with:

$$\begin{array}{ccc}\hfill {\mathsf{\Omega }}_1& \equiv & \sum _{k_A=1}^{L_A}\sum _{1\le i_1<i_2...<i_{k_A}\le L_A}{\mathsf{\Omega }}_{i_1,...,i_{k_A}}\hfill \\ & =& \sum _{k_A=1}^{L_A}\sum _{1\le i_1<i_2...<i_{k_A}\le L_A}\sum _{k_B=1}^{L_B}\sum _{L_A+1\le j_1<j_2...<j_{k_B}\le L_A+L_B}\left|J_{i_1,...,i_{k_A},j_1,...,j_{k_B}}\right|\frac{1}{\sqrt{2(k_A+k_B)W^2}}\hfill \end{array}$$

(55)

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

$$\begin{array}{c}\hfill P_1\left({\mathsf{\Sigma }}_1\right)=L_{\frac{1}{2};{\mathsf{\Omega }}_1}\left({\mathsf{\Sigma }}_1\right)\end{array}$$

(56)

In particular, using the identity:

$$\begin{array}{c}\hfill ln(1+{\mathsf{\Sigma }}_1)={\int }_0^{+\infty }\frac{dt}{t}{e}^{-t}\left(1-e^{-t{\mathsf{\Sigma }}_1}\right)\end{array}$$

(57)

one obtains:

$$\begin{array}{ccc}\hfill \overline{ln(1+{\mathsf{\Sigma }}_1)}& =& {\int }_0^{+\infty }d{\mathsf{\Sigma }}_1{P}_1\left({\mathsf{\Sigma }}_1\right)ln(1+{\mathsf{\Sigma }}_1)\hfill \\ & =& {\int }_0^{+\infty }\frac{dt}{t}{e}^{-t}\left(1-\overline{e^{-t{\mathsf{\Sigma }}_1}}\right)\hfill \\ & =& {\int }_0^{+\infty }\frac{dt}{t}{e}^{-t}{\mathsf{\Omega }}_1{t}^{\frac{1}{2}}+o\left(J\right)\hfill \\ & =& {\mathsf{\Omega }}_1\sqrt{\pi }+o\left(J\right)\hfill \end{array}$$

(58)

3.3. Probability Distribution of the Weight $p_0$

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

$$\begin{array}{c}\hfill p_0=\frac{1}{1+{\mathsf{\Sigma }}_1}\end{array}$$

(59)

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

$$\begin{array}{c}\hfill \mathcal{P}\left(p_0\right)=\frac{1}{p_0^2}{L}_{\frac{1}{2};{\mathsf{\Omega }}_1}\left(\frac{1}{p_0}-1\right)=\frac{{\mathsf{\Omega }}_1{p}_0^{\frac{1}{2}}}{2\sqrt{\pi }{(1-p_0)}^{\frac{3}{2}}}{e}^{-\frac{{\mathsf{\Omega }}_1^2{p}_0}{4(1-p_0)}}\end{array}$$

(60)

In particular, its averaged value:

$$\begin{array}{c}\hfill {\int }_0^{1}d{p}_0{p}_0\mathcal{P}\left(p_0\right)=1-{\mathsf{\Omega }}_1{e}^{\frac{{\mathsf{\Omega }}_1^2}{4}}{\int }_{\frac{{\mathsf{\Omega }}_1}{2}}^{+\infty }dt{e}^{-t^2}\end{array}$$

(61)

is close to unity for small ${\mathsf{\Omega }}_1$:

$$\begin{array}{c}\hfill {\int }_0^{1}d{p}_0{p}_0\mathcal{P}\left(p_0\right)=1-\frac{\sqrt{\pi }}{2}{\mathsf{\Omega }}_1+O\left({\mathsf{\Omega }}_1^2\right)\end{array}$$

(62)

3.4. Probability Distribution of ${\mathsf{\Sigma }}_q$

Let us now consider the variable ${\mathsf{\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,...,i_{k_A}}$ distributed with the Lévy law of Equation (50) is finite and reads:

$$\begin{array}{c}\hfill \overline{D_{i_1,...,i_{k_A}}^q}={\int }_0^{+\infty }dD\frac{{\mathsf{\Omega }}_{i_1,...,i_{k_A}}}{2\sqrt{\pi }}{D}^{q-\frac{3}{2}}{e}^{-\frac{{\mathsf{\Omega }}_{i_1,...,i_{k_A}}^2}{4D}}=\frac{{\mathsf{\Omega }}_{i_1,...,i_{k_A}}^{2q}\mathsf{\Gamma }\left(\frac{1}{2}-q\right)}{4^q\sqrt{\pi }}\end{array}$$

(63)

As a consequence, the average of ${\mathsf{\Sigma }}_q$ is also finite and given by:

$$\begin{array}{c}\hfill \overline{{\mathsf{\Sigma }}_q}=\sum _{k_A=1}^L\sum _{1\le i_1<i_2...<i_{k_A}\le L}\overline{D_{i_1,...,i_{k_A}}^q}=\frac{\mathsf{\Gamma }\left(\frac{1}{2}-q\right)}{4^q\sqrt{\pi }}\sum _{k_A=1}^L\sum _{1\le i_1<i_2...<i_{k_A}\le L}{\mathsf{\Omega }}_{i_1,...,i_{k_A}}^{2q}\end{array}$$

(64)

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

$$\begin{array}{ccc}\hfill \overline{ln(1+{\mathsf{\Sigma }}_q)}& \simeq & \overline{{\mathsf{\Sigma }}_q}=\frac{\mathsf{\Gamma }\left(\frac{1}{2}-q\right)}{4^q\sqrt{\pi }}\sum _{k_A=1}^L\sum _{1\le i_1<i_2...<i_{k_A}\le L}{\mathsf{\Omega }}_{i_1,...,i_{k_A}}^{2q}+{o\left(\right|J|}^{2q})\hfill \end{array}$$

(65)

(ii)

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

$$\begin{array}{ccc}\hfill \overline{e^{-t{\mathsf{\Sigma }}_q}}& \equiv & {\int }_0^{+\infty }d{\mathsf{\Sigma }}_q{P}_q\left({\mathsf{\Sigma }}_q\right)e^{-t{\mathsf{\Sigma }}_q}=\prod _{k_A=1}^L\prod _{1\le i_1<i_2...<i_{k_A}\le L}\left[\overline{e^{-t{D}_{i_1,...,i_{k_A}}^q}}\right]\hfill \\ & =& \prod _{k_A=1}^L\prod _{1\le i_1<i_2...<i_{k_A}\le L}\overline{1-\left(1-e^{-t{D}_{i_1,...,i_{k_A}}^q}\right)}\hfill \\ & =& 1-\sum _{k_A=1}^L\sum _{1\le i_1<i_2...<i_{k_A}\le L}\overline{\left(1-e^{-t{D}_{i_1,...,i_{k_A}}^q}\right)}+o\left(J\right)\hfill \end{array}$$

(66)

Using Equation (50) and the change of variable $x=t{D}^q$, one obtains:

$$\begin{array}{ccc}\hfill \overline{\left(1-e^{-t{D}^q}\right)}& =& {\int }_0^{+\infty }dD\frac{\mathsf{\Omega }}{2\sqrt{\pi }{D}^{\frac{3}{2}}}{e}^{-\frac{{\mathsf{\Omega }}^2}{4D}}\left(1-e^{-t{D}^q}\right)\hfill \\ & =& \frac{\mathsf{\Omega }}{2q\sqrt{\pi }}{t}^{\frac{1}{2q}}{\int }_0^{+\infty }\frac{dx}{x^{1+\frac{1}{2q}}}(1-e^{-x})+o\left(J\right)\hfill \\ & =& \frac{\mathsf{\Omega }\sqrt{\pi }}{sin\left(\frac{\pi }{2q}\right)\mathsf{\Gamma }\left(\frac{1}{2q}\right)}{t}^{\frac{1}{2q}}+o\left(J\right)\hfill \end{array}$$

(67)

Therefore, Equation (66) becomes:

$$\begin{array}{ccc}\hfill \overline{e^{-t{\mathsf{\Sigma }}_q}}& =& 1-\frac{\sqrt{\pi }}{sin\left(\frac{\pi }{2q}\right)\mathsf{\Gamma }\left(\frac{1}{2q}\right)}{t}^{\frac{1}{2q}}\sum _{k_A=1}^L\sum _{1\le i_1<i_2...<i_{k_A}\le L}{\mathsf{\Omega }}_{i_1,...,i_{k_A}}+o\left(J\right)\hfill \\ & =& e^{-\frac{\sqrt{\pi }{\mathsf{\Omega }}_1}{sin\left(\frac{\pi }{2q}\right)\mathsf{\Gamma }\left(\frac{1}{2q}\right)}{t}^{\frac{1}{2q}}}+o\left(J\right)\hfill \end{array}$$

(68)

It is thus convenient to introduce the notation:

$$\begin{array}{c}\hfill {\mathsf{\Omega }}_q\equiv \frac{\sqrt{\pi }{\mathsf{\Omega }}_1}{sin\left(\frac{\pi }{2q}\right)\mathsf{\Gamma }\left(\frac{1}{2q}\right)}\end{array}$$

(69)

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

$$\begin{array}{c}\hfill {\widehat{L}}_{\mu ,\mathsf{\Omega }}\left(t\right)\equiv {\int }_0^{+\infty }dD{L}_{\mu ;\mathsf{\Omega }}\left(D\right)e^{-tD}=e^{-\mathsf{\Omega }{t}^{\mu }}\end{array}$$

(70)

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

$$\begin{array}{c}\hfill P_q\left({\mathsf{\Sigma }}_q\right)=L_{\frac{1}{2q};{\mathsf{\Omega }}_q}\left({\mathsf{\Sigma }}_q\right)\end{array}$$

(71)

That generalizes Equation (56). In particular, using the identity of Equation (58) for ${\mathsf{\Sigma }}_q$, one obtains:

$$\begin{array}{ccc}\hfill \overline{ln(1+{\mathsf{\Sigma }}_q)}& =& {\int }_0^{+\infty }d{\mathsf{\Sigma }}_q{P}_q\left({\mathsf{\Sigma }}_q\right)ln(1+{\mathsf{\Sigma }}_q)\hfill \\ & =& {\int }_0^{+\infty }\frac{dt}{t}{e}^{-t}\left(1-\overline{e^{-t{\mathsf{\Sigma }}_q}}\right)\hfill \\ & =& {\int }_0^{+\infty }\frac{dt}{t}{e}^{-t}{\mathsf{\Omega }}_q{t}^{\frac{1}{2q}}+o\left(J\right)\hfill \\ & =& {\mathsf{\Omega }}_q\mathsf{\Gamma }\left(\frac{1}{2q}\right)+o\left(J\right)\hfill \\ & =& {\mathsf{\Omega }}_1\frac{\sqrt{\pi }}{sin\left(\frac{\pi }{2q}\right)}+o\left(J\right)\hfill \end{array}$$

(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 ${\left|J\right|}^{2q}$ dominates over the contribution of Equation (58), so that the leading order reads:

$$\begin{array}{c}\hfill \overline{S_q}=\frac{\mathsf{\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...<i_{k_A}\le L}{\mathsf{\Omega }}_{i_1,...,i_{k_A}}^{2q}+{o\left(\right|J|}^{2q})\end{array}$$

(73)

(ii)

In the region $q>\frac{1}{2}$, one obtains using Equations (58) and (72):

$$\begin{array}{c}\hfill \overline{S_q}=\frac{{\mathsf{\Omega }}_1\sqrt{\pi }}{1-q}\left[\frac{1}{sin\left(\frac{\pi }{2q}\right)}-q\right]+o\left(J\right)\end{array}$$

(74)

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

$$\begin{array}{c}\hfill \overline{S_1}={\mathsf{\Omega }}_1\sqrt{\pi }+o\left(J\right)\end{array}$$

(75)

In conclusion, the parameter ${\mathsf{\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 ${\mathsf{\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 ${\mathsf{\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,...,i_k}$.

4.1. Statistical Properties of the Couplings $J_{i_1,...,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,...,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,...,i_k}$. To simplify the discussion, let us make the simplest possible choice and assume that the couplings $J_{i_1,...,i_k}$ are independent Gaussian random variables:

$$\begin{array}{c}\hfill P\left(J_{i_1,...,i_k}\right)=\frac{1}{\sqrt{2\pi {\mathsf{\Delta }}^2(i_k-i_1)}}{e}^{-\frac{J_{i_1,...,i_k}^2}{2{\mathsf{\Delta }}^2(i_k-i_1)}}\end{array}$$

(76)

of zero mean:

$$\begin{array}{c}\hfill \overline{J_{i_1,...,i_k}}=0\end{array}$$

(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:

$$\begin{array}{c}\hfill \overline{J_{i_1,...,i_k}^2}={\mathsf{\Delta }}^2(r\equiv i_k-i_1)\end{array}$$

(78)

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

$$\begin{array}{c}\hfill \mathsf{\Delta }\left(r\right)=\frac{v}{r^a}{2}^{-br}\end{array}$$

(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:

$$\begin{array}{c}\hfill B_1^x=-\sum _{k=2}^L\sum _{1<i_2...<i_k\le L}{J}_{1,...,i_k}{S}_{i_2}^x...S_{i_k}^x\end{array}$$

(80)

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

$$\begin{array}{c}\hfill \overline{{\left(B_1^x\right)}^2}\simeq \sum _{k=2}^L\sum _{1<i_2...<i_k\le L}\overline{J_{1,...,i_k}^2}=\sum _{k=2}^L\sum _{1<i_2...<i_k\le L}{\mathsf{\Delta }}^2(i_k-1)\end{array}$$

(81)

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 $\left(\genfrac{}{}{0pt}{}{r-1}{q}\right)$, so that Equation (80) yields:

$$\begin{array}{c}\hfill \overline{{\left(B_1^x\right)}^2}\simeq \sum _{r=1}{\mathsf{\Delta }}^2\left(r\right)\sum _{q=0}^{r-1}\left(\genfrac{}{}{0pt}{}{r-1}{q}\right)=\sum _{r=1}{\mathsf{\Delta }}^2\left(r\right)2^{r-1}\end{array}$$

(82)

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

$$\begin{array}{c}\hfill \sum _{r=1}^{+\infty }{\mathsf{\Delta }}^2\left(r\right)2^r<+\infty \end{array}$$

(83)

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

$$\begin{array}{c}\hfill \sum _{r=1}^{+\infty }{\mathsf{\Delta }}^2\left(r\right)2^r=v^2\sum _{r=1}^{+\infty }\frac{1}{r^{2a}}{2}^{(1-2b)r}\end{array}$$

(84)

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

4.2. Study of the Scale ${\mathsf{\Omega }}_1$

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

$$\begin{array}{c}\hfill \overline{|J_{i_1,...,i_k}|}={\int }_0^{+\infty }dJJ\frac{2}{\sqrt{2\pi {\mathsf{\Delta }}^2(i_k-i_1)}}{e}^{-\frac{J^2}{2{\mathsf{\Delta }}^2(i_k-i_1)}}=\sqrt{\frac{2}{\pi }}\mathsf{\Delta }(i_k-i_1)\end{array}$$

(85)

is also governed by the scale $\mathsf{\Delta }\left(r\right)$ as a function of the spatial range r.

As a consequence, the scale ${\mathsf{\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,...,i_{k_A})$ 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_{k_B-1})$ in the $(j_{k_B}-L_A-1)$ remaining possible positions $(L_A+1,...,j_{k_B}-1)$:

$$\begin{array}{ccc}\hfill {\mathsf{\Omega }}_1& =& \sum _{k_A=1}^{L_A}\sum _{1\le i_1<i_2...<i_{k_A}\le L_A}\sum _{k_B=1}^{L_B}\sum _{L_A+1\le j_1<j_2...<j_{k_B}\le L_A+L_B}\overline{|J_{i_1,...,i_{k_A},j_1,...,j_{k_B}}|}\frac{1}{\sqrt{2(k_A+k_B)W^2}}\hfill \\ & =& \frac{\sqrt{2}}{\pi W}\sum _{i_1=1}^{L_A}\sum _{j_{k_B}=L_A+1}^{L_A+L_B}\mathsf{\Delta }(j_{k_B}-i_1)\sum _{k_A=1}^{L_A-i_1+1}\sum _{k_B=1}^{j_{k_B}-L_A}\frac{1}{\sqrt{(k_A+k_B)}}\left(\genfrac{}{}{0pt}{}{L_A-i_1}{k_A-1}\right)\left(\genfrac{}{}{0pt}{}{j_{k_B}-L_A-1}{k_B-1}\right)\hfill \\ & =& \frac{\sqrt{2}}{\pi W}\sum _{n=0}^{L_A-1}\sum _{m=0}^{L_B-1}\mathsf{\Delta }(1+n+m)\sum _{k_A=1}^{n+1}\sum _{k_B=1}^{m+1}\frac{1}{\sqrt{(k_A+k_B)}}\left(\genfrac{}{}{0pt}{}{n}{k_A-1}\right)\left(\genfrac{}{}{0pt}{}{m}{k_B-1}\right)\hfill \end{array}$$

(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 ${\mathsf{\Omega }}_1$ has always a finite contribution $O\left(1\right)$ coming from the finite distances $(n,m)$ of order $O\left(1\right)$ with respect to the frontier. Now, we wish to evaluate the contribution of large distances $(n,m)$. Using the identity:

$$\begin{array}{c}\hfill \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}\end{array}$$

(87)

one obtains:

$$\begin{array}{ccc}& & \sum _{k_A=1}^{n+1}\sum _{k_B=1}^{m+1}\frac{1}{\sqrt{(k_A+k_B)}}\left(\genfrac{}{}{0pt}{}{n}{k_A-1}\right)\left(\genfrac{}{}{0pt}{}{m}{k_B-1}\right)\hfill \\ & & =\frac{1}{\sqrt{\pi }}{\int }_0^{+\infty }du{u}^{-\frac{1}{2}}\left[\sum _{k_A=1}^{n+1}\left(\genfrac{}{}{0pt}{}{n}{k_A-1}\right)e^{-k_{A}u}\right]\left[\sum _{k_B=1}^{m+1}\left(\genfrac{}{}{0pt}{}{m}{k_B-1}\right)e^{-k_{B}u}\right]\hfill \\ & & =\frac{1}{\sqrt{\pi }}{\int }_0^{+\infty }du{u}^{-\frac{1}{2}}{e}^{-2u}{\left[1+e^{-u}\right]}^{1+n+m}\hfill \\ & & =\frac{1}{\sqrt{\pi }}{\int }_0^{+\infty }du{u}^{-\frac{1}{2}}{e}^{-2u}{\left[2-(1-e^{-u})\right]}^{1+n+m}\hfill \\ & & =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}\hfill \\ & & \underset{n+m\to +\infty }{\simeq }{2}^{1+n+m}\frac{\sqrt{2}}{\sqrt{n+m}}\hfill \end{array}$$

(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\simeq \frac{n+m}{2}$.

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

$$\begin{array}{ccc}\hfill \mathsf{\Delta }(1+n+m)\sum _{k_A=1}^{n+1}\sum _{k_B=1}^{m+1}\frac{1}{\sqrt{(k_A+k_B)}}\left(\genfrac{}{}{0pt}{}{n}{k_A-1}\right)\left(\genfrac{}{}{0pt}{}{m}{k_B-1}\right)& & \underset{n+m\to +\infty }{\propto }\frac{v}{{(n+m)}^{a+\frac{1}{2}}}{2}^{-(1-b)(n+m)}\hfill \end{array}$$

(89)

that we need to integrate over n and m to obtain the scale ${\mathsf{\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 ${\mathsf{\Omega }}_1$ is dominated by the contribution of the long distance:

$$\begin{array}{ccc}\hfill {\mathsf{\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}}}\hfill \\ & \propto & \frac{v}{\left(\frac{3}{2}-a\right)\left(a-\frac{1}{2}\right)}{L}_A^{\frac{3}{2}-a}\hfill \end{array}$$

(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<\frac{3}{2}$, it is possible to construct a critical state with an entanglement growth governed by the exponent:

$$\begin{array}{c}\hfill 0<\alpha =\frac{3}{2}-a<1\end{array}$$

(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 ν 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:

$$\begin{array}{c}\hfill \xi =\frac{1}{b-1}\end{array}$$

(92)

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

$$\begin{array}{c}\hfill \nu =1\end{array}$$

(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 ${\nu }_{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 $\nu =1$ does not satisfy the usual Harris [72] or Chayes–Chayes–Fisher–Spencer [73] inequality:

$$\begin{array}{c}\hfill \nu \ge \frac{2}{d}\end{array}$$

(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^{\frac{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>$ of Equation (22). Then, the full enumeration of possible resonances involves the other $(2^L-1)$ states of energies:

$$\begin{array}{c}\hfill E_{S_1,...,S_L}^{\left(0\right)}=-\sum _{i=1}^L{h}_i{S}_i\end{array}$$

(95)

Therefore, 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 Equation (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 Equation (94), so that there is no theoretical inconsistency in the numerical works where the correlation length exponent ν 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,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 Equation (38), both in the many-body localized phase $b>1$ and at the critical point $b_c=1$.

5.1. Multifractal Exponents $\tau \left(q\right)$

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{array}{ccc}\hfill \overline{S_q}& =& \frac{\mathsf{\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...<i_{k_A}\le L}{\mathsf{\Omega }}_{i_1,...,i_{k_A}}^{2q}+{o\left(\right|J|}^{2q})\hfill \\ & \propto & \sum _{n=0}^{L_A-1}\sum _{m=0}^{L_B-1}{\mathsf{\Delta }}^{2q}(1+n+m)\frac{2^{n+m}}{\sqrt{n+m}}\hfill \\ & \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}}\hfill \\ & \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)}}\hfill \end{array}$$

(96)

For the bipartite case $L_A=L_B=L$, this yields the following power laws in terms of the size $\mathcal{M}=2^{2L}$ of the Hilbert space for the generalized moments:

$$\begin{array}{c}\hfill Y_q^{typ}\left(L\right)=e^{(1-q)\overline{S_q}}\propto v^{2q}{\mathcal{M}}^{-{\tau }_b\left(q\right)}\end{array}$$

(97)

with the multifractal exponents:

$$\begin{array}{ccc}\hfill {\tau }_b\left(q\right)& =& 2qb-1\mathrm{for}q\le \frac{1}{2b}\hfill \\ \hfill {\tau }_b\left(q\right)& =& 0\mathrm{for}\frac{1}{2b}\le q\le \frac{1}{2}\hfill \end{array}$$

(98)

In the region $q>\frac{1}{2}$, Equation (74) yields:

$$\begin{array}{ccc}\hfill {\tau }_b\left(q\right)& =& 0\mathrm{for}q\ge \frac{1}{2}\hfill \end{array}$$

(99)

in continuity with the result of Equation (98).

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

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

$$\begin{array}{ccc}\hfill {\tau }^{criti}\left(q\right)& =& 2q-1\mathrm{for}q\le \frac{1}{2}\hfill \\ \hfill {\tau }^{criti}\left(q\right)& =& 0\mathrm{for}\frac{1}{2}\le q\hfill \end{array}$$

(100)

correspond to the Legendre transform of the multifractal spectrum:

$$\begin{array}{c}\hfill f^{criti}\left(\alpha \right)=\frac{\alpha }{2}\mathrm{for}0\le \alpha \le 2\end{array}$$

(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,49,50,51,52,53,54,55,56,57,58,59,60,61,62,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{array}{ccc}\hfill {\tau }_b^{loc}\left(q\right)& =& 2qb-1\mathrm{for}q\le \frac{1}{2b}\hfill \\ \hfill {\tau }_b^{loc}\left(q\right)& =& 0\mathrm{for}\frac{1}{2b}\le q\hfill \end{array}$$

(102)

is the Legendre transform of the multifractal spectrum:

$$\begin{array}{c}\hfill f_b^{loc}\left(\alpha \right)=\frac{\alpha }{2b}\mathrm{for}0\le \alpha \le 2b\end{array}$$

(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,81,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 $\nu =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 $\nu \ge 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 [38,39].