{'fraction_non_alphanumeric': 0.038503549972692515, 'fraction_numerical': 0.009721463681048608, 'mean_word_length': 4.39178445229682, 'pattern_counts': {'":': 0, '<': 1, '<?xml version=': 0, '>': 2, 'https://': 0, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 5, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
Robust non-adiabatic molecular dynamics for metals and insulators

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


I. INTRODUCTION
Ordinary molecular dynamics (MD) 1 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), 2,3 molecular dynamics with quantum transitions (MDQT), 2,4 , mixed quantum-classical dynamics (MQCD), 5,6 and ab initio multiple spawning (AIMS) 7,8 . They have been employed to simulate quantum dissipative dynamics 9,10 as well as non-adiabatic processes such as photo-chemical reactions, 11,12 polaron formation in conjugated polymers, 13 , and proton transfer in solutions 4 . Although less successful in practical applications, we also mention other dynamical schemes that provide valuable theoretical insights: the frozen Gaussian approximation, 14 which is one of the pillars of AIMS, and the Gauss-Hermite wave-packet expansion, 15 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) 3,16,17 , has been introduced. This method has been mainly applied to study the heat production and dissipation in model metallic nanostructures, 3,18,19 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. 2,8 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 20,21 and conserves the kinematic constraints. On the other hand, timetranslation invariance is valid only approximately in MQCD 20 and numerical instabilities have so far limited the surface-hopping implementation of MQCD to relatively short time simulations. 10 Finally, ED -although it evolves smoothly -has been shown to poorly describe the atomic heating caused by electron-ion interaction. 22 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 23 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. 3 So far, a few schemes to truncate the hierarchy have been proposed 16 in order to simulate the CEID EOM and currently available algorithms are restricted to a meanfield second moment approximation. 17,24 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 25,26 -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 lowdimensional model. 27,28 In general, the Hamiltonian for a system made by electrons and ions (in the absence of external fields) can be written as follows: whereR andP are the quantum operators for the atomic position and momentum, while the electronic dependence of H is collected intoĤ e . 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,Ĥ e , but in this paper we focus only on the following parametrization: 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Ĥ e (R). Eq. (1) and (2) depend on the following parameters: the atomic mass M, the harmonic constant K, the electron-ion coupling constant f c , the surface displacement R 0 , 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 R experiences quantum oscillations and so the system cannot be exactly localized in the minimum of the PES (R = 0 according to our parametrization). Since this initial state is not an eigenstate of the interacting Hamiltonian (i.e. for f c = 0), 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. 22 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 29 describes a 2LS coupled to an environment made by a collections of many quantum harmonic oscillators. This model can be effectively simulated 9,30 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. 22 In order to be as transparent as possible when discussing the dynamical features of our For our purposes, not all the parameters need to be varied during the numerical experiments. First of all, we fixed both the atomic mass (M = M 0 ) and the harmonic constant (K = M 0 ω 2 ) 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. f c ≤ 0.1 ω/a 0 . Finally, we report in this paper the numerical results for two different 2LS geometries only: the unshifted 2LS, with R 0 = 0, and the shifted 2LS, with R 0 = a 0 . 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 R 0 or with other values of K (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 t might be approximated as: where c 0 and c 1 are time-dependent complex coefficients. This state would give a distribution of the atomic position R more or less localized around each of the two classical equilibrium positions, namely R = 0 for the upper PES and R = R 0 for the lower PES.
Eq. (3) might not be a good ansatz for the evolved state as soon as the displacement R 0 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 R = 0 and R = R 0 . Therefore, we will produce a manifest physical inconsistency if we assume that the generic 2LS state can be always well 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. 31 as where R = (R 1 , . . . , R n ), P = (P 1 , . . . , P n ), and the time-dependence has been suppressed.
The latter form in Eq. (4) is more common nowadays. In a straightforward extension of this definition, for an operatorÂ, expandable in basis states {φ a },Â = ab φ a (x)A ab φ * b (x), the WT is given by 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 {|A } and {|a }.
The operator is assumed to be expandable asĈ = ABab |a ⊗ |A C ABab b| ⊗ B|, and its partial WT with respect to the |A -subsystem iŝ where, with Φ A the position representation of |A , Note the important fact thatĈ w,A is still an operator in the |a -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Ĉ w instead ofĈ w,A , and also refer to the |A -subsystem as the atomic and the |a -subsystem as the electronic subsystems.
The time-evolution of the system is generated by a Hamiltonian, whose eigenvalues and eigenvectors shall be E n and |Ψ n , respectively. The basis states {|a } and {|A } of the subsystems are taken to be time-independent from now on. We can expand an eigenstate in the product basis |Ψ n = Aa C Aa n |a ⊗ |A , or vice versa, |a ⊗ |A = n D n Aa |Ψ n , where we have the relations Aa C Aa n D m Aa = δ m n and n C Aa n D n Bb = δ Aa Bb . The expansion coefficients C and D are time-dependent, in any particular case the numerical diagonalization will provide us with the E n and the coefficients C Aa n (0), which make up the eigenvectors, i.e. C Aa n (0) is the Aa-component of eigenvector n.
We now can expand an operatorĜ in the eigenstates or the product states: with The partial WT with respect to the atomic subsystem is then where in turn and 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 ing Eq. (10), the factoring Eq. (12), and Eq. (13), all the components of G ab w can be computed as functions of R, P, t (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 t. 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 N. Then, the diagonalization has to be carried out using a number of product basis states M > N such that the N lowest eigenvalues and corresponding eigenvectors are well converged; typically, we used N ≈ 120 and M ≈ 5N. 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 N ′ (N ′ < N) to states with index larger than N ′ . These first N ′ states, in case of the occupations, and of 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. 32 We start from the well known quantum Liouville equation: which is the EOM for the density matrixρ of the system. 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). 5 By using the standard Dirac's bra and ket notation, this expansion can be expressed as: where the functions {φ n } are a COS in the atomic subspace. As a consequence, in Eq. ( For instance, the total energy of the system is given by: Tr e Ĥ m,nρn,m . [The two traces, Tr and Tr e , apply to different linear spaces, namely the whole Hilbert space and the electronic subspace: see Sec. III.] As a further example, the EOM forρ n,m are obtained by plugging Eq. (17) into Eq. (16): It is easy to prove that E tot 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 N, 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Â, A w (R, P ), is still an operator in the electronic subspace but it explicitly depends on what are now the classical atomic position R and momentum P . In the context of non-adiabatic MD, similar partial WT have been already considered 5, 33 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: where ⋆ is the non-commutative Moyal product: 34,35 [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 |φ n φ m |.
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: It can be shown that from Eq. (22) It is also reasonable to take the phase-space trajectory (R(t), P (t)) as a zero-order approx- This fact suggests that instead of taking the origin as a reference point in the phase-space (R, P ) (i.e a fixed reference frame), one can take advantage of Eq. (23) and use (R(t), P (t)) as reference (i.e. a mobile reference frame). After this mobile reference frame transform, the EOM becomes: Although it is not apparent from Eq. At this stage two paths can be followed. In the first case, the Moyal product is expanded in (usually up to first order 5,33 ) 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, 36,37 and a generalized non-Hamiltonian bracket should be introduced instead. 21  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 T w can be defined as follows: where P n,m (∆R, ∆P ) is the WT of |φ n φ m | 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 P n,m is known, 34 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ρ n,m (in the mobile reference frame): 4i 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 + We now briefly discuss the link between our new formulation and the original CEID.

A. 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.  dRdP ∆R n ∆P mρ w (R, P, t) .
By using the property of the WT, it is easy to find a link between the new and the original notation:μ 0,0 =ρ e ,μ 0,1 =λ 1 , andμ 1,0 =μ 1 . Similar relations for higher CEID moments can be stated, but some extra attention must be payed in the derivation due to the nontrivial commutation relations between positions and momenta. It is worth noting that CEID moments provide valuable information about the system. For instance, the quantities where A n,m r,s = 1 2π dRdP ∆R n ∆P m P r,s (R, P ) .
As usual, a set of recurrence relations for A n,m r,s 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: n(n − 1)ρ n,n−2 − (2n + 1)ρ n,n + n(n − 1)ρ n−2,n , (30d) n(n − 1)ρ n,n−2 + (2n + 1)ρ n,n + n(n − 1)ρ n−2,n .
[At the same CEID orderμ 1,1 = 0 and thatμ 0,2 = (b 2 0 /a 2 0 )μ 2,0 .] The final result is: As explained in appendix C, the CEID evolution of the bare averages defined in Eq. (33) do not give a conserved (bare) total energy, E tot = E kin +E pot , although the error is negligible for large enough CEID order. 38 That is because CEID provides an approximation of the exact evolution (in the truncated Hilbert space) of the observables' averages (see Eqs. C2 and C3).
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 energywhose integrals must be added to Eq. (33)(a) and Eq. (33)(b), respectively -are reported below:Ċ

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 timedependent Schödinger equation and this will be the object of Sec. V A -which contains a discussion of the electronic observable dynamics -and Sec. V B -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 firstorder time-dependent perturbation theory (see appendix D). This comparison is reported in Sec.V C.

A. 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 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ρ e =μ 0,0 (see Sec. IV A) 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 exact integration results. In particular the first order CEID simulation is well converged (i.e. there are no visible differences between the N = 1 and N = 5 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)  This interesting feature is clearly visible in Fig. 2(b) and is also found in the two well converged CEID simulations in Fig. 2 In Fig. 4 we report the results for the variance of the atomic position, ∆R 2 , for both the unshifted and shifted 2LS. This observable can been obtained as the trace of the CEID momentμ 2,0 defined in Sec IV A. Once again, almost perfect agreement has been found between a well converged CEID simulations (here N = 5 and N = 10 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 ∆P 2 and ∆R∆P . However, numerical findings for those cases are nor reported here because they are not qualitatively different from the ∆R 2 case. 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, 16,17 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 2,11 -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 7 and also because a costly diagonalization of the atomic potential energy H e (R) 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 6,39 -cannot thermalize a mixed electron-ion system (the electronic degrees of freedom are too hot with respect to the atomic ones. 39  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 29 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. 40. This scheme is known to fail for ED 41 due to the lack of correct quantum backreaction on the classical bath variables. 21 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 10 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 (N +1) 2Nc , where N is the CEID order and N c is the number of atomic coordinates. 42 Nevertheless, it must be pointed out that the number of relevant atomic coordinates can be effectively much smaller than N c . 43,44 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 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. The P n,m are orthonormal: P n ′ ,m ′ (R, P ), P n,m (R, P ) ≡ 1 2 π dRdP P m ′ ,n ′ (R, P ) P n,m (R, P ) = δ n,n ′ δ m,m ′ , (A1) and form a complete set in the Wigner space. [Also note thatP n,m = P m,n .] The following identities are useful when computing the action of a canonical operator on P n,m : R P n,m (R, P ) = + a 0 2 n 2 P n−1,m + n + 1 2 P n+1,m + m 2 P n,m−1 + m + 1 2 P n,m+1 , Higher order identities can be recursively obtained by using these basic identities.

APPENDIX B: DERIVATION OF THE EQUATIONS OF MOTION
In this appendix we show how Eq. (26) can be derived from Eq. (24) by using the expansion in Eq. (25). First of all, we want to find the action of the position and momentum derivatives on the matrix coefficientρ n,m . We formally introduce the operators D R and D P in the following way: Therefore, by means of the identities reported in Eq. (A2)(b,d), we find that: Finally, previous equations can be also written in a more compact form as follow: For instance, for a quadratic Hamiltonian expanded with respect to the mobile reference one finds that:

4K
R m(m − 1)δ m−2,n + (2m + 1)δ m,n + n(n − 1)δ m,n−2 . (B7) [The first three terms account for the the atomic kinetic energy and the last three for the atomic potential energy.] From Eq. (B7) one can easily obtainĤ n,m whose expression must be substituted in Eq. (B4) to obtain Eq. (26). The requirement of neglecting in the RHS of Eq. (26) all those matrix coefficients whose indices are greater than the CEID order can be directly enforced by constraining the summation index k in Eq. (B4) to be at most equal to the CEID order.
This expression is still true in the mobile reference frame since the extra terms which arise from the implicit time-dependency ofÂ andρ cancel out exactly. The truncated version of Eq. (C1) is naturally given by: where T is the anti-Wigner transform of the truncation operator T w introduced in Sec. IV.
[We recall that this operator depends on the CEID order.] As expected, averages of those observables which commute with the Hamiltonian are constant of motion.
It turns out that, by using the CEID EOM (see Eq. (26)), the bare dynamics of the averages (in the mobile reference frame) is given by: the Eqs. C2 and C3 do not provide the same dynamics, the difference betweenĀ (T ) (t) and A bare (t) -the correcting term CĀ(t) -is expect to be small for large enough CEID order.
This fact has been verified for the atomic kinetic and potential energy, whose analytical expressions of the respective correcting terms are given in Sec. IV B.
Finally, we stress that a general expression for the time-derivative of the correcting term to the average of the observableÂ can be derived analytically. It is given by: O T also depends on the CEID order.] By using this expression, one can implement the correct truncated dynamics of the observable averages (i.e. Eq. (C2)) starting from the CEID EOM and then adding the numerical integral of the appropriate correcting term.

ORY OF A RESONANT 2LS
We start from the non-interacting limit of Eq. (1), i.e. by setting f c = 0 in Eq. (2).
In this easy case, the two PES are dynamically uncoupled. On each surface the atomic configurations can be classified by using the appropriate SHO basis set. We indicate by n+1 are also degenerate, having energy E n+1 = (n + 3/2) ω. The ground-state is not degenerate and its energy is E 0 = 1/2 ω.] This degeneracy is lifted by switching on the interaction (f c > 0). The hybridization energy is: ∆ n+1 = χ (u) n |Ĥ|χ (l) n+1 . As a consequence, a system prepared in the state |χ From Eq. (D2), it is easy to derive the dynamics of the upper electronic population: i |Ψ (0) (t) | 2 = cos 2 |∆ 1 |t .
The correspondent oscillation frequency is given by: We stress that at half period (t = T /2 = π/ω p ) the upper electronic population is exactly zero, i.e. at zero order there is no residual population (see Sec. V A).
To go beyond the zero order, we use the following transformation: |Ψ(t) = e 1 i Ĥ 0 t |φ(t) which leads to an effective Schrödinger equation for |φ : (D7)