Thermalization of dilute impurities in one dimensional spin chains

We analyze a crossover between ergodic and non-ergodic regimes in an interacting spin chain with a dilute density of impurities, defined as spins with a strong local field. The dilute limit allows us to unravel some finite size effects and propose a mechanism for the delocalization of these impurities in the thermodynamic limit. In particular we show that impurities will always relax by exchanging energy with the rest of the chain. The relaxation rate only weakly depends on the impurity density and decays exponentially, up to logarithmic corrections, with the field strength. We connect the relaxation to fast operator spreading and show that the same mechanism destabilizes the recursive construction of local integrals of motion at any impurity density. In the high field limit, impurities will appear to be localized, and the system will be non-ergodic, over a wide range of system sizes. However, this is a transient effect and the eventual delocalization can be understood in terms of a flowing localization length.


I. INTRODUCTION
Understanding, and controlling, the conditions under which dynamical systems thermalize under their own internal dynamics is of fundamental interest and has important technological applications. Avoiding thermalization typically requires careful crafting of the Hamiltonian of the system, but it's been proposed that models with local interactions exhibit non-ergodic behavior, that is stable in the thermodynamic limit, when subject to sufficiently large disorder [1,2]. Following these initial publications there has been very extensive work on understanding this non-ergodic phase (currently going by the name of many-body localization (MBL)) and the nature of transition to the ergodic phase. Existence of a well localized regime has been reported by several state of the art experiments [3,4]. We refer the reader to some recent reviews for further references [5,6]. Nonetheless, several papers have recently questioned the stability of the MBL phase [7][8][9]. In turn the findings of Refs. [7,9] were challenged by some follow up papers [10][11][12][13]. Regardless of where one stands in this debate, one of the key challenges in direct numerical, or experimental, study of the MBL transition is that finite size (time) effects at larger disorder are very strong, making it hard to draw unambiguous conclusions. For example, very recently a series of numerical papers, based on newly developed approaches, moved the lower limit of disorder compatible with the MBL transition to much higher values, by factors of two to five more than was previously believed [14][15][16].
Early analytical approaches to MBL, starting from the pioneering works [1, 2], focused on the stability of the localized phase against the proliferation of resonances at strong enough disorder, in analogy with the noninteracting problem. The resonances are defined as near degeneracies between localized energy states, which are lifted by the hopping of particles and the interaction between them. It was argued that these resonances cannot destabilize the localized phase at sufficiently strong dis-order like in the non-interacting case. A formal mathematical argument for stability of the MBL phase was presented in Ref. [17] for some specific model. In the current work we address the problem from a different angle and come to an opposite conclusion: namely that in the thermodynamic limit the localized phase is always unstable because of off-resonant virtual transitions. While we do not provide a rigorous mathematical proof of this statement, we support our analytical results with a careful numerical analysis. Moreover some of the key numerical results, which were used to demonstrate stability of the localized phase, are consistent with our results.
The approach we develop in this paper is based on first understanding the fate of a single impurity, which is weakly coupled to an ergodic spin chain (bath). At a sufficiently strong local field this impurity undergoes a delocalization crossover as a function of either the bath size L or the impurity observation time t: if L or t is small the impurity is effectively localized, only weakly dressed by the bath spins. In the MBL language, it forms a local integral of motion (LIOM) [6]. However, for large L and at sufficiently long times the LIOM decays. We tie this instability to the Krylov complexity of the bath, which was recently proposed as a generic probe of quantum chaos [18][19][20][21]. Physically this instability manifests itself in a flowing correlation length ξ(x) with the distance: as the LIOM grows in support, its tails decay slower and slower, leading to eventual divergence of ξ(x). Interestingly we find a direct generic connection between the lifetime of the best (slowest decaying) LIOM and the Fermi Golden Rule (FGR) relaxation rate of the impurity spin. Using this approach we avoid the need of making any assumptions about the structure of the bath eigenstates and can work directly in the thermodynamic limit.
Having established the connection between the LIOM instability and the FGR rate for a single impurity, we go on to show that the presence of other impurities does not qualitatively change the situation. That is, any finite impurity density will only lead to a finite renormalization of the LIOM relaxation time. The flow and eventual divergence of the correlation length leading to instability of LIOMs is tied to the operator growth, which is not affected by disorder apart from a finite renormalization. The impurity model allows us to carefully study the effect of finite impurity density on the relaxation rate smoothly connecting decay of a single impurity in the presence of disorder with the clean limit. In this way we are able to make predictions about the thermodynamic limit and test them numerically while only using small systems. Such a study are more difficult in canonical MBL models with large random local fields on every site and much larger finite size effects. It is very hard to imagine that there would be any qualitative difference between our model and canonical MBL models. As we discuss later, our findings for the impurity model are in excellent qualitative agreement with both recent and earlier numerical simulations on fully disordered models. Unlike the previous studies we find that the mechanism of delocalization of the impurity in many body systems is not due to proliferation of the resonances but rather due to virtual off-resonant transitions.
While our analytical constructions are rather general, we use a specific model Hamiltonian to support them numerically, namely Here where S x,y,z j are spin-1/2 operators is the bulk Hamiltonian describing the bath and where {ℓ} is a subset of sites where impurities are located and V ℓ are uniformly distributed in the interval [V /2, 3V /2]. We allow impurities strengths to fluctuate around mean value of V to avoid dealing with any potential resonances. For one impurity this subset consists of a single site with a fixed strength V . The magnetic fields h j in H bulk are small and random uniformly and independently distributed on all sites in the interval [−1/4, 1/4]. These magnetic fields serve a two-fold purpose: firstly they break both integrability and translational symmetry of the Heisenberg chain; secondly averaging over disorder allows us to additionally suppress effects of accidental resonances. We checked that all our results reported in this work are valid for each disorder realization. Open boundary conditions are used unless otherwise stated. We also use the Hamiltonian (1) to explain how earlier analyses of level statistics is affected by our findings. The absence of scale separation between the freezing of impurities and the decoupling of segments of the chain in small systems, creates the illusion of a fixed crossing but this a purely finite size effect. Additionally we analyze the fidelity susceptibility for this system recently proposed by us as a probe of chaos [22,23] and show that its behavior for a single impurity is very similar to that of a fully disordered model [24]. Likewise we find signatures of the inverse frequency scaling of the spectral function (1/f -noise) of the bath spins in the presence of the strong impurity, which are also reminiscent of the results found in the fully disordered model [24].
The paper is structured as follows: In the next two sections we analyze the fate of a single impurity, coupled to a weakly disorder chain which serves as a bath, using Fermi's golden rule and a perturbative Birkhoff construction of the LIOM. We explain how these two apparently different approaches are in fact related through the Krylov complexity, and why they lead to the same criterion for the localization/delocalization crossover at approximately extensive impurity field. The single impurity results not only establish a baseline for understanding how to think about thermalization of the boundaries of rare regions in the putative MBL phase. It also allows us to systematically investigate the affects of additional impurities. We proceed to discuss why adding more impurities to the bath only quantitatively affects the position of this crossover. Finally, we discuss our findings in light of earlier analysis of numerical probes of MBL like level spacing statistics, the fidelity susceptibility and the spectral function of local observables. These probes again highlight qualitative similarity between a single impurity system and fully disordered models.

II. SINGLE IMPURITY
As a first step in understanding the fate of impurity spins in the Hamiltonian (1) we consider a setup where a single impurity is weakly coupled to an ergodic bath such that the Hamiltonian is It is convenient to separate the interaction term of the impurity with the bulk into H int . The small parameter ǫ is introduced to control our analytical results. In the numerical analysis of the model we use ǫ = 1. As it will become clear shortly the longitudinal coupling S z 1 S z 0 between the impurity spin and the bath plays no role in our analysis. Formally this term can be always absorbed in H bulk without affecting any results.

A. FGR relaxation
A standard way to understand relaxation of the impurity coupled to a bath is through the FGR. It is informative to look into the Hamiltonian H bi in the rotating frame defined by the interaction picture of the impurity Hamiltonian H 0 = V S z 0 , which results into mapping of a static Hamiltonian H bi into a Floquet system with no impurity potential but with a periodically driven hopping between the impurity and the boundary spin: The FGR relaxation rate can be extracted from the spectral function of the oscillating spin-spin coupling in the basis of the bath Hamiltonian. Because the matrix elements of S ± 0 are trivial with respect to | ↑ and | ↓ states of the impurity spin it suffices to analyze the spectral function of S x 1 (or equivalently S y 1 ) of the boundary spin A x (ω) defined as: where G n x is the connected correlation function: where {. . . } + stands for the anti-commutator. This spectral function is shown at Fig. 1 for three different system sizes L = 12, 14, 16. At high frequencies the spectral function behaves like where τ ≈ 3.4. Note that with increasing system size the spectral function simply extends to higher frequencies. At frequencies below the high frequency cutoff there are almost no finite size effects. This insensitivity of the high frequency response to the system size is consistent with Ref. [25]. The scaling (7) was predicted earlier as a decay rate of doublons [26,27]. It also saturates the upper bound for the spectral function recently derived in Ref. [18][19][20]. In some closely related models the same exponential form can also be shown to be its lower bound [21]. A slightly weaker bound with no log(ω) correction was derived earlier in Ref.
[28]. According to the fit shown in Fig. 1 this bound is tight and describes the actual spectral function well. In fact this scaling of the spectral function is very easy to understand from simple heuristic considerations. In order to absorb an energy ω ≫ 1, the system it is required to use roughly Cω links as this energy is locally not available. Here C is the constant of the order of one (recall that the spin-spin coupling J on the links is set to unity). Within standard perturbation theory each link will result in 1/ω contribution to the matrix element entering the transition rate, therefore one can estimate the total matrix element as (1/ω) Cω ∝ exp[−Cω log(ω)]. The square of this matrix element defines the spectral function and correspondingly the FGR decay rate of the impurity spin, which agrees with Eq. (7) if we identify τ = 2C. The FGR relaxation rate of a weakly coupled impurity to the boundary inherits the scaling from the spectral function [29]: As we already noted the use of FGR is formally justified if we assume that ǫ ≪ 1, though it is expected that this relationship between Γ and the spectral function holds even when ǫ = O(1). Note that there is an essential singularity in Γ(V ) at 1/V → 0 such that the relaxation rate cannot be captured in any finite order in perturbation theory in 1/V . The FGR relaxation should provide an effective mechanism for the impurity spin relaxation as long as it is much larger than the level spacing: The criterion is equivalent to demanding that the typical unperturbed susceptbility χ for switching on the coupling between the impurity spin and the rest of the chain would be larger than O(1), i.e. that the eigenstates of the impurity and the bath will fully mix with each other, see Appendix A. We thus conclude that the critical impurity potential separating the localized and delocalized regime scales as Up to the logarithmic correction the critical impurity potential separating localized and delocalized regimes scales linearly with the bath size.

B. Asymptotic Birkhoff construction of the LIOM
In small systems, where Γ ≪ ∆ FGR does not apply. Instead it is expected that the boundary spin will only partially relax and form a so called local integral of motion (LIOM) [6]. To test this idea, we will construct the LIOM in the leading order of perturbation theory in the coupling to the bath ǫ and in all orders in 1/V . This is exactly the same order of approximation which is used to derive the FGR. To construct the LIOM we use the so called Birkhoff normal form, where we build a conserved charge iteratively as a series in 1/V : requiring that in each order in 1/V the commutator [Q, H] vanishes to the same order in 1/V . This equation can be solved order by order. Using that [S 0 z , H bulk ] = 0 and {S 0 z , H int } + = 0 it is easy to check (see also Appendix B) that to linear order in ǫ and n-th order in 1/V the LIOM is given by The norm of nested commutators entering the expansion R k ≡ i k Ad k H bulk H int is tied to the parameter τ defining the FGR decay rate. Namely, where L is the system size. At large k this asymptotes to [18,19,21]: (12) Using cyclic properties of the trace it is easy to check that for any integers k and q we have Tr[R k R k+2q+1 ] = 0 and Tr[R k R k+2q ] = Tr[R 2 k+q ]. This observation allows us to exactly account for the interference between different terms in the expansion and express the norm of the conserved operator through the sum of norms of operators R k with positive coefficients: Expression (12) makes clear that the Birkhoff construction is asymptotic. At large V , the decay rate has a non-monotonic dependence on n. It is convenient to introduce the running localization length as This localization length flows with n diverging at At this value of n * the perturbative decay rate Γ 2 n reaches its minimum Apart from an overall prefactor the square of the short time decay rate of the LIOM: Γ 2 min coincides with the FGR rate Γ (8). This situation is not unexpected and a complementary discussion can be found in Ref. [14,15].
Physically for local Hamiltonians the index n represents the spatial range of the approximate LIOM Q n . The flow of the localization length ξ(n) for n < n * indicates that the decay of the tails of this LIOM slows down with the distance. Eventually the decay stops when the localization length diverges.

C. Variational conserved charge
One might wonder if this divergence can be regularized in some way leading to a better conserved charge Q n . To address this question we can use a variational approach using the same commutator ansatz as in perturbation theory as a basis, but allowing for arbitrary coefficients. Instead of computing the coefficients in front of nested commutators in Eq. (11) perturbatively we will assume that they are variational parameters. It is easy to check that in the limit n → ∞ this variational ansatz is exact in the linear order in ǫ. The contributions can be generated recursively as follows. Given an operator O n at order n, define The variational ansatz thus consists of an arbitrary operator in the Krylov subspace of the superoperator B(·) = [S z 0 , [H bulk , ·]]. It's insightful, and indispensable to perform the optimization numerically, to Gram-Schmidt orthogonalize the basis operators as they are generated. The procedure is similar to the familiar Lanczos procedure to generate Krylov space of the Liouvillian [18]. Given a charge q n at order n, define which is proportional to the next order term. To generate an orthonormal set it suffices to orthogonalize it with respect to the last two q n . Hence we define: where which makes the charges obey To generate the same set of operators as the Birkhoff construction from the previous section we use q 0 ∝ H int . Then we can write the variational conserved charge as The best variational solution could be defined as the one which minimizes the residual commutator with the Hamiltonian, i.e.
In the leading order in ǫ this yields the set of linear equations for ψ: Sψ = f with the matrix S having elements and the source vector f defined as Using the basic definitions of q k and the fact that [S z 0 , [S z 0 , q k ]] = q k it's rather straightforward to show that S is a real symmetric pentadiagonal matrix. In addition, the source term f is only non-zero at the first two entries k = 0, 1. Further analytical progress seems possible, but we postpone that to future work. We only point that in the limit of large V we can easily recover the perturbative solution, see Appendix C for more details. To proceed here we solve the problem numerically, and the results are summarized in Fig. 2. We make a number of observations, first at low order n the variational results agree with the perturbative construction from the previous section. At the crossover n * , where the perturbative result yields the minimal relaxation rate, the variational improvement stalls and the residual reaches a plateau. Second, since S is a Hermitian matrix the solution can be decomposed in the eigenmodes of S and we can investigate the stability of the problem by removing the most irrelevant modes one by one. For low orders n < n * , we observe a a drastic increase in the relaxation meaning that the best mode is well defined. At n ≈ n * this gap closes. In Fig. 2, we also observe a clear transition in the ground state of S from localized near the diagonal, to oscillatory with it's mode fixed near n * .
From this analysis we conclude that in finite size systems, as long as L n * ∼ V τ log(V τ ), there is a well defined LIOM, adiabatically connected to the boundary spin operator S z 0 . For larger system sizes this LIOM becomes unstable and delocalizes. In this sense the LIOM is similar to a long lived quasi-particle, which eventually decays. Two seemingly different criteria for localizationdelocalization crossover: i) FGR rate becomes of the order of the level spacing and ii) Birkhoff LIOM construction starts to break down, thus lead to the same estimate of the bath size corresponding to this crossover. This agreement between the two approaches is not accidental as both results are ultimately connected to the universal operator growth of the nested commutators of H int and H bulk . It also substantiates the idea that once the recursive Birkhoff construction breaks down, the system starts thermalizing; eventually becoming ergodic. In that sense the present work is entirely along the lines of seminal works by Abou-Chacra,Thouless and Anderson [31] and Basko, Aleiner and Altshuler [1], where the authors construct a self-consistent theory of localization, solve the equations order by order and interpret the instability of the construction as a sign of delocalization. While our construction is different, we establish a more direct link between both sides of the transition. Finally, in Appendix D we present a brief discussion on how the same construction can be applied to periodically driven systems, where instead of an impurity one can couple a harmonic oscillator to a spin chain. We tie heating to the divergence of the LIOM connected to the "photon" number.

III. FINITE IMPURITY DENSITY
Now let us see how the previous analysis is modified if we consider the full Hamiltonian (1) with a finite density of impurities. We will still use the impurity at the edge as a probe, i.e. analyze the Hamiltonian (4), where Clearly a straightforward application of the Birkhoff construction fails as H bulk contains terms of the order of the impurity potential V such that the expansion (10) becomes much more complicated, e.g. the probe impurity could resonate with some other impurity which could lead to an instability in the naive Birkhoff construction that would not necessarily imply delocalization. To tackle this problem we will first perform a Schrieffer-Wolff (SW) transformation on the bath Hamiltonian to effectively eliminate the impurity spins. In particular, if the bath contains a single impurity at a site ℓ then after the SW transformation we obtain the following effective Hamiltonian describing the bath where H L and H R describe the blocks of the bath Hamiltonian on the left and on the right of the site ℓ. By construction, the Hamiltonian is still diagonal in the bathimpurity spin and one can thus consider the two sectors with S z ℓ = ±1/2 independently. This transformed Hamiltonian can be obtained in two different ways: i) either by performing a standard unitary rotation, which perturbatively removes the coupling between the impurity and the rest of the spins: where or ii) by going to the rotating frame with respect to the impurity potential like in Eq. (5) and performing the leading order van Vleck expansion of the resulting Floquet Hamiltonian [32]. The transformed HamiltonianH SW in Eq. (28) is quite simple, apart from some boundary corrections arising from the bare interaction and virtual coupling with the impurity, it has some effective flip-flop contribution coupling the boundary spins of the two blocks. Since that coupling arises from a virtual process involving the impurity, it's suppressed by 1/V ℓ . Note that the sign of the coupling, determined by the value of S z ℓ , is irrelevant as it can be removed by a simple π-rotation of spins on the right of the impurity around the z-axis. We restrict ourselves to S z ℓ = 1/2 sector in Eq. (28). The other two 1/V ℓ corrections to the transformed Hamiltonian, representing small boundary magnetic fields, are unimportant for the physics of the model and we drop them, defining the effective bulk Hamiltonian: If the bath contains multiple impurities this procedure can be done on each impurity site effectively introducing weak links across them. The new bulk Hamiltonian H ′ bulk contains no large terms of the order of V and one can thus apply both the FGR and the Birkhoff analysis by simply replacing the bulk Hamiltonian H bulk → H ′ bulk in Eq. (4). First, we want to understand how the presence of multiple impurities affects the spectrum of the boundary spins shown in Fig. 1. A sufficient condition for localization is vanishing of the spectral function at ω ≥ V for a sufficiently large V . This would ensure that the FGR rate for impurity relaxation is zero. Conversely finite spectral weight at any (non-extensive) V indicates delocalization of the impurities. Clearly, when focusing on a single impurity, the rest of the system acts as the worst bath when all the other impurities are completely frozen out. It thus suffices to understand the modifications of the boundary spectral function due to the presence of additional weak links in the bulk of the bath. The corresponding spectral functions for the boundary spin S 1 x are shown in  Fig. 1 the spectral functions for different system sizes look identical up to the cutoff scale which increases with the many-body bandwidth. Compared to the case with no impurities, which also corresponds to the top blue line corresponding to V = 1/2, we see two jumps developing in the spectral function at ω ≈ 3.5 and ω ≈ 7. These jumps can be easily explained using the same heuristic argument as before: In order to dump a large amount of energy ω one has to excite τ ω/2 strong links (see discussion after Eq. (7)). However, after each ∆ℓ = 5 strong links in our setup there is a weak link, which almost does not contribute to the energy if J eff ℓ ≪ 1 but leads to an additional 1/(2V ℓ ) suppression to the matrix element and correspondingly 1/(2V ℓ ) 2 suppression to the spectral function and the FGR rate. This simple argument is confirmed numerically in the inset of the top panel Fig. 3 where the two lines show dependence of the drop in the spectral function on V at two values of ω indicated by the arrows in the main plot. The extracted jumps are well described by power laws consistent with the expected 1/(2V ) 2 (after one jump) and 1/(2V ) 4 (after two jumps) scalings. In the bottom panel of Fig. 3 we show similar results for weak links located after every third site as shown in the inset. Now the jumps appear more frequently but the magnitude of each jump is again consistent with V −2 scaling per block. We thus see that the spectral function of the model with weak links is described by This spectral function gives a lower bound on the spectral function of the full model (see Appendix E) and hence defines a lower bound on the FGR relaxation rate of the impurity We thus conclude that the lower bound of the FGR decay rate of the impurity is only weakly affected by the presence of other impurities, which somewhat increase the effective exponent τ → τ ′ . As a consequence, for any impurity at a finite energy V , or more accurately at V which increases slower with system size than L/ log(L), there is a sufficient spectral weight to dissipate energy into the bath.
Like in the single impurity case one can check stability of LIOMs when the FGR relaxation rate becomes smaller than the level spacing. The Birkhoff perturbative construction of the LIOM associated with the boundary impurity looks the same as in the single impurity case with the only difference that we will encounter a finite density of weak links in the nested commutators R k appearing in Eq. (13) such that the norms of such commutators will be suppressed by at most V 2k/∆ℓ if we assume that weak links appear in the rate 1/∆ℓ. Suppression is likely even less as the norm will be dominated by the terms containing fewer than average weak links. In either case this suppression is not enough to counter the factorial growth of the norms nested commutators R k . Moreover Refs.
[33] and [34] argued that for fully disordered models the asymptotic behavior of these norms at large k is not affected by the disorder potential except for finite renormalization of the parameter τ . As a result the LIOM associated with the probe impurity spin remains perturbatively unstable for any density of weak links.
To confirm this, we construct the same variational LI-OMs as for the single impurity problem within the effective weak link model. The results are summarized in Fig. 4, where we observe suppression of the residual commutator with J 2 eff when V = 3 and J 4 eff for V = 4. By increasing the impurity potential from 3 to 4 we increase n * enough so that it encompasses two weak links, substantiating once more that only weak links at a distance less than n * contribute to a suppression of the relaxation and they do so by suppressing the rate by J 2 eff per weak link. These results are again consistent with the steps observed in the FGR rate (see Fig. 3), where the number of active weak links scales with the impurity potential. Finally, the effective weak link couplings can be chosen consistently with the boundary spin V , by fixing them to the SW value J eff = 1/2V . Figure 5 shows the LIOM decay rate, defined as the plateau value of Γ 2 opt , for different impurity configurations and completes the picture.Similarly to the analysis of the FGR rate, we see that finite density of impurities simply shifts localization/delocalization crossover at given V to somewhat larger system sizes. We emphasize again that this instability is associated not with proliferation of resonances but with factorially growing number of virtual transitions encoded in the operator spreading.
As for a single impurity, the breakdown of localization can be understood from the flowing localization length. A careful argument put forward in Ref.
[35], known as an avalanche instability, states that if the correlation length of the LIOMs ξ becomes larger than a constant of the order of the lattice spacing the localized phase becomes unstable to unbounded growth of any ergodic seed. One can thus alternatively interpret delocalization of the impurity spins at any disorder strength as an avalanche induced by a flowing localization length ξ(n) with the distance n. While we only established the flow of ξ(n) in the weak coupling limit to the bath ǫ, it does not look plausible that the situation changes in the higher orders in ǫ. Indeed for ǫ ≪ 1 the shape of the LIOM is ǫ independent, so one would have to imagine very exotic scenarios where ξ(n) is a non-monotonic function of the distance to stabilize LIOMs. The flow of ξ(n) was observed numerically in two recent works which study the decay rate of the slowest operators in fully disordered models (see Fig. 2 in Ref.
[15] and Fig. 5 in Ref. [14]). Moreover, a careful analysis of the numerical data in earlier papers claiming to see the exponential scaling of the LIOMs (constant correlation length) reveals that it actually flows considerably with the system size again in agreement with our results (see for example Fig. 2 in Ref.

IV. LEVEL STATISTICS AND FIDELITY SUSCEPTIBILITY
In the discussion above we established that the localization-delocalization crossover of a single impurity coupled to a bath is only weakly affected by the presence of other impurities, i.e. by the presence of a disorder potential. Let us now look into two other independent measures, both popular probes to study localization in disordered models. The aim of this section is to establish that these measures agree with our previous analysis, i.e. that they show same qualitative behavior for a single impurity as for fully disordered models. Because we will study the system as a whole, i.e. without using a probe spin we consider a setup where a single impurity is added in the middle of the chain such that where L is the chain size which we choose to be odd [38]. This problem was studied earlier in the literature [39][40][41][42] focusing in the regime V ∼ 1. Here we again concentrate in the regime of large V analyzed above. Fig. 6 shows the mean ratio of energy level statistics as a function of the impurity potential for three different short chains of length L = 13, 15, 17. Given subsequent energy level spacings s n = E n+1 − E n , with H = n E n |n n|, this ratio is defined as

A. Level Spacing Statistics
For non-ergodic systems and Poissonian level statistics, the average over eigenstates r ≈ 0.386, whereas for chaotic systems with GOE statistics r ≈ 0.5307 [43]. At sufficiently small impurity potential V , the system The chain has a weak link J eff after every third site and the external impurity field is V = 3; color goes from blue to red with decreasing J eff ranging from 1 to 100. The inset shows the plateau value scales like J 2 eff . (Panel B) Shows the same decay rate for an impurity field of V = 4, with the inset highlighting the plateau value now scales as J 4 eff , ranging from 1 to 20. is observed to be ergodic, as expected. Upon increasing the potential V , ergodicity gets broken in a seemingly two-step way. First, there is a fast drop in r , followed by a much slower further decrease of the level spacing ratio to the Poissonian value. Furthermore, the required V for the initial deviation from the GOE value shifts significantly with system size L. This initial drop is caused precisely by localization of the impurity happening at extensive (up to log corrections) V * ∝ L and agrees with the FGR and Birkhoff predictions for the localization threshold. The further slow decay of r is a consequence of the resulting fragmentation of the chain, which occurs at much larger potential V * * ∝ 2 L/2 . So there is a parametrically large window V * ≪ V ≪ V * * where the impurity is localized and yet the rest of the system is ergodic. Thus the single impurity model is a specific example of a system with a parametrically large difference between the potentials required to localize the impurity spin and to fragment the Fock space into several (three for our setup) disconnected sectors [44,45].
To understand the emergence of the asymptotic behavior of r at large V it is convenient to analyze the effective spin model (31), where the impurity spin is integrated out via a Schrieffer-Wolff transformation. The level spacing statistics of the effective model is illustrated by the dash-dotted lines in Fig. 6. At sufficiently large impurity potential V they asymptote the full model. Note that at large V the full model is better approximated by an unfolded effective Hamiltonian H ′ + V S z ℓ , which consists of two decoupled identical blocks, corresponding to the different values of the conserved magnetization of the impurity spin. The level statistics of this unfolded Hamiltonian is illustrated in Fig. 6  still remains frozen such that the effective model is still accurate but the two blocks start to overlap pushing the level statistics closer to the Poisson value. And indeed we see that the solid lines much better approximate r of the full model. As V decreases further the impurity gets delocalized in the full model such that r approaches the GOE ratio, and so does the effective model. However, the unfolded Hamiltonian always consist of two decoupled blocks and as they overlap more and more with decreasing V , r is pushed down closer and closer to the Poisson value. We thus conclude that the domain of agreement between the data coming from the full model and the unfolded effective model corresponds to the localized impurity regime. The initial drop in r in the full model from the GOE value is therefore associated with localization of the impurity. The remaining physics can be understood within the effective model. The fact that the magnitude of the jump in r decreases with the system size is consistent with the expectation that V * corresponding to the freezing of the impurity scales approximately linearly with L. In this case the energies of two blocks are extensively separated leading to a very small overlap between the corresponding energies and hence a small drop of r .
To highlight the significance of finite size effects on the interpretation of numerical results, we briefly ana-lyze a two-impurity configuration. Figure 7 shows the level spacing ratio for the model with two impurities of opposite strength V and −V located as shown in the inset. Comparing these results to Fig. 6, it becomes immediately clear that the drift in the impurity freezing remains similar, however the drop in level repulsion becomes significantly larger. The latter is easy to understand, as there are now four decoupled blocks once the two impurities have localized. Moreover, because two impurities have been frozen out, the remaining effective model becomes non-ergodic at a smaller impurity potential V , which has the same exponential scaling with the system size, but with considerably enhanced finite size effects. For two smaller system sizes (red and blue) one even observes a crossing at large values of V , which is often interpreted as a signature of the localization transition We move on to analyzing a multiple impurity setup corresponding to a constant spacing between them ∆ℓ = 5. The corresponding dependence of r on V is plotted in Fig. 8. The black lines and dots show the results for the full models of sizes L = 10 and L = 15 as illustrated in the inset. The largest system size corresponding to the green configuration has L = 20 and is outside of reach of exact diagonalization. Nevertheless we can extrapolate the other two lines noting the drift to the right of the departure from the GOE statistics on top and drift to the left of the departure from the folded effective model (full colored lines). This extrapolation would almost certainly lead to a good crossing point with the two other sizes at V ≈ 2.9. However, we see that this feature is an entirely spurious effect. It comes from the real drift of the drop position in statistics to larger values of V with L and simultaneous increase in the drop magnitude with L coming from increasing number of effective blocks corresponding to different frozen impurity arrangements. If we look into the effective model folded or unfolded we see a very clear indication that the effective model is ergodic with no crossing in r developing in the unfolded model and a crossing strongly drifting to the larger values of V with L.

B. Fidelity Susceptibility
The fidelity susceptibility χ, or equivalently the diagonal component of the quantum geometric tensor with respect to some coupling λ can serve as a very sensitive probe of quantum chaos [22][23][24]. Specifically, it has been established that at the crossover from an integrable to an ergodic regime the fidelity susceptibility saturates its upper bound, diverging with the system size as χ ∝ exp[2S(L)], where S(L) is the infinite-temperature entropy of the system. For comparison, in integrable regimes χ diverges at most polynomially with the system size and in the ergodic regime it diverges as exp[S(L)]. For a given eigenstate n the fidelity susceptibility is defined as [46, 47] For concreteness, we use the longitudinal magnetization of the spins in the bulk of the system as a probe, i.e. ∂ λ H = S z 3 . To avoid dealing with large fluctuations due to the broad distribution of χ n in the non-ergodic phase, we look at the typical susceptibility, defined as where the expectation is over all eigenstates and realizations of the weak disorder in the chain. It is convenient to scale χ by the ergodic value corresponding to V = 0 and analyze the ratio χ(V )/χ(0), which should saturate at L independent value in the ergodic regimes and diverge exponentially at the localization transition. This scaled susceptibility is plotted in Fig. 9. In Fig. 9 (A) we illustrate the susceptibility for the full model. At small values of V , i.e. on the ETH side, we identify a good collapse of the data followed by a clear peak in the susceptibility with a height that approximately scales like χ ∼ e 2S(L) . The inset shows the extracted peak position with system size, the latter is linear to good approximation with a numerically extracted slope V * ∝ 0.26L. This expectation up to a log(L) correction fully agrees with the scaling extracted earlier comparing the FGR relaxation rate and the level spacing (see Eq. (9)). The log(L) correction is not visible in numerics due to small system sizes. Further note that the level spacing ratio r (see Fig. 6) at the peak susceptibility is close to the GOE value. The latter is consistent with recent works on MBL [7, 8, 23, 24].
In Fig. 9 (B) we perform the same analysis on the effective model H ′ . Note that folding does not affect χ as the eigenstates in both blocks do not talk to each other. Once more, we observe a peak in the susceptibility, indicating ergodicity breaking in the effective model. However, this time the peak develops much slower and as such appears to drift much faster with the system size. Again, for available system sizes the peak happens at a rather high value of r , where there is still a considerable difference in r between the folded and unfolded models. The inset shows the drift of the peak position on a logscale with the best fit. This drift is well approximated with a linear curve, indicating this time that the critical interaction needed to decouple the effective model with the best linear/exponential fit. Physical system sizes corresponding to the full and effective models shown in the same color are identical, but the as the impurity spin in the effective model is frozen its actual system size is reduced by one.
into the independent left and right blocks scales exponentially with L. The standard expectation, following from many-body perturbation theory, is that the strength of the effective hopping J eff = 1/(2V ) coupling two blocks of length L/2 sufficient for thermalization scales as J eff ∼ exp[−S(L)/2] = 2 −L/2 [23, 42, 48]. Mathematically, this criterion comes from requiring convergence of the leading perturbative correction to eigenstates and an assumption that the spectral function of the perturbation ∂ λ H is flat at small frequencies. The latter assumption is indeed correct (see the inset in Fig. 1). This criterion would predict that V * ∝ exp[L log(2)/2] ≈ exp[0.35L], which gives a somewhat larger slope than that in the inset of Fig. 9 (B). The discrepancy could be due to small system sizes leading to the small dynamical range and/or relevance of various log(L) corrections affecting the observed scaling. Exponential enhancement of the fidelity susceptibility implies an exponential (in L) enhancement of spectral weight at low frequency from O(1) to O(exp(S(L))), accompanied by exponentially slow (in L) relaxation [22,24]. To confirm that this is the case it is thus instructive to look directly at the spectral function of the impurity, which is shown in Fig. 10 for various strengths of the impurity potential V . At intermediate V , before a significant fraction of the magnetization has become conserved and the associated amount of spectral weight has been transferred to ω = 0, we observe a clear 1/ω 2 scaling at low frequencies. The latter was recently observed in other systems with slightly broken integrability [23,49]. This scaling is indicative of Lorentzian line broadening. In turn, the Lorentzian shape of the spectral functions suggests that the relaxation of S z ℓ is simply governed by Fermi's golden rule (FGR) (see the last Appendix in Ref. [23] for a detailed discussion). In passing we note that the spectral function of the bulk spin defining the fidelity susceptibility plotted in Fig. 9 shows slower 1/ω subdiffusive scaling behavior. It is illustrated in Fig. 12 in Appendix F and agrees with the results reported by us earlier in Ref. [24] for a fully disordered model. This 1/ω scaling corresponds to a very slow, logarithmic in time, relaxation which is somewhat surprising for the effective model.

V. CONCLUSION AND OUTLOOK
In this work we have presented a numerical and analytical study of one dimensional Heisenberg chains with dilute sets of defects, being spins with a large external field. We first analyzed the crossover from a localized to a delocalized regime for a single probe impurity weakly coupled to an ergodic bath. We showed that this crossover can be explained from the ergodic side by comparing the FGR decay rate and the mean level spacing of the bath and from the localized side by the divergence of the Birkhoff construction of the LIOM connected to the impurity spin. Interestingly, both approaches give the same criterion for the localization/delocalization crossover.
We tied the divergence of the Birkhoff construction to the Krylov complexity of the bath. In local interacting models (disordered or not) this complexity saturates its upper bound resulting in (almost) factorial growth of norms of nested commutators and as a result to the instability of the LIOM. In this way we avoid any need of making any assumptions about the eigenstates of the bath and can work directly in the thermodynamic limit. Thus we concluded that adding a finite density of disordered sites does not affect the fact that in the thermodyanmic limit there is no localized phase, but does quantitatively affect both the time scales at which impurity delocalizes and the length scale of the crossover between the localized and delocalized regimes. Let us comment that MBL is often argued to be related to localization on graphs like random regular graphs [50] or Caley trees [51]. From the point of view of the Birkhoff construction there is a huge qualitative difference between them and the local models. The nested commutator norms R k on such graphs can only grow exponentially with k such that the Birkhoff construction converges at a sufficiently large impurity potential (see Eq. (14)) even if the bath is not disordered, i.e. ergodic. So a weakly coupled impurity to such a system at a sufficiently large V would be localized (at least in the small ǫ limit). Adding disorder to the system will simply shift the localization transition to a smaller value of V .
We also analyzed numerically various other proxies for egodicity, such as level spacing ratio's, fidelity susceptibilities and spectral functions. All these measures point to the same conclusion that, regardless of the impurity density or the potential, the impurities ultimately relax in the thermodynamic limit by dissipating energy in the remaining bath. Nonetheless the dynamics is exponentially slow in V .
Our conclusions are opposite to previous works which argue for the stability of the MBL phase based on the analysis of the effect of resonances [1, 2, 17, 52]. The physical mechanism of instability, which we found here, is based on virtual non-resonant processes and is ultimately tied to the operator growth which is absent in non-interacting systems. This instability develops at V ∝ L/ log(L), which corresponds to energies, where resonances cannot play a role simply because of a small density of states near the edge of the many-body spectrum. Of course, perturbative divergence of the decay rate does not exclude that there are some other, non-perturbative mechanisms stabilizing LIOMs. But given that our analytical predictions fully agree with all known to us numerical data, as well as with the variational approach (see Appendix C), we find this scenario very unlikely.
We believe that our analysis is fully consistent with most, if not all, numerical results on the MBL transition. Our work suggests that in thermodynamic limit instead of MBL there is a transient glassy-type phase characterized by finite-time subdiffusive, logarithmic in time, spreading of correlation functions. Such slow transport was previously attributed to the MBL phase [60]. Heuristically one can think about this transient regime as a stage of slowly dephasing "quasi-particles" or lbits/LIOMs, which eventually crossovers to their subsequent diffusion. Such a crossover from sub-diffusive to faster transport is not unique to disordered systems and was observed in other setups, see e.g. Refs. [11] Piotr Sierant, Maciej Lewenstein, and Jakub Zakrzewski, "Polynomially filtered exact diagonalization approach to many-body localization," Physical Review Letters 125 (2020), 10.1103/physrevlett.125.15660.
[12] David J. Luitz and Yevgeny Bar Lev, "Absence of slow particle transport in the many-body localized phase," Phys. Rev. B 102, 100202 (2020 where we used that [S z 0 , H bulk ] = 0 which yields the solution We can now continue this construction. In the next order we need to solve the equation Let us point out that the equation only admits a solution for X if the operator A is odd under parity transformation generated by σ z 0 = 2S z 0 : This follows e.g. by multiplying both sides to the equation above by σ z 0 on the left and on the right. If this condition is satisfied then it is easy to check that where B is an arbitrary operator commuting with σ z 0 . Because H int is an even operator and H int is odd, the parity of any nested commutator of these two operators is determined by whether H int appears even or odd number of times. The RHS of Eq. (50) is obviously odd such that Here we set the arbitrary commuting operator B to zero, which as it will become clear shortly is justified in the linear order in ǫ. In general B should be chosen to cancel all even terms appearing in T n . We can now continue this construction iteratively solving the equation (54) This yields the expansion: which as it is easy to see is equivalent to the expansion (10) in the main text if we relabel 2n → n. Alternatively, one could consider a finite system of size L such that the sum in Eq. (10) converges at sufficiently large V and sufficiently small ǫ. Then, the expression above can be resummed to the infinite order, leading to The norm of this operator can be straightforwardly computed in the eigenbasis of the uncoupled Hamiltonian H bulk + V S z 0 : where χ is the eigenstate-average fidelity susceptibility, which we introduced in the previous appendix; ± sign refers to "up" and "down" sectors of the spin S z 0 . This result once again leads to the conclusion that the norm of the conserved charge is related to the fidelity susceptiblity, which as it was already shown in the previous appendix is related to the ratio of the FGR to the level spacing. Because the conserved part of magnetization scales as Tr(S z 0 Q)/Tr(Q 2 ) ∼ 1/ Q 2 we conclude that the condition ǫ 2 χ ≫ 1 ↔ Γ ≫ ∆ implies that this conserved magnetization is small.
While in this paper we focused on quantum systems, let us point that this LIOM construction applies to the classical setup, where the Hamiltonian (4) is expressed not in terms of spin -1/2 operators but in terms of continuous angular momenta satisfying Poisson bracket relations: plus cyclic permutations. Then it is easy to check that the Birkhoff construction for the conserved charge Q, satisfying {Q, H bi }=0 in the linear order in ǫ proceeds as follows: where f (0) y = S y 1 , and for q > 0: The decay rate of this conserved charge is completely analogous to Eq. (14), where R k = Ad k H bulk H int with the "Ad" operator implying the nested Poisson brackets. The norm of the function R k is defined through a phase space average over orientations of all spins: where θ i and φ i are the spherical angles defining spin orientations. As it was argued already in Ref.
[18] the scaling of these norms are expected to have the same factorial scaling for generic local Hamiltonians.
Eq. (14)). When the potential V becomes comparable to b 2n ∼ C n/ log(n) this scaling breaks down and crosses over to much slower decay of the rate with n as shown in the main text. From (65) it also becomes clear why the eigenstates of the quadratic form minimizing the decay rate and shown in Fig. 2 change their structure from positive highly localized states at small n to oscillating delocalized states at large n.

APPENDIX D: BIRKHOF CONSTRUCTION FOR PERIODIC DRIVING
As we discussed analysing the FGR decay, in the rotating frame the impurity problem maps to the Floquet problem. In the Floquet language, the existence of the LIOM is equivalent to the existence of a local Floquet Hamiltonian. While Floquet driving is not the focus of this work, let us briefly show that the Birkhoff's LIOM construction can be applied directly to Floquet systems, where instead of to an impurity, we couple the system to a photon field such that it is described by the Hamiltonian where Ω is the photon frequency, which plays the same role as the impurity potential V , a and a † are the photon creation and annihilation operators and H int can be either a boundary spin of the bath, H int = S x 1 or the total transverse magnetization H int = j S x j or something else. In the rotating frame in the limit of a large photon number this Hamiltonian becomes periodically driven with frequency Ω. We can now construct the LIOM coupled to the photon number using the same spirit as before: This Birkhoff construction has the structure identical to that for the impurity problem. Therefore we can draw identical conclusions about the asymptotic nature of the LIOM. Note that the breakdown of the LIOM in this case indicates that the photon delocalizes and the system heats up. In this context, it's worth to note that this result questions recent claims about MBL stabilizing downconverters (which are sometimes referred to as discrete time crystals) in the thermodynamic limit [65,66] . It should be noted though that if the driving field couples to a Hamiltonian which does not have the factorial growth of the nested commutator norms, then at least in the leading order in ǫ the LIOM converges and there is no heating. This happens, in particular, when one adds a small amplitude drive to the Hamiltonian which consists of a sum of mutually commuting terms. To study heating in those systems one has to extend the Birkhoff construction beyond linear (or possibly any other finite) order in ǫ. We note that in this class of systems it was found numerically that heating is indeed very strongly suppressed [67-69]

APPENDIX E: SPECTRAL FUNCTIONS FOR THE FULL VS. EFFECTIVE MODELS.
In the main text we analyzed the boundary spectral function for a set of coupled blocks, such that the impurities were exactly frozen out. It is intuitively clear that freezing the impurities makes the model more nonergodic transferring the low frequency spectral weight to high frequencies. In this Appendix we will verify that this is indeed the case. In Fig. 11 we show the boundary spectral function log(A x (ω)) for the full model corresponding to the exact same parameters as the effective model spectral function shown in Fig. 3 (B). That is, we add 3 impurities to a chain of 15 spins. The impurity spins will cause resonances at ω ∼ V . In panel C we compare the spectral functions for the full and effective models for two particular values of V : 10.81, 5.85. As we argued here and in the main text while ω ≪ V resonances do not occur and the spectral functions of the two models are indistinguishable. However, at larger frequencies the spectral function for the full model exhibits a non-monotonic behavior due to resonant flipping of the impurity spin. Such resonant process is absent in the effective model where the other impurities are frozen and consequently its spectral function keeps monotonically decreasing with ω.

APPENDIX F: SPECTRAL FUNCTION OF THE BULK SPINS
In the main text, the impurity spectral function is shown to have a 1/ω 2 dependence at low frequency at sufficiently large V , where the impurity starts to decouple. In the context of MBL, similar spectral functions have been analyzed and it's been argued that they should have sub-diffusive scaling on the ergodic side leading up to the transition, i.e. A(ω) ∼ ω 1−1/z , where z goes from z = 2 at weak disorder to z = O(L) at the transition. In general, the sub diffusive behavior is attributed to Griffith's effect, where exponentially rare regions with exponentially slow transport generate anomalous transport behavior. This picture suffers from a number of problems, most notably that the same phenomenology is observed in systems with quasi-periodic potentials in which there are no rare regions. Recently, many-body resonances have been proposed as an alternative explanation [42]. A different phenomenological explanation of this spectral function recently emerged from the work of L. Vidmar et. al. [70], which proposed a scenario of a broad distribution of the FGR relaxation rates. In Fig. 12 we show the spectral function for a spin in the bulk of a block, i.e. for the third spin in the chain, for two weakly coupled blocks described the effective Hamiltonian (31). This is precisely the same spin for which we computed the fidelity susceptibility shown in Fig. 9. One observes a broad region where the spectral function has behavior that is close to 1/ω.