Mapping quantum-classical Liouville equation: projectors and trajectories

# Mapping quantum-classical Liouville equation: projectors and trajectories

Aaron Kelly Chemical Physics Theory Group, Department of Chemistry, University of Toronto, Toronto, ON M5S 3H6, Canada Department of Chemistry, Stanford University, 333 Campus Drive, Stanford, CA 94305, USA    Ramses van Zon Chemical Physics Theory Group, Department of Chemistry, University of Toronto, Toronto, ON M5S 3H6, Canada    Jeremy Schofield Chemical Physics Theory Group, Department of Chemistry, University of Toronto, Toronto, ON M5S 3H6, Canada    Raymond Kapral Chemical Physics Theory Group, Department of Chemistry, University of Toronto, Toronto, ON M5S 3H6, Canada
July 5, 2019
###### Abstract

The evolution of a mixed quantum-classical system is expressed in the mapping formalism where discrete quantum states are mapped onto oscillator states, resulting in a phase space description of the quantum degrees of freedom. By defining projection operators onto the mapping states corresponding to the physical quantum states, it is shown that the mapping quantum-classical Liouville operator commutes with the projection operator so that the dynamics is confined to the physical space. It is also shown that a trajectory-based solution of this equation can be constructed that requires the simulation of an ensemble of entangled trajectories. An approximation to this evolution equation which retains only the Poisson bracket contribution to the evolution operator does admit a solution in an ensemble of independent trajectories but it is shown that this operator does not commute with the projection operators and the dynamics may take the system outside the physical space. The dynamical instabilities, utility and domain of validity of this approximate dynamics are discussed. The effects are illustrated by simulations on several quantum systems.

## I Introduction

Since phenomena such as electron and proton transfer dynamics Bell (1973); Kornyshev et al. (1997), excited state relaxation processes Schryver et al. (2001) and energy transport in light harvesting systems Cheng and Fleming (2009); Scholes (2010) are quantum in nature, the development of theoretical descriptions and simulation methods for quantum many-body systems is a central topic of research. Although various techniques can be used to study such problems, quantum-classical methods Herman (1994); Tully (1998); Billing (2003); Jasper et al. (1981); kap (); Subotnik and Shenvi (2011), where certain degrees of freedom are singled out for a full quantum treatment while other environmental variables are treated classically, permit one to investigate large and complex systems that cannot be studied by other means.

In this article we consider descriptions of the dynamics based on the quantum-classical Liouville equation kap () (QCLE) and, in particular, its representation in the mapping basis Kim et al. (2008); Nassimi and Kapral. (2009); Nassimi et al. (2010). The mapping formalism provides an exact mapping of discrete quantum states onto continuous variables Stock and Thoss. (2005) and in quantum-classical systems leads to phase-space-like evolution equations for both quantum and classical degrees of freedom. The mapping basis has been used in a number of different quantum-classical formulations, often based on semi-classical path integral expressions for the dynamics Miller and McCurdy (1978); Meyer and Miller (1979); Sun et al. (1998); Miller (2001); Ananth et al. (2007); Stock and Thoss (1997); Muller and Stock (1998); Thoss and Stock (1999); Thoss and Wang (2004); Stock and Thoss. (2005); Bonella and Coker (2001, 2003, 2005); Dunkel et al. (2008). The representation of the quantum-classical Liouville equation in the mapping basis leads to an equation of motion whose Liouvillian consists of a Poisson bracket term in the full quantum subsystem-classical bath phase space, and a more complex term involving second derivatives of quantum phase space variables and first derivatives with respect to bath momenta Kim et al. (2008). This latter contribution has been shown to be an excess coupling term related to a portion of the back reaction of the quantum subsystem on the bath Nassimi et al. (2010).

Various aspects of the QCLE in the mapping basis and properties of its full and approximate solutions are discussed in this paper. The solutions of the quantum-classical Liouville equation cannot be obtained from the dynamics of an ensemble of independent classical-like trajectories Kapral and Ciccotti (1999). In the adiabatic basis this equation admits a solution in terms of surface-hopping trajectories Kapral and Ciccotti (1999); MacKernan et al. (2008); Sergi et al. (2003), but other schemes have been used to simulate the dynamics Martens and Fang (1996); Donoso and Martens (1998); Wan and Schofield (2000, 2002). When it is expressed in the mapping basis, we show that a solution can be obtained in terms of an ensemble of entangled trajectories. The excess coupling gives rise to correlations between the dynamics of the quantum mapping degrees of freedom and the bath phase space variables that are responsible for the entanglement of the trajectories in the ensemble. The derivation of the entangled trajectory picture is similar to that for trajectory solutions of the Wigner-Liouville equation Donoso and Martens (2001); Donoso et al. (2003).

If the excess coupling term is dropped and only the Poisson bracket part of the Liouvillian is retained, a very simple equation of motion that admits a solution in terms of characteristics is obtained. Consequently, its solutions can be obtained from simulations of an ensemble of independent trajectories evolving under Newtonian dynamics. The set of ordinary differential equations has appeared earlier in mapping formulations based on semi-classical path integral formulations of the dynamics Stock and Thoss. (2005); Miller (2001); Ananth et al. (2007), indicating a close connection between this approximation to the quantum-classical Liouville equation and those formulations. The solutions of this Poisson bracket approximation to the QCLE, as well as those of other semi-classical schemes that use this set of evolution equations, often provide a quantitatively accurate description of the dynamics Stock and Thoss. (2005); Kim et al. (2008); Nassimi et al. (2010). However for some systems the solutions are not without artifacts and difficulties. Some of these difficulties can be traced to the fact that the independent-ensemble dynamics can take the system out of the physical space and inverted potentials can appear in the evolution equations, which may lead to instabilities Stock and Thoss. (2005); Bonella and Coker (2001, 2005).

The main results of this paper are as follows: We present derivations of expressions for mapping quantum-classical (MQCL) evolution equations and expectation values of operators that explicitly show how projection operators onto the physical mapping eigenstates enter the formulation. We demonstrate that the MQCL operator commutes with this projection operator so that dynamics under this evolution is confined to the physical space. This full quantum-classical dynamics in the mapping basis can be simulated by an ensemble of entangled trajectories. We also show that when the excess coupling term is neglected the resulting Poisson bracket operator no longer commutes with the projection operator so that this approximate dynamics can take the system out of the physical space. Given this context, we revisit the issue of instabilities in the dynamics of the Poisson bracket approximation and discuss the conditions under which such instabilities are likely to arise and lead to inaccuracies in the solutions.

In Sec. II we outline the representation of the quantum-classical Liouville equation in the mapping basis and show how average values of time dependent observables may be computed. We also define a projection operator onto the mapping states and show how this projector enters the expressions for the expectation values and evolution equations. Section III briefly describes the entangled trajectory solution to the QCLE in the mapping basis. This section also shows that when the excess coupling term is neglected, a solution in terms of an ensemble of independent trajectories is possible. In Sec. IV the approximate evolution equation obtained by retaining only the Poisson bracket term in the Liouville operator is considered and the dynamical instabilities that can arise in the course of the evolution are highlighted. Various aspects of the theoretical analysis that concern the approximate solutions and resulting instabilities are illustrated by simulations of a number of model systems. A brief summary of the main results of the study, along with comments, are given in Sec. V. The Appendices provide material to support the text. In particular, we describe an efficient simulation algorithm for the ordinary differential equations that underlie the solutions of the Poisson bracket approximation to the QCLE.

## Ii Quantum-Classical Liouville Equation: Mapping, Projectors and Expectation Values

The quantum-classical Liouville equation (QCLE),

 ∂∂t^ρW(X,t)=−i^L^ρW(X,t), (1)

describes the time evolution of the density matrix , which is a quantum operator that depends on the classical phase space variables of the environment. The quantum-classical Liouville operator is defined by

 i^L⋅=iℏ[^HW,⋅]−12({^HW,⋅}−{⋅,^HW}), (2)

where is the partial Wigner transform of the total Hamiltonian of the system, is the commutator and is the Poisson bracket in the phase space of the classical variables . The total Hamiltonian may be written as the sum of environmental (bath), subsystem and coupling terms, , where is the bath Hamiltonian with the bath potential energy, is the subsystem Hamiltonian with and the subsystem momentum and potential energy operators, and is the coupling potential energy operator. Here and are the masses of the subsystem and bath particles, respectively.

The QCLE may be written in the basis, , that spans the quantum subsystem space with eigenfunctions defined by the eigenvalue problem, . Taking matrix elements of Eq. (1) we obtain

 ∂∂tρλλ′W(X,t)=−iLλλ′,νν′ρνν′W(X,t). (3)

The Einstein summation convention is used here and in subsequent equations although, on occasion, sums will be explicitly written for purposes of clarity. The QCL operator in the subsystem basis is Kapral and Ciccotti (1999)

 iLλλ′,νν′=iωλλ′δλνδλ′ν′−iℏ(δλνVν′λ′c−Vλνcδλ′ν′) +(PM⋅∂∂R+Fe(R)⋅∂∂P)δλνδλ′ν′ −12(δλ′ν′∂Vλνc∂R+δλν∂Vν′λ′c∂R)⋅∂∂P, (4)

where and is the force due to molecules in the environment.

The evolution equation for an observable , analogous to Eq. (1), is

 ddt^BW(X,t)=i^L^BW(X,t), (5)

and its representation in the subsystem basis is analogous to Eq. (3) with a change in sign on the right side.

### ii.1 Representation in Mapping Basis and Projection Operators

In the mapping basis Schwinger (1965); Stock and Thoss. (2005) the eigenfunctions of an -state quantum subsystem can be replaced with eigenfunctions of fictitious harmonic oscillators, , having occupation numbers which are limited to 0 or 1: . Creation and annihilation operators on these states, and , respectively, may be defined. For any operator whose matrix elements in the subsystem basis are , we may associate a mapping basis operator , where

 ^Bm(X)=Bλλ′W(X)^a†λ^aλ′. (6)

It is then evident that the matrix element .

The expression for may also be written in terms of the Wigner transforms in the space of the mapping variables. Inserting complete sets of coordinate states , and making the usual coordinate transformations appropriate for Wigner transforms, , we obtain

 Bλλ′W(X)=⟨mλ|^Bm(X)|mλ′⟩= (7) ∫drdz⟨mλ|r−z2⟩⟨r−z2|^Bm(X)|r+z2⟩⟨r+z2|mλ′⟩

Another form for the matrix element can be obtained by inserting the Wigner transform of an operator and its inverse as

 ⟨r−z2|^Bm(X)|r+z2⟩=1(2πℏ)N∫dpe−ip⋅z/ℏBm(X), Bm(X)=∫dzeip⋅z/ℏ⟨r−z2|^Bm(X)|r+z2⟩. (8)

Here are the extended phase space coordinates for the subsystem mapping variables, , and the environment, . Making these substitutions in Eq. (7) we obtain,

 Bλλ′W(X)=∫dxBm(X)gλλ′(x), (9)

where we have defined foo (a)

 gλλ′(x)=1(2πℏ)N∫dze−ip⋅z/ℏ⟨r+z2|mλ′⟩⟨mλ|r−z2⟩ =1(2πℏ)N∫dzeip⋅z/ℏ⟨r−z2|mλ′⟩⟨mλ|r+z2⟩. (10)

Evaluating the integral we obtain an explicit expression for :

 gλλ′(x)=ϕ(x) (11) ×2ℏ[rλrλ′+pλpλ′−i(rλpλ′−rλ′pλ)−ℏ2δλλ′],

where is a normalized Gaussian function. Here in the Einstein summation convention.

The expression for in Eq. (II.1) can be simplified by evaluating the integral in the Wigner transform. Using the definition of in Eq. (6), Eq. (II.1) may be written as

 Bm(X)=Bλλ′W(X)∫dzeip⋅z/ℏ⟨r−z2|^a†λ^aλ′|r+z2⟩. (12)

Noting that the factor multiplying is the Wigner transform of , , whose explicit value is

 cλλ′(x)=12ℏ[rλrλ′+pλpλ′+i(rλpλ′−rλ′pλ)−ℏδλλ′],

we find

 Bm(X)=Bλλ′W(X)cλλ′(x). (14)

We may deduce a number of other relations given the definitions stated above. A mapping operator acts on mapping functions . In this space we have the completeness relations , where is projector onto the complete set of mapping states. foo (b) Thus, a mapping operator can be written using this projector as

 ^BPm(X)=^P^Bm(X)^P = |mλ⟩⟨mλ|^Bm(X)|mλ′⟩⟨mλ′| (15) = |mλ⟩Bλλ′W(X)⟨mλ′|,

where in the second line we used the equivalence between matrix elements in the subsystem and mapping representations given in Eq. (7). We can make use of the Wigner transforms defined in Eq. (II.1) to write these relations in other forms. Using the first equality in Eq. (15) we have

 BPm(X)=∫dzeip⋅z/ℏ⟨r−z2|^BPm(X)|r+z2⟩= ∫dzeip⋅z/ℏ⟨r−z2|mλ⟩⟨mλ|^Bm(X)|mλ′⟩⟨mλ′|r+z2⟩ =(2πℏ)Ngλ′λ(x)∫dx′gλλ′(x′)Bm(x′,X), ≡PBm(X), (16)

where we used Eqs. (7) and (9). The last line defines the projection operator that projects any function of the mapping phase space coordinates, , onto the mapping states,

 Pf(x)=(2πℏ)Ngλ′λ(x)∫dx′gλλ′(x′)f(x′). (17)

One may verify that since

 (2πℏ)N∫dxgλλ′(x)gν′ν(x)=δλνδλ′ν′. (18)

An equivalent expression for can be obtained by starting with the last equality in Eq. (15) and taking Wigner transforms to find

 BPm(X)=∫dzeip⋅z/ℏ⟨r−z2|mλ⟩Bλλ′W(X)⟨mλ′|r+z2⟩ =(2πℏ)Ngλ′λ(x)Bλλ′W(X). (19)

This result also follows from Eq. (II.1) by substituting Eq. (14) for and using the fact that

 ∫dxgλλ′(x)cνν′(x)=δλνδλ′ν′. (20)

Finally, in view of the definition of the projection operator , in place of Eq. (9) we may write

 Bλλ′W(X)=∫dxBPm(X)gλλ′(x). (21)

An analogous set of relations apply to the matrix elements of the density operator, , where . Taking the Wigner transform of we find

 ρm(X) = 1(2πℏ)N∫dzeip⋅z/ℏ⟨r−z2|^ρm(X)|r+z2⟩ (22) =1(2πℏ)Nρλλ′W(X)cλλ′(x).

Likewise, starting from the expression for the projected density,

 ^ρPm(X) = |mλ⟩⟨mλ|^ρm(X)|mλ′⟩⟨mλ′| (23) = |mλ⟩ρλλ′W(X)⟨mλ′|,

its Wigner transform is

 ρPm(X) = 1(2πℏ)N∫dpeip⋅z/ℏ⟨r−z2|^ρPm(X)|r+z2⟩ (24) =Pρm(X),

which, repeating the steps that gave Eq. (II.1), yields

 ρPm(X)=ρλλ′W(X)gλ′λ(x). (25)

Following the analysis given above that led to Eq. (9) for an operator, and using the relation

 ⟨r−z2|^ρm(X)|r+z2⟩=∫dpe−ip⋅z/ℏρm(X), (26)

 ρλλ′W(X) = (2πℏ)N∫dxgλλ′(x)ρm(X) (27) = (2πℏ)N∫dxgλλ′(x)ρPm(X).

These relations allow one to transform operators expressed in the subsystem basis to Wigner representations of operators in the basis of mapping states. The projected forms of the mapping operators and densities confine these quantities to the physical space and this feature plays an important role in the discussions of the nature of dynamics using the mapping basis. We now show how these relations enter the expressions for expectation values and evolution equations.

### ii.2 Forms of Operators in the Mapping Subspace

We first consider the equivalent forms that operators take, provided they are confined to the physical mapping space. Since

 ⟨mλ|∑ν^a†ν^aν|mλ′⟩=⟨mλ|mλ′⟩, (28)

is an identity operator in the mapping space. (Here we include the explicit summation on mapping states for clarity.) Using the definition of in Eq. (II.1), we may write the right side of Eq. (28) as

 ⟨mλ|mλ′⟩=∫dxgλλ′(x). (29)

The left side of may be evaluated by inserting complete sets of coordinate states and taking Wigner transforms so that an equivalent form for Eq. (28) is

 ∫dxgλλ′(x)∑νcνν(x)=∫dxgλλ′(x). (30)

Thus, we see that

 ∑νcνν(x)=12ℏ∑ν(r2ν+p2ν−ℏ)=1, (31)

provided it lies inside the integral.

This result has implications for the form of operators in the mapping basis. The matrix elements of an operator in the subsystem basis may always be written as a sum of trace and traceless contributions,

 Bλλ′W(X)=δλλ′(TrBW)/N+¯¯¯¯Bλλ′W(X), (32)

where is traceless. Inserting this expression into Eq. (14) for , we obtain

 Bm(X)=(TrBW)/N+¯¯¯¯Bλλ′W(X)¯¯cλλ′(x), (33)

provided appears inside the integral. Note that all subsystem matrix elements are of this form in view of Eq. (9). Here is the traceless form of .

As a special case of these results, we can write the mapping Hamiltonian, in a convenient form. The Hamiltonian matrix elements are given by

 Hλλ′W(X) = He(X)δλλ′+ϵλδλλ′+Vλλ′c(R) (34) ≡ He(X)δλλ′+hλλ′(R),

which can be written as a sum of trace and traceless contributions,

 Hλλ′W(X) = (He(X)+(Trh)/N)δλλ′+¯¯¯hλλ′(R) (35) ≡ H0(X)δλλ′+¯¯¯hλλ′(R).

The Hamiltonian can be written as . From this form for , it follows that

 Hm(X)=P22M+V0(R)+12ℏ¯¯¯hλλ′(R)(rλrλ′+pλpλ′), (36)

again, when it appears inside integrals with . We have used the fact that is symmetric to simplify the expression for in this expression. This form of the mapping Hamiltonian will play a role in the subsequent discussion.

### ii.3 Expectation values

Our interest is in the computation of average values of observables, such as electronic state populations or coherence, as a function of time. The expression for the expectation value of a general observable is

 ¯¯¯¯¯¯¯¯¯¯B(t)=∫dXTr(^BW(X)^ρW(X,t))= (37) ∫dXBλλ′W(X)ρλ′λW(X,t)=∫dXBλλ′W(X,t)ρλ′λW(X),

where the trace is taken in the quantum subsystem space. In the last line the time dependence has been moved from the density matrix to the operator, which also satisfies the QCLE.

The expression for the expectation value can be written in the mapping basis using the results in the previous subsection. For example, using Eq. (9) and the first line of Eq. (27) we find

 ¯¯¯¯¯¯¯¯¯¯B(t)=∫dX[∫dxBm(X,t)gλλ′(x)] (38) ×[(2πℏ)N∫dx′gλ′λ(x′)ρm(x′,X)] =∫dXBm(X,t)ρPm(X)=∫dXBPm(X,t)ρm(X),

where we have made use of the definition of the projection operator in Eq. (17) in writing the second equality. The projection operator can instead be applied to the observable in view of the symmetry in the expression and the resulting form is given in the last equality. We may write other equivalent forms for the expectation value. Starting from the second equality in Eq. (37) involving the time evolved density and the time independent operator, we obtain

 ¯¯¯¯¯¯¯¯¯¯B(t) = ∫dXBm(X)ρPm(X,t) (39) = ∫dXBPm(X)ρm(X,t).

From a computational point of view, the penultimate equality in Eq. (38) is most convenient since its evaluation entails sampling from the initial value of the projected density and time evolution of the operator.

### ii.4 Equations of motion

The most convenient form of the expectation value requires a knowledge of . Of course, if the solution to the QCLE in the subsystem basis, , is known, this definition can be used directly to construct ; however, the utility of the mapping basis representation lies in the fact that one can construct and solve the equation of motion for directly. The derivation of the evolution equation was given earlier. Kim et al. (2008) Here, we derive the evolution equations by taking account of the properties of mapping operators under integrals of in order to make connection with the projected forms of operators and densities. This will allow us to explore the domain of validity of the resulting equations.

The QCLE for an observable is expressed in the subsystem basis by taking matrix elements of the abstract equation with defined in Eq. (2):

 ddt⟨λ|^BW(X,t)|λ′⟩=−iℏ⟨λ|[^HW,^BW(X,t)]|λ′⟩ (40) +12⟨λ|({^HW,^BW(X,t)}−{^BW(X,t),^HW})|λ′⟩.

We may write this equation in terms of mapping variables using Eq. (9) as

 ∫dxgλλ′(x)ddtBm(X,t)= (41) ∫dxgλλ′(x)(−iℏ([^HW,^BW(X,t)])m(X,t) +12({^HW,^BW(X,t)}−{^BW(X,t),^HW})m(X,t)).

The mapping variables occur inside integrals of integral; i.e., they are projected onto the space of mapping states. Since the commutator and Poisson bracket terms in this equation involve products of operators, we must obtain the mapping form for a product of operators . The most direct way to make this transformation is to consider the product of operators as they appear in the subsystem basis and then use Eq. (9) for each matrix element:

 AλνW(X)Bνλ′W(X) = ∫dxAm(x,X)gλν(x) ×∫dx′gνλ′(x′)Bm(x′,X).

This expression does not lead to a useful form for the equations of motion. Instead we may write

 AλνW(X)Bνλ′W(X)=⟨λ|^AW(X)^BW(X)|λ′⟩ =⟨mλ|^Am(X)^Bm(X)|mλ′⟩ (43) =∫dxgλλ′(x)(^Am(X)^Bm(X))W(X)

Given that the Wigner transform of a product of operators is

 (^Am(X)^Bm(X))W=Am(x,X)eℏΛm/2iBm(x,X), (44)

where is the negative of the Poisson bracket operator on the mapping phase space coordinates, we obtain

 AλνW(X)Bνλ′W(X)= (45) ∫dxgλλ′(x)(Am(X)eℏΛm/2iBm(X)).

In Appendix A we establish the equality between this form for the matrix product and that given in Eq. (II.4). Inserting this result into Eq. (41), expanding the exponential operator and noting that the mapping Hamiltonian is a quadratic function of the mapping phase space coordinates, we obtain (details of the derivation are given in Ref. [Kim et al., 2008])

 ∫dxgλλ′(x)(ddtBm(X,t)=iLmBm(X,t)), (46)

where the mapping quantum-classical Liouville (MQCL) operator is given by the sum of two contributions:

 iLm=iLPBm+iL′m. (47)

The Liouville operator has a Poisson bracket form,

 iLPBm=−{Hm,}X=¯¯¯hλλ′ℏ(pλ′∂∂rλ−rλ′∂∂pλ) −(∂Hm∂R⋅∂∂P−PM⋅∂∂R), (48)

where denotes a Poisson bracket in the full mapping-environment phase space of the system, while

 iL′m=ℏ8∂hλλ′∂R(∂2∂rλ′∂rλ+∂2∂pλ′∂pλ)⋅∂∂P. (49)

In writing this form of the mapping Liouville operator we used the expression for the Hamiltonian given in Eq. (36). This is allowed since by Eq. (II.4) the operators appear inside integrals.

The formal solution of the equation of motion for is . The expectation value of this operator is given by (see Eq. (38))

 ¯¯¯¯¯¯¯¯¯¯B(t)=∫dX(eiLmtBm(X))ρPm(X)= (50) ∫dXBm(X)e−iLmtρPm(X)≡∫dXBm(X)ρPm(X,t),

where the evolution operator has been moved to act on the projected density using integration by parts. Thus, we see that the projected density satisfies

 ∂∂tρPm(X,t)=−iLmρPm(t). (51)

Making use of the above results, we can establish relations among the various forms of the expectation values and the dynamics projected onto the physical mapping states. From Eqs. (39) and (50) we have the relation . Differentiating both sides with respect to time and using the MQCLE we may write this equality as

 ∫dXBm(X)iLmPρm(X,t) (52) =∫dXBm(X)PiLmρm(X,t).

This identity, which is confirmed by direct calculation using the explicit form of in Appendix B, shows that commutes with the projection operator. Thus, evolution under the MQCL operator is confined to the physical mapping space.

## Iii Trajectory Description of Dynamics

A variety of simulation schemes have been constructed for the solution of the QCLE, some involving trajectory based solutions Grunwald et al. (2009); Martens and Fang (1996); Donoso and Martens (1998); Santer et al. (2001); Sergi et al. (2003); MacKernan et al. (2008); Wan and Schofield (2000, 2002); Horenko et al. (2002); Hanna and Kapral (2005); Hanna and Geva (2008). These schemes involve either ensembles of surface-hopping trajectories or correlations among the trajectories. A solution in terms of an ensemble of independent trajectories evolving by Netwonian-like equations is not possible Kapral and Ciccotti (1999).

### iii.1 Ensemble of entangled trajectories

A trajectory based solution of the MQCLE can also be constructed but the trajectories comprising the ensemble are not independent. Such entangled trajectory solutions have been discussed by Donoso, Zheng and Martens Donoso and Martens (2001); Donoso et al. (2003) for the Wigner transformed quantum Liouville equation. While our starting equation is very different, a similar strategy can be used to derive a set of equations of motion for an ensemble of entangled trajectories.

The MQCLE (1) can be written as a continuity equation in the full (mapping plus environment) phase space as

 ∂∂tρPm(X,t) = −∂∂X⋅j(X,t) = −∂∂X⋅[v(X;ρPm(X,t))ρPm(X,t)],

where the current has components:

 jrλ′=¯¯¯hλλ′ℏpλρPm,jpλ′=−¯¯¯hλλ′ℏrλρPm,jR=PMρPm, (54) jP=−∂Hm∂RρPm+ℏ8∂¯¯¯hλλ′∂R(∂2∂rλ′∂rλ+∂2∂pλ′∂pλ)ρPm.

The second equality in Eq. (III.1) defines the phase space velocity field through , which is a functional of the full phase space density.

We seek a solution in terms of an ensemble of trajectories, , where is the initial weight of trajectory in the ensemble. To find the equations of motion for the trajectories, consider the phase space average of the product of an arbitrary function with Eq. (III.1):

 ddt∫dXf(X)ρPm(X,t)= (55) ∫dX∂f(X)∂X⋅[v(X;ρPm(X,t))ρPm(X,t)],

where we have carried out an integration by parts to obtain the right side of the equality. Substitution of the ansatz for the phase space density into this equation gives

 (56)

from which it follows that the trajectories satisfy the evolution equations, . More explicitly we have

 ˙rλ = ∂Hm∂pλ