Fully general time-dependent multiconfiguration self-consistent-field method for the electron-nuclear Dynamics

We present the fully general time-dependent multiconfiguration self-consistent-field method to describe the dynamics of a system consisting of arbitrary different kinds and numbers of interacting fermions and bosons. The total wave function is expressed as a superposition of different configurations constructed from time-dependent spin-orbitals prepared for each particle kind. We derive equations of motion followed by configuration-interaction (CI) coefficients and spin-orbitals for general, not restricted to full-CI, configuration spaces. The present method provides a flexible framework for the first-principles theoretical study of, e.g., correlated multielectron and multinucleus quantum dynamics in general molecules induced by intense laser fields and attosecond light pulses.


I. INTRODUCTION
We are now witnessing rapid progress in ultrashort intense light sources in different spectral ranges such as terahertz radiation, optical-parametric-chirpedpulse-amplification mid-infrared lasers, high-harmonic extreme-ultraviolet (XUV) pulses, and XUV/x-ray freeelectron lasers.These technological advances have triggered various research activities, including attosecond science [1][2][3], with a goal to directly measure and, ultimately, control electron and nuclear motion in atoms and molecules.
Ab initio simulations of the electronic and nuclear dynamics in atoms and molecules remain a challenge.The multiconfiguration time-dependent Hartree-Fock (MCT-DHF) method [4,5] has been developed for the investigation of multielectron dynamics in strong and/or ultrashort laser fields [6].In this approach, the timedependent total electronic wave function Ψ(t) is expressed as a superposition of different Slater determinants Φ I I I (t), Ψ(t) = where C I I I (t) is the configuration-interaction (CI) coefficients.Both {C I I I (t)} and the spin-orbitals constituting {Φ I I I (t)} are allowed to vary in time.In the community of high-field phenomena and attosecond physics, the term MCTDHF is conventionally used for the full-CI case, in which the sum in Eq. (1) runs over all the possible ways to distribute the electrons among a given number of spin-orbitals.On the other hand, also under active development are variants without the restriction * anzaki@atto.t.u-tokyo.ac.jp to the full-CI expansion, generically referred to as the time-dependent multiconfiguration self-consistent-field (TD-MCSCF) methods hereafter.The representative examples include the time-dependent complete-activespace self-consistent-field [7,8], the time-dependent restricted-active-space self-consistent-field [9], and the time-dependent occupation-restricted multiple activespace (TD-ORMAS) [10] methods.These allow a compact and computationally less demanding description of the multielectron dynamics, without sacrificing accuracy.
In particular, the TD-ORMAS can treat arbitrary CI expansions of the form Eq. (1) in principle.Among successful approaches for nuclear dynamics is the multiconfiguration time-dependent Hartree (MCTDH) method [11].Developed for systems consisting of distinguishable particles, this method expresses the time-dependent total nuclear wave function as a superposition similar to Eq. ( 1) but that of Hartree products.The other way around, the MCTDHF can be viewed as an extension of the MCTDH to fermions.By hybridizing the MCTDHF for electrons and the MCTDH for nuclei, one can construct a multiconfiguration electron-nuclear dynamics (MCEND) method [12] to describe the non-Born-Oppenheimer coupled dynamics.Nuclei forming molecules are, however, indistinguishable particles, either fermions or bosons.
In this Paper, further stepping forward in this direction, we present a fully general TD-MCSCF method for a system comprising of arbitrary different kinds and numbers of interacting fermions and bosons.Treating all the constituent particles on an equal footing, we expand the total wave function in terms of configurations of the whole system [see Eq. ( 5) below], rather than considering configurations of each particle kind separately as in Ref. [13].Thus, based on the time-dependent variational principle, we derive the equations of motion (EOM) of CI coefficients and spin-orbitals for general configuration spaces, not restricted to full-CI.This paper is organized as follows.Section II introduces our TD-MCSCF ansatz for many-particle systems composed of different kinds of fermions and bosons, and also defines the target Hamiltonian considered in this work.In Sec.III, we derive the general equations of motion, based on the time-dependent variational principle.Explicit working equations for a molecule interacting with an external laser field are shown in Sec.IV.Concluding remarks are given in Sec.V.

II. DEFINITION OF THE PROBLEM A. TD-MCSCF ansatz
We consider a quantum mechanical many-body system with K kinds of fermions or bosons.The subsystem of kind α consists of N α identical particles.Thus, there are N = K α=1 N α particles in whole.For notational brevity, we call such a system an N N N -particle system, where the array of integers N N N = (N 1 N 2 • • • N K ) carries information of both particle kinds and number of particles in each kind.
Let us define, for each kind of particles, the complete orthonormal set of spin-orbitals {χ (α) µα (t) : µ α ∈ Ω α }, which spans the one-particle Hilbert space Ω α , and are time-dependent in general.Then the N N N -particle Hilbert space is spanned by I I Iα (t) is a determinant (or parmanent) of α-kind fermions (or bosons), consisting of N α spin-orbitals chosen from {χ (α) µα }.We call Φ I I I (t) the I I I's configuration, where I I I = I I I 1 I I I 2 • • • I I I K is considered, at the moment, to collectively label the chosen spin-orbitals.The objective of this paper is to formulate the TD-MCSCF theory of the N N N -particle system within the ansatz of total wavefunction analogous to that for electronic system, Eq. ( 1), but using the configurations of Eq. ( 2).
For rigorous and compact presentation of theory, we resort to the second quantization formulation by introducing creation and annihilation operators {ĉ µα }.These operators obey the (anti-)commutation relations of bosons (fermions), for bosons, where [â, b] = âb − bâ, and for fermions, where {â, b} = âb + bâ.Within the TD-MCSCF ansatz, the complete set of spin-orbitals {χ We call the subspace of Ω α spanned by occupied spinorbitals the occupied spin-orbital space Ω occ α , and that spanned by virtual spin-orbitals the virtual spin-orbital space Ω vir α , where Ω α = Ω occ α ⊕ Ω vir α .The total state Ψ(t) is expressed as a superposition of configurations Φ I I I (t) of Eq. ( 2), but constructed from occupied spin-orbitals only.Thus we write where C I I I (t) is the CI coefficient, and |I I I(t) is the occupation number representation of the configuration Φ I I I , Now ), and µ α , ν α , κ α , τ α , ... for general (Ω α ) spin-orbitals of kind α.The indices p α , q α will be used for numbering the coordinates.
It should be noted that we do not restrict the expansion Eq. ( 5) to the full-CI one.It should also be noticed that occupied configurations are specified in terms of the whole system rather than in terms of each particle kind separately as [13], Our approach allows a highly flexible choice of CI space, e.g., including up to double excitation [10] regardless of particle kind, thereby enabling proper account of correlation between different kinds of particles while suppressing computational cost.

B. Target Hamiltonian
In this article, we consider the Hamiltonian of an N N Nparticle system composed of up to M -body terms, The Hamiltonian is explicitly time-dependent in general, but the time argument t is dropped in this section for simplicity.Here, the m-body Hamiltonian is assumed to be given explicitly in terms of the coordinates (and momenta, see below) in a general sense characterizing m particles (or degrees of freedom), and symmetric under exchange of coordinates among particles of the same kind.One-particle Hamiltonian, e.g., is written as where the non-local form allows to describe the momentum dependence of the Hamiltonian, and two-body interaction is generally given by The reasons why we here consider the (non-local) higherthan-two body terms, which will not actually be used in Sec.IV, are (1) that such form is used in multiconfiguration Hartree (MCH) method for distinguishable particles, and (2) their possible appearance upon coordinate transformations, or in the effort of removing translational and rotational degrees of freedom [14,15].The Hamiltonian is equivalently expressed in the second quantization formalism as where the net m-body Hamiltonian is further classified into those contributions Ĥm m m , hereafter called m m mbody Hamiltonian, involving m α particles of the kind α, where indexes the set of spin-orbitals to represent m α particles in the Hamiltonian.Êµ µ µ ν ν ν is the m m m-particle replacement operator Êµ and (H m m m ) µ µ µ ν ν ν is given by where ) is the set of m α coordinates of particle α, and For the later discussion, we define the m m m-body reduced density matrix (RDM) as One-and two-particle RDMs are also denoted as (ρ αβ ) with β = α.

III. EQUATIONS OF MOTION
In this section, we derive the EOMs for the CI coefficients and spin-orbitals by imposing the time-dependent variational principle [16][17][18] on our TD-MCSCF ansatz.We require the action integral to be stationary, δS = 0, with respect to the variation of the total wavefunction δΨ within our TD-MCSCF ansatz Eq. ( 5), subject to the boundary conditions δΨ(t 0 ) = δΨ(t 1 ) = 0. To this end, let us introduce anti-Hermitian matrices ∆ α and X α as, (Recall that indices µ α , ν α refer to both occupied and virtual spin-orbitals.)We also define, with which orthonormality-conserving spin-orbital variations and time derivatives can be written as Then, the variation and time derivative of total state are compactly given by [7,10,19], where Π = I I I |I I I I I I| denotes the projector onto the CI space, i.e., the subspace of N N N -electron Hilbert space spanned by the configurations included in Eq. ( 5).The action functional S should be made stationary with respect to all independent variations; {δC I I I , δC * I I I } for CI coefficiens and {(∆ α ) µα να } for spin-orbitals.First, the EOM for CI coefficients are obtained from δS/δC * I I I = 0, i ĊI I I = J J J Requiring δS/δC I I I = 0 derives the complex conjugate of Eq. ( 27).Next from δS/δ(∆ α ) µα να = 0, one obtains i where Π = 1 − Π. Equation ( 28) is to be solved for (X α ) µα να = χ να , thus determines the time derivative of spin-orbitals.We now take a closer look at Eq. (28) for the following two distinct cases: Case 1: (µ α , ν α ) = (i α , j α ).In this case we focus on the components of the spin-orbital variations within the subspace spanned by the occupied spin-orbitals.Since Π( Êα ) iα jα |I I I = 0 and I I I|( Êα ) iα jα Π = 0 in general, one needs to directly work with Eq. (28) within the occupied spinorbital space i In the full-CI case, where Π( Êα ) iα jα |Ψ = 0, Ψ|( Êα ) iα jα Π = 0, Eq. ( 29) reduces to an identity 0 = 0. Therefore, the corresponding (X α ) iα jα may be arbitrary anti-Hermitian matrix elements, of which the simplest choice is (X α ) iα jα = 0.
Case 2: (µ α , ν α ) = (i α , a α ).In this case we deal with the components of the spin-orbital variations outside the occupied spin-orbital space.Since Ψ|( Êα ) iα aα Π = Ψ|( Êα ) iα aα and ( Êα ) iα aα |Ψ = 0, Eq. ( 28) becomes, i However, the matrix element in the left-hand side of the above equation survives only when β = α and κ α = a α ∈ Ω vir α , namely Ψ|( Êα ) iα aα ( Êβ ) Thus Eq. ( 30) is simplified to Now the RHS of Eq. (31) given by the sum over m m m, where (H α ) jα iα is the effective one-particle operator, and (W m m m ) k k k [α] l l l [α] is given in the coordinate representation as where y y y and x ′ α z z z [α] are defined similarly.Finally, gathering the occupied and virtual components of the time derivative completes the derivation of EOM for spin-orbitals i| χ(α) where Pα = iα |χ iα χ iα | is the spin-orbital projection operator onto the occupied spin-orbital space, with which the virtual space aα |χ (α) aα χ (α) aα | = 1 − Pα is referenced as a whole, thus avoiding explicit use of virtual spinorbitals.(X α ) jα iα in the first term is to be obtained by solving Eq. (29), and, as discussed above, can be set zero in the full-CI case.Equation (27) for CI coefficients and Eq.(36) for spin-orbitals form fully general TD-MCSCF equations of motion, not restricted to full CI, for a system composed of any arbitrary kinds and numbers of fermions and bosons.

IV. MOLECULES INTERACTING WITH AN EXTERNAL LASER FIELD
In this Section we present the working equations for a molecule subject to an external laser field.Let the molecule consist of electrons and K n different kinds of nuclei treated quantum mechanically (the kind does not necessarily corresponds to the nuclear species, see discussion below), and N cl nuclei treated as a classical point charge.For clarity and notational simplicity, we assign the electrons to the first kind of particle (α = 1), and kinds α = 2, 3, • • • , K represent quantum nuclei with K = 1 + K n .The numbers of identical particles are, as before, denoted by {N α }.Then the number of electrons is N 1 , the number of quantum nuclei is N n = K α=2 , and the total number of atoms is N atom = N n + N cl .We use atomic units in this section.
The spin-independent molecular Hamiltonian in the coordinate representation is given by U αβ (|r r r pα − r r r q β |), where U αβ (r) = Z α Z β /r is the Coulomb interaction with Z α being the electric charge, and h α (r r r, r r r ′ , t) = δ(r r r − r r r is the one-particle Hamiltonian composed of the kinetic energy [the first term with m α being the mass (not to be confused with the number of particles)], Coulomb interaction with classical nuclei with the charges {Z A } located at {R R R A } (the second term), and the time-dependent laser-particle interaction V ext α , given, e.g., within the dipole approximation either in the length gauge (LG) or in the velocity gauge (VG), by LG (r r r, r r r ′ , t) = −δ(r r r − r r r ′ )Z α E E E(t) • r r r, V ext α,VG (r r r, r r r ′ , t) = δ(r r r − r r r ′ )i where E E E(t) is the laser electric field, and A A A(t) = − E E E(t)dt is the vector potential.The general formulation of Sec.III is readily applicable to the molecular Hamiltonian of Eq. (37).The CI EOM reads (8) with I I I α , I I I β , I I I γ , • • • being the configuration of particle kind α, β, γ, • • • , respectively, and C I I IαI I I β I I Iγ ••• (t) the CI coefficient.