A criterion for many-body localization-delocalization phase transition

We propose a new approach to probing ergodicity and its breakdown in quantum many-body systems based on their response to a local perturbation. We study the distribution of matrix elements of a local operator between the system's eigenstates, finding a qualitatively different behaviour in the many-body localized (MBL) and ergodic phases. To characterize how strongly a local perturbation modifies the eigenstates, we introduce the parameter ${\cal G}(L)=\langle \ln (V_{nm}/\delta) \rangle$, which represents a disorder-averaged ratio of a typical matrix element of a local operator $V$ to the energy level spacing, $\delta$; this parameter is reminiscent of the Thouless conductance in the single-particle localization. We show that the parameter ${\cal G}(L)$ decreases with system size $L$ in the MBL phase, and grows in the ergodic phase. We surmise that the delocalization transition occurs when ${\cal G}(L)$ is independent of system size, ${\cal G}(L)={\cal G}_c\sim 1$. We illustrate our approach by studying the many-body localization transition and resolving the many-body mobility edge in a disordered 1D XXZ spin-1/2 chain using exact diagonalization and time-evolving block decimation methods. Our criterion for the MBL transition gives insights into microscopic details of transition. Its direct physical consequences, in particular logarithmically slow transport at the transition, and extensive entanglement entropy of the eigenstates, are consistent with recent renormalization group predictions.


I. INTRODUCTION
Experimental advances of the past decade have lead to the realization of isolated quantum many-body systems of cold atoms and trapped ions. These systems have slow intrinsic time scales and are amenable to various quantum-optics experimental probes, thus they are ideally suited for probing far-from-equilibrium quantum dynamics [1]. The studies of cold atoms and trapped ions have inspired a theoretical quest for the universal framework to describe non-equilibrium quantum phenomena. One of the central challenges that has emerged is to understand the microscopic mechanisms of ergodicity and its breakdown in isolated quantum systems [1].
It has been established that there exist two distinct generic classes of many-body systems: the ergodic ones, which thermalize as a result of quantum evolution, and the many-body localized (MBL) [2][3][4][5], which are disordered and avoid thermalization via a mechanism similar to Anderson localization [6] in the Hilbert space. Ergodic and many-body localized systems have drastically different spectral and dynamical properties. Microscopic mechanism of thermalization in ergodic systems is called the Eigenstate Thermalization Hypothesis (ETH) [7,8]. According to the ETH, thermalization in ergodic systems is manifest already in the properties of individual eigenstates, in which physical observables are thermal. In contrast, the eigenstates of MBL systems violate the ETH due to the emergence of extensively many local integrals of motion (LIOMs) [9][10][11][12][13].
One of the invaluable tools for probing ergodicity that has recently emerged is the entanglement entropy. In ergodic systems, the eigenstates are similar to random vectors in the Hilbert space (modulo global conservation laws, e.g. energy); therefore, the eigenstates have an extensive, volume-law entanglement entropy. In contrast, the eigenstates of MBL systems can be obtained from product states by quasi-local unitary transformations, and have low entanglement that obeys the "area law" [9,14]. Further, ergodic and MBL phases are distinguished by the dynamics of entanglement following a quantum quench from initially non-entangled states: in the former case, entanglement spreads ballistically and grows linearly in time [15,16], while in the latter case the spreading is logarithmic in time [17][18][19]. The latter property has been understood from the phenomenological theory of the MBL phase based on the LIOMs [9,10].
Recently, first signatures of MBL have been observed in systems of cold atoms [20]. In this experiment, disorder can be tuned in a broad range, which allows one to probe MBL and ergodic phase, as well as the transition between them. The MBL-delocalization transition is a new kind of a dynamical phase transition, across which the eigenstates undergo a dramatic transformation as their entanglement entropy changes from area-law to volumelaw. Although it is conceivable that certain many-body systems exhibit an intermediate non-ergodic phase, here we focus on the case of a direct transition between the MBL and ergodic phases, which is suggested by numerical studies of 1D systems [5,18,[21][22][23][24][25][26][27], as well as by a recent phenomenological renormalization group study [28], and an effective percolation model [29]. Developing a full theory of this transition remains a major challenge.
In this paper, we propose a new approach to study- . Vertices of the lattice correspond to the unperturbed eigenstates, and matrix elements of the operator V determine the hopping amplitudes. Bottom panels illustrate the possible outcomes of a local perturbation applied to the Hamiltonian: all the eigenstates are mixed in the ergodic phase (panel c), whereas in the MBL phase the perturbation admixes a finite number of states (panel e). In the critical regime, the eigenstates are expected to be multifractal [32].
ing MBL, ergodicity, and the transition between them. We introduce a single dimensionless parameter defined in terms of physically measurable quantities -matrix elements of local operators -which provides a direct and sensitive probe of MBL. We formulate a criterion of MBL-delocalization transition in terms of this parameter, and explore the physical properties at the critical point. The parameter can be viewed as a many-body generalization of the Thouless conductance [30] in the single-particle Anderson problem; the latter plays a central role in the scaling theory of localization [31]. The main idea of our approach, summarized in Fig. 1, is to study the response of the system's eigenstates to a local perturbation V . We consider a lattice many-body system with an unperturbed Hamiltonian H and eigenstates |n , H|n = E n |n . The problem of finding the spectrum and eigenstates of H + V can be viewed as a hopping problem on a lattice where sites correspond to the eigenstates |n of unperturbed Hamiltonian H, with on-site energies E n = E n + V nn and hopping amplitudes V nm between sites n, m (see Fig. 1). Thus, all the information is encoded in eigenenergies E n and in the offdiagonal matrix elements V nm . Note that for local V , V nn ||V ||, and in general the order of energies E n is different from E n .
A key insight is that perturbing the Hamiltonian with a local operator V modifies the eigenstates very differently in the MBL and ergodic phases [see Fig. 1(a-b)].
In the MBL phase, the perturbation V can only affect the degrees of freedom within a localization length ξ away from its support. Thus, the perturbation strongly admixes a given eigenstate with only a finite number of other eigenstates [ Fig. 1(e)], as expected from the existence of LIOMs in the MBL phase. In contrast, in the ergodic phase V hybridizes an extensive number of eigenstates, therefore its effect is highly non-local, and the eigenstates of H + V are completely delocalized in the basis of the unperturbed eigenstates.
To detect the localization/delocalization transition, we ask whether hopping induced by V hybridizes the eigenstates with neighbouring values of E n (the closest states in the many-body spectrum). To characterize this hybridization, we introduce the parameter where we have ordered the eigenstates by their modified energy E n , and ε = E n /L = E n /L+O(1/L) is the energy density. As suggested above, the scaling of G with system size L should be qualitatively different in the MBL and ergodic phases.
In the MBL phase, at L ξ the eigenstates close in energy typically have very different spatial structure due to their different values of LIOMs, and a local operator couples them exponentially weakly in L. Therefore, in the MBL phase G(L) ∝ −L. On the other hand, in the ergodic phase the ETH implies that a local perturbation mixes the energy levels very strongly, and this requires G(L) ∝ +L. The above considerations suggest a natural criterion for the delocalization transition at the energy density ε to be G(ε, L) ∼ 1. This condition implies that a local perturbation significantly changes the structure of the eigenstates, such that the LIOMs become non-local and the effective hopping problem enters the critical regime. These results will be analytically and numerically substantiated in Sections II and III below. It will be shown that G, as a single parameter, characterizes the delocalization transition.
Previous numerical studies of the delocalization transition have relied on the statistics of energy levels [5], as well as the fluctuations of entanglement entropy [22]. In contrast, our approach allows us to directly probe the localization properties of the eigenstates through the matrix elements of local operators. Unlike entanglement entropy and level statistics, these are directly related to correlation functions, and therefore are, in principle, measurable. Moreover, our approach yields further consequences for the properties of the system at the critical point. Those include the extensive entanglement in the eigenstates, and logarithmic in time transport of conserved quantities and entanglement propagation. These properties are qualitatively consistent with the predictions of recent phenomenological renormalization group studies [28,29].
The rest of the paper is organized as follows. In Section II we provide analytical derivation for the scaling of G with energy density and system size in the MBL and ergodic phases, and formulate the delocalization criterion. Using exact diagonalization, in Section III we then test our approach numerically on the random-field XXZ spin chain, where the transition between MBL and ergodic phases can be driven by changing the disorder strength. We demonstrate that the G(L) ≈ 1 criterion gives a location of the transition that is in good agreement with the estimates obtained using other methods [24]. Furthermore, we show that our approach has a good energy-density resolution, and apply it to map out the phase diagram as a function of disorder and energy density, demonstrating the presence of the manybody mobility edge in the random-field XXZ model. In Section IV we address the physical consequences for the MBL transition that follow from our analysis. In particular, we find the logarithmic in time propagation of entanglement and particle number fluctuations, in agreement with recent works [28,29]. We numerically verify those predictions using exact diagonalization and timeevolving block-decimation simulations of a global quench in the random-field XXZ model. Our conclusions and a discussion of future directions is presented in Section V.

II. STATISTICS OF MATRIX ELEMENTS AND THE SCALING OF G: ANALYTIC CONSIDERATIONS
In this Section we discuss the distribution of the parameter G introduced in Eq. (1) and its scaling in the MBL and ergodic phases. In the former case, the behaviour of G can be understood by invoking the ETH. In the latter, our considerations are based on the picture of LIOMs in MBL phase. In the subsequent Section, we present data from exact diagonalization which are fully consistent with these analytical arguments.
A. Ergodic phase: matrix elements from the ETH In the ergodic phase, the eigenstates are expected to be similar to random vectors in the Hilbert space (this analogy is qualitative due to locality of interactions and energy conservation).
Consider two eigenstates |n , |m near the middle of the many-body band. They satisfy the orthogonality constraint m|n = 0, n = m. Intuitively, acting on the state |n with a local perturbation, V |n , gives another random vector in the Hilbert space, which is no longer orthogonal to |m . Thus, it is expected that the off-diagonal matrix element V nm = m|V |n is a scalar product of two random vectors, |m and V |n . Denoting the Hilbert space dimension by D, the matrix element can be estimated as V nm ∼ 1/ √ D. For a quantum many-body system, in general D ∼ exp(sL), where L is the number of lattice sites or spins, and s is statistical entropy density. Since the level spacing scales as the inverse Hilbert space dimension, δ ∼ 1/D, we have δ V nm , and a local perturbation strongly mixes the nearby eigenstates. Our parameter G, which in the leading order is given by G(L) ∼ ln(V nm /δ) ∼ 1/2 ln D ∝ L, thus grows linearly with system size.
The intuitive argument above can be extended to states with an arbitrary energy density using Srednicki's [33] ansatz for the off-diagonal matrix elements of local operators: where S(E) is the statistical entropy at energy E = (E n + E m )/2, R nm is a random number of order one, and f is a smooth function which can be linked to the thermallyaveraged response of the system to the perturbation V .
Since δ ∝ e −S(E) , we obtain that for states with an arbitrary energy density ε = E/L, where s(ε) is the corresponding entropy density. It is worth noting that the off-diagonal matrix elements in the ergodic phase have been studied numerically in Ref. [34], where the scaling V nm ∝ 1/ √ D was confirmed.

B. MBL phase: matrix elements from local integrals of motion
To understand the properties of the off-diagonal matrix elements V nm and G deep in the MBL phase we rely on the picture of LIOMs [9][10][11][12]. It was argued that the key property of a system that shows many-body localization at infinite temperature is that its eigenstates can be connected to product states by a quasi-local unitary transformation U , which entangles the remote degrees of freedom only exponentially weakly.
For concreteness, let us consider an interacting spin-1/2 model [e.g., the XXZ model introduced in Eq. (8) below]. In this case, it is convenient to choose the natural "up-down" basis, i.e., the basis vectors are eigenstates of all σ z i operators. Then, in the MBL phase there exists a quasi-local unitary operator U (quasi-local means that this operator only creates short-range entanglement), which diagonalizes the Hamiltonian in this basis, U † HU = H diag . The unitary U defines a mapping from the physical spin operators, σ α i , to "effective spins" τ α i via The operators τ z i form a complete set of mutually commuting, quasi-local integrals of motion in the MBL phase. In the τ -basis the eigenstates therefore simply become "up-down" states of the effective spins.
In order to analyze the effect of a local perturbation, we expand the corresponding local operator V (which for simplicity is chosen to act on the first few spins in the chain) over the complete basis of operators τ α1 1 τ α2 2 ...τ α L L , where α i ∈ {0, ±, z}. Here τ 0 i = 1 denotes the identity operator on the site i, and ± denote the usual raising/lowering operators. We represent V as a sum over operators V r of range r, which affect only the effective spins 1, . . . , r: where the last α r = 0. The off-diagonal matrix element of V between the eigenstates which differ by ∼ L flips of effective spins is equal to V α1...α L L , where ∼ L coefficients are α i = ± such that the operator τ α1 1 . . . τ αr r maps one configuration into another. For example, if |n 1 =↑↓↑↓, and |n 2 =↓↑↓↑, the matrix element n 2 |V |n 1 = V −+−+ gives us a coefficient in front of τ − 1 τ + 2 τ − 3 τ + 4 term in the expansion (5). Thus the matrix elements of V probe the nature of the spatial decay of the coefficients V α1...α L L , which, in turn, encode the properties of the unitary U (in particular, they control the amount of entanglement generated by U ).
Deep in the MBL phase, the LIOMs are robust under local perturbations. Typically, V hybridizes a given eigenstate only with states which have the same integrals of motion distance x ξ (ξ being the localization length) away from where V acts. Two generic eigenstates typically have different integrals of motion everywhere. Hence, in order for a local perturbation V to couple these states, it has to flip the state of effective spins across the entire volume of the system. This requires tunneling an "excitation" through the entire system, which would leave behind a trace of spin flips [35]. In the localized phase, such tunneling is exponentially suppressed; as a consequence, the probability of V to couple such states decays exponentially with system size. Thus, typically the ratio V nm /δ behaves as e −κL , and the parameter G decays linearly with the system size: Below we study G(L) numerically for the XXZ spin chain, finding the behaviour consistent with the above formula. We also find that G(L) has a broad, nearly normal distribution, similar to the distribution of the wave function tails in the single-particle Anderson insulator [36]. Note, that Eq. (6) indicates that matrix elements of a local perturbation are similar to a random ultrametric matrix [37]. Indeed, the scaling of G in Eq. (6) implies that a local perturbation which changes values of the LIOM up to a distance x away from its center [there are 2 x−1 such eigenstates], has a matrix element V (x) ∝ e −(κ+s)x , where s is the entropy density. For κ > 0 the eigenstates are localized, and the transition occurs at κ = 0 [37], which is similar to our criterion. The problem of Anderson localization for random ultrametric matrices has been studied analytically using strong-disorder renormalization group [38], and, in particular, multifractal exponents have been computed [37].
Note, however, the im-portant difference of the ensemble of ultrametric matrices considered in the literature, is that the matrix elements there obey the gaussian, rather than the log-normal distribution. We reserve the exploration of this analogy to further studies.
We note that the behavior of matrix elements in the MBL phase was used [39] to study the properties of a local spectral function. In addition, the spatial decay of the diagonal part of local operators given by terms with all α i ∈ {0, z} in Eq. (5) was studied in Ref. [13], and was also found to decay exponentially in the MBL phase.

C. A criterion for the transition
Using the predicted linear scaling of the parameter G with system size in the ergodic and MBL phases, it is natural to propose a phenomenological criterion for the delocalization transition as the absence of scaling of G with L: Below, we will use this criterion to determine the location of MBL-ergodic transition as a function of disorder strength and energy density in exact diagonalization. Physically, the above criterion implies that LIOMs, characteristic of the MBL phase, become unstable and are expected to delocalize at the transition. Indeed, any local perturbation strongly couples a given eigenstate with an adjacent eigenstate (in the reshuffled spectrum E n ), which has a different spatial structure, and therefore different LIOMs. In this case, the new eigenstates are superpositions of states with different values of the LIOMs across the entire system. Thus, a local perturbation modifies the integrals of motion even far away from where it is applied, which implies that integrals of motion necessarily become non-local.
We have argued that the above condition implies the non-locality of integrals of motion. It is natural to ask whether the criterion (7) indeed implies delocalization in the effective hopping problem we introduced in Section I above, and what are the properties of the system at this delocalization transition. That is, whether the participation ratio of new eigenstates in the basis of unperturbed eigenstates tends to infinity as L → ∞. One cannot answer this without additional information regarding the distribution of matrix elements. We have studied several quantities which probe the delocalization in the hopping problem directly, and indeed found that the above criterion is consistent with the delocalization in the hopping problem.

III. NUMERICAL RESULTS
In this Section we test the predictions for the distribution and scaling of G using numerical (exact) diagonalization of the one-dimensional random-field XXZ spin-1/2 (a) chain. The Hamiltonian of this model is given by: where we use periodic boundary conditions and identify σ L+1 ≡ σ 1 . Below we set J x = 1. The random z-field h i on each site is drawn from a box distribution [−W ; W ], and σ ± = (σ x ± iσ y )/2. At J z = 0, this model can be mapped onto free spinless fermions moving in a disorder potential via the Jordan-Wigner transformation; in this case, all eigenstates are localized for arbitrary disorder strength W . Introducing J z = 0 is equivalent to turning on the nearest-neighbour interactions between fermions.
Previous work [5] has demonstrated that the model (8) exhibits the MBL phase at strong disorder (W > W c ), and the ergodic phase at weak disorder (W < W c ). The delocalization transition takes place at some critical disorder strength W c . Several numerical studies [5, 21, 23-27, 40, 41] have estimated the location of the transition in this model using various probes, including the statistics of eigenenergies, entanglement entropy and its fluctuations, participation ratios and different dynamical probes. In particular, state-of-the-art exact diagonalization performed on systems up to L = 22 spins [24] has robustly determined the location of the transition at W c ≈ 3.6.
Here we use exact diagonalization to numerically study the statistics of the off-diagonal matrix elements and parameter G for the XXZ spin chain (8) across the MBLdelocalization transition. In order to establish the universality of this approach, we considered the response of the spin chain to two different local perturbations: The interaction strength was fixed to be J z = 1, similar to previous studies. All the data were obtained for chains with periodic boundary conditions constraining the Hilbert space to states with zero total spin, such that The behaviour of the distribution function of G(L) is illustrated in Fig. 2. In this case, the perturbation operator was chosen to beV 1 (we found a qualitatively similar behaviour forV 2 ). We observe that the distribution of G exhibits a qualitative change across the MBL transition. In the ergodic phase [ Fig. 2(a)], the average value of G(L) grows with the system size, signalling that a local perturbation is more and more likely to hybridize the nearby states -a signature of delocalization. If one scales the distribution by the square root of the Hilbert size dimension, G(L) / √ D, the distributions for different system sizes indeed approximately collapse, confirming that G scales as predicted by the ETH.
In contrast, in the MBL phase [ Fig. 2(c)], the averaged value of G decreases as system size is increased; this reflects the fact that a local perturbation is less and less likely to hybridize the nearby states, and therefore the eigenstates and the LIOMs remain robust. Moreover, the distribution of G becomes broader at larger L, and is approximately normal. We found that the distribution becomes very close to normal at large system sizes and strong disorder.
Finally, at intermediate disorder W = 3.6, the average G(L) remains nearly independent of system size. Hence it is natural to identify this point with the MBLdelocalization transition. Below we will use the energy-resolved average value of G to identify the transition point and the mobility edge. Note that although the average value of G(L) remains constant at the transition, the distribution broadens as the system size is increased. Detailed investigations revealed that this effect originates largely from sample-to-sample fluctuations, rather than from the state-to-state fluctuations in the same disorder realization.
B. Average energy-resolved G(ε, L) In order to obtain the location of the MBL transition more precisely, we now study the average value G(ε, L) as a function of the dimensionless energy density and system size. The dimensionless energy density ε (referred to as "energy density" in the following, for simplicity) is defined as ε = (E−E min )/(E max −E min ), where E min/max are the energies of ground state (highest excited state) of our system. Figure 3 shows the system-size dependence of G(ε c , L) for fixed ε c = 0.45, which is the energy density at which the delocalized phase is most robust. Similar to the behaviour already observed for the distribution of G, we see that the behaviour of averaged G(ε c , L) is qualitatively different at weak and strong disorder. At W 3 we have d G(ε c , L) /dL > 0, and the second derivative appears to be positive, signalling that larger systems become more and more thermal. At strong disorder W ≥ 4, G(ε c , L) behaves according to Eq.(6), as expected in the MBL phase. From Fig. 3 we identify the critical value of disorder W c = 3.6 ± 0.15 using our criterion for the MBL transition, Eq. (7). This  agrees with the previous findings of Refs. [24]. Using the same delocalization criteria, we can map out the many-body mobility edge in the random-field XXZ model. We define the many-body mobility edge to be at the energy density ε(W ) where G(ε c , L) is independent of system size. Results of this calculation are shown in Fig. 4. We also notice the intrinsic asymmetry of the mobility edge: it is shifted towards the ground state and centered slightly below ε = 1/2 [24,40,42]. This asymmetry becomes stronger with increasing the interaction strength (not shown), which agrees with the recent arguments of Ref. [42].

IV. ENTANGLEMENT AND DYNAMICS AT THE DELOCALIZATION TRANSITION
Next, we explore the physical properties of the system at the transition, defined using the above criterion. In particular, we will be interested in understanding the entanglement properties of the eigenstates, as well as the dynamical properties -spreading of entanglement and transport of conserved quantities. We discuss qualitative expectations and support them with numerical studies.
In order to qualitatively understand the entanglement properties of the eigenstates, we consider a system of size 2L and initially disconnect it across the middle link into left and right subsystems. The eigenstates of such a system are simply product states of the eigenstates in the left and right parts, |n L ⊗ |m R , n, m = 1, ..., D. Restoring the coupling between L and R systems corresponds to a local perturbation V LR = α V αL ⊗ V αR , where the sum includes a small finite number of terms (three for the case of the XXZ spin chain), and V αL , V αR act only on the degrees of freedom in the left/right part. The W ε Figure 4.
Many-body mobility edge ε(W ) as a function of disorder. Blue (red) color indicates regions where g(ε, L) decreases (grows) with L. Yellow regions correspond to points where we cannot determine the behavior due to errorbars. problem of finding the eigenstates reduces to a hopping problem, where the sites are product states |n L ⊗ |m R , and hopping amplitudes are set by the operator V LR . This problem bears many similarities with the case of a local perturbation acting on a single system, discussed above.
Previously in Section II we established that the MBL transition corresponds to the critical regime of the effective hopping problem. In this case, we expect that the eigenstates of the connected system |I are multifractal when expressed in the product basis |I = n,m A I nm |n L ⊗ |m R , with disorder-averaged participation ratios satisfying a relation P q = n,m |A I nm | 2q ∝ D −2τq . This is consistent with recent work [43] which studied dynamical signatures of the fractality near the MBL transition. We have tested and confirmed the multifractality of the hopping problem numerically (details will be presented elsewhere). Although we are not aware of a direct relation between entanglement entropy and participation ratios in the many-body case (for the singleparticle problem such a relation does exist, see Ref. [32]), it is natural to expect that the entanglement entropy of half the system will be extensive, but sub-thermal: where sL is the maximal (thermal entropy) of the subsystem of size L. This agrees with the predictions of Ref. [28]. We note that the above relation is expected to hold for entanglement between two regions of approximately equal size. If, on the other hand, one of the subsystems is much smaller than the other one, its entanglement entropy should approach the thermal value, in accordance with the general arguments based on the strong subadditivity of entanglement [44].
Next, we discuss the dynamics at the transition. Assuming the analogy with the problem of ultrametric random matrices, described above (although it should be kept in mind that this analogy is almost certainly incomplete -in particular, in our case the off-diagonal matrix elements follow a log-normal distribution rather than a random distribution), we expect that at the transition the typical relaxation time scale for a system of size L is exponentially long in L, and therefore the transport at the transition will be logarithmically slow, consistent with the predictions of Refs. [28,29]. Note that this is a time scale for both dephasing, and the on-shell dynamics (transport) at the MBL transition. For example, let us consider a setup in which the left and right subsystems are initially disconnected, i.e., are prepared in an eigenstate |ψ(0) = |n 0 L ⊗ |m 0 R , and the coupling between them is turned on instantaneously at t = 0. In the MBL phase after a quench, due to dephasing, the entanglement entropy increases logarithmically, and this growth is unbounded in an infinite system; but there is no transport (except for rare resonances). At the transition, however, there is transport of the z−component of spin polarization, which can be characterized by its fluctuations in the left subsystem as a function of time. Hence, both fluctuations of the magnetization and entanglement growth are unbounded in the thermodynamic limit and increase logarithmically in time, where S z,L is the total z−projection of spin in the left part. We expect that these laws will be generic and will hold for a class of initial states, including the case when the system is prepared in a product state (rather than in a product state of eigenstates of the left and right parts).
We have studied the entanglement growth and spin transport numerically to verify Eq. (11). We considered a global quench in which the system is initially prepared in a product state (random up-down configuration of spins), and time evolved with the Hamiltonian in Eq. (8). Using exact diagonalization, we have been able to access long-time dynamics in systems of up to L = 18 spins. In addition, we have performed time-evolving block decimation [45] (TEBD) simulations which allow access to much larger spin chains, but are restricted to moderate times t ∼ 10 2 [46]. For the TEBD algorithm we used the Trotter decomposition with time step ∆t = 0.01 and a maximum bond dimension χ = 2000. The results on the entanglement growth and the fluctuations of the spin polarization, obtained using TEBD for L = 24 spins, are illustrated in Fig. 5. Both quantities spread logarithmically in time. We note that slow transport opens an additional channel for entanglement growth, in addition to dephasing [18], and as a result, the entanglement S ent (t) near the transition does not display a characteristic "knee" which corresponds to the crossover between transport-dominated and dephasing-dominated growth deeper in the MBL phase.

V. SUMMARY
We presented a new approach for studying the MBL and delocalization transition based on the response of system's eigenstates to a local perturbation. We characterized the effect of a local perturbation by a dimensionless parameter G, defined via the ratio of a matrix element of local operator to the level spacing [Eq. (1)]. This parameter is reminiscent of the Thouless conductance in the single-particle localization problem, which describes the response of energy levels to changing the boundary conditions. Similar to the single-particle case, we argued that the transition can be identified by the relation G(L) ≈ G c ∼ 1. However, it should be noted that, unlike the Thouless conductance, the parameter G is not directly related to the physical conductance of a many-body system.
We have verified our approach by studying the MBL and delocalization transition in the random-field XXZ model. We obtained an estimate for the critical disorder strength at the transition which agrees very well with previous numerical studies based on entanglement and energy level statistics. Further, we mapped out the mobility edge in this model, demonstrating that the method has a good energy resolution, while at the same time requiring little beyond exact diagonalization of the Hamiltonian.
In addition to being an efficient tool to detect the MBL transition, our criterion gives insights into its microscopic details. In particular, our criterion implies that at the MBL transition, both transport of conserved quantities (such as spin polarization), as well as spreading of entanglement, are logarithmically slow. Such dynamics at the MBL transition was predicted by recent phenomenological RG studies of the MBL transition [28,29]. Using exact diagonalization and TEBD calculations, we have computed the spreading of entanglement and transport in large systems, confirming the logarithmic in time behaviour for both quantities.
Finally, we mention some further directions opened by this study. First, we expect that this method will be useful for studying MBL and ergodicity breaking in other contexts, for example in the translationallyinvariant models which were conjectured to break ergodicity [47][48][49][50][51][52][53]. Second, our results provide a natural starting point for developing a microscopic renormalization group (RG) procedure for the MBL transition. In spirit, the microscopic parameter G is similar to the phenomenological parameter (ratio of the "entangling rate" to the level spacing), introduced in a recent RG study [28], while having the advantage of being numerically measurable. Hence, one could potentially use the "physical" distribution of G, obtained from exact diagonalization, as an input parameter, and make use of the RG procedure to reach larger length scales. This would complement the existing phenomenological RG studies and is left for future work. Another interesting direction would be to further explore the consequences of the LIOMs in the MBL phase. As we discussed above, our results hint at an underlying ultrametric structure, which could allow one to introduce a random matrix ensemble from which the properties of the MBL phase and the delocalization transition could be analytically computed.