Ergodic inclusions in many body localized systems

We investigate the effect of ergodic inclusions in putative many-body localized systems. To this end, we consider the random field Heisenberg chain, which is many-body localized at strong disorder and we couple it to an ergodic bubble, modeled by a random matrix Hamiltonian. Recent theoretical work suggests that the ergodic bubble destabilizes the apparent localized phase at intermediate disorder strength and finite sizes. We tentatively confirm this by numerically analyzing the response of the local thermality, quantified by one-site purities, to the insertion of the bubble. For a range of intermediate disorder strengths, this response decays very slowly, or not at all, with increasing distance to the bubble. This suggests that at those disorder strengths, the system is delocalized in the thermodynamic limit. However, the numerics is unfortunately not unambiguous and we cannot definitely rule out artefacts.


I. INTRODUCTION
The discovery that noninteracting particles in a disorder potential can become completely immobile by P. W. Anderson in 1958 [1] has created a new field of study, which became enormously active in the last decade after it was predicted that this localized phase could persist in the presence of interactions, leading to a perfect insulator at any temperature [2][3][4][5][6][7].While it is difficult to realize such systems in condensed matter experiments due to the presence of phonons, many-body localized (MBL) systems were realized in synthetic quantum matter [8][9][10].It was first suggested that the interplay of interaction and disorder gives rise to a non-equlibrium phase transition between a thermal phase at weak disorder, which satisfies the eigenstate thermalization hypothesis (ETH) [11][12][13][14] and a many-body localized phase at strong disorder.Such a transition would be unparalleled in equilibrium [15][16][17][18][19][20][21][22].A large body of theoretical work now supports the picture that the many-body localized phase is characterized by an emergent complete set of quasi-local integrals of motion [23,24], which are fully consistent with the observed phenomenology of area law entanglement in all many-body eigenstates [25][26][27] as well as with the logarithmic post-quench entanglement production [28][29][30][31].
Recently, however, the stability of the insulating phase in the thermodynamic limit has been put into question, generating a debate on MBL as a phase of matter [32][33][34][35][36][37][38][39][40][41][42][43][44][45].Some works suggest that the critical disorder might be way higher than what was expected [37,38,[46][47][48], others even predict an infinite critical disorder in the thermodynamic limit [35,39,40,43,49].As a matter of fact, quantum simulations have shown MBL signatures Figure 1.System setup: A spin chain (blue) is coupled to an ergodic bubble (red).The spins at site 0 and L − 1 are both coupled to the bubble with coupling strength g, such that both "ends" of the chain are symmetric.In order to make the bubble perfectly thermal, its Hamiltonian is given by a random GOE matrix, while the spin chain corresponds to a XXZ spin 1/2 model in the presence of a random field coupled to local Ŝz i operators and nearest neighbour interactions.
in systems with few tens of particles [8,10,50,51], thus today's experiments can only access a MBL regime that asymptotically, in time and system size, may or may not thermalize.One central aspect in this debate is the delocalization transition mechanism and the crossover behavior in finite systems, whose understanding is still incomplete.In the vicinity of the transition, anomalously slow dynamics was observed [52][53][54][55][56], which was related to anomalous thermalization behavior [26,[57][58][59][60] and a theoretical description based on rare insulating inclusions was proposed [16,61,62].It remains however unclear how this picture can be reconciled with the observation of slow dynamics in quasiperiodic potentials [19,55,63].
It was also proposed that the many-body resonances are driving the slow dynamics in this regime [64].
It was pointed out in Ref. [65] that many-body localized systems are unstable under certain conditions towards thermal inclusions by a mechanism dubbed "avalanche" [66,67], and such a transition as a function of the localization length was confirmed numeri-arXiv:2308.01350v2[cond-mat.dis-nn]22 Aug 2023 cally in idealized models [68,69], also tailored to address implications for the instability of MBL in higher dimensions [65,70,71].Current activities now focus on the identification of this mechanism driving the transition in more realistic models [40,[72][73][74][75] and experiments [76,77].This is challenging and so far direct evidence for the avalanche mechanism in standard MBL models is still lacking.In this work, we directly address the issue of avalanches in a standard MBL model and consider a thermal inclusion coupled to a disordered spin chain.

II. MODEL
We study the random field XXZ chain, modeled by the Hamiltonian H XXZ ("the chain") coupled to an ergodic (or "thermal") bubble with Hamiltonian R0 .The setup for the model is illustrated in Fig. 1.The total Hamiltonian is of the form Here the first factor of the tensor product refers to the chain and the second factor refers to the thermal bubble.
The XXZ Hamiltonian acts on L spins labelled by i = 0, . . .L − 1; ( where we have taken the random fields h i ∈ [−W, W ] drawn from a box distribution.We consider the isotropic point ∆ = 1.All parameters are measured in J = 1 units.The coupling term reads: where the matrices Rl , l = 0, 1, 2 are independent random matrices from the Gaussian Orthogonal Ensemble (GOE), with a scaled bandwidth defined by ) where norm(0,1) are normal random variables with zero mean and unit variance.The dimension n GOE of the random matrices controls the power of the ergodic bubble, here we use n GOE = 3, 4, 5, 6, 8. β is a real variable that controls the band width of the random matrices and is chosen such that the level mixing (reflected in the overall gap ratio) is maximal for the largest possible range of parameters (see SM Sec.VI A for more details).
Note that the ergodic bubble is coupled to spins 0 and L − 1, thus closing the chain into a ring.Therefore, both ends of the chain are symmetric and the longest distance from the bubble corresponds to the spin situated at the middle of the chain i = ⌊L/2⌋.We also note that the coupling in Eq. (3) of the bubble to the chain is chosen such that all terms preserve the total spin Ŝz =

L−1 i=0
Ŝz i in the chain.This allows us to restrict our study to the largest sector of the Hilbert space, given by S z = 0 for L even and S z = 1 for L odd, with a Hilbert space of dimension dim (H) = n GOE L ⌊L/2⌋ .We perform massively parallel shift-invert diagonalization [78,79] for computing eigenstates close to the exact spectral center of each sample given by (E max + E min )/2, with (E min )E max being the (anti)ground state energy.We confront results for our system with a thermal bubble, Eq. ( 1) with the isolated XXZ chain Eq. ( 2).

Numerical protocol
We consider two types of observables.Firstly, we compute the r n -parameter of adjacent energy gaps in the middle of the spectrum given by r n = min(δ n , δ n+1 )/max(δ n , δ n+1 ) [5], with δ n = E n+1 − E n .Here E n−1 , E n , E n+1 are consecutive energy levels of the XXZ chain, and the system with a bubble, respectively.This parameter is useful to distinguish the ergodic and localized phase in a simple way.
Second, we use the single site purity γ i = Tr(ρ 2 i ), where ρ i = Tr {L−i} (|n⟩⟨n|) is the reduced density matrix on a single site i for a given eigenvector |n⟩.Due to the U (1) symmetry of the model, γ i can be conveniently expressed via the matrix element ⟨n|S z i |n⟩ as (a brief analysis of these matrix elements is presented in SM Sec.VI B ) In an ergodic system, the ETH and random matrix theory predict that γ i − 1/2 ∼ d −1/2 where d is the Hilbert space dimension, at least for states |n⟩ chosen at maximal entropy [14].The matrix elements ⟨n|S z i |n⟩ by themselves show qualitative signatures of thermalization, see SM Sec.VI B for a short discussion.
In models where ergodicity is induced from boundary effects, as it presumably happens in the case investigated here, one should be more careful and write eff,i where d eff,i is the effective dimension, see [26] and SM.In contrast, in a localized system, we expect γ i to depend substantially on the state |n⟩ and the disorder realization, and we expect the average γi to be given by a volume-independent value, tending to 1 as W → ∞.For each disorder realization, we use 50...100 eigenstates |n⟩ and for each set of model parameters at least 2000 disorder realizations of the fields h i and bubble Hamiltonian Rl .

III. SHIFT OF MBL TRANSITION IN THE INTERACTING CHAIN
In an infinitely long disordered chain, it should not matter whether we add a thermal bubble or not, because thermal Griffiths regions acting as ergodic patches are expected to be present anyhow.In a chain of moderate size, we expect that a substantial fraction of samples appears to be localized simply by lack of ergodic regions.This would lead to a shift of the apparent critical disorder value W c for short chains if one compares the isolated chain to our model where we add a bubble by hand.In Fig. 2, we compare the disorder averaged gap ratio r of adjacent energy gaps in the middle of the spectrum given of the isolated XXZ chain (right) to XXZ chain-bubble system (left) as a function of disorder strength W for different system sizes L. Using the crossing of r of size L and L + 2 as a proxy for the apparent critical point at this length scale, the inset shows that it is indeed the case that the crossings appear at slightly larger disorder strengths in the presence of the bubble.However, due to large statistical uncertainty and significant finite size effects, this crossing analysis is not a reliable way to pin down the critical point.

IV. CHANGE OF LOCAL THERMALITY IN INTERACTING CHAINS
We investigate the one-site purities γ i as a measure of local thermality.Throughout this paper, we consistently use the labels "localized" and "thermal" as inferred from the crossing point of r (see Fig. 2), i.e. with the critical point at W ≈ 3.5 − 4, even though recent works locate the transition at larger W . Inspecting the average purities in the presence and absence of the thermal bubble (see Fig. 3) is clear that the effect of the bubble at long distance is rather small and subtle.Therefore, it is useful to first map out what the theory of quantum avalanches predicts.Theoretical background For a simple avalanche model of our setup, we assume the existence of an apparent localization length ξ, describing the system in absence of the bubble, and a number p giving the probability that the bubble can kickstart a thermalization process in its near vicinity, see Ref. [69] for a detailed discussion and arguments why p ≪ 1 in our setup.A simple model, based on an unrealistic dichotomy between ergodicity and localization, leads (see SM) to the following very rough prediction: with γi the average purity at site i and γXXZ the average purity in the bubble-less system.For p = 0, we recover a perfectly localized system, whereas for p = 1, the system is ergodic but the average purity still increases with i, because the effective dimension d eff,i depends on the distance to the bubble.This effective dimension is defined as d eff,i = e −2i/ξ d therm with d therm the total Hilbert space dimension of the thermal region, given by , where ξ * = 1/ log 2 is the critical localization length.According to the above formula Eq. ( 6), there are two regimes in which γXXZ − γi does not decay to zero, or only very slowly.I) ξ > ξ * .Here, the bubble thermalizes a fraction p of samples, resulting in a shift γXXZ − γi which remains finite as i → ∞, even though it decreases due to the decrease of d eff,i .
II) ξ approaches ξ * from below.Then γXXZ − γi → 0 at large i, but the decay is arbitrarily slow when ξ → ξ * , because the decay of d eff,i is arbitrarily slow.
In practice, I) and II) are of course hard to distinguish.Large distance behaviour In Fig. 3, we indeed see a sign of influence of the bubble that does not, or only very slowly, decay with distance from the bubble.We observe that the purity in the presence of the bubble seems to tend to an asymptotic value that is significantly lower than the value for the bubble-less system.This seems to be the case up to disorder strength W = 5.6 − 6, after which the signal, i.e. the difference γXXZ − γi , becomes comparable to the error bars.We will analyze this further below, but let us first note that it is far from clear whether our system sizes are large enough to speculate about the limit i → ∞.Indeed, in Fig. 4, we investigate the dependence of the purity at fixed distance on increasing system size L. We see that for disorder values at and above W = 5.2, the purity γ 5 is actually increasing with increasing system size.This is an effect that is not even accounted for in our model equation and it should be interpreted as a sign that finite-size effects are still important (see SM material Sec.VI C for further analysis).Another way to look for thermalization induced by the bubble is through correlations between different sites of the chain.In that spirit we look at statistical correlations between two purities γ i and γ j .Those are defined as where γ i − γi captures the purity fluctuations around its mean and σ 2 i = (γ i − γ i ) 2 is the variance of the purity at site i.The quantity C(r) is the Pearson correlation coefficient between purities at site next to the bubble (i = 0) and sites at distance r from that spin with r = 0, 1, ..., L/2.This quantity for different disorder and system size is shown in Fig. 5. Interestingly, for W = 4.4 the bubble seems to be enhancing correlations compare to the bubble-less case.This effect is even growing with system size and it persists for W = 5.2 and, up to some extent, for W = 6.0.
Since the signals that we observe, the difference γXXZ − γi and the correlations in Fig. 5, are rather small, there is the concern that they might be caused by some artefact.For example, even in a well-localized system, in the absence of any avalanches, the purity-purity correlation C(r) from Eq. ( 7) has a non-zero limit as r → ∞ of order W 2 in the high-disorder limit W → ∞, see Appendix for more details.While we do not see any clear mechanism how this might pollute our analysis, we find it hard to exclude it.

V. CONCLUSION
We have studied the effect of ergodic inclusions modeled by local random matrices in disordered spin chains of up to 20 sites.We report little or no drift of the thermal-to-localized transition in the average level spacing compare to the bubble less-case, which was also observed in similar settings [72].We also investigate long distance effects of the bubble by looking at functions of the local magnetization expectation value, in this case the purity and its fluctuations.There, we have numerically witnessed a potentially divergent, long wavelength, response of an apparently localized chains (as estimated by ED studies like [78]) to an ergodic bubble.The effect is weak and not unambiguous, but it is compatible with the avalanche theory proposed in earlier works.Importantly, such weak correlations may be influenced by spurious long-range correlations due to the U (1) symmetry (see Sec. VI E in SM), although we can not confirm this from our numerics.
Recent studies also point out weak but persisting correlations produced by many-body resonances in the same disorder regime we study [37,47].Our results suggest that the bubble is effectively enhancing such weak correlations when the system looks well localized in average.However it is yet not clear how to directly relate the observed signal with the many-body resonances.

VI. SUPPLEMENTARY MATERIAL A. Parameters of the model
In this model two systems, the spin chain and the bubble are coupled via H coupling given by Eq. ( 3), where Eq. ( 4) shows the definition of the random matrices that described the bubble Hamiltonian.For coupling efficiently both systems their spectra should have comparable band widths which allows maximum possible mixing of their eigenstates.The parameter β controls the average band width of the random matrices, then an optimal choice of β enhances the mixing of the two systems and therefore maximizes the probabilities of the avalanche to happen.In order to do so, the average gap ratio ⟨r⟩ is computed for different values of band width β at several disorder strength (Fig. 6).We see that ⟨r⟩ always has maxima and they happen to be almost at the same β.Choosing β such that the average gap ratio is always maximum ensures that the spin chain and the bubble are strongly coupled.For 1.0 < W < 4.0 and n GOE = 3, 4, 6, 8 we found that β = 0.55, 0.6, 0.55, 0.5 respectively satisfy this requirement.It could be that at W > 4.0 the band width we chose is de-tuned from the optimal value by a bit, however we do not expect this to change the main results since the maxima blunt and the curves become flatter when disorder is increased.this section we present the distribution of local magnetization in the system with and without bubble (see Fig. 7).There we confirm that at low disorder the bubble has no effect beyond a small shift to "more Guassian" distributions.At disorder W = 4.0, which is supposed to be in the critical region, we see that the localized peaks at ⟨S z ⟩ = ±1 loose weight in favor of the a central "thermal" peak in the presence of the bubble.The latest suggest that an important fraction of samples thermalizes in the presence of the bubble.Unfortunately we could not quantify the fraction of thermalizing samples in a more rigorous manner.A similar but more detailed analysis of distributions was done in Ref. [68] where non-interacting localize bits were coupled to a ergodic bubble.Unfortunately finite size effects are stronger in our current setup, thus avoiding trustworthy conclusions from such distributions.

C. Additional data
In this section we present a further analysis on the purity γ i on localized XXZ chains with and without bubble.Although the system sizes are small and quantitative results are not conclusive, we find instructive to sketch the qualitative signals of thermalization due to the presence of the bubble.The black line corresponds to γXXZ in absence of the bubble (L = 14).The chain length L = 14 is fixed.β = {0.55,0.6, 0.58, 0.55, 0.5} is tuned to be the optimal value for each bubble size nGOE = {3, 4, 5, 6, 8} .In both panels, γ XXZ is computed using L + 2 chain sites compared to the shown system sizes.
In Fig. 8, we study the limit where the bubble size n GOE growth while the length of the chain is kept fixed.The most visible feature in these plots is that the purity γi decreases when a bubble is coupled to the chain, with bigger bubbles having a larger effect.This trend is obvious and expected, because coupling an ergodic bubble can, on average, only increase the thermality.Hence, this trend does not give us any non-trivial information on the nature of the ergodic-to-MBL transition.
To go beyond this, we note first that at very large distance i → ∞, one expects that γi − γ XXZ → 0 independently of whether the systems is localized or ergodic, simply because in a very large system, arbitrarily large bubbles will occur naturally.Numerically, at disorder strengths W < 7, we observe however that the difference γ XXZ − γi seems to decay very slowly, if at all.In Fig. 9, we remark that this function is actually reasonably fitted by an exponential decay law added to a nonzero constant δ 0 .For W ≥ 7, the deviation of δ 0 from 0 is no longer statistically meaningful.It appears hence that for disorder strengths W < 7, there is a effect of the bubble at large distances (see discussion in the main text).This effect is weak, for example 5% at W = 5.2, but; however due to the smallness of the systems considered we can not rule out attributing this signal to finite size effects.

D. Quantum avalanche theory
Let us recall the avalanche theory [65] for a disordered spin chain coupled to a thermal inclusion.If the chain is isolated from the bubble, it is assumed to be manybody localized, forming local integrals of motion (LIOM) [23,24].We assume that the LIOMs can be described by a localization length ξ, whose meaning is that a LIOM centered at some site of the chain is coupled with strength e −ℓ/ξ (operator norm of the coupling) to sites at distance ℓ.The ergodic bubble is assumed to be perfectly thermal and is described by a Hilbert space of dimension n GOE .
The logic of the theory is based on an iterative application of perturbation theory: We start with the ergodic bubble and couple in the first step only the LIOMs which have the largest overlap (and therefore coupling) with the bubble.Those are the LIOMs centered at the sites adjacent to the ergodic bubble.There are two of those sites because we consider a periodic chain.If the coupling is non-resonant, we declare that the adjacent LIOMs have not been thermalized, and the analysis stops there.If the coupling is resonant, then we declare the adjacent spins thermalized and we imagine that the ergodic region has now dimension 4n GOE instead of n GOE , i.e. the LIOMs have been absorbed by the ergodic bubble.This is now repeated, proceeding further and further away from the original bubble, until either the adjacent LIOMs are nonresonant, or until all LIOMs have been thermalized.
Let us make this quantitative.We assume that the bubble has thermalized all LIOM's within a distance ℓ, leading to a delocalized (ergodic) space with dimension d therm = n GOE 4 ℓ .We will now only keep the overall scaling of the relevant terms to obtain a criterion for thermalization: We couple the next LIOM with strength of order e −ℓ/ξ .We use random matrix theory to compute a typical matrix element of such a coupling and the typical level spacing is estimated simply as ∝ d −1 T .We then get for the quantity distinghuishing resonance from nonresonance: The thermalization proceeds until this ratio is of order 1, hence up to ℓ = ℓ c given by where ξ * is the critical localization length.This formula shows that ℓ c is finite if ξ < ξ * and infinite otherwise.Due to the final length of the system, the correct distance up to which LIOMS are thermalized is of course the minimum of ℓ c and L/2.The Hilbert space dimension of the thermal region is then, for ξ < ξ * ,

GOE
Taking also here into account the finite size, we arrive at the expression given in the main text.
Let us now briefly comment on the effective dimension which was mentioned in the main text.This concept was introduced in [68] as follows.The spectral form factor of local observables at site i have a width that is way smaller than of the order of local energy scales, as it is in a normal many-body system.Instead, this width is of order of the mean level spacing times the effective dimension d eff,i = e 2i/ξ d therm .This formula was also verified numerically in [68].It is the smallness of this width, or, equivalently, the fact that d eff,i is way smaller than the total Hilbert space dimension, that gives rise to the extremely slow dynamics in this system, as one goes far away from the ergodic bubble.
Finally, we return to equation ( 6), which can be viewed as a consequence of the above discussion.In the case (occuring with probability 1 − p) the ergodic bubble did not kick-start an avalanche, i.e. when the coupling to the first LIOM was too weak to thermalize it, we naively approximate the influence of the ergodic bubble on the purity as (γ i − 1/2) ≈ (1 − e −|i|/ξ )(γ XXZ−1/2 ).This simply reflects the fact that the purity at site i is influenced by the delocalized states in the ergodic bubble when |i|/ξ is not too large.If an avalanche was indeed started, we can use the theory described above.The purity is computed from local observables and hence it can be computed assuming that one has an ergodic system with dimension d eff,i .Some simple random matrix theory leads then to the conclusion that this purity is given by d −1/2 eff,i .

E. Purity correlations
We compute the purity perturbatively in the highdisorder limit.This will confirm the existence of a residual purity-purity correlation at large distance.

Perturbative expression for purity
We focus on the standard MBL model Let η be classical configurations η i = ±1 and let h i ∈ [−W, W ]. We assume that we are in the MBL phase and the eigenstates connected to the classical configurations.
We will actually do a perturbative expansion in J/W (locator expansion) that disregards resonant regions.It yields an expression for the eigenstates given by where A = i A i,i+1 is a sum of local Hermitian terms, determined by the condition Depending on η, there is either zero or one η ′ such that this matrix element is non-zero.We compute the unique non-zero matrix element as Now we compute +O((J/W ) 3 ).(14) It is easy to see that the first order term [A, σ z i ] does not contribute to the expectation value so we calculate the second order term as well, which will be of order (J/W ) 2 .This suggests that we should have considered also the second order in the expression for A, but this is fortunately not the case.This is because the term linear in the second order A will again be a commutator [A (2) , σ z i ], whose expectation value will vanish, just as it did for the first order A. The result is hence which can be computed to give 2 (16)

Figure 2 .
Figure2.Average gap ratio[5] r for the disordered field XXZ spin chain of size L coupled to a thermal bubble of size nGOE = 4 (left) and the same system without bubble and chain length L + 2 (right).Inset: Pair wise crossing of gap ratio r(W ) as function of 1/L for both systems.

Figure 3 .
Figure 3. Purity γi as function of the distance from the bubble for different disorder strengths W . Left (Right) panel shows system size L = 16 (L = 18) and nGOE = 4. Horizontal shaded areas correspond to the γXXZ in the absence of the bubble with L + 2 chain length compared to the bubble-chain system.

Figure 4 .
Figure 4. Purity γi at fixed distance i = 5 from the bubble plotted as function of system size L for different disorder strengths W . Left panel: XXZ spin chain of length L with bubble of size nGOE = 4. Right panel: XXZ Heisenberg chain of size L. Average is taken over disorder realizations and eigenstates.

Figure 5 .
Figure 5. Average purity correlation C(r) as function of the distance from the bubble for disorder strengths W = 4.4, 5.2, 6.0, 8.0.Dashed lines correspond to bubble-less system of size L + 2 while solid lines to bubble-full system of size L

Figure 6 .
Figure 6.Average gap ratio ⟨r⟩ as function of random matrix band width β for different nGOE = 3, 4, 6, 8.The black dashed line stands for β = 0.55, 0.6, 0.55, 0.5 respectively the chosen values in our simulations that cut off the plateus where the spin chain-bubble mixing is maximal for the relevant parameters.

Figure 7 .
Figure 7. Distribution of local magnetization for XXZ chain with (color lines) and without bubble (black lines) for system size L = 20 and L = 18 respectively (bubble size nGOE = 4).Each color denotes the spin at site i from the bubble (see Fig. 1 for exact layout)

Figure 8 .
Figure 8. Purity γi average over disorder realization and eigenstates as function of bubble size nGOE, several disorder strengths W and different distances i from the bubble.The black line corresponds to γXXZ in absence of the bubble (L = 14).The chain length L = 14 is fixed.β = {0.55,0.6, 0.58, 0.55, 0.5} is tuned to be the optimal value for each bubble size nGOE = {3, 4, 5, 6, 8} .

Figure 9 .
Figure 9. Upper panel: difference γ XXZ − γ i for different disorder strength, L = 18 and nGOE = 4. Dashed lines are fittings γ0 + Ae −i/ξ , with γ0 as the limiting value i → ∞ .Lower panel: γ0 as function of W for system sizes L = 16, 18.In both panels, γ XXZ is computed using L + 2 chain sites compared to the shown system sizes.