{'fraction_non_alphanumeric': 0.04284874114319321, 'fraction_numerical': 0.01215447951382082, 'mean_word_length': 4.359951966376464, 'pattern_counts': {'":': 0, '<': 2, '<?xml version=': 0, '>': 8, 'https://': 0, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 2, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
Spatio-temporal heterogeneity of entanglement in many-body localized systems

We propose a spatio-temporal characterization of the entanglement dynamics in many-body localized (MBL) systems, which exhibits a striking resemblance with dynamical heterogeneity in classical glasses. Specifically, we find that the relaxation times of local entanglement, as measured by the concurrence, are spatially correlated yielding a dynamical length scale for quantum entanglement. As a consequence of this spatio-temporal analysis, we observe that the considered MBL system is made up of dynamically correlated clusters with a size set by this entanglement length scale. The system decomposes into compartments of different activity such as active regions with fast quantum entanglement dynamics and inactive regions where the dynamics is slow. We further find that the relaxation times of the on-site concurrence become broader distributed and more spatially correlated, as disorder increases or the energy of the initial state decreases. Through this spatio-temporal characterization of entanglement, our work unravels a previously unrecognized connection between the behavior of classical glasses and the genuine quantum dynamics of MBL systems.


I. INTRODUCTION
The assumption of local equilibrium is at the core of statistical mechanics: even if isolated from the rest of the universe, a generic many-body system is expected to act as a thermal bath for itself, quickly driving the statistics of local observables to the Gibbs ensemble, by means of classical [1][2][3] or quantum [4,5] chaos. The situations in which ergodization fails, and the system persists in non-thermal states for all relevant times, are therefore of paramount interest, both at the classical and the quantum level [5][6][7].
In this work, we go beyond those approaches, focusing on the local properties of entanglement, therefore aiming at characterizing its spreading in a more detailed way. We focus on the combined temporal and spatial behavior of the local entanglement as measured by the concurrence, which has been used to characterize entanglement in MBL before (but only at the global level) [54], and which can be measured experimentally [55,56]. Our starting point is the finding that the dynamics of the concurrence among couples of spins or l-bits is highly heterogeneous, with a wide range of different relaxation times. To quantify this observation, we investigate the distribution of the relaxation times τ i of the concurrence in the l-bit model, and describe its properties in a wide range of parameters and initial-state energies. Our main findings are the following. First, the strong fluctuations of entanglement at the local level manifest as a powerlaw-tailed probability distribution for τ i which, in turn, is at the origin of the known power-law decay of the average concurrence [54]. Second, we show that the width of such distribution increases as disorder increases or energy arXiv:2108.05594v3 [cond-mat.dis-nn] 16 May 2022 decreases. Third, we illustrate that the local relaxation times τ i are spatially correlated, with the correlations growing as disorder increases or energy decreases: this is a counter-intuitive result, since, naively, one might expect the correlations to increase when approaching the delocalization transition by lowering the disorder. It is worth stressing that this latter result suggests that the associated entanglement correlation length represents a previously unrecognized length scale in the MBL phase.
These results provide a connection with the so-called dynamical heterogeneity observed in classical amorphous materials and spin glasses [10,11,[57][58][59][60][61], adding one more point of contact between the phenomenology of glasses and that of MBL systems [62][63][64][65]. One speaks about dynamical heterogeneity when each local degree of freedom presents an autocorrelation function that decays in time with a different functional form, thus strong spatio-temporal fluctuations are present in the system. It has already been argued that dynamical heterogeneity in classical glass models has a quantum counterpart and, at ultra-low temperatures, can be induced by quantum fluctuations [66][67][68][69]. However, previous studies have investigated quantum glass systems modeled on a classical counterpart. Here, instead, we adopt a different approach: we borrow the theoretical tools of classical glass theory in order to explore the purely quantum MBL phase and study the features of entanglement, which is a genuine quantum feature with no classical analog. In this perspective, the observed similarities in the heterogeneous behavior of local correlations, which are quantum in the MBL case and classical for glasses, appear unexpected. Importantly, entanglement heterogeneity persists also outside the MBL phase. Here, we choose to focus on the deep MBL phase only because, thanks to the l-bit model, our numerics can reach larger system sizes, and we can also get some analytical insight. We defer the detailed study of entanglement heterogeneity out of the MBL phase, and across the MBL-thermal transition, to future studies.
The paper is organized as follows. In Sec. II, we describe the models under study: the XXZ chain with disorder and the l-bit model deep in the MBL phase. In Sec. III, we introduce the quantities we investigate, i.e. the on-site concurrences, and explain how we evaluate their relaxation times and the correlations among such times. In Sec. IV, we show the results on the distribution of the relaxation times: in the MBL phase, we find that they have a power-law tail and that we can analytically predict the dependence of the typical relaxationtime value on the disorder strength. In Sec. V, we discuss the results for the correlation function of the relaxation times, showing that their correlation length increases with the disorder strength. Finally, in Sec. VI we draw our conclusions and discuss future research perspectives.

II. MODEL
We aim at studying the general properties of the spatio-temporal entanglement dynamics of MBL systems. For this purpose, we focus on an effective description in terms of LIOMs, which allows us to access the nonequilibrium real-time dynamics of MBL systems for long times and large system sizes. Deep in the MBL phase, Hamiltonians of short-range interacting quantum spin-1/2 degrees of freedom can be diagonalized through a quasi-local unitary transformation [15,31], yielding a representation of the model in so-called l-bit form: where {σ x i , σ y i , σ z i } are the localized spin-1/2 operators associated with the LIOMs. We neglect further terms in the Hamiltonian which comprise n-body interactions with n ≥ 3, which is a controlled approximation for weakly interacting spins in the original microscopic model. From analytical results [15,28,30,31], it is known that the interactions J ij are exponentially suppressed with the distance r ij between localization centers. To achieve a model-independent effective description, we parameterize the l-bit model as follows. We assume that the h i are independent identically distributed random fields, with a uniform distribution over [−h, h], and that the J ij are uncorrelated Gaussian variables of zero average and standard deviation J 0 e −rij /κ . For numerical purposes, we set h = J 0 = 1.
The particular advantage of the l-bit model (1) is that it allows us to perform analytical estimates of few-body observables, and to efficiently compute them numerically, reaching system sizes up to L = 140 spins for very long times. To give an example, Eq. (D1) in App. D shows how to compute a two-site correlation function in O(N ) steps. We refer the reader to Refs. [36-38, 40, 54, 65] for more general observables.
It is worth stressing that the l-bits become closer and closer to the physical spins as the disorder increases, ultimately coinciding asymptotically at infinite disorder [33,34]. Thus, at small values of κ (i.e. large disorder strength), one can safely consider the l-bits as uniformly spaced on a chain, and compute the distances between them as r ij = |i − j|, i, j = 1, 2, . . . , L. Our numerical results are obtained exactly in this strongly localized regime, deep in the MBL phase. It is important to note that the effective model allows us to tune: i) the interaction decay length κ (equivalent to varying the disorder strength); and ii), the initial condition, i.e. the energy density at which we probe the system's properties. Concerning the latter parameter, we choose as initial state of the dynamics a product state in the effective spin basis: where | ⇑ i ,| ⇓ i are the eigenstates of σ z i , and |A i | 2 + |B i | 2 = 1. Employing Eq. (2), the system is initially prepared in a superposition of eigenstates. Moreover, Eq. (2) provides us with the flexibility to tune the coefficients A i and B i such that we can vary the initial-state energy expectation value E := ψ 0 |H l-bit |ψ 0 . This tuning can be achieved using a classical simulated annealing algorithm (see Appendix A for details), and allows us to explore different regions of the energy spectrum. We measure E in units of the standard deviations of h i and j J ij , defining the dimensionless energy density ε := (E/N )/ h 2 /3 + 2J 2 0 /(e 2/κ − 1). Notice that ε = 0 corresponds to the center of the spectrum, while ε ≈ −1 to the ground state (more details in App. A). Let us remark that the localization properties of MBL systems depend on the energy of the considered state and are stronger near the edges of the spectrum [23,70].
For small system sizes, we will compare the results of the effective model with a full microscopic calculation for the spin-1/2 XXZ chain with random fields:  [33,34] for the relation between W and the effective model parameters h, κ and J 0 ). When employing (3), we probe the centre of the energy spectrum initializing the system in a Néel state |ψ 0 = | ↑↓↑↓ . . . where ↑, ↓ indicates the physical spin basis, and average the results over different disorder realizations. Notice that the XXZ model can be employed to explore the presence of entanglement heterogeneity also in the thermal phase; however, due to the smallness of the accessible system sizes, we will not investigate thoroughly the thermal phase in the present work. We defer such investigation to future studies.

III. METHODS
For the purpose of exploring the spatio-temporal heterogeneity of entanglement in MBL systems, we concentrate on the two-site concurrence, which quantifies the pairwise entanglement between two qubits [42,43,46,54]. For two spins-1/2 located at lattice sites i and j, the concurrence is defined as [42,46] λ 2 a are the eigenvalues of the matrix R ij = ρ ij (σ y ⊗ σ y )ρ * ij (σ y ⊗ σ y ) sorted in descending order, where ρ ij is the two-site reduced density matrix, and the complex conjugation is done in the standard computational basis.
While the general formula (4) must be applied to the l-bit model, for the microscopic XXZ Hamiltonian the concurrence can be computed more easily. Since the dynamics conserves the z component of the total magnetic field, S z tot , and we initialize the system in the Néel state, having S z tot = 0, Eq. (4) can be simplified as [71] and It has been shown that in MBL systems the concurrence averaged over all couples i, j decays in time as a power law [54], while it decays exponentially fast in ergodic systems. Our key goal is to establish a more detailed spatio-temporal analysis of the concurrence beyond its averaged value. Specifically, we wish to investigate, in MBL systems, the relationship between the late-time, power-law behavior of the average concurrence and the relaxation of the concurrence at the local level. Because of the presence of quenched disorder, we expect to observe that the relaxation of the two-site entanglement is heterogeneous, i.e. it has a different functional form and characteristic time scale in distinct spatial regions. We aim at verifying and quantifying such entanglement heterogeneity. Similar questions have been already explored in the framework of classical glasses; in this work, we will adapt some tools and ideas developed in that context to the case of quantum localized systems.
In this perspective, we define a local on-site concurrence as quantifying the total amount of two-qubit entanglement of i with all the other lattice sites. In the case of the XXZ model, we find that C i,j 0 for |i − j| > 1, so that we can trade the sum in Eq. (7) with the nearest-neighbour term: C i (t) := C i,i+1 (t) (see also [54]).
For large systems and for a single disorder realization, we find that the local concurrence defined in Eq. (7) typically decays to zero on a certain time scale, and then definitely remains so (see Fig. 1). This motivates us to define the local relaxation time as where t 0 = J −1 0 . Notice that C i ≥ 0, so the averages above are well-defined and independent of t 0 (in the thermodynamic limit). The definition (8) employs the logarithm ln(t/t 0 ); this ensures that τ i is a good estimator of the typical time scale of the relaxation time of the concurrence even if C i (t) decays very slowly [72]. Notice that typically, for finite systems, C i (∞) O(2 −L ) [65]: thus, the function C i (t) might be interpreted as a probability distribution over R only in the thermodynamic limit.
Our aim is to estimate the distribution function of τ i exactly in the thermodynamic limit. Such limit can be approached by fixing the maximum simulation time of the dynamics t fin and increasing L, until convergence is reached. We find from our numerics that this typically happens for L 30, which is achievable in the l-bit but not for microscopic Hamiltonians by means of exact diagonalization. Indeed, for the microscopic XXZ Hamiltonian in the MBL phase and for finite L, one finds a spurious peak in the distribution of τ i due to those realizations of C i (t) which are still nonzero at the final evolution time t fin . This is precisely due to the fact that, for small L and whatever the choice of t fin , there will be always a finite number of such nonvanishing realizations. Unfortunately, for the XXZ model, we cannot consider system sizes larger than L = 20. This represents a crucial argument for the use of the effective l-bit model.
Let us emphasize again that the average concurrence in the MBL phase decays as a power law, while the on-site concurrence vanishes after a finite time, as depicted in Fig. 1. Thus, the power-law decay of the average is nothing but a consequence of the heterogeneous behavior at the local level. This can be also captured by a simple analytical model, as sketched in the following. Let us schematize the on-site concurrence as a step function f τ (t) = θ(τ − t), where τ is a random variable drawn from a power-law-tailed distribution. Assuming P (τ ) = N θ(τ − τ 0 )τ −γ , with γ > 1 and N a normalization constant, one gets a disorder average f (t) = (τ 0 /t) γ−1 that decays as a power law. This shows that, when f τ (t) has a simple form and vanishes after a finite time in each disorder realization, and when the relaxation times have a long-tailed power-law distribution, the average f (t) decays as a power law. We expect this argument to be at the origin of the powerlaw decay of the average concurrence in MBL systems.
While the τ i 's provide us with temporal information of the entanglement dynamics, we are further interested in the spatial component. For that purpose we quantify the spatial correlations of the local relaxation time via (see also Ref. [57]) where • is denotes the average over different initial states, [•] |i−j|=r the average over all sites i, j separated by a distance r, and • the average over different disorder realizations [73]. In Appendix E, we show that G τ (r) as defined in Eq. (9) is very robust to finite-size effects and disorder fluctuations: it is a self-averaging quantity. From our numerical simulations, we find that G τ (r) experiences in general a stretched-exponential decay as a function of r. This allows us to define a length scale η τ by performing a fit of the form log G τ (r) ∼ a + (r/η τ ) b for some suitable a and b. The length η τ quantifies the distance over which the local entanglement relaxation is spatially correlated, i.e. it gives the size of the typical clusters of fast or slow entangling spins. In App. E, we show also that η τ is almost independent of the system size for L ≥ 40.

IV. DISTRIBUTIONS OF THE RELAXATION TIMES
We show in Fig. 2 the probability distribution function (pdf) of ln τ i , obtained within both the XXZ and the lbit model. We see that within the XXZ model (Fig. 2a) the pdf's show a peak at large relaxation times, corresponding to the final simulation time of the dynamics t fin . In Appendix B we argue that this feature is due to the (typical) asymptotic value C i (∞) O(2 −L ); see also the discussion below Eq. (8). If the time spent in such asymptotic region is too large, the relaxation time is heavily influenced by the final time of the dynamics. This is a finite-size effect, and it does disappear upon considering larger system sizes, as we show for the l-bit model in Appendix B 2 (larger system sizes for the XXZ model cannot be presently considered).
The pdf's obtained considering the l-bit model for L = 80 and for different values of κ, and ε are shown in Figs. 2b-2c. Thanks to the large system size, these plots do not present any peak at large times, and clearly show that the the pdf of ln τ i has a power-law tail; thus the distribution of τ i has a power-law tail as well. We see that the pdf's become broader as the disorder is increased (both in the XXZ and the l-bit model), or the energy is lowered (in the l-bit model).
We define the typical value of τ i as typ[τ i ] := t 0 exp ln(τ i /t 0 ) τi , where • τi is the average over the pdf of τ i . In Fig. 3a, we show the behavior of typ[τ i ] as a function of the parameters κ and ε. Following usual arguments for the l-bit model [36][37][38]40], in Appendix D and ε, averaged over at least 4000 disorder realizations and 20 initial states for each of them. As κ decreases, i.e. disorder increases, the distributions broaden; the same when ε decreases, in analogy with classical amorphous materials approaching the glass transition. We performed power-law fits on the tails of the pdf's, obtaining the exponents β whose behavior is shown in Fig. 3b.
we derive the rough estimate In Fig. 3b we plot the power-law exponents β obtained from the fit of ln P (ln τ i ) ∼ −β(ln τ i ) shown in Fig. 2b. We see that β has a roughly linear dependence on κ, a property which will help us interpreting the behavior of the correlation function in the next Section.
Before moving to the study of the correlation function, let us emphasize that dynamical heterogeneity is not restricted to the MBL phase. In App. C, we present some qualitative results also in the ergodic regime. Further in- vestigations in this direction promise to be fruitful; however, they go beyond the scope of this work and will be discussed in forthcoming publications.

V. SPATIAL CORRELATIONS OF THE
RELAXATION TIMES Fig. 4 shows the spatial correlations between τ i 's. Due to the strong finite-size effects for the XXZ model, we restrict ourselves to the l-bit model. In Fig. 4a we report the spatial distribution of the τ i 's for a disorder realization. As κ decreases, i.e. the disorder increases, the relaxation times of the local entanglement become spatially correlated over longer distances. The correlation function G τ (r), defined in Eq. (9), is shown in Fig. 4b-4c: G τ (r) decays more slowly upon decreasing κ and ε, confirming the pattern observed in Fig. 4a. The same result is also supported by the (qualitative) behavior of the dynamical correlation length η τ as a function of κ. We see in the inset of Fig. 4c that η τ decreases when κ increases, i.e. when disorder decreases.
The implications are twofold. Recall that the local entanglement spreading slows down when κ decreases (i.e. the disorder increases) or the energy decreases (Fig. 2). Therefore, first, increasingly larger clusters of spins emerge, in which the entanglement relaxation is correlated (Fig. 4). Second, since the distribution of relaxation times becomes broader as disorder increases, more clusters are likely to assume an extreme value of the relaxation time in the slow, as well as in the fast tail.
These findings might seem surprising in the quantum case, as a more localized structure might be expected  9), for L = 140, and various κ and ε. We see that the spatial correlation among τi's increases for decreasing κ and ε. Data averaged over at least 1000 disorder realizations, and 20 initial states each. (Inset of (c)) The dynamical correlation length ητ , from stretched exponential fits of Gτ (r), as a function of κ. ητ decreases as κ increases: for larger disorder, larger clusters of dynamically correlated spins emerge.
when disorder increases. Remarkably, a similar growth of dynamical heterogeneous clusters takes place in classical amorphous materials and spin glasses [57], suggesting a connection between the phenomenology close to the glass transition and to the MBL one. We can provide an analytical argument to justify why the correlation length should increase for increasing disorder as follows. An MBL system develops increasing correlations as the time goes on, as witnessed by the entanglement entropy [18,35,36,49]. It is well established that, for finite systems, there is a time t c at which the system becomes fully correlated: therefore, one can roughly assume that relaxation times τ i ≥ t c are more spatially correlated with each other than the times τ i < t c . If, decreasing κ, the fraction of times larger than t c increases, then one can argue that the relaxation times become more correlated, in agreement with the increasing relaxation-time correlation length η τ that we observe.
Let us substantiate the above argument. First, we estimate t c by quantifying the degree of correlation in the system via the half-system entanglement entropy: it increases in time as S(t) ∼ κ ln(t) [18,35,36,49] and, for a finite system of size L, it saturates to a value of order L at a time t c exp(L/κ). For times τ i ≥ t c , the system has reached its maximum correlation, as we have stated before. Now, let us evaluate the fraction f > of relaxation times τ i > t c as i.e. f > quantifies the weight of the relaxation-time distribution P (ln τ i ) corresponding to relaxation times τ i larger than t c . Finally, we can show that f > increases with decreasing κ, meaning that the relaxation times become more correlated for increasing disorder. This simply follows from the results of Sec. IV, for which ln P (ln τ i ) −β ln τ i , with exponent β = Aκ − B (A 2.8 and B 1, see Fig. 3b). Thus, Eq. (11) reads (12) and f > increases with decreasing κ. This means that the weight of the relaxation-time distribution associated with a maximally correlated system increases when increasing the disorder. This conclusion agrees with the behavior found for the relaxation-time correlation length, η τ . Finally, let us notice that the limit κ → 0 is singular (since t c → ∞), and the argument above ceases to be valid.

VI. CONCLUSIONS
In this work, we studied the spatio-temporal spreading of entanglement in MBL systems by monitoring the onsite concurrence. We showed that in the MBL phase the on-site concurrence behaves heterogeneously, with different functional forms and relaxation times for each site of the system. Using the tools developed for dynamical heterogeneity in classical glasses, we quantified such heterogeneous behavior of entanglement by investigating the on-site concurrence relaxation times, τ i , which display a non-trivial spatio-temporal structure.
First, we observed that the local relaxation times τ i increase upon increasing the disorder, or upon lowering the energy of the initial states. Specifically, their distribution broadens significantly, as the exponent β, dictating the decay P (ln τ i ) ∼ (ln τ i ) −β increases as κ → 0. This can be understood in terms of the slowing down of the dynamics, due to the stronger disorder or the vicinity to the edges of the spectrum.
In addition, taking into account the spatial correlations among the τ i 's, we observed that, as disorder increases or energy decreases, increasingly larger dynamically correlated clusters arise. Within a cluster, the relaxation times are close among sites and, due to the broadness of the τ i distribution, are likely to assume an extremely small or large value. We defined a correlation length of the relaxation times, η τ , which quantifies the typical extension of the correlated clusters. It represents a new length scale characterizing the MBL phase.
The emergence of increasingly large clusters of correlated spins is somewhat surprising since one might naively expect the clusters to grow in size upon decreasing the disorder strength, i.e. when the delocalization transition is approached. A possible explanation lies in the fact that the concurrence quantifies only the entanglement shared by two qubits; therefore, at lower disorder, it may miss the tri-partite, or generally multi-partite, entanglement growth.
That the correlation length η τ of the relaxation times increases with increasing disorder is a subtle effect due to the relaxation-time distributions having a power-law tail, with a significant weight in the long-time regime, where the system is already become fully correlated. For this point, we have provided an analytical justification.
Our analysis was mainly focused on the deep, manybody localized phase, where the l-bits are close to the physical spins. Similar properties for the entanglement heterogeneity are expected in all systems that present a long, localized transient before they reach a thermal state. Such systems include MBL systems coupled to a bath [69,[75][76][77][78][79][80][81], MBL systems in d ≥ 2 [32,51,[82][83][84][85][86][87][88], and two-level systems in structural glasses [65]. In addition, it is worth emphasizing that entanglement heterogeneity should be present also in the ergodic phase; however, we leave its characterization to future studies.
Our findings open up future research directions towards the characterization of spatio-temporal entanglement properties. A next crucial step is to explore the spatial correlation of local relaxation times in observables less affected by finite-size effects and disorder fluctuations. It would be desirable to define suitable macroscopic observables, in the way the four-point susceptibility χ 4 (t) is for classical glasses [10,11,60,61]. In the MBL case, such observables need to detect only local entanglement fluctuations and be, possibly, experimentally measurable. A further interesting extension of our contribution could be to consider entanglement heterogeneity for two subsystems consisting of more than one spin, which could provide additional information on the multipartite spatio-temporal structure of quantum entanglement. Finally, we note that the observed entanglement correlated clusters, which grow in size for increasing disorder strength, might be linked to the localized bubbles which have been the subject of numerous recent studies [89,90]. Exploring such connections promises to be an interesting research direction for future studies. In the l-bit model, given a disorder realization {J ij }, we sample the local magnetization configurations {m i } = { σ z i } with probability ∝ e −E/T , T being a fictitious temperature to be gradually lowered. Since m i ∈ [−1, 1] are continuous variables, the annealing procedure has easy access to states down to the edge of the spectrum. From {m i }, we fix the coefficients of the initial states of the dynamics (see Eq. (2), main text) as A i = (1 + m i )/2, and B i = e iφi (1 − m i )/2. This choice guarantees that ψ 0 |H l-bit |ψ 0 = E, i.e. the quantum initial state has an energy expectation value equal to the desired one.
For what concerns our choice of the energy scale (see Sec. II), namely the reasoning goes as follows. The l-bit Hamiltonian (Eq. (1), main text) can be interpreted as a classical spin glass, if one substitutes σ z i −→ s i = ±1. Then, one can compute the (annealed) density of states of the model, finding that with high probability the ground state is at E = −N h 2 /4 + 4J 2 0 /(e 2/κ − 1) (see also Ref. [91]). Changing the spins to continuous variables σ z i −→ m i ∈ [−1, 1] will just modify the prefactors of h 2 and J 2 0 /(e 2/κ − 1), without changing much the scale. For this reason, we have chosen to put in Eq. (A1) simply the sum of the variances of h i and j J ij . The ground state will not be exactly at ε = −1, but close to it. In view of the strong finite-size effects in the results for the XXZ model shown in Fig. 2a, let us better analyze the probability distribution function (pdf) of τ i .
In Fig. 5 we show the distribution of the local relaxation times of the concurrence, computed within the XXZ model. The tail of the distribution is cut away according to the following procedure. We observe that in some instances the nearest-neighbor concurrence C i,i+1 (t) becomes numerically indistinguishable from 0 at a time t * , and then stays equal to 0 definitively. We perform an evolution lasting only a finite time t fin , so for the finite size we consider there will be many sites and realizations with t * > t fin . This is the reason why the full distributions in Fig. 5 and Fig. 2a in the main text show such a huge peak at large times: it is formed by the contributions of C i,i+1 (t) which have not vanished on the finite-time window t fin of our evolution, for the finite system size we consider.  In order to get rid of this peak, we choose a certain truncation time t tr ≤ t fin , and select only the sites and the realizations for which t * < t tr . As we can see in Fig. 5 the huge peak disappears and there is a large-time tail which depends on the chosen value of t tr . The small-time structure is, on the opposite, quite independent of the truncation, so we expect that it has a physical meaning. There is a peak around ln τ i −1, which resembles the one appearing in the l-bit distributions; however, another peak is present around ln τ i 1. The two-peak structure has no equal in the l-bit model results; we argue that this might be due to the n-body interactions with n ≥ 3 missing in the l-bit model. We see that, when L ≤ 20, the probability distribution presents a peak at τ i = O(t fin ), where t fin is the final time used in the numerical simulations for the time evolution. This is the same effect observed in Fig. 2a. In particular, we see that the shape of the pdf's at L = 10, 15, 20 strongly resembles the behavior observed in the XXZ model, confirming that those results are strongly affected by finite-size effects.
In Figs. 6b-6c we reproduce the pdf's at L = 16, ε = 0, and different values of κ an t fin . In the presence of finitesize effects we do not observe the decay of the pdf that is found in larger system sizes. We observe instead a peak in the distribution at τ i = t fin . We see that, when the system size is too small, the pdf's present a spurious peak at large values of τi. Indeed, entanglement cannot spread over many sites, and the concurrence of some strongly interacting couples remains finite even at infinite times (Ci(∞) = O(2 −L )). We collected data from 20 initial states for, at least, 300 disorder realizations.  The peak corresponding to ergodic realizations of the on-site concurrence is very evident and decreases as the disorder increases and the transition to MBL regime is approached. Numerical parameters: J = 2, V = 1 (transition to MBL for W 10), L = 18, at least 300 disorder realizations for each curve.
Appendix C: Spatio-temporal entanglement heterogeneity in the ergodic regime We now present some qualitative results on entanglement heterogeneity in the ergodic regime. As shown in Fig. 7, broad distributions of the relaxation times appear for the XXZ model (3) also in the ergodic regime, and become broader as the disorder W increases and the MBL-ergodicity transition is approached. Such distributions present a clear peak at small values of τ i , which decreases when disorder increases and the transition to the MBL phase is approached. This peak corresponds to realizations of the on-site concurrence that vanish very fast in time, corresponding to situations in which entanglement spreads quickly in the system. The number of such ergodic sites diminishes as the localization transition is approached, as it is quantified by the decrease in the peak height.
Notice that in the thermal phase the distributions do not show the spurious peak for large τ i , which in contrast appears in the MBL phase due to those realizations in which some C i (t) are still nonzero for t = t fin (contrast Fig. 7 with Fig. 2a).

Appendix D: Analytical estimates of local time scales
Let us focus on the l-bit model and on the computation of typical relaxation-time scale of the "one-site concurrence" Eq. (7). The concurrence is a complicated nonlinear function of the two-site reduced density matrix ρ i,j , therefore it is really hard to make analytical predictions for it. However, one can hope to get a rough estimate of its behavior by considering instead the correlation function σ x i (t)σ x j (t) . This type of correlation functions were already considered in previous works (see e.g. Ref. [38]), and are easy to access. Choosing i = 0 and j = 1 without loss of generality, it explicitly reads where ρ 0,0 and ρ 0,1 are the initial density matrices of sites 0 and 1, and θ j is the azimuthal angle on the Bloch sphere for the initial state of site j. We take a further step, and also simplify θ j ≡ π/2, i.e. we choose a particular initial condition at infinite temperature. As a result, we find that σ x 0 (t)σ x 1 (t) is an oscillating function, modulated by envelopes of the form A ± (t) := j =0,1 cos 4J 0j t ± 4J 1j t . (D2) It is clear that, if we want to understand the leading-order behaviour in time, we can reduce to study the simpler function where, we recall, J 0j are Gaussian variables of zero average and standard deviation w j := J 0 e −|j|/κ .
We can estimate the typical value of A(t) by means of typ[A(t)] := exp ln A(t) (we need to average the logarithm of A because, with hindsight, there is a power-law tail in the relaxation-time distribution). Since Appendix E: Self-averaging property of the correlation function As anticipated in the main text, the correlation function G τ (r) (defined in Eq. (9) of the main text) is a selfaveraging quantity, as depicted in Fig. 8. It is indeed only slightly sensitive to finite-size effects (see Fig. 8a), and to disorder fluctuations (see Fig. 8b). As a consequence, the dynamical correlation length η τ is almost independent of the system size as well (see inset of (a)). (1), main text) for: (a) κ = 1, ε = 0, and various system sizes L, averaged over D = 300 disorder realizations and 20 initial states; (b) κ = 1, ε = 0, L = 140, averaged over a different number of disorder realizations D, and 20 initial states. We see that Gτ (r) converges quickly to its thermodynamic value (panel (a)), and is almost independent of the number of disorder realizations (panel (b)). (Inset of (a)) The dynamical correlation length ητ from stretched exponential fits of Gτ (r) reported in (a) (see main text) as a function of the system size L. We see that ητ is almost independent of the system size for L ≥ 40. The error bars are given by the fit errors.