Light cone tensor network and time evolution

The transverse folding algorithm [Phys. Rev. Lett. 102, 240603] is a tensor network method to compute time-dependent local observables in out-of-equilibrium quantum spin chains that can sometimes overcome the limitations of matrix product states. We present a contraction strategy that makes use of the exact light cone structure of the tensor network representing the observables. The strategy can be combined with the hybrid truncation proposed for global quenches in [Phys. Rev. A 91, 032306], which significantly improves the efficiency of the method. We demonstrate the performance of this transverse light cone contraction also for transport coefficients, and discuss how it can be extended to other dynamical quantities.


INTRODUCTION
Tensor networks (TN) [1][2][3] have gained in the last decade a prominent role among numerical methods for quantum many-body systems. Simulating the dynamics of out of equilibrium systems remains nevertheless one of the most challenging open problems for these (and other) techniques.
In one dimensional systems, limitations of TN methods for dynamics are well understood: in global quenches the entanglement may grow fast [4][5][6], and the true state can escape the descriptive power of the TN ansatz. This so-called entanglement barrier limits the applicability of the matrix product state (MPS) [7][8][9][10] description, and makes it difficult to predict the asymptotic long-time behavior, even when local observables in this limit are expected to be well-described by a thermodynamic ensemble, itself well approximated by a matrix product operator (MPO) [11][12][13][14][15][16]. A number of methods have been suggested to try to overcome this issue and extract information about the long-time behavior of local properties [17][18][19][20][21][22][23][24][25][26][27]. While there is no universal solution, understanding the entanglement structures in the evolution TN can be crucial to identify the most adequate one for practical computations.
In particular, the transverse folding strategy [18,28,29] avoids the explicit representation of the evolved state as a MPS and instead focuses on contracting a TN that represents exactly (up to Trotter errors) the timedependent observables. Instead of the standard evolution in time direction, the folding algorithm contracts the TN along space. In some scenarios, this allows local observables to be computed to longer times than other approaches [30], and it is an exact strategy for certain models [31]. Recently, there has been a rekindled interest in this approach, triggered by the interpretation of the network in terms of an influence functional [32][33][34].
In local lattice models, the velocity of propagation of information is upper-bounded [35][36][37] and the exact TN for observables has a light cone structure. While there have been proposals that exploit this fact to reduce the cost of the numerical simulation of the evolved state with TN [38][39][40][41][42][43], and with quantum simulation [44], until now, the potential of combining it with the transverse strategy has not been explored.
Here we propose a strategy to exploit this property, a transverse light cone contraction of the TN (TLCC). As in the original transverse folding, the TLCC does not directly suffer from the entanglement growth in the state, and will be more efficient than standard algorithms when entanglement in the time direction grows slower than in the spatial one. But the TLCC improves the efficiency with respect to the transverse folding in all cases, by reducing the computational effort to that of approximating the minimal network describing the time-dependent ob-arXiv:2201.08402v3 [quant-ph] 6 Dec 2022 servables in a Trotterized evolution. We demonstrate explicitly its performance for global quenches and differenttime thermal correlators at infinite temperature, and investigate how the strategy can make use of the (more efficient) physical light cone determined by the Lieb-Robinson velocity [36]. We discuss possible extensions to other interesting quantities.

LIGHT CONE TENSOR NETWORK FOR GLOBAL QUENCHES
The one-dimensional global quench is a natural test bench for time-evolution TN algorithms. At time t = 0 the system is prepared in a state that can be written as a MPS (e.g. a product state), and then it is let to evolve under a fixed Hamiltonian. For simplicity, we restrict the discussion to a nearest-neighbour model, and a translationally invariant case, but the construction generalizes straightforwardly to any model with local (finite-range) interactions and some non-translationally invariant scenarios.
The transverse folding proposal of [18] starts from a two-dimensional TN whose contraction represents some time-dependent observable, such as a local expectation value. This TN can be constructed from a Suzuki-Trotter approximation of the evolution operator, where the evolution for a discrete step of time δ can be approximated as a matrix product operator (MPO) [11,12] with a small bond dimension, constructed from a product of two-body gates [13]. The TN for the observable at time t = M δ is obtained by applying M copies of this MPO with the initial state, which yields the evolved state, and contracting the operator of interest between this and its adjoint.
While standard TN algorithms as TEBD or tMPS [11,[47][48][49] compute the observable by contracting the network in the time direction, the transverse folding strategy performs the contraction in the spatial direction, after folding the TN in half, such that tensors for the same site and time step in the ket and the bra are grouped together (see figure 1a). After folding, the growth of entanglement in the time direction can be slower than in the spatial one, with the most dramatic difference observed for integrable systems [28,46], but occurring also in generic cases, as the ones shown here. When this difference in growth is present, the transverse strategy allows reaching longer times than standard algorithms.
For a translationally invariant system in the thermodynamic limit, the transverse contraction reduces to an expectation value of the form (L(t)|E O (t)|R(t)), where (L(t)| and |R(t)) are the dominant left and right eigenvectors of the transfer operator and represents the concatenated [50] local tensor of the timedependent state, itself a MPO. In the transverse folding strategy, the boundary vectors (L(t)| and |R(t)) are approximated by MPS. This approximation can be found, for instance, via a power iteration or a Lanczos algorithm, using repeated MPO-MPS contractions. Such strategies do not take into account that the TN has a light cone structure. Because the individual gates are local, outside the causal cone of the operator, each gate cancels with its adjoint. This ensures that each of the required boundary vectors (dominant eigenvectors of the transfer operator) corresponds precisely to the contraction of a triangular network as depicted in fig. 1. We can approximate directly the contraction of such triangle in the space direction by a MPS. This strategy, which we call transverse light cone contraction (TLCC), allows us to obtain (L(t)| and |R(t)) in a fixed number of steps (proportional to M ). Furthermore, once we have found the vectors for M time steps, we can directly obtain them for M + 1 by applying a single MPO (as illustrated in the figure), which increases the length by one, and approximating the result via a single truncation step. This step can be performed using standard MPS truncation algorithms, which reduce the bond dimension by minimizing a distance between the truncated vector and the original one. However, for this particular problem the hybrid truncation algorithm proposed in [29], which effectively evolves the bond of the boundary vector according to the real time dynamics, yields a much more efficient use of the available bond dimension (see also insets of fig. 2).
The TLCC strategy results in a more efficient algorithm than the originally proposed folding, which required iterative MPO-MPS contractions until convergence of the dominant eigenvectors, run independently for each different time step (in particular, for the cases analyzed in this work, we find the power iteration required several tens of MPO-MPS contractions per time step). Notice, nevertheless, that if the bond dimension used is large enough, both the original folding algorithm and the TLCC should result in the same boundary vector. What ultimately determines the applicability of transverse strategies is thus the amount of entanglement present in the transverse network.
To probe the performance of the method, we consider a quantum Ising chain, initialized in a product state We then apply the Hamiltonian, and compute local expectation values after time evolution. In all the following we fix J = 1, and a Trotter step δ = 0.1, and vary the parameters of the model to study integrable (g = {0.5, 1}, h = 0) and non-integrable (g = −1.05, h = 0.5) regimes. Figure 2 shows the results and demonstrates that the TLCC can efficiently simulate the integrable quenches. In the non-integrable regime, the required bond dimension grows much faster with The TLCC contraction has been obtained both with the standard MPS truncation (green squares) and the hybrid truncation of [29] (dark blue circles). For comparison, we also show the results of standard iTEBD (blue diamonds) and Heisenberg picture DMRG (purple triangles). a For the integrable case (a,b) we show also the analytic result (black line). The insets show the scaling of the bond dimension required to keep constant precision in each algorithm [45]. In the integrable case, this is compatible (at the later times) with a polynomial increase D ∼ t α , consistent with observations in [28,46]. In the non-integrable case, the increase is compatible with an exponential growth for both truncation methods, but using the hybrid truncation exhibits a slower rate than standard ones, such that longer times can be reached with the same bond dimension. The left inset in (c) shows a zoom of the main plot to better appreciate the differences.
a Heisenberg picture results are only shown in (c), since, in the integrable case, the operator in (a,b) can be exactly written as an MPO with constant bond dimension at all times [17].
time, but the method is still advantageous as compared to standard evolution, much more so when the truncation is performed as in [29] (see right inset of fig. 2c).

LIGHT CONE TENSOR NETWORK FOR TRANSPORT COEFFICIENTS
The same idea can be adapted to the computation of other dynamical quantities. It is the case of thermal correlators, of the form is the time-evolved operator in Heisenberg picture. Since [ρ β , H] = 0, the thermal state is invariant under the evolution, and using ρ β ∝ ρ β/2 ρ † β/2 we can write (up to normalization), 1 ρ β/2 ). Using a MPO approximation to ρ β/2 (obtained with standard TN methods [11,12,51,52]), and the Trotterized real time evolution as in the previous section, this quantity can be expressed as a two dimensional folded TN, which can be contracted in the temporal [53][54][55] or spatial (transverse) [28] direction.
Due to the invariance of the thermal state, each local observable generates also a light cone structure that can be exploited in the TLCC approach. Now the cancellation of gates outside the causal cone of the operators occurs both at the upper and the lower parts of the network (see figure 3a), and the minimal TN has a rectangular form, resembling a pillow, a structure which was used in [56] to evaluate correlators in random quantum circuits. The TLCC strategy again requires contracting a triangular TN corresponding to the lateral corners of the figure to obtain boundary vectors (L β (t)| and |R β (t)). [57] If both operators act on the same site ( = 0), the time dependent correlators can be expressed as a contraction (L β (t)|T β,O1,O2 (t)|R β (t)), with a single MPO T β,O1,O2 (t) constructed from concatenating the local tensors for the unitaries, the operators and the states (see fig. 3a). For correlators at non-zero distance the minimal TN becomes elongated ( fig. 3a, lower diagrams). To approximate its contraction, the boundary vectors (L β (t)| and |R β (t)) for a certain time t are first grown to incorporate, respectively, O 1 at the bottom of the TN, and O 2 at the top. These extended vectors contain the evolution steps up to time t + 2δ, and can be contracted together to obtain the correlators at = 1 for times t+3δ and t + 4δ. The vectors can be then evolved again, following the TN structure, which does not increase their length, but allows access to correlators at any later time t + (2 + k)δ and distances = k, k + 1. Applying this systematically we can obtain all non-vanishing correlators. This generalizes trivially to operators on more than one site, or with MPO structure.
Here we illustrate the simplest case, infinite temperature, where ρ β=0 ∝ 1 and the contour of the TN becomes uncorrelated. We consider the energy density operator which can be written as a MPO of range 2. Figure 3b shows our results for the correlators C EE (t, , β = 0) as a function of time for several distances in the nonintegrable (g = −1.05, h = 0.5, main plot) and integrable (g = 0.5, h = 0, inset) cases [45].
Specially interesting is the possibility of ab initio calculations of transport properties [58] in non-integrable models. In particular, diffusion constants can be related to the spatial spreading in time of autocorrelations of a density [22,59,60]. Normalizing the correlators asC EE (0, ) := C EE (t, )/ C EE (0, ), a diffusion constant D(t) may be obtained from their spatial variance [59] as ∂W 2 ∂t = 2D(t). Figure 3c shows the (linearly growing) variance W 2 (t) (main plot), and the corresponding diffusion constant (inset) obtained from the correlators for the non-integrable case. The diffusion constant is well fitted by a function D(t) = D E exp(b/t), compatible with saturation to a constant D E ≈ 1.9 in the asymptotic regime [61]. While TEBD (blue diamonds) produces close values for the same quantities, the error is appreciable in the diffusion constant already at short times.

THE PHYSICAL LIGHT CONE
In general, we expect that the physical light cone is much narrower than the trivial one from the Trotterization, used in the previous sections. We could thus approximate the TN by a light cone one in which the slope corresponds to the maximal physical velocity v LR . This can be achieved by implementing a more efficient TLCC growing iteration, in which n T = 1/(v LR δ) time steps are applied at once every time a space site is contracted ( fig. 4a). Notice that this light cone is not exact, but has (exponential) corrections. Thus it is convenient to consider the light cone for a subsystem of size L eff that includes the support of the operator. [62] To probe this reduced light cone we choose an integrable instance, (g = 0.5, h = 0), for which the Lieb-Robinson velocity is known (v LR = 1, corresponding to n T = 10 with our Trotter step), and simulate the global quench of fig. 2a. Compared to TLCC for the full light cone with the same bond dimension, we observe ( fig. 4) that the physical one, determined by v LR , captures indeed the correct evolution: while the narrower light cone deviates from full results, the errors are reduced exponentially (until the level of original truncation error) by considering a small window L eff .

DISCUSSION
We have presented a strategy that builds on the transverse folding [18] to approximate time-dependent observables in a one-dimensional quantum system. Noticing the exact light cone structure of the TN and implementing its transverse contraction, it is possible to compute long time properties in a more efficient manner. Combined with the hybrid truncation [29], this allows us to reach longer times with a smaller bond dimension whenever the temporal entanglement grows slower than the physical one, which, as we have seen, happens not only for integrable systems. It is possible to use the physical upper bound of the Lieb-Robinson velocity to further restrict the width of the relevant TN and define a more efficient iteration.
We have evaluated the performance of the TLCC strategy for integrable and non-integrable global quenches, and for transport properties at infinite temperature. With minimal changes, the method extends to other scenarios, such as finite temperature or non translationally invariant setups including impurities or a contact between two chains. It is furthermore possible to adapt the strategy to other more complex dynamical quantities.
The basic TLCC does not require additional hypothesis to truncate observables or states. Its convergence can be systematically explored as the bond dimension is increased. What ultimately limits the validity of the strategy is the entanglement in the time direction, which strongly depends on the setup and the model [28,46,63]. The behavior of the TLCC can thus provide useful information to determine optimal strategies for different problems. Another parameter in the approximation is the Trotter step, which is known to affect the entanglement growth in standard algorithms [49]. Since simulations with different δ may be necessary to extrapolate the exact results, it is also interesting to study how varying δ affects our observations. Further interesting avenues for future investigation are exploring the TN cut according to different velocities, to explore the propagation of correlations in the TN and effectively measure v LR .
While we were completing this manuscript, an equivalent strategy for global quenches was independently suggested in [64].
We are thankful to J. I. Cirac, M. Hastings and L.  This section shows explicitly how to construct the TN of Fig. 3(a) in the main text, for the case of arbitrary inverse temperature β. We are interested in correlators of the form 1 ), where ρ β = e −βH /Z is the thermal equilibrium state at inverse temperature β, Z = tr(e −βH ) is the partition function, O [ ] k is a (local) operator acting on site , and U (t) = e −iHt is the time-evolution operator for time t. In order to write this quantity as the contraction of the TN shown in the text, we start by finding an MPO approximation of the thermal state. This can be achieved with several algorithms [11,12,51,52]. Here we illustrate (see Fig. 5) the (possibly) most common algorithm, based on a purification and a Trotter expansion of ρ β .  7. TN for the thermal correlators at infinite temperature. In the particular case β = 0, the thermal MPO is exact and proportional to identity, and local unitaries cancel around local O1, resulting in a diamond shape. In the general case, the thermal MPO is not exact, and the diamond-shaped TN will approximate the full one.
The purification approach is equivalent to considering a thermofield double state, i.e. a pure state of the form where |Φ is a maximally entangled state of the system and an ancillary copy of it. Tracing out the ancillary system results (up to normalization) in the Gibbs ensemble ρ β . For any basis {|n } of the system, we can write |Ψ β ∝ e −βH/2 n |n, n .
The most frequently used TNS algorithm for thermal equilibrium states proceeds by approximating |Ψ β by an MPS in a basis in which each system site is grouped with an ancillary one (forming effective sites of dimension d 2 ). To do so, the state is initialized to the maximally mixed one between system and ancilla (equivalent to the vectorized identity operator), i.e. a MPS with bond dimension one. The exponential operator e −βH/2 can be discretized as the product of a finite number M = β/(2δτ ) of imaginary time steps of length δτ . If the Hamiltonian is local, each of them, can be approximated by a MPO (for instance, for nearest-neighbor models, one can use an evenodd Trotter approximation) and successively applied on the MPS, acting on the system degrees of freedom. After each step, a standard truncation can be performed (using any of the algorithms for Trotterized time evolution), such that after a fixed number of steps M = β 2δτ , an MPS approximation to (5) is obtained. As sketched FIG. 9. Alternative construction of the TN for the thermal correlators using the purification structure.
in Fig. 5(a), this procedure is equivalent to approximating the exponential by an MPO (with the Frobenius norm characterizing the quality of the approximation in the operator level). Whereas this procedure could be run for the full inverse temperature, to obtain an MPO approximation of e −βH , the truncation in MPO-MPS products does not preserve positivity. Instead, an MPS approximation of the thermofield state results necessarily in a positive density matrix (with purification structure) when tracing the ancillas, as shown schematically in Fig. 5(b).
The TN for the thermal response functions can be constructed applying the operator O 1 followed by the Trotterized time evolution on this MPO, and finally applying O 2 (possibly on a different site) and taking the trace. This results in a TN that is periodic in the time direction. Folding (or flattening) it results in a doubled TN, similar to the one obtained for pure state evolution, as illustrated in Fig. 6.
The simplest case is that of infinite temperature (β = 0), when the thermal state has an exact MPO representation with bond dimension one, since ρ β=0 ∝ 1. Then the local unitary matrices that represent the real time evolution cancel also around O 1 exactly, and the TN to be contracted has a diamond shape (see fig. 7).
At arbitrary temperature, the cancellation around O 1 is no longer exact, since the MPO representation of the thermal state is only approximate, and local unitaries do not commute exactly with it. Thus we can consider the diamond-shaped TN in this case to be an approximation of the infinite one (see Fig. 9). We expect this to introduce a small error. Notice that in more standard evolution algorithms (i.e. those which evolve the MPS in real time), exploiting this light cone structure has also been shown to be useful in infinite systems, for instance by considering an expanding window embedded in an infinite MPS [41][42][43].
Finally, notice that the TN construction described above is not unique. If we make use of the property which has been previously exploited for the simulation of real time evolution of thermal states [53][54][55], and apply also the cyclic property of the trace, we can write the quantity of interest as (up to a normalization factor) This results in a different TN, illustrated in Fig. 9, in which the upper and lower boundaries are given by the purification tensors. For infinite temperature, both constructions are equivalent, but for the general case, they will give rise to different approximation errors. A systematic analysis of these alternatives, the approximation error and its dependence on temperature will be carried out elsewhere.

ERROR ESTIMATES FOR THE DIFFERENT APPROACHES
In the insets of Figure 2 of the main text, we show the scaling of the bond dimension needed with time to maintain fixed precision in the global quench scenario. Here in this appendix, we provide some extra information on how the scaling is computed and what we mean by fixed precision.
A standard way to bound the error in TN simulations is to sum the squares of the discarded Schmidt weights in every truncation [49]. For the results from TEBD and Heisenberg-picture DMRG, we use that as a measure of the growth in the errors during the simulation. Notice that in the case of TEBD, as we are simulating an infinite system that measure is only an heuristic, as it does not provide an upper bound in the error that one can incur when evaluating some observable. For the TLCC, both when we use the standard and the hybrid truncation from [29] to truncate the boundary vectors, we keep track of the deviation of the expectation value of the identity. As explained in the main text, knowledge of the boundary vectors (L(t)| and |R(t)) gives access to out-of-equilibrium expectation values by computing the expectation value (L(t)|E O (t)|R(t)). That includes the possibility of computing the expectation value of the identity, which should be one for any properly normalized state. The deviation with respect to this value, when using MPS approximations for (L(t)| and |R(t)), gives a good measure of the error of the TLCC.
In order to perform the scaling analysis, we compute as a function of time the quantities mentioned above for the different TN algorithms with simulations with different bond dimensions. Setting a constant value of the precision, that is, of the truncation errors in the TEBD and Heisenberg-picture DMRG cases and of the deviation from the identity for the TLCC, we can keep track of the times when the different bond dimensions exceed the desired precision threshold.

ALTERNATIVE VALUES OF THE PARAMETERS
In the main text, for simplicity, we focused in a particular point in parameter space in the non-integrable case. Here, we show results obtained in a different case, described by the couplings (J = 1, g = 1.4, h = 0.9) and studied in [22,60]. The spatial variance and diffusion constant obtained with TLCC for this case are shown in Fig. 11, along with results obtained with TEBD. As seen in the main text, the introduction of TLCC allows in this case as well to extend the range where it is possible to reliably simulate out-of-equilibrium dynamics by a factor around two.