Non-Abelian Symmetries and Disorder: A Broad Nonergodic Regime and Anomalous Thermalization

Ivan V. Protopopov, Rajat K. Panda, Tommaso Parolini, Antonello Scardicchio, Eugene Demler, and Dmitry A. Abanin Department of Theoretical Physics, University of Geneva, 1211 Geneva, Switzerland L. D. Landau Institute for Theoretical Physics RAS, 119334 Moscow, Russia The Abdus Salam ICTP, Strada Costiera 11, 34151 Trieste, Italy Scuola Internazionale di Studi Superiori Avanzati, Via Bonomea, 265, 34136 Trieste, Italy INFN, Sezione di Trieste, Via Valerio 2, 34127 Trieste, Italy Department of Physics, Harvard University, Cambridge, Massachusetts 02138, USA

Much of the progress in describing MBL and related phenomena is driven by the realization that fully MBL systems of, e.g., spins or fermions on a lattice exhibit a new kind of robust emergent integrability [9][10][11][12]14]. Specifically, it is a complete set of quasilocal integrals of motion (LIOMs) that underlies the ergodicity breaking in MBL phases. The LIOM construction naturally explains the area-law entanglement of the MBL eigenstates [11,32], logarithmic entanglement growth in a quantum quench experiment [7,[11][12][13]27], and a number of other dynamical properties of the MBL phase [33][34][35].
It quickly becomes clear that distinct MBL phases are possible. Much like in the theory of thermodynamic phase transitions, the symmetry of the system plays a central role. For example, systems with a discrete Z 2 symmetry [36][37][38] can exhibit two distinct MBL phases: In one of them, the eigenstates spontaneously break the Z 2 symmetry, and, in the other, the symmetry is preserved. Both phases can be described using the LIOM theory.
Disordered systems with continuous non-Abelian symmetries, which constitute a broad and experimentally relevant class, show a qualitatively different behavior. An example of such a system is a disordered, SUð2Þsymmetric spin chain. Crucially, a non-Abelian symmetry such as SUð2Þ is inconsistent with all eigenstates obeying area-law entanglement entropy. Thus, conventional MBL with a complete set of LIOMs is forbidden by symmetry in this case [39,40]. Some integrals of motion must become nonlocal; accordingly, the entanglement entropy of a subsystem in a typical, highly excited eigenstate must scale at least logarithmically with that subsystem's size l in 1D systems, S ent ðlÞ ≳ c logðlÞ (where c is a coefficient of the order of one).
The incompatibility of conventional MBL with an SUð2Þ symmetry stems from the fact that symmetry imposes exact degeneracies in the spectrum, leading to resonances. The eigenstates of a system of size L come in multiplets with an exact degeneracy 2S 0 þ 1, where the total spin S 0 grows as ffiffiffi ffi L p (see, e.g., Ref. [40]). When two subsystems of size L are connected, even very weakly, the eigenstates of the large system of size 2L are entangled superpositions of the ∼ ffiffiffi ffi L p × ffiffiffi ffi L p original multiplet states; at the same time, by symmetry they must be the eigenstates of the total spin operator of the large system. This situation is qualitatively different from conventional MBL, where two subsystems, when connected, get entangled only weakly, according to the area law.
The fact that an SUð2Þ symmetry enforces a minimum amount of entanglement in the eigenstates raises several fundamental questions. What is the nature of the excited eigenstates and the corresponding dynamical properties of disordered, SUð2Þ-symmetric systems? One exciting possibility hypothesized in Refs. [40,41] is that at sufficiently strong disorder a new kind of dynamical, nonergodic phase may emerge-characterized by an entanglement entropy of excited eigenstates that is subthermal but scales faster than the area law [e.g., as S ent ðlÞ ∼ c logðlÞ]. Such a phase would display only a partial set of LIOMs, being distinct from the conventional MBL phase. Another, equally intriguing possibility is that thermalization may be inevitably enforced by such symmetries in the thermodynamic limit [42]. If this is the case, it would be highly desirable to understand the microscopic processes that govern thermalization, as well as the corresponding time and length scales. This topic has been attracting strong interest, and several works provide valuable complementary insights into the above questions. Reference [42] studies random SUð2Þ k anyonic chains, arguing that the breakdown of the strong-disorder, real-space renormalization group (SDRG) approach as k → ∞ signals self-thermalization of SUð2Þsymmetric spin chains. Reference [41] computes the noise spectrum of random Heisenberg chains using an SDRG approach applied to excited states. Reference [40] introduces a toy model, in which eigenstates of an SUð2Þsymmetric spin chain are described by regular tree tensor networks with S ent ðlÞ ≳ c logðlÞ entanglement entropy scaling; they study the stability of such eigenstates under local perturbations of the Hamiltonian, finding indications of eventual slow delocalization. Furthermore, Refs. [43,44] and Ref. [45] consider spin dynamics in disordered Hubbard and t − J models, respectively. They find that spins are not localized in the conventional sense even at strong disorder and numerically study spin transport, finding indications of subdiffusive behavior. While transport does not imply ergodicity, this behavior is another signal that, in the presence of non-Abelian symmetries, a localized phase cannot have plain-vanilla MBL phenomenology.
In this paper, we study SUð2Þ-symmetric disordered spin chains, focusing on their spectral properties, and the properties of highly excited eigenstates. Our main results can be summarized as follows. First, via a combination of extensive exact-diagonalization studies of system sizes up to L ¼ 26, we establish the existence of a nonergodic regime distinct from MBL. Second, to establish the ultimate fate of such systems, we develop a modified SDRG procedure for excited states. The key innovation is the inclusion of multispin processes. We analyze the number and structure of resonances, finding that they become relevant at large system sizes. This result signals that at L → ∞ SUð2Þ-symmetric systems thermalize but via unconventional, multispin processes.
The starting point of our analysis is the SDRG which is used to approximately construct excited eigenstates. This procedure, originally introduced to describe low-energy properties of random spin chains [46][47][48], has been recently applied to the highly excited eigenstates in a range of systems [9,36,41,42]. Applied to the random Heisenberg chains, the SDRG yields a caricature of an eigenstate in the form of an (irregular) tree tensor network, with the structure that depends on the disorder realization; at each step of this construction, two spins which are strongly interacting with each other (relative to their interactions with their other neighbors) are added to form some other total spin. This caricature is illustrated in Fig. 1. Naturally, such states are strongly nonergodic, although distinct from the conventional MBL eigenstates, e.g., in their entanglement properties (see below). So, if the SDRG procedure remains accurate, the system is in a novel nonergodic phase. However, typically, the SDRG procedure allows one only to test for "local" resonances involving a small number of nearby spins, and therefore it is an open question when and whether SDRG is reliable and gives a good approximation of the system's eigenstates at large system sizes.
Below, we investigate how well the SDRG procedure approximates the system's eigenstates. To that end, we first perform extensive numerical simulations of spectral statistics and eigenstates. In particular, we test the eigenstate thermalization hypothesis (ETH), which is believed to underlie thermalization in ergodic systems [49]. We find that, at strong disorder, there is a broad nonergodic regime in which the SDRG accurately captures the eigenstates. At weak disorder, above a certain length scale, we find evidence for thermalization and breakdown of the SDRG. We investigate how this length scale depends on the strength of the effective disorder, in the regime where it is smaller than the largest system size accessible numerically (L ¼ 26). We note that previous work [44] focuses on exact diagonalization (ED) studies of transport in small systems, where, for a particular distribution of spin exchange couplings, subdiffusive behavior is found and interpreted as a signature of ergodicity. Our conclusions differ from that work in that we focus on eigenstates, which instead provide a clear indication of ergodicity breaking.
To describe the behavior of large chains, far beyond those accessible via conventional numerical techniques, we develop a novel approach to describe nonlocal, multispin processes that are not captured by the conventional SDRG.
For that purpose, we analyze the relevance of terms in the Hamiltonian that are responsible for the processes that are usually neglected in the SDRG. These terms mix different states in the SDRG, and if the mixing is sufficiently strong, they cannot be neglected and give rise to resonances. We study how the number of resonances grows with the system's size and describe their properties, such as energy scales and the number of physical spins involved.
We find that at strong disorder resonances are absent in a surprisingly broad range of length scales, signaling a regime in which the SDRG describes eigenstates accurately. Eventually, in sufficiently large systems, resonances proliferate, the perturbation theory in the terms neglected by the SDRG does not converge, and the system thermalizes. We expect this process to give rise to full ergodicity, in an unconventional way that we describe. Thus, our conclusion favors the scenario of "non-Abelian-symmetryprotected thermalization." Our work shows that this thermalization proceeds via long-range resonances that involve many spins; we extract the corresponding timescales and find them to be extremely long at strong disorder. Thus, for all practical purposes, the strongly disordered system appears nonergodic, for reasonably short experimental observation times.
The rest of the paper is organized as follows: In Sec. II, we introduce the model and describe the SDRG procedure which is used to find approximate eigenstates. In Sec. III, we first use ED and a measure of participation ratios to check how well the approximate eigenstates given by the SDRG agree with the exact one. Then, we investigate the onset of the ETH and its breakdown at strong disorder using various measures (level statistics, statistics of matrix elements, and entanglement entropy). In Sec. IV, we develop our SDRGbased approach to the analysis of resonances. We show how the terms neglected in the SDRG give rise to resonances which eventually proliferate, leading to the thermalization of very large systems. Section V closes the paper with a recapitulation and suggestions for future work.

II. STRONG-DISORDER RENORMALIZATION GROUP AND TREE STATES
We start this section by introducing the model of a disordered Heisenberg chain in Sec. II A. We refresh the well-studied example of a random-field Heisenberg model which lacks the SUð2Þ symmetry and compare it to the symmetric Heisenberg chain. We then review the construction of the approximate eigenstates based on the The strong-disorder renormalization group aims to construct approximate eigenstates. It yields a tree state, characterized by its geometry and the choice of total block spins at each node (see the main text). The Heisenberg Hamiltonian, written in this basis, gives rise to processes which can change the spins along the "causal" path connecting two neighboring spins, to one of the block spins. (c) A schematic dynamical phase diagram of the random Heisenberg model. There are three regimes: (I) At short length scales L < L 1 ðαÞ, the SDRG tree states are accurate approximations of the eigenstates; (II) at intermediate length scales L 1 ðαÞ < L < L erg ðαÞ, there are resonances but the system remains nonergodic; (III) above some large length scale L > L erg ðαÞ, the resonances proliferate and the system becomes thermalizing (see Secs. IV B and IV C for a definition of these scales).
SDRG [46][47][48] paradigm. Basic properties of tree states obtained by the SDRG are discussed (Sec. II B). Finally, in Sec. II C, we qualitatively describe our approach to probing the stability of tree states obtained by the SDRG. The detailed numerical studies are presented in the subsequent sections.

A. The model and preliminary remarks
The model we study is the disordered, Heisenberg spin-1=2 chain with the Hamiltonian where S i ¼ ðS x i ; S y i ; S z i Þ are the standard spin operators and the couplings J i are independent random variables with a probability distribution PðJÞ specified by two parameters, η and α. The parameter 0 < η < 1 gives the fraction of antiferromagnetic (positive) couplings in the system. Throughout this work, we assume η ¼ 0.5. We do not expect the properties of highly excited eigenstates in the middle of the many-body band to exhibit a significant dependence on the choice of η. Note that the ground state properties do not depend strongly on the value of η (except for the extremal points η ¼ 0, 1-see Ref. [50]).
The parameter α > 0 controls the distribution of jJj. We assume that the probability density function of this distribution has a power-law form with a cutoff at jJj ¼ 1: where ΘðxÞ stands for the Heaviside function. This distribution of couplings emerges naturally in a wide range of low temperatures, as shown in the seminal papers [46,47].
In that context, under the assumptions of what would now be called the ETH, it can explain the anomalous exponent of the specific heat observed in early experiments [51]. The exponent α effectively controls the strength of disorder, with smaller α corresponding to stronger disorder. Indeed, for the distribution (2), the ratio of two neighboring couplings in the system has a typical value maxðjJ 1 j; jJ 2 jÞ minðjJ 1 j; jJ 2 jÞ typ ≡ exp ðhj ln jJ 1 j=jJ 2 jjiÞ ¼ e 1=α : ð3Þ This ratio increases exponentially when α → 0. Therefore, at small α it becomes more and more likely to find exchange constants in the system that are much larger than the two neighboring ones. This ratio is exactly the condition that enables SDRG, as we discuss below.
Another quantity of interest is the smallest coupling J (in absolute value) in the whole system, representing the "weakest link." We find that For α ¼ 0.3 and L ≃ 20, this coupling can be as small as 10 −3 hJi (here, hJi is the mean value of J i ). Throughout the paper, it is helpful to contrast our findings with the properties of the random-field XXZ model, which is studied extensively in the literature (see Refs. [28,52,53] for recent reviews): The model (5) can be mapped, via Jordan-Wigner transform, onto an interacting fermionic problem with t representing the hopping amplitude, U the nearest-neighbor interaction, and h i the random on-site potential with a variance that we denote by W. In the following, for concreteness, we assume that t ∼ U, such that the disorder strength is described by a single dimensionless parameter L W ¼ W=t. Then, the XXZ model is known to have a diffusive-subdiffusive dynamical transition [54] at L W ≃ 0.55 and an MBL-thermal transition at L W ≃ 3.5 [15]. Note that, in the limit of strong disorder W ≫ t, the parameter L W can be interpreted as a typical distance between (rare) pairs of "resonant" sites in the model that happen to have close enough values of the magnetic field to enable resonant spin exchange (or, equivalently, hopping in the fermionic model). A resonance between spins 1 and 2 appears, e.g., if jh 1 − h 2 j ≲ t. Starting at a very large disorder, these resonant sites are typically well separated by distances of OðL W Þ, and one can show that they will not mix at any order of perturbation theory [55,56]. By mixing, we mean that the resonant pairs can exchange energy and become strongly entangled in the eigenstates. The fact that resonances are rare and isolated at L W ≥ 3.5 is intimately related to the low, area-law entanglement scaling of eigenstates and the existence of a complete set of LIOMs [9][10][11][12]14]. As the disorder strength is decreased, eventually the resonant pairs of spins become mixed, forming a connected network; then, LIOMs are destroyed, becoming nonlocal, and the system exits the MBL phase.
Below, we find that SUð2Þ-symmetric spin chains exhibit a qualitatively different delocalization scenario. First, in Sec. IV, we compute the system size L 1 ðαÞ at which the system typically has just one resonance. We find the scaling L 1 ðαÞ ∝ α −0.4 , such that L 1 ðαÞ diverges at strong disorder. An MBL-type scenario, if realized, would mean that at sufficiently small α, when L 1 ðαÞ is very large, resonances remain isolated and the system is in a nonergodic phase. We show below that, in contrast to this naive expectation, multispin resonances, which involve changing spins at different levels of the RG (see Fig. 1), play an important role. In contrast to the MBL case, the density of resonances grows nonlinearly with L, which can be attributed to the nontrivial entanglement patterns of the SDRG tree states. As the RG proceeds, the probability of composite spins being involved in such resonances grows, eventually leading to the restoration of ergodicity. The length scale that controls it, denoted by L erg ðαÞ > L 1 ðαÞ, is estimated below. While at weak disorder L erg ðαÞ can be extracted from ED studies using standard indicators of ergodicity (such as the ETH tests for matrix elements), at stronger disorder it is well beyond the system sizes reachable by ED, and we estimate it by performing a detailed analysis of resonances (see Sec. IV).

B. SDRG and excited eigenstates of the Heisenberg chain
In this subsection, we qualitatively describe the SDRG approach to the disordered Heisenberg chains and discuss the properties of tree tensor-network states that it yields. A detailed description of the method is provided in the Appendixes. We emphasize that such states differ from the conventional MBL ones in two crucial aspects: First, they have a parametrically larger entanglement entropy, and, second, one cannot define a complete set of LIOMs for them.
A very large typical ratio of two neighboring couplings found for small α [Eq. (3)] suggests that the properties of the system can be described using the SDRG framework. The idea of SDRG is to identify a local "grain" in the system that is strongly coupled inside but, due to strong disorder, only weakly coupled to the rest of the system. The state of the grain is then approximated by one of the eigenstates of its Hamiltonian, with the rest of the system decoupled. If one is looking for the ground state, the eigenstate of the grain is chosen to be its ground state. Alternatively, if one aims to construct a random highly excited eigenstate that is effectively at an infinite temperature, as we do in this paper, some eigenstate of the grain is randomly chosen. Then, the effective Hamiltonian of the system in which the grain is in the chosen eigenstate (or, more generally, a multiplet of states if symmetries dictate degeneracies in the spectrum of the grain's Hamiltonian) is calculated by the perturbation theory in the grain-system coupling.
One can continue this procedure, assuming that the disorder in the effective Hamiltonian remains strong. Such is indeed the case for, e.g., ground states of random antiferromagnetic (AFM) Heisenberg chains [48]. Then, a repeated application of the SDRG rules results in an approximate wave function of the whole system, obtained by "patching" together the wave functions of the grains.
A detailed discussion of the SDRG rules for excited states of the Heisenberg chain can be found in Ref. [41], and we provide it in Appendix A. Qualitatively, for this system a grain is a pair of neighboring spins coupled by a strong bond; its eigenstates [which come in SUð2Þ multiplets] are labeled by the total spin of the grain. The SDRG procedure replaces such spin pairs by effective (typically larger) spins; i.e., it assigns some total spin to larger and larger blocks of contiguous spins in the system. The resulting approximation for an eigenstate (more precisely, for a degenerate symmetry-enforced multiplet [57]) is a kind of a tree tensor network, illustrated in Fig. 2. The nodes of the tree represent the block spins identified in the SDRG process. The structure of the tree reflects the order in which the elementary spins of the system should be added up to give an (approximate) eigenstate.
The fusion of spins in the course of the SDRG must be supplemented by a perturbative account of the interaction of merging spins with the rest of the system. In the present setting of the infinite-temperature SDRG, where spins typically fuse into nonsinglet states, a first-order perturbation theory (that simply amounts to the projection of the fusing spins onto the direction of the total spin) suffices in most cases. The resulting renormalization of couplings is weaker than the one that occurs in the low-temperature SDRG for AFM spin chains, where the spins always fuse into singlets, and, therefore, a second-order perturbative treatment is required to find new renormalized couplings (see Appendix A). Still, the distribution of couplings developed in the course of SDRG turns out to be broad (see Ref. [41] and below).
Within the SDRG approximation, the values of the block spins (the numbers associated to the nodes of the tree in Fig. 2) label the eigenstates of the Hamiltonian and bear a similarity to the LIOMs of the conventional MBL phase. An eigenstate of the Hamiltonian is also an eigenstate of a sequence of these operators, just as an eigenstate of an MBL Hamiltonian is simultaneously an eigenstate of FIG. 2. A multiplet of eigenstates predicted by the SDRG for a system of 12 spins 1=2. The leaves of the tree represent elementary spins in the system. The tree describes the way the elementary spins are fused into larger block spins in the course of the SDRG. The numbers in the nodes indicate the resulting spins of the blocks. The value in the top node (marked red) is the total spin S 0 of the system (S 0 ¼ 1 in the present example). S 0 is an exact integral of motion. ð2S 0 þ 1Þ different states in the multiplet can be distinguished by additionally specifying the projection of the spin in the top node to the z axis. each LIOM. However, there are two major differences between these quantum numbers and LIOMs.
First, in the MBL phase, the eigenstates of H are at the same time eigenstates of a fixed set of LIOMs. Total spins of the blocks in our problem would form conserved operators if different eigenstates are represented by geometrically identical trees, which differ only in the values of the block spins. In reality, the order in which the spins are merged in the course of the SDRG depends not only on the particular disorder realization, but also on the eigenstate of the grain, which is randomly picked at any given step of the SDRG (see Appendix A for details). Thus, the values of the block spins, in general, cannot be promoted from labels of a particular eigenstate to operators acting in the full Hilbert space. The structure of larger blocks depends on the history of choosing total spins at the earlier steps of the SDRG.
Second, LIOMs in an MBL system are quasilocal, exponentially localized in space operators [10][11][12]. In contrast, the block spins of the strongly disordered Heisenberg chain have a hierarchical structure. While some of them (living near the bottom of the tree) can be expressed in terms of an Oð1Þ number of the original spin operators S i , the other ones, found at the higher levels of the tree, are highly nonlocal in terms of the original spins. Thus, SUð2Þ symmetry forces some integrals of motion to become nonlocal. Therefore, the SDRG (in the regime of its validity) describes a nonergodic phase of a new kind, with a partial rather than complete set of LIOMs. Our goal is to investigate the stability of this putative phase.
The novel nonergodic character of tree eigenstates manifests itself in the scaling of entanglement entropy. For simplicity, we consider the entanglement entropy of an eigenstate with respect to a half-chain cut: where ρ L=2 is the reduced density matrix of half the chain in the chosen eigenstate and the trace is taken over the degrees of freedom (d.o.f.) in the other half of the system. A bound for the entanglement entropy depends on the tree structure describing a given state, in particular, on the tree depth d (the number of levels between the very top node of the tree and the original physical spins). We find, via numerical simulations, that typical states produced by the SDRG procedure have a logarithmic depth d ∝ ln L. It is then possible to show (see Appendix B) that the entanglement entropy of a single typical [58] tree satisfies where c 1 and c 2 are numerical constants of the order of unity that depend on the statistical properties of the tree. Thus, the entanglement of the tree states scales faster than the area law found in MBL but significantly slower compared to the thermal entanglement of an infinitetemperature state, S th ðL=2Þ ≈ ðL=2Þ (measured in bits). The upper bound on the entanglement entropy in Eq. (7) can also be generalized (see Appendix B) to the case when the state in question is not a single tree state but rather a linear combination of n T tree states: Although this bound might seem weak, it has an important implication, which is used below: If the system's eigenstates become ergodic, they must be represented by an exponentially large number of tree states.

C. Validity of SDRG and the (in)stability of tree states
The SDRG is a heuristic procedure relying on strong disorder. The tree states generated by the SDRG are not exact eigenstates of the Heisenberg spin chain, but how accurate are they? Historically, at each step of the SDRG, one checks that the disorder in the effective Hamiltonian remains strong, such that strong couplings can be found; one can then check for the absence of resonances involving a small number of spins, to make sure that the neglected processes do not destroy the tree structure. While for the analysis of ground states this process is often sufficient, it is unclear whether such tests can guarantee the accuracy of SDRG for the excited states.
Below, we check the validity of the SDRG for excited states using several approaches. First, we compare SDRG tree states to the exact eigenstates for system sizes up to L ¼ 26, obtained numerically. We use a number of measures, such as level statistics, and the ETH and its breakdown. Second, to describe large system sizes, we develop an approach to account for many-body processes that are usually neglected in the SDRG and to test their relevance. We introduce this approach qualitatively now, and we apply it in what follows. Suppose that the SDRG yields some tree state jΨ 0 RG i, specified by the tree geometry and the choice of total spins in each node. Instead of considering the effective Hamiltonian at every step, we can write the original Hamiltonian exactly in the basis of tree states with the geometry identical to that of jΨ 0 RG i. The first key observation is that the selection rules imposed by symmetry facilitate the analysis of relevant processes; more specifically, the block spins along a path connecting a pair of contiguous physical spins to the top of a tree can change by ΔS ¼ 0; AE1 [59]. The second observation is that, given that the typical spins of larger blocks grow (as the square root of the block size), the tree states connected to jΨ 0 RG i by the Hamiltonian are expected to have the same geometrical structure. This expectation is because, for large spins, strong bonds remain strong when a value of some spins is changed by ΔS ≪ S.
We search for resonances between different tree states and characterize their properties. Solving the full eigenvalue problem for large L is hopelessly complicated; thus, we focus on low-order resonances. Effectively, we check whether the Hamiltonian hybridizes a given tree state with its neighbor, say, jΨ 1 RG i (a neighbor is a state such that hΨ 1 RG jHjΨ 0 RG i ≠ 0). As long as the probability of finding resonances is sufficiently low, we expect that true eigenstates are localized in the tree basis, which corresponds to a nonergodic phase or regime (if it occurs only for sufficiently small system sizes). Alternatively, if there are many resonances which proliferate, it is natural to expect the SDRG to break down and the system to become ergodic.
It is instructive to draw parallels with the conventional MBL phase of the strongly disordered XXZ spin chain in a random magnetic field. The caricature of MBL eigenstates is just product states with a well-defined S z i projection for each spin. While corrections to this picture certainly exist (e.g., LIOMs are not strictly equal to S z i operators), we know that MBL is stable if the disorder is sufficiently strong. Our aim is to understand whether for SUð2Þ-symmetric chains tree states, with their built-in correlations and unusual entanglement properties, can be stable, representing a dynamical phase distinct from both MBL and ergodic phase.
Below, we use the above approach to reveal a broad nonergodic regime where tree states are, in fact, stable. We also provide evidence that trees eventually become unstable above a certain system size (dependent on α) for all values of α that we study. We therefore propose the picture that, while for finite systems the dynamics is nonergodic at strong disorder, in the thermodynamic limit the ETH should be recovered (see Fig. 1).

III. EXACT DIAGONALIZATION STUDIES
In this section, we present our numerical results from exact diagonalization.

A. Probing the stability of tree states
To analyze the accuracy of the SDRG procedure, we first study the participation ratios of exact eigenstates of the system (1) in the basis of the tree states generated by the SDRG. More precisely, for a given disorder realization fJ i g, we first run the SDRG to generate some tree state jΨ 0 RG i with a total spin S 0 . A complete basis of states in the sector with a given total spin S 0 (and some fixed z projection of the total spin) can be built out of jΨ 0 RG i by fixing the geometry of the underlying tree but allowing the block spins in the tree (apart from the top one, S 0 ) to take all possible values consistent with the angular momentum addition rules. We denote the basis obtained in this manner by where L is the length of the chain and D S 0 ;L is the Hilbert space dimension of the sector with a total spin S 0 and a fixed projection S z ¼ 0: The state with index a ¼ 0 is the original SDRG state jΨ 0 RG i. In general, due to the correlations between the geometric structure of the tree and the values of the block spins discussed in Sec. II B, many of the states in the basis (with indices a > 0) would not be approximate eigenstates constructed in the real-space renormalization group. We expect, however, that at strong disorder the geometry of the tree that corresponds to the state jΨ 0 RG i is also appropriate for a number of other SDRG states that do not differ too much from jΨ 0 RG i in the values of block spins. Therefore, a significant fraction of jΨ a>0 RG i are, in fact, "SDRG eigenstates." We then perform an exact diagonalization of the Hamiltonian in the basis (9) [60]. Among all the eigenstates of the Hamiltonian, we focus on a single one, denoted by jEi, that has the maximum overlap with a given jΨ 0 RG i. The quality of jΨ 0 RG i as an approximation to jEi can be quantified by the inverse participation ratio (IPR) of the state jEi in the SDRG basis (9): and its inverse N E ≡ 1=I E , which can be viewed as the number of tree states jΨ a RG i (of a given topology) that one needs to accurately represent the eigenstate jEi. Thus, small values of N E ∼ 1 indicate that the SDRG is accurate, while very large N E ≫ 1 signals an instability of tree states.
Computing the participation ratio N E for 10 3 disorder realizations fJ i g (and a single random SDRG state jΨ 0 RG i for each fJ i g), we investigate the statistical properties of this quantity. We perform numerical simulations for the disorder parameter α ranging from α ¼ 1.2 (weak disorder) to α ¼ 0.3 (strong disorder). The results are summarized in Figs. 3 and 4. Figure 3 shows several examples of the distributions of log 10 N E for different system sizes and two different disorder strengths. We observe that in short systems, L ¼ 10, jΨ 0 RG i is very close to an exact eigenstate even for weak disorder, α ¼ 1, in the sense that N E ∼ 1. Upon increasing the system size, N E grows, signaling that approximating the eigenstate jEi with a tree state jΨ 0 RG i becomes less accurate.
The evolution of the typical value of N E (defined as e hln N E i ) with the system size is illustrated in the top in Fig. 4. Interestingly, even in the weak disorder regime α ¼ 1 and for the largest system size L ¼ 20, the typical N E ∼ 25 remains small compared to the dimension of the Hilbert space D S 0 ;L . The latter depends on the spin sector S 0 , which is chosen at random in the present analysis.
The SDRG procedure we use generates states with different S 0 in accordance with their probability in the infinitetemperature ensemble, PðS 0 Þ ∝ ð2S 0 þ 1ÞD S 0 ;L . For L ¼ 20, the most frequently encountered value of S 0 is 3, corresponding to the Hilbert space dimension D 3;20 ¼ 38760. Moreover, for 90% of the SDRG states S 0 ≤ 5 and D S 0 ;20 ≥ 10659. The length dependence of the typical Hilbert space fraction occupied by the energy eigenstate jEi (in the tree basis), e hln N E =D S 0 ;L i , is shown in the bottom in Fig. 4.
It is instructive to compare the above findings to the behavior of the IPR in the product-state basis for the conventional MBL phase. Viewing MBL as a kind of Anderson localization in the Hilbert space, one might naively expect that in the MBL regime the eigenstates would exhibit a system-size-independent IPR, N E ≳ 1. It is known [15,32,52,61], however, that in reality MBL eigenstates are rather fractal when viewed in the product-state basis: the participation ratio N E scaling as N E ∝ D γ ∝ 2 γL with an exponent that depends on disorder strength. The fractal behavior stems from perturbative corrections and resonances discussed at the end of Sec. II A (or, equivalently, it is due to the fact that local integrals of motion have support over more than one lattice site). In the strong-disorder limit, γ ∝ t=W ¼ 1=L W ≪ 1. The MBL transition is thus marked not by the emergence of the growth of N E with the system size but rather by a jump of the exponent γ to its thermodynamic value γ ¼ 1 (at an infinite temperature).
The behavior shown in Fig. 4 for the Heisenberg chain is qualitatively similar. At strong disorder α ≤ 0.8, the dependence N E ðLÞ for the available system sizes can be approximated by an exponential fit, N E ∝ 2 L=L 1 ðαÞ [see dashed lines in the top in Fig. 4; the corresponding values of the fitting parameterL 1 ðαÞ are indicated in the legend]. The lengthL 1 ðαÞ grows as the disorder strength is increased. By analogy with the conventional MBL, we can expect that the length scaleL 1 ðαÞ characterizes the density 1=L 1 ðαÞ of the rare local resonant d.o.f. in the system; see also the discussion at the end of Sec. IV B.
At weaker disorder α ¼ 1, 1.2, the naive exponential fit produces a very smallL 1 ðαÞ < 2.5. Moreover, the slope d ln N E =dL of the corresponding lines shows a clear increase as the system size grows. Accordingly, the fraction N E =D L;S 0 (bottom in Fig. 4) displays a tendency toward saturation, suggesting that, ultimately, the scaling of N E in long systems becomes ergodic, N E ∝ 2 L .
In view of the results above, it may be tempting to conclude that the strongly disordered Heisenberg spin chains do indeed display a nonergodic, non-MBL phase with unusual treelike eigenstates that are only slightly dressed by perturbative corrections and occasional resonances (similar to how in the conventional MBL phase the eigenstates are perturbatively dressed product states). At weaker disorder, one would then expect a transition into an ergodic phase. However, the crucial question concerns the ultimate fate of the putative nonergodic behavior in the thermodynamic limit. In particular, does the observed fractal scaling N E ∝ 2 L=L 1 ðαÞ persist, or does it eventually cross over to the ergodic scaling, as for the weakly disordered case? In order to answer these questions, in the next sections we subject the hypothetical nonergodic phase to several stringent tests.

B. Level statistics
Our main goal in this subsection is to further characterize the nonergodic behavior found above and its dependence on the system size. We employ the standard diagnostics of ergodicity and its breakdown: the level statistics in the center of the many-body band. An extensive use of the constraints imposed by the SUð2Þ symmetry allows us to perform exact diagonalization on spin chains of up L ¼ 26 spins. Larger system sizes that we can achieve here compared to Sec. III A are due to the use of massively parallel algorithms together with the possibility to focus on a small number of eigenstates near the band center (recall that the identification of the eigenstate jEi studied in Sec. III A requires the knowledge of the full set of eigenstates). In most of our studies, we concentrate on the S 0 ¼ 0 sector, and data for S 0 ¼ 1, 2 do not show any qualitative differences. For each L, α, and each disorder realization fJ i g i¼1;…;L , up to 50 eigenstates around the middle of the spectrum (fewer for L ¼ 10, 12 and L ¼ 26) are obtained, and a total of at least 1000 disorder realizations (except for L ¼ 26) are considered.
We characterize the level statistics by the r parameter, defined as follows [6]: with Δ n and Δ nþ1 being two consecutive level spacings. The distribution of the parameter r and its dependence on the system size and disorder strength are shown in Figs. 5 and 6. The distributions of r change qualitatively as α is decreased at a fixed L: For largest α ¼ 1.6 (very weak disorder), r is described by the standard Wigner-Dyson distribution, while for small α ¼ 0.3 (strongest disorder considered), one observes the Poisson distribution, with virtually no level repulsion. This result supports the existence of a nonergodic regime at accessible system sizes. For α ∈ ½0.6; 1, the level statistics is intermediate between the Wigner-Dyson and Poisson distributions. We also illustrate the dependence of the distribution PðrÞ on the system size for weak disorder α ¼ 1. 3. It is evident that the distribution flows toward Wigner-Dyson, albeit relatively slowly.
Furthermore, we study the flow of the average value, hri, with the system size, in an attempt to extract relevant length scales. hri as a function of L for different values of α is illustrated in Fig. 6. For weak disorder α ≥ 0.8, the dependence of hri on L is nonmonotonic. Our data show a tendency toward the Poisson statistics for small system sizes L < L Ã ðαÞ, but for L > L Ã ðαÞ the value of hri starts growing, moving toward the Wigner-Dyson (WD) value. Upon decreasing α to the value of 0.8, the length scale L Ã ðαÞ increases, while its value r Ã ðαÞ ≡ r½L Ã ðαÞ decreases. The ultimate flow of hri toward the WD value is consistent with the expectation that at weak disorder the system becomes ergodic for modest system sizes. One can estimate the scale where the system becomes ergodic, L erg , by extrapolating the hrðLÞi dependence up to the crossing with the WD line. The length scale extracted in this way is larger than the maximum system sizes accessible numerically for α < 1. The extrapolation procedure suffers from a large uncertainty. Therefore, we choose instead to characterize the delocalization crossover by the length L Ã ðαÞ, and we expect that L erg ðαÞ ∝ L Ã ðαÞ. We extract the dependence L ⋆ ðαÞ by fitting the data in Fig. 6 through the quadratic dependence rðLÞ ¼ r ⋆ ðαÞ þ a½L − L ⋆ ðαÞ 2 (with a being an additional fitting parameter). The resulting behavior is illustrated in Fig. 7.
The data at stronger disorder α ∈ ½0.3; 0.6 show prima facie a qualitatively different behavior. For the strongest disorder α ¼ 0.3, the parameter hri slowly increases for small L, in stark contrast with the behavior found for α ≥ 0.8. Interestingly, at small L this parameter is below the Poisson value of hri P ≈ 0.39. We attribute this result to strong disorder leading to the appearance of very small couplings in a typical disorder realization (smaller than the level spacing at small L). The chain is then effectively broken into smaller, almost noninteracting, spin chains. This breaking leads to level clustering, and the r parameter becomes sub-Poissonian. However, since the level spacing decreases exponentially with the system size, while the weakest coupling decreases only as a power law [see Eq. (4)], the level clustering is eventually washed out, and for L > 18 the parameter hri rapidly approaches the standard Poisson value. For disorder strengths α ¼ 0.45, 0.6, hri is initially slightly above the Poisson value, but it decreases as the system size is increased; no flow toward WD is seen. For the system sizes analyzed, it is evident that ergodicity has not developed and a single SDRG tree state provides a good approximation to the eigenstates, as we also demonstrate in the previous subsection.
The exact diagonalization results for strong-disorder values α ∈ ½0.3; 0.6 may be consistent with two scenarios. One scenario is that (much like in the usual MBL) the system experiences a phase transition at some critical disorder strength. Another scenario is that, even at strong disorder, the system eventually flows to ergodicity, similar to what we find for weaker disorder values. Assuming that this second scenario is realized, in large enough systems the curves for α ¼ 0.45, 0.6 first develop a minimum and then flow to the WD value at yet larger system sizes. The corresponding scale L Ã can be heuristically extracted by fitting the ED data shown in Fig. 6 against a quadratic polynomial, as mentioned above, and then performing an extrapolation to the numerically inaccessible system sizes. The dependence of the length L Ã on disorder, as extracted by the analysis outlined above, is illustrated in Fig. 7(c). It is consistent with a power-law scaling L Ã ðαÞ ∝ α −1. 4 . We note that the curves hrðLÞi for α ≥ 0.6 (including extrapolated data for α ¼ 0.45, 0.6) can be collapsed (in the vicinity of L ¼ L Ã ) into a single one by simultaneously rescaling r → r=r Ã and L → L=L Ã ; see Fig. 7.
To sum up, the length L Ã beyond which the spectral parameter starts flowing toward the WD value (but the system of size L Ã is still nonergodic, because r Ã is closer to the Poisson value) grows rapidly with the increase of disorder. Although the trend is clear, we are extrapolating significantly away from the accessible system sizes L ≤ 26. Thus, the law governing L Ã ðαÞ which we propose should be taken with a grain of salt. In the next subsection, we proceed to test the eigenstate thermalization hypothesis.

C. Eigenstate thermalization hypothesis and its breakdown
We characterize the eigenstates of random Heisenberg chains by testing the ETH and its breakdown. The ETH  provides a microscopic picture of thermalization in ergodic quantum systems [62][63][64]. Specifically, it states that individual ergodic eigenstates appear to be thermal, from the point of view of simple physical observables (e.g., fewbody operators). The ETH formalizes and extends the intuition that the eigenstates of an ergodic system should be "as random as possible," up to a small set of global constraints (in our case, energy and total spin). For our purposes, the ETH can be formulated in terms of the expectation values of local observables. LetÔ be an operator representing some physical observable. Then, for every pair of eigenstates jai; jbi of an ergodic system, the ETH yields an ansatz for matrix elements ofÔ [62,63]: where E and ΔE are, respectively, the average and the difference between the energies of the two eigenstates, and SðEÞ is the microcanonical entropy. Function fðE; ΔEÞ is a smooth function of its arguments, which reflects dynamical properties of observableÔ and is system specific. Finally, R ab is a normally distributed random variable with unit variance. Notably, in this formula the diagonal partŌ is assumed to be a smooth function of E alone and is equal to the microcanonical average ofÔ. This assumption reflects the fact that observables in eigenstates are equal to their microcanonical ensemble values.
According to Eq. (13), in a thermalizing system the distribution of values hajÔjai for eigenstates jai that are sufficiently close in energy should display a reasonably smooth dependence on E, with only small, normal fluctuations about the average, suppressed exponentially in the system size by the factor e −SðEÞ=2 so as to reproduce the microcanonical ensemble in the infinite-size limit.
We focus our attention on the following two local observables:Ô where ði ⋆ ; i ⋆ þ 1Þ is the pair of spins coupled the most strongly (jJ i ⋆ j ¼ max jJ i j) and ðj ⋆ ; j ⋆ þ 1Þ is its antipodal pair (j ⋆ ¼ i ⋆ þ L=2 modulo L). Since couplings are independent, the latter pair is coupled by an interaction J j ⋆ of typical (or "random") strength, hence the nameÔ rand .
Let us discuss our expectations for the averages of these operators over eigenstates, depending on whether the SDRG is accurate. First, suppose that jai is exactly an SDRG tree state. Then the spins ði ⋆ ; i ⋆ þ 1Þ are going to be paired in either a S ¼ 0 or a S ¼ 1 state, and the value of hajÔ max jai is going to be either −3=4 or 1=4, respectively. Even for hajÔ rand jai, these two values are going to be likely, although in many cases the pair ðj ⋆ ; j ⋆ þ 1Þ will not be coupled directly by the SDRG procedure but rather at a higher level, resulting in some intermediate value.
However, in the ergodic regime-when the SDRG breaks down-local thermalization implies that the local state of any pair of spins will be a uniform (at T ¼ ∞ ) mixture of the four possible above-mentioned states, resulting in a thermal average of zero for both observables.
The distributions of the expectation values ofÔ max =rand over eigenstates at system size L ¼ 20 are shown in Fig. 8. It is clear that the system is perfectly compliant with the ETH at sufficiently high values of α, whereas at smaller values of α the behavior is consistent with the eigenstates being close to tree SDRG states. This phenomenology, which we interpret as a finite-size crossover between ergodic and nonergodic structure of the system's eigenstates, is compatible with the observed behavior for the level statistics (cf. Sec. III B).
In order to validate our interpretation, we characterize the finite-size flow to ergodicity by looking at the percentage of eigenstates whose corresponding values ofÔ max =rand fall within some fixed window centered at zero. Figure 9 confirms that the "ergodic fraction" of infinite-T eigenstates is increasing with L for bothÔ max andÔ rand , though much more slowly for strong disorder. Crucially, at disorder α ¼ 0.6, the ETH is still strongly violated, which is consistent with the nonergodic behavior observed in level statistics above.

D. Entanglement entropy
Another witness of nonergodic behavior can be found in the scaling of the half-chain entanglement entropy with the system size, which is known to obey an area law for MBL systems and a volume law for ergodic ones. More precisely, in a system that thermalizes, generic eigenstates are expected to be similar to random states; their entanglement entropy equals the thermodynamic one, yielding for states in the middle of the band a value of S ent ðL=2Þ ¼ L=2 þ oðLÞ, when measured in bits [65,66].
The numerical results are reported in Fig. 10. The median entanglement entropy of the infinite-temperature eigenstates exhibits linear scaling for all considered values of α, but the linear coefficient observed at L ≤ 20 deviates substantially from the ergodic prediction at strong disorder (although significant curvature is present). This result is once more consistent with the results of Secs. III B and III C.
To summarize the results of this section, ED data show a clear trend toward the ETH for moderate-to-weak disorder (i.e., α ≳ 0.6) while indicating a novel nonergodic regime for the case of strong disorder. To determine the behavior of the system in the thermodynamic limit, we have to resort to a completely different approach, presented in the next section, which surpasses ED.

IV. RESONANCE COUNTING: FROM A SINGLE TREE TO A FOREST
As we showed in the previous section, at strong disorder, finite-size random Heisenberg chains exhibit a nonergodic regime, in which their eigenstates are well approximated by tree states. Here, to determine the eventual fate of these systems in the thermodynamic limit L → ∞ , we develop an approach to analyze resonances between different tree states. We are able to capture long-range, multispin processes, which are beyond the conventional SDRG. We obtain the asymptotic behavior of the resonance number and their spatial structure. We find that the resonance density grows for all studied disorder strengths, leading to an eventual delocalization at very large length scales, which we estimate. Beyond this length scale, the system presumably becomes ergodic.
Given a tree state generated by the SDRG, we can construct a complete basis [Eq. (9)] in the Hilbert space (with the total spin of the system fixed) by allowing the values of the block spins identified by the SDRG to take all possible values consistent with the rules of angular momentum addition. The Hamiltonian (1), written in this basis, will then connect the initial SDRG state to a certain number of other tree states. We consider the eigenvalue problem in this basis. The localization in this problem corresponds to true eigenstates being close to the tree states; in contrast, delocalization signals the breakdown of the SDRG approximation, suggesting ergodicity. The criteria for delocalization is studied below.

A. Connectivity of the hopping problem
First, we investigate the connectivity of this eigenvalue problem. That is, we analyze how many matrix elements of the Hamiltonian between a given tree state and other ones are nonzero. The SUð2Þ symmetry of the model imposes stringent constraints on the matrix elements of the Hamiltonian [40]. Specifically, let us consider one of the terms in the Hamiltonian, J i S i S iþ1 . It can be shown that the action of such an operator on a tree state can affect only the block spins that lie on the path in the tree connecting spins i and i þ 1; see Fig. 1(c).
Moreover, each of those block spins on the path, if affected, can change only by 0 or AE1. It then follows that FIG. 9. Fraction of eigenstates with a value of hÔ max =rand i between −1=8 and 1=8, for α ∈ f0.6; 0.8; 1.0g and S 0 ¼ 0. The ETH predicts this fraction to become 1 in the infinite-size limit. the number of states connected to a given one by the operator J i S i S iþ1 is given by where l i;iþ1 is the length of the path in the tree connecting physical spins i and i þ 1. The factor 3 arises from the selection rules: The operator can change the value of the representation at a node by ΔS ¼ −1; 0; þ1. The ≃ sign is due to the constraint that the new values of block spins in the tree must still be consistent with the rules of angular momentum addition (in particular, they cannot be negative). Sufficiently far from the bottom of the tree, the typical values of block spins are large, and the latter constraint can influence only the prefactor in Eq. (14). Taking into account that the Hamiltonian (1) is just a sum of local terms of the form discussed above, we conclude that the total connectivity in the Hilbert space induced by the Hamiltonian (1) is We are now left with the task of computing the distribution PðlÞ of these lengths l i;iþ1 , for the SDRG trees. The SDRG fuses spins that are most strongly coupled. Neglecting the correlations between the (renormalized) couplings at any step of the SDRG as well as the dependence of those on the couplings at earlier stages of the SDRG, we can assume that the pair of spins to be fused is just randomly chosen among all possibilities (with the only requirement that the fusing spins are nearest neighbors so that locality is respected). In such an ensemble of maximally random SDRG trees, the distribution PðlÞ can be computed analytically. As we show in Appendix D 1, it turns out that PðlÞ falls down exponentially with l, and in the limit L → ∞ it becomes [the normalization is the correct one considering l ≥ 2, so P l≥2 PðlÞ ¼ 1]. With this distribution PðlÞ, the sum (15) is dominated by the maximum l M over the L terms. To the leading order in L, the value of l M can be estimated from the condition LPðl M Þ ∼ 1. This estimate follows from the distribution of the largest of L random variables with the distribution (16), which is given by Plugging back into K ∼ 3 l M , we find with κ ¼ ln 3= lnð3=2Þ ≃ 2.71. The power-law scaling (18) and the value of the exponent κ are in good agreement with the numerical simulations of the SDRG trees (with the full set of SDRG rules taken into account); see Fig. 11.

B. Local resonances
Our next goal is to find resonances among the K ∝ L κ "hopping" processes generated by the Heisenberg Hamiltonian for a given tree state. We first focus on investigating relatively small system sizes L ≲ 30, comparable to those accessible by exact diagonalization.
To study the number of resonances, we first use the SDRG procedure to generate a random tree state jΨ 0 RG i and identify resonant neighbors for that state-that is, the ones for which the ratio of matrix element connecting them to jΨ 0 RG i and the energy difference is larger than one. These resonances invalidate the perturbative expansion around the "infinite disorder" eigenstates (the SDRG states). Their proliferation signals the instability of tree states, strongly suggesting that ergodicity is restored. The SDRG is essentially a local optimization procedure that aims to construct basis states free of such resonances. Based on the results presented above, we expect that, at strong disorder and in relatively short systems, these resonances should be few in number, because the SDRG is accurate.
The average number of resonant neighbors, hK res i, of an SDRG tree state is shown in Fig. 12. We observe that, for relatively small systems discussed here, hK res i scales linearly with the system size L. As expected, the slope of this linear growth becomes smaller for stronger disorder. The condition hK res i ¼ 1 defines an important (disorderdependent) length scale in the problem, L 1 ðαÞ, at which resonances start appearing. Naively, this length scale plays the same role as the length scale L W introduced in Sec. II A to characterize the resonances in a random-field XXZ chain. We find that the scale L 1 ðαÞ grows at stronger disorder, crudely following a power-law dependence, L 1 ðαÞ ∝ α −0.4 . Figure 12 shows that at relatively strong disorder values α ≤ 0.6, the average number of resonant neighbors for an SDRG tree state is 1 or less for all system sizes available in ED. This number agrees with the observation that such chains display a nonergodic behavior in all of the ED studies of Sec. III, with eigenstates being well approximated by the tree states.
In particular, the low number of resonances is in agreement with the slow growth of N E (the participation ratio of eigenstates in the tree basis) found in Sec. III A. Drawing parallels to conventional MBL systems, it is tempting to identify the length scaleL 1 ðαÞ that controls the exponential growth of N E with the system size (see Sec. III A) with L 1 ðαÞ. However, the comparison of the values of L 1 ðαÞ and L 1 ðαÞ reveals that the latter is several times shorter. We attribute this difference to the effect of the second-order perturbative corrections that contribute to the spreading of the exact eigenstate jEi over the SDRG tree state. Such higher-order perturbative corrections lie beyond the firstorder resonance counting that underlies the scale L 1 ðαÞ. The perturbative corrections are expected to be more significant at weak disorder; in accordance with this intuition, we find a more significant difference between L 1 ðαÞ andL 1 ðαÞ for such disorder strengths.

C. Longer systems and the proliferation of resonances
Does the linear scaling of K res with the system size discussed in the previous subsection persist in the thermodynamic limit? Such behavior closely resembles that of the strongly disordered XXZ model. It implies that the resonant neighbors can be attributed to the existence of local subsystems with resonating levels which, if sufficiently separated in space, remain isolated and do not cross talk (in the sense that there is no significant entanglement in the eigenstates between such "local" resonances). If true, this argument would be strongly in favor of the SDRG tree states surviving in an infinitely long system, up to corrections due to local, isolated resonances. We now perform a detailed analysis of resonances in large systems, up to L ∼ 2 × 10 3 , and find that Heisenberg chains actually behave qualitatively differently compared to the plainvanilla MBL systems: The number of resonances grows faster than linearly with the system size.

Number of resonances and their structure
The probability for an SDRG tree state to have no resonant neighbors vanishes in sufficiently long systems (see the top in Fig. 13). Then, a typical tree state has a large number of resonances attached to it. The bottom in Fig. 13 shows (in log-log scale) the dependence of the typical number of resonant neighbors (defined as e hln K res i ) for an SDRG tree state. The dashed lines represent power-law fits: with the scale L 1 ðαÞ determined from the short-scale behavior of K res ; see Sec. IVA. The "anomalous" exponent μ is approximately disorder independent, μ ≈ 0.38. The details of the numerical procedure employed to find and characterize resonances are given in Appendix C. The power-law scaling of the number of resonant neighbors, K res ∝ L 1þμ , implies that, in stark contrast to the random-field XXZ model, the density of the resonating d.o.f. grows with the system size. Accordingly, at least some of the resonant transitions must originate not from the rearrangement of a few local spins but rather involve a growing number of spins. To support this conclusion, we analyze the structure of typical resonances. The top in Fig. 14 shows the average number of block spins, N bs , that are changed in the course of a resonant transition. We observe that N bs grows (albeit rather slowly) with the length of the chain. The decrease of N bs with increasing disorder can be understood as follows: At weak disorder, the possibility to change a block spin in a resonant manner is often accompanied by an "instability" (with respect to resonances) of the block spins higher up in the hierarchy. The more complicated (involving flipping of more than one block spin) resonant neighbors appear and contribute to the increase of the average N bs . On the other hand, at strong disorder an instability of a single block spin is more likely to remain "localized" and not to "propagate" upwards in the tree. The size of a typical resonance in real space also grows as the system size is increased; see the bottom in Fig. 14. An elementary physical spin is affected by a resonant transition if at least one of its descendant block spins changes its state. The physical size of a resonance is then defined as the total number of the elementary spins involved in it. Essentially, it is the level (as counted from the bottom of the tree) of the highest block spin affected by the resonance that sets the size of the resonance.
We observe that at moderate system sizes the typical spatial size of a resonance at strong disorder exceeds that at weak disorder. This observation is in accord with our intuition: At strong disorder, a large number of spins need to rearrange collectively in order for a transition to be resonant. In terms of the SDRG, this need means that many SDRG steps can be performed before the resonances start to play any role. On the other hand, in sufficiently long systems, we see the opposite tendency: Weakly disordered chains typically exhibit resonances of larger size. This result is the manifestation of the propagation of an instability of block spins upwards in the tree; cf. the discussion above.
The growing length scale characterizing the resonances comes together with a decreasing energy scale. The latter is  given by a typical matrix element for a resonant transition, V typ . Its system size dependence is shown in Fig. 15. In the following subsection, we use V typ to estimate the energy scale associated with the crossover to ergodicity.

Breakdown of SDRG and delocalization
The results presented in the previous subsection (most importantly, the power-law growth of the resonance density) strongly suggest that, even in the strongly disordered chains with α ≤ 0.6, where ED studies of Sec. III reveal little (if any) sign of ergodicity, the resonant transitions missed by the SDRG eventually proliferate. In this subsection, we estimate the corresponding thermalization scale L erg ðαÞ.
Given an SDRG tree state and a set of resonant transitions associated with it, one can identify a set of block spins that can be changed via at least one resonant process. We refer to those block spins as resonant, or unstable, ones. For a chain of L spins, there are 2L − 1 nodes in the SDRG tree (L of them are leaves corresponding to the physical spins). At each stage of the SDRG procedure, some number L RG of the block spins play the role of the physical spins of the system. For example, in the initial state of the SDRG, L RG ¼ L and the SDRG spins are just the physical ones. The final stage of the SDRG corresponds to L RG ¼ 1 and the top node of the SDRG tree being the only remaining spin. The ratio L=L RG is nothing but the average size of the spin clusters in the system.
At any given moment in the course of the SDRG, only the unstable block spins that are among the L RG spins currently comprising the system are relevant for potential delocalization. The others either have not yet formed or have already been decimated by the SDRG; they are not expected to contribute directly to the physics at the current energy scale. It is thus natural to ask how the number and density of the resonant block spins evolve in the course of the SDRG.
This result is illustrated in Fig. 16, which shows the dependence of the number (top) and the density (bottom) of unstable block spins for two different values of α and several values of the physical chain length L. These quantities are plotted as a function of the running RG length L RG , normalized by L. For not too small L RG =L, we observe that, for a fixed disorder strength, the density of resonant spins exhibits a universal (L-independent) behavior, ρ res ðL RG =L; LÞ ¼ ρ res ðL RG =LÞ. The density ρ res ðL RG =LÞ is higher at weaker disorder and also grows in the course of the SDRG. In contrast, at small L=L RG (corresponding to the final stages of the SDRG in a finite chain), a rather pronounced dependence on L is observed.
Next, let us denote by ρ max ≡ ρ max ðα; LÞ the maximum density of unstable spins developed during the SDRG process. Small ρ max means that ρ res remains small at all steps of the SDRG. We then expect the resonances to be of little importance for our system. On the contrary, ρ max ∼ 1, indicates that at some stage the SDRG inevitably runs into a state where almost all spins participate in resonances. Then, the basic assumptions of SDRG are violated and we expect it to break down-this breakdown means that block spins are no longer well defined and start resonating. Presumably, this resonation signals the onset of ergodicity.
It is natural to assume that there exists a critical value ρ max;c < 1 at which the crossover between nonergodic (SDRG valid) and ergodic (breakdown of the SDRG) regimes occurs. We can then identify the length of the system for which ρ max ¼ ρ max;c , as the ergodicity scale: The scale L erg ðαÞ along with the typical matrix element for resonant transitions V typ gives an estimate for the ergodicity time and energy scales: Figure 17 shows the dependence of ρ max ðLÞ for different disorder strengths. Estimating L erg ðαÞ requires fixing the critical density ρ max;c . While we have no general theory for ρ max;c , we observe (see Fig. 17) that the value   Fig. 17.
It is evident that, at strong disorder α ¼ 0.3, resonances start to proliferate only at very large length scales L erg ≈ 300, and, moreover, the corresponding timescales are extremely long. Such timescales are beyond the limitations of the synthetic platforms where ergodicity and its breakdown are actively investigated (see Ref. [28] for a review). Thus, in experiments, strongly disordered, SUð2Þ-symmetric systems are expected to display the novel nonergodic regime described above.
Systems of size L ≫ L erg are slowly thermalizing and presumably display slow diffusive transport at low frequencies. An interesting open question concerns the eventual fate of the integrals of motion obtained in the first steps of the SDRG (when the typical cluster size is much smaller than L erg ). Such nearly conserved operators arise due to strongly coupled cluster spins, and therefore destroying them would typically involve a relaxation process with a large energy scale ΔE. In very large systems, slow thermalizing processes eventually destroy the conservation of these operators. However, since thermalization processes typically occur on a much smaller energy scale, E erg ≪ ΔE, we expect that the decay time of such operators will be parametrically large in ΔE=E erg . An instructive example is that of a narrowbandwidth thermal bath with energy scale E 0 ; there, the relaxation of excitations with energy ω ≫ E 0 is exponentially slow in ω=E 0 [67]. We expect that the integrals of motion obtained within the SDRG before its breakdown will be similarly long-lived (but we leave a detailed investigation of this issue for a followup work). Thus, we propose a picture that the dynamical properties of the strongly disordered Heisenberg chains are captured by the SDRG at frequencies ω ≳ E erg (in particular, they have nontrivial noise properties, described in Ref. [41]). At lower frequencies, ω ≲ E erg , a crossover to a diffusive behavior is expected.

V. CONCLUSIONS AND OUTLOOK
To sum up, the goal of this paper is to investigate the effects of continuous non-Abelian symmetries on dynamical properties of disordered systems. We have considered a concrete example of disordered Heisenberg spin chains, which are characterized by an SUð2Þ symmetry. To describe the properties of this model, we combined state-of-the-art exact diagonalization studies with a new approach that allows us to include long-range resonances into the strong-disorder renormalization group.
We have found that, in a broad range of disorder strengths and system sizes, Heisenberg chains exhibit a new kind of nonergodic behavior. In this regime, the highly excited eigenstates have a scaling of the entanglement entropy that is intermediate between the area law characteristic of MBL states and the volume law found in thermalizing systems. This behavior stems from the tree tensor network structure of the eigenstates obtained within the SDRG. Simultaneously, in this regime the system exhibits a different kind of integrability, with integrals of motions having a varying degree of locality: Some of them act on a small number of neighboring spins, while others act on larger and larger spin clusters.
Furthermore, we found that, for weak disorder, the behavior crosses over from nonergodic to ergodic as the system size is increased. For stronger disorder, all system sizes accessible numerically exhibited nonergodic behavior. To address the eventual fate of the nonergodic phase in this case, we have extended the SDRG approach, characterizing resonances that endanger the stability of tree states. Our results strongly suggest eventual delocalization and ergodicity, albeit at very large system sizes; delocalization occurs via unconventional, multispin processes, which is yet another unique feature of disordered systems with non-Abelian symmetries. In future work, we plan to describe the transition between nonergodic and ergodic regimes as a function of the system size. A promising starting point seems to be to formulate an effective model in terms of resonant d.o.f., with parameters extracted using the methods described above.
Another interesting direction is to better understand dynamical signatures of the new nonergodic regime uncovered here. One natural experiment would be to probe the dynamics of the most local integrals of motion (e.g., total spin of a pair of strongly coupled physical spins) and to observe that, for system sizes L < L erg , it is conserved to a good precision and for arbitrarily large times. Another interesting open question concerns spin transport in N-species, disordered Hubbard models [21,43]. In case of flavor SUðN ¼ 2Þ symmetry, our work suggests that a sufficiently large system should show thermalizing behavior. Further work is required to establish the details of the dynamics (e.g., diffusion vs subdiffusion).
More broadly, this work sets the stage for the future discovery of new nonergodic regimes and true dynamical phases that survive in the thermodynamic limit. The approach introduced here can be naturally extended to other symmetry groups, for example, SUðNÞ spins. We expect that adapting our method to disordered systems with discrete non-Abelian symmetries should be possible as well, but it will require a modified renormalization scheme. Another possible application is to establish the fate of the transition between distinct MBL phases in systems with Z 2 symmetry [36][37][38]. We leave studies of these questions for future work. Even more generally, it would be interesting to investigate the stability of other tree tensor network structures with intermediate entanglement scaling as possible good approximations of eigenstates in physical systems.

APPENDIX A: SDRG FOR HEISENBERG SPIN CHAINS
The SDRG procedure for Heisenberg spin chains is formulated and discussed comprehensively in a number of publications [41,[46][47][48]50]. As is common for the RG studies, the aforementioned works focus on the flow of the system parameters under RG transformation and the consequences of this flow for the thermodynamic properties of the system. The interpretation of the SDRG approach from the perspective of the many-body eigenstates (including highly excited ones) is put forward in Ref. [36]. In this Appendix, we briefly review the SDRG protocol for SUð2Þ-symmetric Heisenberg spin chains with the emphasis on this latter aspect of the problem. We also discuss several subtle points of the procedure.
The SDRG protocol we design deals with the spin Hamiltonian of the form In the initial state of the SDRG, the spin operators S i represent the elementary spins 1=2 that constitute the system. The summation runs over nearest-neighbor links in a 1D lattice, i ¼ 1; …; L − 1. The SDRG procedure merges individual spins into clusters. Correspondingly, at later stages of the SDRG, S i represent the total angular momentum of clusters of elementary spins (block spins in the terminology of Sec. II B). The corresponding quantum number S i can take arbitrary integer or half-integer values limited from above by half of the size of the cluster. The eigenstates and eigenvalues of each of the "link" Hamiltonians H i are completely fixed by symmetry. Its spectrum consists of n i ≡ 2 minðS i ; S i þ 1Þ þ 1 levels with energies where by jSj we denote the absolute value of the spin, jSj ≡ ffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi ffi SðS þ 1Þ p , andS i ¼ jS i − S iþ1 j; jS i − S iþ1 j þ 1; …; S i þ S iþ1 stands for the total spin of the cluster formed by the spins S i and S iþ1 .
Every link i in the system is thus associated with a set of energy gaps in the link Hamiltonian H i : The gaps Δ AE i;S i have the physical meaning of the precession frequencies for the vector S i − S iþ1 in the state of the ith link characterized by the total spinS i .
The SDRG procedure aims to eliminate from the system the fastest d.o.f. Therefore [41], it looks for the link i 0 with maximal value of min AE Δ AE i;S i and approximates the state of the link i 0 by the one with definite total spinS i . Thereby, we eliminate from the consideration the rapidly oscillating vector S i 0 − S i 0 þ1 .
Note a subtle point here: At any stage of the SDRG, each link in the system is generically characterized by a set of n i − 1 > 1 energy gaps, and the judgement on which link represents the strongest-coupled subsystem requires a guess about the total spinS i associated to each link. Our present situation is to be compared with the SDRG for the ground state of the Heisenberg spin chains [50] or the SDRG for the highly excited states of less symmetric systems [36]. In both these cases, the relevant energy gap for each of the links is uniquely defined either as the gap between the ground state and the first excited state of the system or just due to the fact that each link is associated with a two-dimensional Hilbert space and is characterized by a single energy gap to begin with. This fact allows one, in particular, to apply the SDRG for the construction of a full basis of (approximate) eigenstates in, e.g., the generalized quantum Ising model of Ref. [36].
The dependence of the definition of the strongestcoupled subsystem of the spinsS i makes it difficult to generate the full set of eigenstates by the SDRG procedure [68]. We do not attempt to solve this problem here. Instead, in our numerical analysis, we resort to the probabilistic sampling of the SDRG tree states in the middle of the many-body band. To this end, we argue that in the infinitetemperature ensemble the probability for a couple of spins to have the total spin S is dominated by the entropic factor (2S þ 1). Therefore, starting from the initial SDRG state with S i representing the elementary spins for each link in the system, we generate randomlyS i with probability pðS i Þ ∝ ð2S i þ 1Þ. We then use Eqs. (A2) and (A3) to find the gaps associated to the links and identify the link i 0 to be removed by the SDRG. The link i 0 is removed, and the pair of spins S i 0 and S i 0 þ1 is replaced by a new spinS i 0 (see Fig. 18). Before the removal, the spins S i 0 and S i 0 þ1 have spins S i 0 −1 and S i 0 þ2 as their left and right nearest neighbors, respectively. The corresponding links i 0 − 1 and i 0 þ 1 have the spinsS i max −1 andS i max þ1 associated to them. After the removal of the link i 0 , the spins S i 0 −1 and S i 0 þ2 become the neighbors of the spinS i 0 . We check at this point if the values ofS i 0 −1 andS i 0 þ1 are still consistent with the rules of the angular momentum addition for the new spin configuration. If this is the case, we keep them as the spin values associated to the newly created links (see Fig. 18).
Otherwise, the newly created links receive new randomly generated values of the associated spins.
An important ingredient of the SDRG procedure is the renormalization of spin-spin couplings. In the zeroth order of the perturbation theory, the two strongly interacting spins S i 0 and S i 0 þ1 are treated as decoupled from the rest of the system. To establish the coupling of the newly created block spinS i 0 to the outside world, one needs to take into account the higher-order terms of the perturbation theory in the interactions J i 0 −1 and J i 0 þ1 (see Fig. 18). Specifically, we consider the 4-spin Hamiltonian and integrate out fast fluctuations of S i 0 − S i 0 þ1 . In a generic case, the first-order treatment suffices, and we end up with the effective Hamiltonian after an SDRG step: Here and below, to simplify our notations, we denote the spinS i 0 simply by S; the shorthand notation v stands for On going over from the Hamiltonian (A4) to the Hamiltonian (A6), we simply project the spin vectors S i 0 and S i 0 þ1 on the direction of the (approximately) conserved spin S. Such an approximation is not sufficient, however, if the spins S i 0 , S i 0 þ1 and S form a "quantum Pythagorean triangle," i.e., satisfy one of the two conditions One of the couplings J 0 i 0 AE1 turns then to zero cutting the chain into two independent pieces, and the perturbation theory should be developed further.
A particular case of Eq. (A8) is the singlet formation: In that situation, the second-order perturbative Hamiltonian takes the form [50] H ð2Þ eff ¼JS i 0 −1 · S i 0 þ2 ; ðA9Þ Note that the singlet formation is the only instance of Eq. (A8) relevant in the context of the SDRG near the ground state. In a more general case of S ≠ 0, straightforward but lengthy algebra leads to where [69] Here, we denote by ϵ αβγ the Levi-Civita tensor and Note that, while the explicit expression for g αβ is complicated, its tensor structure is fully determined by the SUð2Þ symmetry. Equation (A11) shows that in the general case of S ≠ 0 the form of the Hamiltonian (A1) is not preserved under the SDRG transformation if the second-order terms are taken into account. We argue, however, that, while being of central importance for the SDRG flow for the ground states of antiferromagnetic chains, the second-order renormalizations play only a minor role for the physics at an infinite temperature. The reason for this result is the growth of spins under the SDRG transformation and the fact that in the set of all triples ðS; S i 0 ; S i 0 þ1 Þ the "Pythagorean" ones have measure zero. Of course, in the spin chains consisting of spins 1=2 that we consider in the present work, the application of the second-order perturbation theory is quite frequent in the early stages of the SDRG, since about 1=4 of the microscopic spins may fuse into singlets. However, as mentioned above, the second-order treatment is controlled in this case. In contrast, the nontrivial Pythagorean triples ðS; S i 0 ; S i 0 þ1 Þ with S ≠ 0 appear only at the later stages of the SDRG when the typical spins are already large; such triples are extremely rare. Numerically, we find 2-5 occurrences of nontrivial Pythagorean triples in a typical SDRG tree state for a system of 10 4 spins.
Correspondingly, instead of treating Eq. (A11) in its full form, we apply to it several (generically uncontrolled) approximations. First, we focus on the coupling of spin S to that of the spins S i 0 −1 and S i 0 þ2 for which the corresponding first-order coupling in Eq. (A6) vanishes. For example, in the case of the plus sign in Eq. (A8), leading to the vanishing of J 0 i 0 −1 in Eq. (A6), we make a replacement Explicitly, using Eq. (A12) (and taking into account that v ¼ S 2 in the present case), we find Finally, we use the expected value of the spinS i 0 −1 (see Fig. 18) to estimate the last term in Eq. (A15) in a kind of mean-field approximation according to After the manipulations outlined above, the Hamiltonian H ð2Þ eff reduces back to the Heisenberg model expected by the SDRG. We stress that, despite being uncontrolled, the approximations we employ are expected to produce a correct order-of-magnitude estimate for the coupling of spins S i 0 −1 and S (that vanishes in the first-order perturbation theory). This estimate should be enough to capture the physics at an infinite temperature; see the discussion after Eq. (A13).
Before closing this section, let us stress once again that in the present work we are primarily interested in the properties of wave functions generated by the SDRG: SDRG tree states. Each tree generated by the SDRG describes the way the elementary spins in the system fuse to organize an approximate eigenstate of the Hamiltonian [or rather an SUð2Þ multiplet thereof]. Given an SDRG tree, one can use the Clebsch-Gordan coefficients to write down the corresponding wave function in terms of the elementary spin d.o.f. We illustrate this process for the two SDRG trees shown in Fig. 19.
In a system of two spins S 1 and S 2 , a wave function with total spin S 0 and the z projection of the total spin M 0 , Here, jS i ; M i i, i ¼ 1, 2, is the states of the spin S i with the z-axes projection M i . In a similar manner for a system of free spins S 1 , S 2 , and S 3 , the wave function corresponding to the tree of Fig. 19

APPENDIX B: ENTANGLEMENT ENTROPY OF TREE STATES
In this Appendix, we discuss the entanglement properties of the tree states.
Let us consider a single tree state jΨi in a system of L spins 1=2; see Fig. 20. We are interested in the entanglement entropy where ρ L=2 stands for the density matrix of, e.g., the left half of the system.
To estimate S ent ðL=2Þ, we observe that the Schmidt cut in the middle of the chain naturally gives rise to a cut of the tree representing the state into a "forest" and a decomposition of the chain into a collection of clusters in the manner exemplified in Fig. 20. We denote by L i (R i ) the clusters to the left (right) from the cut. It can be readily seen that with the whole system in the state jΨi the quantum state of each of the clusters described above lies in the multiplet specified by the subtree build above that cluster. In particular, all the clusters have well-defined total spin. The only d.o.f. for each cluster that is not locked by the state jΨi is the projection of its total spin. It follows then that the rank of the density matrix ρ L=2 is limited by where the product runs over all the clusters to the left of the cut and S L i are corresponding total spins. Each of the spins S L i is limited by L=2, while the number of clusters cannot exceed the depth d of the tree. Correspondingly, the entanglement entropy of the tree state jΨi satisfies In the case of a logarithmic tree, d ∼ log 2 L, Eq. (B3) implies the estimate with some numerical constant c of the order of 1 that depends on the statistical properties of the tree. This result proves the upper bound for the entanglement stated in Eq. (7) in the main text. From the consideration above, we see that the smallest value of entanglement is to be expected when the Schmidt cut at the middle of the chain cuts the tree into just two subtrees, so that there is only one left and one right cluster (L 1 and R 1 ) in Fig. 20. Let us denote the spins of the left FIG. 20. Tree state generated by the SDRG and its entanglement properties. A Schmidt cut at the middle of the system gives rise to a cut of the tree into a forest and prescribes a view of the chain as a collection of clusters lying to the left (L i ) and to the right (R i ) of the cut. With the full chain in the quantum state described by the tree, each cluster has the projection of the total momentum as the only d.o.f. and right cluster by S L and S R , respectively. The density matrix ρ ≡ ρ L=2 depends on the total spin S and its projection M in the state jΨi. It is of dimension ð2S L þ 1Þ and can be written explicitly in terms of Clebsch-Gordan coefficients C: The conservation of the projection of the angular momentum forces then the density matrix to be diagonal: where A particularly simple case is that of S ¼ 0 (which implies S L ¼ S R and M ¼ 0), where Eq. (B7) reduces to and gives the entanglement entropy For a typical tree state, S L ∝ ffiffiffi ffi L p and the entanglement entropy in agreement with Ref. [40]. While we have no proof of the logarithmic scaling of entanglement for arbitrary values of S, M, S L , and S R in Eq. (B7) (and this scaling certainly does not hold in some specific cases, e.g., S ¼ S L þ S R , M ¼ S), we expect that the lower bound on entanglement and the entanglement entropy satisfies S ent ðL=2Þ < clog 2 2 L þ log 2 n T : ðB13Þ We conclude that the entanglement entropy grows logarithmically with the number of tree states involved and of the order of 2 L=2 of them are required to recover the volume-law scaling of ergodic eigenstates.

APPENDIX C: SEARCHING FOR RESONANCES
In this Appendix, we briefly review our numerical procedure for searching resonances.
Let us consider a tree state jΨ 0 RG i generated by the SDRG. Fixing the tree geometry but allowing the values of the block spins in the nonleaf nodes of the tree to take arbitrary values consistent with the rules of the angular momentum addition provides us with the basis in the Hilbert space. In Sec. III A, such a basis is denoted by jΨ a RG i (with a ¼ 1; …; D S 0 ;L ). We are interested in the matrix elements of the Hamiltonian H 0a between the original RG state jΨ 0 RG i and other members of the basis. Of particular importance for us are the resonant situations when jH 0a j > jH aa − H 00 j.
The dimension of the Hilbert space D S 0 ;L scales exponentially with the length L. Fortunately, most of the matrix elements of the Hamiltonian are, in fact, identically zero due to the SUð2Þ symmetry. For an arbitrary pair of spins i and j, the operator S i · S j acting on jΨ 0 RG i can change only those block spins that lie on the path in the tree connecting the spins i and j. Moreover, the selection rules analogous to the ones in the optics limit the possible change in each block spin S to ΔS ¼ AE1 or 0. In addition, ΔS ¼ 0 is forbidden in the case of S ¼ 0.
Using this selection rule together with the fact that the Hamiltonian is just a linear combination of operators S i · S iþ1 , we are able to count and index all the states jΨ a RG i such that H a0 ≠ 0 (we call them the neighbors of the state jΨ 0 RG i) without actually generating them. We denote by K the number of available neighbors.
We then start the random search of resonances among the neighbors. To pick a random neighbor, we generate a random integer from an interval [1; K] and recompute the corresponding neighbor. We then evaluate the matrix element H 0a and the energy difference H aa − H 00 [70]. If the resonance condition is met, we record the information about the resonant neighbor. The random search runs over some number n samp of neighbors. Depending on the size of the system, n samp can reach values up to 2 × 10 7 . Given the number n res discovered during the random sampling, we estimate the total number of resonant neighbors of the given SDRG state by Note that for system sizes L < 2000 the sampling number n samp actually exceeds the total number of available neighbors K so that the random sampling could be replaced by an exhaustive search. For larger system sizes, an exhaustive search becomes, however, unfeasible. We repeat the procedure outlined above for 5000 disorder realizations, generating for each disorder realization a single SDRG state. In long systems, K res (as well as K) fluctuates strongly from sample to sample. The ln K res develops, however, a well-behaved distribution exemplified in Fig. 21. Therefore, in the main text we characterize the proliferation of resonances in long systems by the typical value of K res , e hln K res i , whose dependence on the size of the system is discussed in Sec. IV C.

APPENDIX D: STATISTICAL PROPERTIES OF SDRG TREES
In this Appendix, we study the properties of random SDRG trees. These are obtained by randomly picking spins to fuse together, retaining only their spatial arrangement, ignoring the values of the (bare or renormalized) J's. This simplification allows us to get some analytical results.

Distribution of nearest-neighbor graph distances
We want to prove the claim (16) in the main text, where PðlÞ is the distribution of the random variable l i;iþ1 , namely, the graph distance of two neighboring spins in a generic SDRG tree.
To this end, we consider the ensemble of trees constructed by taking a chain of L spins and fusing them all together, two neighbors at a time. After each fusion, the chain effectively shrinks by one site, and the neighbor structure gets updated accordingly. This process is an approximation of the SDRG procedure where we completely neglect the detailed structure of the J couplings.
More precisely, let a tree be described by a sequence of fusions ði 1 ; …; i L−1 Þ, where i k means that we are fusing, at the kth step, the pair ði k ; i k þ 1Þ (with periodic boundary conditions). In order to emulate the SDRG algorithm, we sample the sequence of fusions uniformly randomly among the L! possible (L − 1) permutations of ð1; …; LÞ, which results in a biased distribution on the set of all binary trees, with "taller" trees being less likely. Now take a generic pair ði; i þ 1Þ in a given tree, and suppose that their common block spin descendant is created at the (k þ 1)th step of the tree construction. All fusions taking place after that step are irrelevant for determining l i;iþ1 , whereas each of the k previous ones may contribute either 0 or 1 to such a distance. In fact, the distance contributed by the jth fusion is a Bernoulli random variable with success probability p j ¼ ð 2 L−j Þ, because the distance between i and (i þ 1) increases only if either one of their descendants is picked out of the L − j possible spins at that step. Moreover, the contributions are uncorrelated, since all free indices are sampled with equal probability regardless of the previous history of the tree construction. Therefore, we have where and the (k) superscript serves as a reminder that our random variable is now being conditioned on k.
Let us compute the cumulant generating function for l ðkÞ − 2: the logarithm of which is By defining x ¼ j=L and α ¼ k=L and taking L → ∞ , we have lnhe −sðl ðkÞ −2Þ i ∼ position of index i in the tuple ði 1 ; …; i L−1 Þ. We can then get rid of the k conditioning by averaging over α ∈ ½0; 1. This averaging gives and, by expanding the denominator in a geometric series, we get Inverting the Laplace transform results in Eq. (16).

Size of the block spins
We now set out to determine the average size of the support of a randomly chosen block spin operator for a random SDRG tree state, which amounts to estimating the number of leaves which connect to a node picked uniformly randomly from the set of nonleaf nodes in a generic fusion tree.
To this end, it is convenient to introduce an alternative (but equivalent) construction for our random ensemble. Consider a single node, and start by attaching two children nodes to it, one to the left and one to the right. We can see this step as a "splitting" for the original node. Now pick with equal probability either one of the resulting leaves and perform the same kind of splitting. Iterate the procedure for a total of (L − 1) times, such that the final number of leaves is L. The leaves are spatially ordered by the order relation induced in an obvious way by the distinction of left and right children. In this way, we obtain a binary tree whose geometry is compatible with an SDRG tree. We can call this tree the "fission tree" ensemble.
We are now going to prove by induction that the fission tree and fusion tree ensembles are equivalent [71].
Suppose that the above claim holds after the (k − 1)th splitting, that is to say, for the ensembles of k-leaved fission and fusion trees. Now, when constructing a fusion tree on (k þ 1) leaves, after the first fusion we are left with an effective k-leaved tree. In order to prove the claim, it is then enough to show that the first fusion does not spoil the ensemble equivalence. By definition of the fusion tree ensemble, it is the case that each one of the initial (k þ 1) leaf pairs has the same probability of being fused at the first step, which means that every one of the k effective leaves after the first step has the same likelihood of being the one resulting from the fusion. Therefore, upon reversing the "time direction," we see that if we allow all the k leaves to split with the same probability, both fission and fusion trees on (k þ 1) leaves are sampled with the same distribution, and the inductive step is completed. It also holds trivially that the two ensembles coincide when k ¼ 1, providing the basis of the induction.
In light of this result, it is possible to assign to each node of a tree the step at which it is split. For instance, the root is always labeled by 1, and the maximum label is L − 1 (note that, similarly to the case of the "fusion labeling" in Appendix D 1, this labeling is not uniquely defined). Now fix k ∈ f1; …; L − 1g and consider the node labeled by k. Introduce the variable t to measure the number of fissions occurring after the kth one, t ∈ f0; …; L − k − 1g, and call NðtÞ the total number of leaves which affect the state of the initial node at "time" t. Since every fission can increment N by 1 at time (t þ 1) only if one of the NðtÞ leaves is picked for the fission, we have the stochastic recursion equation where B½p is a Bernoulli variable with success probability p and p ðkÞ ðtÞ ¼ f½NðtÞ=ðk þ t þ 1Þg. The initial condition must be set to N ðkÞ ð0Þ ¼ 2. This equation is hard to treat due to the NðtÞ dependence hidden inside p k ðtÞ, but it is linear and therefore easily solved in the expectation values: where we use B½p ¼ p. By iterating and simplifying the product on the right-hand side and then looking at the final time, we get This result is the average number of ancestors of a node that is split at the kth fission step. In order to answer our initial question-what is the average number of ancestor elementary spins of a random nonleaf node-we simply take the average on all possible values of k. This process yields showing that the block spins have on average unbounded support in space. [