Quantum Dynamics in Open Quantum-Classical Systems

Quantum Dynamics in Open Quantum-Classical Systems

Raymond Kapral Chemical Physics Theory Group, Department of Chemistry, University of Toronto, Toronto, ON, M5S 3H6 Canada

Often quantum systems are not isolated and interactions with their environments must be taken into account. In such open quantum systems these environmental interactions can lead to decoherence and dissipation, which have a marked influence on the properties of the quantum system. In many instances the environment is well-approximated by classical mechanics, so that one is led to consider the dynamics of open quantum-classical systems. Since a full quantum dynamical description of large many-body systems is not currently feasible, mixed quantum-classical methods can provide accurate and computationally tractable ways to follow the dynamics of both the system and its environment. This review focuses on quantum-classical Liouville dynamics, one of several quantum-classical descriptions, and discusses the problems that arise when one attempts to combine quantum and classical mechanics, coherence and decoherence in quantum-classical systems, nonadiabatic dynamics, surface-hopping and mean-field theories and their relation to quantum-classical Liouville dynamics, as well as methods for simulating the dynamics.

I Introduction

It is difficult to follow the dynamics of quantum processes that occur in large and complex systems. Yet, often the quantum phenomena we wish to understand and study take place in such systems. Both naturally-occurring and man-made systems provide examples: excitation energy transfer from light harvesting antenna molecules to the reaction center in photosynthetic bacteria and plants, electronic energy transfer processes in the semiconductor materials used in solar cells, proton transfer processes in some molecular machines that operate in the cell, and the interactions of the qbits in quantum computers with their environment. Although the systems in which these processes take place are complicated and large, it is often the properties that pertain to only a small part of the entire system that are of interest; for example, the electrons or protons that are transferred in a biomolecule. This subsystem of the entire system can then be viewed as an open quantum system that interacts with its environment. In open quantum systems the dynamics of the environment can influence the behavior of the quantum subsystem in significant ways. In particular, it can lead to decoherence and dissipation which can play central roles in the rates and mechanisms of physical processes. This partition of the entire system into two parts has motivated the standard system-bath picture where one of these subsystems (henceforth called the subsystem) is of primary interest while the remainder of the degrees of freedom constitute the environment or bath.

Most system-bath descriptions focus on the dynamics of the subsystem density matrix, which is obtained by tracing over the bath degrees of freedom: . If such a program were carried out fully an exact equation of motion for could be derived and no information about the bath would be lost in this process. Of course, for problems of most interest where the bath is very large with complicated interactions this is not feasible and would defeat the motivation behind the system-bath partition. Consequently, the influence of the bath on the dynamics of the subsystem is embodied in dissipative and other coupling terms in the subsystem evolution equation.

There are many instances where more detailed information about the bath dynamics and its coupling to the subsystem is important. Examples are provided by proton and electron transfer processes in condensed phases or biological systems. As a specific example, consider the proton transfer reaction in a phenol-amine complex, , when the complex is solvated by polar molecules (see Fig. 1). The proton transfer events are strongly correlated with local solvent collective polarization changes. Subtle changes in the orientations of neighboring solvent molecules can induce proton transfers within the complex, which, in turn, influence the polarization of the solvent. The treatment of the dynamics in such cases requires detailed information about the dynamics of the environment and its coupling to the quantum process. It is difficult to capture such subtle effects without fully accounting for dynamics of individual solvent molecules in the bath.

Figure 1: Schematic representation showing the local solvent structure around the phenol-triethylamine complex. The covalent form of the phenol-amine complex (left) is unfavorably solvated by the polar solvent molecules. This induces a proton transfer giving rise to the ionic form (right). Subsequent solvent dynamics can lead to solvent polarization that favors the covalent form and the reverse proton transfer.

When investigating the dynamics of a quantum system it is often useful and appropriate to take into account the characteristics of the different degrees of freedom that comprise the system. The fact that electronic and nuclear motions occur on very different time scales, as a result of the disparity in their masses, forms the basis for the Born-Oppenheimer approximation where the nuclear-configuration-dependent electronic energy is used as the potential energy for the evolution of the nuclear degrees of freedom. This distinction between electronic and nuclear degrees of freedom is an example of the more general partition of a quantum system into subsystems with different characteristics.

Since the scale separation in the Born-Oppenheimer approximation is approximate, it can break down and its breakdown leads to coupling among many electronic energy surfaces. When this occurs, the evolution is no longer described by adiabatic dynamics on a single potential energy surface and nonadiabatic effects become important. Nonadiabatic dynamics plays an essential role in the description of many physical phenomena, such as photochemical processes where transitions among various electronic states occur as a result of avoided crossings of adiabatic states or conical intersections between potential energy surfaces.

In the examples presented above the molecules comprising the bath are often much more massive than those in the subsystem (). This fact motivates the construction of a quantum-classical description where the bath, in the absence of interactions with the quantum subsystem, is described by classical mechanics. Mixed quantum-classical methods provide a means to investigate quantum dynamics in large complex systems, since fully quantum treatments of the dynamics of such systems are not feasible. The study of such open quantum-classical systems is the main topic of this review. Since quantum and classical mechanics do not easily mix, one must consider the properties of schemes that combine these two types of mechanics. One such scheme, quantum-classical Liouville dynamics, will be discussed in detail and its features will be compared to other quantum-classical and full quantum methods.

Ii Open Quantum Systems

Since the quantum systems we study are rarely isolated and interact with the environments within which they reside, the investigation of the dynamics of such open quantum systems is a worthy endeavor. The full description of the time evolution of a composite quantum system comprising a subsystem and bath is given by the quantum Liouville equation,


where is the density matrix at time , is the total Hamiltonian, and the square brackets denote the commutator.

Introducing some of the notation that will be used in this paper, we denote by the coordinate operators for the subsystem degrees of freedom with mass , while the remaining bath degrees of freedom with mass have coordinate operators . (The formalism is easily generalized to situations where the masses and depend on the particle index.) The total Hamiltonian takes the form


where the momentum operators for the subsystem and bath are and , respectively. It is convenient to write the potential energy operator, as a sum of subsystem, bath and coupling contributions: . In this case the Hamiltonian operator can be written as a sum of contributions,


where is the quantum subsystem Hamiltonian, is the quantum bath Hamiltonian and is the coupling between these two subsystems.

Most often in considering the dynamics of such open quantum systems one traces over the bath since it is the dynamics of the subsystem that is of interest. As noted in the Introduction, a considerable research effort has been devoted to the construction of equations of motion for the reduced density matrix, . The Redfield equation Redfield (1965) describes the dynamics of a subsystem weakly coupled to a bath with suitably fast bath relaxation time scales, since a Born-Markov approximation is made in its derivation. In a basis of eigenstates of , , it has the form,


where the summation convention has been used. This convention will be used throughout the paper when confusion is unlikely. Here , while the second term on the right accounts for dissipative effects due to the bath. Remaining within the Born-Markov approximation, the general form of the equation of motion for a reduced density matrix that guarantees its positivity is given by the Lindblad equation Lindblad (1976),

where the are operators that account for interactions with the bath. In addition to these equations, a number of other expressions for the evolution of the reduced density matrix have been derived. These include various master equations and generalized quantum master equations. There is a large literature dealing with open quantum systems, which is described and surveyed in books on this topic. Davies (1976); Weiss (1999); Breuer and Petruccione (2006) In such reduced descriptions information about the bath is contained in parameters that enter in the operators that describe the coupling between the subsystem and bath. Also, quantum-classical versions of the Redfield Toutounji (2005) and Lindblad Toutounji and Kapral (2001) equations have been derived.

There are many applications, such as those mentioned in the Introduction, where a more detailed treatment of the bath dynamics and its interactions with the subsystem is required, even though one’s primary interest is in the dynamics of the subsystem. If, as we suppose here, the systems we wish to study are large and may involve complex molecular constituents, a full quantum mechanical treatment is beyond the scope of existing computational power and algorithms. Currently, the only viable way to simulate the dynamics of such systems is by using mixed quantum-classical schemes. Quantum-classical methods in a variety of forms and derived in a variety of ways have been used to simulate the dynamics. Herman (1994); Tully (1998); Ben-Nun and Martińez (1998); Kapral (2006); Tully (2012) Mean field and surface-hopping methods are widely employed and will be discussed in some detail below. Mixed quantum-classical dynamics Agostini et al. (2013) based the on the exact time-dependent potential energy surfaces derived from the exact decomposition of electronic and nuclear motions Abedi, Maitra, and Gross (2010) has been constructed. In addition, semiclassical path integral formulations of quantum mechanics Sun and Miller (1997); Sun, Wang, and Miller (1998); Makri and Thompson (1998); Miller (2009); Lambert and Makri (2012) and ring polymer dynamics methods Habershon et al. (2013) have been developed to approximate the dynamics of open quantum systems.

In the next section the specific version of mixed quantum-classical dynamics that is the subject of this review, quantum-classical Liouville dynamics, will be described. The passage from quantum to classical dynamics is itself a difficult problem with an extensive literature, and decoherence is often invoked to effect this passage. Joos et al. (2003); Zurek (1991) Considerations based on decoherence can also be used motivate the use of mixed quantum-classical descriptions. Shiokawa and Kapral (2002) Mean-field and surface-hopping methods suffer from difficulties related to the treatment of coherence and decoherence, and these methods will be discussed in the context of the quantum-classical Liouville equation, which is derived and discussed in the next two sections. Some applications of the theory to specific systems will be presented in order to test the accuracy of this equation description and the algorithms used to simulate its dynamics.

Iii Quantum-Classical Liouville Dynamics

The first step in constructing a quantum-classical Liouville description is to introduce a phase space representation of the bath degrees of freedom in preparation for the passage to the classical bath limit. This is conveniently done by introducing a partial Wigner transform Wigner (1932) over the bath degrees of freedom defined by


with an analogous expression for the partial Wigner transform of an operator in which the prefactor is absent. We let to simplify the notation. The quantum Liouville equation then takes the form,

To obtain this equation the formula for the Wigner transform of a product of operators Imre et al. (1967),


was used. Here the operator , where the arrows denote the directions in which the derivatives act, is the negative of the Poisson bracket operator,


The partial Wigner transform of the total Hamiltonian is,


We have dropped the subscript W on the potential energy operator to simplify the notation; when the argument contains the partial Wigner transform is implied.

Derivation of the QCLE

The quantum-classical Liouville equation (QCLE) can be derived by formally expanding the exponential operators on the right side of Eq. (III) to Aleksandrov (1981); Gerasimenko (1982) The truncation of the series expansion can be justified for systems where the masses of particles in the environment are much greater than those of the subsystem, Kapral and Ciccotti (1999) Scaling similar to that in the microscopic derivation of the Langevin equation for Brownian motion from the classical Liouville equation Mazur and Oppenheim (1970) can be used for this purpose, and we may write the equations in terms of the reduced bath momenta, where . In this variable the kinetic energies of the light and heavy particle systems are comparable so that is of order . To see this more explicitly we introduce scaled units where energy is expressed in the unit , time in and length in units of . Using these length and time units, the scaling factor for the momentum is . Thus, in terms of the scaled variables and we have


where the prime on indicates that it is expressed in the primed variables. Note that for a system characterized by a temperature the small parameter can be written as the ratio of the thermal de Broglie wavelengths and of the heavy bath and light subsystem particles, respectively, , and truncation of the dynamics to terms of effectively averages out the quantum bath oscillations on the longer quantum length scale of the light subsystem.

Inserting the expression for the exponential Poisson bracket operator, valid to , into the scaled version of Eq. (III) and returning to unscaled variables we obtain the quantum-classical Liouville equation Kapral and Ciccotti (1999),


Additional discussion of this equation can be found in the literature Kapral (2006); Kapral and Ciccotti (1999); Aleksandrov (1981); Gerasimenko (1982); Boucher and Traschen (1988); Zhang and Balescu (1988); Donoso and Martens (1998); Horenko et al. (2002); Shi and Geva (2004); Thorndyke and Micha (2005); Bousquet et al. (2011). Comparison of the second and third equalities in this equation defines the QCL operator , and given this definition the formal solution of the QCLE is


where we let here and in the following to simplify the notation. The QCLE (12) may also be written in the form Nielsen, Kapral, and Ciccotti (2001),


which resembles the quantum Liouville equation but the quantum Hamiltonian operator is replaced by the forward and backward operators,


This form of the evolution equation has been used to discuss the statistical mechanical properties of QCL dynamics Nielsen, Kapral, and Ciccotti (2001), and will be used later to derive approximate solutions to the QCLE.

In applications it is often more convenient to evolve an operator rather than the density matrix and we may easily write the evolution equations for operators. Starting from the Heisenberg equation of motion for an operator ,


one can carry out an analogous calculation to find the QCLE for the partial Wigner transform of this operator:


whose formal solution can be written as


QCLE from linearization

The QCLE, when expressed in the adiabatic or subsystem bases, has been derived from linearization of the path integral expression for the density matrix by Shi and Geva Shi and Geva (2004). It can also be derived in a basis-free form by linearization Bonella, Ciccotti, and Kapral. (2010) and it is instructive to sketch this derivation here to see how the QCLE can be obtained from a perspective that differs from that discussed in the previous subsection.

The time evolution of the quantum density operator from time to a short later time is given by


Writing the Hamiltonian in the form , for this short time interval, a Trotter factorization of the propagators can be made:


For simplicity, we have suppressed the dependence in but kept the dependence since it is required in the derivation. Working in the representation for the bath, inserting resolutions of the identity, and evaluating the contributions coming from the kinetic energy operators that appear in the resulting expression, we obtain


Next, we make the change of variables and , along with similar variable changes for the momenta and . In the new variables, the density matrix element is


We may now make use of the definition of the partial Wigner transform (see Eq. (6)) in the expression for the matrix element of the density operator in Eq. (22) to derive an equation of motion for . To do this we first expand the exponentials that depend on to first order in this parameter; e.g., . We may then use this expansion to compute the finite difference expression . Finally we multiply the equation by , integrate the result over and take the limit . The result of these operations is


where . This integro-differential equation describes the full quantum evolution of the density matrix element; however, it is not a closed equation for because of the dependence of on . If we make use of the expansion of this operator to linear order in , when performing the integrals in the right side of Eq. (23), we obtain the QCLE in Eq. (12). The linearization approximation can be justified for systems where Bonella, Ciccotti, and Kapral. (2010) The same scaled variables introduced above in the first derivation may also be used to re-express Eq. (23) in scaled form. In this scaled form one may show that the expansion in is equivalent to an expansion in the mass ratio parameter .

QCLE in a dissipative environment

At times it may be convenient to further partition the bath into two subsets of degrees of freedom, , where the variables are directly coupled to the quantum subsystem and the remainder of the (usually large number of) degrees of freedom denoted by only participate in the subsystem dynamics indirectly through their coupling to . In such a case we can project these degrees of freedom out of the QCLE to derive a dissipative evolution equation for the quantum subsystem and the directly coupled variables Kapral (2001). For example, such a description could be useful in studies of proton or electron transfer in biomolecules where remote portions of the biomolecule and solvent need not be treated in detail but, nevertheless, these remote degrees of freedom do provide a source of decoherence and dissipation on the relevant degrees of freedom.

For a system of this type the partially Wigner transformed total Hamiltonian of the system is,


The potential energy operator, , includes all of the coupling contributions discussed above, namely, the potential energy operator for the quantum subsystem and directly coupled degrees of freedom, the potential energy of the outer bath and the coupling between the two bath subsystems, .

An evolution equation for the reduced density matrix of the quantum subsystem and directly coupled degrees of freedom,


can be obtained by using projection operator methods. Nakajima (1958); Zwanzig (1961) The result of this calculation is a dissipative QCLE, which takes the form Kapral (2001),


where is the QCL operator introduced earlier but now only for the quantum subsystem and bath degrees of freedom. The effects of the less relevant bath degrees of freedom are accounted for by the mean force defined by , where the average is over a canonical equilibrium distribution involving the Hamiltonian . The Fokker-Planck-like operator in Eq. (26) depends on the fixed particle friction tensor, , defined by


where , and its time evolution is given by the classical dynamics of the degrees of freedom in the field of the fixed coordinates. The quantum-classical limit of the multi-state Fokker-Planck equation introduced by Tanimura and Mukamel Tanimura and Mukamel (1994) is similar to the dissipative QCLE (26) when expressed in the subsystem basis.

Iv Some properties of the QCLE

The QCLE specifies the time evolution the density matrix of the entire system comprising the subsystem and bath and conserves the energy of the system. If the coupling potential in the Hamiltonian is zero, the density matrix factors into a product of subsystem and bath density matrices, . In this limit the subsystem density matrix satisfies the quantum Liouville equation,


and bath phase space density satisfies the classical Liouville equation,


While the bath evolves by classical mechanics when it is not coupled to the quantum subsystem, its evolution is no longer classical when coupling is present. As we shall see in more detail below, not only does the bath serve to account for the effects of decoherence and dissipation in the subsystem, it is also responsible for the creation of coherence. Conversely, the subsystem can interact with the bath to modify its dynamics. This leads to a very complicated evolution, but one which incorporates many of the features that are essential for the description of physical systems.

Often, when considering the dynamics of a quantum system coupled to a bath, the bath is modeled by a collection of harmonic oscillators which are bilinearly coupled to the quantum subsystem. In this case we may write the coupling potential as . The partially Wigner transformed Hamiltonian then takes the form,


where is the harmonic oscillator bath Hamiltonian. When the Hamiltonian has this form one may show easily that . Consequently, when the exponential Poisson bracket operators in Eq. (III) are expanded in a power series, the series truncates at linear order and we obtain the QCLE in the form given in Eq. (14); thus, the QCLE is exact for general quantum subsystems which are bilinearly coupled to harmonic baths. For more general Hamiltonian operators the series does not truncate and QCL dynamics is an approximation to full quantum dynamics.

Quantum and classical mechanics do not like to mix. The coupling between the smooth classical phase space evolution of the bath and the quantum subsystem dynamics with quantum fluctuations on small scales presents challenges for any quantum-classical description. The QCLE, being an approximation to full quantum dynamics, is not without defects. One of its features that requires consideration is its lack of a Lie algebraic structure. The quantum commutator bracket and Poisson bracket for quantum and classical mechanics, respectively, are bilinear, skew symmetric, and satisfy the Jacobi identity, so that these brackets have Lie algebraic structures. The quantum-classical bracket, , which is the combination of the commutator and the Poisson bracket terms does not have such a Lie algebraic structure. While this bracket is bilinear and skew symmetric, it does not exactly satisfy the Jacobi identity. Instead, the Jacobi identity is satisfied only to order (or if scaled variables are considered): . The lack of a Lie algebraic structure, its implications for the dynamics, and the construction of the statistical mechanics of quantum-classical systems were discussed earlier Nielsen, Kapral, and Ciccotti (2001); Kapral and Ciccotti (2002) where full details may be found. For example, the standard linear response derivations of quantum transport properties have to be modified, and in quantum-classical dynamics the evolution of a product of operators is not the product of the evolved operators; this is true only to order . These feature are not unique to QCL dynamics and almost all mixed quantum-classical methods used in simulations suffer from such defects, although they are rarely discussed. Mixed quantum-classical dynamics and its algebraic structure continue to attract the attention of researchers. Salcedo (1996, 2012); Prezhdo and Kisil (1997); Sergi (2005); Prezhdo (2006); Salcedo (2007); F. Agostini and Ciccotti (2007); Hall and Reginatto (2005); Hall (2008)

One way to bypass some of the difficulties in the formulation of the statistical mechanics of quantum-classical systems that are associated with a lack of a Lie algebraic structure is to derive expressions for average values and transport property using full quantum statistical mechanics. Then, starting with these exact quantum expressions, one may approximate the quantum dynamics by quantum-classical dynamics. Sergi and Kapral (2004); Kim and Kapral (2005a, b); Hsieh and Kapral (2014) In this framework the expectation value of an observable is given by,


where is the partial Wigner transform of the initial quantum density operator and the evolution of is given by the QCLE. Similarly, the expressions for transport coefficients involve time integrals of correlation functions of the form,


where is the partial Wigner transform of the product of the quantum canonical density operator and the operator , and the time evolution of is again given by the QCLE. Such formulations preserve the full quantum equilibrium structure which, while difficult to compute, is computationally much more tractable than full quantum dynamics. Poulsen, Nyman, and Rossky (2003); Ananth and Miller (2010) The importance of quantum versus classical equilibrium sampling on reactive-flux correlation functions, whose time integrals are reaction rate coefficients, has been investigated in the context of quantum-classical Liouville dynamics. Kim and Kapral (2005b) In this review we shall focus on dynamics but, when applications are considered, the above equations that contain the quantum initial or equilibrium density matrices will be used.

V Surface hopping, coherence and decoherence

Surface-hopping methods are commonly used to simulate the nonadiabatic dynamics of quantum-classical systems. In such schemes the bath phase space variables follow Newtonian trajectories on single adiabatic surfaces. Nonadiabatic effects are taken into account by hops between different adiabatic surfaces that are governed by probabilistic rules.

One of the most widely used schemes is Tully’s fewest-switches surface hopping. Tully (1990, 1991, 1998) In this method one assumes that the electronic wave function depends on the time-dependent nuclear positions , whose evolution is governed by a stochastic algorithm. More specifically, choosing to work in a basis of the instantaneous adiabatic eigenfunctions of the Hamiltonian , , we may expand the wave function as . An expression for the time evolution of the subsystem density matrix can be obtained by substitution into the Schrödinger equation. The equations of motion for its matrix elements, are given by,


In this equation is the nonadiabatic coupling matrix element, . From this expression the rate of change of the population in state may be written as


where, for simplicity, we have suppressed the time dependence in the variables, and stands for the real part. This rate has contributions from transitions to and from all other states . Consider a single specific state . Then transitions into from and out of to will determine the rate of change of the population due to transitions involving this state. In fewest-switches surface hopping the transitions are dropped and the transition rate for , , is adjusted to give the correct weighting of populations:


This transition rate is used to construct surface-hopping trajectories that specify the evolution of the phase space variables as follows: When the system is in state , the coordinates evolve by Newtonian trajectories on the adiabatic surface. Transitions to other states occur with probabilities per unit time, . Since the rates may take negative values, the Heaviside function sets the probability zero for negative values of the rate. If the transition to state occurs, the momentum of the system is adjusted to conserve energy and the system then propagates on the adiabatic surface. The momentum adjustment is taken to occur along the direction of the nonadiabatic coupling vector and is given by , with


The form that the stochastic evolution takes can be seen from an examination of Fig. 2, which schematically shows the evolution of a wave packet that starts on the upper adiabatic surface of a two level system with a simple avoided crossing. (This is Tully’s simple avoided crossing model. Tully (1990)) When the system enters the region of strong nonadiabatic coupling near the avoided crossing, nonadiabatic transitions to the lower state are likely, a surface hop occurs and the system then continues to evolve on the lower surface after momentum adjustment.

Figure 2: Schematic representation of the evolution of a wave packet in a two-level system with a simple avoided crossing. The diabatic (crossing curves) and adiabatic (avoided crossing curves) are shown. Following the nonadiabatic transition from the upper to lower adiabatic surfaces, the system continues to evolve on the lower surface until the next nonadiabatic transition.

For upward transitions it may happen that there is insufficient energy in the environment to insure energy conservation. In this case the transition rule needs to be modified, usually by setting the transition probability to zero. This scheme is very easy to simulate and captures much of the essential physics of the nonadiabatic dynamics.

Fewest-switches surface hopping does suffer from some defects associated with the fact that decoherence is not properly treated. The transition probability depends on the off-diagonal elements of the density matrix but no mechanism for their decay is included in the model. As a result, the fewest-switches surface hopping model overestimates coherence effects and retains memory which can influence the probabilities of subsequent hops. Several methods have been proposed to incorporate the effects of decoherence in mixed quantum-classical theories and, in particular, in surface-hopping schemes. Neria and Nitzan (1993); Hammes-Schiffer and Tully (1994); Bittner and Rossky (1995); Bittner, Schwartz, and Rossky (1997); Schwartz et al. (1996); Bedard-Hearn, Larsen, and Schwartz (2005); Subotnik and Shenvi (2011); Shenvi, Subotnik, and Yang (2011); Landry, J.Falk, and Subotnik (2013); Subotnik, Ouyang, and Landry (2013); Subotnik (2011); Jaeger, Fischer, and Prezhdo (2012) In many of these methods a term of the form, , is appended to the equation of motion for the off-diagonal elements of the subsystem density matrix to account for the decay of coherence. The decoherence rate is estimated using perturbation theory or from physical considerations involving the overlap of nuclear wave functions. In the remainder of this section we discuss how the QCLE accounts for decoherence and comment on its links to surface-hopping methods.

QCL dynamics in the adiabatic basis and decoherence

Since surface-hopping methods are often formulated in the adiabatic basis, it is instructive to discuss the dynamical picture that emerges when the QCLE is expressed in this basis. Adopting an Eulerian description, the adiabatic energies, , and the adiabatic states, , depend parametrically on the coordinates of the bath. We may then take matrix elements of Eq. (12),


to find an evolution equation for the density matrix elements, . Evaluation of the matrix elements on the right side of this equation yields an expression for the QCL superoperator Kapral and Ciccotti (1999),


Here the frequency (now in the adiabatic basis), and is the classical Liouville operator


and involves the Hellmann-Feynman forces, . The superoperator, , whose matrix elements are


couples the dynamics on the individual and mean adiabatic surfaces so that the evolution is no longer described by Newtonian dynamics.

The resulting QCLE in the adiabatic representation reads,


To simplify we shall often use a formal notation and write Eq. (41) as


where and (without “hats”) are understood to be a matrix and superoperator, respectively, in the adiabatic basis.

Insight into the nature of QCL dynamics can be obtained as follows. If the operator is dropped the resulting equation of motion for the diagonal elements of the density matrix is


which implies that the phase space density is constant along trajectories on the adiabatic surface,




with the notation . The off-diagonal density matrix elements satisfy


whose solution is


where and the evolution of the phase space coordinates of the bath is given by


The off-diagonal elements accumulate a phase in the course of their evolution on the mean of the two and adiabatic surfaces.

The momentum derivative terms in are responsible for the energy transfers that occur to and from the bath when the subsystem density matrix changes its quantum state. Consequently the subsystem and bath interact with each other and the dynamics of both the subsystem and bath are modified in the course of the evolution. Further, we can see from the structure of the QCLE that there are continuous changes to the subsystem quantum state and bath momenta during the evolution, as opposed to the jumps that appear in surface-hopping schemes. Nonetheless, links to surface-hopping methods can be made.

Subotnik, Ouyang and Landry Subotnik, Ouyang, and Landry (2013) established a connection between fewest-switches surface-hopping and the QCLE. They investigated what must be done to the equations describing fewest-switches surface hopping in order to obtain the QCL dynamics. Since there are continuous bath momentum changes in QCL dynamics and discontinuous changes in fewest-switches surface hopping, there are limitations on the nuclear momenta. An important element in their analysis is the fact that terms of the form, , that account for decoherence must be added to the fewest-switches approach. The specific form of the decoherence rate in their analysis is


The superscript indicates that evolution is on the adiabatic surface and all quantities on the right are taken to evolve on this surface. An analogous expression can be written for .

Recall that surface-hopping schemes assume that the dynamics occurs on single adiabatic surfaces between hops. Given this fact, we can understand the need for such a term by viewing QCL dynamics in a frame of reference corresponding to motion along single adiabatic surfaces. To see this consider the equation of motion for an off-diagonal element of the density matrix as given by the QCLE. From Eqs. (38)-(41) we have


Defining the material derivative for the flow on the adiabatic surface as


we obtain


We see that the second term on the right side of this equation is just the decoherence factor that appears in Eq. (49). The fact that decoherence depends on the difference between the forces is a common factor in many of the models for decoherence mentioned above. The decoherence contribution is difficult to compute in its current form because of the bath momentum derivative and it is usually approximated in applications. Subotnik, Ouyang, and Landry (2013)

Surface-hopping solution of the QCLE

As discussed above, the dynamics prescribed by the QCLE is not in the form of surface hopping since quantum state and bath momentum changes as embodied in the superoperator occur continuously throughout the evolution. The effects of can be seen by considering the formal solution of Eq. (41),