The discretized path integral expression for the reduced density matrix (RDM) of a system interacting with a dissipative harmonic bath is fully entangled because of influence functional terms that couple the variables at different time points. The iterative decomposition of the path integral, which exploits the finite length of influence functional memory, involves a tensor propagator whose size grows exponentially with the memory length. The present Communication disentangles the path integral by recursively spreading the temporal entanglement over longer path segments, while decreasing its contribution. Eventually, the entangled term becomes sufficiently small and may be neglected, leading to iterative propagation of the RDM through simple multiplication of matrices whose size is equal to that of the bare system. It is found that the temporal entanglement length is practically equal to the bath-induced memory length. The small matrix decomposition of the path integral (SMatPI) is stable and very efficient, extending the applicability of numerically exact real-time path integral methods to multi-state systems.

Besides serving as the paradigm for quantum dissipation, the dynamics of a quantum system interacting with a harmonic bath offers a realistic and often very accurate simplification of processes in polyatomic molecular or condensed phase environments. Many analytical treatments have provided qualitative pictures and much insight into system-bath dynamics.^{1} While some wavefunction-based methods^{2,3} are able to treat system-bath dynamics in specific parameter regimes, these methods are limited by the growth of the full system-bath Hilbert space and, thus, are not practical for situations involving strongly dissipative environments with many low-frequency bath modes at moderate or high temperatures. Feynman’s path integral formulation^{4,5} is particularly attractive for treating system-bath dynamics,^{6} as it directly targets the system’s reduced density matrix (RDM) without having to deal explicitly with the huge Hilbert space of the bath. This is achieved by analytically capturing the effects of the harmonic bath on the dynamics of the system through the Feynman–Vernon influence functional.^{7}

The main practical issue with the path integral is that the influence functional renders the RDM dynamics nonlocal in time. In particular, the path coordinates at each time point are coupled to the coordinates at all the other time points, leading to a fully entangled structure [see Fig. 1(a)]. As a result, it appears that evaluation of the RDM which develops in an *n*-state system from a given initial condition after *N* time steps requires the full evaluation of a sum of (*n*^{2})^{N} terms. Unfortunately, Monte Carlo methods^{8} converge extremely slowly in the case of real-time dynamics because of the oscillatory phase factors arising from quantum time evolution. In practice, the exponential scaling with propagation time restricts the full path sum to short times, although the blip decomposition^{9} can sometimes extend calculations to ∼10^{2} time steps in incoherent regimes. The iterative quasi-adiabatic propagator path integral (i-QuAPI) methodology^{10–15} takes advantage of the finite length of influence functional interactions (a consequence of decoherence) to overcome the exponential scaling with the propagation length, leading to an algorithm that scales as $N(n2)\Delta kmax+1$, where Δ*k*_{max} is the length of influence functional couplings (in units of the path integral time step Δ*t*). The iterative decomposition employs a propagator tensor of rank Δ*k*_{max}, i.e., a matrix of path segments spanning Δ*k*_{max} time points; thus, the basic i-QuAPI algorithm requires the storage of an array of $(n2)\Delta kmax$ elements. A number of techniques have been developed^{16–21} for reducing the size of this array, and the resulting size reduction can be dramatic in some regimes. Still, iterative path integral calculations often become prohibitively costly as the number of system states increases. An efficient alternative is offered by the quantum-classical path integral (QCPI) methodology,^{22–24} which treats the dynamics of general polyatomic environments via classical trajectories within the semiclassical approximation and which gives the exact, fully quantum mechanical result in the special case of a harmonic bath. The main appeal of the QCPI approach is the removal of the classical memory,^{25} which leads to smaller values of Δ*k*_{max}. An efficient approach is offered by the hierarchical equations of the motion^{26} scheme, which is restricted to bath spectral densities of the simple Debye–Drude form; unfortunately, this method is impractical for simulating the dynamics of systems with typical molecular and phonon modes whose frequencies and coupling parameters cannot be expressed in this simple form.

The present Communication introduces a new decomposition of the path integral which removes the temporal entanglement of path variables even within the memory length. This is achieved through recursive shifting of the entanglement to longer path segments, while simultaneously decreasing its magnitude. This process is repeated until the entanglement terms practically vanish, and when the memory length has also been fully accounted for, the result converges to the full path sum. The disentangled path sum involves multiplication of small *n*^{2} × *n*^{2} matrices (i.e., of size equal to that required to store the density matrix of the bare system), and the construction of these matrices requires a single sum of $(n2)rmax$ terms, where *r*_{max} is the temporal entanglement length, which in practice is equal to Δ*k*_{max}. Thus, the small matrix path integral (SMatPI) decomposition circumvents the exponential scaling of array storage required by earlier iterative path integral methods, while remaining a numerically exact path integral algorithm.

The Hamiltonian describing a system interacting with a dissipative harmonic bath is given by

where *s* is the discrete coordinate of the *n*-state system,

and {*x*_{j}, *p*_{j}} denote the coordinate and momentum variables of the bath degrees of freedom. The focus is on observables pertaining to the system, which may be obtained from the RDM. It is commonly assumed that the initial RDM is a product of system and bath components, *ρ*(0) = *ρ*_{sys}(0)*ρ*_{b}(0).

Consider the reduced propagator of the system at the time *N*Δ*t* (where Δ*t* is the path integral time step), defined as the matrix **U**^{(N0)} with elements

The RDM is easily obtained from Eq. (3),

Provided the path integral time step is sufficiently short, Eq. (3) may be computed exactly from a discretized path integral expression. With the QuAPI splitting,^{27} the discretized path integral has the following form:

where

is the forward–backward propagator of the shifted system Hamiltonian along the adiabatic path, $\u01240=\u0124sys\u2212\u2211jcj2\u015d2/2mj\omega j2$, and *F* is the QuAPI-discretized Feynman–Vernon influence functional, which depends on the coordinates of the system’s discrete forward and backward paths. The discretized influence functional has the following form:^{13}

where the value of the system coordinate $sk+$ at the time *k*Δ*t* on the forward path is $\sigma ik+$, etc., and $\eta k\u2032k\u2033$ are coefficients^{11} related to integrals of the spectral density function,^{28}

The influence functional coefficients may also be obtained directly from time correlation functions computed either quantum mechanically or via molecular dynamics simulations.^{29} If the bath is initially in equilibrium with a particular system site, the influence functional contains additional terms.^{30}

The time nonlocality encoded in the double sum in Eq. (7) [see Fig. 1(a)] entangles the path integral variables, and Eq. (5) requires the full evaluation of a 2^{N−2}-term sum. However, the strength of these nonlocal interactions (i.e., the magnitude of the $k\u2032k\u2032\u2032$ coefficients) decays with increasing time separation Δ*k* = $k\u2032$ − $k\u2032\u2032$ and becomes negligible beyond some Δ*k*_{max}. Taking advantage of this finite-length entanglement, the path sum may be evaluated iteratively^{12} by connecting a rank-Δ*k*_{max} tensor that contains the path segments $sk+\Delta kmax\u22121\xb1,\u2026,sk\xb1$ to $sk+\Delta kmax\xb1,\u2026,sk+1\xb1$. Thus, the i-QuAPI decomposition reduces the exponential scaling with the propagation length *N* to linear scaling. The algorithm requires the storage of $n2\Delta kmax$ path amplitudes at each step, and each iteration involves $n2\Delta kmax+2$ operations; thus, the cost of the algorithm in its simplest form is exponential with the memory length Δ*k*_{max}. Path filtering procedures,^{16–18} singular value decomposition,^{21} and the blip sum^{19} can dramatically reduce the size of these arrays, in particular when Δ*k*_{max} is large. Convergence of the i-QuAPI methodology for two-state systems is fairly straightforward in many regimes of interest, and several applications with *n* > 2 systems have been reported. However, calculations can become very challenging as the number of states increases.

Equation (7) is a product of two-time factors,

The reduced propagator of the system at the first time point is given by the matrix **U**^{(10)} with elements

**U**^{(20)} is given by

where

Equation (11) involves a simple product of two *n*^{2} × *n*^{2} matrices and, thus, requires *n*^{2} operations for each element. Proceeding to the third time step,

This expression contains a double sum that cannot be factored into single sums. Equation (13) shows that the temporal entanglement of path variables is present at *N* = 3 and above (unless only the *η*_{k+1,k} influence functional coefficients are nonzero, in which case the RDM obeys Markovian dynamics). Similarly, the expression for **U**^{(40)} involves a triple sum, which is entangled through Δ*k* = 2 and Δ*k* = 3 influence functional couplings. It will be shown below that the path sum can be disentangled by successively diminishing the magnitude of problematic influence functional connections, while shifting them to larger Δ*k*, until they approach zero.

To illustrate the basic idea, consider first the case of Δ*k*_{max} = 2, which implies that $Fi3\xb1i0\xb1(30)=0$, and thus,

Observing that

Eq. (14) can be rewritten as

where $Ti3\xb1i2\xb1(32)=Ai3\xb1i2\xb1(32)$,

and

The influence functional couplings that are included in these terms are illustrated in Fig. 2. In the first term, the two-step coupling is effectively included through **U**^{(20)}. In the second term, the two-step coupling is in the propagator matrix **T**^{(31)}. These two terms involve only *n*^{2} × *n*^{2} matrix multiplications. The influence functional entanglement is in the last “remainder” term [Eq. (18)], which contains both two-step influence functional couplings. While this term cannot be disentangled, its magnitude is now much reduced compared to Eq. (14). This is so because the influence functional coefficients are small, $\eta k\u2032k\u2033\u226a1$, and, consequently, $Fik\u2032\xb1ik\u2033\xb1(k\u2032k\u2033)$ do not differ significantly from unity.

The goal is to evaluate the RDM for times significantly exceeding the memory length. Continuing to the fourth time step, the full path sum with Δ*k*_{max} = 2 influence functional terms has the following form:

Proceeding as before and using a factorization similar to Eq. (15), it is straightforward to show that Eq. (19) can be written in the following form:

where **T**^{(43)}, **T**^{(42)}, and **T**^{(41)} are one-, two-, and three-step *n*^{2} × *n*^{2} propagator matrices analogous to those defined in Eqs. (17) and (18) and

Again, the decomposition of Eq. (20) has shifted the entanglement to this last term, **T**^{(40)}. Most importantly, the presence of a third small factor in Eq. (21) implies that the magnitude of this term is now further reduced. Continuing this process to the fifth time point, one finds

Here, **T**^{(52)} is a three-step propagator similar to **T**^{(41)}, **T**^{(53)} is a two-step propagator similar to **T**^{(42)}, etc. In fact, the translational invariance of the influence functional coefficients implies that like matrices are identical as long as neither of their indices corresponds to endpoints, and thus, such matrices need to be computed only once. The remainder in Eq. (22) is a five-step propagator, which is even smaller than the four-step remainder given by Eq. (21).

Equations (20) and (22) are just different forms of the full path integral expression for Δ*k*_{max} = 2, i.e., no terms have been left out. This hierarchy continues indefinitely, giving **U**^{(r0)} in terms of small matrix products with a remainder **T**^{(r0)}, which contains the influence functional entanglement. Because the remainder matrices decrease with increasing *r*, there exists a value *r*_{max}, the “temporal entanglement length,” after which they may be dropped with negligible error. This gives the following SMatPI decomposition:

Note that, unlike the tensor decomposition which with Δ*k*_{max} = 2 is based on a two-step object, the SMatPI decomposition has spread the entanglement over objects that span *r*_{max} > 2 time steps.

The most time consuming part in Eq. (23) is the calculation of $T(N,N\u2212rmax)=T(rmax+1,1)$, which involves *a single* (*r*_{max} − 1)-dimensional sum for each of the *n*^{4} elements. Once this matrix has been computed, propagation according to Eq. (23) scales linearly with time and is extremely fast as it involves just ordinary multiplications of *n*^{2} × *n*^{2} matrices. Note that the memory requirements of Eq. (23) are similar to those for propagation of the bare system.

Next, consider the inclusion of Δ*k*_{max} = 3 terms. Returning to Eq. (13), it is clear that the three-step coupling will affect only the entangled remainder in the SMatPI decomposition of Eq. (16). It is not hard to show that the remainder is given by the following expression:

Equation (24) is the most general form for the entangled three-step remainder, even if the influence functional includes terms with Δ*k*_{max} > 3. The dominant term in Eq. (24) is the three-step influence functional memory. Evaluation of the *r*_{max} = 3 remainder involves exactly the same effort as with Δ*k*_{max} = 2. Proceeding as before, one can calculate the contribution of three-step influence functional memory at the *r*_{max} = 4 level. As observed previously, the magnitude of the remainder matrix decreases with increasing *r*_{max}. This trend holds for any value of Δ*k*_{max} > 2.

As already explained, the *r*_{max}-level SMatPI decomposition shifts the entanglement to a “remainder” matrix $T(rmax0)$, which requires the calculation of an expression involving (*r*_{max} − 1) nested sums and whose elements can be made arbitrarily small by increasing the value of *r*_{max}. Clearly, to include influence functional terms with a particular Δ*k*_{max}, one must use *r*_{max} ≥ Δ*k*_{max}. To see what values of *r*_{max} are likely to be needed, consider the following error scaling consideration. Suppose that *F*^{(k+1,k)} influence functional factors are of order 1 − *λ*, that *F*^{(k+2,k)} are of order 1 − *λ*^{2}, etc. According to Eq. (24), the **T**^{(30)} remainder has a contribution of order *λ*^{4} from *F*^{(k+2,k)} terms and another contribution of order *λ*^{3} from the *F*^{(k+3,k)} terms. Truncating the influence functional at Δ*k*_{max} = 2 implies dropping terms of order *λ*^{3}. In that case, the dominant error from neglecting the **T**^{(30)} entangled remainder would be of order *λ*^{4}, thus less significant than the omitted memory terms. Similarly, truncating the influence functional terms at Δ*k*_{max} = 3 implies omitting contributions of order *λ*^{4} from the full path sum. It is easy to show that terminating the SMatPI decomposition at the *r*_{max} = 3 level, i.e., dropping **T**^{(40)}, would lead to error of order *λ*^{5}, which is considerably smaller. One sees that truncation of the SMatPI decomposition at *r*_{max} = Δ*k*_{max} leads to error that is less prominent than that arising from memory truncation. Note that the matrix elements fluctuate considerably, so these scaling criteria indicate a trend rather than strict bounds. If the error arising from the neglect of the remainder with *r*_{max} = Δ*k*_{max} is not as small as desired, one could use *r*_{max} = Δ*k*_{max} + 1. A number of numerical tests indicate that *r*_{max} = Δ*k*_{max} is adequate for accurate RDM propagation.

The numerical calculations presented in Figs. 3–6 corroborate this analysis. Shown in Fig. 3 is the real part of $T12,11(r0)$, one of the largest elements of the **T**_{(r0)} matrices, for a symmetric TLS (*ε* = 0) described by the Hamiltonian,

(which, thus, includes the “counterterms”) with coordinates *σ*_{1} = 1, *σ*_{2} = −1, which is coupled to a dissipative harmonic bath with an Ohmic spectral density,

where Δ*σ* = *σ*_{1} − *σ*_{2}, *ω*_{c} = 5Ω, and *ξ* = 0.3, at a low temperature corresponding to *ℏ*Ω*β* = 5. Before the memory is saturated, $T12,11(r0)$ is seen to be dominated by the magnitude of the influence functional factor for *r* ≤ Δ*k*_{max} and drops abruptly to much smaller values when *r* > Δ*k*_{max}. It is thus clear that the error arising from the neglect of the remainder matrix with *r*_{max} ≃ Δ*k*_{max} is negligible compared to the error arising from premature truncation of influence functional memory. When the system-bath coupling is strong and the memory is truncated prematurely, the value of *r*_{max} that must be used to reproduce the full path integral results for a given Δ*k*_{max} value may slightly exceed Δ*k*_{max}. As the value of Δ*k*_{max} increases to convergence, the errors from the neglected remainder and the neglected memory become comparable and approach zero.

The time evolution of the diagonal and off-diagonal RDM elements, *U*_{11,11} and *U*_{12,11}, with the initial condition $\rho ^sys=11$ (and the bath initially at thermal equilibrium) is displayed in Fig. 4 with ΩΔ*t* = 0.2 and *r*_{max} = Δ*k*_{max} = 6. The SMatPI results are indistinguishable from those obtained via i-QuAPI calculations with the same time step and memory length. In fact, the SMatPI remainder error is too small to be discerned on the population curve. A magnified graph is given in the right panel, which shows the effect of increasing *r*_{max} with full memory and also for unconverged calculations with Δ*k*_{max} = 3. The remainder error with *r*_{max} = Δ*k*_{max} is about 0.001 in both cases and drops to about 0.0001 upon using *r*_{max} = Δ*k*_{max} + 1. If the acceptable error threshold arising from various factors (spectral density integrals, time step, and memory) in these calculations is set to about 0.01, one concludes that the *r*_{max} = Δ*k*_{max} SMatPI results are identical to the full path internal results.

The convergence behavior discussed above and illustrated in Fig. 4 does not seem to depend on the spectral characteristics of the bath. Figure 5 displays the diagonal and off-diagonal elements of the RDM with an initially localized state, i.e., *U*_{11,11} and *U*_{12,11}, for a strongly damped symmetric TLS with *ω*_{c} = Ω and *ξ* = 2 at a high temperature *ℏ*Ω*β* = 0.125, as well as an asymmetric TLS with *ε* = 1, *ω*_{c} = 1.5 Ω and *ξ* = 0.75 at a low temperature, *ℏ*Ω*β* = 5. The TLS level splittings of these two TLSs are 2Ω and $22\Omega $, respectively; thus, both of these models involve sluggish baths which induce long memory. SMatPI calculations were run to Ω*t* = 50 in order to verify stability and norm conservation. The results converged with time steps ΩΔ*t* = 0.125 and 0.15, respectively, and *r*_{max} = Δ*k*_{max} = 20 in both cases. The path sums in the SMatPI calculations did not take advantage of the blip decomposition. An i-QuAPI calculation with 400 time steps would require 400 times the central processing unit (CPU) time of the SMatPI run and storage of an array of 2^{42} ≃ 4.4 × 10^{12} complex-valued elements, an impractical if not impossible task. The population dynamics of the asymmetric TLS exhibits an interesting pattern, which results from the interplay among the bare TLS tunneling frequency, the asymmetry, and the slow-modulating bath at a low temperature and moderate dissipation.

Figure 6 shows the population dynamics of a 10-site dissipative tight binding model specified by the system Hamiltonian,

with *n* = 10, which is again coupled linearly to a harmonic bath with the spectral density given in Eq. (26), where Δ*σ* = *σ*_{1} − *σ*_{n} is the outermost site distance, *ξ* = 0.5, and *ω*_{c} = 8Ω. The bath is initially in equilibrium with the first site at *ℏ*Ω*β* = 1. Converged SMatPI results are presented with a time step ΩΔ*t* = 0.25 and *r*_{max} = Δ*k*_{max} = 5. To test the stability of the SMatPI method again, the site populations were followed for 500 iteration steps, up to Ω*t* = 125. The sum of the populations remained equal to 1 within 8 decimal places over the entire propagation. Figure 5 shows results only up to Ω*t* = 60 as populations have almost reached equilibrium by this time. The populations exhibit complex oscillatory dynamics, with multiple reflections off the edges, interference, and damping due to the dissipative bath. It is seen that the populations of the two edge states approach slightly different values compared to nonterminal sites. This subtle feature is a consequence of the edge-site components in the eigenstates, and the ability of the SMatPI calculation to capture it reflects the high accuracy of the algorithm.

With these parameters, direct computation of the full forward-backward path sum for the same final time has 10^{1002} terms. The i-QuAPI algorithm would require the storage of a tensor with 10^{10} complex-valued amplitudes and 500 × 10^{12} operations and is impractical. The most demanding step of the SMatPI calculation required a single sum of 10^{12} terms and minimal storage (of just 10 100 × 100 matrices) and was performed on the author’s laptop.

In summary, the SMatPI decomposition offers an extremely efficient, numerically exact alternative to iterative path integral simulation of quantum dissipative systems. The i-QuAPI iterative decomposition exploits the finite length of bath-induced memory to replace the full path sum of (*n*^{2})^{N+1} terms by $N(n2)\Delta kmax+1$ operations, along with array storage of size $n2\Delta kmax$. The reduction in the exponent from *N* to Δ*k*_{max}, and the resulting linear scaling with propagation length, makes long-time calculations feasible and has found a plethora of applications. However, the i-QuAPI algorithm maintains the full influence functional entanglement within the memory length, leading to exponential scaling with Δ*k*_{max}. The SMatPI algorithm requires only the storage of a few *n*^{2} × *n*^{2} matrices and a single path sum that involves $(n2)rmax+1$ terms. This dramatic gain is achieved by recursively removing the influence functional entanglement, shifting it to longer paths while simultaneously reducing the magnitude of its contribution, until it can be neglected. The analysis presented here, along with several calculations, suggests that the temporal entanglement length *r*_{max} does not exceed the required Δ*k*_{max}, and thus, the cost of SMatPI is $(n2)\Delta kmax+1$. The SMatPI effort is, thus, equal to the effort required for a single i-QuAPI iteration, and most importantly, no storage of path arrays spanning the memory length is required.

The SMatPI circumvents the exponential scaling of array storage with the memory length of iterative path integral methods and drastically reduces the required CPU time for propagation to long times. This advantage will allow fully quantum mechanical simulations of dissipative dynamics in systems of many states and/or more than one degree of freedom. Finally, the conceptual insights gained by the new picture of the disentangled path integral may prove even more valuable in the quest for understanding the quantum dynamics of systems in dissipative environments.

This material was based on the work supported by the National Science Foundation under Award No. CHE-1665281.