Electronic and nuclear fluxes induced by quantum interference in the adiabatic and nonadiabatic dynamics in the Born-Huang representation

We perform an electronic and nuclear ﬂux analysis for nonadiabatic dynamics and its corresponding adiabatic counterpart, both of the wavefunctions of which are represented in the Born-Huang expansion. It is well known that the electronic-nuclear con-ﬁgurations (terms) in the expansion of the total wavefunction interfere each other through the nonadiabatic interactions and give birth to electronic and nuclear ﬂuxes. Interestingly, even in the adiabatic dynamics without such nonadiabatic interactions, a wavefunction composed of more than one adiabatic state can undergo interference among the components and give the electronic and nuclear ﬂuxes. That is, the individual pieces of the wavepacket components associated with the electronic wave-functions in the adiabatic representation can propagate in time independently with no nonadiabatic interaction, and yet they can interfere among themselves to generate the speciﬁc types of electronic and nuclear ﬂuxes. We refer to the dynamics of this class of total wavefunction as multiple-conﬁguration adiabatic Born-Huang dynamics. A systematic way to distinguish the electronic and nuclear ﬂuxes arising from nonadiabatic and the corresponding adiabatic dynamics is discussed, which leads to the deeper insight about the nonadiabatic dynamics and quantum interference in molecular processes. The so-called adiabatic ﬂux will also be discussed.


I. INTRODUCTION
Nonadiabatic transition is one of the key phenomena in molecular science in which chemical and physical properties of molecules can undergo sudden change during the dynamical processes.2][3][4][5] (We refer to Refs. 2, 6, and 7 for full quantum mechanical calculations of the nonadiabatic coupling elements and to Refs. 8 and 9 for recent progress in quantum mechanical wavepacket dynamics.)2][13][14][15][16][17][18] This is because nonadiabatic dynamics can be sensitively reflected in the electronic flux induced by the dynamics.
The quantum mechanical flux, electronic and nuclear ones, are quite sensitive to the dynamical change of the associated wavefunctions, and therefore they appear to be vital and characteristic means to extract insights about chemical dynamics.The graphical representation of the motion of radical electrons, negative charges, and rearrangement of chemical bonds are widely adopted in the textbook of chemical reactions, which is very useful to facilitate our intuitive understanding of complicated chemical reactions.Interestingly, however, there have not been many quantum mechanical studies until recently which are devoted to The Journal of Chemical Physics ARTICLE scitation.org/journal/jcpjustifying or correcting such an intuitive view of electronic flow associated with chemical reactions.Among others, the group of Manz has extensively performed a flux analysis on the pericyclic reactions 19 and on intramolecular concerted reactions in terms of electronic and nuclear fluxes. 20They have also shown electronic flow induced by crossing the transition state. 21Refer also to Ref. 22 for the simultaneous electronic and nuclear fluxes.Flux analyses have been studied for nonadiabatic processes 23 and charge migration dynamics in systems such as HCCl, 24 benzene, 25 and LiH. 26katsuka and his co-workers have revealed the electrondynamical mechanism of proton transfer 11 and double proton transfer 12 in terms of electronic flux, which was taken from the complex-valued electron wavepackets.In these studies, they have shown that the flux analysis is quite useful, particularly for the study of nonadiabatic chemical reactions and dynamics.However, the flux analysis made by them has been based mostly on the mixed quantum (electronic) and classical (nuclear) representation of the total molecular wavefunctions, in which the electronic wavepackets are described as a complex-valued function but the nuclei are treated "classically" (but not necessarily Newtonian).2][13][14][15][16][17][18] Therefore we herein concentrate on a canonical flux analysis on the nonadiabatic dynamics in terms of the wavefunctions in the Born-Huang expansion, which is the theoretically standard representation of the total wavefunctions.[In the Appendix, we briefly summarize the flux arising from the Ab Initio Molecular Dynamics (ABMD) and the semiclassical Ehrenfest theory (SET) since these are among the most frequently applied methods in the analysis of chemical reactions.]Matsuzaki and Takatsuka 27 have recently applied a full flux analysis for a wavepacket bifurcation process due to the nonadiabatic interactions in the Born-Huang representation by calculating both the electronic and nuclear fluxes.They first prepared a single nuclear Gaussian wavepacket at a Franck-Condon region on the first excited potential energy curve of LiF in the adiabatic representation and let it proceed towards the avoided crossing numerically.This single-configuration wavepacket splits into two pieces, one remaining on the second potential curve (representing Li + F − ) and the other running away to the dissociation channel as Li and F atoms (ground state).The characteristic features of the relevant fluxes such as the spatiotemporal oscillation of the electronic flux have been observed, which facilitates a deeper understanding of the nonadiabatic dynamics.They also studied the Rabi-like oscillation in the flux based on the two-state interaction model to account for a part of the temporal oscillation of the flux.A correlation function between the nuclear and electronic fluxes has been numerically investigated for the first time, along with the physical interpretation behind.9][30][31][32][33][34] This nonadiabatic dynamics has been compared with the corresponding adiabatic dynamics (with no nonadiabatic coupling elements on the exactly the same potential surfaces), in which the single-configuration wavepacket remains on the same adiabatic potential curve.In this adiabatic dynamics, no electronic flux is generated throughout since the electronic wavefunction is kept real-valued.However, the above study is never all about the comparison between the adiabatic and nonadiabatic dynamics, since the wavefunction they studied as a Born-Oppenheimer adiabatic wavefunction is only of a single-configuration form χ I (R, t)Φ I ( q ; R) in the expansion of Eq. (5).However, the Born-Oppenheimer approximation is not limited to the form of the single-configuration, and the multiple-configuration adiabatic wavefunction, which is composed of plural electronic and nuclear wavefunctions as in Eq. ( 5), can bear the totally different characteristics in the flux behavior even under the absence of nonadiabatic interactions.Such a multiple-configuration adiabatic wavefunction has been studied by Bandrauk et al., 35 who promoted the ground state of H + 2 to the plural electronically excited states coherently with an ultrashort pulse laser.Their finding is remarkable; a class of charge migration in this molecule has been induced due to the quantum coherence despite the absence of any mutual mechanical (nonadiabatic) interactions.This suggests a novel aspect of the Born-Oppenheimer adiabatic dynamics.
As will be shown later, the flux property of the multipleconfiguration adiabatic state is quite different from that of the single-configuration adiabatic counterpart. 27Therefore, when we discuss quantum coherence and decoherence in the nonadiabatic dynamics, we need to explicitly specify whether an adiabatic wavefunction to be compared is of the singleconfiguration or multiple one.This paper therefore complements our former study on the nonadiabatic dynamics 27 in that the scope of the adiabatic dynamics is extended to the multiple-configuration Born-Huang representation.
We begin with the formal study in Sec.II and define the various fluxes based on the adiabatic Born-Huang representation, deriving the relevant identities of them.We will show that there are flux components that can distinguish singleconfiguration adiabatic, multiple-configuration adiabatic, and nonadiabatic wavefunctions.In Sec.III, we prepare a LiF system to verify the theoretical objects.The model system is designed so that the magnitude of the nonadiabatic interaction can be artificially varied from the genuine value.We numerically present in Sec.IV the three major fluxes (electronic, direct nuclear, and indirect nuclear fluxes) and their reduced counterparts.The correlation function between the nuclear and electronic flux will also be presented.The paper concludes in Sec.V.

II. ELECTRONIC AND NUCLEAR FLUXES FOR MOLECULAR WAVEFUNCTIONS IN THE BORN-HUANG EXPANSION
We first present the formal expressions of molecular flux and its components of three types, based on the Born-Huang J. Chem.Phys.150, 014103 (2019); doi: 10.1063/1.5066571150, 014103-2 representation.This section is complementary to the literature, for example, Refs.20, 22, 23, and 36, and a particular emphasis is placed on the distinction between the fluxes arising from a nonadiabatic wavefunction and the corresponding multiple-configuration adiabatic function.

A. Flux in the Born-Huang expansion
As usual, we begin with the time-dependent Schrödinger equation with the canonical Hamiltonian Here q = {r, ω } with q i = (r i , ω i ), where r i is the coordinate of the ith electron and ω i being the spin variable.{r} and R in V({r}, R), representing the Coulombic potentials among electrons and nuclei, collectively denote the electronic and nuclear coordinates, respectively.m e is the electron mass and M denotes a collective notation of the nuclear masses M k .The kinetic energy operators are assumed to have the canonical forms such that and In solving Eq. ( 1), we resort to the Born-Huang expansion We do not have to specify Φ I ( q ; R) at this point whether they are in the adiabatic representation or else, and the total electronic and nuclear flux as a physical quantity should be invariant with respect to the choice of the representation of wavefunctions.However, the components, such as nuclear flux and electronic flux, can individually depend on the functional form of the total wavefunction used.Therefore we consistently adopt the adiabatic representation for the flux quantities (except in the practical calculations of nuclear nonadiabatic wavepackets with the use of "diabatic representation"), and the electronic wavefunctions are all supposed to be realvalued.Equation (1) then gives rise to the coupled equations of motion for the nuclear wavepackets where Pk indicates the nuclear momentum operator in the kth direction, Pk = −i ∂/∂R k , and the nuclear kinematic coupling elements are defined as and H el IJ are the matrix elements of the electronic Hamiltonian The terms Y k IJ are associated with 2 and are supposed to be small.Therefore the terms X k IJ in the adiabatic representation are supposed to be responsible for the nonadiabatic transition throughout this paper.
As usual, we consider the conservation of the quantum probability where denotes the integration over the electronic coordinates from q 2 to q N (except q ≡ q 1 ) in their entire spaces, with N being the total number of electrons.Ω r and Ω R are designated to be the spaces for the electron of r 1 and the nuclei to be contained, respectively.∂Ω r and ∂Ω R are the surfaces surrounding Ω r and Ω R , respectively, and dS r and dS R are the volume (surface) elements of the integration, which are normal to ∂Ω r and ∂Ω R , respectively.With no loss of generality, one can rigorously define the flux where the electronic flux is defined as with and The significance of this difference should manifest itself in the chemical reaction dynamics of open-shell electronic states.In this paper, we study only the spin singlet case, and we refer to the electronic flux to that defined in Eq. ( 12), which use only r without q.The nuclear flux J nu (r, R, t) + J el nu (r, R, t) is defined as The total flux is thus given as which is supposed to be integrated only individually as in Eq. (11).We assign the lower case j and upper case J to the electronic and nuclear fluxes, respectively.In both J nu (r, R, t) and J el nu (r, R, t), the nuclear derivative operator ∇ R is operated on the total wavefunctions Ψ({r}, R, t), and therefore they are collectively called the nuclear flux.However, it should be noted that while ∇ R in J nu (r, R, t) is operated on the nuclear wavefunctions, that in J el nu (r, R, t) is operated on the electronic wavefunctions.[The definitions of J nu (r, R, t) and J el nu (r, R, t) will be described below.]Hence the "nuclear flux" J el nu (r, R, t) emerges from the spatial dynamics of electronic wavefunctions in the similar way to the matrix elements X k IJ of Eq. (7).We survey the individual flux components below.

Electronic flux j el (r, R, t)
The electronic flux j el (r, R, t) is explicitly written by inserting the total wavefunction in the form of the Born-Huang expansion into Eq.(12) as where we have defined an electronic flux operator ˆ el such that with ∇ r being supposed to be operated on the function of the coordinate r .The off-diagonal electronic density is given as and ρ el JI (r, R) = ρ el JI (r, R; r, R).Note in the double summation of Eq. ( 18), the diagonal terms are not included since ˆ el ρ el JI (r, R) takes a large value only when the off-diagonal density (transition density in the terminology of quantum chemistry) ρ el JI (r, R) contains the asymmetric momentum information with respect to the electronic motion.j el (r, R, t) can take a significant value when not only the component ˆ el ρ el JI (r, R) but also the nuclear differential overlaps χ * J (R, t)χ I (R, t) are large enough.
To facilitate the intuitive understanding of the physical properties of j el (r, R, t), we will use the following reduced fluxes later in the numerical analyses.These are f R (R, t) = j el (r, R, t)dr and f r (r, t) = j el (r, R, t)dR.
Their physical meanings are rather obvious.

Direct nuclear flux J nu (r, R, t)
The direct nuclear flux is defined and explicitly written as where we have Similarly to j el (r, R, t), J nu (r, R, t) can take a large value when not only the pure nuclear flux component ˆ nu χ * I (R, t)χ I (R, t) but also ρ el II (r 1 , R) happens to have a significant value simultaneously.Note that the diagonal terms in the representation of Eq. (25) do not vanish in J nu .
We define a reduced flux of J nu , by integrating J nu (r, R, t) over r as The Journal of Chemical Physics ARTICLE scitation.org/journal/jcp

Indirect nuclear flux J el nu (r, R, t)
The other component of the "nuclear flux" is defined as where the operator ˆ nu is defined as with ∇ R and ∇ R being supposed to be operated on ρ el JI (r, R ; r, R), respectively.Due to the identity we have no sum over the diagonal elements, and consequently it becomes This flux emerges characteristically from the Born-Huang expansion, in which the electronic wavefunctions parametrically contain the nuclear coordinates.As the nuclear geometry changes a little, the electronic state can dramatically change as in the nonadiabatic transition, which in turn drives (accelerates or decelerates) the nuclear flux.(Note that increase in the nuclear flux does not necessarily imply that the corresponding classical nuclear velocity becomes larger.)Again χ * J (R, t)χ I (R, t) should be simultaneously large, too.Notice that the smaller coefficient /2iM, rather than /2im e , is associated with J el nu (r, R, t).We also note that j el (r, R, t) and J nu (r, R, t) + J el nu (r, R, t) are individually invariant with respect to the correct transformation among the basis functions Φ I ({r}; R) as long as the total wavefunction Ψ({r}, R, t) is kept invariant.However, each of J nu (r, R, t) and J el nu (r, R, t) can depend on the choice of the representation.
We define the reduced flux of J el nu (r, R, t) by integrating over r such that where with its kth component being X k IJ in Eq. ( 7), and we thus have Both F nu (R, t) and F el nu (R, t) indicate the nuclear flux without respect to the electronic position.

B. Flux in the Born-Oppenheimer approximation
The total wavefunctions Ψ({r}, R, t) undergoing nonadiabatic transitions can have finite values for all the flux components j el (r, R, t), J nu (r, R, t), and J el nu (r, R, t), whereas in the Born-Oppenheimer approximation, some of these quantities become identically zero, depending on how the total wavefunction is represented mathematically.Thus the nonadiabatic dynamics can be characterized in terms of the flux.

Single-configuration BO representation
Suppose we have a Born-Oppenheimer state of a single configuration As far as the ground state of a stable molecule is concerned, this functional form usually gives an accurate approximation to the full Born-Huang expansion. 37Since Φ I ({r}; R) is an eigenfunction of the electronic Hamiltonian and is usually real-valued (even the degenerate wavefunctions can be transformed to be real), it simply gives Also, Eq. ( 30) indicates that Thus only the direct nuclear flux can have a finite value such that which gives a clear interpretation that J nu (r, R, t) carries the electronic density at R with the pure nuclear flux ˆ nu χ * I (R, t)χ I (R, t) .Since only J nu (r, R, t) is not identically zero, the so-called adiabatic "electronic" flux, if any, must be extracted from J nu (r, R, t).
If a chemical reaction proceeds in a concerted manner in that the total wavefunction could be represented ideally in an adiabatic single-configuration Born-Huang representation, the net "electronic flux" to be induced should be zero all the way through the reaction, as in Eq. ( 35).This is the physical background of the minimum flux principle 12 asserting that the concerted reactions such as the Woodward-Hoffmann allowed reactions will proceed inducing the electronic flux to the minimum extent possible.

Multiple-configuration BO representation
Another but less frequently studied types of the Born-Oppenheimer wavefunctions are those composed of multiple configurations, which is under our current attention.Shining a ultrashort pulse laser (in the attosecond scale) on a molecular ground state χ 0 (R, t)Φ 0 ({r}; R), one can coherently excite it with a large energy-time uncertainty to a multipleconfiguration state 35 such that If there is no nonadiabatic couplings among the electronic states, those nuclear states χ I (R, t) thus created should run on the individual potential energy surface offered by Φ I ({r}; R).This is exactly what we mean by "multiple-configuration BO state."The time-propagation of the entire wavefunction of Eq. ( 38) is readily performed by propagating each of χ I (R, t)Φ I ({r}; R) independently.
In clear contrast to the single-configuration BO state, however, its flux calculation is more involved, since it holds that j el (r, R, t) 0 (39) besides Thus, when at least two of the components, say, χ I (R, t)Φ I ({r}; R) and χ J (R, t)Φ J ({r}; R) come close in space, they can interfere with each other and may generate a new flux component for j el (r, R, t) and/or J nu (r, R, t) even in the absence of the nonadiabatic coupling elements X k IJ .
Nevertheless, J el nu (r, R, t) is zero or extremely small.This identity is crucial since J el nu (r, R, t) 0 in the corresponding nonadiabatic dynamics.

C. Quantum beat in the flux
Vitanov and Garraway have found an oscillatory behavior of the nonadiabatic transition probabilities through their studies of the so-called transition time in nonadiabatic dynamics in the Landau-Zener model. 38,39We have reported spatiotemporal oscillations in the fluxes before. 27Such oscillatory behaviors are not always identified, but if found, they may reflect a specific dynamical feature behind.We here formulate two types of oscillatory patters in the fluxes, one being the simple quantum beat and the other arising from the Rabi-like oscillation among the states that are coupled by the nonadiabatic interactions.

Effect of the electronic phases
The flux can exhibit a quantum mechanical oscillatory pattern even in the absence of nonadiabatic transitions.Suppose a ground state vibronic wavefunction is initially photoexcited to various adiabatic states as as in a multiple-configuration adiabatic Born-Huang representation.The time-evolution of this wavefunction, with or without nonadiabatic coupling elements, may be represented as A comparison of this total wavefunction with that of Eq. ( 5) shows in which the electronic phases are explicitly taken account.
Then the flux components are accordingly rewritten as follows: The electronic flux is and f R becomes where Thus both j el (r, R, t) and f R (R, t) can have a quantum beat due to the presence of the factor On the other hand, possible oscillatory behavior in f r (r, t) is a little complicated.If the product of the nuclear wave packets χ 0 * J (R 0 (t), t)χ 0 I (R 0 (t), t) happen to be well localized at the position R 0 (t) for a sufficient time, f r (r, t) can be approximately represented as which suggests that f r (r, t) may retain a possible oscillatory property.However, this should be rarely the case in practical cases.
The direct nuclear flux reads The Journal of Chemical Physics

ARTICLE scitation.org/journal/jcp
and its reduced version is where we have used It is interesting to note that F nu (R, t) bears the factor which is proportional to the "acceleration" of the nuclear motion or the "force" working on the nuclei running on the potential energy surface E I (R).
The indirect nuclear flux looks as along with its reduced flux This beat, and that in j el (r, R, t) as well, fades away as Since the origin of the oscillatory feature in j el (r, R, t), J nu (r, R, t), and J el nu (r, R, t) comes from exp 1 i (E I (R) − E J (R))t , it has been integrated out from the reduced flux.No beat is expected of course from the single-valued Born-Huang representation.

Rabi-like oscillation in the nonadiabatic dynamics
In case where nonadiabatic interactions are present, the electronic states involved are given additional quantum phases, which are reflected in the electronic flux j el (r, R, t).To see this in the simplest way, we approximate the total wavefunction [Eq.(42)] in the form of the semiclassical Ehrenfest theory only in a narrow area nearby an avoided crossing where G(R, t; R(t)) is supposed to be the relevant nuclear wavefunction like a sharp Gaussian function, the center of which is located at R(t), and c I (t) are the coefficients to represent the electron dynamics as above.This approximation is justified since the nuclear wavefunction varies much slower than the electronic dynamics represented by c I (t) .In the case of the two level problem, N st = 2, we have where v(t) is the velocity of the center of the nuclear wavepacket G(R, t, R(t)).These coupled equations of motion are isomorphic to those of the two-state model quantum dynamics in an intense laser field, the most well-known phenomenon of which is the Rabi oscillation.Similarly, the nonadiabatic dynamics can be accompanied by oscillations in the populations and phases.Matsuzaki and Takatsuka have studied this aspect in detail and gave the angular frequency of the oscillation. 27This type of Rabi-like oscillation can compete with the simple quantum beat discussed above.When the nonadiabatic coupling element is small enough, the frequency of the beat is to be determined mainly by the energy gap.The significance of the Rabi-like oscillation in the study of nonadiabatic dynamics has been discussed in Ref. 27.Note that for the nonadiabatic dynamics, it is practically impossible to separate the effects of the simple beat and the Rabi-type oscillation in the oscillatory structure of the fluxes.

D. Summary
The above properties of the components of electronic and nuclear fluxes are summarized in Table I.

III. ELECTRONIC AND NUCLEAR WAVEPACKET DYNAMICS OF LIF MOLECULE IN TWO-STATE MODEL
We perform numerical studies to see how we can or cannot distinguish the nonadiabatic dynamics and the multipleconfiguration adiabatic dynamics in terms of the quantum flux components.To do so, we perform the full quantum mechanical calculations over the electronic and nuclear degrees of freedom for the two-state model of the spin-singlet LiF molecule in the molecular frame only, surveying only the intramolecular flux as a case study.LiF molecules in the excited state are well studied in many aspects. 7,457][48][49][50] Although the flux in the corresponding laboratory frame is directly relevant to experimental observation such as photoelectron spectroscopy and cross sections of chemical reactions, we do not consider these aspects in this paper.No molecular rotation is considered either, and the nuclear coordinates so far generally denoted as R are reduced to only the internuclear distance Q in what follows.Thus the total wavefunction in the Born-Huang expansion to be numerically studied is where Φ ad I ({r} ; Q) are the adiabatic electronic wavefunctions.
In our previous flux analysis on the nonadiabatic dynamics compared with the single-configuration adiabatic dynamics, we considered only a wavepacket bifurcation. 27Here in this paper, we study the dynamics of a wavefunction starting in the multiple-configuration Born-Huang representation at the outset.A pair of these wavefunctions is compared; one undergoes nonadiabatic transition and the other does not.To see the effect of the magnitude of the nonadiabatic coupling element on the dynamics and the resultant flux, three sets of adiabatic potential curves and the nonadiabatic coupling elements are prepared, one is the native set attained with the quantum chemical calculations and the other two counterparts are rather artificial in that the nonadiabatic coupling elements are deliberately made larger and smaller than the native one.
While the calculations of quantum wavepackets and the potential curves along with their coupling elements have been consistently carried out in the diabatic representation, we represent all the resultant fluxes in the adiabatic representation throughout this paper.

A. Potential surfaces and coupling elements
The adiabatic ground state (Φ ad 1 ) and the first excited state (Φ ad 2 ) of the LiF molecule in spin-singlet are first obtained with the GAMESS quantum chemistry package, 51,52 using the two-state averaged complete active space self-consistentfield (SA-CASSCF) calculation of the level of a 6-electron 6orbital active space with the use of aug-cc-pVDZ basis sets at every 0.1 Å internuclear distances ranging from 0.5 Å to 11.0 Å.We have numerically calculated the nonadiabatic coupling elements within the SA-CASSCF scheme by means of the finite difference method, where • | • r means integration over the electron coordinates.Since all the quantum wavepacket dynamics calculations are performed in the diabatic representation, we adopt the standard method 2 using the adiabatic-to-diabatic transformation angle α(Q) Then the approximate diabatic electronic states Φ I ({r} ; Q) (I = 1, 2 and with no superscript) are attained by the transformation such that ARTICLE scitation.org/journal/jcpFIG. 1. Schematic representation of the potential curves in the diabatic representation and the associated electronic Hamiltonian matrix elements H el IJ (Q).Two wavepackets (A) and (B) constitute an initial total wavefunction both for the adiabatic and nonadiabatic dynamics.In the nonadiabatic dynamics, the packets (A) and (B) merge into a single piece (C), while in the adiabatic dynamics, each of (A) and (B) is propagated in time independently on the individual "adiabatic potential curves" without the nonadiabatic interaction.
Figure 1 shows the diabatic potential energy curves and the electronic matrix elements H el 11 (Q) represents the electronic energy of the covalent state, while H el 22 (Q) of the ionic character.These states cross at Q = 6.1 Å, with the coupling elements being referred to as H el 12 (Q) in Fig. 1.
The nuclear wavepackets in this diabatic representation χ I (I = 1, 2) are obtained by integrating the coupled Schrödinger equations We integrate the coupled equations ( 67) and (68) numerically with the split-operator FFT method. 40The spatial grid points as many as 2048 have been taken between 0.0 Å and 20.0 Å and the time step has been set to 25 attoseconds.Since the ab initio potential curves have been generated only in the interval from 0.5 Å to 11.0 Å, we have installed an absorbinglike potential beyond the interval.Hence, no information on the fluxes beyond the interval is given.The nuclear wavepackets χ I (Q, t) thus attained are transformed into the adiabatic counterparts χ ad I (Q, t) to calculate the relevant fluxes.The total wavefunction Ψ({r}, Q, t) in Eq. ( 62) is kept invariant throughout the transformation.

B. Modification of nonadiabatic coupling
To investigate the relationship between the magnitude of the nonadiabatic couplings and the resultant fluxes, we introduce the following method that systematically scales the nonadiabatic coupling elements.The diabatic Hamiltonian matrix elements H IJ and adiabatic potential curves E I can be transferred to each other Since the integral of X 12 in Eq. ( 63) must be constant, we use a transformation that makes the integral of X 12 invariant.There can exist many such transformations, we adopt a scaling of the coordinates such that where Q 0 is the location satisfying α(Q 0 ) = π/4.With this transformation, we define the corresponding scaled potential curves Ēad I and the nonadiabatic coupling element X12 from Eq. ( 70) along with the off diagonal Hamiltonian H12 as The value of α(Q 0 ) = π/4, which gives H 11 − H 22 = 0, has been chosen to avoid the singularity due to cos 2 ᾱ = 0 in Eqs. ( 73) and ( 74) at Q = Q 0 .We use σ = 5 Å and A = 0.5, 1.0, and 2.0 in the below calculations.Figure 2 shows X12 and H12 as a function of the internuclear distance Q.
By varying the above parameter A in Eq. ( 72) as A = 0.5, 1.0, and 2.0, we prepare three sets of the potential curves Ēad 1 (Q) and Ēad 2 (Q) with and without the associated coupling elements X12 (Q).The choice of A = 1.0 among them gives the original potential set for LiF.

C. Wavepackets to be surveyed
The main subject of the present study is the quantum interference effect arising in two dynamics: one undergoing nonadiabatic transition and the other being adiabatic all the way through, but they commonly have the multipleconfiguration form [Eq. ( 62 ).We terminate the dynamics at this time.As a pair of two initial wavepackets, we let them move back numerically to the reversed direction with the inverse momenta, as indicated by the arrows in Fig. 1

on (A) and (B).
The nonadiabatic version of this total wavefunction consisting of (A) and (B) proceeds towards the asymptotic region, and in passing across the avoided crossing, (A) and (B) merge into a single piece like the packet (C).On the other hand, each of the components of the wavepackets (A) and (B) in the adiabatic counterpart proceeds individually on their adiabatic potential curves without being merged retaining each population.

FIG. 4.
Space-time distribution of the probability density Γ(r, Q, t) in (z, Q) space (at x = y = 0) for t = 20 fs (before penetrating into the avoided crossing region), 30 fs (just passing across the avoided crossing), and 40 fs (almost after the passage).The panels are arrayed in the same manner as in Fig. 3.

ARTICLE scitation.org/journal/jcp
The resultant nuclear paths are depicted in Fig. 3, which are actually represented in terms of the space-time trajectories of the nuclear wave packet densities running on the adiabatic potential curves for the ground state Ēad 1 (Q) (of blue back-ground) and the first excited state Ēad 2 (Q) (of brown back-ground).The sets in the top, middle, and bottom row of the panels exhibit the results in the cases of the weak coupling [A = 0.5 in Eq. ( 72)], the native one for LiF (A = 1.0), and the strong coupling (A = 2.0).In the left and right columns, the results of nonadiabatic (with the given X12 , see Fig.  Let us track the nuclear wavepacket dynamics in Fig. 3 on the native potential curves and nonadiabatic interaction for LiF (A = 1.0) (see the middle row and the left column).The initial wavepackets are about to run on both the potential curves (t = 0), and the nonadiabatic interaction merges them into a single-configuration wavefunction to run on the Ēad 1 (Q) alone.The two adiabatic counterparts (the middle row and the right column) start in the same manner at t = 0, but in the absence of the nonadiabatic interaction, these pieces do not merge and keep running on the individual adiabatic potential curves.The similar dynamics are seen in the bottom row for the case of the large value (A = 2.0).The figures at the top row (A = 0.5) show different features.The nonadiabatic dynamics in the left FIG. 5. Space-time distributions of the electronic flux e z • j el (r, Q, t) in (z, Q) space (at x = y = 0) at t = 20.0 fs, 40.0 fs, and 60.0 fs.The panels are arrayed in the same manner as in Fig. 3.
column is closer to the adiabatic dynamics due to the weak nonadiabatic interaction.The majority of the wavepacket component is visually confirmed only on Ēad 1 (Q).After the launch of the packets at t = 0, they merge here again on the Ēad 1 (Q).This global feature is very similar to the adiabatic dynamics shown in the right column.
Figure 4 displays the snapshots of the simultaneous electronic and nuclear distributions defined as at selected times for both the nonadiabatic dynamics (left panels) and the adiabatic counterparts (right panel).In Fig. 4, both the F and Li atoms are placed on the z-axis and located at the nuclear coordinates (x, y, z) = (0, 0, 0) and (0, 0, Q), respectively.These six panels are arranged in the same manner as in Fig. 3. Γ(r, Q, t) is not a particularly interesting quantity in this study.Nonetheless, we show Fig. 4 as a reference because (i) we want to stress graphically that Γ(r, Q, t) looks much more smooth than the corresponding flux to be extensively presented below, (ii) we need to confirm in Fig. 8 that the space-time distributions of J nu (r, Q, t) are indeed quite similar to those of Γ(r, Q, t), and (iii) despite the graphical similarity between Γ(r, Q, t) arising from the nonadiabatic dynamics (the left column) and adiabatic counterpart (the right column), the fluxes arising from them, J el nu (r, Q, t) and F el nu (Q, t), in particular, can definitely distinguish from each other after all.

IV. FLUX FROM THE ADIABATIC MULTIPLE-CONFIGURATION BORN-HUANG REPRESENTATION AND THE NONADIABATIC WAVEFUNCTION
We next present numerical studies on the various components of the fluxes of the wavefunctions described above, FIG. 6. Distribution of the reduced electronic flux f Q (Q, t) defined in Eq. (78).The panels are arrayed in the same manner as in Fig. 3. which are nonadiabatic, adiabatic of multiple-configuration, and adiabatic of single-configuration in the Born-Huang representation.This section is devoted to the numerical realization of the critical properties summarized in Table I, with a particular attention paid to the distinction between the nonadiabatic and the multiple-configuration adiabatic wavefunctions.
A. Electronic flux j el (r, Q, t) We first show in Fig. 5 the electronic flux j el (r, Q, t), arising from the six different dynamics corresponding to Fig. 4 at time 20 fs (before penetrating into the avoided crossing region), 30 fs (passing across the avoided crossing), and 40 fs (after the passage).The fluxes given by the nonadiabatic and adiabatic wavefunctions at t = 20 fs are essentially the same as each other, since they simply run on the same potential curves with no nonadiabatic interactions.In contrast, it is immediately noticed that the patterns of the two fluxes in the top row (A = 0.5) at t = 30 fs and t = 40 fs are significantly different from each other, despite their very close similarity in the space-time distribution of the densities (see Figs. 3 and 4).Hence it is confirmed that even if the nonadiabatic interaction is small, the flux can show somewhat a large difference.The flux quantities can thus magnify a small difference in the wavefunctions or densities, simply because the flux reflects the derivative properties and is sensitive enough to the spatial variation of quantum phases.
A major difference between the nonadiabatic and adiabatic dynamics at t = 40 fs for A = 1.0 and A = 2.0 is that the flux in the adiabatic dynamics extends further to the long domain in the z coordinate (the longer internuclear distance) than the nonadiabatic flux does.This aspect is better confirmed in Fig. 6, which shows a reduced quantity from j el (r, Q, t) as defined as (78) FIG. 7. Distribution of the reduced electron flux f z (z, t) defined in Eq. ( 79).The panels are arrayed in the same manner as in Fig. 3.The quantum beats are apparent in the two panels for A = 1.0 with X 12 = 0 (the middle row in the right column) and A = 2.0 with X 12 = 0 (the bottom row in the right column).
This quantity represents the electronic flux to be observed at the internuclear distance Q at time t, no matter where the electronic components lie.Striking patterns of the oscillatory features in f Q (Q, t) are first observed.The flux in red indicates that the overall electron flux at given Q and t is directed towards the positive direction from F towards Li and blue in the opposite direction.This temporal oscillatory feature in the adiabatic dynamics should partly reflect the quantum beat as discussed above in Eq. ( 45).
Most importantly, the electron flux in the adiabatic dynamics keep appearing in the longer time and longer distance Q.First, f Q (Q, t) for the nonadiabatic dynamics commonly disappear around t = 40 fs.This is because the two wavepacket components initially prepared merge to a single-configuration real-valued function due to the nonadiabatic interaction (recall Fig. 3), and the electronic flux for such a wavefunction is nullified as shown in Eq. ( 35) and Table I.Remarkably though, the multiple-configuration FIG. 8. Space time distribution of the direct nuclear flux J nu (r, Q, t) in (z, Q) space (at x = y = 0) at time t = 20.0 fs, 40.0 fs, and 60 .0fs.The panels are arrayed in the same manner as in Fig. 3. adiabatic wavefunction gives a finite value in the electron flux, even if it undergoes perfect adiabatic dynamics (under X12 = 0).Yet, the reduced electron flux f Q (Q, t) in the multipleconfiguration adiabatic wavefunctions also fade away around t = 80 for A = 1.0 and 2.0, and before t = 40 for A = 0.5.This indicates that the coherence between χ ad 1 (Q, t)Φ ad 1 ({r}; Q) and χ ad 2 (Q, t)Φ ad 2 ({r}; Q) there has become too weak to be monitored by the flux quantity.Our studies 5,27 and the literature [30][31][32][33][34] discussed the coherence-decoherence problem in the context of the nonadiabatic wavepacket bifurcation.However, it is now apparent that the coherence and decoherence among the wavefunctions should be considered also from the view point of the multiple-configuration Born-Huang representation, irrespective of nonadiabaticity.
Figure 7 shows another reduced electron flux, which is defined as (79) Figure 7 shows even more clearly the difference in the dynamics between the nonadiabatic and multi-configuration adiabatic time-propagation.The electronic flux diminishes beyond 40 fs in the left panels (nonadiabatic dynamics), since the two packets merge into a single piece due to the nonadiabatic interaction, while the adiabatic dynamics continues propagating beyond the avoided crossing region (without nonadiabatic interaction).Thus a new feature (the flux directed inward represented in blue color) appears in the adiabatic dynamics as long as the coherent interaction continues in the two wavepackets.Also, the graphs demonstrate that the fluxes are more intensively induced as the parameter A of Eq. ( 72) becomes larger, which indicates that the diabatic character of the potential curves becomes larger with X12 becoming larger.In other words, the coherence becomes smaller as the adiabaticity between the two potential curves is larger.

FIG. 9.
Spatiotemporal distribution of the reduced direct nuclear flux F nu (Q, t).The panels are arrayed in the same manner as in Fig. 3.The direct nuclear flux can arise from any type of dynamics in the Born-Huang representation.Figure 8 exhibits the direct nuclear flux J nu (r, Q, t) in (z, Q) space (with x = y = 0) at time t = 20.0 fs, 40.0 fs, and 60.0 fs.A quick comparison between Fig. 8 and Fig. 4 shows that J nu (r, Q, t) is quite similar in shape of space-time distribution to that of the density Γ(r, Q, t) of Eq. (77).

The Journal of
For the single-configuration adiabatic wavefunctions Φ ad 1 ({r}; Q) is usually a real-valued function, and since there is no nonadiabatic interaction applied, both j el (r, Q, t) and J el nu (r, Q, t) are identically zero, and only J nu (r, Q, t) can have a finite value (see Table I).These results are mathematically correct, but are rather counter intuitive in that the electrons can flow synchronously with the nuclear motion.Therefore, several papers have been published to define the electronic flux from an adiabatic electronic state within the Born-Oppenheimer approximation, which is referred to as adiabatic electronic flux: Diestler and his coworkers have developed their coupled-channel theory to detect the electronic-state variation synchronously transported by the nuclear motion in general coordinate systems including those of the laboratory-frame. 36,41Pohl and Tremblay 42 and Engel and his co-workers 43 have given extensive discussions from symmetry-breaking and a flux-flux FIG. 10.Space time distribution of the indirect nuclear flux J el nu (r, Q, t) in (z, Q) space (with x = y = 0) for t = 20.0 fs, 40.0 fs, and 60.0 fs.The panels are arrayed in the same manner as in Fig. 3.The panels in the right column are all void representing J el nu (r, Q, t) = 0 for the Born-Oppenheimer dynamics with X 12 = 0. reflection principle, respectively.Okuyama and Takatsuka proposed time-shift flux, in which time information is externally introduced in terms of molecular deformation as in the ab initio MD. 44 Despite the intuition that the adiabatic electronic wavefunction should explicitly bear the information about electronic flow that will lead chemical reactions, Eqs. ( 35) and ( 36) are rigorously correct as far as the adiabatic wavefunction of Eq. ( 80) with the real-valued Φ ad 1 ({r}; Q) is concerned.The above seen similarity between J nu (r, Q, t) and Γ(r, Q, t) along with Eq. ( 37) suggests that the information about the adiabatic electron flux should be embedded in J nu (r, Q, t).Indeed, it is shown in the Appendix that the adiabatic electronic flux in the Ab Initio Molecular Dynamics (ABMD) can be taken by the "nuclear derivative" over the "electronic wavefunction." For the multiple-configuration adiabatic wavefunctions, the discussion about the adiabatic (electron) flux is more complicated, since j el (r, Q, t) is not zero as shown above.Besides, the J nu (r, Q, t) cannot be simply approximated by the product of a single nuclear velocity and a single electron density [see Eq. ( 37)], since the interference among the different configurations is expected as seen in Eq. (25).Therefore the notion of "adiabatic flux" of the multiple-configuration wavefunctions, if any, should be reconsidered.
The reduced nuclear flux F nu (Q, t) of Eq. ( 26) is presented in Fig. 9.All the panels there are not particularly interesting except that the multiple-configuration adiabatic wavefunctions are seen to give the bifurcation pattern in F nu (Q, t).As a matter of fact, this is not the wavepacket bifurcation.Two belts representing F nu (Q, t) come close to each other in the avoided crossing region and leave from each other asymptotically.On the other hand, F nu (Q, t) arising from the nonadiabatic wavefunction merge into a single piece after passing across the nonadiabatic region.The adiabatic wavefunction of A = 0.5 (the first row, right column) is composed of a major component and a very small part, and thereby the resultant F nu (Q, t) looks like the nonadiabatic counterpart.It turns out that F nu (Q, t) is far less sensitive than j el (r, Q, t) of Fig. 5 in detecting the effect of nonadiabaticity.

FIG. 11.
Space time distribution of the reduced indirect nuclear flux F el nu (Q, t).The panels are arrayed in the same manner as in Fig. 3.A close look at the left panels of A = 1.0 and 2.0 shows a spacetime oscillation of the reduced flux.Again all the panels in the right column are completely blank.
C. Indirect nuclear flux J el nu (r, Q, t), F el nu (Q, t) Among the flux components so far considered, only J el nu (r, Q, t) and F el nu (Q, t) can distinguish the nonadiabatic dynamics from the dynamics of multiple-configuration adiabatic wavefunction in a very clear cut manner.As shown in Table I and confirmed in Figs. 10 and 11, both J el nu (r, Q, t) and F el nu (Q, t) are identically zero only for the multipleconfiguration adiabatic wavefunctions.Therefore J el nu (r, Q, t) and F el nu (Q, t) can distinguish the nonadiabatic dynamics and adiabatic dynamics qualitatively.Up to this stage of the flux calculations, the Bon-Oppenheimer (adiabatic) dynamics and the nonadiabatic dynamics could not be differentiated.This is one of the highlights of the present study.
Unfortunately, however, the magnitudes of J el nu (r, Q, t) and F el nu (Q, t) arising from the nonadiabatic dynamics are about 10-100 times smaller than those of J nu (r, Q, t) and F nu (Q, t).Nevertheless, a close look at the panels of A = 1.0 and 2.0 with X 12 0 identifies a space-time oscillation of the reduced flux.

D. Correlation function
As the final part of the present numerical survey, we study a correlation between the electronic and nuclear fluxes with an expectation that this quantity should give a deeper insight about the correlated motion between electrons and nuclei.A significant difference between the nonadiabatic dynamics and multiple-configuration adiabatic dynamics will be identified.Incidentally, if the adiabatic dynamics is represented in a single-configuration wavefunction, the corresponding correlation function is identically zero. 27The flux-flux correlation function between electronic and nuclear dynamics is defined as where the operators fz (z ) and FQ (Q ) are defined as follows:

FIG. 12. The correlation function C(t)
between the electronic and nuclear fluxes.The panels are arrayed in the same manner as in Fig. 3. and where the arrow → (←−) on top of the differential operators indicates that they should be applied to the ket (bra) vector, and µ is the relevant reduced mass.Since the flux-flux correlation functions appear in the study of the rate constant for chemical reactions, [53][54][55] the difference of the adiabatic and nonadiabatic dynamics may send a signal to the studies of the relevant physical quantities.
Figure 12 exhibits the correlation functions in the same arrangement of the fluxes as before.The left column shows C(t) for the nonadiabatic dynamics with growing order of A = 0.5, 1.0, and 2.0 from the top to the row.The correlation functions in the right column are for the multipleconfiguration adiabatic dynamics, which shows a significant difference from the nonadiabatic counterparts.
The qualitative interpretation on the global feature of C(t) for the nonadiabatic dynamics is rather simple.First we notice that the correlation function beyond about t = 60 fs is zero since the electronic fluxes both of the nonadiabatic and adiabatic dynamics become diminished as studied above.On the way of the nonadiabatic dynamics, the electronic flux tends to be induced in the vicinity of the avoided crossing.After the two components χ ad 1 (Q, t)Φ ad 1 ({r}; Q) and χ ad 2 (Q, t)Φ ad 2 ({r}; Q) merge into a single piece, the electronic flux becomes zero, which accordingly nullify C(t) as well.This gives the flat feature at C(t) already about t = 50.The large single peak around t = 40 is conceived to be given as follows: Before the nonadiabatic transition, χ ad 1 (Q, t)Φ ad 1 ({r}; Q) represents a continuous dynamics of F − Li + → F + Li (F is at z = 0, while Li is at z that becomes larger).Therefore electrons jump in the positive direction.On the other hand, χ ad 2 (Q, t)Φ ad 2 ({r}; Q) is supposed to represent the process FLi → F − + Li + , which should give the electron flow in the negative direction.However, due to the nonadiabatic interaction, the component χ ad 2 (Q, t)Φ ad 2 ({r}; Q) becomes quickly smaller, and eventually χ ad 2 (Q, t)Φ ad 2 ({r}; Q) becomes zero in the asymptotic region, and the flux arising from F − Li + → F + Li given by χ ad 1 (Q, t)Φ ad 1 ({r}; Q) becomes predominant in the entire correlation function, which is roughly a product of the two positive quantities.
In the multiple-configuration adiabatic dynamics, the role of χ ad 1 (Q, t)Φ ad 1 ({r}; Q); F − Li + → F + Li is not different much from the nonadiabatic dynamics, but that of χ ad 2 (Q, t)Φ ad 2 ({r}; Q); FLi → F − + Li + , gives an electron flux in the negative direction and makes a dramatic difference.Unfortunately, the effect of the interference between them cannot be simply broken down to allow for a simpler explanation.The numerical calculation claims that the factor χ ad 2 (Q, t)Φ ad 2 ({r}; Q); FLi → F − + Li + , which is absent from the nonadiabatic dynamics should result in the negativity of the correlation function.

V. CONCLUSIONS
We have presented an electronic and nuclear flux analysis on the nonadiabatic dynamics in comparison with the corresponding adiabatic dynamics of a wavefunction that is of multiple configuration in the Born-Huang representation.We have shown both theoretically and numerically that only the indirect nuclear flux J el nu (r, R, t) and its reduced fluxes give a clear-cut distinction between the two dynamics, since J el nu (r, R, t) is simply zero for any kind of adiabatic dynamics.The role of the nonadiabatic dynamics in chemical dynamics needs no additional comment.On the other hand, the experimental and theoretical studies on the multipleconfiguration adiabatic dynamics are rare but will become more and more increasing, since the excitation dynamics triggered by ultrashort pulse lasers of the order of tens of attoseconds may produce the multiple configuration states due to a large energy-time uncertainty. 35As explicitly shown in the present paper, the quantum mechanical interference can induce all the flux components except for J el nu (r, R, t) even without nonadiabatic interactions.This is rather counter-intuitive since each electronic-nuclear configuration χ I (R, t)Φ I ( q ; R) involved in the Born-Huang expansion in Eq. ( 5) can be propagated in time in a completely independent manner.Therefore it is not appropriate to claim that the adiabatic wavefunctions running on the different potential curves should not result in any physical consequence.On the contrary, a theoretical or experimental manipulation that can respond to such induced fluxes can make the otherwise independent dynamics couple each other, resulting in a modulation of the Born-Huang total wavefunction in the asymptotic regions.Yet, one notices a theoretical difficulty: If the induced electron current in the adiabatic multiple-configuration dynamics happens to be large enough to emit light, the system should lose the energy to the corresponding amount despite the absence of any external perturbation.Nonetheless, the present multipleconfiguration adiabatic dynamics should behave as though nothing happened by definition.This apparent inconsistency should be a natural consequence of the naive application of the primitive Schrödinger equation without taking account of the radiation field.
Nonetheless, the present analysis thus comprehended should give a theoretical foundation to study the chemical and/or physical phenomena that emerge from the induced electronic and nuclear fluxes due to the interference among the different configurations.

ARTICLE
scitation.org/journal/jcpsemiclassical Ehrenfest theory (SET), not only to show that the flux are much dependent on the total wavefunction chosen but also because ABMD and SET are widely applied in the quantum chemical analysis of reactions.

Adiabatic flux in ABMD
The so-called Ab Initio MD (ABMD) electronic wavefunction Φ(r; R(t)) (A1) contains the time direction through the nuclear MD path R(t).
Here again, the electron flux j el (r, R, t) is identically zero, since Φ itself is real-valued.However, the amount of the electron flow can be estimated as follows: First consider the electron density ρ(r, R(t)) = Φ(r; R(t)) 2 .
(A2) The slight shift of the nuclear configuration may be associated with the time variation of the electronic state through bond formation, cleavage, and so on, and the resultant change of the electron density is Theoretically (not numerically), the right-hand side of this expression is similar to that of Eq. ( 28) in the indirect nuclear flux J el nu defined in Eq. ( 27).One can define the "electronic flux" j ABMD (r, R(t)) in such a manner that ∂ρ ∂t = −∇ r • j ABMD (r, R(t)), (A4) which may be regarded as a class of "adiabatic electronic flux." The ABMD electronic wavefunction is exactly the same as the adiabatic electronic wavefunction at each nuclear configuration R = R(t).However no nuclear wavefunction has to be calculated and is tractable for the simple ABMD scheme.
Incidentally, ∇ R Φ(r; R(t)) can be estimated by the finite difference method or a scheme via the Hellmann-Feynman theorem along with the resolution of identity.where the coefficients C I (t) are generally complex and R(t) represents a single path averaged over the electronic states.We here consider only the SET, and the extension to the more general case as in the path-branching representation such as Eq.(A6) 4,5,[13][14][15][16] is rather straightforward.The main difference of the "electronic" wavefunction of ϕ(r, t; R(t)) of Eq. (A7) from the ABMD electronic wavefunction Φ(r; R(t)) of Eq. (A1) is obviously that ϕ(r, t; R(t)) is an explicit function of time.Therefore there are two time coordinates in ϕ(r, t; R(t)).Hence j el (r, t, R(t)) arises from the electronic state change, which may be symbolically written as In the earlier studies by the Takatsuka group, 5,11,12 J SET nu (r, t, R(t)) has not been considered explicitly.
FIG. 3. Space-time distribution of the nuclear wave packet density.Left (right): Pairs of the panels show non-adiabatic (adiabatic) dynamics.The top, middle, and bottom rows show the case of A = 0.5, 1.0, and 2.0, where A is the amplitude of the artificial nonadiabatic term defined in Eq. (72).

2
) and adiabatic dynamics ( X12 = 0) on thus created sets of the potential curves, respectively, are shown.The figures following from Fig. 4 up to Fig. 11 are presented in this manner.

TABLE I .
Components of the electronic and nuclear fluxes and their properties.