Loschmidt echo singularities as dynamical signatures of strongly localized phases

Quantum localization (single-body or many-body) comes with the emergence of local conserved quantities -- whose conservation is precisely at the heart of the absence of transport through the system. In the case of fermionic systems and $S=1/2$ spin models, such conserved quantities take the form of effective two-level systems, called $l$-bits. While their existence is the defining feature of localized phases, their direct experimental observation remains elusive. Here we show that strongly localized $l$-bits bear a dramatic universal signature, accessible to state-of-the-art quantum simulators, in the form of periodic cusp singularities in the Loschmidt echo following a quantum quench from a N\'eel/charge-density-wave state. Such singularities are perfectly captured by a simple model of Rabi oscillations of an ensemble of independent two-level systems, which also reproduces the short-time behavior of the entanglement entropy and the imbalance dynamics. In the case of interacting localized phases, the dynamics at longer times shows a crossover in the decay of the Loschmidt echo singularities, offering an experimentally accessible signature of the interactions between $l$-bits.

Introduction. Localization phenomena in systems of quantum particles offer striking evidence of the role of interference in quantum dynamics. Constructive interference of paths bringing a particle back to its initial location in real space is at the heart of single-particle (or Anderson) localization (AL) [1,2]; more recently a similar phenomenon occurring in Hilbert space (dubbed many-body localization -MBL) [3][4][5][6] has been shown to prevent many-body quantum systems from relaxing to thermal equilibrium, undermining the ergodic hypothesis in a large class of models of interacting quantum particles. Localized phases are generally characterized in the negative (absence of transport, of long-range order, of spectral gaps, etc.), while positive characterizations are generally elusive. A crucial aspect of localization is the persistence of initial conditions, which, in the case of AL of non-interacting particles, is related to the conservation of populations in the localized single-particle eigenstates of the Hamiltonian. In the case of MBL, the analog phenomenon would be the appearance of local conserved quantities (called local integrals of motion or l-bits [7][8][9][10]) which are obtained by unitary transformations of local operators; and which, if extensive in number, constrain the dynamics of the system to the point of preventing relaxation. The existence of l-bits in disordered spin chains can be mathematically proven under the assumption of limited level attraction [9], and approximate l-bits for many-body systems can be constructed with a variety of analytical as well as numerical methods [11][12][13][14][15][16][17][18][19][20]. Much of the phenomenology of MBL dynamics (persistence of traits of the initial state, logarithmic growth of entanglement entropies, etc.) -observed in numerical studies as well as in experiments [4,5] -can be directly explained in terms of the existence of l-bits and interactions between them. Yet observing l-bits directly is an arduous task, given that their expression is highly disorder-dependent (and generally unknown even in theory), and it would require high-precision measurements of local observables in different local bases. An even more arduous task is the one of probing directly the existence of interactions among l-bits, which is a defining feature distinguishing MBL from AL.
The purpose of this Letter is to show that, in the case of strongly localized phases, the existence of l-bits can offer striking signatures in the dynamics of the Loschmidt echo (LE), namely in the logarithm of the return probability to the initial state |ψ 0 Here H is the system's Hamiltonian and L its size; [...] av indicates the disorder average. When |ψ 0 has a simple factorized form and in the case of strong disorder, we find that the LE displays periodic singularities, decaying very slowly in amplitude -as illustrated using a model of disordered spinless fermions in 1d (corresponding to the S = 1/2 XXZ model in a fully random or quasiperiodic field) initialized in a charge density-wave (CDW) state. The singularities in the LE are fully explained quantitatively by a simple model of a collection of localized 2-level systems (2LS) undergoing independent Rabi oscillations, and approximating strongly localized l-bits. The same minimal model captures quantitatively the dynamics of the entanglement entropy at short times as well as of the number entropy at longer times; and the dynamics of the density imbalance characterizing the initial state. At longer times the deviation of the exact results for the MBL dynamics from the predictions of the 2LS ensemble offers direct evidence of the interactions among the l-bits in the form of a faster decay of the LE singularities and imbalance oscillations. As the LE is generally accessible to quantum simulators measuring individual degrees of freedom [21,22], our results show that strong signatures of l-bits dynamics are within the immediate reach of state-of-the-art experiments on disordered quantum systems.
Model. Our platform for the investigation of LE dynamics is given by a paradigmatic model, namely the S = 1/2 XXZ chain in an inhomogeneous magnetic field [23], corresponding to a model of spinless fermions with nearest-neighbor interactions in an inhomogeneous local chemical potential where S α i (α = x, y, z) are spin operators and c i , c † i and n i = c † i c i are fermionic operators; the equality between the two Hamiltonians is true up to an additive constant [24]. In the following the external field/potential h i is taken to be either quasi-periodic (QP), namely h i = ∆ cos(2πκi + φ) with κ = 0.721 (inspired by experiments on bichromatic optical lattices [25,26]) and φ a random phase; or to be fully random (FR) and uniformly distributed in the interval [−∆, ∆]. We consider chains of length L (up to L = 22) with open boundaries, and we average our results over ∼ 10 3 realizations of the random phase (QP) or of the full random potential (FR). All the unitary evolutions considered in this study are obtained using exact diagonalization (ED), and they start from the charge-density wave state |ψ 0 = |1010101... , corresponding to a Néel state for the spins.
We shall focus on the case of interacting fermions J z = J (corresponding to an SU(2) invariant spin-spin interaction) and contrast it with the limit of free fermions J z = 0. In the latter case, the QP potential leads to a transition to fully localized single-particle eigenstates for ∆ ≥ J, with an energy independent localization length ξ = 1/ log(∆/J); while the FR potential leads to AL of the whole spectrum at any infinitesimal value of disorder. In the interacting case, instead, a QP potential of strength ∆ 4J [27] and a FR potential of strength ∆ 3.5J [28] are numerically found to lead to MBL, although the exact critical value ∆ c has been recently the object of further investigation [28][29][30][31][32][33][34][35]. In the following we will conduct our discussion starting from the case of the QP potential, which has a simpler spatial structure devoid of rare regions, leading to stronger localization effects; and we shall later discuss how the picture should be enriched in the case of the FR potential to account for the existence of rare regions.
LE singularities. Fig. 1 shows the dynamics of the LE, λ(t), along with that of the imbalance I(t) = i (−1) i (2[ n i ] av − 1)/L -the latter saturates to its maximum value of 1 in the initial state and probes the persistence of the initial density/spin pattern [36]. We observe that for both the QP and FR potentials, and for disorder strengths compatible with the onset of the MBL regime, the LE displays a sequence of periodic cusp-like peaks at times t n = (2n + 1)π/J (n = 0, 1, 2, ...). These times correspond to minima in the imbalance, as the system reaches instantaneous configurations which are the farthest from the initial spin/density pattern. A closer inspection shows that, for sufficiently strong disorder, all the peaks become sharp cusps, namely they represent genuine non-analyticities of the LE. They are rather remarkable given that they survive disorder averaging, and they seemingly appear in a finite fraction of disorder realizations (see [37] for further details); and in particular they decay very slowly in time, as we shall discuss in detail in the following.
2LS model and l-bits. All the essential details of the short-time evolution of the LE can be captured with a surprisingly simple, yet rather insightful model. This model is best understood (and justified) in the case of the QP potential, as illustrated in Fig. 2. In the case of strong disorder, the fastest dynamics in the system starting from a Fock state will be offered by those particles that sit on a site i which is nearly resonant with its unoccupied neighbor (say i + 1), because the hopping J/2 is either larger than the energy offset δ i = h i+1 − h i (in the non-interacting case) or larger than the screened offset δ i − J z (in the presence of nearest-neighbor repulsion). These 2-site clusters, representing nearly resonant two-level systems (2LS), have the property of being spatially isolated in the QP potential, because of the strong  (1)  (1) presents a pair of quasi-resonant sites for the particle in orange; in the case of interacting particles, region (2) shows two quasi-resonant sites for the orange particle, thanks to the partial screening of disorder offered by the interaction with the red particle.
anticorrelation among two consecutive energy offsets (δ i and δ i±1 -see [37]). As a consequence, a nearly resonant 2-site system will be generally surrounded by highly nonresonant pairs of sites, which can be considered as nearly frozen to the initial state. This invites us to write for the evolved state the following 2LS Ansatz where the first tensor product ⊗ n runs over the nearly resonant 2LS, while the second tensor product ⊗ i runs over the leftover sites (we have taken the freedom of reordering the sites arbitrarily in the tensor product). |ψ (n) 2LS (t) is the evolved state of the n-th (isolated) 2LS system, corresponding to two states split by an energy difference δ n = δ n − J z and connected by a Rabi coupling J [37]; while |ψ 0,i is the (persistent) initial state of the site i belonging to the remainder of the system.
The LE for such a system is readily calculated as with p(δ, Ω, t) = (Ω/Ω ) 2 sin 2 (Ω t/2) (and Ω = √ Ω 2 + δ 2 ) the well-known probability of finding the 2LS in the state orthogonal to the initial one while performing Rabi oscillations [38]. When averaging Eq. (4) over disorder, it is immediate to obtain the following simple expression Here P (x) is the probability that the energy offset between two neighboring sites takes the value x [39]. Eq. (5) is an analytical integral formula which depends uniquely on the (known) statistics of the disorder potential via the P distribution. In the case of the QP potential [40]; while for the FR potential P (x) is the normalized triangular distribution defined on the [−2∆, 2∆] interval. Fig. 3(a-c) shows that, for the case of the QP potential, Eq. (5) is able to predict with high accuracy the ED results deep in the MBL phase without any adjustable parameter. In particular the cusp singularities of the ED results are easily explained as descending from the divergent singularity of the LE for a fully resonant 2LS with Ω = Ω = J, reaching a state orthogonal to the initial one after odd multiples of half a Rabi oscillation t n = (2n + 1)π/Ω. These divergences are smoothened into cusp singularities due to the fact that such resonant 2LSs are a set of zero measure in the disorder statistics. This result has important consequences. Indeed the nearly resonant 2LSs captured by the model are clearly an ensemble of approximate l-bits with Hamiltonian where , built from the original spin operators for the pair n = (i, i + 1); and K n = δ 2 n + J 2 is the l-bit splitting. Hence the LE singularities are a striking manifestation of the existence of such (nearly free) l-bits, to be found in the short-time dynamics of the system.
It is worthwhile to mention at this point that the existence of singularities in the quench dynamics of the LE is currently the subject of several theoretical and experimental investigations, as they represent the main signature of so-called dynamical quantum phase transitions, studied both in non-random systems [22,[41][42][43][44] as well  as in disordered quantum systems [45][46][47]. Nonetheless our observation of LE singularities is fully explained by a model of individual 2LS, without the need of any manybody effect; therefore we shall refrain from associating them to any form of time-dependent transition. From 2LS to 3LS. Fig. 3(d-f) shows that, in the case of the FR potential, the 2LS model of Eq. (5) still predicts the correct frequency of the LE singularities, but not the correct height and it also misses a global offset. This is not surprising, as in the case of the FR potential the assumption of anticorrelation between the energy offset of contiguous pairs is no longer valid, namely the potential can host "rare" regions in which contiguous pairs of sites -(i − 1, i) and (i, i + 1) -are nearly resonant at the same time. To take those regions into account (at least partially) one can easily promote the 2LS model to a model of 3-site systems (amounting to effective threelevel systems -3LS), and approximate the evolved state as that of a collection of independent 3LS. Unlike the case of the 2LS model, the 3LS model requires a fully numerical treatment (detailed in [37]), which amounts to numerically scanning the ensemble of 3-site systems in a single realization of the disorder potential in an arbitrarily large chain. As shown by Fig. 3(d-f), the improvement offered by the 3LS model for the FR potential is substantial; these results can further be improved by moving to 4-site clusters etc. albeit at an exponential cost.
Slow dephasing and l-bit interactions. A significant feature of the LE singularities is their slow decay in time -which is remarkable given that they result from the Rabi oscillations of a collection of 2LS with a distribution of frequencies that can be a priori expected to lead to fast dephasing. The reason behind the slow decay is also captured by the 2LS model, Eq. (5), namely by the fact that the integral expressing LE takes contributions from a small window of detunings δ around zero, the smaller the longer the time. When looking at the singularity times t = t n , a direct inspection of the function log(1 − p(δ , Ω, t n )) seen as a function of δ shows that it has a large peak centered on δ = 0 with a width depending on time as t −1/2 n [37]. The singularity in the average LE comes from the integral of this peak, while the rest of the integral contributes essentially to the regular part of the LE; hence it is immediate to predict that the height of the cusp singularity should decay as the peak width, (namely as t −1/2 n ). Fig. 4(a) shows the time evolution of singularity peaks in the LE for free as well as interacting fermions in the QP potential, compared to the prediction of the 2LS model (for the interacting case): we observe that the t −1/2 decay is indeed confirmed by the ED data for free fermions, as well as by the ED data for interacting fermions at sufficiently short times (tJ t * ≈ 100 for ∆ = 8J). On the other hand, at longer times the interacting data are found to display a strong deviation from the 2LS model prediction, exhibiting a much faster decay. This crossover to an interaction-induced dephasing (IID) regime clearly shows the limits of the 2LS model as a model of free l-bits expressed by Eq. (6), and it marks a fundamental difference between AL and MBL in the QP system. Indeed the faster decay of the LE must be related to the effect of l-bit interactions, which are a defining feature of MBL, and which add terms of the kind nm U nm τ n τ m + nml V nml τ n τ m τ l + ... to the effective l-bit Hamiltonian. Such terms are responsible for the persistent growth of entanglement entropy in the system [7] as the logarithm of time, and indeed the onset of the log t growth of entanglement occurs at a time compatible with t * (see [37]). A similar crossover from a slow power-law decay of the LE peak height to a faster decay, dictated by the presence of interactions, is also exhibited by the comparison between the ED data for interacting fermions in the FR potential with the same data for noninteracting fermions and for the 3LS model -as shown in Fig. 4(b).
Further predictions. The 2LS and 3LS models allow us not only to predict the dynamics of the LE, but also that of the imbalance and of the entanglement entropy -which in the case of the 2LS model acquire an explicit integral formula similar to Eq. (5) for the LE dynamics. The results (detailed in [37]) for the imbalance dynamics show that the oscillations exhibited by the imbalance (see Fig. 1(a,c)) decay in amplitude as t −1/2 similarly to what is seen for the LE singular peaks, crossing over to a faster decay due to the l-bit interactions in the case of interacting fermions. The entanglement entropy is also correctly captured by the 2LS (3LS) models (as the average entropy across a bipartition of the 2-site (3-site) cluster) at all times for free fermions, and at short times for interacting ones. In the case of interacting fermions at longer times, the models in question only capture the number entropy component of the entanglement entropy [48,49], something which is readily understood in that the residual entropy component (the so-called configurational entropy) precisely originates from the l-bit interactions.
Conclusions. In this work we have shown that sharp cusp-like singularities in the Loschmidt echo (LE) are a generic feature of the localized dynamics of an extended quantum system initialized in a factorized state. These features can be fully explained by the dynamics of a simple model, describing an ensemble of effective independent two-level (or even three-level) systems, offering an explicit approximation to the conserved l-bits in the MBL regime. Such a model predicts very accurately the LE singularities for strongly disordered systems as well as their decay; the faster decay in the dynamics compared to that predicted by the model is a direct manifestation of the interactions between the l-bits, namely the defining feature of many-body localization (MBL) with respect to Anderson localization (AL). Based on our results, we can conclude that experimental evidence of l-bit dynamics and of their interactions is readily accessible to state-of-the-art quantum simulators which have direct access to the Loschmidt echo, such as e.g. trapped ions [21,22], quantum-gas microscopes [48] or superconducting circuits [50].
Acknowledgements. L.B. gratefully acknowledges hospitality and financial support from the Laboratoire de Physique of the ENS Lyon. R.A.R. and L.B. acknowledge funding from the CY Initiative of Excellence (grant"Investissements d'Avenir" ANR-16-IDEX-0008) where this work developed during R.A.R.'s stay at the CY Advanced Studies. Part of the exact diagonalization simulations were performed using routines contained in the QuSpin [51,52] and quimb [53]

SUPPLEMENTAL MATERIAL
Loschmidt echo singularities as dynamical signatures of strongly localized phases

I. TWO-SITE SYSTEM AS TWO-LEVEL SYSTEM AND ITS RABI OSCILLATIONS
Let us isolate a two-site system (i, i + 1) hosting one particle in the fermionic chain, with Hamiltonian where we assume that the site i + 2 is occupied by a (pinned) particle, while size i − 1 is empty (or occupied by a pinned hole). Introducing the spin operators the Hamiltonian becomes simply namely, a two-level system (2LS) with splitting δ and Rabi frequency J. If the system starts from the |10 state, the return probability is given by the well-known formula for the probability of persistence in the initial state during Rabi oscillations [38], namely 1 − p(δ; t), where p(δ; t) = 1 1 + (δ/J) 2 sin 2 1 + (δ/J) 2 2 tJ . (S4)

II. SPATIAL CORRELATIONS IN THE QUASI-PERIODIC VS. FULLY RANDOM POTENTIAL
A fundamental assumption of the 2LS model described in the main text is that quasi-resonant two-site systems are spatially isolated in a (quasi)-disordered chain -namely, if a pair of sites (i, i + 1) is quasi-resonant for the motion of a particle, the two adjacent pairs of sites (i − 1, i) and (i + 1, i + 2) are not resonant. Defining δ 1 = h i+1 − h i as the energy difference of the two sites in question, and δ 2 = h i+2 − h i+1 as that of the following pairs of sites, in the case of non-interacting fermions, the above condition requires that the two energy differences do not vanish simultaneously.
Such a form of correlation is indeed observed in the case of the quasiperiodic (QP) potential: Fig. S1(a) shows the joint probability P (δ 1 , δ 2 ) for two adjacent energy differences, displaying a dip for δ 1 = δ 2 = 0 -an aspect which prevents two successive pairs of sites from being resonant simultaneously. In the case of interacting fermions, on the other hand, the above condition requires that if, e.g. J z ± δ 1 ≈ 0, then J z ∓ δ 2 is non-zero, or vice versa -this prevents a state of the type |101 on the sites (i, i+1, i+2) from being simultaneously (quasi-)resonant with |110 and |011 , or, similarly, the state |010 from being quasi-resonant with |100 and |001 . This is indeed guaranteed by the fact that P (δ, −δ) is nearly vanishing for any finite δ, except for δ ≈ 1.25∆ -but the latter situation does not lead to consecutive resonances when ∆ > 1.25J, which is always the case in our study.
On the other hand the uniform potential has no correlations between two consecutive energy differences, and the P (δ 1 , δ 2 ) distribution is the product of two triangular distributions for δ 1 and δ 2 -shown in Fig. S1(b). This implies that a fundamental assumption behind the 2LS model description is not guaranteed to be satisfied -while it is more likely to have two adjacent pairs of sites with different energy offsets than with similar ones, one cannot exclude the existence of "rare" regions with consecutive nearly resonant pairs. This requires to improve the 2LS model to a three-site (three-level) one, as detailed in Sec. IV.°2°1

III. DECAY OF THE LOSCHMIDT-ECHO SINGULARITIES AND IMBALANCE OSCILLATIONS FROM THE 2LS PREDICTION
As seen in the main text, the 2LS predicts the LE in the form of the integral and p(δ; t) as given in Eq. (S4). The cusp singularities in λ(t) at times t n J = (2n+1)π, n = 0, 1, 2, . . ., descend from the fact that the function f (δ; t n ), seen as a function of δ, develops a logarithmic singularity at δ = 0, as shown in Fig. S2, while it is fully regular at any other time. The singular peak centered at δ = 0 has a support shrinking with t n as t −1/2 n -as seen in Fig. S2 when plotting the function f (δ; t n ) vs δ √ t n , which leads to a collapse of the peak widths at different times (when n 1). The integral of the f function outside the peak contributes to the regular part of the LE, while the integral of the peak dictates fundamentally the height of the cusps above the regular background (estimated as the long-time averageλ), namely the quantity λ P (t n ) = λ(t n ) −λ. Therefore one can expect that for n 1, the height of the cusps in the LE decay with time as λ P (t n ) ∼ t −1/2 n . This prediction is confirmed in Fig. 4 of the main text.
The 2LS model prediction for the imbalance is very similar to that of the LE, as the imbalance is simply related to the persistence probability of the initial state (|10 or |01 ) on the 2-site cluster -given that the orthogonal state contributes zero to the imbalance. Therefore the 2LS expression for the imbalance simply reads The times t n giving cusp singularities in the LE correspond to dips in the imbalance, and these dips come from local dips in the g(δ; t) = 1 − p(δ, t) function centered around δ = 0 and touching zero for t = t n . The width of these dips is also shrinking in time as t −1/2 n . Therefore one expects the depth of the minima in the fluctuations of the imbalance to decay to the long-time average as t −1/2 n as well -this prediction will be verified in Sec. V.

IV. THREE-SITE MODEL
The Hamiltonian of a three-site system (i, i + 1, i + 2) containing two interacting fermions in an initial |101 state is explicitly given by Here all single-site energies are referred to the energy of site i, and δ i = h i+1 − h i . The above Hamiltonian assumes that the sites i−1 and i+3 remain empty during the time evolution. The Hilbert space of the 3-site system is restricted to the three states |101 , |110 and |011 , making of it a three-level system (3LS), with a generic time-dependent wavefunction Its explicit form can be easily calculated numerically for any specific choice of the energy differences δ i .
We can then model a chain in a QP or FR potential as an ensemble of independent 3LSs by generating sequences of energy offsets δ i , δ i+1 between adjacent site pairs according to the distribution P (δ i , δ i+1 ). The 3LS prediction for the Loschmidt echo of the ensemble is λ 3LS (t) (S12) In practice, the above integral can be sampled numerically by simply averaging over a large number of different realizations of the potential on 3-site systems, such as those offered by a very long chain, namely where α i = 101 if i is odd and 010 if i is even, and L 1. Similarly the imbalance can be calculated as with α i , β i , γ i = α(t), β(t), γ(t) orα(t),β(t),γ(t) depending on whether i is odd or even. The above expression for the LE (Eq. (S13)) has the apparent drawback of triple-counting each site. Nonetheless, similarly to what is done in the main text for the 2LS case, it is fair to assume (and it can be numerically tested) that, out of the three clusters containing each site, only one at most will contribute significantly to the LE. As a consequence the triple counting has a mild effect on the final result. One could avoid triple counting by thoughtfully decomposing a chain into non-overlapping clusters of up to 3 sites, in such a way as to maximize the LE; yet this procedure introduces significant complications which are not justified a posteriori, given the quality of the results offered already by the naive ensemble average (see Fig. 3 of the main text and further results in Sec. VI). Fig. S3 shows the comparison between the ED results for the imbalance dynamics of interacting fermions immersed in a QP and fully random potentials of variable strength, compared with the predictions of the 2LS and the 3LS models. For sufficiently strong disorder (∆ 6J) the 2LS predictions are already rather accurate in the case of the QP potential, and the 3LS model offers further improvement. On the other hand in the case of the FR potential the 3LS model offers a more substantial improvement, fixing an overall offset (for sufficiently strong disorder) which is seen in the 2LS predictions. Fig. S4 shows the evolution of the depth of the minima of the imbalance at times t = t n , taken with respect to the longtime average, namely the quantity I M (t n ) =Ī − I(t n ). We observe that the predictions of the 2LS system for the fermionic chain immersed in the QP potential shows a clear power-law decay at long times, compatible with t −1/2 , which is indeed reproduced in the case of noninteracting fermions. In the case of interacting fermions,  on the other hand, a crossover is observed at long times (t t * ≈ 100) to a faster decay (interaction-induced dephasing -IID -regime); this crossover is compatible with what is observed in the decay of the LE (see main text) as well as with the evolution of the entanglement entropy (see below -Sec. VI). As discussed in the main text, this crossover is to be attributed to the interactions between the l-bits.

V. IMBALANCE DYNAMICS
A similar picture is offered by the case of the FR potential. There the ED results are compared with the predictions from the 3LS model; the latter model predicts correctly the decay of the minima depth in the noninteracting case at all times, while the exact results for the interacting system show a clear crossover towards a faster decay for times t t * ≈ 50; this crossover to an IID regime is again fully compatible with that observed in the decay of the cusp maxima of the LE (see main text), as well as with the onset of the logarithmic growth in the entanglement entropy (see Sec. VI).

VI. ENTANGLEMENT DYNAMICS: FULL ENTROPY VS. NUMBER ENTROPY
A. Entanglement entropy from the 2LS and 3LS model The 2LS and 3LS models allow for a simple calculation of the entanglement entropy of a A/B bipartition of the system into two adjacent chains, defined as the von Neumann entropy of the reduced density matrix where ρ A (t) = Tr B |ψ(t) ψ(t)| is the partial trace (over the degrees of freedom in B) of the instantaneous purestate density matrix associated with the evolved state |ψ(t) . Within those models, the entanglement associated with such a bipartition simply comes from the entanglement inside the 2-site or 3-site cluster which contains the cut defining the bipartition. In the case of the 2LS model the disorder-averaged entanglement entropy of a bipartition is simply predicted as the entropy of the reduced state of one site in the 2-site cluster, namely and . Notice that, unlike for the formulas of the LE and of the imbalance, no-double counting is implied in the above formula, since the entanglement is referred to a cut of the chain, and there is one unique cut per 2-site cluster.
On the other hand a 3-site cluster can be cut in two different ways, that we will indicate as • | • • and • • | • in the following (where • stands for a site and | stands for the cut). The reduced density matrices for the two cuts are readily obtained from the cluster wavefunctions described in Sec. IV; e.g. for a 101 cluster the reduced density matrix associated with the • | • • cut reads with associated entanglement entropy (S18) the one associated with the • • | • cut reads The density matrices ρ •|•• 010 and ρ ••|• 010 and related entropies associated with a 010 cluster can be calculated similarly.
The disorder-averaged entanglement entropy of the whole system within the 3LS model is then given by where The factor 1/2 in (S20) comes from the double counting of each cut (which is contained both in a 101 cluster as well as in a 010 cluster).

B. Entanglement entropy vs. number entropy
The 2LS and 3LS models picture the entanglement between two adjacent subsystems as arising uniquely from the coherent motion of particles within the restricted size of the clusters they describe. When starting from a factorized state, this picture is certainly valid at short times. At long times it remains valid only if particles remain localized within the size of the clusters (namely if the localization length is smaller than the cluster size), and if this is a sufficient condition for entanglement not to spread any further. The latter aspect is true in the case of non-interacting fermions, for which the only mechanism behind entanglement of different spatial partitions is particle motion between them. On the other hand, in the case of interacting fermions in the MBL regime, entanglement keeps growing due to the interactions between l-bits, and distant degrees of freedom can become entangled even without any net particle exchange.
In this context it is useful to decompose the entanglement entropy of a subsystem A into a number entropy contribution, and a remainder part (called the configurational entropy), S A = S A,N + S A,c [48,49]. The number entropy is given by where p N A is the probability of having N A particles in subsystem A, and it accounts for the particle number uncertainty appearing in subsystem A because of the coherent exchange of particles with its complement B. On the other hand the configurational entropy accounts for correlations establishing between the particle arrangements in A and B once the partitioning of the particles between A and B has been fixed. The 2LS and 3LS models, lacking completely any form of correlations among the clusters, can only capture the number entropy contribution in systems with a localization length smaller than the cluster size; as we shall see in the next section, this Jt FR, J z = 1 limited picture still offers a faithful description of entanglement in the MBL regime for short times (the longer the stronger disorder is), while it can describe entanglement at all times for strongly localized non-interacting particles.
C. 2LS/3LS entropy vs. exact entanglement and number entropy Fig. S5 shows a comparison between the entanglement entropy of interacting fermions in a QP potential and the 2LS prediction. We observe that at moderate disorder in the MBL phase (∆ = 8J) the 2LS and 3LS models only capture the initial rise of the entanglement entropy and (partly) the first maximum; in particular the very existence of a maximum is explained by the models as the result of nearly resonant small clusters returning close to the initially factorized state -albeit at different times due to the inhomogeneously broadened local frequencies, which explains why the entanglement entropy does not come back to (nearly) zero. The 2LS and 3LS models on the other hand completely miss the long-time logarithmic growth of the entanglement entropy -something which is fully expected, given that such a growth is the consequence of interactions between l-bits, not included in the 2LS and 3LS models by construction. On the other hand, at stronger disorder (∆ = 20J) the interactions between l-bits are parametrically suppressed, and the 2LS and 3LS description of entanglement becomes accurate up to very long times.
As suggested in the previous section, a more appropriate comparison with the entanglement entropies of the 2LS and 3LS models would involve the number entropy from the ED data -shown in Fig. S6. For non-interacting fermions in a QP potential of strength ∆ = 8J the number entropy is found to nearly coincide with the full entanglement entropy, and to be very well described by the 2LS prediction -see Fig. S6(a). When adding the interactions, the agreement between the number entropy and the 2LS entropy deteriorates mostly at long times, seemingly due to the ∼ log log t growth of the number entropy observed in the MBL phase [49]. Similar considerations can be made in the case of the FR potential -the 3LS models describe well the entanglement and number entropy in the non-interacting case, and they miss the slow long-time growth of the number entropy in the interacting case.

VII. LOSCHMIDT-ECHO DYNAMICS FOR DIFFERENT DISORDER REALIZATIONS
Figs. S7 and S8 show the disorder average of the LE for a chain of L = 22 sites, along with all the disorder realizations (> 10 3 ) contributing the average, for various strengths (∆/J = 1, 2, ..., 10) of the QP and FR potential, respectively. We observe that sharp cusp singularities are exhibited by a signification portion of the disorder realizations, and that for sufficiently strong disorder these realizations are a finite fraction of the disorder statistics (in the asymptotic limit), so that cusp singularities persist in the disorder-averaged results as well. These plots also suggest the fact that cusp singularities can be observed with a limited disorder statistics, under realistic experimental conditions.