Universal properties of many-body delocalization transitions

We study the dynamical melting of"hot"one-dimensional many-body localized systems. As disorder is weakened below a critical value these non-thermal quantum glasses melt via a continuous dynamical phase transition into classical thermal liquids. By accounting for collective resonant tunneling processes, we derive and numerically solve an effective model for such quantum-to-classical transitions and compute their universal critical properties. Notably, the classical thermal liquid exhibits a broad regime of anomalously slow sub-diffusive equilibration dynamics and energy transport. The subdiffusive regime is characterized by a continuously evolving dynamical critical exponent that diverges with a universal power at the transition. Our approach elucidates the universal long-distance, low-energy scaling structure of many-body delocalization transitions in one dimension, in a way that is transparently connected to the underlying microscopic physics.


I. INTRODUCTION
The laws of thermodynamics can break down in disordered quantum systems that are isolated from external heat sources, due to the localization of excitations that would ordinarily transport energy between distant regions to establish thermal equilibrium [1][2][3][4]. Such many-body localized (MBL) systems have the remarkable property that all high-energy excited states behave exactly like zero-temperature quantum ground states. This raises the intriguing possibility that quantum coherent phenomena, typically associated with zero-temperature systems, can occur in arbitrarily "hot" matter. Examples of such quantum coherent phenomena include symmetry breaking below the equilibrium lower critical dimension [5], topological edge states [5,6], and quantum criticality [7][8][9]. MBL systems also host novel out-ofequilibrium dynamical phase transitions, where thermodynamics breaks down sharply at a critical point [10][11][12]. Cold atomic and trapped ion systems offer a promising experimental platform to explore these theoretical ideas. Indeed, issues of thermalization and excited state dynamics necessarily arise in such systems, as they are inherently well isolated from their surroundings, and typically cannot be cooled to low temperatures (compared to their characteristic energy scales).
Theoretical investigations into these questions, however, must tackle a daunting combination of out-ofequilibrium quantum dynamics, interactions, and disorder. Consequently, most existing theoretical work on MBL systems has been confined to small-scale numerics and analysis of phenomenological models. Remark-ably, because of the short localization length and shortrange entanglement structure of MBL eigenstates deep within the localized phase, these approaches have met with considerable success in gleaning properties of the MBL phase itself. In contrast, these same techniques are poorly suited to access the universal properties of disordered criticality and dynamical phase transitions out of the localized phase. Such transitions are driven by long distance properties whose characteristic length scale diverges as the critical point is approached, making numerical studies increasingly unreliable in its vicinity. New analytic approaches are clearly needed to understand the nature of this new class of dynamical transitions. To this end, we expand upon a real-space renormalization group (RSRG) procedure introduced in Ref. 7 and 8, that extends previously developed RSRG techniques for zero-temperature ground states [13][14][15][16][17][18][19][20][21], in order to treat properties of highly excited states. Whereas prior RSRG studies of excited states required (often costly) numerical sampling procedures [8], a dramatic simplification occurs for very high energy states near the center of the many-body spectrum, that enables asymptotically exact, closed-form analytic solutions [9].
Using this method, we compute the properties of highly excited states of a broad family of random quantum spin chains in the limit of strong disorder. We find a sequence of non-ergodic phases in which every excited eigenstate behaves like a zero-temperature quantum critical system -including algebraic decay of correlation functions and critical scaling of entanglement entropy S(L) ∼ log L for a subregion of size L. We dub such systems quantum critical glasses (QCG). QCGs are characterized by an absence of energy transport and arXiv:1501.03501v1 [cond-mat.dis-nn] 14 Jan 2015

QCG Transition
FIG. 1. Generic phase diagram and finite-size crossovers for one-dimensional delocalization transitions for (a) quantum critical glasses (QCG) (ψ < 1) and (b) non-critical many-body localized (MBL) glasses (ψ = 1). W parameterizes disorder strength. For strong disorder W > Wc, a non-ergodic glass is obtained, which lacks thermal transport but exhibits slow dephasing and entanglement dynamics with polynomial scaling of length (L) with log-time τdeph (where the subscript "deph" emphasizes that this is the dephasing time associated with virtual quantum fluctuations rather than the energy transport time scale, which is strictly infinite in the strong disorder glass). Below a critical disorder strength, Wc, the quantum glass melts into a thermal, energy-conducting liquid. (a) Melted QCGs for W Wc exhibit a sequence of three crossover length scales, ξ ∼ |W −Wc| −ν , L * ∼ e ξ ψ , and L diff ∼ exp ξ ψ 1−ψ that delineate three regimes with distinct scaling relations between transport time τ and length L (as indicated in the figure). (b) Many-body localized phases (ψ = 1) have a similar phase diagram except that L diff = ∞, and there is no change in the transport properties across L ∼ L * -hence for MBL systems the thermal liquid in this case shows subdiffusion at all scales with power law scaling between time and length, τ ∼ L z with dynamical exponent z ∼ ξ that diverges continuously upon approaching the transition.
slow dynamics of dephasing and entanglement, that display a power-law relation between length and the logarithm of time: log τ (L) ∼ L ψ , with universal exponent ψ < 1 [22]. Concrete examples include critical random Ising and Potts models. The universal critical properties of these models, as well as disordered versions of an infinite family of one dimensional quantum critical points described by the minimal models and parafermionic models of 1+1d conformal field theories (CFTs) are neatly unified by a description in terms of SU(2) k anyon chains. We then analyze the stability of these quantum critical glass phases by modifying the RSRG procedure to account for the formation of rare resonant tunneling links, which become increasingly important for weak disorder. From these results, we derive an effective model for the melting of these quantum glasses via the formation of a percolating network of resonantly coupled links. Numerical simulation of this model reveal a continuous dynamical phase transition characterized by a diverging length scale ξ, in which quantum glasses melt into classical thermal liquids. The predicted phase diagram and crossover structure of this effective model are summarized in Fig. 1a. Crucially, we find that the critical resonant network is very sparse, containing a vanishing density of spins (corresponding to fractal dimension less than one) in the limit of infinite system size. We argue that the sparsity of the critical cluster enables a controlled calculation of the critical properties within our effective model.
We also study non-critical MBL glasses (Fig. 1b) that, like QCGs, are characterized by an absence of thermal transport and logarithmically slow dephasing and entanglement dynamics (corresponding to effective ψ = 1). Unlike QCGs, the spectrum of MBL systems can be labelled by an extensive set of exponentially-well localized conserved quantities [23,24], such that eigenstates exhibit boundary law scaling [25] of entanglement entropy S(L) ∼ constant, characteristic of gapped quantum ground states. We argue that a similar effective model of percolating resonance networks also accurately describes the disorder-tuned melting of such non-critical MBL glasses. This MBL-to-classical transition was recently studied by exact diagonalization [10,11] and a phenomenological renormalization group procedure [12]. The effective model also provides insight into the dynamics of thermal transport and entanglement of the delocalized thermal liquid, near the melting transition (see Fig. 1). For quantum critical glasses (ψ < 1), our approach predicts that the delocalized phase exhibits ordinary diffusive thermal transport at the very longest time scales. However, we also find an extremely broad intermediate regime of subdiffusive transport characterized by a variable dynamical critical exponent that diverges continuously upon approaching the transition. For ψ < 1, the subdiffusive behavior persists up to a crossover length scale L diff that diverges faster than exponentially upon approaching the transition, with ordinary diffusion recovered on longer length scales. In contrast, we find that L diff becomes infinite (even away from the critical point) for melted MBL glasses, which consequently exhibit subdiffusive dynamics on all length scales in agreement with a recent phenomenological renormalization group approach [12].

II. MODELS AND STRONG DISORDER GLASS PHASES
A. Critical Spin Chains The starting point of our analysis is a family of strongly disordered quantum spin-chains, exemplified by the random quantum Ising and Potts models: where σ µ i in (1) are 2 × 2 Pauli matrices [26], and iσi . Here, the bond-and transverse field-strengths J i > 0 and h i > 0 are independent random variables whose logarithms have averages log J and log h respectively and standard deviation For log J = log h, the zero-temperature ground state of these models exhibits a quantum phase transition between a symmetry-preserving localized paramagnetic phase (log h > log J) and a symmetry-breaking localized magnetic phase (log J < log h). In this paper, we are interested in whether localization and quantum criticality persist in highly excited states.

B. Anyonic Description
The universal critical properties of the critical random Ising (1) and Potts (2) chains, as well as disordered analogs of an infinite family of one-dimensional quantum critical points described by the minimal models and parafermionic models of 1+1d CFT are neatly unified by a description in terms of chains of SU(2) k anyonic spins, S i . Though not fundamentally necessary, this anyonic description provides a convenient language for analyzing excited state dynamics of strongly random critical points.
The detailed properties of anyons are detailed in [19-21, 27, 28], and summarized in Appendix A. For most purposes, SU(2) k anyons can simply be thought of as Heisenberg spins, except that their maximum spin size is bounded by S ≤ k 2 . Specifically, the anyonic spins take values S ∈ {0, 1 2 , 1, . . . , k 2 } and obey the modified spin addition (fusion) relations S 1 ⊗ S 2 = |S 1 − S 2 | ⊕ · · · ⊕ min (S 1 + S 2 , k − (S 1 + S 2 )), rather than the usual angular momentum addition rules familiar from elementary quantum mechanics. A corollary of this modified algebraic structure is that the underlying Hilbert space is nonlocal. Hence, unlike spins, anyonic particles do not possess an 'on-site' degeneracy in a strict sense, though we may define a local quantum dimension d S = sin π(2S+1) k+2 / sin π k+2 , that describes the multiplicative growth of the Hilbert space dimension upon adding a single anyon of spin S.
The dynamics of these chains are governed by the Heisenberg-like Hamiltonian: The truncated dot product S i · S i+1 k generalizes the usual Heisenberg Hamiltonian to truncated SU(2) k spins [9], assigning monotonically increasing energies to larger spin sizes resulting from fusing S i and S i+1 (see Appendix A). This energetically favors S = 0 singlet bonds between neighboring spins, just as for conventional Heisenberg chains. Absent disorder, chains of S = 1/2 anyons described by (3) are gapless and described by the k th minimal model (resp. Z k parafermionic model) of conformal field theory for antiferromagnetic (resp. ferromagnetic) couplings [27,28]. Here, we will be interested in the high energy density properties of strongly disordered chains. In this limit, the physics of such chains can be understood in terms of the fixed-point behavior of a strong disorder renormalization group method [9] (see Appendix A for details), and we will refer to these disordered phases as quantum critical glasses (QCGs). These dynamical fixed points also describe the universal properties of the critical points of the random Ising (1) and Potts (2) chains for k = 2 and k = 4, respectively [29].

C. Non-critical MBL glasses
We also consider non-critical MBL glasses, exemplified by interacting spinless fermions subject to a random chemical potential, where J ⊥ is the hopping amplitude, V J ⊥ is the interaction strength, and µ i is the chemical potential whose log-standard deviation, W , serves as a measure of disorder. This fermionic problem maps (via a standard Jordan-Wigner transformation) onto an XXZ spin chain with random fields µ i , with anisotropy parametrized by the interaction V . This Hamiltonian has been considered as a prototypical example of a many-body localized (MBL) system (see e.g. [30] for a recent review), and when V = 0, the system is believed to exhibit a MBLto-thermal dynamical phase transition as a function of disorder strength [10][11][12].

D. Strong disorder glass phases
For strong disorder (W 1), all of these models possess a non-ergodic glass phase characterized by the absence of d.c. energy transport and anomalously slow dynamics, reflected in a power-law relation between length and the logarithm of time: | log t| ∼ L ψ . Physically, ψ ≤ 1, since processes with ψ > 1 are shorted out by quantum tunneling that decays at most exponentially with distance. The specific case of a quantum glass with ψ = 1 corresponds to the non-critical MBL phase (e.g., the strong-disorder phase of (4)) that is characterized by exponential decay of disorder-averaged correlation functions and short-range boundary-law scaling of entanglement, S(L) ∼ constant for a subregion of size L [25]. These results are consistent with the existence of an extensive set of exponentially well-localized integrals of motion [7,23,25]. In contrast, phases with ψ < 1 -exemplified by the family of critical SU(2) k anyon chains, with ψ = 1 2 for k = 2 (Ising), ψ = 1 6 for k = 4 (Potts), and ψ ∼ 4π 2 k 3 for k 1 (see Appendix A for derivation and general formula) -have quantum critical structure in every eigenstate, with algebraic decay of disorderaveraged correlations, and critical entanglement scaling, S(L) ∼ log L.
Though there is currently no known analytic approach for solving non-critical MBL models, the properties of these quantum critical glasses can be exactly computed [9] in the strong disorder limit using an analytically tractable real space renormalization group (RSRG) approach. We refer the reader to Appendix A for details, and simply quote the results where needed. In the balance of this paper, we will construct and analyze an effective model of the resonant tunneling processes that determine the fate of both MBL and QCG phases as the strength of disorder is weakened.

III. DELOCALIZATION VIA A FRACTAL RESONANT NETWORK
To motivate the problem of delocalization, imagine creating an energy wavepacket initially localized near a single site (bond) in the MBL (QCG) [31]. How does such a packet evolve with time? Deep in the localized phase, it is extremely unlikely that it can tunnel from a given site (bond) to another, due to the wide variance of site (bond) energies. More precisely, the amplitude to tunnel through a distance x of the localized phase is Here, J 0 ≈ V is the interaction strength for MBL systems (ψ = 1) while for QCG systems we compute J 0 = e −W using the RSRG method (in units where the maximal bond strength is set to 1, see Appendix A). If the tunneling strength J(x ij ) between two sites (bonds) i and j separated by distance x ij is much smaller than their energy difference δ ij , then the true many-body eigenstate is very close to a product state of independent configurations of each bond. On the other hand, if J(x ij ) δ ij , then the manybody eigenstate contains an entangled superposition of the degrees of freedom on the two sites (bonds), and we say that the two sites (bonds) are resonantly linked. We estimate the probability that a given site or bond is resonantly linked to at least one other is ( is the Gamma function and ρ(0) is the density of sites/bonds per unit energy and length: ρ(0) MBL ≈ 1 J0 and ρ(0) QCG ≈ 1 W . For very strong disorder (x 0 1), the density of resonantly linked pairs is very small, ρ res ≈ λ R 1, so that the resonant links are well isolated from each other, and they do not disrupt the properties of the surrounding localized phase. As disorder is weakened, these resonant links occur more frequently, and eventually disrupt local- ization below a critical disorder strength W c . Crucially, however, we will see that the delocalization transition occurs in a regime where resonances are still dilute and rare. We now construct an effective model to analyze this delocalization transition.
A cluster of m resonantly linked sites is characterized by non-overlapping bands of N ≈ 2 m energy levels, each spanning a bandwidth Λ m = m i=1 J(x i,i+1 ), corresponding to typical level spacing δ m ≈ Λ m 2 −m . Roughly speaking, excitations can tunnel between two such clusters containing m 1,2 sites respectively and separated by distance This motivates an effective model for the delocalization transition, in which we examine a disordered anyon chain to identify resonantly linked pairs of sites (bonds). Once we form all such clusters, we account for the reduction of level spacing, and examine whether the any of newly formed clusters interact sufficiently strongly to merge. This process can then be repeated iteratively until all possible resonant clusters are formed, which we do numerically.
Typical results are shown in Fig. 3 for a QCG with ψ = 1/(2 + ϕ), corresponding to SU(2) 3 Fibonacci anyons (ϕ is the golden ratio). The results for this specific model exemplify all of the generic features of the one-dimensional dynamical melting transitions. Similar results for other ψ are shown in Appendix C. We plot the disorder-averaged length of the longest cluster, denoted ξ loc -a convenient measure of the localization length -against the disorder strength. For strong disorder (W > W c ) the growth of resonant clusters ceases after a finite number of iterations, leaving small, well-separated resonant clusters embedded in a localized phase. For weaker disorder (W ≤ W c ), a connected resonant cluster traverses the entire system, enabling thermal transport over arbitrary distances, resulting in thermalization and breakdown of localization. The data for different system sizes collapses when scaled by |W − W c |L 1/ν , indicating the existence of a length scale ξ that diverges when the transition is approached from either side, ξ ∼ |W − W c | −ν , with correlation length exponent ν = 3.3 ± 0.02. This indicates a sharp continuous delocalization phase transition, with the correlation length ξ identified with the localization length [32] and the size of the typical insulating region for W > W c and W < W c respectively. Notably, related finite size scaling data for different models with various values of ψ collapse onto the same universal curve (see Appendix C) -indicating a super-universality of the correlation length exponent, ν.

A. Fractality of the critical cluster
Before turning to a more detailed numerical analysis of this transition, a few comments are in order. First, owing to the exponential relation between level spacing and cluster size, very few bonds are required to form a resonance that spans a distance L. To see this, let us estimate the minimal number of sites (bonds), N (L) required to span a system of length L. In order to form a connected resonance chain, the level spacing δ N ∼ e −N must be comparable to or larger than the tunneling between adjacent links, which are typically a distance L/N apart. Setting δ N ≈ J(L/N ) yields N (L) ≈ L ψ 1+ψ . This corresponds to a density ρ(L) = N/L ≈ L − 1 1+ψ , which vanishes in the limit of large system size.
At criticality, a single resonant cluster just barely spans the system, and contains only the minimal number of bonds, N (L). Hence, we see that the critical resonant cluster actually has fractal dimension ψ 1+ψ < 1. Related fractal structure also arises in equilibrium delocalization transitions, including at the mobility edge in a non-interacting 3D Anderson insulator, and at quantum Hall plateau transitions [33], and is indicative of the self-similar scaling structure of the delocalization critical point.
The fractality of the critical resonant cluster has several important consequences. First, in the limit of large system size the distance between adjacent sites (bonds), ∼ ρ(L) −1 in the critical cluster diverges, and the bandwidth, Λ crit ≈ N (L)J ρ(L) −1 vanishes. Hence, while there is a finite density ≈ λ R of resonantly linked pairs, the critical cluster skips over nearly all of these and involves only a vanishingly small fraction of atypically wellcoupled pairs. Therefore, near criticality, we are well justified in using the tunneling amplitude J(x) (valid at Resonant backbone (a) For an infinite system on the delocalized side of the melting transition (W < Wc), every spin (black dots) is resonantly connected (dashed link) to every other spin, forming a dense cluster with infinite bandwidth and zero level spacing. (b) Performing the resonance merging algorithm only on a finite size subregion, A, with length L L * ∼ e ξ ψ identifies a dilute subset of spins -the "resonant backbone" -that are resonantly connected within a narrow bandwidth, self-thermalize faster than the transport time to cross A, and hence dominate transport. asymptotically long distances), since the critical cluster forms form divergingly large segments of localized phase. For this reason, we expect that our model accurately captures the universal properties (e.g., critical exponents) of the transition, though likely not non-universal properties (e.g., the precise value of W c ).

IV. NEAR-CRITICAL DYNAMICS IN THE DELOCALIZED PHASE
We have seen that the delocalization transition of the strong disorder quantum glass phase is driven by the formation of an infinitely sparse cluster of resonances consisting of a fractal subset of spins. The critical resonant cluster mediates coherent transport of energy with characteristic scaling of length and time log τ ∼ L ψ . This scaling relation is reminiscent of the one governing the dynamics of dephasing and entanglement growth in the localized glass phase -with the important distinction that dynamics in the strong disorder glass phase describe only virtual ("off-shell") transitions rather than real ("on-shell", or resonant) processes required for thermal transport. Hence, whereas the glass is non-ergodic, the energy eigenstates at the glass melting critical point are ergodic and thermal, with volume law scaling of entanglement in every eigenstate (consistent with general entanglement monotonicity requirements [34]). However, while the pure eigenstates at the delocalization transition exhibit thermal behavior, starting from a superposition of energy eigenstates (the only initial conditions that can be prepared experimentally in finite time), the system will take super-polynomially long time (in system size) to equilibrate and thermalize.
Detuning slightly away from criticality into the delocalized phase, W W c , the critical cluster is bolstered by additional resonant links, and develops a finite density, ρ(W ), of spins (independent of system size). In a truly infinite system with W < W c , this finite density of resonantly coupled bonds corresponds to infinite bandwidth and zero level spacing, and hence all degrees of freedom in the system will inevitably be subsumed into the resonant cluster.
To analyze the delocalized and near-critical dynamics within our effective model, consider a sub-region A, of length L, within an infinite chain with W W c . Applied to the full system (Fig. 4a), the iterative resonance merging procedure used to analyze the transition from the localized side will eventually subsume all degrees of freedom in the system. On the other hand, implementing the iterative procedure only on degrees of freedom within A (Fig. 4b) identifies a connected cluster that spans A (with probability exponentially close to one for L ξ), but which contains only a subset of the degrees of freedom in A. We refer to this subset as the resonant backbone of A. Note that there is an exponentially small probability ∼ e −αL/ξ that the backbone cluster does not span A, but as we will argue below, these exponentially rare events can be largely ignored (except in a window of length scales L * L L diff , see below) for quantum critical glasses (ψ < 1).
It is natural to expect that this resonant backbone dominates the dynamics of thermal transport and entanglement spread across A, since peripheral spins not included in this backbone must interact with spins far outside of A in order to thermalize and participate in the resonant cluster. Specifically, we assume that dynamics across A occur more slowly than the time for the resonant backbone of A to self-thermalize, but much faster than the time, τ therm to thermalize peripheral degrees of freedom outside the backbone. We will later critically assess the validity of this approximation, and show that it is self-consistent over a broad range of time scales e ξ ψ τ τ therm ≈ e ξ 2ψ [35]. In this regime, transport occurs via the thermally incoherent transfer of energy among bonds in the resonant backbone. Hence, we model the spread of excess energy (initially localized near position x as a function of time) as a classical random walk across the resonances, with a timescale 1/J typ (L AB ) to transfer excitations between resonantly coupled bonds A and B separated by distance L AB . The iterative merging procedure yields detailed information about the connectivity structure and the coupling strengths for each participating link in the resonant backbone. Using standard Green's function methods, we compute the time evolution of an initially well-localized energy wave packet spreading across the resonances network via a random walk.
Averaging over initial position x, and over disorder configurations, we find that the mean-square displacement of an excitation grows as a power law in time (Fig. 5), For ordinary classical diffusion, z = 2. In contrast, we find that z(W ) diverges in the glassy phase, indicating an absence of energy transport and a breakdown of thermal equilibration. On the delocalized side of the critical point, 1/z(W ) increases continuously from zero as a function of detuning from the critical disorder strength: Though seemingly a distinct universal exponent, ζ is related by a general scaling relation to the log-dynamical exponent ψ of the localized phase, and the correlation length exponent ν. To see this, note that at the critical point, transport along the critical resonant chain occurs by tunneling through long (diverging as a power law in system size) localized regions, hence energy scales with length as ∼ e −L ψ at the critical point, just as in the localized phase. This must cross over smoothly to the power law scaling of (6) on length scales L ≈ ξ, implying the scaling relation: Using this scaling relation and the definition of the correlation length exponent we may rewrite (7) as z ∼ ξ ψ .

A. Structure of resonant backbone
Subdiffusive transport arises from the broad distribution of effective tunneling links, J ij = J typ (L ij ), in the resonant backbone, corresponding to a power-law distribution of timescales [10,36], τ ij ∼ 1/J ij : with 1 < α ≤ 2 such that the mean waiting time ∞ 0 τ p(τ )dτ is divergent. To see the relation between this power-law distribution and subdiffusive transport properties, note that energy transport occurs as a random walk on a one dimensional chain whose "sites" are twospin bonds in the resonant cluster, and whose links are weighted by a waiting time τ corresponding to the inverse effective coupling between "sites" [37]. In N steps, a random walker traverses ∼ O( √ N ) different bonds, revisiting each approximately √ N times, and moving total distance L ≈ √ N (in units of the average step size). For such a broad distribution of waiting times, the time taken for N steps is dominated by the longest waiting time, T , encountered. The probability of encountering a bond with waiting time of at least T is We can reasonably expect to find a link with waiting time τ T among L ≈ √ N different bonds if LP (τ ≥ T ) ≈ 1, and therefore the slowest link in region of size L has waiting time T (L) ≈ L 1 α−1 . Since this weak link is revisited order √ N ≈ L times, the total time to move distance L scales as t energy (L) ≈ L α α−1 . Comparing to (6), we identify the dynamical exponent for energy transport as: implying that α approaches one near the transition as: Note also that the broad distribution of waiting times (9) corresponds to a probability of encountering a link of length l in the transport path along the resonant backbone, near the transition, given by Hence the typical spacing between "vertebrae" in the resonant backbone is ≈ ξ, corresponding to energy scale e −ξ ψ . However, in a region of length L ξ, it is extremely unlikely to avoid encountering an atypically long gap of length , defined by L ξ P ( ) ∼ 1, i.e. ∼ ξ log L ξ 1/ψ . Such long rare links involve waiting time τ ∼ e l ψ ∼ L ξ ψ , which dominates the time required to traverse a segment of size L.
We note that the dynamics of energy transport scale differently than those of entanglement since energy is a conserved quantity that can only be transferred among subregions whereas entanglement can be freely generated [12,38]. Hence, whereas energy transport can be viewed as a random walk of conserved particle-like excitations, entropy spreads deterministically in all directions simultaneously [38]. In order to entangle regions separated by distance L , entanglement must spread across order L links, only visiting each once (in the same spirit as the second law of thermodynamics, once a bond is entangled with many others, it is extremely unlikely to later disentangle itself). Therefore the timescale for entanglement spread is dominated by the longest typical waiting time encountered, t ent ≈ T (L) ≈ L 1 α−1 . Comparing, we find a different effective dynamical exponent for the spread of entanglement: The power-law distribution of time-scales implies an anomalous power-law frequency dependence for energy conductivity, σ E (ω) ∼ ω 3−z −1 (obtained by scaling analysis of the high-temperature Kubo formula, neglecting subleading logarithmic corrections to scaling). This power-law behavior crosses over continuously from σ E (ω) ∼ ω 3 dependence in the strong-disorder glass phase [8,39] and at the critical point towards σ E (ω) ∼ ω 5/2 as disorder is weakened.

B. Consistency and limitations of the effective model for transport
The crucial assumption in our transport calculations was that transport across a region A of length L occurs with time τ (L) longer than the time required for segments of the resonant backbone to self-thermalize, but much shorter than the time required for peripheral spins outside the resonant backbone to thermalize. Having obtained detailed information about the structure of the resonant backbone, we are now in a position to critically assess the range of validity of this approximation.
Let us first estimate the time required for the resonant backbone of a region of length L to thermalize a bond or small resonant cluster, C, by computing the time required for C to resonantly exchange an excitation with the rest of the backbone. Since links in the backbone have typical size ξ and corresponding coupling J(ξ) ∼ e −ξ ψ , N C ∼ δ C e ξ ψ bonds of the backbone must be excited in order to resonantly absorb an excitation of C with energy ∼ δ C . This implies that it takes characteristic time τ therm (C) ≈ e ξ ψ log N C for C to thermalize (dominated by the largest bottleneck encountered in N C links in the resonant cluster).
If C is a sub-cluster of the resonant backbone, then, by definition, its level spacing δ C , is less than the total bandwidth of the backbone: δ C < L ξ e −ξ ψ , i.e. N C < L ξ . Hence, in this case, the thermalization time for C is faster than the transport time across the region L, τ therm (C) < τ (L) ∼ e ξ ψ log L ξ . This justifies our assumption that the resonant backbone self-thermalizes on faster time scales than transport through the backbone.
On the other hand, if C is a peripheral spin, or short segment of resonant cluster outside the resonant backbone, then δ C > L ξ e −ξ ψ (otherwise C would be subsumed by the resonance cluster), and the thermalization time for C exceeds the transport time across L. However, since the maximal level spacing is bounded by δ C,max = J 0 e W for QCGs and δ C,max = max i µ i for MBL glasses, in sufficiently long segments, namely those of length L L * ≈ δ C,max ξe ξ ψ , every spin thermalizes faster than transport through the resonant backbone, invalidating our assumption that transport is dominated by a sparse resonant backbone. In the regime L L * , the backbone has a very large bandwidth (compared to δ C,max ) and contains essentially all bonds. Here, our picture of a sparse backbone with resonating bonds separated by insulating (random-singlet like) regions distributed according to P ( ) ∼ e −( /ξ) ψ (which produces subdiffusive transport τ ∼ L ξ ψ ) breaks down. Instead, for length scales L L * , transport occurs not via a rarified resonant backbone, but via sequentially hopping along adjacent spins, suggesting diffusive transport.
However, even for L L * this thermal liquid is not diffusive as it contains exponentially rare regions of length l ξ where our resonance-merging process would yield a non-percolating backbone. To describe the effects of these rare regions, we adopt a coarse-grained picture in the spirit of [12] where we think of these regions of length as insulating (on the time scale for transport across them -though they will be eventually thermalized on some much longer time scale), distributed according to P ( ) ∼ e − /ξ . In a region of length L L * , the largest such insulating region of length l max that we can expect to encounter satisfies (L/ξ) P (l max ) ∼ 1, so that l max ∼ ξ log L ξ . It will thus take τ max ∼ e ξ ψ (log L) ψ for this rare insulating region to thermalize. If L L diff ≡ e ξ ψ/(1−ψ) , this thermalization time is negligible compared to the time τ ∼ e γ log L it takes for entanglement (resp. energy) to spread ballistically (resp. diffusively) with γ = 1 (resp. γ = 2) through the region A, and such rare events can be ignored compared to diffusion. If on the other hand L * L L diff , these exponentially rare events actually dominate the dynamics, and we expect anomalous transport τ ∼ e ξ ψ (log L) ψ through the region A.

C. Phase diagram
Let us summarize our results in the thermal phase (W < W c ). Near the critical point, the system exhibits anomalously slow subdiffusive dynamics τ ∼ L ξ ψ over a parametrically large range of lengthscales (ξ L L * ∼ e ξ ψ ). In this regime, energy is transported along a sparse backbone with small bandwidth. We expect this regime to be the only relevant one for future numerical studies of the critical point. For L L * , this picture then crosses over to that of a dense backbone with very large bandwidth containing almost all spins, with exponentially rare bottlenecks leading to anomalously slow transport τ ∼ e ξ ψ (log L) ψ for L * ∼ e ξ ψ L L diff ∼ e ξ ψ/(1−ψ) . For MBL systems (ψ = 1), there is no change in the scaling of dynamics between L L * and L L * as both regimes exhibit subdiffusion with dynamical critical exponent z ∼ ξ. Moreover, L diff is strictly infinite, indicating that subdiffusive transport persists z ∼ ξ on all scales, in agreement with [12]. For critical glasses (ψ < 1) however, these exponentially rare bottlenecks are eventually negligible and we expect diffusion at scales L L diff . The resulting phase diagram and crossover scales are shown in Fig. 1.

V. DISCUSSION
To summarize, we have explored the high-energy excited state dynamics for a family of one-dimensional spinchains that include Ising, Potts, and Heisenberg models, as well as an infinite class of anyonic-chains. At strong disorder, exact analytic results for the properties of the highly excited spectrum can be obtained using a realspace renormalization group approach. This identifies a family of non-equilibrium quantum critical glass phases in which every eigenstate has the properties of a zero temperature quantum critical system.
Upon weakening disorder, these glass phases melt into thermal liquids. We have derived an effective model for this melting transition from the scaling properties of the strong disorder phase. This model can be efficiently simulated for very large system sizes (with computational cost scaling linearly in system size), and we argue that it gives an accurate description of its long-distance, lowenergy, universal critical properties. Our effective model can also be applied to study the delocalization transition driven by disorder in many-body localized systems.
This analysis reveals an unusual scaling structure characterized by an extremely broad regime with anomalously slow subdiffusive dynamics, which crosses over to ordinary diffusion on very long length scales for critical glasses, but not for MBL systems. The subdiffusive regime is characterized by a continuously evolving dynamical critical exponent that diverges in a universal fashion upon approaching the melting transition from the thermal liquid side. Lastly, we describe the set of universal critical exponents that characterize this melting transition, which are connected by scaling relations.
We anticipate that these methods can be extended to treat a variety of other out-of-equilibrium dynamical phase transitions in disordered quantum systems. For instance, our analysis in terms of quantum percolation of resonant clusters can be straightforwardly adapted to treat higher dimensional delocalization transitions, which are so far inaccessible by other methods.  (3)

is a free-fermion
Hamiltonian with a fully integrable spectrum. Strictly in this non-interacting limit, we do not expect a delocalized phase for any disorder strength. To move away from this fine-tuned point, we imagine adding weak interactions that break integrability and produce a nonintegrable thermal liquid phase at high temperature and weak disorder. Such perturbations are irrelevant at strong disorder, and hence do not disturb the scaling properties of the quantum critical glass, which forms the basis of our effective model for the transition. In this appendix, we review the salient properties of SU(2) k anyons and derive the RSRG rules for renormalized couplings. We then analyze the main features of the resulting RG flow.

Hamiltonian and renormalization
For the purposes of this paper, SU(2) k may be viewed as a "quantum deformation" of SU(2) that serves to bound the spin size in a controlled way. The irreducible representations of SU(2) k are labeled by their "spin" j = 0, 1 2 , 1, . . . , k 2 , and obey modified fusion rules 'truncated at level k': j 1 ⊗ j 2 = |j 1 − j 2 | ⊕ · · · ⊕ min(j 1 + j 2 , k − (j 1 + j 2 )). Such SU(2) k algebras arise naturally in different contexts, but the one most relevant for our purposes is in describing the fusion of anyonic quasiparticles proposed as a basis for topological quantum computation [40]. We study one-dimensional chains of such 'anyons' [27,28] with random couplings [19][20][21], each with topological charge S = 1 2 . The Hilbert space H of the chain is then spanned by the basis states 1 2 , 1, . . . , k 2 } and j i is contained in the fusion product of j i−1 with S = 1 2 . Owing to the truncated fusion rules, H for N anyons cannot be written as a tensor product of local Hilbert spaces, in contrast to conventional SU(2) spins.
Since the RG rules can generate arbitrarily high spins, is necessary to work with generalizations of this truncated spin-1 2 chain that allow the spin S i on each site to take any value in { 1 2 , 1, . . . , k−1 2 }, and consider configurations in which j i ∈ j i−1 ⊗ S i . We work with the Hamiltonian introduced in Ref. 21, whereP S is the projector onto the fusion channel S in the fusion S i ⊗ S i+1 , that can be expressed in the basis introduced above using "F -moves" (to be defined shortly). The operatorQ is the same as the truncated dot product introduced in the main text. The different fusion channels are weighted by (A1) recovers the Heisenberg chain as k → ∞, since S 1 · S 2 = 1 2 S∈S1⊗S2 [S(S + 1) − S 1 (S 1 + 1) − S 2 (S 2 + 1)]P S . It is convenient to introduce a graphical depiction of the Hilbert space in which the fusion rules are simply represented by the merging of particle lines labeled by a particular spin representation. For SU(2) k , in contrast to more complicated theories, there is a single way to fuse an anyon into each of its outcomes. We may represent the 'link basis' states that span the Hilbert space pictorially by a trivalent graph whose 'legs' are the site spins, and whose links are labeled by the possible fusion outcomes of the spins to the left (Fig. A1). The projection operator onto a particular fusion channel of a pair of anyons has a simple graphical representation, also shown in Fig. A1, upto an undetermined constant that is fixed by normalization of the projector. Such fusion diagrams may be simplified by using a set of simple identities: the 'no tadpole' rule that dictates the vanishing of any graph with a closed loop that can be detached by a single bond; the 'loop rule' that allows us to replace a closed loops of anyon m by a factor equal to the quantum dimension d m ; and the 'F -move' that expresses how to relate the two independent fusion paths in combining three particles into a fourth, that involves an F -matrix (also called a 6j-symbol, see e.g. [21]). These rules are depicted graphically in Fig. A2.
The bond strengths J i are strongly random. The precise distribution of J is unimportant, but for concreteness we take the logarithm of couplings to be distributed ac- . The parameter W indicates the standard deviation in the order of magnitude of the couplings, J i , and is a convenient measure for the strength of disorder. We start with all spins- 1 2 , but will soon see that spins of arbitrary size are generated in the RSRG procedure.
For strong disorder (W 1), each eigenstate of the chain can be accurately constructed via a real-space renormalization group (RSRG) procedure that begins by identifying the strongest bond with strength J S and choosing a fusion channel for its two spins (Fig. A3a). For strong disorder, W 1, the neighboring bonds on the left (L) and right (R) are typically much weaker than the strong bond (J L/R J S ) and are therefore unable to dynamically alter the chosen fusion channel. If the strong bond spins fuse to a singlet (S = 0, k 2 ), then they drop out of future stages of the RG. Virtual excitations of this frozen singlet then mediate effective cou-plingJ ≈ A(S, S L , S R ) J L J R J S , between the neighboring spins. Here A is a J-independent function of spin sizes whose precise form will be unimportant at strong disorder. Another possibility is that the strong bond spins fuse to a "superspin" of size S = 0, k 2 , which continues to participate in later stages of the RG and interacts with its neighboring spins via renormalized coupling J L/R ≈ B(S, S L/R )J L/R . The coefficient B can take ei-  ther sign, but this will turn out to be unimportant for dynamics in highly excited states. Crucially, apart from the generating different spin-sizes, each RSRG step preserves the form of the Hamiltonian (A1).

Decimation rules
We emphasize that the general form of the decimation rules is entirely dictated by the structure of perturbation theory, and that the precise coefficients appearing in the renormalized couplings are not important at large disorder. Nevertheless, we give the explicit form of the decimation rules below for concreteness. This first part of our discussion very closely parallels that of [21], but we reproduce the details here for completeness; we caution the reader that we use a different notation, in particular differing on the definition of . . . k and {. . . } k .
In deriving these decimation rules, we consider the case k odd, where we can restrict ourselves to integer spins to make the calculation a bit simpler, but the final results hold for even k and half-integer spins as well. The easiest way to phrase the decimation rules is in terms of projection operators; in essence, we project the spins on a strong bond into a specific fusion channel and then compute the coupling of the fused spins with the remaining ones using perturbation theory. This procedure is most conveniently implemented using the graphical representation of the fusion procedure introduced above, and is especially transparent once one realizes that (A1) may be rewritten by inverting an F -move, so that the op-eratorQ i,i+1 represents the exchange of the '1' anyon between sites i, i + 1 (indeed, this is how the form of H was originally inferred [21].) The graphical depiction of the operatorQ i,i+1 ( rewritten as an anyon exchange) is given in Fig. A4.
We first explain how to handle first-order decimations. Consider the case where two spins S 1 and S 2 fuse to a new effective spin S; we wish to compute, for instance the effective couplingJ R between the spin S R on the right of S 2 and S in terms of J R , the old coupling between S R and S 2 (the results forJ L in terms of J L are similar and so we do not show them explicitly.) The fusion of S 1 and S 2 is captured by the projection operatorP 12 S , and the corresponding renormalized coupling is given bŷ P 12

SQ 23P 12
S (see Fig. A5(a)). By performing a (rather tedious) sequence of F -moves, we find that the coupling between the composite spin S and S R takes the form J RQSS R , with A4. TheQ-operator that defines the truncated dot product. The operator acts on nearest-neighbor sites and can be rewritten as a sum of projectors using an F -move. Note that it is nonzero only when the fusion rules are obeyed at its vertices.
In the k → ∞ limit, we recover which concurs with the result for conventional SU(2) Heisenberg spins, as expected.
Let us now imagine that the two spins S 1 = S 2 = j fuse to a singlet with topological charge (spin) 0. If we denote the adjacent spins on either side of the strongly coupled bond to be S L , S R , then we may write the 4spin Hamiltonian H 4 in suggestive form: H 4 = H 0 + H , where the 'bare' Hamiltonian has a strong bond coupling S 1 , S 2 , and H denotes the weak residual couplings, where Q 23 denotes the expectation value ofQ 23 in the state when the two spins fuse to the trivial state. Using second-order perturbation theory we find that the effec-  [21]. (a.) First-order decimation when spins S1, S2 fuse to a new spin S. In order to determine the effective Hamiltonian between the composite spin S and SR, we need to computeP 12

12
S ; using a a sequence of F -moves within the red box, we can show that this isJRQSS R , i.e. has the same form as the original Hamiltonian. The explicit form of JR is given in (A2). (b.) Second-order decimation, where spins S1 = S2 = S fuse to a singlet. Once again, a series of F -moves may be employed and after a tedious computation we can show that the effective Hamiltonian between SL, SR is J effQS L S R , with J eff given by (A6).
tive Hamiltonian between S L , S R takes the form as depicted graphically in Fig. A5(b). After a tedious calculation, we find that once again we may write the Hamiltonian between the residual spins in the original form, H eff = J effQS L S R , with In the Heisenberg limit, this becomes

Infinite temperature strong disorder fixed point
Iterating the renormalization procedure generates a tree of fusion choices (Fig. A3b), each branch of which corresponds to a different many-body eigenstate. Sampling these eigenstates equiprobably results in enormous simplifications that enable analytic progress. This averaging is analogous to working at infinite temperature, and for large systems is overwhelmingly dominated by states in the middle of the many-body spectrum (as there are exponentially more states there than at any other energy for entropic reasons). Such equiprobable sampling of eigenstates (tree branches) is accomplished by weighting each fusion choice (tree node) by the fraction of branches descending from this node. Specifically, two spins S 1 , S 2 are fused to spin S with probability proportional to the quantum dimension d S of the new spin.
This RSRG procedure generates superspins of various sizes, characterized by a distribution P(S), which quickly asymptotes to a steady-state form: . From P (S) it is straightforward to compute the probability of decimating a strong bond into a singlet (fusing to spin 0 or k 2 ), which asymptotes to Π s = 1 . This quantity plays a key role in determining properties of the strong disorder phase.
One can deduce most of the important properties of the strong disorder phase from the spin distribution P(S) alone. While spins initially clump into large super-spins, these super-spins eventually form pairwise singlet bonds. The resulting eigenstate has a strongly random pattern of singlets formed between super-spins, with a broad distribution of singlet bond sizes (see Fig. A3c). This random-singlet state has all the properties of a zero-temperature random quantum critical point: for example, it exhibits power law scaling of (disorder averaged) correlation functions [14,15], and critical scaling of (disorder averaged) entanglement entropy for sub-regions of size L: S(L) ∼ log L [9,41]. Remarkably, however, we find that every eigenstate in the spectrum of the anyon chains exhibits these quantum critical properties typically associated only with zero-temperature quantum phases. Hence, we see that the strong disorder phase of these random anyonic spin chains does not reach thermal equilibrium, and is not governed by the laws of thermodynamics. We refer to such a non-ergodic system as a "quantum critical glass".
A more detailed understanding of the strong disorder phase can be obtained by analyzing how the distribution of bond strengths evolves along the RSRG flow. For this purpose, it is convenient to define an energy scale for the RG, Ω (we choose units in which Ω ≡ 1 at the start of the RG), and to work with logarithmic variables Γ = log 1 Ω and β = log Ω J . Denoting the probability of having bond strength β at energy scale Γ by ρ(β, Γ) one can derive the evolution of ρ with energy scale Γ: The first term simply maintains the normalization of the bond probability distribution. The second, bracketed term in (A8) represents the probability that the strongly-bonded spins fuse to a singlet (∼ Π s ) or residual super-spin (∼ (1 − Π s )), thereby producing a renormalized bond-strengthJ = e −β/Γ . In (A8), we have neglected the spin-dependent prefactors A, B, as these are unimportant in the limit of strong disorder. Note that we have also assumed that the spin distribution is already converged to its fixed point value P (S) so that Π s does not depend on Γ.
Despite its complicated integro-differential character, (A8) is readily solved by the simple function where W is the initial disorder strength (log-standard deviation of J, or standard deviation of β). As the flow progresses, the effective disorder strength increases as W (Γ) = W + Π s Γ due to the bond renormalizations.
Hence we see that as long as the initial disorder is sufficiently strong, the system flows to an infinite-randomness fixed point where the RSRG procedure becomes asymptotically exact.
The quantities P * (S) and ρ(β, Γ) completely characterize the long-distance universal properties of the strong disorder phase. For instance, the singlet bonds between superspins are formed from clumps of lattice-scale spins of characteristic size ξ k ≈ e Π −1 s ; this defines an effective localization length for the random singlet phase, and characterizes the length scale at which the system's dynamics cross over from delocalized thermal behavior (L ξ k ) to quantum critical glassiness (L ξ k ). For large k, Π s ≈ 4π 2 k 3 indicating a very long localization length ξ k ≈ e k 3 /4π 2 . Since ordinary Heisenberg chains are recovered as k → ∞, this implies that Heisenberg chains are inevitably delocalized even for arbitrarily strong disorder. The dynamical scaling properties of the strong disorder fixed point are readily extracted from ρ(β, Γ). For example, if we wish to compute the residual coupling between two bonds separated by distance L, we simply perform the RSRG until these bonds are adjacent, and compute the typical coupling at this scale Γ. To estimate the scale at which two distant bonds become adjacent, we first compute the density, n, of surviving spins. The total number of spins, N , decreases by 2 (1) each time a strong bond is decimated to a singlet (new spin), so that implying a power-law decay of the number of spins with RG scale: n(Γ) = n(Γ 0 ) W is the exponent that characterizes scaling within the quantum critical glass phase. The typical distance between superspins at scale Γ is then L ≈ n(Γ) −1 . Inverting this relation, we find the typical energy scale associated with distance L 1: From this, we deduce that the typical coupling between bonds separated by length L 1 is as claimed in the main text (5). Hence the critical properties of the strong-disorder phase are characterized by the logarithm of energy scaling like a power law in distance, with characteristic exponent ψ. This power law scaling of log-energy is very different from the usual power law relation between energy and distance found in clean or weakly random critical points, but is typical for infiniterandomness fixed points [14,15,18,20,21]. While our RSRG procedure simplifies dramatically for highly excited states near the center of the many-body spectrum, results for lower-energy eigenstates in the tails of the spectrum with energy density E may be approximately obtained by following the ground-state RSRG procedure to scale Ω ≈ E, and then continuing with the highly-excited state RSRG procedure for lower energy scales. This argument suggests that energy density is a relevant perturbation to the ground-state fixed point, and that the dynamical scaling properties exhibit a crossover, with the infinite-randomness exponent evolving from the ground state value ψ GS = 1 2 to the excited state result (A11) on characteristic length scales L ∼ log 1/E. A more quantitative analysis of the intermediate energy spectrum requires numerical Monte Carlo sampling of the RSRG spectral tree, and is left for future work. In this Appendix we give a more detailed account of the prescription for numerically simulating the merging of resonances to analyze the effective model for the transition. We will focus on the example of critical glasses with ψ < 1 (in the language of SU(2) k anyon chains) for concreteness, but the results of this appendix hold for any value of ψ, and extend straightforwardly to MBL systems by setting ψ = 1. We begin by taking L anyons with bonds of strength J whose logarithm, β = log 1 J , is drawn from the distribution ρ(β, Γ 0 ) = 1 W e −β/W with disorder strength W . Since resonant links predominately form at the highest energy scales (since bonds at lower energy scales experience stronger effective disorder due to the RSRG flow), we restrict our attention to the energy range β ∈ [0, β max ], with cutoff scale β max ≈ O(1). We expect that changing the details of this cutoff procedure, or altering the precise initial form of the bond-strength distribution will not alter the results for universal scaling exponents, but rather only effect nonuniversal physics such as the precise value of the critical disorder strength.
In the first stage of resonance formation, we check for pairwise resonances between all pairs of bonds i and j for  which J(L ij ) > |J i − J j |, where J(L) is the typical effective coupling through the intervening (infinite temperature) random singlet phase given by (5), see also (A13), and L ij is the distance separating bonds i, j. Bonds that can resonate with multiple partners are lumped together into a single cluster.
After this initial stage, we iteratively check whether excitations can resonantly tunnel between any two newly formed clusters. Specifically, at any stage in the iterative process, we have a collection of N resonant chains labeled by j = 1 . . . N , each comprised of (m j + 1) bonds (chosen so that isolated bonds have m j = 0). Each resonant cluster has ≈ 2 mj energy levels, spread over bandwidth Λ j equal to the sum of couplings J for the resonant links inside the cluster, and hence have typical level spacing δ a ≈ Λj 2 m j . The energy levels within a given resonant chain are all strongly admixed and exhibit level repulsion.
If the residual coupling J(L ij ) between two clusters i, j can resonantly drive transitions among the levels, then we merge these two resonances into a new, larger one. The precise rule for whether or not to merge two resonances depends on the relative size of the bandwidths and level spacings for each. Consider two resonant clusters, labeled j = 1, 2. Without loss of generality, assume that cluster 1 has larger bandwidth Λ 1 > Λ 2 . There are two distinct cases to consider. If the level spacing of the first cluster exceeds the bandwidth of the second, δ 1 > Λ 2 (Fig. B1a), then the minimum energy difference associated with changing the states of both resonances is: On the other hand, if δ 1 < Λ 2 (Fig. B1b), then the minimum energy cost for transitions is of order ∆ 12 ≈ δ1δ2 Λ2 = Λ1 2 m 1 +m 2 . If J(L 12 ) ∆ 12 , then the residual coupling between the two resonant clusters merely weakly perturbs the energy levels of each. In this case, the combined energy levels of both clusters are well-described by the tensor product of the independent levels associated with each cluster indicating that energy excitations are well localized independently in each cluster. On the other hand, if J(L 12 ) ∆ 12 , then the energy levels of both resonances are strongly admixed such that excitations are spread throughout both clusters, and energy can resonantly shuttle between them. In our numerical procedure, we take ∆ as a sharp cutoff, merging resonances coupled by J(L ij ) > ∆ ij . Another significant feature is that the value of the critical disorder W c increases with increasing k. The numerical results for k ≤ 5 show a logarithmic growth of W c with k (though we cannot rule out more complicated large-k behavior as the numerics become increasingly demanding due to stronger finite-size corrections from the increasingly long range interactions at large k). Together with the k-dependence of ζ, this shows that Heisenberg spin-chains (recovered in the limit of k → ∞) are delocalized thermal liquids with diffusive thermal transport even for arbitrarily strong disorder, in agreement with previous analysis [9].
In terms of transport, we also find subdiffusive transport for other values of k with continuously evolving dynamical critical exponent z ∼ |W −W c | −ζ for all k. However, unlike ν, ζ exhibits strong k-dependence decreasing with larger k compatible with ζ = ψν, and tending to zero in the large k limit (Fig. C2). This shows that for large k, the delocalized side of the transition becomes increasingly conventional, with the dynamical critical exponent z tending rapidly towards 2 (ordinary diffusion).