The TEMPO network can be made independent of the system size

Abstract

The time-evolving matrix product operator (TEMPO) method efficiently simulates non-Markovian open quantum system dynamics by contracting a two-dimensional tensor network representing the time-discretized Feynman-Vernon influence functional. The standard form of this network has an unfavorable scaling with the size of the system. We show that, within the usual second-order time-step error, the dimensions of all internal network bonds can be reduced to $(2m+1)(m+1)$, where $m$ is the number of independent coupling channels. This dimension does not depend on the number of distinct eigenvalues of the coupling operators. The system size only enters in the final boundary layer of the network, making large system sizes accessible to TEMPO methods. This structure was discovered by Claude Opus 5.5 during guided algorithm development.

TEMPO represents Gaussian-bath influence functionals as two-dimensional tensor networks, but conventional internal bond dimensions can scale as $d^{2m}$ for system dimension $d$ and $m$ coupling channels. This paper asks whether that dependence can be removed without changing the usual second-order time-step accuracy. The construction expands each influence-functional gate through second order in the discretized bath correlations, then projects network bonds onto the span of coupling-eigenvalue monomials of degree at most two. This gives internal bond dimensions bounded by $(2m+1)(m+1)$, including dimension $6$ for one channel; linear dependencies or bath structure can reduce the dimension further. System-size dependence remains in the final boundary layer. Although the individual exponential expansion has error $\mathcal O(\delta t^6)$, projecting the full network introduces an estimated total error $\mathcal O(\delta t^2)$ at fixed physical time, matching the order of symmetric Trotter splitting. A damped-oscillator benchmark compares the reduced network with SVD filtering of the full network against an exact Gaussian reference. Influence-functional contraction runtime becomes independent of the Fock cutoff, with comparable accuracy at fixed MPO bond dimension. A second example simulates three and four emitters coupled to a shared two-dimensional bosonic reservoir, using bond dimensions around $200$ and roughly two days per simulation on one CPU. The reduction is approximate, not an exact finite-step identity, and does not remove system-size dependence from the entire simulation. Time-dependent coupling operators and continuous coupling spectra are identified as directions for further generalization.