High-order geometric integrators for representation-free Ehrenfest dynamics

Ehrenfest dynamics is a useful approximation for ab initio mixed quantum-classical molecular dynamics that can treat electronically nonadiabatic effects. Although a severe approximation to the exact solution of the molecular time-dependent Schr\"odinger equation, Ehrenfest dynamics is symplectic, time-reversible, and conserves exactly the total molecular energy as well as the norm of the electronic wavefunction. Here, we surpass apparent complications due to the coupling of classical nuclear and quantum electronic motions and present efficient geometric integrators for"representation-free"Ehrenfest dynamics, which do not rely on a diabatic or adiabatic representation of electronic states and are of arbitrary even orders of accuracy in the time step. These numerical integrators, obtained by symmetrically composing the second-order splitting method and exactly solving the kinetic and potential propagation steps, are norm-conserving, symplectic, and time-reversible regardless of the time step used. Using a nonadiabatic simulation in the region of a conical intersection as an example, we demonstrate that these integrators preserve the geometric properties exactly and, if highly accurate solutions are desired, can be even more efficient than the most popular non-geometric integrators.


I. INTRODUCTION
Mixed quantum-classical methods, such as the surface hopping, 1-5 mean-field Ehrenfest dynamics, [6][7][8][9][10][11][12][13][14][15] and methods based on the mixed quantum-classical Liouville equation [16][17][18] or the Meyer-Miller-Stock-Thoss mapping Hamiltonian, [19][20][21][22][23][24][25] remedy one of the shortcomings of classical molecular dynamics: its inability to describe electronically nonadiabatic processes [26][27][28] involving significantly coupled [29][30][31] states. Although a severe approximation to the exact quantum solution, 8,32,33 Ehrenfest dynamics can provide a useful first picture of nonadiabatic dynamics in some, especially strongly coupled systems. Indeed, Ehrenfest dynamics was successfully used to describe electron transfer, [34][35][36][37][38] nonadiabatic processes at metal surfaces, [39][40][41][42] and photochemical processes. [43][44][45] The mean-field theory was also employed to simplify the evaluation of the memory kernel in the generalized master equation formalism. 46 In addition, Ehrenfest dynamics provides a starting point for various refined methods. For example, a multi-trajectory, locally mean-field generalization of Ehrenfest dynamics was used to evaluate vibronic spectra 47 and, when combined with the semiclassical initial value representation, can describe even wavepacket splitting. 48 A further generalization, the multiconfigurational Ehrenfest method, [49][50][51] includes correlations between Ehrenfest trajectories. In what follows, we shall only consider the basic, mean-field Ehrenfest method, whose validity conditions were formulated by Bornemann et al. 52 The coupling between nuclear and electronic dynamics complicates the numerical integration in Ehrenfest dynamics. The widely-used two and three time step methods [11][12][13]53,54 improve the efficiency by using different integration time steps that account for the different time scales of nuclear and electronic motions (see Appendix A). However, such integration schemes violate the geometric properties of the exact solution: the simpler, two time step method is irreversible and neither method is symplectic (see Fig. 6 in Appendix A).
Almost every geometric property [55][56][57] can, however, be preserved exactly by employing the symplectic integrators 58 based on the splitting method. 59,60 This splitting method is widely applicable-so long as the Hamiltonian can be decomposed into exactly solvable parts-and was employed to obtain symplectic integrators in many well-known applications, including molecular quantum 61 and classical 62 dynamics, Schrödinger-Liouville-Ehrenfest dynamics, 63 and the Meyer-Miller-Stock-Thoss mapping approach. [64][65][66] In particular, because the Ehrenfest method in either the adiabatic or diabatic representation can be formu-lated as a special case of the mapping method, 19,67 the integrators developed for one of these two methods should also be applicable to the other. Motivated by on-the-fly ab initio applications that employ increasingly practical real-time time-dependent electronic structure methods, 68,69 here we present integrators that-in contrast to the integrators 64-66 formulated in the mapping approach-do not rely on any particular representation of electronic states and thus avoid the expensive construction of a truncated diabatic or adiabatic electronic basis.
Typically, to reach the same accuracy, geometric integrators need greater computational effort than their non-geometric counterpart. 55,56 Yet, the efficiency of geometric integrators can be improved significantly by employing various composition methods. 55,56,[70][71][72][73][74] Thus obtained integrators of high orders of convergence in the time step offer the best of both worlds: they are efficient while conserving the relevant geometric structure exactly. 56 After showing analytically, in Sec. II, that the high-order geometric integrators preserve almost all of the geometric properties of Ehrenfest dynamics, in Sec. III, we numerically demonstrate the efficiency and geometric properties of these integrators on a fourdimensional extension 75-77 of the Shin-Metiu model. 78,79 In this system, the first and second excited adiabatic states are coupled significantly due to a conical intersection. Section IV concludes the paper.

II. THEORY
A. Time-dependent Hartree approximation for the molecular wavefunction Quantum evolution of a molecule is governed by the time-dependent Schrödinger equation where Ψ t denotes the molecular state at time t and H is the molecular Hamiltonian. In general, we will denote operators acting on both nuclei and electrons by a calligraphic font, whereas the operators acting either only on nuclei or only on electrons will have a hat. The molecular Hamiltonian is equal to the sum of the nuclear kinetic energy operator electronic kinetic energy operatorT and potential energy operator V(q,Q). We assume that the nuclear position Q and momentum P are D-dimensional vectors, whereas the electronic position q and momentum p are d-dimensional vectors. The nuclear and electronic mass matrices, M and m, can, in general, be real symmetric D × D and d × d matrices, respectively.
The time-dependent Hartree (TDH) 57,80-84 approximation is an optimal approximate solution to the molecular TDSE (1) among those in which the molecular state can be written as the Hartree product of the nuclear wavepacket χ t and electronic wavepacket ψ t ; the complex number a t is inserted for convenience. In the TDH approximation, obtained by applying the Dirac-Frenkel timedependent variational principle 57,85,86 to ansatz (5), the prefactor evolves as 57 and the nuclear and electronic states satisfy the system of coupled nonlinear Schrödinger equations with mean-field nuclear and electronic Hamilto- The mean-field operators satisfy the obvious identity Ĥ nu χt = Ĥ el ψt = E. Note that the solution expressed by Eqs. (6)-(8) is unique except for an obvious gauge freedom in redistributing the phase among a t , χ t , and ψ t .
B. Mixed quantum-classical limit: Ehrenfest dynamics In the classical limit for nuclei, the nuclear position and momentum operatorsQ andP are replaced with classical variables Q and P . Then, the mean-field nuclear Hamiltonian (9) is no longer an operator, but a phase space function where H(Q, P ) =T el + T nu (P ) +V (Q), and the mean-field electronic Hamiltonian (10) becomes the molecular Hamiltonian evaluated at the current nuclear positions and momenta: We thus obtain the mixed quantum-classical Ehrenfest dynamics, in which the nuclear positions and momenta evolve according to classical Hamilton's equations of motion with Hamiltonian H nu (Q, P ), and the electronic state evolves according to the TDSE with a time-dependent HamiltonianĤ el (Q t , P t ): Note that these three differential equations are coupled and, moreover, that the electronic TDSE is nonlinear due to this coupling.
Equations (13)- (15) can be re-expressed as more compact Hamilton's equationṡ associated with an effective mixed quantum-classical Hamiltonian acting on an extended, effective mixed quantum-classical phase space with coordinates x eff = (q eff , p eff ) = (Q, q ψ , P, p ψ ). The "quantum" Darboux coordinates (q ψ , p ψ ) consist of the real and imaginary part of the electronic wavefunction in position representation: q ψ := √ 2 Reψ(q) and p ψ := √ 2 Imψ(q) (we omit the dependence of q ψ and p ψ on q for brevity).
In practical calculations, q ψ and p ψ are usually represented in a finite basis or on a grid as N -dimensional vectors, where N is either the size of the basis or number of grid points.
In such cases, functional derivatives (21) and (22) reduce to partial derivatives where H el (Q, P ) is an N × N matrix representation of operatorĤ el (Q, P ).

Norm conservation
Ehrenfest dynamics conserves the norm of the electronic wavefunction because where we used Eq. (15) and the hermiticity ofĤ el (Q t , P t ).

Energy conservation
The total energy E = H nu (Q t , P t ) = Ĥ el (Q t , P t ) ψt of the system is conserved, in general, by the time-dependent variational principle and, in particular, by the TDH approximation.
However, because we have also taken the mixed quantum-classical limit, let us verify the conservation of energy explicitly: where we used the hermiticity ofĤ el (Q t , P t ) and Eqs. (13)- (15). The energy conservation also follows directly from the effective Hamiltonian structure: where we used Eqs. (16) and (17).

Time reversibility
An involution is a mapping S that is its own inverse, i.e., S(S(x)) = x. We will consider the involution that changes the sign of the nuclear momenta and conjugates the electronic wavefunction in position representation (i.e., changes the sign of p ψ ). Following Ref. 56 Because effective Hamiltonian H eff is an even function in p eff , i.e., H eff (x eff ) = H eff (Sx eff ), its Hamiltonian flow satisfies 56 Since S −1 = S and, by definition, any flow is symmetric (i.e., Φ −t = Φ −1 t ), 55,56 the satisfaction of Eq. (32) implies the satisfaction of the time reversibility condition (31).

D. Geometric Integrators
As in the split-operator algorithm 61,90-92 for the TDSE or in the Verlet algorithm 62 for Hamilton's equations of motion, we can obtain a symmetric potential-kinetic-potential (VTV) algorithm of the second order in the time step ∆t by using the Strang splitting 59 and performing, in sequence, potential propagation for time ∆t/2, kinetic propagation for ∆t, and potential propagation for ∆t/2. The second-order kinetic-potential-kinetic (TVT) algorithm is obtained similarly, by exchanging the potential and kinetic propagations. Either of the two second-order algorithms can be symmetrically composed 55,93 to obtain an algorithm of an arbitrary even order of accuracy in ∆t. This is achieved by using the recursive triplejump 70,71 or Suzuki-fractal 71 composition schemes, or with a more efficient scheme specific to each order (which we shall call "optimal"). 71,73,74 We, therefore, only need to present the analytical solutions of the kinetic and potential propagation steps for arbitrary times t.
During the kinetic propagation, the Hamiltonian reduces tô and the equations of motion (13)-(15) becomė which are equivalent to Hamilton's equations (16) and (17) with Because nuclear momenta P t do not evolve during the kinetic propagation, Eqs. (34)- (36) can be solved analytically to obtain AsT el = T (p), Eq. (39) is easily evaluated in momentum representation.
During the potential propagation, the Hamiltonian reduces tô and the equations of motion (13)-(15) becomė which are equivalent to Hamilton's equations (16) and (17) with H eff = V (Q) ψ . Because nuclear positions Q t do not evolve during the potential propagation, one can replace Q t with Q 0 in Eqs. (42) and (43). Even after the substitution of Q t = Q 0 , Eq. (42) seems hard to solve due to an apparent coupling to Eq. (43). However, this coupling can be removed by noting that As a result, Eqs. (41)- (43) can be solved analytically to obtain (47) is easily evaluated in position representation.

III. NUMERICAL EXAMPLE
We use a low-dimensional model that is solvable "numerically" exactly to demonstrate the geometric and convergence properties of the presented integrators. Since the high efficiency and geometric properties of these integrators are most meaningful when the mean-field Ehrenfest approximation is valid, we have chosen the system and initial state carefully so that numerically converged Ehrenfest and exact quantum simulations yield similar results.
where Q a = (−L/2, 0) and Q b = (L/2, 0) are the positions of the two fixed protons, and are attractive and repulsive regularized Coulomb potentials; following Ref. 75, we take L = with Q c = 3.5 a.u. ensures that the system remains bound.
For the dynamics, we considered the initial state where ϕ i is the ith excited adiabatic electronic state, and is the ground vibrational eigenstate of a harmonic fit to the ground electronic state; here, σ = 0.24 a.u. The displacement of the initial wavepacket by Q 0 = (0.5 a.u., 1.5 a.u.) from the ground state equilibrium is motivated by the displaced excitation of molecules: Suppose state ϕ 2 is dark; then a wavepacket may reach it at a nuclear geometry that is not the ground state equilibrium via an intersection with a bright state. To obtain ϕ 2 (q; Q), we solved the electronic time-independent Schrödinger equation whereT el = T (p) andV (Q) = V (q, Q) are operators acting on electrons; in position representation, Eq. (53) takes a more familiar form Section S1 of the supplementary material describes the method we employed to solve this equation.
Because our approach does not rely on a specific electronic basis (such as the basis of adiabatic or diabatic electronic states) to represent the molecular wavepacket, the initial state can be a general function of q and Q. To be specific, however, we chose to start the dynamics from a single excited adiabatic electronic state (here ϕ 2 ), which is the most common choice in the literature studying nonadiabatic dynamics following a light-induced excitation. 27,[94][95][96][97] In the model described by Eqs. (48)-(50), the second excited adiabatic state ϕ 2 is, indeed, significantly coupled to the first excited state ϕ 1 by a conical intersection depicted in Fig. 1.
In Fig. 2, we compare the exact quantum dynamics Ψ t = exp (−itH/ )Ψ 0 with Ehrenfest dynamics x eff,t = Φ H eff ,t (x eff,0 ). The initial state of the system is (Q t , P t , ψ t ) | t=0 = (Q 0 , 0, ϕ 2 (Q 0 )) and the corresponding initial mixed quantum-classical phase space point can be thought of as state (51) with an infinitesimally narrow Gaussian wavepacket; the fourth component in Eq. (55) is zero because the state ϕ 2 (q; Q 0 ), in position representation, is purely real: Imϕ 2 (q; Q 0 ) = 0. We compare three observables: nuclear position Q(t), adiabatic population P i (t), and electronic density ρ el (q, t). In quantum dynamics, they are obtained from the full wavefunction Ψ t as 75,77 To find population P i (t) from Eq. (57), we computed the expectation value of the population operatorP i := |ϕ i ϕ i | in position representation: 77 In Ehrenfest dynamics, the nuclear position Q(t) is simply the current position Q t of the trajectory, whereas the adiabatic population P i (t) = | ϕ i (Q t )|ψ t | 2 and electronic density ρ el (q, t) = |ψ t (q)| 2 depend on the electronic wavefunction ψ t .  eff,t is measured by the distance of the forward-backward propagated state x fb eff,t := SΦ appr,t (x eff,0 )] from the initial state x eff,0 . The norm x eff := x eff , x eff 1/2 of an effective phase space point x eff is defined using the scalar product x eff,1 , x eff,2 := Q T 1 · Q 2 + P T 1 · P 2 + ψ 1 |ψ 2 of x eff,1 and x eff,2 . The corresponding squared "distance" x eff,1 − x eff,2 2 between points x eff,1 and x eff,2 is simply the sum Q 1 − Q 2 2 + P 1 − P 2 2 + ψ 1 − ψ 2 2 of the squared distances between Q 1 and Q 2 , between P 1 and P 2 , and between ψ 1 and ψ 2 . Fig. 3 show that the energy is only conserved approximately, to the same order as the order of convergence of the integrator. The loss of exact energy conservation is standard for any splitting method 55,91 and is due to alternating kinetic and potential propagations: the effective Hamiltonian alternates between H eff = T el + T nu (P ) ψ and H eff = V (Q) ψ , and its time-dependent nature breaks the conservation of energy. Figure 4 confirms the predicted asymptotic order of convergence of the geometric integrators. However, panel (c) may wrongly suggest that the Suzuki-fractal scheme leads to the most efficient method, as it has the smallest error for each time step size. What Fig. 4 does not show is that the sixth-order Suzuki-fractal scheme has a factor of 25/9 more substeps per time step than either the triple-jump or optimal scheme. If we instead consider the dependence of the convergence error on the computational cost [measured by the central processing unit (CPU) time], the optimal composition scheme indeed yields the most efficient method for each order of accuracy (see Fig. S1 in Sec. S2 of the supplementary material).

Panels (g) and (h) of
To reach a modest convergence error of 10 −3 , the most efficient geometric integrator (obtained using the optimal fourth-order composition scheme) is 15 times faster than the two time step method and roughly twice slower than the three time step method (see Fig. 5).
Yet, a clear advantage of the geometric integrators over the other methods is the exact conservation of geometric properties. Both the two and three time step methods violate symplecticity; the two time step method, in addition, violates time reversibility (see Fig. 6  in Appendix A). Moreover, due to its higher order of convergence in ∆t, the fourth-order geometric integrator becomes more efficient than even the three time step method to reach convergence errors below 10 −4 (see Fig. 5). 10 -4

IV. CONCLUSION
We have demonstrated that the high-order geometric integrators for Ehrenfest dynamics can be obtained by simultaneously employing the splitting and composition methods.
Since Ehrenfest dynamics already involves a rather severe approximation, one is often not interested in numerically converged solutions. In such cases, geometric integrators become much more relevant because only they guarantee the exact conservation of the geometric invariants regardless of the accuracy of the solution. That is not to say that the high-order geometric integrators are inefficient for high accuracy simulations. On the contrary, to reach an error of, e.g., 10 −6 , using the eighth-order geometric integrator yields a four-fold speedup over the three time step method (see Fig. 5) and a 5000-fold speedup over the two time step method. High-accuracy results with negligible numerical errors may be desirable when the error introduced by the mean-field Ehrenfest approximation is either very small or unknown.

SUPPLEMENTARY MATERIAL
See the supplementary material for the computational details (Sec. S1), efficiency of the high-order geometric integrators obtained using the triple-jump, Suzuki-fractal, optimal composition schemes (Sec. S2), and detailed algorithms of the two and three time step methods (Sec. S3).
agreement No. 683069 -MOLEQULE) and thank Tomislav Begušić and Nikolay Golubev for useful discussions.

Conflict of interest
The authors have no conflicts to disclose.

DATA AVAILABILITY
The data that support the findings of this study are openly available in Zenodo at https://doi.org/10.5281/zenodo.5167211.
Appendix A: Two and three time step methods

Two time step method
Unlike the geometric integrators, which propagate Q t , P t , and ψ t simultaneously, the two time step method 11,12,53 consists in alternately propagating the classical nuclear phase space point and electronic wavefunction. The time step ∆t el = ∆t/n el for the electronic propagation is typically much smaller than the time step ∆t for the nuclear propagation (we used n el = 100).
In the two time step method, we employed the second-order Verlet algorithm 62 to propagate the nuclear phase space point and the second-order VTV split-operator algorithm 61 to propagate the electronic wavefunction. However, because the nuclear phase space point and electronic wavefunction are propagated separately and alternately, the overall two time step method is only first-order accurate in the time step. Moreover, the method is neither time-reversible nor symplectic (see Fig. 6). See Sec. S3 of the supplementary material for the detailed algorithm of the two time step method.

Three time step method
The three time step method, 13,54 owing to its symmetry, is both time-reversible and second-order accurate in the time step. 55 However, the method is still not symplectic (see Fig. 6).
In because dq ψ ∧ dp ψ (ψ 1 , ψ 2 ) = dq ψ (ψ 1 )|dp ψ (ψ 2 ) − dp ψ (ψ 1 )|dq ψ (ψ 2 ) where we have used that the tangent space of a vector space can be identified with the vector space itself, 112 i.e., d(M T t JM t )/dt = 0, which follows easily from the calculation where we have used the fact that the time derivative of the stability matrix satisfies 56 Although we did not need the explicit form of the Hessian of H eff (x eff ) to prove the symplecticity of M t , this expression will be useful for the numerical propagation of the stability matrix; we defined and used During the kinetic propagation, the equation of motioṅ for stability matrix M t is obtained by reducing the effective Hamiltonian to H eff (x eff ) = T nu+el (P ) ψ = T nu (P ) +T el ψ in Eq. (C6). The explicit form of the Hessian is and Eq. (D1) can be solved analytically to yield where and the propagation of quantum Darboux coordinates p ψ,t = −s T (P 0 )q ψ,0 + c T (P 0 )p ψ,0 for time t is equivalent to the standard propagation of electronic wavefunction During the potential propagation, the equation of motioṅ for stability matrix M t is obtained by reducing the effective Hamiltonian to H eff (x eff ) = V (Q) ψ in Eq. (C6). The explicit form of the Hessian is and Eq. (D8) can be solved analytically to yield where and the propagation of quantum Darboux coordinates for time t is equivalent to the propagation of the electronic wavefunction ψ t = exp[−itV (Q 0 )/ ]ψ 0 since ψ t = (q ψ,t + ip ψ,t )/ √ 2 .

Supplementary material for: High-order geometric integrators for representationfree Ehrenfest dynamics
This document provides information supporting the main text. It contains the computational details (Sec. S1), efficiency of the high-order geometric integrators obtained using the triple-jump, Suzuki-fractal, and optimal composition schemes (Sec. S2), and detailed algorithms of the two and three time step methods (Sec. S3).

S1. COMPUTATIONAL DETAILS
In the exact quantum simulation, the full wavefunction Ψ t was represented on a uniform four-dimensional grid that is a tensor product of two different two-dimensional grids for the nuclear and electronic degrees of freedom. A uniform grid of 101 × 101 points defined between Q n = −3 a.u. and Q n = 3 a.u. was used for the nuclear degrees of freedom.
As for the electronic degrees of freedom, a uniform grid of 64 × 64 points defined between q n = −10 a.u. and q n = 10 a.u. was used; this grid was also used to represent the electronic wavefunction in Ehrenfest simulations. We employed the second-order VTV split-operator algorithm with time step ∆t = 0.1 a.u. for the numerical propagation of Ψ t .
To obtain the initial states, we solved the electronic time-independent Schrödinger equa- here,T nu+el (P ) = T nu (P ) +T el .