On iterative path integral calculations for a system interacting with a shifted dissipative bath

System-bath models, where the effects of a condensed phase or molecular environment on the dynamics of a system are described in terms of a bath of harmonic oscillators linearly coupled to the system,^1–4 continue to find widespread use. Apart from inviting analytical treatments, the harmonic bath model allows calculation of system properties using numerically exact quantum mechanical methods. In particular, the iterative quasi-adiabatic propagator path integral (i-QuAPI) methodology^5–26 offers an efficient approach for treating continuous or discrete dissipative systems, which may interact with external time-dependent fields,^27–30 over extremely long time lengths. Because of its iterative nature, the i-QuAPI algorithm circumvents the exponential proliferation of quantum paths with propagation time, as well as the Monte Carlo sign problem associated with the phase of the real-time path integral. The methodology has been employed for generating accurate benchmark results and has found numerous applications in chemical and condensed matter physics (for example, see Refs. 31–43); it has also been extended to fermionic environments.^44–46 

The i-QuAPI methodology employs a path integral^47,48 representation of the system’s reduced density matrix in terms of localized system states.⁷ Many generic studies (e.g., investigations of “spin-boson” dynamics^1,2,4) assume that the bath is initially isolated from the system of interest. However, in many chemical processes, the system is described in terms of such a localized state (the “donor” in the case of charge transfer), and the dissipative environment is initially in equilibrium with this localized state of the system. Such an initial condition requires a modification of the influence functional⁴⁹ employed in the i-QuAPI methodology.

In this paper, we discuss two equivalent procedures for generalizing the i-QuAPI algorithm to situations where the harmonic bath is initially in equilibrium with the localized state of the system. The first approach⁵⁰ involves evaluation of the influence functional with a shifted bath. This shift augments the system amplitude via a time-local phase, which is given in terms of the spectral density. The alternative, second approach consists of shifting the coordinate of the system to bring its initial state in equilibrium with the unshifted bath. This approach requires no modification of the influence functional and has already been used in earlier work by our group.³² We show that the two approaches are entirely equivalent in the full path integral expression.

However, we point out that the first approach must be implemented rather cautiously in iterative path integral calculations. Specifically, the phase supplied by the shifted bath must be subject to memory truncation consistent with that employed in the iterative decomposition of the influence functional, otherwise spurious memory effects lead to very slow convergence. The second approach does not require particular caution. At the same time, that treatment suggests the most convenient form for the memory-truncated phase of the shifted bath approach, in terms of QuAPI-discretized influence functional coefficients. Numerical calculations demonstrate the equivalence of the two procedures and the importance of consistent memory truncation.

In Sec. II, we derive the influence functional for a shifted bath using a simple classical procedure. We also consider the alternative approach of shifting the system coordinate while retaining the influence functional in terms of an unshifted bath and show that it leads to the same modification of the phase of the system. In Sec. III, we consider the use of the shifted bath influence functional in the iterative decomposition of the path integral. We show that failure to implement proper memory truncation of this phase leads to slow convergence and suggest a simple and accurate treatment in the context of the standard i-QuAPI formulation. In Sec. IV, we present some concluding remarks.

In the QuAPI formulation, the Hamiltonian describing the system is written in the general form

where $σ n$ represent the localized discrete states or a discrete variable representation (DVR)^51,52 obtained from continuous coordinates, i.e., the system position operator is⁷ 

It is assumed that the system interacts with a bath of harmonic oscillators 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¹ 

The reduced density matrix of the system, $ρ ̃$⁠, is given by the expression

where F is the Feynman-Vernon influence functional. It is often assumed that the bath is initially at equilibrium, i.e., its initial density matrix is e^−βH_b/Tr e^−βH_b with

With this initial condition, which is illustrated in Fig. 1(a), the influence functional is given by the expression

where

is the response function of the bath. The last exponential in Eq. (2.7) arises from the counter terms¹ included in the system-bath Hamiltonian. (This particular form, Eq. (2.3), is necessary in order for the equilibrated donor-acceptor energy difference to be equal to 2ε.) Eq. (2.7) can also be written in the form

where Δs = s⁺ − s⁻ and $s ̄ = 1 2 ( s + + s − )$⁠. With the QuAPI propagator factorization,⁵ where the system path coordinates are specified by

the influence functional takes the form

where the coefficients η_k′k″ have been given in earlier publications.^12,13

However, when one is modeling a particular chemical process, such as electron or proton transfer, the bath is initially in equilibrium with the populated system state (the “donor”). This can be achieved either by shifting the bath initial state and calculating the influence functional that arises from that change, or by shifting the Hamiltonian of the system, to bring the initially populated right state in equilibrium with the unshifted bath. These two arrangements are illustrated in Figures 1(b) and 1(c) for a two-level system (TLS).

There are several ways to proceed in order to calculate the influence functional for this shifted initial condition of the bath. For example, one may follow the original procedure of Feynman and Vernon.⁴⁹ It is instructive to follow a simpler procedure, motivated by our analysis⁵³ in the context of the quantum-classical path integral (QCPI) formulation,^53–55 where the influence functional from a complex polyatomic environment is obtained in terms of classical trajectories within a forward-backward semiclassical^56,57 or quasiclassical⁵⁸ approximation. The starting point for the latter is the expression of the influence functional in the form

where Φ_j is the forward-backward classical action along a bath trajectory that experiences the average of the forces exerted by the system along its forward and backward quantum paths, and W_j is the Wigner function⁵⁹ of the bath. Eq. (2.12) derives from the linearized initial value representation⁶⁰ or linearized path integral approximation.⁶¹ In the case of a harmonic bath experiencing the averaged system force, the forward-backward action is

The Wigner function is given by the expression

Use of these expressions in Eq. (2.12) leads to the exact Feynman-Vernon influence functional.

For clarity, we focus below on a system consisting of just two (M = 2) “right” and “left” states, denoted $R$ and $L$⁠, respectively. The resulting TLS Hamiltonian has the form

where −ħΩ is the coupling (tunneling) matrix element, the left-right energy difference is 2ε, and σ_x, σ_z are the usual Pauli spin matrices.

Setting the position operator as $s ˆ = σ z$⁠, i.e., with σ_1,2 = ± 1 in Eq. (2.2), one can bring the initial state of the bath in equilibrium with the donor by shifting the Hamiltonian that defines the initial condition of the bath, i.e., redefining the initial bath Hamiltonian as

The corresponding Wigner function can be obtained from Eq. (2.14) via the substitution

Note, however, that the trajectories are not affected by the shift of the equilibrium density. Changing integration variables to y_j,0 = x_j,0 − λ_j, it is easy to see that the influence functional for this shifted initial condition is given by the expression

Noting that when $s ˆ = σ z$⁠, the last term in Equations (2.7) vanishes because $σ 1 2 = σ 2 2 =1$⁠, Eq. (2.18) becomes

Finally, Eq. (2.18) can be rewritten in terms of the spectral density,

Eq. (2.20) is the result obtained in Ref. 50.

To arrive at the discretized form of Eq. (2.19), we set t = NΔt and evaluate the time integral along each time segment of constant system coordinate. Following the QuAPI propagator splitting, Eq. (2.10), the time integral in Eq. (2.20) becomes

Performing the time integrals analytically, we obtain the following expression for the discretized influence functional that corresponds to the shifted bath,

where the new coefficients are given by the expressions

(Note that these coefficients are not identical to those in Ref. 50, which did not employ the QuAPI discretization of the phase in Eq. (2.20).) Thus, to bring the initial state of the bath in equilibrium with the initial state of the system, one should use the modified influence functional given by Equations (2.22) and (2.23).

An attractive alternative³² that does not require any additional calculations is to shift the initial state of the system, to place it in equilibrium with the unshifted harmonic oscillator bath, as illustrated in Figure 1(c). This is achieved by shifting the coordinates of the system to σ₁ = 0, σ₂ = − 2, i.e., setting $s ˆ = σ ̃ z$⁠, where

Note that the phase from the counter terms in Eq. (2.7) no longer vanishes. Because the bath Hamiltonian and initial state are unchanged, the influence functional is still given by Eq. (2.9) by replacing the system position s by s − 1. This substitution leaves Δs unchanged but replaces $s ̄$ by $s ̄ −1$⁠. Thus, the influence functional becomes

Comparison to Eq. (2.9) shows that shifting the system coordinate introduces the factor

Performing the inner integral of the first term leads to the result

which agrees with the shift factor in Eq. (2.19).

Alternatively, we can apply the system coordinate shift to the discretized form, Eq. (2.11),

It follows that the shift-induced change in the influence functional is given by the expression

with

Eq. (2.29) has the form of the shift-induced factor in Eq. (2.22). It is straightforward to show that $γ k = γ ̃ k$⁠.

To summarize this approach, by shifting the coordinates of the system, one can bring the initial state of the system in equilibrium with an unshifted bath. This way one does not need to alter the expression for the influence functional.

The expressions derived in Sec. II may be used in a full path integral calculation with a bath initially in equilibrium with a localized state of the system, such as the donor in a charge transfer reaction. However, long time calculations require an iterative decomposition of the path integral.^9,11,12 The iterative algorithm exploits the decay of the bath response function with time, which implies the system coordinates s^±(t′) and s^±(t″) become uncoupled when t′ > > t″. In its discretized form, the coefficients η_k′k″ tend to zero as k′ − k″ ≡ Δk increases. Truncating the double sum in Eq. (2.11) at some Δk_max restricts the range of nonlocality of the influence functional, allowing iterative propagation via a matrix-vector procedure which employs path segments that span Δk_max time steps.

The second procedure discussed in Sec. II, which leaves the influence functional unchanged, allows application of the iterative algorithm without modification. On the other hand, the first approach⁵⁰ results in a modified influence functional, so one needs to reformulate the iterative decomposition. It is clear that the shift factor in Eq. (2.22) does not involve nonlocal interactions, and thus this factor may be included in the system propagator. Thus, it appears best to include the full phase (without performing memory truncation).⁵⁰ Using the expression for γ_k given in Eq. (2.23), adaptation of the i-QuAPI algorithm to account for a shifted bath seems straightforward.

However, this procedure leads to very slow convergence with increasing Δk_max. It appears that shifting the bath initial condition introduces long-lived memory, which can span a time length comparable to the time scale of the system’s decay to equilibrium. This behavior is illustrated in Figure 2, which shows the population of the right state for a TLS with ħΩ = 1, ε = 0 coupled to a dissipative bath characterized by the Ohmic spectral density,

with ω_c = Ω, at a temperature ħ ω_cβ = 1. The bath is initially in equilibrium with the right-localized TLS state. Full path integral results, obtained via the blip decomposition,^62,63 are shown for 25 path integral time steps of length ΩΔt = 0.125 and compared to results obtained via iterative evaluation. Iterative evaluation of the path integral with Δk_max = 10, using the shifted bath procedure with the shift coefficients given by Equations (2.23), is seen to fail after about 12 path integral steps. By contrast, iterative calculations with the same memory length remain accurate for 20 path integral steps if the shifted system procedure is employed. The reason for this discrepancy is that retention of full memory in the evaluation of the shift term is inconsistent when the memory has been truncated to Δk_max steps in the rest of the calculation. In fact, repeating the shifted system coordinate procedure, Equations (2.28)-(2.30), with the influence functional memory truncated at Δk_max time steps, one finds

i.e., the shift phase should be evaluated by retaining only the interval Δk_maxΔt preceding the current time.

Modeling condensed phase environments via a harmonic bath coupled to the system of interest is a common approach that allows numerically exact calculation of dynamical properties. In many such situations, the bath is initially in equilibrium with a localized state of the system, such as the state of the electron donor in the case of charge transfer. We have discussed two straightforward procedures for generalizing the QuAPI methodology to these shifted bath situations. The first approach involves an additional phase, which augments the action (or propagators) of the system, and which is given in terms of time integrals of the spectral density, or, equivalently, in terms of the coefficients that enter the standard influence functional. The second approach consists of shifting the coordinate of the system, to bring the donor state in equilibrium with the unshifted bath. This latter approach requires no additional phase factors.

Iterative decomposition of the path integral with the first approach requires consistent memory truncation in the evaluation of the shift phase, otherwise the spurious memory effects can lead to very slow convergence. The second approach (shifting the system coordinate) automatically employs consistent memory truncation and requires no modification of the i-QuAPI methodology.

While the i-QuAPI methodology with an analytic influence functional is applicable to systems in contact with harmonic baths, the ideas and subtleties pointed out in Sec. II and III have implications for the simulation of quantum mechanical processes in complex polyatomic environments using the QCPI formulation, where the influence functional is obtained from information extracted via classical molecular dynamics calculations.

This material is based upon work supported by the National Science Foundation under Award No. CHE 13-62826.