Relativistic theory of electron-nucleus- radiation coupled dynamics in molecules: Wavepacket approach

We propose a general theoretical scheme of relativistic electron-nucleus coupled dynamics of molecules in radiation fields, which is derived from quantum electrodynamical formalism. Aiming at applications to field-induced dynamics in ultrastrong laser pulses to the magnitude of 10 W/cm or even larger, we derive a nonperturbative formulation of relativistic dynamics using the Tamm-Dancoff expansion scheme, which results in, within the lowest order expansion, a time-dependent Schrödinger equation with the Coulombic and retarded transversal photon-exchange interactions. We also discuss a wavepacket type nuclear dynamics adapted for such dynamics. © 2019 Author(s). All article content, except where otherwise noted, is licensed under a Creative Commons Attribution (CC BY) license (http://creativecommons.org/licenses/by/4.0/). https://doi.org/10.1063/1.5109272., s


I. INTRODUCTION
Electrons accelerated close to the speed of light exhibit relativistic effects. [1][2][3][4] In molecular science, relativistic electrons appear in orbitals of heavy atoms where it has been established that not only the inner-shell electrons [5][6][7] but also valence electrons are subject to relativistic effects such as the orbital contraction/expansion effects. 8,9 Another important type of relativistic electrons in molecular science is those accelerated in ultraintense optical fields. 10,11 The former aspect has been discussed separately in our recent work 12 within the path integral formalism. We herein develop a relativistic "wavepacket" theory for coupled electronic and nuclear dynamics of molecules, which is applicable to dynamics in external fields of nonperturbative intensity regime.
Rapid progress in experimental techniques 10,11 has realized ultraintense laser pulses of petawatt intensity and femtosecond width, 13 and a number of advanced light sources are being developed worldwide. 14 The strongest laser field ever reported 15 exceeds 10 22 W/cm 2 . Nonrelativistic quantum mechanics breaks down in such ultrastrong fields as the vector field amplitude in the Hamiltonian compares mec 2 /|e|, with me and |e| being the electronic mass and charge magnitude, respectively. Introducing nondimensional relativistic scale for a given vector field A by a 0 ≡ |e| |A|/(mec 2 ), a typical infrared (IR) laser pulse of wavelength λ = 800 nm ( ̵ hω = 1.55 eV) of intensity 2.14 × 10 18 W/cm 2 corresponds to a 0 = 1. Indeed, with pulse intensity of order ∼10 16 W/cm 2 or more, experimental studies have demonstrated relativistic behaviors of electrons including relativistic trajectories of photoelectrons, 16 radiation reactions, 17,18 and multiphoton Thomson scatterings. 19 Molecular dissociations in laser fields of relativistic intensity have also been reported in Refs. 20 and 21. In Ref. 20, CH 4 molecules irradiated by IR laser fields of intensity exceeding 10 16 W/cm 2 were observed to undergo multiple tunnel ionizations, followed by rapid fragmentation. 20 Such dynamics is particularly challenging in that strong electron-radiation coupling does not allow perturbative approach which is the most standard tool in the quantum electrodynamics (QED). Another interest is in the interplay of such relativistic electron-radiation coupled dynamics and the nuclear dynamics. Taking into account the typical nucleus-to-electron mass ratio of order 10 3 or more, we can assume that the nuclear dynamics is nonrelativistic even in ultrastrong fields with a 0 ≲ 100. Relativistic electronic dynamics, however, should induce unconventional nuclear dynamics such as Coulomb explosion of heavily ionized molecules. 20 Relativistic theory of stationary electronic many-body states within the Born-Oppenheimer (BO) approximation has been intensively developed in the so-called quantum chemistry. [22][23][24] Numerical studies have revealed that relativistic electrons, which exist in the inner-core levels of heavy element atoms, play several important roles in the chemical properties of molecules. 22 Static properties of those molecules containing heavy elements, including their bonding properties [25][26][27] and spectroscopic observables, 28 have been successfully reproduced by numerical ab initio calculations. Formulations of numerical calculations include the original four-component formulation and more tractable two-component formulations based on the Foldy-Wouthuysen transformation [29][30][31] or even more sophisticated theories. [32][33][34][35] Such two-component reduction techniques are sometimes indispensable for the study of large systems, whereas its extension to dynamical problems is not straightforward. We therefore have the present formulation remain in the original fourcomponent theory unless specified explicitly. A remarkable merit in the quantum-chemical approach is its large variety of methodologies of constructing (correlated) many-body wavefunctions; techniques for constructing four-component wavefunctions include the variational theory 36 and the exact solution technique. 37 Multiconfiguration theories to include correlation effects to the self-consistent field orbitals have also been developed. 38,39 Along with the methodological sophistication, there are a large number of "tools" available for achieving efficient calculations; studies have established a number of efficient basis sets 40 and calculation techniques. 41,42 A number of computational packages for general purposes have also been developed. 43,44 However, in the BO framework, the involved nuclei are made fixed in space and no time variable is considered in the electronic wavefunctions, let alone in the radiation fields. Despite its great success in the calculations of static molecular properties, the main assumptions presupposed in the BO approximation (and thereby the conventional quantum chemistry) are not valid when the kinematic coupling between electrons and nuclei are strong. Indeed, the physical and chemical significance of the so-called nonadiabatic reactions has long been recognized since the early stage of quantum mechanics [45][46][47] to date. [48][49][50] Overcoming such breakdown has become even more important in the regime of laser chemistry, 51 which makes it possible to observe 52 and control (see Ref. 53) the kinematic coupling between electrons and nuclei. In particular, it has been found that the notion of ultrafast nonadiabatic electron wavepacket dynamics is particularly crucial. [54][55][56][57][58][59] Therefore, one of the main objectives in this paper is to extend such non-BO wavepacket dynamics to the relativistic regime. Another shortcoming of the conventional type of relativistic quantum chemistry is that it often shows lack of theoretical rigor; predominantly many of those studies have been made based on the Dirac-Coulomb or the Dirac-Breit theory, which are the low-energy static reductions of the exact QED. Many features in the exact QED treatment, including dynamical coupling to the radiation-field and/or energy-dependence of the electron-electron interactions, are missing in those types of quantum chemical approaches, whereas it can easily be expected that those missing components are essential in high-energy dynamics of our interest. We therefore start from the exact QED formulation to construct a framework of relativistic nonadiabatic electron-nucleus coupled dynamics. The most standard QED formulation would then take a perturbative approach to include electron-radiation coupling effects. We indeed have taken such an approach in our recent study. 12 It is, however, obvious that such an approach is not applicable to dynamics with nonperturbative couplings to the (external) radiation fields, which is our main target in this paper. In order to develop a nonperturbative approach, we consider the Tamm-Dancoff (TD) 60,61 expansion of the electron-radiation coupled states. We also note that there have been the development of advanced formulations of quantum chemistry that resolves such discrepancy from the exact QED treatment. For example, the effective quantum electrodynamics (eQED) approach 23,24 constructs an effective Hamiltonian that includes static QED effects such as the vacuum polarization effect and the one-loop electronic self-energy arising from virtual interactions with the radiation field. It is of theoretical interest to make comparison between these static QED theories and our formulation.
The original idea of the Tamm-Dancoff approximation (TDA) was independently developed by Tamm 60 and Dancoff 61 earlier than the establishment of the covariant perturbation theory. [62][63][64][65][66] The original idea of Dancoff 61 was to expand the Heisenberg state vector, or the time-independent energy eigenstate, by a well-defined set of basis vectors in various "sectors" to derive a set of Schrödingerlike equations. A problem in this original approach, pointed out by Dyson,67 is that the number of interacting particles adds up in a less-controlled manner than what is in the covariant perturbation theory, making the theory nonrenormalizable. Several modifications to overcome such a problem were proposed including Dyson's formulation. [67][68][69] Among them, Cini 68 and Visscher 69 developed an alternative formulation closer in idea to the covariant perturbation theory, which shows that divergences arising from lower order expansions (Cini 68 considered up to what corresponds to the twoloop expansion in diagrammatic language) can be removed using the renormalization technique. Even with those modifications, TDA has several fundamental problems 70 including apparent violation of covariance by its naive truncation scheme. Wilson and co-workers formulated TDA in the light-cone frame 70 to overcome these inconveniences. Their light-front TDA approach 71 has then been applied to quantum color dynamics (QCD). To the best of our knowledge, however, it is not clear that such a drastic change in the frame is applicable or suitable for molecular dynamics. We therefore do not use the light-cone frame but stick to more familiar types of space-time frame (the laboratory or the center-of-mass frame) which may sacrifice the formal rigor. Nevertheless, it appears that the Tamm-Dancoff approach provides a method to formulate electronradiation coupled dynamics of relativistic energy scale based on "equal-time" amplitudes, which describes a many-body system on a spacelike surface with a single timelike variable. We consider that, in our future study, the equal-time representation would work favorably in importing many-body techniques developed in quantum chemistry in which the same type of space-time frame is (implicitly) assumed.
This paper is organized as follows. In Sec. II, we first prepare the fundamental quantities relevant in the present relativistic dynamics. In Sec. III, we formulate a theory based on the Tamm-Dancoff 60,61 type expansion of the electron-radiation coupled states. This paper concludes in Sec. IV.

II. DEFINITIONS OF FUNDAMENTAL QUANTITIES
We first define fundamental quantities that appear in our discussion. Although this paper is self-contained, a more comprehensive discussion on some of those quantities is shown in our recent work. 12 Throughout this paper, we use the sign convention (1, −1, −1, −1), and the metric tensor ημν is a diagonal tensor with η 00 = 1, η 11 = η 22 = η 33 = −1. Symbol qe represents the electronic charge, which takes a negative value; qe = −|e|. We use the Gauss unit for electromagnetic field, and the fine-structure constant is e 2 / ̵ hc ≈ 1/137. We use the Hermitian Dirac matrices, α and β. Four-dimensional notation α μ stands for α μ ≐ (1, α), i.e., α 0 is the four-dimensional unit matrix. Other notations, which follow the standard convention in QED, 2 are described explicitly in the main text and also summarized at the end of the main text.
We consider an electron-nucleus-radiation coupled system in the Coulomb gauge, whose full Hamiltonian reads where α and β are the Dirac matrices, Pa and Ma represent the threedimensional momentum and mass of the ath nucleus, whereas ψ and A tr are the field operators of the electronic and transversal radiation fields, respectively. The canonical conjugate of A tr is represented by Π tr =Ȧ tr /(4πc 2 ), which is negative of the transversal part of the electric field. Symbols ρmat and Jmat represent the matter field charge density and current, respectively, as and with qe representing the electronic charge (qe = −|e|), Ra being the ath nuclear (spatial) coordinate, and Za|e|f a (R) being the charge distribution of the ath nucleus, which could simply be a delta function Za|e|δ 3 (R − Ra) or one of the existing model functions. 72 Possible nontrivial nuclear charge distributions are taken into account only for calculation of the electronic wavefunction, whereas we approximate them by the pointlike model for description of the field-nucleus and the nucleus-nucleus Coulombic interactions. Such an approximation should be appropriate in the energy range where the typical wavelength of the radiation field is much larger than the length scale of the nuclear charge distribution. We then introduce the nuclear potential field and a set of electronic "mean-fields" ρ(x, t) and W HF (x, t) and introduce a Fermionic eigenvalue equation with 0 ℓ representing the ℓth orbital energy. The eigenfunctions of Eq. (4) are hereafter referred to as molecular orbitals (MOs). We assume that we are interested in the state with a given configuration Iocc of Ne positive energy MOs (if we are interested in the ground state, for example, the lowest Ne are occupied and occupied orbitals are Iocc = {1, 2, . . . , Ne}), The most natural choice of the mean-field potential W HF must be the Hartree-Fock potential, which is However, we find it more convenient to use its local approximation and hereafter consider the localized Hartree-Fock field, W loc HF , which is to be obtained by one of existing techniques. [73][74][75][76] For later convenience, we also introduce MO creation/annihilation operators asĉ ℓ,m,n,... = {ĉ r,s,t,... positive energy stateŝ b a,b,c,... negative energy states, (7) and the whole set of positive (negative) MO indices are denoted by I + (I − ).
Our Hamiltonian now reads H tot = Hnuc +H el mf +H el int +H rad with expansions, we have to rewrite each term in a symmetric form 77 when we calculate its vacuum expectation value. The current operator J μ (x, t) = cψ † (x, t)α μ ψ(x, t) is therefore to be interpreted as when we calculate its vacuum expectation value.

III. THEORY OF ELECTRON-NUCLEUS-RADIATION COUPLED DYNAMICS
Having formulated fundamental quantities, we now develop the Tamm-Dancoff expansion of the electron-nucleus-radiation coupled dynamics.
A. The lowest order expansion in electron-radiation coupled dynamics We first discuss electron-radiation coupled dynamics for a given nuclear configuration. Couplings to dynamical nuclei are restored in Subsection III B. Following the original idea of Tamm 60 and Dancoff, 61 we expand the electron-radiation coupled state vector as a superposition of zero, one, two, . . . photon states. We formally introduce notation (Ne, Ne, N ph ), where Ne, Ne, and N ph represent the number of electrons, positrons, and photons. Each triplet represents a "sector." We assume that the physical state of our interest is dominated by a small number of "active" sectors around (Ne, 0, 0) and truncate other sectors. We first consider the electron-radiation dynamics for a given nuclear coordinate. We also assume that the mean-field equation [Eq. (4)] is already solved and the relevant MOs are known. We then set the zeroth order Hamiltonian as and the residual interaction for the correction is given as Eq. (8c). We first consider the smallest nontrivial set of "sectors," which consists of (Ne, 0, 0), (Ne, 0, 1), and (Ne + 1, 1, 1). For simplicity of notation, we also use more intuitive labels for those sectors: 000, 001, and 111, respectively.
A brief note is made here on the (Ne + 1, 1, 1) or 111 sector; it consists of states with one electron-positron pair plus a photon, wherein we are not necessarily assuming existence of "real" electronpositron pair, which would require energy of order ∼1 MeV. In the lower energy dynamics, those states are regarded as "virtual" states, which are, nevertheless, not simply negligible for its possible roles in four-component dynamics (see Ref. 78, for an illustrative example) and for achieving mathematically consistent description of dynamics.

Coordinate representation
Following Ref. 68, we define an "amplitude" [hereafter referred to as the Tamm-Dancoff (TD) amplitude] of the electronic state vector in (Ne, 0, 0) sector as where the electronic annihilation operatorψ, the vacuum state |0⟩, and the time-dependent electronic state vector Ψ int t are all in the interaction representation. The time evolution equation reads We do not have to take account of all terms in H el int of Eq. (8c) to obtain the correct equation of F 000 t . Since there is a strict oneto-one correspondence in the transversal photon exchange and the Coulombic interaction (plus its associated mean-field subtraction terms), we can concentrate on the transversal interaction part , are restored when we evaluate the obtained diagrams by re-interpreting the photon lines as a summation of the Coulombic and transversal interactions. Except for those representing the self-interactions (i.e., the photon lines emitted and absorbed by the same electron), we also attach the mean-field subtraction terms to those photon lines.
For later convenience, we introduce the mean-field electronic propagator S F by with T being the time-ordering operator. We now expand the last term in Eq. (12) using H el int tr given as Eq. (13). Details of the expansion procedure is given in Appendix A, but in short, we move all creation operators appearing in H el int tr to left by application of the equal-time (anti-)commutation rule. We then obtain where Λ + (Λ − ) is the projection operator to positive (negative) energy states. The index kλ represents the photon wavevector k and the polarization index λ, whose associated frequency being ω k = c|k|.
For notational convenience, we combine the polarization vector e k,λ μ and all numerical factors associated with the photon annihilation operator A tr and further contract with the electronic Dirac spinors α μ to define that the corresponding spinor matrix operates on the ℓth component of the rank-Ne spinor F as J μ k symbolizes the Fourier transformation of the vacuum expectation value of the current where J μ (x, t) represents the vacuum expectation value of the current operator with t − 0 in the last side indicating an infinitesimally small time before t. While the corresponding quantity vanishes in the free-particle system, 64 it generally has a nonzero value in our model because of the existence of the mean-field potential. It can be evaluated either by direct numerical calculation using a large set of MOs or by perturbation expansion of the mean-field propagator S F by the free-particle propagator We also note that the term "vacuum" comes from the vacuum state in the interaction picture |0⟩, which is actually affected by the mean-field Hamiltonian. All related ideas in our following discussion, including the vacuum polarization and vacuum energy shift, reflect the mean-field potential as their physical origin.
New amplitudes belonging to the (Ne, 0, 1) and (Ne + 1, 1, 1) sectors appearing in Eq. (15) are defined, respectively, as and The derivation of Eq. (15) is given in detail in Appendix A. The time evolution equation for F 001 reads and that for F 111 becomes + (terms that do not belong to the (Ne, 0, 0) or (Ne, 1, 1) sector) + (terms arising from H el int ′ ), (23) where the expression F| i ℓ →j p indicates the spinor F whose ℓth index i ℓ is replaced by jp.
Concentrating on the terms proportional to the (Ne, 0, 0) sector amplitudes only, a formal solution of Eq. (22) can be symbolically written as whereas Eq. (23) is solved as Substituting these expressions into Eq. (15) and restoring the Coulombic and mean-field subtraction terms, we obtain whereĥ mf j and W loc HFj represent the mean-field Hamiltonian in the operator representation and the mean-field potential acting on the jth electronic coordinate in F 000 t , whereas VC and V SI represent the Coulombic (inter-electronic) interaction and the Coulombic self-interaction, 79 respectively. Symbols C > and C < represent the zeroth order propagators defined as with ψ> (ψ<) representing the positive (negative) energy part of the electronic field operator. Using those propagators, the for- appearing in Eqs. (24) and (25) is expanded as with ϕt representing an arbitrary physical quantity and {x} ({x}) representing electronic (positronic) coordinates. Equation (26) describes the time evolution of the electronic system that interacts through the Coulombic interaction and the transversal photon exchange, and yet Eq. (26) also includes apparently divergent self-interaction terms whose renormalization is better formulated in the MO-representation as we discuss below.

Molecular orbital representation
We reformulate the above procedure in the MO-representation. The MO-representation of Eq. (16) is defined as We also use a collective notation I, which represents a set of MO indices of the Ne (positive energy) electronic orbitals as I ≡ {I 1 , I 2 , . . . , IN e }, and I r j is given by replacing the jth component Ij by the rth MO as I r j ≡ {I 1 , I 2 , . . . , I j−1 , r, I j+1 , . . . , IN e }. We then define amplitudes in 000, 001, and 111 sectors by Substitution of solutions of Eqs. (31b) and (31c) into Eq. (31a) yields the following equation of f 000 t : The Journal of Chemical Physics ARTICLE scitation.org/journal/jcp with ⟨Ij I ℓ ∥r s⟩ representing the instantaneous Coulomb-exchange interaction among four MOs Ij, I ℓ , r, and s (in the physicists' notation).
By the symbolΔ vp f 000 t (I), we separated out contributions that purely arise from vacuum polarization and does not have direct effects on the dynamics of molecular electrons. It readŝ which will be discussed later; below we concentrate on the terms of dynamical relevance. Equation (32) contains a number of self-interaction terms, the Coulombic as well as the transversal photon exchange interactions. These are apparently divergent, and we need appropriate renormalization to obtain a meaningful result. We first combine two transversal self-interaction terms in Eq. (45) to define a quantity Σ tr 1LI j ,s (t − t ′ ), and its Fourier transformation where η represents an infinitesimal positive constant and the meanfield propagator S F in the MO representation expands as We then combine the Coulombic self-interaction contribution as well, to obtain Self-energy renormalization using the Coulomb gauge is formulated in Refs. 80 and 81, and we may well use these results. For simplicity of discussion, however, we here assume that the final result should be invariant with respect to the choice of the gauge, and we replace the photon propagator by the Feynman gauge propagator, The Journal of Chemical Physics

ARTICLE scitation.org/journal/jcp
where Σ 1L (x, y; ω) is the one-loop electronic self-energy constructed from the Hartree-Fock propagator. Apparent divergence in Σ 1L (x, y; ω) can be renormalized using the standard procedure by subtracting the counter terms, which leaves a finite contribution Σ fin. 1L (x, y; ω). Formal expression of the finite part of the single-loop self-energy can be written as where S F 0 is the free-electron propagator and δmc 2 and Z 2 are the mass and field-renormalization constants. 2 Detailed discussion of the self-energy, including its numerical calculation procedure, is found in a number of excellent reviews including Refs. 7 and 82.
We next consider terms including J μ , which arises from the vacuum polarization. The two terms with (J ν −k ε kν ) or its conjugate but not the both describe scattering of molecular electrons by the current J ν , as are schematically shown in Figs. 1(d) and 1(e). In order to get better insight, we introduce the following potential field arising from J μ : (40) Indeed, the two scattering terms in Eq. (32) essentially reduce to scattering by V vp μ. Apparent complication arises from formal distinction of the time ordering of the source term (J μ ) and the scattering event ( ε kλ ), as well as existence of the upper limit of the integration range in Eq. (32). Assuming approximate stationarity of J μ , we therefore should be able to replace those two terms by V vp μ. We further combine the Coulombic contribution to Eq. (40) as  32) and (33).
Applying the same argument as Eq. (38), we can replace the square bracket in Eq. (41) by the corresponding expression in the Feynman gauge, −4πq 2 e ημν/(k 2 + i0). We can then combine it into the selfenergy as We finally discuss the vacuum polarization term we separated out as Eq. (33). It consists of two types of interaction, the Coulombic and the transversal photon exchange, and each has direct and exchangelike terms [see Figs. 1(f) and 1(g)]. A tricky point is that some of those terms have apparent dependence on the zeroth electronic energy, and it includes convolution with f 000 t (I). We can, however, separate it by taking account of time dependence of the amplitude, ft(I) ∝ e −i ∑ k 0 I k and extending the time integration to ∞ to rewrite it asΔ vp f 000 where ΔE vp is the vacuum energy shift,  (44) is a formally correct expression for the one-loop vacuum energy shift in QED. Furthermore, as such terms only contribute to a phase factor of the system's wavefunction, we hereafter drop this term for simplicity.
Thus, we can now rewrite Eq. (32) as where the self-energy term is defined as and the delta functions in the curly bracket indicate that either τ ′ 1 or τ ′ 2 is to be fixed to t (with t − 0 being an infinitesimal time before t) and τ ′ < represents the smaller one of τ ′ 1 and τ ′ 2 . We can also rewrite the real-space expression Eq. (26) as The real-time expression of the self-energy operator Eq. (39) represents "fluctuation" (i.e., emission and absorption of virtual photons) induced by the radiation field. Assuming such effects can essentially be included in the form of energy correction, we can approximately rewrite Eq. (45) into a simpler form, where Σ ℓ,m is a time-independent matrix, defined by the Fourier transformation of the self-energy operator as which is seen to converge to an appropriate expression 83 in the static limit. We thus have removed the divergent terms in Eq. (47). Although our results Eq. (48) or its real-space representation Eq. (47) are aimed at relativistic dynamics, they have interesting applications in the nonrelativistic limit as is discussed in Appendix B.
In this paper, we have limited our discussion to the lowest order expansion. The higher order expansions [i.e., formulation taking account of higher orders of H el int appearing in the RHS of Eq. (12)] can be in principle obtained in the same framework although it is not straightforward. We show a sketch of the higher order extensions in Appendix C, and the full expressions will be worked out in future study.

B. Coupling to nuclear dynamics
We next consider inclusion of nuclear degrees of freedom. We start with the formal second-quantized representation. We disregard internal structures of the atomic nuclei and assume that all nuclei in the system are distinguishable nonrelativistic particles. The nuclear Hamiltonian Eq. (8a) can formally be rewritten as ]χa(R) withχa (χ † a ) representing a formal annihilation (creation) operator of the ath nucleus. We then extend our definition of the TD amplitudes as A tricky aspect in the present procedure is that our definition of the TD amplitudes requires the zeroth order Hamiltonian in order to specify the interaction representation. We therefore first assume a fictitious potential function U(r, t) which mimics the true electron-nucleus Coulombic interaction potential Ûnuc(r) defined as an operator,Û Assuming such a space-time function U(r, t), the "mean-field" Hamiltonian is redefined as and we can set the zeroth order Hamiltonian H 0 along with the correction term H int as with and where the latter one H ′ el−nuc represents the difference between the fictitious field and the true electron-nucleus potential. Having thus defined the zeroth order Hamiltonian, the extended TD amplitude given by Eq. (51) is now well-defined. We can then derive an equation of motion for F 000 We find that commutation of the nuclear annihilation operators with Hnuc formally yields additional couplings to the (Ne, 0, 1) sector although such radiation-nucleus coupling terms in general should be much smaller in amplitude compared to its electronic counterpart because of large masses of the atomic nuclei and/or slowness of the nuclear motion. Below we only keep the first order coupling terms, which are proportional to ̵ h∇a/i, whereas we neglect the second order terms, which are proportional to (A tr ) 2 . Terms arising from H el int can be expanded in the same manner as we did in Subsection III A. The remaining terms are expressed in terms of the derivatives of the (Ne, 0, 0) sector amplitude.
We therefore find that the amplitude F 000 t ({x}; R) behaves as a "first-quantized nuclear wavefunction" on which Hnuc operates as a differential operator. We therefore rewrite the definition of the TD amplitude as in which the electronic mean-field Hamiltonian is constructed with the use of the true electron-nucleus interaction, H mf (R) = H mf [Ûnuc] (with the operator Ûnuc rewritten in the first-quantized operator form), so that we no longer need the fictitious field U. The above definition can be naturally extended to other sectors as and and so on. The equation of motion for F 000 Since we have switched from the fictitious U to the true Ûnuc, the electronic wavefunctions depend explicitly on the nuclear coordinate R.
We therefore have included the derivative coupling operator Xa, which acts on electronic variables in the manner (Xa) ij = ⟨φi| ∂/∂Ra|φj⟩ and possible mass-polarization terms, K. The latter one K represents a set of coordinate-frame dependent differential operators arising from the coupling to the nuclear dynamics; in the center of the mass of the nuclei (CMN) frame, for example, there should be mass-polarization terms as we discuss in Appendix D. For simplicity, however, we here use the Cartesian coordinates in the laboratory frame, where K vanishes.
We next consider the equation of motion of F 001 t ({x }; k λ; R), which reads + (terms arising from H el int ′ ) + (terms that do not belong to the (Ne, 0, 0) or (Ne, 0, 1) sector), and that for F 111 becomes ).
In addition to those belonging to other sectors than (Ne, 0, 0), (Ne, 0, 1), or (Ne + 1, 1, 1), we neglected all terms arising from nucleusradiation couplings as these describe the emission of radiation by the nuclear charge current, which should be negligible due to slowness of the nuclear motion (|Ṙa| ≪ c). Switching to the MO representation and solving Eqs. (62) and (63), we obtain a closed equation of (Ne, 0, 0) amplitude, where we neglected terms that contain nucleus-radiation coupling to the second order. Equation (64) is to be compared to the previous result [Eq. (45)], where the nuclear dynamics is absent. In addition to the nuclear kinematic energy term, which is in the first term in Eq. (64), we have additional terms that arise from electron-nucleus interaction by the transversal photon exchange. Although we admit those terms should be small in general by the factor |Ṙ|/c, it gives the lowest order relativistic correction to the electron-nucleus interaction due to relativistic electronic motion. We also note that nuclear derivative operators in Eq. (64) accompany the derivative couplings i ̵ hX. Since the amplitude is explicitly dependent on the nuclear coordinates, appearance of the derivative coupling terms is a manifestation of the "gauge invariance," or covariance of the expression with electronic basis set transformations. Equation (64) represents the main conclusion of this subsection. We also derive the corresponding equation in an external field in Appendix E, where we show that the external field should formally be put into H 0 so that one can treat it in a nonperturbative manner. Our formulation thus defines a Schrödinger-like equation of motion in the equal-time representation (i.e., all coordinates sharing the same timelike coordinate) of a many-body system; hence, the space-time frame is the same as that of the conventional nonrelativistic many-electron theory. Also taking account of the fact we are assuming the nuclear motion as nonrelativistic, we expect that there should be no fundamental difficulty in integration of our equation of motion, Eq. (64) using one of existing calculation techniques. Yet, we propose in Appendix F a wavepacket implementation of the nuclear dynamics in an attempt to seek for a best suited numerical algorithm for relativistic electron-nucleus coupled dynamics.
For clarity of later discussion, we summarize the assumptions we made in this formulation. First, we assumed (i) nonrelativistic nuclear motion. The assumption is essential not only for describing the nuclear kinematic term in the Schrödinger-like form in Eq. (64) but also for treating nucleus-radiation coupling as a weak perturbation. For laser-induced dynamics, this assumption puts a formal upper limit of the laser intensity around ∼10 22 W/cm 2 , assuming λ = 800 nm and nuclear mass of order 10 3 me.
Second, we truncated the higher order terms in the Tamm-Dancoff expansion beyond the lowest order, which in turn implicitly assumes (ii) irrelevance of the higher order electron-radiation coupled effects. Such assumption should be valid in models where smallness of electron-radiation coupling suppresses the higher order effects, but it would break down in unconventional types of dynamics such as those discussed in Refs. 17-19, where higher order QED effects play vital roles. This assumption (ii), in practice, may therefore impose much lower limit than (i) does to the external field intensity acceptable in our formulation. Relevance of higher order effects is, however, not clear just from an order estimate. A desirable strategy is therefore to start from a parameter range ∼10 16 W/cm 2 to ∼10 18 W/cm 2 and critically examine the validity of the lowest order treatment.
We, however, note that the absence of higher order effects in assumption (ii) refers to the "dynamical" electron-radiation coupling such as those arising from physical processes represented by more complex Feynman diagrams. It should not be confused with coupling with "external" radiation fields such as ultrastrong laser fields, which can be included in We have also mentioned about the gauge invariance, which we used to simplify the expression of the photon propagator, Coulomb plus transversal photon to the Feynman gauge propagator [see, for example, Eqs. (37) and (38) and related discussions]. We do not consider it as an essential assumption since numerical studies indeed show gauge independence of the final results, 80,81 let alone difference in the contribution from each individual diagram, and there is also an established technique to directly evaluate the Coulomb gauge propagator. 80,81 Another feature of our present approach is that our working equation [Eq. (64)] is a direct extension of the conventional time-dependent Schrödinger equation (TDSE) approach in the sense it reduces to TDSE in the limit of vanishing dynamical electronradiation coupling. We therefore expect that we can apply it to those models where the conventional TDSE approach provide at least qualitatively reasonable results in order to obtain an insight on how the relativistic effects affect the dynamics.

IV. SUMMARY AND DISCUSSIONS
In this report, we have developed a nonperturbative approach for relativistic electron-nucleus coupled dynamics. We consider that this approach has potential applications to the electron dynamics induced by ultraintense IR laser fields.

ARTICLE scitation.org/journal/jcp
Our approach essentially reduces to a set of effective Schrödinger-like equations of equal-time amplitudes of the system; hence, it can be combined with nonperturbative calculation techniques including grid-type time-dependent wavepacket approaches.
Possible limitations yet encountered in this approach are rooted to the problems inherent to the Tamm-Dancoff formulation. 67,70 Nonetheless, the present theory has resolved difficulties faced in application of QED to a marked extent.
The formulation presented here only uses the lowest order expansion, which reproduces two of the essential features in the quantum electrodynamics, retarded and/or energy-dependent interactions among electrons and self-energies arising from the selfinteraction of electrons. However, we need higher order expansions to include other types of QED effects including vertex corrections and multiple scatterings. Although we included a sketch of such a formulation in Appendix C, full formulation is to be discussed in another paper.
As we discussed in our closely related publication, 12 the performance in practical calculations should depend strongly on the quality of electronic wavefunctions to be used, which are supposed to take into account an electron-electron correlation. Compared to the path-integral formulation, 12 the present formulation, which derives a Schrödinger-like equation for a many-body wavefunction in the equal-time representation, has better similarity and affinity to the conventional Hamiltonian formalism of nonrelativistic theory, which might work favorably. We therefore can rest on many computational tools 40 and techniques 41,42 well developed in quantum chemistry.
Another theoretical interest is to make a close comparison between our obtained equation [Eq. (45)] in the static limit and existing static QED theory such as the eQED approach. 23, 24 We can indeed find rough correspondence to their theory by performing the Fourier transformation to those integral kernels, as we did in Eq. (49). Nevertheless, partly because it is not very straightforward to define a strict static limit of Eq. (45), a detailed comparison is beyond our current scope and left for future study.
• Symbol ≷ attached as a subscript or superscript of an electronic annihilation or creation operator specifies the energy sector; for example, ψ> and ψ< represent a positive-energy electron annihilation operator and a negative-energy electron (positron) creation operator, respectively. • Symbol T represents the time-ordering operator unless specified otherwise. • Symbol R represents the nuclear coordinate, which is a 3Nn dimensional vector. The mass and charge of the ath nucleus are denoted by Ma and Za|e|, respectively. • Symbol a kλ represents the annihilation operator of a (transversal) photon with the wavevector k and polarization λ. The associated normalized polarization vector is denoted by e μ k λ .

ACKNOWLEDGMENTS
The work has been supported by JSPS KAKENHI, Grant No. JP15H05752.

APPENDIX A: DETAILS OF THE TAMM-DANCOFF EXPANSION
We here describe some details of the Tamm-Dancoff expansion. As a simple example, we show the derivation of Eq. (15), the time evolution equation of an (Ne, 0, 0) sector amplitude.
We first recall the definition of F 000 t in Eq. (11). Since each operator is in the interaction representation,ψ( its time-derivative yields commutator with H 0 , as shown in the first term in the RHS of Eq. (12). It then contributes to the mean-field single-particle Dirac Hamiltonian term H 0 F 000 t ({x}) in the RHS of Eq. (15). On the other hand, the state vector appearing in the definition of F 000 t [Eq. (11)] is in the interaction representation. Hence, its time derivative yields i ̵ h ∂ ∂t |Ψ int t ⟩ = H int (t)|Ψ int t ⟩, which appears in the second term in the RHS of Eq. (12). All Fermionic and Bosonic creation operators appearing in H el int (t) are then moved leftward until they come next to the bra-vector ⟨0|, where they eliminate the vacuum state. During this operation, the Fermionic (Bosonic) anticommutation (commutation) rules between equal-time operators apply.
Following the main text, we consider the transversal part, H el int tr , in the interaction Hamiltonian. We see that H el int tr is decomposed as It then follows that, in the equation of F 000 , all the "photon creation" operators trivially move to the left end and vanish, whereas all the photon annihilation operators survive. We next see that, in the electronic operators, ψ † > (ξ, t) and ψ<(ξ, t) contain the creation operators. We see that the positron creation operator ψ<(ξ, t) moves to the left end and vanishes, whereas the electron creation operator ψ † > (ξ, t) anticommutes withψ(xj, t) to yield Λ+δ 3 (xj − ξ); hence, there appears ε kλ (xj)ψ>(xj, t) in the position of the jth operator, yielding a vector in the 001 sector.

APPENDIX B: SPIN-ORBIT AND SPIN-SPIN COUPLING IN THE NONRELATIVISTIC LIMIT
Although it is not our major task, the low-energy effective theory with nonrelativistic two-component approximation gives by-products which are certainly useful in many applications.
As is suggested in the textbook of Bethe and Salpeter, 84 an effective Hamiltonian for interacting particles can be derived from the Tamm-Dancoff expansion. First we note that, in contrast to the remark made in Ref. 84 that (standard) TDA only reproduces the Coulombic term, our Tamm-Dancoff-like expansion does indeed recover those equations that include transversal interaction terms.
Following Ref. 84, we study an Ne-particle system in an external vector field potential V ext which is a slight modification of our real-space equation (47), in which we use, following Ref. 84 Sec. 39, the free-particle Hamiltonian as the zeroth order one; The corresponding positive-energy propagator is denoted by K 0 > , which reads with upσe ipx/ ̵ h being the positive energy solution of h 0 with spinprojection σ and εp being its energy, εp ≡ √ c 2 p 2 + m 2 c 4 . Assuming the existence of its stationary solutions, we perform Fourier transformation of Eq. (B1) to obtain where we introduced frequencies ω can be seen as an implicit eigenvalue problem, where one solves an eigenvalue ω by ω-dependent matrix in the RHS. We also see that by neglecting the frequency dependent term (ω (2) i (k)/c) 2 , the transversal photon contribution reduces to  (2) i (k)/c is not negligible in comparison with a typical value of k, which is the inverse of the typical length scale of the electronic wavefunction. Assuming inner core electron, such an inverse length scale is of order ∼ (aB/Z) −1 , with aB being the Bohr radius and Z being the (effective) atomic number. This implies the frequency condition ̵ hω (2) i (k) ∼ Z ̵ hc/(aB), which is usually much smaller than the "relativistic" energy scale ̵ hω (2) i (k) ∼ mc 2 .
As for Fω, we apply, following Ref. 84 where the indices of type aj take 0 or 1, whereas those of type Cj take + or −. Symbol θ C 1 ,C 2 ,...,C Ne a 1 ,a 2 ,...,a Ne represents a 2Ne-component spinor, which is labeled by an index C 1 , C 2 , . . . , CN e and has 2 N e -dimensional spinor index a 1 , a 2 , . . . , aN e . We also recall that, in our convention, ij represents the jth (four-component) spinor index of the amplitude F 000 ω . Four-by-two spinor matrices ΞC are given by with 1 2 representing the two-dimensional unit matrix. Positive (negative) energy projection Λ + p (Λ − p ) is defined as with which we have where the explicit expression of α μ CC ′ (p, p ′ ) are found in Ref. 84 Eq. (16.12) [α 0 CC ′ (p, p ′ ) in our notation corresponds to ICC(p, p ′ ) in Ref. 84]. Here, we concentrate on the lowest order expansion with respect to the |p|/mec term of α μ +,+ (p, p ′ ), which are, according to Eq. (16.14) in Ref. 84, Substitution of expansion Eq. (B5) into Eq. (B3) and application of positive-energy projection operator ∏ j Λ + p j retaining only a single component C 1 , C 2 , . . . , CN e = + + ⋯ +, we obtain where we have neglected the self-energy term whose contribution is negligible in nonrelativistic limit and summation over spinor indices is made implicit, that is, spinor products are to be understood as ∑ a 1 ⋅ ⋅ ⋅ ∑ a Ne Ξ a 1 ,...,a Ne ++⋯+ θ ++⋯+ a 1 ,...,a Ne . Equation (B10) only differs by the retarded interaction from Ref. 84.
We next expand the retarded interactions. Substitution of Eqs. (B9a) and (B9b) into the last term in the RHS of Eq. (B10) yields the following three terms: which are then transformed into real-space representation as where q ≡ ω (2) i (k)/c, rij ≡ri − rj, and nij ≡ rij/|rij|, and the primed bracket [⋯] ′ indicates that the singularity of the expression in the square bracket at the origin has been removed. 84 We note that these expressions reduce to static expressions in the limit ω → 0 as The Journal of Chemical Physics ARTICLE scitation.org/journal/jcp A qualitative difference in our results [Eqs. (B12b) and (B12c)] in contrast to its well-known static counterpart [Eqs. (B13b) and (B13c)] is that the present expressions bear the longer-range terms (with respect to spatial separation rij) arising from higher order terms of qrij.

APPENDIX C: HIGHER ORDER EXPANSIONS IN THE ELECTRON-RADIATION COUPLED DYNAMICS
Here, we discuss higher order expansions in the electronradiation coupled dynamics to be appended to the lowest order expansion in Subsection III A. We find that a straightforward extension of our discussion made on the lowest order expansion in Subsection III A does not work to a larger number of sectors because of divergences arising from the summation over photon k vectors. Since such k vector summation arises from all loop structures in the diagrammatic representation, we first need to separate out these divergent summations. We restart from a formal perturbation expansion of the state vector |Ψ int t ⟩, which is assumed to correspond to a reference state (eigenstate of H 0 ) |Φ 0 ⟩ in the limit of far past t → −∞ as where η is an infinitesimal positive constant, which is to be later taken the limit η → 0, and the associated normalization factor is defined as Nη ≡ √ ⟨Φ 0 |Uη(t, −∞)|Φ 0 ⟩. We then find that the amplitude is written as whose time derivative leads to an integrodifferential equation, We herein concentrate on the formal theory for simplicity and assume that each integral can be extended to t → −∞ and that the state in the far past t → −∞ is a radiation-free reference state of Ne electrons with no photons or antiparticles. We also extend our definition of the amplitudes. Only in this section, we allow appearance of multiple electronic timelike coordinates in the amplitude, where we require that all timelike variables should not be smaller than t; x 0 j /c ≥ t. Those amplitudes can then be calculated by spatial convolution with function C > , with nlm being either of 000, 001, or 111, and F nlm t ({x}, ⋯) in the RHS being the corresponding amplitude with a single timelike variable t.
As we did in the main text, we first replace H el int by H el int tr and expand the interaction terms. We then find that, in Eq. (C3), the first order correction vanishes, whereas the second order survives. After restoring Coulombic contributions, the second order contribution reads The remaining part in Eq. (C3) is expanded as and where the four-dimensional integral with upper limit is a shorthand notation defined as whereas the symbol v ′ C represents the Coulombic interaction minus the mean-field interaction, We see that Eqs. (C10b)-(C10d) are divergent. It is, however, clear that summation over all the different time-ordering with the same topology as Fig. 2(b) makes an expression that is obtained from Eq. (C10b) by reordering the time. We can then apply the standard renormalization procedure for the electronic self-energy Σ to remove formal divergences and obtain a finite result. Application of the same procedure to diagrams Figs. 2(c) and 2(d) [Eqs. (C10c) and (C10d)] should yield the polarization and vertex correction functions, respectively. We then need expressions for the amplitudes of the (Ne, 0, 1) and (Ne + 1, 1, 1) sectors. Applying the simplest truncation scheme to Eqs. (22) and (23), to truncate sectors other than (Ne, 0, 0), (Ne, 0, 1), and (Ne + 1, 1, 1), we obtain and Substitution back of the formal solutions of Eqs. (C11) and (C12) then leads to an effective equation of the (Ne, 0, 0) sector such that Here, an additional divergence arises from the summation over kλ if the corresponding photon line makes a loop. But it can be treated in the same manner as Eq. (48). Although Eq. (C13) only shows an outline of the procedure, it is clear that summation over all the diagrams arising from expansion of Eq. (C3) corresponds to the standard Feynman diagrams of the two loop expansion in the standard perturbation theory, in which all formal divergences can be removed.
We therefore conceive that there should be no divergence problem although we have yet to calculate all the possible diagrams in MO representation carefully examining whether this expansion works.

APPENDIX D: MASS-POLARIZATION TERMS
Here, we discuss the mass-polarization terms K we encounter in the derivation of Eq. (61). In the main text, we adopted the lab frame and simply set it zero. On the other hand, if we use the coordinate representation relative to the center of the mass of the nuclei (CMN), Nonzero mass polarization emerges. Hence, we need to consider a general transformation from the lab frame to the CMN frame to derive those terms. Below, the lab frame electronic coordinates are denoted by {rj} (j = 1, 2, . . ., Ne) and the nuclear coordinates in the mass-weighted representation are denoted by {Qa} (a = 1, 2, . . ., Nn), whereas in the CMN frame, electronic coordinates are denoted by {rj} and the mass-weighted nuclear coordinates are denoted by {Ξ λ }, among which the CMN coordinate is denoted by We then consider the following generalized transformation with an orthogonal matrix R: The associated differential operators then transform as Substitution of these results to the electronic Dirac equation makes no apparent change, whereas that to the kinematic part of the nuclear Schrödinger equation yields with {A, B} + ≡ AB + BA representing the anticommutator, and Ra in the RHS are to be understood as a function of new coordinate vectors {Ξ λ } in the sense Ra = ∑ λ (R −1 ) aλ Ξ λ . We can then find the mass-polarization term which couples among all the electronic momenta and also charge-polarization terms that add to the electron-radiation coupling term although both terms accompany small factor me/(∑ M b ). We also see that the CMN mode formally couple with the electronic momentum operators reflecting the Galilean covariance although we can simply assume that the CMN is fixed at the origin in many of applications.
In dynamical models of molecules, however, the derivative couplings, which approximately scales as √ me/M λ for the λth (λ ≠ ω) nuclear displacement modes around the equilibrium configuration, can be even larger in magnitude than these mass-polarizations.

APPENDIX E: EXTERNAL FIELDS
We here consider chemical dynamics (CD) in an external vector field ext (r, t). We keep discussions in the main text where we set that A tr in Eqs. (8a)-(8d) as "internal" or "dynamical" field while adding new terms In principle, one can treat these terms in an analogous manner as what we did in Subsection III A. We then need to expand the electronic part in Eq. (E1), H el int ext ≡ −qe ∫ d 3 xψ † (x)αψ(x)⋅ ext (x, t), which introduces additional couplings between (Ne, 0, 0) and (Ne + 1, 1, 0) sectors, (Ne, 0, 1) and (Ne + 1, 1, 1) sectors, etc. It then appears, in the lowest order, we need to solve a set of coupled equations of the (Ne, 0, 0), (Ne, 1, 0), and (Ne, 1, 1) sector amplitudes, which is not a trivial task, since there appears potentially divergent summation over photon wavevectors.
A better insight is obtained by first solving the singleparticle propagator in the external field and then by solving the radiation coupling. In such a scheme, one can directly reach the final conclusion just by replacing propagators in Eq. (64) as The Journal of Chemical Physics ARTICLE scitation.org/journal/jcp state). In such a formulation, it is, in general, not appropriate to use a single wavepacket, but one needs to introduce bifurcation at a point where the character of the electronic state changes drastically. We therefore develop a general (multiple) wavepacket expansion of the nuclear dynamics. We then augment our theory with an appropriate "branching" algorithm in order to take account of bifurcations.
In this section, we use the mass-weighted representation of the nuclear coordinates, and hence, R here represents mass-weighted coordinates. We start with a formal expansion of the electronnucleus coupled state, where |ΦA : R⟩ represents a time-independent electronic state at nuclear coordinate R and χA(R, t) represents its associated nuclear wavefunction. Recalling Eq. (64) in Sec. III, a formal equation of motion of χ can be rewritten in the following form: where Tnuc is the (nonrelativistic) kinetic operator for nucleus, where ∇ = ∂/∂R is the 3Nn-dimensional (mass-weighted) gradient, Qnuc represents a diagonal matrix which consists of reduced charges (the nuclear charges divided by the square-root of nuclear masses), and ext represents the external field supervector ext ≐ (A ext (R 1 ) t , A ext (R 2 ) t , . . .). Since we neglect the radiation from the nuclear charge current, the resultant vector fields are regarded as an external field. Matrices XAB and YAB represent the derivative couplings of the first and second order; XAB ≡ ⟨A|∇|B⟩ and YAB ≡ ⟨A|∇ 2 |B⟩. We then expand χA(R) as with normalized Gaussian g, defined as where Qσ, ̵ hKσ Γσ, and Θσ represent the central coordinate, momentum, inverse width, and phase of the σth Gaussian. The inverse width matrix Γσ is symmetric. For convenience, we also introduce a shorthand notation for nuclear coordinate integral in the bracket form, where O represents an arbitrary function or operator in the nuclear coordinate representation. The overlap matrix is then defined as In the RHS of Eq. (F1), coefficients C σ A are related to the total population of the Ath adiabatic state nA by where the unitarity condition requires ∑AnA = 1. We then apply the time-dependent variational equation, with H tot representing the total (electron-nucleus) Hamiltonian and bracket ⟨⋯⟩ representing the integration of electronic as well as nuclear coordinates, to derive equations of motion for C σ A (t), and the parameters in the Gaussian functions, Qσ, Kσ, and Γσ.

A
We first decompose the total Hamiltonian H tot by the electronic basis as where T is the kinetic operator acting only on the nuclear wavefunction, whereasKAB is the remaining part of the Hamiltonian, where "hat" on K is meant to emphasize that it is an operator working on the nuclear wavefunction. We also use notation H eff AB to denote the "electronic part," We now write down the variational stationary condition for C σ where we introduced the notation and f (Γτ,Γτ; τ) represents with the operator Tr representing the trace over coordinatelike indices. We then obtain the time evolution equation for C σ A as where∂t represents time derivative whose operation is limited to inside the bracket and bra-vectors with tilde and ⟨σ| are defined as ⟨σ| ≡ ∑ τ (S −1 ) στ ⟨τ|.

Variation of the Gaussian parameters
We first note that variation of the Gaussian parameters ξ, which is one of the four parameters K ℓ σ , Q ℓ σ , Γ ℓm σ , and Θσ, can be written in the following form: where the superscript (ξ) represents the variational parameter and the corresponding O (ξ) x is in general a function of nuclear coordinates, We then write the variational stationary condition for those parameters, where the matrices with Gaussian superscripts σ, τ, . . . represent and we also define related quantities for the unit matrix Although we have to find such a solution that eliminates all variation in principle, it would lead to complex expressions for numerical calculations. Here, we rather derive a simpler approximate solution and later confirm that major contributions to variation do vanish.
We first set ̵ hKσ =Q σ (F20) so that all terms containingQ σ vanish at all order. We then concentrate on large contributions, arising from σ = τ, or diagonal terms in Eq. (F17), and eliminate variation of type O = 1 and O = y σ . It can then be shown that all the other residual terms arise from the higher order correlation. Diagonal contribution in the RHS of Eq. (F17) then reads where we have taken account of the fact X σσ = 0, which follows from the antisymmetry of XAB. Requiring Eq. (F21) for O = 1, O = y ℓ σ and O = y ℓ σ y m σ to vanish, we obtain In deriving these equations, we took account of the fact that any odd momenta of Gaussian ⟨σ|y ℓ 1 σ ⋅ ⋅ ⋅ y ℓ 2n+1 σ |σ⟩ = 0 and that ⟨σ|y j σ y ℓ σ |σ⟩ =  We now fix the time derivatives of the first two set of parameters,Θσ and ̵ hKσ by Eqs. (F22a) and (F22b), respectively, together withQ σ by Eq. (F20). On the other hand, we formally leaveΓ in the following expressions as we later employ the frozen Gaussian approximation to fix Γσ at a reasonable value, such as that satisfies Eq. (F25c) which shows that the variation indeed vanishes for Oσ = 1 and the residual part consists only of the higher correlation for the other Oσ.

A practical calculation scheme
Having established formal equations of motion, we now consider a practical calculation scheme with branching. We first apply the frozen-Gaussian approximation, to avoid complexity and possible unphysical deformation of Gaussian wavepackets. We set Γσ of the initial wavepacket to some reasonable value, in such a manner that Γ 2 equals the Hessian of the electronic Hamiltonian at the initial geometry and setΓ = 0 in subsequent time evolution. We thus start from a single wavepacket, but the final state should be in general a superposition of multiple wavepackets that asymptote to mutually independent adiabatic states. Since our equations of motion do not automatically increase the number of Gaussian wavepackets, we need an additional algorithm that introduces new wavepackets. Following Ref. 55, we call it a "branching" algorithm.

ARTICLE scitation.org/journal/jcp
There are several clues to judge path-branching that directly follows from our formulation: (i) the variation residual of Eq. (F17), which does not exactly vanish except for O = 1, or (ii) K (2) στ AA and Δ ofd στ in Eq. (F30), which should be small by physical requirements. One can therefore introduce a new wavepacket such that it reduces those inconvenient quantities. Such algorithm should be most desirable for fully consistent realization of our Gaussian wavepacket formulation; however, we here seek for a simpler algorithm that directly realizes branchings.
An alternative approach is to judge from the nuclear kinematics. In case we know the asymptotic adiabatic states, we can (a) introduce branching so that each wavepacket should consist of a single adiabatic state in the final asymptotic region. Even if we do not know those states a priori, we can also use more intuitive quantity, "force matrix," introduced in Refs. 54 and 55. In the present formulation, it appears as the first (nonradiative) term in the RHS of Eq. (F26), where we find a state-dependent force, or a variant of the force matrix, One can then (b) use the eigenvectors of the force matrix to decompose the wavepacket into most rapidly departing parts. 55 In either scheme, (a) or (b), we can construct branching algorithm following the idea proposed in Ref. 55: (1) Set a relative population threshold n rel thr and relative adiabatic coupling thresholds ξ rel thr .
(2) Start checking wavepackets which have a mixture of distinct adiabatic states more than n rel thr when the strength of the nonadiabatic coupling decreases and passes through the thresholds ξ rel thr (relative to the peak) downward, and divide the wavepacket by projecting onto either (a) the adiabatic state vectors or (b) the eigenstates of the force matrix.
At this point, our discussion on the branching algorithm cannot be further extended since branching does not directly follow from the variational principle [Eq. (F8)]. The best algorithm is therefore only be judged through numerical applications by comparing accuracy, convergence, computational costs, etc., which we leave to our future study.