A Pairwise Connected Tensor Network Representation of Path Integrals

It has been recently shown how the tensorial nature of real-time path integrals involving the Feynman-Vernon influence functional can be utilized using matrix product states, taking advantage of the finite length of the non-Markovian memory. Tensor networks promise to provide a new, unified language to express the structure of path integral. Here, a generalized tensor network is derived and implemented specifically incorporating the pairwise interaction structure of the influence functional, allowing for a compact representation and efficient evaluation. This pairwise connected tensor network path integral (PCTNPI) is illustrated through applications to typical spin-boson problems and explorations of the differences caused by the exact form of the spectral density. The storage requirements and performance are compared with iterative quasi-adiabatic propagator path integral and iterative blip-summed path integral. Finally, the viability of using PCTNPI for simulating multistate problems is demonstrated taking advantage of the compressed representation.


I. INTRODUCTION
Tensor networks (TN) are designed to be compact "factorized" representations of high-ranked tensors. Probably the most common use of TN in physics is related to representations of the quantum many-body wavefunction which, in general, is also a high-ranked tensor. This use has been widely demonstrated in a multitude of methods such as the density matrix renormalization group (DMRG) [1,2] which uses a Matrix Product State (MPS) [3,4] representation, and multiconfiguration time-dependent Hartree (MCTDH) [5] and its multi-layer version (ML-MCTDH) [6][7][8] which use tree tensor networks. For multidimensional systems, an "extension" of MPS to multiple dimensions called projected entanglement pair states (PEPS) [9] is used. For systems at critical points, an MPS representation does not work because of long-range correlations necessitating the use of the so-called multi-scale entanglement renormalization ansatz (MERA) [10,11]. Tensor networks, since its introduction, have proliferated in various diverse fields requiring the use of compact representations of multidimensional data like machine learning and deep neural networks.
While quantum dynamics at zero temperature can often be simulated using wave-function based methods like time-dependent DMRG [2,12,13] or MCTDH, at finite temperatures, owing to the involvement of a manifold of vibrational and low frequency ro-translational states in the dynamics, they suffer from an exponentially growing computational requirements. Feynman's path integral provides a very convenient alternative for simulating the time-dependent reduced density matrix (RDM) for the system. The vibrational states of the "solvent" introduced as harmonic phonon modes under linear response [14] are integrated out leading to the Feynman-Vernon influence functional [15]. Identical influence functional also arises in dealing with light-matter interaction through the integration of the photonic field.
The primary challenge in using influence functionals and path integrals is the presence of the non-local history-dependent memory that leads to an exponential growth of system paths. While many recent developments have helped improve the efficiency of simulations [16][17][18][19][20], each of them utilize very different and deep insights into the structure of path integrals. It has recently been shown that the MPS representation can be very effectively utilized to reformulate real-time path integrals involving the influence functional leveraging the finite nature of the non-local memory [21][22][23][24]. While the MPS structure is the simplest tensor network that can be used, the 1D topology is probably not optimal when the non-Markovian memory spans a large number of time-steps and suffers from growing bond dimensions. In this paper, an alternate generalized tensor network that directly captures the pairwise interaction structure of the Feynman-Vernon influence functional, is introduced. This pairwise connected tensor network path integral (PCTNPI) has an extremely compact representation, that can be efficiently evaluated, allowing us to go to much longer non-Markovian memories without resorting to various techniques of path filtration. Tensor networks show great promise in being a unifying language for formulating and thinking about path integral methods.
The construction and evaluation of the tensor network is discussed in Sec. II. In Sec. III, we illustrate some typical applications of the algorithm. The memory usage is also reported for various parameters. The implementation of this method utilized the open-source ITensor [25] library for tensor contractions allowing for extremely efficient tensor contractions using highly efficient BLAS and LAPACK libraries. We end the paper in Sec. IV with some concluding remarks and outlook on future explorations.

II. METHODOLOGY
Consider a quantum system coupled to a dissipative environment described by a Caldeira-Leggett model [26][27][28] whereĤ 0 is the Hamiltonian of the D-dimensional system of interest shifted along the adiabatic path [29]. If the quantum system can be described by a two-level Hamiltonian, thenĤ 0 = σ z − Ωσ x , whereσ z andσ x are the Pauli matrices.Ĥ env represents the Hamiltonian of the reservoir or environment modes which are coupled to some system operatorŝ. The strength of the j th oscillator is c j . While we are using a time independent Hamiltonian for simplicity, time-dependence from an external field in the system Hamiltonian can be captured through the corresponding system propagator in a straightforward manner. For a problem where the environment is in thermal equilibrium at an inverse temperature β = 1 k B T , and its final states are traced out, the interactions between the system and the environment is characterized by the spectral density [14,26] In fact, the spectral function, S(ω) corresponding to the collective bath operator X = − j c j x j is related to the spectral density as follows [30]: For environments defined by atomic force fields or ab initio calculations, it is often possible to evaluate the spectral density from classical trajectory simulations. The dynamics of the RDM of the system after N time steps, if the initial state is a direct product of the system RDM and the bath thermal density is given as: where Here,Û is the short-time system propagator for ∆t and {s + j } and {s − j } are the forward-backward system paths. The Feynman-Vernon influence functional [15], F [s ± j ], is dependent on the system path s ± j and the bath response function that is discretized as the η kkcoefficients [31,32]. The influence functional depends upon the history of the system path, leading to the wellknown non-Markovian nature of system-environment decomposed quantum dynamics. Notice that it can be factorized based on the "range" of interaction in the following manner: The influence functional creates pairwise interactions between points that are temporally separated. As it has been shown, if MPS and MPO are used to model the influence functional, the fact that these interactions can spread across long temporal spans leads to an increase in the effective bond dimension. Here, the goal is to create a structure that naturally and efficiently accounts for the pairwise interactions that span long temporal separations while not being associated with any one particular representation.
To motivate the tensor network representation, first consider the Markovian part of Eq. (5), involving just the propagators and the terms of the influence functional coupling consecutive time points. These terms can be simply rearranged as: Here, we are implicitly summing over repeated indices that do not appear on both sides of the equation. The ± labels on the site indices of the tensors are omitted for convenience of notation. The superscript, 1, on P is there to denote the maximum distance of interaction that we have incorporated. Equation (9) is already a tensor network; more specifically it is series of matrix multiplication as shown in Fig. 1. Let us now bring the "next-nearest neighbor" interactions I (2) sj−2,sj . Clearly, it is not possible to directly contract the I (2) sj−2,sj tensor to the P (1) tensor because the internal s j 's have already been traced over. To make it possible to incorporate the I (2) tensors, we augment the K tensors as follows: It is convenient to think of lower indices as the "input" indices and the upper indices as the "output" indices, though there is no other mathematical significance to the positioning of the indices. With this input-output convention in mind, it is easy to see that the internal augmented K tensors duplicate and flip the order of the input indices, s j . This ensures that indices that differ by two time steps are now placed adjacent in the output layer. Now, the Markovian terms and the I (2) terms can be combined and we get: which is depicted in Fig. 2. Notice that the index s j still connectsK sj−1,sj andK sj ,sj+1 , as in P (1) , but now there is another connection that goes through the I (2) tensor in a "triangular" form. This feature of an I (α) with a higher α acting as a bridge betweenK or I tensors with smaller values of α would become a recurring motif in this tensor network. The pattern for inclusion of the rest of the non-local interactions is quite similar. Note that in Fig. 2, if we did the same "trick" of duplicating and flipping the order of the inputs, in the next layer indices that differ by three time points, like s 0 and s 3 , s 1 and s 4 , are going to be adjacent. Hence, this can now be multiplied sj−3,sj . Continuing like this, we can complete the network. The diagram is shown in Fig. 3 [33], these augmented tensors are going to be written as K, I (2) , . . . . The tensor network shown in Fig. 3, which we will conventionally denote by P (∞) , represents the final Green's function for the propagation of the system RDM having incorporated the non-local influence from the environment. So, ρ(t) = P (∞) ρ(0).
If the system is defined to have D states, then in Fig. 3, all the indices have D 2 dimensionality corresponding to each of the possible combination of forward-backward states. However, this is not optimal. Notice that the influence functional tensors, I (α) s k−α ,s k for a time difference of α, depends only on the "difference" coordinate, ∆s k = s + k − s − k of the latter time point. So, currently, we are carrying over more information than we need to.
To take care of this redundancy, we need to redefine the K tensors to not just duplicate the input indices, but to project the "latter" index onto its difference coordinates as follows: Notice that the upper-right output indices in the diagrams remain exactly the same. Only the upper-left output index of K changes. Therefore, the tensor K r1 s0,s1 remains unchanged. The dimensionality of the "l" indices is the number of unique values of ∆s = s + − s − that the system can have. For a general D-level system, this value is B = D 2 − D + 1 instead of D 2 , however the actual symmetries present in the system might reduce this even further. Finally, the influence functional tensors have to be changed to be consistent, viz.Ĩ (α) these changes, the basic topology of the network remains the same. The new network with the different dimensions is shown in Fig. 4.
Having discussed the tensor network, now let us turn to the job of contracting it. Typically, many tensor networks are constructed using singular value decomposition (SVD) and evaluated via the truncation of the singular values [21,24]. The PCTNPI network is constructed without resorting to any SVD calculations and consequently "exact." The goal now is to find an optimal contraction scheme that preserves this "exactness." The storage cost, S, is also evaluated at the end of every step. The canonical contraction order that we discuss below has been marked out in cyan arrows in Fig. 4. For a simulation with N time steps: s0,s N and contract it withĨ  9. Finally, the last tensor,K s N −1 ,s N is contracted.

Multiply allĨ
In the above contraction scheme, we multiply the initial condition, ρ s0 , and get the final RDM. While this leads to a more efficient algorithm in terms of the storage and computational cost, it is possible to reformulate the scheme in terms of the Green's function by not involving the initial condition in the contractions and evaluating P (∞) . An in-depth analysis of the memory and computational cost is given in Appendix A. Of course, the storage requirement grows to a maximum of (D 2 ) 2 B N −2 before decreasing continuously. This naïve contraction scheme does not solve the problem of storage. Still, as would be illustrated in Sec. III, PCTNPI outperforms both traditional iterative quasi-adiabatic propagator path integrals (QuAPI) [31,32], and iterative blip summed path integral (BSPI) [17,18], when used without path filtration, in the memory lengths that can be accessed without any sort of filtration. In a future work, filtration schemes on top of PCTNPI would be introduced that can not only deal with this problem, but would also avoid the construction and storage of the full tensor. The focus of this paper is however on the tensor network and its performance in the most naïve implementation.
It is well-known that the non-local memory of the influence functional dies away with the distance between the points, allowing for a truncation of memory. This idea is commonly used both in Nakajima-Zwanzig generalized quantum master equations [34][35][36] and iterative QuAPI [31,32]. In the framework of PCTNPI, the length of the non-Markovian memory is equal to the depth of the resultant network. The topmost tensor encodes the interaction between the most distant points, while the bottom most tensor captures the Markovian interactions coming through the propagator and the I (1) terms.
At two time-steps of memory, that is L = 2, we basically get Fig. 2. In Fig. 5, we show the structure of the network for a 5-step propagation with L = 3. Because s 0 does not interact with s 4 or s 5 , it is not necessary to store and evaluate the full diagram at once, but it can be built iteratively. The first edge, corresponding to interactions with s 0 is contracted, and multiplied by the second edge, using the canonical contraction scheme discussed previously. As soon as this is done, the storage of the first edge can be freed, and the third edge can be contracted. This iteration scheme turns out to be identical to the iteration scheme in iterative QuAPI. The first steps of the iteration algorithm is pictorally outlined in Fig. 6. Makri [17] has shown that it is possible to think of the memory as arising from two different causes. The influence functional F can be rewritten in terms of the real and imaginary parts of the η-coefficients as: where ∆s k = s + k − s − k ands k = 1 2 (s + k + s − k ). The part of the influence functional that arises from Re η kk is called the classical decoherence factor. It corresponds to stimulated phonon absorption and emission [37]. This can also be obtained through classical trajectory-simulations and reference propagators [38] in a Markovian manner. All effects of temperature is captured in the classical de-coherence term. The term with the Im η kk is the backreaction that leads to quantum decoherence. This part of the memory is truly non-local and temperature independent.
As a cheap approximation to the dynamics, it is possible to do a simulation with classical decoherence, that would become increasing accurate as the temperature of the simulation rises. In this, the full η kk coefficients are used only when k = k or k = k + 1, and otherwise the imaginary part of η kk is ignored. (Actually, the true expressions for classical decoherence would include the full η kk coefficients only when k = k and the real part otherwise. In PCTNPI, we can include the case of k = k + 1 as well at the same storage and computational cost.) Effectively, we are modifying the I (α) s k ,s k operators to be exp (−1/ Re(η kk )∆s k ∆s k ) when α = k − k ≥ 2. Just like before when the s k lines carried unnecessary information, now the s k lines carry more information than they need to. We only need to know about ∆s k . Thus we can make the required changes to the dimensionality of the indices by putting in the corresponding projector operators in the K tensors, thereby reducing the cost of computation even further. The network for the classical decoherence simulations would have exactly the same structure as Fig. 3 with all edges except the base ones being B dimensional. This approximation is especially accurate at short times.

III. RESULTS
As illustrative examples, we apply PCTNPI to a twolevel systems (TLS) coupled bilinearly to a dissipative environment:Ĥ The dissipative environment is chosen to be defined by Ohmic model spectral densities, which are especially useful in modeling the low frequency ro-translational modes. We use the very common Ohmic form with an exponential cutoff, where ξ is the dimensionless Kondo parameter and ω c is the characteristic cutoff frequency and the Ohmic form with a Drude cutoff, where κ is a measure of the coupling strength. Generally these model spectral densities are often thought to be fully characterized by a reorganization energy (22) and the cutoff frequency, ω c . The reorganization energies for the exponential and the drude cutoff spectral densities are as listed below: As we demonstrate through the examples, though the reorganization energy and the cutoff frequency are same, the exact dynamics of the reduced density matrix is highly dependent on the form of the "decay function." Consider a symmetric TLS ( = 0) and Ω = 1 interacting strongly (ξ = 2) with a sluggish bath (ω c = Ω) initially localized on the populated system state 1. The bath has a reorganization energy of λ = 4 and is held at an inverse temperature of Ωβ = 1. The dynamics was converged at ∆t = 0.125, and a memory length L = 16. The convergence is shown in Fig. 7 (a) for an Ohmic bath with an exponential decay. Full quantum-classical simulations for this parameter is available [39]. If the Drude form of decay is used, the dynamics changes quite signif-icantly. The comparison between the dynamics arising from the two spectral densities is shown in Fig. 7 (b). Next, consider a case where not only is the dynamics different between the two different decay functions, but the converged non-Markovian memory length is different as well. The dynamics of the same TLS as above ( = 0, Ω = 1) is now simulated in a bath with the reorganization energy λ = 8 and a characteristics cutoff frequency ω c = 5. The bath is equilibrated at an inverse temperature of Ωβ = 5. The time-step is converged at Ω∆t = 0.125. The convergence of the dynamics of the reduced density matrix on changing the memory length, L, is shown in Fig. 8. While the memory length for the exponential decay function spectral density is quite close to convergence at L = 14, for the Drude spectral function, it converges at L = 10.
In Fig. 9, we consider a TLS coupled to a strongly coupled Ohmic bath with an exponential cutoff (ξ = 1.2, ω c = 2.5Ω) equilibrated at a high temperature Ωβ = 0.2. The converged time step is ∆t = 0.125. The classical memory calculations converge at a comparatively lower memory length, L and agree quite well with the full simulations at short times. Though at intermediate and long times, the classical decoherence dynamics differs from the true dynamics, this can often be enough for estimating timescales of processes, especially using rate theory [40][41][42][43].
Next, the storage requirements of PCTNPI is compared with that of iterative QuAPI and iterative BSPI in Fig. 10 (a). To keep the comparisons fair, the iQuAPI and iBSPI methods were run without any path filtering. It is quite clear from the plot that the scaling of PCT-NPI is essentially "like" that of iBSPI, i.e. for a TLS, 3 L scaling for iBSPI and PCTNPI vs 4 L scaling of iQuAPI. However, the prefactor is much smaller, allowing us to access much longer memories with limited resources. In fact, this difference in the prefactor would grow with the dimensionality of the quantum system. In iBSPI, there would be 3 L paths for a memory length of L, but the storage is more than just a number corresponding to each path. It stores a small dimensional matrix for each path. This is the cause of the larger prefactor.
A comparison of the run times of PCTNPI with respect to iBSPI without any filtration for a simulation of 100 time steps is presented in Fig. 10 (b). A laptop with Intel® Core TM i5-4200U CPU with a clock speed of 1.60GHz was used for these benchmark calculations. These measurements are not going to be consistent with similar benchmarks run on other machines, but the basic trends would continue to hold. The PCT-NPI algorithm is built on top of ITensor and automatically uses parallel BLAS and LAPACK wherever possible. There is no standard iBSPI code. The iBSPI program used for these benchmarks was manually parallelized with OpenMP loop parallelization.
As a final example, consider a molecular wire described by the tight-binding Hamiltonian involving N sites:   The site energy of the j th site is j and the nearest neighbor couplings are V . The sites are separated by unit distance such that |σ j are eigenstates of the position operator,ŝ |σ j = (j − 1) |σ j . The site energy of all but the first site is chosen to be zero j = 0 for j = 1 and 1 = 1. The intersite coupling is chosen to be V = 0.025 [44]. The computational cost grows exponentially with the number of sites. To test the efficiency of the basic contraction scheme outlined here, we use a system with N = 4 sites. The bath is characterized by an Ohmic spectral density with an exponential cutoff, Eq. (20) with ω c = 4 and ξ = 0.12 [24] equilibrated at an inverse temperature of β = 0.1. As discussed in Sec. II, the scaling of the algorithm would go as B < D 2 . The symmetry of the Hamiltonian in this case ensures that the number of unique values of ∆s, B = 7 for this 4 state system, which is even less than the D 2 − D + 1 for a completely general Hamiltonian. The population dynamics of all the states is shown in Fig. 11. An initial state with only the first site populated was used. Because of the high temper- ature of the bath, the classical decoherence simulation produces practically identical dynamics but converges at a smaller memory length L.

IV. CONCLUSION
A novel tensor network is introduced to perform path integral calculations involving the Feynman-Vernon influence functional. This pairwise connected tensor network path integral (PCTNPI) captures the pairwise interaction structure of influence functional. PCTNPI can be contracted efficiently, and minimizes the storage requirements as far as possible without resorting to various path filtration algorithms. Iterative decomposition of the memory is also possible in an elegant manner. Comparisons between PCTNPI and iQuAPI and iBSPI show the scaling of memory requirements of PCTNPI to be similar to iBSPI, but much smaller.
PCTNPI provides an alternative to the MPS representation [21,24], serving as a small step in further elucidating the deep relation between tensor networks and path integrals. While no path filtration scheme has been developed, PCTNPI is already quite useable. It can easily incorporate classical trajectories through harmonic backreaction quantum-classical path integrals [37] thereby making it possible to include anharmonic effects of the environment in an approximate manner without any additional cost. Additionally, harmonic backreaction also leads to an increase in the converged time-step and a decrease in the effective memory length such that some ultrafast reactions can be simulated directly. Taking advantage of the extended memories that are accessible with PCTNPI, the combined method would be able to simulate systems with strongly coupled sluggish realistic solvents with high reorganization energy. This promises to be a fruitful avenue of research in terms of applications (a) Contraction of internal tensor along the left edge (j ≥ 2). to electron and proton transfer reactions.
Algorithms based on MPS representations of the augmented reduced density tensor [21] or of the pathdependent Green's function [24] can be thought of as particular optimized re-factorizations of the PCTNPI network. We have demonstrated the viability of evaluating PCTNPI in a brute force manner compared to other methods. This suggests that using the PCTNPI network directly to generate other optimized representations might also lead to novel methods.
While ideas of path filtration were not a consideration of the present paper, schemes based on the singular value decomposition (SVD) can be incorporated with PCTNPI, leading to a method that significantly reduces the storage, since the full tensor would not need to be computed and stored. This development would be discussed in a future publication. cost of contraction, the left "edge" of the triangular network is first considered. Consider contractingĨ (j) s0,sj , for j ≥ 2, with two D 2 indices and one B index, as schematically indicated in Fig. 12 (a). The part that has already been contracted has one D 2 index and (N − j) B indices. Therefore, the cost of contraction is (D 2 ) 2 B N −j+1 . The space requirement at this stage is D 2 B N −j+1 . To finish the contraction of the left-most edge of the triangle, we need to multiply byK s0,s1 ρ s0 leading to the tensor network shown in Fig. 12 (b). The resultant tensor does not have a index corresponding to s 0 because that has been traced over. The computational cost of this step is C = (D 2 ) 2 B N −1 and the storage becomes S = D 2 B N −1 . Now, the second parallel edge is to be contracted. This step however is started from the bottom, i.e. fromK s1,s2 .
The first contraction, shown in Fig. 13 (a), is the most costly step in the entire algorithm. The computational cost of this step is C = (D 2 ) 3 B N −1 and the storage requirement increases to S = (D 2 ) 2 B N −2 . Continuing with the other intermediate tensors of the first parallel edge, notice that the cost of contraction remains constant at C = (D 2 ) 3 B N −1 and the space required remains constant at S = (D 2 ) 2 B N −2 . Finally, the last, top-most tensor of this edge is to be contracted. This is illustrated in Fig. 13 (b). The computational cost is C = (D 2 ) 2 B N −1 . The storage cost now drops to S = D 2 B N −2 . Now, consider contracting a general diagonal edge, say the j th one. The resultant tensor from the previous contraction has one D 2 index and N − j + 1 B indices. Contracting theK tensor leads to a tensor with two D 2 indices and (N − j) B indices. The cost of this contraction is C = (D 2 ) 3 B N −j+1 and the storage is S = (D 2 ) 2 B N −j . For all the intermediate tensors at this stage, once again both the computational costs and the storage costs remain the same. On contracting the last tensor of this diagonal, the storage drops to S = D 2 B N −j .
Below we list the total computational cost for contracting each of the "parallel" edges. The edge number is given as the subscript.
The prefactor of the computational and storage costs is lower for classical decoherence simulations: It goes from a power of D 2 to the corresponding power of B. It is clear that the complexity of the entire contraction goes as O B N −1 and the peak storage requirement is O B N −2 .