Many-body localization transition in large quantum spin chains: The mobility edge

Thermalization of random-field Heisenberg spin chain is probed by time evolution of density correlation functions. Studying the impacts of average energies of initial product states on dynamics of the system, we provide arguments in favor of the existence of a mobility edge in the large system-size limit.

Alternatively, time evolution of large [30] (or even infinite [31]) disordered many-body systems can be simulated with tensor network algorithms. Reaching large time scales, necessary to assess thermalization properties [32] is challenging especially in the vicinity of transition to ergodic phase. Nevertheless, such an approach allows to obtain estimates for critical disorder strength for large system sizes [30,32], in quasiperiodic systems [33], or even beyond one spatial dimension [34,35]. An advantage of such an approach is that it directly mimics experimental observations of MBL [36][37][38][39][40][41][42][43].
Exact diagonalization study [12] of the Heisenberg chain showed that there exists, for a fixed disorder strength W , a threshold energy such that eigenstates above (below) that energy are extended (localized). This means that the system possesses the so-called many-body mobility edge. The question of precise location of mobility edges, even in a well understood case of single particle Anderson localization [44,45] is not fully resolved [46,47]. Expectedly, the situation is even more complex for interacting systems. The existence of many-body mobility edges was questioned in the work [48] which argues that local fluctuations in presence of a mobility edge would induce a global delocalization. That would leave out only the two possibilities: either the states at all energies are delocalized or localized. The mobility edges found in earlier works [12,49,50], were argued in [48] to be indistinguishable from finite-size effects.
The aim of our work is to study many-body mobility edges at much larger system sizes than those available to exact diagonalization in an attempt to verify conclusions of [48]. To that end, we employ Chebyshev polynomial expansion of the evolution operator [51][52][53] and the time-dependent variational principle (TDVP) applied to matrix product states (MPS) [54][55][56][57] to simulate time dynamics of random-field Heisenberg spin chain. Our approach is, in spirit, similar to that of [58] (and used for bosons in [59]). However, instead of considering an injection of controllable amount of energy into ground state of the system, we consider time evolution of initial product states with specified average energies, exactly similar to what was done recently in spin quantum simulator [60]. Probing time decay of density correlation functions allows us to estimate the critical disorder strength as a arXiv:2006.02860v1 [cond-mat.dis-nn] 4 Jun 2020 function of energy of the initial state. Studying systems of size up to L = 100, we perform a finite size scaling of our results which provides arguments in favor of existence of mobility edge even in large systems.
The paper is organized as follows: Sec. II briefly introduces the disordered model as well numerical techniques considered in this study. We present our main results regarding the estimation of many-body mobility edge via the time evolution of large systems in Sec. III . Finally, we conclude in Sec. IV.

II. THE MODEL AND METHOD
We consider 1D random-field Heisenberg (XXZ) spin chain with the Hamiltonian given by where S α i , α = x, y, z, are spin-1/2 matrices, J = 1 is fixed to be the unit of energy, and h i ∈ [−W, W ] are independent, uniformly distributed random variables. In this work, we consider open boundary conditions in the Hamiltonian (1). The random-field Heisenberg spin chain has been widely studied in the MBL context [12,31,[61][62][63][64][65][66][67], which has made it the de facto standard model of MBL studies.
The transition between ergodic and MBL phases is reflected in change of statistical properties of energy levels of the system. A common approach is to consider the gap ratio r i = min{Ei+2−Ei+1,Ei+1−Ei} max{Ei+2−Ei+1,Ei+1−Ei} , where E i are the energy eigenvalues of the system. Averaging the gap ratio over part of the spectrum of the system and over disorder realizations, one obtains an average gap ratio r, which differentiates between level statistics of ergodic system [10,68], well described by Gaussian orthogonal ensemble of random matrices, for which r ≈ 0.53 and between Poissonian statistics of eigenvalues in MBL phase (for which r ≈ 0.39). The later arises due to emergent integrability resulting from the presence of local integrals of motion [69][70][71][72][73][74][75].
To reveal the dependence of ergodic-MBL transition on energy, the gap ratios r i are averaged over only a certain number of eigenvalues with energies close to a rescaled en- it the energy of the ground (highest excited) state. Such a calculation of average gap ratio (supported with results for other probes of localization) for random-field Heisenberg spin chain reveals that the ergodic region has shape of a characteristic lobe on the phase diagram in variables of the rescaled energy and disorder strength W [12] . The average gap ratio, obtained in exact diagonalization of random field Heisenberg spin chain of size L = 16, is plotted as a function of and W in the background of Fig. 1.
To probe the transition between ergodic and MBL phases with time evolution, we propose the following pro- tocol. We consider an initial state |ψ = |σ 1 , . . . , σ L , where σ i =↑, ↓ are chosen randomly with constraint that the average rescaled energy ψ = ( ψ|H|ψ − E min )/(E max − E min ) of this state lies withing the range [ − δ , + δ ] corresponding to a given rescaled energy , where δ is a small tolerance (we take δ = 0.01).
To calculate ψ for L 26 we find E max , E min with the standard Lanczos algorithm [76]. For larger system sizes, E min and E max are calculated using density matrix renormalization group (DMRG) algorithm [77][78][79][80][81](see Appendix A).
Our quantity of interest is the density correlation function where the constant D assures that C(0) = 1 and l 0 > 0 diminishes the influence of boundaries (in our calculations, we take l 0 = 2). The standard deviation of the rescaled energy (3) is smaller than 0.1 for disorder strengths that we consider in this work as shown in Fig. 2. Those relatively small fluctuations of energy suggest that the properties of eigenstates at the rescaled energy can be well probed by time evolution of the state |ψ and reflected, in particular, by the density correlation function C(t).

III. QUENCH DYNAMICS
A. Disorder strength dependence Fig. 3 shows the density correlation functions C(t) obtained for the random-field Heisenberg spin chain of a fixed size L = 20, for rescaled energy = 0.5 of the initial state. The correlation function decreases in time, with some oscillations superimposed [84]. For small disorder strength, e.g. W = 2.8 the eigenstate thermalization hypothesis [5,85] is valid for the system, and in the long time limit the correlation function vanishes C(t) t→∞ → 0 as system loses the memory of the initial state. In contrast, for large disorder strength, e.g., W = 5, a non-zero stationary value of the correlation function C(t) t→∞ → c 0 > 0 is admitted showing that the system is non-ergodic. The first experimental signatures of MBL were obtained in study of time evolution of imbalance [36], quantity analogous to the density correlation function -for quantitative comparison of the two quantities see Appendix B.
At large times (t > 100), the decay of the correlation function is well described by a power law, C(t) ∝ t −β . Griffiths rare regions are one possible explanation of this behavior [86]. However, it was shown experimentally and numerically that time dynamics in quasiperiodic potentials, where Griffiths regions are necessarily absent, have analogous features [64,87,88]. Regardless of the origin of the power law decay of the correlation function, the disorder strength dependence of the exponent β can be used to locate the onset of ergodicity breaking in the system.
The exponent β governing the decay of the density correlation function is shown in Fig. 4. Let us first concentrate on the results in the middle of the spectrum ( = 0.5). In the considered interval of disorder strength W , the exponent decreases exponentially with W with a good approximation β ∝ e −W/Ω . The large number of disorder realizations (10000) used in calculation of C(t) allows us to see that even at the large disorder strength W = 5 the exponent β = 4.1(4)·10 −3 is non-vanishing. If the power-law decay C(t) = a 0 t −β prevailed for t → ∞, the density correlation function would vanish in the longtime limit and the system would be ergodic. This, however, does not happen for L = 20, as after the so-called Heisenberg time t H discreteness of spectrum manifests itself in saturation of the correlation functions [89][90][91][92], so that one would observe C(t) t→∞ → c 0 > 0 for W = 5 and L = 20. The Heisenberg time t H increases exponentially with the system size L. This illustrates a difficulty in locating the MBL transition using time dynamics of large systems on time scales of few hundred J −1 accessible to tensor network methods (or to current experiments with e.g., ultra-cold atoms): one cannot predict whether a slow decay of correlation functions governed by an exponent β 1 observed, for example t ∈ [100, 500], will eventually lead to C(t) t→∞ → 0 or not. To resolve the difficulties, the work of [30] assumes that the value of the exponent β must be vanishing within error bars to be compatible with saturation of correlation functions in the long time limit. The drawback of this criterion is that the error bar of β depends on the number of disorder realizations used in calculation of the correlation function. Therefore, we introduce a cut-off β 0 : disorder strength W C (L) for which β = β 0 is regarded as disorder strength for transition to MBL phase at system size L. Exact diagonalization results show that: i) collapse of data for L 22 gives a critical disorder strength W C ≈ 3.7; ii) the similar values W C ≈ 3.8 or W C ≈ 4.2 are obtained in asymmetric scaling on the two sides of the transition; iii) the breakdown of the volume-law scaling of entanglement entropy gives an estimate W C = 3.75 at system size L = 20 [29]. The obtained results for β at L = 20 and = 0.5 show that the cut-off value β 0 = 0.014 is consistent with the above estimates for the critical disorder strength obtained from exact diagonalizations. Consequently, throughout this work, we use β 0 = 0.014 as a threshold value which separates ergodic and MBL regimes for all system sizes and energies of the initial state we consider.
The values of β presented in Fig. 4 show that the increase of disorder strength W slows down the dynamics more severely for rescaled energies of initial state different than = 0.5. Notably, the exponent β decreases exponentially with W : β ∝ e −W/Ω (where Ω is a constant) in a wide regime of disorder strengths. This resembles the scaling of Thouless time t T h ∝ e −W/W0 observed in exact diagonalization data in [25].

B. System-size dependence
Density correlation function C(t) for larger system sizes are shown for two exemplary pairs of disorder strength W and initial rescaled energy in Fig. 5. The decay of C(t) at large times is well fitted by an algebraic dependence C(t) ∝ t −β . The exponents β obtained in the fitting of power-law decay to C(t) are shown for two exemplary values of the rescaled energy of the initial states in Fig. 6. For a given disorder strength W , we observe a clear increase of β with increasing system size. Interestingly, the shift is, to a good approximation, uniform for all disorder strengths so that the exponential decrease β ∝ e −W/Ω (at sufficiently large W ) is observed for all considered system sizes. Let us mention here that we consider 400 realizations of disorder for L = 34, 50, and 200 realizations for L = 100 for each values of and W. For small system sizes (L = 20, 22, 26) we consider between 10000 and 500 disorder realizations.
We obtain estimates for disorder strength W C (L) for transition to MBL phase by finding the crossings of β(W ) curve for given system size L with the β = β 0 line. Results of this procedure are shown in Fig. 7. The disorder strength W C (L) depend, within the estimated error bars, linearly on the inverse of the system size L with clear growth of W C (L) as the system size increases. On one hand, this trend allows us, by means of a linear fit W C (1/L) = A/L + W C (∞), to extrapolate the results to L → ∞ and to obtain the estimate of critical disorder strength W C (∞) for transition to MBL phase.
On the other hand, we observe that the slopes A are similar for all of the considered rescaled energies of the initial state. Thus, the shape of the boundary between ergodic and MBL regimes observed for L = 20 does not change considerably when the system size is increased. This is visible in Fig. 1. The points for various system sizes L are precisely the values of W C (L) obtained from the condition β = β 0 . The characteristic shape of the lobe does not change when the system sizes increases from L = 20 to L = 100 and is preserved even after the extrapolation to L → ∞. Therefore, there exists a certain range of disorder strengths such that the states for < L M E are localized, states for L M E < < U M E are extended and states for > U M E are again localized. Thus, our results indicate that the system indeed possesses a many-body mobility edge in the thermodynamic limit.

IV. DISCUSSION AND OUTLOOK
Chebyshev polynomial expansion of the time evolution operator and the TDVP method applied to MPS allowed us to study the problem of energy dependence of the transition between ergodic and MBL phases in large disordered quantum spin chains. Introducing a cut-off value of exponent β of power-law decay in time of density correlation function, we were able to probe the transition for different rescaled energies of the initial state. For small system-sizes (e.g., L = 20), our approach gives results consistent with exact diagonalization. Importantly, our method allows to consider much larger system sizes (L = 100) for which it predicts an existence of a mobility edge.
The disorder strength W C (L) is a lower bound on the transition to MBL phase: the residual decay of density correlations with exponent β 0 is insufficient to restore the uniform density profile for system size L = 20, but it is possible that it leads to an eventual decay of correlation function for larger system sizes.
The protocol we considered is, in principle, experimentally realizable. Our results can be verified experimentally if the setup of [60] was scaled to larger system sizes. Many-body mobility edge arises also in disordered Bose-Hubbard models [93]. It can be probed by a quench protocol analogous to the one considered in this work. Since the bosonic models allow for occupations in each site larger than unity, density wave-like states that are easier to obtain in an experiment with ultra-cold atoms can be use to probe the many-body mobility edge [59,93].  For large system sizes, i.e., L > 26, we use matrix product states (MPS) ansatz [80,81] to represent the quantum state of the disordered Heisenberg chain. The energies of the ground state (E min ) and highest excited state (E max ) are obtained using standard density matrix renormalization group method (DMRG) [77][78][79][80][81] . E min can simply be found from the ground-state search of H of Eq. (1) using DMRG, while that of −H gives us −E max . In these calculations, we keep truncation error due to singular value decomposition (SVD) below 10 −12 , i.e., we discard states corresponding to smallest singular values λ n that satisfy n∈discarded λ 2 n n∈all λ 2 n < 10 −12 . For time evolution large systems, i.e., L > 26, we employ MPS based time-dependent variational principle (TDVP) method [54][55][56][57]. Specifically, we use 'hybrid' scheme of TDVP [32,57,82,83] with step-size δt = 0.1, where two-site variant of TDVP is used initially to dynamically grow the bond dimension of the MPS, until the bond dimension in the bulk of the MPS is saturated to a predetermined value χ max . One-site variant of TDVP is employed in the subsequent dynamics. In our calculations, we fix χ max to be 384, which gives reliable results up to t = 500. For a more detailed analysis of the convergence of TDVP in the disordered Heisenberg chain, see [32].
was used in previous works [30,32] to locate the transition to MBL phase. It is worth to mention that density wave states, analogous to the Néel state were used in experiment with ultra-cold fermions [36] that demonstrated signatures of MBL transition. Imbalance I(t) is the common name for the density correlation function (2) calculated for the Néel state. The imbalance behaves similarly to the density correlation function: initially it decreases and oscillates, subsequently it decays algebraically in time. Fig. 8 shows comparison of exponents β C and β Imb that govern the power-law decay of C(t) correlation function for rescaled energy of initial state = 0.5 (which is the rescaled en-ergy corresponding to the Néel state), and power-law decay of imbalance of Néel state. The exponents β Imb are consistently bigger than β C . This fact may be traced back to the number of states coupled by the Hamiltonian (1) to state |ψ Néel which is higher than the average number of states coupled to a typical Fock state with = 0.5. The absence of aligned spins on neighboring sites makes the thermalization process more efficient for Néel state. The inequality β C < β Imb justifies the cut-off beta value β Imb 0 = 0.02 for decay of imbalance in Néel state assumed in [32].