Many-body localization and thermalization in the full probability distribution function of observables

We investigate the relation between thermalization following a quantum quench and many-body localization in quasiparticle space in terms of the long-time full distribution function of physical observables. In particular, expanding on our recent work [E. Canovi {\em et al.}, Phys. Rev. B {\bf 83}, 094431 (2011)], we focus on the long-time behavior of an integrable XXZ chain subject to an integrability-breaking perturbation. After a characterization of the breaking of integrability and the associated localization/delocalization transition using the level spacing statistics and the properties of the eigenstates, we study the effect of integrability-breaking on the asymptotic state after a quantum quench of the anisotropy parameter, looking at the behavior of the full probability distribution of the transverse and longitudinal magnetization of a subsystem. We compare the resulting distributions with those obtained in equilibrium at an effective temperature set by the initial energy. We find that, while the long time distribution functions appear to always agree {\it qualitatively} with the equilibrium ones, {\it quantitative} agreement is obtained only when integrability is fully broken and the relevant eigenstates are diffusive in quasi-particle space.


Introduction
The physics of thermalization in isolated quantum systems has a long and debated history. Recent groundbreaking experiments on the non-equilibrium dynamics of lowdimensional condensates [1,2,3] triggered a great deal of attention on this topic, which, up to them, was mostly addressed as an academic question, in connection with the notion of quantum chaos [4,5,6,7,8,9]. The demonstration of the lack of thermalization in two colliding bosonic clouds confined in a cigar-shaped potential [3], and the attribution of this observation to quantum integrability, generated a lot of theoretical activity devoted to the study of the connections between integrability, ergodicity and thermalization in strongly correlated quantum systems [10,11]. The main focus efforts has been on the characterization of thermalization resulting from the simplest possible non-equilibrium protocol: an abrupt change in time of one Hamiltonian control parameter, that is a quantum quench.
At long times after the quench, the lack of thermalization in an integrable system can be seen as a consequence of the sensitivity to the specifics of the initial state, that are encoded in the values of the constants of motion along the whole time evolution. Following the prescriptions of Jaynes [12], this qualitative understanding was made rigorous by the proposal of describing the steady state after a quench by means of a generalized Gibbs ensemble, which keeps track of the initial value of all the non trivial constants of motion [13]. The conditions of applicability and drawbacks of this approach have been extensively tested (see Ref. [10] and references therein). If the system is in turn far enough from the integrable limit, thermalization is expected in general to occur. This expectation is based on the eigenstate thermalization hypothesis, stating that expectation values of few-body observables in a given eigenstate equal thermal averages with the corresponding mean energy [4,5], and it has been tested by means of several numerical techniques (see Refs. [14,10] and references therein). The issue is still under debate [15,16,17].
A natural scenario to describe the effects of integrability and its breaking on thermalization is that of many-body localization. Building on a series of seminal papers in disordered electron systems [18,19,20], the interplay between integrability breaking and many-body localization has been recently studied in the context of thermalization [21,22,23,24], In analogy to a construction originally conceived for disordered electron systems, the quasi-particle space can be thought of as a multidimensional lattice where each point is identified by the occupations of the various quasi-particle modes. As long as states are localized in quasi-particle space, the system behaves as integrable: any initial condition spreads into few sites maintaining strong memory of the initial state. Thermalization will not occur. At the same time, qualitative behavior of local and non-local operators in the quasi-particle is naturally going to be different: locality in quasi-particle space implies sensitivity to the localization/delocalization of states, while non-local operators display always an effective asymptotic thermal behavior. Once a strong enough integrability-breaking perturbation hybridizing the various states is applied, the consequent delocalization in quasi-particle space will lead to thermalization. An initial state is allowed to diffuse into all states in a micro-canonical energy shell generating a cascade of all possible lower energy excitations.
While this scenario appears physically sound, it does not give information on what is the degree of sensitivity of the various physical quantities of a many-body quantum system to integrability breaking and thermalization. This question is particularly important in view of the fact that recent studies on the dynamics of quantum field theories lead to the proposal of a dynamics of thermalization comprising two stages [25]: first the system decays to a so-called pre-thermalized state, where the expectation value of certain macroscopic observables is to a good approximation "thermal", while the distribution function of the elementary degrees of freedom is not [25,26,27]. At a second stage, when the energy is efficiently redistributed by scattering processes, real thermalization eventually occurs. Pre-thermalization has been shown to occur theoretically for quantum quenches in a variety of systems [26,27,28,29,30,31,32]). Moreover, the study of pre-thermalization in weakly perturbed integrable systems shed light on the nature of the pre-thermalized state, which is nothing but a close relative of the non-thermal steady state attained asymptotically by its integrable counterpart [29]. Signatures of pre-thermalization have been observed in split one dimensional condensates [33], which have been shown to be characterized by an intermediate, pre-thermalized stationary state. The latter has been investigated by studying the full probability distribution of the interference contrast [33,34] which turns out to have a closely thermal behavior, even though the distribution of quasi-particles is non-thermal.
The purpose of this paper is to provide a detailed characterization of the effect of integrability breaking on the asymptotic state after a quantum quench, not only by studying the average expectation values of certain selected observables, but also focusing on their full probability distribution function (FPDF), following the suggestions that were recently put forward in Ref. [34]. We stress that the latter quantity is experimentally accessible by studying shot-to-shot variations of a physically measurable observable, as it has been performed for the FPDF of matter-wave interference in a coherently split Bose gas [33]. In the first part we introduce our model and discuss the long time limit of two kind of observables, thus summarizing the results contained in Ref. [22]; in the second part we take a considerable step forward, by extending the discussion to quantities which are closer to actual experiments [33], and clarifying more precisely to what extent thermalization takes place. We will show that, in the presence of manybody localization, the entire probability distribution function describes a canonical distribution of the degrees of freedom.
In the specific, we consider a quantum XXZ spin-1/2 chain undergoing a sudden quench of the anisotropy, in the presence of an integrability breaking term in the form of a random transverse field. As the strength of the integrability breaking term is cranked up, the many-body level statistics has a well defined transition from Poisson (Integrable) to Wigner-Dyson (non-Integrable), closely associated to the localized/diffusive character of eigenstates in quasi-particle space [22]. Focusing on the asymptotic state attained after a quench from the antiferromagnetic to the critical phase, we compute the FPDF of both transverse and longitudinal magnetization densities in a given spatial interval by means of exact diagonalization techniques. We compare the results with the FPDF of the same observables obtained in a canonical ensemble, with an effective temperature fixed by the initial energy. We show that, while for both weak and strong integrability breaking the FPDFs attained after a quench agree qualitatively with those obtained in the corresponding canonical ensemble, a full quantitative agreement for the whole distribution is only obtained whenever the level statistics in the bulk of the spectrum is of Wigner-Dyson type, corresponding to diffusive eigenstates in quasi-particle space.
This paper is organized as follows. In Sec. 2 we define the models and the details of the quantum quench protocol. In the following two sections we first discuss the spectrum of the Hamiltonian, which provides an insight in the localization/delocalization transition (Sec. 3), and then we briefly discuss how the long-time asymptotics of different correlation functions is affected by the localization/delocalization in many-body space (Sec. 4), thus summarizing the results of Ref. [22]. In Sec. 5 we take a step further: we define the FPDF of the transverse and longitudinal spin and discuss their behavior in the asymptotic long-time state attained after a quench (diagonal ensemble). Finally, in Sec. 6 we draw our conclusions.

Model
Throughout this paper we will consider a quantum quench described by a timedependent Hamiltonian: where: The Hamiltonian H(t) is composed of an integrable part H 0 [g(t)], and an integrabilitybreaking term given by H ib . Concerning the integrable part, we will consider the anisotropic spin-1/2 Heisenberg chain of length L (also called the XXZ model) with open boundary conditions: where σ α i (α = x, y, z) are the spin-1/2 Pauli matrices on site i, J is the planar xycoupling, while J z is the nearest-neighbor anisotropy parameter in the z direction, which coincides with the parameter g that will be quenched at time t = 0. In what follows we take = k B = 1, we adopt J = 1 as the energy scale and work in the zero total magnetization sector along the z axis. This model is integrable by Bethe Ansatz [35]. The zero-temperature phase diagram is characterized by three regions: a gapped ferromagnetic phase (J z < −1), a gapped antiferromagnetic phase (J z > 1), and a gapless critical phase for −1 ≤ J z ≤ 1. In this gapless region the critical exponents depend on J z and the system is characterized by a quasi-long-range order in the xy plane [36]. In the following we break the integrability of the model by applying a random magnetic field: where the quantities h i are randomly chosen in the interval [−1, 1]. We point out that the model described in Eqs. 3,4 is equivalent to a system of hardcore bosons, as one can show by applying the Jordan-Wigner transformation. Therefore it is also interesting from an experimental perspective, since cold-atom gases in one dimension and at low densities behave like impenetrable bosons [37].
After investigating the transition from integrability to non-integrability, we will address the long-time behavior of the system following a sudden quench of J z . In this context, the knowledge of a substantial part of the spectrum is required, making it necessary to resort to a standard exact diagonalization technique. We will diagonalize Hamiltonian systems with up to 14 sites, and only consider the zero magnetization sector, thus working in Hilbert spaces with up to 3432 states.
The zero-temperature phase-diagram of the XXZ model in presence of disorder is well established [38]. Much less is known at infinite temperature: the phase-diagram has been conjectured to be composed of two phases, a non-ergodic many-body localized phase (in real space) at ∆ > ∆ crit , and an ergodic one at ∆ < ∆ crit , where ∆ crit ∼ 6 ÷ 8 at J z = 1 [21,39]. Our results indicate the presence of a second non-ergodic localized phase (in quasi-particle space) for ∆ ≡ ∆ ⋆ close to zero that crosses over to the ergodic phase upon increasing ∆. The fate of this crossover in the thermodynamic limit and the eventual value of the critical ∆ ⋆ are yet to be determined ‡.
As the strength of ∆ is varied, the system deviates from integrability, as discussed more in detail in Sec. 3, where the spectrum and the properties of the eigenstates are quantitatively investigated. Indeed, as the level statistics changes from Poissonian (typical of integrable systems) to Wigner-Dyson (characterized by the level repulsion of non-integrable systems), the eigenstates are also modified. More specifically, when the system is integrable they are localized in quasi-particle space, while they become diffusive as the system becomes non-integrable. This transition affects the dynamics of the system. As summarized in Sec. 4, there is a connection between the onset of thermalization and the many-body localization transition of the eigenstates (for further details, we refer the reader to Ref. [22]). Later we will go one step forward and investigate how the localization/delocalization transition emerges in the FPDF of operators [25,26,27].

Spectral characterization of the integrability-breaking crossover
In this section we follow the approach of Ref. [22] and consider the properties of the eigenvalues and eigenvectors of the Hamiltonian in Eq. 1 for a fixed value of the anisotropy J z = 0.5, and show that the integrability breaking is associated with a localization/delocalization transition in quasi-particle space.

Statistics of the energy level spacings
In finite-size systems, it is commonly believed that the statistical distribution of the energy spacings of the quantum energy levels directly reflects the integrability properties of the model [40]. In particular an integrable quantum system is typically signaled by the presence of a Poissonian statistics in the distribution of its level spacings {s n }, s n ≡ E n+1 −E n being the spacing between two adjacent levels normalized to the average level spacing, Physically this means that the eigenvalues of the Hamiltonian within a given symmetry sector are allowed to cluster. On the contrary, for non-integrable systems level crossing is inhibited. The level statistics has a Wigner-Dyson (WD) distribution, which embeds the level repulsion in a power-law: lim s→0 P WD (s) ∼ s γ . More precisely, when antiunitary symmetry is preserved, as in the present case, the statistics is described by a Gaussian Orthogonal Ensemble [41]: By tuning the parameter ∆, the XXZ chain undergoes a transition from Poissonian to WD statistics. In a finite-size system, this transition takes the form of a smooth crossover which can be studied within a specific energy shell, by means of the following level spacing indicator (LSI) [8]: where P [E,E+W ] (s) is the level statistics computed in the window [E, E + W ], while s 0 is the first intersection point of P P (s) and P WD (s).
In Fig. 1 we show the level spacing distribution function P (s) for the levels with an excitation energy less than a given cutoff , i.e. E < E c , (see caption of Fig. 1). We see that for small ∆ the distribution is closer to an exponential, while for ∆ = 1 it almost coincides with a WD distribution, and it is shifted towards larger values of the spacings. In the inset of Fig. 1 we show η w for different intensities of the integrabilitybreaking perturbation. As the strength of the integrability-breaking term increases, the LSI approaches values close to unity for ∆ ∼ 1. At large values of ∆ ≫ J z the system tends to another integrable limit [42,43], indeed we can see that η w decreases again for ∆ 1. We note that only the states in the middle of the band display level repulsion (i.e. values of η w closer to unity), while states in extreme regions of the spectrum do not. This comes as a consequence of the two-body interaction of the system, as opposed to the behavior that is typically observed in full random matrices [44].

Inverse participation ratio vs. number of available states
Along with the Poisson-to-WD transition of the level spacing statistics, the eigenvectors also undergo a transition in their statistical properties. Indeed, when the level statistics is Poissonian and ∆ is small, the eigenstates of the Hamiltonian are close to those of the integrable XXZ chain. In other words, they are localized in quasi-particle space. On the contrary, when the statistics is WD, the eigenstates are delocalized in quasi-particle space. This picture can be made quantitative by using the inverse participation ratio (IPR) [45,46,47,48,49,50]. Given a pure state |ψ and an arbitrary basis {|n } with N elements, the IPR is defined by [41]: § (8) § With this definition we fix the typo in Ref. [22] where there was a wrong factor 1/N in the definition. We remark that the results shown in Ref. [22] are not affected by that mistake.  If a state is a superposition of n st basis states, the corresponding contribution to ξ is of order n st . Interesting results emerge if one considers two different bases: (i) the "site basis" or "computational basis" of the eigenvectors of σ z i with the constraint of zero total magnetization, and (ii) the "integrable basis" of the eigenstates of the integrable model under investigation (the XXZ chain with J z = 0.5 and ∆ = 0). Given an eigenstate of H(J z ), we found that the IPR in the integrable basis is of the order of unity if ∆ ≪ J z , and it grows with increasing ∆. Eventually, for large values of ∆, the eigenstates of H(J z ) become diffusive superpositions, with random phases and amplitudes, of the eigenstates of the integrable model. On the contrary, in the site basis, we found an opposite behavior, since the states of this basis approach the eigenstates of the system at ∆ ≫ J z [22].
We can draw an intuitive picture of the localization/delocalization transition in quasi-particle space by comparing the IPR in the integrable basis, in a given microcanonical shell, with the number of available eigenstates N [E,E+W ] in the same shell. This is shown in Fig. 2. The microcanonical shell has a width W = 2∆ ∼ V , where V is the typical matrix element of the integrablity-breaking perturbation. In quasi-integrable situations (∆ ≪ 1, left panel) the IPR is much lower than the available microcanonical states, thus meaning that the degree of delocalzation of the system is very low. On the contrary, in a chaotic situation (∆ ∼ 1, right panel), the perturbation is able to hybridize nearly all the states in the microcanonical energy shell. We point out that our approach has been also recently adopted in order to identify the emergence of chaos in a very similar spin chain model [50,51].

Energy scales involved in the quench
Let us now consider a sudden quench of the parameter g ≡ J z , from J z0 at time t ≤ 0 to J z at time t > 0, as in Eq. 2. We assume that the initial state |ψ 0 of the system is the ground state of H(J z0 ). Since after the quench the Hamiltonian is time-independent, the energy is conserved and is given by E 0 = ψ 0 |H(J z )|ψ 0 . For large values of J z0 the ground state of the Hamiltonian is the classical antiferromagnetic Néel state and the energy density E 0 /L converges to a constant value, below the middle of the spectral band of the final Hamiltonian. This means that the eigenstates of H(J z ) involved in the time evolution of the system are only those lying in the first half of the band. This is shown in Fig. 3, where the energy E 0 after the quench from an initial antiferromagnetic state (vertical lines) are compared with the number of available states in a given microcanonical energy shell. Following Rossini et al. [52,53], we then proceed to define an effective temperature associated to the quench as the solution of the equality: where ρ(T eff ) is the equilibrium density matrix at temperature T eff : Equation 9 can be solved numerically for each realization of disorder and then averaged. As shown in Ref. [22], the effective temperature increases with the quench strength |J z − J z0 |, and tends to saturate for large values of J z , since the ground state approaches the antiferromagnetic Néel state. As an example, for a system of L = 12 sites and a quench from J z0 = 10 to J z = 0.5 (that are the typical parameters we considered hereafter) the effective temperature at ∆ = 1.0 is T eff = 3.4 ± 0.4 (units of J/k B ), after averaging over 200 disorder instances.

Asymptotics of observables
The connection between the localization/delocalization transition and the onset of thermalization deeply affects the dynamics following the quantum quench. A possible way to test this interplay is to compare the long-time behavior of the system with the thermal behavior at the effective temperature T eff . This can be done at different levels.
The first possibility is to take a traditional point of view and study the correlation functions of a given observable, as we illustrate now. A second more general (to be discussed in the following) approach consists in investigating directly the full probability distribution function of selected observables. Let us start by considering the two-spin correlation functions constructed as expectation values of the following operators: The average value predicted by the canonical ensemble at temperature T eff is given by: The asymptotic value after the quench is found from the diagonal ensemble [54,14,55]: where |ψ(t) = e −iH(Jz)t |ψ o is the state of the system at time t, while c i = ψ o |φ i is the i-th component of the initial state |ψ 0 in the basis of the eigenstates {|φ i } of the final Hamiltonian H(J z ). In Fig. 4, we compare the expectation values in the diagonal ensemble (black data) of these operators with those predicted by the canonical ensemble at the temperature T eff (red data). We have chosen ∆ = 0.4, so that the system is still close to integrability and one can appreciate the discrepancies with the thermal behavior only for the correlator of σ z k (right panel). On the other hand for σ x k (left panel) a quantitative agreement with the thermal ensemble is observed. Note that the different behavior of the two observables is most visible at k = π, where the system is less sensitive to boundary effects. As discussed in Ref. [22], the different behaviors of these correlators at long times could be related to the different nature of the operators involved in the low energy limit [52,53]. Indeed, in the low energy limit [35], where the XXZ chain critical phase maps onto a Luttinger liquid and quasi-particles are approximately free bosons, σ z k turns out to be a local operator, coupling a finite number of quasi-particle states, while σ x k couples all the states and is thus non-local. To be more quantitative, we define the absolute "residual" of the operator n α k between the diagonal and the canonical ensemble: In Fig. 5 we show such quantity at k = π, as a function of ∆. The residual for the non-local operator, δn x π , does not depend substantially on the strength of ∆ nor on the size of the system. On the contrary, the residual of the local operator, δn z π , decreases significantly as the system departs from integrability and shows a minimum for∆ = 1.
For larger values of ∆ the residual δn z π starts to grow. This is consistent with the fact that for large values of ∆ the system approaches another integrable limit. Moreover, the quantity δn z π , and in particular the position of the minimum∆ are size-dependent. Given the numerical limitations of exact diagonalization, we cannot guess what is the limit of∆(L) when L → ∞. The value∆ of the perturbation strength at which the system is closer to a thermal behavior, at a given size, corresponds to the situation in which delocalization in Fock space is most pronounced.

Full probability distribution functions
The study of FPDFs of observables, instead of only their average value, has attracted a lot of interest, both theoretically [34,56,57,58] and, very recently, experimentally [33]. Residuals δn x π (black curves) and δn z π (red curves) between the diagonal and the canonical ensemble predictions. Different symbols refer to different system sizes, as depicted in the caption. Averages over 200 instances are performed.
Here we want to study the probability distribution function of the transverse and of the longitudinal spin of a subsystem S of size R ≤ L.
The xy transverse-spin magnitude of a subsystem S made by R spins is given by: where S α i = 1 2 σ α i is a spin operator on site i and direction α, while spins inside S are chosen contiguously in the central portion of the total system and are identified by the indexes j = L 2 − R 2 +1, . . . , L 2 + R 2 . This is exactly the same quantity considered in Ref. [34], where the full distribution of quantum noise in Ramsey interference experiments is studied in detail.
The z component of the spin for a subsystem of size R is given by: Since the eigenvectors of σ z R coincide with the states of the computational basis, the computational cost of evaluating the probability distribution of this quantity is less than that for (Ŝ ⊥ R ) 2 .

Asymptotics of the FPDFs
Similarly to what we have done in Sec. 4.2 for the correlation functions, we now consider the behavior of the FPDFs in the diagonal [54,14] and in the canonical ensemble at the effective temperature T eff , and compare the results. The FPDF can be evaluated in the canonical ensemble at the effective temperature T eff according to where p α T eff (s R,n ) ≡ Tr ρ(T eff ) P α R,n denotes the thermal expectation value of the projector P α R,n on a given eigenstate of σ α R (α =⊥, z) corresponding to the eigenvalue s R,n (see the correspondence with Eq. 12). On the other side, analogously to what we did in Eq. 13, we define the FPDF in the diagonal ensemble as The comparison between the two ensembles can be made quantitative by computing the absolute difference (see Eq.14): The latter quantity depends on R and on the values of s R specific to each R. We may be interested in comparing the FPDF's at different distances: in this case a useful quantity is the integrated difference, which is given by: where ′ s R,n denotes the sum over the eigenvalues of σ α R compatible with the zero total magnetization condition, and ν(R) is their number. In the following, we will average all these quantities over disorder, and identify their uncertainty with the standard deviation. We are going to show our results for a quench from J z0 = 10 to J z = 0.5 in a system with L = 12 spins.
In Fig. 6, 7 and 8 we show the probability distribution of the transverse spin (Ŝ ⊥ R ) 2 at ∆ = 0.1, ∆ = 0.5 and ∆ = 1.0, respectively. In all the cases, the distribution functions have qualitatively the same pattern in both the diagonal and the canonical ensemble. But when disorder is small, e.g. for ∆ = 0.1 in Fig. 6, there is a visible discrepancy. Increasing the intensity of the disorder a quantitative agreement is established, until for ∆ = 1, that is, where the system is furthest from integrability, the two distributions perfectly overlap (Fig. 8). The discrepancy between the two distributions is quantified by means of the integrated difference ∆p ⊥ defined in Eq. 20, which also allows to compare the distributions for different values of R. The results are shown in Fig. 9. The common feature for all values of R is that, when the disorder intensity is small, i.e. ∆ ≪ J z , the discrepancy ∆p ⊥ is larger and more sensitive to the size of R than for bigger ∆. However it is not clear from our results what the dependence on R should look like in  the thermodynamic limit. Unfortunately with our simulations we can only work with very small sizes L, such that the smallest R is only one order of magnitude less that L.
In Fig. 10 we show the FPDF of the longitudinal spin σ z R in the diagonal and the canonical ensemble for a system of L = 12 sites and R = 6. The constraint j σ z j = 0 allows only some eigenvalues of σ z R , so when the size of the subsystem increases from R = L/2 to R = L progressively less values of s R are permitted. For this reason we choose R = 6, which is the case with the maximum number of eigenvalues allowed by the symmetries. We observe in all the cases that there is a maximum at s R = 0, which is the most sensitive point to the variation of ∆. For odd R (data not shown) there are instead two symmetrical eigenvalues, s R = ±1, where the function has two maxima and which are most suitable to check the dependence against ∆. From Fig. 10 we see that at small values of ∆ there is a difference at s R = 0 between the two ensembles, while the FPDFs almost coincide for ∆ ∼ 1. The scenario is similar to that of (S ⊥ R ) 2 . In Fig. 11 we show the integrated difference ∆p z . We see a slight dependence on whether R is even or odd. Namely, for even values of R the discrepancy between the ensembles at small values of ∆ is more pronounced than for odd R. This effect is more visible for small values of ∆. For even values of R the difference ∆p z has its smallest value for ∆ = 1, while for odd values of R the minimum is at ∆ ∼ 0.5. Looking at the error bars, we see that both p z and ∆p z are generally characterized by larger fluctuations with respect to the transverse spin, perhaps because the spectrum of σ z R has a higher degeneracy with respect to (S ⊥ R ) 2 . However the error bars, both for the transverse and longitudinal spin, are proportional to the disorder intensity. As we already pointed out for the transverse spin, we cannot fully characterize the dependence of ∆p z on R, because we have too small systems. Intuitively, what we can expect is that the constraint on the total magnetization j σ z j = 0 should relate the situations in which the subsystem has R and L − R sites respectively. On a qualitative level this is also suggested by the panels of Fig. 11.

Conclusions
In this paper we analyzed the connection between thermalization and many-body localization. We studied the long time behavior of a XXZ spin-1/2 chain following a quantum quench, when integrability is broken by a random magnetic field. In particular we addressed the issue of pre-thermalization, not only looking at the behavior of average values of observables, but also at their FPDF. We found a qualitative agreement between the FPDF in the asymptotic state and that predicted by a thermal distribution in the whole range of disorder intensity we considered. Nevertheless, the situation in which the two FPDFs are quantitatively indistinguishable occurs only when the system is furthest from integrability, that is when the eigenstates are diffusive superpositions in quasiparticle space.