A constructive theory of the numerically accessible many-body localized to thermal crossover

The many-body localised (MBL) to thermal crossover observed in exact diagonalisation studies remains poorly understood as the accessible system sizes are too small to be in an asymptotic scaling regime. We develop a model of the crossover in short 1D chains in which the MBL phase is destabilised by the formation of many-body resonances. The model reproduces several properties of the numerically observed crossover, including an apparent correlation length exponent $\nu=1$, exponential growth of the Thouless time with disorder strength, linear drift of the critical disorder strength with system size, scale-free resonances, apparent $1/\omega$ dependence of disorder-averaged spectral functions, and sub-thermal entanglement entropy of small subsystems. In the crossover, resonances induced by a local perturbation are rare at numerically accessible system sizes $L$ which are smaller than a \emph{resonance length} $\lambda$. For $L \gg \sqrt{\lambda}$, resonances typically overlap, and this model does not describe the asymptotic transition. The model further reproduces controversial numerical observations which Refs. [\v{S}untajs et al, 2019] and [Sels&Polkovnikov, 2020] claimed to be inconsistent with MBL. We thus argue that the numerics to date is consistent with a MBL phase in the thermodynamic limit.

However, theoretical descriptions and numerical observations of the MBL-thermal transition remain at odds with one another. Phenomenological models suggest that the transition has Kosterlitz-Thouless-type scaling [45][46][47], and occurs when the localised phase is destabilised by rare thermal regions which seed "thermalisation avalanches" [48][49][50][51][52][53][54]. Numerical studies, which are limited to small systems, do not find any evidence of rare thermal regions [55,56], but are known to be plagued by unexplained finite-size effects [57][58][59][60]. The absence of a theory of the finite-size crossover leaves unclear which features of the numerical data may survive in the thermodynamic limit, and has led Refs. [1,2] to claim that the numerical data precludes the possibility of an MBL phase altogether.
We develop a microscopically motivated resonance * philip.jd.crowley@gmail.com model for the one-dimensional MBL-thermal crossover at finite sizes. In this model the MBL phase is not destabilised by rare thermal regions, but by many-body resonances involving macroscopically distinct l-bit states. Although this mode of instability was previously identified [61] and observed in finite size numerics [62], it has received little attention in the literature. Specifically, we consider a presumptively many-body localised chain, analyse the statistics of resonances induced by local perturbations, and establish when these resonances destabilise MBL. The detailed analysis is different in the Floquet (Sec. II) and Hamiltonian (Sec. III) settings. However, in both cases, the same set of nontrivial length scales emerge which control the physics. The first of these is the bare localisation length ζ, which governs the exponential decay of off-diagonal matrix elements of local operators in the l-bit basis. A site-local perturbation introduces many-body resonances between eigenstates. The probability that a given eigenstate finds a first-order resonance involving l-bits within a range r (in the Floquet case) is given by Here, two additional lengths emerge: the correlation length ξ sets the typical range of resonances, while the resonance length λ determines their density. The RM predicts that ξ diverges as the localisation length approaches the critical value ζ c . This marks the transition between a localised phase in which the number of resonances is finite and a delocalised phase (dubbed thermal in Fig. 1a) in which the number of resonances grows exponentially with range. The finite-size behaviour near the transition depends crucially on the resonance length λ which is much larger than the lattice scale. For system size L λ (region I, Fig. 1a), typical eigenstates have no resonances and non-thermal expectation values. For  a) The resonance model (RM) predicts a continuous transition (orange point) between a localised (blue) and a thermal (red) phase, and an inverse correlation length |ξ| −1 (orange lines) that vanishes with exponent ν = 1 at the transition. In region I at system sizes smaller than a resonance length λ (purple), typical eigenstates have no resonances and spectrally averaged properties resemble those of the localised phase. b) The MBL-thermal finite-size crossover : At large L in the vicinity of the RM transition (hatched region), localisation is inconsistent due to overlapping resonances. The RM is however self-consistent in the blue regions. The RM thus describes the MBL-thermal crossover in small system numerics (horizontal line), even though it does not describe the asymptotic transition (black point).
system sizes L λ (region II), typical states participate in L/λ 1 resonances even at first-order [63]. The first-order analysis is clearly incomplete in regimes where the number of resonances induced by a single local perturbation grows with L (region II and thermal). In fact, the region of instability is somewhat larger if we consider locally perturbing the system at every site. In this case, a typical eigenstate develops a density ∼ ξ/λ of resonances each of which rearranges a region of size ξ (here and henceforth we measure lengths in units of the lattice constant). For ξ √ λ , the resonances typically spatially overlap and we expect them to lead to l-bit rearrangements on the scale of the system. The hatched region in Fig. 1b indicates the parameter regime and finite sizes where localisation in the RM is inconsistent due to this instability.
Nevertheless, we present analytical arguments in Sec. IV that the RM is self-consistent outside of the hatched region -i.e. at small enough L in region I and at any L for large enough disorder (i.e. 1/ζ). Rough estimates of the resonance length in Floquet and Hamiltonian disordered chains suggest 15 λ 50 for models numerically studied to date (see Sec. IV) . Thus, we believe that numerically accessible system sizes correspond to the horizontal dashed line in Fig. 1b, so that the observed crossovers in spectral quantities, spectral functions, finite-size drifts, etc. can all be predicted within the region of validity of the RM. Summarising the more detailed results in Sec. V, the RM reproduces many features of numerically exact data: • Localised region I: As typical eigenstates do not find a resonance for L √ λ , the RM predicts that region I displays the phenomenology of the localised phase: long-time local memory, a logarithmically growing light cone, sub-thermal eigenstate entanglement entropy of small sub-systems etc.. Spectrally averaged quantities are thus insensitive to the boundary between the MBL phase and region I (ξ = L, Fig. 1b), in agreement with Ref. [58].
• Drift of the critical disorder strength W c with L: The RM predicts the controversial observation of Refs. [1,2] that W c ∝ L at small L. . As the corrections are small (θ c 1), the RM explains the apparent 1/ω behaviour reported in Refs. [2,67].
• Scale-free resonances: Within regions I and II, q(r) is scale-invariant and resonances form at all ranges, in agreement with a numerically exact calculation of q(r) [62].
• Apparent sub-diffusion: On the thermal side of the transition (0 < −ξ < L), the dynamics at short times t < ω −1 ξ is critical. The RM describes this dynamics, and predicts a continuously varying exponent z in spectral functions ∼ 1/ω 1−1/z (see Figs 2a-b for Floquet, and Figs 2d-e for the Hamiltonian case). The RM thus explains the apparent sub-diffusion (as measured by z) reported in several studies [2,[67][68][69][70][71], without invoking rare region effects, which Ref. [55] finds are absent in numerically accessible systems.
• Exponential increase of Thouless time at weak disorder W W c : This numerical observation of Refs. [1,2] follows from the logarithmic growth of the light cone until time t ≈ ω −1 ξ in the thermal phase of the RM.
As the resonance model of the finite-size crossover assumes the existence of MBL, and reproduces the numerical observations of Refs [1,2], we conclude to their contrary, that the numerics to date appears consistent with a stable MBL in the thermodynamic limit.
We additionally predict three interesting features of the dynamical phase diagram that could be tested numerically in the near future.
• The exponents controlling the strongest lowfrequency divergence of [S(ω)] ∼ 1/ω 1−θc in region I: We predict that the exponent θ c is a nonzero non-universal value in the Floquet setting, while θ c → 0 + (corresponding to log corrections) in the Hamiltonian setting with energy conservation. That is, the existence and number of conservation laws affects the scaling theory of the finitesize MBL-thermal crossover.
• An empirical criterion for MBL: In localised systems, the distribution (v) of matrix elements of a local operator V that couple eigenstates in two small non-overlapping mid-spectrum energy (or quasi-energy) windows takes the form, with 0 < θ 0 ≤ 1 (see Fig. 2c). A simple numerical criterion follows: with ρ denoting the mid-spectrum many-body density of states. This criterion generalises the avalanche stability criterion of Ref. [48] to a setting without l-bits or rare thermalising regions.
• Detecting the crossover between MBL and region I: In region I, scale free resonances form, but remain rare. Thus eigenstate averaged observables are largely insensitive to the formation of resonances. However, by analysing the distribution of an observable over eigenstates, or conditioning on the formation of resonances, it is possible to numerically detect the crossover between MBL and region I. Such an analysis is performed in Ref. [62].
We proceed as follows. In Section II, we describe the Floquet resonance model, couple the RM to a probe spin, compute the statistics of many-body resonances that a reference l-bit state is involved in, and thus derive the disorder-averaged spectral function of a local operator. In Sec. III we repeat the analysis for a Hamiltonian system. In Sec. IV we establish the regime in which the RM is self-consistent, showing it to apply to small and strongly disordered systems (small L in region I in Fig. 1). In Sec. V we discuss the implications of this analysis for interpreting finite-size numerical data, before concluding in Sec. VI.

II. FLOQUET RESONANCE MODEL
After a brief overview of the set-up of the Floquet RM (Fig. 3) and the definition of the localisation length ζ, we detail a careful counting of resonances induced by a probe spin in Sec. II B. Panels (a), (b) and (f) in Fig. 2 summarise the results for the spectral function of the probe spin in the Floquet RM.
Resonances do not span the system for 1/ζ > 1/ζ c := log 2; this is the MBL phase of the Floquet RM. The RM MBL phase has infinite time memory of initial conditions, and a power-law divergence of the spectral function at small frequency (53).
The point 1/ζ = 1/ζ c marks the transition out of the RM MBL phase, at which resonances occur on all length scales. The statistics of the strongest resonances determine the low-frequency scaling of [S(ω)] in regions I and II within Fig. 1a. The exponent θ characterising the lowfrequency divergence of [S(ω)] in region II jumps at the transition (57).
Although typical states find increasingly many resonances at long ranges for 1/ζ < 1/ζ c , they remain rare on the scale of the correlation length ξ. Consequently, the RM predicts the behaviour of [S(ω)] at intermediate frequencies (59) in the thermal phase.
A. Set-up

Chain Hamiltonian
Consider a generic strongly disordered and interacting quantum spin chain with periodic boundary conditions, and subject to a periodic (Floquet) drive. For example, the Heisenberg model with random O(3) fields: where W , J and Ω set the disorder strength, interaction strength and fundamental frequency of the drive respectively, σ n = (σ x n , σ y n , σ z n ) is the usual vector of Pauli matrices acting on the nth site, and σ L+1 = σ 1 enforces periodic boundary conditions. The v n are independent and identically distributed (iid) random vectors with zero mean [v n ] = 0 and unit variance [v n · v n ] = 1, with, for example, iid Gaussian distributed entries. Here [·] denotes disorder averaging.
We assume two key properties of H(t): (i) it has no global conservation laws, and (ii) for some finite Ω, W J, the model is Floquet many-body localised, as per Ref. [72]. The specific form of H(t) is otherwise unimportant.
The dynamics of the chain is characterised by the Flo- In the vicinity of the RM transition, the correlation length |ξ| sets the cross-over frequency scale ω ξ ∼ exp(−1/|θ0|). The low-frequency behaviour (ω ω ξ ) is determined by the phase, while the intermediate frequency behaviour ω ω ξ J −1 is determined by the transition. The two other frequency scales are set by the system size: the Heisenberg scale ωH is the inverse level spacing, while ωc is the scale of the smallest off-diagonal matrix elements. The thermal-region I crossover occurs when ω ξ ∼ ωH ∼ ωc. In region I, only the exponent controlling the ω > ω ξ decay is visible. This exponent is continuously varying and is significantly corrected from its value at the transition in region I (as quantified by the O(θ0) term). The smallest value of the exponent is however set by θc. Panel (c): The exponent θ0 may be directly extracted from (v), the distribution of off-diagonal matrix elements of a local operator. In the localised phase, there are exponentially many off-diagonal matrix elements which are exponentially small in range, so (v) diverges as a power-law at small v. The exponent defines θ0. Panel (f ): The time averaged correlator [Czz] serves as an order parameter for the MBL phase.
[Czz] goes to zero smoothly as 1/ζ → 1/ζc is approached from the MBL side, faster than any power law in both the Hamiltonian and Floquet cases. quet operator where T = 2π/Ω and T is the usual time ordering operator. The associated Floquet states | a , and quasienergies a are defined by

Localisation in the l-bit basis
At sufficiently strong disorder in the MBL phase, we assume that the Floquet states | a may be identified with configurations of quasi-local integrals of motion, or l-bits [7,38,39] (in Sec. IV A, we discuss how this assumption may be relaxed). Each l-bit τ z n is traceless tr (τ z n ) = 0, squares to the identity (τ z n ) 2 = 1, is exponentially localised around the physical site n, and commutes with the Floquet operator Each Floquet state | a can be identified with an l-bit configuration τ a ∈ {−1, 1} L . The scalar element τ an = ±1 of τ a specifies the state of the nth l-bit: A quasi-local operator U diagonalises the Floquet unitary, and maps the physical spin operators to l-bits, Thus the σ α n are similarly exponentially localised operators in the l-bit basis.

Probe spin
Probe-chain coupling Probe spin Probe-chain coupling Figure 3. Set-up in the physical and l-bit bases respectively: a) A "probe" spin-1 2 (orange) couples to a strongly disordered chain (blue) at the site n = 0 (magenta). b) Transforming to the l-bit basis renders the Floquet unitary of the chain diagonal and the probe-chain coupling quasi-local. The coupling strength decays exponentially with distance from n = 0.
Consider two eigenstates | a , | b . We say two states differ at range r ab if the furthest flipped l-bit is at distance r ab from the site n = 0.
r ab := max{|n| : τ an = τ bn } The range is depicted in Fig. 4b. If the matrix element V ab := a |V | b of an operator V is non-zero, then V ab is also said to have range r ab . The length scale on which a physical spin operator is localised in the l-bit basis defines the localisation length ζ. Consider a local operator V acting on the physical site of index n = 0. The operator V can be decomposed into a sum of terms of increasing range where all the non-zero matrix elements of V r have range r. The asymptotic decay of the norm of V r defines ζ: We use the re-scaled Frobenius norm as it is simple to calculate analytically, and captures the typical expectation value of an arbitrary vector | ψ|V r |ψ | ≈ |V r |.

Coupling a probe spin to the disordered chain
To probe the dynamical phase of the disordered chain, we introduce a probe spin-1 2 σ P subject to a z-field of strength W . The combined Hamiltonian of the probe spin and disordered chain, is periodic with fundamental frequency Ω. Here H 0 encodes the part of the Hamiltonian in which the probe spin and disordered chain are decoupled and H 1 (t) encodes their coupling. Throughout we use cursive letters to denote properties of the combined Hilbert space of the disordered chain and the probe spin, and roman letters to denote properties of the reduced Hilbert spaces. The spin and chain are coupled an inter- Here V is some local operator which acts only on the n = 0 site of the chain, and which is assumed to have norm |V | = J, e.g. V = Jσ x 0 . The Floquet operator of the combined system is given by where U 0 is the Floquet unitary for H 1 = 0, and U 1 encodes the interaction Each eigenstate of the unperturbed Floquet unitary U 0 |ε 0 α = e −iε 0 α T |ε 0 α is a tensor product of a quasi-energy state of the disordered chain | a and a z-polarised state of the probe spin |σ , where α = (a, σ) is a composite label.
B. Spectral function of σ z P in the RM MBL phase ζ < ζc Our aim is to calculate the disorder averaged infinite temperature zz spin correlator, in the RM. Here the normalization by D, the Hilbert space dimension, ensures that [C zz (0)] = 1. For simplicity, we restrict to stroboscopic observations at the drive Figure 4. Organising resonances by range: a) The many-body spectrum of H0 in a small quasi-energy window is divided into two sectors labelled by the state of the probe spin σ =↑, ↓. | a ↑ labels a specific reference state. b) The l-bit configuration corresponding to reference state (red spectral line) is shown. The states | d ↓ in the opposite sector (green lines) can be grouped according to their range r from the reference state (ranges r = 0, 1, 2 shown); states at range r differ only on the l-bits with index |n| ≤ r (highlighted in orange). A state | d ↓ at range r is resonant with | a ↑ if its quasi-energy separation is less than the matrix element size v(r) (i.e. if it lies within the magenta region). In the plot, the first resonance occurs at range r = 2.
period t ∈ T N. The Heisenberg operator σ z P (t) at integer periods is given by The spectral function [S(ω)] is obtained by Fourier transformation of (21), The basic steps in the calculation are as follows. We resolve the trace in the correlator (21) over the eigenstates | a σ of H 0 , and argue in Sec. II B 1 that each term is well approximated by either unity or a pure tone: cos |V ab |t (resonance) (24) Above, |V ab | is the largest matrix element that couples | a σ to a resonant state | bσ whereσ is the opposite z-spin projection as compared to σ. Taking the matrix elements at range r to have a characteristic scale v(r), we obtain where p(r) is the probability (upon varying the initial state, and disorder realisation) that the resonant process with the largest matrix element is at range r, and is the probability of no resonances. As p(r) and v(r) are exponentially decaying in r, we find that the spectral Figure 5. Cartoon of approximate eigenstates: For the purposes of calculating the spectral function, the resonant eigenstates may replaced with cat states. Here the resonance is of range r = 2, so that only l-bits with indices n ∈ {−2, −1, 0, 1, 2} (red box) are reconfigured.
function is a power law at low frequencies, The exponent θ approaches zero as ζ → ζ − c from the localised side, but jumps to a non-zero θ c precisely at the critical point ζ = ζ c . Ref. [61] gave a similar resonance counting argument for the low frequency properties of the spectral function in the localised phase.
We now detail how these results are obtained. The final expressions for the spin-spin correlator are given in Secs. II B 4, II B 5.

Contribution of a resonance to the spectral function
Let us define a resonance. Consider a Floquet state |ε α of combined system, Expanding these Floquet states to leading order in V , we obtain where α = (a, ↑) [73]. We define the two states | a ↑ and | b ↓ to be resonant if the first-order correction is large, that is, if If g ba < 1 for all b, then we approximate |ε α by the unperturbed eigenstate | a ↑ .
If g ba > 1 for a single b, then degenerate perturbation theory yields 'cat' Floquet states to good approximation (Fig. 5). These states are depicted in Fig. 5. The two cat states (31) are split in quasi-energy by the matrix element |V ba |, Ignoring the sub-leading corrections, we thus obtain The corresponding contribution to the spectral function is two delta function peaks at ω = ±|V ba |. The absence of weight at zero frequency is a consequence of the equal amplitudes in the RHS of (31). We argue in Appendix A that extending this calculation to include a small nonzero weight at ω = 0 does not alter the low frequency behaviour of the disorder-averaged spectral function.
If g ba > 1 for multiple indices b, the eigenstates do not have the simple form in (31). Nevertheless, we argue in Appendix A that the strongest resonance, corresponding to the largest matrix element, sets the frequency of oscillation if ζ < ζ c . That is, In other words, for an initial state | a ↑ , the probe spin oscillates at a frequency ω a↑ for a window of time t ω −1 a↑ , and thus the Fourier transform of (34) is sharply peaked at ±ω a↑ . Analogous expressions for an initial state in the down sector are easily obtained.
2. The probability q(r) of resonance at range r We take all the matrix elements at range r to have a single characteristic value v(r) which is a monotonically decreasing function of the range r. This recasts the problem of finding the resonance with the largest matrix element as the problem of finding the resonance with the smallest range r. We now calculate v(r), and subsequently the probability q(r) of finding a resonance at range r.
As described in Sec. II A, V = L/2 r=0 V r may be decomposed into terms of increasing range r in the MBL phase. V r couples a given state | a to N r other states | b at range r, where The characteristic scale v(r) of each matrix element is determined by, where the correlation length ξ is defined by The omission of the unimportant pre-factor of 3/2 makes (38) approximate.
Two properties of ξ are noteworthy. First, ξ has the interpretation of a length only in the MBL phase of the RM, in which it is positive. Second, ξ diverges as ζ → ζ − c . When we use results of the RM to discuss the shorttime dynamics as ζ → ζ + c , we will be careful to use the absolute value of ξ.
Let ρ(r)dr denote the density of states per unit quasienergy with range in the interval [r, r + dr]; from here on we will coarse grain and treat the range r as a continuous variable. As the states are uniformly distributed in quasienergy b ∈ [0, Ω], and the total number of states within range r is given by 2 2r+1 we have Consider the dn = Ωρ(r)dr states with ranges in the interval [r, r +dr]. As they are uniformly distributed over the quasi-energy interval [0, Ω], the probability that an arbitrarily selected one of them has a quasi-energy in the interval ε 0 , and is thus resonant, is given by 2v(r)/Ω. It follows that the probability that at least one of these states is resonant with | a ↑ is given by where higher-order corrections in v(r)/Ω can be dropped for q(r) 1. Combining (38), (40) and (41) q(r) = e −r/ξ λ (42) with ξ as in (39) and the resonance length λ defined as We expect that λ 1 as MBL in the RM requires Ω J. Put another way, deep in the MBL phase where ξ 1, the probe spin will typically induce resonances of range r = 0 (i.e involving only the l-bit n = 0 to which it is directly coupled). For stable MBL, the probability of such resonances q(0) = 1/λ should be small so that nearest-neighbour resonances are atypical.
3. The probability p(r) that the strongest resonance is at range r The fraction F (r) of states that have not resonated up to range r satisfies the differential equation with solution The probability p(r)dr that the strongest resonance with the largest matrix element has range in the interval [r, r+ dr] is then determined by

The time domain correlator [Czz(t)] and the logarithmically growing light cone front
We now have all the pieces in place to write down the spin-spin correlation function. The strongest resonance for each state is mediated by a matrix element of size v(r) with probability p(r). Plugging this into the pure tone ansatz (34), and treating the disorder average as simply sampling the distribution p(r), we obtain the Floquet RM spectral function where the integral runs over all possible ranges 0 ≤ r ≤ L/2, and the infinite time average is simply the probability that a state of the uncoupled system is not resonant with any other state. We were unable to exactly perform the integral (47). However a crude approximation allows us to extract the asymptotic behaviour in the time domain. Specifically where the Iverson bracket takes values [[P ]] = 1, (0) when the proposition P = true, (false). Within this approximation we obtain where r(t) is obtained by solving v(r)t = 1, where we have defined The position r(t) has a simple interpretation as the front of a logarithmically growing light cone. Only the cat states formed from l-bits states with r < r(t) contribute to the correlation function at time t.

The spectral function [S(ω)]
From (47) it is straightforward to obtain the spectral function. For brevity we first recast the matrix element (38) as v(r) = Je −r/(θ0ξ) using (51). Then by inverse Fourier transform of (47) Inserting the calculated form of p(r) (46) into (52) yields in the MBL phase of the Floquet RM. The cutoff scale ω c is set by the smallest matrix elements at distance L/2, where the Heisenberg frequency ω H := Ω2 −L is set by the typical many-body level spacing. The high-frequency (ω ≈ J) behaviour of [S(ω)] depends on the microscopic Hamiltonian in the immediate vicinity of the spin, and is thus non-universal. In contrast, the exponent θ characterising the power-law at low frequency: is a consequence of distant resonances which reconfigure large regions of the chain. Thus, as ζ → ζ − c , we expect θ to have a universal functional dependence on |ζ − ζ c |.
For L λ (region II in Fig. 1a), it follows from (53) That is, θ vanishes linearly with |ζ − ζ c | as ζ → ζ − c , but jumps to a non-universal non-zero value at the transition.
For L λ (region I in Fig. 1a), θ = θ c + O(θ 0 ), so that the exponent is continuously varying. The low-frequency divergence in [S(ω)] is strongest when θ = θ c , we return to this in Sec. V B [74].
Eq. (56) implies that that disorder-averaged correlators exhibit a power-law decay at long times t J −1 in the RM MBL phase: The decay persists until time ∼ ω −1 c , which is exponentially larger than the Heisenberg time ∼ ω −1 H . The dynamics at these long time scales are due to the exponentially small (in L) fraction of cat states involving re-configurations of l-bits on the scale of the system size L.
A fraction of the eigenstates | a σ do not hybridise with any other states despite the coupling with the probe spin to the chain, even as L → ∞. As the probe spin has a well-defined orientation in these states (even upon including perturbative corrections), these states contribute to the infinite-time memory [C zz ] of the MBL phase.
We defer more detailed discussion of the finite-size behaviour of [S(ω)] to Sec. V.
C. Spectral function of σ z P in the RM thermal phase ζ < ζc In the thermal phase, we expect that the off-diagonal matrix elements obey the eigenstate thermalization hypothesis. In particular, the off diagonal matrix elements they do not decay exponentially with range r at large r, as assumed by the RM in (12). Consequently, the RM does not apply in this regime.
Despite being generally inapplicable, the early time predictions of the RM are found to hold even in the thermal regime. Specifically, as the probability of resonance q(r) is small for r |ξ|, [S(ω)] exhibits power-law decay (as in (56)) for J ω > ω ξ where, That is, the correlator's dynamics are critical until a time-scale ∼ ω −1 ξ . This result is obtained exactly as in the MBL case, with the refinement that, instead of working in a basis of l-bits (which do not exist in the thermal regime), it is necessary to work in a basis of "almost-lbits"τ z n . These operators have the same properties as l-bits (mutually commuting exponentially localised etc.), but only "almost commute" with the Hamiltonian

III. HAMILTONIAN RESONANCE MODEL
We describe the computation of the spectral function of the RM with Hamiltonian dynamics. Despite the Hamiltonian case appearing superficially simpler than the Floquet case (as it lacks the additional "ingredient" of a drive frequency) the analysis is more complicated due to the conservation of energy. The associated hydrodynamic mode constrains the late time dynamics, and hence the low frequency behaviour of the spectral function.
For simplicity, we assume that the chain has a single hydrodynamic mode. The analysis is easily generalised to accommodate further conservation laws, such as the spin conservation present in the "standard model of MBL" the Heisenberg model with random z-fields.
A. Set-up

Chain Hamiltonian
Consider a strongly disordered static chain with disorder strength W and interaction strength J. For specificity, consider the Ω → ∞ limit of the Floquet model in (4), that is, the Heisenberg model with O(3) random fields As before, the details of this model will be unimportant except for two key properties: (i) energy is the only conserved extensive quantity at any W, J, and (ii) the model is many-body localised for some finite W J.

The local energy a
In addition to its energy eigenvalue E a , each eigenstate |E a of H can be assigned a local energy a (r) which can loosely be understood as the expectation value of the Hamiltonian restricted to the sites n ∈ [−r, r]: Here H [−r,r] is the Hamiltonian (61) with the summation restricted to terms acting on the sites n ∈ [−r, r].
We make this notion sharp with the following definition where the energy shift E 0 (a, r) is obtained by averaging the energies of the 2 2r+1 states within range r of |E a The local energy has two useful properties. First, for two states |E a , |E b within range r, energy differences are preserved exactly Second, given a state |E a , the distribution of the local energies b (r) of the states within range r is Gaussian and centred at = 0. Specifically, where ∼ denotes convergence in distribution at large r.
Neglecting sub-leading corrections in J/W , the width of the Gaussian is given by

Coupling a probe spin to the disordered chain
The Hamiltonian of the chain coupled to a probe spin is given by The eigenvectors of H , H 0 and H are denoted |E α , |E 0 α and |E a respectively. These vectors play roles in direct analogy with |ε α , |ε 0 α and | a from the Floquet case in Sec. II. The eigenvectors and corresponding eigenvalues of H and H 0 are related by Each eigenstate |E a , σ of H 0 is assigned a local energy e (a,σ) (r) = a (r) + σh/2.
B. Spectral function of σ z P in the RM MBL phase ζ < ζc Our aim is to calculate the disorder averaged infinite temperature spin-spin correlator for time evolution generated by the Hamiltonian As in Sec. II B, states with resonant partners contribute a pure tone, while states with no resonant partners contribute unity (see (24)), and hence [S(ω)] follows.
The key difference between the Floquet and Hamiltonian cases stems from the energy dependence of the density of states at range r. In the Floquet case, at sufficiently large range r, the density of states at range r is independent of quasi-energy, thus all states states have an equal probability of finding a resonance at range r. In contrast, in the energy conserving case, states with unusually high/low local energy e α (r) couple to an atypically small density of states at range r. As such these atypical states find resonances at a significantly lower rate (see Fig. 6). We thus adapt the calculation to keep track of the local energy e α (r) of the states. This leads to a slower decay of F (r), and hence a slower than power law decay of correlations.

Identifying resonances
Recall the resonance condition: two states |E a ↑ and |E b ↓ that differ at range r are said to be resonant if Using (65), this condition is recast as 2. The probability q(e, r) of finding a resonance at range r, and local energy e Define q ↑ (e, r)| e=e (a,↑) (r) , the probability that a state |E a ↑ with finds a resonant partner state |E b ↓ at range r. Analogous to the Floquet case, q ↑ (e, r) is given by where ρ ↓ (e, r) is the density of states in the down sector (i.e. the opposite spin sector) at local energy e = e (b,↓) (r) and range r, and v(e, r), the characteristic size of matrix elements, coupling states from the two spin sectors at local energies e, and range r.
Consider the characteristic matrix element v(e, r). To begin with, we neglect the energy dependence of v and assume that the matrix element have the same form as in (38), v(e, r) = Je −r/ξ−2r/ζc .
We later discuss refinements to this approximation. Next, the density of states ρ σ (e, r) follows from (66), and the variance s e (r) is set by (67). Differentiating (78) and taking the asymptotically dominant behaviour we obtain (80) Equivalently stated, the asymptotic behaviour of ρ σ (e, r) is dictated by the growth-diffusion equation where the boundary condition is obtained by matching the solutions with (80).
3. The probability p(r) that the strongest resonance is at range r The growth diffusion equation (81), which describes the total density of states at local energy e and range r, is easily modified to describe the density of states which have not found a resonant partner by range r. At each range r, the hybridisation probability is set by q σ (e, r).
We thus obtain: Here the superscript 'u' (for unhybridised) distinguishes ρ u σ from the total density of states ρ σ . We now extract the probability p(r) that a state |E a σ finds its strongest resonance at a range r. Observe that the second term in (84) leads to exponential growth with r. Define a distribution that scales out this exponential growth: Substituting in (84), we obtain The substitution (85) has a simple interpretation: is the fraction of states which have not hybridised by range r. Eq. (86) is invariant under the replacements (e, σ) → (−e, −σ), by this symmetry F (r) is independent of σ. It follows that the probability p(r)dr that the strongest resonance of a given state is in the interval [r, r + dr] is given by Eq. (88) is the generalisation of the Floquet result (46) to the energy conserving case. Here it is necessary to solve the two-dimensional partial differential equation (86) rather than the simpler one-dimensional ordinary differential equation (44). What do the solutions of (86) and (88) look like? We discuss two regimes. The first regime in Sec. III B 4 is most relevant for the numerically accessible MBLthermal crossover in Fig. 1b. The second regime of L, |ξ| λ determines properties of the Hamiltonian RM in the vicinity of ζ = ζ c as L → ∞ and is discussed in Appendix B.
4. Far from criticality |ξ| < λ, or small critical systems L < λ < |ξ| Neglecting the energy dependence of q σ (e, r), Substituting (89) into (88), we obtain an approximate equation for F (r), which we denote as F 1 (r) to distinguish it from a true solution to the growth diffusion equations (86) and (87). Let us justify the approximation above a posteriori. For ξ L, the solution F 1 (r) of (90) decays exponentially on the length scale set by λ. Thus for r < λ, the bulk of the weight of the distribution of unhybridised states f σ (e, r) is at typical energies |e| < s e (r), where the energy dependence of q σ (e, r) can be neglected by making the replacement q σ (e, r) → q σ (0, r) in (88) to obtain (90). The approximation is thus valid for small critical systems L < λ < |ξ| (region I of Fig 1). Far from the crossover on the MBL side |ξ| < λ, few resonances form after the length scale ξ and f σ (e, r) does not becomes small at typical energies |e| < s e (r). The bulk of the weight of the distribution of unhybridised states f σ (e, r) is thus at typical energies and the approximation is justified.
On longer length scales r λ at 1/ξ = 0, the weight of f σ (e, r) at typical energies is depleted by the exponential decay. The weight of the distribution is instead concentrated at atypical energies |e| > s e (r) where the resonance probability q σ (e, r) is much smaller. Appendix B discusses the behaviour at r λ in detail. The solution to the approximated equation (90) is where Erf(·) and Erfi(·) are the usual error function and imaginary error function respectively. The correlator then immediately follows Using (88) and (52), we obtain the desired result: The spectral function exhibits the same ω −1+θ0 low frequency behaviour as (53) in the Hamiltonian RM MBL phase and at intermediate frequencies in the thermal phase. However, as the localisation length approaches the critical value ζ → ζ c , the correlation length diverges 1/ξ → 0, the correlation decay exponent θ 0 → 0 + , and the correction to the low-frequency ω −1 behaviour of the spectral function is logarithmic rather than power law. We further discuss the logarithmic corrections in Sec. V B 2.

IV. REGIME OF SELF CONSISTENCY OF THE RESONANCE MODEL
The RM assumes a characteristic range-dependence for the matrix elements v(r) of a local operator V acting at site n = 0 (see (38)). The coupling to the probe spin induces hybridisation between the eigenstates of H 0 . The reader might thus worry that the off-diagonal matrix elements of a local operator between the hybridised eigenstates is not consistent with the RM assumption in (38). In other words, the distribution of matrix elements after having introduced the probe spin is inconsistent with the distribution we assumed at the beginning.
We address this question in two parts. First, we show that [S(ω)] ∼ ω −1+θ0 at low frequencies even if the matrix elements at range r have a generic distribution p(v|r), as opposed to a single value v(r), so long as the aggregate distribution of off-diagonal matrix elements is distributed as a power-law in v at small v. Thus, we can relax the assumption in (38) to allow for a pre-existing population of resonant cat pairs states, as the matrix elements between such cat pairs and the reference state can differ from v(r). Next, we imagine perturbing a MBL RM chain, with a given p(v|r), weakly at every site. The local perturbations induce local resonances. When these resonances do not overlap, we argue that the distribution p(v|r) is unaffected at large r, and thus that the perturbed chain presents the same statistics of off-diagonal matrix elements v as the unperturbed chain at small v. Consequently, the exponent θ 0 that sets the low-frequency divergence of [S(ω)] is stable to local perturbations.
Specifically, we argue that the resonance model is perturbatively stable, and consequently our conclusions hold, in the regime in which resonances do not typically overlap. Eq. (95) holds deep in the RM MBL phase as L → ∞ and in region I (see Fig. 1) for sufficiently small systems. Three important conclusions follow: 1. As the RM is self-consistent deep in the MBL phase, the RM predicts and describes a stable MBL phase in the thermodynamic limit.
2. The RM describes the MBL-thermal crossover in short chains, despite being an inapplicable at large L.
3. The RM describes dynamics in the MBL-thermal crossover at short times as L → ∞, or equivalently on frequency scales: A. Generalised RM with p(v|r) Define the aggregated distribution of off diagonal matrix elements (v) as the distribution of matrix elements |V ba | that couple two narrow energy windows E a ∈ [E, E + ∆] and E b ∈ [E , E + ∆] at maximum entropy: where (v) and the distribution of matrix elements p(v|r) are related by (94). In Secs. II and III, we took the matrix elements at range r to be single valued p(v |r) = δ(v − v(r)). In the Floquet case the corresponding aggregated distribution of off diagonal matrix elements at small v is (36), 0 < θ 0 ≤ 1 is defined in (51), and the power law is obtained by coarse-graining over the scale separating the delta functions.
Eq. (53) follows from (98), independent of the precise model p(v|r) for the matrix elements at range r. Consider the Floquet RM. A change of variables in (44) yields The solution is the fraction of states which do not have a resonance induced by a matrix element of size v or larger. Note that F (v = ∞) = 1. Similarly we may define so that p(v)dv is the fraction of eigenstates of H 0 whose strongest resonance is due to a matrix element in the range [v, v + dv]. The spectral function is then given by, Substituting (98), we recover the previously calculated spectral function (53). The calculation presented in Sec. III for the Hamiltonian RM can be similarly generalised.
Note that a general model for the matrix elements alters the simple relationship between the localisation length ζ and the exponent θ 0 , and thus leads to an altered critical value of the localisation length ζ c := ζ| θ0=0 .

B. Self-consistent and stable localisation
To be self-consistent, the RM must have the same statistical distribution of resonances before and after a local perturbation.
Consider a perturbation V of strength |V | ≈ J applied at a single site n = 0 (as in Sec. (II B)). The effect of this perturbation is straightforward: first the eigenstate energies are corrected by the diagonal elements of V (i.e. E a → E a + V aa ) and second, each state |E a finds a resonance at range r (i.e. |V ab | > |E a − E b |, where r ab = r) with probability q(r) = e −r/ξ /λ. This leads to a pair of resonant 'cat' states with corresponding energies E a , E b and splitting |E a − E b | ≈ |V ab |. We now apply a second perturbation U , also of strength |U | ≈ J, at a site m a finite distance from n = 0. Naively, the arguments of Sec. (II B) imply each such subsequent perturbation causes more long range resonances to develop. However, this is not the case. The matrix element E a |U |E b ≈ Je −s/ξ where s = max(0, m − r ab ) acts to disentangle cat state pairs (103) whose splitting is small |V ab | Je −s/ξ . This removes all resonances due to V which are of long range r ab m/2. This disentangling of resonances is counterbalanced by the formation of new long range resonances due to the combined action of U and V . Their distribution is statistically identical to that induced by a single local perturbation. Specifically, the range of typical resonances remains O(ξ).
Short range (r ab m/2) resonances induced by V survive the second perturbation. When the surviving resonances overlap with those induced separately by U , the eigenstate entanglement further increases. Specifically, two cat pairs |E a , |E b (103) and |E c , |E d with respective level splittings |V ab | and |V cd | survive if In the state |E a , a small subsystem in the vicinity of n = 0 has entanglement entropy S ≈ 2 log 2. Similarly two "cats of cats" |E a , and |E e may be hybridised by a third perturbation W to form |E a ≈ (|E a + |E e )/ √ 2 , with entropy S ≈ 4 log 2. Here we have illustrated the increase of entanglement entropy due to overlapping resonances for the case The general case is more complex. However we suspect similar increases of the entanglement entropy when resonances overlap.
The merging of local resonances into larger resonant clusters with larger entanglement entropies represents an instability of the "l-bits + local resonances" picture assumed by the RM unless the localisation length is sufficiently short ξ √ λ . Consider perturbing the RM at every site. At each site, the probability of inducing at least one resonance between the reference state |E a and a second state |E b is 1 − F (r = ∞) ≈ ξ/λ. If the typical spacing between these resonances λ/ξ exceeds their typical size ξ, then they remain spatially separated. We conclude that for ξ 2 /λ 1 resonances do not merge, and do not alter the asymptotic distribution of matrix elements at low frequencies. The RM is thus self consistent and stable to local perturbations in this regime. This case is depicted in Fig. 7a where the extent of each resonance is indicated by the red arcs. We note that rare states participate in long range resonances r ξ; however these do not destabilise the localisation.
Finally, we note that the RM describes dynamics in the thermodynamically large thermalising phase at short times, or equivalently at frequencies satisfying (96). At these short times, resonances are rare and thus the RM is controlled. As noted in Sec. II C, the derivation of [S(ω)] proceeds through "almost-l-bits" that almost commute with the Hamiltonian.

V. RM PREDICTIONS FOR FINITE-SIZE NUMERICS
The RM is self-consistent in short chains in region I and provides a simple model for the MBLthermal crossover. Could the RM describe the numerically accessible MBL-thermal finite size crossover? A naive estimate of the resonance length λ comes from Eqs. (83) and (43) using numerical and experimentally reported values for the critical frequency or critical disorder strength [6,64,75]. This gives 15 λ 50. Physically, λ has to far exceed the lattice scale, as q(0) = 1/λ is the probability of a nearest neighbour resonance in the MBL phase. We thus reason that numerically accessible chain lengths L are smaller than or comparable to 2 √ λ , and that the RM is an analytically tractable model for the numerics.
In what follows, we describe several properties of the RM in short chains that explain numerical observations about the finite-size MBL-thermal crossover. The crossover occurs around the line |ξ| = L/2 separating the thermal phase from region I in Fig. 1a). We also explain the numerical observations of Refs. [1] and [2] within the RM. As the RM has a stable MBL phase, we weigh in on the controversy of the existence of MBL in favour of MBL.

A. Correlation length exponent ν = 1
The thermal-MBL crossover in the resonance model is characterised by a correlation length |ξ|: which diverges with exponent ν = 1. This value is close to the numerically reported values of 0.77 ≤ ν ≤ 1.02 reported for data collapses of different quantities in Ref [64]. Note that the RM exponent, as well as the numerically reported ones, violate the Harris bound for randomly disordered systems ν ≥ 2 [57,65,66], as they only capture the pre-asymptotic in L scaling.

B. Apparent 1/ω divergence of the spectral function
The RM predicts a power-law divergence in [S(ω)] at low frequencies in the MBL phase and in region I: Above ∼ indicates asymptotic equality up to a constant factor, and θ > 0. Deep in MBL phase, the following hierarchy of frequency scales hold: and [S(ω)] takes the form in (107) for ω > ω c with the exponent θ given by θ 0 > 1 in (57). In region I in Fig. (1), |ξ| L/2, and the frequency scales are arranged as: Below, we show that the low-frequency divergence of [S(ω)] is strongest in the middle of region I and is given by [S(ω)] ∝ ω −1 up to logarithmic corrections. Ref. [2] interpreted the apparent ω −1 behaviour as inconsistent with MBL. The RM however predicts this behaviour near the finite-size MBL-thermal crossover in region I and allows for a stable MBL phase.

Floquet systems
The exponent θ 0 in (57) vanishes as |ξ| → ∞ in the RM. The strongest low-frequency divergence [S(ω)] is however not ∼ 1/ω (indeed, as noted in [2] such a strong divergence would violate an elementary sum rule) because the exponential term in (53) modifies the exponent. The RM instead predicts the following spectral function in the middle of region I: with θ c = ζ c /2λ, as given by (57).
As λ 1 and ζ c is on the lattice scale, we conclude θ c = ζ c /2λ 1. The strongest low-frequency divergence in (110) is thus close to 1/ω.
Note that (110) implies a power law decay of correlations at late times. Such decay can only be consistent with a logarithmically spreading light cone (50) in the absence of any conserved quantities, such as in a Floquet system.

Hamiltonian systems
Hamiltonian systems conserve energy, which results in a logarithmic, rather than power law, correction to 1/ω scaling of [S(ω)]. Specifically, for |ξ| λ L, we simplify (93) to obtain: Here ∼ indicates equivalence up to an ω independent prefactor. Observe that this decay is not asymptotically consistent with hydrodynamics. The light-cone only grows logarithmically in time in the RM (see Fig. 8), but (111) implies critical correlations that decay faster than 1/ log(Jt) as t → ∞, r ∝ log(Jt) Figure 8. Logarithmically growing light cone: the Heisenberg operator σ z P (t) (in (22)) is localised to the probe spin site time t = 0. Under time evolution, the support spreads and defines a light cone. After a time t, this light-cone has width r(t) ∝ log Jt (green).
More careful analysis of the Hamiltonian resonance model finds that below a frequency timescale ω λ := v(λ), the decay of [C zz (t)] is dictated by a form consistent with hydrodynamics. We note this corresponds to a time averaged value which goes to zero as [C zz ] ∼ ξ −1/2 . However, as (113) applies outside of the regime of self-consistency of the resonance model, we relegate further discussion to Appendix B.

C. Localised finite-size crossover
As the resonance probability is small for L 2 √ λ , the RM predicts a localised finite-size crossover (i.e. a localised region I).
First, the time-averaged correlator [C zz ] is close to unity in both the Floquet and energy conserving cases, and thus retains long-time memory: (Energy conserving) (114) Next, the late-time memory implies that small subsystems of the chain have sub-thermal entanglement entropy. This prediction is in agreement with numerical observations in Ref. [58].
Finally, dynamics in the finite size crossover is characterised by a dynamical exponent z = ∞ as per the logarithmically growing light cone (see (50) and Fig. 8). The length-energy relationship set by the matrix elements t ∼ v(r) −1 determines the light cone; any l-bits outside the light cone are not entangled with the probe spin. In the thermal phase, we expect that the logarithmic expansion of the light cone crosses over to ballistic or diffusive expansion for t > ω −1 ξ in Floquet and Hamiltonian systems respectively. Ref. [62] numerically observed stretched exponential decay of typical spatial correlations in eigenstates in the MBL-thermal crossover region and noted the similarity of their numerical results to that near an infiniterandomness fixed point. Although we do not flesh out the connection between the RM transition and the infiniterandomness transition here, we note that both theories predict z = ∞ and logarithmically growing light cones.

D. Scale-free resonances near the finite-size crossover
In region I (and II), the probability of resonance at range r is scale free resulting in the formation of resonances on all length scales. This feature of the thermal-MBL crossover in small systems has been observed numerically in Ref. [62].

E. Linear drift of critical disorder strength with L
The RM predicts a ubiquitous feature of small system numerics on disordered chains: that the critical disorder strength increases approximately linearly with L.
Refs. [1] and [2] argued this drift to be inconsistent with the existence of MBL; the RM however provides an alternative explanation.
The origin of the drift lies in the localised nature of region I. On increasing 1/ζ at small sizes, the chain crosses over from thermal to localised behaviour when the correlation length first exceeds the system size |ξ| ≈ L (see Fig. 1). The critical 1/ζ (and equivalently the critical disorder strength) thus increase with L.
This drift can be quantified: let W δ (L) denote the disorder strength at which the time-averaged correlator [C zz ] deviates from its value in the infinite temperature Gibbs ensemble by some small amount δ, For the Hamiltonian RM, algebraic manipulation of (91) with 1/ξ = log(W/W c ) yields: for some δ-dependent constant δ . Over a regime of sufficiently small L, this function is approximately linearly increasing with L (see Appendix D for derivation).
More generally the linear growth of W δ follows from Taylor expanding ξ near W = W δ . Precisely, if we identify ξ(W δ (L)) ∝ L, (for some δ-dependent constant of proportionality), and consider the taylor expansion .
Eq. (119) and the linear-in-L drift of the critical point follow provided W is sufficiently far from the transition that i) the Taylor expansion is valid (i.e. |W − W δ (L)| < |W δ (L) − W c |) and ii) that ξ (W δ (L)) is slowly varying in L. Fig. 9 plots the probe spin entanglement entropy [S P ] averaged over all eigenstates, in the Hamiltonian RM vs the re-scaled disorder strength (using 1/ξ = log(W/W c )). The probe spin entropy is maximal in the cat states, and is zero is the fraction [C zz ] = F (L/2) of states that do not resonate. The inset confirms that the deviation (W δ − W c ) increases linearly with L at small L, before converging to zero from below at large L. A similar analysis in the Floquet RM predicts a linear drift of the critical frequency at which localisation sets in with L for fixed disorder strength. Refs. [1] and [2] numerically studied the scaling of the Thouless time with disorder strength in the thermalising phase. The Thouless time is defined as the time-scale above which random matrices govern quantum dynamics in chaotic systems, or equivalently as the inverse of the energy scale below which the random matrices govern eigenstate properties. Through a detailed study of the spectral form factor and [S(ω)], Refs. [1] and [2] argued that the inverse of the Thouless time ω Th. exponentially decreases with disorder strength: Should this behaviour continue asymptotically as L → ∞, then the numerically observed MBL-thermal crossover is simply a finite-size effect caused by ω Th. becoming smaller than the Heisenberg time ω −1 H . That is, the observed localisation is simply a consequence of the small sizes accessible to exact numerics.
The RM provides an alternate explanation for (121) while allowing for a MBL phase. In a diffusive system, the Thouless time is set by the time taken by a localised packet of energy to spread over the system. For diffusion constant D, thus ω th. = D/L 2 . As the packet takes time ω −1 ξ to spread a distance ξ, D = ω ξ ξ 2 . Combining these estimates where ≈ indicates the dropping of an O(1) factor. Next, consider the correlation length ξ(W ). It is a smooth function of the disorder strength W and diverges at the critical disorder W c defined by ζ = ζ c . As discussed in Sec. V E, the crossover from spectrally averaged statistics being close to their thermal values, to close to their localised values occurs at disorder strength W δ , a much weaker disorder strength than W c in small systems sizes. We may thus Taylor expand ξ near W = W δ (as in (118)) from which the exponential dependence of the Thouless time on the disorder strength W of (121) follows.

G. Apparent sub-diffusion in the RM thermal phase
Eqs. (53) and (91) predict a continuously varying exponent for the spectral function [S(ω)] ∼ ω −1+θ above a threshold frequency scale ω ξ in the thermal phase. The RM thus explains the apparent sub-diffusion (as measured by the dynamic exponent 1/θ) reported in several studies [67][68][69][70][71] without any reference to rare regions, and indeed predicts such apparent sub-diffusive behaviour even in Floquet systems without any conservation laws. This prediction of the RM may resolve a mystery about the absence of broad distributions of the conductivity (across disorder realisations) that are expected in a sub-diffusive regime characterised by weak links [55,56].
We note that Ref. [76] (in the supplementary material) previously speculated that rare resonances may lead to apparent sub-diffusive behaviour in the thermal phase.
H. Exponentially enhanced sensitivity to eigenstates or 'maximal chaos' The fidelity susceptibility χ a measures the sensitivity of an eigenstate |E a to perturbation by a local operator U . It is defined as (123) The mean of the logarithm of χ (defined as the average of log χ a across infinite temperature eigenstates and disorder realisations) shows the following scaling with L: Ref. [2] made two observations about the distribution of log χ a at numerically accessible sizes. First, there is a regime of maximal chaos separating the thermalising and MBL regimes in which [log χ] ∼ L · 2 log 2, ("maximal chaos".) Second, the tails of the distribution in the putative MBL regime (in which [log χ] saturates) are fatter than expected from a Poisson distribution. The authors explained both observations through the exponential enhancement of matrix elements between eigenstates with energy differences comparable to the many-body level spacing, and concluded that such enhancement is inconsistent with MBL. The RM explains both observations in Ref. [2] assuming a thermodynamic MBL phase.
Consider a pair of resonant cat states |E a,b = (|E a ± |E b )/ √ 2 involving the re-arrangement of l-bits at range r = L/2 and splitting comparable to or less than the many-body level spacing. A generic local perturbation U will couple these states as E a |U |E b = O(|U |) [77]. Consequently, their fidelity susceptibility is very large, increasing as ∼ 2 2L .
In the numerically accessible MBL-thermal crossover, a finite fraction q(L/2)∆L of the eigenstates are involved in resonances with range between L and L + ∆L and splitting comparable to the many-body level spacing. The RM thus predicts maximum chaos (125) at the finite-size crossover. More precisely, in regions I and II of the Floquet RM where ρ(r) sets the typical inverse level spacing for a resonance at range r, and s(r) = q(r) exp(− L/2 r q(r )dr ), is the probability that the longest range resonance for a given state is at range r. Thu, maximum chaos is approached as L becomes closer to λ.
In the RM MBL phase, the fraction of states involved in system-wide resonances q(L/2) is exponentially small in L. These states thus do not contribute to [log χ], which is independent of L. Nevertheless, these rare states lead to increased weight in the tail of the distribution of log χ. This explains the second observation of Ref. [2]. Note that (130) makes no reference to a l-bit basis. colourWhen (v) ∝ v −2+θ0 at small v, the stability criterion implies that 0 ≤ θ 0 < 1 for MBL. Eq. (130) generalises the stability criterion to thermalising avalanches introduced in Ref. [48]. Ref. [48] studied the stability of a MBL system composed of l-bits to a thermalising inclusion, and argued that ζ (the length scale controlling the localisation of a physical spin operator in the l-bit basis) must be smaller than ζ c = 1/ log 2. Re-writing the avalanche criterion in terms of properties of off-diagonal matrix elements, we obtain (130) with no reference to either rare regions or to l-bits.

VI. DISCUSSION
We have presented the RM, a model of the finite-size MBL-thermal crossover in which the localised phase is destabilised by many-body resonances, rather than rare low-disorder regions. The RM is consistent with a stable MBL phase, and reproduces several numerically observed features of the MBL-thermal finite-size crossover, including the controversial observations of Refs. [1,2]. Fig. 10 re-plots the [S(ω)] data in Fig. 2 of Ref. [2]. The plot shows the frequency dependence of [S(ω)] at several disorder strengths 0.5 ≤ W ≤ 2.5 in the putative thermalising phase of the disordered spin-1 2 XXZ chain. Ref. [2] argued that the data is consistent with the scaling law [S(ω)] ∼ C/(W ω) (black horizontal line) over an increasing range of frequencies. We instead argue that the data is consistent with the scaling law predicted by the Hamiltonian RM with a logarithmic correction (dashed black line). Indeed, the curves for W 1 align with the RM prediction over ≈ 1.4 decades in frequency, while evidence of the plateau predicted by Ref. [2] is visible only in two of the curves with W ≈ 1.5, 1.75, and over less than a decade in frequency. The behaviour of the curves with W ≈ 1.5, 1.75 is however noteworthy, and not immediately explained by the RM. To settle the debate between the two scaling predictions requires more systematic numerical investigation of the effects of system size on the curves in Fig. 10. Specifically, numerics at larger L should reveal which of the two regimes (the linear growth or the plateau) expands with increasing L.
The RM makes several numerically testable predictions about Floquet and quasi-periodically modulated spin chains. First, Sec. V applies without alteration to the quasi-periodic case. Second, the exponent θ c controlling the strongest low-frequency divergence of the spectral function in region I in the Floquet case is nonuniversal and non-zero, in contrast to the Hamiltonian RM with θ c → 0 + . Third, Floquet systems on the thermalising side of the finite-size crossover would also exhibit apparent sub-diffusive scaling in their spectral functions. The origin of this apparent sub-diffusion is the formation of many-body resonances on length scales shorter than ξ. Fourth, irrespective of the type of disorder or the number of conservation laws, we predict logarithmically growing light cones in the thermalising phase for t ω −1 ξ . Finally, observables conditioned on the formation of resonances could detect the MBL-region I crossover in Fig. 1a.
Eq. (131) offers a new numerical criterion to differentiate localised and thermalising systems. Analogous to the G parameter in Ref. [78] and the typical fidelity susceptibility [2], ρv is exponentially larger in L in the thermalising phase as compared to the MBL phase. Preliminary work on a disordered Ising model suggests that (131) bounds the transition out of the localised phase to larger disorder strengths than other standard criteria based on energy level statistics or eigenstate entanglement entropies.
Future work could explore the RM along several axes. The first is to establish whether the distribution of sample conductivities (across disorder realisations) predicted by the RM is consistent with the observations of Ref. [55]. This would add further evidence to the claim that manybody resonances, and not rare regions, give rise to the apparent sub-diffusion observed numerically.
The second is to compare the eigenstate correlations predicted by the Hamiltonian RM to those from the Anderson model on the random regular graph (RRG) [79,80]. The RRG Anderson transition is believed to model the MBL-thermal transition if one identifies each site of the RRG with a computational basis state of a disordered spin chain [81]. Using Mott-type resonance arguments similar to those of Sec. III, Ref. [79] recently argued that in the RRG localized phase, the correlator [tr (Π n (t)Π n (0))] (where Π n (t) is the time evolved single site projector onto the site n) has a Fourier spectrum β(ω) which diverges as a power law as ω → 0. Identifying each Π n with |E a σ E a σ|, a product state of the probe spin and the disordered chain, the RM predicts that β(ω) diverges exactly as [S(ω)] (27). The reconcilation of the RM with the RRG is however less apparent in the thermal phase, where the latter predicts a correlation length that diverges with a different exponent than in the RM.
(A4) Accounting for the distribution of g ba in (25) corrects λ, the weight at zero frequency and the exact form of [S(ω)]. However, it does change universal features, such as the vanishing of the exponent θ with |ζ − ζ c | and the exponential decay in r of F (r), the weight at zero frequency after all range r ≤ r processes have been accounted for, as per (45).

Multiple resonances
Suppose an eigenstate | a , ↑ is resonant with multiple other eigenstates of H 0 . Here we argue that the strongest resonance (defined by (35)) sets the frequency of oscillation of a ↑ |σ z P (t)σ z P (0)| a ↑ . Consider the case of two resonances at different ranges. Let |ε α = 1 √ 2 (| a ↑ + | b ↓ ) denote the cat state resulting from the strongest resonance (at the shorter range). Suppose that |ε α is now resonant with another state | c ↓ at larger range with some matrix element This matrix element is much smaller than |V ba | as |V ac |, |V bc | |V ba |. Treating this resonance within degenerate perturbation theory splits the peak at ω = |V ba | into two peaks at ω = |V ba | ± |V αc |. As this further splitting is small, we neglect it and assume that the spectral weight remains sharply peaked around ω = |V ba |.
In the time domain this statement is as follows: an initial state | a ↑ oscillates between | a ↑ and | b ↓ on a time scale |V ba | −1 and tunnels into the state | c ↓ on the much longer timescale |V αc | −1 .
We generalise the above argument to many-resonance case. Suppose |ε α has a resonance meditated by a matrix element |V αc |, which leads to hybridised states Take one of these states |ε α+ . Suppose this state has a longer-range resonance mediated by a matrix element |V αd |. We obtain two new cat states. Suppose one of these two cat states |ε α+ has an even longer-range resonance mediated by |V αe | and so on. The initial peak at ω a↑ = |V ba | splits into several peaks at An analogous procedure splits each of the peaks with a minus sign in the RHS above into many sub-peaks.
To show that such shift ∆ω remain unimportant we calculate the root-mean-square size shift ∆ω 2 as show that ∆ω 2 ω a↑ . To do this we first note that the matrix elements v(r) connecting an already hybridised state to other unhybridised states at range r are a factor where as the density of states is twice as large yielding a probability of hybridising at range r of Thus, supposing that the initial resonance is at a range r (i.e. that ω a↑ = v(r)) we find where X(r) is a random variable which takes values X(r) = 1, −1, 0 with probabilities q (r)/2, q (r)/2, 1 − q (r) respectively. Thus ∆ω has mean ∆ω = 0 and, measured in units of the initial resonant frequency ω a↑ , has variance On the localised half of the phase diagram (ξ > 0) this quantity is exponentially decaying in r, indicating this approximation scheme is asymptotically improving at low frequencies.
In the crossover region it is bounded by its critical value, which is much smaller than unity and so does not alter the asymptotic form of the spectral function [S(ω)], whereas on the thermal this approximation breaks down only for r > ξ, outside the regime of validity of our calculation. In this regime hydrodynamic constraints become important. These constraints highlight the limitations of the approximation made in (90), as F 1 predicts unphysical behaviour. Specifically which using [C zz (t)] = F (r(t)) (49), and the logarithmically growing light cone r(t) ∝ log t implies that the correlations decay as a stretched exponential in log t. This decay is slower than any power law, but much faster than the maximum possible decay rate permitted by energy conservation of This maximum rate follows as the z-field on the probe spin σ z P has overlap with the Hamiltonian tr σ z P H = W , and any initial energy on the probe spin cannot have spread further than the light cone front r(t).
In order to address this inconsistency we turn to a more careful treatment of Eqs. (86) and (88). By direct numerical integration (see Appendix C 1) we find that the stretched exponential decay is cut-off at r λ by an asymptotic decay F (r) ∼ r −2 , implying a decay [C zz (t)] ∼ log −2 t. This decay is still too fast to be consistent with hydrodynamics, however, the weakness of this violation means there are many small corrections which yield a late time dynamical regime consistent with hydrodynamics. For example, a sub leading power law in r on the matrix elements v(e, r) will suffice. However, here we explore the effect of energy dependency of the matrix elements.
Instead of the energy independent form for the matrix elements (77), we now consider v(e, r) = J exp − r ζ(e/r) − r ζ c .
where we now allow the localisation length to vary as a function of the energy density e/r of the patch of the system which must be rearranged to relate the two states |E a ↑ and |E b ↓ (As we are interested only in behaviour at asymptotically large r, we consider these states to be at the same energy density, despite their energy difference of ±W due to the probe spin). We consider only the leading order dependence on energy density of the localisation length 1 ζ(e/r) where ζ is the localisation length at maximum entropy, the constant energy densities µ, η determine scales over which ζ varies, and we have suppressed higher powers of e/r. We will assume η = ∞ as the statistical symmetry of the model impliesζ should be an even function, and µ positive and finite. This corresponds to a localisation length which is shorter away from maximum entropy. The energy dependence of the matrix elements then alters the form of q σ (e, r): For µ positive and finite q σ (e, r) is asymptotically narrower than ρ σ (e, r) at large r, we can extract the asymptotic behaviour of f σ by replacing q σ (e, r) with a delta function where γ = deq σ (e, r) is an r independent constant at the critical point. Solving (B6) (see Appendix (C 2)) we find asymptotic decay where here ∼ indicates asymptotic equality up to an overall constant. This yields consistent with hydrodynamic restrictions.
Appendix C: Solutions to the loss-diffusion (86) In this appendix we consider the loss-diffusion equation (86) We study two regimes:  Figure 11. Decay in F (r) for energy independent matrix elements: Values of λ|dF/dr| are plotted versus r/λ, these are obtained by numerically solving (C2) and (C4). The point r/λ = 1 is marked with a vertical grey line. For r/λ < 1, the behaviour is consistent with F (r) = exp(− r/λ ) (dotted line). For r/λ > 1, the decay is slower F (r) ∝ (λ/r) 2 (dashed). Different series correspond to different values of λ (legend inset).
• We first study the critical dynamics (ζ = ζ c ) with energy independent matrix elements (v a function of r only). We show that the asymptotic decay of F (r) = def σ (e, r) is given by F (r) ∝ r −2 as quoted in the main text. This behaviour is not permitted asymptotically due to hydrodynamic restrictions.
• We then study the asymptotic critical dynamics for energy dependent matrix elements (B3) with η = ∞, and 0 < µ < ∞. We show that in this case F (r) ∼ r −1/2 , behaviour consistent with hydrodynamics.

Critical point with energy independent matrix elements
Here we study the equation defined in the main text, specifically for the loss function q ↑ (e, r) ∼ 1 √ 4λr exp − (e + 1 2 W ) 2 2s 2 e (r) .
where s e (r) = W √ 2r + 1 . We numerically solve these equations by stochastic sampling of trajectories. In Fig 11 we plot dF/dr for different values of the parameter λ where as before F (r) = def ↑ (e, r).
We see that for all trajectories the initial decay at small r λ is consistent with the approximate solution F (r) = exp(− r/λ ) (grey vertical line marks r = λ) at which there is a crossover to F (r) ∝ r −2 behaviour. For these equations this latter behaviour continues asymptotically.
In Fig. 12 we show the variation of f ↑ (e, r) with e, specifically we plot f ↑ (e, r) for a series of fixed log-spaced values of r. For clarity we also re-scale e by the width of the distribution s e (r) = W √ 2r + 1 (i.e. so that for λ = ∞ the plots would collapse for all r). From these plot it is clear that the centre of the distribution is depleted faster than the mean, that is f ↑ (0, r) decays asymptotically faster than F (r). This behaviour is exhibited for r λ and violates the approximation scheme of Sec. III B 4.

Critical point with energy dependent matrix elements
We now study the same loss-diffusion equation (C1) for dynamics in the crossover region with energy dependent matrix elements. Specifically we now set for some finite µ in the range 0 < µ < ∞.
To simplify the problem we make several approximations which do not alter the asymptotic behaviour of these equations. First, as the width of q σ is asymptotically smaller (in r) than s e (r), for r λ we can approximate q ↑ (e, r) with a delta function placed at the origin with weight γ = deq ↑ (e, r) = W µ λ ζcπ (4W 2 + ζ c µ 2 ) + O(r −1 ).
(C6) Second, we neglect the sub-leading r-dependent correction to γ, and thirdly we neglect the initial energy offset of f ↑ . This yields the equation with boundary condition f ↑ (e, r = 0) = δ(e).
With this f 0 is straightforwardly identified and it further follows that for n > 0 f n (e, r) = −γ and recognising the resulting summation as a Taylor series, this yields F (r) = e r/ Erfc r/ where is the usual complementary error function. From (C17) it follows that F (r) decays asymptotically as as quoted in the main text. The constant pre-factor here is liable to be altered by the simplifications we made earlier in the calculation, however the asymptotic behaviour F (r) ∝ r −1/2 is robust.
Appendix D: Linear drift of the deviation from thermal behaviour In this appendix we derive (117) from the main text where δ is some δ dependent constant, and W δ (L) is defined as the disorder strength at which the time averaged correlator [C zz ] deviates from thermal behaviour by some small amount δ Recalling that [C zz ] = F (L/2) and using the form (91) for F (r) on the thermal side where we have explicitly labelled disorder dependence of the correlation length ξ and the resonance length λ. We use whereas λ is given by (83). Let us extract from (D3) how W δ varies with L. Away from the crossover region the imaginary error function can be written in terms of more familiar functions Substituting both (D5) and λ(W δ ) = (W δ /W c ) 2 λ(W c ) into (D3) and rearranging we obtain Consider the RHS of (D6): for sufficiently small δ we are far from the crossover L |ξ| and the corrections may be neglected. Now consider the leading term on the RHS of (D6): this term exhibits weak logarithmic dependence of L, and, recalling that ξ(W δ ) ≈ 1/ log(W δ /W c ), doubly logarithmic dependence on W δ , thus to first approximation the RHS may be replaced by a (negative) constant − δ : Then, again using ξ(W δ ) ≈ 1/ log(W δ /W c ), by rearranging we obtain the desired result (D1). This function is approximately linear for sufficiently small L. To see this, note that the RHS of (117) has an inflection point at L = δ /2 − 1, and thus has zero curvature at this point. Taylor expanding about the inflection point and demanding that the cubic term is not larger than the linear term reveals the approximate linearity to persist for L + 1 δ (1/2 + 3/4 ).