Correlation Functions in Open Quantum-Classical Systems

Abstract

Quantum time correlation functions are often the principal objects of interest in experimental investigations of the dynamics of quantum systems. For instance, transport properties, such as diffusion and reaction rate coefficients, can be obtained by integrating these functions. The evaluation of such correlation functions entails sampling from quantum equilibrium density operators and quantum time evolution of operators. For condensed phase and complex systems, where quantum dynamics is difficult to carry out, approximations must often be made to compute these functions. We present a general scheme for the computation of correlation functions, which preserves the full quantum equilibrium structure of the system and approximates the time evolution with quantum-classical Liouville dynamics. Several aspects of the scheme are discussed, including a practical and general approach to sample the quantum equilibrium density, the properties of the quantum-classical Liouville equation in the context of correlation function computations, simulation schemes for the approximate dynamics and their interpretation and connections to other approximate quantum dynamical methods.

1. Introduction

The dynamical properties of condensed-phase or complex systems are often investigated experimentally by applying external fields to weakly perturb a system and observe its relaxation back to the thermal equilibrium state. In such experiments, measurable quantities can be related to equilibrium time correlation functions via linear response theory [1,2]:

$$C_{\mathit{AB}}(t)=\frac{1}{Z_Q}\text{Tr}\hspace{0.17em}\hspace{0.17em}\left[e^{-\beta \hat{H}}\hat{A}(0)\hat{B}(t)\right]=\frac{1}{Z_Q}\text{Tr}\hspace{0.17em}\hspace{0.17em}\left[e^{-\beta \hat{H}}\hat{A}{e}^{\frac{i}{\hslash }\hat{H}t}\hat{B}{e}^{-\frac{i}{\hslash }\hat{H}t}\right]$$

(1)

where  and B̂ are operators corresponding to some specific dynamical variables under investigation, Ĥ is the unperturbed Hamiltonian and Z_Q is the quantum canonical partition function associated with Ĥ. Many experiments employing spectroscopic methods directly probe such time correlation functions.

Exact numerical evaluation of Equation (1) for real condensed phase quantum systems is prohibitive, since the computational cost scales exponentially with respect to the number of degrees of freedom (DOF). Various approaches have been developed to address this challenging problem. A common approach shared by many methods is to partition the entire system into a subsystem (whose dynamical properties are of interest) and an environment (or bath) in which the subsystem resides. Other recently developed schemes for computing quantum correlation functions do not rely on such a partition and instead utilize approximations to treat the quantum evolution of the entire system in conjunction with quantum equilibrium sampling [3–5]. In this paper, we focus on schemes based on the system-bath partition, and using this partition, the Hamiltonian reads: Ĥ= Ĥ_b+ĥ_s+V̂_c(R̂); where ${\hat{H}}_b=\frac{{\hat{P}}^2}{2M}+{\hat{V}}_b(\hat{R})$ and ĥ_s represent the pure bath and subsystem Hamiltonians, respectively. The last term in Ĥ is a coupling potential that depends on the spatial coordinates of the bath wave functions. We shall always take the bath part of the Hamiltonian in the coordinate representation; however, we can represent ${\hat{h}}_s=\frac{{\hat{p}}^2}{2m}+{\hat{V}}_s(\hat{r})$ in some quantum basis: ${\hat{h}}_s={\sum }_{ij}|i\rangle \langle i\left|{\hat{h}}_s\right|j\rangle |j\rangle$.

Several methods based on various master equations [6–10] and path integral influence functional methods [11,12] provide approximate schemes, often in the weak coupling limit, to systematically project out the environmental DOF and yield a subsystem dynamics that incorporates dissipation and decoherence, due to coupling to the environment. However, for many applications, such as proton and electron transfer in condensed phases, it is desirable to explicitly simulate, even approximately, the bath dynamics, since specific local bath DOF may be crucial for a description of the dynamics of the quantum subsystem. For this purpose, several semiclassical [13–15] and mixed quantum-classical [16,17] (MQC) methods, which either treat the entire dynamics semiclassically or simulate the dynamics of the bath and subsystem with different levels of rigor (e.g., classical versus quantum mechanical), have been formulated. Many semiclassical and mixed quantum-classical approaches, adopting powerful classical simulation techniques, evaluate Equation (1) by combined Monte Carlo-molecular dynamics (MC-MD) techniques.

In this paper, we formulate MC-MD schemes to evaluate Equation (1) within the framework of the quantum-classical Liouville equation (QCLE) [18]. The QCLE employs a partial Wigner representation of the environmental (bath) DOF and may be derived from full quantum dynamics by truncating the quantum evolution operator to the first order in a small parameter related to the ratio of the characteristic masses of quantum and bath DOF [18]. In particular, we suppose that the quantum subsystem has a finite-dimension Hilbert space. Under this assumption, Equation (1) is cast in the following form [19,20]:

$$C_{\mathit{AB}}(t)=\frac{1}{Z_Q}\sum _{n_1,n_2}\int dX\left[{\left(e^{-\beta \hat{H}}A\right)}_W^{n_1{n}_2}\left(X\right)B_W^{n_2{n}_1}\left(X,t\right)\right]$$

(2)

where the n_j indices label the basis states (in some chosen quantum basis), X = (R, P) represents the Wigner-transformed phase space point for the bath, N_B is the number of bath DOF and the subscript, W, on an operator indicates a partial Wigner transform on the bath DOF; e.g., an operator is partially Wigner transformed as ${\hat{B}}_W(X)=\int dZ\langle R-\frac{Z}{2}\left|\hat{B}\right|R+\frac{Z}{2}\rangle e^{\frac{i}{\hslash }P\cdot Z}$.

Two main tasks are involved in evaluating Equation (2) with an MC-MD algorithm. First, one needs to sample initial conditions (for an ensemble of trajectories) from the partially Wigner-transformed quantum density, ${\left({\hat{\rho }}_{eq}\hat{A}\right)}_W(X)$ with ρ_eq = e^−βĤ/Z_Q. There exist numerical algorithms to accomplish such a task [21,22]. Second, one needs to propagate the initial points in the phase space. These time-evolved trajectories may then be used to construct the matrix elements, $B_W^{nm}(X,t)$, needed to compute the correlation function. Various simulation methods, whose structure depends on the basis chosen to represent the quantum degrees of freedom in the QCLE, have been devised to simulate the mixed quantum-classical dynamics [23–31]. Simulation methods that utilize an adiabatic basis can be cast into the form of surface-hopping dynamics, but in a way that includes coherent evolution segments that account for creation and destruction of coherence in a proper manner. More recently, as in some semiclassical approaches [32], the mapping basis [33] was used to describe the quantum degrees of freedom in the QCLE in a continuous classical-like manner, leading to a trajectory description in the full system phase space [30,31,34–36].

The goals and outline of the paper are as follows: We first consider how the two ingredients, quantum equilibrium sampling and evolution of quantum operators, which are needed to compute quantum correlation functions, may be carried out. In Section 2, we describe a path-integral scheme to perform MC sampling from the partially Wigner transformed quantum density. In the Appendix, we also discuss a simplified, but approximate sampling scheme that is useful in the high-temperature limit. Another aim of this paper is to demonstrate how a recently-developed simulation method for the QCLE, the forward-backward trajectory solution (FBTS), can be used to efficiently obtain quantum correlation functions. To place these results in proper context, in Section 3, we sketch the important features and properties of the QCLE and discuss both the adiabatic Trotter-based surface-hopping (TBSH) algorithm and the FBTS, which is formulated in the mapping basis. In this section, we also present the explicit form of the N-level generalization of the TBSH algorithm. Comparisons of the trajectories that underlie these algorithms allow us to investigate how completely different ensembles of trajectories can be used to simulate the same observable correlation function. The implementation and utility of the simulation algorithms are illustrated on the dynamics in a two-level system coupled to a quartic oscillator embedded in a bath of independent harmonic oscillators, described in Section 4. Finally, in Section 5, we comment on the advantages, challenges and potential problems in adopting an approximate mixed quantum-classical dynamics for the computation of quantum time-correlation functions.

2. Sampling from the Partially Wigner-Transformed Density

In general, analytical expressions for the Wigner transform of the density matrix cannot be determined easily. In this section, we present a path-integral-based scheme to perform MC sampling from the Wigner-transformed density, ${\left({\widehat{\rho }}_{\mathit{eq}}\hat{A}\right)}_W^{n_1{n}_2}(X)$, in Equation (2).

First, we recall the definition of partial Wigner transform:

$${\left({\widehat{\rho }}_{\mathit{eq}}\hat{A}\right)}_W^{n_1{n}_2}(X)=\frac{1}{Z_Q}\int dZ\langle n_1,R-\frac{Z}{2}\left|e^{-\beta \hat{H}}\hat{A}\right|n_2,R+\frac{Z}{2}\rangle e^{\frac{i}{\hslash }P\cdot Z}$$

(3)

where R represents the vector of bath coordinates, n denotes a basis state for the subsystem and $\hat{H}=\frac{{\hat{P}}^2}{2M}+{\hat{V}}_b(\hat{R})+\hat{h}(\hat{R})$ with ĥ(R̂) ≡ ĥ_s + V̂_c(R̂). One way to compute the integral on the right side of Equation (3) is to first factorize $e^{-\beta \hat{H}}=\prod e^{-{\beta }_L\hat{H}}$ into L − 1 pieces with β_L = β/(L − 1). Following the standard procedures for path integral calculations, we then insert resolutions of the identity, $ℐ=\int d{R}_i{\sum }_{m_i}|m_i,R_i\rangle \langle m_i,R_i|$, between every pair of factorized operators and apply the approximation, $e^{-{\beta }_L\hat{H}}\approx e^{-{\beta }_L\frac{P^2}{2M}}{e}^{-{\beta }_L\left({\hat{V}}_b(\hat{R})+\hat{h}(\hat{R})\right)}$. The integrand on the right side of Equation (3) can then be written as follows:

$$\begin{array}{l}\langle n_1,R-\frac{Z}{2}\left|e^{-\beta \hat{H}}\hat{A}\right|n_2,R+\frac{Z}{2}\rangle \\ \hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}=\int \prod _{i=1}^{L-1}d{R}_i\sum _{\left\{m_i\right\}}\langle n_1,R-\frac{Z}{2}\left|e^{-{\beta }_L\hat{H}}\right|m_1,R_1\rangle \langle m_1,R_1\left|e^{-{\beta }_L\hat{H}}\right|m_2,R_2\rangle \dots \\ \hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\times \langle m_{L-1},R_{L-1}\left|\hat{A}\right|n_2,R+\frac{Z}{2}\rangle ,\\ \hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}={\left(\frac{M}{2\pi {\beta }_L{\hslash }^2}\right)}^{\frac{N_B(L-1)}{2}}\int \prod _{i=1}^{L-1}d{R}_i\sum _{\left\{m_i\right\}}\left\{\prod _{i=1}^{L-2}{ℳ}_{i,i+1}(R_i)e^{-{\beta }_L{V}_b(R_i)}{e}^{-\frac{M}{2{\beta }_L{\hslash }^2}{\left(R_i-R_{i+1}\right)}^2}\right\}\\ \hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\times \langle n_1,R-\frac{Z}{2}\left|e^{-{\beta }_L\hat{H}}\right|m_1,R_1\rangle \langle m_{L-1},R_{L-1}\left|\hat{A})\right|n_2,R+\frac{Z}{2}\rangle \end{array}$$

(4)

where:

$${ℳ}_{i,j}=\langle m_i\left|e^{-{\beta }_L\hat{h}(R_i)}\right|m_j\rangle =\{\begin{array}{ll}{e}^{-{\beta }_L}{h}^{\mathit{ij}}(R_i),\hfill & i=j\hfill \\ -{\beta }_L{h}^{\mathit{ij}}(R_i)e^{-{\beta }_L{h}^{\mathit{ij}}(R)},\hfill & i\ne j\hfill \end{array}$$

(5)

which is correct to order $𝒪\left({\beta }_L^2\right)$. Substituting Equation (4) into Equation (3), the new integrand of the Wigner transform becomes $\widehat{𝒜}=\left(e^{-{\beta }_L\hat{H}}|m_1,R_1\rangle \langle m_{L-1},R_{L-1}|\hat{A}\right)$, as shown in the last line of Equation (4). An analytical approximation for the Wigner transform of 𝒜̂ can be obtained easily in most cases when  is a pure observable subsystem or if it depends on just one of the conjugate variables: R or P. Since β_L ≪ 1, it is possible to replace the term, e^−β_LĤ, inside 𝒜̂ with its high-temperature approximation (discussed in the Appendix). Letting 𝒜̂_W (X) be the partial Wigner transform of 𝒜̂, Equation (3) reads:

$${\left({\widehat{\rho }}_{\mathit{eq}}\hat{A}\right)}_W^{n_1{n}_2}(X)=\frac{{𝒢}^{N_B(L-1)/2}}{Z_Q}\int \prod _{i=1}^{L-1}d{R}_i\sum _{\left\{m_i\right\}}\left\{\prod _{i=1}^{L-2}{ℳ}_{i,i+1}(R_i)e^{-{\beta }_L{V}_b(R_i)}{e}^{-\pi 𝒢{\left(R_i-R_{i+1}\right)}^2}\right\}{𝒜}_W^{n_1{n}_2}(X)$$

(6)

where $𝒢=\left(\frac{M}{2\pi {\beta }_L{\hslash }^2}\right)$. Substituting Equation (6) into Equation (2), the time correlation function becomes:

$$\begin{array}{ll}{C}_{\mathit{AB}}(t)=\hfill & \frac{{𝒢}^{N_B(L-1)/2}}{Z_Q}\sum _{n_1,n_2}\sum _{\left\{m_i\right\}}\int \prod _{i=1}^{L-1}d{R}_i\left\{\prod _{i=1}^{L-2}{ℳ}_{i,i+1}(R_i)e^{-{\beta }_L{V}_b(R_i)}{e}^{-\pi 𝒢{\left(R_i-R_{i+1}\right)}^2}\right\}\hfill \\ \hfill & \times \int dX{\left(\widehat{𝒜}\right)}_W^{n_1{n}_2}(X)B_W^{n_2{n}_1}(X,t)\hfill \end{array}$$

(7)

Following [37], we remark that the initial phase space coordinate X = (R, P) and auxiliary variables, {R_i}, can be sampled from probability densities constructed from 𝒜̂_W(X) and |$ℳ$_i,i+1(R_i)|e^{−β_LV_b(R_i)}e^{−π$𝒢$(R_i−R_i+1)2}, respectively.

3. Quantum-Classical Liouville Equation

In this section, we discuss how one can simulate the time-evolved matrix elements, $B_W^{n_2{n}_1}(X,t)$, in Equation (2) using the QCLE:

$$\begin{array}{ll}\frac{\partial {\hat{B}}_W(X,t)}{\partial t}\hfill & =\frac{i}{\hslash }[{\hat{H}}_W,{\hat{B}}_W]-\frac{1}{2}\left(\{{\hat{H}}_W,{\hat{B}}_W\}-\{{\hat{B}}_W,{\hat{H}}_W\}\right)\hfill \\ \hfill & =i\widehat{ℒ}{\hat{B}}_W(X,t)=\frac{i}{\hslash }\left({\overrightarrow{\mathscr{H}}}_{\Lambda }{\hat{B}}_W-{\hat{B}}_W{\overleftarrow{\mathscr{H}}}_{\Lambda }\right)\hfill \end{array}$$

(8)

where $\Lambda ={\overleftarrow{\nabla }}_P{\overrightarrow{\nabla }}_R-{\overleftarrow{\nabla }}_R{\overrightarrow{\nabla }}_P$. The arrow on top of a differential operator indicates the direction in which it acts. In the first line, the square bracket and the curly brackets denote the quantum commutator and classical Poisson brackets, respectively. The two kinds of Lie bracket act together as the generator of the mixed quantum-classical dynamics. Due to the fact that Ĥ_W(X) and B̂_W(X, t) are quantum operators with respect to the subsystem DOF, two differently ordered Poisson brackets are needed to properly account for the mixed dynamics. However, in general, the dynamics described by the QCLE does not have a Lie algebraic structure, a feature that is common to mixed quantum-classical approaches [38]. In the second line, we introduce the abstract, quantum-classical Liouville (QCL) superoperator, $ℒ$̂. Finally, the third equality is another equivalent representation of QCLE in terms of the forward and backward mixed quantum-classical Hamiltonians:

$${\overrightarrow{\mathscr{H}}}_{\Lambda }={\hat{H}}_W\left(1+\frac{\hslash \Lambda }{2i}\right),\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}{\overleftarrow{\mathscr{H}}}_{\Lambda }=\left(1+\frac{\hslash \Lambda }{2i}\right){\hat{H}}_W$$

(9)

The QCLE has many desirable features, such as the conservation of energy, momentum and phase space volumes. Furthermore, the QCLE is equivalent to full quantum dynamics for arbitrary quantum subsystems, which are bilinearly coupled to a harmonic bath. For instance, commonly used spin boson models are of this type. In this circumstance, the combination of quantum and classical brackets in the QCLE does have a Lie algebraic structure. For the more general bath and coupling potentials, the QCLE provides an approximate description of the quantum dynamics. In this case, comparisons of simulations of QCL dynamics with exact quantum results have indicated that it is quantitatively accurate for a wide range of systems [36,39–48]

The QCLE equation can be simulated using ensembles of trajectories, which, in combination with the quantum initial condition sampling discussed above, provides a way to compute quantum correlation functions. As we shall see, the nature of the trajectories that enter in the simulations depends on the algorithm and should not be ascribed physical significance. It is only the observable, in this case, the correlation function, that has physical meaning and is independent of the manner in which it is simulated, provided the simulation algorithm is capable of exactly solving the QCLE, which is not always the case. One of the goals of this paper is to illustrate how a recently-developed FBTS [31] can be used to easily compute quantum correlation functions. For this purpose, it is interesting to contrast the solution using this scheme, and the trajectory description that underlies it, with the previously-developed and frequently-used TBSH algorithm [26]. Taking the adiabatic representation of the QCL superoperator is the key step in implementing the TBSH algorithm. The last representation of QCLE in Equation (8) resembles the quantum Liouville equation and forms the starting point of the FBTS.

3.1. Adiabatic Trotter-Based Surface Hopping

In order to discuss the nature of the trajectory description involved in the TBSH algorithm, we briefly describe how it is implemented and, in particular, present the explicit generalization to an N-level quantum subsystem, which was only outlined in [26]. We first consider the adiabatic representation of the QCLE, since the TBSH algorithm is cast in this basis. The adiabatic basis is defined by ĥ_W(R) | α; R〉 = E_α(R) | α; R〉, where ĥ_W(R) = Ĥ_W(R) − P²/2M is taken to be the adiabatic Hamiltonian for a static configuration of R in this section. In the adiabatic basis, the QCLE reads:

$$\frac{\partial B_W^{\alpha {\alpha }^{\prime }}}{\partial t}=i{ℒ}_{\alpha {\alpha }^{\prime },\beta {\beta }^{\prime }}{B}_W^{\beta {\beta }^{\prime }}(X,t)$$

(10)

where the matrix elements of the QCL superoperator are given by:

$$i{ℒ}_{\alpha {\alpha }^{\prime },\beta {\beta }^{\prime }}=\left(i{\omega }_{\alpha {\alpha }^{\prime }}+i{L}_{\alpha {\alpha }^{\prime }}\right){\delta }_{\alpha \beta }{\delta }_{{\alpha }^{\prime }{\beta }^{\prime }}-{𝒥}_{\alpha {\alpha }^{\prime },\beta {\beta }^{\prime }}=i{ℒ}_{\alpha {\alpha }^{\prime }}^0{\delta }_{\alpha \beta }{\delta }_{{\alpha }^{\prime }{\beta }^{\prime }}-{𝒥}_{\alpha {\alpha }^{\prime },\beta {\beta }^{\prime }}$$

(11)

with ω_αα′ = (E_α − E_α′)/ℏ. (The Einstein summation convention will be used throughout the following sections, although sometimes, sums will be explicitly written if there is the possibility of confusion.) The Liouville operator, i$ℒ$, may be separated into two contributions: The classical propagator is defined as:

$$i{L}_{\alpha {\alpha }^{\prime }}=\frac{P}{M}\cdot \frac{\partial }{\partial P}+\frac{1}{2}\left(F_{\alpha }+F_{\alpha \prime }\right)\cdot \frac{\partial }{\partial R}$$

(12)

where $F^{\alpha }=\langle \alpha ;R\left|\frac{\partial {\hat{h}}_W(R)}{\partial R}\right|\alpha ;R\rangle$ is the Hellmann-Feynman force. The superoperator, $𝒥$_{αα′,ββ′} is responsible for nonadiabatic transitions and associated momentum changes in the bath. For an N-level system, there exist N(N − 1)/2 unique transitions. In the following, we define $𝒥$ as a sum of J_λλ′, which introduces transitions only between the specific pair of λ and λ′ adiabatic states:

$$\begin{array}{ll}{𝒥}_{\alpha {\alpha }^{\prime },\beta {\beta }^{\prime }}\hfill & =\sum _{\lambda >{\lambda }^{\prime }}{\left({𝒥}_{\lambda {\lambda }^{\prime }}\right)}_{\alpha {\alpha }^{\prime },\beta {\beta }^{\prime }}\hfill \\ \hfill & =\sum _{\lambda >{\lambda }^{\prime }}\{-d_{\lambda {\lambda }^{\prime }}\cdot \frac{P}{M}(({\delta }_{\lambda \alpha }{\delta }_{{\lambda }^{\prime }\beta }-{\delta }_{{\lambda }^{\prime }\alpha }{\delta }_{\lambda \beta }){\delta }_{{\alpha }^{\prime }{\beta }^{\prime }}+(({\delta }_{\lambda {\alpha }^{\prime }}{\delta }_{{\lambda }^{\prime }{\beta }^{\prime }}-{\delta }_{{\lambda }^{\prime }{\alpha }^{\prime }}{\delta }_{\lambda {\beta }^{\prime }}){\delta }_{\alpha \beta })\hfill \\ \hfill & \hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}-\frac{1}{2}\hslash {\omega }_{\lambda {\lambda }^{\prime }}{d}_{\lambda {\lambda }^{\prime }}\cdot \frac{\partial }{\partial P}(({\delta }_{\lambda \alpha }{\delta }_{{\lambda }^{\prime }\beta }-{\delta }_{{\lambda }^{\prime }\alpha }{\delta }_{\lambda \beta }){\delta }_{{\alpha }^{\prime }{\beta }^{\prime }}+(({\delta }_{\lambda {\alpha }^{\prime }}{\delta }_{{\lambda }^{\prime }{\beta }^{\prime }}-{\delta }_{{\lambda }^{\prime }{\alpha }^{\prime }}{\delta }_{\lambda {\beta }^{\prime }}){\delta }_{\alpha \beta })\}\hfill \\ \hfill & =-\frac{P}{M}\cdot d_{\alpha \beta }\left(1+\frac{1}{2}{S}_{\alpha \beta }\cdot \frac{\partial }{\partial P}\right){\delta }_{{\alpha }^{\prime }{\beta }^{\prime }}+\frac{P}{M}\cdot d_{{\beta }^{\prime }{\alpha }^{\prime }}\left(1-\frac{1}{2}{S}_{{\beta }^{\prime }{\alpha }^{\prime }}\cdot \frac{\partial }{\partial P}\right){\delta }_{\alpha \beta }\hfill \end{array}$$

(13)

where d_αβ = 〈α; R| ∂/∂R | β; R〉 and $S_{\alpha \beta }=\hslash {\omega }_{\alpha \beta }{d}_{\alpha \beta }{\left(\frac{P}{M}\cdot d_{\alpha \beta }\right)}^{-1}$. The second equality gives the adiabatic representation of $𝒥$_λλ′. We remark that it is difficult to exactly simulate the term, $𝒥$, involving bath momentum derivatives within the context of a trajectory-based algorithm. Using the identity that $\frac{1}{2}{S}_{\alpha \beta }\cdot \frac{\partial }{\partial P}=\hslash {\omega }_{\alpha \beta }M\cdot \partial /\partial {\left({\hat{d}}_{\alpha \beta }\cdot P\right)}^2$, where M is a diagonal matrix of the masses of the bath particles and d̂_αβ is the unit vector along d_αβ, allows us to employ the momentum-jump approximation:

$$\left(1+\frac{c}{2}{S}_{\alpha \beta }\cdot \frac{\partial }{\partial P}\right)f(P)\approx e^{\frac{c}{2}{S}_{\alpha \beta }\cdot \frac{\partial }{\partial P}}f(P)=e^{c\hslash {\omega }_{\alpha \beta }M\cdot \partial /\partial {\left({\hat{d}}_{\alpha \beta }\cdot P\right)}^2}f(P)=f(P+\Delta P_c)$$

(14)

where c = 1, 2 corresponding to single and double hops, respectively, and $\Delta P_c={\hat{d}}_{\alpha \beta }\text{sgn}\left(\hat{d}\cdot P\right)\sqrt{{\left({\hat{d}}_{\alpha \beta }\cdot P\right)}^2+c\hslash {\omega }_{\alpha \beta }M}-\hat{d}\left(\hat{d}\cdot P\right)$. We have a translation operator with respect to the variable, (d̂_αβ· P)², in the above equation. Decomposing $P=P_{\perp }+P_{\parallel }=P_{\perp }+{\hat{d}}_{\alpha \beta }\text{sgn}\left({\hat{d}}_{\alpha \beta }\cdot P\right)\sqrt{{\left({\hat{d}}_{\alpha \beta }\cdot P\right)}^2}$, it becomes obvious that the translation operator updates P_‖ components by ΔP_c, as presented in Equation (14). This momentum update conserves the energy of surface-hopping trajectories. Apart from technical issues associated with sampling when the algorithm is implemented, this is the only approximation made to QCL evolution. In fact, it is this approximation that gives this algorithm a surface-hopping structure that has some features in common with Tully’s surface-hopping method; however, coherence and decoherence are automatically incorporated in the evolution. The QCLE does not have such sudden momentum changes, and its evolution is described by continuous momentum changes in the course of the evolution. Comparisons of results using this algorithm with exact quantum solutions indicate that the momentum-jump is rarely the source of problems.

Equation (10) admits a formal solution:

$${\hat{B}}_W^{\alpha {\alpha }^{\prime }}(X,t)={\left(e^{iℒt}\right)}_{\alpha {\alpha }^{\prime }\beta {\beta }^{\prime }}{B}_W^{\beta {\beta }^{\prime }}(X)$$

(15)

Thus, our following discussion focuses on evaluating:

$${\left(e^{i\widehat{ℒ}t}\right)}_{\alpha {\alpha }^{\prime },{\alpha }_K{\alpha }_K^{\prime }}=\sum _{\left({\alpha }_1{\alpha }_1^{\prime }\right)\dots \left({\alpha }_K{\alpha }_K^{\prime }\right)}\prod _{j=1}^K{\left(e^{i\widehat{ℒ}\Delta t_j}\right)}_{{\alpha }_{j-1}{\alpha }_{j-1}^{\prime },{\alpha }_j{\alpha }_j^{\prime }}$$

(16)

In the above equation, we simply factorize the propagator into K pieces with Δt_j = t_j − t_j−1 = Δt. In each small time slice, we perform the symmetric Trotter decomposition:

$${\left(e^{i\widehat{ℒ}\Delta t_j}\right)}_{\alpha {\alpha }^{\prime },\beta {\beta }^{\prime }}\approx {𝒲}_{\beta {\beta }^{\prime }}\left(t_{j-1},t_{j-1}+\frac{\Delta t}{2}\right)e^{i{L}_{\beta {\beta }^{\prime }}\Delta t/2}{𝒬}_{\beta {\beta }^{\prime },\alpha {\alpha }^{\prime }}{𝒲}_{\alpha {\alpha }^{\prime }}\left(t_{j-1}+\frac{\Delta t}{2},t_j\right)e^{i{L}_{\alpha {\alpha }^{\prime }}\Delta t/2}$$

(17)

where: $𝒲$_αα′ (t₁, t₂) = e^{iω_αα′(t₂−t₁)}, and:

$${𝒬}_{\alpha {\alpha }^{\prime },\beta {\beta }^{\prime }}={\left(e^{𝒥\Delta t}\right)}_{\alpha {\alpha }^{\prime },\beta {\beta }^{\prime }}={\left(e^{{\sum }_{\lambda >{\lambda }^{\prime }}{𝒥}_{\lambda {\lambda }^{\prime }}\Delta t}\right)}_{\alpha {\alpha }^{\prime },\beta {\beta }^{\prime }}\approx {\left(\prod _{\lambda >{\lambda }^{\prime }}{e}^{{𝒥}_{\lambda {\lambda }^{\prime }}\Delta t}\right)}_{\alpha {\alpha }^{\prime },\beta {\beta }^{\prime }}$$

(18)

We observe that it is possible to express e^{$𝒥$_λλ′Δt} in the following block-diagonal matrix form:

$$e^{{𝒥}_{\lambda {\lambda }^{\prime }}\Delta t}={ℳ}^{\lambda {\lambda }^{\prime }}\oplus {𝒦}_{{\xi }_1}^{\lambda {\lambda }^{\prime }}\dots \oplus {𝒦}_{{\xi }_{N-2}}^{\lambda {\lambda }^{\prime }}\oplus {𝒩}^{\lambda {\lambda }^{\prime }}$$

(19)

where ξ_i is one of the N − 2 adiabatic states other than λ and λ′. In the above equation, $ℳ$ is a four by four matrix, defined with respect to the basis, {(λ, λ), (λ, λ′), (λ′, λ), (λ′,λ′)}:

$${ℳ}^{\lambda {\lambda }^{\prime }}=\left(\begin{array}{cccc}\hfill {\text{cos}}^2(a)\hfill & \hfill -\text{cos}(a)\text{sin}(a){\hat{j}}_{\lambda {\lambda }^{\prime }}\hfill & \hfill -\text{cos}(a)\text{sin}(a){\hat{j}}_{\lambda {\lambda }^{\prime }}\hfill & \hfill {\text{sin}}^2(a){\hat{j}}_{\lambda \to {\lambda }^{\prime }}\hfill \\ \hfill \text{cos}(a)\text{sin}(a){\hat{j}}_{\lambda {\lambda }^{\prime }}\hfill & \hfill {\text{cos}}^2(a)\hfill & \hfill -{\text{sin}}^2(a)\hfill & \hfill -\text{sin}(a)\text{cos}(a){\hat{j}}_{\lambda {\lambda }^{\prime }}\hfill \\ \hfill \text{cos}(a)\text{sin}(a){\hat{j}}_{\lambda {\lambda }^{\prime }}\hfill & \hfill -{\text{sin}}^2(a)\hfill & \hfill {\text{cos}}^2(a)\hfill & \hfill -\text{sin}(a)\text{cos}(a){\hat{j}}_{\lambda {\lambda }^{\prime }}\hfill \\ \hfill {\text{sin}}^2(a){\hat{j}}_{\lambda \to {\lambda }^{\prime }}\hfill & \hfill \text{cos}(a)\text{sin}(a){\hat{j}}_{\lambda {\lambda }^{\prime }}\hfill & \hfill \text{cos}(a)\text{sin}(a){\hat{j}}_{\lambda {\lambda }^{\prime }}\hfill & \hfill {\text{cos}}^2(a)\hfill \end{array}\right)$$

(20)

with a = (P/M) · d_λλ′ Δt, and ĵ _λλ′ and ĵ _λ→λ′ are the momentum-jump operators, $e^{\frac{1}{2}{S}_{\lambda {\lambda }^{\prime }}\frac{\partial }{\partial P}}$ and $e^{S_{\lambda {\lambda }^{\prime }}\frac{\partial }{\partial P}}$, defined in Equation (14) with c = 1, 2, respectively. In Equation (19), there exists another set of four by four matrices, ${𝒦}_{{\xi }_i}^{\lambda {\lambda }^{\prime }}$, with i = 1, …, N − 2. Each of these matrices is defined with respect to a basis of the form, {(λ, ξ_i), (λ′, ξ_i)} ⊕ {(ξ_i, λ), (ξ_i, λ′)}:

$${𝒦}_{\xi }^{\lambda {\lambda }^{\prime }}=\left(\begin{array}{cc}\hfill \text{cos}(a)\hfill & \hfill -\text{sin}(a){\hat{j}}_{\lambda {\lambda }^{\prime }}\hfill \\ \hfill \text{sin}(a){\hat{j}}_{\lambda {\lambda }^{\prime }}\hfill & \hfill \text{cos}(a)\hfill \end{array}\right)\oplus \left(\begin{array}{cc}\hfill \text{cos}(a)\hfill & \hfill -\text{sin}(a){\hat{j}}_{\lambda {\lambda }^{\prime }}\hfill \\ \hfill \text{sin}(a){\hat{j}}_{\lambda {\lambda }^{\prime }}\hfill & \hfill \text{cos}(a)\hfill \end{array}\right)$$

(21)

Finally, there is a null matrix, $𝒩$^λλ′, of a size of (N − 2)², and the associated null space is spanned by basis vectors, (ξ₁, ξ₂), where ξ_i ≠ λ^(′). We remark that one has to permute the basis vectors in order to construct these block-diagonal matrices [26].

At this point, we have specified all the necessary details in order to simulate the QCL dynamics in the adiabatic basis:

$$B_W^{\alpha {\alpha }^{\prime }}(X,t)=\sum _{({\alpha }_1{\alpha }_1^{\prime }),\dots ,({\alpha }_K{\alpha }_K^{\prime })}\left[\prod _{j=1}^K{𝒲}_{{\alpha }_{j-1}{\alpha }_{j-1}^{\prime }}{e}^{i{L}_{{\alpha }_{j-1}{\alpha }_{j-1}^{\prime }}\Delta t}{Q}_{{\alpha }_{j-1}{\alpha }_{j-1}^{\prime }.{\alpha }_j{\alpha }_j^{\prime }}{𝒲}_{{\alpha }_j{\alpha }_j^{\prime }}{e}^{i{L}_{{\alpha }_j{\alpha }_j^{\prime }}\Delta t}\right]B_W^{{\alpha }_K{\alpha }_K^{\prime }}(X)$$

(22)

where ${\alpha }_0^{(\text{'})}={\alpha }^{(\text{'})}$. The explicit summation over all quantum indices, ( ${\alpha }_1{\alpha }_1^{\prime }$) ⋯ ( ${\alpha }_K{\alpha }_K^{\prime }$), can also be evaluated stochastically. For instance, given a pair of indices, (α_jα_j−1), one can determine the next pair at the time slice, j + 1, by drawing an MC sample from the transition probability:

$$P({\alpha }_{j+1},{\alpha }_{j+1}^{\prime }|{\alpha }_j,{\alpha }_j^{\prime })=\frac{\left|{𝒬}_{{\alpha }_j{\alpha }_j^{\prime },{\alpha }_{j+1},{\alpha }_{j+1}^{\prime }}\right|}{{\sum }_{{\beta }_{j+1},{\beta }_{j+1}^{\prime }}\left|{𝒬}_{{\alpha }_j{\alpha }_j^{\prime },{\beta }_{j+1},{\beta }_{j+1}^{\prime }}\right|}$$

(23)

If the sampled new pair of indices differs from the starting pair, then the sampled $𝒬$ matrix element must contain the proper momentum-jump operators to update the energy of the trajectory after the jump. In any actual implementation of this algorithm, it is desirable to restrict to nonadiabatic transitions between one pair of states in every time slice. Under this assumption, one can then approximate:

$${𝒬}_{\alpha {\alpha }^{\prime },\beta {\beta }^{\prime }}\approx \{\begin{array}{ll}{\delta }_{\alpha \beta }{\delta }_{{\alpha }^{\prime }{\beta }^{\prime }}\hfill & \text{if}\hspace{0.17em}\text{no}\hspace{0.17em}\text{hop}\hspace{0.17em}\text{happens},\hfill \\ {\left(e^{{𝒥}_{\mu \gamma }}\right)}_{\alpha {\alpha }^{\prime },\beta {\beta }^{\prime }}\hfill & \text{if\hspace{0.17em}}(\alpha ,{\alpha }^{\prime })\to (\beta ,{\beta }^{\prime })\hspace{0.17em}\text{involves}\hspace{0.17em}\text{transition}\hspace{0.17em}\text{between}\hspace{0.17em}(\mu ,\gamma )\hspace{0.17em}\text{states},\hfill \\ 0\hfill & \text{if}\hspace{0.17em}(\alpha ,{\alpha }^{\prime })\to (\beta ,{\beta }^{\prime })\hspace{0.17em}\text{involves}\hspace{0.17em}\text{transition}\hspace{0.17em}\text{between}\hspace{0.17em}\text{two}\hspace{0.17em}\text{or}\hspace{0.17em}\text{more}\hspace{0.17em}\text{pair}\hspace{0.17em}\text{of}\hspace{0.17em}\text{states}.\hfill \end{array}$$

(24)

In this algorithm, we see that the trajectories in the ensemble that are used to simulate the time evolution are non-Newtonian in character, consisting of Newtonian segments where the system evolves on adiabatic surfaces, or the mean of two adiabatic surfaces, interspersed with quantum transitions and momentum changes.

3.2. Forward-Backward Trajectory Solution

This scheme is motivated by another way of writing the formally exact solution [38] of the QCLE using the last line of Equation (8):

$${\hat{B}}_W(X,t)=𝒮\left(e^{i{\overleftarrow{\mathscr{H}}}_{\Lambda }t/\hslash }{\hat{B}}_W(X)e^{-i{\overleftarrow{\mathscr{H}}}_{\Lambda }t/\hslash }\right)$$

(25)

The $𝒮$ operator [31,38] specifies the order in which the forward and backward evolution operators act on B̂_W (X). The ordering of evolution operators is critical because of the lack of an underlying Lie algebraic structure [38] of the QCLE.

One approach to solve Equation (25) is to apply the mapping transformation in which N discrete quantum states of the subsystem are represented by the continuous position and momenta of N fictitious harmonic oscillators. The properties of the original subsystem are then obtained via an ensemble average involving trajectories in the phase space of the fictitious oscillators. More precisely, in the mapping representation, a subsystem state, |λ〉, is replaced by |m_λ〉 = |0₁, …, 1_λ, … 0_N〉, a product state specifying the occupation numbers (limited to zero or one) of N fictitious harmonic oscillators [33,49]. Creation and annihilation operators, ${\hat{a}}_{\lambda }^{†}$ and â_λ, satisfy the commutation relation $\left[{\hat{a}}_{\lambda },{\hat{a}}_{{\lambda }^{\prime }}^{†}\right]={\delta }_{\lambda {\lambda }^{\prime }}$ for harmonic oscillators. The actions of these operators on the single-excitation mapping states are ${\hat{a}}_{\lambda }^{†}|0\rangle =|m_{\lambda }\rangle$ and â_λ|m_λ〉 = |0〉, where |0〉 = ||0₁ . . . 0_N 〉 is the ground state of the mapping basis.

Next, we define the mapping version of operators, ${\hat{B}}_m(X)=B_W^{\lambda {\lambda }^{\prime }}(X){\hat{a}}_{\lambda }^{†}{\hat{a}}_{{\lambda }^{\prime }}$, such that matrix elements of B̂_W in the subsystem basis are equal to the matrix elements of the corresponding mapping operator: $B_W^{\lambda {\lambda }^{\prime }}(X)=\langle \lambda |{\hat{B}}_W(X)|{\lambda }^{\prime }\rangle =\langle m_{\lambda }|{\hat{B}}_m(X)|m_{{\lambda }^{\prime }}\rangle$. In particular, the mapping Hamiltonian is:

$${\hat{H}}_m=H_b(X)+h^{\lambda {\lambda }^{\prime }}(R){\hat{a}}_{{}^{{}_{\lambda }}}^{†}{\hat{a}}_{{\lambda }^{\prime }}\equiv H_b(X)+{\hat{h}}_m$$

(26)

where we applied the mapping transformation only on the part of the Hamiltonian that involves the subsystem DOF in Equation (26). The mapping Hamiltonian,ĥ_m, is always a quadratic Hamiltonian with respect to the quantum DOF. The pure bath term,Ĥ_b(X), acts as an identity operator in the subsystem basis and is mapped onto the identity operator of the mapping space directly. The mapped formal solution of QCLE now reads:

$${\hat{B}}_W(X,t)=𝒮\left(e^{i\overrightarrow{{\mathscr{H}}_{\Lambda }^m}t/\hslash }{\hat{B}}_W(X)e^{-i\overleftarrow{{\mathscr{H}}_{\Lambda }^m}t/\hslash }\right)$$

(27)

where $\overrightarrow{{\mathscr{H}}_{\Lambda }^m}$ is given by $\overrightarrow{{\mathscr{H}}_{\Lambda }^m}={\hat{H}}_m\left(1+\hslash \Lambda /2i\right)$, with an analogous definition for $\overleftarrow{{\mathscr{H}}_{\Lambda }^m}$.

We now introduce the coherent states, |z〉, in the mapping space, â_λ|z〉 = z_λ|z〉 and $\langle z|{\hat{a}}_{\lambda }^{†}=z_{\lambda }^{*}\langle z|$, where |z〉 = |z₁, . . . , z_N 〉, and the eigenvalue is $z_{\lambda }=\left(q_{\lambda }+i{p}_{\lambda }\right)/\sqrt{\hslash }$. The variables q = (q₁, . . . , q_N) and p = (p₁, . . . , p_N) are mean coordinates and momenta of the harmonic oscillators encoded in the coherent state, |z〉, respectively. The coherent states form an overcomplete basis with the inner product between any two such states, 〈z| z′〉 = e^{−(|z−z′|2)−i(z·z′*−z* ·z′)}. Finally, we remark that the coherent states provide the resolution of identity:

$$ℐ=\int \frac{d^{2}z}{{\pi }^N}|z\rangle \langle z|$$

(28)

where d²z = d(ℜ(z))d($𝒯$(z)) = dqdp/(2ℏ)^N.

Similar to the path integral approach for solving the quantum dynamics, we decompose the forward and backward evolution operators in Equation (27) into a concatenation of M short-time evolutions with Δt_i = τ and Mτ = t. In each short-time interval, Δt_i, we introduce two sets of coherent states, |z_i〉 and $|z_i^{\text{'}}\rangle$, via Equation (28) to expand the forward and backward time evolution operators, respectively. The time evolution (generated by a quadratic Hamiltonian) of coherent states can be represented by trajectory evolution in the phase space of (q, p). After some algebra, the matrix elements of Equation (27) can be approximated by:

$$\begin{array}{l}{B}_W^{\lambda {\lambda }^{\prime }}(X,t)=\sum _{\mu {\mu }^{\prime }}\int dxd{x}^{\prime }\varphi (x)\varphi (x^{\prime })\frac{1}{\hslash }\left(q_{\lambda }+i{p}_{\lambda }\right)\left({q^{\prime }}_{{\lambda }^{\prime }}+i{p^{\prime }}_{{\lambda }^{\prime }}\right)B_W^{\mu {\mu }^{\prime }}(X_t)\\ \hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\times \frac{1}{\hslash }\left(q_{\mu }(t)+i{p}_{\mu }(t)\right)\left({q^{\prime }}_{{\mu }^{\prime }}(t)+i{p^{\prime }}_{{\mu }^{\prime }}(t)\right)\end{array}$$

(29)

where x = (q, p) gives the real and imaginary parts of z, dx = dqdp and $\varphi (x)={\left(\hslash \right)}^{-N}{e}^{-{\sum }_{\nu }\left(q_{\nu }^2+p_{\nu }^2\right)/\hslash }$ is the normalized Gaussian distribution function. In deriving Equation (29), we have invoked an orthogonality approximation on the inner product between subsequent coherent state variables, $\langle z_i\left|e^{\frac{i}{\hslash }\hat{h}t}\right|z_{i+1}\rangle =\langle z_i(t)|z_{i+1}\rangle \approx {\pi }^N\delta \left(z_{i+1}-z_i(t_i)\right)$, with i being the time step index. This approximation is necessary to construct a continuous trajectory of z(t). In the extended phase space of (X(t), z(t), z′(t)), the trajectories follow Hamiltonian dynamics:

$$\frac{d{\chi }_{\mu }}{dt}=\frac{\partial H_e(\chi ,\pi )}{\partial {\pi }_{\mu }},\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\frac{d{\pi }_{\mu }}{dt}=-\frac{\partial H_e(\chi ,\pi )}{\partial {\chi }_{\mu }}$$

(30)

where $H_e\left(\chi ,\pi \right)=P^2/2M+V_0(R)+\frac{1}{2\hslash }{h}^{\lambda {\lambda }^{\prime }}(R)(q_{\lambda }{q}_{{\lambda }^{\prime }}+p_{\lambda }{p}_{{\lambda }^{\prime }}+{q^{\prime }}_{\lambda }{q^{\prime }}_{{\lambda }^{\prime }}+{p^{\prime }}_{\lambda }{p^{\prime }}_{{\lambda }^{\prime }})$ with V₀(R) = V_b(R) − T rĥ(R), χ = (R, q, q′) and π = (P, p, p′). We remark that the FBTS trajectories manifestly conserve energy. Furthermore, simulating the dynamics with a standard velocity Verlet type of symplectic integrator has a stationary solution proportional to H_pseudo = H_e(χ, π) + Δt²δH, as discussed in [35].

The main approximation introduced in the derivation of the FBTS, Equation (29), is the orthogonality approximation. The simplest improvement to the algorithm is to refrain from applying this approximation at every time step. In [36], we outlined a practical approach to evaluate the set of selected integrals of z_i and $z_i^{\text{'}}$ (which could be evaluated analytically if the orthogonality approximation were applied). We termed this extension of FBTS as the jump FBTS (JFBTS). Since the computational cost grows quickly with respect to the number of jumps inserted, one needs to make a trade-off between numerical efficiency and accuracy.

In the simplest approach, one selects every (M/K) time step from a total of M steps to fully evaluate the coherent state integrals:

$$\begin{array}{l}{B}_W^{\lambda {\lambda }^{\prime }}(X,t)=\sum _{\mu {\mu }^{\prime }}\sum _{\begin{array}{l}{s}_0{s^{\prime }}_0\dots \\ s_{K-1}{s^{\prime }}_{K-1}\end{array}}\int \prod _{v=0}^{K}dxd{x}^{\prime }\varphi (x_v)\varphi ({x^{\prime }}_v)\\ \times \frac{1}{\hslash }(q_{0\lambda }+i{p}_{0\lambda })({q^{\prime }}_{0\lambda }+i{p^{\prime }}_{0\lambda })B_W^{\mu {\mu }^{\prime }}(X_t)\\ \times \frac{1}{\hslash }\left\{\prod _{v=1}^K\left(q_{(v-1)s_{v-1}}({\tau }_v)-i{p}_{(v-1)s_{v-1}}({\tau }_v)\right)\left(q_{v{s}_v}+i{p}_{v{s}_v}\right)\right\}\\ \times \frac{1}{\hslash }\left\{\prod _{v=1}^K\left({q^{\prime }}_{(v-1)s_{v-1}}({\tau }_v)+i{p^{\prime }}_{(v-1)s_{v-1}}({\tau }_v)\right)\left({q^{\prime }}_{v{s}_v}+i{p^{\prime }}_{v{s}_v}\right)\right\}\\ \times \frac{1}{\hslash }\left(q_{K\mu }\left({\tau }_{K+1}\right)-i{p}_{K\mu }\left({\tau }_{K+1}\right)\right)\left({q^{\prime }}_{K{\mu }^{\prime }}\left({\tau }_{K+1}\right)-i{p^{\prime }}_{K{\mu }^{\prime }}\left({\tau }_{K+1}\right)\right)\end{array}$$

(31)

where the subscripts, v and s, refer to the v-th time step and the s-th component of the q and p vectors, respectively, and τ_v = t_iv − t_{iv−1 with ti0} = 0 and t_{iK+1 = t}. According to this prescription, the continuous FB trajectories experience K discontinuous jumps in the (x, x′) phase space. Between subsequent jumps, the evolution of the FB trajectory is governed by Equation (30). Simulations show that with a sufficient number of jumps, numerically exact solutions of the QCLE can be obtained [36].

3.3. Comparisons between Algorithms

The differences between the two QCLE simulation algorithms can be traced to the quantum basis that is used and the way that feedback between quantum and classical systems is treated. In the case of the TBSH algorithm, the trajectories are propagated through a Hellmann-Feynman force, or the mean of two Hellmann-Feynman forces [Equation (12)], with intermittent surface hops that switch the adiabatic surfaces on which the trajectories propagate. In the case of FBTS, one not only propagates the bath dynamical variables as trajectories, but also the quantum dynamical variables, which are associated with fictitious harmonic oscillators. In this extended phase space, we have exact Hamiltonian dynamics. In particular, the force acting on the bath particles simultaneously involves all N adiabatic surfaces, which is similar to, but different from, the Ehrenfest mean-field approach. The very different characteristics of the trajectories in two algorithms manifest the artificial character of the trajectory dynamics. Thus, one should not attach physical significance to single trajectories in the computation. All physical properties of the system can only be extracted from a proper ensemble average of a large set of trajectories, as implied in Equation (2). Nevertheless, insight into the trajectory dynamics of each algorithm will help to judge the simulation efficiency for various classes of models.

For certain problems, such as proton transfer reactions, where the time scales of the bath and subsystem are well-separated, even during nonadiabatic transitions, the TBSH algorithm can yield quantitatively accurate results with a few hops. There are also dynamical problems in which distinct bath motions can be explicitly correlated with the subsystem’s quantum states. For instance, in the simple Tully I model [35,50], trajectories populated on the excited state will cross the avoided crossing point, while the ground state trajectories will eventually be reflected and retrace their paths in the opposite direction. This kind of behavior is, however, completely missed when one propagates trajectories in a single effective mean field. Again, the inherent multi-configuration nature of surface-hopping-like algorithms is a more appropriate choice for this case. However, a recent study [51] has indicated that the “jump” version of mean-field-like algorithms can improve the simulation results in cases of this type.

Alternatively, there are also many examples where one would expect FBTS to be the preferred simulation method. In general, the TBSH algorithm has convergence issues, as the MC weights associated with nonadiabatic hops grows rapidly. Even for the simple spin boson model, one can identify parameter regimes where this numerical instability is clearly observed. In these cases, the FBTS and JFBTS are certainly the alternatives that one should adopt for efficient simulations.

4. An Example: Quartic Oscillator in a Harmonic Bath

As a specific example to illustrate the formalism outlined above, we consider a two-level system coupled to a quartic bistable oscillator with a single pair of phase space coordinates X₀ = (R₀, P₀). The quartic oscillator is, in turn, coupled to an Ohmic heat bath of N_b independent harmonic oscillators with phase space coordinates X_i = (R_i, P_i) and i = 1 . . . N_b. The partially Wigner transformed Hamiltonian, expressed in the diabatic basis, {|R〉, |L〉}, reads:

$${\hat{H}}_W=\left(\begin{array}{ll}\hslash {\gamma }_0{R}_0\hfill & -\hslash \Omega \hfill \\ -\hslash \Omega \hfill & -\hslash {\gamma }_0{R}_0\hfill \end{array}\right)+\left(\frac{P_0^2}{2M_0}+V_n(R_0)+\sum _{j=1}^{N_b}\frac{P_j^2}{2M_j}+\frac{M_j{\omega }_j^2}{2}{\left(R_j-\frac{{\gamma }_b{c}_j}{M_j{\omega }_j^2}{R}_0\right)}^2\right)\mathbf{I}$$

(32)

where $V_n(R_0)=-M_0{\omega }_0^2{R}_0^2/2+A{R}_0^4/4$ and I is an identity matrix. We take N_b = 40 harmonic oscillators for the discretization of the Ohmic heat bath. Following the discretization scheme introduced in [52], we set ω_j = ω_c ln (1−jω_c/δω) and c_j = (ξℏωM_j)^1/2ω_j with δω = (1−exp(ω_max/ω_c))/N_b. The parameters, ω_c and ω_max, are the characteristic and cut-off frequencies for the Ohmic bath, respectively. The Kondo parameter is ξ.

The adiabatic states for the subsystem are:

$$\begin{array}{l}|+;R_0\rangle =\frac{1}{𝒩(R_0)}[(1-G)|R\rangle -(1+G)|L\rangle ]\\ |-;R_0\rangle =\frac{1}{𝒩(R_0)}[(1+G)|R\rangle +(1-G)|L\rangle ]\end{array}$$

(33)

where $𝒩(R_0)=\sqrt{2\left(1+G^2(R_0)\right)}$ and $G(R_0)={\left({\gamma }_0{R}_0\right)}^{-1}\left[-\Omega +\sqrt{{\Omega }^2+{\gamma }_0^2{R}_0^2}\right]$. The adiabatic energies are given by $E_{\pm }=V_n(R_0)\pm \hslash \sqrt{{\Omega }^2+{\gamma }_0^2{R}_0^2}=V_n(R_0)\pm {∊}_{\pm }(R_0)$.

We shall study the autocorrelation functions, C_LL, with  = B̂ = |L〉 〈L|. The entire system is assumed to be in thermal equilibrium initially. Using the high-temperature approximation presented in the Appendix, the correlation function of interest can be given in a compact form:

$$\begin{array}{ll}{C}_{LL}(t)=\hfill & \int d{X}_{0}d{X}_b𝒲(R_0)𝒢\left(P_0;\frac{M_0}{\beta }\right)\prod _{j=1}^{N_b}𝒢\left(P_j;\frac{M_j}{\beta }\right)𝒢\left(R_j-\frac{{\gamma }_b{c}_j}{M_j{\omega }_j^2}{R}_0;\frac{1}{\beta M_j{\omega }_j^2}\right)\hfill \\ \hfill & \times \sum _{n=L,R}\sum _{\alpha ,{\alpha }^{\prime }}{F}_{\alpha {\alpha }^{\prime }}(X_0)\langle n|\alpha ;R_0\rangle \langle {\alpha }^{\prime };R_0|L\rangle B_W^{\mathit{Ln}}(X_t)\hfill \end{array}$$

(34)

where $𝒢$(x; σ²) = (2πσ²)^−1/2 e^−x2/2σ2, and:

$$𝒲(R_0)=\frac{e^{-\beta \left(\frac{A}{4}{R}_0^4-\frac{1}{2}{M}_0{\omega }_0^2{R}_0^2\right)}\left(e^{-\beta {∊}_{+}(R_0)}+e^{\beta {∊}_{-}(R_0)}\right)}{\int d{R}_0{e}^{-\beta \left(\frac{A}{4}{R}_0^4-\frac{1}{2}{M}_0{\omega }_0^2{R}_0^2\right)}\left(e^{-\beta {∊}_{+}(R_0)}+e^{\beta {∊}_{-}(R_0)}\right)}$$

(35)

An MC evaluation of the integrals can be done by sampling P₀, R_b, P_b from the Gaussian distributions and sampling R₀ from $𝒲$(R₀), respectively. The time-evolved matrix element, $B_W^{nm}\left(X_t\right)$, will be computed using both the TBSH and the FBTS algorithms. Finally, we note that the path-integral-based sampling scheme introduced in Section 2 should be adopted to sample phase-space points from (ρ̂_eq)_W(X) for more generalized situations, including cases of low-temperature, arbitrary subsystem-bath divisions of a composite system, strong subsystem-bath couplings and an arbitrary potential energy profile.

In this study, we report numerical results in the energy unit, ℏω_c, and distance unit, $\sqrt{\hslash /M_j{\omega }_c}$, for each environmental DOF. We consider two sets of parameters. In the first case, we use the following parameter values, a = 1.0, ω₀ = 1.2, γ₀ = 0.05 γ_b = 1.0, Ω = 0.3, ξ = 0.1, ω_max = 3 and β = 0.2, in the dimensionless units. Figure 1a presents the potential surface profiles [53], W_α(R₀). The two diabatic surfaces, W_L,R(R₀), remain close to each other, and the two adiabatic surfaces, W_±(R₀), share essentially the same characteristics. In this case, a mean-field-based algorithm, like FBTS, should be accurate and efficient. This problem can also be handled easily in the adiabatic basis, since the surface-hopping trajectories will be initialized in both the adiabatic ground and excited states, because the system is in a thermal equilibrium state at t = 0. Furthermore, the coupling parameter, γ_b, was purposely chosen to be small in order to minimize the number of nonadiabatic transitions (or hops) encountered in the TBSH algorithm. In panel (b), C_LL(t) is computed using both algorithms. The agreement between these results is good.

Next, we consider the following parameter set, a = 0.8, ω₀ = 0.6, γ₀ = 0.3 γ_b = 1.0, Ω = 0.1, ξ = 0.1, ω_max = 3, and β = 0.2 in the dimensionless units. Figure 2a shows the potential surface profiles, W_α(R₀), obtained from this set of parameters. In this case, the adiabatic, W_±(R₀), and diabatic surfaces, W_L,R(R₀), only differ markedly near the region of the barrier top, where an avoided crossing point indicates significant mixing of the two diabatic states. Nonadiabatic effects should be most prominent near this barrier top. A stronger coupling, γ₀, is also chosen in this case. Figure 2b presents the autocorrelation functions. In the main figure of panel (b), the blue curves (C_LL(t) computed by the FBTS) start with the full correlation at one, then gradually reduce to 1/2, which implies that the subsystem is in an equal admixture of the two diabatic states in the asymptotic limit. The TBSH simulation results are only valid for very short times (as shown in the inset of the Figure 2b), due to instability arising from the accumulation of weights, even with filtering [54]. The thermal equilibrium distribution, $𝒲$(R₀), has a bimodal distribution profile, as illustrated in Figure 2a; however, for the (inverse) temperature, β = 0.2, the double-peaked structure is very broad. The $𝒲$(R₀) distribution profile (blue curve in Figure 2a) suggests that the thermal equilibrium state has a non-trivial contribution from the excited surface. Sampling from $𝒲$(R₀) yields many R₀ values near the barrier top, where several hops immediately take place for this strong-coupling case, and the instability sets in early in the simulation. Lowering β will produce a more pronounced double-peak structure for $𝒲$(R₀), but the quartic oscillator’s momentum, P₀, will fluctuate with a larger variance in the presence of the heat bath in this case. Since nonadiabatic transitions depend non-trivially on a = P₀· d₁₂(R₀)Δt in the TBSH algorithm, large momentum fluctuations will eventually affect the long-time result. This case shows some of the practical limitations of the TBSH algorithm for the computation of this correlation function.

5. Conclusions

The scheme for computing the quantum correlation function in Equation (2) combines a numerically exact quantum initial sampling method with dynamics described by the QCLE; thus, the approximations in the simulation method reside in the dynamics. It is easier to compute the equilibrium properties of a quantum system, for instance, by using the imaginary-time Feynman path integral method, than to obtain dynamical properties by using similar real-time Feynman path integrals without adopting further approximations. Since we approximate the quantum dynamics of the entire system, quantum subsystem plus bath, by QCL dynamics, it is appropriate to comment on some of its features.

It is known that the quantum-classical bracket, defined in terms of the commutator and Poisson brackets in Equation (8), does not possess a Lie algebraic structure, since it fails to satisfy the Jacobi identity [2,38]. This lack of a proper algebraic structure is shared by all known MQC methods and simply reflects the inconsistency in mixing classical and quantum mechanical dynamics. One consequence of this inconsistency is that the partial Wigner transform, ρ̂W_e(R, P), of the full canonical equilibrium density function, ρ̂_eq = e^−βĤ/Z_Q, is not stationary under the QCLE; however, ρ̂W_e(R, P) can be written as an expansion in the mass ratio (or ℏ), and it has been shown that the full quantum equilibrium density is conserved under the QCL dynamics up to $𝒪$(ℏ). Therefore, the detailed balance relation is also satisfied to this order. The violation of a detailed balance is a common problem that affects all major MQC methods, including the two most popular approaches, Ehrenfest mean-field [55] and Tully’s fewest switching surface hopping [56] (FSSH), to various degrees. Of course, as noted earlier, for the class of models where an arbitrary quantum system is bilinearly coupled to a harmonic bath, the dynamics is exact, and a detailed balance is exactly satisfied.

The dynamics described by the QCLE can be related to that prescribed by other methods. In [57], it was shown that one could derive both Ehrenfest mean-field dynamics and a version of surface-hopping dynamics starting from the QCLE. In the former case, one simply drops all the “correlations” (including entanglement) between the subsystem and bath densities in the QCLE [58]. In the later case, one projects out all the off-diagonal matrix elements of the density in the QCLE to obtain a generalized master equation for the subsystem alone. Then, one considers decoherence to suppress the coherences in order to recover a simple “surface hopping” dynamics [59] similar to that prescribed in the FSSH algorithm. Furthermore, it had been proven [60] that the QCLE and the partially linearized path integral (PLPI) method [61–64] share the same starting mathematical foundation. In particular, the most recent PLPI algorithm, called PLDM (Partially Linearized Density Matrix) method [64], is very similar to the FBTS presented in this paper [31]. One can also draw comparisons between methods based on the QCLE and semiclassical initial value representations. For instance, numerical schemes based on the Poisson bracket mapping equation (PBME) [30], an approximate equation derived from the mapping-transformed QCLE, and the linearized semiclassical initial value representations [65] share the same set of equations of motion for the trajectories.

Mixed quantum-classical methods are often the only feasible approach to explore the dynamics of large complex systems, such as condensed phase or biochemical systems, where only a few light-mass DOF need be treated quantum mechanically. In many rate processes of interest, such as electron transfer or proton transfer, the local polar solvent motions are responsible for important features of the reaction mechanism. As a result, it is essential that the dynamics of these environmental degrees of freedom be treated in detail. Open quantum system methods that trace out all bath details cannot capture important aspects of such dynamics.

Some recent work [48,66] has suggested interesting ways to combine the QCLE and the generalized master equation [67–69] approach. Simulation tests on spin boson models [48] and a two-level system coupled to an anharmonic bath [68] indicate that accurate, long-time dynamical properties of such systems can be efficiently calculated with an improved memory kernel (which takes the short-time QCLE computation of some bath correlation functions as the input) for the general master equation. This type of hybrid approach may eventually prove to be useful for studies of more complex systems.

Finally, we provide comments that may help in choosing between the two algorithms for simulations. The TBSH algorithm, without filtering, provides a very accurate QCL dynamics before the onset of the sign problem associated with its heavy reliance on Monte Carlo sampling. While filtering can be used to extend simulations to much longer times, the problems related to Monte Carlo sampling limits its usefulness in performing long-time simulations, as vividly illustrated in Section 4. However, the TBSH is found to be the preferred simulation method (in comparison to the FBTS) when one investigates bath dynamical properties of systems in the vicinity of conical intersections and avoided crossings. For instance, the TBSH results accurately capture the intricate geometric phases [46] and the bimodal structure in the momentum distribution [35] in the Tully 1 model (a single avoided crossing model), while the FBTS fails to reproduce these delicate features, even though it provides fairly accurate population dynamics, as reported in [36]. Since the FBTS trajectory dynamics is based on a mean-field description, one finds that the results are usually very accurate (even in the long-time limit) when the energy gap between diabatic energy surfaces is small in comparison to the typical subsystem-bath coupling strength. Another advantage of the FBTS is the availability of the JFBTS [36] algorithm, which implements systematic correction of FBTS results towards the exact QCL dynamics and provides a simple method to gauge the sufficiency of the FBTS results.