(4)

and the time-dependent Berry connection

,

(5)

are defined as rigorous concepts. These two quantities mediate the coupling between the nuclear and the electronic degrees of freedom in a formally exact way. The vector potential can be expressed as

.

(6)

This equation is interesting in several respects. First, writing , the last term on the RHS of Eq. (6) can be represented as , so it can be gauged away. Consequently, any true Berry connection must come from the first term. If the exact Ψ(t) is real-valued (e.g. for a non-current-carrying ground state) then the first term on the RHS of Eq. 6 vanishes and hence the exact Berry connection vanishes. Second, since , is the true nuclear (many-body) current density, Eq. 6 implies that the gauge-invariant current density, , that follows from Eq. (3) does indeed reproduce the exact nuclear current density. Hence, the solution of Eq. (3) is, in every respect, the proper nuclear many-body wavefunction: Its absolute-value squared gives the exact nuclear (N-body) density while its phase yields the correct nuclear (N-body) current density.

Fig. 1: Snapshots of the TDPES (blue lines) and nuclear density (black) at times indicated,for the H₂⁺ molecule subject to the laser-field (see text), I₁ = 10¹⁴W/cm² (dashed line) and I₂ = 2.5 ×10¹³W/cm² (solid line). The circles indicate the position and energy of the classical particle in the exact-Ehrenfest calculation (I₁: open, I₂: solid). For reference, the ground-state BO surface is shown as the thin red line.

To demonstrate the usefulness of the approach we have calculated the TDPES's for a numerically exactly solvable model: the H₂⁺ molecular ion in 1D, subject to a linearly polarized laser field. TDPES's along with the corresponding nuclear density, |χ(R,t)|², are plotted in Fig.1 at six snapshots of time. For the stronger field the dissociation of the molecule is dramatically reflected in the exact TDPES. Also in the "exact Ehrenfest" approximation where the nuclei are treated as classical particles moving on the exact TDPES, the molecule dissociates (open circles in Fig.1). However, for the weaker field only exact quantum calculation leads to dissociation while in the "exact Ehrenfest" calculation the system gets stuck in a local minimum of the TDPES (solid circles), suggesting that tunneling is the leading dissociation mechanism. This reveals that the TDPES is a powerful interpretive tool to analyze and interpret different types of dissociation processes such as direct vs. tunneling.

The concept of non-adiabatic transition has been widely used to describe various dynamical processes. The root of this concept is the adiabatic treatment of dynamical processes in which the complete system is approximately decomposed to two parts, based on the assumption that part of the complete system usually changes in a much shorter time-scale than the rest and can be adjusted instantaneously to the adiabatic changes of the rest. The "fast" part of the system, within the adiabatic approximation, depends on the rest only via an "environmental" parameter that changes adiabatically compare to the time-scale in which the fast part changes. This implies that, the fast part of the system is treated independently for each environmental parameter that represents a specific configuration of the rest. However, this ideal picture may break down in different situations. This is when the concepts such as "non-adiabatic coupling" and "non- adiabatic transition" between the adiabatic states come to remedy the adiabatic approximation and describe the dynamical processes. For example, in the Born-Oppenheimer (BO) approximation, the "faster" electronic motion is separated from the "slower" nuclear motion and the Hamiltonian that describes the electronic motion has a parametric dependence on the nuclear configuration. However, this approximation may break down, for example when two or more BOPESs come close or even cross for some nuclear configurations. Some of the most fascinating and most challenging molecular processes occur in the regime where the BO approximation is not valid, e.g. ultrafast nuclear motion through conical intersections, radiationless relaxation of excited electronic states, intra- and inter-molecular electron and proton transfer, to name a few. On the other hand, an adiabatic description of dynamical processes, such as the BO approximation, provides a great deal of intuition that is fundamental to our understanding of the processes. Hence, the non-adiabatic treatments are usually based on the essential pictures that the adiabatic approximation provides. The standard way of studying and interpreting "non-adiabatic" molecular processes is to expand the full molecular wave-function in terms of the BO electronic states. Within this expansion, non-adiabatic processes can be viewed as a nuclear wave- packet with contributions on several BOPESs, coupled through the non -adiabatic coupling (NAC) terms which in turn induce transitions between the BOPESs. While this provides a formally exact description, one may nevertheless ask: Is it also possible to study the molecular process using a single PES? This question is particularly relevant if one thinks of a classical or semi-classical treatment of the nuclei where a well-defined single classical force would be highly desirable.

In our previous works, we have introduced an exact time-dependent potential energy surface (TDPES) that, together with an exact time-dependent vector potential , governs the nuclear motion. These concepts emerge from a novel way to approach the coupled electron-nuclear dynamics via an exact factorization, , of the electron-nuclear wave function^[1]. The crucial point of this representation of the correlated electron-nuclear many-body problem is that the wave-function that satisfies the exact nuclear equation of motion leads to an N- body density and an N-body current density that reproduce the true nuclear N-body density and current density obtained from the full wave -function^[2]. In this sense, , can be viewed as the proper nuclear wave-function. The time evolution of , on the other hand, is completely determined by the TDPES and the vector potential. Moreover, these potentials are unique up to within a gauge transformation. In other words, if one wants a TDSE whose solution yields the true nuclear N-body density and current density, then the potentials appearing in this TDSE are (up to within a gauge transformation) uniquely given by and ; there is no other choice. This also implies, that the gradient of this exact TDPES is the only correct force on the nuclei in the classical limit (plus terms arising from the vector potential, if those cannot be gauged away).

In our recent work, we have investigated the generic features of the exact TDPES without external laser but in the presence of strong non-adiabatic couplings. As a major result we observe that the exact TDPES exhibits nearly discontinuous steps connecting different static BOPES, reminiscent of Tully's surface hopping [R.K. Preston and J.C. Tully, J. Chem. Phys. 54, 4297 (1971)] in the classical limit.

To investigate the TDPES in detail, we write it as a sum of two parts: , defined as

,

(7)

is form-invariant under the gauge-transformation, whereas , defined as

,

(8)

is the part that depends on the choice of the gauge.

We have calculated the TDPES for a numerically exactly solvable model of Shin and Metiu [S. Shin and H. Metiu, J. Chem. Phys. 102, 23 (1995)], whose schematic representation in given in Fig.2.

Fig. 2: Schematic representation of the Shin-Metiu model system.R and r indicate the coordinates of the moving ion and electron, respectively, in one dimension. L is the distance between the fixed ions.

The parameters of the model were chosen such that there is a very strong coupling between the first two BOPESs. The first four BOPESs are shown in Fig.3 (left panel), along with the initial nuclear density. The same figure (right panel) presents the time-evolution of the populations of the BO states. We have shown that the exact TDPES exhibits nearly discontinuous steps connecting different static BOPESs, reminiscent of Tully's surface hopping in the classical limit.

Fig. 3: Left: lowest four BO surfaces, as functions of the nuclear coordinate. The first (red line) and second (green line) surfaces will be considered in the actual calculations that follow, the third and forth (dashed black lines) are shown as reference. The squared modulus (reduced by ten times and rigidly shifted in order to superimpose it on the energy curves) of the initial nuclear wave-packet is also shown (black line). Right: populations of the BO states along the time evolution. The strong non-adiabatic nature of the model is underlined by the population exchange at the crossing of the coupling region.

For the 1D model system studied in our work, the TDPES is the only potential that governs the dynamics of the nuclear wave-function (the vector potential can be gauged away) and provides us with an alternative way of visualizing and interpreting the non-adiabatic processes. We have shown that (see Fig. 4) the gauge-invariant part of the TDPES,ϵ_gi(R,t) , is characterized by two generic features: (i) in the vicinity of the avoided crossing,ϵ_gi(R,t) becomes identical with a diabatic PES in the direction of the wave-packet motion, (ii) far from the avoided crossing, ϵ_gi(R,t) , as a function of R, is piecewise identical with different BOPESs and exhibits nearly discontinuous steps in between. The latter feature holds after the wave-packet branches and leaves the avoided crossing. The gauge-dependent part, ϵ_gd(R,t) , on the other hand, is piecewise constant in the region where ϵ_gi(R,t) coincides with different BOPESs. Hence ϵ_gd(R,t) has little effect on the gradient of the total TDPES, but may shift the BOPES-pieces of ϵ_gi(R,t) by different constants causing the exact TDPES to be piecewise parallel to the BOPESs.

Fig. 4: TDPES and nuclear densities at different time-steps, namely t=0 fs, t=10.88 fs and t=26.61 fs. The different panels show: (top) GI part of the TDPES (black dots) and the two lowest BOPESs (first, dashed red line, and second, dashed green line) as reference; (center) the GD part of the TDPES (green dots); (bottom) nuclear density (dashed black line) and its components, |F_l(R,t)|² (l=1 red line and l=2 green line), on different BO surfaces. The gray boxes define the regions in R-space where the energies have been calculated, since the nuclear density is (numerically) not zero.

The diabatic feature (i) of the TDPES supports the use of diabatic surfaces as the driving potential when a wave-packet approaches a region of strong non-adiabatic coupling. The step feature (ii) is in agreement with the semi-classical picture of non-adiabatic nuclear dynamics provided by the Tully's surface hopping scheme, that suggests to calculate the classical forces acting on the nuclei according to the gradient of only one of the BOPESs.

The exact TDPES represented in Fig. 4 can beviewed from a different perspective. The nuclear wave-packet from a semi-classical point of view can be represented as an ensemble of classical trajectories, along which point-particles evolve under the action of a classical force which is the gradient of ϵ_gi. According to our observations, on different sides of a step such a force is calculated from different BOPESs. This is reminiscent of Tully's surface hopping approach, that deals with the problem of coupled electron-nuclear dynamics semi-classically. The method introduces stochastic jumps between BOPESs to select the adiabatic surface that, at each point in time, governs the classical nuclear dynamics. The nuclear density is reconstructed from bundles of classical trajectories. Such bundles evolve independently from one another on different adiabatic surfaces and are a semi-classical approximation of the components, labeled as |F_l(R,t) |² in Fig. 4, of the exact nuclear density. The step feature of the TDPES, following from the exact solution of the full TDSE, makes clear that, after the wave-packet splits at the avoided crossing, the motion of its components (the bundles in surface hopping language) is driven by single adiabatic surfaces and not (like, e.g., in Ehrenfest dynamics) by an average electronic potential.

We have started our journey towards a full ab initio treatment of the coupled electron-nuclear dynamics by introducing an exact separation of electronic and nuclear motions ^[1]. We have derived the equations of motion that govern the dynamics of each subsystem and studied the potentials that mediate the couplings between the two subsystems in a formally exact way ^[2]. In particular, we have studied features of the TDPES in two topically demanding situations: molecules in strong fields ^{[1, 2]} and splitting of a nuclear wave-packet at avoided crossings ^[3] of BO potential energy surfaces. These studies provide us with the essential elements (fundamental equations of motion and insights of the coupling potentials) for making approximations, especially for the systematic development of semiclassical approximations.

Here we report on the recent progress towards developing a mixed quantum-classical scheme based on the exact separation of the electronic and nuclear motions ^[3]. In order to examine the classical approximation of the nuclear dynamics, we first study the classical nuclear dynamics on the exact TDPES that are obtained from the exact solution of the time-dependent Schrödinger equation for an exactly solvable model, developed by Shin and Metiu (Fig. 2). We have studied the classical nuclear motion using the full TDPES and the gauge-invariant component of the exact TDPES by integrating Hamilton's equations

(9)

using the classical force calculated as the gradient of the potential (full or GI) at the position of the classical particle. We present some of the results in Fig.5, where we follow classical dynamics along the exact potential (EX) and along its GI part.

Fig. 5: The figure shows classical positions (dots) at different times, as indicated in the plots, with the corresponding potentials, ϵ_GI(R,t) (orange lines) and ϵ(R,t) (blue lines). The nuclear density (dashed black line) is plotted as reference.

As expected, the single trajectory approach implemented here is not able to capture all dynamical details of quantum evolution. However, some observables can be adequately reproduced by the evolution on the exact TDPES, as shown in Fig.6.

Fig. 6: Classical position (left panel) and velocity (right panel) and mean nuclear position and velocity as functions of time. The dashed black line represents the average nuclear values from quantum calculation, the blue and orange lines are the positions and velocities of the classical particle when it evolves on the exact potential and on the GI part of the potential, respectively.

Our results show the importance of the gauge-dependent component of the TDPES as the results obtained from the classical evolution on the GI part of the full potential only deviate largely from the exact results.

This procedure has been refereed to as "exact Ehrenfest" in the study of the dissociation of H₂⁺ in^[1]. The Ehrenfest theorem relates the time-derivative of the expectation value of a quantum-mechanical operator to the expectation value of the commutator of that operator with the Hamiltonian, i.e.,

,

(10)

where the average is calculated over the full wave-function, . When the electron-nuclear wave-function is represented in the factorized form, , Ehrenfest theorem can be extended to the nuclear wavefunction

	(11)
,	(12)

provided that the nuclear momentum operator is replaced with , and the full Hamiltonian is replaced by the nuclear Hamiltonian

,

(13)

where, the average operations in Eqs. (11) and (12) are performed over the nuclear density only. In 1D cases (A(R,t)=0 according to the gauge condition) and for a nuclear density that is infinitely localized at the classical position R_c(t) (the trajectory), |χ(R,t)|²→δ(R−R_c(t)), Eqs. (11) and (12) lead to Hamilton's equations (9).

After studying various aspects of the classical nuclear dynamics on the exact TDPES and its gauge-invariant component, we have developed a mixed quantum-classical method based on the exact factorization of the electronic and nuclear motions. Starting from the exact equations of motion for electrons and nuclei, we rigorously take the classical limit of the nuclear motion and derive a mixed quantum-classical scheme . In order to take the classical limit, the nuclear wave-function is expanded in an asymptotic series in powers of ħ following the van Vleck proposal [J. H. van Vleck, Proc. Natl. Acad. Sci. 14, 178 (1928)]. If only the zero-th order terms of the expansion are considered, the nuclear time-dependent Schrödinger equation leads to the Hamilton-Jacobi equation for the classical action associated to the classical nuclear dynamics. The nuclear Hamiltonian generating the evolution can be obtained from the exact nuclear Hamiltonian (13) by replacing with P_ν, the nuclear momentum evaluated along the classical trajectory. Newton's equation can be easily derived from the Hamilton-Jacobi equation. The classical evolution obtained according to the described procedure is coupled to the electronic equation obtained by expanding the exact electronic equation (2) on the adiabatic basis, , formed by the eigenstate, , of the BO Hamiltonian , and neglecting, as a first approximation, the spatial dependence of the coefficients, . The set of ordinary differential equations for these coefficients is coupled to the classical nuclear evolution equation. The performance of the mixed quantum-classical scheme in comparison with the exact solution of the TDSE is examined by calculating the populations of the BO states and the nuclear kinetic energy.

Fig. 7: Left panel: populations of the BO states as functions of time determined by quantum (blue) and mixed quantum-classical (MQC, dashed red) propagation schemes. Right panel: nuclear kinetic energy as function of time.

As it is seen in Fig. 7, the results of the mixed quantum-classical approach, rigorously derived from the exact Eqs. (2-3), are in a very good agreement with results obtained from the exact solution of the time-dependent Schrödinger equation.