The QuAPI Method for System-Bath Dynamics

Our work in the 1990s focused on processes that can be described in terms of a low-dimensional system interacting with a bath of harmonic oscillators . The system-harmonic bath model offers a realistic picture of a variety of processes. For example, it provides an accurate description of crystalline solids since lattice vibrations are harmonic to a very good approximation. If the coupling between system and bath is bilinear, the dynamics is equivalent classically to that of a one-dimensional system obeying a generalized Langevin equation. Further, the system-bath model is often useful in reaction path Hamiltonian descriptions of polyatomic processes in the gas phase, where the bath consists in molecular vibrations orthogonal to the reaction coordinate.

We construct accurate propagators for a system coordinate coupled to a harmonic bath by choosing as the reference the one-dimensional Hamiltonian with the potential along the adiabatic path. The latter is defined as the path of minimum potential energy for fixed values of the system coordinate and is designed to offer an accurate description of high-frequency bath degrees of freedom. This partitioning of the system-bath Hamiltonian is similar in spirit to a small polaron transformation. Employing the adiabatic reference to split the short-time evolution operator and integrating out the harmonic bath degrees of freedom leads to a quasi-adiabatic propagator path integral (QuAPI) that involves one-dimensional system propagators which describe the exact dynamics along the adiabatic path along with an influence functional which incorporates non-adiabatic corrections. The latter enter in the form of Franck-Condon factors that arise from displaced harmonic oscillator modes. It can be shown that the error in this adiabatic splitting of the time evolution operator depends sensitively on the time step as well as the non-adiabaticity of the Hamiltonian and for this reason the quasi-adiabatic propagator allows fairly large time steps if the bath consists of high-frequency degrees of freedom. In addition, the quasi-adiabatic propagator is also accurate for slow bath degrees of freedom whose characteristic periods of motion exceed significantly the time step.

For short to intermediate time, the dimension of the QuAPI expression is not prohibitively large, allowing use of multidimensional integration methods. For condensed phase processes characterized by strong dissipative forces Monte Carlo methods can often be used successfully. In most other cases importance sampling fails due to the phase oscillation problem and one must resort to multidimensional quadrature. We have adopted system-specific discrete variable representations (DVR) for discretization of the path integral coordinates, as these lead to compact grids and allow a uniform description of continuous potentials as well as systems defined in terms of a finite number of levels. Nevertheless, there are several situations where conventional discrete variable representations do not lead to compact representations of time-dependent quantities. Typical examples include motion in unbound potentials as well as relaxation of highly excited states in the condensed phase. We have recently extended these schemes by developing time-dependent discrete variable representations which employ a moving basis. Such a time-dependent basis samples at any instant in time only part of the available energy spectrum and leads to significant savings in problems where the wavefunction (or density matrix) visits many time-independent states during the course of evolution.

While the above techniques enable evaluation of the QuAPI expression typically up to one or two periods of motion, propagation for long times must be addressed differently. The need for use of multidimensional integration methods arises from the nature of the influence functional, which is nonlocal in time.  Although the density matrix of all system and bath degrees of freedom obeys Markovian dynamics and the corresponding path integral is equivalent to the Schrödinger equation, the path integral expression for the reduced density matrix where the bath has been traced out involves memory effects and cannot be mapped onto an ordinary differential equation. The presence of memory prohibits iterative evaluation of the path integral by means of matrix-vector multiplication schemes routinely used for wavefunction propagation.

When the bath corresponds to a macroscopic environment it is effectively described by a continuous spectral density function. In that case the influence functional in the QuAPI expression can be expressed in the form given by Feynman and Vernon with the time-dependent force calculated along a modified contour and the time nonlocality is characterized by the decay time of the bath response function. The latter always decays within a finite time interval if the bath is described by a broad spectral density, implying that the nonlocality in the path integral spans a finite range. This is a consequence of destructive phase interference among a continuum of frequencies characteristic of condensed phase environments. Due to the finite range of nonlocal interactions, the dissipative influence functional resembles the partition function for an Ising model with non-nearest-neighbor interactions, which can be evaluated by transfer matrix techniques. Thus, by discarding negligible long-range influence functional connections, one can break up the path integral into multiple integrals of lower dimension and evaluate it iteratively. This is accomplished by defining a higher-dimensional object, a tensor of rank equal to the number of time steps necessary to span the memory length, which obeys Markovian dynamics.  This object can also be thought of as an array of path segments that span the bath correlation length.  This array (or tensor) is propagated in time through multiplication with a propagator matrix (or higher-rank tensor) which contains all relevant couplings. Projection after each iteration to eliminate the auxiliary variables yields the system reduced density matrix whose diagonal elements represent state populations and which can be used to extract observables.


By avoiding global summation over paths, the iterative quasi-adiabatic propagator path integral (i-QuAPI) scheme requires effort that scales linearly with the total number of propagation steps. However, the number of path segments that must be stored may be unrealistically large, in particular if the time nonlocality is long. This is the case with sluggish environments, such as solvents containing heavy atoms and/or low frequency skeletal vibrations of biological molecules, which give rise to correlation functions that can decay very slowly. It is thus necessary to restrict the propagated array to include only statistically significant path segments. These are selected either by Monte Carlo filtering, or by a deterministic procedure based on propagator values.  In incoherent regimes, a vast reduction of the number of necessary path segments is achieved through the blip decomposition of the forward-backward path sum.

Further, we have shown that even in situations where the bath is characterized by high-frequency modes (which necessitate a small time step) and slow oscillators (which give rise to very long memory), efficient propagation that takes into account the memory on all time scales is possible by renormalizing the propagator matrix to include successively larger-length influence functional correlations. At its crudest level, the renormalized propagator is an ordinary matrix, thus the computational demands of the method are very modest, similar to those for propagating the bare system.  By including memory on successively increasing time scales, the renormalized propagator is capable of capturing the correct nonexponential kinetics induced by sluggish solvents.

Related articles:


Back to Research