Robust non-adiabatic molecular dynamics for metals and insulators

# Robust non-adiabatic molecular dynamics for metals and insulators

L. Stella Department of Physics and Astronomy, University College London, Gower Street, London WC1E 6BT, UK London Centre for Nanotechnology, 17–19 Gordon Street, London WC1H OAH, UK    M. Meister Department of Physics and Astronomy, Queen’s University Belfast, Belfast BT7 1NN, UK    A.J. Fisher Department of Physics and Astronomy, University College London, Gower Street, London WC1E 6BT, UK London Centre for Nanotechnology, 17–19 Gordon Street, London WC1H OAH, UK    A.P. Horsfield Department of Materials, Imperial College London, South Kensington Campus, London SW7 2AZ, UK
###### Abstract

We present a new formulation of the correlated electron-ion dynamics (CEID) scheme, which systematically improves Ehrenfest dynamics by including quantum fluctuations around the mean-field atomic trajectories. We show that the method can simulate models of non-adiabatic electronic transitions, and test it against exact integration of the time-dependent Schrödinger equation. Unlike previous formulations of CEID, the accuracy of this scheme depends on a single tunable parameter which sets the level of atomic fluctuations included. The convergence to the exact dynamics by increasing the tunable parameter is demonstrated for a model two level system. This algorithm provides a smooth description of the non-adiabatic electronic transitions which satisfies the kinematic constraints (energy and momentum conservation) and preserves quantum coherence. The applicability of this algorithm to more complex atomic systems is discussed.

###### pacs:
31.15.Qg, 31.50.Gh

## I Introduction

Ordinary molecular dynamics (MD) Tuckerman and Martyna (2000) uses Hamiltonian equations of motion (EOM) to describe the evolution of an atomic system. The validity of that approach relies on the Born-Oppenheimer (BO) approximation which states that — in most cases — the atomic motion induces a slow (adiabatic) perturbation of the electronic dynamics. Therefore, if the system is originally prepared in an electronic eigenstate (e.g. the ground-state) for a given atomic configuration, it will evolve into the same electronic eigenstate for the evolved atomic configuration. Furthermore, according to the BO approximation, the ground-state electronic energy for a given atomic configuration can be employed as an effective (i.e. low energy) potential energy surface (PES) for the atomic motion.

The BO approximation is not longer applicable whenever electronic transitions between surfaces are relevant: for instance, where there is a non-radiative transition following an earlier photo-excitation. However, it is still possible to extend the scope of MD by allowing more than a single PES — one for every electronic state — and by providing a meaningful way to make transitions among PES i.e. non-adiabatic electronic transitions.

The extension of ordinary MD to deal with processes involving many PES has been extensively pursued over recent decades and some effective algorithms are available for atomistic simulations. The most used are: Ehrenfest dynamics (ED), Tully (1998); Horsfield et al. (2006) molecular dynamics with quantum transitions (MDQT), Hammes-Schiffer and Tully (1994); Tully (1998), mixed quantum-classical dynamics (MQCD), Kapral and Ciccotti (1999); Kapral (2006) and ab initio multiple spawning (AIMS) Ben-Nun et al. (2000); Ben-Nun and Martinez (2002). They have been employed to simulate quantum dissipative dynamics Mac Kernan et al. (2002); Sergi (2007a) as well as non-adiabatic processes such as photo-chemical reactions, Müller and Stock (1997); Kin et al. (2007) polaron formation in conjugated polymers, An et al. (2004), and proton transfer in solutions Hammes-Schiffer and Tully (1994). Although less successful in practical applications, we also mention other dynamical schemes that provide valuable theoretical insights: the frozen Gaussian approximation, Heller (1981) which is one of the pillars of AIMS, and the Gauss-Hermite wave-packet expansion, Adhikari and Billing (1999) which shares similarities to the method presented in this paper.

In order to extend the scope and the accuracy of the aforementioned methods, a new approach called correlated electron-ion dynamics (CEID) Horsfield et al. (2004a); Horsfield et al. (2005, 2006), has been introduced. This method has been mainly applied to study the heat production and dissipation in model metallic nanostructures, Bowler et al. (2005); Horsfield et al. (2006); McEniry et al. (2007) a problem relevant for nanotechnology. Indeed, other algorithms based on non-smooth dynamics, that is, those which allow for either sudden surface hopping (like MDQT) or sudden wave-packet spawning (like AIMS), are expected to be less efficient in the simulation of systems with a dense, gapless electronic spectrum, like a metal. That is a consequence of the slight adjustments of the average atomic positions and/or velocities that might need to be imposed after a sudden transition in order to ensure total energy and momentum are conserved. Tully (1998); Ben-Nun and Martinez (2002) If the number of crossings is large, the computational overhead due to these adjustments might be non-negligible. MQCD, although it is often implemented by using both surface hopping and mean-field evolutions, is based on a sound theory Nielsen et al. (2001); Sergi (2005) and conserves the kinematic constraints. On the other hand, time-translation invariance is valid only approximately in MQCD Nielsen et al. (2001) and numerical instabilities have so far limited the surface-hopping implementation of MQCD to relatively short time simulations. Sergi (2007a) Finally, ED — although it evolves smoothly — has been shown to poorly describe the atomic heating caused by electron-ion interaction. Horsfield et al. (2004b)

A further complication arises when a quantum sub-system is coupled with large quantum reservoirs (which are not treated in detail), as for a nanostructure connected to macroscopic leads. In this case, it is hard to define meaningful PES in terms of the sub-system degrees of freedom only, and so algorithms based on the surface hopping paradigm are expected to be less accurate. This problem is absent if the quantum sub-system is coupled to a classical dissipative environment; in this case MQCD or its latest variant Sergi (2007b) can give reliable results.

In principle, the CEID EOM form an exact, yet infinite, kinetic hierarchy which corrects ED by means of the so-called small amplitude moment expansion (SAME) of the Liouville equation. Horsfield et al. (2006) So far, a few schemes to truncate the hierarchy have been proposed Horsfield et al. (2004a) in order to simulate the CEID EOM and currently available algorithms are restricted to a mean-field second moment approximation. foo (a); Horsfield et al. (2005) Although those algorithms are accurate enough to describe — at least qualitatively — nanostructure heating, the existence of a practical truncation scheme which converges to the exact quantum dynamics has not been demonstrated until now: in this paper we show that a convergent truncation scheme actually exists and we provide a new practical CEID algorithm whose accuracy depends on a single tunable parameter. We provide a validation of our theory by comparing CEID results against exact integration of the time-dependent Schrödinger equation for a model two level system (2LS) — i.e. the simplest system which displays non-adiabatic transitions Landau (1932); Zener (1932) — in two contrasting parameter regimes.

The rest of the paper is organized as follows: In Sec. II we introduce the physics of the model 2LS employed to test our new CEID algorithm. In Sec. III we describe the exact algorithm by which we produced the benchmark calculations, while in Sec. IV the new CEID scheme is derived in detail. Finally, numerical results from simulations of the model 2LS are collected in Sec. V and the conclusions and perspectives of this work are discussed in Sec. VI.

## Ii A case study: the two level system

In this section we introduce the model 2LS that we employed to test the convergence properties of our CEID scheme. [Numerical findings are reported in Sec. V.] Here we discuss the physics we expect to address before explaining the details of the algorithms we use.

It is widely recognized that a one-dimensional 2LS illustrates many of the fundamental features of a non-adiabatic system and, at the same time, it retains the simplicity of a low-dimensional model. Tully (1990); Martens and Fang (1997) In general, the Hamiltonian for a system made by electrons and ions (in the absence of external fields) can be written as follows:

 ^H=^P2/2M+^He(^R), (1)

where and are the quantum operators for the atomic position and momentum, while the electronic dependence of is collected into . In particular, the first term in the RHS of Eq. (1) accounts for the atomic kinetic operator while a sort of atomic potential operator is described by the second term.

There is a lot of freedom in constructing the potential term, , but in this paper we focus only on the following parametrization:

 ^He(R)=(12K(R−R0)2−fcR−fcR12KR2+Δε) (2)

which describes two parabolic PES (diagonal entries) linearly coupled through a kind of dipolar interaction (off-diagonal entries). This is a non-adiabatic representation of the electronic PES in which the electronic basis is independent of the atomic coordinate. The adiabatic representation can be obtained as usual by diagonalizing . Eq. (1) and (2) depend on the following parameters: the atomic mass , the harmonic constant , the electron-ion coupling constant , the surface displacement , and the PES energetic offset .

Since both the PES are confining, we expect to see periodic electronic transitions between the two PES driven by the electron-ion interaction. For instance, a state can be prepared as the atomic ground-state on the upper electronic PES. The atomic position experiences quantum oscillations and so the system cannot be exactly localized in the minimum of the PES ( according to our parametrization). Since this initial state is not an eigenstate of the interacting Hamiltonian (i.e. for ), the system will eventually make a transition into the lower PES. We stress that this process must conserve the total energy so that an atomic transition must accompany the electronic transition. For instance, the electronic process described above can be viewed as an initial decay since the atomic potential energy is effectively decreased. Therefore, an increase of the atomic kinetic energy is expected as a consequence of this decay in order to conserve the total energy. This, in a nutshell, is the heating of an atomic degree of freedom caused by the electron-ion interaction. Horsfield et al. (2004b) However, we stress that the aim of the present work is not the study of a quantum decay process, but the illustration of the convergence properties of a new CEID algorithm. Other models must be used to address the physics of quantum thermalization. For instance, the spin-boson model Leggett et al. (1987) describes a 2LS coupled to an environment made by a collections of many quantum harmonic oscillators. This model can be effectively simulated Makri and Thompson (1998); Mac Kernan et al. (2002) and it might provide a future test case for our new CEID algorithm.

There is an interesting general feature displayed by the electronic dynamics of our model 2LS. According to the initial condition described above, the electronic transition is due to the quantum fluctuations only, because the dipolar interaction is exactly zero for a classical atom perfectly localized in the minimum of the upper PES. On the other hand, quantum fluctuations are completely neglected in the sort of mean-field description of the atomic motion employed in ED. As a consequence, we do not expect ED to reproduce the initial transition from the upper PES to the lower PES, which can be thought of as spontaneous phonon emission. If confirmed, this behavior will provide further evidence that the exchange of energy from the electronic to the atomic degrees of freedom cannot be properly addressed by ED. Horsfield et al. (2004b)

In order to be as transparent as possible when discussing the dynamical features of our model 2LS, we choose to measure the values of the parameters in the Hamiltonian in terms of natural units. The natural energy scale is given by the harmonic quantum, , where , and so the time can be measured in units of the harmonic period, . Introducing a mass scale, , (its actual value is not important here), a length scale and a linear momentum scale are immediately obtained: and , respectively.

For our purposes, not all the parameters need to be varied during the numerical experiments. First of all, we fixed both the atomic mass () and the harmonic constant () and then we took the PES offset to be equal to one harmonic quantum (). This is equivalent to saying that the atomic ground-state on the upper PES, , has got exactly the same energy as the first harmonic excitation on the lower PES, . [Here we are neglecting the electron-ion coupling and using the notation introduced in appendix D.] If the coupling constant is small enough (see appendix D), the time-evolution of our 2LS can be understood starting from these two states. We confined this study to this weak coupling regime, i.e. . Finally, we report in this paper the numerical results for two different 2LS geometries only: the unshifted 2LS, with , and the shifted 2LS, with . The PES for these two cases are shown in Fig. 1(a) and Fig. 2(a), respectively. Although many other 2LS geometries have been studied, e.g. with larger or with other values of (including the possibility of different harmonic force constants for the lower and the upper PES), they gave numerical outcomes qualitatively similar to either the shifted or unshifted 2LS, and so the details are not reported here.

A linear combination of the two low-lying resonant states might be employed to give a qualitative account of our model 2LS dynamics i.e. the wave-function at time might be approximated as:

 ψ(t)≃c0(t)|χ(u)0⟩+c1(t)|χ(l)1⟩, (3)

where and are time-dependent complex coefficients. This state would give a distribution of the atomic position more or less localized around each of the two classical equilibrium positions, namely for the upper PES and for the lower PES. Eq. (3) might not be a good ansatz for the evolved state as soon as the displacement is large enough. Indeed, even a classical atom can be found far from its equilibrium position whenever this is energetically allowed. In this case, we should be able to describe an atomic wave-packet almost localized far from both and . Therefore, we will produce a manifest physical inconsistency if we assume that the generic 2LS state can be always well approximated by Eq. (3). This physical inconsistency can be easily fixed by taking a longer expansion of the exact evolved state in terms of the harmonic excitations of the atomic degrees of freedom. Unfortunately, there is a cost to pay: the longer — and so the more accurate — the expansion, the more time-consuming will be the simulation. [See Sec. VI for a further discussion of this point.]

In the next two sections, two different ways to compute non-adiabatic dynamics are considered in detail. In Sec. III a method based on the numerical diagonalization of the Hamiltonian of the full system is presented, while a new formulation of CEID is introduced in Sec. IV.

## Iii Exact integration

The ‘exact’ method of integration is based on numerical diagonalization of the full Hamiltonian matrix and subsequent time evolution exploiting the eigenvectors and eigenvalues obtained in the diagonalization step. While it is highly accurate, in practice, such an approach naturally has to be limited to systems with a small number of degrees of freedom. However, for our purposes it provides an ideal scheme for producing benchmark results. A central quantity in this paper is the Wigner transform (WT) of an operator. Originally, the WT of a wave-function was defined in Ref. Wigner, 1932 as

 W(R,P) = 1(ℏπ)n∫ψ∗(R+s)ψ(R−s)exp(2iℏP⋅s)dns = (4)

where , , and the time-dependence has been suppressed. The latter form in Eq. (III) is more common nowadays. In a straightforward extension of this definition, for an operator , expandable in basis states , , the WT is given by

 (5)

Of particular interest to us here, however, are WT with respect to the degrees of freedom of a subsystem, i.e. partial WT. These appear naturally when the system can be divided into subsystems on physical grounds. In this paper, where the system is a molecule, the system can be split into an electronic subsystem and the subsystem of the ions/nuclei. Consider an operator acting on a system divisible into subsystems with basis sets and . The operator is assumed to be expandable as , and its partial WT with respect to the -subsystem is

 ^Cw,A(R,P)=∑ab|a⟩Cabw,A(R,P)⟨b|, (6)

where, with the position representation of ,

 Cabw,A(R,P)=∑ABCABab1(2πℏ)n∫exp(iℏP⋅s)Φ∗B(R+12s)ΦA(R−12s)dns. (7)

Note the important fact that is still an operator in the -subsystem. As in the rest of the paper the only WT used are partial WT with respect to ‘ionic’ or ‘atomic’ degrees of freedom, we shall henceforth write instead of , and also refer to the -subsystem as the atomic and the -subsystem as the electronic subsystems.

The time-evolution of the system is generated by a Hamiltonian, whose eigenvalues and eigenvectors shall be and , respectively. The basis states and of the subsystems are taken to be time-independent from now on. We can expand an eigenstate in the product basis , or vice versa, , where we have the relations and . The expansion coefficients and are time-dependent,

 CAan(t)=CAan(0)exp(−iℏEnt); (8)

in any particular case the numerical diagonalization will provide us with the and the coefficients , which make up the eigenvectors, i.e. is the -component of eigenvector .

We now can expand an operator in the eigenstates or the product states:

 (9)

with

 GABab(t)=∑nmCAan(0)Gnm(CBbm)∗(0)exp(−iℏ(En−Em)t). (10)

The partial WT with respect to the atomic subsystem is then

 ^Gw(R,P,t)=∑ab|a⟩Gabw(R,P,t)⟨b|, (11)

where in turn

 Gabw(R,P,t)=∑ABFBA(R,P)GABab(t), (12)

and

 FBA(R,P)=1(2πℏ)n∫exp(iℏP⋅s)Φ∗B(R+12s)ΦA(R−12s)dns. (13)

Specializing to the particular case of the 2LS and its actual implementation on a computer, we introduce dimensionless quantities using the scale factors mentioned in Sec. II. In particular we obtain the dimensionless wave-function , the dimensionless version of , and dimensionless eigenvalues . Using Eq. (10), the factoring Eq. (12), and Eq. (13), all the components of can be computed as functions of (or their dimensionless counterparts). Note that this does not involve the numerical solution of a (potentially partial) differential equation, so we are not required to advance over many small time-steps in order to reach a given value of . Limitations are, however, introduced by the need to truncate the oscillator basis to a finite number of states. As the purpose of the exact approach is to provide a benchmark for the CEID method, some care has to be taken to avoid truncation errors here. First one has to decide how many product and eigenbasis states to use in expansions like Eq. (10), say . Then, the diagonalization has to be carried out using a number of product basis states such that the lowest eigenvalues and corresponding eigenvectors are well converged; typically, we used and . The operator one is considering (in what follows it will be the density operator) should only have negligible coupling between states with index equal to or less than () to states with index larger than . These first states, which in the diagonalization step are produced as -component quantities, should not have significant contributions from product basis components with index greater than . Only if these conditions are met, and the initial conditions for the operator are chosen to involve only the first states, can we consider the numerical results for the time-evolution reliable and an adequate benchmark for CEID. The quality of the choice of can therefore only be assessed after diagonalization.

From the results produced for the purpose of comparison we show the occupations , , of the electronic levels, and expectation values of position, momentum and the variance of the position, as functions of time. These have been calculated via the WT of the density operator of the system, by numerical evaluation of the integral

 Na(t)=∬ρaaw(R,P,t)dRdP, (14)

in case of the occupations, and of

 ⟨f(R,P,t)⟩=2∑a=1∬f(R,P,t)ρaaw(R,P,t)% dRdP (15)

for , , , respectively, in the rest of the cases.

## Iv Correlated electron-ion dynamics

In this section we describe a new formulation of the correlated electron-ion dynamics (CEID) while the original formalism can be found in Ref. Horsfield et al., 2004a. [See also Ref. Horsfield et al., 2005 and Ref. Horsfield et al., 2006 for further details.] For the sake of simplicity, the new CEID EOM has been derived here only for the one-dimensional case although multi-dimensional EOM are also known. foo (b)

We start from the well known quantum Liouville equation:

 ˙^ρ=1iℏ[^H,^ρ], (16)

which is the EOM for the density matrix of the system. [We use a dot to indicate time-derivative.] Unfortunately, a direct integration of Eq. (16) is exceedingly time-consuming because it scales approximately as the cube of the Hilbert space dimension which is very large in most cases of interest. On the other hand, since atoms are much heavier than electrons, an expansion of their motion around the classical trajectories is often justified. It turns out that this kind of expansion cuts off the quantum fluctuations of the atomic degrees of freedom and so it effectively reduces the Hilbert space dimension. For instance, simulations of the semi-classical limit of Eq. (16) Martens and Fang (1997) have been shown to reproduce — at least qualitatively — the correct non-adiabatic dynamics of a few interesting test-cases.

More generally, the density matrix can be partially expanded with respect to the atomic degrees of freedom by means of a complete orthonormal system (COS). Kapral and Ciccotti (1999) By using the standard Dirac’s bra and ket notation, this expansion can be expressed as:

 ^ρ=∞∑n=0∞∑m=0|ϕn⟩^ρn,m⟨ϕm|, (17)

where the functions are a COS in the atomic subspace. As a consequence, in Eq. (17) the electronic degrees of freedom are included in the matrix coefficients . A natural choice dictated by the kind of 2LS physics introduced in Sec. II (i.e. confining PES) is to use the simple harmonic oscillator (SHO) eigenfunctions as atomic COS. We stress here that, although these functions are usually centered around a classical equilibrium point, any other reference point can be taken instead (see below).

Since all the observables can be expanded as in Eq. (17), all the operations involving observables (e.g. averages) can be worked out by means of the observable matrix coefficients only. [In general, the matrix coefficients of the observable are given by: .] For instance, the total energy of the system is given by:

 Etot=Tr{^H^ρ}=∞∑n=0∞∑m=0Tre{^Hm,n^ρn,m}. (18)

[The two traces, and , apply to different linear spaces, namely the whole Hilbert space and the electronic subspace: see Sec. III.] As a further example, the EOM for are obtained by plugging Eq. (17) into Eq. (16):

 ˙^ρn,m=1iℏ∞∑k=0[^Hn,k^ρk,m−^ρn,k^Hk,m]. (19)

It is easy to prove that is in fact a constant of motion by taking the time-derivative of Eq. (18) and then by using Eq. (19).

The complete set of EOM, Eq. (19), cannot be directly simulated because it is not finite. Therefore, we must make an approximation and, quite naturally, we set to zero (as a matrix) every matrix coefficient with indices greater than , the CEID order. Nevertheless, after this truncation, the EOM are still fully quantized (and expressed by a proper Lie bracket), but they are restricted to a smaller Hilbert subspace.

As for the exact scheme described in Sec. III, in order to make the classical limit of Eq. (19) more manifest, we use a partial Wigner transform (WT) i.e. a WT taken only with respect to the atomic degrees of freedom. Therefore, the partial WT of the operator , , is still an operator in the electronic subspace but it explicitly depends on what are now the classical atomic position and momentum . In the context of non-adiabatic MD, similar partial WT have been already considered Kapral and Ciccotti (1999); Thorndyke and Micha (2005) and they have been shown to provide correct numerical results. In order to avoid confusion, we stress that, although WT seems to map a quantum operator into a classical distribution, the dynamics remains non-classical because the WT of the product of two operators is in general not the product of the WT of the two operators:

 (^A^B)w=^Aw⋆^Bw, (20)

where is the non-commutative Moyal product: Grönewald (1946); Moyal (1949)

 ⋆=exp[iℏ2(←∂R→∂P−←∂P→∂R)]. (21)

[The arrows indicate the directions in which the derivative operators act.] The WT of any operator can be expanded as in Eq. (17) by using the WT of the basis operators . As a consequence, the matrix coefficients of an operator are the same both in the original Hilbert space and the in transformed one.

The WT of the Liouville equation can be formally stated as:

 ˙^ρw(t)=1iℏ(^Hw⋆^ρw−^ρw⋆^Hw). (22)

It can be shown that from Eq. (22) the usual Hamilton-Ehrenfest equations (i.e. the EOM for and ) can be derived: Horsfield et al. (2006)

 ⎧⎪ ⎪⎨⎪ ⎪⎩˙¯R=¯P/M,˙¯P=¯F=−Tr{(∂^He∂R)^ρ}. (23)

It is also reasonable to take the phase-space trajectory as a zero-order approximation of the true atomic dynamics if the quantum fluctuations are not too large. [The semi-classical limit of the WT is explained more extensively in Ref. Wigner, 1932 (see also Ref. Grönewald, 1946).] This fact suggests that instead of taking the origin as a reference point in the phase-space (i.e a fixed reference frame), one can take advantage of Eq. (23) and use as reference (i.e. a mobile reference frame). After this mobile reference frame transform, the EOM becomes:

 ˙^ρw=1iℏ[^Hw⋆^ρw−^ρw⋆^Hw]+(∂^ρw∂¯R)¯PM+(∂^ρw∂¯P)¯F. (24)

Although it is not apparent from Eq. (24), the EOM can be still expressed by a proper Lie bracket — as in Eq. (22) — by means of a different time-translation generator i.e. by a different Hamiltonian. [Mathematical details can be found in appendix B.]

At this stage two paths can be followed. In the first case, the Moyal product is expanded in (usually up to first order Thorndyke and Micha (2005); Kapral and Ciccotti (1999)) and a quantum-classical extension of the Liouville equation is obtained. However, although this approach is physically appealing, one must be aware that the EOM obtained this way cannot be formulated through a proper Lie bracket, Caro and Salcedo (1999); Prezhdo (2006) and a generalized non-Hamiltonian bracket should be introduced instead. Sergi (2005) As a consequence, the evolution of a composite operator, , might be different from the composition of the separated evolutions, , Caro and Salcedo (1999); Nielsen et al. (2001) (because the non-Hamiltonian bracket does not define a proper derivative), and the conservation of dynamical symmetries might be a problem Caro and Salcedo (1999) (due to the violation of the Jacobi identity). It must be also recalled that the difference between the non-Hamiltonian and Lie brackets is only of order Nielsen et al. (2001) and that the non-Hamiltonian structure arises in a quite natural way for open systems, whether classical or quantum. Sergi (2005)

By following the other route, one uses the exact expression for the Moyal product and takes the WT of Eq. (17) as a natural way to truncate . [We have found a direct expansion in the transformed space — e.g. by using weighted orthogonal polynomials — to cause dangerous instabilities in the truncated dynamics.] Therefore, the action of the truncation super-operator can be defined as follows:

 Tw[^ρw(R,P,t)] =Tw[^ρw(¯R+ΔR,¯P+ΔP,t)] ≡N∑n=0N∑m=0^ρn,m(¯R,¯P,t)Pn,m(ΔR,ΔP), (25)

where is the WT of in the mobile reference frame. [Properties of these functions are given in appendix A.] We opted for this approach because it is still a full quantum scheme — but in a truncated Hilbert space — and the EOM for the density matrix can still be formulated by means of a proper Lie bracket (see appendix B).

Although an analytical expression of is known, Grönewald (1946) it is far more convenient to state a set of recurrence relations by taking advantage of the well-known SHO algebra. Details can be found in appendix B; here we report only the final form of the CEID EOM for the matrix coefficients (in the mobile reference frame):

 ˙^ρn,m=−b204iℏM(√(n+2)(n+1)^ρn+2,m−(2n+1)^ρn,m+√n(n−1)^ρn−2,m+−√m(m−1)^ρn,m−2+(2m+1)^ρn,m−√(m+2)(m+1)^ρn,m+2)++1iℏ[^He(¯R),^ρn,m]−a0iℏ(Δ^F(¯R)√n+12^ρn+1,m+Δ^F(¯R)√n2^ρn−1,m+−√m2^ρn,m−1Δ^F(¯R)−√m+12^ρn,m+1Δ^F(¯R))+a204iℏ(^K(¯R)√(n+2)(n+1)^ρn+2,m++^K(¯R)(2n+1)^ρn,m+^K(¯R)√n(n−1)^ρn−2,m−√m(m−1)^ρn,m−2^K(¯R)+−(2m+1)^ρn,m^K(¯R)−√(m+2)(m+1)^ρn,m+2^K(¯R)), (26)

where , , and . [Terms involving higher derivatives of should also appear in Eq. (26), but vanish in this case since the 2LS Hamiltonian we want to study is quadratic — see Eq. (2).] We recall that, according to our truncation scheme, one must neglect in the RHS of Eq. (26) those matrix coefficients whose indices are greater than the CEID order. Those equations, along with Eq. (23), have been used to simulate the 2LS dynamics described in Sec. II. In particular, the current implementation uses a second order Runge-Kutta non-adaptive algorithm to integrate Eq. (26) and the standard velocity-Verlet algorithm to integrate Eq. (23). [A time-step in our natural unit (see Sec. II) has been found to be appropriate for the precision required by the comparison between CEID an exact approaches reported in Sec. V.] At every integration step, the averaged coordinates, and , are evolved according to Eq. (23) for half a time-step, then the matrix coefficients are propagated through Eq. (26) for a whole time-step, and finally the averaged coordinates are evolved by another half time-step. We verified that the accuracy achieved by this kind of symmetric Trotter decomposition is greater than what is obtained by means of a single evolution of the averaged coordinates for a whole time-step followed (or preceded) by a matrix coefficients propagation for a whole time-step. Numerical results can be found in Sec. V.

We now briefly discuss the link between our new formulation and the original CEID.

### iv.1 Comparison with former CEID integration schemes

At variance with the scheme described so far, the original formulation of CEID makes use of a completely different expansion which directly provides EOM for the moments of the density matrix. Horsfield et al. (2005, 2006) The most relevant CEID moments are: , , and , where is the partial trace with respect to the atomic degrees of freedom, , and . [Higher order moments must be carefully defined because and do not commute.] On the other hand, analogous objects can be also introduced in the new formulation:

 ^μn,m(t)=12πℏ∫dRdPΔRnΔPm^ρw(R,P,t). (27)

By using the property of the WT, it is easy to find a link between the new and the original notation: , , and . Similar relations for higher CEID moments can be stated, but some extra attention must be payed in the derivation due to the non-trivial commutation relations between positions and momenta. It is worth noting that CEID moments provide valuable information about the system. For instance, the quantities give the probability of observing the system on the -th PES and the average force (see Eq. (23)) can be easily computed from . Higher moments can be used to study electron-ion correlations. Horsfield et al. (2005)

Moments defined in Eq. (27) can be expressed in terms of the matrix coefficients by means of the following linear transform:

 ^μn,m=∑r,sAn,mr,s^ρr,s, (28)

where

 An,mr,s=12πℏ∫dRdPΔRnΔPmPr,s(R,P). (29)

As usual, a set of recurrence relations for can be found and — at least in theory — CEID moments of any order can be computed. In practice, only the low lying moments are relevant and here we give a short selection of them:

 ^μ0,0 = N∑n=0^ρn,n, (30a) ^μ0,1 = −ib0N∑n=0√n2[^ρn,n−1−^ρn−1,n], (30b) ^μ1,0 = +a0N∑n=0√n2[^ρn,n−1+^ρn−1,n], (30c) ^μ0,2 = −b202N∑n=0[√n(n−1)^ρn,n−2−(2n+1)^ρn,n+√n(n−1)^ρn−2,n], (30d) ^μ1,1 = −ia0b02N∑n=0[√n(n−1)^ρn,n−2−√n(n−1)^ρn−2,n], (30e) ^μ2,0 = +a202N∑n=0[√n(n−1)^ρn,n−2+(2n+1)^ρn,n+√n(n−1)^ρn−2,n]. (30f)

As anticipated in Sec. I, the zero order CEID (i.e. for ) is equivalent to the ED: Horsfield et al. (2004a)

 ˙^μ0,0=1iℏ[^He(¯R),^μ0,0]. (31)

[We have used Eq. (30)(a) and Eq. (26).] We also stress that, in this case, and that both and are proportional to .

To clarify the link between the new and the original formalism, it is helpful to write down the first order () EOM in terms of the CEID moments. This can be done by inverting Eqs. (30)(a-c,f) to express , , , and as functions of , , , and . [At the same CEID order and that .] The final result is:

 ˙^μ0,0 = 1iℏ[^He(¯R),^μ0,0]−1iℏ[^F(¯R),^μ1,0]+12iℏ[^K(¯R),^μ2,0], (32a) ˙^μ0,1 = {Δ^F(¯R),^μ0,0}−14{^K(¯R),^μ1,0}+1iℏ[^He(¯R),^μ0,1]−ia02b0[^K(¯R),^μ0,1]+ −1a20{Δ^F(¯R),^μ2,0}−b202a20M^μ1,0, (32b) ˙^μ1,0 = 12M^μ0,1+1iℏ[^He(¯R),^μ1,0]−ia02b0[^K(¯R),^μ1,0]+ia02b0[^F(¯R),^μ0,0]+ +a204b20{^K(¯R),^μ0,1}, (32c) ˙^μ2,0 = 1iℏ[^He(¯R),^μ2,0]+ia0b0[^F(¯R),^μ1,0]−ia0b0[^K(¯R),^μ2,0]+3iℏa208b20[^K(¯R),^μ0,0]+ +a202b20{Δ^F(¯R),^μ0,1}. (32d)

These equations might be compared with the Eq. (8) of Ref. Horsfield et al., 2005 (the mean-field second moment approximation) keeping in mind that in that paper the following ansatz as been made: , , and , where , , and are time-dependent quantities. According to our initial conditions (see Sec. II), the initial values of these variables are: , and . In order to find the EOM for , one can trace Eq. (32)(d) and it turns out that, to the first order in the coupling constant , . Remarkably, by assuming that the matrix is proportional to the unit matrix and by substituting into Eq. (32)(a-c), we obtain EOM for , , and which are equal to the ones stated in Eq. (8) of Ref. Horsfield et al., 2005 up to the first order in the coupling constant. Although there is no reason to believe that this agreement must be restricted only to a given set of initial conditions, it is not clear yet how it might be proved right for higher CEID order or different basis set expansion.

### iv.2 Energy conservation

By using the original CEID scheme, it is possible to write the total energy, Eq. (18), in terms of the density matrix moments. Horsfield et al. (2004a) A similar expansion is also obtained by computing Eq. (18) explicitly. [The matrix coefficients of the Hamiltonian can be found in appendix B.] During this computation, it is quite useful to distinguish between atomic kinetic energy, , and atomic potential energy, (see Eq. (1)). In terms of the CEID moments (see Eq. (30)), those two quantities are given by:

 Ekin = ¯P22M+¯PMTr{^μ0,1}+12MTre{^μ0,2}, (33a) Epot = Tre{^He(¯R)^μ0,0}−Tre{^F(¯R)^μ1,0}+12Tre{^K(¯R)^μ2,0}, (33b)

As explained in appendix C, the CEID evolution of the bare averages defined in Eq. (33) do not give a conserved (bare) total energy, , although the error is negligible for large enough CEID order. foo (c) That is because CEID provides an approximation of the exact evolution (in the truncated Hilbert space) of the observables’ averages (see Eqs. 46 and C). On the other hand, the exact evolution (in the truncated Hilbert space) of every observable can be retrieved starting from the CEID EOM and then adding a correcting term whose general analytical expression is reported at the end of appendix C

The time-derivatives of the corrections for the bare atomic kinetic and potential energy — whose integrals must be added to Eq. (33)(a) and Eq. (33)(b), respectively — are reported below:

 ˙C(N)Ekin= +¯PM(N+1)Tr{Δ^F(¯R)^ρN,N}+ −ib04M(N+1)√N2Tr{Δ^F(¯R)(^ρN,N−1−^ρN−1,N)}+ +¯Pb204a0M2(N+1)√N2Tr{^ρN,N−1+^ρN−1,N}+ −¯Pa04M(N+1)√N2Tr{^K(¯R)(^ρN,N−1+^ρN−1,N)}, (34a) ˙C(N)Epot= −ia20¯F4b0(N+1)√N2Tr{^K(¯R)(^ρN,N−1−^ρN−1,N)}+ +ib04M(N+1)√N2Tr{^F(¯R)(^ρN,N−1−^ρN−1,N)}. (34b)

## V Comparison between exact integration and CEID

In this section we present the main results of this work. They were obtained by means of the two numerical algorithms described in the previous sections, namely the exact integration scheme of Sec. III and the CEID scheme of Sec. IV. We recall that our main goal is to attest the convergence of CEID (by increasing its order) and to verify that the converged results agree with the exact dynamics of the 2LS geometries introduced in Sec. II. That can be safely done by a direct comparison between CEID and exact integration of the time-dependent Schödinger equation and this will be the object of Sec. V.1 — which contains a discussion of the electronic observable dynamics — and Sec. V.2 — which contains a discussion of the atomic observable dynamics.

The agreement of CEID with exact integration is clearly a fundamental numerical achievement, but it does not directly help in the interpretation of the simulation findings. Further insights can be obtained by comparing CEID against analytical results derived through first-order time-dependent perturbation theory (see appendix D). This comparison is reported in Sec.V.3.

### v.1 Electronic observables

Here we present the results of the electronic dynamics. As initial condition, we always choose the atomic vibrational ground-state on the upper PES and than we let the system to evolve according to either the exact Schrödinger evolution or the CEID equations. The WT of the initial uncorrelated density matrix is: , where , according to the definition given in Sec. IV, is the WT (in the mobile reference frame) of the atomic vibrational ground-state (which is centered in and ) and

 ^ρe(0)=(0001) (35)

describes a pure excited electronic state in the non-adiabatic representation introduced in Sec. II.

The most informative electronic observables are the probabilities to find the system in the upper or lower electronic state. Those are obtained as the diagonal entries of the electronic density matrix (see Sec. IV.1) and we shall call them electronic populations.

In Fig. 1 we collect the numerical results for the unshifted case.

A sketch of the PES is plotted in the first panel. [In this particular kind of 2LS, the difference between adiabatic and non-adiabatic surfaces is negligible.] In Fig. 1(b) the exact evolution of electronic populations is reported. We stress that the oscillatory population transfer between the two PES clearly confirms the crude picture we guessed in Sec. II.

In Fig. 1(c) the time-evolution of the electronic populations from CEID is reported. Different CEID orders, namely , , and , have been considered. As we expected, the (which is equivalent to ED, see Sec. IV.1) does not display any electronic transition. On the other hand, the outcomes of higher order CEID simulations present oscillations of the electronic populations which are in almost perfect agreement which the exact integration results. In particular the first order CEID simulation is well converged (i.e. there are no visible differences between the and findings). On the other hand, this is not very surprising; the unshifted 2LS is the easiest case, since the symmetries involved ensure that the evolved state is well described by the simple ansatz stated in Eq. (3) (see also appendix D).

In Fig. 2 we show the results for the shifted case following the same scheme employed for the previous figure. For this kind of 2LS, adiabatic and non-adiabatic PES are qualitatively different, but since in Fig. 2(a) the difference can be appreciated only close to the crossing, we provided a magnified plot of that region in a small inset. Once again, almost perfect agreement is seen between a well converged CEID simulation (here for at least ) and the exact integration of the time-dependent Schrödinger equation, while at the level of ED (i.e. ) the system is stuck in the upper PES. The fact that a first order CEID simulation is not yet well converged is not surprising and is a confirmation of the general trend predicted in Sec. II: the larger the surface displacement, , the higher will be the CEID order required to obtain a well converged simulation.

We see in Fig. 2(b) that the period of the electronic oscillations is larger for the shifted case than for the unshifted 2LS. [A perturbative account of that effect can be found in appendix D.] Moreover, in this shifted case the population exchange between the two PES is not complete: the minimum of the electronic population on the upper surface (corresponding to the maximum of the electronic population on the lower surface) is not exactly zero (one). This interesting feature is clearly visible in Fig. 2(b) and is also found in the two well converged CEID simulations in Fig. 2(c) so it is not a numerical feature. This is instead a non-trivial fingerprint of otherwise elementary dynamics which is caused by the virtual transitions — a clear quantum effect — between the low lying resonant states and more energetic atomic vibrational states. Further details can be found in Sec. V.3, in which we study the dependence of such residual population on the coupling constant .

### v.2 Atomic observables

Since the atomic motion is actually non-classical, we expect to find quantum fluctuations of the atomic observables around their average values. We study the time-evolution of the average atomic position and momentum, and , because they provide a sort of effective trajectory in the phase space which represent an important link with classical MD. Obviously, the concept of trajectory is not well defined in quantum mechanics and is useful as an approximation only if the fluctuations are not too large. So, we also consider the variance of the atomic position, , in order to test the accuracy of CEID in describing possibly non-classical atomic dynamics.

We start by reporting results for the average atomic position and momentum, and . They are evolved by means of the Hamilton-Ehrenfest equations (see Sec. IV) according to CEID while in the exact integration scheme they are obtained by means of Eq. (15). For the unshifted 2LS, and for all time due to the inversion symmetry displayed by the system. On the other hand, the findings reported in Fig. 3 for the shifted case once again show almost perfect agreement between CEID and exact integration in a completely non-trivial case. In particular, CEID not only reproduce the general trend of both and (large period oscillations), but also gives the short time scale details (rapid oscillations).

In Fig. 4 we report the results for the variance of the atomic position, , for both the unshifted and shifted 2LS. This observable can been obtained as the trace of the CEID moment defined in Sec IV.1. Once again, almost perfect agreement has been found between a well converged CEID simulations (here and for the unshifted and shifted 2LS, respectively) and the exact integration of the time-dependent Schrödinger equation. We stress that such fluctuations are quite significant and so the atomic dynamics is only poorly approximated by its average trajectory in the classical phase space. CEID is working properly even in those highly non-classical cases.

Finally, we have also verified the agreement between CEID and exact integration for the other entries of the covariance matrix, namely and . However, numerical findings for those cases are nor reported here because they are not qualitatively different from the case.

### v.3 Comparison between time-dependent perturbation theory and CEID

In this last section we briefly compare the CEID outcomes against time-dependent perturbation theory results. [Mathematical details are collected in appendix D.]

First of all, from a well converged CEID simulation (e.g. and for the unshifted and shifted case, respectively), the values of the electronic oscillation frequency and the residual electronic population can be obtained by means of a straightforward numerical interpolation. Then, this procedure can be repeated for the same 2LS geometries, but different electron-ion coupling constant, (see Eq. (1)). It is instructive to study the effect of the atomic motion on the electronic transitions because it might cause non-classical phenomena, like quantum interference between different transition paths. Those effects are usually hard to interpret without a model which can — at least qualitatively — describe the physics involved. Fortunately, for the kind of model 2LS we considered in this paper, a simple model can be obtained by means of time-dependent perturbation theory, whose prediction for the electronic transition frequency and the residual population are summarized in Eq. (52) and Eq. (55), respectively.

In Fig. 5(a) numerical values of electronic population frequencies are reported against several values of the coupling constant, . A clear linear trend is manifest in all the 2LS geometries. Moreover, numerical values are in almost perfect agreement with the analytical results, Eq. (52).

In Fig. 5(b) the residual populations are plotted against the same coupling constant values. Although the analytical trend, (see Eq. (55)) is confirmed, in general is not easy to give an estimate of the prefactor, . On the other hand, for the unshifted 2LS case, only one term in Eq. (55) is non-zero due to the SHO selection rules. As a consequence, a numerical estimate can be obtained and it gives , while a direct numerical interpolation gives . We stress that such disagreement might depend on the kind of approximation we made in order to derive Eq. (55) and it does not effect the general scaling trend of the residual population with the coupling constant.

## Vi Discussion and conclusions

We have presented a new formulation of correlated electron-ion dynamics (CEID). It is based on a suitable expansion of the quantum fluctuations around the mean-field atomic trajectories and its lowest accuracy limit has been proved to be equivalent to the well-known Ehrenfest dynamics (ED). This new formulation has been obtained by a combined use of: 1) an expansion of the density matrix in terms of atomic harmonic states centered around the average instantaneous atomic positions; 2) an exact Wigner transform with respect to the atomic degrees of freedom of the expanded density matrix. The validity of this scheme has been successfully tested by simulating the non-adiabatic time-evolution of a model two level system (2LS). The accuracy of our simulations is determined by a single parameter which is related to the order of the density matrix expansion and is called the CEID order. We then verified that, for all the considered 2LS geometries, the exact quantum dynamics — obtained by exact integration of the time-dependent Schrödinger equation — is eventually retrieved by increasing the CEID order. We think that this is a crucial property of our new CEID scheme which allows us to estimate the convergence of a numerical simulation even when reliable benchmarks are not available.

As for the other proposed CEID schemes, Horsfield et al. (2004a); Horsfield et al. (2005) our algorithm only needs the Hamiltonian and the initial conditions to start and the subsequent evolution is computed smoothly, without resorting to any kind of surface hopping or wave-function spawning. No a posteriori position, velocity, or density matrix adjustment is needed. The exact evolution (in the truncated Hilbert space) of every observable average can be obtained starting from the CEID EOM by adding a correction term (whose analytical expression is known) which is anyhow negligibly small for large CEID order. Moreover — and at variance with other available algorithms Müller and Stock (1997); Tully (1998) — it works perfectly well within a non-adiabatic representation of the electronic PES. This is desirable because non-adiabatic PES may be smoother than the adiabatic ones Ben-Nun et al. (2000) and also because a costly diagonalization of the atomic potential energy at each step is avoided. All these dynamical properties make our new CEID algorithm a good candidate for simulating atomic systems in which quantum coherence is relevant.

The advantages of a coherent quantum scheme might be relevant even when macroscopic quantum coherence is not shown. It is well known that ED — at variance with other quantum-classical methods Parandekar and Tully (2005); Kapral (2006) — cannot thermalize a mixed electron-ion system (the electronic degrees of freedom are too hot with respect to the atomic ones. Parandekar and Tully (2005)) This failure of the mean-field approximation depends on the absence of quantum fluctuations which cause spontaneous phonon emission from an excited electronic state. [This drawback is apparent also in the 2LS simulations considered in this paper (see Sec. V): the ED is always stuck in the initial excited state.] As a consequence, ED does not satisfy microreversibility. Müller and Stock (1997); Tully (1998) On the other hand, a CEID simulation beyond ED can describe quantum fluctuations and meet the coherence requirements for microreversibility in a very natural way. [Once again, see the results of Sec. V.]

Although it is not the main concern of this paper, our group is considering a viable way to approach quantum thermalization physics by means of CEID. A first possibility is given by the spin-boson model Leggett et al. (1987) in which the bath degrees of freedom are treated explicitly by means of a collection of many quantum harmonic oscillators. On a more speculative ground, one can think to implement the generalization of the Nosè-Hoover thermostat introduced in Ref. Grilli and Tosatti, 1989. This scheme is known to fail for ED Mauri et al. (1993) due to the lack of correct quantum back-reaction on the classical bath variables. Sergi (2005) On the other hand, as we have shown again in this paper, CEID corrects this ED drawback and it might be better suited for that sort of thermostat. Moreover, a successful attempt to couple the Nosè-Hoover thermostat to the spin-boson model is known in literature Sergi (2007a) and it can provide an interesting test case for future CEID simulations.

Our CEID algorithm is computationally demanding and is expected — in the worst case scenario — to scale as , where is the CEID order and is the number of atomic coordinates. foo (d) Nevertheless, it must be pointed out that the number of relevant atomic coordinates can be effectively much smaller than . Ness and Fisher (1999); Tamura et al. (2007) In this case, one might accelerate a CEID simulation by allowing for quantum atomic fluctuations along the relevant directions only. We also stress that the CEID algorithm is still faster than the exact integration scheme employed to produce benchmark calculations in this paper (see Sec. III) which should scale as since a numerical diagonalization of the Hamiltonian in the truncated Hilbert space is implied. We are also considering alternative truncations of the Hilbert space in order to restore the polynomial scaling with atomic degrees of freedom of the early CEID algorithms.

We see another possible advantage of this CEID scheme over exact integration: the former expands the quantum fluctuations around mean-field atomic trajectories, while the latter expands with respect to a fixed reference frame. Now, consider a quantum motion in which there are fluctuations about the mean-field atomic trajectories that are very tightly confined along a given direction. With our CEID formulation such fluctuations can be treated accurately with a low order expansion. However, schemes that employ basis functions which are not centered around the atomic trajectories could require a very high order expansion to reproduce that confined behavior if the trajectories are remote from the center of the basis functions.

Other algorithms, such as molecular dynamics with quantum transitions or ab initio multiple spawning, might have a lower computational complexity, especially if the region of the configuration space where non-adiabatic effects are relevant is small and crossed only few times during the time-evolution. This is the case, for instance, for many chemical reactions in a gaseous or diluted phase. On the other hand, we recall here that CEID was explicitly devised to deal with electron-ion correlations in metals, a kind of systems in which the aforementioned algorithms are expected to be less efficient.

Needless to say, a reliable algorithm to simulate microscopic electro-mechanical effects, including joule heating, will find important application in nanostructure design. Our CEID algorithm is a good candidate because its accuracy can be systematically increased by tuning a single parameter that allows us to approximate the quantum atomic fluctuations in a physically transparent way. Moreover, since quantum coherence is well addressed by CEID, subtle photo-physical effect like luminescence in conjugated polymers might be addressed by this method. Applications of our algorithm to larger atomic systems and different thermodynamical ensembles are subjects of ongoing study.

###### Acknowledgements.
LS is supported by EPSRC under grant EP/C524381/1 and MM is supported by EPSRC under grant GR/S80165. The authors like to thank Tchavdar Todorov for illuminating suggestions and critically reading of this paper. MM also thanks A.T. Paxton for helpful comments and LS acknowledges useful discussions with R. Peixoto Miranda, D.B. Bowler, P. Delaney, and A.M. Stoneham.

## Appendix A Properties of the Pn,m(R,P) functions

All the properties of the functions introduced in Sec.