Small matrix disentanglement of the path integral: Overcoming the exponential tensor scaling with memory length

Makri, Nancy

doi:10.1063/1.5139473

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.¹ 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,⁶ 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.⁷

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²)^N terms. Unfortunately, Monte Carlo methods⁸ 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⁹ can sometimes extend calculations to ∼10² 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 {(n^{2})}^{Δ k_{max} + 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 ${(n^{2})}^{Δ k_{max}}$ 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,²⁵ which leads to smaller values of Δk_max. An efficient approach is offered by the hierarchical equations of the motion²⁶ 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.

FIG. 1.

View large Download slide

Left: influence functional connections in a path of length 4Δt. Right: influence functional terms in the rank-3 tensor employed in i-QuAPI with Δk_max = 3.

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² × n² 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 ${(n^{2})}^{r_{max}}$ 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

Ĥ = Ĥ_{sys} + \sum_{j} (\frac{{\hat{p}}_{j}^{2}}{2 m_{j}} + \frac{1}{2} m_{j} ω_{j}^{2} {\hat{x}}_{j}^{2} - c_{j} ŝ {\hat{x}}_{j}),

(1)

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

ŝ = \sum_{i = 1}^{n} |σ_{i}⟩ σ_{i} ⟨σ_{i}|, Ĥ_{sys} = \sum_{i = 1}^{n} \sum_{j = 1}^{n} |σ_{i}⟩ h_{i j} ⟨σ_{j}|,

(2)

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

U_{i_{N}^{\pm} i_{0}^{\pm}}^{(N 0)} = T r_{b} ⟨σ_{i_{N}^{+}}| e^{- i Ĥ N Δ t / ℏ} |σ_{i_{0}^{+}}⟩ {\hat{ρ}}_{b} (0) ⟨σ_{i_{0}^{-}}| e^{i Ĥ N Δ t / ℏ} |σ_{i_{N}^{-}}⟩ .

(3)

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

{\tilde{ρ}}_{i_{N}^{\pm}}^{(N)} \equiv T r_{b} ⟨σ_{i_{N}^{+}}| \hat{ρ} (N Δ t) |σ_{i_{N}^{-}}⟩ = \sum_{i_{0}^{\pm} = 1}^{n} U_{i_{N}^{\pm} i_{0}^{\pm}}^{(N 0)} {\tilde{ρ}}_{i_{0}^{\pm}}^{(0)} .

(4)

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,²⁷ the discretized path integral has the following form:

U_{i_{N}^{\pm} i_{0}^{\pm}}^{(N 0)} = \sum_{i_{N - 1}^{\pm} = 1}^{n} \dots \sum_{i_{1}^{\pm} = 1}^{n} K_{i_{N}^{\pm}, i_{N - 1}^{\pm}} \dots K_{i_{1}^{\pm}, i_{0}^{\pm}} F_{i_{N}^{\pm}, i_{N - 1}^{\pm}, \dots, i_{0}^{\pm}},

(5)

where

K_{i_{k + 1}^{\pm}, i_{k}^{\pm}} = ⟨σ_{i_{k + 1}^{+}}| e^{- i Ĥ_{0} Δ t / ℏ} |σ_{i_{k}^{+}}⟩ ⟨σ_{i_{k}^{-}}| e^{i Ĥ_{0} Δ t / ℏ} |σ_{i_{k + 1}^{-}}⟩

(6)

is the forward–backward propagator of the shifted system Hamiltonian along the adiabatic path, $Ĥ_{0} = Ĥ_{sys} - \sum_{j} {c_{j}}^{2} ŝ^{2} / 2 m_{j} ω_{j}^{2}$ ⁠, 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:¹³

\begin{matrix} F_{i_{N}^{\pm}, i_{N - 1}^{\pm}, \dots, i_{0}^{\pm}} & = & F (s_{N}^{\pm}, \dots, s_{0}^{\pm}) \\ = & exp [- \frac{1}{ℏ} \sum_{k^{'} = 0}^{N} \sum_{k^{″} = 0}^{k^{'}} (s_{k^{'}}^{+} - s_{k^{'}}^{-}) (η_{k^{'} k^{″}} s_{k^{″}}^{+} - η_{k^{'} k^{″}}^{*} s_{k^{″}}^{-})], \end{matrix}

(7)

where the value of the system coordinate $s_{k}^{+}$ at the time kΔt on the forward path is $σ_{i_{k}^{+}}$ ⁠, etc., and $η_{k^{'} k^{″}}$ are coefficients¹¹ related to integrals of the spectral density function,²⁸

J (ω) \equiv \frac{π}{2} \sum_{j} \frac{c_{j}^{2}}{m_{j} ω_{j}} δ (ω - ω_{j}) .

(8)

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

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^{'} k^{''}$ coefficients) decays with increasing time separation Δk = $k^{'}$ − $k^{''}$ and becomes negligible beyond some Δk_max. Taking advantage of this finite-length entanglement, the path sum may be evaluated iteratively¹² by connecting a rank-Δk_max tensor that contains the path segments $s_{k + Δ k_{max} - 1}^{\pm}, \dots, s_{k}^{\pm}$ to $s_{k + Δ k_{max}}^{\pm}, \dots, s_{k + 1}^{\pm}$ ⁠. Thus, the i-QuAPI decomposition reduces the exponential scaling with the propagation length N to linear scaling. The algorithm requires the storage of $n^{2 Δ k_{max}}$ path amplitudes at each step, and each iteration involves $n^{2 Δ k_{max} + 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,²¹ and the blip sum¹⁹ 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,

\begin{matrix} F_{i_{N}^{\pm}, i_{N - 1}^{\pm}, \dots, i_{0}^{\pm}} & = & \prod_{k^{'} = 1}^{N} \prod_{k^{″} = 1}^{k^{'}} F_{i_{k^{'}}^{\pm} i_{k^{″}}^{\pm}}^{(k^{'} k^{″})}, \\ F_{i_{k^{'}}^{\pm} i_{k^{″}}^{\pm}}^{(k^{'} k^{″})} & = & exp [- \frac{1}{ℏ} (σ_{i_{k^{'}}^{+}} - σ_{i_{k^{'}}^{-}}) (η_{k^{'} k^{″}} σ_{i_{k^{″}}^{+}} - η_{k^{'} k^{″}}^{*} σ_{i_{k^{″}}^{-}})] . \end{matrix}

(9)

The reduced propagator of the system at the first time point is given by the matrix U⁽¹⁰⁾ with elements

U_{i_{1}^{\pm} i_{0}^{\pm}}^{(10)} = F_{i_{1}^{\pm}}^{(11)} F_{i_{1}^{\pm} i_{0}^{\pm}}^{(10)} F_{i_{0}^{\pm}}^{(00)} K_{i_{1}^{\pm} i_{0}^{\pm}} .

(10)

U⁽²⁰⁾ is given by

U_{i_{2}^{\pm} i_{0}^{\pm}}^{(20)} = \sum_{i_{1}^{\pm} = 1}^{n} F_{i_{2}^{\pm}}^{(22)} F_{i_{2}^{\pm} i_{1}^{\pm}}^{(21)} F_{i_{2}^{\pm} i_{0}^{\pm}}^{(20)} K_{i_{2}^{\pm} i_{1}^{\pm}} U_{i_{1}^{\pm} i_{0}^{\pm}}^{(10)} = F_{i_{2}^{\pm} i_{0}^{\pm}}^{(20)} \sum_{i_{1}^{\pm} = 1}^{n} A_{i_{2}^{\pm} i_{1}^{\pm}}^{(21)} A_{i_{1}^{\pm} i_{0}^{\pm}}^{(10)},

(11)

where

A_{i_{1}^{\pm} i_{0}^{\pm}}^{(10)} \equiv U_{i_{1}^{\pm} i_{0}^{\pm}}^{(10)}, A_{i_{k + 1}^{\pm} i_{k}^{\pm}}^{(k + 1, k)} \equiv F_{i_{k + 1}^{\pm}}^{(k + 1, k + 1)} F_{i_{k + 1}^{\pm} i_{k}^{\pm}}^{(k + 1, k)} K_{i_{k + 1}^{\pm} i_{k}^{\pm}}, k = 1, \dots .

(12)

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

U_{i_{3}^{\pm} i_{0}^{\pm}}^{(30)} = F_{i_{3}^{\pm} i_{0}^{\pm}}^{(30)} \sum_{i_{2}^{\pm} = 1}^{n} \sum_{i_{1}^{\pm} = 1}^{n} F_{i_{3}^{\pm} i_{1}^{\pm}}^{(31)} F_{i_{2}^{\pm} i_{0}^{\pm}}^{(20)} A_{i_{3}^{\pm} i_{2}^{\pm}}^{(32)} A_{i_{2}^{\pm} i_{1}^{\pm}}^{(21)} A_{i_{1}^{\pm} i_{0}^{\pm}}^{(10)} .

(13)

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⁽⁴⁰⁾ 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 $F_{i_{3}^{\pm} i_{0}^{\pm}}^{(30)} = 0$ ⁠, and thus,

U_{i_{3}^{\pm} i_{0}^{\pm}}^{(30)} = \sum_{i_{2}^{\pm} = 1}^{n} \sum_{i_{1}^{\pm} = 1}^{n} F_{i_{3}^{\pm} i_{1}^{\pm}}^{(31)} F_{i_{2}^{\pm} i_{0}^{\pm}}^{(20)} A_{i_{3}^{\pm} i_{2}^{\pm}}^{(32)} A_{i_{2}^{\pm} i_{1}^{\pm}}^{(21)} A_{i_{1}^{\pm} i_{0}^{\pm}}^{(10)} .

(14)

Observing that

F_{i_{3}^{\pm} i_{1}^{\pm}}^{(31)} F_{i_{2}^{\pm} i_{0}^{\pm}}^{(20)} = F_{i_{2}^{\pm} i_{0}^{\pm}}^{(20)} + (F_{i_{3}^{\pm} i_{1}^{\pm}}^{(31)} - 1) + (F_{i_{3}^{\pm} i_{1}^{\pm}}^{(31)} - 1) (F_{i_{2}^{\pm} i_{0}^{\pm}}^{(20)} - 1),

(15)

Eq. (14) can be rewritten as

U_{i_{3}^{\pm} i_{0}^{\pm}}^{(30)} = \sum_{i_{2}^{\pm} = 1}^{n} T_{i_{3}^{\pm} i_{2}^{\pm}}^{(32)} U_{i_{2}^{\pm} i_{0}^{\pm}}^{(20)} + \sum_{i_{1}^{\pm} = 1}^{n} T_{i_{3}^{\pm} i_{1}^{\pm}}^{(31)} U_{i_{1}^{\pm} i_{0}^{\pm}}^{(10)} + T_{i_{3}^{\pm} i_{0}^{\pm}}^{(30)},

(16)

where $T_{i_{3}^{\pm} i_{2}^{\pm}}^{(32)} = A_{i_{3}^{\pm} i_{2}^{\pm}}^{(32)}$ ⁠,

T_{i_{3}^{\pm} i_{1}^{\pm}}^{(31)} = (F_{i_{3}^{\pm} i_{1}^{\pm}}^{(31)} - 1) \sum_{i_{2}^{\pm} = 1}^{n} A_{i_{3}^{\pm} i_{2}^{\pm}}^{(32)} A_{i_{2}^{\pm} i_{1}^{\pm}}^{(21)},

(17)

and

T_{i_{3}^{\pm} i_{0}^{\pm}}^{(30)} = \sum_{i_{2}^{\pm} = 1}^{n} \sum_{i_{1}^{\pm} = 1}^{n} (F_{i_{3}^{\pm} i_{1}^{\pm}}^{(31)} - 1) (F_{i_{2}^{\pm} i_{0}^{\pm}}^{(20)} - 1) A_{i_{3}^{\pm} i_{2}^{\pm}}^{(32)} A_{i_{2}^{\pm} i_{1}^{\pm}}^{(21)} A_{i_{1}^{\pm} i_{0}^{\pm}}^{(10)} .

(18)

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⁽²⁰⁾. In the second term, the two-step coupling is in the propagator matrix T⁽³¹⁾. These two terms involve only n² × n² 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, $|η_{k^{'} k^{″}}| ≪ 1$ ⁠, and, consequently, $|F_{i_{k^{'}}^{\pm} i_{k^{″}}^{\pm}}^{(k^{'} k^{″})}|$ do not differ significantly from unity.

FIG. 2.

View large Download slide

Path integral couplings in the SMatPI decomposition for N = 3. Top: T⁽³²⁾ · U⁽²⁰⁾. Middle: T⁽³¹⁾ · U⁽¹⁰⁾. Bottom: T⁽³⁰⁾. Solid lines indicate the full path integral couplings given by $A^{(k k^{'})}$ or $F^{(k k^{'})}$ ⁠, while dashed lines indicate reduced couplings $F^{(k k^{'})} - 1$ ⁠. Terms included in T matrices are shown in red, while those in black correspond to U matrices.

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:

U_{i_{4}^{\pm} i_{0}^{\pm}}^{(40)} = \sum_{i_{3}^{\pm} = 1}^{n} \sum_{i_{2}^{\pm} = 1}^{n} \sum_{i_{1}^{\pm} = 1}^{n} F_{i_{4}^{\pm} i_{2}^{\pm}}^{(42)} F_{i_{3}^{\pm} i_{1}^{\pm}}^{(31)} F_{i_{2}^{\pm} i_{0}^{\pm}}^{(20)} A_{i_{3}^{\pm} i_{2}^{\pm}}^{(43)} A_{i_{3}^{\pm} i_{2}^{\pm}}^{(32)} A_{i_{2}^{\pm} i_{1}^{\pm}}^{(21)} A_{i_{1}^{\pm} i_{0}^{\pm}}^{(10)} .

(19)

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:

\begin{matrix} U_{i_{4}^{\pm} i_{0}^{\pm}}^{(40)} = \sum_{i_{3}^{\pm} = 1}^{n} T_{i_{4}^{\pm} i_{3}^{\pm}}^{(43)} U_{i_{3}^{\pm} i_{0}^{\pm}}^{(30)} + \sum_{i_{2}^{\pm} = 1}^{n} T_{i_{4}^{\pm} i_{2}^{\pm}}^{(42)} U_{i_{2}^{\pm} i_{0}^{\pm}}^{(20)} + \sum_{i_{1}^{\pm} = 1}^{n} T_{i_{4}^{\pm} i_{1}^{\pm}}^{(41)} U_{i_{1}^{\pm} i_{0}^{\pm}}^{(10)} + T_{i_{4}^{\pm} i_{0}^{\pm}}^{(40)}, \end{matrix}

(20)

where T⁽⁴³⁾, T⁽⁴²⁾, and T⁽⁴¹⁾ are one-, two-, and three-step n² × n² propagator matrices analogous to those defined in Eqs. (17) and (18) and

\begin{matrix} T_{i_{4}^{\pm} i_{0}^{\pm}}^{(40)} & = & \sum_{i_{3}^{\pm} = 1}^{n} \sum_{i_{2}^{\pm} = 1}^{n} \sum_{i_{1}^{\pm} = 1}^{n} (F_{i_{4}^{\pm} i_{2}^{\pm}}^{(42)} - 1) (F_{i_{3}^{\pm} i_{1}^{\pm}}^{(31)} - 1) (F_{i_{2}^{\pm} i_{0}^{\pm}}^{(20)} - 1) \\ \times A_{i_{4}^{\pm} i_{3}^{\pm}}^{(43)} A_{i_{3}^{\pm} i_{2}^{\pm}}^{(32)} A_{i_{2}^{\pm} i_{1}^{\pm}}^{(21)} A_{i_{1}^{\pm} i_{0}^{\pm}}^{(10)} . \end{matrix}

(21)

Again, the decomposition of Eq. (20) has shifted the entanglement to this last term, T⁽⁴⁰⁾. 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

\begin{matrix} U_{i_{5}^{\pm} i_{0}^{\pm}}^{(50)} & = & \sum_{i_{3}^{\pm} = 1}^{n} T_{i_{5}^{\pm} i_{4}^{\pm}}^{(54)} U_{i_{4}^{\pm} i_{0}^{\pm}}^{(40)} + \sum_{i_{3}^{\pm} = 1}^{n} T_{i_{5}^{\pm} i_{3}^{\pm}}^{(53)} U_{i_{3}^{\pm} i_{0}^{\pm}}^{(30)} + \sum_{i_{2}^{\pm} = 1}^{n} T_{i_{5}^{\pm} i_{2}^{\pm}}^{(52)} U_{i_{2}^{\pm} i_{0}^{\pm}}^{(20)} \\ + \sum_{i_{1}^{\pm} = 1}^{n} T_{i_{5}^{\pm} i_{1}^{\pm}}^{(51)} U_{i_{1}^{\pm} i_{0}^{\pm}}^{(10)} + T_{i_{5}^{\pm} i_{0}^{\pm}}^{(50)} . \end{matrix}

(22)

Here, T⁽⁵²⁾ is a three-step propagator similar to T⁽⁴¹⁾, T⁽⁵³⁾ is a two-step propagator similar to T⁽⁴²⁾, 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:

U^{(N 0)} = \sum_{m = N - r_{max}}^{N - 1} T^{(N m)} \cdot U^{(m 0)}, N = r_{max} + 1, \dots .

(23)

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 - r_{max})} = T^{(r_{max} + 1,1)}$ ⁠, which involves a single (r_max − 1)-dimensional sum for each of the n⁴ 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² × n² 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:

\begin{matrix} T_{i_{3}^{\pm} i_{0}^{\pm}}^{(30)} & = & \sum_{i_{2}^{\pm} = 1}^{n} \sum_{i_{1}^{\pm} = 1}^{n} [(F_{i_{3}^{\pm} i_{1}^{\pm}}^{(31)} - 1) (F_{i_{2}^{\pm} i_{0}^{\pm}}^{(20)} - 1) + (F_{i_{3}^{\pm} i_{0}^{\pm}}^{(30)} - 1) F_{i_{3}^{\pm} i_{1}^{\pm}}^{(31)} F_{i_{2}^{\pm} i_{0}^{\pm}}^{(20)}] \\ \times A_{i_{3}^{\pm} i_{2}^{\pm}}^{(32)} A_{i_{2}^{\pm} i_{1}^{\pm}}^{(21)} A_{i_{1}^{\pm} i_{0}^{\pm}}^{(10)} . \end{matrix}

(24)

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^{(r_{max} 0)}$ ⁠, 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 − λ², etc. According to Eq. (24), the T⁽³⁰⁾ remainder has a contribution of order λ⁴ from F^(k+2,k) terms and another contribution of order λ³ from the F^(k+3,k) terms. Truncating the influence functional at Δk_max = 2 implies dropping terms of order λ³. In that case, the dominant error from neglecting the T⁽³⁰⁾ entangled remainder would be of order λ⁴, thus less significant than the omitted memory terms. Similarly, truncating the influence functional terms at Δk_max = 3 implies omitting contributions of order λ⁴ 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⁽⁴⁰⁾, would lead to error of order λ⁵, 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 $T_{12,11}^{(r 0)}$ ⁠, one of the largest elements of the T_(r0) matrices, for a symmetric TLS (ε = 0) described by the Hamiltonian,

Ĥ_{0} = - ℏ Ω (|σ_{1}⟩ ⟨σ_{2}| + |σ_{2}⟩ ⟨σ_{1}|) + ε (|σ_{1}⟩ ⟨σ_{1}| - |σ_{2}⟩ ⟨σ_{2}|),

(25)

(which, thus, includes the “counterterms”) with coordinates σ₁ = 1, σ₂ = −1, which is coupled to a dissipative harmonic bath with an Ohmic spectral density,

J (ω) = \frac{2 π ℏ ω}{Δ σ^{2}} ξ e^{- ω / ω_{c}},

(26)

where Δσ = σ₁ − σ₂, ω_c = 5Ω, and ξ = 0.3, at a low temperature corresponding to ℏΩβ = 5. Before the memory is saturated, $T_{12,11}^{(r 0)}$ 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.

FIG. 3.

View large Download slide

Real part of $T_{12,11}^{(r 0)}$ as a function of r for Δk_max = 2, 3, …, 8. The matrix element for r = 2 has the same value for all Δk_max ≥ 2. With Δk_max = 2 (brown circles), the remainder matrix element drops abruptly by two orders of magnitude as r is increased from 2 to 3, in accordance with Eq. (18). Inclusion of the Δk_max = 3 influence functional factor (orange circles) leads to a large increase in the remainder matrix element at r = 3, but a similar abrupt decrease is observed again as r ≥ 4. The same trend is observed with Δk_max = 4–8 (blue, violet, light green, dark green, and red circles).

FIG. 4.

View large Download slide

Density matrix elements for a symmetric TLS with parameters described in the text, obtained via SMatPI (markers) and also through i-QuAPI (lines). Left: Diagonal (red) and off-diagonal (green and blue) elements with Δk_max = 6 and i-QuAPI results with the same time step and memory length. Red solid circles show the diagonal element with $r_{max}$ = Δ $k_{max}$ = 6. Green and blue hollow squares show the real and imaginary parts of the off-diagonal element with Δk_max = 6 and r_max = 6, while green and blue solid circles show the off-diagonal element with Δk_max = 6 and r_max = 7. The right panel shows magnified results of the diagonal element using SMatPI with Δk_max = 6 and r_max = 6 and 8 (hollow green and solid red squares) vs full i-QuAPI results (red line), as well as unconverged results using SMatPI with Δk_max = 3 and r_max = 3, 4, 5, 6 (hollow orange, hollow violet, hollow green, and solid blue circles) vs full i-QuAPI results with the same memory length (blue line).

FIG. 5.

View large Download slide

Reduced density matrix elements for a symmetric and an asymmetric TLS with parameters described in the text. Converged SMatPI results obtained with r_max = Δk_max = 20. The lines show results for the symmetric TLS, while the lines with markers show results for the asymmetric model. Black: U_11,11. Red and blue: real and imaginary parts of U_12,11.

FIG. 6.

View large Download slide

Site populations for the tight binding Hamiltonian of Eq. (27) with parameters described in the text, computed via SMatPI with r_max = Δk_max = 5. Initially only the first site is populated.

The time evolution of the diagonal and off-diagonal RDM elements, U_11,11 and U_12,11, with the initial condition ${\hat{ρ}}_{sys} = |1⟩ ⟨1|$ (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 $2 \sqrt{2} Ω$ ⁠, 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⁴² ≃ 4.4 × 10¹² 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,

Ĥ_{0} = - ℏ Ω \sum_{i = 1}^{n - 1} (|σ_{i}⟩ ⟨σ_{i + 1}| + |σ_{i + 1}⟩ ⟨σ_{i}|),

(27)

with n = 10, which is again coupled linearly to a harmonic bath with the spectral density given in Eq. (26), where Δσ = σ₁ − σ_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¹⁰⁰² terms. The i-QuAPI algorithm would require the storage of a tensor with 10¹⁰ complex-valued amplitudes and 500 × 10¹² operations and is impractical. The most demanding step of the SMatPI calculation required a single sum of 10¹² 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²)^N+1 terms by $N {(n^{2})}^{Δ k_{max} + 1}$ operations, along with array storage of size ${n^{2}}^{Δ k_{max}}$ ⁠. 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² × n² matrices and a single path sum that involves ${(n^{2})}^{r_{max} + 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 ${(n^{2})}^{Δ k_{max} + 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.