Converting long-range entanglement into mixture: tensor-network approach to local equilibration

In the out-of-equilibrium evolution induced by a quench, fast degrees of freedom generate long-range entanglement that is hard to encode with standard tensor networks. However, local observables only sense such long-range correlations through their contribution to the reduced local state as a mixture. We present a tensor network method that identifies such long-range entanglement and efficiently transforms it into mixture, much easier to represent. In this way, we obtain an effective description of the time-evolved state as a density matrix that captures the long-time behavior of local operators with finite computational resources.

Reconciling the time-reversal invariant unitary evolution of closed quantum many-body systems with the emergence of statistical mechanics and its well defined arrow of time is still an open question, hindered by the exponential complexity of simulating this problem in the generic case.
In contrast to the dynamical scenario, the equilibrium wave functions of a large class of physically relevant systems live in a small corner of the exponentially large Hilbert space.Such corner is characterized by a bounded amount of entanglement [1] and states therein admit an efficient approximation as a tensor network (TN) [2].This property implies we can perform very precise numerical simulations with only polynomial resources [3][4][5][6][7][8][9].The most notable example of successful TN simulations are those in one-dimension performed with matrix product states (MPS) [4,[10][11][12][13].
Out of equilibrium, instead, initially localized correlations can propagate over arbitrarily large distances.As a result, the entanglement in the system increases rapidly with time [14][15][16][17][18][19] and simple TN ansatzes such as MPS have limited applicability.The exponential complexity of the full quantum state evolution originates from an extremely non-local pattern of correlations, and might be circumvented by focusing instead on a local description of the state.Significant simplifications using this approach have already been observed [20][21][22][23][24][25][26][27][28][29][30].
Here we develop this idea and propose a new algorithm that explicitly identifies long-range entanglement in the system and trades it for mixture (see also [26]).In particular, our method focuses on separating fast and slow propagating degrees of freedom.Already from an early time we here observe that the fast degrees of freedom mediate non-local correlations whose effect on local observables is similar to that of a statistical mixture.
By trading such correlations for the corresponding mixture we devise a TN algorithm that allows to simulate the evolution of local observables for considerably longer times than what was achievable before, with finite computational resources.In fact, our algorithm provides an effective description of the time-evolved state as a density matrix, represented by a matrix product operator (MPO) with bounded bond dimension, that accurately captures the long-time behavior of local operators.
We support these claims by providing benchmark numerical results for quenches in the transverse field Ising model and its non-integrable generalization.We find that our algorithm predicts accurately the long-time values of local observables as dictated by, respectively, the analytical solutions and the prediction of the diagonal ensemble [31][32][33].
Quenches and quasiparticles.We focus on the outof-equilibrium dynamics induced by quantum quenches [14,35], in which a one-dimensional system in the thermodynamic limit is prepared at t = 0 in a product state |ψ 0 ⟩, and later evolved with an entangling Hamiltonian H. Our aim is to compute the time-dependent expectation value of a local observable O.
The initial state has a finite energy density above the ground state, and thus a large occupation of excited states, often described as quasiparticles states (QP).The subsequent dynamics can be described as the radiation of entangled pairs of QPs [14].In translational invariant systems indeed QPs possess a well-defined energymomentum dispersion relation ϵ(k), and a QP wavepacket centered at k 0 propagates with group velocity v k0 = ∂ k ϵ(k)| k0 , while momentum conservation enforces equal occupation of k and −k states that propagate in opposite directions.Identifying long-distance entanglement.In order to identify the long-distance contributions to the entanglement we need to focus on a subsystem S of ℓ neighboring spins, and label L and R the remaining left and right regions of the system (see Fig. 1).Entanglement across a bipartition (e.g.L vs SR) is created by correlated pairs of QPs with support on both sides of the cut.We define as long-distance entanglement the contribution to the entanglement generated by any QP pair supported on L and R, but not on S.
We can visualize this in the cartoon of Fig. 1(a), where we consider a quench that excites a single pair of entangled QPs, initially located inside S, that later travel with opposite velocities ±v.[36].In this simple example, S becomes entangled with the rest when one member of the pair leaves the region, but when both QPs have left, after t ∼ ℓ/v, S is again in a product state with respect to LR.At this stage we can factorize the total state into the product of an entangled state of LR and a state of S.
In a more general scenario such complete factorization is not possible, since S and LR might contain other sources of entanglement.The simplest generalization depicted in Fig. 1(b) involves two pairs of QPs, one fast (orange) and one slow (blue).When the former has left S the latter is still partially in S. At that stage, it is still possible to disentangle a part of L and R from S. We can indeed identify a factorization of the Hilbert spaces, In the simple scenario presented above, L f ⊗ R f is defined by the degrees of freedom that describe the fast pair, and we can factorize the state as |ψ (slow) where the second factor captures the pair of fast modes.
Trading long-distance entanglement for mixture.In a translationally invariant system, we can repeat the cartoon picture above for each block of ℓ sites.When computing the expectation value of a local observable, we need to trace out all the degrees of freedom but the ones where the observable is defined.Given that fast constituents of entangled pairs of QP are separated by at least ℓ sites, one of the two partners will always be traced out in this procedure leaving the other in a mixed state.
Thus, when focusing on local observables, we can describe the system as a mixed state, originated by having discarded the coherence between each pair of fast QP.Their joint pure entangled state is substituted by a mixed state built as the product of the two mixed states obtained by tracing out the partner in the pair, where ρ are the reduced density matrices of each QP.This provides a more efficient local description in terms of entanglement, as the long-distance components have been removed.
From the QP intuition to a TN algorithm.
We can use the above observation in an actual TN algorithm.Standard MPS algorithms for this setup attempt to represent the time-evolved state as a MPS, parametrized by one or few tensors A m of small size d × D × D, where d is the physical dimension of the chain sites and D the bond dimension [37][38][39].Since the half-chain entropy of the MPS is upper-bounded by log D [4,40,41], the typical linear growth in time of entropy after a quench [14,16] requires an exponential increase of the tensors dimensions to maintains a constant precision.As a result, given finite computational resources, standard TN algorithms [42] only give reliable predictions for relatively short times After identifying the optimal unitaries UL and VR we find a residual entropy shown as a function of time for the Ising model ( 4), (a) for several integrable quenches (g ̸ = 0, J2 = 0) and (b) non-integrable ones (g = 2, J2 ̸ = 0), always starting from the state |X+⟩.The amount of entanglement between L and R we are able to isolate with the identified unitaries is quantified by the logarithmic negativity [34] between L f and R f and is reported as EN = • • • for each of the lines.The inset of (a) shows the effect of the block size ℓ for the case (g = 2, J2 = 0) (from darker to lighter blue, ℓ = 2, 4, 6, 8).Last, the inset of (b) shows the logarithmic negativity of the reduced density matrix for the fast degrees of freedom, according to the partition L f vs R f .
(typically of the order Jt ≃ 10 with J the relevant energy scale).
Translating the QP intuition above to the TN setting we can however obtain a more efficient TN description for the local observables.Let us, for simplicity, assume that each MPS tensor represents precisely ℓ sites [43].
Singling out one subsystem, the whole state can thus be written as where the sum is over orthonormal bases of the block and its left and right environments {|Φ L/R ⟩}, meaning that we use the gauge in which the tensor C s ℓ αβ for the subsystem {|s⟩} is the orthogonality center of the MPS [13,37,44,45].
If there are long-range entangled degrees of freedom that decouple from S, the state will have a product structure, with one component completely disentangled from the physical degrees of freedom.There will thus exist basis transformations (disentanglers [46,47]) on L and R identifying the decomposition as depicted in Fig. 2(I.a).Assuming such U L and V R exist, they need to be determined variationally, for example by minimizing the Euclidean distance between the left and right-hand side of (3).Namely, given the evolved state in the form (2), and for fixed dimension d fast of the long-range component, we iteratively optimize each of the disentanglers U L and V R and the vectors |ψ LsSRs ⟩ and |ϕ L f R f ⟩ until we reach the minimal possible Euclidean distance among the right and left hand sides of Eq. ( 3) [48].Since the dimension d fast of the fast factors L f and R f is not known, the procedure needs to be repeated for different trial values.
Assuming the above procedure is successful, we now can transform the identified long-range entanglement into mixture by applying the substitution (1) to each consecutive subsystem in the chain.As schematically shown in Fig. 2(I.b), this exchanges the pure MPS description of the system by a mixed (purified) MPO with smaller bond dimension D/d fast , at the expense of additional (purification) indices.
When the decomposition (3) is exact, this mixed state has the property that the reduced density matrices for S ℓ R (resp.LS ℓ ) are unchanged [49].As a consequence, also the reduced density matrix for two adjacent blocks, and thus all observables with support on up to ℓ sites are preserved.We can then continue the evolution with the original Hamiltonian acting on the physical indices.Fast degrees of freedom will continue propagating entanglement through the subsystem S. Thus, the step of finding and mixing the long-range entangled components is repeated periodically.
An improved TN algorithm.In practice, the optimal decomposition we are able to identify using the strategy just outlined still retains some residual entropy between fast and slow degrees of freedom S(ρ fast ).We apply the truncation when this entropy falls below a given threshold η S .This results into small errors when building the The purification tensors are found variationally by minimizing the distance between the reduced density matrices ρ SR (and ρ LS ) obtained from the original state (2) and from the purification.As initial guess we use the tensors obtained from solving (3), and we use gradient descent for the minimization problem (see [49]) .
Notice that the number of purification legs increases after each mixing step.To handle this in practice, we group together all the purification indices between two physical sites and apply a cutoff to the total purification dimension [50].
Numerical results.We benchmark the algorithm using the transverse field Ising model with an additional integrability-breaking next-to-nearest neighbor interaction, Initially in a product state |X+⟩ ≡ ⊗ 1 √ 2 (|0⟩ + |1⟩), corresponding to the ground state deep in the paramagnetic phase (g → ∞), the system is quenched to a finite value of the transverse field g (and potentially J 2 ).
To explore systematically the structure of the timeevolved MPS, we simulate accurately the short-time evolution using the infinite time-evolving block decimation (iTEBD) algorithm [37] with a large enough bond dimension.We then probe the disentangling of fast and slow degrees of freedom at different times solving (3) for a subsystem of ℓ = 2 sites for several quenches.
As shown in Fig. 3, after a certain time we can identify subsystems of L and R that practically disentangle from the local region and carry long-range entanglement between the environments, measured by the logarithmic negativity [34] of the reduced density matrix for the fast degrees of freedom tr LsSRs |Ψ⟩⟨Ψ|.For all quenches, the residual entanglement entropy between the fast degrees of freedom and the rest decays fast with time, at a slower rate when the quench is into the ordered phase [51].Also as expected, the time at which fast degrees of freedom start decoupling from the local system depends linearly on its size ℓ, as illustrated in the inset of Fig. 3(a)].
We now use the TN algorithm described above and check its performance simulating the long-time evolution of local observables.Fig. 4(b) shows the results for the integrable quench g = 2, J 2 = 0 and the only nonvanishing (due to the Z 2 symmetry) single-site observable ⟨σ x ⟩.Solving (3) and trading entanglement by mixture captures the qualitatively correct dynamics, but introduces discontinuous jumps (orange line and lower-left inset) that we attribute to the decomposition not being exact.The improved heuristic algorithm achieves much better results (blue line and lower-right inset).While the approximation induces some residual oscillations, these are very small and close to the exact equilibration results even after considerably long times.In contrast, the iTEBD results with a maximal bond dimension D = 500 (pink line) start to severely deviate from the exact solution at around t = 10.
Our algorithm uses instead a much smaller bond dimension: during the time evolution steps, the required bond dimension of the purified MPS grows exponentially, but each time we trade some of the entanglement for mixture, the bond dimension gets halved [52] As we iterate the procedure, the maximum required bond dimension tends to a constant D ∼ 100 (upper inset).
In Fig. 4(b) we repeat a similar analysis for the nonintegrable quench g = 2, J 2 = 0.1.The different colors show the results for different thresholds for the residual entropy η S .Smaller values require longer evolution between mixing steps, and thus larger bond dimension, but improve the results systematically.
Our results systematically approach the long-time limit predicted by the diagonal ensemble ⟨σ x ⟩ DE = 0.852 [53].Even for the largest threshold, the relative error in the time-averaged value after t ∼ 50 is below 1% (0.848).Also the time-averaged magnetization exhibits a similar precision, as shown in the lower inset.Also in this case the largest bond dimension saturates with t (upper right inset).
In the supplementary material, we also repeat the calculations in the integrable case directly using the freefermionic formalism where we can track the coherence we discard in our simulation and confirm it does not play a relevant role in the long-time dynamics of local observables [49].
Discussion.By identifying the long-range entanglement produced in the out-of-equilibrium dynamics after a quantum quench and converting it into mixture, we have proposed an explicit approach to avoid the entanglement barrier and simulate the out-of-equilibrium dynamics of an infinite quantum chain with MPS using finite computational resources.
Our approach is inspired by the intuitive understanding of entanglement dynamics in terms of the radiation of QP.It relies on the hypothesis that fast degrees of freedom propagate correlations to steadily growing distances, contributing to the linear growth of entanglement across the system, but only as statistical mixture to sufficiently local observables.Our numerical results for the Ising chain show that the intuition is accurate in the freefermionic case, which best fits the QP picture.
We have generalized our intuition to a heuristic algorithm that goes beyond the QP picture and also performs well in non-integrable regimes of the model.It is thus important to pursue further characterization of the algorithm in order to chart its potential and limitations.
LT would like to thank Jacopo Surace who has contributed to developing the original ideas that have motivated this work.This work was partially supported by the Deutsche Forschungsgemeinschaft (DFG, German Research Foundation) under Germany's Excellence Strategy -EXC-2111 -390814868; 499180199 and and by the EU-QUANTERA project TNiSQ (BA 6059/1-1).LT acknowledges support from the Proyecto Sinérgico CAM 2020 Y2020/TCS-6545 (NanoQuCo-CM), the CSIC Research Platform on Quantum Technologies PTI-001 and from Spanish projects PID2021-127968NB-I00 and TED2021-130552B-C22 funded by MCIN/AEI/10.13039/501100011033/FEDER,UE and MCIN/AEI/10.13039/501100011033,respectively.In the main text we discussed a way to transform the long-range entanglement generated in a quench into mixture using tensor networks.The transformation is such that it reduces the bond dimension of the TN, while aiming to preserve the local observables and their evolution.In this section of the supplementary material we provide a complete description of the algorithms we use for that purpose.First, in order to set the notation and introduce some useful concepts, we briefly describe the uniform MPS formalism.For a more complete description, we refer the reader to the excellent review [1], whose notation we follow here.

Uniform MPS
Uniform MPS (uMPS) are TN states that represent one-dimensional quantum many-body states in the thermodynamic limit.These states are constructed by repeating infinitely many times a tensor along a line (or repeating in a cyclic manner a finite set of them, if we work with a unit cell larger than one).In graphical TN notation, Where A is a d×D×D tensor, d being the physical dimension of the sites of the chain and D the bond dimension.
A central object in uMPS simulations is the transfer operator, defined as the following four-legged tensor: We can interpret this object as a D 2 ×D 2 matrix (left vs. right indices).Many of the important properties of the state are captured by the dominant left and right eigenvectors of this matrix.For instance, they determine the reduced density matrix of any (connected) subsystem and thus appear in the calculation of any local expectation value.In the generic case, these dominant eigenvectors are non-degenerate [2].If the state is normalized, the largest eigenvalue is 1, and the leading eigenvectors correspond to the fixed points of the transfer operator seen as a map from the left to the right indices and viceversa.It is possible to exploit the gauge freedom in the TN to transform these eigenvectors into a canonical form, which is particularly simple.In particular, in the so-called leftcanonical form, the MPS tensor for a normalized state satisfies the conditions where λ 2 is a diagonal matrix with positive entries.Cutting (1) in half across one cut defines D states for the left half-chain |Φ L α ⟩, and the same number for the right half-chain, indexed by the (now open) virtual leg.The first condition in (3) ensures that the D states on the left half-chain are orthonormal.Alternatively, one can impose a right-canonical form A R in which the conditions are the symmetrical ones and orthonormality is satisfied by the states to the right of the cut.
In practice, it is instead common to work in a mixed gauge, in which one can define [1] a single central tensor C s αβ = γ (A L ) s αγ λ γβ , where α (γ) is the left (right) virtual index of the tensor A L , and s denotes its physical index.The state can then be written as a half-chain of left-canonical tensors to the left of C and a half-chain of right-canonical ones to the right, such that both left and right virtual bonds of C are in orthonormal bases.Taking the tensor product of these two bases with the single-site physical basis, the state can be written as, We can rewrite the full uMPS in terms of the single tensor C, as .

arXiv:2308.04291v1 [quant-ph] 8 Aug 2023
In this work, we are interested in the long-range entanglement between distant regions of the system L and R, separated by a block S of finite size ℓ.Thus, it is convenient to use the central canonical form with respect to the whole block S, such that we can study the longrange entanglement at the level of the virtual indices of the MPS.After blocking the ℓ sites of the subsystem S, the wave function of the system can be written as As we will see in the next three subsections, our study of the long-range entanglement will allow us to construct a low-rank decomposition of the tensor C when partitioned between L and SR (or LS versus R), at the cost of introducing indices that couple the tensors in the bra and the ket.This decomposition will be such that the local expectation values and their evolution will be preserved, while the long-range coherences will be discarded.

Detection of long range entanglement
In this subsection we describe the algorithm that we use to identify the long-range entanglement in our MPS.
As discussed in the main text, if there are long-range entangled degrees of freedom that decouple from S, the state will factorize in a tensor product with one component completely disentangled from the physical degrees of freedom in S.There will thus exist basis transformations on L and R identifying the decomposition L = L f ⊗ L s and R = R f ⊗ R s , such that the degrees of freedom in the Hilbert spaces L f and R f are in a tensor product structure with the rest of the system.Graphically, .
To look for such a decomposition we numerically minimize the square of the Euclidean distance between the left and the right-hand side of the previous diagram (7).
The variables in the optimization problem are the unitary matrices U L and V R , which respectively determine the left and right basis transformations, and the two resulting vectors, |ψ LsSRs ⟩ and |ϕ L f R f ⟩.The full optimization is a complicated non-linear problem.However, as a function of each of the individual pieces, the cost function is quadratic and can be solved using basic linear algebra.Hence, to find the decomposition we use an alternating scheme, first used in [3]: we fix all the variational tensors but one, and substitute the free piece by the one that minimizes the value of the cost function.When the unitaries are fixed, the optimal for the two vectors are just the leading left and right singular vectors across the partition L s SR s versus L f R f .For the case of the unitaries, after imposing the unitarity constraint, the cost function is linear in each of them.That kind of problems can also be easily solved with a singular value decomposition of the environment [4].
The algorithm is efficient, with cost scaling as O(D 3 ) with the bond dimension, and we find that it shows a good performance.In particular, we find components of the wave function which mediate long-range entanglement and decouple from the slower degrees of freedom, as expected according to the intuition drawn from the quasiparticle picture.

Simple truncation algorithm
In the previous section we have discussed the TN algorithm we use to identify the degrees of freedom that contribute to the long-range entanglement in our MPS representation of the time-evolved state.In practice, we combine the previous algorithm with an iTEBD simulation that provides us with a faithful representation of the state at short times and later allows us to evolve the local tensors in the ansatz.
Starting from a pure state (uMPS) description of the initial state, we evolve unitarily with iTEBD.At periodic time intervals, we look for the above decomposition of the simulated wave function.The extent to which our time-evolved state has the sought structure can be characterized in several ways.We choose to use the entanglement entropy between the fast and slow degrees of freedom, which can be identified after we have transformed the basis of the bond dimension indices according to the found unitaries.We thus need to compute the entropy of the reduced density matrix ρ (fast) which can be graphically represented as A vanishing entanglement entropy means that our state has, in the found virtual bases, exactly the desired tensor product structure, |Ψ⟩ = |ψ In case, it is possible to modify the state in such a way that local expectation values are preserved while discarding some long-range correlations, as described in the main text.More concretely, we substitute the density operator of the full system |Ψ⟩⟨Ψ|, represented graphically in Fig. S1 (left), by a mixed state (right) that, in the basis of the virtual degrees of freedom in which the state factorizes, corresponds to with ρ as defined in the main text.
The state on Fig. S1(a) is an exact representation of the tensor C in the case in which the time-evolved state presents the exact sought structure.The state on the right (Fig. S1(b)) has a different structure: it represents a mixed state.Despite the two states being globally different, the two have the same reduced density matrices on the subsystems LS, as graphically shown in Fig. S2.An analogous computation shows that the reduced density matrix for SR is also preserved by this substitution.Furthermore, as can be seen from Fig. S1, the rank of the state (or bond dimension) across the cut L versus SR is reduced by a factor which is given by the dimension of the fast degrees of freedom.L f R f ⟩ is normalized, and thus, the tensor contraction in the righthand side of the equation in (b) gives 1.This same diagram also shows that the fixed points are unchanged by the mixing step: notice that by tracing subsystem R, we get a condition equivalent to the second equality in 3 (namely, that Both the state before and after the truncation give the same resulting tensor.A similar argument can be made for the left fixed point.S3.Reduced density matrices on two consecutive blocks of the system, before (a) and after (b) the truncation of the fast degrees of freedom.As explained in the text, the fixed points are also unchanged by the transformation.Using the unitarity of the two basis transformations in the virtual degrees of freedom, one can show that the expectation value for both is given by the same contraction (c).
Because the state of the chain can also be explicitly written as a sequence of blocks, using the form (5), we can perform the same transformation in all blocks of the chain at the same time [i.e., we substitute all the tensors C in (5) as shown in Fig. S4].After this substitution, it can be shown that the leading left and right eigenvectors of the transfer operator for the block are unchanged (see Fig. S2).This can be used to show graphically that the reduced density matrices of up to two blocks (2ℓ consecutive sites) are preserved (Fig. S3).In practice, as can be seen from Fig. 3 of the main text, the decomposition we obtain variationally never reaches full separation between fast and slow degrees of freedom and as a result the substitution is never exact.However, the fast-slow entropy seems to decay fast with time.Our simple truncation consists just on performing an iTEBD simulation, while monitoring the "disentangling" entanglement entropy.When it falls below a certain threshold, we substitute all the C tensors on the present state by its tensor product approximation, and mix the fast degrees of freedom according to the diagrams above.After that, because the last step is not exact, we normalize again the state.The results from this evolution are the orange line in Fig. 4(a) of the main text.In them, it can be appreciated that at the times where the truncation is performed, we incur in an instantaneous error in the local observable.However, as we evolve the state again the dynamics seems to be correctly captured modulo the instantaneous shift introduced by the truncation.The same phenomenon repeats itself as we perform more and more truncations in the same simulation.A closer-up figure of the results for this truncation in the main text can be seen in Fig. S5.The qualitative agreement with the true dynamics of this simple algorithm corroborates the intuitive picture discussed in the main text, and justifies using these tensors as starting point for an improved method that corrects for those discontinuous jumps in a heuristic manner.We discuss the improved technique in the next section.

Heuristic truncation
Inspired by the results of our initial truncation, we present in this section a small modification of the previous algorithm that attempts to correct the shifts induced by the fact that we do not find an exact factorization between fast and slow degrees of freedom.
Our improved heuristic algorithm takes advantage of this fact and tries to find an alternative purification form that generalizes the one of the simple algorithm and guarantees the correct marginals for LS and SR.The ansatz we use is composed by three tensors M L , B ℓ and N R .We initialize their values from the pieces of the decomposition recovered by the simple algorithm and modify them to ensure the desired property (see Fig. S6 and Fig. S7).
In order to improve the purification, we minimize the sum of the Euclidean distances between the new and previous reduced density matrices, namely we are interested in FIG.S7.In order to improve our purification, we minimize the squares of the euclidean distances between the left and the right hand side of the two approximate equations above.We minimize the sum of the two distances, as we need a resulting state satisfies both.
where ρ denotes the new reduced density matrices obtained from the mixed state defined by the tensors M L , B ℓ and N R , which are the variables in the optimization.To minimize the cost function we use a gradientdescent scheme, as the optimization contains terms which are of fourth order in the entries of the tensors.The dimensions of the tensors, and the fact that they have purification legs guarantee that the original pure state will be approximated by a mixed state with a lower bond dimension.
The same as in the simple algorithm, we combine this decomposition with an iTEBD simulation.As the evolution unfolds, we monitor the behavior of the entanglement entropy between fast and slow degrees of freedom.Once it falls below a pre-established threshold η S , we perform the truncation described by the direct decomposition and group the pieces obtained in it into the tensors M L , B ℓ and N R , as can be seen in Fig. S6.Lastly, we use these tensors as a starting guess for the heuristic optimization described in the paragraph above, which modifies them slightly to correct for the small deviation in the reduced density matrices.Once the optimization is stopped, we substitute the pieces for all the tensors C in our chain, normalize again the state and continue the evolution.After the tensors C are substituted by their lower-rank counterparts for the first time, the structure of the uMPS describing the purification contains a sequence of tensors B ℓ that carry the physical indices, separated by two purification tensors N R and M L , themselves connected by a diagonal of inverse Schmidt values.Notice that the subsequent time evolution is applied only to the physical indices of each B ℓ tensor, but not to the purification ones.Nevertheless, these tensors mediate the entanglement between physical sites, which in general will continue building up in time.We thus need to continue monitoring the fast-slow entanglement entropy of the B ℓ tensors in order to apply further truncations as required.
Two extra aspects of the algorithm deserve a comment.The first one is that with each truncation/mixing step, the number of purification sites between a pair of blocks increases by two.In principle, since there is a large freedom in the representation of the purification indices, it would be possible to look for a more efficient description of this set of tensors (e.g.minimally entangled, in the spirit of [5]).In this work, for simplicity, we group all the purification legs together and work with a unit cell that has two tensors, one with the physical dimension and another one with the purification dimension.The dimension of the purification index would thus in principle grow exponentially with the number of truncations.In order to efficiently deal with it, we perform a SVD of the purification tensor, splitting virtual vs purification legs, and implement a hard cutoff on this dimension (with a value of a d purification = 1000, for our simulations).Then we continue the evolution while keeping the purification dimension fixed.The second aspect that deserves a mention refers to the performance of the heuristic truncation without using the input from the entanglement decomposition.We observe that whenever we do not make use of it, the heuristic optimization results (i.e. the distance between the correct marginals and the density matrices of the varitional mixed state) are not good, and thus, we incur in errors that accumulate during the simulation.

FREE FERMION CALCULATIONS
In the integrable case (J 2 = 0), the model can be mapped to an exactly solvable free fermionic system.We can use this solution to perform the same calculations as with the TN algorithm, and thus corroborate and gain a deeper understanding of the observed numerical results.
Disentangling the long-range components Fig. 3(a) of the main text shows that, for quenches from the ground state at g → ∞, the residual entropy between fast and slow degrees of freedom decays slower for smaller values of g, even though for all the cases studied the decay is compatible with an exponentially fast decrease in time.At the same time, the saturation value of the logarithmic negativity carried by the fast degrees of freedom seems to grow as we scan the final value of the transverse field from g = 2 to g = 0.5.In fact, at g = 1, the logarithmic negativity seems to saturate to 1, the maximal possible value for a bipartite state of two qubits [6].
We can relate this behavior to the distribution of quasiparticles in the initial state, illustrated in figure S8 for the studied quenches.As we scan from g = 2 to g = 0.5 in the phase diagram, the occupation number of the freely propagating modes of the quench Hamiltonian grows.Something similar happens to the coherences between opposite momenta.This is reminiscent of the behavior of the logarithmic negativities in the TN computations and of the fact that as we quench to smaller values of g, the energy of the initial state according to the final Hamiltonian gets closer to the middle of the spectrum.Furthermore, by analysing the group velocities one can get some intuition of the rates of decay of the residual entropy.For a quench inside the paramagnetic phase, the maximum of the occupation coincides with the maximal group velocity.In contrast, when we quench to the ferromagnetic phase, the maximally occupied mode has vanishing group velocity.This could be an explanation of the slower decay rate observed when quenching to the ferromagnetic phase.

Eliminating the long-range coherences
The free fermionic chain can also provide insight on the transformation performed by the TN algorithm.A Jordan-Wigner transformation maps the spin chain to a system of N fermionic modes, created by operators a † i , i = 1, . . .N , that fulfill standard commutation relations a † i , a j = δ i,j , with δ the Kronecker delta.Because the resulting Hamiltonian involves only quadratic terms in a i , a † j , the fermionic system is free, and ground states and their evolution are described by Gaussian states.By virtue of Wick's theorem, these are completely characterised by two-point correlation functions, which can be organized in the 2N × 2N correlation matrix [7] Γ i,j = ⟨α i α † j ⟩, where the operator-valued vector N holds all creation and annihilation operators of the chain.Other sets of fermions c = U α can be defined through symplectic transformations, which conserve the anticommutation relations.
We can now divide the system in three parts, namely a central region S containing few fermions and playing the role of the local block in the above description, and two regions L and R, containing respectively all the fermions to the left or right of S. Because the state is Gaussian, all correlations of operators supported in one connected region are determined by the restriction of the correlation matrix to the modes in that region.The symplectic transformation representing a unitary that acts on the left (right) virtual leg of the MPS would in this setting act only on the modes defined in L (R).FIG.S9.We reproduce the TN calculations for the integrable quench g = 2 using the free fermionic formalism.For a system of size N = 800 we consider a block S of 4 sites and identify a pair of modes, respectively in L and R, that disentangle from S, as explained in the text.The residual fast-slow entropy (a) shows the same behavior as in the TN calculation, as does the logarithmic negativity between the L and R fast components (inset).(b) shows the contribution to the correlation matrix of the coherences of this long-range pair when it is identified at t = 5 (inset) and after an additional evolution for ∆t = 20.
The free-fermionic equivalent of the TN algorithm should try to single out, after some evolution time, at least one fermionic mode localized in L and entangled only with modes in R. Namely we want to identify modes l = j U L 1,j α L j and r = j U R 1,j α R j such that ⟨λ • α S † ⟩ = ⟨ρ • α S † ⟩ = 0, where ρ † = (r, r † ), λ † = (l, l † ) and α S(L,R) is the restriction of α to the modes in S (L, R).To this end, after discarding modes that are pure in either L or R, we optimize numerically the linear combinations above that produce the lowest overlap with modes in S.
In order to compare to the TN results, we have applied this procedure to the quench g = 2 for a system of N = 800 fermionic modes, in which we singled out a central region of four sites.After time t = 5, we solve the optimization described above to identify a pair of fermionic modes carrying long-range entanglement.Fig. S9(b) shows (in absolute value) the contribution of this pair to the (a † a part of the) correlation matrix (inset) as well as how this has evolved after a time ∆t = 20 has elapsed.The empty region around (0, 0) corresponds to our block S.
In panel (a) of the same figure we show how the participation of those modes changes in time.In particular we observe that, despite the fact that these mode spread and eventually enter S, their participation on the S modes is extremely small compared to their participation elsewhere.Therefore at the level of the local system S, this pair can be substituted by a density operator, completely disentangled from the local degrees of freedom, starting at time t = 5 and for the rest of the evolution.Fig. S9(a) shows the residual entanglement between the l, r pair and the modes in S (main panel) for several values of the quench parameter, together with the negativity between the two modes l and r (inset), to be directly compared to the results in Fig. 3 in the main text.Both quantities show good agreement with the TN results and confirm that in all scenarios, as time passes, the identification of l, r becomes more and more accurate, and that both modes keep being robustly entangled with each other.Notice that the agreement is not only qualitative, but the magnitude of the entanglement quantities is very similar to those found with the TN algorithm, which supports our interpretation of the TN transformation as correctly identifying the quasiparticles that carry away long-range entanglement.

FIG. 1 .
FIG. 1.(a) A correlated QP pair initially located inside region S creates long-distance (pure LR) entanglement when both QPs propagate outside of S. (b) At that time, slower QPs (blue) can still contribute to entanglement between S and its surroundings.

FIG. 2 .
FIG. 2. First iteration in the transformation of long-range entanglement into mixture in TN. (I.a) Decomposition (3) for a block S, represented by a single tensor; (I.b)The density matrix of the full state (left) can be written in terms of the new tensors (middle) when we apply the decomposition to each block [black circles represent the inverse Schmidt values matrix, inserted to compensate the central gauge in (I.a)].We then substitute each long-range component by the product of its marginals, giving rise to a mixed description (right).(II.a)The heuristic algorithm directly searches for tensors that (approximately) preserve the reduced density matrices for LS and SR.(II.b)The density matrix for the full state is replaced by a purification defined by the solution.The structure is analogous to that in (I.b) (in each case we have indicated the relevant dimensions).

FIG. 3 .
FIG. 3. The entropy of the reduced density matrix for the fast degrees of freedom (c), computed from the decomposition in fig.2(I.a), measures the residual fast-slow entanglement.After identifying the optimal unitaries UL and VR we find a residual entropy shown as a function of time for the Ising model (4), (a) for several integrable quenches (g ̸ = 0, J2 = 0) and (b) non-integrable ones (g = 2, J2 ̸ = 0), always starting from the state |X+⟩.The amount of entanglement between L and R we are able to isolate with the identified unitaries is quantified by the logarithmic negativity[34] between L f and R f and is reported as EN = • • • for each of the lines.The inset of (a) shows the effect of the block size ℓ for the case (g = 2, J2 = 0) (from darker to lighter blue, ℓ = 2, 4, 6, 8).Last, the inset of (b) shows the logarithmic negativity of the reduced density matrix for the fast degrees of freedom, according to the partition L f vs R f .

FIG. 4 .
FIG. 4. (a) Evolution of the transverse magnetization for the integrable quench (g = 2, J2 = 0) as obtained with the direct decomposition shown in fig.2(I) (orange) and the heuristic algorithm of fig.2(II) (blue).The latter shows very good agreement with the exact result (dashed black line), and converges towards the equilibration value (dashed grey).The iTEBD simulation with D = 500 (pink), in contrast, deviates at times O(10).The lower insets show close-ups of the shortand long-time regions.In the left one, the vertical dashed lines indicate mixing steps, inducing jumps in the simple but not the heuristic algorithm.(b) Same observable for the nonintegrable quench (g = 2, J2 = 0.1) using the heuristic algorithm.Different colors indicate different residual entropy thresholds ηS.Our simulations converge towards the diagonal ensemble value ⟨σx⟩DE (dashed gray line).This can be more clearly seen in the time average of the observable (lower inset) or the value of the time average at time t as a function of ηS (upper-left inset).For both quenches, the upper (right) inset shows the bond dimension as a function of time.

1
Max-Planck-Institut für Quantenoptik, Hans-Kopfermann-Str. 1, D-85748 Garching, Germany 2 Munich Center for Quantum Science and Technology (MCQST), Schellingstr.4, D-80799 München 3 Institute of Fundamental Physics IFF-CSIC, Calle Serrano 113b, 28006 Madrid, Spain DESCRIPTION OF THE ALGORITHMS FIG. S1.Graphical representation of the density matrices of the entire state before (a) and after the truncation of the fast degrees of freedom (b).
FIG. S2.Reduced density matrix on the subsystem LS before (a) and after (b) the truncation.Using the unitarity of VR, one can check that the reduced density matrices are identical (right-hand side of the equalities).The state |ϕ (fast) FIG.S3.Reduced density matrices on two consecutive blocks of the system, before (a) and after (b) the truncation of the fast degrees of freedom.As explained in the text, the fixed points are also unchanged by the transformation.Using the unitarity of the two basis transformations in the virtual degrees of freedom, one can show that the expectation value for both is given by the same contraction (c).
FIG.S4.After identifying the long range entanglement, we truncate the fast degrees of freedom in all block simultaneously.
FIG.S5.Evolution of the transverse magnetization for the integrable quench (g = 2, J2 = 0) as obtained with the algorithm discussed in the section above, in orange.The exact result is shown by the black dashed line.The inset shows a close-up of the evolution.In the inset, the vertical gray dashed lines show the times in which the truncation is performed.

FIG
FIG. S6.(a) Tensor C obtained from the uMPS dynamical simulation in the mixed gauge.(b) Using the entanglement decomposition, we find an approximation to the tensor C in terms of the different pieces that mediate the fast and slow degrees of freedom.(c) We group the pieces according to the dotted lines to construct the new tensors ML, B ℓ and NR, which will be the starting point for the heuristic optimization variables.
FIG.S8.Occupation of the eigenmodes of the quench Hamiltonian (solid lines) and their group velocity (dotted lines) as a function of their momentum in the integrable quenches shown in figure3(a) in the main text.Dashed lines show the initial coherence between opposite momenta.In the leftmost plot (g = 0.5), the coherence between opposite momenta and the group velocity are on top of each other.