Blip decomposition of the path integral: Exponential acceleration of real-time calculations on quantum dissipative systems

Dissipative processes are ubiquitous in chemistry and are known to affect the dynamics of condensed phase systems in profound ways. While isolated systems often exhibit pronounced quantum mechanical behavior such as rephasing or tunneling, interaction with a dissipative medium generally tends to quench these effects, eventually leading the system to the state of thermal equilibrium (for example, see Refs. 1–10 and references therein.) (However, driven dissipative systems display more complex behaviors; quantum stochastic resonance^11–13 is a fascinating phenomenon where dissipation can enhance coherent quantum oscillations.). Generic effects in the dynamics of quantum dissipative processes are often studied in terms of system-bath models, which consist of a low-dimensional system interacting with a “bath” of harmonic oscillators. Specifically, the “spin-boson” model, which consists of a two-level system (TLS) bilinearly coupled to a harmonic bath, is the simplest model of quantum dissipation and has received much attention. Apart from their significance as generic models for elucidating the effects of condensed phase environments, system-bath Hamiltonians provide usefully accurate descriptions of phenomena induced by lattice phonon vibrations, as well as molecular reactions^14,15 and charge transfer in solution or in biological molecules.^16,17 Even though no adequate harmonic modes can be identified at the molecular level in the latter case, the validity of a harmonic bath model¹⁸ is commonly attributed to linear response; a formal proof based on the path integral formulation is available,¹⁹ and recent mean-field calculations in crystals led to consistent findings.²⁰ 

The dissipative effects of a harmonic bath on tunneling dynamics have been studied via a variety of analytical and numerical treatments, many of which are based on the path integral formulation of quantum mechanics.^21,22 The first comprehensive study of TLS tunneling used the non-interacting blip approximation⁴ to perform an approximate summation of the path integral, with the dissipative effects of the harmonic bath included via the Feynman-Vernon influence functional.¹ Early numerical simulations using Monte Carlo sampling^23,24 to evaluate the multidimensional path integral were able to capture the short-time behavior and relaxation dynamics of dissipative tunneling systems.^25–29 The development of quasi-adiabatic propagators^30–33 and numerically exact iterative decomposition schemes,^34–45 which do not suffer from the Monte Carlo sign problem,^46,47 extended the feasibility of real-time path integral calculations to long times. Several other path integral algorithms and hybrid methods have also been reported in recent years (see, for example, Refs. 48–54).

The present paper focuses on processes characterized by strong coupling to sluggish dissipative media (e.g., solvent or low-frequency protein vibrations), which leads to quenching of coherence along with a near-monotonic (although not necessarily exponential) decay to equilibrium. These conditions pose a challenge to path integral methods, because they dictate a small time step, while coupling to a slow bath implies³⁴ long-lived “memory,” leading to a substantial increase in computation cost. In particular, when the length of the bath-induced memory is comparable to the entire time length of interest, iterative decompositions are no longer meaningful. When this time spans several tens or hundreds of time steps, evaluation of the path integral poses an extremely challenging task.

Under the stated conditions, it is clear that the structure of conventional path integral schemes, where the main building block is the full path sum (which accounts for all coherence effects of the isolated system) may not offer the optimal starting point for such strongly dissipative, near-incoherent processes. Instead, one expects that a fully incoherent approximation to the dynamics should provide a better framework that should lead to faster convergence. Such a fully incoherent limit is easy to construct by collapsing the forward and backward branches of the path sum for the reduced density matrix into a single path. Within the framework of semiclassical theory, such a path forms the basis of forward-backward semiclassical and linearized semiclassical treatments,^55–66 which have been shown to lead to excellent approximations of the dynamics of liquids where quantum coherence effects are drastically quenched.^67–72 In the present case though, the path sum is not restricted to classical trajectories, but includes all discretized Feynman paths along the forward-backward time contour.

The connection between this single-path incoherent limit and the full path integral expression is easily established by examining the decoherence-inducing role of the influence functional,¹⁰ which is manifested through blips, i.e., path segments where the forward and backward branches do not coincide.⁴ Thus, the fully incoherent result is the zero-blip limit of the path integral expression. Starting from this limit, one gradually includes coherence effects by adding path pairs with one, two,^73–75 or more blips. For processes that display mostly incoherent behavior, one intuitively expects this strategy to be much more efficient than the conventional approach which starts with all distinct forward and backward paths (whose phases account for the coherence of the isolated quantum system in its entirety) and subsequently lets the influence functional wipe out this coherence by annihilating the contribution of the vast majority of these paths.

Yet, the anticipated rapid convergence is not the most notable advantage of the blip decomposition. Section IV shows that expression of the forward-backward path integral in terms of sum and difference coordinates maps entire segments of the inner sum on the time-local dynamics of an externally driven Hamiltonian, which allows these portions to be evaluated by ordinary matrix procedures, reducing the computational effort dramatically. Even in the limit where all possible path pairs are included, the blip decomposition leads to an exponential decrease of the number of operations compared to the conventional procedure.

Section II reviews the discretized path integral formulation for a system coupled to a harmonic dissipative bath. The blip decomposition of the full path integral, along with the mapping that allows iterative evaluation of blip-free segments, is described in Sec. III. Section IV discusses the rapid convergence of the blip series, which leads to dramatic savings in path integral calculations of strongly damped processes. Representative calculations on a dissipative TLS are shown in Sec. V. Section VI presents some concluding remarks.

Grid representations of discretized path integrals³⁰ offer the possibility of deterministic (quadrature-based) evaluation, thus potentially circumventing the “sign problem” that plagues Monte Carlo methods in real time.⁴⁷ Earlier work has shown³³ that discrete variable representations^76–79 (DVR) provide the most efficient discretization of the path integral because they require the same number of terms as the number of populated system eigenstates, yet they are local in position space, allowing straightforward evaluation of the influence functional. Further, truncation of the system's Hilbert space eliminates the rapid oscillations of the real-time propagator,^32,80 leading to well-behaved integrands. The Hamiltonian describing the system is thus written in the form

where |σ_n⟩ represent spin-type states or a DVR obtained from continuous coordinates, i.e., the system position operator is

It is assumed that the system interacts with a harmonic oscillator bath via bilinear coordinate coupling terms, such that the combined bath and interaction Hamiltonian has the form

The characteristics of the bath pertaining to the dynamics of the system can be specified collectively via the spectral density function³ 

System observables can be obtained from knowledge of the reduced density matrix

where |$\hat H = \hat H_0 + \hat H_{{\rm sb}} $|$Ĥ=Ĥ0+Ĥ sb$ is the total Hamiltonian and the trace is with respect to the bath. It is commonly assumed⁴ (though not necessary in the quasi-adiabatic propagator formulation^41,42) that the initial density operator factorizes and that the isolated bath is initially in equilibrium at a temperature 1/k_Bβ, where k_B is the Boltzmann constant, i.e.,

where H_b is the Hamiltonian for the uncoupled bath (obtained from Eq. (2.3) by setting c_j = 0) and Z_b is the canonical partition function of the bath. Expressing each of the forward and backward time evolution operators in discretized path integral form²² with time step Δt = t/N and integrating out the harmonic bath degrees of freedom leads to the following form for the reduced density matrix:

Here F is the DVR-discretized³³ Feynman-Vernon influence functional,¹ 

The short time propagators in Eq. (2.7) are obtained using the quasiadiabatic splitting.³⁰ The coefficients η_k′k″ have been given in previous publications^36,81 in terms of integrals of the spectral density function.

The double sum in Eq. (2.8), which couples non-adjacent time points, introduces interactions between path integral coordinates separated by more than a single time step. These nonlocal in time effects prevent stepwise evaluation of the path integral by means of matrix-vector multiplication,⁸² a procedure routinely applied (in various forms⁸³) to the propagation of low-dimensional systems in the absence of an influence functional. The coefficients η_kk′ responsible for these nonlocal influence functional couplings are discretized forms of the bath force autocorrelation function.¹ The critical observation that enables an iterative evaluation of Eq. (2.7) is that decoherence leads to irreversible decay in the correlation function from a bath with a continuous spectral density.³⁴ When the time interval (k′ − k″)Δt separating the path integral variables |$\sigma _{k^{\prime} }^ \pm $|$σk′±$ and |$\sigma _{k^{\prime \prime} }^ \pm $|$σk′′±$ exceeds the time over which the bath correlation function is nonzero, these variables become decoupled and need not be summed simultaneously. Inclusion of the time points that span the bath correlation (or “memory”) time leads to multi-time propagators,^34–36 which can be thought of as tensors³⁵ (with respect to system coordinates) or ordinary matrices³⁷ (with respect to discrete path segments).

The iterative quasi-adiabatic propagator path integral (i-QuAPI) methodology summarized above is a deterministic (i.e., free of statistical error) procedure which converges to the full discretized path integral result. The only input required is the system Hamiltonian (either in matrix form or in terms of a continuous potential that can converted to matrix form via a DVR transformation) and the bath spectral density,³ which may be specified numerically on a frequency grid. Thus, i-QuAPI offers a simple and efficient way of propagating the reduced density matrix, provided the bath induced memory is not very long (i.e., Δk_max ∼ 10). This is often the case when the bath includes a continuum of modes with frequencies that are high compared to the frequency components of the system. On the other hand, slow baths strongly coupled to the system often lead to memory that spans several tens of time steps. Direct implementation of the procedure described above becomes impractical in these cases. Progress is made by realizing that the vast majority of multi-time path segments that may span the memory length make exponentially small contributions. Filtering procedures^37,38,43,44 discard the vast majority of path segments, retaining only those that are statistically significant. Such filtering techniques extend significantly the applicability of the i-QuAPI methodology to multi-level systems coupled to relatively slow solvents. The i-QuAPI methodology has produced useful benchmark results for spin-boson dynamics and barrier crossing and been applied successfully to test the validity of analytic and numerical approximations, to investigate the interplay between tunneling, coherence and dissipation, and to study a host of proton, charge and energy transfer processes (for example, see Refs. 84–100). Its extension to include time-dependent fields is straightforward and has been used to investigate tunneling control.^101–105 The methodology has also been generalized to equilibrium (non-product) initial conditions and thus can be used to calculate a variety of time correlation functions.^41,42 More recently, the methodology was extended to include multiple bosonic and fermionic reservoirs.^106,107

Baths composed of strongly coupled fast and slow degrees of freedom present more challenging situations, as the high frequency oscillators enforce a small time step while the slow modes produce long-lived memory. An attractive way to proceed is renormalization of the propagator.¹⁰⁸ The main idea is to keep increasing the path integral time step successively, while including longer memory in each step of the process. While the renormalization procedure is not an exact decomposition of the path sum, even its crudest two-time version often leads to nearly quantitative results. The algorithm can be made arbitrarily accurate by augmenting the scheme to a multitime propagator.

The focus of the present paper is on situations where the influence functional memory is comparable to the entire time length of interesting dynamics, prohibiting iterative decompositions of the path integral. Section III presents an exact decomposition of the path sum that reduces the number of paths that need to be included explicitly without involving any additional approximations. Further, Sec. IV shows that paths can be classified into groups of decreasing weight based on the number of blips they contain, and that discarding exponentially small multi-blip contributions allows a dramatic increase of efficiency.

For the diagonal element of the reduced density matrix (⁠|$\sigma _N^ + = \sigma _N^ - \equiv \sigma _N $|$σN+=σN−≡σN$⁠) the number of terms in the double path sum of Eq. (2.7) is M^2N. By transforming to sum and difference coordinates,

the double path sum is re-expressed as

As a result of this change of integration variables, the inner summation index |$\bar \sigma _k $|$σ¯k$ at each time point in the sum-difference coordinate representation of the path integral depends on the outer summation index Δσ_k. For example, in the case of a TLS with coordinate values σ = ±1, the sum/difference coordinate pair takes on the values |$\Delta \sigma = \pm 2,\;\bar \sigma = 0$|$Δσ=±2,σ¯=0$ and |$\Delta \sigma = 0,\;\bar \sigma = \pm 1$|$Δσ=0,σ¯=±1$⁠. Time points kΔt with Δσ_k ≠ 0 are the “blips,” while those with Δσ_k = 0 are the “sojourns.”⁴ Thus the outer sum in Eq. (3.2) contains all forward-backward path pairs with b = 0, 1, 2, …, N blips. For each path pair, sojourn time points have identical forward and backward coordinates. Eq. (3.2) shows that sojourns do not interact with other sojourns (except with their immediate neighbors through the short-time system propagators). Yet, explicit summation with respect to the coordinates of all sojourn points gives rise to a very large number of terms. The absence of nonlocal couplings for the sojourn coordinates suggests that these terms can be summed more efficiently.

To illustrate the idea, consider the diagonal element of the reduced density matrix, |$\left\langle 1 \right|\tilde \rho (N\Delta t)\left| 1 \right\rangle $|$1ρ̃(NΔt)1$⁠, for a TLS with site coordinates σ = ±1 prepared initially in a diagonal state. The zero-blip (b = 0), fully incoherent limit is given by the sum of all sojourn paths, which with the given endpoints has the form

where the second equality follows by setting |$\sigma _k^ + = \sigma _k^ - \equiv \sigma _k $|$σk+=σk−≡σk$⁠. (Note that in the case of an off-diagonal initial condition, the blip at the initial coordinate |$\sigma _0^ \pm $|$σ0±$ must be allowed.) Defining the incoherent propagator

which in the case of a TLS is a 2 × 2 matrix, allows Eq. (3.3) to be evaluated via a sequence of N matrix-vector multiplications. Note that the incoherent propagator does not satisfy the group property, therefore the result of the fully incoherent path sum depends on the chosen time step. For a symmetric TLS described by

the incoherent propagator matrix is

The eigenvalues of this matrix are 1 and λ = cos ²ΩΔt − sin ²ΩΔt. Since 0 < λ < 1 for |$\Omega \Delta t < \pi {\kern 1pt} /{\kern 1pt} 4$|$ΩΔt<π/4$⁠, repeated operation with this propagator leads to exponential decay of the initially localized state to equilibrium. Note, however, that the incoherent propagator does not depend on the temperature, thus the final state reached generally is not the proper thermal equilibrium state. Path pairs with blips must be included to correct the time step-dependent, temperature-independent fully incoherent limit. As argued in Sec. IV, the zero-blip approximation provides a good starting point under conditions that lead to near-monotonic relaxation.

The one-blip (b = 1) correction to the fully incoherent approximation corresponds to terms in Eq. (3.2) that are made up of two sojourn segments connected to the single blip. The latter may be located at any of the time points k₁ = 0, …, N − 1. Since the forward and backward paths are identical at all but one of the path integral time points, the contribution to the b = 1 correction from each of the two configurations (⁠|$\sigma _{k_1 }^ + = 1,\,\;\sigma _{k_1 }^ - = - 1$|$σk1+=1,σk1−=−1$ and |$\sigma _{k_1 }^ + = - 1,\,\;\sigma _{k_1 }^ - = 1$|$σk1+=−1,σk1−=1$⁠) that have a blip at time point k₁ is given by the expression

Figure 1(a) shows a particular realization of the sojourn configurations for a forward-backward path pair with N = 20 time steps and a single blip located at k₁ = 3 with Δσ₃ = 2. The solid line shows the single blip-blip self-interaction (which corresponds to the real part of the exponential in Eq. (3.7)), while the dashed lines show blip-sojourn interactions (which correspond to the imaginary terms). In addition to the influence functional interactions, each time point is coupled to its nearest neighbors via the short-time system propagators (not shown in the figure). According to Eq. (3.7), one must sum the contributions from all sojourn paths. For the given blip configuration, the sojourn coordinates can take two values at each of the time points 0, 1, …, k₁ − 1, k₁ + 1, …, N − 1, giving rise to 2¹⁹ terms for the chosen path length N = 20. Fortunately, this sum can be evaluated with very little effort. Because influence functional interactions are zero unless the later time point is a blip, the sum with respect to all sojourn paths that span the points k₁ + 1, …, N − 1 may be evaluated by straightforward matrix-vector multiplication as in the zero-blip case. Next, the sojourn paths that span the time points 0, 1, …, k₁ − 1 contain in addition an influence functional phase that couples each of these time points to the coordinates of the blip point. This sum is of the form

where |$w_{k_1 k^{\prime \prime} } = - 2\Delta \sigma _{k_1 } {\mathop{\rm Im}\nolimits} \eta _{k_1 k^{\prime \prime} }$|$wk1k′′=−2Δσk1 Im ηk1k′′$⁠. If the standard short time propagator is replaced by Eq. (3.4), the sojourn sum at hand becomes isomorphic to the path integral expression for the evolution of the wavefunction (not density matrix) of a system augmented by an external driving term. Thus, defining the single-blip propagator matrix,

the sum with respect to all sojourn segments can be obtained by sequential matrix-vector multiplications.

Next, consider the particular term in Eq. (3.2) where the outer sum contains two blips at the time points k₁ and k₂. A particular realization of the sojourn path segments for k₁ = 3, |$\Delta \sigma _{k_1 } = 2$|$Δσk1=2$ and k₂ = 15, |$\Delta \sigma _{k_2 } = - 2$|$Δσk2=−2$ is shown in Figure 1(b), along with all influence functional interactions for this path pair. For the given blip configuration, the inner sum corresponds to all possible coordinates of the remaining 18 sojourn time points, i.e., 2¹⁸ terms. The goal is to sum the contributions from all sojourn paths (⁠|$\sigma _k^ + = \sigma _k^ - $|$σk+=σk−$⁠) that can be constructed with the chosen blip configuration. Again, the sum with respect to the path coordinates of time points k₂ + 1, …, N − 1 contains no influence functional terms and can be evaluated iteratively using the incoherent propagator matrix. For the case shown in Fig. 1(b), this procedure replaces a sum over 2⁴ terms by four 2 × 2 matrix-vector multiplications. The sum with respect to the sojourn paths that span the points k₁ + 1, …, k₂ − 1 contains influence functional phase factors that couple these points to the blip at k₂, thus that segment can be evaluated iteratively by using the single-blip propagator |$Q_1^{k_2 } $|$Q1k2$⁠, replacing 2¹¹ terms by eleven sequential 2 × 2 matrix-vector multiplications.

Finally, the coordinates of the sojourn segments that correspond to time points 0, 1, …, k₁ − 1 are coupled to both blips. The relevant sum has the form

which can be evaluated iteratively as well, using a two-blip propagator with an external phase from two source points

This last step is performed with three matrix-vector operations rather than the full sum of 2³ terms. Thus, mapping blip-sojourn interactions to a local Hamiltonian with external time-dependent interactions from the blips reduces the number of terms for this blip configuration from 2¹⁸ to 18 multiplications of a 2 × 2 matrix and a 2-component vector, a total of only 18 × 2² terms.

The above ideas are easily extended to all possible arrangements of the outer sum of Eq. (3.2) (blip configurations). For each such arrangement, the iterative evaluation of the inner (sojourn segment) sum involves on average |${\textstyle{1 \over 2}}N$|$12N$ matrix-vector multiplications. Since there are three possible values of Δσ for each of the N time points in the case of a TLS, the total number of outer sum arrangements is 3^N. Thus the blip decomposition of the path sum described in this section reduces the number of operations to |${\textstyle{1 \over 2}}3^N N$|$123NN$⁠. While the scaling is still exponential with the number of time slices, it is clear that the blip decomposition represents very substantial savings when compared to the full evaluation of the path sum, which requires evaluation of 2^2N terms. It is worth stressing that the blip decomposition is an exact transformation which involves no approximations. The procedure can be straightforwardly extended to multistate Hamiltonians.

In the absence of a dissipative environment, a system with a quantized level structure characterized by two or more distinct frequencies exhibits effects that are collectively referred to as “quantum coherence.” The interference fringes in Young's two-slit experiment, along with rephasing of time correlation functions, are most evident examples, but quantum tunneling is also intimately connected to quantum coherence. In time-dependent semiclassical theory, it is well known^109,110 that these processes arise from phase interference between distinct forward and backward trajectories of the system. This picture applies also within the path integral framework. Pairs of identical forward-backward quantum paths cannot contribute to quantum coherence because their complex conjugate phases multiply to unity.

As discussed in Sec. III, the quenching of coherence induced through interaction with a dissipative environment involves the damping of the quantum phase of forward-backward path pairs. This mechanism is present only if the forward and backward paths differ at least somewhat, as identical forward-backward paths make purely incoherent contributions. In the case of a harmonic bath, Eq. (3.2) shows that decoherence involves two-time interactions, where at least one of the two time points must be a blip, i.e., the forward and backward paths must differ over at least one time slice. In addition, the decoherence quenching is a real exponential, and thus much more efficient, when the forward and backward paths differ at both time points (blip-blip interactions). Blip-sojourn interactions are indirectly associated with decoherence, which manifests itself as a medium-induced phase (the imaginary term in the exponent of Eq. (3.2)). It has been shown^10,111 that blip-blip interactions are associated with a classical decoherence mechanism that involves primarily induced phonon transitions, while blip-sojourn interactions are of a purely quantum mechanical origin and occur via spontaneous phonon emission. (Interestingly, enforcing zero-point energy via quantization of the initial density allows inclusion of a portion of the spontaneous phonon contribution via the classical decoherence process.¹⁰) These observations are exploited in our quantum-classical path integral (QCPI) formulation,^73–75 which shows much promise as a rigorous tool for simulating complex systems.

Naturally, if the forward and backward paths differ substantially, the added classical decoherence interactions in the influence functional can practically annihilate the given path pair by making its weight negligible. Thus, forward-backward path pairs with many blips make exponentially small contributions to the path sum. Under classical-like conditions and strong system-bath coupling, the dynamics will approach the fully incoherent limit where forward and backward paths can differ only infinitesimally. Since the damping effects of multiple blips are cumulative, the contribution of forward-backward path pairs generally decreases monotonically as the number of blips is increased. As a consequence, the number of blips provides an ideal convergence parameter.

To illustrate these ideas, consider again the TLS example of Sec. III and a propagation time equal to NΔt. For diagonal elements of the reduced density matrix (assuming a general initial condition), blips can occur at any of the N time points with k = 0, …, N − 1. There are N!/b!(N − b)! possible configurations of b blips, and

possible coordinate arrangements of these blip configurations. Thus, if the maximum number of blips is restricted to b_max, the total number of blip coordinate arrangements will be equal to

Each of these blip coordinate arrangements requires an iterative propagation of (N − b)-point sojourn segments. If b_max = N, i.e., if all blips are included, the sum in Eq. (4.2) simplifies to 3^N, which is the expected result for distributing three difference coordinate values (Δσ = 0, ±2) in N time slots. Table I shows the resulting number of terms in the blip sum for various values of N and b_max. As explained in Secs. III and IV, even in the full blip sum limit (which in the case of a TLS requires the evaluation of 3^N terms), the blip decomposition results in far fewer terms than a direct summation with respect to all forward-backward path coordinates (which amounts to summing 2^2N terms). Since the required CPU time is directly proportional to the number of terms, the exact blip decomposition reduces the computational effort to ∼(3 / 4)^N of that required by the conventional, direct forward-backward path sum. For example, the CPU time is reduced by a factor of ∼300 for N = 20, making exact summation of the path integral easily feasible. Perhaps more importantly, as seen from the representative numbers shown in Table I, the savings realized through the combined blip decomposition and blip truncation of the path sum are even more dramatic.

The ideas described in Secs. II–IV are illustrated using a model symmetric TLS described by the Hamiltonian of Eq. (3.5). The TLS is coupled bilinearly to a harmonic dissipative bath described by the commonly employed Ohmic spectral density,⁴ 

with ω_c = 0.25Ω and Kondo parameter ξ = 10. The ideas described in this paper are not restricted to a particular form of the spectral density; the Ohmic form is often chosen because it leads to rich dynamical behaviors. The TLS is initially in a localized state, |$\tilde \rho (0) = \left| 1 \right\rangle \left\langle 1 \right|$|$ρ̃(0)=11$⁠. Since the peak of the spectral density lies at ω = ω_c and the tunneling frequency of the TLS is 2Ω, the most strongly coupled bath modes are slower than the system by a factor of 8. This sluggish bath leads to very long influence functional memory, posing a challenge to iterative procedures. The calculations presented in this section proceed via direct (non-iterative) evaluation of the forward-backward path sum through blip decomposition and truncation.

Figure 2 shows converged results for the population of the initially localized state as a function of time at two temperatures which correspond to ℏΩβ = 0.5 and ℏΩβ = 0.125. These results were obtained with a path integral time step Δt = 0.375 and N = 32 time slices. The populations are shown for two different initial preparations of the bath. The first of these corresponds to a bath that is initially at thermal equilibrium in isolation. This is traditionally the choice in theoretical investigations of dissipative tunneling dynamics. In the second case the bath is at t = 0 in equilibrium with respect to the populated localized TLS state. This choice arises naturally in simulations of charge transfer, where the bath is initially in equilibrium with the electron donor. Figure 2 shows that the effects of initial bath preparation on the population dynamics are substantial at the lower temperature for this sluggish bath; however, these effects are very minor at the higher temperature. At long times the population approaches the thermodynamic value ρ₁₁ = 0.5.

In spite of the large bath reorganization energy (λ = 5ℏΩ), one observes a substantial degree of coherent oscillation in the TLS population at the lower temperature. This coherence is more pronounced in the case of an initially unshifted bath, whose role in quenching the TLS oscillations is severely reduced compared to a bath initially in equilibrium with the populated TLS state. Thermal excitation quenches the TLS coherence significantly but not completely, leading to a local flattening of the population curve.

Also shown in Figure 2 is the fully incoherent limit, which corresponds to identical forward-backward path pairs. The fully incoherent approximation captures correctly the early time dynamics, but soon blip corrections become sizable. The contributions of paths with 1–8 blips are shown in Figure 3 in the case of an initially isolated bath at thermal equilibrium at the two chosen temperatures. At short times, path pairs with small numbers of blips tend to make the largest contribution to the population.

In general, the contribution from even numbers of blips tends to be mostly positive, while path pairs with an odd number of blips tend to contribute negative corrections. This can be understood by considering the number of kinks in the TLS paths. The presence of a single blip in a path pair with identical endpoints implies a kink, i.e., a transition from one TLS site to the other and back. This introduces two off-diagonal propagator factors in the amplitude of the path pair. Since the off-diagonal propagator for a symmetric TLS is a purely imaginary number, a path pair with a single kink has a negative contribution. Consider now path pairs with two blips. For large N, these blips are non-adjacent in the majority of these path pairs, thus they correspond to two kinks. Since the negative contributions from each of the kinks multiply, the overall contribution will be positive. Similarly, one concludes that a path pair with three separated blips (thus three kinks) will make a negative contribution, etc. (However, a path pair with two adjacent and a separated blip has only two kinks, thus it will make a positive contribution.) As seen in Fig. 3, these trends are violated mostly at short times (when the number of blips is not very small compared to the number of path integral steps) because the number of paths with adjacent blips can then be significant, leading to different kink counts.

As the temperature is increased, the damping from blip-blip interactions in the influence functional becomes more effective, especially for the path pair combinations associated with negative corrections. This effect, which is seen clearly in Fig. 3(b), leads to quenching of the oscillatory behavior of the TLS population dynamics.

The magnitude of the total contribution from all path pairs with a given number of blips peaks at a time that increases with increasing blip number and then decays, as the given paths become outnumbered by paths with more blips. As the number of blips is increased, the largest contribution of the blip set decreases in magnitude. In practice, one needs to increase the number of blips until convergence is reached over the time interval of interest. For the TLS parameters examined in this section, the number of blips required for convergence is smaller than 10. For a fixed path integral time step, the convergence becomes faster with increasing temperature and/or dissipation strength.

In addition to the dramatic savings attained through the blip decomposition and truncation, further acceleration of the path sum is possible by eliminating path configurations with very small weights. These small weights arise from the large number of kinks associated with particular blip arrangements. The simplest way of implementing this acceleration is by choosing a threshold (which will serve as a convergence parameter) and terminating the blip insertion process in the generation of blip configurations once a particular blip configuration drops below this threshold.

The blip decomposition of the path integral leads to a dramatic acceleration of real-time path integral calculations on quantum dissipative systems. This decomposition takes advantage of the absence of dissipative influence functional interactions between sojourn points, such that the influence functional terms for those can be mapped on an incoherent propagator with an external phase. This mapping allows iterative integration of those segments, replacing the number of discrete Feynman paths by the number of possible blip arrangements. This replacement leads to an exponential reduction in the number of path integral terms without introducing any approximations.

Further, it was argued that the contribution of forward-backward path pairs decreases exponentially with increasing number of blips. Thus, the number of blips can be used as a convergence parameter. At moderate to high temperature and dissipation strength, where the evolution of the quantum system tends to be mostly monotonic, the fully incoherent limit (which corresponds to the zero-blip term) provides a good approximation to the dynamics, and inclusion of correction terms from the first few terms in the blip series leads to rapid convergence to the exact result. Finally, a very simple procedure for filtering out blip configurations that make negligible contributions allows additional savings.

The dramatic savings attainable through the blip decomposition will be important in the simulation of processes where the timescale of bath correlations is comparable to the propagation time of interest, thus obviating iterative decompositions of the path integral. Under conditions of long bath-induced memory with a much longer relaxation time, iterative decompositions offer the only viable approach. However, the blip decomposition can be exploited to reduce the number of path segments that span the bath-induced correlation time very substantially compared to the reduction achievable via filtering procedures.^37,38,43,44 Progress along these lines will be reported elsewhere.

This material is based upon work supported by the National Science Foundation under Award Nos. CHE 11-12422 and CHE 13-62826.