Avalanche induced coexisting localized and thermal regions in disordered chains

We investigate the stability of an Anderson localized chain to the inclusion of a single ﬁnite interacting thermal seed. This system models the effects of rare low-disorder regions on many-body localized chains. Above a threshold value of the mean localization length, the seed causes runaway thermalization in which a ﬁnite fraction of the orbitals are absorbed into a thermal bubble. This “partially avalanched” regime provides a simple example of a delocalized, nonergodic dynamical phase. We derive the hierarchy of length scales necessary for typical samples to exhibit the avalanche instability, and show that the required seed size diverges at the avalanche threshold. We introduce a dimensionless statistic that measures the effective size of the thermal bubble, and use it to numerically conﬁrm the predictions of avalanche theory in the Anderson chain at inﬁnite temperature.


I. INTRODUCTION
Sufficiently strong disorder localizes noninteracting particles in any dimension [1].Each single-particle orbital is then exponentially localized with an (energy-dependent) localization length.Rare low-disorder regions have no deleterious effects on such systems as the asymptotic decay of the associated orbitals is set by the disorder strength of the typical regions.
What of the effects of rare regions in interacting localized systems?Interacting rare regions are much better baths than their noninteracting cousins as their densities of states are exponentially larger [44][45][46][47].They seed thermal bubbles which can grow by thermalizing nearby l-bits.As the bubble grows, it becomes a more effective bath that can absorb more distant l-bits.If this process runs away, then the system undergoes a thermalization "avalanche." De Roeck and Huveneers (DRH) argued that avalanches destabilize MBL at any disorder strength in d > 1, so that MBL is an unstable dynamical phase [44].In d = 1 however, the instability is present only if the l-bit localization length exceeds a threshold value.Avalanche theory thus predicts the existence of stable MBL and thermalizing phases, and the properties of the phase transition between them [48].
Despite the far-reaching implications of the avalanche instability, it has been tested in relatively few studies.References [49,50] numerically found the instability in toy models defined in the l-bit basis, while Ref. [51] showed that a central spin can thermalize many l-bits at a vanishing interaction scale.However, a recent numerical study reports no evidence of avalanches in MBL spin chains [52].Furthermore, experiments with ultracold atoms suggest a surprising robustness of the MBL phase in two dimensions [39].Using two-component mixtures in an optical lattice with only one of the components experiencing a disorder potential, Ref. [39] reported persistent local memory in the disordered component despite interactions with the clean component.
In this paper, we test the avalanche instability for an Anderson chain coupled to a single interacting thermal seed (Fig. 2).We focus on the Anderson chain as the l-bits are uniquely defined and numerically accessible.The Anderson l-bits have a distribution of localization lengths, whereas DRH assumed that the MBL l-bits are characterized by a unique localization length [53].
Our first result is that introducing a single thermal seed into a system with a distribution of localization lengths will result in a partial avalanche.This occurs when the mean localization length exceeds the critical value When the system avalanches the resulting thermal bubble absorbs a noncontiguous finite fraction f∞ of the chain's l-bits, and coexists with the remaining localized l-bits.The partially avalanched Anderson-seed model thus provides an example of a delocalized nonergodic phase that is robust to all perturbations which preserve the free-fermion structure of the chain Hamiltonian.For the Anderson chain with box disorder, we compute f∞ as a function of the disorder strength h (Fig. 3).Next, for [ξ ] > ξ c (or equivalently for disorder strengths h < h c ), we identify the length scale up to which the thermal bubble must grow before the instability argument ensures FIG. 1. Phase diagram of avalanches.Density plot illustrating the fraction of samples which avalanche close to the threshold.Typical samples avalanche (region I) only when h < h c and the initial thermal seed is sufficiently large, or more precisely, when the hierarchy of scales Eq. ( 2) is observed.As h approaches h c from below, the necessary seed size for typical samples to avalanche diverges as |h − h c | −2 after preasymptotic scaling of |h − h c | −1 .For h < h c with small seeds, rare samples may still avalanche (region II).For h > h c no samples avalanche (region III).
that typical samples exhibit avalanches.As this length scale diverges as h → h c from below (with exponent ν = 2), the probability that a thermal seed of fixed finite size sets up a runaway avalanche vanishes in the same limit.Figure 1 summarizes the behavior of the ensemble averaged fraction of l-bits absorbed into the thermal bubble f as a function of h and the thermal seed size.
In more detail, we identify three length scales: l dir , l FS , and l .The scale l dir sets the number of l-bits with significant direct coupling to the bare thermal seed, and is proportional to the number of degrees of freedom in the thermal seed [Eq.(4)].The second scale l FS is the Harris-Luck scale; this sets the number of l-bits required to accurately determine whether the disorder strength h is above or below the threshold value h c [54][55][56][57].Thermal bubbles of spatial extent greater than the third scale l typically grow indefinitely if h can be determined accurately on the scale l , i.e., if l l FS .We present analytical arguments that typical samples avalanche if As l and l FS diverge near h = h c , while the seed size and thus l dir remain finite, typical samples partially avalanche only in region I of Fig. 1.
In contrast, though rare samples may avalanche in region II, in typical samples the thermal bubble is of finite extent.Third, we derive the conditions on the coupling strength between a single l-bit and a thermal seed which determine when the two systems are either very strongly or very weakly hybridized [Eqs.(30), (31)].On applying these conditions to the Anderson-seed model near the avalanche threshold, we obtain the remarkable result that there are approximately ten l-bits with intermediate values of couplings to the bare thermal seed.This provides conservative bounds on l dir : l S l dir l S + 9.9, where l S is proportional to the number of degrees of freedom in the thermal seed.This implies that most l-bits in numerically accessible finite-size chains are in the intermediate regime in which the extent of hybridization between the lbit and the seed is partial with large eigenstate-to-eigenstate variations (Appendix A).
As we find significant seed sizes to be necessary, the chain lengths L accessible to exact diagonalization are L ≈ 10 ∼ l dir .As Eq. ( 2) is not satisfied, it is difficult to conclusively identify avalanches.Furthermore, the partially avalanched Anderson-seed system is delocalized, but nonergodic.These phases are also notoriously hard to characterize numerically [58][59][60][61].
To address the challenge of identifying the nonergodic delocalized phase, we present our fourth result: a dimensionless statistic [v] which measures the size of a noncontiguous thermal bubble.This statistic characterizes the seed connected to the Anderson chain by its effectiveness as a bath for a probe l-bit.Specifically, we use a modified ratio of the matrix elements between eigenstates to their energy level spacing to determine the number of degrees of freedom in the total system that couples to the probe l-bit (previously introduced in Ref. [48]).The measure [v] grows linearly with L with slope f /ξ c if and only if a finite fraction f of the l-bits in the chain have been absorbed into a thermal bubble.We see robust evidence for linear growth for [ξ ] > ξ c in the Anderson-seed system.
As an aside, we also confirm that [v] detects the MBL transition in the interacting disordered fermionic chain (corresponding to the disordered Heisenberg model).We find that the finite-size estimate of the critical disorder strength agrees with that extracted from standard spectral measures, despite the inclusion of an external thermal seed (Sec.VII).
The outline of the paper is as follows.In Sec.II we introduce the model of interest, a disordered fermionic or spin chain coupled to a thermal seed.In Sec.III, we derive the conditions under which the thermal seed absorbs the first few l-bits.We subsequently generalize the theory of DRH in Sec.IV, and show that due to fluctuations on the localization length the avalanche instability only leads to at most partial thermalization.We then turn to the length scales that control the probability that a finite thermal seed can set up a runaway (partial) avalanche in Sec.V.After deriving further conditions on the Anderson-seed model to see the runaway avalanche instability in Sec.VI, we turn to finite-size numerics in Sec.VII.We define the dimensionless measure [v] and present numerical evidence in favor of partial avalanches in the Anderson-seed model.Finally, we discuss the implications of this work for MBL systems in Sec.VIII.

II. MODEL
The random-XXZ chain provides a canonical model of both single-particle and many-body localization: Here s ν n = 1 2 σ ν n are the usual spin-1 2 operators, the local fields are drawn independently from the box distribution h n ∈ [−h, h], and we use open boundary conditions.Under the Jordan-Wigner transformation, the XXZ chain maps to an interacting disordered fermionic model [62,63]: up to an additive constant.Above, and the fermionic creation and annihilation operators satisfy the usual anticommutation relations {a n , a † m } = δ mn .In this article, we focus on two values of in the random-XXZ chain: (1) The random-XX chain ( = 0), which maps to the non-interacting Anderson model Eq. ( 6) [1].
To probe the avalanche instability, we couple one end of the disordered chain to an external thermal seed with even Hilbert space dimension d 0 .The total Hamiltonian acts on the Hilbert space obtained by taking the tensor product of the Hilbert spaces of the thermal seed and the chain and is given by The Hamiltonian is depicted in Fig. 2(a) in the fermionic representation.We discuss the first and third term in Eq. ( 8) in turn below.The Hamiltonian of the thermal seed H T is given by a d 0 × d 0 random matrix of bandwidth 4W T drawn from the Gaussian orthogonal ensemble (GOE).We take d 0 to be an even integer; we explain why below.The eigenvalues w a of H T have mean [w a ] = 0 and variance [w 2 a ] = W 2 T .The square brackets [•] denote ensemble averaging with respect to the GOE matrix H T and the on-site disorder in the spin chain.We refer to a single instance of this ensemble as a sample.An important quantity in subsequent analysis is the density of states of the thermal seed at maximal entropy: For even d 0 , we define spin-1 2 operators on the thermal seed as follows.Let d 0 ∈ 2N.Such a thermal seed may be obtained by coupling a d 0 2 -level system to a spin-1 2 particle (with label T).The required seed spin- 1  2 operators are then s ν T for The decomposition of the thermal seed into a d 0 2 -level system and spin- 1  2 particle is useful to define the seed-chain interaction term V in Eq. ( 8) to have the same form as the term in H D that couples neighboring sites: See Fig. 2(a).

A. l-bit basis
In the Anderson model ( = 0), H D is bilinear in the fermionic operators, and hence easily diagonalized: where f † α = n φ α n a † n creates a fermion in the αth diagonal orbital, and the single-particle energies ˜ α and orbitals φ n α are the solutions to the equation The diagonal orbitals define the localized bits (or l-bits) of the disordered chain in the absence of coupling to the seed (V = 0) [26][27][28][29]31].The single-particle energies have mean and variance given by They are further bounded by the following inequality (which is saturated as L → ∞): The l-bit orbitals φ α n are exponentially localized [1] with localization length ξ α and localization centers nα We index the orbitals according to their distance from the thermal seed so that α < β ⇐⇒ nα < nβ .This indexing scheme implies that [|n α − α|] = O(ξ α ).
Rewriting the seed-chain coupling V in the diagonal basis: The αth l-bit is coupled to the seed with strength that decays exponentially with α: The Hamiltonian of the Anderson-seed system in the l-bit basis is shown in Fig. 2(b).The interacting fermion system with = 0 is more complex as H D cannot be simply diagonalized.Nevertheless, full localization persists above a nonzero critical disorder strength .This many-body localized (MBL) phase is characterized by an extensive set of quasilocal integrals of motion (the l-bits) which generalize the diagonal orbitals of the Anderson model to the interacting case [26][27][28][29][30][31][32][33].The expansion of H in terms of the l-bit operators contains higher-order terms as compared to the Anderson case: The higher-order terms in H D encode diagonal interaction energies between the l-bits.These interaction energies decay exponentially with the distance between the l-bits [11,31,33,64].We ignore these diagonal interaction energies henceforth.The higher-order terms in V allow several l-bits to reconfigure simultaneously through their interaction with the thermal seed.The typical strength of interaction terms that reconfigure the αth l-bit (i.e., terms J(n) α in V which do not commute with the l-bit occupation number f † α f α ) decays exponentially in the order n; however as there are exponentially many-in-n such terms they cannot be discarded by any naive argument.From the Heisenberg equation of motion we see that the energy scale associated with l-bit flip processes is given by Thus the effective coupling Jα between the αth l-bit is given by the norm of the commutator: The coupling Jα determines the interaction energy scale and may be understood as a resummation of the terms J(n) α .Thus, the qualitative features of the Hamiltonian in the l-bit basis in the interacting case is the same as in Fig. 2(b) with ξ α in Eq. ( 20) defining the effective localization length of the αth l-bit.
We note that the effective localization length defined through Eq. ( 20) may be longer than that controlling the decay of any single bare coupling, e.g., J (1) α in Eq. ( 18), as multiple interaction terms contribute to the decay.This discrepancy has been observed numerically in Ref. [33].

B. The localization length distribution in Anderson chains
The envelope of each single-particle orbital decays exponentially, with a characteristic localization length ξ α .This localization length varies as a smooth function of the single-particle energy as L → ∞, Thus, the distribution of localization lengths converges to a simple analytic form set by the single-particle density of states g( ): The distribution p(ξ ) is easily interpreted: consider the orbitals 1 < α < l which are centered within a region of length l 1 of the Anderson chain.The expected number of orbitals in this region with localization lengths ξ α ∈ [ξ, ξ + dξ ] is given by l p(ξ )dξ .

III. STARTING AN AVALANCHE: ABSORBING THE FIRST FEW L-BITS
In Ref. [44], DRH studied the effects of coupling a finite thermal seed to a localized chain.They assumed that the avalanche could start, that is, that the seed hybridized strongly with nearby l-bits and forms a small thermal bubble.They then derived conditions that the thermal bubble grows asymptotically, so that l-bits that are very distant from the seed eventually hybridize with the thermal bubble and reach thermal equilibrium.If the avalanche does not start, then the asymptotic thermalization instability is immaterial.
Section V discusses several length scales that control the probability to start an avalanche and their requisite hierarchy for a typical sample to avalanche.In this section, we have a more modest aim: we identify parameter regimes of H in which the bare thermal seed strongly hybridizes with a nonzero number of nearby l-bits l S due to the direct interactions in V : The above condition is necessary but not sufficient for starting an avalanche with high probability.We return to the sufficient conditions in Sec.VI.

A. Bandwidth condition: There exist resonances between the thermal seed and nearby l-bits
The first l-bit can flip only if it can exchange an amount of energy ˜ 1 with the thermal seed.For small J, we thus require the seed bandwidth 4W T to be larger than ˜ 1 .Using the bound on single-particle energies in Eq. ( 14), we obtain the condition in Eq. ( 28) for the Anderson model.When = 0, and indeed in general MBL models, there is no known general bound on |˜ α |.Nevertheless, at strong disorder strengths, we expect that, up to small corrections, The first relationship implies that the variance of the l-bit energies is set by the variance of Using a property of the box distribution that , we obtain the following relation at strong disorder: In summary, enforcing the condition that 4W In numerical experiments, this bound must be satisfied by a reasonable margin as the density of states of the thermal seed vanishes at the edge of its spectrum.Let the closest l-bit with index α = 1 be unoccupied, i.e., f 1 |c b = 0.At small V , the l-bit hybridizes with the seed if the typical matrix element V i j between |ψ i and a second eigenstate |ψ j = |w c ⊗ |c d exceeds the typical energy level spacing of the seed δ i j ∼ ρ −1 0 .Here δ i j = |E i − E j | is the separation of the energies of the state |ψ i , |ψ j .The typical matrix element is given by Hybridization requires that typically V i j ρ −1 0 .In the Heisenberg case the interaction in the l-bit basis is more highly structured.However the structure of matrix elements remains simple; specifically we take the matrix elements of local operators to follow the eigenstate thermalization hypothesis (ETH) [65][66][67].From the ETH it follows that the distribution of matrix elements is characterized by a single scale, here Jα , and Eq. ( 29) follows.
How much larger should V i j be as compared to ρ −1 0 in order for the l-bit to be strongly hybridized?This is tricky to quantify as the distribution of |V i j /δ i j | is broad even within a sample.In Appendix A, we use a second probe l-bit to define a statistic [v] (detailed in Sec.VII) which depends on the |V i j /δ i j |.Using [v] we determine that a thermal seed hybridizes strongly with a single l-bit if When Eq. ( 30) is satisfied by the first l S l-bits, the seed strongly hybridizes with these l-bits and the thermal bubble grows by l S l-bits.
Similarly, the typical configurations of the thermal seed and l-bit do not hybridize significantly if When Eq. ( 31) is satisfied by the first l-bit, the seed and lbit remain very weakly coupled.Consequently, the thermal bubble does not grow larger than the seed and the avalanche does not start.We observe that c W and c S are approximately two orders of magnitude apart.At values of Jα ρ 0 / √ d 0 between c W and c S , the l-bit is partially absorbed into the thermal bubble and there is significant eigenstate-to-eigenstate variation in the degree of absorption.As the numerically accessible chain lengths are short, most l-bits are in this intermediate coupling regime in our numerical study in Sec.VII.

IV. THEORY OF THE RUNAWAY INSTABILITY
We first review the avalanche theory of DRH [44] and derive the condition for the runaway thermalization instability for a unique l-bit localization length.We then extend the avalanche theory to the case of finite localization length distributions and apply it to the Anderson-seed model.We identify the critical disorder strength h c below which the instability is present, and discuss the properties of the partially avalanched regime.
In order to differentiate between complete and partial avalanches, we define f ∈ [0, 1] to be the fraction of l-bits that are a part of the thermal bubble when the avalanche ceases.In infinitely long chains, the fraction f depends on the disorder strength h and the number of l-bits that are directly thermalized by the seed l dir .
In this section, we assume that the avalanches start with probability 1.That is, we study the functional dependence of f on h in the following limit: We emphasize the order of limits L → ∞ first so that l dir L always.DRH predict that f∞ (h) is a step function of unit height.Our extended avalanche theory predicts a nontrivial function f∞ (h) for h < h c in the Anderson-seed model (horizontal axis of Fig. 1), with an intermediate critical value

A. Review of DRH avalanche theory
The setup in Ref. [44] is identical to that in Fig. 2(b) with all ξ α equal.Suppose the starting seed has strongly hybridized with the first (α − 1) l-bits with The thermal bubble then has Hilbert space dimension d = d 0 2 α−1 and density of states ρ = ρ 0 2 α−1 .That is, the thermal bubble has an exponentially larger density of states as compared to the bare seed.From Eq. ( 30), the l-bit with index α strongly hybridizes with this thermal bubble if its direct coupling is large enough: Using Eq. ( 17), we find The above condition is satisfied by a greater and greater margin with increasing α provided that the localization length of each l-bit is above the threshold value DRH thus predict complete thermalization à la the ETH with f∞ = 1 when Eq. ( 37) holds, and localization with f∞ = 0 otherwise.Numerical studies on toy models in the l-bit basis confirm this prediction [49,50].

Spectral function of the seed operator
The DRH avalanche theory assumes that the ETH holds in the thermal bubble at every stage of the avalanche process.Precisely, the theory requires that after α l-bits have been absorbed into the thermal bubble, the matrix elements of the seed operator are given by the ansatz where |b a , |b b are eigenstates of the thermal bubble with respective energies b a and b b corresponding to infinite temperature, is the spectral function of a † T in the thermal bubble, and η ab are random numbers with zero mean and unit variance.One might be concerned that this ansatz is inaccurate because the larger the α, the (exponentially) longer it takes for the αth l-bit to be absorbed into the thermal bubble, and thus the smaller the scale of structure in the spectral function F B (ω).Specifically, the αth l-bit decays at a rate set by Fermi's golden rule, implying that the peaks of the spectral function have width α .
The ETH ansatz in Eq. ( 38) is valid in a small energy window if this scale far exceeds level spacing, i.e., α ρ 1.As exponentially grows with α when ξ α > ξ c , Eq. ( 38) holds with increasing accuracy as the size of the thermal bubble increases in the avalanching regime.References [50] and [44] also discuss this "back-action" of the absorbed l-bits on the spectral function.As both studies conclude that more detailed treatments only provide subleading corrections to the DRH theory, we do not discuss the issue of back-action further.

B. Partial avalanches in the Anderson model
When the distribution p(ξ ) of l-bit localization lengths has finite width (as in the Anderson chain), only a fraction of the l-bits may have long enough localization lengths by Eq. ( 37) to be absorbed into the thermal bubble.The avalanche ceases when the l-bits that are not a part of the thermal bubble are too weakly coupled to the bubble to be absorbed.The stopping condition is that the inequality is satisfied by a fraction (1 − f∞ ) of the total number of the l-bits.We derive a self-consistency equation for f∞ below.The fraction f∞ is determined by the properties of the disorder distribution {h n } alone and is generally between 0 and 1 below a critical value of the disorder strength.
A few notational points before we proceed.Let ξ min (ξ max ) denote the minimum (maximum) localization length (i.e., p(ξ ) > 0 only if ξ ∈ [ξ min , ξ max ]).The mean localization length is given by where p(ξ ) is the ensemble distribution of localization lengths.Consider a small segment α/L ∈ [s, s + ds] of the disordered chain for s ∈ [0, 1] for α l dir .Let f ∞ (s) be the fraction of l-bits absorbed by the thermal bubble in this segment.This "local" fraction is related to the global fraction by integrating over s: In the limit of large L, it follows from Eqs. ( 17) and ( 41) that the l-bit with index α does not hybridize with the bubble if The local fraction of l-bits which satisfy the above condition is thus Integrating over s, we obtain an equation which determines f∞ implicitly: Solutions to Eq. ( 46) determine the fraction of l-bits absorbed into the thermal bubble.These solutions define three distinct regimes.We describe these regimes below, and delegate some mathematical details to Appendix B. (i) No avalanches for sufficiently localized chains.A trivial solution to Eq. ( 46) is f∞ = 0.The thermal bubble then fails to absorb a nonzero fraction of the l-bits in the chain.The stability of the f∞ = 0 solution determines whether or not the l-bits are asymptotically localized for arbitrarily large thermal seeds.For [ξ ] < ξ c , f∞ = 0 is the only solution and is stable: Thus, when [ξ ] < ξ c , there is no avalanche instability.The condition [ξ ] = ξ c determines the avalanche threshold.Below threshold h < h c ([ξ ] > ξ c ), f∞ = 0 remains a solution, but is unstable.
(ii) Complete avalanches for sufficiently delocalized chains.When the localization lengths ξ α are all larger than the threshold value ξ c , all l-bits are absorbed into the thermal bubble f∞ = 1 for ξ min > ξ c .(48) This is the DRH result of Sec.IV A.
(iii) Partial avalanches for generic distributions.When ξ min < ξ c < [ξ ], there is a stable solution to Eq. ( 46) with a nonzero f∞ satisfying f∞ fc := ξ c /ξ max : As f∞ < 1, the resulting avalanche is only partial.At late times, a noncontiguous thermal bubble comprising a fraction f∞ of the l-bits coexists with the remaining localized l-bits.This intermediate partially avalanched regime is generic: if p(ξ ) has any finite width as [ξ ] is tuned through the critical value ξ c , then Eq. ( 49) will be satisfied.Figure 3 shows the dependence of f∞ on the disorder strength h for a model with the same p(ξ ) as the Anderson model in Eq. (11).The p(ξ ) for the Anderson model with box disorder was obtained numerically; see Appendix C. Numerically, the disorder strength and thermal bubble fraction at threshold take the values h c ≈ 1.37 and fc ≈ 0.67. ( The inset of Fig. 3 schematically depicts the spatial structure of the thermal bubble below threshold. A number of comments are in order.First, Fig. 3 shows the unique stable solution to Eq. (46).Although f∞ = 0 is a solution at any h, its instability below threshold indicates that if an avalanche starts, it forms a thermal bubble with ∼ f∞ L fc L l-bits.Second, the f∞ vs h curve is discontinuous at threshold if and only if the disorder distribution has bounded support.That is, fc is nonzero if and only if ξ max < ∞, as FIG. 3. Thermal fraction of l-bits in the Anderson-seed model.The fraction f∞ of l-bits absorbed into a thermal bubble as a function of disorder strength h (solid green line).Insets schematically depict the resulting system.The leftmost red rectangle represents the thermal seed, while the circles represent the l-bits ordered by their distance from the seed.The red l-bits participate in the thermal bubble; the red lines remind the reader that thermalization is mediated by interactions with the seed.For h > h c , the thermal seed only strongly hybridizes with finitely many l-bits in its vicinity.For h < h c however, a thermalization avalanche forms a noncontiguous thermal bubble of ∼ f∞ L l-bits.
is the case for the Anderson model with box disorder.Third, the fraction of l-bits absorbed into the thermal bubble within a segment α/L ∈ [s, s + ds] [given by f ∞ (s)] can exhibit strong dependence on the distance from the seed.For example, Eq. ( 45) implies that f ∞ (s = 0) = 1 and f ∞ (s = 1) 1 for h < h c .Most dramatically, f ∞ (s = 1) → 0 as h → h − c ; i.e., the l-bits at the end of the chain are always localized at threshold.
The partially avalanched phase is thus an example of a delocalized nonergodic phase.A number of spectral and dynamical properties immediately follow from the simple picture of a thermal bubble with ∼ f∞ L l-bits coexisting with localized l-bits.These are summarized in Table I.We return to these spectral properties in Sec.VII, when we numerically test the predictions of the extended avalanche theory for an Anderson chain coupled to a thermal seed.

C. Avalanches in the random-Heisenberg model
The key difference between an MBL system and an Anderson chain coupled to a single thermal seed is that, in the MBL chain, for [ξ ] > ξ c , the avalanche instability leads to complete thermalization irrespective of the distribution of l-bit localization lengths.This is because an MBL system has a finite density of rare low-disorder regions of any given size, and thus a finite density of thermal seeds.As the coupling strength between any l-bit and the closest seed participating in the bubble Jα is finite as L → ∞, the left-hand side of the inequality is exponentially larger than the right-hand side and every l-bit is absorbed into the thermal bubble.Nevertheless, it may be that the distribution of localization lengths plays a role in analytical predictions about the MBL transition based on avalanches [48].The numerical study of Ref. [33] suggests that deep in the MBL phase, the lbits of the random-Heisenberg model are characterized by a unique localization length.Ascertaining whether the localization length is necessarily unique at intermediate strengths near the putative MBL transition is an interesting topic for future study.

V. AVALANCHES AT FINITE l dir
The number of l-bits that are directly thermalized by the bare thermal seed, l dir , controls the probability of an avalanche starting below threshold.If l dir is too small, then only rare samples will avalanche and f ≈ 0 irrespective of the value of [ξ ] (see Fig. 1).
In this section, we derive the minimum value of l dir such that typical samples avalanche below threshold and f (h, l dir ) ≈ f∞ (h) for h < h c .We find that the minimum value is set by the maximum of two length scales, l and l FS , so that l dir > max(l , l FS ): typical samples avalanche; l dir < max(l , l FS ): typically finite thermal bubble.
We show that l and l FS diverge with different exponents as h → h c .Consequently, the probability to start an avalanche with fixed l dir rapidly (exponentially) goes to zero as threshold is approached h → h − c , and f (h, l dir ) rapidly approaches zero as h → h − c [see Fig. 5(b)].In comparison, f∞ (h) jumps at h = h c (see Fig. 3).The different length scales l dir , l FS and l * are illustrated in Fig. 4, which shows a sample in which the thermal bubble fails before growing beyond l * .The resulting phase diagram is summarised in Fig. 1.
We note that multiple diverging length scales near a dynamical transition in a disordered system is reminiscent of infinite-randomness physics.We return to this point in the discussion section.
A. l dir : Locus of direct thermalization due to the seed L-bits with index α l dir are sufficiently strongly coupled to the seed so as to directly hybridize with it.Thus, l dir sets the scale to which the thermal bubble is guaranteed to grow.As the weak and strong coupling thresholds, c W and c S , are separated by over an order of magnitude, it is not possible to define this length precisely.However, by the arguments in Sec.IV, the length is bounded by the relation FIG. 4. Schematic avalanche.The thermal seed (red rectangle) and the red l-bits form a noncontiguous thermal bubble.The blue l-bits have short localization length and are weakly coupled to the thermal bubble.Direct coupling to the seed adds l dir l-bits to the thermal bubble (red shaded area).The boundary to the direct coupling region is spread over ≈9.9 sites.Below the Harris-Luck scale l FS , fluctuations in the window averaged disorder make it impossible to determine whether the system is above or below the avalanche threshold.The thermal bubble has to grow to a size comparable to l to cause runaway thermalization.In the figure, l dir < max(l , l FS ) and the thermal bubble is typically of finite extent.

Using |φ α
1 | ∼ e −α/ξ α and considering a typical sample (i.e., taking ξ α ≈ [ξ ] in the vicinity of the seed) we obtain where the length scale l S sets the range of l-bits strongly coupled to the seed, As the avalanche threshold is approached [ξ ] → ξ c the length scale l dir remains finite and specified by l S l dir l S + 9.9.

B. l : A length scale controlling runaway thermalization
Define p stop (l ) to be the probability that a thermal bubble with ∼ f∞ l l-bits is stable to the inclusion of more l-bits, i.e., that the thermal bubble "stops" at length l.For h < h c , this probability is exponentially decaying.We define this decay length to be l A thermal bubble of length l > l has a finite probability of absorbing O(L) l-bits.Below, we derive an expression for l , and discuss the implications of its divergence near threshold.Consider a partially avalanched Anderson chain below threshold (h < h c ) of length l 1.By Eq. ( 46), the thermal bubble has ∼ f∞ l l-bits.Suppose we increase the length of chain to some L l.If the (L − l ) l-bits in the extended part of the chain are too weakly coupled to the bubble, then the thermal bubble's size remains ∼ f∞ l.The condition for weak coupling is By taking Jα = 1 2 Jφ α 1 ≈ 1 2 Je −α/ξ α and rearranging, we find that the thermal bubble's size is self-consistently ∼ f∞ l if Here l W [Eq. ( 54)] sets the number of l-bits with strong or intermediate coupling to the seed.The probability that Eq. ( 57) is satisfied for all α > l is given by the product of the probabilities that each l-bit satisfies Eq. ( 57) [68].Thus, the stopping probability for a thermal bubble with ∼ f∞ l l-bits is The number of terms in the product which are less than unity [i.e., those with α < ξ max ( f∞ l/ξ c + l W /[ξ ])] scales with l, and hence the probability that the thermal bubble stops at length l is exponentially small for l > l [Eq.( 55)].
In Appendix B 1, we show that l diverges as the avalanche threshold is approached from the avalanching phase: Above, ∼ indicates asymptotic equality.The proportionality to |h − h c | −1 follows as [ξ ] is a continuous function of h with finite gradient at the avalanche threshold.The diverging length scale l is shown in green in Fig. 5(a).
For h < h c , the chain typically avalanches if l < l dir .The resulting thermal bubble has ∼ f∞ L l-bits and is described by the partial avalanche theory of Sec.IV B. This regime is depicted to the left of the dashed red line in Fig. 5

(b).
For h < h c and l > l dir , the thermal seed typically stops after absorbing O(l dir ) l-bits.This regime occurs between the red and blue dashed lines in Fig. 5(b).As runaway thermalization occurs in atypical samples where the thermal seed absorbs of order l l-bits, the mean fraction of l-bits in the thermal bubble satisfies the following inequality: As h → h − c , stopping probability in Eq. ( 58) becomes nondecaying with l and consequently f approaches zero [magenta line in Fig. 5(b)].
Let us turn to the effects of a finite chain length L. At finite L, there is a regime near threshold where l > L > l dir .In this regime, the stopping probability p stop is effectively constant over the length of the chain.The thermal bubble then typically grows to a size of order 1/| ln p stop | before stopping.Consequently, the fraction of the l-bits in the chain that typically participate in the thermal bubble decreases with L. Using the relationship between l and disorder strength near threshold [Eq.( 59)], See the cyan curve in Fig. 5(b).
In this argument we have treated the localization lengths of orbitals as iid drawn from the marginal distribution p(ξ ).In reality there are correlations: nearby levels repel in energy, and as the localization length is a function of the orbital energy only, there are corresponding correlations in the localization lengths of nearby levels.However these correlations decay exponentially on a characteristic length set by the localization lengths of the orbitals in question; as length scale is asymptotically smaller than l as threshold is approached, the short range correlations are irrelevant to the argument presented here.
Thus far, we have ignored fluctuations of the disorder strength on the scale of the thermal bubble; i.e., we have implicitly assumed that l, l l FS .This assumption is manifest in the appearance of the global mean localization length [ξ ] in the definition of l .Sufficiently close to threshold, the windowed average of the disorder shows significant spatial fluctuations, and the assumption that l, l l FS is violated.We derive l FS next.

C. l FS : The Harris-Luck scale
The Harris-Luck scale [54][55][56][57] l FS sets the minimum size of the window required to accurately estimate the average disorder strength in the chain.As the average disorder strength controls runaway thermalization (see Fig. 3), quantities such as f which can differentiate between the avalanched and localized phases cannot be accurately determined on length scales much smaller than l FS .
A precise definition of l FS is as follows.Consider the localization length averaged over the l-bits within a spatial window of length l: The windowed average ξl varies from one sample to the other (or equivalently, with the position of the window within an infinite sample).The deviation from the ensemble mean is set by the central limit theorem: where If the deviation of the windowed average from the ensemble mean is greater than the difference between [ξ ] and ξ c , no window averaged "order parameter" can consistently determine whether the chain is above or below threshold.l FS is defined as the threshold window size at which the deviation is comparable to Close to the avalanche threshold, l FS sets the minimum size of a thermal bubble that can cause runaway thermalization.This follows from the corrected version of Eq. ( 58) with [ξ ] replaced by ξal for l < l FS .Asymptotically, we then find where a is an O(1) number.By the central limit theorem where η l is a standard normal random variable [η l ] = 0, [η 2 l ] = 1.Substituting ( 67) into (66) we obtain When the second term in the exponent is comparable to the first term, the two terms may cancel and produce significant stopping probability.Such cessation of the avalanche due to such local fluctuations in the disorder cannot happen once the first term dominates: Close to threshold, when l l l FS , the growing thermal bubble will absorb typical regions.However, it will also encounter rare regions where the Anderson chain appears locally to be more localized.Such rare regions can stop the growth of the bubble, and hence stop the avalanche.This mode of failure due to rare regions becomes exponentially unlikely once the thermal bubble grows to a length scale l > l FS .As l FS diverges near threshold, the probability of an avalanche starting goes to zero (magenta curve in Fig. 5), and the thermal bubble size is typically l dir .
Finally, we emphasize that l FS controls finite-size scaling near the avalanche threshold for large enough L [54][55][56][57]. Figure 5 depicts a preasymptotic regime in which the finitesize corrections are set by l .We depict this preasymptotic regime as it is likely that numerically accessible chains are in this regime.

VI. OBSERVING THE AVALANCHE INSTABILITY
The requisite hierarchy of length scales to conclusively observe the avalanche instability in typical samples below the avalanche threshold is l , l FS l dir L. In this section, we identify parameter regimes of the Anderson-seed chain in Fig. 2 in which this hierarchy holds to test our extended avalanche theory.Numerically however, the chain lengths are so short that l dir is comparable to L.
These conditions supplement those of Sec.III which derived conditions on the thermal seed so that the first l S l-bit strongly hybridizes with the seed.That is, Sec.III asked that l dir > l S > 0.

A. L > l dir condition: Farthest l-bits are not thermalized by direct processes
To numerically determine whether l-bits are thermalized by avalanches or by direct coupling with the seed requires sufficiently long Anderson chains.That is, L should exceed l dir .Using the inequality derived in Sec.V A, we find Above, we assumed that the chain is in the avalanching regime ([ξ ] > ξ c ), and that l S 1 (see Sec. III B).

B. Partial avalanche condition: Dominant fraction of l-bits not at intermediate couplings
Localized l-bits coexist with a macroscopically large thermal bubble in the partially avalanched regime.The coexistence results in a bi-modal distribution of single orbital entanglement entropies.The single orbital entanglement entropies are obtained by tracing out all but one l-bit in each eigenstate.The bimodality quantifies the nonergodic nature of the partially avalanched phase (see Table I).
At finite L, the bimodal signature can be washed out by the l-bits that have intermediate coupling to the bubble and hence intermediate orbital entropy values.The condition for intermediate coupling is where c 0 = 1 2 Jρ 0 / √ d 0 .This condition follows from the discussion around Eq. ( 41).
As φ α 1 decays exponentially with α, the number of intermediately coupled l-bits is O(1), as compared to the O(L) number of localized l-bits.Nevertheless, at finite L, the former can exceed the latter.Below, we numerically compute the L at which the two are equal in the Anderson chain.
Figure 6(a) shows the disorder averaged distributions of The orbital index is shown in the legend.The violet region marks the intermediate coupling regime in Eq. (72).
Orbitals for which x falls to the right of the violet region in the thermal region (x > ln c S ) are strongly coupled to the thermal bubble and avalanche theory predicts that they will be absorbed into the thermal bubble.Those to the left in the localized region (x < ln c W ) are only perturbatively corrected by the thermal bubble.As the distributions of x widen and decrease in height with increasing L, the fraction of l-bits in the intermediate regime decrease with increasing L.
The disorder averaged total fraction of l-bits in the thermal, intermediate, and localized regions are plotted in Fig. 6(b).As expected, the fraction of intermediately coupled l-bits falls as O(L −1 ).However, the fraction of intermediate l-bits falls below the fraction of localized l-bits only at L ≈ 26.Thus, to see clear signatures of the predicted bimodality of single-site eigenstate entanglement entropy requires

C. Other conditions
As l dir could be as large as l S + 11 sites, the numerics described below with L ≈ 10 is likely in the regime l dir L. We therefore do not identify further parameter regimes based on the hierarchy between l dir and l or l FS .

VII. NUMERICAL EVIDENCE FOR PARTIAL AVALANCHES
We present numerical tests which support the partial avalanche theory in the Anderson-seed model of Eq. ( 8).Specifically, we find that below a threshold disorder strength (i) l-bits with exponentially small in L bare couplings to the seed have thermal eigenstate expectation values, (ii) the matrix-element-to-level-spacing ratio of an operator on the thermal seed is exponentially large in L, and (iii) the nearestneighbor level spacing ratio r flows toward the Poisson value r = 2 ln 2 − 1 with increasing L in all cases.These results are explained by the presence of a thermal bubble with ∼ f L < L l-bits.
We introduce a statistical indicator (denoted by v) in Sec.VII C which directly probes the effective size of the thermal bubble, and thus detects the nonergodic nature of the delocalization in the avalanched regime.It is also insensitive to the spatial inhomogeneity of the thermal bubble.In contrast, standard spectral signatures such as the half-cut entanglement entropy (at late times or in eigenstates) and the energy level spacing that detect violations of the ETH are difficult to interpret conclusively.The behavior of standard spectral signatures is discussed in Table I.Finally, we present numerical tests of the Heisenberg-seed model in Sec.VII D. As discussed in Sec.IV C, interacting fermionic chains have a finite density of low-disorder regions, and thus a finite density of thermal seeds.We therefore expect that (i) MBL is stable if and only if [ξ ] < ξ c , (ii) the avalanche instability when [ξ ] > ξ c leads to complete thermalization and an ETH phase, and (iii) the L → ∞ dynamical phase diagram is unaffected by the inclusion of a single finite thermal seed.Nevertheless, our numerical study in Sec.VII D has three aims.The first is to demonstrate that the v statistic introduced in Sec.VII C detects the MBL-thermal transition as a function of disorder strength.The second is to show that thermalization is complete in the delocalized phase by measuring the effective size of the thermal bubble to be ∼L.The third is to test whether the critical disorder strength shifts to larger values due to the inclusion of a finite thermal seed in the Heisenberg model at numerically accessible system sizes.We however find no evidence of such a shift within the accuracy of our numerics.

A. Parameter regime
We set J = 1, = 0, d 0 = 16, and W T = 0.8, and present data on either side of the avalanche threshold h c ≈ 1.37.These parameters satisfy the conditions laid out in Secs.III and VI.
Precisely, we satisfy the bandwidth condition Eq. ( 28), and the coupling strength condition Eq. ( 30), Above, we have used and the GOE density of states at maximum entropy ρ 0 = d 0 /(πW T ).We have numerically checked that the first l-bit is strongly hybridized with the thermal seed at these parameter values (data not shown).We note the chain length condition, Eq. ( 71), is L > 10.9 for the observation of avalanches.Furthermore, we expect significant bimodality in the distribution of eigenstate expectation values only for L 26 (see Sec. VI B).Consequently, for numerically accessible L, we indeed see no evidence for bimodality.

B. Thermalization of l-bits with exponentially weak coupling to the thermal seed
Extended avalanche theory (Sec.IV B) predicts that the thermal bubble has Hilbert space dimension d = d 0 2 f∞ L and density of states ρ = ρ 0 2 f∞ L for h sufficiently smaller than h c (so that the typical sample avalanches).Thus, an l-bit hybridizes with the bubble if the coupling to the seed Jα exceeds a threshold which is exponentially small in L: We test Eq.( 76) using the von Neumann entropy of the reduced density matrix obtained by tracing out all l-bits except α in each eigenstate, S α,i = tr( ρα,i ln ρα,i ), (77) where ρα,i = tr ᾱ (|ψ i ψ i |) is the reduced density matrix of the αth l-bit in an eigenstate |ψ i taken from the middle third of the energy spectrum of the Anderson-seed chain.The eigenstate averaged quantity Sα is defined by averaging S α,i over eigenstates taken from the middle third of the spectrum.Heuristically we expect Sα = 0 (localized l-bit), ln 2 (thermal l-bit). ( With each value of Sα , we also record the orbital overlap of the l-bit onto the first site of the chain φ α 1 .Our data thus consist of a list of pairs ( Sα , φ α 1 ) from all l-bits α and different samples.As we seek to measure the average coupling strength necessary to thermalize an l-bit, we further average Sα over a window of values of φ α 1 : where the normalization factor In the plots, we use r = 1.5.
In Fig. 7, we plot [S(φ α 1 )] as a function of φ α 1 (solid lines).The points (faded colors) are a random subsample of the data prior to window averaging [Eq.( 79)].Data are shown for chain lengths L = 3 . . . 10 (legend inset) for disorder strengths h = 0.7, 1.3, 1.9 which are respectively well below, close to, and well above the theoretically predicted avalanche threshold h c = 1.37.The spread on the unaveraged data is due to the sample-to-sample variation in the final size of the thermal bubble.At weak disorder, the leftward drift in the curve [S(φ α 1 )] with increasing L indicates the exponentially decreasing value of the direct coupling Jα ∼ φ α 1 necessary to thermalize an l-bit.This qualitatively confirms the avalanche FIG. 7. L-bit entanglement entropy in eigenstates as a function of the orbital matrix element.The eigenstate averaged l-bit entanglement entropy Sα is plotted as a function of |φ α 1 | (points), the orbital overlap which sets the coupling Jα to the thermal seed.Each point corresponds to a sample.The window average of Sα is additionally shown (solid lines).For h < h c (top panel), the leftward drift with increasing L (legend) indicates that the matrix element required to thermalize l-bits is exponentially decreasing with L. For h > h c (bottom panel), this drift saturates at sufficiently large L. For h ≈ h c (center panel), a small drift is visible.
prediction.This drift disappears at strong disorder where the instability is not present, and thus the thermal bubble is not growing with L.
To quantitatively study the drift with L identified in Fig. 7, for each L, h we extract the threshold value φ * at which [S(φ * )] = 0.7 ln 2. ( For φ α 1 φ * , the coupling is sufficiently small so that the l-bit is not fully absorbed into the thermal bubble.Equations ( 76) and ( 49) predict the following behavior for φ * as L → ∞: The numerically extracted values of φ * are plotted in Fig. 8.The figure confirms that φ * saturates at large h (h h c = 1.37) and exponentially decreases with L at smaller values of h (h h c = 1.37).We note that the rate of exponential decrease of φ * decreases with h for h < h c , and approaches fc near h = h c (dotted line), as predicted by extended avalanche theory.We also observe that the saturation value of φ * as L → ∞ increases with h above threshold.This is consistent with a finite thermal bubble whose size decreases with increasing disorder strength.
To confirm that the observed thermalization cannot be explained by direct hybridization with the bare seed, we note FIG. 8.The threshold coupling strength for l-bit thermalization vs L. At small disorder, the avalanche is almost complete and the threshold coupling strength decays as φ * ∼ 2 −L/2 (dashed line).For 0 < h < h c = 1.37, we see that φ * ∼ 2 −aL/2 with a continuously decreasing from 1 (dashed line) to fc (dotted line).In contrast, for h > h c , φ * saturates to a constant value as L → ∞.
that φ * = 0.016 for h = 0.7, L = 11.This value corresponds to a bare matrix element of level spacing ratio which is more than an order of magnitude below c S , the ratio necessary to induce thermalization in the absence of the thermal bubble, and approximately equal to c W = 0.01, where the effect of the bare coupling is typically perturbative.Near h = h c , it is unclear whether there is a runaway thermalization instability in finite-size numerics.As discussed in Sec.V, the stopping probability is significant for l dir < max(l , l FS ), and the probability of starting an avalanche rapidly goes to zero as h → h − c .Moreover, the fraction of l-bits absorbed into the thermal bubble will appear to be nonzero in a crossover region around h = h c at finite L (see Fig. 5).As the L range is limited, we do not attempt finite-size scaling here, and leave a detailed numerical study of the threshold to future work.

C. Exponential enhancement of the bath
Extracting φ * requires direct access to the l-bits in the chain.The method in Sec.VII B thus cannot be easily applied to interacting chains.In this section, we introduce a statistic which measures the thermal bubble size and is easily generalizable to interacting chains.
Consider the full Hamiltonian H [Eq. (8)] with eigenstates |ψ i and eigenenergies E i .Let us introduce a hypothetical probe spin with fixed energy splitting h P .This spin is coupled to the thermal seed by the same operators as the disordered chain, so that the part of the Hamiltonian involving the probe spin is The two states |ψ i ⊗ |↑ and |ψ j ⊗ |↓ are related by the matrix element and separated by the level spacing By the arguments in Sec.III B, the two states hybridize if V i j is much larger than δ i j .Thus, the state |ψ i ⊗ |↑ will hybridize with at least one other state if the quantity is sufficiently large.We take the logarithm in Eq. ( 86), to capture the typical value as the ratio V i j /δ i j is broadly distributed.Finally, we calculate the mean [v] by averaging over the ensemble of Hamiltonians H and over states i from the middle third of the spectrum.
What does [v] measure?We argue that as the statistic [v] thus probes the number of l-bits absorbed into the thermal bubble, Here [v 0 ] is the value of this statistic in the absence of the coupling between the thermal seed and the disordered chain (i.e., for V = 0).Suppose first that the full Hamiltonian H satisfies the ETH so that f∞ = 1.The ETH implies that ln |V i j /δ i j | ∼ 1 2 L ln 2 is maximal for j close to i. Consequently, [v] approaches 1  2 L ln 2 = L/ξ c as L → ∞.Next, suppose that the full Hamiltonian is many-body localized so that f∞ = 0. Here, ln |V i j /δ i j | is largest for states j at finite energy above or below i, so that [v] approaches a finite O(L 0 ) value as L → ∞.Finally, in a delocalized nonergodic phase with 0 < f∞ < 1, ln |V i j /δ i j | is maximal when the state j only involves reconfigurations of the thermal bubble.Such states are typically separated by exponentially many (in L) intervening states.These intervening states differ in the state of the thermal bubble as well as the localized l-bits .Assuming that the thermal bubble satisfies the ETH, we obtain Eq. ( 88).
The quantity exp(v i ) is identical to the quantity G i defined in Ref. [48].However, it has not been numerically studied before.The quantity v i is however distinct from the G statistic introduced in Ref. [14] because the latter parameter is defined with j = i + 1 rather than with max j in Eq. (86).As the matrix elements between nearby states of a delocalized nonergodic phase are typically exponentially small in L 2 , the G statistic studied in Ref. [14] will not detect partial avalanches.growth rates of the thermal bubbles.Indeed, the curves at h = 0.7, 0.9, and 1.1 appear to grow linearly with L with slopes between fc /ξ c and 1/ξ c .Above threshold (h > h c ), there is no runaway thermalization instability, and we observe that [v] indeed saturates.The crossover between linearly increasing and saturating behavior is approximately at h = h c .
In Fig. 10 we show spectral statistics for the Anderson-seed model with the same parameters as Fig. 9(a).Specifically, we plot the mean level spacing ratio r = [min(r i , r −1 i )] where and averaging is performed over the middle third of spectrum, and over the ensemble of Hamiltonians.r measures level repulsion and takes the values r = 0.531 . . . in an ergodic phase, and r = 0.386 in a localized phase [4,69].We see that, in the Anderson-seed model, for all disorder strengths, r flows toward the localized value with increasing chain length L. Level repulsion is absent in this model as even below threshold, where the model avalanches, .For all disorder strengths (legend inset) the flow is toward the localized value, with no indication of flow reversal.This is due to the coexistence of l-bits with the thermal bubble even below avalanche threshold.
there are many localized l-bits which coexist with the thermal bubble.

D. Heisenberg
We now turn to the numerical study of the Heisenberg-seed model using the v statistic introduced in Sec.VII C.

Parameters
We set J = 1, = 1, W T = 1.6, and d = 32.We collect data in the vicinity of the MBL transition in the random Heisenberg model at h = h MBL c ≈ 3.7 [5,[70][71][72].These parameters ensure that the bandwidth condition is satisfied, as is the condition to thermalize the first l-bit, J1 Here we have used that , and the GOE density of states at maximum entropy ρ 0 = d 0 /(πW T ).We numerically check that our choice of parameters leads to strong hybridization between the seed and the first physical bit (data not shown).
The length condition, Eq. ( 71), L > 10.9 is model independent and the same as in the Anderson-seed model.(86)] as a function of the chain length L for varying disorder strengths h (legend inset).The dashed line shows the theoretical upper bound of [v] = L/ξ c , when the entire chain satisfies the ETH.At large disorder, we see that [v] saturates, indicating that the interacting chain is stably localized.At low disorder, [v] grows with L; the slope shows a slight increase with L toward the ETH value (dashed line) for h < h MBL c ≈ 3.7.Our finite-size numerics is thus consistent with a delocalized ETH phase with a thermal bubble size ∼L.Remarkably, the value of h at which the scaling behavior of [v] changes is close to the critical value reported in previous studies h MBL c ≈ 3.7 [5,[70][71][72], and does not seem significantly altered by the presence of the seed.

Exponential enhancement of the bath
If avalanches underlie the transition out of the MBL phase, we expect the inclusion of the thermal seed to enhance the delocalization tendency at finite chain lengths, and thus shift the location of the transition to larger h.However, we do not see clear evidence of this shift at accessible system sizes.

VIII. DISCUSSION
We have extended and tested the theory of quantum avalanches in several respects.
(1) Assuming that a large enough thermal bubble forms, a single seed partially thermalizes a system with a distribution of l-bit localization lengths with [ξ ] < ξ c .The resulting avalanched phase violates the ETH and is delocalized but nonergodic.
(2) Thermal bubbles need to reach a certain size in a typical sample before they cause runway thermalization.The required bubble size diverges near the avalanche threshold and is set by the larger of the two lengths scales l and l FS (Fig. 1).
(3) For the Anderson chain coupled to a single thermal seed, we derived the critical disorder strength h c below which the chain partially avalanches, the fraction of l-bits absorbed into the thermal bubble for h < h c , and the system parameters to typically observe avalanches.
(4) We introduced a dimensionless measure [v] of the effective size of the thermal bubble that the seed is a part of, and numerically confirmed that the Anderson-seed model partially avalanches.
We describe a few broader implications of our results below.
The thermal bubble typically absorbs ∼l dir l-bits as h → h − c .However, rarely, it grows past the length scale l asy ∼ max(l FS , l ), and absorbs ∼ fc L l-bits.On blocking the chain into regions of size l dir , the probability of the bubble growing to size l asy is seen to be exponentially small in the ratio l asy /l dir .We therefore predict that the mean fraction of l-bits in the thermal bubble as Thus, the growth of the probability in region II of Fig. 1 with |h − h c | is slower than any power law.This low probability may explain the absence of avalanches reported in Ref. [52], as we discuss further below.Avalanches are suppressed near the threshold even if the localization length is unique (as is believed to be the case in MBL systems).A thermal bubble of size l now has significant stopping probability for l l because of the broad distribution of the couplings between the l-bits and the seed [33,73].Taking the log couplings to be independently normally distributed with mean [ln Jα ] ∼ α and variance Var(ln Jα ) ∼ α, we obtain the following expression following the steps in Sec.V: p stop (l ) ∼ e −(l/l asy ) 2 (92) with l asy ∼ max(l FS , l ).Although the exponent controlling the divergence of l ∼ |h − h c | −3/2 is larger than that in the Anderson-seed case, the asymptotic behavior of l asy sufficiently close to h = h c is unaffected as l asy ∼ l FS ∼ |h − h c | −2 .Following the discussion in the previous paragraph, the growth of the probability with |h − h c | in region II of Fig. 1 is slower than any power law in this case as well.
The suppression of avalanches near threshold explains a number of contradictory results in the literature.Reference [49] took the coupling between the l-bit α and the seed to be equal to J exp(−α/ξ ) and reported an ETH avalanched phase induced by a thermal seed of dimension d = 8 in an L = 12 chain.By neglecting the broad distribution of the couplings, we believe that Ref. [49] significantly overestimated the probability of realistic models to avalanche close to threshold [74].Reference [52] studied a realistic Heisenberg chain and reported that a thermal seed of dimension d = 8 has no destabilizing effect on the MBL phase.This result contradicts Fig. 1.As short Heisenberg chains typically do not generate their own large thermal seeds, Fig. 1 predicts a shift of the MBL-ETH transition to larger disorder strengths on coupling the chain to an external seed.The apparent contradiction is explained by the d = 8 seed, which our analysis suggests is too small to start avalanches in typical samples.
The transition between the localized and avalanched regime is reminiscent of the infinite-randomness transition [75][76][77][78].Specifically, for h < h c , the transition is characterized by two length scales which diverge with different exponents, l and l FS .The more divergent length scale controls the universal scaling function of the fraction of absorbed l-bits near the avalanche threshold: where F is specified by Eq. ( 91) and δh = |h − h c |.For uncorrelated disorder, l FS ∼ δh −2 while l ∼ δh −1 , and thus ν = 2.As in the infinite-randomness case, hyperuniform correlations in the disorder can reduce the exponent controlling the divergence of l FS while leaving the exponent of l unchanged [79].When the disorder is sufficiently hyperuniform that l is the most divergent length scale, the exponent is thus ν = 1.The quasiperiodic potentials used in previous studies provide access to this case [7,[80][81][82][83][84][85].Finally, we comment on the implications of this work for the MBL transition [19,[46][47][48]57,[70][71][72][73]80,[86][87][88][89][90][91][92][93].Numerically accessible systems of length L ∼ 10 near the MBL-ETH transition are likely to be in the intermediate coupling regime.Assuming that typical samples contain small thermal seeds, the effective coupling Jα between the αth l-bit and the seed takes values between the weak and strong limits for almost all α: In this regime, our single l-bit study suggests that l-bits partially hybridize with the seed with significant eigenstate-to-eigenstate variation within a single sample.
Remarkably, numerical studies of the MBL transition have reported significant eigenstate-to-eigenstate variation in entanglement entropy and other spectral quantities [80,91].
Interesting directions for future work include the quantitative description of the finite-size numerics in these terms and the modification of phenomenological renormalization group theories of the MBL-ETH transition [19,48,[86][87][88]90,92,93] to account for the intermediate coupling regime.Note added.Recently the authors became aware of a related forthcoming work by Colmenarez, Luitz, and De Roeck [94] that studies avalanches in MBL systems.We are grateful to W. De Roeck for discussing these preliminary results.A prominent feature of these plots is that the distributions of values across eigenstates is significantly wider than the shift induced by absorbing a single spin.Recalling that max j |V i j /δ i j | = e v i sets the extent of hybridization of the ith eigenstate, the width of these distributions is understood as the origin of the significant spread between the constants c W and c S .Specifically, the width of the distribution entails that for spins (or l-bits) with coupling strengths spanning about two orders of magnitude, there are significant numbers of resonant and nonresonant eigenstates.The peak in the standard deviation of S in the lower panel of Fig. 12 also reflects the wide distribution of the extent of hybridization of the spin across eigenstates.

APPENDIX B: VARIOUS PROOFS REGARDING PARTIAL AVALANCHES
In Sec.IV B, we introduced the self-consistency condition Eq. ( 46) for the fraction of thermalized l-bits in the disordered chain.In this Appendix, we derive several useful results which are stated in the main text.To study Eq. ( 46) in more detail, it FIG.14. Plots of ( f ).( f ) is plotted for various disorder strengths (see legend).For strong disorder ( f ) → f and for weak disorder ( f ) → f − 1 (dashed black lines).The solid black line shows = 0. ( f ) is asymptotically increasing, convex, and always satisfies (0) = 0.These properties are visually clear.From these properties various useful results regarding the roots ( f ) = 0 follow.These forms were computed numerically by exactly diagonalizing the Anderson model at L = 3000 and 4000 samples. is useful to define the quantity whose roots ( f ) = 0 are the solutions to the self-consistency equation (46).( f ) has several useful properties: (a) f = 0 is always a root.Specifically, lim f →0 + ( f ) = 0.
(b) ( f ) is differentiable for f 0; ( f ) has the first derivative given by This is seen by further differentiating to obtain This is easily seen by noting that the integrand in the second term on the right-hand side of Eq. (B1) is strictly non-negative over the domain of integration.
(e) ( f ) is continuous provided we make only the assumption that the cumulative distribution function ξ 0 dξ p(ξ ) is continuous.
These properties are all evident in Fig. 14.Using these properties, we prove the following propositions regarding the partial avalanche theory: Proposition 1.If [ξ ] < ξ c , then f∞ = 0 is the only solution to Eq. (46).
Proof.We prove the equivalent claim, that ( f ) has exactly one non-negative root, given by f = 0.It was previously noted that f = 0 is always a root.To see that this is the only root for [ξ ] < ξ c , we note that the first derivative of , evaluated at f = 0, is positive for [ξ ] < ξ c .As ( f ) is convex for f > 0, Eq. (B5) further implies there are no further non-negative roots.Proposition 2. If [ξ ] > ξ c , there is exactly one other solution f∞ to Eq. (46).
Proof.We prove the equivalent claim, that ( f ) has exactly two non-negative roots, one at f = 0, and one in f ∈ (0, 1].It was previously noted that f = 0 is always a root.The existence of the second root in the interval f ∈ (0, 1] follows from (i) the first derivative ∂ f | f =0 < 0, which implies that ( f ) < 0 for sufficiently small f > 0; (ii) (1) 0, which follows from ( f ) f − 1; (iii) the continuity of ( f ).That there are no further roots follows from the convexity of ( f ) for f > 0.

Evaluation of p stop
In this section, we derive the result quoted in Sec.V, that the probability p stop [Eq.( 58)] of a thermal bubble of size l ceasing to grow is exponentially small in l/l , where the Next, we replace the sum with an integral.The Riemann sum over α in the above equation is O(l 1 ), whereas the error incurred from replacing it with an integral is O(l 0 ).Thus, for sufficiently large l we are at liberty to make this replacement.We further neglect the subleading in l corrections to the upper bound of the integral over ξ .We then obtain p stop (l ) = exp   This sets the length scale to which the thermal bubble must grow to set off the thermalization instability.
As the avalanche threshold is approached, this length scale diverges.This can be seen by rewriting Eq. (B10), In the second step, we switched the order of integration, and performed one of the subsequent integrals.The higher-order terms may be neglected as they are asymptotically smaller than the leading order terms in the limit of interest, where 1/l → 0. Rearranging the above equation and substituting in Eq. ( 46) we then obtain This is the result given in the main text as Eq. ( 59).Sufficiently close to the threshold, we may expand in powers of h − h c , to obtain

APPENDIX C: NUMERICAL EVALUATION OF THE ANDERSON CRITICAL FIELD h c AND THRESHOLD THERMAL FRACTION fc
In this Appendix, we numerically calculate the critical disorder strength h c ≈ 1.37 and fc ≈ 0.67 [as quoted in the main text (50)] for the Anderson model

Calculation of h c
The mean localization length ξ max is a monotonically decreasing function of the disorder strength.From these distributions, [ξ ] is easily calculated to equal ξ c at h c = 1.37 . . .(see Fig. 16).

Calculation of the thermal fraction of l-bits f∞
To obtain f∞ as a function of h, we numerically solve for the roots of ( f ) in Eq. (B1).The form of f∞ extracted in this way is plotted in Fig. 3 in the main text (green solid line).

Calculation of the critical thermal fraction fc
As the avalanche threshold is approached from the avalanching regime, the fraction of thermal l-bits approaches the threshold value of fc : To obtain ξ max at the avalanche threshold, we use that the localization length is maximal in the center of the singleparticle spectrum (see, e.g., Fig. 17): Numerically, the average is taken over orbitals with energies in the window ˜ α ∈ [− J+h 100 , J+h 100 ], a range significantly smaller than that over which ξ ( ) varies (see Fig. 17).

FIG. 2 .
FIG. 2. The Anderson-seed model.An Anderson chain coupled to a thermal seed in (a) the physical (p-bit) basis, and (b) the diagonal (l-bit) basis.Under the Jordan-Wigner transformation, the random-XX chain coupled to a thermal seed maps to this model.

B.
Coupling strength condition: The thermal seed is sufficiently strongly coupled to nearby l-bits For V = 0 the eigenstates of H are product states of the thermal seed and a configuration of the l-bits in the disordered chain |ψ i = |w a ⊗ |c b .Here i = (a, b) where a = 1 . . .d enumerates the states of the bath and b = 1 . . . 2 L enumerates the states of the chain.Consider the fate of the product eigenstate |ψ i = |w a ⊗ |c b when interactions are reintroduced.

FIG. 5 .
FIG. 5. Length scales near the avalanche threshold.The upper plot shows the instability scale l (green), and the Harris-Luck scale l FS (yellow) diverging as δh → 0, while the length of the chain L (blue) and the locus of direct thermalization l dir (red) remain finite.The lower plot depicts f as a function of δh for fixed l dir in the thermodynamic (magenta) and finite chain (cyan).The mean thermodynamic fraction f deviates from fc when l becomes comparable to l dir and approaches zero as δh → 0 − .At finite values of L (cyan), fL saturates to an O(L −1 ) value when l > L or δh ∼ L −1 .This O(L −1 ) value decreases with increasing δh.At larger L however, the blue dashed line is to the right of the yellow dashed line, and fL saturates to an O(L −1 ) value when l FS > L.

FIG. 6 .√ 2 f L ) for |φ α 1 |
FIG. 6. Finite-size effect due to intermediate couplings between the l-bits and the thermal bubble.Panel (a) shows the distributions of the quantity x = log 2 (c 0 |φ α 1 | √ 2 f L ) for |φ α 1 | obtained from exact diagonalization of an open Anderson chain.Data are shown for different α (legend inset) and different chain lengths.L-bits with intermediate coupling to the thermal bubble fall in the violet region.(b) The total fraction of l-bits lying in the intermediate region (violet line) scales as L −1 .This fraction is less than the total localized fraction (blue line) at L ≈ 26.Parameters: h = h c , f = fc from Eq. (50).

Figure 9 (
a) shows [v] − [v 0 ] as a function of L for the Anderson-seed model.Here we take h P = 2πW T /d 0 .At sufficiently small disorder strengths (h < h c ), typical samples are predicted to avalanche, and the lines [v] − [v 0 ] = fc L/ξ c (black dashed line) and [v] − [v 0 ] = L/ξ c (black dotted line) respectively provide the lower and upper bounds for the

FIG. 9 .
FIG. 9. Scaling of [v] with L. Upper panel: The matrix element to level spacing ratio [v] ∼ f L/ξ c is plotted as a function of chain length L in the Anderson-seed model.Data are shown for various disorder strengths (legend inset) above and below the transition h c = 1.37.For a fully thermalizing model [v] ∼ L/ξ c (dashed), at threshold in the Anderson-seed model [v] ∼ fc L/ξ c (dotted), whereas below threshold [v] saturates.These different behaviors are reproduced for h far above and far below the transition.Lower panel: Analogous data are shown for the Heisenberg-seed model in the vicinity of the MBL transition (h MBL c ≈ 3.7).The crossover from growing to saturating behavior is consistent with this value.

FIG. 10 .
FIG. 10.Variation of level spacing ratio r with L in the Anderson-seed model.The mean level spacing ratio is plotted over the middle third of the spectrum for various disorder strengths in the Anderson-seed model.The thermal (GOE) value r = 0.531 . . .and the localized (Poisson) values r = 0.386 . . .are shown (horizontal dashed lines).For all disorder strengths (legend inset) the flow is toward the localized value, with no indication of flow reversal.This is due to the coexistence of l-bits with the thermal bubble even below avalanche threshold.

Figure 9 (
Figure 9(b) shows the ensemble and spectrally averaged mean [v] [Eq.(86)] as a function of the chain length L for varying disorder strengths h (legend inset).The dashed line shows the theoretical upper bound of [v] = L/ξ c , when the entire chain satisfies the ETH.At large disorder, we see that

FIG. 13 .
FIG. 13.Distributions of v i .The distributions of v i for different values of Jρ/ √ d (legend inset) on a linear scale (upper panel) and a logarithmic scale (lower panel).The widths of the distributions are significantly larger than the drift with increasing Jρ/ √ d.The exponential tail (dashed line) indicates the broadly distributed nature of the ratio max j |V i j /δ i j |.Parameters: d 0 = 2048, and those given in the caption of Fig. 11.

FIG. 15 .
FIG. 15.The distribution of rescaled localization lengths ξ/[ξ ] in the Anderson model.At large disorder strengths, the distribution is narrow and peaked about 1.However, with decreasing disorder, the distribution develops a double peak, with the taller peak moving to smaller values of ξ/[ξ ].The distributions were computed numerically by exactly diagonalizing the Anderson model at L = 3000 and 4000 samples.

ξ c / f∞ dx ln x 0
dξ p(ξ ) + O(l 0 ) .(B9)In the second line of the equation above we have simply made the replacement x = (ξ c α)/( f∞ l ).Thus, we see that for

FIG. 17 .
FIG. 17.Localization length as a function of single-particle energy ˜ .At any disorder strength (see legend), the localization length is maximal at ˜ = 0. We extract ξ max at each h by averaging the localization length over a small energy window around ˜ = 0 [see Eq. (C4)].Data for L = 3000 and 4000 samples.

l 1 ,
p stop ∼ e −l/l with l given by that (i) x ∈ [ξ c / f∞ , ξ max ], and (ii) f∞ → ξ c /ξ max as h → h − c .Thus, the argument of the logarithm is close to unity, and may be Taylor expanded:

1 =FIG. 18 .
FIG.18.Extracting the threshold thermal fraction fc .As the avalanche threshold is approached from the avalanching regime, f∞ → ξ c /ξ max .The plot shows the ratio ξ c /ξ max as a function of h.The vertical dashed line marks h = h c , at which we find fc ≈ 0.67.Data for 4000 samples.
At h = h c , ξ max equals the DRH value ξ c .The only ingredient necessary to find [ξ ] is the distribution of localization lengths p(ξ ).To obtain p(ξ ) we numerically obtain the eigenorbitals of the Anderson model (C1) with periodic boundary conditions.The localization length of each orbital is then determined by a least squares fit to the relationship − |r| ξ α + const.= ln L n=1 φ α n φ α n+r , (C2) where −L/2 < r < L/2.The distribution of values of ξ α obtained from diagonalizing many such Anderson Hamiltonians provides a numerical estimate of p(ξ ).Example distributions of p(ξ ) obtained this way are shown in Fig. 15.

TABLE I .
Properties of avalanched systems.Partially avalanched systems are delocalized, but nonergodic, and do not satisfy the predictions of the ETH.The table summarizes the behavior of various spectral properties in the localized, partially avalanched, and completely avalanched regimes.The symbol s th denotes the thermal entropy density at the appropriate energy density from which eigenstates are drawn.