Time-dependent variational dynamics for nonadiabatically coupled nuclear and electronic quantum wavepackets in molecules

We propose a methodology to unify electronic and nuclear quantum wavepacket dynamics in molecular processes including nonadiabatic chemical reactions. The canonical and traditional approach in the full quantum treatment both for electrons and nuclei rests on the Born–Oppenheimer fixed nuclei strategy, the total wavefunction of which is described in terms of the Born–Huang expansion. This approach is already realized numerically but only for small molecules with several number of coupled electronic states for extremely hard technical reasons. Besides, the stationary-state view of the relevant electronic states based on the Born–Oppenheimer approximation is not always realistic in tracking real-time electron dynamics in attosecond scale. We therefore incorporate nuclear wavepacket dynamics into the scheme of nonadiabatic electron wavepacket theory, which we have been studying for a long time. In this scheme thus far, electron wavepackets are quantum mechanically propagated in time along nuclear paths that can naturally bifurcate due to nonadiabatic interactions. The nuclear paths are in turn generated simultaneously by the so-called matrix force given by the electronic states involved, the off-diagonal elements of which represent the force arising from nonadiabatic interactions. Here we advance so that the nuclear wavepackets are directly taken into account in place of path (trajectory) approximation. The nuclear wavefunctions are represented in terms of the Cartesian Gaussians multiplied by plane waves, which allows for feasible calculations of atomic and molecular integrals together with the electronic counterparts in a unified manner. The Schrödinger dynamics of the simultaneous electronic and nuclear wavepackets are to be integrated by means of the dual least action principle of quantum mechanics [K. Takatsuka, J. Phys. Commun. 4, 035007 (2020)], which is a time-dependent variational principle. Great contributions of Vincent McKoy in the electron dynamics in the fixed nuclei approximation and development in time-resolved photoelectron spectroscopy are briefly outlined as a guide to the present work.


Introduction
To address one of the ultimate goals of theoretical molecular science, we here propose to incorporate nuclear quantum wavepackets represented in the Cartesian Gaussian functions multiplied plane waves into nonadiabatic electron wavepacket dynamics we have long been developing [1][2][3][4][5][6][7] with use of our proposed quantum least action principle (time-dependent variational principle) [8]. In order to describe why the formulation of simultaneous nuclear and electronic dynamics is needed in current scientific status and to clearly distinguish the present theory from canonical and the traditional methods based on the Born-Huang expansion [10], we track our pathways to the goal starting from two early works with the Born-Huang dynamics; one is electron scattering by molecules as an electron dynamics in the fixed nuclei approximation and the a e-mail: kaztak@fukui.kyoto-u.ac.jp (corresponding author) other being direct observation of nuclear wavepacket dynamics through time-resolved photoelectron spectroscopy. Both of them have been studied together with Vincent McKoy.
Theoretical foundation of molecular science was established far long ago in 1927 by Born and Oppenheimer [9], shortly after the birth of quantum mechanics, and is later followed by the so-called Born-Huang expansion [10] to describe a total electronicnuclear molecular wavefunction. This view naturally brought about the Born-Oppenheimer (BO) approximation, or fixed nuclei approximation, in which nuclei are supposed to undergo their dynamics on stationary electronic state energy hypersurfaces (potential energy hypersurfaces abbreviated as PES) [11][12][13][14][15][16]. Quantum chemistry responsible mainly to analyze electronicstructures has been making tremendous contributions to molecular science and chemistry, mostly in the aspect of energetics such as the estimate of energy barrier for chemical reactions and so on [17].
Such fixed nuclei approximation, or the notion of "instantaneous" electronic-state adjustment to moving nuclei, has been known to be quite useful in the calculations of electronic bound states. Interestingly though, the fixed nuclei approximation turns out to work well for "slow" electron collision with molecules and lowenergy photoionization [18]. Quantum chemistry has thus successfully extended its wing so as to cover electron dynamics. Vincent McKoy was a pioneer and had been the central figure in these fields.
Another interesting feature of the BO approximation is that it is valid even for fast nuclear wavepacket dynamics, in the time scale of as fast as femtoseconds. It had been predicted by Domcke [19], Engel [20], and their coworkers that such nuclear wavepacket dynamics can be observed with time-resolved photoelectron spectroscopy (TRPES). However, in order to carry out ab initio calculations of TRPES for realistic molecules, one needs accurate photoelectron scattering amplitude and nuclear wavepackets running on accurate (sometimes nonadiabatically coupled) PES. There is no surprise that McKoy, who was leading in the calculations of photoionization, led and guided the field of TRPES too. Indeed, a computational scheme to attain angleresolved femtosecond TRPES was proposed for the first time in 1999 by the joint group of McKoy and Takatsuka [21][22][23]. Now in the current age of attosecond laser chemistry, the approximation of stationary electronic state is not realistic enough and the BO approximation is not convenient either. The cutting-edge electronic state theory should be capable of describing dynamical electrons as quantum wavepackets that are kinematically (nonadiabatically) coupled with nuclear motions. Therefore, we have been developing a theory of such nonadiabatic electron dynamics, in which electron wavepackets propagate in time along simultaneously generated nuclear "paths" [1][2][3][4][5][6][7]. These paths can naturally and smoothly branch into pieces at each significant nonadiabatic transition region as many as the number of adiabatic potential energy surfaces that commit the nonadiabatic avoided crossings and conical intersections. Therefore, those paths are not classical (Newtonian) in general and driven by the so-called force matrix, in which nonadiabatic interactions are taken into account besides the adiabatic nuclear forces [2]. The advantages starting from the nonadiabatic electron wavepacket representation are (1) electron dynamics can be tracked directly and vividly as chemical reactions proceed, and (2) one can treat large molecules in size and in the number of involved electronic states, since there is no need to prepare global PES because of the on-the-fly nuclear dynamics [7]. More essential virtue of the nonadiabatic electron wavepacket dynamics is that it can cope with (1) highly quasi-degenerate electronic states, the fluctuation of which is huge and extremely frequent due to their mutual mixing by nonadiabatic coupling among them, and (2) the real-time dynamics of electron wavepackets that are nonadiabatically coupled with nuclei like protons as in charge separation in water splitting in Photo System II of plants, and so on. These materials and phenomena are practically inaccessible by the Born-Huang representation.
On the other hand, an obvious shortage of our nonadiabatic electron wavepacket approach is that the nuclear quantum wavefunctions are not yet fully incorporated. With these experience, progress, and difficulty in mind, we here propose a method to quantize nuclear dynamics in terms of nuclear Gaussian wavepackets (Cartesian Gaussians multiplied by plane waves as in the squeezed states), replacing the nuclear path dynamics. The nuclear wavepackets are determined simultaneously with electron wavepackets in full quantum mechanical manner by means of the dual least action principle, a time-dependent variational principle, developed by the present author [8]. Not only the shapes and spatial distributions but also the positions of the Gaussian centers are determined quantum mechanically.
This paper is organized as follows. In Sect. 2, we briefly review the theoretical framework of the Born-Huang expansion, highlighting its success in terms of a theory of electron scattering by molecules with fixed-nuclei approximation, and a photoelectronspectroscopic real-time observation of nonadiabatic nuclear wavepacket dynamics. We move on from the Born-Huang representation to what we call nonadiabatic electron wavepacket representation with nuclear path approximation in Sect. 3. Section 4 achieves the aim of this paper to incorporate nuclear wavepackets into the nonadiabatic wavepacket representation, by treating electronic and nuclear wavepackets uniformly on equal footing by means of the time-dependent variational principle. This paper concludes in Sect. 5 with some remarks.

Dynamics in the Born-Huang representation
Prior to presentation of nonadiabatic electron wavepacket approach, we briefly summarize the success and limitation of the Born-Huang representation based on the stationary-state electronic state theory.

Born-Huang representation
Our nonrelativistic molecular Hamiltonian has a form where many-body electronic Hamiltonian is defined as with the symbols r and R being the electronic and nuclear coordinates, respectively, whilep j andP k denoting the operators for their conjugate momenta of the jth and kth component of r (specifically written as r j ) and R (R k ), respectively. r and R are independent variables at this stage. m and M k are the relevant masses. V c (r, R) is the collective representation of the Coulombic interactions among electrons and nuclei where r a and R A are the ath electron and the Ath nuclei, respectively, and Z A indicate the nuclear charge on the Ath nuclei.
Since the nuclei move far more slowly than electrons to the order of the square root of the mass ratio or are typically 10 2 times slower, the notion of instantaneous adjustment of electrons to nuclear reconfiguration seems rather natural. This obviously breaks the relativity yet is acceptable mathematically. Therefore, the total wavefunction Ψ(r, R, t) subject to the time-dependent Schrödinger equation with H(r, R) is to be expanded in stationary electronic basis functions {Φ I (r; R)}, as [10] Ψ(r, R in Φ I (r; R) are regarded as a parameter, while those are to be treated as independent variables for χ I (R, t). The electronic basis functions {|Φ I } are naturally orthonormalized at each nuclear configuration as The bra-ket inner products in what follows represent integration over the electronic coordinates only. χ I (R, t) in Eq. (4) are supposed to describe the nuclear wavepacket subject to the coupled equations of motion as where and These matrix elements are the main players in the theory of nonadiabaticity. In addition to these, H el may contain spin-orbit couplings to take a partial account of the relativistic effects. The expansion of Eq. (4), called the Born-Huang expansion [10], is mathematically rigorous in itself.

Born-Oppenheimer approximation for energetics
Φ I (r; R) in Eq. (4) are quite often chosen to be eigenfunctions of H (el) (r; R) as at each R, and V K (R) are denoted as potential energy (hyper)surface (PES). The notion of the fixed nuclei based on Eq. (9) has given a theoretical ground of molecular science and has made huge contributions to various aspects in chemical science. In particular, a single term truncation of the expansion of Eq. (4), which is with Eq. (9) is widely called the Born-Oppenheimer (BO) approximation and has been very successfully applied to the study of the molecular ground states. Then, the electronic energy V K (R), the nuclear repulsion energy included, serves as a potential function (K = 1, 2, · · · ) for the nuclear wavefunction in such a manner that The BO approximation is expected to be very good as long as the ground state under study is well separated from the excited state manifold. The error in the total molecular energy up to the rotational and vibrational energies on a given potential energy surface for the stable ground state is expected to be very small down to the order [24] m M where m and M represent the electron and involved nuclear masses. When this ratio is 10 −4 , the error will be in the order of 10 −6 .

Electron dynamics in the fixed nuclei approximation with stationary-state theory for electron-molecule scattering
The notion of fixed nuclei approximation has been widely applied to stationary-state scattering theory for electron-molecule scattering and direct photo-ionization dynamics as well. It was virtually impossible before 1980s to calculate the differential cross sections of electron scattering by anisotropic polyatomic molecules even for elastic scattering, let alone the inelastic processes involving electronic-state excitation [25,26]. The stationary scattering theory for molecular photoionization within the fixed nuclei framework has also been extensively developed by Lucchese and McKoy [18].  [25,26]. As usual the total electronic Hamiltonian H (el) we consider is written as where H N is the target electronic Hamiltonian at a given molecular geometry and T N +1 is the kinetic energy operator for the incident electron. The Coulombic interaction potential V between an incident electron and the target is given as usual Under this situation we want to solve the stationary Schrödinger equation, instead of the eigenvalue problem like Eq. (9), for the continuum state, where E is a parameter that specifies the experimental collision energy and is not an eigenvalue as in Eq. (9). ψ (+) has a nonzero component in the asymptotic region (|r| → ∞) satisfying an appropriate boundary condition for a given collision status. Therefore, the scattering wavefunction ψ (+) is generally not square-integrable and thereby is not a member of the standard Hilbert space. However, what we really need is not the entire wavefunction ψ (+) but the scattering amplitudes, which give rise to the differential cross sections as physical observables. Therefore, ab initio calculations of the scattering amplitude is the goal in this study. Electron-impact vibrational excitation can be approximately evaluated with use of the resultant scattering amplitudes.
To this goal, we first define a projection operator P which specifies the open target eigenstates. By open channels we mean the states that are accessible within the collision energy [25,26]. P is an N -body operator unlike the (N + 1)-body projection operator of Feshbach [27]. Practically, P is simply represented in terms of target eigenstates Φ I such that P = open channels With this operator, we first project the formal Lippmann-Schwinger equation for the scattering state P ψ where ψ (+) I has the asymptotic form of an incoming plane wave and a scattered component, and S I is hence a product of the target wavefunction Φ (ad) I and an incident plane wave, i.e., The projected Green function in Eq. (17) is formally written with the one-body Green function g , respectively. We multiply the potential V on both sides of Eq. (17) such that In contrast to the original Lippmann-Schwinger equation, Eq. (20) is N + 1-body equation, and moreover, it is not symmetric. The asymmetry of Eq. (20) comes from the fact that we have simply projected only the open channel components of the entire dynamics. To remedy, we need to find projected counterparts from the closed channels. Such components are supposed to be well represented by the Schrödinger equation without need of an asymptotic form. We thus project the total Schrödinger Eq. (15) in such a way that Ĥ − a PĤ +ĤP ψ (el) and an arbitrary parameter a remains to be chosen later. Combining Eq. (20) and Eq. (21), we obtain an integro-differential equation in a symmetric form Closer examination of the operators demands that the parameter a must be chosen to be (N +1)/2 to make the operatorĤ − a PĤ +ĤP in the left-hand side Hermitian and the correct permutation symmetry among incident and molecular electrons be satisfied. (For an extension to positron-molecule scattering, see ref. [28].) Rewriting Eq. (22) with an operatorÂ aŝ a variational functional is given in a fractional form [29] such that Then the variational expressions for the inelastic scattering amplitudes naturally result as in Eq. (25) are always multiplied effectively by the potential function V (in the denominator too), which extends only in the finite molecular region [25,26]. This stationary-state scattering theory turned out to be quite powerful, and McKoy and his coworkers have indispensable contributions to science and technology of molecular electron scattering [30][31][32][33][34].

Atomic and molecular integrals with Cartesian Gaussians and plane waves
A technical aspect in carrying out the integration of the functional in Eq. (26) is highly relevant to the present study. In the evaluation of the matrix elements of Eq. (26), one needs the atomic integrals (kinetic, electronnucleus attraction, and electron-electron repulsion integrals) at each nuclear position, in which the Cartesian Gaussian functions multiplied by plane waves are involved as seen in Eq. (18). Fortunately, such techniques have been practically developed in electronscattering works. Besides, Pulay and his group have developed a very fast and accurate algorithms of those Gaussian and plane-wave atomic integrals [35]. These are particularly vital in the present work, since later in Sect. 4, we will represent nuclear wavefunctions in Cartesian Gaussians multiplied by plane waves and thereby treat both the electronic and nuclear wavepackets in a unified manner.

Nonadiabatic quantum dynamics
In rising from the polyatomic dynamics based on the Born-Oppenheimer approximation, Eq. (11), to the coupled equations, Eq. (6), one faces big technical difficulties. Review articles about chemical dynamics for polyatomic molecules in this line are given in Refs. [11][12][13]. To integrate Eq. (6), one needs to prepare the potential energy surfaces V I (R) in a wide spatial range prior to the computation of the {χ I (R, t)}, since the spatial extent of χ I (R, t) s are not known beforehand. Moreover, the nuclear wavefunctions χ I (R, t) are very oscillatory in space due to the heavy masses of nuclei. Besides, the physical dimension of R increases by 3 as one atom is added. Hence, 5 to 6 atomic systems already hit the practical limitation of the direct applications.
Historically the theory of nonadiabatic transition was introduced by Landau, Zener, and Stuekelberg as early as in 1932 as a one-dimensional two-state coupled Schrödinger equations with simple model potential curves. Such one-dimensional semiclassical theory has been mathematically finalized by Nakamura and Zhu [36]. For further progress, we refer to reviews about the rich physical implications emerging from the realistic molecular nonadiabatic chemistry, including effects arising from the genuine multidimensional nonadiabatic interactions [5][6][7][36][37][38][39][40].
Below are among the modern theories of nonadiabatic dynamics relevant to the present work in that quantum wavepackets dynamics are directly considered. Martínez and his colleague have formulated the so-called ab initio multiple spawning method and applied very successfully to nonadiabatic and excitedstate dynamics of chemically and biologically interesting molecules [41,42]. The very basic gradient is to propagate Gaussian nuclear wavepackets along on-thefly paths [43] running on nonadiabatically coupled ab initio potential surfaces, incorporating the surface hopping algorithm [44]. (Yet, the theory itself has been deepened so as to allow the mixing of the nuclear wavepackets on the way of reactions [45]). Burghardt and her coworkers proceeded to incorporate a diabatic two-state linear vibronic coupling model into the onthe-fly multi-configuration (nuclear) Gaussian wavepacket dynamics [46]. The theory is sophisticated in that it address the two-(or more-) layer Gaussian multiconfiguration time-dependent Hartree method [13] in the context of on-the-fly scheme, which are expected to have a vast application realm of large molecular systems of many vibronic modes. Worth and his group have developed the improved direct dynamics variational multi-configurational Gaussian method for nonadiabatic dynamics [47].
A quite unique approach, which determinedly abandons the Born-Huang representation, is the so-called exact-factorization approach by Gross and coworkers [48][49][50]. Both electronic and nuclear states are to be determined starting from their newly established coupled equations of motion. This theory is very unique and brings about a new way of studying nonadiabatic dynamics. It is so revolutionary that even the very basic concepts historically well established based on the Born-Oppenheimer approximation are to be revisited from scratch. Finally, we slightly touch upon the works of nuclear-electronic orbital methods, mainly developed in quantum chemistry as a direct extension of the molecular orbital model. We simply refer to the review article by Pavošević, Culpitt, and Hammes-Schiffer [51], and extensive citation in it, including Thomas. [52] and Nakai et al. [53].
The present paper also aims at a full quantum mechanical construction of electronic and nuclear wavepackets beyond the Born-Huang representation, yet with a deep appreciation of the great heritage and methodologies that have been developed in the long history of quantum chemistry.

Theory of time-resolved photoelectron spectroscopy in the Born-Huang representation:
Monitoring the nonadiabatic nuclear wavepacket bifurcation

Theory of time-resolved photoelectron spectroscopy (TRPES)
Before proceeding to the study of electron wavepacket dynamics, it would be fair to show an example of the remarkable success of the Born-Huang representation in full quantum nonadiabatic dynamics. In this particular study, we have made a combination of the nuclear wavepacket dynamics on nonadiabatically coupled potential surfaces and photoelectron scattering amplitudes to realize energy-, angle-, and time-resolved photoelectron spectroscopy. This turns out to provide an extremely powerful means both theoretically and experimentally to monitor femtosecond real-time chemical dynamics [21,22].
a. Theoretical outline of pump-probe photoelectron spectroscopy Let us generally consider a nuclear wavepacket running on the first electronic excited state, say Φ 2 (r; R), that is nonadiabatically coupled with the ground state, Φ 1 (r; R). These electronic states may be represented adiabatically or diabatically, subject to the two state equivalence. The initial state prepared on Φ 1 (r; R) with a nuclear wavefunction χ 1 (R, t) is to be pumped up by a short pulse laser to partially create χ 2 (R, t)Φ 2 (r; R), and up there it begins to move spontaneously after the pump laser is switched off. Eventually χ 2 (R, t)Φ 2 (r; R) are supposed to be bifurcated to χ 2 (R, t)Φ 2 (r; R) and χ 1 (R, t)Φ 1 (r; R) and/or mixed up together in nonadiabatic regions. To track and probe such wavepacket dynamics, we shine another pulse laser with delay times, the energy of which is high enough to ionize the dynamical molecular state. Since photoionization does not have optically dark states, and since photoelectrons are very well captured with good resolutions of energy and angle by means of the imaging technique, time-resolved photoelectron spectroscopy is now one of indispensable techniques for ultra-fast chemical dynamics [54][55][56].
Theoretical task to estimate photoelectron distributions in the space of the wave vector k, measured with respect to the frames of molecule and polarization vector of probe laser, is as follows [21,22]. The molecular system is expanded in the Born-Huang expansion, as where Φ (−) k (ion state) is a photoionizing stationarystate electronic wavefunction with an appropriate boundary condition. The calculations of Φ (−) k to represent ionization from both Φ 1 (r; R) and Φ 2 (r; R) at each nuclear configuration R and each wave vector k are among difficult parts in this work and the solving the nuclear dynamics for the continuum nuclear wavefunction χ k (R, t) is another one. The coupled timedependent Schrödinger equations for χ 1 , χ 2 , and χ k are of the form (in atomic units) and (29) where T N is the kinetic energy operator; V k,d is the ionization interaction matrix element, including dependence on molecular geometry, photoelectron energy and angle, probe pulse shape, and polarization (θ P ). Molecular rotation has been neglected. The equation for χ 2 is analogous to that for χ 1 and omitted here. The photoelectron spectra are calculated with χ k after the probe interaction is well over as

Wavepacket bifurcation as observed with TRPES
In a series of papers on time-resolved photoelectron spectroscopy (TRPES), Arasaki, Takatsuka, Wang, and McKoy showed that the precise instance of the wavepacket bifurcation of NaI at an incidence of intramolecular electron transfer NaI ↔ Na + I − as schematically drawn in Fig. 1(a) can be tracked in terms of the femtosecond time-energy-and angleresolved photoelectron spectra [57,58]. In these potential curves in the diabatic representation, the character of Φ , remains as is in the region of R > 7Ȧ. And they cross each other at a point near R ∼ 7Ȧ. As illustrated in the figure, nuclear wavepackets are supposed to bifurcate. To make the bifurcation directly observable, a linearly polarized laser pulse of frequency ω 1 pumps a groundstate wavepacket onto the lowest covalent state (A 1 Σ + , Ω = 0 + ), and then another linearly polarized pulse of frequency ω 2 is shot to ionize the pumped states with a variable delay time. The mutual angle between the pump and probe laser vectors is varied parametrically. Figure 1(b) and (c) shows the photoelectron kinetic energy spectra and the angular distributions, respectively, for pump-probe at delay times 190, 220, and 245 fs with the probe laser of parallel polarization to that of the pump laser. The absolute square of the wavepackets on the diabatic potentials at selected delay time is shown in the panel (a). The red and blue ones represent χ (R), respectively. These wavepackets are seen to pass across the crossing point from the left to the right of the crossing point R X after pumping. It is clearly seen that reshaping of the peaks of photoelectron signals, from one to three peaks, well reflect the wavepacket bifurcation. These peaks can be comprehended primarily with the Condon principle: The wavepackets running with the lower potential energy make the higher photoelectron signals. As for the interpretation of the photoelectron angular distribution, we refer to the original papers [57,58].
The above NaI study is just one of illuminating examples of the ab initio calculations of TRPES. Yet, experimental observation is of course not as easy. Nevertheless, direct experimental methodologies to observe the passage of wavefunctions through nonadiabatic regions have been developed since then. For instance, Wörner and his group for the first time observed the passage of NO 2 molecule across the conical intersection between the lowest A and B states by means of high harmonic spectroscopy [59] and more recently time-resolve photoelectron spectroscopy [60]. More of cutting-edge experimental methods have been introduced to identify the passage of nonadiabatic transitions, such as attosecond stimulated X-ray Raman spectroscopy [61], timeresolved fluorescent spectroscopy [62], ultrafast electron diffraction technique [63], attosecond XUV transient absorption spectroscopy [64], and so on.

Nonadiabatic electron wavepacket representation with nuclear path approximation
In the preceding section, we have seen a typical example of wavepacket bifurcation in NaI dynamics and shown how it can be directly observed with TRPES. In this example, only two adiabatic states are involved in nonadiabatic transitions, which is very well approximated as χ 1 (R, t)Φ 1 (r; R) + χ 2 (R, t)Φ 2 (r; R) in the Born-Huang representation, Eq. (4), where the view of nuclear wavepacket bifurcation seems natural.
Meanwhile, the Born-Huang (B-H) expansion can be also regarded as a representation of electronicstate fluctuation in terms of time-dependent coefficients χ I (R, t) over time-independent electronic-state basis functions Φ I (r; R). This can be comprehended clearly if we consider a fixed nuclei situation in Eq. (4) by setting R = R 0 (a constant), since the total wavefunction looks Ψ(r, t) = which is a linear combination of the electronic basis functions with time dependent coefficients. However, as stressed several times, obtaining χ I (R, t) numerically is extremely hard in general. Therefore, in order to survey the effects of large electronic fluctuation and/or ultrafast events for electrons in attosecond-scale, the formalism of nonadiabatic electron wavepacket representation is far more realistic and appealing to our intuition than the B-H one. To convince that this is certainly the case, we begin with a couple of such examples below.

Some characteristics of electron wavepacket dynamics
To represent a single bunch (before bifurcation) of electronic-nuclear wavepacket state, the following representation of such total wavefunction should be convenient in that the electronic wavefunction is also an explicit function of time as the nuclear counterpart is. Thus, both Φ elec and G nuc are wavepacket states. R t is the quantum mechanical average of the nuclear coordinate R over G nuc (R, t) and serves as a molecular frame for the description of the electronic states. The primary consequence of this representation is that nonadiabatic electron wavepackets Φ elec (r, t; R t ) are in general complex-valued having the dynamical phases, in a marked contrast to the real-valued feature of the stationary state representation of Φ I (r; R) of Eq. (4). Therefore, we can calculate electron flux (current of probability density), which is defined as [65] with ρ(r , r, t) being the first-order spinless reduced density matrix taken from the electron wavepacket Φ elec (r, t; R t ). (For molecular flux, see ref. [66][67][68].) With the flux one can directly track the electron flow in molecules as chemical dynamics proceeds. [7] Moreover, even real-time ionization dynamics such as photoionization and autoionization can be tracked by making appropriate use of the electron flux. For instance, Matsuoka and Takatsuka have shown how electrons are pumped up to an ionization manifold by the assistance of nonadiabatic couplings in autoionization. [69,70] The second illuminating example is a huge electronicstate fluctuation in the manifold of densely quasidegenerate electronic states, such as those in the excited states of boron clusters. [71][72][73] Another example is often seen in electron capture by the so-called Rydberg state on atoms in molecules. Since the electronic states are so close to their neighbors in energy, and since the strong nonadiabatic couplings mix them up frequently everywhere in the clusters, the notion of "isolated electronic state" gets lost within several femtoseconds. Even in such extremely tough situations, the clusters keep to exist without breaking apart and, moreover, specific chemical reactions are expected to occur in them. It is obvious in these cases that the notion of isolated adiabatic electronic states having an isolated potential energy surface (PES) based on the Born-Huang representation loses the sense. We refer to this situation as chemistry without potential energy surfaces. [71][72][73] The third one is taken from our study about the mechanism of charge separation dynamics in water splitting catalyzed by Mn 4 CaO 5 cluster in Photo System II. The most critical aspect in this process is how electrons and protons are separated and transmitted from water molecules via Mn 4 CaO 5 to the relevant protein residues. As a matter of fact, the water oxidation catalyst is not merely a bare Mn 4 CaO 5 and charge separation reactions proceed (surprisingly) in the electronic ground states. [74] The actual cluster we treated is as large as Mn 3 Ca(H 2 O) 2 (OH) 4 (HCOO) 5 -OH-Mn(H 2 O) 2 , which is further tied with protein residues through water clusters. These examples and much more suggest that the nonadiabatic electron wavepacket dynamics play inevitable roles in the studies of molecular systems of large electronic fluctuation.

Entanglement representation of electronic and nuclear Hamiltonian
To overcome the difficulties inherent to the Born-Huang expansion, and moreover, to cope with experimental studies on electron dynamics in modern chemical physics, we resume our basic discussion with writing down the following total Hamiltonian operator [1,2] H(R, elec) which is represented in a basis set {|Φ I (R) |R } andP k is the nuclear canonical momentum operator. Note that P k are to be operated only on the nuclear wavefunctions but not on the parameter R in the electronic wavefunctions Φ I (r; R), since such kinematic interactions have been already taken into account in Eq. (34) through X k IJ , which was defined in Eq. (7). (See Appendix, too.) The second-order derivatives Y k appear after the square operation of the kinetic energy operator. In this representation, the nuclear motion is described in coordinate space (R-space), while the electron dynamics is represented in the electronic Hilbert space. The electronic state summation over H

Nuclear path approximation
It is still not easy to directly solve the Schrödinger dynamics with Eq. (34) numerically in the current stage. We therefore first take an approximate approach to "electron wavepacket dynamics along nuclear (nonclassical) trajectories", by approximating the nuclear momentum operatorP k with its classical counterpart P k in such a way as This is an ultimate form of the so-called mixed quantumclassical representation. We often adopt the massweighted coordinates in which the nuclear masses are all set to unity. The advantage of this quantum-classical representation manifests when the quantum Hamiltonian of Eq. (1) is classicalized as above by replacingP k with P k , where we only have in which no kinematic couplings X k IJ , the key quantities in the study of the beyond-Born-Oppenheimer chemistry, are left behind.
One of the systematic methods to treat electron dynamics on the basis ofH(R, P, elec) of Eq. (35) is to resort to the path-branching representation theory: Electron wavepackets are full quantum mechanically propagated in time along branching nuclear paths. [2,[5][6][7]76] As usual, the electron wavepacket Φ elec (r, t; R(t)) (with R t = R(t) and G nuc (R, t) = δ (R − R) 1/2 in Eq. (32), R(t) being a nuclear path to be determined later) is expanded in a set of time-independent wavefunctions {Φ I (r; R)} at each nuclear configuration R such that Φ elec (r, t; R(t)) = I C I (t)Φ I (r; R)

R=R(t)
, (37) with the time-dependent coefficients C I (t) to be evaluated along R(t). Then the coupled equations of motion for electron wavepackets are reduced to Here again the bra-ket notation used demands integration over the electronic coordinates r. The secondorder derivative terms Y k IJ in Eq. (7) are quite often neglected because they are always accompanied by the small quantity 2 , although it is not always negligible in general [5,6].
The nuclear paths R(t) are driven in turn by the force matrix F k IJ (see ref. [2]) expressed as which arises from the formal Hamilton canonical equations of motion [2]. If the force matrix is diagonal, that is, F k IJ = 0 for I = J, each of F k II gives a Newtonian force on each potential function I, as in the socalled ab initio dynamics, molecular dynamics, and so on. This has given the theoretical foundation of ab initio molecular dynamics. On the other hand, nonzero offdiagonal elements of F k IJ can induce branching of the nuclear paths, realizing a nonadiabatic dynamics. From such path branching at each nonadiabatic region, infinite number of branching-paths are proliferated in the exact solutions of Eq. (39). To avoid such infinite generation of paths, we have proposed a systematic approximation of phase-space average of naturally branching paths (PSANB) [76]. We do not track this route in this paper, since the goal of this paper is to replace the branching paths with quantum nuclear wavepackets.
Incidentally, the mean-field path approximation, or the so-called semiclassical Ehrenfest theory (SET) is the simplest approximation to avoid the difficulty of path branching, in which the force matrix of Eq. (39) is to be averaged over the electron wavepacket in such a way that This gives a single scalar force to drive a single path. If the basis set {Φ I (r; R)} happens to be complete, Eq. (40) is reduced to the form of Hellmann-Feynman force Precisely speaking, the conventional SET does not include the second-order derivative terms Y k IJ in Eq. (38), which is widely applied in the literature [77][78][79][80][81].

Putting semiclassical nuclear wavepackets along the branching paths
We here show that those branching paths driven by the matrix force can be quantized my means of our developed semiclassical nuclear quantum wavepacket method. [82,83] The packets are referred to as the normalized variable Gaussian (NVG) wavepacket. [84] The NVGs vary in time their amplitudes, complex exponents, and phases subject to the dynamical equations for action decomposed function (ADF), to be defined below, along "predetermined paths". This procedure is just one-step before the simultaneous electronic and nuclear full quantum wavepacket dynamics.

Dynamics of Action Decomposed Function
We first present a semiclassical theory to attach quantum wavepackets on given classical paths [82,83]. Let us begin with a time-dependent wavefunction of the Maslov form [85,86] on a coordinate R in configuration space where S cl is assumed to satisfy the classical Hamilton-Jacobi (HJ) equation We refer to F (R, t) of Eq. (42) as Action Decomposed Function (ADF). By this factorization, the purely quantum factors are all involved in the function F (R, t), [86] and the Schrödinger equation for Ψ (R, t) is transformed to a linear equation of motion for the complex valued amplitude function F (R, t) as where p is a (nuclear) momentum at (R, t), which is The Trotter decomposition for a very short time allows for F (R − R (t) , t) here represents a solution, which is carried by the classical flow in the Lagrange view. [82,83] The term represented by exp − 1 2 (∇ · P )Δt is referred to as momentum gradient, while that by exp i 2 Δt∇ 2 to quantum diffusion. We can track F (R − R (t) , t) according to ref. [82][83][84].

Gaussian approximation
So far the formal theory is rigorous in the limit of Δt → 0. However, to step further for numerical realization, we approximate F (R − R (t) , t) with a Gaussian function of the form, N (t) being a normalization factor. Within the single Gaussian approximation both the dynamics of momentum gradient and quantum diffusion can be performed almost rigorously. It has been found [82,83] that the inverse exponents c (t) and d (t), both being real number (matrix), satisfẏ where the dot above the symbols indicates the firstorder time derivative. The so-called deviation determinant σ (t) is expressed as which is an N -dimensional orientable tiny volume surrounding the point R (t) in configuration space, representing how the classical flow nearby R (t) behaves, where R i (t) is the ith classical trajectory running nearby R (t). Equations (48) and (49) take account both of momentum gradient and quantum diffusion collectively. We refer to refs. [82,83] for multi-dimensional calculation of σ (t), its meaning of the absolute invariance, the associated Maslov index, and classical chaos. [87] We call the function of Eq. (47) as ADF-NVG, with NVG being for normalized variable Gaussian.

Example of branching nuclear packets along branching paths
We next show an example in which the NVG are propagated in time along the non-Born-Oppenheimer branching paths (PSANB) created by electron wavepackets. The total wavefunction is of the general form where k and K, respectively, indicate that the kth path is running on the Kth potential energy surface.
In what follows, we use the ADF-NVG described above as G Kk (R − R Kk (t) , t) in Eq. (1) with integrating the classical action S Kk along the path. In the following calculations, the momentum gradient in Eq. (46) has been treated by another form of approximation other than one using the deviation determinant of Eq. (50), which is (using the mass-weighted coordinates, in which since this expression can be applied even to non-Newtonian paths, provided that the information (q(t), v(t)) is available. [84] The adiabatic potential functions along with the nonadiabatic coupling element are drawn in each column of Fig. 2 in red solid curves (see ref. [84] for details). These curves model the essential character of the ab initio potential curves of LiH molecule and alkali halides. The bell-shape function also in a red solid curve superposed on these potential curves represents the nonadiabatic coupling element. A path-branching calculation based on PSANB is carried out first. In this example system the semiclassical Ehrenfest theory fails, since any Ehrenfest path starting from the basin area can never proceed to the asymptotic region of the dissociation channel. [76] Next, G(R − R (t) , t) of Eq. (47) and the relevant (formal) action integral are associated so as to run along the branching paths. As for the initial ADF-NVG in the PSANB approximation, we prepare the following function , all in atomic (or absolute) units. S 0 is chosen to be the same as that of Eq. (55). The corresponding full quantum nuclear wavepackets have been also generated for comparison. The initial Gaussian packet is prepared at a position R 0 , with the initial wave number k 0 . The width is set to δ 0 = 10/k 0 so as to be in the same order of its corresponding de Broglie wave length. The explicit form is where S 0 is the function mimicking the classical action with V 0 being a potential gradient at R = R 0 on the initial electronic state energy surface. S 0 is chosen so as to make a consistent comparison. We here exemplify two cases: One is a vibrational decay through the nonadiabatic coupling, the initial wavepacket of which is prepared at the left cliff of the lower curve (R 0 = 2, with the initial momentum k = 25) in panel (A) of Fig. 2. The other one, panel (B), is a collision event, with an initial packet coming in from the dissociation channel (from R 0 = 7 , with k = −45). These total wavefunctions of the full quantum and ADF-NVG are projected onto the nuclear configuration space and compared at selected times. In both cases the real parts of Full Quantum (in red solid curve) and PSANB-ADF-NVG (blue broken curves) are compared. The agreement up to the fine oscillatory structures is seen to be quite good even for this lowest level approximation in the construction of the PSANB-ADF-NVG wavefunctions. [84] Meanwhile, we also see the limitation or the validity range of the present theory. Among others, there is no mathematical mechanism involved to feed the effect of such quantization of nuclear dynamics back to the electron wavepacket dynamics.

Time-dependent variational determination of simultaneous electronic and nuclear quantum wavepackets
As shown above, the semiclassical incorporation of the Gaussian wavepackets (ADF-NVG) into nonadiabatic electron wavepacket dynamics has marked an excellent achievement as a nuclear and electronic quantum wavepacket method. Thus far, however, the molecular nuclei have been treated within the scheme of path dynamics, that is, the nuclear Gaussian functions thus running have no freedom for themselves to determine their own (quantum) pathways. We therefore extend the scheme so that not only the electronic wavepackets but the nuclear wavepackets can determine their dynamics subject to full quantum dynamics. In doing so, we apply our developed time-dependent variational principle (TDVP), with which we track in time the set of parameters that characterize trials electronic and nuclear wavefunctions. The bottom line of the accuracy of the designed method is warranted by the above pathbranching method using ADF-NVG.

Dual least action principle to determine the total wavefunctions
We begin with the equations of motion in a general parameter space first with which to track the Schrödinger dynamics. TDVP is a practical methodology in that it transforms the Schrödinger partial differential equation to a set of coupled ordinary differential equations over the space of variational parameters. There have been proposed various elaborated formalisms in the literature, such as those of the so-called Dirac-Frenkel, McLachlan [88], Kan [89], Kramer and Saraceno [90]. A unified account of these seemingly dif-ferent theories was given by Broeckhove, Lathouwers, Kesteloot and van Leuven [91]. These methods and their variants are widely used in many fields of physics and chemistry. (For more recent progress in other general theories of TDVP, we refer to refs. [92][93][94][95].) Yet, it is widely known that the existing TDVPs generally bear a divergence problem, commonly arising from an inversion procedure of singular matrices. To overcome the technical matters fundamentally, and moreover, in order to explore the theoretical structure of quantum mechanics from the scope of axiomatic variational principles, we have proposed a TDVP that is based on the dual quantum mechanical Maupertuis least actions. [8] Below we apply this quantum mechanical Maupertuis-Hamilton least action principle having a set of dual variational functionals.

Dual variational functionals
To represent the variational principle generally, we adopt general trial wavefunctions in a form φ (u, v) in this particular subsection, where real-valued variational parameter vectors u and v are supposed to have pairwise components (u i , v i ). Any TDVP generally starts from the stationarity of the following single variational functional On the other hand, the Hamilton action principle plays an axiomatic role in classical mechanics, which is with q and p being a coordinate point and its associated momentum, respectively. (The mass-weighted coordinates are used here also.) The canonical equations of motion resulting from Eq. (57) are which is essentially all about classical mechanics. Therefore, it is attractive to take an analogy to Eq. (57) by defining the following variational functional where C indicates arbitrary time interval under study.
In what follows, φ is always assumed to be normalized, and therefore the associated normalization factor may be subject to the variational operation. An application of the variational principle to S H alone gives under the fixed boundary condition Since both δu and δv are individually arbitrary, Eq. (60) results in the canonical equations of motion in the parameter space and It is obvious that along the flow lines (u(t), v(t)) in parameter space thus determined, the energy conservation is ensured. This seems fine, yet is only partially satisfactory from the view point of Eq. (56), in which the Maupertuis (reduced) action v ·u of Eq. (59) does not appear there. The Maupertuis action can be actually retrieved in such a way that This subtraction is mathematically valid since the variation is a linear operation. Therefore, the Maupertuis action has been cancelled in Eq. (56) and is formally hidden. Thus, the variational functional appearing as a counterpart in Eq. (64) is defined as Obviously, S W and S H are not independent from each other, and indeed they are required to satisfy the simultaneous set of variational principles , v (t)) = 0 (66) or more simply where (u (t) , v (t)) in S H and S W should be common. Equation (66) or (67) is referred to as a dual least action principle.
The variational counterpart δS W = 0 has a different physical meaning from that of δS H = 0. The formal application of δS W = 0 gives rise to and as Eqs. (62) and (63). Equations of motion (68) and (69) represent a "flow conservation" (see below) along the path The duality in the present least action principle is a reflection of the particle-wave duality of quantum mechanics, that is, δS H (u, v) = 0 gives a dynamical (or geometrical) restriction over the particle nature, while δS W (u, v) = 0 provides a constraint over the flow of matter wave. In particular, we may emphasize δS W (u, v) = 0 is responsible for the correct description of the quantum phase associated with φ. This is understandable if we recall that is widely denoted as the Berry phase. Incidentally, we here confirm that v ·u is regarded not merely as a factor to induce the Legendre transformation but as a machinery to make resultant paths stationary and invariant in the parameter spaces, realizing v ·udt = dv ∧ du = − du ∧ dv = − u ·vdt. (72) v ·u in this context is identified as a quantity related to an absolute invariance in the parameter space [96].

Quantum flux in the parameter space for the wave dynamics
We further survey the physical meaning of δS W = 0. Rewrite the main part in S W as which is real-valued. Note in this expression for instance, corresponds to the gradient in u-direction of the field of fluid. This is regarded as a quantum mechanical flux [65] in parameter space. (Recall Eq. (33) for the quantum flux in configuration space.) Thus, by analogy, we may define fluxes as which is the probability current in the u direction, and likewise the flux in the v-direction is We have after all which gives a natural interpretation of φ|φ in terms of the flow in the parameter space. Equations (68)-(69) are hence rewritten, respectively, as and Note, however, that since du/dt and dv/dt appear in the right-hand sides, a nonlinear nature of the dynamics seem to make the situation complicated.
Note that we see the following factor in Eqs. (79) and (80), which is a symplectic inner product between the Hamiltonian derivative vector and the flux vector. This quantity demonstrates a unique coupling between the Hamilton dynamics (particle dynamics) and the motion of the spatial redistribution of probability density (wave dynamics).

Simultaneous dynamics of particles and waves
The two sets of dynamical equations, Eqs. , v(t)) along the correct trajectories. More precisely, we regard them as a short time approximate solution in such a way that and likewise from Eqs. (77)-(78) The suffices H and W indicate, respectively, the Hamilton dynamics and the wave-flow dynamics. In practical approximations, these two solutions at a finite time t will deviate from each other after some finite time.
Likewise, we need to be careful about the choice of initial conditions (u(0), v(0)). Practically, we need to impose the conditions on each pair of the initial components (u i (0), v i (0)), in a manner consistent with the conditions for other parameters.
Thus, the equations of motion to determine variational trial functions for simultaneous full quantum electronic and nuclear wavepacket dynamics have been set up. Our next task is to prepare the trial functions.

Total wavefunctions in Gaussian representation
We next expand the entire electronic-nuclear wavefunctions in terms of the atomic-orbital-like Cartesian Gaussian functions, and only the nuclear basis functions are multiplied by plane waves. The variational parameters appearing below in those functions are to be determined variationally with the above dual least action principle.

The simplest case: nuclear packets running on an adiabatic potential
The starting total wavefunction is chosen to be in the simplest form as which fits in a physical situation before the electron wavepacket undergoes bifurcation. Φ (ad) K (r; R t ) is the Kth adiabatic electronic wavefunction at molecular geometry R t , which indicates the position averaged over the nuclear wavefunction G (R, t) such that where G (t) is normalized. G (R, t) in Eq. (90) is can be primarily written as in which the so-called contracted Cartesian Gaussians are adopted with plane wave components in such a way that where R = (R A , R B , ...) with R A = (x A , y A , z A ), and the vector Λ A indicates a set of (l A , m A , n A ) as in the atomic orbital representation.
ΛA requires a summation with respect to possible combinations of Λ A within the coordinate R A . The function g A (R A , t) for the nuclear coordinate R A is supposed to be localized around a nuclear Gaussian center Q A (t), and Q A (t) are supposed to be moving with its effective associated momentum P A (t), as in the sense of the squeezed state representation. The (complex-valued) exponents matrix α A ΛA (t) give the width of the Gaussian, each of which is a little skewed by pre-exponential factor with (l A , m A , n A ) and coefficients d A ΛA (t) . Accordingly, the elementary setup of the Gaussian functions as basis functions is almost exactly the same as those in the electronic state calculations, except that the nuclear part is characterized by the time dependence of P A (t) and α A ΛA (t). Thus, the molecular frame to define electron dynamics is specified by the quantum mechanical average of the nuclear positions with the average three-dimensional position for the nucleus A being under conditions g A | g A = 1. Practically, a convenient approximation for the averages must be to choose which is practically identified as the central position of the Ath nuclei at time t.
The squeezed-state-type Gaussian functions, without the polynomial functions in the pre-exponential factors, have been extensively used mainly in semiclassical context in the Born-Huang representation (see a series of pioneering works by Heller [97,98]). As discussed in the preceding section, we have taken the normalized variable Gaussians as ADF in the nonadiabatic electron wavepacket dynamics. [82,83,86] However, we hereafter proceed in a full quantum manner without resorting to any classical and/or semiclassical dynamics. The greatest advantage of the functions of Eq. (93) in use of the molecular problems is that all the atomic integrals with respect to the nuclear kinetic energy, electron-nuclear Coulomb interaction, nuclear-nuclear Coulomb interaction are basically available. These integrals are direct extensions of the ab initio electronic state counterparts; the energies of electronic kinetic, electron-nuclear interaction, and electron-electron interactions. Besides, as stressed above (in Sect. 2.3.2), much experience about the atomic and molecular Gaussian integrals including plane waves has been extensively accumulated in the studies of electron scattering by molecules. [18,30,35] Incidentally, the electronic part of Ψ(r, R, t) in Eq. (90) can be expanded in general basis functions other than the adiabatic wavefunctions such that Ψ(r, R, t) = I C I (t)Φ I (r; R t ) G (R, t) , (97) where Φ I (r; R) are those functions like configuration state functions, Slater determinants composed of any orbitals, and so on. By variationally propagating C I (t) in time simultaneously with G (R, t), one can track Ψ(r, R, t) without diagonalization for obtaining Φ (ad) K (r; R t ) at each time. [70] As for the electronic counterparts, we also adopt the standard Cartesian Gaussian functions χ A i (r; R A ), which is positioned at R A . Relevant molecular orbitals φ m (r) are expanded as with The total electronic and nuclear wavefunctions are expanded in terms of the configuration functions where A is the antisymmetrizer for the electronic coordinates only. Note again that Φ I (r; R t ) is just one of the configuration state functions or the Slater determinants adopted.

Nuclear molecular orbitals and beyond
It is not hard to extend the localized nuclear Gaussian functions so as to extend in space to be referred to as "nuclear molecular orbital" such that By assigning g 1 (R, t) to R A and g 2 (R, t) to R B , and so on we obtain a Hartree product in terms of time-dependent nuclear molecular orbitals. Compared with Eq. (100), with this straightforward extension the electronic and nuclear parts are treated on equal footing, except for the antisymmetrization over electronic coordinates. Generalization to the multi-reference form further consisting of the direct products of those nuclear molecular orbitals is rather straightforward.

Configurations after wavepacket bifurcation
A general form of a electronic and nuclear electronic wavepacket can be represented such that (103) We here choose the adiabatic functions Φ (ad) K just for simpler interpretation. One may use more general basis functions like the Slater determinants as done in Eq. (97). G K (R, t) indicates the nuclear wavefunction in the Kth nuclear frame defined as with R K t being the average over the Kth frame. Each component such as g K A (R A , t) is further expanded as where Λ K,A = (l K A , m K A , n K A ). These are rather a straightforward extension of the definition of Eq. (93). One may want to use the nuclear molecular orbital of Eq. (101) for the individual Kth state. Ψ(r, R, t) in Eq. (103) is the most general form of the electronic and nuclear wavepacket dynamics in this stage.

Dynamical parameters to be determined
Almost all the matrix elements needed for the atomic integrals for electronic and nuclear components are essentially available in the standard quantum chemical programs, even though some necessary modifications are demanded. The parameters to be tracked in TDVP are summarized as and for nuclei (and D M A (t) for nuclear molecular orbitals), and for electrons. Note that the variational parameters for electronic configurations are just linear and thereby can be determined with the Dirac-Frenkel variational principle in the manner similar to Eq. (38).

Mean field approximation
The wavefunction in Eq. (103) can be determined variationally for all the time throughout reactions, which represents a coherent set of electronic and nuclear wavepackets, each running on each (the Kth) adiabatic potential surface at the center of the nuclear positions Q K A (t), Q K B (t), · · · . Suppose the total wavepacket Ψ(r, R, t) lies in an area of strong nonadiabatic coupling among the component states (say the Kth adiabatic states). In such a case, it should be impractical to calculate the Hamilton matrix elements between electronic wavefunctions sitting at Q K A (t), Q K B (t), · · · and that at Q L A (t), Q L B (t), · · · for K = L even with the current technology of quantum chemistry.
To cope with this situation, we consider a meanfield approximation with respect to the centers of the nuclear wavepackets, which shares the same spirit with the semiclassical Ehrenfest theory, in which the nuclear force is averaged over the relevant electron wavepackets (recall Eq. (41)). We below describe how to propagate in time the mean-field wavefunction Ψ MF (r, R, t) in a recursive manner. Suppose we have already had Ψ MF (r, R, t) such that All the electronic states have a common molecular frame R AV t , which is With the initial condition the Kth states. Then we take average only for Q and P parameters in such a way that and We thus can proceed to determine the parameters at time step next to t + Δt. Note here again that the mean-field approximation is generally not valid for an exceedingly long-time.
Incidentally, one can track back to the semiclassical Ehrenfest theory (SET) [77][78][79][80][81] from the form of Eq. (109), which is shown in Appendix, justifying the mathematical procedure of the present framework as well.

Wavepacket bifurcation
We next consider the quantum mechanical extension of the idea of naturally branching paths [76], which was devised for the nonadiabatic electron wavepacket dynamics of Eq. (37) in the path approximation with the force matrix. Let us consider a wavepacket bifurcation Ψ K (r, R, s) = Φ in which a wavepacket Ψ K (r, R, s) passes across a nonadiabatic region on the way from time s to t.
Suppose Ψ K enters into the nonadiabatic region. Then we switch to the mean-field approximation within the adiabatic states involved in this nonadiabatic region with an appropriate set of the conditions to resume the mean-field wavepacket dynamics. are approximately satisfied, we terminate the meanfield approximation, whereĤ is the total electronic and nuclear Hamiltonian. After time t B , each wavepacket component resumes to run on each adiabatic state such that of Eq. (105). C K (t B ) is the final coefficient at t = t B , and thereafter C K (t) is fixed as giving the transition amplitude. By the collection of those bifurcated wavepackets, each being represented as in Eq. (116), into a single grand wavefunction, the form of Eq. (103) is retrieved.
Another physically important aspect of wavepacket dynamics is the coherent merging between the electronic and nuclear packets running on different adiabatic potential energy surfaces in nonadiabatic interaction regions. In this case the total Hamiltonian matrix elements in Eq. (115) grows to a significant level. Here at this point the simple path approximation discussed in Sect. 3.3 faces a fundamental difficulty due to the lack of the information about nuclear wavefunction. The present approach of full quantum wavepacket treatment has lifted this fundamental adversity. We will discuss the relevant technical issues somewhere else.
We have just formulated a nonadiabatically coupled dynamics in terms of simultaneous electronic and nuclear quantum wavepackets.

Concluding remarks
We have partially accomplished the aim to formulate simultaneous electronic and nuclear quantum wavepacket dynamics in molecules by incorporating nuclear wavepackets into the previously formulated nonadiabatic electron wavepacket dynamics. A logical trail to reach there has been sketched together.
In the proposed scheme, the nuclear wavepackets are to be expanded in terms of the (Cartesian) Gaussian functions, exactly as the electronic wavefunctions are. A formal difference is that the nuclear wavefunctions need plane-wave components as in the so-called squeezed-state (or vaguely coherent-state) representation. The accumulated experience in the theory of electron-molecule scattering and others, the atomic and molecular integrals including those plane waves, has already been customarily performed. Therefore, it is anticipated that the extension of the molecular integrals is of no technical difficulty.
The positions of the center of the nuclear Gaussians are determined through the phase-space-like parameters involved in the squeezed-states, and the spatial distributions and anisotropy of the shapes are to be represented by the Gaussian exponents and the preexponential polynomial terms as well. Those parameters, along with the relevant coefficients of electronic molecular orbitals and electronic configurations, are to be uniformly and systematically determined by the quantum mechanical time-dependent variational principle. Among many variants of the time-dependent variational principle (TDVP), we propose to apply our developed one, the dual least action principle based on the quantum mechanical Maupertuis-Hamilton least action principle, which has resolved the mathematical difficulties inherent to the existing TDVP.
If the nuclear wavepacket dynamics is mimicked with the path-branching dynamics, the theory comes back to the mixed quantum-classical representation as presented in Sect. 3.4. We have shown there that the variable Gaussian functions (ADF-NVG) that are designed to run on the branching paths could approximate the corresponding full quantum nuclear wavefunctions very accurately. Thus, the proposed full electronic-nuclear simultaneous wavepacket dynamics is well guaranteed in accuracy. We anticipate much development in new realm of molecular science, in which quantum nature of nuclei and entanglement between electrons and nuclei are significant.
Finally, the present work does never ever underestimate the roles of the Born-Huang representation based on the molecular view of Born and Oppenheimer. Indeed, we have emphasized this aspect through our own contributions in Sect. 2. Nonetheless, it is worth refraining that the simultaneous nonadiabatic electron and nuclear wavepacket representation should be indispensable to analyze the molecular systems of very large electronic fluctuation and/or to track real-time electronic and nuclear coupled dynamics of significant quantum effects.
where both the nuclear and electronic parts are normalized such that (A7) Then it appears as Ψ(t) Ĥ (R, elec) Ψ(t) Next, let us classicalize them such that Then we have Ψ(t) Ĥ (R, elec) Ψ(t) Applying to the Dirac-Frenkel variational principle with respect to the set of coefficients {C I (t)}, we directly attain the coupled equations of motion for the electronic states. Note that the second-order terms 2 2 k Y k IJ exist in the dynamics [1][2][3][4][5][6][7], in contrast to the previous theories commonly referred to as SET [77][78][79][80][81].