Correlation Functions in Open Quantum-Classical Systems

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.


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]: 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][4][5].In this paper, we focus on schemes based on the system-bath partition, and using this partition, the Hamiltonian reads: Ĥ = Ĥb + ĥs + Vc ( R); where Ĥb = P 2 2M + Vb ( 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 ĥs = p2 2m + Vs (r) in some quantum basis: ĥs = ij |i i| ĥs |j j|.
Several methods based on various master equations [6][7][8][9][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][14][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]: 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 i h P •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, ρeq Â 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 nm W (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][24][25][26][27][28][29][30][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][35][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.

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, ρeq Â First, we recall the definition of partial Wigner transform: where R represents the vector of bath coordinates, n denotes a basis state for the subsystem and Ĥ = 2M + Vb ( R) + ĥ( R) with ĥ( R) ≡ ĥs + Vc ( R).One way to compute the integral on the right side of Equation ( 3) is to first factorize e −β Ĥ = e −β L Ĥ into L − 1 pieces with β L = β/(L − 1).Following the standard procedures for path integral calculations, we then insert resolutions of the identity, I = dR i m i |m i , R i m i , R i |, between every pair of factorized operators and apply the approximation, e −β L Ĥ ≈ e −β L P 2 2M e −β L( Vb ( R)+ ĥ( R)) .The integrand on the right side of Equation ( 3) can then be written as follows: where: which is correct to order O(β 2 L ).Substituting Equation (4) into Equation (3), the new integrand of the Wigner transform becomes 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: 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

Quantum-Classical Liouville Equation
In this section, we discuss how one can simulate the time-evolved matrix elements, B n 2 n 1 W (X, t), in Equation (2) using the QCLE: where 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 BW (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, L. Finally, the third equality is another equivalent representation of QCLE in terms of the forward and backward mixed quantum-classical Hamiltonians: 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][40][41][42][43][44][45][46][47][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.

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 2M is taken to be the adiabatic Hamiltonian for a static configuration of R in this section.In the adiabatic basis, the QCLE reads: where the matrix elements of the QCL superoperator are given by: with ω αα = (E α − E α )/h.(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, iL, may be separated into two contributions: The classical propagator is defined as: where ∂R |α; R is the Hellmann-Feynman force.The superoperator, J αα ,ββ , 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 J as a sum of J λλ , which introduces transitions only between the specific pair of λ and λ adiabatic states: where The second equality gives the adiabatic representation of J λλ .We remark that it is difficult to exactly simulate the term, J , involving bath momentum derivatives within the context of a trajectory-based algorithm.Using the identity that , 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: where c = 1, 2 corresponding to single and double hops, respectively, and We have a translation operator with respect to the variable, ( dαβ • P ) 2 , in the above equation.Decomposing P = P ⊥ + P = , 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: Thus, our following discussion focuses on evaluating: 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: , t j e iL αα ∆t/2 (17) where: W αα (t 1 , t 2 ) = e iω αα (t 2 −t 1 ) , and: We observe that it is possible to express e J λλ ∆t in the following block-diagonal matrix form: where ξ i is one of the N − 2 adiabatic states other than λ and λ .In the above equation, M is a four by four matrix, defined with respect to the basis, {(λ, λ), (λ, λ ), (λ , λ), (λ , λ )}: ∂P and e S λλ ∂ ∂P , defined in Equation ( 14) with c = 1, 2, respectively.In Equation ( 19), there exists another set of four by four matrices, K λλ ξ i , with i = 1, . . ., N − 2. Each of these matrices is defined with respect to a basis of the form, {(λ, ξ i ), (λ , ξ i )} ⊕ {(ξ i , λ), (ξ i , λ )}: Finally, there is a null matrix, N λλ , of a size of (N − 2) 2 , and the associated null space is spanned by basis vectors, (ξ 1 , ξ 2 ), 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: where α The explicit summation over all quantum indices, (α 1 α 1 ) . . .(α K α K ), 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: If the sampled new pair of indices differs from the starting pair, then the sampled Q 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: ) involves transitions between two or more pairs of states. ( 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.

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): The S operator [31,38] specifies the order in which the forward and backward evolution operators act on BW (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 Next, we define the mapping version of operators, Bm (X) = B λλ W (X)â † λ âλ , such that matrix elements of BW in the subsystem basis are equal to the matrix elements of the corresponding mapping operator: In particular, the mapping Hamiltonian is: 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: where ), with an analogous definition for ← H m Λ .We now introduce the coherent states, |z , in the mapping space, âλ |z = z λ |z and z| â † λ = z * λ z|, where |z = |z 1 , . . ., z N , and the eigenvalue is z λ = (q λ + ip λ )/ √ h.The variables q = (q 1 , . . ., q N ) and p = (p 1 , . . ., 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: where 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 , 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: where x = (q, p) gives the real and imaginary parts of z, dx = dqdp and φ 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, , 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: where 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 2 δ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 (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: 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 i v−1 with t i 0 = 0 and t i K+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].

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 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.

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 0 = (R 0 , P 0 ).
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: where 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 = (ξhδω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: where The adiabatic energies are given by ).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: where G(x; σ 2 ) = (2πσ 2 ) −1/2 e −x 2 /2σ 2 , and: An MC evaluation of the integrals can be done by sampling P 0 , R b , P b from the Gaussian distributions and sampling R 0 from W(R 0 ), respectively.The time-evolved matrix element, B nm W (X t ), 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, hω c , and distance unit, h/M j ω c , for each environmental DOF.We consider two sets of parameters.In the first case, we use the following parameter values, a = 1.0, ω 0 = 1.2, γ 0 = 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 0 ).The two diabatic surfaces, W L,R (R 0 ), remain close to each other, and the two adiabatic surfaces, W ± (R 0 ), 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 = 0.6, γ 0 = 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 0 ), obtained from this set of parameters.In this case, the adiabatic, W ± (R 0 ), and diabatic surfaces, W L,R (R 0 ), 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, γ 0 , 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, W(R 0 ), 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 W(R 0 ) distribution profile (blue curve in Figure 2a) suggests that the thermal equilibrium state has a non-trivial contribution from the excited surface.Sampling from W(R 0 ) yields many R 0 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 W(R 0 ), but the quartic oscillator's momentum, P 0 , 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 0 • d 12 (R 0 )∆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.

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, ρWe (R, P ), of the full canonical equilibrium density function, ρeq = e −β Ĥ /Z Q , is not stationary under the QCLE; however, ρWe (R, P ) can be written as an expansion in the mass ratio (or h), and it has been shown that the full quantum equilibrium density is conserved under the QCL dynamics up to O(h).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][62][63][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][68][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.
We next apply a symmetric Trotter decomposition to the matrix element of Equation ( 36 In this equation, the symmetric Trotter decomposition separates the subsystem potential in ĥW (R n ) into harmonic, Vho = 1 2 ω 2 R 2 n , and anharmonic, ∆ Ĥ(R n ) = ĥW (R n ) − V ho (R n ), contributions; furthermore, we define Ĥho = Kn + Vho .
The anharmonic term in Equation ( 37) can be approximated as follows: where |n is the subsystem basis and |α; R n is the real-valued adiabatic state with adiabatic energy E α (R n ) with respect to the Hamiltonian, ĥW (R n ).The adjusted energy is The O function in Equation ( 38) reads: where: Rn)   λ e −βE λ (Rn) δ αα − i Now, the canonical partition function is determined by: where Z b is defined by the expression in the second bracket on the second line and Z sn is defined by the expression behind the second bracket on the second line.Z b and Z sn are the bath and subsystem (with its immediate environment) canonical partition functions, respectively.In summary, the time correlation function takes the following simple form: where ρW (X) and Z Q are given by Equations ( 40) and (42), respectively.

Figure 1 .
Figure 1.(a) Potential surface profiles, W α (R 0 ), for the ground adiabatic state (black, dotted), excited adiabatic state (black, dotted) and for the diabatic states, L (green) and R (red).(b) C LL (t) correlation function.These results are associated with the first set of parameters.

Figure 2 .
Figure 2. (a) Potential surface profiles, W α (R 0 ), for the ground adiabatic state (black, dotted), excited adiabatic state (black, dotted) and for the diabatic states, L (green) and R (red).The blue curve is a plot of the un-normalized distribution function, W(R 0 ), Equation (35).(b) C LL (t) correlation functions.(Inset) Short-time C LL (t) computed by the FBTS (blue) and TBSH (red) algorithms.These results are associated with the second set of parameters.