The random quantum comb: from compact localized states to many-body scars

In this work we investigate the effects of configurational disorder on the eigenstates and dynamical properties of a tight-binding model on a quasi-one-dimensional comb lattice, consisting of a backbone decorated with linear offshoots of randomly distributed lengths. We show that all eigenstates are exponentially localized along the backbone of the comb. Moreover, we demonstrate the presence of an extensive number of compact localized states with precisely zero localization length. We provide an analytical understanding of these states and show that they survive in the presence of density-density interactions along the backbone of the system where, for sufficiently low but finite particle densities, they form many-body scar states. Finally, we discuss the implications of these compact localized states on the dynamical properties of systems with configurational disorder, and the corresponding appearance of long-lived transient behaviour in the time evolution of physically relevant product states.


I. INTRODUCTION
The role of disorder in quantum systems with itinerant degrees of freedom has been the subject of much interest over several decades since the pioneering work of Anderson [1] and collaborators [2]. The field has recently received renewed interest with a specific focus on the interplay between disorder and interactions in many-body localization [3][4][5][6][7], particularly with reference to questions about ergodicity breaking in closed quantum systems and the eigenstate thermalisation hypothesis [8][9][10][11].
Disorder can take many forms, and in some contexts it can be due to the structure or configuration of the system rather than random potential energy or random interaction terms in the Hamiltonian. For example, an imperfect (quasi-) one-dimensional system may have dendritic offshoots [12], as in biological systems [13], percolation clusters [14], and spin chains [15][16][17][18][19][20]. Another example is that of quasi-particles in dimer, vertex and ice models -spin ice materials being a case in point [21]; in these systems, spin correlations impose local constraints on the quasi-particles' dynamics that can result in quasi-1D structures with a backbone and a distribution of dangling ends or offshoots [22,23].
Here we take a closer look at this phenomenon by investigating a model system consisting of a 1D tight-binding chain with linear offshoots whose lengths are distributed randomly (i.e., a random quantum comb model). Classical variants of comb models have been introduced as a stepping stone to understanding diffusion in percolation clusters and other fractal systems [24]; the offshoots lead to the trapping of particles, which inhibits their motion along the backbone, and can result in subdiffusive behaviour [25][26][27][28][29][30][31][32][33].
In the quantum model, we demonstrate that integrating out the offshoots of random lengths results in a localizing disorder for the quasi-particles moving along the backbone. Remarkably, the configurational nature of the disorder produces resonances between the total energy of the particle and the energy levels of the offshoots. These in turn give rise to an extensive number of states with vanishing localization length along the backbone, namely compact localized (CL) states [34,35]. We provide an analytical understanding of these states and a numerical study of their effect on the dynamics along the backbone. Furthermore, we give conditions for their existence in comb-like structures with generic offshoots. Importantly, these states appear also in specific translationally-invariant configurations of the comb structure, where they form flat bands.
We demonstrate that these CL states survive the addition of non-integrable interactions along the backbone. They correspond to atypical, area law states in an otherwise thermal spectrum, i.e., quantum many-body scars [36][37][38][39][40][41][42][43][44][45] that give rise to weak ergodicity breaking. We argue that these scar states have a strong overlap with several physical states of interest in these systems, and are likely to give rise to long-lived transients in their dynamical properties before thermalization eventually sets in, as is expected asymptotically.
Our results were obtained by studying a model system in order to achieve an analytical understanding as well as allowing for numerical simulations of reasonably large systems and times. However, the key properties that we uncovered ultimately hinge on a simple phenomenon: Integrating out the offshoots gives rise to an energy-dependent effective disorder that leads to resonances between the particle energy and the energy of the available states at each site. We expect this behaviour to feature in a broader class of systems. Indeed, we find that the density of states in our model exhibits sharp features reminiscent of the ones observed in a recent study of quantum spin ice [22], alluding to the possibility that such compact localized states may indeed play a role in the properties of that system.
The manuscript is structured as follows. We first introduce the random quantum comb model in Sec. II. Then, in Sec. III, we show that all eigenstates of the model are exponentially localized along the backbone, and look at the implications of this on the dynamics of the system. We discuss the compact localized states in Sec. IV and finally, in Sec. V, we add density-density interactions to the backbone of the comb and show that the aforementioned CL states form many-body scars. We draw our conclusions and outlook in Sec. VI.

II. MODEL
We consider a single quantum particle hopping on a random, comb-like structure (as shown in Fig. 1) defined by the Hamiltonian whereĤ describes hopping along the one-dimensional backbone andĤ is the Hamiltonian on the offshoots, which take the form of one-dimensional chains of varying lengths. The coupling between the two Hamiltonians is given byĤ B-O = −t 2 x (ĉ † x,0ĉx,1 + H.c.). The index x labels the sites on the backbone, which satisfy periodic boundary conditions, and the indices {i x } label the sites on the offshoots. L and x are the lengths of the backbone and the offshoot on site x, respectively, in units of the lattice spacing. The randomness in the structure derives from the lengths of the offshoots, which are drawn from a probability distribution p( ) with ∈ N. The hopping amplitudes t 1 and t 2 , corresponding to the backbone and the offshoots, respectively, are real and positive. For t 1 = t 2 , Eq. (1) represents the quantum version of the random comb model studied in Refs. [26,27,33].

III. LOCALIZATION PROPERTIES
In this section, we prove analytically that all the eigenstates ofĤ 0 in Eq. (1) are exponentially localized along Localization length ξ loc as a function of energy E and the ratio of the hopping parameters t2/t1 for the power law probability distribution of lengths p( ) ∼ −γ with γ = 2.5, calculated using the transfer matrix technique, for systems of size L ∈ [2 × 10 5 , 2 × 10 6 ] sites. The minima in ξ loc correspond to the discrete energy levels of the offshoots, E ∈ {E cl n }, which lead to a resonance in the magnitude of the effective disorder. The dashed region, 2t2 < |E| < t1 + 1 2 t 2 2 /t1, contains no resonances.
the direction of the backbone by investigating the localization length ξ loc (E) at energy E using transfer matrix techniques. We then study numerically the effects of these localized states on the dynamics of the system along the backbone.

A. Transfer matrix results
Denoting the projection of the wavefunction onto the sites with indices (x, i x ) by ψ x,ix , the discrete form of the Schrodinger equation according to the nearest neighbour tight-binding Hamiltonian in Eq. (1) is given by for sites belonging to the backbone. Further, for sites belonging to the offshoots, we can write down where the length of the offshoot on site x is x . This set of equations, determining the wavefunction on the offshoots, can be solved recursively to give where we defined with ≡ E/2t 2 the energy of the particle in units of half the bandwidth of the offshoots. We further define W (0) = 0 for consistency of notation. Equation (6) relates the wave function on the backbone ψ x,0 to the offshoot ψ x,1 connected to the same site. Substituting Eq. (6) into the Schrodinger equation on the backbone in Eq. (4), we obtain Equation (8) describes an effective one-dimensional Anderson model with a random on-site potential energy term x µ xĉ † x,0ĉx,0 , with µ x = −W ( x ), the magnitude of which depends explicitly on the energy E of the particle. Indeed, this process may be thought of as integrating out the offshoots in order to provide the sites on the backbone with a random self-energy. The explicit energydependence induced by the offshoots can result in curious dynamics along the backbone. In the case of infinite offshoots attached to each site, the backbone dynamics can be described by a fractional time Schrodinger equation [46,47] (see also Appendix C).
Crucially, the function W ( ) has simple poles. As a result, the magnitude of the effective disorder diverges at these resonant energies, leading to a vanishing localization length as a function of energy, and the existence of CL states [48]. Such states are also known to occur in percolation clusters in higher dimensions [49].
In our case, the poles of W ( ) are at the energies of a chain of length with open boundary conditions, which has energy levels E cl n ( ) = −2t 2 cos nπ +1 , where n = 1, . . . , . More generally, in the case of an arbitrary (non-interacting) offshoot with HamiltonianĤ x connected to the backbone at site x, the function W (x) is simply proportional to the diagonal element of the offshoot Green's functionĜ x ( ) = (Ĥ x − ) −1 on the site connecting the offshoot to the backbone (see Appendix A for further details). The Green's function will then exhibit poles at the single-particle energies of that offshoot, i.e., the eigenvalues of the offshoot HamiltonianĤ x .
For E = {E cl n }, we introduce the 2 × 2 transfer matrices which allow the Schrodinger equation on the backbone, Eq. (8), to be rewritten as The localization length at energy E is then given by is the largest Lyapunov exponent of the product of transfer matricesT L = . In all panels, the probability distribution for the lengths of the offshoots is indicates an average over the distribution of lengths of the offshoots. Since det T n = 1, the product of transfer matricesT L also has unit determinant and hence its eigenvalues are reciprocals of one another. For a nontrivial [50] probability distribution of the offshoots, p( ), it is possible to use Furstenberg's theorem [51] to show that the Lyapunov exponent is strictly positive and thus the localization length is finite ξ loc (E) < ∞. Figure 2 shows a colour plot of the localization length ξ loc as a function of the particle energy E and the ratio t 2 /t 1 , where ξ loc is computed numerically using standard transfer matrix techniques [52,53]. We checked numerically that the behaviour of ξ loc (E) changes only quantitatively, but not qualitatively, if we use a different probability distribution p( ). As already discussed, for energies satisfying the resonance condition E ∈ {E cl n } the localization length drops to zero, as one can see from Fig. 2. It is interesting to note that ξ loc (E) is largest at the edges of the energy spectrum, while in a standard Anderson model the localization length is instead largest for energies belonging to the middle of the spectrum [52]. This behaviour may be understood by noting that the effective disordered potential W ( x ) is smallest at the edges of the spectrum, implying a larger localization length. For example, in the simplest case of a binomial distribution of lengths, p( ) = 1 2 [δ ,0 + δ ,1 ], the fluctua- . Further, the offshoot Green's function decays with energy outside of the band −2t 2 < E < 2t 2 , leading to a suppression of the disorder: The localization length therefore becomes exponentially large in k within this region.

B. Dynamical properties
We focus on the dynamics of the probability density Π(x, t) starting from a wave packet localized on a single site of the chain, |ψ(0) =ĉ † 0,0 |0 . In particular, we study the probability density marginalised over the offshoot indices that is, the probability of finding the particle on any site of the comb with backbone index x. To quantify the spread of Π(x, t), we define its return probability by R(t) = Π(x = 0, t) and its mean-square displacement As expected, and in agreement with the analysis of ξ loc (E), both the return probability R(t) and the meansquare displacement X 2 (t) exhibit behaviour typical of a localized system. Figures 3 (a)-(b) show R(t) and X 2 (t) for several system sizes (i.e., lengths of the backbone L), where the dynamics has been computed using exactdiagonalization. The lengths of the offshoot are power law distributed according to p( ) ∝ 1/ 2.5 [54], and t 2 /t 1 = φ FIG. 5. Translationally-invariant comb structure that hosts flat bands at E = ±t2 corresponding to the CL states discussed in Sec. IV. The support of the CL states is depicted by the blue sites. t1/t2 = 1 is chosen for the plot of the band structure, although the flat bands persist irrespective of the value of t1/t2.
Both quantities, after some initial transient dynamics, saturate to an L-independent value, implying that the particle cannot propagate through the system beyond the maximal localization length. Finally, Fig. 3 (c) shows the density profile averaged over both disorder realisations and time: Interestingly, for t 1 = t 2 , where the maximum localization length becomes large, ξ loc 1 (see Appendix B), we found a transient dynamics that is consistent with an algebraic propagation X 2 (t) ∼ t α , where α ≈ 2 − 1/γ, with γ > 2 the decay rate of the power law probability distribution of the offshoots p( ) ∼ 1/ γ .

IV. COMPACT LOCALIZED STATES
Having shown that all the eigenstates ofĤ 0 are exponentially localized along the backbone, we now turn our focus to the eigenstates ofĤ 0 with ξ loc (E) = 0, referred to previously as 'compact localized' (CL) states. As we already discussed, these states may be found at the zeros of the function W −1 ( ), which may be thought of as the inverse of the disordered onsite potential.
We are able to construct families of exact eigenstates of the full Hamiltonian, Eq. (1), for all values of t 1 /t 2 at energies E ∈ {E cl n } by considering symmetric clusters containing few sites. For example, has zero mode where the site labels correspond to their position in the state vector, and the blue sites correspond to the sites on which the wavefunction has a nonzero projection. These states are in fact exact eigenstates of the full Hamiltonian in Eq. (1) by virtue of having precisely zero projection onto the sites that connect this cluster to the remainder of the lattice. Further, placing any chain of even length on the intervening site will give rise to an eigenvector of the form (t 1 /t 2 , 0, −1, 0, t 1 /t 2 ) ⊕ (0, 1, 0, −1, . . .) with zero energy.
Similarly, we are able to find clusters of sites that give rise to ±t 2 eigenvalue pairs. For example, 3 4 5 1 2 6 7 6 has two eigenstates with E = ±t 2 Once again, similar eigenstates with energy ±t 2 may be written down for any intervening chain of length = 3k−2 with k ∈ N.
In the case of generic offshoots connected to the backbone, the existence of such CL states is not guaranteed. Indeed, using the above construction, if the two ends of the cluster at sites x − 1 and x + 1 have an identical offshoot with an eigenstate at energy E cl n , then the intervening offshoot at site x must satisfy E cl n + W E cl n ( x ) = 0 for there to exist a CL state with energy E cl n (i.e., an eigenstate satisfying ψ x−1,0 = ψ x+1,0 = 0, which does not connect to the remainder of the chain). This general condition, twinned with the exact expression for W ( ) in Eq. (7), can then be used to deduce the rule for the length of the intervening chain required to give rise to a CL state at energy E cl n : with k ⊂ N such that x ∈ N.
It is important to note that these CL states are distributed with a finite density throughout the bandwidth of the offshoots |E| < 2t 2 , and the corresponding eigenvalues, E cl n ( ) = −2t 2 cos nπ +1 , are extensively degenerate. It is easy to estimate the average number of such states at E ∈ {E cl n } by counting the expected number of occurrences of the above structures. For example, for a given probability distribution p( ), the average number of zero-modes (E = 0) is simply given by N 0 ∼ Lp(1) 2 ∞ k=1 p(2k + 1). A consequence of these macroscopically degenerate energies is that the configuration-averaged density of states is the dimension of the Hilbert space, will exhibit non-analytic points at E ∈ {E cl n }, as shown in Fig. 4. These CL states can also be found in special translationally invariant comb structures, where they form flat bands. Figure 5 shows a comb lattice that hosts 2L/3 such CL states at energies E = ±t 2 [see Eq. (14)].

V. MANY-BODY SCARS
We now tackle the question of adding interactions on the quantum comb model in Eq. (1). We focus specifically on the fate of the CL states introduced in the previous section once density-density interactions are added between adjacent sites on the backbone. Specifically, we consider the following deformation of the free HamiltonianĤ 0 defined in Eq. (1) whereĤ int = xn x,0nx+1,0 , withn x,0 =ĉ † x,0ĉx,0 , corresponds to density-density interactions [55] on the backbone of magnitude V [56]. It is easy to see that due to the spatial structure of the CL states, a Slater determinant of CL states belongs to the kernel of the interaction operatorĤ int . Specifically, let us definê x,ix (E cl n )ĉ † x,ix as the creation operator for the single-particle CL state at energy E cl n , where the index s labels its macroscopic degeneracy. Now, noninteracting eigenstates of the form |ψ cl = n,sη † s (E cl n )|0 remain eigenstates ofĤ since they satisfyĤ int |ψ cl = 0. One can only construct such states as long as the total number of particles does not exceed the total number of CL states. Importantly, these eigenstates are highly nonthermal and violate the eigenstate thermalization hypothesis [11], since they are exact eigenstates of a quadratic Hamiltonian. Moreover, as a result of their strictly localized nature, they satisfy exact area law scaling of the entanglement entropy.
It is important to point out that the existence of these states does not depend on the value of the interaction strength V or the hopping parameters t 1 and t 2 . Thus, in general, the integrability of the model is broken (see Appendix D), and hence these CL states are surrounded by thermal (volume law entangled) eigenstates. These special states therefore constitute an example of exact many-body scar states [37]. Indeed, their construction is reminiscent of Refs. [58][59][60].
As we already discussed for the non-interacting problem, these special states, located at some of the noninteracting energies, appear with probability one in random comb structures, as well as in special translationally invariant models. For example, in the comb structure in Fig. 5 with N = L/3 number of particles, we will have 2L 3 L 3 ∼ 2 2L/3 / πL/3 many-body scar states, which appear at energies ±nt 2 , with n ∈ Z. Since these many-body scars are Slater determinants of CL states, it is easy to find physical (product) states in the computational basis, x,ixĉ † x,ix |0 , that have a large overlap with them. As a result, the dynamics starting from these special initial conditions will be strongly dominated by the existence of the many-body scars. Importantly, such initial conditions are relevant both theoretically as well as experimentally, e.g., in cold atom setups.
For concreteness, we focus onĤ defined on the translationally invariant comb structure shown in Fig. 5. In this case, the charge density wave state |ψ Néel = L/3 x=1ĉ † 3x−2,0 |0 maximises the overlap with the sub-Hilbert space spanned by the scar states. Figure 6 (a) shows the overlap | ψ Néel |E | 2 , between |ψ Néel and energy-resolved eigenstates ofĤ for L = 12, N = L/3, t 2 /t 1 = 2 and V = 1. This shows that |ψ Néel is pre- dominantly supported by the scar states, located at energies E = 0, ±2t 2 and ±4t 2 . It is important to note that |ψ Néel belongs to the middle of the spectrum of H. Conversely, starting from a random product state, i.e., infinite-temperature in the computational basis, the overlap | ψ random |E | 2 is homogeneously spread over the entire energy spectrum ofĤ, as shown in Fig. 6 (b).
The presence of the scar states can be probed dynamically by studying the return probability R(t) = | ψ Néel |e −itĤ |ψ Néel | 2 . We compute the time evolution of R(t) using Chebyshev polynomial techniques [61], which is shown in Fig. 6 (c). The return probability R(t) exhibits coherent oscillations, implying that the time-evolved state of the system is confined within the sub-space spanned by the scar states, whose energies are commensurate. Nevertheless, for asymptotically large times in the thermodynamic limit (L → ∞), we expect the breakdown of this oscillatory behaviour and a decay of R(t) to its equipartition value in which |ψ Néel (t) is spread over the entire available Hilbert space, implying thermalization [62].

VI. CONCLUSIONS
In this work we studied tight-binding Hamiltonians on comb-like structures as a model system to investigate the effects of configurational disorder on localization and transport properties. The model is composed of a one-dimensional backbone decorated with offshoots attached to each site of the backbone. These models are experimentally accessible in quantum synthetic platforms, such as cold atoms and other artificial systems. Moreover, we argue that it may be relevant to the motion of quasi-particles in dimer and vertex models, quantum spin liquids and fractonic systems, supported by recent results on quantum spin ice [21,22].
We focused primarily on the case in which the offshoots assume the form of one-dimensional chains whose lengths are randomly distributed. Pictorially, this model represents a quasi-one-dimensional system in which a particle is able to escape from the main chain (backbone) due to the presence of the offshoots (see Fig. 1).
We considered first the non-interacting limit of the model. We showed analytically and numerically that all the eigenstates are exponentially localized along the direction of the backbone for any amount of disorder in the lengths of the offshoots. Using transfer matrix techniques, we mapped the problem onto a one-dimensional Anderson model with an effective, energy-dependent, onsite disorder. Analysing the behaviour of this effective onsite disorder, we identified special energies for which the onsite disorder diverges. As a result, at these energies, which coincide with the energy levels of the Hamiltonian of the offshoots, the eigenstates are compact localized (CL), characterised by a vanishing localization length. Moreover, the energy degeneracy of these states is extensive, leading to non-analytic points in the density of states. These CL states are also present in translationally-invariant comb structures, where they form flat bands.
Finally, we considered the interacting case, where the interactions act between adjacent sites on the backbone only. Analytically, we proved that Slater-determinants of CL states are exact eigenstates of the interacting model, as long as the particle number is less than the number of such single-particle states. These states, which are highly non-thermal, have area law entanglement scaling while belonging to the bulk of the energy spectrum, providing an example of exactly solvable many-body scars.
Importantly, these states have a large overlap with experimentally relevant states. Therefore, we argue that for a range of physical initialisations of the system, the compact localized states alter dramatically the real-time evolution for experimentally accessible time scales. As a result, we were able to provide a quench protocol to probe dynamically the existence of scar states in the system. Numerically, we tested it by investigating the dynamics of the interacting model in a translationally-invariant comb structure, which hosts perfect quantum many-body scars.
Interesting directions for future work include investigating the possibility of a many-body localization [6] tran-sition in the interacting random comb model, or testing the robustness of the many-body scars uncovered in this work with respect to more generic types of interactions.
Further, comb-like structures can be used to study the dynamics of one-dimensional open quantum systems, where the offshoots play the role of a reservoir/sink for the particles. This connection is related to the paradigm of fractional quantum mechanics [63,64]. For example, the quantum dynamics on the backbone of the non-interacting tight-binding model with infinite offshoots is described by a fractional time Schrodinger equation that exhibits nonunitary time evolution [46,47]. Whether the fractional time Schrodinger equation could be used to investigate interacting, open many-body systems has not studied extensively so far and it is an interesting question that deserves further attention in the future.
Whereas we considered here a model system in order to facilitate in-depth analytical and numerical understanding, the underlying mechanism at play is generic: Integrating out the offshoots produces an energy-dependent effective disorder that leads to resonances between the particle energy and the energy of the available states at each site. For example, the density of states in our model exhibits sharp features reminiscent of the ones observed in quantum spin ice [22], suggesting that compact localized states may play a significant role in the dynamics of those systems and materials.

ACKNOWLEDGMENTS
We are deeply indebted to Tom Gray for several insightful discussions throughout the course of this work from its inception. We would also like to thank Pasquale Calabrese and Sarang Gopalakrishnan for useful discussions. This work was supported in part by the Engineering and Physical Sciences Research Council (EPSRC) Grants No. EP/P034616/1 and No. EP/M007065/1. Some of the numerical simulations were performed using resources provided by the Cambridge Service for Data Driven Discovery (CSD3) operated by the University of Cambridge Research Computing Service (www.csd3.cam.ac.uk), provided by Dell EMC and Intel using Tier-2 funding from the Engineering and Physical Sciences Research Council (capital grant EP/P020259/1), and DiRAC funding from the Science and Technology Facilities Council (www.dirac.ac.uk).

Appendix A: Generic offshoots
In this appendix, we consider the case of generic offshoots emanating from-and possibly connectingvarious sites on the backbone. Let us begin by considering an offshoot which, when isolated from the backbone, is described by a generic HamiltonianĤ x , quadratic in fermionic operators. If the offshoot is then connected to the backbone via a nearest-neighbour hopping term with amplitude −t 2 , the relationship between the projection An offshoot that connects sites x and y on the backbone, leading to disorder in both the potential on sites x and y, and a hopping term between the two sites.
of the wavefunction onto |x, 0 and |x, 1 is given by where we have defined the Green's function of the offshoot G x (ω) = (Ĥ x − ω) −1 . Hence, the value of the effective disordered potential is proportional to the diagonal element of the Green's function on the site that connects the offshoot to the backbone. If instead the offshoots connect to the backbone in multiple locations, then the offshoots not only induce effective on-site disorder, but also induce an effective disordered hopping term between the sites that are connected by the offshoot. For the case shown in Fig. 7 (b), where there exists an offshoot with (isolated) HamiltonianĤ xy that connects to the backbone at sites x and y, then where G αβ (ω) ≡ α, 1|Ĝ xy |β, 1 . The generalisation to offshoots that connect more than two sites on the backbone is simple: The matrix that connects the projection of the wavefunction onto the backbone and the connected offshoot sites is given by the appropriate sub-matrix of the offshoot Green's function.

Compact localized states
If there exist two identical offshoots with Hamiltonian H x±1 =Ĥ (1) on sites x ± 1, and an offshoot on the central site with HamiltonianĤ x =Ĥ (0) , we look for the conditions imposed on the HamiltoniansĤ (i) for there to exist a compact localized (CL) state on this cluster. Again H (i) correspond to generic noninteracting Hamiltonians, quadratic in fermionic operators. The CL state (if it exists) corresponds to an eigenstate of the Hamiltonian H (1) with energy ω n , which requires that ψ x±1,0 = 0 (also required for the state to be an exact eigenstate of the full chain). We then require that the remaining sites satisfy ψ x+1,0 = ψ x−1,0 = −(t 1 /t 2 )ψ x,0 . The final requirement places restrictions on the Hamiltonian of the central offshoot. We find that ω n ψ x,0 = −t 2 ψ x,1 and, from (A1), ψ x,1 = t 2 G (0) (ω n )ψ x,0 , where the Green's function G (0) corresponds to the HamiltonianĤ (0) . We therefore require the consistency condition to be satisfied in order for the CL state at ω n to exist. This equation determines which pairings of offshoots allow for the existence of CL states at the eigenvalues ofĤ (1) . For the case of linear offshoots, the appropriate diagonal element of the Green's function is given by Eq. (7) in the main text. If the offshoots on sites x ± 1 have length x+1 , and the central offshoot has length x , then (7) evaluates to for the energy E n = −2t 2 cos k n . Solving the consistency relation (A3) gives rise to the length constraint on the central offshoot stated in the main text: x = ( x+1 + 1)k/n − 2, for k ∈ N.
If the offshoots are instead given, for example, by Bethe branches with branching ratio z − 1 (i.e., coordination number z), then the Green's function defined by (7) is re- where is now interpreted as the depth of the tree (i.e., the maximum recursion depth). If we repeat the above analysis to find the required depth of the interstitial offshoot, we find in general that the depth required to satisfy (A3) would be non-integer, x / ∈ N. The one exception is for zero modes, E n = 0, in which case the consistency condition can be trivially satisfied [G(0) = 0] by having no offshoot on the intervening site.
These restrictions can be relaxed somewhat if we allow each site on the backbone to be connected to two (or more) offshoots. For example, if two offshoots with identical HamiltoniansĤ x are connected to the backbone at site x, then we can construct a CL state satisfying ψ x,0 = 0 for each eigenstate ofĤ x , which therefore does not connect to the rest of the backbone. The CL state is formed by antisymmetrising an eigenstate |ϕ ofĤ x so that the eigenstate of the full Hamiltonian is of the form |ϕ ⊕ 0 ⊕ |−ϕ . in terms of a fractional time Schrodinger equation of the form given in Ref. [46], i.e., of the form i α ∂ α t ψ = Hψ, where α = 1/2. In the lattice model, and in the thermodynamic limit (in the offshoot direction, N → ∞), one may convert the sum in Eq. (C3) to an integral and hence arrive at the expression [for Re(s) > 0] x h x,x Ψ x ,0 (s) + iψ x,0 (0) s 2 + 4t 2 2 = iΨ x,0 (s) , which has the formal solution ψ x,0 (t) = L −1 −iG x,x i s 2 + 4t 2 2 ψ x ,0 (0) , (C5) where the summation over x is implicit, and the Green's function on the backbone is defined as G x,x ( ) = [(ĥ − E) −1 ] x,x .

Bandwidth
In order to bound the bandwidth, we compute the spectrum of the translationally-invariant model with infinite offshoots. We will begin with offshoots of uniform length N , and take the N → ∞ limit at the end of the calculation. Taking the Fourier transform over the backbone direction, we arrive at the Hamiltonian H 0 = −2t 1 kx cos k xâ † kx,0âkx,0 − t 2 kx, i,j â † kx,iâkx,j .
(C6) In particular, to find the bandwidth, we seek the extremal eigenvalues ofĤ 0 above. Parametrising the energy as E = −2t 2 cosh η, the quantisation condition for offshoots of length N is found to be 2t 1 cos k x − t 2 e η 2t 1 cos k x − t 2 e −η = e −2N η . (C7) In the thermodynamic limit N → ∞, the solution of this equation for real η is given by 2t 1 cos k x = t 2 e η . This solution remains finite for t 2 → 0 + , whilst the other N − 1 eigenvalues vanish. Therefore, in the thermodynamic limit, we find the extremal eigenvalue for t 2 < 2t 1 : This result defines the white dashed region in Fig. 2.

Appendix D: Level statistics
To provide evidence that the system is non-integrable, we analyse the level statistics of the interacting Hamiltonian (16). In particular, we study the distribution of (unfolded) level spacings P (s), where s n = E n+1 − E n [65], and the level spacing statistics value r n = min(s n , s n−1 )/ max(s n , s n−1 ) [66]. It is crucial to resolve all symmetries of the system and calculate the statistics within each symmetry sector separately. For the system shown in Fig. 5, the Hamiltonian possesses both translational invariance and inversion symmetry. In Fig. 8, we plot the level statistics for a system of size L = 12, with N = L/3 particles averaged over momentum sectors (excluding k = 0, π). We find that the distribution P (s) is in good agreement with the Gaussian Orthogonal Ensemble (GOE) from random matrix theory, as expected for non-integrable models [11]. Further, the rvalue is r = 0.534, to be compared with the value of the GOE, r GOE = 0.5359, and of uncorrelated energy levels r Poisson = 2 log 2 − 1 ≈ 0.3863 for integrable models.