Measuring out quasi-local integrals of motion from entanglement

Quasi-local integrals of motion are a key concept underpinning the modern understanding of many-body localisation, an intriguing phenomenon in which interactions and disorder come together. Despite the existence of several numerical ways to compute them - and astoundingly in the light of the observation that much of the phenomenology of many properties can be derived from them - it is not obvious how to directly measure aspects of them in real quantum simulations; in fact, the smoking gun of their experimental observation is arguably still missing. In this work, we propose a way to extract the real-space properties of such quasi-local integrals of motion based on a spatially-resolved entanglement probe able to distinguish Anderson from many-body localisation from non-equilibrium dynamics. We complement these findings with a new rigorous entanglement bound and compute the relevant quantities using tensor networks. We demonstrate that the entanglement gives rise to a well-defined length scale that can be measured in experiments.


INTRODUCTION
It is widely believed that generic quantum systems isolated from their environments will evolve under their own dynamics until they reach an apparent equilibrium state that locally resembles the expectations of a thermal equilibrium state [1,2].This expectation is seen as a stepping stone to reconcile predictions from statistical mechanics and those of basic quantum mechanics.One major exception to this rule is the case of low-dimensional quantum systems in the presence of random disorder.Non-interacting quantum systems in one dimension will entirely fail to thermalise due to any finite concentration of disorder [3], and in recent decades it has been shown that interacting many-body systems appear to suffer the same fate [4,5], leading to the phenomenon now known as many-body localisation (MBL) [6][7][8][9][10][11][12].From a theoretical standpoint, MBL is now fairly well understood in terms of the emergence of an extensive number of conserved quantities known as (quasi-)local integrals of motion (LIOMs, also known as localised bits or l-bits) which can prevent many-body systems from reaching thermal equilibrium [7,13].While phenomenological models based around the concept of l-bits have seen great success [14,15], and there are several approaches that can map microscopic mod-FIG.1. Division into subsystems and computation of negativity: a) A sketch showing how a one-dimensional spin chain is partitioned into three subsystems.We are interested in computing the entanglement between subsystems A and B after subsystem C has been traced out, giving rise to a spatially-resolved entanglement measure.b) Sketch of the initial quantum state in matrix product operator (MPO) form, made by taking the outer product of two matrix product state vectors.c) Sketch of how the negativity is computed: the partial transpose of subsystem A corresponds to 'twisting' the MPO legs while tracing out subsystem C corresponds to contracting the relevant MPO indices.els onto effective l-bit models [16][17][18][19][20][21][22][23][24][25][26], the l-bits themselves remain a strictly theoretical construct, inaccessible to any experimental probes.This is in contrast with the case of Anderson localised systems, where the exponentially localised l-bits can be straightforwardly related to the real-space decay of the single-particle states, which has been experimentally observed [27].
In this work, we propose an experimentally feasible approach to measuring the actual real-space properties of local integrals of motion in many-body quantum systems using the entanglement negativity, a sensitive entanglement monotone that allows for the recovery of spatially resolved entanglement information.In this way, we accommodate the above missing link.Various quantities capturing correlations and entanglement, including the negativity, have been measured in recent experiments with ultra-cold bosons: Ref. [28] has measured single-site and half chain number and configurational entanglement for a system subject to a quasi-periodic potential.These can be seen as a witness detecting the absence of thermalization, but they do not provide a length scale.The authors have also measured classical density-density correlationsakin to the proposal of Ref. [29]-showing exponentially decaying correlations, but this is a two-point classical measure, in contrast to the genuine entanglement between two halfregions considered here.Reference [30] has studied a disordered system and has shown that the entanglement negativity can be directly measured.In a first experiment, the authors of arXiv:2301.01787v4[cond-mat.dis-nn]21 Nov 2023 the latter work prepared the system in a product state and measured the two-qubit entanglement of formation as they vary the separation between the qubits.While this setting is close in spirit to our approach, they have chosen a two-qubit setting, which is the only setting in which one can compute this quantity, so that the diagnostic time scale that allows observation of any spatial dependence is short.In a second experiment, the preservation of entanglement has been studied, departing from the approach taken here.Here, we demonstrate that the negativity itself gives direct access to a unique length scale that characterises the l-bits.

RESULTS AND DISCUSSION
Quasi-local operators.The question of whether manybody localisation is a well-defined stable phase in the thermodynamic limit remains unsettled, nevertheless systems showing MBL-like phenomenology for experimentally accessible times and system sizes appear to be well described by l-bit models, and their length scale is physically meaningful regardless of whether this description keeps holding for very long times and very large system sizes.Since the goal of the present work is to characterize the l-bits, we will now give precise definitions of what they are and the role they play in the dynamics.
Definition 1 (Quasi-local operators).An operator O on the lattice Λ is said to be quasi-local around a region R with localisation length ξ if for any region where K > 0 is a universal constant, X C denotes the complement of X, d(X c , R) = min x∈X c ,r∈R d(x, r) is the length of the shortest path from X to R, and Intuitively, this definition says that most of the support of O is concentrated in and immediately around the region R, in the sense that if we truncate O to an operator only supported on a sphere centered at R, the resulting operator differs from O only by an error decaying exponentially in the radius of the sphere.Interestingly, the relatively loose sense of decay in 2norm seems crucial, an insight that is often under-appreciated [31,32].It is important to note that this is an abstract definition: it does not give operational advice on how to find those l-bits.What is more, even if they exist, they are by no means unique [33].There could be "more local" l-bits than those given that still give rise to a complete set of l-bits.Either way, as is common, such l-bits serve as our definition for manybody localisation.

Definition 2 (Many-body localisation). A Hamiltonian
with real weights {ω j } and {ω (2) j,k }, is called many-body localised if it can be written as a sum of mutually commuting ([h j , h k ] = 0 for all j, k) quasi-local terms h j , each centred around site j, and if ω i1,...,in ≤ ωe −|i1−in|/κ for some constant κ > 0, where i In other words, a many-body localized Hamiltonian can be written as an effective model which is classical, in the sense that all Hamiltonian terms commute, and quasi-local, in the sense that the total coupling strength between two sites decays exponentially with distance.
Premise of the approach.When written in the basis that diagonalises the Hamiltonian, as in Eq. ( 2), these l-bits are strictly local objects, but in real-space they are quasi-local with exponentially decaying tails.In order to extract properties of l-bits from experiments, we shall consider the evolution of an arbitrary initial state under the following Hamiltonian dynamics as for times t ≥ 0. To simplify the notation, we will suppress the time argument for time t = 0. How can this time evolution be exploited to measure out real-space properties of the l-bits?Some intuition can be attained in the situation when the {h j } are strictly local.The terms that do not overlap do not contribute to the entanglement evolution at all.So in the end, it is the overlapping tails that will lead to entanglement growth.
Model.We will demonstrate our scheme numerically using the 'standard model' of MBL, namely the XXZ spin-1/2 chain with random on-site fields, while it should be clear that the approach taken would be applicable to any many-body localised model.Its Hamiltonian is given by with We shall set J 0 = 1 as the unit of energy throughout, with ∆ = 1.0 unless otherwise stated, and will use open boundary conditions.This model has been thoroughly studied and shown to exhibit a phase with anomalous thermalisation properties above a disorder strength of d ≳ 3.7 [9], although recent work has suggested that the true phase transition in the thermodynamic limit could be at much larger values of d if it exists at all [34][35][36][37][38].
The characteristic growth in time of the von Neumann entanglement entropy [39,40] (or its correlation-based analogues [29]) mentioned above has been shown to be a good indicator of many-body localisation, able to distinguish it from single-particle Anderson localisation via the late-time logarithmic growth.Motivated by this, our aim in this work is to show that other entanglement measures which provide spatially-resolved information can not only distinguish manybody localisation from Anderson localisation, but can also allow direct quantitative measurement of the properties of many-body local integrals of motion.
Diagnostic entanglement quantity.The main quantity of interest in this work is the logarithmic negativity, a measure of the entanglement between two subsystems of the spin chain, denoted A and B, separated by a distance r, which together with C constitutes the entire system (sketched in Fig. 1).It is defined as [41][42][43][44] where ∥O∥ 1 = tr|O| denotes the trace norm, ρ A,B (t) = tr \{A,B} [ρ(t)] is the time-dependent quantum state of subsystems A and B after tracing out all other lattice sites, and the superscript T A indicates the partial transpose with respect to subsystem A. This has been shown to be an entanglement monotone meaningfully quantifying entanglement [43][44][45].In the following, we shall refer to this quantity simply as 'negativity'.By contrast to the more commonly studied bipartite von Neumann or Rènyi entanglement entropies which consider a single bi-partition between two connected subsystems, the entanglement negativity allows for a meaningful spatially resolved measure of mixed-state entanglement, as the two subsystems can be separated by an arbitrary distance r := dist(A, B), a feature the von Neumann entropy cannot capture as a pure state entanglement measure.This measure can also be used to study the entanglement between subsystems of arbitrary size.However, for conceptual clarity, we shall mainly consider A and B to cover the entire chain except for a piece C with |C| = r+1 separating A and B, as shown in Fig. 1.That said, the concept works as well for small regions A and B, as they are accessible in experiments and are discussed in the rigorous bounds.Numerical evidence is shown in Supplementary Note 5.The negativity has previously been investigated in the context of ground states of disordered spin chains [46], quenches in random spin chains [47], the manybody localisation transition [48], and quench dynamics in the presence of a defect [49].For clarity, in the following, we shall drop the explicit dependence of E N on the quantum state and instead use the notation E N (r, t) to represent the negativity associated to two subsystems separated by a distance r and a time t following a quench from an initial product state, emphasizing that this is indeed a spatially resolving entanglement measure.
A heuristic argument for why this quantity is relevant in our case can be given in the following manner.Ref. [13] has shown that the von Neumann entanglement entropy grows in time following a quench according to S ent ∝ ln(J 0 t/ℏ), once the system enters the late-time equilibration regime.If we wish to consider the entanglement negativity between two subsystems separated by a distance r, a reasonable starting assumption is that the negativity will vary in time according to the same ∼ ln(t) growth but will be exponentially suppressed in magnitude due to the spatial separation of the two subsystems, leading to an overall behaviour of E N ∝ exp(−r/ξ) ln(J 0 t/ℏ).We shall show that this ansatz is a good match for the numerical results.We also wish to emphasise that this logarithmic growth is characteristic of the interacting system and is entirely absent from Anderson-localised systems, meaning that the existence of this length scale is a distinct fingerprint of a many-body localised system.
Corroborating the reasoning with rigorous bounds.We see that Hamiltonians that are many-body localised in the sense of Definition 2 create entanglement at a rate that decays exponentially in the distance r = dist(A, B) between parts A and B, reflecting the exponential decay of the tails in quasi-local l-bits.In fact, not only this intuition can be made entirely rigorous, but, at the cost of slightly weakening the definition of quasi-locality, we are in the position to state precise upper bounds for the negativity for all times and distances.
Theorem 1 (Rigorous entanglement bounds).Let ρ be an initial product state.Let H be a many-body localised Hamiltonian as per Definition 2 with localisation length ξ < 1/(4 log(2)) and 2(1/κ − log(2)) > 1/ξ, consider three blocks A, C, B such that C divides A from B, with |C| = r + 1.The growth of the negativity of the state ρ(t) = e −itH ρe itH restricted to the regions A, B is bounded as 1), (6) for times t ≥ e r/(4ξ) , while for t < e r/(4ξ) , We hence find a short time behaviour signifying a linear growth in time, a cross-over regime governed by the correlation length, and a logarithmic growth for long times.These bounds-interesting in their own right and complementing and refining those of Ref. [50]-are perfectly compatible with the above numerical assessment.In Supplementary Note 8, we state details of the proof of the bound that makes extensive use of the precise form of the tails of the l-bits.Based on our numerical results, we expect that our assumptions on the localisation length and the definition of quasi-locality can be relaxed without affecting the result.We also note that the observed ξ−dependence of the late time decay of the entanglement with the size of C is not visible in this bound, though we expect that it can be refined to show this.
Numerical results.We first discuss qualitatively the results for the growth of the entanglement negativity with time for various different distances r, as shown in Fig. 2a) for a disorder strength d = 8.0 (deep in the localised phase), where we find that indeed the negativity grows logarithmically with time at late times.Results for further disorder strengths, system sizes, and subsystem sizes are available in Supplementary Notes 2,3,4 and 5.At short times, the negativity is dominated by diffusive transport on length scales shorter than the localisation length.At large distances r, the negativity remains close to zero until a time exponentially large in r, which can be used to define a 'light cone' that characterises the spreading of the entanglement negativity, shown in Fig. 2b).The three lines indicate when the negativity grows above a threshold ε ∈ {0.0001, 0.001, 0.005}, mapping out an approximately logarithmic light cone.As the negativity outside of this length cone is exponentially small, in the following analysis, we restrict ourselves to space-time coordinates (r, t), which are within the light cone.The existence of this light cone means that we gain only diminishing returns by going to larger system sizes: although we are able to separate the subsystems by a larger value of r, the evolution time required to obtain meaningful entanglement scales exponentially in r, which incurs a large computational cost for large systems and quickly becomes prohibitive.
In the late-time logarithmic growth regime, where the dynamics are dominated by the quasi-local nature of the l-bits, we extract the value of the negativity at a given time t * following the quench from an initial Néel state and plot it versus the subsystem separation r.We show this in Fig. 2c) for several different choices of time t * [indicated by the dashed lines in Fig. 2a)].The data points form a straight line (on a logarithmic scale), and at late times the gradient of the line does not strongly change with the choice of time t * , appearing to saturate at a fixed value (although the y-axis offset will, of course, Orange lines mark the localisation length obtained through exact diagonalisation following Ref.[31].For further details on calculating the localisation length using exact diagonalisation, see Supplementary Note 7 [51].The black line indicates the localisation length of the corresponding Anderson-localised system, obtained by directly diagonalising the Hamiltonian in the non-interacting limit (∆ = 0), for a system of size L = 128 with Ns = 10000.continue to increase in time).Further details are available in Ref.Under the assumption that the negativity decays exponentially with distance like E N (r, t * ) ∝ exp(−r/ξ), we can perform a linear fit to the data shown in Fig. 2c) and extract a well-defined length scale ξ which characterises the spatial extent of the l-bits.The results are shown in Fig. 3, where we find that the length scale ξ exhibits monotonic decay with increasing disorder strength, as expected.Note that no assumptions are involved other than the exponential decay of the negativity with distance at some fixed time t * : the resulting length scale is an emergent property of the many-body system.This assumption does not hold in the delocalised phase, where the entanglement does not enter a regime of logarithmic growth.We can further compare the length scale extracted from our procedure with the l-bit decay lengths computed using the established numerically exact method of Ref. [31], using the definition of quasi-locality from Eq. ( 1).We find excellent agreement between the entanglement-based length scale and the l-bit localisation length obtained independently from this method, confirming that the length scale probed by the negativity is the localisation length of the l-bits.
For comparison, we also indicate the corresponding localisation length of an Anderson localised system, here obtained by directly diagonalising the Hamiltonian with ∆ = 0 (following a Jordan-Wigner transform into the fermionic representation).We compute the eigenvectors of the Hamiltonian in the non-interacting setting, which decay in real space as exp(−r/ξ) [52], average over disorder realisations and extract the localisation length ξ from a least-squares fit.The length scale extracted from the TEBD data behaves in a qualitatively similar manner to the single-particle localisation length but is always larger, confirming that we are not measuring single-particle properties but are indeed extracting a gen-uinely many-body feature of the system.In the delocalised phase, our assumed form of the negativity is no longer valid, and as such, the method cannot extract a reliable length scale.
We also note that the entanglement negativity is not the only entanglement measure which may be used in this way: any spatially-resolved entanglement probe should behave similarly.In Supplementary Note 6, we demonstrate that the mutual information also gives consistent results.

CONCLUSION
In this work, we have outlined an experimentally feasible procedure for measuring local integrals of motion based on their contribution to the slow growth of the negativity at long times following a quench from an arbitrary initial state.We have demonstrated that the length scale which we obtain from this procedure, which characterises the l-bits, is in good agreement with that obtained using other theoretical methods in the literature.The crucial advantage is that our scheme is experimentally tractable, unlike other purely theoretical/numerical methods, which cannot be verified in real experiments.It would be extremely interesting to apply this method to other scenarios where many-body localisation is believed to exist, such as in disorder-free systems and two-dimensional models, in order to see if well-defined length scales based on the spreading of entanglement may still be extracted in these situations.This work paves the way for the application of spatially-resolved entanglement probes to phenomena in quantum simulation beyond many-body localisation, where such methods may be able to provide valuable insight into emergent length scales associated with other types of quasiparticles.

METHODS
We compute the negativity using time-dependent matrix product state simulations -an instance of a tensor network method [53] -implemented in the Quimb package [54] using the time-evolving block decimation (TEBD) algorithm to perform the evolution [55,56], with the system initially prepared in a Néel state.We use system size L = 24 with a maximum bond dimension of χ = 192.We perform the time evolution using a maximum time step dt = 0.05, at each step discarding singular values smaller than ϵ = 10 −10 .We have checked that the results are well-converged.Detailed benchmarks are shown in Supplementary Note 1.Our TEBD results are compared against l-bit length scales obtained using exact diagonalisation, following Ref.[31].
The negativity can be computed straightforwardly from a matrix product state (MPS) representation [57].The state vector can be turned into a matrix product operator (MPO) (sketched in Fig. 1) representing the quantum state by considering vectors and dual vectors represented as MPS.The partial transpose can be computed by 'twisting' the legs of the MPO tensors, while the partial trace over the subsystem C can be performed by contracting the free indices of the MPO tensors in this subsystem.At long times, the negativity should saturate at a value controlled by the size of the subsystems, and at any time t < t s (where t s is the saturation time), the negativity should satisfy the hierarchy E N (r 1 , t) < E N (r 2 , t) for any two distances r 1 > r 2 .10 0 10 2 0.000 0.002 0.004

Supplementary note 2: Comparison of disorder strengths
In addition to the data presented in the main text, here we show the behaviour of the entanglement negativity over a range of different disorder strengths, demonstrating the logarithmic growth at late times in the localised phase, and the qualitatively different behaviour seen in the delocalised phase.The results are shown in Fig. S2, for a system of size L = 14, bond dimension χ = 128, and averaged over N s = 240 disorder realisations.The circle markers represent the data points, while the solid lines are smoothed guides to the eye.Deep in the delocalised phase (d = 1.0), we see that the negativity saturates for all values of r at relatively early times, making it difficult to pinpoint a regime where the growth of the negativity can be associated with a length scale.In contrast, the negativity in the localised phase increases much more slowly with time, and the spacing of the curves is consistent with an exponential suppression of the negativity with distance, as demonstrated in the main text.In the delocalised phase, the negativity saturates to a value determined by the size of the subsystems A and B, while in the localised phase the negativity displays a slow ∝ log(t) growth even at late times.In this dephasing regime, we are able to use the data shown here to extract a length scale that characterises the localised phase, as detailed in the main text.

Supplementary note 3: Bond dimension
In Fig. S3, we show the negativity dynamics for system size L = 24 and varying bond dimension χ, demonstrating that for χ = 192 (the choice used in the main text) the results are well-converged.The crucial factor for our work is the rate of growth of the negativity, which appears largely unaffected by the choice of bond dimension, although deviations can be seen for the smallest value shown in Fig. S3.

Supplementary note 4:Comparison of different measurement times
In the main text, all results for the length scale ξ ≡ ξ(t * ) are taken with t * = 500, i.e., the maximum evolution time of our simulations, however, it is clear from the negativity dynamics that there should be a weak dependence of ξ on the measurement time t * .Here we demonstrate this effect.Figs.S4 and S5 show how the linear fit used to extract the decay of E N (r, t) against r depends on the choice of time t * for a variety of different disorder strengths, and Fig. S6 shows how the resulting values of ξ(t * ) change as t * and d are varied.At short times there is a visible change in the length scale ξ(t), however, at longer times we see that it appears to saturate towards a well-defined length scale with only a weak dependence on time.It is possible that simulations which extend to longer times may be able to improve upon the results presented here, but our results suggest this will be a meagre quantitative improvement in exchange for a great deal of computational effort.The circular markers show the data, while the solid lines indicate the linear fits used to extract the l-bit localisation length.Error bars showing the standard error in the mean are typically smaller than the marker size.At large disorder strengths (i.e., in the localised phase), the linear fit is very good, confirming that the negativity does indeed decrease exponentially with distance in this phase.After a time of t * ≈ 100, the gradient of the decay (and hence the corresponding l-bit localisation length) does not strongly change with time in the localised phase.

Supplementary note 5: Negativity between subsystems of fixed size
It is also possible to compute the entanglement negativity on a more general interval between two subsystems of fixed size R b separated by a distance r, as sketched in Fig. S7.In this case, a very similar procedure to that proposed in the main text is possible, with the caveat that one must carefully choose both the block size R b and the distance r to ensure that the extraction of the l-bit length scale is done during the dephasing regime.
To be specific, if the block size R b is small and the subsystems are close together (i.e., r is small), then the entanglement negativity will rapidly saturate (as the maximum entanglement is controlled by the size of the subsystems under consideration).S7.We can see that for small values of R b , the negativity at small separations r saturates quickly and that only for larger values of R b does the negativity increase in a manner that allows extraction of the relevant l-bit length scale.Error bars show the standard error in the mean, and the solid lines are a smoothed guide to the eye.

Supplementary note 6: Mutual information
In the main text and in our analytical results, we specified the logarithmic negativity as our chosen spatially-resolved entanglement probe.Here we briefly demonstrate that another spatially-resolved entanglement measure, the mutual information, also exhibits qualitatively similar behaviour.The mutual information between subsystems A and B is defined as where S(ρ) := −trρ log(ρ) represents the von Neumann entropy of a quantum state ρ, and ρ A is the reduced quantum state of subsystem A. Fig. S9 shows the results for a small system of size L = 12 with bond dimension χ = 128, averaged over N s = 100 disorder realisations.The mutual information is qualitatively -and even quantitatively, in many cases -similar to the entanglement negativity, strongly suggesting that it would be a more than acceptable substitute and that much of the intuition developed in the main text should also apply to the mutual information.These findings are compatible with those using the multi-partite mutual information as a probe for many-body localization [60,61].The l-bits in this work are calculated using the method described in Ref. [31].For completeness, we offer a brief explanation of the method: A system in the fully many-body localised phase can be fully characterized by a complete set of l-bits.Let H be the MBL Hamiltonian of a system of size N and τ i , i = 1, . . ., N , be the l-bits, then they need to satisfy the following properties: 1. [H, τ z i ] = 0, 2. [τ z i , τ z j ] = 0 for all i and j, 3. τ i is quasi-local in real space (will be elaborated in the next section).
A simple prescription in Ref. [50] has been given to construct the l-bits out of infinite time averages of terms that the Hamiltonian contains.The infinite-time average of a term h j is assuming a non-degenerate spectrum {E k } of H. Therefore, E(h j ), being sums of projectors onto the eigenstates of H, automatically fulfill properties 1 and 2. The authors have also demonstrated the quasi-locality of the resulting operator, showing that the non-local contributions from the off-diagonal elements, i.e., contributions from τ x and τ y will be removed through the infinite-time averaging procedure.The downside with this approach is that the resulting operator has a degenerate spectrum that is distinct from what the Pauli algebra mandates.Thus, we can no longer associate this operator with the picture of having ladder operators that help us transverse through the different modes/eigenstates.The authors in Ref. [31] have computed E(σ z j ) from ( 10) and re-arranged the order of the eigenstates in U d to minimize the pair-wise differences between the spectrum of E(σ z j ) and that of starting with j = 1 and sequentially optimising the eigenstate order with larger j while keeping the already optimized partial ordering from smaller j intact.The resulting operators will preserve the Pauli algebra by construction while becoming quasi-local in real space.For simplicity, we only deal with trace-less and Hermitian operators for the moment.Operators with non-vanishing traces require special procedures to meet the orthonormality of the Hilbert-Schmidt inner product.
Definition 3 (Quasi-locality).Let τ be a trace-less Hermitian operator, normalized with respect to the Frobenius norm, and L i its orthogonal trace-less Hermitian basis.Consider the decomposition Let S(L i ) be the support of L i in real-space.τ is quasi-local around site j, if and only if for any connected region B containing j, we have This definition, although rigorous, can be cumbersome, and it is not easy to compute the relevant quantities.Let us consider the partial trace of the operator τ , where Li are now defined on a smaller Hilbert space contained in B and the extra identity operators outside B become the pre-factor 2 |B C | after the partial trace.Then, we can compute the square of the Frobenius norm of this truncated operator to get due to the orthogonality of the basis operators L i and their restrictions to any region B. Finally, we leverage the normalisation of τ to observe that Therefore, This is exactly the quasi-locality measure proposed in Refs.[31,33].The spatial decay of the l-bits can also be computed using the weight measure in Eq. ( 6) from Ref. [62].This notion of quasi-locality is agnostic to operator content and where it is quasi-local around because one only needs to compute the weight of the operator filtered at each site.

Supplementary note 8: Entanglement growth bound
In this section, we prove an entanglement growth bound that has been tailored for the case at hand.Generic entanglement growth bounds for the entanglement entropy have been proven for local Hamiltonians [50,[63][64][65].Here, we prove novel bounds for the negativity, ones that are tailored to be applicable to many-body localised systems and the geometry at hand.We now look at the specific setting of having two regions, A and B separated by a region C that has been traced out, and ask how much entanglement can be generated by a many-body localised Hamiltonian.We consider states evolving in time as ρ(t) := e −itH ρe itH .For the purpose of this proof, we will consider a slightly strengthened definition of quasi-locality with respect to the definition in the main text, namely, instead of the normalized Frobenius norm, we will assume the l-bits are localised in the operator norm, i.e., the following.

Definition 4 (Strong quasi-locality
).An operator is said to be quasi-local around a region R, if for any region X ⊂ R, it holds that for some constant K > 0.
This definition is stronger than the definition in the main text, in the sense that if an operator is quasi-local in the sense above, it is also quasi-local in the sense of the main text.For convenience, we repeat here the definition of a many-body localized Hamiltonian:

Definition 5 (Many-body localisation). A Hamiltonian
with real weights {ω (1) j } and {ω (2) j,k }, is called many-body localised if it can be written as a sum of mutually commuting ([h j , h k ] = 0 for all j, k) quasi-local terms h j , each centred around site j, and if ω i1,...,in ≤ ωe −|i1−in|/κ , where The following is the central statement from which the entanglement bound easily follows: Theorem 2 (Entanglement growth bound for sums of quasi-local operators).Let ρ be any initial state.Let H be a many-body localised Hamiltonian as per Definition 5 with localisation length ξ ≤ 1/(4 log(2)) and 2(1/κ − log(2)) > 1/ξ, consider three blocks A, C, B such that C divides A from B the growth of the negativity of the state ρ(t) = e −itH ρe itH restricted to the regions A, B for times t ≥ 0 and for any r ≥ |C|/2 Notice the the r in the above is a parameter one can choose freely.We obtain the following statement, morally, by picking r = 2ξ log(t) in the above bound if t > e |C|/(4ξ) and r = |C|/2 otherwise.Some additional details are required in order not to consider a time dependent r in the proof of Theorem 2.
Corollary 1. [Logarithmic growth of the negativity] If t ≥ e |C|/4ξ , under the assumptions of Theorem 2, we have while for t < e |C|/4ξ , Proof.For the t < e |C|/4ξ case, one need only pick r = |C|/2.Otherwise, for any time t, let r * be the positive integer such that e r * /(2ξ) ≤ t ≤ e (r * +1)/(2ξ) .Then r * ≤ 2ξ log(t) ≤ r * + 1. Picking r = r * we get since the bound with r = |C|/2 still holds, one can pick the minimum between these two bounds.
To turn this result into a proper bound for many-body localised Hamiltonians, we need to reduce such Hamiltonians to the form considered in Theorem 2. For this purpose, we will assume a slightly stronger definition of the many-body localised Hamiltonian.We call a unitary quasi-local with localisation length ξ if it maps any local operator on a region R to a quasi-local operator around R with localisation length ξ.We will assume that the Hamiltonian is diagonalised by a quasi-local unitary, which means that the l-bits are simply dressed Pauli Z operators, h i = U σ i z U † .In particular this implies that a product of l-bits For the purposes of the subsequent discussion, it will be useful to define the projector for any region X of the lattice.It is easy to verify that P X is a projector and that it is self-adjoint in the Hilbert Schmidt inner product.In what follows, ∥ • ∥ denotes the operator norm.
Lemma 1 (Quasi-local sums).Let H be a many-body localised Hamiltonian in the sense described above with localisation length ξ, and in addition, assume 2(1/κ − log(2)) > 1/ξ, where κ has been defined in Definition 5. Then the Hamiltonian can be written as where H i is quasi-local around the site i with localisation length 2ξ.
We note that, a quasi-local operator with localisation length ξ is quasi-local for any localisation length ξ ′ > ξ, if 2(1/κ − log(2)) ≤ 1/ξ, it suffices to relax the quasi-locality of the l-bits by increasing the localisation length until this condition is satisfied, i.e., choose for some constant ϵ > 0. This will result in faster (∼ r/ξ ′ ), but still exponentially slow, growth of entanglement in the bound.
Proof.Define H i as where I k = {i 1 , . . ., i k } are collections of k sites contained in i − l . . .i + l, containing the left boundary and either the right boundary or the site immediately to the left of it.We denote the set of all such I k as B k (i, l) (c.f. Figure S10).In the above, we have defined h I k = h i1 . . .h i k .We have then ω 2 I k ≤ ω 2 e −(2l−1)/κ and that h I k is quasi-local around [i − l, i + l] with localisation length ξ.We have H = i H i , it remains to show that the H i are quasi-local as claimed.Let R = [i − r, i + r] be a stretch of sites of radius r around i. We have, by the triangle inequality, . . .where we have used quasi-locality and in the l ≥ r part of the sum, that ∥h i − P R (h i )∥ ≤ 2. Now notice that

< l a t e x i t s h a 1 _ b a s e 6 4 = " 6 T q o t s f I 6 p N G 5 D h t 8 n A y K 1 W A p Y 4 = " >
we then have that ωe −2l(1/κ−log( 2)) e −(r−l)/(2ξ) + 2ω l≥r e −2l(1/κ−log(2)) . ( The first term is already bounded as needed, and for the third term, we use 2(1/κ − log( 2)) > 1/(2ξ), to get We note in passing that if δ = 0, the bound does not diverge, but an additional linear term is added to get a bound of the form re −rξ .In conclusion, with that is, which is the definition of quasi-locality noticing The rest of this section is dedicated to proving Theorem 2. The technique used for the proof is analogous to that of Ref. [50], where an analogous bound was derived for the entanglement entropy in the case |C| = 0, the properties of the negativity, which is a meaningful entanglement measure for mixed states, allow the generalisation to disconnected regions.The following two properties of the 1-norm and the partial transpose will be used repeatedly: • The 1−norm contracts under the partial trace, i.e., ||ρ A || 1 ≤ ||ρ A,B || 1 for any bipartite state ρ A,B [66], [67].
We will now prove Theorem 2. From now on we will assume that |C| is even, the odd case is analogous.We will label the two central sites of C as ±1, all the sites to the right of the center of C have positive integer labels, and to the left negative integer labels (notice that there is no site 0).In addition, we will use the standard notation [a, b) denoting real intervals to denote intervals in the chain, so from now on [a, b) is understood to mean [a, b) ∩ Z/{0}.
Divide the chain into the following regions (see Fig. S11), recall that r is an arbitrary integer greater than |C|/2.
and divide the terms in the Hamiltonian accordingly as where h S = i∈S h i for S ∈ L, ∂, R. We will now split the term h S into a main local part acting on S and a tail whose norm is exponentially small in r.More specifically, we can isolate the tail evolution as follows Lemma 2. Let H S = P S (h S ) and define the corresponding tail as T S = h S − H S , then e ih S t = e iH S t U T S (t) (40) where the tail evolution is given by where TS (τ ) = e −iH S τ T S e iH S τ .
Proof.We have U T S (t) = e −iH S t e ih S t , hence d dt U T S (t) = −iH S e −iH S t e ih S t + ie −iH S t h S e ih S t = ie −iH S t T e ih S t = TS (t)U T S (t) (43) this is the differential equation satisfied by the time evolution operator with the time-dependent Hamiltonian TS (t), it is well known that its solution is given by the Dyson series in the statement.
Notice that the tail evolution is an operator acting globally on the lattice, hence it could still, in principle, create an amount of entanglement extensive in the system size.The intuition is that it is generated by TS , which has exponentially small norm.We need to formalize the idea that this leftover tail evolution cannot create much entanglement.For each S, define S + q as the region S extended to the left and the right (if possible) by q sites.We can then write the tail T S as a telescopic sum, defining H S+q = P S+q (h S ) where ZS+q (t) = e −iH S t Z S+q e iH S t acts only on S + q.We then have, using Eq. ( 50) and the triangle inequality  (53) before plugging this in in the initial expression, notice that we then have We can conclude by using that by Lemma 3 that this expression has a norm bounded by ∥ ZS+q (t)∥ = ∥Z S+q ∥ ≤ O e −(r+q)/(2ξ) .
This allows us to prove the first part of the bound: We can write Given a generator h, define the evolution operator U h (ρ) := e ih ρe −ih , furthermore define the tail evolution operator as U T,t S (ρ) = U T (t) S ρ(U T (t) S ) † .Consider a state σ, notice that the overlap between B and ∂ is of length 2r − |C|, in particular, Here, we simply assume sufficient time has passed that H ∂ has had time to saturate the entanglement in the region ∂, that is, create a maximally entangled state between the sites on either side of C. We have then Next, by Lemma 4, we have ≤ exp tO e −r/(2ξ) ∞ q=0 2 2q e −q/(2ξ) ≤ exp tO e −r/(2ξ) , where we have used the assumption ξ < 1/(4 log(2)), so that the terms in the sum above are exponentially decaying and sum to a constant.Since L has no overlap with (2r − |C|/2, ∞), we can just bring the partial transpose inside the L evolution and eliminate it by unitarity.Afterwards, we eliminate U T,t L , exactly as we have eliminated U and in the same way one can show that E L,t ≤ exp t O e −r/(2ξ) .Finally, notice that for any region X and a state σ, ||σ T X || 1 = ||σ T X c || 1 , since σ is Hermitian.Then where we have used that the initial state is a product state, E R,t can be bounded as in the previous cases.Overall, we have then proven that ||ρ(t)

FIG. 2 .
FIG. 2. Behaviour of entanglement negativity in time and space: Results showing the growth of the negativity EN (r, t) with time for different distances r.Data is shown for a system size L = 24 and a disorder strength d = 8.0, averaged over Ns = 100 disorder realisations.a) The dynamics of EN (r, t) following a quench from a Nèel state, showing the logarithmic growth at late times.The circular markers are the raw data points, while the solid lines are a smoothed guide to the eye.The error bars indicate the standard error in the mean.We note that these error bars show agreement on average between the various disorder realizations, but they are not fully statistically independent errors, as would be expected in an experiment where each data point would come from a different run.b) The full dynamics of EN (r, t), reflecting the logarithmic 'light cone'.Each circle maps the point where the negativity grows beyond the corresponding threshold ε and the lines are linear fits.c) By extracting the behaviour of EN (r, t * ) ∝ exp(−r/ξ) at fixed times t * [dashed vertical lines in panel (a), horizontal lines in panel (b)], we can extract a well-defined length scale ξ(t), which depends only weakly on time.The solid lines indicate the fits to the data points which are used to extract the l-bit length scale, demonstrating convergence at late times.

FIG. 3 .
FIG. 3. Extracted l-bit length scale: The characteristic l-bit length scale ξ extracted from the entanglement negativity at time t * = 500, shown in blue for L = 24 with Ns ∈ [50, 100] disorder realisations and various values of the disorder strength d.Error bars indicate the fit error and are roughly the same size as the plot markers.Orange lines mark the localisation length obtained through exact diagonalisation following Ref.[31].For further details on calculating the localisation length using exact diagonalisation, see Supplementary Note 7[51].The black line indicates the localisation length of the corresponding Anderson-localised system, obtained by directly diagonalising the Hamiltonian in the non-interacting limit (∆ = 0), for a system of size L = 128 with Ns = 10000.

Fig
Fig. S1.Relative error of the simulation: A comparison of the relative error in the energy of the time-evolved state for different values of the disorder strength d, shown for system sizes L = 14 (χ = 128, averaged over Ns = 240 disorder realisations) and -at strong disorder only -L = 24 (χ = 192, averaged over Ns = 100 disorder realisations).The relative error remains below 1% for all disorder strengths.Error bars indicate the variance over disorder distributions and in most cases, are of comparable size to the plot markers.

Fig. S2 .
Fig.S2.Entanglement negativity at different disorder strengths: A comparison of the dynamics of the entanglement negativity EN (r, t) for different values of the disorder strength d, shown for L = 14 with bond dimension χ = 128 and averaged over Ns = 240 disorder realisations.In the delocalised phase, the negativity saturates to a value determined by the size of the subsystems A and B, while in the localised phase the negativity displays a slow ∝ log(t) growth even at late times.In this dephasing regime, we are able to use the data shown here to extract a length scale that characterises the localised phase, as detailed in the main text.

Fig
Fig. S3.Entanglement negativity for different bond dimensions and disorder strengths: The dynamics of the negativity EN (r, t) with r = 0, for different bond dimensions and disorder strengths.Data shown is for L = 24, averaged over Ns = 100 disorder realisations.Error bars show the standard error.

t
Fig.S4.Fits of the spatial decay of entanglement negativity for different times and disorder strengths: A comparison of the fits to the entanglement negativity EN (r, t * ) at different times t * and for different values of the disorder strength d, shown for L = 14 and averaged over Ns = 240 disorder realisations.The circular markers show the data, while the solid lines indicate the linear fits used to extract the l-bit localisation length.Error bars showing the standard error in the mean are typically smaller than the marker size.At large disorder strengths (i.e., in the localised phase), the linear fit is very good, confirming that the negativity does indeed decrease exponentially with distance in this phase.After a time of t * ≈ 100, the gradient of the decay (and hence the corresponding l-bit localisation length) does not strongly change with time in the localised phase.

t
Fig.S5.Fits of the spatial decay of entanglement negativity for different times and disorder strengths: The same as in Fig.S4, but for L = 24, χ = 192 and averaged over Ns = 100 disorder realisations, shown only in the localised phase.The solid lines indicate the range of points over which the exponential fits were performed in order to extract the length scale shown in Fig.3of the main text.

Fig. S6 .
Fig. S6.Extracted length scales at different times: A comparison of how the length scale ξ(t * ) changes as the measurement time t * is varied.The results here are extracted from the fits shown in Fig. S5.At short times, the value of ξ(t * ) changes rapidly, however, at longer times, the dependence of t * weakens significantly.The black line is the Anderson localisation length for comparison, as discussed in the main text.For clarity, error bars are not shown except on the Anderson localisation length.(Note that error bars are shown in the data presented in the main text.)

Fig. S8 .
Fig.S8.Growth of negativity between subsystems of fixed size: Entanglement negativity following a quench from a Néel state, shown for disorder strength d = 6.0 with system size L = 16, averaged over Ns = 96 disorder realisations.Here we compute the negativity for a subsystem of fixed size R b := |A| = |B|, as sketched in Fig.S7.We can see that for small values of R b , the negativity at small separations r saturates quickly and that only for larger values of R b does the negativity increase in a manner that allows extraction of the relevant l-bit length scale.Error bars show the standard error in the mean, and the solid lines are a smoothed guide to the eye.

Fig. S9 .
Fig.S9.Mutual information between separated regions over time: A comparison of mutual information (solid lines) and entanglement negativity (dashed lines) between subsystems A and B separated by a distance r, for a system of size L = 12 with bond dimension χ = 128, averaged over Ns = 100 disorder realisations.
Fig. S10.Examples of sets of sites contained in B k (i, l): In this example, k = 4 and l = 3.The sites i, i + l, i − l are highlighted in red.

@
< l a t e x i t s h a 1 _ b a s e 6 4 = " 8 G k L H F t Z d P E N P u L x e t D I b 2 + W B P U = " > A A A C D X i c d V D L S g N B E J y N r x h f U Y 9 e B o P g a Z l N o i a 3 g B e P E Y w J J C H 0 T j o 6 O P t g p l e Q x W 8 Q v O p v e B O v f o N / 4 S e 4 G y O o a J 2 K q u 6 p n v J j r S w J 8 e Y U 5 u Y X F p e K y 6 W g 7 V 9 T 8 Q 0 T P P X p M b b 0 i C x G I O 8 g g v s Z z S E A O 0 w n d 5 7 y / c S C x T x G A 1 X m k 9 F / L 6 R Q m D t T e B n k w H Q p f 3 t 5 e K f n u 9 H e v w r n C a N Y a r C O C E M Z Z 5 N S u M 0 2 0 q j s m 6 Q j 5 V B I s g / g 1 y F X I I B I j S K g 5 S Z m G R l l b K K v n r g / 5 P z q u v V 3 O p p v d K q z 8 o q s h 2 2 y / a Z x 4 5 Y i 5 2 w N u s w y T S 7 Z w / s 0 b l z n p x n 5 + V z t O D M d r b Z D z i v H 1 L D n M E = < /l a t e x i t > @ < l a t e x i t s h a 1 _ b a s e 6 4 = " e X P 0 O y 5 a R W a t e m Z S 2 x H b b L 9 p n L D l i b H b F j 1 m G S 3 b B 7 9 s A e n T v n y X l 2 X j 5 H Z 5 z p z j b 7 A e f 1 A + C o n q Q = < / l a t e x i t > L < l a t e x i t s h a 1 _ b a s e 6 4 = " y F d R 2 B h m W J L x t j 4 D h 4 S a X b 0 W K i o = " > A A A C B n i c d V D L S g N B E J y N r x h f U Y 9 e B o P g a d n d x C T e A l 4 8 e E j A P C B Z w u y k E w d n H 8 z 0 y B P b I n 6 9 5 6 t l 6 s 1 6 / R l D X f 2 W e / Y L 1 9 A u A b m 1 4 = < / l a t e x i t > r < l a t e x i t s h a 1 _ b a s e 6 4 = " i J U l Z t A i s u f x 1 e 9 s o t Z M m O U 8 T z M = " > A A A C B n i c b V D L S g N B E J y N r x h f U Y 9 e B o P g K e x G Q Y 8 B L x 4 T M A 9 I l t A 7 6 c Q h s w 9 m e o W w 5 C 5 4 1 d / w J l 7 9 D f / C T 3 B 3 3 Y M m 1 q m o 6 q a 6 y 4 u U N G T b n 1 Z p b X 1 j c 6 u 8 X d n Z 3 d s / q B 4 e d U 0 Y a 4 E d E a p Q 9 z 0 w q G S A H Z K k s B 9 p B N 9 T 2 P N m N 5 n f e 0 B t Z B j c 0 T x C 1 4 d p I C d S A K 2 y s 6 Z w 6 5 Y k 9 2 y F u s w w Z A 9 s W f 2 Y j 1 a r 9 a b 9 f 4z W r K K n W P 2 B 9 b H N z B K m V k = < / l a t e x i t > r < l a t e x i t s h a 1 _ b a s e 6 4 = " i J U l Z t A i s u f x 1 e 9 s o t Z M m O U 8 T z M = " > A A A C B n i c b V D L S g N B E J y N r x h f U Y9 e B o P g K e x G Q Y 8 B L x 4 T M A 9 I l t A 7 6 c Q h s w 9 m e o W w 5 C 5 4 1 d / w J l 7 9 D f / C T 3 B 3 3 Y M m 1 q m o 6 q a 6 y 4 u U N G T b n 1 Z p b X 1 j c 6 u 8 X d n Z 3 d s / q B 4 e d U 0 Y a 4 E d E a p Q 9 z 0 w q G S A H Z K k s B 9 p B N 9 T 2 P N m N 5 n f e 0 B t Z B j c 0 T x C 1 4 d p I C d S A K Fig. S11.Division of the chain into regions for the proof of Theorem 2: The central region bounded by the vertical lines is C, represented by the dashed red and blue line.Note that C is contained in both L and R.