# Extended nonergodic states in disordered many-body quantum systems

## Abstract

This work supports the existence of extended nonergodic states in the intermediate region between the chaotic (thermal) and the many-body localized phases. These states are identified through an extensive analysis of static and dynamical properties of a finite one-dimensional system with onsite random disorder. The long-time dynamics is particularly sensitive to changes in the spectrum and in the structures of the eigenstates. The study of the evolution of the survival probability, Shannon information entropy, and von Neumann entanglement entropy enables the distinction between the chaotic and the intermediate region.

## 1 Introduction

The Anderson localization in real space of a particle in a disordered medium (1) is theoretically well understood (2) and experimentally confirmed in various settings (3); (4); (5); (6); (7). In systems with many interacting particles, despite the consensus that many-body localization should still take place in the limit of strong disorder, the details of the metal-insulator transition are richer than in noninteracting systems and not yet completely understood.

The seeds for the analysis of many-body localization were sown in the 80’s (8); (9); (10) and 90’s (11); (12); (13); (14); (15). The general expectation was that interaction might hamper localization, but not completely prevent it. The study of two interacting particles in a one-dimensional (1D) disordered chain, for instance, showed that the localization length was larger than in a one-particle system (12). A number of theoretical studies followed in the 2000’s (16); (17); (18); (19); (20); (21); (22); (23); (24); (25); (26), but the current decade has witnessed a burst of interest in the subject, with not only theoretical (27); (28); (29); (30); (31); (32); (33); (34); (35); (36); (37); (38); (39); (40); (41); (42); (43); (44); (45); (46); (47); (48); (49); (50); (51); (52); (53); (54); (55); (56); (57); (58); (59); (60); (61); (62); (63); (64); (65), but also experimental works (66); (67); (68); (69).

Among the several questions that have been raised, the existence or not of an intermediate phase between the chaotic and the many-body localized phase (14); (36); (55); (56); (70); (58); (71); (72) is the one that mainly motivates the present work. Our numerical results for a finite 1D spin-1/2 system with nearest-neighbor interaction and onsite random magnetic fields suggest a positive answer. A key aspect of our approach is the shift of the emphasis from spectral properties to dynamical properties of the system. This establishes a strong connection between our studies and current experiments where dynamics is typically investigated. We show that detailed information about the spectrum, eigenstates, and initial state can be obtained from analyzing the system’s time evolution.

As the strength of the disorder increases from zero, where the system is integrable, the following regions are covered, before localization is finally reached: (i) transition from the integrable to the chaotic domain, (ii) chaotic regime, and (iii) intermediate region between the chaotic limit and the many-body localized phase. Chaotic states prevail in the chaotic domain, while in the intermediate region before localization, the eigenstates are delocalized (extended) in the configuration space but nonchaotic (nonergodic).

Our definition of a chaotic state is as follows. Given an eigenstate , the inverse participation ratio

(1) |

measures how much delocalized the eigenstate is in the basis vectors , with being the dimension of the Hamiltonian matrix. The eigenstate is chaotic if it samples most of the Hilbert space without any preference, which gives . The paradigmatic example of chaotic states are the eigenstates of full random matrices (FRM), which are (pseudo-)random vectors and therefore maximally delocalized. When these matrices are the real and symmetric ones from Gaussian orthogonal ensembles (GOE) (73); (74); (15), the components of their eigenstates are real uncorrelated random numbers from a Gaussian distribution satisfying the normalization condition, so (75); (76).

In FRM, the notion of basis is not well defined, but in realistic systems, they certainly affect the results for . In this case, the choice of basis is done according to the question we are addressing. In studies of spatial localization, the focus is on the basis that represents the configuration space. Another important aspect is the energy dependence of . In realistic finite chaotic systems with two-body interactions, since the density of states is Gaussian (77), it is away from the edges of the spectrum that we expect to find eigenstates that are close to random vectors and have , while near the borders they tend to be more localized. Strict ergodicity in the sense of full random matrices does not hold in these systems.

A delocalized but nonchaotic eigenstate has with . For a localized eigenstate, . Our focus is on the middle of the spectrum of regions (ii) and (iii), where we find respectively:

Chaotic states: | (2) | ||||

Delocalized nonchaotic states: | (3) |

The three different regions between the integrable point and the localized phase are distinguished through the analysis of static and dynamical properties. We do not aim at locating exactly the point that separates one region from the other, but discuss the best quantities to identify and characterize each one. System sizes larger than the ones treated here are necessary to delineate the borders and precisely determine how they depend on parameters of the spin-1/2 systems, such as anisotropy, and on the energy of the states.

Our study of the statistic properties involves level statistics, structure of the eigenstates and structure of the initial states. They are used for comparisons with our results for the dynamics. The latter is the actual center of our attention.

We analyze the spectrum using nearest-neighbor level spacing distributions and level number variance. The structures of the eigenstates are investigated with the inverse participation ratio, the Shannon information entropy, and the von Neumann entanglement entropy. The distribution of the components of the initial state written in the energy eigenbasis are contrasted with the Porter-Thomas distribution (77). In regions (i) and (iii), the level statistics is intermediate (78); (18); (19) between that of uncorrelated and repelling eigenvalues. The results indicate that region (i) may disappear in the thermodynamic limit. In region (iii), the eigenstates are nonchaotic but delocalized, and there is no agreement with the Porter-Thomas distribution. This region appears to shrink with system size and may turn into a critical point in the thermodynamic limit.

For the dynamics, we investigate the time evolution of both the Shannon and entanglement entropies and of the survival probability. The latter is defined as

(4) |

where is the initial state, is the Hamiltonian describing the system, and . The decay of the survival probability at long times is necessarily power law, , even in the chaotic regime (79); (80). The value of the power-law exponent is very sensitive to the properties of both the initial state and the Hamiltonian. In region (ii), where there is level repulsion following the Wigner-Dyson distribution and the eigenstates are chaotic, we have that . In region (iii), and it coincides with from Eq. (3)(42); (43).

The presence of level repulsion is manifested also in the so-called correlation hole (81); (82); (83); (84). This corresponds to a time interval where drops below its infinite time average. It happens before saturation takes place and is very evident in region (ii). In region (iii), the correlation hole becomes less deep and fades aways at the approach of the localized phase.

The dynamical behavior of both entropies is very similar. In the chaotic region, they grow linearly in time and quickly saturate, while in region (iii), the linear increase is followed by a logarithmic behavior, before saturation. Since the results for both entropies are comparable, either one can be used in the analysis. The advantage of the Shannon entropy is to be computationally less expensive than the entanglement entropy, because contrary to the latter, it does not require performing a partial trace over degrees of freedom.

The main characteristics of regions (i), (ii), and (iii) are summarized in Fig. 1.

## 2 Model and Static Properties

The 1D disordered spin-1/2 system that we consider has two-body nearest-neighbor couplings, sites, and periodic boundary conditions. The Hamiltonian is given by

(5) |

where are the spin operators of each site . We set . The disorder corresponds to random static magnetic fields of amplitude , where are random numbers from a uniform distribution and is the disorder strength. The total spin in the -direction, , is conserved. We study the largest subspace, , which has dimension .

### 2.1 Level statistics

The observation of intermediate level statistics in the vicinity of the metal-insulator transition dates back to (85); (86). In the particular case of many-body quantum systems described by 1D disordered spin-1/2 models, level statistics was studied in (78); (18); (19). More recent works include the detailed analyses developed in (50); (60). Here, the brief discussion of spectral statistics is used for comparison with our results for delocalization measures and dynamics. In addition to the frequently investigated nearest-neighbor level spacing distribution, we consider also the level number variance.

Level repulsion is a main signature of quantum chaos. After unfolding the spectrum (15); (87) and separating the rescaled eigenvalues by symmetry sectors (87); (88), level repulsion can be detected, for example, by computing the nearest-neighbor level spacing distribution , where is the spacing between two neighboring levels. For FRM from GOE or any chaotic real and symmetric Hamiltonian matrix, is the Wigner-Dyson distribution of shape . In systems where the levels are not prohibited from crossing and where the number of degeneracies are not excessive (89), the typical distribution is Poissonian, . To quantify the proximity to the Wigner-Dyson distribution, we use the chaos indicator (90)

(6) |

where is the first intersection point of and . For a Poisson distribution, , and for the Wigner-Dyson, .

In Fig. 2 (a), we show averaged over several disorder realizations. The average is denoted with . As increases from zero, the level spacing distribution for (5) first transitions abruptly from Poisson () to Wigner-Dyson (small ) in region (i). It is Wigner-Dyson in region (ii), where plateaus to a small value. It then transitions from Wigner-Dyson back to Poisson in region (iii). The boundary of region (iii) after which localization emerges occurs approximately where the curves for different system sizes cross, seen in the figure for . Beyond this point, the larger the system is the closer the distribution is to .

The level spacing distribution provides information about the short-range fluctuations of the spectrum. Information about long-range correlations, which better captures how rigid the spectrum is, can be obtained with quantities such as the level number variance , where is the number of unfolded energy levels in the interval and the bar denotes here the average over different initial values of (15). For FRM from GOE, , where is the Euler’s constant. Uncorrelated levels lead to . A careful study of the level number variance in disordered interacting systems and how it connects with the existence of the Thouless energy was done in (60).

In Fig. 2 (b), we show how varies with the disorder strength. The results are very similar to those in Fig. 2 (a). To draw attention to region (i), we show in a logarithmic scale. In this region, the disorder strength for which and reach the smallest values decreases as the system size increases. This suggests that in the thermodynamic limit, region (i) may disappear and an infinitesimally small should take the system into the chaotic domain (91); (92). The chaotic region (ii) stretches with not only in the direction of smaller ’s, but also for larger disorder strengths. Whether region (iii) will also disappear or persist in the thermodynamic limit, possibly as a critical point, is an open question.

We stress that the values of and especially of are very sensitive to the unfolding procedure used (60); (93). The results suffer also from finite size effects and, in the case of , from additional symmetries. Hence, the purpose of our illustrations is to support the existence of different regions associated with different values of , but not to exactly delineate their boundaries.

### 2.2 Delocalization measures

The eigenstates provide much information about the system. One way to study their structures is with the moments . One finds, for example, that at the Anderson transition critical point, the scaling analysis of with the dimension (that is, ) leads to a nonlinear dependence of the generalized fractal dimension on . This indicates that the eigenstates become multifractal (94). Here, we focus on and drop the subscript . This corresponds to the inverse participation ratio defined in Eq. (1).

Fractal dimensions can also be obtained from , where is the Rényi entropy (95). When , the Rényi entropy coincides with the Shannon information entropy,

(7) |

Like the inverse participation ratio, is extensively used to measure the level of delocalization of quantum states (75). In this work, we compute from an from .

#### Inverse Participation Ratio

Figure 3 (a) shows the ratio , where . We average over 10% of the eigenstates that have the closest energies to the middle of the spectrum, where they are more delocalized, and over several disorder realizations, so that the total amount of data is . Table 1 gives the specific values for each system size.

System | Number of States | Number of |

Size | from the | Disorder |

Middle of the Spectrum | Realizations | |

8 | 8 | 14 285 |

10 | 26 | 3968 |

12 | 93 | 1082 |

14 | 344 | 291 |

16 | 1288 | 77 |

The value of depends, of course, on the basis selected. To study localization in real space, the natural basis corresponds to the states where on each site the spin either points up or down in the direction, as for example . We refer to these states as site-basis vectors; in quantum information theory they are called computational basis vectors.

The behavior of from region (i) to (iii) is nonmonotonic with . This is better seen in a linear plot (26); (96) than in the log-plot of Fig. 3 (a). The highest level of delocalization occurs in the chaotic region (ii). Notice also that for all ’s, as expected for sparse banded Hamiltonian matrices (76).

The curves of the ratio for different ’s cross at a value of . The fact that after this point the ratio decreases with implies that (equivalently, ) with . The eigenstates are no longer chaotic, although they are still delocalized. This is the beginning of region (iii), while the crossings in Fig. 2 mark its end.

The values of the fractal dimension as a function of the disorder strength are shown in Fig. 3 (b). is the slope of the plot of vs . It is obtained from a scaling analysis with . Some illustrations can be found in Ref. (42). The values of decrease significantly in region (iii), while in region (i) they are not far from 1.

Close to the integrable point, the eigenstates remain delocalized in configuration space, although for the system sizes considered . A question that has been debated due to its importance for the subject of thermalization (97); (98); (99); (100); (101) is whether for the eigenstates right in the middle of the spectrum are chaotic or not.

#### Shannon and Entanglement Entropies

With Figs. 3 (c) and (e), we compare the Shannon information entropy with the von Neumann entanglement entropy, . The latter is vastly employed in quantum information theory to quantify the amount of entanglement in a state (102); (103). It is obtained by dividing the system in subsystems A and B and tracing out the degrees of freedom of one of the two, that is,

(8) |

where and is the density matrix of the whole system. The dimension of the Hilbert space of the remaining chain of length is . We consider .

The values of the entropies are maximum for the random vectors of FRM. In this case, (75); (76). For the entanglement entropy, it was shown that when and (104). In analogy with the Shannon entropy, the alternative approximation was considered in (76) for and is used here.

Figures 3 (c) and (e) display, respectively, the ratios and as a function of the disorder strength for different system sizes. The results for both entropies are very similar and the curves intersect approximately at the same point. This crossing point has been used to detect the transition to localization (29); (105). It occurs at a value of similar to the one in Fig. 2 (a). Since the results for both entropies are similar, either one can be used in the studies of many-body localization. The computational advantage of the Shannon entropy is that it does not involve the trace operation.

Figures 3 (d) and (f) give and as a function of the disorder strength. They are, respectively, the slopes of the plots of vs and vs . The scaling analysis is done with . The behaviors of and are comparable. In the chaotic region (ii), they are 1. Values slightly above 1 are sometimes seen, but they must be caused by finite size effects, the entanglement entropy being more sensitive to them. Both and become smaller than 1 for close to the point where is also , that is at the beginning of region (iii).

We notice that in region (iii), the plots for vs and vs are still linear, yielding the volume-law scaling of the entropies with system size, but since the slopes are smaller than 1, the eigenstates are no longer chaotic. It is in this region that we expect to find the “two-component” structure of the entanglement spectrum described in (44).

## 3 Dynamical Properties

To study the dynamics, we consider as initial state a site-basis vector, . Our analysis is performed for the 10% site-basis vectors that have energy closest to the middle of the spectrum. The averages are done over these states and over several disorder realizations (see Table 1). The total amount of data for each is approximately .

We investigate the behavior of the survival probability given by Eq. (4). It can also be written as

(9) |

where is the energy distribution of weighted by the components . This distribution is referred to as the local density of states (LDOS). The survival probability is the absolute square of the Fourier transform of the LDOS.

The perturbation that takes the system out of equilibrium is very strong. The Hamiltonian initially has , so , and the coupling strength is abruptly changed to the finite value . Since the strength of the perturbation is strong, the envelope of the LDOS is Gaussian. This mirrors the shape of the density of states of systems with two-body interactions (see Refs. (106); (107); (108) and references therein for more details). The LDOS may, however, be very sparse. This is what happens in region (iii), where the disorder is strong (42) and the energy eigenbasis are no longer chaotic.

At very short times, the decay of the survival probability is quadratic, , where is the width of the LDOS. Subsequently, for intermediate times, the behavior depends on the shape of the LDOS. The decay remains Gaussian when the envelope of the LDOS is also Gaussian or it switches into an exponential decay when is Lorentzian (106); (107); (108). For long times, however, it is the filling of the energy distribution of the initial state that determines the behavior of (79); (80). Independently of how fast the initial evolution may be, the long-time decay necessarily slows down and becomes power law (79); (80), . The long-time evolution is the focus of Secs. 3.2 and 3.3.

We also study the evolution of the entanglement entropy and how the Shannon entropy spreads in time over the other site-basis vectors, that is,

(10) |

where is the probability for the initial state to be found in state at time . For both entropies, the evolution at very short times is (107); (76). Our interest is in what happens afterwards.

### 3.1 Infinite time averages and Porter-Thomas distribution

We start by analyzing how the structure of the initial state written in the energy eigenbasis depends on . The inverse of the participation ratio of the initial state, , coincides with the survival probability after its saturation, that is after it reaches its infinite time average,

(11) | |||||

The results for are shown in Fig. 4 (a). The scaling analysis gives , depicted in Fig. 4 (b). There is great similarity between Figs. 4 (a), (b) and Figs. 3 (a), (b). In fact, .

Initial states for which are chaotic, which implies the uniform and unbiased sampling of the Hilbert space. In a chaotic initial state, are approximately zero-centered Gaussian random variables. For a large system, the weights follow the Porter-Thomas distribution (109); (77); (75),

(12) |

Notice that the weights appear also in the denominator of the first term. The Porter-Thomas distribution is a result from random matrix theory. It is an additional tool to determine whether the (initial) states are ergodic or not. As the states localize, the components fluctuate more and deviate from the Porter-Thomas. Disagreement with this distribution occurs whenever .

In Figs. 4 (c) and (d), we show the distribution of the weights for two different disorder strengths. The logarithmic of is used to better visualize the distribution of the small components, which appear in large amounts. For , the system is in the chaotic regime () and consequently the distribution of the components follows very well the Porter-Thomas distribution. As the disorder strength increases and the system approaches the localized phase, the distribution of the components broadens. In Fig. 4 (d) we depict the case of , where ; the deviation from the Porter-Thomas distribution is evident.

### 3.2 Time evolution under strong disorder

As mentioned above, the long-time decay of the survival probability becomes power law, so the average over initial states and disorder realizations gives . There are different causes for this algebraic behavior. In region (iii), where the LDOS is sparse, and the eigenstates are correlated. In this case, the power-law exponent reflects the level of correlations of the eigenstates and of the components . We have that (42); (43). This is similar to what one finds at the metal-insulator transition of noninteracting models (110); (111).

Some illustrations of the behavior of survival probability in region (iii) are provided in Figs. 5 (a), (b), and (c) (more examples can be found in Ref. (42)). The initial decay of is Gaussian up to and then it slows down. For the disorder strengths of the three panels, . One sees that for , the numerical results for (solid lines) agree very well the power-law decay (dashed lines).

For both entropies, after the quadratic behavior for very short times (107); (76), the evolution becomes linear in approximately up to where the decay of is still Gaussian. The linear increase of the Shannon entropy was studied in (76); (98); (99); (112); (113), where semi-analytical expressions were also provided. In region (iii), the linear increase is followed by a logarithmic behavior in time. This occurs for the entanglement entropy [see Figs. 5 (g), (h), (i)], as has been discussed before (24); (28), but also for the Shannon entropy, as seen in Figs. 5 (d), (e), (f). The analogous behavior reinforces our statement that either one of the entropies can be used in these studies.

Motivated by the power-law decay of the survival probability, , and by the fact that, from Eq. (10), , we fit the logarithmic growth of with

(13) |

Our results indicate a very good agreement between and . This comparison is done in Fig. 4 (b) for . We should expect improvement for larger system sizes. In fact, in Figs. 5 (d), (e), (f), the dashed lines correspond to , with only being fitted. Similarly to the power-law decay of with , the onset of the logarithmic behavior for must also be associated with the emergence of correlated eigenstates. This connection deserves further investigation.

### 3.3 Time evolution under weak disorder

All the results for region (ii) indicate chaos. The level spacing distribution has the Wigner-Dyson shape, and are , and the distribution of the components of the initial state agrees very well with the Porter-Thomas distribution. The quantities that measure the level of chaoticity saturate: reaches its smallest values and . Yet, the power-law exponent still varies with . Throughout region (ii), , so it no longer coincides with . As decreases below 1, increases monotonically above 1 up to its maximum value (79); (80).

The value is a consequence of the Khalfin effect, which refers to the onset of the power-law decay of the survival probability caused by the unavoidable presence of energy bounds in the spectrum (114). The phenomenon is well understood in the case of continuous models (115); (116); (117). We have extended these studies to finite many-body quantum systems, where the spectrum is discrete (79); (80). In this case, when the initial state is chaotic and the LDOS is very well filled, the same sort of analysis used for continuous remains valid. The Fourier transform of a Gaussian LDOS that takes into account the bounds in the spectrum leads asymptotically to . This is the behavior seen in Fig. 6 (a). This figure is obtained for , where approaches the largest value for . As increases above up to , decreases from and approaches . This progressive reduction of is shown in Figs. 6 (b), (c), and (d).

We expect a chaotic initial state, which implies , to thermalize. A necessary condition for thermalization (that is, for the coincidence between the infinite-time averages of few-body observables and their thermodynamic averages) is the presence of chaotic states, which leads to the unbiased sampling of the Hilbert space (118); (119). As we saw above, when , the power-law exponent ranges between 1 and 2, while for , the exponent is smaller than . This indicates that from the value of we can anticipate whether or not the initial state will thermalize. This point was carefully studied in Refs. (79); (80).

Combining the results from Sec. 3.2 and Sec. 3.3, one sees that by varying the strength of the disorder in (5), all values of that are reachable by realistic models with two-body interactions, that is , can be covered. These same values can be achieved also with banded random matrices (79); (80). Examples of these matrices include the power-law banded random matrices that have been widely used in studies of the Anderson metal-insulator transition (94) and the Wigner banded random matrices (120). The variances of their random numbers are large within a bandwidth around the diagonal elements and very small or zero outside. As increases from one, all values of from 0 to 2 are found. In addition, with a very broad we can go beyond these realistic values of . The maximum power-law exponent is reached in the limit of full random matrices when (79); (80); (76). Banded random matrices provide a general picture of the available power-law exponents for the quench dynamics of isolated many-body quantum systems without being restricted to a particular model. This emphasizes that the results presented here are general and not specific to the spin-1/2 model considered.

In contrast to the survival probability, the entropy growth in time does not provide a clear signature of the Khalfin effect for the system sizes examined here. After the linear increase, where one might expect it to show up, the Shannon entropy and the entanglement entropy simply saturate, as shown for the Shannon entropy in Figs. 6 (e) and (f). The survival probability is more sensitive to changes in the spectrum and in the structures of the eigenstates than the entropies.

### 3.4 Correlation hole

The results above show that details about the Hamiltonian and the initial state can be extracted from the dynamics. This is very useful for experiments that routinely study dynamics. The evolution of the survival probability, in particular, offers direct information about level statistics. The presence of level repulsion can be inferred from the difference between the minimum value of the survival probability, , and its infinite time average, . This difference is known as correlation hole (81); (82); (83); (84). It does not exist in systems with uncorrelated eigenvalues, while in FRM the hole is very deep. The notion of correlation hole was introduced as a method to obtain information about level statistics from systems where one has only partial access to the spectrum (81); (82); (83).

In Figs. 7 (a), (b) and (c), we show for long times, up to saturation. The disorder strength increases from (a) to (c). The dot-dashed line at the bottom of each panel indicates and the dashed line marks the infinite time average . In Fig. 7 (a), where the system is still chaotic (), the correlation hole is deep. As increases, the hole fades away and practically disappears in Fig. 7 (c), where .

To measure the depth of the correlation hole, we compute

(14) |

For FRM of GOE, and (84), which leads to .

In Fig. 7 (d), we show the dependence of on the disorder strength. It peaks in the chaotic region, where it approaches . It decreases in both directions, that of small , where the system gets closer to integrability, and that of large , where the system approaches localization.

## 4 Conclusions

From the analysis of static and dynamical properties of a finite 1D spin-1/2 system with onsite random disorder, we identified three regions between the integrable limit, where the disorder strength is zero, and the many-body localized phase. They are: (i) the transition region between the integrable and the chaotic domain, (ii) the chaotic region, and (iii) the intermediate region between chaos and localization. We argued that the eigenstates in region (iii) are delocalized (extended) but nonchaotic (nonergodic).

Special attention was given to the dynamics, which makes a strong connection with experiments that routinely study dynamics. Information about the spectrum, eigenstates, and initial state can be obtained from the evolution of different quantities. The survival probability, in particular, reveals details that the Shannon information entropy and the von Neumann entanglement entropy cannot capture for the system sizes considered. These include the correlation hole, which is a way to directly detect level repulsion from dynamics, and the Khalfin effect caused by the unavoidable energy bounds of the spectrum.

We also showed that the results for both entropies are comparable, thus either one can be used for the studies of many-body localization. The calculation of the Shannon entropy is more straightforward, because it does not involve any partial trace.

The main characteristics of regions (i), (ii), and (iii) are summarized below.

Region (i) has intermediate level statistics (level spacing distribution between Poisson and Wigner-Dyson). This region may disappear in the thermodynamic limit.

Region (ii) shows level repulsion. The eigenstates away from the borders of the spectrum are chaotic, that is . The components of the initial states satisfy the Porter-Thomas distribution and . The power-law decay of the survival probability at long times has exponent , where indicates ergodic filling of the energy distribution of the initial state. The correlation hole is deep and close to the limits established by random matrix theory. The linear growth of the Shannon entropy and of the von Neumann entanglement entropy is followed by saturation. This is the region where thermalization is expected.

Region (iii) has intermediate level statistics. It seems to become smaller as the system size increases, so it may reduce to a critical point in the thermodynamic limit. In this region, the eigenstates and initial states are delocalized, but nonchaotic, that is , with . From the analysis of the entropies, we also find that . Fractal dimensions smaller than indicate that the states are fractal and therefore nonergodic. This lack of ergodicity is reflected into the dynamics. The power-law exponent of the long-time decay of is smaller than and agrees with the fractal dimension, . The linear increase in time of the entropies is followed by a logarithmic growth before saturation. Both and the logarithmic behavior are caused by correlations in the eigenstates. This is the region where one expects subdiffusive transport.

## Acknowledgments

EJTH acknowledges funding from CONACyT, PRODEP-SEP and Proyectos VIEP-BUAP 2016, Mexico. EJTH is also grateful to LNS-BUAP for allowing use of their supercomputing facility. This work was supported by the NSF grants No. DMR-1147430 and No. DMR-1603418.

### References

- P. W. Anderson Phys. Rev. 109, 1492 (1958).
- P. A. Lee and T. V. Ramakhrishnan Rev. Mod. Phys. 57, 287 (1985).
- B. Kramer and A. MacKinnon Rep. Prog. Phys. 56, 1469 (1993).
- J. Billy, V. Josse, Z. Zuo, A. Bernard, B. Hambrecht, P. Lugan, D. Clement, L. Sanchez-Palencia, P. Bouyer, and A. Aspect Nature 453, 891 (2008).
- G. Roati, C. D’Errico, L. Fallani, M. Fattori, C. Fort, M. Zaccanti, G. Modugno, M. Modugno, and M. Inguscio Nature 453, 895 (2008).
- G. Modugno Rep. Progr. Phys. 73, 102401 (2010).
- W. R. McGehee, S. S. Kondov, W. Xu, J. J. Zirbel, and B. DeMarco Phys. Rev. Lett. 111, 145303 (2013).
- L. Fleishman and P. W. Anderson Phys. Rev. B 21, 2366 (1980).
- B. L. Altshuler, I. K. Zharakeshev, S. A. Kotochigova, and B. I. Shklovskii Sov. Phys. JETP 67, 625 (1988).
- T. Giamarchi and H. J. Schulz Phys. Rev. B 37, 325(1988).
- O. N. Dorokhov Sov. Phys. JETP 71, 360 (1990).
- D. L. Shepelyansky Phys. Rev. Lett. 73, 2607 (1994).
- Y. Imry Europhys. Lett. 30, 405 (1995).
- B. L. Altshuler, Y. Gefen, A. Kamenev, and L. S. Levitov Phys. Rev. Lett. 78, 2803(1997).
- T. Guhr, A. Mueller-Gröeling, and H. A. Weidenmüller Phys. Rep. 299, 189 (1998).
- B. Georgeot and D. L. Shepelyansky Phys. Rev. E 62, 6366 (2000).
- G. P. Berman, F. Borgonovi, F. M. Izrailev, and V. I. Tsifrinovich Phys. Rev. E. 64, 056226 (2001).
- L. F. Santos J. Phys. A 37, 4723 (2004).
- L. F. Santos, G. Rigolin, and C. O. Escobar Phys. Rev. A 69, 042304 (2004).
- L. F. Santos, M. I. Dykman, M. Shapiro, and F. M. Izrailev Phys. Rev. A 71, 012317 (2005).
- I. V. Gornyi, A. D. Mirlin, and D. G. Polyakov Phys. Rev. Lett. 95, 206603 (2005).
- D. M. Basko, I. L. Aleiner, and B. L. Altshuler Ann. Phys. 321, 1126 (2006).
- V. Oganesyan and D. A. Huse Phys. Rev. B 75, 155111 (2007).
- M. Žnidarič, T. Prosen, and P. Prelovšek Phys. Rev. B 77, 064426 (2008).
- W. G. Brown, L. F. Santos, D. Starling, and L. Viola Phys. Rev. E 77, 021106 (2008).
- F. Dukesz, M. Zilbergerts, and L. F. Santos New J. Phys. 11, 043026 (2009).
- A. Pal and D. A. Huse Phys. Rev. B 82, 174411 (2010).
- J. H. Bardarson, F. Pollmann, and J. E. Moore Phys. Rev. Lett. 109, 017202 (2012).
- A. D. Luca and A. Scardicchio Europhys. Lett. 101, 37003 (2013).
- J. Z. Imbrie J. Stat. Phys. 163, 998 (2016).
- J. Z. Imbrie Phys. Rev. Lett. 117, 027201 (2016).
- D. A. Huse, R. Nandkishore, and V. Oganesyan Phys. Rev. B 90, 174202 (2014).
- J. A. Kjäll, J. H. Bardarson, and F. Pollmann Phys. Rev. Lett. 113, 107204 (2014).
- Y. Bar Lev and D. R. Reichman Phys. Rev. B 89, 220201 (2014).
- Y. Bar Lev, G. Cohen, and D. R. Reichman Phys. Rev. Lett. 114, 100601 (2015).
- T. Grover, Certain general constraints on the many-body localization transition, arXiv:1405.1471.
- R. Vosk, D. A. Huse, and E. Altman Phys. Rev. X 5, 031032 (2015).
- A. Chandran and C. R. Laumann Phys. Rev. B 92, 024301 (2015).
- D. J. Luitz, N. Laflorencie, and F. Alet Phys. Rev. B 91, 081103 (2015).
- D. J. Luitz, N. Laflorencie, and F. Alet Phys. Rev. B 93, 060201 (2016).
- V. Ros, M. Müller, and A. Scardicchio Nucl. Phys. B 891, 420 (2015).
- E. J. Torres-Herrera and L. F. Santos Phys. Rev. B 92, 014208 (2015).
- E. J. Torres-Herrera, M. Távora, and L. F. Santos Braz. J. Phys. 46, 239 (2016).
- Z. C. Yang, C. Chamon, A. Hamma, and E. R. Mucciolo Phys. Rev. Lett. 115, 267206 (2015).
- E. Altman and R. Vosk Ann. Rev. Cond. Mat. Phys. 6, 383 (2015).
- R. Nandkishore and D. Huse Annu. Rev. Condens. Matter Phys. 6, 15 (2015).
- M. Serbyn, Z. Papić, and D. A. Abanin Phys. Rev. Lett. 111, 127201 (2013).
- M. Serbyn, Z. Papić, and D. A. Abanin Phys. Rev. B 90, 174302 (2014).
- M. Serbyn, Z. Papić, and D. A. Abanin Phys. Rev. X 5, 041047 (2015).
- M. Serbyn and J. E. Moore Phys. Rev. B 93, 041424 (2016).
- R. Singh, J. H. Bardarson, and F. Pollmann New J. Phys. 18, 023046 (2016).
- K. Agarwal, S. Gopalakrishnan, M. Knap, M. Müller, and E. Demler Phys. Rev. Lett. 114, 160401 (2015).
- S. Gopalakrishnan, K. Agarwal, E. A. Demler, D. A. Huse, and M. Knap Phys. Rev. B 93, 134206 (2016).
- C. Monthus Entropy 18, 122 (2016).
- X. Li, S. Ganeshan, J. H. Pixley, and S. Das Sarma Phys. Rev. Lett. 115, 186601 (2015).
- X. Li, J. H. Pixley, D. L. Deng, S. Ganeshan, and S. Das Sarma Phys. Rev. B 93, 184204 (2016).
- T. Devakul and R. R. P. Singh Phys. Rev. Lett. 115, 187201 (2015).
- J. Goold, C. Gogolin, S. R. Clark, J. Eisert, A. Scardicchio, and A. Silva Phys. Rev. B 92, 180202 (2015).
- L. F. Santos, F. Borgonovi, and G. L. Celardo Phys. Rev. Lett. 116, 250402 (2016).
- C. L. Bertrand and A. M. García-García Phys. Rev. B 94, 144201 (2016).
- C. Gogolin and J. Eisert Rep. Prog. Phys. 79, 056001 (2016).
- D. Cohen, V. I. Yukalov, and K. Ziegler Phys. Rev. A 93, 042101 (2016).
- O. S. Bari šić, J. Kokalj, I. Balog, and P. Prelovšek Phys. Rev. B 94, 045126 (2016).
- V. Khemani, F. Pollmann, and S. L. Sondhi Phys. Rev. Lett. 116, 247204 (2016).
- T. Enss, F. Andraschko, and J. Sirker, Many-body localization in infinite chains, arXiv:1608.05733.
- M. Schreiber, S. S. Hodgman, P. Bordia, H. P. Lüschen, M. H. Fischer, R. Vosk, E. Altman, U. Schneider, and I. Bloch Science 349, 842 (2015).
- S. S. Kondov, W. R. McGehee, W. Xu, and B. DeMarco Phys. Rev. Lett. 114, 083002 (2015).
- P. Bordia, H. P. Lüschen, S. S. Hodgman, M. Schreiber, I. Bloch, and U. Schneider Phys. Rev. Lett. 116, 140401 (2016).
- J. Smith, A. Lee, P. Richerme, B. Neyenhuis, P. W. Hess, P. Hauke, M. Heyl, D. A. Huse, and C. Monroe Nat. Phys. 12, 907 (2016).
- A. De Luca, B. L. Altshuler, V. E. Kravtsov, and A. Scardicchio Phys. Rev. Lett. 113, 046806 (2014).
- V. E. Kravtsov, I. M. Khaymovich, E. Cuevas, and M. Amini New J. Phys. 17, 122002 (2015).
- M. Žnidarič, A. Scardicchio, and V. K. Varma Phys. Rev. Lett. 117, 040601 (2016).
- E. P. Wigner Ann. Math. 67, 325 (1958).
- M. L. Mehta, Random Matrices (Academic Press, Boston, 1991).
- V. Zelevinsky, B. A. Brown, N. Frazier, and M. Horoi Phys. Rep. 276, 85 (1996).
- E. J. Torres-Herrera, J. Karp, M. Távora, and L. F. Santos Entropy 18(10), 359 (2016).
- T. A. Brody, J. Flores, J. B. French, P. A. Mello, A. Pandey, and S. S. M. Wong Rev. Mod. Phys 53, 385 (1981).
- Y. Avishai, J. Richert, and R. Berkovitz Phys. Rev. B 66, 052416 (2002).
- M. Távora, E. J. Torres-Herrera, and L. F. Santos Phys. Rev. A 94, 041603R (2016).
- M. Távora, E. J. Torres-Herrera, and L. F. Santos, Power-law decay exponents: a dynamical criterion for predicting thermalization, arXiv:1610.04240.
- P. Pechukas J. Phys. Chem. 88, 4823 (1984).
- L. Leviandier, M. Lombardi, R. Jost, and J. P. Pique Phys. Rev. Lett. 56, 2449 (1986).
- A. Delon, R. Jost, and M. Lombardi J. Chem. Phys. 95, 5701 (1991).
- Y. Alhassid and R. D. Levine Phys. Rev. A 46, 4650 (1992).
- F. M. Izrailev J. Phys. A 22, 865 (1989).
- B. I. Shklovskii, B. Shapiro, B. R. Sears, P. Lambrianides, and H. B. Shore Phys. Rev. B 47, 11487 (1993).
- A. Gubin and L. F. Santos Am. J. Phys. 80, 246 (2012).
- L. F. Santos J. Math. Phys 50, 095211 (2009).
- P. R. Zangara, A. D. Dente, E. J. Torres-Herrera, H. M. Pastawski, A. Iucci, and L. F. Santos Phys. Rev. E 88, 032913 (2013).
- P. Jacquod and D. L. Shepelyansky Phys. Rev. Lett. 79, 1837 (1997).
- L. F. Santos and M. Rigol Phys. Rev. E 81, 036206 (2010).
- E. J. Torres-Herrera and L. F. Santos Phys. Rev. E 89, 062110 (2014).
- J. M. G. Gómez, R. A. Molina, A. Relaño, and J. Retamosa Phys. Rev. E 66, 036209 (2002).
- F. Evers and A. D. Mirlin Rev. Mod. Phys. 80, 1355 (2008).
- Y. Y. Atas and E. Bogomolny Phil. Trans. R. Soc. A 372 (2014).
- For fixed , decreases monotonically as the ratio between the Ising interaction strength and the flip-flop term strength increases (26).
- E. J. Torres-Herrera and L. F. Santos Phys. Rev. E 88, 042121 (2013).
- L. F. Santos, F. Borgonovi, and F. M. Izrailev Phys. Rev. Lett. 108, 094102 (2012).
- L. F. Santos, F. Borgonovi, and F. M. Izrailev Phys. Rev. E 85, 036209 (2012).
- K. He and M. Rigol Phys. Rev. A 87, 043615 (2013).
- M. Rigol Phys. Rev. Lett. 112, 170601 (2014).
- L. Amico, R. Fazio, A. Osterloh, and V. Vedral Rev. Mod. Phys. 80, 517 (2008).
- N. Laflorencie Phys. Rep. 646, 1 (2016).
- D. N. Page Phys. Rev. Lett. 71, 1291 (1993).
- The study in (29) actually considered the diagonal entropy, which is the Shannon entropy of the initial state written in the energy eigenbasis (121); (122).
- E. J. Torres-Herrera and L. F. Santos Phys. Rev. A 89, 043620 (2014).
- E. J. Torres-Herrera, M. Vyas, and L. F. Santos New J. Phys. 16, 063010 (2014).
- E. J. Torres-Herrera and L. F. Santos Phys. Rev. A 90, 033623 (2014).
- C. E. Porter, Statistical Theories of Spectra: Fluctuations (Academic Press, New York, 1965).
- R. Ketzmerick, G. Petschel, and T. Geisel Phys. Rev. Lett. 69, 695 (1992).
- B. Huckestein and L. Schweitzer Phys. Rev. Lett. 72, 713 (1994).
- G. P. Berman, F. Borgonovi, F. M. Izrailev, and A. Smerzi Phys. Rev. Lett. 92, 030404 (2004).
- V. V. Flambaum and F. M. Izrailev Phys. Rev. E 64, 036220 (2001).
- L. A. Khalfin Sov. Phys. JETP 6, 1053 (1958).
- J. G. Muga, A. Ruschhaupt, and A. del Campo, Time in Quantum Mechanics, vol. 2 (Springer, London, 2009).
- K. Urbanowski Eur. Phys. J. D 54, 25 (2009).
- A. del Campo New J. Phy. 18, 015014 (2016).
- F. Borgonovi, F. M. Izrailev, L. F. Santos, and V. G. Zelevinsky Phys. Rep. 626, 1 (2016).
- L. D’Alessio, Y. Kafri, A. Polkovnikov, and M. Rigol Adv. Phys. 65, 239 (2016).
- E. P. Wigner Ann. Math. 62, 548 (1955).
- A. Polkovnikov Ann. Phys. (N.Y.) 326, 486 (2011).
- L. F. Santos, A. Polkovnikov, and M. Rigol Phys. Rev. Lett. 107, 040601 (2011).