Constructing Tensor Network Influence Functionals for General Quantum Dynamics

We describe an iterative formalism to compute influence functionals that describe the general quantum dynamics of a subsystem beyond the assumption of linear coupling to a quadratic bath. We use a space-time tensor network representation of the influence functional and investigate its approximability in terms of the bond dimensions and time-like entanglement in the tensor network description. We study two numerical models, the spin-boson model and a model of interacting hard-core bosons in a 1D harmonic trap. We find that the influence functional and the intermediates involved in its construction can be efficiently approximated by low bond dimension tensor networks in certain dynamical regimes, which allows the quantum dynamics to be accurately computed for longer times than with direct time evolution methods. However, as one iteratively integrates out the bath, the correlations in the influence functional can first increase before decreasing, indicating that the final compressibility of the influence functional is achieved via non-trivial cancellation.


I. INTRODUCTION
Obtaining the long-time dynamics of a large quantum system is in general intractable due to the exponential scaling of Hilbert space with respect to system size and the associated exponential growth of spatial entanglement with time. Fortunately, in many cases, one is most interested in the dynamics of observables defined on a small subset of the full system, and can thus reframe the dynamics of the observable from the viewpoint of a subsystem coupled to a bath 1 . If one is able to determine the dynamics of the bath and its influence on the subsystem, then the original dynamics problem is reduced to obtaining the dynamics of the subsystem [2][3][4][5] .
The influence functional (IF) 2 method provides an exact framework for computing the dynamics of an arbitrary bath and its interactions with the subsystem. However, the cost of computing the IF without approximation is comparable to determining the dynamics of the full system and the size of the IF scales exponentially with the number of time steps. Thus it is not usually possible to use the IF method without additional approximations.
The IF can be viewed as reweighting the path integral of the subsystem. In the case of linear coupling to a harmonic bath, Feynman and Vernon derived an analytical form of the IF 6,7 , which takes the form of the Boltzmann weight of a complex valued Hamiltonian defined in the time direction with pairwise interactions between time points.
For many physical bath spectral densities, it is natural to assume that the pairwise time interaction is shorta) Electronic mail: eyerik.a@gmail.com b) Electronic mail: gkc1000@gmail.com ranged in time, corresponding to a finite "memory" in the influence of the bath, and many numerical approximations have successfully taken advantage of this shortrange temporal nature [8][9][10][11][12][13][14][15][16][17][18][19][20][21] . For IF methods, the assumption of limited memory allows one to remove the exponential growth of cost of the quantum dynamics with simulation time, thus making long time-scale quantum dynamical simulations possible. For example, in QuAPI 8 , one approximates the analytical harmonic bath IF by only including terms acting on time steps within a finite time window. Alternatively one can construct an ansatz for the IF; a natural choice is a matrix product state (MPS) in the time direction, which compactly encodes shortrange time correlations. This approach has been used in a number of recent works, which, although they do not necessarily use the language of IFs, all proceed by constructing a compressed version of the IF 22 or a closely related object such as the auxiliary density operator as defined in QuAPI 23,24 , the process tensor 25,26 , or other variants [27][28][29][30] .
In this work, we are interested in using the IF method for computing the dynamics of a subsystem within a general quantum system. Such a subsystem may arise as part of a larger interacting problem (in which case, the subsystem might not be different from other parts of the system) or it may arise from a system-bath model. In either case, the couplings and bath cannot be assumed to be linear and quadratic respectively and thus the analytical form of the IF is not known. Instead, the IF is simply a particular integral of the space-time dynamics that must be obtained numerically. To do this concretely, we can use a tensor network description of the space-time dynamics. For a 1-dimensional representation of the system and bath, the tensor network is thus defined in 1+1 dimensions. Traditional time evolution corresponds to contracting this network first along the time axis, and if done exactly leads to exponential cost with time. Computing the IF corresponds to the contraction of the network in the spatial direction to yield a final matrix product state IF defined along the time direction, and if done exactly has an exponential cost with the spatial length.
A similar framework appears in the recently proposed modular path integral (MPI) procedure [31][32][33][34][35][36] . However, to avoid exponential computational cost, the contractions must be performed approximately. There has been some discussion regarding using filtering methods to reduce the cost of the MPI algorithm 31 . In contrast, we utilize tensor network compression and propose a tensor network contraction procedure similar to that presented by Banuls et al. 27,28 and Lerose et al. 22 , assuming a translationally invariant infinite problem. However, we will describe a general tensor network procedure to construct the IF without such assumptions and explore the numerical feasibility of doing so to compute quantum dynamics in different regimes beyond the commonly considered linear coupling to quadratic baths.
The paper is organized as follows. We first translate the IF into a space-time tensor network language and describe an iterative algorithm to compute it. We then investigate the compressibility of the IF and its ability to produce long-time dynamics, first for the canonical spinboson model where the analytical IF is known, and then an interacting hard-core boson model where there is no analytical expression, which corresponds to the case of general quantum dynamics. We analyze the time-like entanglement both in the IF itself as well as the intermediates that arise as the bath is numerically integrated out. We end with a brief discussion of the implications of this work for future studies.

II. THEORY
A. Definition of the Influence Functional To introduce notation, we first recall the definition of the influence functional (IF) 2 . The influence functional describes how the path integral of a subsystem is reweighted, under the influence of dynamical coupling to a bath. To obtain an explicit form, we define a full system composed of the subsystem of interest and the coupled bath. At time , we denote the subsystem density matrix by ( ), where is a basis for the density matrix, and the bath density matrix is analogously written as ( ). The basis of the full system is spanned by the product space { }⊗{ }. The evolution of the density matrix is given by a linear operator, the Liouville operator , which we partition as = + where contains the component operating only on the subsystem and contains the component on the bath and interactions between the subsystem and bath. If we further assume the system dynamics obeys Hamiltonian evolution, then the Liouville action can be written as ⋅ = [ , ⋅].
Formally, ( ) is obtained by time evolving the entire system and tracing out the bath degrees of freedom. The path integral expression, assuming a second-order Trotter decomposition of the time evolution operator into timesteps of length , is where ( 0 , 0 ) is the initial state of the system, and the double bra/ket notation indicates we are working in Liouville space, with the density matrix being a vector in this space. For simplicity, we assume there are no correlations between the subsystem and bath initially such that Furthermore, is typically assumed to be diagonal in the basis { } (we lift both these restrictions below). Then Eq. (1) becomes where ( 1 , 2 , ..., ) is the influence functional The IF assigns a complex weight to each configuration of the system path integral. Consequently, the storage of the IF grows exponentially with number of time steps .

Influence Functional Structure
Translating Eq. (1) and Eq. (2) into the tensor network language is straightforward, and is shown in diagrams Fig. 1(a) and (b), respectively. The time evolution operators of the system and system-bath dynamics respectively appear as boxes with two and four legs in Fig. 1(a), with the legs labelled by the bra and ket basis states. The elements in the tensors are ⟨⟨ | − | ′ ⟩⟩ and ⟨⟨ , | − | ′ , ′ ⟩⟩, respectively. In the case that is diagonal with respect to the subsystem basis states { }, then the elements of the system-bath time evolution operator are ⟨⟨ , | − | , ′ ⟩⟩ and can be depicted as a box with three legs labelled by , , and ′ , as in Fig. 1(b). The final trace operation over the bath degrees of freedom at the last time step can be written as a vector that is contracted with the corresponding leg of the time evolution operator tensor.
The influence functional element ( 1 , 2 , … , ) is the object within the blue rectangle of Fig. 1(b), obtained after performing tensor contractions over all bath tensors as denoted by the connected lines. Within this diagrammatic picture it is easy to depict the generalization of the influence functional to a correlated initial state. In this case, the dotted line indicates a correlated initial state with entanglement between the subsystem and bath, and the influence functional is defined with an additional index, ( 1 , 2 , … , ; ). Similarly, if cannot be diagonalized in the subsystem basis { } then we can generalize the influence functional to contain two subsystem indices at each intermediate time, , and it is the object within the blue rectangle of Fig. 1(a). Given the influence functional, arbitrary time-correlation functions can be computed as shown in Fig. 1(c).
Because the influence functional has a onedimensional structure along the time axis, it is natural to rewrite it as a matrix product state of tensors (see Fig. 1 where ( ) denotes a matrix of dimension × for each element of the basis , except for (1) ( 1 ) and ( ) ( ) which are dimensional row and column vectors respectively. In MPS language, is referred to as the virtual bond dimension. Although any IF can be represented as a MPS for sufficiently large , the MPS of small bond dimension naturally capture sums of exponentially decaying time-correlations along the time axis. The key system-specific questions to understand are thus (i) is the IF itself representable by an MPS of low bond dimension, in physically relevant dynamical and interaction regimes, and (ii) can the IF be constructed with manageable cost in those regimes. It is important to note that an affirmative answer to (i) does not imply an affirmative answer to (ii).

Space-Time Tensor Network Representation
To define an approximate procedure to construct the influence functional for complex bath dynamics, we first write down a space-time representation of the full system dynamics. We first assume that the bath Hilbert space is a product space over modes, We can then formally express the system density matrix at any time as a matrix product state (MPS) with a bond dimension denoted . Similarly, the Liouville evolution operator can be written as a matrix product operator (MPO) with the elements and has a bond dimension . Note that since the Liouville operator is assumed time-independent, is fixed. At last the time step, the bath degrees of freedom are traced out so the bath sites are now in MPS form, wherẽ = Tr ( ) with as defined in Eq. (8). If the Hamiltonian consists of only nearest neighbor interactions, one can obtain the operator using a Trotter-Suzuki decomposition of nearest neighbour gates and then directly map them onto a matrix product operator. Otherwise, for more general interactions one can use a 4th-order Runge-Kutta expansion 37 . In this case, the matrix product operator can have large bond dimension, but can be compressed by allowing for truncation errors of ( 5 ). In the cases studied here, the subsystem is small enough such that − can be obtained exactly. The full time evolution of the system with bath modes and time steps thus corresponds to the twodimensional tensor network diagram shown in Fig. 2(a). Correspondingly, the space-time representation of the influence functional is shown in Fig. 2

Transverse Contraction Scheme
The most common way to contract a 2D space-time tensor network is from bottom to top, i.e. in the direction of increasing time [38][39][40][41][42][43][44][45] . We refer to this as direct time evolution. For example, contracting the network in Fig. 2(a) row by row yields the system density matrix at each time step as an MPS. The cost of contracting two rows is ( 2 2 2 ) where is the dimension of the density matrices represented at each MPS site, and is the bond dimension of the density matrix MPS. For exact time evolution grows by a factor of at each time step. This means that in the worst case, grows exponentially with time.
However, if the time-correlations in the influence functional decay with long time, then this implies that the influence functional ultimately can be represented by a matrix product with low bond dimension along the time axis. This suggests that a more efficient contraction strategy is to contract column wise (in the transverse direction to time). The cost of contracting columns together in the process of construct the final influence functional is ( 2 2 2 ) where is the bond dimension of the column bath MPS (defined along the time direction). For exact contraction, will grow by at each contraction step and thus in the worst case scales exponentially with . But, it need not have a dependence on the total simulation time.
In practice, due to the exponential growth of the bond dimension, exact contraction of the 2D tensor network (in either direction) is often too expensive. In Banuls et al 27,28 , the explicit transverse contraction of the 2D tensor network was avoided by assuming that the system is infinite and translationally invariant, in which case the result of the infinite contraction of columns is proportional to the maximal eigenvector of the column transfer operator.
Alternatively, one may use standard matrix product state techniques to compress the intermediates that arise during the contraction to restrict bond dimensions or to some constant value [38][39][40] . In this case, the cost of the algorithm is dominated by the cost of the MPS compression (which requires performing a series of singular initialize → contract + canonicalize at each time step are the rows of the grid and are each represented as an MPO, and the subsystem of interest is at the left-most site. Before contraction, we first canonicalize each row into left canonical form as indicated by the right pointing triangles along the rows. The rightmost two columns are then contracted and compressed to fixed bond dimension using the standard MPS compression algorithm, where the column is first converted into a canonical form (here, top canonical form) and then compressed by singular value decomposition in the reverse direction (leaving it in bottom canonical form). The canonical form implies that the tensors satisfy an isometric condition (see diagram on the right); e.g. the right pointing arrow implies contraction of a tensor with its complex conjugate over the left, up, and down indices yields the identity matrix. The procedure is repeated until all columns have been contracted. value decompositions), which scales like ( ) for the IF time evolution procedure (sideways contraction by column). In this paper, we use such an approximate transverse contraction scheme (compressing to a fixed bond dimension of the contraction intermediates) to compute the IF (Fig. 2(b)) for systems with arbitrary baths. The algorithm is to iteratively contract and compress the columns from the edges of the bath inwards to the sites connected to the subsystem (Fig. 3). Assuming the subsystem is defined as the leftmost site, we start from the rightmost boundary column, and then the column is absorbed leftward to make a new boundary column, which is compressed using standard MPS compression to a pre-specified maximum bond dimension . This is the standard "boundary contraction" algorithm of 2D tensor networks 46 . Key to the success of the algorithm and the quality of the compression is the choice in gauge (a redundant degree of freedom in all tensor networks). We choose the gauge as shown in Fig. 3. In the case where the bath extends both to the left and right of the site of interest, we compute the IFs corresponding to the left and right bath sites separately (see Appendix for details). This contraction scheme was implemented within the Quimb tensor network library 47 .
For purposes of comparison, we will also present reference dynamics generated by standard MPS time evolution (ie. contracting the space-time tensor network in the usual time direction) 37,38,40,48 . Because the underlying full system dynamics is governed by Hamiltonian evolution in the problems that we study, we have the option to apply − as a commutator to the square root of the density matrix ([ − 1∕2 ][ 1∕2 † ]) or via − directly. We refer to the former as Hilbert time evolution (HTE) and the latter as Liouville time evolution (LTE). In HTE, the compressed tensor network dynamics is carried out for the pseudowavefunction = − 1∕2 49-51 . HTE has the advantage that the compressed density matrix is always positive definite, although correlations between the bra and ket sides of the density matrix are less compressible. In existing literature, this is sometimes referred to as purification-based time evolution. In the case that is a pure state, this method is equivalent to traditional MPS Hilbert space time evolution.

A. Spin Boson Model
First, we consider the well-studied spin-boson model, in which a single spin is linearly coupled to a bath of noninteracting harmonic oscillators, where Δ is the tunneling strength between the two subsystem states, and the system-bath coupling strength ( ) is determined from the bath spectral density function ( ) by In the case of an Ohmic bath with exponential cut-off, where is the Kondo parameter and is the cut-off frequency. Typically one computes the dynamics from a factorized initial state | ( 0 )⟩⟩ | ( )⟩⟩, where | ( )⟩⟩ is the Gibbs thermal state of the isolated bath at finite temperature . In this paper, we set Δ = 1.0, = 7.5, and = 5.0. The spin-boson model exhibits a dynamical phase transition from thermalizing to localizing behavior at = 1.0 + (Δ∕ ) 23 , and is often cited as an example of physically relevant non-Markovian dynamics 1,52 . Because of the linear coupling and harmonic bath, the IF may be computed via an analytical expression. There already exist several methods of obtaining accurate dynamics for various bath coupling strengths and spectral densities 8,23,[53][54][55] . We thus use this model as a benchmark to understand the properties of the influence functional, its compressibility, and the accuracy of the tensor network contraction approximation.

Compressibility of Analytical IF
We first investigate the compressibility of the analytical expression for the influence functional for the spinboson model. Denoting the system basis , |−1⟩} and |1⟩ , |−1⟩ are the eigenstates of the operator, then the influence functional can be written as This explicitly shows the form of the influence functional as the Boltzmann weight of a complex spin Hamiltonian with the spins interacting along the time axis via the long-range pairwise "interaction" ′ . We can further factorize the weights into contributions for times ( 1 ),  8 . Then, one can evaluate the influence functional (or its effect on the dynamics 54,56 ) with a computational cost exponential in max . Since QuAPI often converges rapidly with max , one also expects the matrix product state representation of the IF (Eq. (5)) to be compressible to small bond dimension. One way to verify this would be to construct the large influence functional object as an exact tensor, and then compress it into a matrix product state. Because of the exponential storage of the tensor with time, this is possible only for a small number of time points . Alternatively, one could build the influence functional iteratively (i.e. piece by piece in Eq. (14)) and compress at each step. This is the idea behind TEMPO and related methods 23,24 which exploit the compressibility of the augmented density matrix, the influence functional applied to the subsystem density matrix, i.e.
Note that because is composed of commuting pieces there are many possible decompositions and thus sequences of iterative constructions.
To verify the compressibility of the IF itself, we approximate with no memory cutoff ( max = ∞) as an MPS of bond dimension using an iterative scheme (see Appendix) and determine the error in the resulting on-site dynamics, using the = 256 result as reference. Here and throughout the paper, the error is computed as the r.m.s. deviation of the dynamics of an observable with respect to some reference over the time interval of the plot. The = 256 results are used as the reference because most of the = 128 dynamics calculations are already converged to within an r.m.s. error of 0.001 with respect to it. Fig. 4 compares convergence with respect to for infinite max , as well as convergence with respect to max for fixed = 256. As expected, the influence functional generally becomes less compressible with increased Kondo parameter . However, while convergence of the IF with respect to max and seem comparable for baths with larger cut-off frequencies ( = 7.5), the IF converges more quickly with respect to for smaller . The key advantage of using a compressed matrix representation, as opposed to truncating , + at some max as in QUAPI, is that this does not eliminate the effects of long-range memory 23 . Overall, this result confirms that for certain spectral densities, the influence functional can be efficiently written as a low-rank MPS, consistent with the findings from TEMPO 23,24 .

Finite Size Harmonic Bath
We now use the spin-boson model to examine if the IF can be constructed efficiently from our tensor network contraction scheme. To do so, we consider a finite size harmonic bath with sites, where the bosons are capped to some finite number of states. For small baths and small boson cap, the transverse contraction can be performed without compression, allowing a numerical test of the compression procedure. The analytical IF, assuming bosons with infinite boson cap, can also be computed. The bath is characterized by the discretized spectral density, We use a linear discretization of the bath sites, such that the bath density is ( ) = ∕ , where is the maximum boson frequency used. Here, we set = 10. First, we consider a system with only 2 bath modes (in the number basis) each with a maximum boson number of 2, for = 100 time steps of size = 0.05. For this small system, we use the exact time evolution operator of and compute the IF by exact transverse contraction, applying compression only to the final IF object. Fig. 5(a) shows the error of the exact IF compressed to bond dimension . As expected, the IF is much less compressible than with the continuous bath density in the last section, due to the small bath size. The error decreases only slightly until it drops suddenly once the bond dimension is large enough to capture the IF exactly, and further, the compression error increases with . We then perform the same analysis but for IFs computed by transverse contraction with compression, for both exact and RK4 time evolution. As seen in Fig. 5, the errors are comparable to those obtained when compressing the final exact IF. This indicates that for this small problem, there is little additional error added by the iterative contraction, and that the time-step error is negligible: the error is dominated by the compressibility of the final IF itself (which is low when the bath size is small).
Next we investigate systems with larger bath sizes in Fig. 6(a). We first examine the time-dynamics of the analytical IF (i.e. without any boson cap, and without transverse contraction) for a discretized spectral density with 11 bath sites, as well as the IF computed by transverse contraction, using a boson cap of 2. Because the system size is so small, the exact reference dynamics for a boson cap of 2 can be generated by direct MPS time evolution (here we use = 64 and Hamiltonian time evolution). From this comparison we observe two things. First, compared to using the continuous bath density, the error of the analytical IF dynamics is increased, although it is still somewhat compressible. For the same , the errors us- ing transverse contraction are larger, suggesting that at intermediate points in the transverse contraction, there is more time-like entanglement than in the final IF itself. In Fig. 6(b) we show the time-averaged error of the IF dynamics as a function of the number of bath sites. We see that this error decreases as the number of bath sites increases, both for the analytical IF and the transverse contraction. This is consistent with the idea that smoother bath densities are more "compressible".

B. 1D Hard-core Boson Model
We next study dynamics of a 1D hard-core boson (HCB) lattice model, in which each lattice site is either unoccupied (|0⟩) or occupied (|1⟩) by a single bosonic particle. We consider the Hamiltonian = |0⟩ ⟨1| are hard-core boson creation and annihilation operators at the th lattice site, and = † is the number operator. This Hamiltonian is intended to mimic the Bose-Hubbard model dynamics often simulated by cold-atom experiments, which was shown to be difficult to compute using direct time evolution methods 57 , where the on-site interaction term is replaced by a nearest-neighbor interaction term. We also include a harmonic potential term to emulate a cold atom trap. We assume a pure initial state |0, 1, 0, ...⟩ such that there is one particle at every other lattice site, and set the parameters = 1 and = 10 −2 while varying . For non-zero interaction term , there is no analytical form for the IF.
We compute the dynamics of ⟨ ( )⟩ at lattice sites = {0, ∕4, ∕2}, where is the length of the 1-D chain. For longer chains, the rapid growth of entanglement means that direct MPS time evolution (either using Hamiltonian evolution, denoted HTE or Liouvillian evolution, denoted LTE) with a finite can only obtain converged dynamics up to a finite time. We consider two chain lengths: = 9 where converged HTE MPS dynamics can be used as a reference, and = 43, where the HTE MPS dynamics appear to be not fully converged (no longer within 0.03 of = 512 results) for the full simulated time. To obtain the dynamics using the influence functional method, we partition the lattice such that The = 0 dynamics for = 9 and = 43 is shown in Fig. 7. Direct MPS LTE shows unphysical behaviour for large system sizes, presumably because of the loss of positivity at some point in the dynamics. In contrast, the dynamics obtained using the iteratively contracted IF are more stable, highlighting the innate compressibilty of time evolution tensor network along the time axis as opposed to the spatial axis. For = 9, the IF dynamics only appears to begin to converge by = 128 with respect to the (exact) HTE dynamics, having less than 0.03 r.m.s error in ⟨ ( )⟩ over the simulated time interval and deviations within 0.08. Thus, there appears to be no significant advantage to using the the IF method over direct HTE for small system sizes. Conversely, for the = 43 system, the IF dynamics are converged with respect to = 256 results by = 64, with less than 0.02 r.m.s error and a maximum deviation of 0.04 for = 10 and less than 0.001 r.m.s error and a maximum deviation of 0.004 for = 21. Note that the = 10 dynamics converge more slowly because effectively the site is coupled to two separate baths, one of which is small. However, the IF method still outperforms direct HTE.
For > 0, as shown in Fig. 8, the = 64 results appear less converged than the non-interacting case, but the = 128 results are converged (r.m.s. errors of the = 128 observable dynamics with respect to the = 256 results are less than 0.03), for times longer than that accessible by direct time evolution. This shows that the IF-based dynamics can produce the correct oscillatory behaviour of the density as a function of time, which is not captured by the direct MPS time evolution despite using a larger bond dimension (this difficulty with the long time oscillatory tail has previously been noted in other cold atom simulations 57 ). However, while the IF method notably outperforms direct HTE at small , the two methods become comparable at larger ≈ 1.5 once the oscillatory tail is sufficiently dampened.

C. Entanglement Spectrum
The accuracy of the transverse contraction scheme depends on the entanglement in the time-like direction. Recall that our contraction algorithm starts with the farthest column (an MPS), and at each iteration another column is contracted into this boundary. Thus, as the iteration number increases, the boundary column represents more of the bath. For both the spin-boson model and HCB  Fig. 7. Compared to the = 0 case, a larger bond dimension is needed, particularly at around = 3.0 where the r.m.s. error peaks. In contrast, the HTE dynamics converge more quickly with increasing . model, we measure the singular values at the middle of the boundary "bath" MPS during each step of the iterative contraction scheme. The entanglement entropy (EE) and spectrum of the singular values (normalized so that ∑ | | 2 = 1) are plotted in Fig. 9. For the spin-boson model, only results for = 1.0 are shown; the only notable difference for other is that the EE increases with . Consistent with the observations in our simulations above, the EE decreases as one increases bath size. For the SB model, the EE decreases with increased number of bath sites in the discretization until convergence. For the HCB model, only if sufficiently large enough bond dimension is used does the EE decrease with increasing iteration number. Otherwise, the EE stays at a large value throughout the contraction scheme and the gap between the dominant and non-dominant singular values decreases; this makes the EE of the smaller approximation larger than that of the larger approximation. Overall, this suggests that the final compressibility of the IF emerges from the cancellation of many different correlations as one iteratively contracts out the bath.

IV. CONCLUSIONS
In this work, similar to some other recent contributions 22,26 , we have used the representation of the influence functional within the tensor network language, motivated by the limitations of modeling spatial entanglement growth in quantum dynamics. We have discussed a transverse tensor network contraction algorithm that allows us to compute the influence functional in cases where the analytical form is not known. We have applied this algorithm to study both the canonical spin-boson model as well as an interacting hard-core boson chain where the bath is not quadratic (i.e. interacting). We find that the compressibility of the influence functional is controlled by several factors, principally the size of the bath, as well as the nature of the interactions. In addition, although the time-like correlations may ultimately be short-ranged in the final influence functional, during the transverse tensor network contraction to construct it, it is possible to proceed through intermediate quantities with larger timelike entanglement. This suggests a complicated picture where time-like correlations first accumulate as the bath is integrated out before finally cancelling in the influence functional itself. In the regimes where the influence functional and all intermediate quantities are compressible, as in some interaction regimes in the interacting hardcore boson model we have studied, it is possible to outperform conventional tensor network time evolution methods at longer times.
There are many possible directions for further investigation. For example, there are natural extensions to higher-dimensional interacting problems and fermionic systems, as well as more complicated correlation functions. Also, a better theoretical understanding of how correlations grow and cancel out in the transverse contraction scheme may lead to a deeper understanding of the generation of memory in quantum dynamics, the master equation formalism 17 , improved contraction schemes, and ultimately new algorithms to carry out longer time dynamical simulations.

VI. DATA AND CODE AVAILABILITY
Data and code is available upon request.

Appendix A: Transverse Contraction Procedure
We describe a single step of the transverse contraction procedure, which in essence is performing an MPO-MPS contraction and then compressing the resulting MPS to the desired bond dimension , described in detail in Ref. 40. We assume that all tensors have been prepared in the appropriate gauge (discussed later) as shown in Fig. 3. This algorithm is also the standard boundary contraction algorithm for 2D tensor networks, as described in Ref. 46. For simplicity, consider the influence functional tensor network for bath sites, as discussed in the main text. We start with the farthest boundary column MPS (the column for the th bath site) which is given by where ( ) and ( ) are the right-most tensors of the density matrix MPS (Eq. (7)) and Liouville evolution operator (Eq. (8)), respectively, and̃ ( ) denotes the Liouville evolution operator with traced out bath degrees of freedom. The next column can be interpreted as an MPO, given by where the sum over indicates sums over all ( , ), and  We are able to define canonical forms because tensor networks have a gauge degree of freedom. This means that the choice of tensors in the network is not unique.
One can see this by introducing a set of matrices , −1 , which clearly satisfy −1 = , along any line connecting two tensors.
Canonicalization can be performed using singular value decomposition (SVD). Suppose that we are interested in written the MPS in left canonical form, and that all tensors left of the th tensor are already in left canonical form. We then take the SVD of the th tensor, where Σ are the singular values from the decomposition of the tensor. Note that by definition, is left canonical, as desired, and thus will be used as the new th tensor. The remaining matrices are then pushed into the ( + 1) th tensor, By iteratively performing this operation starting from the left-most tensor all the way to the right end of the MPS, the MPS is put into left canonical form. The procedure for expressing the MPS in right canonical form is analogous. MPS compression is performed in the same way, except only the largest singular values are retained, generating some error. For minimal compression errors, the MPS must be in left (right) canonical form prior to performing the iterative compression procedure starting from the right (left) end.
In Fig. 3 we depict the contraction of the columns of the (1+1)D influence functional tensor network. The rows are initialized in left canonical form. After the two rightmost columns are contracted, the product is canonicalized and the compressed using the procedure discussed above. Because of the vertical orientation, the left and right canonical forms are depicted by triangles pointing upwards and downwards along the column.

Appendix B: Matrix Product Form of Analytical IF
As discussed in the text, the expression for the analytical IF in discretized time steps is given by Eq. (13) and we wish to write it in the MPS form with physical bonds that index the states of the density matrix at each timestep, ( 1 ) 1 ( 2 ) 1 , 2 … ( ) −1 (B1) One possible way to construct the IF is to take the product of factors in Eq. (13) in the order = 0 1 … −1 . We start by using the 0 terms which are in the form of a product state (MPS with bond dimension 1). We then multiply by each of the subsequent and compress into an MPS after each is applied. Multiplying by can be viewed as multiplication by an MPO where the tensors are very sparse. The two-body terms for > 1 are long-range operators, and must be padded with identities to skip over the times in the middle. More explicitly, the MPOs are where we first decompose into two tensors (eg. via SVD or QR decomposition) and then define the tensors in the MPO as In our contraction scheme, we start from the bottom row and contract upwards. However, because each factor commutes with the rest, other choices of ordering are possible and are the basis of algorithms such as TEMPO 23,24 .
Appendix C: Influence functional transverse contraction around an arbitrary site Sometimes the site whose dynamics we are interested in may be at the middle of the MPS representation of the system (e.g. in the hard core boson model). Thus, we need to generalize the tensor network diagrams presented in the main text to consider IFs for subsystems at arbitrary lattice site .  Fig. 1). It is cheapest to contract this network vertically from the row at one end and continuing to the other end. Note that in using this method one does not explicitly compute the IF itself.
If the Hamiltonian only consists of nearest-neighbor interactions, the IFs from the two sides of site are separable and can be computed independently. Otherwise, the tensor network can be initialized as shown in Fig. 11, and one contracts inwards from the outer columns separately. Once only the column corresponding to the isolated site is left, one can now include on-site terms (initial state, onsite time evolution operators, observable) such that the network now corresponds to the expectation value of the observable at the desired time step (a scalar). The cost of contracting this network scales like ( 3 ), which is much cheaper than explicitly computing the full IF first and then computing the observable expectation values.