Quantum-classical path integral with a harmonic treatment of the back-reaction

Understanding the intricate dynamics of many processes in condensed phase and biological systems often hinges on one’s ability to perform all-atom simulations. Molecular dynamics (MD) tools are often the method of choice, because of the favorable, practically linear scaling of Newton’s equations with the number of degrees of freedom. However, it is often necessary to account for at least some quantum mechanical effects, such as nonadiabatic transitions (as in the case of electron transfer), zero-point energy (ZPE), and sometimes tunneling. Since a fully quantum mechanical treatment of the time-dependent Schrödinger equation for hundreds or thousands of atoms remains computationally prohibitive, there have been intense efforts to develop approximate simulation tools that can correctly capture the most critical aspects of the dynamics.

Over the last few years, our group has pursued a rigorous approach to mixed quantum-classical dynamics, i.e., situations where the most important quantum dynamical effects are confined to a small subsystem (typically an electron transfer pair or a proton coordinate), while the remaining atoms (those of the environment, such as a solvent or biological molecule) may be treated by classical trajectories. In the traditional Schrödinger description of quantum mechanics, specifying the contribution to a classical trajectory from a delocalized wavefunction leads to Ehrenfest’s model,¹ which involves a mean-field treatment that cannot correctly account for wavepacket splitting. A variety of approximate schemes have been developed for circumventing the most serious shortcomings of the Ehrenfest model (for example, see Refs. 2–6). Another possibility involves mapping discrete electronic states on continuous degrees of freedom,^7,8 circumventing the need for trajectory hops. Another approach that retains the discrete character of the quantum system proceeds by solving the quantum-classical Liouville equation via an expansion in the mass ratio.^9–11 

Our group has utilized the path integral formulation^12,13 of quantum mechanics to remove the incompatibility of the wavefunction-based formulation with the local description of classical trajectories. Expressing the time-dependent density matrix in the path integral form and taking the classical limit with respect to the degrees of freedom composing the system’s environment led to a rigorous quantum-classical path integral (QCPI) formulation.^14,15 The QCPI methodology treats the interaction between system and environment in full atomistic complexity and is free of any assumptions. This is possible because the path integral allows a local description of the system dynamics, which (in contrast to wavefunction treatments) leads to unambiguous determination of the quantum-classical force. The QCPI expression faithfully accounts for quantum interference and its destruction (i.e., the dephasing and decoherence effects from the environment) through quantum mechanical phases that interfere constructively or destructively, without the need for ad hoc decoherence factors or adjustable parameters. While the computational cost of QCPI certainly exceeds significantly that of classical MD, charge transfer simulations with thousands of solvent atoms¹⁶ have demonstrated that it is not prohibitive.

The present work aims at increasing the efficiency of QCPI calculations by reducing the number of required trajectories. To this end, we revisit the phase that contains the interplay between the quantum system and the classical environment. The interaction phase contains time-local as well as time-nonlocal contributions. This can be seen by separating each trajectory into a path-independent component (the “classical path” approximation) and the deviations from that trajectory caused by the force supplied by the quantum system (the “back-reaction”). The path-independent component can be added to the system Hamiltonian, giving rise to an effective propagator which accounts for the system-environment interaction along the particular reference trajectory.¹⁷ It has been shown¹⁸ that upon performing the phase space integration, the contributions from the path-independent trajectories amount to all decoherence effects arising from stimulated phonon absorption and emission, which constitute the classical decoherence. While these effects appear as time-nonlocal terms in the integrated form,¹⁹ the dynamics in the reference propagator is entirely local and thus requires minimal computational effort, and the automatic inclusion of the bulk of decoherence in the reference propagator leads to large path integral time steps. The remaining effects arise from the “kicks” by the system, which change the course of a classical trajectory. These effects, i.e., the back-reaction, lead to additional quantum decoherence via spontaneous phonon emission,¹⁸ which is responsible for detailed balance. This mechanism is a truly nonlocal quantum memory.

The majority of effort in the evaluation of the QCPI expression is related to these time-nonlocal, quantum memory effects. As the unique sequence of forces along each path gives rise to a different classical trajectory (from the same initial condition), the number of trajectories at first appears computationally prohibitive. The proliferation of classical trajectories (whose number grows exponentially with propagation time) is the quantum-classical manifestation of time nonlocality.²⁰ Fortunately, the decoherence induced by the environment effectively leads to the loss of memory of a trajectory’s distant past, reducing the necessary trajectories to those that span this (typically short) quantum memory and allowing a tensor decomposition similar to that in the iterative quasiadiabatic propagator path integral^21–26 (i-QuAPI) methodology. Furthermore, judicious choices of reference trajectories²⁷ can further shorten the quantum memory, typically to a few path integral time steps, making full QCPI calculations manageable, at least for two-state systems such as donor-acceptor pairs.

Interestingly, all the effects arising from the back-reaction may be obtained in terms of a single (per initial condition) trajectory in the case of a harmonic environment.²⁸ This is so because the state-hopping of harmonic oscillator trajectories amounts to a phase independent of the initial phase space location. While this feature is unique to quadratic potentials, we consider, in this paper, the possibility of exploiting it as an approximation to QCPI dynamics in cases of complex, anharmonic environments. There are several justifications for such a treatment.²⁸ First, as discussed above, the bulk of the environment-induced decoherence is fully accounted for through the system-independent component of classical trajectories, which evolve on the full, anharmonic potential, and quantum decoherence effects play a more subtle role. Second, quantum memory is associated with spontaneous emission, which is a “vacuum” effect, and potentials are nearly quadratic near the minimum. Third, complex, strongly anharmonic environments, such as liquids and biomolecules, often collectively exhibit Gaussian-like effects on the dynamics of a quantum system; thus, their dynamical behavior is not very different from that induced by a harmonic bath.²⁹ 

In this paper, we pursue these ideas by exploring the harmonic treatment of the back-reaction, arriving at a methodology that requires a single classical trajectory from each initial condition. In Sec. II, we give an overview of the QCPI methodology. In Sec. III, we analyze the system-environment interplay and describe the harmonic back-reaction (HBR) approximation to QCPI dynamics, along with two possible schemes for identifying a proper harmonic bath. The accuracy of the HBR approximation is assessed through model calculations in Sec. IV, and Sec. V presents some conclusions.

The QCPI formulation employs a discrete state representation of the quantum system. The Hamiltonian for the isolated M-state system is written as

Here, |σ_n⟩ represent discrete charge transfer states or system-specific discrete variable representations³⁰ (DVRs) of a continuous coordinate, i.e., the system position operator has the form

where σ_n are the DVR eigenvalues. The total Hamiltonian is given by

where $q^$ and $p^$ are the coordinates and momenta of the particles comprising the system’s environment. The properties of the system at the time t = NΔt can be obtained from the reduced density matrix,

where Δt is the path integral time step. The QCPI expression of the reduced density matrix has the form¹⁴ 

where P(q₀, p₀) is the initial phase space distribution function of the environment, which specifies the weights of the classical trajectories, and Q is the function that contains all the dynamical effects in the quantum-classical approximation, i.e., the sum of amplitudes along all QuAPI-discretized²¹ system paths, along with phases arising from the interaction of the quantum system with the classical trajectory of the system’s environment along each quantum path. The classical trajectories of the environment are subject to a Newtonian force obtained from the gradient of H_env, with the system coordinate set equal to the average of the instantaneous forward and backward system coordinates. The dynamical function Q is the path integral expression (in the potential-optimized³¹ DVR path integral discretization³²) of the time-dependent reference Hamiltonian¹⁷ 

which augments the system by the external time-varying potential V_env along a classical trajectory q^ref(t), p^ref(t) (obtained from a particular initial condition q₀, p₀) that evolves on a single state or a predefined sequence of states. This function is given by the expression

where $Ûref$ is the time evolution operator for the time-dependent Hamiltonian defined in Eq. (2.6). Each path integral variable $sk±$ specifies the system coordinate along a particular QuAPI-discretized forward or backward path and takes on all possible DVR values. The function

(with ΔV_env = H − H_ref) is the residual action that accumulates from the hops of the trajectory to states that are different from that of the reference Hamiltonian.¹⁷ It has been shown¹⁸ that upon integrating with respect to the trajectory initial conditions, the dynamics generated by the reference Hamiltonian, i.e., the single-step QCPI propagator, captures correctly the energy exchange between system and environment that corresponds to stimulated absorption and emission of phonons, which amounts to the classical mechanism of decoherence, even with the simplest choice of a reference trajectory, and that some of the remaining quantum decoherence process (which is associated with spontaneous phonon emission) can be recovered within the single-step propagator through improved choices of the reference trajectory.²⁷ The remaining spontaneous phonon emission events, which are crucial for satisfying detailed balance, are accounted for through the “back-reaction,” i.e., the dependence of classical trajectories on the forward-backward system paths resulting from time slicing. The inclusion of the system-environment phase in the reference Hamiltonian leads to convergence with larger path integral time steps. The reference propagator can be constructed by numerically solving the time-dependent Schrödinger equation, a first order differential equation, in the basis of system states,

Since the sequence of system states along each system path specifies a different classical trajectory, the path sum with respect to system coordinates $s0±,s1±,...,sN±$ gives rise to exponential proliferation of classical trajectories with time. By taking advantage of the memory-quenching effects of dissipative environments, it is possible to drop trajectory branches corresponding to the distant past, maintaining a constant number of trajectories. This allows an iterative decomposition of the QCPI expression, which leads to linear scaling with propagation time.¹⁵ We emphasize that the memory arising from the classical decoherence mechanism is not truncated at all, as the corresponding action is captured entirely through the reference propagator. Judicious choices of the reference trajectories²⁷ can account for some quantum memory as well, allowing convergence with a shorter memory length. Furthermore, path filtering schemes^33–35 may be implemented to drastically reduce the number of paths that must be included, and blip decompositions^36,37 should lead to an exponential acceleration of the QCPI sum.

Finally, we note that (unlike the sum with respect to system paths, which is evaluated by quadrature) the multidimensional phase space integral in the QCPI expression is evaluated by Monte Carlo³⁸ using P(q₀, p₀) as the sampling function. This phase space function may be given by the classical approximation (e.g., the Boltzmann factor for finite-temperature calculations) or by its quantum mechanical analog, the Wigner function.³⁹ The classical form would be the natural choice if the dynamics of the system’s environment is described by classical force fields, which are designed to account for quantum effects to some extent. However, if the trajectory forces are obtained through ab initio electronic structure calculations, initial conditions should be sampled from a quantized distribution. The Wigner function is available analytically for quadratic Hamiltonians but must be evaluated numerically for other cases. This is not an easy task because the multidimensional Fourier-type integral leads to extensive cancellation that appears like a Monte Carlo sign problem. Because of the utility of the Wigner function in quasiclassical trajectory calculations, a number of methods have focused on this task. These include local⁴⁰ or variationally optimized⁴¹ Gaussian wavepacket approaches and the thermal Gaussian approximation^42,43 (which employ frozen Gaussian dynamics⁴⁴ in imaginary time), and extensions of these methods can account for quantum corrections.⁴⁵ Another trajectory-based approach^46,47 exploits the classical adiabatic theorem. Path integral methods have also been developed for constructing the Wigner function for system-bath models⁴⁸ and also for Hamiltonians with unrestricted many-particle interactions as path integral methods.⁴⁹ 

For a quantum memory length equal to LΔt, there are M^2L forward-backward system paths and (in general) [M(M + 1)/2]^L classical trajectories from each initial conditions. In the special case of a harmonic environment, the residual action, Eq. (2.8), is independent of trajectory initial conditions and can be evaluated analytically, reducing the effort to a single trajectory per initial condition.²⁸ Furthermore, if the reference Hamiltonian corresponds to a fixed system state, this action is that of a forced harmonic oscillator and amounts to the imaginary part of the exponent in the Feynman-Vernon influence functional.⁵⁰ Thus, the contribution of the back-reaction is easily obtained by analytical expressions similar to those in the QuAPI-discretized influence functional coefficients,²³ with minor differences arising from the QCPI discretization.

The ability to capture the system-bath interaction along [M(M + 1)/2]^L trajectories from a single one is a special feature of the harmonic potential, where the potential at different values of the system coordinate (i.e., different system states) corresponds to displaced parabolas of the same frequency. The resulting computational simplification of the QCPI algorithm for a system-bath Hamiltonian is dramatic.²⁸ Since the residual action (i.e., the imaginary part of the influence functional) is a “vacuum” effect related to spontaneous emission, one wonders if a quadratic treatment of this action can adequately capture the back-reaction effects also in the case of anharmonic environments, eliminating the need for the explicit numerical integration of a large number of trajectories. Thus, we consider replacing Eq. (2.8) by its harmonic back-reaction (HBR) approximation. For this purpose, we define the HBR Hamiltonian,

To proceed, one must first decide how to choose the harmonic bath frequencies and coupling coefficients. We consider two ways: (i) using the quadratic expansion of the potential, i.e., the normal modes of the environment at the potential minimum and (ii) using the effective harmonic bath that corresponds to Gaussian response (GR), as suggested by the central limit theorem and obtained through a cumulant expansion of the influence functional.²⁹ This approach requires the calculation of the appropriate time correlation function of the environment, from which harmonic bath parameters (e.g., a spectral density) may be obtained. With the fixed-state reference, the HBR effects are equivalent to those in the imaginary part of the influence functional and may be computed either by the standard procedure²⁴ (in the QCPI discretization¹⁵) or directly from the time correlation function.⁵¹ The details and the extension to the dynamically consistent state-hopping reference will be described in another paper.⁵² 

In Sec. IV, we assess the accuracy of these treatments.

We use a model two-level system (TLS) described by the Hamiltonian

where ℏΩ₀ is the tunneling matrix element and ε is the asymmetry parameter. We choose the two-level system (TLS) tunneling frequency Ω₀ = 1.

We test the HBR approximation on an environment composed of independent quartic oscillators,

To avoid confusion, we emphasize that the parameters Ω_j and C_j specify the model Hamiltonian and are not necessarily equal to the HBR parameters ω_j and c_j. [The HBR mass parameters are arbitrary, and thus, we do not use different symbols for the masses in Eqs. (3.1) and (4.2).] The parameters b_j quantify the potential anharmonicity. We choose values such that the ZPE of each quartic oscillator, as obtained by first-order perturbation theory, is 10% larger than the harmonic ZPE, i.e.,

where Ψ_0,j is the harmonic oscillator ground state wavefunction. This leads to the identification

We choose m_j = 1 for all oscillators. Figure 1 compares this quartic potential to its harmonic counterpart, along with the two lowest energy levels. It is seen that the anharmonicity is quite large.

The harmonic frequencies and coupling coefficients of the model Hamiltonian are specified by the relations

which are motivated by the logarithmic discretization^29,53 of a harmonic bath with cutoff frequency Ω_c and dimensionless dissipation parameter ζ.⁵⁴ The number of oscillators is chosen to be n = 20 − 30. With these choices, the dissipative effects are sufficient to prevent recurrences over the times of interest, but the number of oscillators is not yet sufficiently large for Gaussian response.

The bath initial conditions are taken to be the Boltzmann distribution,

To have the oscillators in equilibrium with the initial localized state of the system,⁵⁵ we choose the two eigenvalues of the TLS as σ₁ = 0 and σ₂ = 2.

Choice (i) of the harmonic bath employed in the HBR treatment consists of setting ω_j = Ω_j. We note that this identification is obvious for the given model, but would generally require knowledge of the Hessian matrix and a normal mode analysis, and may not be meaningful at all if the environment does not possess a well-defined equilibrium geometry, as in the case of a liquid where a normal mode analysis would give rise to imaginary frequencies.

For choice (ii), we need to identify the effective harmonic bath that matches the force-force autocorrelation function of the system’s environment,

where the force is given by the derivative of the potential with respect to the system coordinate,

The discretized effective bath modes are given by⁵⁵ 

The anharmonic oscillator trajectories were obtained by integrating the classical equations of motion with the Verlet algorithm.⁵⁶ 

Figures 2 and 3 show results for a symmetric TLS (ε = 0) coupled to an anharmonic bath of n = 25 oscillators at various temperatures and coupling strengths. In each case, we present the population of the initially localized TLS state as obtained from the full QCPI expression and also from QCPI calculations with the HBR treatment, with the two harmonic bath choices discussed in Sec. III. In addition, we show results for the equivalent harmonic bath obtained through the mapping Eqs. (4.7)–(4.9).

With the HBR bath frequencies chosen from the harmonic frequencies at the potential minimum, the HBR-QCPI results show considerable deviations from the exact QCPI results. However, the HBR-QCPI treatment with the GR effective harmonic bath produces results in excellent agreement with the full QCPI results across a wide range of temperatures. The full GR harmonic bath results, which include no explicit anharmonicity in the bath potential, exhibit small deviations from the full QCPI results, implying that the collective response of the environment has not reached the full Gaussian response limit with these parameters.

Similar trends are observed in Fig. 4, which shows results on an asymmetric TLS. Again, the HBR-QCPI results with the GR harmonic bath parameters are in excellent agreement with the full QCPI values.

Many condensed phase processes are described reasonably well by system-bath models. Such a behavior is generally referred to as linear or Gaussian response and is a consequence of the central limit theorem. Electron transfer processes often fall in this category, as the long-range nature of the Coulomb potential implies the involvement of a large number of atoms in the solvent reorganization.⁵⁷ On the other hand, many processes exhibit significant deviations from the GR model. It is thus important to allow classical trajectories in a QCPI simulation to experience the nonlinear forces of the system’s environment, without invoking a harmonic bath assumption.

The HBR treatment described in this paper leads to a dramatic acceleration of QCPI calculations, reducing the computational cost to a single classical trajectory per initial condition, thus achieving MD efficiency (with the relatively small additional effort associated with evaluating the path sum within the quantum memory length). As discussed in Sec. III, the gain scales exponentially with the length of quantum memory induced by the environment, which is more pronounced under low-temperature conditions in solvents with a large reorganization energy. The HBR calculations presented in Sec. IV, which employ a model bath of strongly anharmonic oscillators, led to results practically identical to the full QCPI dynamics. Furthermore, application of the HBR treatment to an electron transfer reaction in solution quantitatively reproduced the full QCPI results.⁵² We thus envision that the HBR-QCPI methodology will find wide application as an accurate, yet efficient tool for molecular simulation of condensed phase and biological processes.

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