On temporal entropy and the complexity of computing the expectation value of local operators after a quench.

We study the computational complexity of simulating the time-dependent expectation value of a local operator in a one-dimensional quantum system by using temporal matrix product states. We argue that such cost is intimately related to that of encoding temporal transition matrices and their partial traces. In particular, we show that we can upper-bound the rank of these reduced transition matrices by the one of the Heisenberg evolution of local operators, thus making connection between two apparently diﬀerent quantities, the temporal entanglement and the local operator entanglement. As a result, whenever the local operator entanglement grows slower than linearly in time, we show that computing time-dependent expectation values of local operators using temporal matrix product states is likely advantageous with respect to computing the same quantities using standard matrix product states techniques.

The complexity of simulating quantum-many body systems increases exponentially with the number of its constituents.Over the last decades, the development of tensor networks (TN) techniques has however helped in gaining better insight on the equilibrium properties of many-body quantum systems.It is now understood that, at equilibrium, quantum complexity is mostly related to the amount of entanglement, and we have been able to design TN Ansätze that can encode the structure of typical equilibrium states.
Out of equilibrium, the situation is different: even in the simplest protocol, such as e.g.quantum quenches, correlations spread over large distances, quickly producing robustly entangled states [1].Standard TN techniques such as those based on matrix product states (MPS) thus struggle to cope with the fast growth of entanglement and, as a result, their cost increases exponentially with the duration of the evolution [2]: this is often referred to as "entanglement barrier" in the literature [3].
This fact is a consequence of trying to represent the full quantum state during the evolution.The situation might change if we focus on a local description of the state, by trying to describe the evolution of the expectation value of local operators.In principle, this is a much simpler task, and several approaches along these lines have been proposed [4][5][6][7][8][9][10][11][12].However, except for few studies on integrable systems [13,14], there has been no systematic understanding of the real computational cost of such approaches and thus no concrete understanding on the complexity of simulating the evolution of local observables.
In this work we make a step in this direction.We consider the evolution under a local Hamiltonian of the expectation value of a local observable and follow the inspirational papers [4-6, 10, 12].These algorithms rely still on a matrix product state, which is however now defined in time, making it a temporal MPS (tMPS) along a Keldysh contour [6,15].It has recently been noticed that such tMPS encodes the influence matrix of the system [16] that drives the evolution of a region of the full system, and thus provides the systematic way of translating the global evolution into a local one [11,17,18].In specific cases, such tMPS can be described with small bond dimension, though this is not always the case [4,5].
Here we provide some theoretical backup to these numerical observations.Our main result in this direction is to bound the complexity of the tMPS with the one of encoding the time evolution of operators in the Heisenberg picture.Such complexity is encoded in the so called "operator entanglement" (OE) [19,20], which is expected to grow at most logarithmically in time for integrable systems, whereas it increases linearly for non-integrable systems [3,[21][22][23].The key to obtaining this result is a slightly modified version of the algorithms introduced in [4,6], which, as shown in the following, allows to build a direct connection between the two quantities.
Whenever the OE grows logarithmically in time, we can thus show that the tMPS bond dimension is bounded by a polynomial growth in time.On the other hand, when the OE grows linearly in time, our bound on the tMPS bond dimension is exponential.While this does not rule out that there might be specific cases in which the tMPS might still have a small bond dimension and thus provide a computational advantage, see e.g.[22][23][24], we believe that these are but fortuitous exceptions to the general rule, and indeed in the scenarios we presented for this case as well as others in the literature [10] we find that such bound is saturated.

II. SETUP
Given a lattice system, we consider the complexity of computing the evolution of the expectation value of a local operator after a quench.We start from a product state |ψ 0 ⟩ and a local Hamiltonian H.The system starts to evolve and it is described by the state |ψ(T )⟩ = exp (−iHT ) |ψ 0 ⟩ ≡ U (T ) |ψ 0 ⟩.In general, the entanglement entropy of the state increases linearly in time [1] and standard time-dependent MPS simulations become exponentially expensive [25,26] (for recent reviews see also [27,28]).Here, rather than evolving the state, we focus on the evolution of the expectation value of a local operator acting on two neighbouring sites 1 , Approximating U (T ) by a sequence of short evolutions U (δt) N T with N T = T /δt and using a Trotter expansion, this quantity is encoded in the contraction of a two dimensional TN containing order N x × N T tensors [29,30], see Fig. 1(a).Furthermore, one can fold the network following [4,5], leading to a double-layer structure as shown in Fig. 1 (b).Given the Trotter approximation, the time evolution of a local operator has an exact causal cone obtained by cancelling all the unitary gates that are contracted with their respective Hermitian conjugates.As a result, the network can be simplified and acquires the triangular shape shown in Fig. 1 (c) [10,12].

A. Compressing the tMPS
Being two-dimensional, the best contraction path for this TN is not a priori obvious.a possible strategy is to identify two tMPS defining the contraction of the left and the right half of the system (see Fig. 2), and interpret the triangular network as the scalar product between the two, ⟨OQ⟩ (t) = ⟨L O |R Q ⟩2 .The identification is purely formal, since the bond dimension of the individual tMPS tensors can grow exponentially with the number of time steps.The construction of the tMPS thus only makes sense if one can show that, at least for specific scenarios, its bond dimension increases mildly with the number of time steps.
The task is therefore to identify the relevant rank of the tMPS matrices and compress them on their support.We know that the tMPS allows to compute correlation functions of the type ⟨OP(t) • • • QR(t ′ )⟩ for an arbitrary number of insertions of local operators at different times between 0 and T .Following the standard DMRG recipe [26], this requires having a faithful representation of all reduced transition matrices (RTM) T which can be seen as generalizations of reduced density matrices.
The graphical representation of these objects is found in Fig. 2(b).Interestingly, such transition matrices have also been considered in the context of holography: the properties of the equivalent in the field theory of the transition matrices described here have a geometrical interpretation in the bulk [31][32][33][34].It is then the rank of such RTMs which dictates the rank of the tMPS matrices, rather than the temporal entropy that was previously considered [4,6,10,12].
Generalizing the DMRG prescription, the bond dimension of the tMPS tensors should allow to obtain a low-rank approximation of the reduced transition matrix T with the desired precision.The reduced transition matrices are complex-valued and, for reflection invariant Hamiltonians, whenever Q = O they are symmetric.As a result, their low rank approximation is better defined in terms of their singular values, so we will project the tMPS bond dimension at time t on the largest singular values of T Following [6], we now define Λ s (t) and Λ o (t) respectively as the contractions of the overlap ⟨L O |R Q ⟩ until t (the top part of the network contraction, including the initial state) and below t (the bottom part, containing the operator), respectively, see Fig. 2. We also define Λs L (t), Λs R (t) as the contraction up to t of the network obtained by contracting With these definitions, we have that where we use the similarity symbol to indicate that the two operators share the same singular values.For reflection invariant systems, such as those we analyze here, one has Λs where R denotes the rank of a given matrix.This equation is at the basis of our analysis.Since for each choice of O and Q we obtain a different tMPS representation or ⟨L O (t)| and |R Q (t)⟩, we define the cost of the algorithm as the cost of simulating the operators that require the highest bond dimension.

B. The rank of the reduced transition matrices
Whenever O = Q = I, the transition matrices T L I |R I t for all t are projectors.As a result the states |L I ⟩ and |R I ⟩ consist of trivial singlets at the virtual level, as sketched in Fig. 2 (c), implying that the tMPS matrices have always bond dimension 1, for all t.In order to obtain non-trivial tMPS we thus include the operators in the construction of the states.
From the previous definitions it follows that meaning that Λ o (t) exactly encodes the Heisenberg evolution of the initially localized operator where ρ(t) is the time-evolved density matrix of the initial state.Finally, for Hamiltonians that are invariant under reflections with respect to the center of the chain, where T R stands for the partial transpose on the semiinfinite right part of the system.Since we are evolving a pure state, there are known relations between ρ T R (t) and ρ(t) [35], and in particular we know that R Λs [36,37].Having identified the elements composing the RTMs in Eq. ( 3), we can now use their physical properties together with Eq. ( 4) to identify some useful bound on the rank of the tMPS.A first scenario one can consider is that of a local quench.In this case, it is well known that the entropy of the evolved state only increases logarithmically with time, and as a result the rank of Λ s (t) increases at most polynomially [38,39], thus using Eq. ( 4) we have Such a case is perhaps not particularly interesting, since a similar polynomial cost is obtained also by using standard MPS algorithms.We thus turn to the more interesting scenario of a global quench.Here we know that R (Λ s (t)) increases exponentially with T , as already mentioned.If now R (Λ o (t)) only increases polynomially with T , we have and the temporal MPS strategy provides a polynomial algorithm to compute the out-of-equilibrium dynamics of local observables, as already anticipated in [10,12].In the literature, it is conjectured (and checked in several scenarios) that in the case of integrable systems the entanglement of local operators only grows logarithmically [3,19,20,22,23], implying that the rank of Λ o (t) only increases polynomially with T as required in Eq. ( 9).As a result, if the conjecture is correct, the tMPS provides an efficient method to characterise the time evolution of local operators for integrable systems.
Since generically the rank of Λ o (t) is expected to grow exponentially with t, we expect that with α and β > 0 model-dependent constants.In this case, there is no guarantee that the tMPS provides an efficient compression for the problem at hand.

C. Relation to other definitions of temporal entanglement
The basic ingredient in our formalism are the transition matrices defined in Eq. 2, and the generalized entanglement related to their rank.It is worth noting that these quantities differ from the definition of temporal entanglement recently proposed in the literature [4,6,10,12,13], which is built through density matrices obtained individually from either the left or right vectors (ie.ρ R ∼ |R⟩ ⟨R| , ρ L ∼ |L⟩ ⟨L|).
This can be most easily seen in the case of a reflectionsymmetric problem: here one has |R⟩ = | L⟩, so that the two definitions are connected by a nontrivial complex conjugation on the left state, which can be seen as a partial transpose operation acting on the legs of the folded transition matrix (see Fig. 2 (d)).This global unitary operation (the tensor product of swap gates on every physical leg, ie. at each time site) acting on one side of the transition matrix can change the entanglement properties of the objects involved, so that a clear connection between the two quantities is not evident and would definitely require a deeper investigation in the future.
Another aspect which has been recently pointed out is that the traditional definition of temporal entanglement can depend on the gauge chosen to define the transverse transfer matrix, and can be made arbitrarily small by using such gauge freedom.As a result, being able to accurately represent |R⟩ and ⟨L| does not guarantee obtaining a good precision on an expectation value ⟨O(T )⟩ = ⟨L| O |R⟩, since the leading Schmidt vectors can be orthogonal to such overlap [40].Thus, as we argued, for our intent of determining the complexity of simulating expectation value of local operators, the use of the generalized entanglement proposed here is the appropriate object to consider, close in spirit to the concept of biorthonormal truncation basis, which is already wellestablished for non-hermitian DMRG problems [40][41][42].

III. NUMERICAL RESULTS.
We now report the numerical evidence to support our claims.We consider a transverse field Ising Hamiltonian where σ x,z are Pauli matrices.We perform a secondorder Trotter expansion of the time evolution operator and cast it into a matrix product operator (MPO) following [30], with a timestep δt = 0.05.For our examples we consider two cases: an integrable one for g = 0.5, h = 0, and a non-integrable one g = −1.05,h = 0.5.We consider different initial product states, namely . The time evolution is performed by the variant of the standard folding algorithms we have described, exploiting the causal-cone of the network [4,6,10,12] and relying on a low-rank approximation of the RTMs defined in Eq. ( 2).
As in previous approaches [10,12]  for every bi-partition of the system into t and T − t.Further details about the algorithm and accurate comparison with other prescriptions in the literature are presented in the accompanying Supplementary Material.
In the case of the Ising model, we found the operator with the largest bond dimension defined at most on two sites to be a single-site σ x .In the following we will focus on optimizing with respect to it 4 .
We start by considering the bond dimension necessary to keep the truncation error above a threshold ϵ = 10 −4 in the SVD spectrum of the reduced transition matrices, and compare them with those obtained imposing the same truncation on the Heisenberg evolution for the operator O(t), which entails the evolution of the vectorized local operator under the Hamiltonian H ⊗ I − I ⊗ H, resulting in the MPS, |ψ O (t)⟩.The results are shown in Fig. 3, where we can observe that, as expected: 1) the behavior is different in the integrable and non-integrable cases, 2) the bond dimension necessary to correctly describe the RTM are always below the ones necessary to describe the Heisenberg evolution of the operator.
We note however that in order to estimate the computational cost of simulating the expectation value of local operators with a finite precision using a tMPS we should determine the bond dimension that is required to (swapping the bra with the ket legs).Such an operation does not modify the rank of each state, but it reduces the overlap between the two, so that a shift in the resulting singular values can be expected.
As a result, the computational cost required for simulating the expectation value of a local operator with a finite precision is constantly offset from that of computing ⟨L O |R O ⟩, the latter being upper-bounded by the cost of the Heisenberg evolution of the operator.The scaling in time of the cost to maintain a constant error in local observables is the same as the one of ⟨L O |R O ⟩ (see also the discussion in the Supplementary Material).
We verify this explicitly by following a procedure similar to the one proposed in [19] to estimate the truncation error in our tMPS compared to an exact result.The resulting bond dimensions for a fixed truncation error are reported in Fig. 4. For the integrable case, we find a polynomial increase, consistent with the behavior of the OE.Different initial states require different bond dimensions for a faithful description, though the growth follows the same power law.For the non-integrable case, on the other hand, the OE requires an exponentially growing bond dimension for a faithful description, and the tMPS bond dimension follows the same behavior, consistent with our result.

IV. CONCLUSIONS.
We have provided a new connection between two separate concepts, temporal entanglement and the operator entanglement.In particular, we have shown that the rank of the RTM necessary to compute the temporal entanglement is upper-bounded by the rank of the OE of a bi-partition.As a result, the scaling in time of the computational cost of simulating the evolution of a local observable with constant error is upper-bounded by the scaling of the OE, even though the exact value of the cost can be offset by a constant value whose origin has been discussed in detail.Our presented algorithm yields the most efficient performance observed thus far, with a marginal advantage over previous approaches.
As a result, we can claim that whenever the OE only grows logarithmically, using a tMPS for simulating the evolution of a local observable is the best choice in terms of computational cost and scales only polynomially with time, whereas in the generic case the entanglement barrier cannot be circumvented by using tMPS.On the other hand, a linear growth of the OE does not necessarily imply that the tMPS cannot provide an efficient compression, see eg. [22][23][24].The transverse contraction methods discussed here would then be useful even in those cases.
Our results are a step towards understanding the cost of simulating the out-of-equilibrium dynamics of quantum many-body systems.It would be interesting to check the bounds we have obtained on a larger class of integrable and non-integrable models, and more importantly to explore if the transverse contraction can help in simulating the dynamics of higher dimensional systems.
While working on this paper, we found out about Ref. [44] discussing similar issues and suggesting a connection between the OE and the temporal entanglement, a connection we hope to have elucidated with our work.MCIN/AEI/10.13039/501100011033,respectively.CR acknowledges support from a "la Caixa" Foundation fellowship (ID 100010434, code LCF/BQ/DI21/11860031).This work has been financially supported by the Ministry of Economic Affairs and Digital Transformation of the Spanish Government through the QUANTUM ENIA project call -Quantum Spain project, and by the European Union through the Recovery, Transformation and Resilience Plan -NextGenerationEU within the framework of the Digital Spain 2026 Agenda.
V. APPENDIX

The algorithm
In this Appendix we describe the algorithms we employ to obtain the left and right tMPS encoding the contraction of the TN associated with the expectation value ⟨O(T )⟩, where O is a local (or a set of local) operator acting on one or few neighbouring sites.In practice, for the iterative methods we employ it is convenient to incorporate the operator(s) on one side of the network, as will become clear in the following.If we choose eg. to incorporate O on the right side of the network, our algorithms allow to build ⟨L I | and |R O ⟩ for a given time T .
The basic building blocks of the TN are the MPO tensors representing the trotterized time evolution (here we employ the parametrization from [30]).For a given T which determines the extension of the network in the temporal direction (N T = T /δt, δt being the Trotter step), we contract them in columns to generate the transfer matrices encoding the rotated (space-like) evolution E T for the folded system of which ⟨L| and |R⟩ are the leading left and right eigenvectors in the thermodynamic limit.Depending on what we want to construct, we can also include the local operator at the bottom of the transfer matrix, in that case we label it as E T O .We have implemented both a power-method for extracting the leading eigenvectors of E T for a fixed time T , as well as an algorithm based on the strict causal cone of the network when dealing with local operators.Both our algorithms are small variations of those presented in [4,10,12]), the main difference being the cost-function we use to perform the truncation required in the various iterations: we project the tMPS matrices on the support of the largest singular values of the respective RTM, as will be explained in the following.
The power method (see Fig. 5(a)) works by repeatedly applying the transfer matrix to an initial guess tMPS of length N T until convergence is reached.More specifically, we start from the left ⟨L I | and right |R I ⟩ vectors, and apply a column E T to the left and a E T O to the right.At each step, the bond dimension of the tMPS increases by a factor d 2 , where d is the physical dimension of the system constituents, so in order to proceed we truncate following our prescription and take the updated ⟨L I | as input for the next step.With our parametrization of the MPO tensors, E T is symmetric in left-right legs, so we can use ⟨L I | also as the new |R I ⟩.Alternatively, the optimization for |R I ⟩ simply requires an analogous step involving ⟨L O |.
In order to determine whether the power method has converged, we calculate several entropies associated with the ⟨L I | and |R O ⟩ vectors: most notably, the Von Neumann entropy S V N t (L I ) and the generalized Rényi 2, which we define as S r2 t (R O ) = − log n λ 2 n , where λ 2 n are the eigenvalues of the RTM T Convergence is reached at a given step i of the power method if, for the aforementioned entropies, one has ϵ, which in our case we take to be ϵ = 10 −6 .
The light-cone method works in a slightly different way, as it allows to build the tMPS for a time T starting from the one for T −1 with a single optimization (see Fig.As for the power method, after each iteration the bond dimension of the tMPS grows by a factor d 2 , so that truncating is required to avoid an exponential computational cost.The truncated ⟨L I | and |R I ⟩ are then used as starting point for the next iteration.The cone is thus built from the center moving outwards, assuming that the truncation of its edges does not spoil the causal structure induced by the local operator.The algorithm seems in any case to be quite stable, as we checked by piling up periodically a few layers of E T before truncating and verifying that the final result is the same as the one obtained when we truncate after each iteration.
In spite of their differences, both methods end up giving comparable results, as they share the same truncation procedure, which goes as follows: We focus on optimizing the overlap between ⟨L I | and |R O ⟩, which includes the operator O whose expectation value we are interested in.Since our aim is to find a low-rank approximation for the RTM T L I |R O t at any bi-partition t and T − t, we focus on computing the eigenvalues of the matrix (T depicted in Fig. 5(c), which can be associated with the squares of the relevant singular values.As usual in TN calculations, we can transform a global optimization problem into a local one by making use of gauge transformation [45,46].In particular, we start by bringing both L and R vectors individually to a standard orthogonal gauge, starting from the side of the initial state (Fig. 5(c)).This gauge transformation allows to incorporate the Λs L,R (t) factors towards the operator side, so that we don't need to compute them explicitly as their contraction up to t < T now reduces to an identity.We now only need to contract the bottom environments and project them over their largest D ′ singular values at each t, thus optimizing the overlap between the two tMPS.The advantage of keeping the local operators on either the left or the right side of the network is that, while we optimize with respect to ⟨O⟩, we always end up with a new ⟨L I | (|R I ⟩) which does not include the operator itself.This tMPS, which is related rather to the time evolution of the state, thus encodes the influence functional of the system, and can be used as starting point for the next iteration of our algorithms.

Estimate of the bond dimension for a given fidelity
In order to estimate the bond dimension required to obtain our tMPS with a given fidelity, we follow a procedure similar to the one employed in [19] to estimate the truncation error in our tMPS compared to an exact result, which in our case corresponds to the evolution obtained with the maximum bond dimension we can afford, D max (≃ 1000).The largest T we consider for this estimation are thus restricted to times for which D max does not induce any sizeable truncation error with respect to the exact (Trotterized) dynamics.O ⟩ |.Once more by fixing such an error to a given threshold ϵ we can identify the corresponding D * (t) for the tMPS.The typical behavior of the resulting curves for F can be seen in Fig. 6.

Comparison with other approaches
As previously mentioned, our optimization prescription is a relatively small modification with respect to already existing algorithms, which allows nevertheless to make a clearer connection with the underlying structures and estimate the computational complexity associated with the calculation of time-dependent expectation values.It might nevertheless be interesting to check how this method compares with other prescriptions employed in the literature.
The most straightforward optimization procedure when dealing with (t)MPS after applying a column of MPO is to truncate using the standard canonical forms for the left and right vectors separately.In our framework, with left-right symmetry of the Hamiltonian and  translational invariance, this would amount to an optimization of ⟨L| with its conjugate |L⟩ instead of |R⟩.The operation of complex conjugation here can be seen as a partial transpose involving the forwards and backwards leg of the tMPS (basically an insertion of a series of swap operators), see Fig. 7 (a), since (not considering additional complications in case the initial state is complex) it basically amounts to an exchange U (t) ↔ U † (t) of the time evolution operators.Due to the non-trivial structure induced by this, the projector structure of the TN is lost even if no operators are present, so that one obtains a non-trivial bond dimension even when O = I.
Instead of optimizing with respect to a single operator inserted at the edge of the MPO column, the conjugation would imply acting on the whole column, ie. to perform these transpositions at each timestep.
Another possible strategy, which has been suggested in [6,10] involves bringing the tMPS to canonical form starting from the bottom side and then doing a sweep from the initial state optimizing ⟨L I |R I ⟩, Fig. 7 (b).In this case, while the spectrum of the RTM is trivial since no operator is present, the singular value decomposition is not: this can be seen again as due to an insertion of swap operators which generate a non-trivial spectrum, which however cannot be directly related to the causal structure generated by a local operator.In this case, one is unable to give a priori a statement on the maximum rank required by the algorithm.
In spite of this difference, remarkably this latter method returns comparable bond dimensions to the ones obtained with our method when O = σ x , see Fig. 8. On the other hand, the optimization using the standard canonical forms requires a larger bond dimension, although the behavior of the various algorithms is comparable: this confirms our expectation that the scaling of our truncation based on the overlap ⟨L|R⟩ is the same as the one based on ⟨L|L⟩.8. Bond dimensions of the tMPS as function of time obtained with different methods (as described in the text): our method (blue line), the optimization starting from the side of the initial state (orange) and the canonical form (green), imposing the same truncation error (ϵtrunc = 10 −6 ).Our method leads to the smallest bond dimension.

FIG. 1 .
FIG. 1.(a) The time-dependent expectation value of a local operator acting on a one-dimensional system in the Keldysh representation (b) can be described by a double sheet 2D TN contraction.(c) Upon using a Trotter approximation and the locality of the operator, it simplifies to a triangular TN.

FIG. 2 .
FIG. 2. (a): The triangular TN can be contracted from the sides, identifying a left ⟨LO| and a right |RQ⟩ tMPS.The upper-most tensors of the temporal MPS are dictated by the initial state, while the lower-most ones by the choice of the operator.(b) The reduced transition matrix.(c) In the absence of operators, the folded tensors of the time evolution resolve to identities.(d) The partial transpose of the time-evolved left-right density matrix of the system, where forwards and backwards legs are swapped.
we iteratively construct the left and right vectors ⟨L O (T )| and |R Q (T )⟩ starting from ⟨L O (T − 1)| and |R Q (T − 1)⟩ by absorbing a new column of MPOs into them.These MPOs are obtained by contracting a column of the original TN and they thus represent transfer matrices E T evolving states of a time-slice in space, a rotated version of the standard evolution in time.The bond dimension of the corresponding MPS increases, and we then keep it under control by projecting onto the largest D singular values of T L O |R Q t keep the distance in norm || ⟨L O | − ⟨L D O | || 2 constant, where ⟨L D O | is the truncation of the tMPS ⟨L O | to a given bond dimension D. Using the reasoning in Ref. [39], the term to consider would then be the overlap ⟨L O |L D O ⟩, which is upper bounded by the sum of the norm of the residuals discarded in the truncation process.So far we have just provided a bound on the overlap ⟨L O |R O ⟩ which, even for reflection-invariant systems and considering the symmetric case O = Q, differs from ⟨L O |L D O ⟩ by the absence of a complex conjugation.However, the two quantities are related, since ⟨L O |L D O ⟩ can be obtained from ⟨L O |R O ⟩ by evolving |R O ⟩ with one layer of local unitaries, representing the partial transpose

FIG. 3 .
FIG. 3.Left: Bond dimensions D obtained by imposing a truncation error of 10 −4 in the singular value decomposition of the RTMs, as function of time (see main text).The shaded areas are delimited by the minimum and maximum D we found by varying the initial state.The thick solid lines represent the bond dimension curves for the operator entanglement (OE): they lie consistently above the corresponding ones for the tMPS.Blue colors denotes the integrable case (int), orange-red the non-integrable one (non-int).
5(b)): Starting from ⟨L I | and |R I ⟩ of a given length N T , we apply a transfer matrix E T +1 of length N T + 1 to the left and a E T O to the right, extending the network both in the time and space direction.

FIG. 5 .
FIG. 5. Illustrations of the methods used for building the tMPS.(a) Power method (b) Light cone method (c) The matrix involved in our low-rank approximation of T and the role of the gauge transformation in our calculation.

FIG. 6 .
FIG. 6. Fidelity for overlap ⟨L Dmax O |L D O ⟩ for different bond dimensions, as function of time, for the integrable case with starting state |+⟩.After fixing a threshold ϵ we can extract the required bond dimension to give a faithful approximation of |L D O ⟩.
(a) A possible method for truncating the tMPS relies on the optimization of the overlap ⟨L|L⟩, where in the symmetric case |L⟩ can be related by a partial transpose to the vector |R⟩.(b) Alternatively, it is also possible to consider the RTM starting from the side of the initial state, a procedure which can give a nontrivial SVD spectrum even in the absence of a local operator.
FIG.8.Bond dimensions of the tMPS as function of time obtained with different methods (as described in the text): our method (blue line), the optimization starting from the side of the initial state (orange) and the canonical form (green), imposing the same truncation error (ϵtrunc = 10 −6 ).Our method leads to the smallest bond dimension.
If we work with normalized MPS, the error induced by truncating the operator MPS to bond dimension D is encoded in the fidelity F = | ⟨ψ Dmax O (t)|ψD O (t)⟩ |, where |ψ D O (t)⟩ represents the time-evolved state up to time t truncated at every time