Real-time path integral calculations for the propagation of a system in contact with a harmonic dissipative environment often employ the iterative quasi-adiabatic propagator path integral (i-QuAPI) methodology. We compare two simple ways of applying this methodology to a bath initially in equilibrium with the localized state of the system (e.g., the donor in the case of charge transfer). The first way involves modifying the phase of the system via a time-local phase given in terms of integrals of the spectral density or in terms of the coefficients entering the QuAPI-discretized influence functional. In the iterative decomposition of the path integral, this approach requires consistent memory truncation to avoid extremely slow convergence. The second, alternative approach involves shifting the coordinate of the system, to bring the donor state in equilibrium with the bath, and requires no further modification of the i-QuAPI algorithm.
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) methodology5–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 integral47,48 representation of the system’s reduced density matrix in terms of localized system states.7 Many generic studies (e.g., investigations of “spin-boson” dynamics1,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 functional49 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 approach50 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.32 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.
II. INFLUENCE FUNCTIONAL FOR BATH INITIALLY IN EQUILIBRIUM WITH DONOR
In the QuAPI formulation, the Hamiltonian describing the system is written in the general form
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 function1
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−βHb/Tr e−βHb with
With this initial condition, which is illustrated in Fig. 1(a), the influence functional is given by the expression
is the response function of the bath. The last exponential in Eq. (2.7) arises from the counter terms1 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 . With the QuAPI propagator factorization,5 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.49 It is instructive to follow a simpler procedure, motivated by our analysis53 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 semiclassical56,57 or quasiclassical58 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 Wj is the Wigner function59 of the bath. Eq. (2.12) derives from the linearized initial value representation60 or linearized path integral approximation.61 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 and , 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.
A. Influence functional for shifted bath
Setting the position operator as , 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 yj,0 = xj,0 − λj, it is easy to see that the influence functional for this shifted initial condition is given by the expression
Finally, Eq. (2.18) can be rewritten in terms of the spectral density,
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).
B. Shifting the system Hamiltonian
An attractive alternative32 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 σ1 = 0, σ2 = − 2, i.e., setting , 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 by . 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
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.
III. MEMORY TRUNCATION FOR ITERATIVE DECOMPOSITION
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 Δkmax restricts the range of nonlocality of the influence functional, allowing iterative propagation via a matrix-vector procedure which employs path segments that span Δkmax 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 approach50 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).50 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 Δkmax. 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 Δkmax = 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 Δkmax 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 Δkmax time steps, one finds
i.e., the shift phase should be evaluated by retaining only the interval ΔkmaxΔt preceding the current time.
IV. CONCLUDING REMARKS
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.