Theory Department
Max Planck Institute of Microstructure Physics
Theory Department   >   Research   >   Correlated electron-nuclear motion

Exact factorization of the time-dependent electron-nuclear wavefunction

A. Abedi, F. Agostini, S. K. Min, C. Proetto, Y. Suzuki, F. Tandetzky
The Born-Oppenheimer (BO) approximation is among the most basic approximations in the quantum theory of molecules and solids. It is based on the fact that electrons usually move much faster than the nuclei. This allows us to visualize a molecule or solid as a set of nuclei moving on the potential energy surface generated by the electrons in a specific electronic state. The total wavefunction is then a product of this electronic state, , and a nuclear wavefunction satisfying the Schrödinger equation

Here, =(R1 ...RNn) , denotes the nuclear configuration and = (r1 ...rNe) represents the set of electronic positions. The concept of the potential energy surface, given in the BO approximation by is enormously important in the interpretation of all experiments involving nuclear motion. Likewise, the vector potential and the Berry phase associated with it, provide an intuitive understanding of the behavior of a system near conical intersections. Here and in the following, denotes the inner product over all electronic coordinates. Berry- Pancharatnam phases are usually interpreted as arising from an approximate decoupling of a system from "the rest of the world", thereby making the system Hamiltonian dependent on some "environmental" parameters. The best example is the BO approximation, where the electronic Hamiltonian depends parametrically on the nuclear positions; i.e., the stationary electronic Schrödinger equation is solved for each fixed nuclear configuration , yielding -dependent eigenvalues and eigenfunctions . Hence one has to acknowledge the fact that in the traditional treatment of molecules and solids the concepts of the potential energy surface and the Berry phase arise as consequence of the BO approximation. Yet, this is "just" an approximation, and some of the most fascinating phenomena of condensed-matter physics, like superconductivity, appear in the regime where the BO approximation is not valid. This raises the question: If one were to solve the Schrödinger equation of the full electron-nuclear Hamiltonian exactly (i.e. beyond the BO approximation) does the Berry phase and the potential energy surface survive, and if so, how and where does it show up? Moreover, many interesting phenomena occur when molecules or solids are exposed to time-dependent external fields, such as lasers. One may ask: Can one give a precise meaning to a time-dependent potential energy surface and a time-dependent Berry phase?

In a recent Letter [1] we are able to answer all of the above questions. We prove that the exact solution of the TDSE, , can be written as a single product

where satisfies the partial normalization condition, , for any fixed nuclear configuration, , at any time t. An immediate consequence of the identity (1) is that, is the probability density of finding the nuclear configuration at time t. The electronic wavefunction, , in the gauge where satisfies the equation




The nuclear wavefunction obeys the equation
. (3)

Via these exact equations the time-dependent potential energy surface (TDPES) :

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 H2+ molecule subject to the laser-field (see text), I1 = 1014W/cm2 (dashed line) and I2 = 2.5 ×1013W/cm2 (solid line). The circles indicate the position and energy of the classical particle in the exact-Ehrenfest calculation (I1: open, I2: 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 H2+ molecular ion in 1D, subject to a linearly polarized laser field. TDPES's along with the corresponding nuclear density, |χ(R,t)|2, 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.

Steps in the exact time-dependent potential energy surface

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, |Fl(R,t)|2 (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 |Fl(R,t) |2 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.

Non-adiabatic processes via mixed quantum-classical dynamics: A novel perspective from the exact time-dependent potential energy surface

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

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 H2+ 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

, (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 Rc(t) (the trajectory), |χ(R,t)|2→δ(R−Rc(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.


A. Abedi, N. T. Maitra and E. K. U. Gross, Phys. Rev. Lett., 105, 123002, (2010).
A. Abedi, N. T. Maitra and E. K. U. Gross, J. Chem. Phys., 137, 22A530, (2012).
A. Abedi, F. Agostini, Y. Suzuki and E. K. U. Gross, arXiv:1302.1453v1. Accepted for publication in Phys. Rev. Lett.

To top