Exact factorization of the timedependent electronnuclear wavefunction
A. Abedi, F. Agostini, S. K. Min, C. Proetto, Y. Suzuki, F. Tandetzky
The BornOppenheimer (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,
=(R_{1} ...R_{Nn})
, denotes the nuclear configuration and
= (r_{1} ...r_{Ne})
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 condensedmatter 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 electronnuclear 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 timedependent external fields, such as lasers.
One may ask: Can one give a precise meaning to a timedependent
potential energy surface and a timedependent 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
with
and
The nuclear wavefunction obeys the equation
Via these exact equations the timedependent potential energy surface (TDPES) :
and the timedependent Berry connection
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
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
realvalued (e.g. for a noncurrentcarrying 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 (manybody) current density, Eq.
6 implies that the gaugeinvariant 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 manybody wavefunction: Its absolutevalue squared
gives the exact nuclear (Nbody) density while its phase yields the correct
nuclear (Nbody) current density.


Fig. 1: Snapshots of the TDPES (blue lines) and nuclear density (black) at
times indicated,for the H_{2}^{+} molecule subject to the
laserfield (see text), I_{1} = 10^{14}W/cm^{2}
(dashed line) and I_{2} = 2.5 ×10^{13}W/cm^{2}
(solid line). The circles indicate the position and energy of the classical
particle in the exactEhrenfest calculation (I_{1}: open,
I_{2}: solid). For reference, the groundstate 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_{2}^{+}
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 timedependent potential energy surface
The concept of nonadiabatic 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 timescale 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 timescale 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 "nonadiabatic coupling" and "non
adiabatic transition" between the adiabatic states come to remedy the
adiabatic approximation and describe the dynamical processes. For example,
in the BornOppenheimer (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 intermolecular 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 nonadiabatic treatments are usually based on the essential pictures
that the adiabatic approximation provides. The standard way of studying and
interpreting "nonadiabatic" molecular processes is to expand the full
molecular wavefunction in terms of the BO electronic states. Within
this expansion, nonadiabatic 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 semiclassical treatment of the nuclei where a welldefined
single classical force would be highly desirable.
In our previous works, we have introduced an exact timedependent potential
energy surface (TDPES)
that, together with an exact timedependent vector
potential , governs the nuclear motion. These concepts emerge from a
novel way to approach the coupled electronnuclear dynamics via an exact
factorization,
, of the electronnuclear wave function^{[1]}.
The crucial point of this representation of the correlated electronnuclear
manybody problem is that the wavefunction
that satisfies the exact nuclear equation of motion leads to an N
body density and an Nbody current density that reproduce the
true nuclear Nbody density and current density obtained from the full wave
function^{[2]}.
In this sense, , can be viewed as the proper nuclear wavefunction. 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 Nbody 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 nonadiabatic 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
is forminvariant under the gaugetransformation, whereas
, defined as
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 ShinMetiu 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 timeevolution 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 wavepacket is also shown (black line). Right: populations of the BO states along the time evolution. The strong nonadiabatic 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 wavefunction (the vector potential
can be gauged away) and provides us with an alternative way of visualizing
and interpreting the nonadiabatic processes. We have shown that (see Fig. 4)
the gaugeinvariant 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 wavepacket 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 wavepacket branches and leaves the avoided crossing. The gaugedependent 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 BOPESpieces 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 timesteps, 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)^{2} (l=1 red line and l=2
green line), on different BO surfaces. The gray boxes define the regions in Rspace 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 wavepacket approaches a region of strong nonadiabatic coupling. The step feature (ii) is in
agreement with the semiclassical picture of nonadiabatic 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 wavepacket from a semiclassical point of view can be represented as an ensemble of classical
trajectories, along which pointparticles 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
electronnuclear dynamics semiclassically. 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 semiclassical approximation of the components, labeled as F_{l}(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 wavepacket 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.
Nonadiabatic processes via mixed quantumclassical dynamics: A novel perspective from the exact timedependent potential energy surface
We have started our journey towards a full ab initio treatment of the coupled electronnuclear 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 wavepacket 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 quantumclassical 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 timedependent 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 gaugeinvariant 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 gaugedependent 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_{2}^{+} in^{[1]}.
The Ehrenfest theorem relates the timederivative of the expectation value of a quantummechanical operator
to the expectation value of the commutator of that operator with the Hamiltonian, i.e.,
where the average is calculated over the full wavefunction,
.
When the electronnuclear
wavefunction is represented in the factorized form,
,
Ehrenfest theorem can be extended to the nuclear wavefunction
provided that the nuclear momentum operator is replaced with
,
and the full Hamiltonian
is replaced by the nuclear Hamiltonian
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)^{2}→δ(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 gaugeinvariant component, we have developed a
mixed quantumclassical 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 quantumclassical scheme
.
In order to take the classical limit, the nuclear wavefunction 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 zeroth order terms of the expansion
are considered, the nuclear timedependent Schrödinger equation leads to the HamiltonJacobi 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 HamiltonJacobi 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 quantumclassical 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 quantumclassical (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 quantumclassical approach, rigorously derived from the exact Eqs. (23), are in a very good agreement with results obtained from the exact solution of the timedependent Schrödinger equation.