## Abstract

Many-body localized (MBL) systems lie outside the framework of statistical mechanics, as they fail to equilibrate under their own quantum dynamics. Even basic features of MBL systems, such as their stability to thermal inclusions and the nature of the dynamical transition to thermalizing behaviour, remain poorly understood. We study a simple central spin model to address these questions: a two-level system interacting with strength *J* with *N*≫1 localized bits subject to random fields. On increasing *J*, the system transitions from an MBL to a delocalized phase on the *vanishing* scale *J*_{c}(*N*)∼1/*N*, up to logarithmic corrections. In the transition region, the single-site eigenstate entanglement entropies exhibit bimodal distributions, so that localized bits are either ‘on’ (strongly entangled) or ‘off’ (weakly entangled) in eigenstates. The clusters of ‘on’ bits vary significantly between eigenstates of the *same* sample, which provides evidence for a heterogeneous discontinuous transition out of the localized phase in single-site observables. We obtain these results by perturbative mapping to bond percolation on the hypercube at small *J* and by numerical exact diagonalization of the full many-body system. Our results support the arguments that the MBL phase is unstable in systems with short-range interactions and quenched randomness in dimensions *d* that are high but finite.

This article is part of the themed issue ‘Breakdown of ergodicity in quantum systems: from solids to synthetic matter’.

## 1. Introduction

There has been considerable recent theoretical and experimental activity in discovering and understanding the dynamical phases of generic isolated quantum systems. The many-body localized (MBL) phase has been of particular interest, as it is a dynamical phase in which an isolated many-body quantum system is not able to act as a reservoir and bring its subsystems to thermal equilibrium [1–7]. MBL phases exhibit a complex of intriguing properties absent in conventional thermalizing phases, including an extensive set of localized conserved quantities (l-bits) [8–14], slow entanglement dynamics [15–18], and dynamically protected long-range orders in highly excited states that are forbidden in equilibrium [19–23] among other interesting features [24–30]. MBL phases (or at least very slow dynamics approximating MBL) have been experimentally observed in cold atoms [31–34], in trapped ions [35] and in diamond with nitrogen-vacancy centres [36], in various regimes (one and two dimensions, with long-range interactions, with quasi-periodic potentials, etc.) and there are many more exciting experiments in the pipeline.

While MBL phases have been essentially proven to exist in one-dimensional systems [10], their existence and stability in higher dimensions remains controversial. Chandran *et al.* [37] studied the effect of a thermalizing boundary on a *d*>1 dimensional MBL phase and argued that the l-bits and other eigenstate measures of localization are unstable; instead, approximately conserved ‘l*-bits’ underlie what remains of localization in such systems. Another study of the effect of thermal inclusions was presented in [38], in which the authors argue that any non-zero density of large enough thermalizing inclusions will destroy the MBL phase in *d*>1, converting it to an extremely slow, but ultimately thermalizing, regime. The outline of their argument is as follows. Consider a single thermal inclusion—for example, a rare finite region where the disorder is atypically small. If the inclusion is sufficiently coupled to the MBL phase in which it is embedded, then it will thermalize the localized degrees of freedom (l-bits) bordering the inclusion, and this allows the thermal inclusion to grow in size. This then makes the inclusion a more effective ‘bath’ for l-bits that are further away.

In *d*=1, there is a critical value for the exponential decay length of the interaction between the inclusion and the l-bits beyond which this process runs away. This was recently demonstrated numerically in [39]. In *d*>1 systems with randomness, this process *always* runs away so that there is a finite, but possibly extremely large, time scale over which a putative MBL phase is thermalized by the non-zero density of rare thermal inclusions which arise for any strength of disorder [38]. This argument leaves open the possibility that MBL phases can remain stable in higher *d* in non-random systems with quasi-periodic fields or potentials that do not produce such rare thermalizing regions.

One of the motivations of this article is to study the key step in the destabilization process described above in *d*>1; that is, how a single inclusion interacting with the same strength *J* with many l-bits on its border leads to their thermalization. We show that a single two-level inclusion is sufficient to thermalize *N* neighbouring l-bits in the large-*N* limit at any non-zero *J*/*W*, where *W* is the strength of the random fields on the l-bits. In order for such a ‘central spin’ system to remain in the MBL phase in this large-*N* limit, the interaction with the inclusion must be scaled down with increasing *N*. Specifically, we find that the critical *J*_{c}(*N*) separating the MBL phase for *J*<*J*_{c}(*N*) from the thermalizing phase at *J*>*J*_{c}(*N*) scales as *J*_{c}(*N*)∼*W*/*N* up to multiplicative corrections (figure 1).

At the transition for our central spin model, the interactions are crucial to the dynamics of the system but do not contribute to the system’s equilibrium thermodynamics, as the coupling *J* would need to be much larger (of order *N*^{0}) to do so. Like other fully connected quantum spin glass models, the dynamical phase diagram simply does not coincide with the thermodynamic phase diagram [40–44]. These models have been extensively studied in the context of quantum dots, ion traps and other long-range interacting spin systems; they provide a tractable analytic setting for studying many-body quantum dynamics. In all of these models, there is a polynomially small in *N* coupling strength *J*_{c}(*N*) separating an MBL phase from the phases that fully or partially thermalize. At finite temperature, these models can also have glass phases where the system fails to fully thermalize due to divergent free energy barriers in the limit. It is important to make the distinction between arrested dynamics due to (i) detuning, which is the cause of Anderson and many-body localization,^{1} and due to (ii) divergent free energy barriers. In the latter case, the system can remain delocalized and function as a thermal bath for itself, only failing to thermalize a small number of collective degrees of freedom [43]. Finally, for these fully connected models there are no destabilizing rare region effects, allowing for a stable MBL phase at *J*<*J*_{c}(*N*).

In the sequel, we analyse the infinite-temperature dynamical phase diagram of our central spin model using exact diagonalization with up to *N*=13 l-bits, and using small *J* perturbation theory at large *N*. Numerically, the single-site eigenstate entanglement entropy, energy level repulsion and many-body eigenstate participation ratios all support the existence of an MBL phase for *J*<*J*_{c}(*N*)∼*W*/*N* and a thermalizing (ETH)^{2} phase for *J*>*J*_{c}(*N*). They also reveal many interesting features about the localized phase and the crossover region at *J*≈*J*_{c}(*N*). In the localized phase, the mean single-site eigenstate entanglement entropy [*S*_{1}] decreases as 1/*N*, while the participation ratio distributions are *N*-independent.

In the crossover region, the l-bits are either ‘on’ (strongly entangled) or ‘off’ (weakly entangled) in eigenstates and the pattern of ‘on’ l-bits varies significantly between states of the same sample (and, of course, between samples). Thus, single-site observables are very heterogeneous in real space, in energy space and across disorder realizations in the crossover region, suggesting that they change discontinuously as , in line with recent proposals [50,51] that few-body observables are similarly discontinuous across the MBL transition in one dimension.

The perturbative analysis at small *J* on the classical hypercube explains many of these numerical observations. At second order, a typical initial configuration is resonant with *K*∼*J*^{2}*N*^{2}/*W*^{2} states in which two l-bits are flipped and the central spin is flipped. These states are in turn resonant with ∼*K* *completely new* states at the next order. The resulting ‘resonant subgraph’ is therefore locally tree-like and we argue that the statistical properties of the resulting eigenstates in the localized phase can be understood via an associated bond percolation problem on the hypercube.

This mapping, however, does not capture the transition region, in particular, the heterogeneity in single-site observables. This is not particularly surprising, as we have neglected higher-order processes. Altshuler *et al.* [52] treat these processes within the tree approximation in a related model. Their arguments seem to predict an intermediate delocalized non-ergodic phase between and *J**(*N*)∼*W*/*N*, that is, between the MBL phase and the thermal phase in the central spin model. The numerical data at accessible system sizes are not conclusive about the existence of this possible intermediate phase or the logarithmic suppression of the critical coupling. We speculate on the possible phase diagrams in the thermodynamic limit in §4f.

One difference between our central spin model in the limit of large *N* and a short-range model in the limit of large *d* is due to the order of how these two limits are taken. Suppose we are in large but finite dimension *d*, with a system of size *N* spins and with coordination number *z*. If we keep the coordination number *z* fixed and take the thermodynamic limit of , then our sample is infinite and it has an infinite number of chances to produce a rare thermalizing inclusion, so for any non-zero spin-flip interaction *J* these inclusions are present at some non-zero density and they destabilize the MBL phase with probability one. In our central spin model, on the other hand, we take *z* to infinity along with *N*, and scale *J* accordingly, which yields a stable MBL phase. The reason this MBL phase of the central spin model can be stable is that for any finite *N* each sample is finite, so the probability that it has a destabilizing thermal inclusion (a rare set of mutually resonant l-bits) is less than one. And this probability goes to zero in the limit when the scaled *J* is within the MBL phase.

In what follows, we first describe the model in §2 and present the numerical diagonalization results in §3. We then turn to the perturbative analysis at low orders in §4b,c and compare their quantitative predictions for the MBL phase and the crossover region to the numerical results in §4d. Finally, we discuss the role of higher-order processes in §4e.

## 2. Model

Our Hamiltonian for a thermal inclusion connected to *N* l-bits is
2.1where are independently sampled random variables drawn from the box distribution [−*W*,*W*], for *α*=*x*,*y*,*z* are the Pauli operators of l-bit *i*, and *A*_{i} and *B*_{i} are real Gaussian symmetric random matrices (Gaussian orthogonal ensemble, GOE) acting on the thermal inclusion Hilbert space of dimension 2^{M}. We caution the reader that, while the *τ*_{i} operators act on l-bit *i*, the *A*_{i} and *B*_{i} operators *act* on the central spin; the subscript *i* indicates which l-bit they are coupled to in *H*. The normalization of the GOE matrices is such that the off-diagonal elements have variance 1/2^{M}, while diagonal elements have twice the variance. This normalization guarantees that the operator norm of *A*_{i},*B*_{i} is order 1, which is the appropriate scaling for a local operator acting on the inclusion.

In this paper, we focus on the simplest case of *a central spin model* with *M*=1, as it captures most of the physics of the general *M*>1 model in equation (2.1). At *M*=1, it is convenient to expand the *A*_{i} matrices in the Pauli basis *σ* of the central spin:
2.2As *A*_{i} is a GOE matrix, , while for *α*=0,*x*,*z* are independent Gaussian random variables with zero mean and variance 1/2. The Hamiltonian for the central spin model is therefore
2.3where is a renormalized field.

There are two dimensionless scales in the model: the temperature and the coupling constant *J* in units of *W*. As we focus on infinite temperature in this article, *J*/*W* determines the entire dynamical phase diagram. At *J*=0, none of the spins are coupled and the eigenstates are trivially localized product states with . Each *τ*^{z} state of the l-bits is doubly degenerate as the energy is independent of the state of the central spin *σ*. In anticipation of the leading interaction term at non-zero *J*, we choose to work in a basis in which *σ* points along or against the effective -dependent field:
2.4Thus, for *s*=±1 is an eigenstate of equation (2.3) with *J*=0, and, in fact, even for *J*≠0 for *B*_{i}=0. This set of ‘Fock states’ are naturally viewed as the vertices of an (*N*+1)-dimensional hypercube. We use the shorthand and drop the explicit dependence of ** h** on

*τ*where it is clear from context.

The statistics of the effective field ** h** play an important role in the dynamics of the model. As the sum of

*N*independent Gaussian random vectors of

*O*(

*J*) in the

*xz*-plane,

**has mean zero and typical (root-mean square) length . Its distribution is Gaussian with respect to disorder fluctuations at fixed**

*h**τ*

^{z}(by construction) and with respect to varying in a fixed sample (by the central limiting theorem).

In the opposite limit of *J*≫*W*, the bare random fields become negligible. Each l-bit feels a random field of order *J* in the *xz*-plane and interacts with the central spin with strength order *J*. As the disorder and the interactions are of comparable strength, we expect the system to thermalize as . Both the perturbative (in *B*) arguments and numerical exact diagonalization study discussed below confirm this expectation. However, as *W* drops out for large *J*, we cannot simply increase *J* in order to approach a more robust thermalizing limit at finite size. This unfortunately results in large finite-size effects in the thermalizing regime.

## 3. Numerical results

We study equation (2.3) using full diagonalization for an *M*=1 central spin coupled to *N*=7 to *N*=13 l-bits. The number of samples at each (*N*,*J*) is 2500, except for *N*=13 where this number is 600, and within each sample, we restrict our analysis to the eigenstates within the central half of the energy spectrum in order to study the properties of infinite temperature. The mean of a quantity *q* is denoted by [*q*], while the standard deviation is denoted by Δ*q*. Unless specified otherwise, the mean and standard deviation are taken with respect to all the l-bits (*l*), the eigenstates in the central half of the spectrum (*E*) and all samples (*s*). When the statistical operation is restricted to a particular subset of these, we include the subset in the subscript. We measure energy in units where *W*=1.

### (a) Three regimes

Numerically, we find three distinct regimes as a function of rescaled coupling *JN*: localized (MBL) at small *JN*, thermal (ETH) at sufficiently large *JN* and a wide crossover at intermediate *JN*. This is conveniently summarized by the behaviour of the mean single l-bit entanglement entropy [*S*_{1}] (figure 2). The entropy, , measures the degree to which l-bit *i* is thermalized within an eigenstate |*E*〉:
3.1Here, is the reduced density matrix for site *i*. If the state |*E*〉 is thermal (at infinite temperature) for l-bit *i*, then obtains its maximal value. In a localized state, on the other hand, can be <1.

The entropy shows the three qualitatively distinct regimes in figure 2. At sufficiently small coupling, we find that [*S*_{1}]*N* collapses onto a single curve to excellent precision. This strongly localized behaviour ([*S*_{1}]∼1/*N*→0 as ) also follows in the perturbative treatment of the localized phase in §4. This is a particularly strong form of localization, seen also in other long-range models [41], as [*S*_{1}]∼*O*(1) in one-dimensional MBL phases [53,54]. We consider this regime the finite-size precursor to the MBL phase. At sufficiently large coupling , [*S*_{1}]*N* develops plateaus which increase with *N*, consistent with thermalizing behaviour in the thermodynamic limit ([*S*_{1}]→1 as ). Finally, the wide crossover between these two behaviours exhibits growth of [*S*_{1}]*N* with *N*, which suggests at least partial delocalization.

In order to focus in on the intermediate transition regime, we look at the finite-size behaviour of [*S*_{1}] (unscaled) along with that of the mean level spacing ratio [*r*] (figure 3).^{3} The ratio measures the level repulsion in a system and is commonly used to diagnose (de)localization. It flows to *r*_{GOE}≈0.53 in systems with random matrix level statistics and to *r*_{Poisson}≈0.39 for systems with Poisson level statistics [4,55]. Figure 3*a* indicates that [*S*_{1}] has a sharpening crossover from 0 to 1 near *JN*∼0.5 while figure 3*b* shows a similar sharpening crossover in [*r*] at *JN*∼1. This suggests that the crossover between the localized and thermal phases sharpens into a phase transition on the scale *JN*. However, as the finite-size effects are large, it is difficult to separate two possible scenarios. The first scenario posits that the two crossover points coalesce as and there is a direct transition from MBL to a fully delocalized thermal phase on the scale *JN* (up to multiplicative logarithmic corrections). In the second scenario, the crossover points remain separate as , sharpening into two phase transitions surrounding an intervening partially delocalized phase [52], where [*S*_{1}]=1, the level statistics are Poisson and the off-diagonal matrix elements of local operators would presumably not satisfy ETH. This is often referred to as a delocalized non-ergodic phase.

We comment that [*S*_{1}] and [*r*] both attain plateaus at large *JN* in the thermal phase. We have checked that the plateau values approach their limiting ETH values as 2^{−N/2}.

### (b) Heterogeneity in observables

The most striking feature of the crossover regime is the extreme heterogeneity of single-site observables across eigenstates, samples and l-bits. As a coarse measure, consider the standard deviation of the single-site entanglement entropy Δ*S*_{1} shown in figure 4. In both the MBL and ETH regimes of figure 4, the ‘flow’ of Δ*S*_{1} is to 0 as *N* grows. This is to be expected in the infinite-temperature ETH phase, where the fluctuations in single-site observables across eigenstates are exponentially small in *N* (Δ*S*_{1}∼2^{−N/2} to be precise). On the MBL side, the perturbative picture of §4 also indicates that the entire distribution of *S*_{1} scales to 0 as 1/*N* so that Δ*S*_{1}∼1/*N*.^{4} The data in figure 4 are consistent with these predictions.

In the crossover regime, however, Δ*S*_{1} shows a peak whose height increases with *N*. As *S*_{1} is a bounded variable, the growth of must saturate at larger sizes. Nonetheless, should the peak persist to the thermodynamic limit, it must converge to a critical point (see the discussion below). Previous works in one-dimensional models have likewise used the peak as a sensitive proxy for the critical point [22,51]. We accordingly define *J*_{p}(*N*) by the location of the peak in Δ*S*_{1} and study the properties of the critical region by following this coupling. The inset to figure 4 confirms that it scales as *J*_{p}(*N*)∼*W*/*N*.

The relative width of the peak at *J*=*J*_{p}(*N*) narrows as *N* increases (figure 5). This trend is consistent with the existence of a sharp phase transition in *S*_{1} at *J*_{p}(*N*) in the thermodynamic limit. Indeed, *J*=*J*_{p}(*N*) approaches the location of the crossing point in [*S*_{1}] in figure 3*a* at *JN*≈0.5 with increasing *N*. As *S*_{1} detects if l-bits are entangled in eigenstates, we expect that it is sensitive to whether the system is localized or delocalized on the classical hypercube of ‘Fock states’, but not necessarily to whether the delocalized phase satisfies the off-diagonal criteria for ETH. Figure 5 is, therefore, consistent with both scenarios for the phase diagram discussed in §3a. In the first scenario, there is a direct phase transition from an MBL to an ETH phase at *J*=*J*_{p}(*N*) as . In the second scenario, in which there are two phase transitions with an intervening partially delocalized phase, *J*=*J*_{p}(*N*) is the position of the first transition out of the localized phase.

We note that the narrowing of this peak on the localized side is slower than on the delocalized side. This makes the extraction of a finite-size scaling exponent problematic and indicates that the numerics are not yet in the asymptotic scaling regime, as we discuss below. Moreover, the narrowing on the delocalized side suggests exponents which violate finite-size scaling bounds [56,57].

We now characterize the heterogeneity in *S*_{1} in the crossover region in more detail. Figure 6*a* shows the wide variation of *S*_{1} across both l-bits and eigenstates within a typical sample. More quantitatively, figure 6*b* shows that the distribution of *S*_{1} across l-bits, eigenstates and samples at *J*=*J*_{p}(*N*) is increasingly bimodal with increasing *N*, with increasing weight near zero and one bit and decreasing weight at intermediate values of *S*_{1}. Thus, l-bits are mostly either ‘on’ (highly entangled) or ‘off’ in any given eigenstate in the crossover regime. The bimodal distribution leads to the large peak in Δ*S*_{1} seen in figure 4.

As l-bits in weak (strong) fields might strongly (weakly) entangle with the central spin and with one another, the reader may not be surprised by the bimodality seen at *J*_{p}. This explanation, however, misses the remarkable feature that *different* l-bits are active in different eigenstates. To demonstrate this effect, we investigate the sample-averaged variation of across eigenstates, , for fixed l-bit *i* (figure 7*a*). If l-bit *i* were entangled/unentangled across all the states, would be small. Instead, we see a robust peak in this quantity at *J*=*J*_{p}(*N*) at each *N*, with a magnitude comparable to the total variation of *S*_{1} across all samples, states and l-bits (Δ*S*_{1}).

In the remainder of figure 7, we parse the contributions to Δ*S*_{1} coming from sample, eigenstate and l-bit fluctuations, all of which contribute significantly to the total variation at these sizes [51]. We find that for our range of *N* the largest contributions to Δ*S*_{1} near the transition come from variations between l-bits in the same eigenstate, i.e. [Δ_{l}*S*_{1}]_{E,s} (figure 7*b*), and across different eigenstates for the same l-bit, i.e. (figure 7*a*). This is in contrast with previously studied one-dimensional models with quenched disorder exhibiting an MBL–thermal transition: in Khemani *et al.* [51], the sample-to-sample variation is the dominant contribution at the largest sizes accessible to diagonalization.

The sample-to-sample variation shown in figure 7*d* is the only component which shows an *accelerating* trend to stronger and sharper peaks as *N* increases. This increasing trend cannot continue with increasing *N* as Δ*S*_{1} is bounded; this indicates that we are not in the asymptotic scaling regime at the numerically accessible system sizes. Nevertheless, should the sample-to-sample variation become the dominant contribution to Δ*S* in the thermodynamic limit, the transition is first-order in the sense that we discuss in §5.

The inhomogeneity between eigenstates in the l-bit-averaged entanglement at the transition is shown in figure 7*c*. The peak shows a slight decrease with *N*. If this persists to larger *N*, then all of the eigenstates in one sample would become similar in this respect, although it seems likely that they will still differ in which l-bits are less versus more entangled.

Finally, we note that this central spin model exhibits the strongest single-site bimodality in the MBL–ETH crossover of any numerically studied model to date, suggesting that the transition may indeed be ‘first-order’ [50]. Furthermore, figure 7*a* shows that fluctuates strongly between eigenstates of a single sample, even at infinite temperature. No renormalization group treatment to date takes these intra-sample variations between eigenstates into account, so might be missing some important physics of this transition [58,59].

### (c) Eigenstate distributions

As discussed in §2, the eigenstates at *JN*=0 are the product states . It is therefore natural to study the support of the eigenstates at *JN*>0 on this basis (which forms a hypercube). This is measured by the participation ratio PR of an eigenstate |*E*_{m}〉:
3.2

Figure 8 shows the evolution of the distribution of the participation ratios with coupling. In the localized phase (illustrated by *JN*=0.2), the distribution strongly decays with PR and is independent of *N*. This implies that the eigenstates are essentially localized to finite regions of the hypercube even in the thermodynamic limit. This is consistent with the perturbative picture of the localized phase in §4.

In the crossover regime (at *JN*=0.8), the distribution remains peaked at small PR, but the weight in the tail grows with *N*. Further, the scale of the cut-off grows exponentially with *N*. In the second-order perturbative treatment developed in §4, the transition coincides with the high-dimensional percolative transition on the hypercube. In this picture, the PR distribution should decay as a power law with exponent 3/2 which is cut off by *n*_{ξ}∼2^{2N/3}. This scaling of the cut-off is in agreement with the numerical data at *JN*=0.8, although the curvature in the bulk of the distribution is clearly inconsistent with a simple power law.

Finally, deep in the thermalizing phase (at *JN*=4), the PR at the peak of the distribution increases exponentially with *N*, consistent with ETH. The probability density at finite PR appears to decrease exponentially with *N*, consistent with exponentially rare samples probing the localized side of the transition.

## 4. Perturbative analysis

The central spin model is amenable to perturbative study in small *B* relative to the ‘classical’ configurations defined in §2. In this approach, delocalization takes place when a typical starting configuration resonates with a divergent subgraph of degenerate configurations on the hypercube.

In §4b,c, we consider first- and second-order processes and show that a delocalization transition takes place at *J*_{c}∼*W*/*N*. At this coupling, only a small (finite) number of spins can resonate by first-order spin-flip processes. On the other hand, second-order ‘flip-flop’ processes, in which two l-bits and the central spin are flipped, produce a divergent subgraph of resonant configurations for *J*>*J*_{c}. Thus, the low-order perturbative analysis captures a delocalization transition on a scale *J*_{c}∼*W*/*N*, consistent with the numerical observations.

As the resonant subgraphs produced at second order are locally tree-like with an expected branching number *K*∼*J*^{2}*N*^{2}/*W*^{2}, their statistical properties can be determined from independent second-neighbour bond percolation on the hypercube with probability of placing a bond, *p*∼*J*^{2}/*W*^{2}. We summarize the relevant properties of this percolation problem in §4d. In §4e, assuming that actual many-body eigenstates are simply delocalized over these resonant subgraphs, we explain the numerical observations that [*S*_{1}]∼1/*N* in the localized phase, as well as the distribution of the participation ratios in the classical basis. Thus, the low-order perturbative analysis also quantitatively describes the localized phase for *J*<*J*_{c}.

However, the low-order perturbative analysis does not capture all of the features of the crossover region at *J*≈*J*_{c}. This is not surprising as we expect higher-order processes to be important in the vicinity of the delocalization transition. In §4f, we adapt the arguments of Altshuler *et al.* [52] that suggest a transition from localized to a delocalized non-ergodic phase on the scale . As discussed in §3, the numerical evidence for this intervening phase is inconclusive.

### (a) Preliminaries

At *B*=0, the state has classical energy
4.1These 2^{(N+1)} configurations are naturally viewed as the corners of a hypercube. The interaction term, , defines a short-range hopping model on the hypercube as it connects the corners related by single l-bit flips either with or without flipping the central spin *s*. As the perturbative arguments will predict a delocalization transition at *J*_{c}∼*W*/*N*≪*W*, we assume *J*≪*W* henceforth, since here we are only considering the large-*N* regime.

To perturb around a typical infinite-temperature configuration , we need the classical energies of nearby configurations. We denote the l-bit configuration obtained by flipping *k* l-bits *i*_{1},…,*i*_{k} by . From equation (4.1), the energy of such a state relative to the initial state is
4.2For small numbers of flipped l-bits (*k*≪*N*), the deviation in the effective field on the central spin,
4.3is much smaller than the field itself, . Thus, expanding equation (4.2) in *J*** δ**/|

**(**

*h**τ*)| (and suppressing the indices of the

*k*flipped l-bits), 4.4The first term is the dominant energy change on flipping the central spin relative to its local field. It is of order but independent of the choice of flipped l-bits. The second term is the total field energy of the flipped l-bits. The bare fields Δ

_{i}∼

*W*while the correction due to the central spin is only ∼

*J*≪

*W*. Finally, the third- and higher-order terms represent the interactions between the flipped l-bits induced by their interaction with the central spin. In particular, the leading interaction is of order .

### (b) First-order processes

At first order in *B*, there are 2*N* neighbouring configurations with a single l-bit flipped with or without flipping the central spin. The matrix element
4.5while the energy differences Δ*E*_{s′τ′j} are distributed on a band of width *W*. As the level spacing is ∼*W*/2*N*, there are ∼2*NJ*/*W* resonant neighbours. Thus, first-order processes begin to find *O*(1) resonances when *J*∼*W*/*N*. We call these ‘step 1’ first-order resonances.

Diagonalizing the perturbation exactly on these resonant configurations produces new basis states delocalized over the resonant cluster. These in turn connect at first order to configurations on the hypercube with two l-bit flips. For *J*∼*W*/*N*, there are again *O*(1) next-neighbour states which satisfy the resonance condition (‘step 2’ resonances). However, as the interaction energy in equation (4.4) is much smaller than the level spacing *W*/*N*, approximately the same l-bit flips that were resonant at step 1 are resonant at step 2.

More precisely, if l-bit flip *i*_{1} is resonant at step 1, Δ*E*_{s′τ′i1}<*O*(*J*), then l-bit *i*_{2} can resonantly flip at step 2 only if Δ*E*_{s′′τ′i1i2}<*O*(*J*). This requires that the effective field on l-bit *i*_{2}, , itself should be close to 0 (if *s*′′=*s*) or |** h**| (if

*s*′′=−

*s*) to an accuracy of

*O*(

*J*). The interaction energy between l-bits

*i*

_{1}and

*i*

_{2}does not change the above resonance condition for l-bit

*i*

_{2}at step 2 because it is negligible on the energy scales

*J*≈

*W*/

*N*. As there are only

*O*(1) l-bits whose field energy is close to 0 or |

**| to accuracy**

*h**O*(

*J*), even multi-step first-order processes only produce a small cluster of resonant configurations (figure 9).

We have neglected the *O*(*J*) shifts in the reference energy *E*_{sτ} which arise due to the diagonalization over the resonant configurations. This is analogous to neglecting the real part of the self-energy in a locator expansion. As these corrections effectively shift the resonance condition rigidly for all l-bits, they do not modify the statistics of the resonant clusters generated.

The structure of the resonant clusters suggests that there will be eigenstate dependence of which l-bits resonantly flip at *J*∼*W*/*N*. There is an *O*(1) subset of the l-bits with fields Δ_{i}∼*O*(*J*) which resonate across all eigenstates without flipping the central spin. Meanwhile, those l-bits with larger fields resonate (in *O*(1) groups) only in those eigenstates whose central spin field matches Δ_{i}.

### (c) Second-order processes

As all but *O*(1) of the 2*N* nearest-neighbour states are off-resonant when *J*∼*W*/*N*, in the following we neglect the first-order resonances entirely. At second order in *B*, there are neighbouring configurations with two l-bits flipped with or without flipping the central spin. In appendix A, we show that the effective matrix elements are largest to the subset of states in which the central spin is flipped from *s* to −*s* relative to the central field. This is physically sensible, as delocalization in this model proceeds due to the interactions with the central spin, so that the divergent resonant subgraph should include classical configurations in which the central spin is flipped. The effective matrix elements to the states are typically ∼*J*^{2}/*W*. As there are ∼*N*^{2} such final states with energies Δ*E*_{−sτ′ij} spread over a bandwidth *W*, the level spacing is ∼*W*/*N*^{2}. This leads to (*J*^{2}/*W*)/(*W*/*N*^{2}) resonances, so that second-order processes begin to find *O*(1) resonances when *J*∼*W*/*N* (‘step 1 second-order’ resonances).

Now consider the ∼*N*^{2} states accessible from a step 1 resonant state with two more l-bits flipped and the central spin returned to *s* (as the central spin flips at each step). In order for the l-bit pair *kl* to be resonant, the state must have energy |Δ*E*_{sτ′ijkl}|<*J*^{2}/*W*. This is on the scale of the level spacing *W*/*N*^{2}. As Δ*E*_{sτ′ijkl} does not contain the energy change on flipping the central spin , and is much greater than the level spacing *W*/*N*^{2}, the pairs *kl* that satisfy the resonance condition at step 2 are typically not the same as the pairs *ij* that satisfy the resonance condition at step 1. Thus, the resonant subgraph includes new l-bit flips at step 2 as compared to step 1.

At step 3, the ∼*N*^{2} accessible states involve three l-bit pair flip-flops with the central spin again flipped to −*s*. The resonance condition |Δ*E*_{−sτ′ijklmn}|<*J*^{2}/*W* picks *different* pairs as compared to step 1 because the interaction energy,
4.6has cross terms between the pairs *ij*, *kl* and *mn* which cannot be ignored. The cross terms are typically ∼ *J*^{2}/|** h**|∼

*W*/

*N*

^{3/2}which is much larger than the level spacing

*W*/

*N*

^{2}. Thus, the pairs

*mn*that satisfy the resonance condition |Δ

*E*

_{−sτ′ijklmn}|<

*J*

^{2}/

*W*are not the same as the pairs

*ij*that satisfy the resonance condition |Δ

*E*

_{−sτ′ij}|<

*J*

^{2}/

*W*at step 1. When

*J*∼

*W*/

*N*, we therefore find

*O*(1)

*completely new*resonant states among the configurations accessible by flip-flops at

*each*step.

From the above construction, we see that the resonant subgraph of the hypercube generated by second-order processes is tree-like for small enough *J* (figure 10). The expected branching number of the tree is *K*∼*J*^{2}*N*^{2}/*W*^{2}. As is well known, such random trees undergo a continuous percolation transition at *K*=*K*_{c}∼1 at which the probability that the resonant subgraph is infinite starts growing from zero.^{5} This provides a beyond which second-order resonances guarantee that typical eigenstates cannot remain localized on the hypercube.

### (d) Properties of resonant subgraphs

In this section, we model the resonant subgraphs generated by second-order processes by independent bond percolation on the hypercube with second-neighbour bonds placed with probability *p*∼*J*^{2}/*W*^{2}. The following properties of resonant subgraphs then follow [60,61]:

(i) For

*J*<*J*_{c}, the number of resonant subgraphs of size*n*per site in Hilbert space is given by 4.7with 4.8*n*_{ξ}is the characteristic number of sites in a resonant subgraph in the localized phase, which diverges at the transition. At the transition and for*n*≪*n*_{ξ}in the vicinity of the transition, the subgraphs are ‘critical’ and are distributed according to a power law. A randomly chosen site in the hypercube lies in a connected cluster of size*n*with probability 4.9(i) At the critical point

*J*=*J*_{c}at finite*N*, we expect equation (4.9) to describe the statistics of cluster sizes with a cut-off size*n*_{ξ}that diverges with increasing*N*. Above the upper critical dimension for percolation, the largest cluster scales with the volume of space to the 2/3 power, 4.10This is well known for the emergence of the giant component in Erdös–Renyi random graphs and likewise for percolation on Euclidean lattices of finite size*L*in dimensions*d*>*d*_{u}=6 [60,62]. Numerical investigation of independent bond percolation on the hypercube (up to*N*=20, not shown) is consistent with this scaling.Thus, at the critical point at size

*N*, the largest resonant subgraphs, while exponentially large in*N*, are still an exponentially small fraction of the hypercube. We expect that such large subgraphs are equally likely to include states in which l-bit*i*is up or down.(i) For

*J*>*J*_{c}, there is a unique giant component which absorbs a fraction*P*∝|*J*−*J*_{c}|^{1}of the total volume 2^{N}of the hypercube. The remaining clusters are of finite size with a distribution ∝*n*^{−3/2}e^{−n/nξ}. In the giant component, all l-bits are equally likely to be up or down, while the remaining clusters have only a finite number of flipped l-bits.

### (e) Consequences

We now assume that the actual many-body eigenstates are simply delocalized over the resonant subgraphs generated by second-order processes. This leads to several quantitative predictions about the localized phase, which agree with numerical observations in §3. On the other hand, the predictions for the crossover region do not explain the numerical data very well and suggest that higher-order processes are important.

(i) The transition from localized to delocalized takes place at

*J*=*J*_{c}(*N*)∼*W*/*N*. This is consistent with our numerical observations in §3. However, as we discuss below, higher-order processes may decrease*J*_{c}logarithmically in*N*.(ii) In a cluster with

*m*edges, the number of flipped l-bits is at most 2*m*(as at most two distinct l-bits can flip on every edge of the cluster, figure 10). Thus, the entanglement entropy of only*O*(*m*) l-bits can be non-zero in a many-body eigenstate delocalized over this cluster. In the localized phase, the clusters are finite and thus 4.11This explains the remarkable collapse in figure 2 in the localized phase. Indeed, the entire distribution of has weight only at zero as , which implies that Δ*S*_{1}∼1/*N*in the localized phase.(iii) Within this model, the participation ratio PR of an eigenstate |

*E*〉 is simply given by the size*n*of the associated cluster. Thus, the PR distribution is given by equation (4.9). This explains several aspects of the numerical PR distributions in figure 8: the PR distribution is nearly*N*-independent for small*JN*and has a sharp cut-off*n*_{ξ}; and the cut-off scales exponentially with*N*as predicted by equation (4.10) in the crossover region. However, the power-law regimes suggested by this model are not clearly seen in the numerics.(iv) At finite

*N*, the localized phase in the vicinity of the transition is predicted to be heterogeneous across eigenstates, l-bits and samples as each initial configuration*τ*provides a different landscape for building the resonant subgraph. However, this heterogeneity vanishes as , as almost all the l-bits are unentangled in eigenstates (see the discussion of [*S*_{1}] above), so this does not explain the numerical observations in the crossover region.(v) Another possible model of the crossover region is that it is described by the bond percolation problem at

*J*>*J*_{c}with a giant component. For eigenstates delocalized on the giant component, , while for eigenstates delocalized on finite clusters, is non-zero only for a finite number of l-bits. This model therefore predicts that the distribution for*S*_{1}becomes bimodal as with*p*(*S*_{1}=0)=(1−*P*) and*p*(*S*_{1}=1)=*P*, where*P*∼|*J*−*J*_{c}| is the fraction of sites belonging to the giant component.However, the numerics presented in §3 do not support this picture. First, there is no evidence of a giant component in the PR distributions in the crossover region in figure 8

*b*. Furthermore, this percolation model does not predict a significant variation of between l-bits*i*in the same state. Numerically, this is the biggest source of the variation in in the crossover region (figure 11).

### (f) Higher-order processes

It is known that the Anderson localization transition for a single particle hopping on a Bethe lattice with large coordination *Z* takes places at hopping strength , where *W* is the bandwidth of the independently sampled disorder on each site [63]. Altshuler *et al.* [52] argued that this result applies directly to the ‘hopping’ problem in Fock space induced by two-body interactions in a disordered quantum dot. They further conjectured a phase diagram in which the eigenstates are localized for , delocalized sparsely in Fock space for and fully delocalized for *W*/*Z*<*J*.

Physically, the transition at arises due to long-distance, high-order resonances in the Fock space. The associated states appear mixed from the point of view of few-body observables, just as the GHZ (Greenberger–Horne–Zeilinger) state does, even though their participation ratios remain a vanishingly small fraction of the total volume of the hypercube. They further would not exhibit level repulsion, as states neighbouring in energy typically do not overlap.

There are a number of assumptions that go into mapping the interacting dot model onto the Bethe lattice localization problem. The most crucial of these are (i) that interference due to the presence of loops on the hypercube is unimportant at long distances and (ii) that the energy correlations are likewise unimportant. These are plausible but not rigorously justified. While there is evidence from several detailed studies of related fully connected spin glass models for the enhancement of the critical coupling by higher-order processes [41–43], the delocalized non-ergodic phase remains controversial. Indeed, even in the Bethe lattice problem, whether this regime is a phase or a slow crossover is still debated [64–66].

The central spin model can be mapped to the Bethe lattice problem under similar assumptions as the quantum dot model of Altshuler et al. [52] with *Z*=*N*. This mapping predicts that higher-order processes destabilize MBL as at asymptotically smaller *J* than the second-order processes described in §4c. If this is true, then it seems that there are two possibilities:

(i) There is a direct transition from the localized to the ETH phase at . If this is the case, then it must be that, at the small sizes

*N*that we can reach numerically, whatever non-perturbative processes cause the sparsely delocalized states to fully thermalize in the limit of large*N*are ineffective at these small*N*, and we primarily observe behaviour governed by the second-order picture up to the scale*J*_{c}∼*W*/*N*. This is why we develop a detailed picture of the second-order percolation model in §4e.(ii) There is an intermediate delocalized ergodic phase between a transition at and another transition on the scale

*J**∼*W*/*N*. This would imply, for example, that [*S*_{1}] transitions at asymptotically smaller coupling than [*r*] by this factor of .

We have attempted to interpret our numerical results from both of these viewpoints in §3, but the results are inconclusive, so we leave this difficult question about the possibility of such an intermediate phase in our central spin model for future work.

## 5. Discussion

We have explored a central spin model for many-body localization and thermalization. We show that a single central spin can serve as a ‘seed’ for interactions between otherwise non-interacting l-bits and thus produce a thermal bath. From the point of view of finite-dimensional disordered systems, our model captures the interaction between a small inclusion and the layer of its immediate neighbours. As the number *N* of immediate neighbours grows with dimensionality *d*, the coupling necessary for this layer to thermalize tends to zero. This bath, ‘nucleated’ by just one central spin, can then thermalize further layers following a ‘layer-by-layer’ version of the arguments presented in [38] and accordingly destabilize MBL.

We show, based both on numerics for small systems and on perturbative arguments, that a central spin system with *N* l-bits exhibits an MBL phase and phase transition in the limit of large *N* if we scale the interaction down with *N*. This dynamical transition occurs at an interaction strength that is thermodynamically insignificant in the large-*N* limit. We argue that the MBL phase can be understood as localization on a hypercube of Fock states [52]. We discuss and explore the possibility of an intermediate ‘delocalized but non-ergodic’ phase between the MBL and thermal phases. Our numerical results on small systems do not allow us to draw a conclusion about whether or not such an intermediate phase exists in this model. Our results also apply to more realistic models of the l-bits that include diagonal interactions between them.^{6}

We have also examined the finite-size behaviour near the phase transition, particularly of the single-l-bit entanglement within eigenstates. We find that the probability distribution of this quantity is bimodal near the transition, with peaks near full entanglement and near zero entanglement. This heterogeneity is shown to occur both across l-bits within single eigenstates as well as across eigenstates for a single l-bit. This is a feature of the finite-size data that has also been seen in one-dimensional models and is not yet well understood [50,51]. The finite-size data for our central spin model do not scale well, indicating that the sizes we can access numerically are well short of any asymptotic large-*N* scaling regime. One indication of this is that the sample-to-sample variations are small, but appear to be accelerating in their growth with increasing *N*.

Let us speculate on possible scenarios for the phase transition out of the localized phase into a delocalized phase (which may be non-ergodic or obey ETH). In the limit of large systems, each sample *s* has its own transition point *J*_{c}(*s*). Owing to the quenched disorder, we expect the probability distribution of these transition points to have a width [56,57]. Within a single finite-size sample, the transition may be sharper than this, of similar width, or broader than this. The apparent accelerating increase of the sample-to-sample differences seen numerically (figure 7*d*) argues against the last of these possibilities. This leaves the two other possibilities. At the accessible system sizes, the inter-sample and within-sample variations are comparable; if this persists to large *N*, we would expect the rounding of the transition to be of width both within samples and between samples. The decreasing trend of the peak with *N* in figure 7*c* suggests that the inhomogeneity between eigenstates in the l-bit-averaged entanglement at the transition might go away at large *N*, thus making all eigenstates in one sample similar in this respect (although it seems likely they will still differ in which l-bits are less versus more entangled).

The other possible scenario is that the transition in a single sample is sharper by a power of *N* than the variation in *J*_{c}(*s*) between samples. A limiting case of this scenario in equilibrium is provided by first-order transitions with quenched randomness [67,68]. In this scenario, if we scale *J* by the sample-averaged *J*_{p}, then at *J*_{p} almost all samples will be away from their transition and be in their MBL or delocalized phases, respectively. Almost all the variation in *S*_{1} then comes from sample-to-sample variations. We are clearly not in this regime for the sample sizes we study numerically here, but it remains possible that this happens at much larger *N*. To study the single-sample transition in this regime, one should instead scale *J* by the sample-specific transition point *J*_{c}(*s*). It is an interesting question to ask how sharp this transition could possibly be. What would be the nature of the sharpest possible, thus most discontinuous, MBL-to-delocalized phase transition? We leave these questions for future research.

## Data accessibility

The code to generate the data reported here is available at https://bitbucket.org/ppontem/centralspinmbl.

## Authors' contributions

P.P. carried out numerical simulations and data analysis. All authors contributed to the writing of the manuscript.

## Competing interests

The authors declare that they have no competing interests.

## Funding

P.P. acknowledges financial support from Fundação para a Ciência e a Tecnologia (Portugal) through grant SFRH/BD/84875/2012 and the Perimeter Institute for Theoretical Physics. Research at Perimeter Institute is supported through Industry Canada and by the Province of Ontario through the Ministry of Research & Innovation. C.R.L. acknowledges support from the Sloan Foundation through a Sloan Research Fellowship and the NSF through grant PHY-1520535.

## Acknowledgements

We thank Wojciech de Roeck, Sarang Gopalakrishnan, Francois Huveneers and Vedika Khemani for helpful discussions.

## Appendix A. Second-order processes

At second order in perturbation theory in *B*, there are configurations, connected to an initial state . These correspond to flipping a pair of l-bits *i*<*j* and possibly flipping the central spin (*s*′) relative to its local field. There are four channels contributing to the second-order amplitude:
A 1where Δ*E*_{s′τ′j}=*E*_{s′τ′j}−*E*_{sτ} are the changes in energy relative to the initial state.

Naively, the final energy denominator, Δ*E*_{s′τ′ij}, takes ∼2*N*^{2} values on a band of width *W* and thus has density of states ∼2*N*^{2}/*W*. The intermediate energy differences, , are typically *O*(*W*) and the numerators are *O*(*J*^{2}), so large amplitudes (resonances) arise in equation (A 1) for ∼*J*^{2}*N*^{2}/*W*^{2} final states. This predicts *O*(1) second-order resonances when
A 2

It turns out that this argument is essentially correct for those final states in which the central spin is flipped, *s*′=−*s*. However, destructive interference between the four channels reduces the probability of resonance for final states in which the central spin is not flipped, *s*′=*s*. Indeed, such cancellations must prevent second-order resonance if the two l-bits involved were non-interacting and thus in product states. Since the interaction terms in this model are *O*(*J*), care must be taken to check the interference for *J*∼*W*/*N*.

More precisely, consider the interference between two second-order channels with intermediate states 1 and 2, respectively, to reach a final state 12,
A 3where Δ*E*_{12}=Δ*E*_{1}+Δ*E*_{2}+Δ^{(12)} and are the numerators of the respective processes. The amplitude can be rewritten as
A 4The first term does not involve the final energy denominator Δ*E*_{12} at all and thus does not cause second-order resonance (assuming that the first-order differences Δ*E*_{1},Δ*E*_{2} are themselves not resonant). Resonances only arise due to interaction energy Δ^{(12)} and matrix element deviations *η*. If both of these are parametrically small in *N*, then resonances are suppressed.

Let us apply equation (A 4) to the determination of resonances in equation (A 1). There are two types of final states depending on whether the central spin is flipped (*s*′=−*s*) or not (*s*′=*s*) with respect to the initial configuration.

*Central spin unflipped ( s′=s)—*In this case, the four channels destructively interfere in pairs according to the intermediate central spin state .

For the pair of channels with ,
A 5so that
A 6To leading order in 1/*N*, the matrix elements
A 7and
A 8where we have used that for *τ*′=*τ*′_{i},*τ*′_{j},*τ*′_{ij}. We have also introduced the shorthand notation for the matrix elements of *B*_{i} with respect to the initial central spin field ,
A 9

For the pair of channels with ,
A 10which provides a larger interaction energy as compared to equation (A 6):
A 11The matrix elements are:
A 12and
A 13For the case of GOE random matrices *B* treated in this paper, this simplifies to as .

With reference to equation (A 4), the largest amplitudes arise from the Δ^{(12)} interaction term for the pair of channels. Its contribution to the amplitude exceeds 1 when
A 14and
A 15Thus, at *J*∼*W*/*N*≪*W*/*N*^{5/6}, there are no second-order resonances for final states with the central spin unflipped.

*Central spin flipped ( s′=−s)—*In this case, rather than being paired by the intermediate state of the central spin , it is natural to pair the channels by whether the central spin flips simultaneously with the

*i*’th l-bit or with the

*j*’th. This is because these pairings lead to the smallest interaction energy Δ

^{(ij)}and correspondingly most destructive interference: A 16and A 17The interaction energy is only Δ

^{(ij)}=

*O*(

*J*), which is again small for

*J*∼

*W*/

*N*.

However, unlike the previous cases, the matrix elements for the four channels are quite different. For example, in the first pair (central spin flips with *i*):
A 18and
A 19where the last step follows because and are independent Gaussian random numbers of *O*(1).

This suggests that the amplitude in equation (A 4) can become large because of the *η* term already if
A 20that is, for
A 21as expected by the naive argument.

Let us check that no further cancellations arise from summing all four channels. For simplicity, we may assume that the interaction energies Δ^{(ij)} are zero, as this only increases the destructive interference. In this case, we have that Δ*E*_{sτ′i}≈−Δ*E*_{−sτ′j} and Δ*E*_{−sτ′i}≈−Δ*E*_{sτ′j}. Thus, the sum of the two most divergent (*η*) terms in equation (A 4) is (to leading order)
A 22The two terms in the square brackets are independent random variables of *O*(1). Thus, they typically add in quadrature and there is no cancellation. This term is therefore ∼*J*^{2}/(*W*/*N*^{2})×1/*W*, which reproduces the estimate that *J*∼*W*/*N* produces *O*(1) resonances.

## Appendix B. Central spin entanglement

The focus of this appendix is the behaviour of the central spin entanglement in eigenstates *S*_{cs}. At *B*_{i}=0, the eigenstates are product states between the central spin and the l-bits with *S*_{cs}=0 in all eigenstates and samples. At small *B*_{i}, the central spin flips in the resonant processes that hybridize the classical states (see §4). As the central spin is substantially entangled in eigenstates that are spread over a few classical states, we expect that *S*_{cs} is close to one in eigenstates with even a few resonant l-bits, while it is close to zero in eigenstates with no resonant l-bits. This has two important consequences: (1) the moments of the distribution of *S*_{cs} are non-zero in the MBL phase even in the thermodynamic limit and (2) the central spin crosses over to high entanglement at *J*^{cs}_{p}(*N*)<*J*_{p}(*N*) within the MBL phase. Note, however, that the scaling of *J*^{cs}_{p}(*N*) with *N* is the same as that of *J*_{p}(*N*); that is, *J*^{cs}_{p}(*N*)∼*W*/*N*.

Feature (1) is already visible in the mean [*S*_{cs}] of the central spin entanglement (figure 11*a*) where the *N* dependence only emerges for *J*∼*J*_{p}(*N*). The total standard deviation Δ*S*_{cs} brings out feature (2) more clearly (figure 11*b*) as the peak lies at *J*∼*J*^{cs}_{p}(*N*)<*J*_{p}(*N*) to the left of the global transition point. Furthermore, the peak value of Δ*S*_{cs} exceeds 0.28 (the standard deviation of a uniform distribution), showing that the distribution is bimodal in the vicinity of *J*^{cs}_{p}(*N*). Finally, in the delocalized phase at *J*/*J*_{p}(*N*)>1, Δ*S*_{cs} decreases with increasing *N*, in agreement with ETH.

In order to bring out the sharp contrast between Δ*S*_{1} of the l-bits (figure 5) and Δ*S*_{cs} of the central spin, we plot Δ*S*_{cs} normalized by its peak value versus *J*/*J*^{cs}_{p}(*N*) in figure 11*b*. Unlike the l-bit entanglement entropy variation, Δ*S*_{cs} does not narrow on its own scale with increasing *N* in the MBL phase. Hence the central spin’s entanglement undergoes a crossover within the MBL phase, as opposed to the l-bits for which the peaks narrow as *N* is increased as seen in figure 5.

## Footnotes

One contribution of 10 to a discussion meeting issue ‘Breakdown of ergodicity in quantum systems: from solids to synthetic matter’.

↵1 Here detuning corresponds to the

*O*(1) mismatch of energy configurations accessible with local moves.↵2 We refer to the thermalizing phase as an ETH phase because it satisfies the eigenstate thermalization hypothesis [45–49].

↵3 The level spacing ratio

*r*is defined as the ratio of consecutive level spacings*r*(*n*)=min(*δ*(*n*),*δ*(*n*+1))/*max*(*δ*(*n*),*δ*(*n*+1)) with*δ*(*n*)=*E*_{n}−*E*_{n−1}when the energies*E*_{i}are enumerated in increasing order.↵4 By contrast, in the one-dimensional MBL phase, Δ

*S*_{1}is non-zero as as the entanglement entropy of any given l-bit is an*O*(1) function of its local environment and the choice of state. This can be seen explicitly in, for example, the one-dimensional phenomenological models in [8,9].↵5 When the expected branching number

*K*becomes too large, we expect the tree approximation to break down (as short loops in the interaction graph become very important). As the critical branching number*K*_{c}is of order one, our treatment is self-consistent within the MBL phase.↵6 With diagonal interactions between the l-bits, first-order processes can already produce a divergent subgraph of resonant configurations for

*J*>*J*_{c}(*N*).

- Accepted September 11, 2017.

- © 2017 The Author(s)

Published by the Royal Society. All rights reserved.