Analysis of the quantum-classical Liouville equation in the mapping basis
aa r X i v : . [ phy s i c s . c h e m - ph ] J un Analysis of the quantum-classical Liouville equation in the mapping basis
Ali Nassimi , Sara Bonella and Raymond Kapral Chemical Physics Theory Group, Department of Chemistry,University of Toronto, Toronto, Ontario M5S 3H6, Canada Dipartimento di Fisica and CNISM Unit 1, Universit`a La Sapienza, P.le A. Moro 2, 00185 Rome ∗ The quantum-classical Liouville equation provides a description of the dynamics of a quantumsubsystem coupled to a classical environment. Representing this equation in the mapping basisleads to a continuous description of discrete quantum states of the subsystem and may providean alternate route to the construction of simulation schemes. In the mapping basis the quantum-classical Liouville equation consists of a Poisson bracket contribution and a more complex term. Bytransforming the evolution equation, term-by-term, back to the subsystem basis, the complex term(excess coupling term) is identified as being due to a fraction of the back reaction of the quantumsubsystem on its environment. A simple approximation to quantum-classical Liouville dynamics inthe mapping basis is obtained by retaining only the Poisson bracket contribution. This approximatemapping form of the quantum-classical Liouville equation can be simulated easily by Newtoniantrajectories. We provide an analysis of the effects of neglecting the presence of the excess couplingterm on the expectation values of various types of observables. Calculations are carried out onnonadiabatic population and quantum coherence dynamics for curve crossing models. For theseobservables, the effects of the excess coupling term enter indirectly in the computation and goodestimates are obtained with the simplified propagation.
I. INTRODUCTION
The construction of algorithms for simulating thequantum dynamics of an arbitrary many-body systemis a long standing problem. Since the computationalcost of a quantum simulation scales exponentially withthe system size, a full quantum calculation is impracti-cal except for small systems. This fact has promptedthe construction of a variety of approximate methods forthe simulation of quantum dynamics. One class of ap-proximate schemes singles out a portion of the systemfor a full quantum treatment while its environment—therest of the system—is treated classically. A number ofsuch mixed quantum-classical schemes have been pro-posed and references to this literature can be found inreviews on this topic .We focus on one scheme of this type based on thequantum-classical Liouville equation (QCLE), ∂∂t ˆ ρ W ( X, t ) = − i ¯ h [ ˆ H W , ˆ ρ W ( t )] (1)+ 12 (cid:0) { ˆ H W , ˆ ρ W ( t ) } − { ˆ ρ W ( t ) , ˆ H W } (cid:1) , where [ · , · ] is the commutator and {· , ·} is the Poissonbracket. Here X = ( X , X , . . . , X N e ) = ( R, P ) =( R , R , . . . , R N e , P , P , . . . , P N e ) are the positions andmomenta of the 2 N e environmental degrees of freedomand the index W stands for a partial Wigner transform over the environmental degrees of freedom, so that thedensity matrix, ˆ ρ W ( X ), and Hamiltonian, ˆ H W ( X ), areoperators in the Hilbert space of the subsystem whilefunctions of the phase space variables for the bath de-grees of freedom. (The dependence on the phase spacecoordinates will be omitted when confusion is unlikely toarise.) The nature of QCL dynamics and the statisticalmechanics of systems following this dynamics have been described . For a review with references to the literatureon this topic, see Ref. [5].While the QCLE has a number of attractive featuresand has been shown to provide an accurate descriptionof the dynamics in many instances , it is difficult tosolve. A direct numerical integration of Eq. (1) requiresvery fine spatial grids and can be used for the study ofsmall systems . The QCLE can be cast in any basiswhich spans the Hilbert space of the quantum subsys-tem . When written in the quantum subsystem basis,the QCLE can be solved using a trajectory-based algo-rithm where trajectories are not independent of one an-other . Representation of QCLE in the adiabatic basisgives a more intuitive picture in terms of classical tra-jectories moving on single or mean Born-Oppenheimersurfaces . Also, diagonalizing the Hellman-Feynmanforce derived from the Born-Oppenheimer potential leadsto the force basis which yields yet a different route tosimulation of the QCLE .Another representation of QCL dynamics is in termsof the mapping basis . This basis, which provides adescription of a discrete quantum system in continuousvariables, has been used in a number of different appli-cations to quantum dynamical problems . When theQCLE is expressed in this basis one obtains an evolu-tion equation where the evolution operator consists ofa Poisson bracket contribution which can be solved bycharacteristics, and a more complex term that involvesspecific correlations of the quantum subsystem and clas-sical environment. Thus far, simulations of the QCLEin the mapping basis have neglected this more complexcontribution. In this paper we provide an analysis of themapping form of the QCLE that elucidates the natureof this neglected term. By transforming the QCLE inthe mapping basis back to the subsystem basis, wherean interpretation of the physical meaning of the differentterms is easier, we show that the neglected term (calledthe excess coupling term) accounts for a fraction of backreaction of the quantum subsystem on its environment.We also describe how the contribution of the excess cou-pling term to the average values of certain observablescan be estimated numerically.The outline of the paper is as follows: In Sec. II wesummarize the equations to be analyzed and introducesome key definitions, while in Sec. III we perform theterm-by-term analysis of the mapping equation based onthe back transformation to the subsystem basis. Sec-tion IV presents an analysis of the approximate formof the mapping QCLE that includes only the Poissonbracket term, called the Poisson bracket mapping equa-tion (PBME). This section also describes how neglectingthe excess coupling term affects average values of differ-ent types of observables. In Sec. V, the performance ofthe PBME is tested on two often-studied curve crossingmodels, and the effects of the excess coupling term inthe evolution operator are discussed. The last sectioncontains the conclusions of our investigation. II. SUBSYSTEM AND MAPPING BASISREPRESENTATIONS
For a system with partially Wigner transformed Hamil-tonian ˆ H W = P M + ˆ p m + ˆ V s (ˆ r ) + V e ( R ) + ˆ V c ( R, ˆ r ), where V e , ˆ V s and ˆ V c are, respectively, the environment, subsys-tem and coupling potentials, R and P are again the N e -dimensional coordinates and momenta of particles of theenvironment with mass M , and ˆ r = (ˆ r , ˆ r , . . . , ˆ r N s ) andˆ p = (ˆ p , ˆ p , . . . , ˆ p N s ) are the N s coordinate and momen-tum operators of the particles of the quantum subsystemwith mass m , the QCLE in the subsystem basis takes theform , ∂∂t ρ αα ′ W ( X, t ) = − i ˜ ω αα ′ ρ αα ′ W ( t )+ i ¯ h (cid:0) ρ αα ′′ W ( t ) V α ′′ α ′ c − V αα ′′ c ρ α ′′ α ′ W ( t ) (cid:1) + (cid:16) ∂V e ∂R · ∂∂P − PM · ∂∂R (cid:17) ρ αα ′ W ( t )+ 12 (cid:16) ∂V αα ′′ c ∂R · ∂ρ α ′′ α ′ W ( t ) ∂P + ∂ρ αα ′′ W ( t ) ∂P · ∂V α ′′ α ′ c ∂R (cid:17) = − i ˜ L αα ′ ,µµ ′ ρ µµ ′ W ( X, t ) . (2)Here and in the following we use the Einstein sum-mation convention where repeated indices are summed.The subsystem basis is defined by the eigenvalue prob-lem, ˆ h s | α i = ǫ α | α i , where ˆ h s = ˆ p m + ˆ V s . In Eq. (2),˜ ω αα ′ = ( ǫ α − ǫ α ′ ) / ¯ h and V αα ′ c = h α | ˆ V c | α ′ i . The lastline in Eq. (2) defines ˜ L αα ′ ,µµ ′ , the QCL operator in thesubsystem basis.Starting from this subsystem representation of theQCLE, we may transform it to the mapping basis in the following way. Suppose that there are N subsys-tem quantum states | λ i , λ = 1 , . . . , N . These subsystemstates may be mapped onto harmonic oscillator states viathe transformation, | λ i → | m λ i = | , · · · , λ , · · · , N i ,where h q | m λ i = h q , q , · · · , q N | , · · · , λ , · · · , N i = φ ( q ) · · · φ ( q λ − ) φ ( q λ ) · · · φ ( q N ) , (3)with φ and φ , respectively, being the ground and thefirst excited state wave functions of a harmonic oscillator.The creation and annihilation operators on the mappingstates, ˆ a λ = (ˆ q λ + i ˆ p λ ) / √ h and ˆ a † λ = (ˆ q λ − i ˆ p λ ) / √ h ,are used to define a mapping of operators as ˆ A = A λλ ′ | λ ih λ ′ | → ˆ A m = A λλ ′ ˆ a † λ ˆ a λ ′ . (4)For example, the partially Wigner transformed mappingform of the system Hamiltonian isˆ H m = P M + V e ( R ) + 12¯ h h λλ ′ ( R ) (cid:0) ˆ q λ ˆ q λ ′ + ˆ p λ ˆ p λ ′ − ¯ hδ λλ ′ (cid:1) ≡ P M + V e ( R ) + ˆ h m ( R ) . (5)Here we used the expression for the creation and annihi-lation operators in terms of the coordinates and momentaof the harmonic oscillators, together with the appropri-ate commutation relationship. Also, h λλ ′ ( R ) = h λ | ˆ h | λ ′ i with ˆ h = ˆ h s + ˆ V c ( R, ˆ r ) so that, more explicitly, h λλ ′ ( R ) = ǫ λ δ λλ ′ + V λλ ′ c ( R ) . (6)In writing Eq. (5), we assumed that h λλ ′ = h λ ′ λ . Follow-ing this prescription, after a further Wigner transforma-tion over the mapping variables, ρ m ( X, x ) = 1(2 π ¯ h ) N Z dz e ip · z/ ¯ h h r − z | ˆ ρ m ( X ) | r + z i , (7)where x = ( x , x , . . . , x N ) = ( r, p ) =( r , r , . . . , r N , p , p , . . . , p N ), the QCLE takes theform ∂∂t ρ m ( X, x, t ) = − h λλ ′ ¯ h (cid:18) p λ ′ ∂∂r λ − r λ ′ ∂∂p λ (cid:19) ρ m ( t ) (8)+ (cid:16) ∂H m ∂R · ∂∂P − PM · ∂∂R (cid:17) ρ m ( t ) − ¯ h " ∂h λλ ′ ∂R (cid:16) ∂ ∂r λ ∂r λ ′ + ∂ ∂p λ ∂p λ ′ (cid:17) · ∂ρ m ( t ) ∂P . Note that the Wigner transform variables r , p and z havethe same dimension as the number of quantum states, N .The Wigner transform of ˆ H m over mapping variables is H m ( X, r, p ) = Z dz e ip · z/ ¯ h h r − z | ˆ H m | r + z i = P M + V e ( R ) + 12¯ h h λλ ′ ( R ) (cid:0) r λ r λ ′ + p λ p λ ′ − ¯ hδ λλ ′ (cid:1) ≡ P M + V e ( R ) + h m ( R ) . (9)We can also write the mapping form of the QCLE as ∂∂t ρ m ( X, x, t ) = { H m , ρ m ( t ) } X,x (10) − ¯ h " ∂h λλ ′ ∂R (cid:16) ∂ ∂r λ ∂r λ ′ + ∂ ∂p λ ∂p λ ′ (cid:17) · ∂ρ m ( t ) ∂P . In Eq. (10), { H m , ρ m } X,x denotes a Poisson bracket withrespect to both the bath X and mapping x variables. Ne-glecting the last term of the dynamics in Eq. (10) yieldsa Hamiltonian system of equations which can be easilysolved using Newtonian trajectories. The last term hasa complex form involving both bath and mapping dif-ferential operators which make its interpretation difficultand precludes the implementation of simple algorithmsfor the simulation of the full mapping form of the QCLE. III. TRANSFORMING MAPPING DYNAMICSBACK TO THE SUBSYSTEM BASIS
In order to understand the nature of the last term inthe mapping QCLE, we shall transform each term in thisequation back to the quantum subsystem basis. Natu-rally, combining all the back-transformed terms, the orig-inal subsystem basis equation (Eq. (2)) is recovered; how-ever, the term-by-term transformation highlights howeach contribution in the mapping representation of theQCL operator is associated with a specific contributionin the subsystem basis form of the operator. This proce-dure leads to a simple physical interpretation of the lastterm in Eq. (10).We begin by recalling the expression for the Wignertransformed density matrix in the mapping basis. Us-ing the definition of a mapping operator in Eq. (4) andthe expression for the Wigner transform of the densityoperator in Eq. (7), one has ρ m ( X, x ) = 1(2 π ¯ h ) N Z dze ip · z/ ¯ h h r − z | ˆ ρ m ( X ) | r + z i = 1(2 π ¯ h ) N Z dze ip · z/ ¯ h ρ λλ ′ W ( X ) h r − z | ˆ a † λ ˆ a λ ′ | r + z i = ρ λλ ′ W ( X ) c λλ ′ ( x ) , (11)where c λλ ′ ( x ) = 12¯ h (2 π ¯ h ) N (12) × h r λ r λ ′ + i ( r λ p λ ′ − r λ ′ p λ ) + p λ p λ ′ − ¯ hδ λλ ′ i . Details of the calculation of c λλ ′ ( x ) are given in the Ap-pendix.Next, given ρ m ( X, x ) we consider how to recover theexpression for ρ λλ ′ W ( X ). The mapping relationship inEq. (4) ensures that the matrix elements of an opera-tor are the same in the subsystem and mapping bases,thus, ρ λλ ′ W ( X ) ≡ h λ | ˆ ρ W | λ ′ i = h m λ | ˆ ρ m | m λ ′ i . (13) Inserting resolutions of the identity in the coordinate rep-resentation, the last term in the equality above can bewritten as h m λ | ˆ ρ m | m λ ′ i = Z dydy ′ h m λ | y ih y | ˆ ρ m | y ′ ih y ′ | m λ ′ i (14)= Z drdz h m λ | r − z ih r − z | ˆ ρ m | r + z ih r + z | m λ ′ i Mean r = ( y + y ′ ) / z = y ′ − y coordi-nates were introduced in the last line to pave the way forthe introduction of the Wigner transform of the densityoperator in the mapping coordinates. The off-diagonalelement of the operator in the equation above can in factbe expressed as h r − z | ˆ ρ m | r + z i = R dpe − ip · z/ ¯ h ρ m ( X, x )(the inverse transform of the expression in Eq. (7)). Com-bining this expression with the identity in Eq. (13), weget ρ λλ ′ W ( X ) = (15) Z drdzdp e − ip · z/ ¯ h ρ m ( X, x ) h m λ | r − z ih r + z | m λ ′ i In the Appendix we show that the integral over z in theequation above can be performed analytically to obtain ρ λλ ′ W ( X ) = Z dx ρ m ( X, x ) g λλ ′ ( x ) , (16)where g λλ ′ ( x ) = 2 N +1 ¯ h e − x / ¯ h (17) × h r λ r λ ′ − i ( r λ p λ ′ − r λ ′ p λ ) + p λ p λ ′ − ¯ h δ λλ ′ i . (In the notation of this paper x = x · x = r + p = P λ ( r λ + p λ ). For clarity, in some places we shall use themore expanded forms of this condensed notation.) Rela-tionships analogous to those in Eqs. (11) and (16) holdfor the representation of any operator in the two bases,except for the fact that the multiplicative prefactors inEqs. (12) and (17) are interchanged.These results can be used to establish the relation-ship between the mapping and the subsystem bases. Toset the stage for our analysis, we multiply each term ofEq. (10) by g αα ′ ( x ) and integrate over the mapping co-ordinates to obtain ∂∂t ρ αα ′ W ( X, t ) = Z dx g αα ′ ( x ) { H m , ρ m ( t ) } X,x − ¯ h Z dx g αα ′ ( x ) " ∂h λλ ′ ∂R (cid:16) ∂ ∂r λ ∂r λ ′ + ∂ ∂p λ ∂p λ ′ (cid:17) · ∂∂P × ρ m ( t ) . (18)In the expression above, the integrand contains ρ m ( X, x, t ). We can use Eq. (11) to express this quantityin terms of the density matrix in the subsystem basis, ρ µµ ′ W ( t ), to get ∂∂t ρ αα ′ W ( X, t ) = Z dx g αα ′ ( x ) { H m , ρ µµ ′ W ( t ) c µµ ′ ( x ) } X,x − ¯ h Z dx g αα ′ ( x ) " ∂h λλ ′ ∂R (cid:16) ∂ ∂r λ ∂r λ ′ + ∂ ∂p λ ∂p λ ′ (cid:17) · ∂∂P × ρ µµ ′ W ( t ) c µµ ′ ( x ) . (19)In order to proceed with our analysis, we consider sepa-rately the Poisson bracket and second complex terms onthe right hand side. A. Transformation of the Poisson bracket term
The first contribution in the Poisson bracket, { H m , ρ m ( t ) } X,x , arising from derivatives with respectto the mapping variables is − h λλ ′ ¯ h (cid:16) p λ ′ ∂∂r λ − r λ ′ ∂∂p λ (cid:17) ρ m (see Eq. (8)). Its back transformation is given by P B Z drdp g αα ′ ( − h λλ ′ ¯ h ) (cid:18) p λ ′ ∂∂r λ − r λ ′ ∂∂p λ (cid:19) ρ µµ ′ W c µµ ′ (20)Direct substitution of the expressions for g αα ′ and f µµ ′ results in a sum of integrals that can be evaluated ana-lytically (they all are in the form of a polynomial timesa Gaussian function) to obtain P B − i ˜ ω αα ′ ρ αα ′ W + i ¯ h (cid:0) ρ αα ′′ W V α ′′ αc − V αα ′′ c ρ α ′′ αW (cid:1) (21)after substitution of h αα ′ from Eq. (6). This contributionrepresents the evolution of the quantum subsystem andincludes the influence of the environmental degrees offreedom through the coupling potential ˆ V c .Next, we consider the back transformation of the sec-ond contribution in Eq. (8) coming from derivatives withrespect to environmental coordinates, (cid:0) ∂H m ∂R · ∂∂P − PM · ∂∂R (cid:1) ρ m ( t ). Given the form of H m in Eq. (9), one canwrite this contribution as ∂H m ∂R · ∂ρ m ∂P − PM · ∂ρ m ∂R = ∂V e ∂R · ∂ρ m ∂P − PM · ∂ρ m ∂R + ∂h m ∂R · ∂ρ m ∂P . (22)Because ∂V e /∂R = − F e ( R ) and P/M are independentof the mapping variables, the back transformation of thefirst two terms in Eq. (22) is simple:
P B Z drdp g αα ′ (cid:16) ∂V e ∂R · ∂∂P − PM · ∂∂R (cid:17) ρ λλ ′ W c λλ ′ = − (cid:16) F e · ∂∂P + PM · ∂∂R (cid:17) ρ αα ′ W . (23)The last term in Eq. (22), arising from ∂h m ∂R , gives P B Z drdp g αα ′ h ( r λ r λ ′ + p λ p λ ′ − ¯ hδ λλ ′ ) × ∂h λλ ′ ∂R · ∂ρ µµ ′ W ∂P c µµ ′ ( x ) (24) This integral too can be performed after substitution ofthe g αα ′ and c µµ ′ functions to get P B (cid:16) ∂V αα ′′ c ∂R · ∂ρ α ′′ α ′ W ∂P + ∂ρ αα ′′ W ∂P · ∂V α ′′ α ′ c ∂R (cid:17) + 14 Tr ′ (cid:16) ∂ ˆ V c ∂R · ∂ ˆ ρ W ∂P (cid:17) δ αα ′ , (25)where Tr ′ indicates a trace over the quantum subsystemstates.Combining the results above, we obtain Z dx g αα ′ ( x ) (cid:8) H m , ρ m ( t ) (cid:9) X,x = − i ˜ L αα ′ ,µµ ′ ρ µµ ′ W ( t )+ 14 Tr ′ (cid:16) ∂ ˆ V c ∂R · ∂ ˆ ρ W ∂P (cid:17) δ αα ′ . (26)Thus, we find that the back transformation of the Pois-son bracket term yields the QCL operator in the subsys-tem basis defined in Eq. (2), plus an extra contribution Tr ′ (cid:16) ∂ ˆ V c ∂R · ∂ ˆ ρ W ∂P (cid:17) δ αα ′ . It is this term that is responsi-ble for any errors in the approximate simulations of theQCLE that include only the Poisson bracket term.An understanding of the physical meaning of this extracontribution can be obtained from the following consid-erations: First, we observe that the trace of the partiallyWigner transformed density matrix over the quantumsubsystem states is just the phase space density of theenvironment, ρ b ( X, t ) = Tr ′ ˆ ρ W ( X, t ). Taking the traceof Eq (2) yields, ∂∂t ρ b ( X, t ) = − (cid:16) F e ( R ) · ∂∂P + PM · ∂∂R (cid:17) ρ b ( t )+Tr ′ (cid:16) ∂ ˆ V c ∂R · ∂ ˆ ρ W ( t ) ∂P (cid:17) , (27)or, equivalently,Tr ′ (cid:16) ∂ ˆ V c ∂R · ∂ ˆ ρ W ( t ) ∂P (cid:17) = (cid:16) ∂∂t + F e ( R ) · ∂∂P + PM · ∂∂R (cid:17) ρ b ( t ) ≡ (cid:16) ∂∂t + iL e (cid:17) ρ b ( t ) , (28)where iL e is the classical Liouville operator for the envi-ronment in isolation from the quantum subsystem. Con-sequently, the trace term in Eq. (28) can be interpretedas the time variation of the phase space density of the en-vironment along the flow lines generated by the classicalevolution of the environment in isolation from the quan-tum subsystem. In a system where there is no couplingbetween the environment and the quantum subsystem,this term is zero by Liouville’s theorem. As a result,a non-zero value of the trace contribution can be inter-preted as an effect arising from the back reaction of thequantum subsystem on its environment. B. Excess coupling term
To complete the back transformation of the mappingQCLE to the subsystem basis, we must back transformthe complex last term in Eq. (10): − ¯ h Z drdp g αα ′ " ∂h λλ ′ ∂R (cid:0) ∂ ∂r λ ∂r λ ′ + ∂ ∂p λ ∂p λ ′ (cid:1) · ∂∂P × ρ µµ ′ W c µµ ′ = − (cid:16) ∂h λλ ′ ∂R · ∂ρ λλ ′ W ∂P + ∂h λλ ′ ∂R · ∂ρ λ ′ λW ∂P (cid:17) δ αα ′ = −
14 Tr ′ (cid:0) ∂ ˆ V c ∂R · ∂ ˆ ρ W ∂P (cid:1) δ αα ′ . (29)Since this result is again proportional to Tr ′ (cid:16) ∂ ˆ V c ∂R · ∂ ˆ ρ W ( t ) ∂P (cid:17) , we call this contribution a excess coupling term.Inserting these results into Eq. (19), the QCLE in thesubsystem basis (Eq. (2)) is obtained as expected. How-ever, from this derivation, we learn that the excess cou-pling term is equal and opposite in sign to the last termin Eq. (26). As a result, this contribution exactly can-cels the analogous term in the back-transformed Poissonbracket expression to yield the desired result. IV. APPROXIMATE MAPPING QCLE
The results of the previous section can be used to as-sess the utility of the approximate mapping QCLE whereonly the Poisson bracket contribution is retained in theLiouville operator. This Poisson bracket mapping equa-tion (PBME) is ∂∂t ρ m ( t ) ≈ { H m , ρ m ( t ) } X,x = − h λλ ′ ¯ h (cid:18) p λ ′ ∂∂r λ − r λ ′ ∂∂p λ (cid:19) ρ m ( t )+ ∂H m ∂R · ∂ρ m ( t ) ∂P − PM · ∂ρ m ( t ) ∂R . (30)It can be solved by characteristics and the resulting setof ordinary differential equations is dr λ ( t ) dt = 1¯ h h λλ ′ ( R ( t )) p λ ′ ( t ) ,dp λ ( t ) dt = − h h λλ ′ ( R ( t )) r λ ′ ( t ) ,dR ( t ) dt = P ( t ) M , dP ( t ) dt = − ∂H m ∂R ( t ) . (31)Consequently, the simulation of the dynamics describedby Eq. (30) is an easy task. From the results in Sec. III it follows that Eq. (30)is equivalent to the following equation in the subsystembasis: ∂∂t ρ αα ′ W ( X, t ) ≈ − i ˜ L αα ′ ,µµ ′ ρ µµ ′ W ( t )+ 14 Tr ′ (cid:16) ∂ ˆ V c ∂R · ∂ ˆ ρ W ( t ) ∂P (cid:17) δ αα ′ . (32) Because of the presence of the second term on the rightside, this form of the equation shows that the back reac-tion of the quantum subsystem on the dynamics of theenvironmental degrees of freedom is incorrectly describedand the solution by characteristics is only an approxima-tion to full QCL dynamics. Two features should be keptin mind when considering the error incurred in the use ofthe PBME. First, the effects of the environment on thequantum subsystem are fully accounted for in the Pois-son bracket, and the effects of the quantum subsystem onthe environment are also included in this term throughthe first term in PB3 in Eq. (25). Second, note the factorof 1 / B W ( X ) can bewritten as B ( t ) = Z dX B λλ ′ W ( X ) ρ λ ′ λW ( X, t )= Z dX B λλ ′ W ( X, t ) ρ λ ′ λW ( X )= Z dXdx B m ( x, X, t )˜ ρ m ( x, X ) , (33)In the first line of Eq. (33) the expectation value is com-puted in the subsystem basis, in the second line the timeevolution is moved from the density matrix to the oper-ator, while the last line gives the expectation value com-puted in the mapping basis. The quantities entering inthe mapping form of the expectation value are the timedependent observable, B m ( x, X, t ) = (2¯ h ) − B λλ ′ W ( X ( t )) (cid:16) r λ ( t ) r λ ′ ( t ) (34)+ p λ ( t ) p λ ′ ( t ) + i ( r λ ( t ) p λ ′ ( t ) − p λ ( t ) r λ ′ ( t )) − ¯ hδ λλ ′ (cid:17) , and the initial density ˜ ρ m ( x, X ) = R dx ′ f ( x, x ′ ) ρ m ( x ′ , X )with f ( x, x ′ ) = 1(2 π ¯ h ) N Z dzdz ′ h m λ | r − z ih r + z | m λ ′ i×h m λ ′ | r ′ − z ′ ih r ′ + z ′ | m λ i e − i ( p · z + p ′ · z ′ ) / ¯ h . (35)The integrals involved in the definition of the function f ( x, x ′ ) in Eq. (35) can be performed analytically. For atwo-level system the result is f ( x, x ′ ) = 8 π ¯ h h r · r ′ + p · p ′ ) + 2( r · p ′ − r ′ · p ) − ¯ h ( x + x ′ ) + ¯ h i e − ( x + x ′ ) / ¯ h , where r · r ′ = r r ′ + r r ′ , with similar expressions forother scalar products. Thus, these formulas allow one touse either the mapping or subsystem bases to computethe expectation values of any observable.While Eq. (33), along with Eqs. (34) and (31), providea simple means to compute the expectation value of anoperator using the PBME, it is convenient to considerthe computation of the general operator ˆ B W ( X ) usingthe PBME written in the subsystem basis in order togain further insight into the effects of the excess couplingterm on such average values. Taking the time derivativeof Eq. (33) and using Eq. (32) we find dB ( t ) dt ≈ − i Z dX B α ′ αW ( X ) ˜ L αα ′ ,µµ ′ ρ µµ ′ W ( t ) (36)+ 14 Z dX B α ′ αW ( X )Tr ′ (cid:16) ∂ ˆ V c ∂R · ∂ ˆ ρ W ( t ) ∂P (cid:17) δ αα ′ , = − i Z dX B α ′ αW ( X ) ˜ L αα ′ ,µµ ′ ρ µµ ′ W ( t ) − Z dX ∂B α ′ αW ( X ) ∂P · ∂V λλ ′ c ∂R ρ λ ′ λW ( t ) δ αα ′ , where, as usual, the Einstein summation convention isused. In the last line of this equation an integration byparts with respect to P was performed so that the mo-mentum derivative of the density matrix no longer ap-pears. Therefore, the effect of the excess coupling termon the evolution of the expectation value of ˆ B W ( X ) canbe estimated by computing the average values on theright side of this equation and determining their relativemagnitudes.If the initial value of the operator ˆ B W ( X ) is a functiononly of the configuration space coordinates, or indepen-dent of environmental phase space coordinates, then thelast excess coupling term in Eq. (36) will vanish becausethe momentum derivative is zero, and for this case wehave dB ( t ) dt = − i Z dX B α ′ αW ( X ) ˜ L αα ′ ,µµ ′ ρ µµ ′ W ( t ) . (37)For such observables the excess coupling term will enterthe computation of the average value indirectly throughthe density matrix on the right side. This can be seen bywriting the formal solution of Eq. (32) as ρ αα ′ W ( t ) = (cid:16) e − i ˜ L t (cid:17) αα ′ ,νν ′ ρ νν ′ W (0) (38)+ Z t dt ′ (cid:16) e − i ˜ L ( t − t ′ ) (cid:17) αα ′ ,νν ′ ∂V µµ ′ c ∂R · ∂ρ µ ′ µW ( t ) ∂P δ νν ′ , and inserting the result on the right of Eq. (37). The eval-uation of this higher order correction is difficult becauseof the convolution in the above equation. In particularthese considerations apply to the computation of impor-tant observables such as the quantum subsystem popu-lations and quantum coherence, since the initial values of these observables are independent of the phase spacecoordinates.If instead the observable is a function G ( P ) only ofthe environmental momenta, then the excess couplingterm enters the computation of the time rate of change ofits expectation value directly. Consider the expectationvalue of B α ′ αW ( X ) = G ( P ) δ αα ′ , G ( t ) = X α Z dX G ( P ) ρ ααW ( X, t ) . (39)If we compute the time derivative of this quantity usingEq. (36) we obtain, dG ( t ) dt = Z dX ∂G ( P ) ∂P · F αα ′ ρ α ′ αW ( X, t ) (40)+ N Z dX ∂G ( P ) ∂P · F αα ′ c ρ α ′ αW ( X, t ) , where F αα ′ = − ∂ ( V e ( R )+ V αα ′ c ( R )) /∂R is the total forceon the environment and F αα ′ c = − ∂V αα ′ c ( R ) /∂R is theforce contribution arising from the quantum subsystem-environment coupling. The first contribution on theright side of Eq. (40) is obtained from the explicitcomputation of − i R dX G ( P ) δ α ′ α ˜ L αα ′ ,µµ ′ ρ µµ ′ W ( t ), whilethe second contribution comes from the evaluation of R dX G ( P ) δ α ′ α Tr ′ (cid:16) ∂ ˆ V c ∂R · ∂ ˆ ρ W ( t ) ∂P (cid:17) δ αα ′ . These formu-las provide an indication of how the time evolution ofthe environmental momenta are influenced by the excesscoupling term. In particular, if the coupling is weak or ifthe environment is large and only a portion of the envi-ronment is directly coupled to the quantum subsystem,the total force on the environment will dominate the cou-pling force and this effect will be small. V. CURVE CROSSING MODELS
In this section we consider the simulation of quantum-classical Liouville dynamics in the mapping basis for twocurve crossing models studied earlier by Tully . Inparticular we carry out simulations using the PBMEand compare the results with numerically exact quan-tum dynamics using a discrete Fourier transformationmethod for these systems. Previous studies of the spin-boson model showed that the PBME yields results whichare in excellent agreement with the numerically exactquantum simulations for the entire range of system pa-rameters that were studied . In the spin-boson model,the coupling between the environment and the quantumsubsystem is linear in the environmental coordinates andis non-zero for all values of this quantity. The curvecrossing models studied here provide a further test ofthe PBME for a situation where the coupling among thedifferent degrees of freedom is nonlinear in the environ-mental coordinates and is localized in the configurationspace of the environment. A: Simple avoided crossing
The simple two-state avoided crossing is defined by thefollowing Hamiltonian matrix in the diabatic basis : h = A [1 − e − B | R | ] R | R | Ce − DR Ce − DR − A [1 − e − B | R | ] R | R | ! . (41)The corresponding total Hamiltonian in the mapping ba-sis is H m = P M + 12¯ h A (cid:2) − e − B | R | (cid:3) R | R | ( r + p − r − p )+ C ¯ h e − DR ( r r + p p ) . (42)The diagonal elements of the Hamiltonian in the dia-batic basis and the adiabatic energies for this model aresketched in Fig. 1. Our simulations are carried out usingatomic units (a.u.) and the parameters in these units aretaken to be A = 0 . B = 1 . C = 0 .
005 and D = 1 . ρ m ( X, x ) =˜ ρ s ( x ) ρ e ( X ), where ˜ ρ s ( x ) is the subsystem density matrixin the mapping basis and ρ e ( X ) is the distribution func-tion for the environment. The density ρ e ( X ) is chosento be the Wigner transform of a Gaussian wave packetcentered at R with momentum P : ρ e ( X ) = σ √ π ¯ h e − ( P − P ) σ / ¯ h √ πσ e − ( R − R ) /σ . (43)The mass M of the environmental degree of freedom istaken to be 2000 atomic mass units, while R = − . σ = 1.We assume the subsystem is initially in state 1, so thatin the mapping basis˜ ρ s ( x ) = 2 π ¯ h ( r + p − ¯ h e − x / ¯ h . (44)We consider the time evolution of the subsystem popu-lations, ρ ααs ( t ) , α = 1 ,
2, and coherence as measured bythe off-diagonal element of the subsystem density matrix, ρ s ( t ). These quantities are conveniently computed us-ing the formulas in Eq. (33) for a general operator byselecting B λλ ′ W = δ λα δ λ ′ α and B λλ ′ W = δ λ δ λ ′ , respec-tively. The asymptotic values of the populations in states1 and 2 after the system has passed through the inter-action region are shown in Fig. 1 as a function of theinitial momentum of the wave packet. We see that thePBME and exact quantum results are in excellent agree-ment for values of P above 10. Under this threshold,nuclear quantum effects play an important role so themost likely cause for the discrepancies in the figure is thebreakdown of the classical dynamics approximation forthe nuclei (note that usually, mixed quantum-classicalnon-adiabatic tests do not explore the range below thisvalue of the momentum for this model). -0.01-0.005 0 0.005 0.01 -3 -2 -1 0 1 2 3 E R ρ s α α P FIG. 1: (top) Diagonal elements of the Hamiltonian in the di-abatic basis (dashed lines) and adiabatic (solid lines) energiesfor the simple avoided crossing model. (bottom) Asymptoticpopulations of state 1 (solid squares with error bars) and state2 (solid circles with error bars) as a function of the initial mo-mentum of the wave packet, P , for the simple avoided cross-ing. The corresponding numerically exact quantum resultsare indicated by open squares for state 1 and open circles forstate 2. The exact results are connected with lines as a guidefor the eye. The time evolution of the real and imaginary partsof ρ s ( t ) are displayed in Fig. 2. These results for thequantum coherence are also in good agreement with thenumerically exact full quantum simulations. The com-parisons shown in these figures test two effects at thesame time: the validity of the QCLE and its approxima-tion by the PBME. Recall that from the considerationsin Sec. IV the effects of the excess coupling term enterthe populations and quantum coherence only as higherorder contributions.For the simple avoided crossing (and the dual avoidedcrossing considered below), V e ( R ) = 0 so that F αα ′ = F αα ′ c and the two integrals on the right of Eq. (40) areequal, and the two contributions differ only in their pref-actors. For this two-level model N = 2 and the PBMEprediction for the time rate of change of G ( t ) is 3 / V e ( R ) = 0 for this model. Thus, there is -0.3-0.2-0.1 0 0.1 0.2 0.3 0 150 300 450 ρ s t -0.1 0 0.1 0.2 0 400 800 1200 1600 ρ s t FIG. 2: The real (thick lines) and imaginary (thin lines) partsof ρ s ( t ) versus time for the simple avoided crossing model.The dashed lines denote the full quantum results while thesolid lines are the results of computations using the PBME.(top) P = 50, (bottom) P = 10. a substantial influence on the environmental momentabut, as noted above, these enter only in higher ordercorrections to the quantum subsystem populations andcoherences, contributing to the good agreement with theexact results seen in the figures. B: Dual avoided crossing
The dual avoided crossing model provides a further testof the theory. The Hamiltonian matrix in the diabaticbasis takes the form , h = (cid:18) Ce − DR Ce − DR − Ae − BR + E (cid:19) , (45)resulting in the mapping Hamiltonian, H m = P M + 12¯ h ( − Ae − BR + E )( r + p − C ¯ h e − DR ( r r + p p ) . (46)The system parameters were taken to be A = 0 . B = 0 . C = 0 . E = 0 .
05 and D = 0 .
06, again in -0.04-0.02 0 0.02 0.04 -4 -3 -2 -1 0 1 2 3 4 E R ρ s α α P FIG. 3: (top) Diabatic diagonal elements of the Hamiltonian(dashed lines) and adiabatic energies (solid lines) for the dualavoided crossing model. (bottom) Populations of the diabaticstates with same symbols as those in Fig. 1. atomic units. The initial conditions have the same formas for the simple avoided crossing, with R = −
10. Thediagonal matrix elements in the diabatic representationand the adiabatic energies for this model are plotted inFig. 3 (top) and the populations as function of P in thebottom panel of this figure. One can see that the agree-ment with exact quantum dynamics is very good overthe entire range of P , although there are small discrep-ancies in the magnitudes of maxima and minima in thepopulation curves.The time evolution of the real and imaginary partsof ρ s ( t ) for the dual avoided crossing model are shownin Fig. 4. The results for the quantum coherence areagain in good agreement with the numerically exact fullquantum simulations for this model. VI. CONCLUSION
By transforming the mapping form of the QCLE term-by-term back to the subsystem basis we were able toidentify the nature of the complex term in this equationthat makes its simulation difficult. This observation then -0.3-0.15 0 0.15 0.3 0 200 400 600 800 ρ s t -0.3-0.2-0.1 0 0 1000 2000 3000 4000 5000 ρ s t FIG. 4: The real and imaginary parts of ρ s ( t ) versus timefor the dual avoided crossing model. Same line types as inFig. 2. (top) P = 50, (bottom) P = 10. led to an analysis of the PBME, the approximate form ofthe mapping QCLE that retains only the Poisson bracketterm in the Liouville operator. The PBME is easily simu-lated by classical trajectories. However, when comparedto the full QCLE, it contains an extra term that canbe interpreted as being proportional to the effect of thequantum subsystem on the time evolution of the densitymatrix of the environment along the flow lines generatedby the classical evolution of the environment in isolation from the quantum subsystem. The excess coupling termin the full mapping QCLE cancels this contribution toyield the original QCLE.Earlier simulations on spin-boson model with bi-linear coupling to the bath, along with the present simu-lations on two curve crossing models with nonlinear bathcoupling, indicate that the correction terms are small andthe approximate PBME yields good agreement with theexact quantum results for these models. These resultssuggest that this approximate evolution equation will beuseful in many applications; however, the validity of thePBME must be tested for the specific problem under con-sideration and it may not always be a good approxima-tion to full QCL dynamics . The effects of the excesscoupling term on the time evolution of the observablescan be estimated by computing the expectation values inEq. (36) and these values provide one way of gauging theimportance of the deviations from the QCLE. Althoughcare must be exercised in the use of the PBME when ap-plied to a given system, the ease with which it can besimulated and its accuracy for many systems of interestmake it a powerful simulation scheme that should findincreased use in future applications. Acknowledgement
This work was supported in part by a grant from theNatural Sciences and Engineering Council of Canada.SB acknowledges funding from a SEED grant from IIT-Genova.
Appendix
In this Appendix we show how c λλ ′ ( x ) in Eq. (12) and g λλ ′ ( x ) in Eq. (17) can be computed. Calculation of c λλ ′ ( x ) Starting from its definition we have,(2 π ¯ h ) N c λλ ′ ( x ) = (ˆ a † λ ˆ a λ ′ ) W = 12¯ h Z dze ip · z/ ¯ h h r − z | [ˆ q λ ˆ q λ ′ + ˆ p λ ˆ p λ ′ + i (ˆ p λ ˆ q λ ′ − ˆ q λ ˆ p λ ′ )] | r + z i = 12¯ h Z dze ip · z/ ¯ h h ( r − z λ ( r + z λ ′ − ¯ h ∂∂z λ ∂∂z λ ′ + ¯ h ( r + z λ ′ ∂∂z λ − ¯ h ( r − z λ ∂∂z λ ′ i δ ( z )= 12¯ h h r λ r λ ′ − ¯ h Z dzδ ( z ) ∂∂z λ ∂∂z λ ′ e ip · z/ ¯ h − ¯ h Z dzδ ( z ) ∂∂z λ (cid:2) ( r + z λ ′ e ip · z/ ¯ h (cid:3) + ¯ h Z dzδ ( z ) ∂∂z λ ′ (cid:2) ( r − z λ e ip · z/ ¯ h (cid:3)i = 12¯ h h r λ r λ ′ − ¯ h ( i ¯ h p λ )( i ¯ h p λ ′ ) − ¯ h ( δ λλ ′ i ¯ h p λ r λ ′ ) + ¯ h ( − δ λλ ′ i ¯ h p λ ′ r λ ) i = 12¯ h h r λ r λ ′ + i ( r λ p λ ′ − r λ ′ p λ ) + p λ p λ ′ − ¯ hδ λλ ′ i . (47)0 Calculation of g λλ ′ ( x ) Similarly for g λλ ′ ( x ) we have g λλ ′ ( x ) = Z dze − ip · z/ ¯ h h m λ | r − z ih r + z | m λ ′ i = Z dze − ip · z/ ¯ h φ ( r − z · · · φ ( r λ − z λ · · · φ ( r N − z N φ ( r + z · · · φ ( r λ ′ + z λ ′ · · · φ ( r N + z N (cid:16) π ¯ h (cid:17) N/ h Z dze − ip · z/ ¯ h (cid:26) δ λλ ′ e [ − h (( r − z ) +( r + z ) )] · · · ( r λ − z λ e [ − h (( r λ − zλ ) +( r λ + zλ ) )] × · · · e [ − h (( r N − zN ) +( r N + zN ) )] + (1 − δ λλ ′ ) e [ − h (cid:0) r + z )] · · · ( r λ − z λ e [ − h (cid:0) r λ + z λ )] · · · ( r λ ′ + z λ ′ × e [ − h ( r λ ′ + z λ ′ )] · · · e [ − h ( r N + z N )] (cid:27) = 2 N +1 ¯ h e − x / ¯ h h δ λλ ′ ( r λ −
14 (2¯ h − p λ )) + (1 − δ λλ ′ )( r λ + ip λ )( r λ ′ − ip λ ′ ) i (48)= 2 N +1 ¯ h e − x / ¯ h h r λ r λ ′ − i ( r λ p λ ′ − r λ ′ p λ ) + p λ p λ ′ − δ λλ ′ ¯ h i , where we have used Eq. (3) together with φ ( x ) = ( π ¯ h ) / e − h x , φ ( x ) = ( π ¯ h ) / q h xe − h x . In this expressionrecall that z = ( z , z , . . . , z N ), p = ( p , p , . . . , p N ) and p · z = p z + p z + · · · + p N z N , with similar vector notationfor other unscripted quantities. ∗ Electronic address: [email protected]; Elec-tronic address: [email protected]; Electronic ad-dress: [email protected] G. D. Billing,
The Quantum Classical Theory (Oxford Uni-versity Press, 2003). J. C. Tully, in
Modern Methods for Multidimensional Dy-namics, Computations in Chemistry , edited by D. Thomp-son (World Scientific, 1998), p. 34. M. F. Herman, Annu. Rev. Phys. Chem. , 83 (1994). A. W. Jasper, C. Zhu, S. Nangia, and D. G. Truhlar, Fara-day Discuss. , 1 (2004). R. Kapral, Annu. Rev. Phys. Chem. , 129 (2006). K. Imre, E. Ozizmir, M. Rosenbaum, and P. Zweifel, J.Math. Phys. , 1097 (1967). M. Hillery, R. O’Connell, M. Scully, and E. Wigner, Phys.Rep. , 121 (1984). S. Nielsen, R. Kapral, and G. Ciccotti, J. Chem. Phys. , 5805 (2001). A. Donoso and C. Martens, J. Phys. Chem. A , 4291(1998). A. Donoso, Y. Zheng, and C. Martens, J. Chem. Phys. , 5010 (2003). C. Martens and J. Fang, J. Chem. Phys. , 4918 (1997). J. M. Riga, E. Fredj, and C. C. Martens, J. Chem. Phys. , 064506 (2006). M. Santer, U. Manthe, and G. Stock, J. Chem. Phys. ,2001 (2001). I. Horenko, C. Salzmann, B. Schmidt, and C. Schutte, J.Chem. Phys. , 11075 (2002). Y. Tanimura and S. Mukamel, J. Chem. Phys. , 3049 (1994). M. Toutounji, J. Chem. Phys. , 194520 (2006). G. Hanna and R. Kapral, J. Chem. Phys. , 244505(2005). G. Hanna and R. Kapral, J. Chem. Phys. , 164520(2008). G. Hanna and R. Kapral, Acc. Chem. Res. , 21 (2006). H. Kim, G. Hanna, and R. Kapral, J. Chem. Phys. ,084509 (2006). G. Hanna and E. Geva, J. Phys. Chem. B , 4048 (2008). R. Grunwald, A. Kelly, and R. Kapral, in
Energy TransferDynamics in Biomaterial Systems (Springer, 2009). R. Kapral and G. Ciccotti, J. Chem. Phys. , 8919(1999). A. Sergi, D. MacKernan, G. Ciccotti, and R. Kapral,Theor. Chem. Acc , 49 (2003). D. MacKernan, G. Ciccotti, and R. Kapral, J. Phys. Chem.B , 424 (2008). C. Wan and J. Schofield, J. Chem. Phys. , 7047 (2000). C. Wan and J. Schofield, J. Chem. Phys. , 494 (2002). H. Kim, A. Nassimi, and R. Kapral, J. Chem. Phys. ,084102 (2008). A. Nassimi and R. Kapral, Can. J. Chem. , 880 (2009). J. Schwinger, in
Quantum Theory of Angular Momentum (Academic press, 1966), p. 229. T. Holstein and H. Primakoff, Phys. Rev. , 1098 (1940). W. Miller and C. McCurdy, J. Chem. Phys. , 5163(1978). X. Sun, H. Wang, and W. Miller, J. Chem. Phys. ,7064 (1998). W. Miller, J. Phys. Chem. A , 2942 (2001). G. Stock and M. Thoss, in
Advances in Chemical Physics (Wiley-IEEE, 2005), vol. 131, p. 560. M. Thoss and G. Stock, Phys. Rev. A , 64 (1999). U. Muller and G. Stock, J. Chem. Phys. , 7516 (1998). G. Stock and M. Thoss, Phys. Rev. Let. , 578 (1997). S. Bonella and D. Coker, J. Chem. Phys. , 4370 (2003). S. Bonella and D. Coker, J. Chem. Phys. , 194102(2005). E. Dunkel, S. Bonella, and D. F. Coker, J. Chem. Phys. , 114106 (2008). Since the mapping oscillators are fictitious, we set mω = 1so that the mapping coordinates and momenta are definedto have dimensions of √ ¯ h . This set of evolution equations appears in Ref. [36] wherea mapping description of the dynamics was investigated in the context of a semiclassical formulation. A. Kelly, R. van Zon, J. Schofield, and R. Kapral, unpub-lished (2010). Efficient algorithms can be constructed to simulate theseequations that account for the time scale disparity be-tween the classical and mapping quantum degrees of free-dom. J. Tully, J. Chem. Phys. , 1061 (1990). D. Kosloff and R. Kosloff, J. Comp. Phys. , 35 (1983). Depending on the nature of the potential and the choice ofinitial conditions, there are situations where the dynam-ics described by this set of equations will yield unphysi-cal results unless modifications to the dynamics are intro-duced . S. Bonella and D. Coker, J. Chem. Phys.114