Signatures of a critical point in the many-body localization transition

Disordered interacting spin chains that undergo a many-body localization transition are characterized by two limiting behaviors where the dynamics are chaotic and integrable. However, the transition region between them is not fully understood yet. We propose here a signature that unambiguously identifies a possible finite-size precursor of a critical point, and distinguishes between two different stages of the transition. The kurtosis excess of the diagonal fluctuations of the full one-dimensional momentum distribution from its microcanonical average is maximum at this singular point in the paradigmatic disordered J1-J2 model. Both the particular value of this maximum and the disorder strength at which it is reached increase with the system size, as expected for a typical finite-size scaling. We completely characterize the short and long-range spectral statistics of the model and find that their behavior perfectly correlates with the properties of the diagonal fluctuations. For lower values of the disorder, we find a chaotic region in which the Thouless energy diminishes up to the transition point, at which it becomes equal to the Heisenberg energy. For larger values of disorder, spectral statistics are very well described by a generalized semi-Poissonian model, eventually leading to the integrable Poissonian behavior.


Introduction
In the classical world the connection between the ergodic properties of a system and its class of regularity (chaotic and integrable) is well understood both mathematically and phenomenologically [1,2]. The classical phase space is completely covered by the erratic trajectories followed by classical chaotic systems while integrability means that only a certain subset defined by the constraints of the constants of motion is accessible. In the first case the system will thermalize after being driven away from thermal equilibrium while in the second case it will generally not. Thus ergodicity and chaotic behavior are generally tied together and one is expected if the other is present.
In the quantum realm, however, the situation is much more involved. Our understanding of thermalization lies in the eigenstate thermalization hypothesis (ETH) [3][4][5][6][7][8]. According to this theory, a quantum system will thermalize if the fluctuations of the long-time averages of local observables around their standard microcanonical equilibrium value decrease fast enough with system size. Plainly put, the ETH states that the diagonal matrix elements of observables in the eigenstates of the system Hamiltonian must change with energy smooth enough, coinciding with the micranonical average up to some random and sufficiently small fluctuation. The connection between quantum thermalization and quantum chaotic behavior [9] is still an open question. The conjecture put forward by Bohigas, Giannoni and Schmit [10] indicates that the spectral fluctuations of quantum systems with an ergodic classical analogue follow exactly the predictions of random matrix theory (RMT) [11]. Hence the statistical analysis of the eigenlevel distribution of quantum spectra has been established as a very valuable tool to investigate quantum ergodicity. The ETH has been long thought to be a consequence of RMT, and thus the thermal properties of quantum systems seem to depend on the onset of quantum chaos in a non-trivial way. As in the classical case, quantum thermalization and chaos are two terms almost invariably associated with each other. However, the limits of applicability of RMT are more rigid than those of the ETH, and thus the latter can still apply beyond the capabilities of the former.
Originally, the RMT provided a solid framework to identify chaos in isolated quantum systems. Perhaps unexpectedly, the last decade has seen an enormous revival of the interest in spectral analysis due to the very exotic nature of quantum many-body systems. In this direction, disordered lattice models from condensed matter have gifted us with an unprecedented playground to test quantum ergodicity [12][13][14]. It seems by now well established that this is a property that clean (i.e., not disordered) quantum many-body systems generally have [4]. Introducing sufficiently strong disorder gives rise to one of the most striking exceptions to quantum thermal behavior: many-body localization (MBL) [3,15], which is now taken as the prototypical nonergodic phenomenon. MBL is observed in some disordered lattice models with local interactions such as the Heinsenberg spin-1/2 chain or its generalization to next to nearest-neighbor interactions, the J 1 -J 2 model, for high enough disorder [16][17][18][19][20][21][22][23][24][25][26][27][28]. Systems in this class often display two completely opposed regimes: an ergodic phase for small disorder strength, in which both the ETH and RMT provide conclusive results, and the MBL phase, where the spectrum does not follow RMT and the system does not equilibrate at all. The region in between these two phases reveals highly unusual properties that have attracted very active investigation [29][30][31][32][33][34]. One of the most pressing questions arguably concerns the nature of the transition between the metallic, thermal phase and the insulating, nonthermal one, i.e., whether there is an actual phase transition in which a critical point can be observed or if what we are dealing with is a smooth crossover instead. Although exploring finite systems rather clearly leads to the second option, the answer in the thermodynamic limit (i.e., at macroscopic scales) is a lot more difficult to obtain, mainly due to the exponential increase of the dimensionality of the Hilbert space with system size in many-body quantum systems. A variety of ergodicity breaking indicators [15,23,24,28,[35][36][37][38][39][40][41][42] have been used in the past to try and identify a hypothetical value of disorder strength that, in the thermodynamic limit, completely separates the ergodic and the nonergodic phases. The search for such a critical point implicitly assumes a nature of real phase transition in MBL systems, although this has not been proved and some apparent contradictions seem to exist. To obtain this value, different kinds of involved finite size scalings of such indicators have been employed. However, the problem is not simple at all and a strong dependence on the particular indicator and the number of sites means that some of the previously obtained values for the critical point are often incompatible with each other.
The aim of this paper is to introduce a new finite-size precursor of the critical point separating the transition between the ergodic phase and the many-body localized phase in disordered manybody quantum systems. The key quantity is the probability of extreme events in the fluctuations from the equilibrium microcanonical value of the diagonal matrix elements of the momentum distribution. These fluctuations are the basic quantities in the ETH. This probability of extreme events is quantified by the kurtosis excess of their ensemble probability distribution. The kurtosis excess is shown to be highest at the transition between the ergodic phase and the MBL phase in the J 1 -J 2 disordered spin chain, coinciding with the result provided by a full analysis of short and long-range correlations of the spectral statistics. The transition point can be determined with precision for each size of the spin chain, allowing a transparent separation of phases irrespective of the Hilbert space dimension. For disorder strengths smaller than the critical one, the system is in its chaotic phase. However, the Thouless energy scale varies greatly depending on the value of disorder until it disappears at the transition. For larger disorder strengths the region can be described by a family of semi-Poissonian statistics only reaching the fully Poissonian limit for infinitely large values of disorder.
The remainder of this paper is structured as follows. In Section 2 we introduce the J 1 -J 2 spin chain used in the computations. In Section 3, we present the concepts of quantum thermalization and ETH in subsection 3.1 which is at the core of our research. In the subsection 3.2 we introduce the diagonal fluctuations of the momentum distribution from their microcanonical averages and compute the results as a function of disorder, which is the main result of our paper. In Section 4 we study the spectral statistics across the transition introducing a full description including longrange correlations. In the first subsection 4.1 we present the semi-Poisson spectral statistics that is followed by chains in the MBL phase; then, in the second subsection 4.2 we complement the previous results with the study of the Thouless energy E Th and long-range correlations on the ergodic side of the transition; finally, in subsection 4.3 we provide a general landscape of the MBL and ergodic phases in terms of their spectral properties. In Section 5 we relate the results of the spectral statistics with the behavior of the diagonal fluctuations presenting a coherent interpretation of our results. We gather the main conclusions of our work in Section 6.
The J 1 -J 2 model consists on a one-dimensional chain with L sites, on-site magnetic fields ω , and coupling parameters J 1 and J 2 where next to nearest-neighbors interactions are additionally considered. Its Hamiltonian can be cast in the form whereŜ x,y,z are the total spin-1/2 operators at site ∈ {1, . . . , L}. The Hamiltonian is simulated on a lattice with periodic boundary conditions as these minimize finite-size effects; thus, S x,y,z =Ŝ x,y,z +L . Disorder enters the chain via the uniformly, independently and randomly distributed magnetic fields ω ∈ [−ω, ω]. The coupling constants J 1 and J 2 are associated to nearest and next to nearest-neighbors interactions, respectively, and they can be used to set the unit of energy. In our simulations, we fix J 1 = J 2 = 1 throughout -since J 1 , J 2 ≥ 0, we are dealing with the anti-ferromagnetic variant of the model . Note that H(J 1 = 0, J 2 = 0) is simply the famous XXZ Heisenberg chain. The terms λ 1 and λ 2 quantify the intensity of the interactions, which we choose λ 1 = λ 2 = 0.55 (this choice is also made, e.g., in [28]). Spin chains such as this one can be mapped exactly onto a spinless fermionic chain with an extra boundary term via the Jordan-Wigner transformation [44].
The crux of the matter of systems that reveal a MBL transition is that both their dynamical and static physical properties critically depend on the intensity of the disorder strength, ω. Although the precise boundaries have not been completely delimited and in fact they vary depending on each particular model, the general scenario is more or less shared by all of them. While in the absence of disorder the XXZ chain is solvable by Bethe ansatz and thus integrable, the J 1 -J 2 chain is considered as a paradigmatic model to exhibit quantum chaos and quantum thermalization [4].
For intermediate values of ω, the chain shows an ergodic phase where most initial conditions are expected to thermalize. This region is generally said to be quantum chaotic in the sense that energy levels are characterized by level repulsion and the eigenvalue distribution is close to the Gaussian orthogonal ensemble (GOE) as in the predictions of RMT. This picture can be easily checked by statistical measures such as the nearest-neighbor spacing distribution (NNSD), P (s), or the adjacent eigenlevel gap ratio, P (r). However, as we will show later, this metallic region hides long-range deviations from RMT universal results, which can be analyzed by computing the Thouless energy scale as was done in Refs. [23,45]. Dynamically, this region is dominated by sub-diffusive processes, multifractal scalings and other unexpected behavior [29][30][31]. The region in between the ergodic and the localized phases has been extensively studied. Recent results indicate that it may be characterized by a Griffith-like phase in which anomalously different disorder regions seem to dominate the dynamics [32][33][34]46]; nonetheless, the debate is still open [47]. To describe the flow of intermediate statistics observed in this region, mean-field plasma models with effective power-law interactions between energy levels [24,48], the Rosenzweig-Porter ensemble with multifractal eigenvectors [49,50], a family of short-range plasma models [51] and generalizations [41,42] have been used. Finally, for disorder strengths larger than a critical value that is dependent on the dimension of the Hilbert space, the chain gradually reaches the MBL phase. Here, the ETH is violated, so generic initial conditions do not relax to their microcanonical average value. This region shows Poissonian spectral statistics instead, from which it follows that the spectrum behaves as a set of uncorrelated random numbers where level crossings can potentially take place [52]. This is explained by the identification of a complete set of local integrals of motion in the MBL phase [53,54]. Whether the MBL phase is reached in the thermodynamic limit from the ergodic phase in the form of an actual phase transition (i.e., by means of a critical point), or as a dynamical crossover, exhibiting a finite transition region (dubbed bad metal) with surviving non-ergodic but extended states is an open question [38,[55][56][57]. The answer seems to be forever eluding the community due to the (strictly speaking) impossibility to access the thermodynamic limit and the importance of finite size effects [58,59], which are strong in these chains.
As the operatorŜ z := L =1Ŝ z commutes with the Hamiltonian, in this work we restrict our attention to the sector S z = 0. The total dimension of the Hilbert space is then d = L L/2 which grows asymptotically when L → ∞ as d ∼ 2 L / πL/2. Thus, exact diagonalization with full calculation of all the eigenstates becomes realistic only up to chain lengths 16 L 18. The spin chains that are frequently used to investigate the MBL phenomenon are known to exhibit strong border effects, meaning that the eigenstates are more localized as they get closer to the boundaries of the band, regardless of the disorder strength [36,38,60,61]. The question of whether true mobility edges survive in the thermodynamic limit in MBL systems is still open as they have also been argued to be indistinguishable from finite size effects [62]. In any case, to avoid border effects we will only consider the central N = d/4 eigenstates {|E n } N n=1 . The number of states we use range from N = 63 for L = 10 to N = 3217 for L = 16.
3 Extreme events across the transition 3.1 Quantum thermalization: the concept Quantum thermalization refers to the equilibrium state that a quantum system reaches for sufficiently long times. Let us consider an isolated quantum system with Hamiltonian H with eigenstates {|E n } n satisfying H |E n = E n |E n , and an arbitrary initial condition |ψ(0) . The latter evolves in time under H and after an interval of t ≥ 0 its wavefunction can be written |ψ(t) = exp(iHt/ ) |ψ(0) . For a typical observableÔ, one might consider the long-time average Ô t in the eigenbasis of H. If, for simplicity, the spectrum is non-degenerate, this is given by where |C n | := | ψ(0)|E n | is a c-number representing the probability of finding the system in the eigenstate |E n . Under very relaxed assumptions [63,64], long-time averages of the kind of Eq.
(2) are known to reach a certain equilibrium value and remain close to it at all times. Nonetheless, the system is often said to thermalize only if that value corresponds to the particular case of the microcanonical average, where E is the macroscopic energy of the system and ∆E is a small energy window, ∆E/E 1, containing a large but finite number of levels, 1 N < ∞. The connection between Eqs. (2) and (3) is well understood in classical dynamics; in particular, for classical chaotic systems, long-time averages are equivalent to phase averages if one fixes the right energy E, whereas integrable systems generically do not relax to Eq. (3) at all. In quantum mechanics, the eigenstate thermalization hypothesis [3][4][5][6][7][8] underlies the equivalence between these two averages. It concerns the diagonal matrix elements of a typical observable O nn := E n |Ô |E n , which can always [8] be written where N denotes the Hilbert space size. The values ∆ n are called diagonal fluctuations, and describe the fluctuations of O nn around the microcanonical average Ô ME . We can consider them as random numbers verifying ∆ n = 0 and ∆ 2 n = 0. The ETH states that a quantum system thermalizes if ∆ n (or, rather, its standard deviation, σ ∆n = ∆ 2 n ) decreases fast enough with N . In other words, the ETH is nothing more than a statement about how smoothly the diagonal terms must change with energy for a quantum system to reach the equilibrium value Eq. (3) for sufficiently long (but not exponentially long in the system size) times.
The standard tool to determine whether a quantum system thermalizes or not consists in studying if σ ∆n decreases fast enough with the system size. However, besides this fact, it has been recently shown that the presence of correlations in the ∆ n constitutes a signature of the existence of anomalous non-thermalizing initial conditions [45]. To better study these correlations with different observables, it is convenient to make the diagonal fluctuations in Eq. (4) dimensionless, by normalizing by its standard deviation, which yields The consequence is that the new quantity in Eq. (5) is exactly a standard Gaussian random variable G(0, 1), that is, with expected value 0 and variance 1, for any observable in any quantum system obeying RMT. Hence, as previous evidence strongly suggests that the breakdown of the ETH close to integrable regions may be a generic result [65][66][67][68][69][70][71][72], whereas thermalization is widely associated with ergodicity and and RMT, the quantity ∆ n sharply discriminates between these two limiting regularity classes: for quantum chaotic systems it behaves as an uncorrelated, white random noise and is Gaussian, whereas for quantum integrable ones an emerging structure in its Fourier modes results from the existence of integrals of motion (i.e., the noise is no longer featureless). Thus, the (normalized) diagonal fluctuations ∆ n provide a powerful tool to identify small deviations from thermalizing behavior, even in otherwise rather chaotic regions. In this direction, the power spectrum of ∆ n has been recently used in conjunction with long-range spectral statistics to connect thermalization and Thouless energy in the prototypical disordered XXZ spin-1/2 chain [45].

Analysis of the diagonal fluctuations
Here, we study the diagonal fluctuations Eq. (5) across the MBL transition in the J 1 -J 2 model, Eq. (1). As representative physical observables, we choose the full momentum distribution on a one-dimensional lattice with lattice constant set to unity, i.e., where := 1 andŜ ± are the usual ladder spin operators, which are related to those in Eq. (1) by the well-known expressionsŜ ± =Ŝ x ± iŜ y . Thus, to calculate ∆ n , we first need to evaluate the diagonal matrix elements O nn , for which the complete set of eigenstates {|E n } n is needed. This is obtained by full diagonalization of Eq. (1). Next, the microcanonical average Ô ME is obtained by fitting a polynomial of degree 4 to the previous matrix elements O nn to avoid the spurious effects originating in averages over finite energy windows [73] (this is the same procedure that was used in Ref. [45]). For reference, a summary of the number of central states N = d/4 and the number of realizations for each value of the number of sites L and the disorder strength ω can be found in Table  1. These quantities have been chosen explicitly so that their product remains approximately constant. As not only the eigenlevels {E n } n but also the eigenvectors {|E n } n are needed, we cannot realistically increase the value of L beyond 16. On the other hand, system sizes below L = 10 are indeed accessible but not very representative as the samples become statistically poor.
We can now compute the probability distribution of such quantities ∆ n , P ( ∆ n ). It is then straightforward to obtain the kurtosis excess which is directly related to the probability of extreme events. For a random variable X, the kurtosis excess γ 2 (X) is defined as L Levels in the central region Realizations Here, µ := X is the expected value of X while σ 2 := X 2 − X 2 is its variance. The kurtosis of any univariate Gaussian distribution is 3, and thus we may use this quantity to compare how important is the presence of extreme events in a certain probability distribution with respect to a Gaussian. A positive kurtosis excess means that the distribution of X is more tailed than a Gaussian while a negative kurtosis excess indicates the opposite. In what follows, the random variable of interest is X = ∆ n .
In the main panel of Fig. 1 we show the kurtosis excess as a function of disorder for L = 10, 12, 14, 16. For small values of ω, those corresponding to the ergodic phase of the model, the probability of extreme events roughly corresponds to the expected for a Gaussian distribution, γ 2 ( ∆ n ) ≈ γ 2 (G) = 0; the larger the chain size, L, the better the agreement. This means that the observables in Eq. (6) behave as expected in RMT in this region, and therefore ETH is fulfilled and thermalization is expected. Note that ω = 0 has been excluded from our analysis as the clean system, a spin chain with periodic boundary conditions, is translationally invariant. The influence of these extra symmetries survives up to larger values of ω for smaller L and γ 2 ( ∆ n ) increasingly deviates from 0 as a consequence. Then, the probability of extreme events as measured by γ 2 ( ∆ n ) shows a neat increase with disorder. However, this increase is not monotonic. On the contrary, γ 2 ( ∆ n ) has a maximum at a certain value of the disorder, which depends on L; this disorder strength will be denoted by ω c (L) from now onward. In fact, we take the preceding statement as the definition of ω c (L). After reaching this peak, γ 2 starts decreasing as a function of disorder, becoming negative for the largest values of disorder considered.
A precise calculation of ω c (L) is a complicated task, requiring a huge number of realizations. To illustrate this fact, in the inset of Fig. 1 we plot the same results for L = 10, obtained with 10 3 and 10 5 realizations. The last curve is smooth and provides a pretty good value for ω c (L = 10). On the contrary, the curve obtained with 10 3 realizations shows large fluctuations. As collecting so many realizations is too expensive for larger values of L, we provide here just an estimate of the singular disorder strength, ω c (L), for the cases displayed in the main panel of Fig. 1: Despite these uncertainties, it is clearly seen that both the disorder strength for the maximum in γ 2 ( ∆ n ), ω c , and its magnitude increase with the system size. Furthermore, as the dimension of the Hilbert space grows, the probability of finding extreme events is increasingly higher for intermediate and large ω but also smaller for small ω.
At this stage, we can draw the following picture for the transition from the ergodic to the MBL phases in the J 1 -J 2 model. This transition starts with an increase of extreme events for the diagonal fluctuations, ∆ n , which favor the existence of non-thermalizing initial conditions [45]. However, the integrable MBL phase is not characterized by a large probability for such extreme events, but by a large value of σ ∆n with a negative kurtosis excess. Hence, we can write a preliminary  Figure 1: Main panel: kurtosis excess, γ 2 ( ∆ n ), Eq. (7), as a function of the disorder strength ω for the number of sites L ∈ {10, 12, 14, 16}. The black, dashed line represents the kurtosis excess for a standard Gaussian distribution G(0, 1), γ 2 (G) = 0. Inset: kurtosis excess for L = 10 obtained with 10 5 realizations (red squares) and 10 3 realizations (yellow circles).
statement for the main conclusion of the paper: ω c (L) is a singular point in the transition from the ergodic to the MBL phase in the J 1 -J 2 model, characterized by a maximum probability of extreme events for the diagonal fluctuations. Furthermore, its behavior when the chain size L is increased resembles the finite-size scaling of observables that diverge at the critical point of a standard phase transition, like, say, the magnetic susceptibility.
From this result, one may ask about the actual shape of the distribution P ( ∆ n ). This matter is addressed in Fig. 2, which consists of two panels. Here we fix the number of sites at L = 16. On the right-hand side we show P ( ∆ n ) for three representative values of disorder. We also plot with black, dashed lines the probability density function of a Gaussian G(0, 1), P (x) Gaussian := exp −x 2 /2 / √ 2π, x ∈ R, which underlies this quantity in the case of thermalizing systems.
For ω = 1, we can see that this is indeed the case: P ( ∆ n ) matches perfectly such a Gaussian distribution. It is worth to mention that no fitting has been performed here; we merely plot the Gaussian on top of the distributions. For ω = 4.7, which lies inside the interval 1 for ω c (16) inferred from Fig. 1, we observe that P ( ∆ n ) no longer agrees with a Gaussian, and its tails display a slower decay, in concert with the high value of γ 2 ( ∆ n ) at this disorder. Instances of long tailed distributions and the breakdown of the ETH near the MBL transition are known, often associated with Griffith effects [33,34,40,45]. Finally, for a value of disorder well within the localized phase, ω = 100, the distribution has been completely distorted and the tails decay faster than those of a Gaussian. To get a clearer picture, on the left-hand side of Fig. 2 we have chosen to plot the tails of the same distributions for ∆ n ∈ [2,6]. This panel further confirms that the case ω = 4.7 decays much more slowly than the others; as ω → ∞, the decay is faster and faster; and ω = 1, matching very approximately the tails of the G(0, 1) distribution, seems to be in an intermediate situation.
This discussion reinforces our previous conclusion. ω c (L) may be associated with a precursor of a critical point in the ergodic-MBL transition. Unfortunately, much larger systems, far above the reach of current computers, are required to determine if this is an actual transition, or just a smooth crossover. In any case, from the results discussed in this section, we can conjecture that, if there exists a critical point, ω c (∞), then it must be the one in which γ 2 ( ∆ n ) is maximum . In the next section we will show that the singular character of ω c (L) is fully compatible with the behavior of spectral statistics across the transition, widely used to characterize the static properties of the chain.

Spectral statistics across the transition
The statistical analysis of eigenlevel fluctuations is by now well established as one of the main tools to study quantum complex systems. This is because the spectral properties of quantum systems with a chaotic classical counterpart do not depend on the particular features of the Hamiltonian but only on its more general global symmetries. As a result, the level statistics of a variety of apparently completely unrelated physical systems universally coincide with the predictions of RMT for the random matrix ensemble associated to their particular symmetry [10]. The eigenlevel distribution of the kind of spin chains associated to the MBL transition has been abundantly studied by many previous works. The goal of this section is to show that the maximum of the probability of extreme events as measured by the kurtosis in Fig. 1 acts as an indicator of the value of disorder beyond which the chaotic regime is completely abandoned and the transition to the MBL phase begins. We will show that the spectral statistics are not only quantitatively but also qualitatively different on both sides of this point. For ω < ω c (L), the system belongs in the chaotic region where the Thouless energy, the scale associated with RMT correlations, is larger than the Heisenberg energy, the mean energy distance between levels. To model the spectral statistics for larger disorder strengths, ω > ω c (L), we employ a family of generalized semi-Poisson statistics.
Here we make use of two different measures: the NNSD, P (s), which is concerned with shortrange level correlations, and the δ n , which is a long-range spectral statistic. For the universal predictions of RMT to hold it is necessary to obtain the cumulative spectral function, where Θ is the unit step function, representing the number of levels with energy less than or equal to a certain value E [9][10][11]. This function can be separated into a smooth part G and a fluctuating part, G, as G(E) = G(E) + G(E). Then, the original eigenlevels {E n } N n=1 are mapped onto the dimensionless quantities {ε n } N n=1 as E n → G(E n ) =: ε n . As a consequence, the mean level density is unity, and RMT results can be applied. This procedure is termed unfolding [73]. Since there is no statistical theory for Eq. 1 that provides a theoretical expression for the smooth G(E), in this work we numerically obtain it by fitting a polynomial of degree 10 to the original energies {E n } N n=1 , i.e., G(x) = 10 k=0 a k x k . The famous NNSD, P (s), is the distribution of the (unfolded) consecutive level spacings, For systems that are invariant under orthogonal transformations, the eigenlevel distribution follows the results for GOE random matrices, and quantum chaos manifests via level repulsion in the Wigner-Dyson surmise P (s) = π 2 s exp −πs 2 /4 . By contrast, the spectrum of quantum integrable systems is equivalent to a set of independent, identically distributed, Poisson random variables [52] where level repulsion is no longer present, and thus P (s) = exp(−s).
Long-range statistics describe the spectral properties of eigenlevels separated by large energy index distances (i.e., levels E i , E j ∈ {E n } N n=1 with |i − j| N ), as opposed to level distances of one or two units (or in general |i − j| N ). Contrary to short-range statistics, these allow to obtain the so-called Thouless energy scale, E Th : the energy scale beyond which universal RMT results break down [74,75]. It has been recently shown that, at least in the disordered XXZ spin chain, this scale determines to what extent a given system thermalizes [45]. A decreasing value of E Th was linked to an increase of the probability of extreme events from the ergodic region. RMT results have been shown to describe well the statistics of eigenvalues separated by less than E Th ; however, for level distances larger than E Th , long-range deviations towards the corresponding integrable result can be detected even in the region interpreted as chaotic by short-range statistics such as the P (s). In the case of disordered spin chains, the number variance [23], the spectral form factor [58], and very recently the δ n [45] have allowed to calculate the value of E Th .
In this work we will use the δ n spectral statistic as a long-range measure of level correlations [76][77][78][79][80][81]. It is defined as the difference between the nth unfolded level and the corresponding energy value in an equiespaced spectrum (i.e., with ε n = n), that is, δ n := ε n − n, n ∈ {1, . . . , N }.
The quantity δ n can be understood as a discrete time series where the level order index n plays the role of a discrete time. Thus a discrete Fourier transform can be applied to δ n , one of the most common techniques in time series analysis. Taking its square modulus gives the power spectrum P δ k of the signal. The quantity of interest is the averaged power spectrum, which reads where k Ny := N/2 is the Nyquist frequency. Here the angular brackets · denote average over realizations for a fixed disorder strength. Analytical expressions are known for the main random matrix ensembles [77]. In the limit k/N 1 and N 1, the power spectrum of δ n in quantum integrable systems exhibits the neat power-law decay P δ k 1/k 2 , whereas for quantum chaotic ones this is P δ k 1/k [76]. One says that this is a universal feature because it is only dependent on the regularity class (integrable or chaotic) of the Hamiltonian matrix; however, it does not depend on the underlying symmetry.
In this section we show results for the number of sites L = 16 in all cases. Smaller sizes reveal the same features. However, we have focused on the larger size that shows best statistical significance, especially when dealing with long-range correlations.

Semi-Poisson model for short-and long-range spectral statistics
While a quantum chaotic spectrum exhibits correlations between levels, these are completely absent from quantum integrable ones (i.e., levels are uncorrelated). Wigner-Dyson and Poisson statistics are valid, respectively, in the limiting metallic and insulator regimes of the metalinsulator transition (MIT) of the Anderson model [82]. At the MIT critical point a new universal model of spectral statistics is believed to be applicable, the semi-Poisson [75], which is connected to the fractal nature of the critical wavefunctions, intermediate between the extended and localized cases [48,83]. Semi-Poisson spectra are intermediate between the two opposed regimes of Wigner-Dyson and Poisson in the sense that they reveal some degree of level repulsion but for asymptotically large level distances the NNSD decreases much more slowly than for GOE spectra, producing P (s) = 4s exp(−2s) [84,85].
It has been shown [86] that a simple way to generate semi-Poisson statistics from a given Poissonian spectrum is to keep every two eigenvalues from it so that the new spacing is simply s i = (s i + s i+1 )/2. If instead one keeps every η ∈ N eigenvalues from a Poissonian spectrum, one gets the new sequence of spacings {s i } i where From now onward we continue to denote the spacings defined in Eq. (10) by s i . The sum of η (unitary) exponential random variables (i.e., Poisson in the sense of Berry and Tabor) has the Erlang probability density, P (s) = s η−1 exp(−s)/(η − 1)!. In particular, the case η = 2 corresponds to semi-Poisson, while η = 1 gives the original Poissonian statistics. This can be further generalized by extending the domain, η ∈ N → [1, +∞). The corresponding NNSD is simply the marginal probability density of any spacing in Eq. (10), where Γ(z) = ∞ 0 dt t z−1 e −t is the gamma function. In Ref. [84], an equivalent family of distributions was obtained from a phenomenological short-range plasma model in order to fit the spectrum of a variety of systems displaying dynamics between chaotic and integrable [87][88][89]. In fact, Eq. (11) is the exact same as Eq. (5) in [84] with η := β + 1 and n = 1, the NNSD of the well-known short-range plasma model based on a one-dimensional Dyson gas introduced by Bogomolny and co-workers. As soon as the continuous parameter η of Eq. (11) departs from its corresponding value in a Poissonian spectrum, the NNSD reveals level repulsion proportional to P (s) ∝ s η−1 , but the fall-off exp(−ηs) is always slower than exp −s 2 in the Wigner-Dyson surmise of chaotic systems. In addition, the property of level spacing independence is preserved in a so-defined spectrum, as the covariance between any two spacings s i and s j is where δ ij is the Kronecker delta. This property of independent level spacings can be used to derive the power spectrum of δ n exactly for the family of distributions in Eq. (11). The exact result [90] is given by for unfolded sequences each one having s = 1 and in the case where the total number of levels is N (i.e., there are N − 1 level spacings; we are of course concerned with the case N 1). Equation (13) depends on the single continuous parameter η ∈ [1, +∞), with η = 1 again corresponding to Poissonian statistics, and η = 2 being semi-Poisson in the strictest sense of the term. From Eq. (13) it follows immediately that the crossover between these two limits takes the form of a vertical translation of P δ k but its overall structure remains unchanged. Similarly, Eq. (12) implies linear behavior of the number variance Σ 2 (L) with a proportionality constant η −1 , i.e., Σ 2 (L) ∼ L/η when L → ∞. In the MIT transition this generalized semi-Poisson model seems to describe the spectral statistics of the critical region with a value of η that changes from 2 to 1 as the dimensionality is increased as has been numerically investigated up to six spatial dimensions [91].
Therefore, Eq. (11) and Eq. (13) allow to completely characterize a crossover from semi-Poissonian to Poissonian statistics. Although the semi-Poissonian limit strictly refers to the case η = 2 alone, one usually still uses the term to allude to the whole family of distributions. As the value of η is reduced the family of distributions exhibit a decreasing intensity of level repulsion, which completely vanishes only at the Poissonian limit.
Here we show that the semi-Poisson model provides a characterization of the level statistics for ω > ω c (L) in Eq. (1), for both short-and long-range correlations. To that end, first we have computed the NNSD for the eigenlevels of our model, Eq. (1). We have unfolded the sequence of original energies {E n } N n=1 with a polynomial of degree 10. Starting from the central d/3 levels, we have removed the 2 d/48 levels closest to the edges before and after unfolding, which gives N = d/4 for our analysis (see Table 1). To avoid the spurious effects introduced by spectral unfolding in systems close to integrable dynamics, we have divided each level spacing by its mean value for each realization, so that the mean level spacing becomes s = 1 exactly [90]. Then, we have constructed the histograms of P (s) with a bin size ds = 0.1. After that, we performed a single-parameter fit of Eq. (11) to the P (s) to obtain the value of η. Results are shown in Fig.  3 for a set of disorder values and L = 16. As can be seen, P (s) can be very approximately characterized by Eq. (11) for disorder values larger than ω c (16). Incidentally, for ω = 5.0, we obtain almost fully semi-Poissonian level spacings, η ≈ 2.079. For ω = 4.7, we obtain η ≈ 2.274, implying that the level repulsion is stronger that in the ergodic case. However, it is difficult to infer whether this is a spurious result due to the proximity to the singular point, ω c (L), or if it is just a finite-size effect. As ω is increased into the localized phase, the parameter η decreases, until it reaches a value of η ≈ 1.058 for ω = 12. Well within the localized phase, the NNSD agrees well with η ≈ 1 (not shown). For ω smaller than ω = 4.7, the system enters the chaotic region. Equation (11) is explicitly derived on the assumption of statistically independent spacings as in the Berry-Tabor result [52]. Then, since quantum chaotic eigenlevels are statistically correlated, Eq. (11) cannot satisfactorily describe this side of the transition. However, for ω 4.7, level spacings approximately become independent random variables, but as we have seen they differ from (generic) integrable systems in that their distribution is not Poissonian.
This picture is also obtained in long-range measures of level statistics. This means that this side of the transition can be correctly described by spectral statistics ranging from semi-Poisson to Poisson for eigenergies separated by not only short but also large level distances.
In Fig. 4   GOE results are also shown for reference. For brevity we omit the analytic expression of the GOE but it can be found in Ref. [77]. For ω 5.0, the agreement between the two curves is almost perfect. Below this disorder strength, for ω = 4.7, the power spectrum is slightly touching the GOE result. This reflects that this value of ω can indeed be taken as the limit beyond which the model of statistically independent spacings is valid. As ω is increased, η decreases towards the Poisson result, and the power spectrum undergoes a smooth crossover, approaching the theoretical Poisson curve vertically. We stress that no fitting has been performed in this case; the black solid lines merely correspond to Eq. (13) for η as obtained in Fig. 3, for each value of ω there shown. The first frequencies are spoiled in all six cases, but this is an expected consequence from the unfolding procedure [73], not a shortcoming of our model.

Long-range statistics on the ergodic region
Short-range spectral statistics in the chaotic region show results consistent with GOE random matrices as has been extensively shown before (see, e.g., [38]), so we will not consider them further. Additional investigation of the long-range spectral correlations across the transition is provided in Fig. 5, where we focus on the chaotic side of the model Eq. (1) instead. There we show the numerical P δ k , obtained following the same steps as before. For ω = 0.5 and ω = 1.0, the differences between the chaotic theoretical curve and the numerics are minimal. However, as ω is increased, the power spectrum gradually deviates from the ergodic towards the integrable curve. This effect is blurred in panels (a) and (b), due to the unfolding procedure [73], but it is clearly seen in the rest of the panels. This deviation is linked to the Thouless energy E Th , the energy scale beyond which energy levels are no longer correlated like in RMT. Following the same procedure that in [45], the power spectrum of δ n gives us direct access to the Thouless frequency  k Th , which determines a characteristic length Th = N/k Th : two energy levels, E n and E m are correlated like in RMT if their level index distance satisfies |n − m| < Th . In this sense, a fully chaotic RMT spectrum is one that has the highest possible value of Th (or the lowest possible value of k Th ). This would indicate that GOE correlations are shared by levels separated by any distance within the spectrum boundaries. As k Th increases, the spectrum is thus 'less chaotic' in this particular sense 2 . A good estimation for k Th is to choose the lowest possible frequency for which P δ k fluctuates below the GOE curve, which gives the approximation k min ≈ k Th . This point is identified in all the panels of Fig. 5 by means of a vertical arrow. We can see that k min monotonically increases with ω, showing that the system becomes less chaotic as the singular point, ω c (L), is approached. It is worth to note that this degree of detail cannot be achieved by analyzing short-range spectral statistics, even in the case where k min is very large. Because they are by definition insensitive to the spectral properties of distant levels, short-range statistics would still produce a rather chaotic result (this can be seen, e.g., in the mean value of the adjacent level gap ratio, sometimes employed for finite-size scaling considerations).

Transition landscape
Results plotted in Figs. 4 and 5 suggest a scenario summarized in Fig. 6. In panel (a) of this figure, we display k min as a function of ω. We observe that the value of k min for ω ω c (L) is very small compared to the number of levels, k min /N 1, which gives a large value of the Thouless energy. It is seen that k min then grows quite fast with ω, explaining the subsequent separation from the GOE curve of the power spectrum for those values of disorder. The limiting value k min = k Ny = N/2 is also shown with a dashed line in panel (a) of Fig. 6. Reaching the Nyquist frequency indicates that the power spectrum has completely separated from the GOE curve, and hence the quantum correlations of RMT are destroyed at all scales (i.e., there is no level index distance such that levels separated by that distance are correlated). Importantly, for ω = 4.5 the characteristic frequency is still smaller than the Nyquist frequency, k min < k Ny , but for ω = 5.0 it has already reached it. This means that P δ k completely separates from the RMT result at a value of the disorder compatible with ω c (L).
From this point onward, the power spectrum is characterized by Eq. (13) instead because level correlations are no longer present. This is shown in panel (b) of Fig. 6, where we display η as a function of ω. In concert with previous Fig. 1, η smoothly decreases as the MBL phase is approached. Furthermore, a value close to η = 2, which implies a level repulsion equal to the GOE, is found around the singular point ω c (L). As we have discussed before, our results are compatible with a strange singular point ω c (L) characterized by semi-Poisson statistics with level repulsion larger than in the GOE. Notwithstanding, this may be also a finite-size effect. In any case, we find a neat transition in the spectral statistics between full or partial correlations with level repulsion [the chaotic region, 0.5 ω ω c (L)], to no correlations but still level repulsion [semi-Poisson region, ω c (L) ω < ∞], to finally no correlations and no level repulsion whatsoever (Poisson limit, ω → ∞).
All these results reinforce our previous conclusion regarding the singular character of ω c (L). Furthermore, the disorder value at which the kurtosis of the diagonal fluctuations, γ 2 ( ∆ n ), is maximum , ω c (L), is also a singular point regarding spectral statistics.

Discussion of results
The results from previous Secs. 3 and 4 provide an interesting insight into the mechanism for the transition between the chaotic and the localized phases in the J 1 -J 2 model. It can be summarized as follows: If there exists a critical point in the transition from the ergodic to the MBL phases, ω c (L → ∞), then it must correspond to the disorder strength at which the kurtosis of the distribution of the (normalized) diagonal fluctuations γ 2 ( ∆ n ) is maximum . For finite systems, ω c (L < ∞) unambiguously splits the chaotic and localized phases.
In Fig. 7 we illustrate this scenario, exemplified with L = 16 but that we conjecture to hold with generality (at least for finite sizes). There, we find two different regions: • Chaotic region. This corresponds to small values of disorder, ω < ω c (L). Here spectral statistics coincide with RMT up to a certain characteristic length max , beyond which RMTlike correlations are lost [45]. For very small values of ω, the ETH is fulfilled, and generic observables relax to their microcanonical equilibrium value. As can be seen in Fig. 1, the probability of extreme events is very approximately the same than for a Gaussian distribution, which underlies the diagonal matrix elements for thermalizing systems. Even more, panels (a) and (b) of Fig. 2 show that ∆ n agrees almost perfectly with such a distribution, including in the tails. As ω is increased, the kurtosis excess γ 2 ( ∆ n ) increases too, significantly separating from the Gaussian expectation. As was shown in Ref. [45], this means that, although generic initial conditions will thermalize to the microcanonical average for these disorder strengths, it is also a lot more probable to find anomalous initial conditions that do not. In terms of spectral statistics, even though for small ω the eigenlevel distribution is very close to the GOE predictions of RMT, Fig. 5 shows that as ω increases (but still remains on the chaotic side) long-range deviations can be attributed to the anomalous behavior of the Thouless energy. The increasing values of the probability of extreme events are thus connected to the gradual loss of chaos in the spectrum. For small values of ω, ∆ n is Gaussian, the system thermalizes, and GOE level correlations are maintained between levels separated by distances comparable to the total size of the spectrum. As ω increases, extreme events take place with more and more probability and the Thouless energy starts decreasing, meaning that level correlations between levels very far apart from each other are being destroyed. For disorder values close to ω ≈ ω c (L), the Thouless energy is minimal and the model departs from its chaotic phase.
• Semi-Poisson region. This corresponds to ω > ω c (L). At this point the probability of extreme events of Fig. 1 starts diminishing; at around ω ≈ 11, it crosses the corresponding value for a Gaussian distribution. As can be seen in panels (a) and (d) of Fig. 2, well within the localized phase, at ω = 100, the distribution of the diagonal fluctuations has been completely distorted and for asymptotically large values it decreases faster than a Gaussian distribution, in concert with the finding that the probability of extreme events is smaller than that of a Gaussian for large disorder. On this side of the transition, Eqs. (11) and (13) account for both short and long-range spectral statistics, as can be seen in Figs. 3 and 4. This means that the spectrum is here approximately composed of independent, identically distributed random numbers that still show level repulsion, so they are intermediate between GOE and Poisson. For L < ∞, the Poisson limit is only strictly reached when ω → ∞.
The singular point separating these two regions, ω c (L), shows the following features. First, it is the disorder strength for which the maximum probability of extreme events in the diagonal fluctuations occur. Here ∆ n is no longer well described by a Gaussian, and the decay of its tails is much slower, almost exponential as panel (a) of Fig. 2 suggests. And second, it indicates certain singularity in the spectral statistics in the sense that below it level correlations exists between levels separated by certain distances, but beyond it no such feature can be found, even though some degree of level repulsion is still preserved.
Thus, we observe that the kurtosis excess γ 2 ( ∆ n ) together with the spectral analysis provided in the previous sections strongly suggest that the flow from Wigner-Dyson to Poisson statistics is a two-stage process, in concert with previous numerical findings [24,30]. In Ref. [30] the transition from the ergodic to the MBL phase was identified by means of a nonuniversal jump of the multifractal dimensions (both in Fock and spin configuration basis). We find that a similar effect gives rise to a maximum value of γ 2 ( ∆ n ), which also hints towards the existence of a critical transition point in the J 1 -J 2 model as we have shown. The multifractal dimensions vanish only in the infinite disorder limit, coinciding with full Poissonian statistics where level repulsion completely vanishes as well since η = 1 [see Eq. (11)]. Our results are also consistent with those presented in Ref. [24], where the effective interaction between eigenlevels in the disordered XXZ spin chain was analyzed. In the ergodic phase, level statistics were characterized by a longrange plasma model. However, upon reaching the MBL transition, a power-law local interaction between levels means that these are intermediate between Wigner-Dyson and Poisson, leading to the family of semi-Poisson distributions as we have seen. The locality of interactions on this side of the transition is also consistent with the level repulsion P (s) ∝ s η−1 , which becomes increasingly weaker as we approach the Poisson limit and hence the interaction between level spacings also diminishes up to this point. Other more common signatures of the transition from ergodicity to MBL, like the mean value of the adjacent level gap ratio or the family of Rényi entropies (of which the Shannon entropy is a particular case) [15,23,24,[35][36][37][38][39][40][41][42], change monotonically with ω, and therefore do not give rise to a neat singular point. For the previous indicators usually some form of scaling is involved in order to identify the transition point, and its value is generally largely influenced by several factors among which the most important is the number of simulated sites, L. By contrast, as can be seen in Fig. 1, the probability of extreme diagonal fluctuations allows to separate the dynamical sides transparently and is valid irrespective of L. The transition point, located at the maximum of γ 2 ( ∆ n ), grows with L, in concert with previous numerical findings. The fact that this maximum value itself also increases with size indicates that the transition is more pronounced for larger systems, but it still does not rule out the possibility that there might exist an extended transition region between the completely chaotic and integrable regimes even in the thermodynamic limit. Experiments in one-dimensional interacting bosons seem also to point to the critical point scenario by studying spatial correlations at long distances after some time evolution. However, sizes are not large enough to make a sensible extrapolation to the thermodynamic limit [92]. Thus, although our findings seem more compatible with a critical point leading to an actual phase transition, they could signify a change between two distinct, extended regimes. In this sense it is not clear whether both the semi-Poisson and the chaotic regions are finite size effects, disappearing altogether at macroscopic scales and giving rise to an abrupt change from ergodicity to MBL, or if it is a robust characteristic of disordered interacting spin chains that, like the J 1 -J 2 model, undergo a MBL transition. The answer to this question, however, lies out of the scope of this manuscript.

Conclusion
We have studied the probability of extreme events of the (normalized) fluctuations of the diagonal matrix elements of physical observables around its microcanonical equilibrium value for the J 1 -J 2 disordered quantum spin chain. For intermediate values of disorder this probability exhibits a maximum value that increases and whose exact location shifts with the system size. We interpret this result as a possible finite-size precursor of the critical point of the ergodic-MBL transition. Below this value of disorder the model is in its chaotic phase, characterized by GOE as in RMT spectra but with long-range deviations due to the Thouless energy. Beyond this value of disorder, an extended region can be identified whose spectral statistics can be described by a family of generalized semi-Poissonian statistics which show level repulsion but not chaotic correlations. Both short and long-range spectral measures can be accurately taken into account by this model. For very large values of disorder, the standard Poissonian statistics associated to the integrability of the localized phase are recovered, where both level repulsion and correlations are lost. Contrary to other ergodicity indicators such as the adjacent level gap ratio or the family of Rényi entropies, this probability is not a monotonous quantity and allows to distinguish these two regimes for any value of the number of sites unambiguously. The main conclusion of our work is that the maximum of the probability of extreme events as represented by the kurtosis excess is an indicator of the hypothetical critical point of the transition. In other words, if the ergodic and MBL phases are indeed connected by an actual phase transition and not by a smooth crossover, then the critical point must correspond to the value of disorder strength that yields the maximum of the kurtosis excess. It is interesting to note that the behavior of γ 2 ( ∆ n ) closely resembles that of a magnetic susceptibility, a robust indicator of a phase transition.
We believe our results are an important contribution for a better understanding of the ergodic-MBL transition and the MBL phase itself, opening up several new avenues of research. The behavior of the fluctuations of matrix elements from the microcanonical result has been studied before but not with our focus of computing the probability of extreme events. These ideas immediately call for a study of their relationship with Griffiths effects and their generality in other many-body systems with MBL transition. It should be mentioned that the momentum distribution is a quantity experimentally accessible in time-of-flight experiments with cold atoms [93]. By performing many copies of the same experiment with different disorder strengths, our results could be tested experimentally. One of the main questions in the field is the relationship between the standard one-body Anderson localization transition and the Anderson localized phase and the many-body localization transition and phase. Our description of the spectral statistics of the MBL phase with a generalized semi-Poisson model which is known to describe the critical behavior of the Anderson model at high dimensions should be an important motivation for the current effort of understanding the relationship between the MBL transition and the Anderson localization transition in high dimensional lattices [94]. The main problem in the study of many-body disordered systems, which is also the main limitation of the results presented in this paper, is the scaling to the thermodynamic limit, L → ∞. However, we believe that our results could serve as guidance for the search of a theoretical framework capable of overcoming this limitation.