Implementation of a time-dependent multiconfiguration self-consistent-field method for coupled electron-nuclear dynamics in diatomic molecules driven by intense laser pulses

We present an implementation of a time-dependent multiconfiguration self-consistent-field (TD-MCSCF) method [R. Anzaki et al., Phys. Chem. Chem. Phys. 19, 22008 (2017)] with the full configuration interaction expansion for coupled electron-nuclear dynamics in diatomic molecules subject to a strong laser field. In this method, the total wave function is expressed as a superposition of different configurations constructed from time-dependent electronic Slater determinants and time-dependent orthonormal nuclear basis functions. The primitive basis functions of nuclei and electrons are strictly independent of each other without invoking the Born-Oppenheimer approximation. Our implementation treats the electronic motion in its full dimensionality on curvilinear coordinates, while the nuclear wave function is propagated on a one-dimensional stretching coordinate with rotational nuclear motion neglected. We apply the present implementation to high-harmonic generation and dissociative ionization of a hydrogen molecule and discuss the role of electron-nuclear correlation.


I. INTRODUCTION
The rapid development of laser technology enables the generation of coherent light sources that cover a broad frequency range, from the mid-infrared, to the x-ray regime [1][2][3][4][5]. The emergence of these new light sources has opened up new opportunities for many scientific fields, such as ultrafast atomic and molecular physics [6][7][8][9], attosecond chemistry [10,11], and material science [12][13][14][15], with a goal to visualize, understand and ultimately control electron and nuclear dynamics in atoms, molecules, and solids. In particular, when molecules are subject to ultrashort intense laser pulses, the induced molecular dynamics is the simultaneous motion of atomic nuclei and electrons forming a molecular entity [16,17]. An accurate description of such processes requires treating both electronic and nuclear degrees of freedom on an equal footing. At present, ab initio theoretical description of coupled electron-nuclear dynamics in molecules remains a challenge.
For the investigation of multielectron dynamics in intense laser fields, the multiconfiguration time-dependent Hartree-Fock (MCTDHF) method has been developed [18][19][20][21]. In this approach, the time-dependent total electronic wave function Ψ(t) is given in the configurationinteraction (CI) expansion, Ψ(x 1 , · · · , x n , t) = I C I (t)Φ I (x 1 , · · · , x n , t), (1) * yangli@atto.t.u-tokyo.ac.jp † sato@atto.t.u-tokyo.ac.jp ‡ ishiken@n.t.u-tokyo.ac.jp with x = {r, σ} being the spatial-spin variable of the electron. The electronic configuration Φ I (x 1 , · · · , x n , t) is a Slater determinant constructed from spin orbital functions {ψ(r, t) × s(σ)}, where {ψ(r, t)} and {s(σ)} are electronic spatial orbitals and spin functions, respectively. Both CI coefficients {C I } and orbitals are simultaneously propagated in time, which provides an adequate representation of the temporal change of the wave function with a reduced number of determinants compared to fixed orbital treatment. In the MCTDHF method, the summation in Eq. (1) runs over all possible configurations for a given set of spin orbitals, which are conventionally referred to as full-CI expansion. In the full-CI case, the computational cost due to CI coefficients is proportional to the number of electron distributions among all spin orbitals, which increases factorially with the number of electrons, thus hindering its application beyond small molecular systems. To subjugate this difficulty, variants that go beyond the restriction to the full-CI expansion are now under active development, which we refer to as timedependent multiconfiguration self-consistent-field (TD-MCSCF) methods [22,23] in the general case. Representative examples include the time-dependent completeactive-space self-consistent-field (TD-CASSCF) method [24][25][26], the time-dependent occupation-restricted multiple active-space (TD-ORMAS) method [27,28], and the time-dependent restricted-active-space self-consistentfield (TD-RASSCF) method [29][30][31]. These methods are based on the truncated CI expansion within the chosen active orbital space, which are computationally more efficient and at the same time provide compact representation of multielectron dynamics in a given physical condition without sacrificing accuracy.
Efforts have been made to describe the coupled elec-tronic and nuclear dynamics of molecules within the framework of the TD-MCSCF method by multiple groups and from diverse perspectives. Nest developed the multiconfiguration electron-nuclear dynamics (MCEND) method [32,33], where the molecular wave function is represented as a sum over products of determinants for the electrons and Hartree products for the nuclei. By introducing a set of atomic orbitals on "ghost" centers along the internuclear axis of diatomic systems, this method is computationally effective to investigate both electronic and vibrational photo-excitations. However, the MCEND method fails to describe photoionization and dissociation processes due to the use of Gaussian bases functions as electronic basis and limited grid range of nuclear bases [34]. These two processes are keys in understanding the strong coupling of electronic and nuclear dynamics to strong laser fields. Kato and Yamanouchi developed an extended MCTDHF method [35] by using a basis of Slater determinants with time-dependent nuclear orbitals for the fermionic protons while keeping heavy nuclei fixed. This method has been applied to calculate the ground-state properties of a model one-dimensional hydrogen molecule [36], which shows good agreement with the direct solution of the time-dependent Schrödinger equation (TDSE). Very recently, this method has also been adopted for real-time simulation of hydrogen molecular ions (with one electron and one nuclear degree of freedom) irradiated by an intense laser pulse [37]. In another work, Haxton and McCurdy developed a promising way to include quantum treatment of nuclei into MCT-DHF for the special case of diatomic molecules. By introducing prolate spheroidal coordinates for the electronic wave function [38], the total wave function is expanded in terms of configurations of orbitals which depend on the internuclear distance. Thus, the trivial but strong correlation between the positions of the nuclei and electrons due to their Coulomb attraction is naturally included and an efficient procedure for evaluating overlap and Hamiltonian matrix elements is formulated. On one hand, geometry-dependent electron orbitals satisfy the key physical condition for the inclusion of strong electronnuclear coupling. It is demonstrated that such treatment provides a rapidly convergent representation of the vibrational states of diatomics such as HD + , H 2 , HD, and LiH. On the other hand, the parametric dependence of electron orbitals on the nuclear geometry causes qualitative failure in calculating the dissociative photoionization cross section of H + 2 [39]. Some alternatives of the prolate spheroidal coordinates have been suggested to resolve this issue. In a series of papers [40][41][42][43], Alon and Cederbaum have formulated a multiconfiguration method to describe systems consisting of different types of identical particles with particle conversion taken into account.
Stepping forward in this direction, a fully general TD-MCSCF method [44] has been developed in our group, which can describe the dynamics of a many-body system comprising any arbitrary kinds and numbers of fermions and bosons subject to an external field. The equations of motion (EOMs) of CI-coefficients and spin orbitals are derived for general configuration spaces, based on the time-dependent variational principle. In particular, working equations for molecules subject to intense laser fields are formulated in a very flexible manner, without restriction to the full-CI case, and nuclei can been treated as either classical or quantum particles.
In this paper, as the first step towards the numerical implementation of the general TD-MCSCF method [44], we report our implementation specializing for diatomic molecules in the full-CI case. By transforming to Jacobian coordinates and neglecting the rotational motion of the nuclei, we have only one degree of freedom for the nuclear motion. Thus nuclear orbitals can be defined as functions of internuclear distance. This treatment highly simplifies the underlying numerical procedures and at the same time retains the electron-nuclear correlation. For the electron motion, we use three-dimensional curvilinear coordinates to effectively reduce the computational cost. We apply the present implementation to a hydrogen molecule in intense laser fields. By comparing with the fixed-nuclei case, we address the importance of electronnuclear correlation.
This paper is organized as follows. In Sec. II, we present EOMs of the TD-MCSCF method specialized to diatomic molecules. Numerical implementation is described in Sec. III, and its applications to H 2 and D 2 are described in Sec. IV. Section V concludes this paper. Atomic units are used throughout unless otherwise stated.

A. The system Hamiltonian
We consider a neutral diatomic molecule with N electrons exposed to a laser field linearly polarized parallel with the molecular axis in the z direction. The motions of the two nuclei [with mass and charge (M 1 , Z 1 ) and (M 2 , Z 2 )] and N electrons are characterized by their position vectors R 1 , R 2 , and {r ′ i : i = 1, 2, · · · , N }, respectively, in the laboratory frame. Hiskes [45] has shown that it is possible to separate the center-of-mass motion by choosing the center of mass of the two nuclei as the origin (Jacobian coordinates [46]) and introducing N + 2 new variables: a center-of-mass coordinate R c , a relative nuclear coordinate R and N additional coordinates {r i : i = 1, 2, · · · , N } representing the position vectors of the i-th electron from the center of mass of the two nuclei. For aligned molecules, we neglect the nuclear rotation and retain only vibrational motion. In this case, the internuclear distance vector R can be simply expressed as R ≡ Re z , where e z is the unit vector in the z direction. After several algebraic manipulations (see e.g., Ref. [35]), the Hamiltonian for the internal motion readŝ whereĤ (n) (t),Ĥ (e) (t), andĤ ne) are the nuclear part, electronic part of the Hamiltonian, and electron-nuclear Coulomb interaction, respectively.Ĥ (n) (t) describes the relative motion of the two nuclei with the repulsive Coulomb potential and the coupling with the laser field. Within the dipole approximation, it is given bŷ in the velocity gauge (VG), where µ n = M1M2 M1+M2 is the reduced nuclear mass, and E(t) and A(t) = − with the one-body electronic Hamiltonian either in LG or in VGĥ LG (r, t) = − the electron-electron Coulomb interaction and the mass-polarization term where µ e = M1+M2 M1+M2+1 is the reduced electron mass. The mass-polarization term Eq. (9) describes the deviation of the center of mass of the nuclei from the center of mass of the nuclei plus any N − 1 subset of the N electrons, which is a minor term [47] and omitted in the present paper. Finally, the last term in Eq. (2) represents the electron-nuclear interaction, which has the form It is worth noting that in the fixed-nulcei model, the electron-nuclear Coulomb interaction is trivially a timeindependent quantity which enters to the one-body electron Hamiltonian. However, in the quantum-nuclei case, this term is a two-body operator and needs to be propagated explicitly.
B. The TD-MCSCF method for the diatomic molecule In this subsection, we present the TD-MCSCF method [44] specializing for diatomics in the full-CI case. For a compact representation, the formulation of the theory is given in the second quantization formalism. Hereafter, we use orbital indices {µ, ν, λ, γ, δ} for electronic orbitals, {p, q, r} for nuclear orbitals, and {i, j} for the general (electronic and nuclear) occupied orbitals. For electrons, the Hamiltonian, Eq. (5), can be rewritten as (the masspolarization term is neglected) being the fermionic electron creation and annihilation operators, which change the occupation of spin orbitals {ψ µ × s(σ)} with M e orthonormal electron spatial functions {ψ µ }. The one-and two-electron Hamiltonian matrix elements are defined as and Similarly, the operators, Eq. (3), (4) and (10), in the second quantization notation are expressed aŝ andĤ where (Ê n ) p q =b † pbq andb † p (b p ) is the nuclear creation (annihilation) operator for the set of M n nuclear spatial orbitals {χ p }. The matrix elements, (h n ) p q and (g ne ) pµ qν , are defined as and For the convenience of later discussion, we define the oneand two-body reduced density matrix (RDM) elements as where Ψ is the total molecular wave function defined below in Eq. (19).
In the TD-MCSCF method limited to full-CI expansion, the total electron-nuclear wave function of a diatomic molecule is expanded as where {Φ I } denotes the N -electron Slater determinants constructed from electronic orbitals {ψ µ }, and {χ p } are nuclear orbitals.
The EOMs for the TD-MCSCF method have been derived based on the time-dependent variational principle [48,49], where the action integral S, is made stationary with respect to any possible variations of orbitals {ψ µ } and {χ p } as well as the CI coefficients {C I,p }. The detailed derivation has been given in Ref. [44]. Here in our notation, the resulting EOMs of the CI coefficients read The EOMs of the nuclear and electronic orbitals are given by In Eq. (21)-(23), the operatorsQ (n) = 1 − p |χ p χ p | andQ (e) = 1 − µ |ψ µ ψ µ | are projectors onto the orthogonal complement of the occupied orbital space, which guarantee the orthonormality of orbitals during time propagation. (Ŵ e ) µ ν , (Ŵ ne ) µ ν , and (Ŵ en ) p q are meanfield potentials, given, in the coordinate space, by and The operators,X (n) andX (e) , the matrix elements of which are defined by and determine the components of the time derivation of the orbitals within the occupied orbital space. These two operators impose constraints to the EOMs, which can be arbitrary Hermitian matrices. Depending on the choice of the constraint operators, the EOMs have various equivalent forms. In the present paper, we choose the naturalorbital constraint, i.e., the form of the constraint opera-torsX (n) andX (e) that satisfies the condition to diagonalize the one-body RDMs for both electrons and nuclei [50]. In other words, electronic and nuclear natural orbitals [51] are propagated directly during the interaction with the laser pulse. This form is particularly suited for our propagation scheme applied below. Explicit expressions for the natural-orbital constraint operators together with their variants in imaginary-time propagation for ground-state calculation are given in Appendix A.
For completeness, we also briefly summarize the working equations of TD-MCSCF method for a pure multielectron system within fixed nuclei in Appendix B.

III. IMPLEMENTATION
One of the important aspects for the implementation of TD-MCSCF method using real-space grid-based techniques is the choice of coordinate system. Many existing implementations of TD-MCSCF method are optimized to efficiently study the strong-field phenomena in atoms and diatomic molecules in the presence of either linearly or elliptically polarized laser pulse, exploiting the underlying symmetries of the system Hamiltonian with, for instance, the spherical [25,26,28,52], cylindrical [37], or prolate spheroidal coordinates [38]. To retain the possibility of carrying out first-principle simulations of the response of molecules to intense laser pulses without assuming symmetry, we do not pursue this direction here.
The most straightforward way is using Cartesian coordinates with an equal grid mesh. While this approach offers better opportunities for parallel computing due to the highly sparse and structured representation of the Hamiltonian, it suffers from drawbacks of requiring small grid spacings to accurately represent the electron-nucleus interaction and a large number of grid points to sustain photoelectrons (or dissociated nuclei). One way to overcome this problem is using multiresolution Cartesian grids [53], which has previously been demonstrated for the implementation of MCTDHF method. Alternatively, we can use curvilinear coordinates [54][55][56]. By choosing an appropriate coordinate transformation, it is possible to achieve a high spatial resolution in the vicinity of nuclei and, at the same time, a sufficiently large simulation box. Furthermore, the coordinate transformation can be given in analytical form. In addition, the EOMs can be analytically transformed. The transformed EOMs are discretized using higher-order finite difference methods on an equal grid mesh in curvilinear coordinates. Such a treatment overcomes the shortcomings of the equal grid mesh while retaining its advantages. This scenario is demonstrated for the implementation of the time-dependent Hartree-Fock (TDHF) method elsewhere. In Sec. III A, we describe the coordinate transformation for the present case.
Another important issue is the time propagation method. The orbital EOMs Eq. (22) and (23) consist of stiff linear terms and nonstiff nonlinear terms. The stiffness of the linear part comes from the one-body kinetic-energy operators while the nonlinear terms contain dependencies on the CI coefficients and the orbitals. For an efficient propagation of the EOMs, we adopt the exponential time differencing second-order Runge-Kutta scheme (ETDRK2) [57,58], as described in Sec. III B.

A. Curvilinear Coordinate
Here we mainly focus on the electron coordinate transformation and briefly mention the nuclear transformation at the end of this subsection. Let us consider the physical space to be described by Cartesian coordinates r = (x, y, z) = (r 1 , r 2 , r 3 ) and an analytical transformation r = f (ξ) to curvilinear coordinates ξ = (ξ 1 , ξ 2 , ξ 3 ) that act as the computational space. This transformation is single-valued and at least twice continuously differentiable on a three-dimensional domain. The Jacobian of the transformation J is defined as a 3 × 3 matrix of which the (i, j) matrix element is J i j = ∂r i /∂ξ j , with |J| = detJ being its determinants. The metric tensor g i j in the curvilinear coordinates is and the Laplacian operator is given by It is beneficial to transform the electron orbitals as such that the inner product obeys Inserting Eq. (31) into Eq. (23), any electron operatorĥ in the EOMs of transformed electron orbitals will change according toĥ −→ |J|ĥ 1 √ |J| . The transformed Laplacian reads Special care must be taken to keep the resulting Laplacian symmetric after applying the finite difference discretization [56]. In our implementation, we rewrite Eq. (33) as where All the terms in Eq. (33) involving the first derivatives are now rewritten so that they only involve the second derivatives in Eq. (34). When applying the finite difference method, the symmetry of the resulting Laplacian is preserved. Other operators are formulated in a similar way.
The z component of the transformed momentum operator is expressed as The electron mean-field potential Eq. (24) is evaluated by solving the corresponding Poisson equation , in the same way as in Eq. (31), the Poisson equation is transformed as With all the transformed operators available, the transformed EOMs can be discretized on the computational space ξ = (ξ 1 , ξ 2 , ξ 3 ) by equally spaced grids. The fifth-order central finite difference is applied for spatial differential operators. To incorporate the exterior complex scaling (ECS) absorbing boundary [26,[59][60][61][62], we use an orthogonal transformation which smoothly scales the real ξ i to complex r i . The explicit expression for this transformation reads where the computational space along one direction is divided into five regions: the inner region (−ξ (0) ≤ ξ ≤ ξ (0) ), two intermediate regions (−ξ (1) ≤ ξ < −ξ (0) and ξ (0) < ξ ≤ ξ (1) ), and two outer regions (ξ < −ξ (1) and ξ (1) < ξ). The parameters ∆ s and ∆ t are the scaling factors in the inner and outer regions, respectively. The polynomials are used in the intermediate regions, where parameters a p and b p = (−1) p a p are determined to satisfy the conditions of continuity of f (ξ) and its first and second derivatives at the boundaries. θ is the ECS scaling angle.
Special care should be taken for solving the Poisson equation when analytic continuation is applied for the orbitals. In this case, solving the Poisson equation in the complex-scaled region (|ξ| > ξ (0) ) turns out to be numerically unstable [26]. In the present paper, we approximately neglect the Coulomb interaction between electrons in the scaled region and adopt octupole expansion for the calculation of the electron mean-field potential in the scaled region as well as the boundary condition of solving the Poisson equation in the unscaled region.
The nuclear coordinate is transformed in a similar way as that of electrons, i.e., R = f (ζ), together with orbital transformationχ p (ζ) = √ f ′ χ p (ζ). The transformed EOMs of nuclear orbitals are numerically discretized on equally spaced points in ζ and readily propagated. The same transformation function in Eq. (40) is adopted except for the condition ζ > 0.

B. Time propagation scheme
We rewrite the electronic orbital EOMs Eq. (23) as a nonlinear differential equation whereĥ is a stiff linear operator and N [t, u(t)] a nonstiff nonlinear reminder. We first consider the stiff part to be time-independent, i.e.,ĥ = − 1 2µe ∇ 2 , and put the laser interaction term into the nonlinear part. Eq. (41) can be integrated over a small time interval [t n , t n+1 = t n + ∆t] as which is exact for time-independentĥ. This scheme is called the exponential integrator [58]. By approximating N [t ′ , u(t ′ )] in the integrand using second-order Runge-Kutta (RK2) method, we obtain the expressions of ET-DRK2 method and where the ϕ-functions are defined as We use the EXPOKIT package [63] to efficiently calculate the action of ϕ functions on the operand vectors in Eq. (43) and (44). The modification for time-dependent stiff operatorĥ(t) is given in Appendix C. The nuclear EOM Eq. (22) is propagated in the same manner. For CI EOM Eq. (21), we treat the totality of the right-hand side as a nonlinear term and time propagation of CI reduces to the RK2 scheme.

IV. NUMERICAL RESULTS AND DISCUSSIONS
A. Ground-state properties of H2 To assess the performance of the method and our implementation described in Sec. II and III, we first do a series of ground-state calculations taking the H 2 molecule as an example, either with fixed-nuclei approximation or fully quantum mechanically with nuclear motion (hereafter refer to as quantum nuclei). We assume a H 2 molecule is aligned along the z axis. The imaginarytime propagation method is adopted to obtain the ground state. The grid resolution is optimized for ground-state calculations in order to have a better comparison with the existing literature values of ground-state properties such as ground-state energy and natural occupation number. In the quantum-nuclei case, we use ∆ξ = 0.08 a.u. and box size of ξ max = 8 a.u. for all the three directions of the electronic coordinates, and ∆ζ = 0.05 a.u. and total length of ζ max = 6 a.u. for nuclear coordinates. Curvilinear transformation is not used in this subsection. In the fixed-nuclei case, the same spatial grid is used for electrons, and nuclei are fixed at equilibrium internuclear distance R = 1.4 a.u.
Let us first consider the case of the fixed-nuclei H 2 molecule. In Table I, we present the ground-state energies for different numbers M of electron orbitals. Up to eight orbitals with six σ and two degenerate π orbitals are used. Clearly, with increasing M , the ground-state energy consistently converges to the exact value [64] obtained by the direct solution of two-electron Schrödinger equation, and only a few orbitals are needed to obtain good agreement. The energy values in this paper (third column) have a higher precision than those from Ref. [53] (fourth column) due to finer spatial resolution. For instance, the Hartree-Fock energy presented here is nearly identical to that calculated using the finite-element discrete variable representation basis set in Ref. [38]. In the fifth column, we also give the fraction of the correlation energy defined as  The total wave function Ψ(t) of the fixed-nuclei H 2 molecule can be expressed by using different types of electron orbitals. In the present paper, we use naturalorbital representation as explained in Appendix A, where the one-body RDM is kept diagonal during the propagation. In this way, natural orbitals and natural occupation numbers, defined as eigenvectors and eigenvalues of the one-body RDM, are directly obtained after the propagation. Taking M = 8 for example, we present the natural occupation numbers in Table II. The natural occupation numbers represent the distribution of all electrons among natural orbitals. Their sum is equal to the number of electrons (2 for a hydrogen molecule). As can be seen, the first natural orbital is highly occupied with an occupation number close to 2 while for higher orbitals the occupation numbers decrease rapidly. The two degenerate π orbitals have exactly the same occupation numbers, as expected. Our results are in excellent agreement with those in Ref. [38] calculated by prolate spheroidal coordinate code, validating our implementation. The  two-dimensional slices of natural orbitals are displayed in Fig 1, where orbital symmetries can be clearly identified. The information of orbital symmetries can guide us to study the time-evolution behavior of the orbitals in the response to the laser field. For instance, the two π orbitals in Fig 1(d) and 1(e) will respond in the same fashion to a laser field polarized parallel to the molecular axis.
Next, we examine the ground-state properties of the quantum-nuclei H 2 molecule. In the present paper, we  Table IV, with reference values from Ref. [38] for comparison. The corresponding electronic natural orbitals are shown in Fig. 2. We see in Fig. 1 and Fig. 2 that the overall shapes of the two sets of natural orbitals are remarkably similar, although small differences can be viewed, for instance, in the 4σ and 6σ orbitals. Unlike the treatment in Ref. [38], where the electron orbitals have parametric dependence on nuclear coordinate R through the prolate spheroidal coordinate system, in our electronnuclear wave function expansion Eq. (19), the electronic and nuclear primitive orbital functions are strictly independent without parametric dependence. The similarity of the natural orbitals in Fig. 1 and Fig. 2 thus suggests that the strong correlation between the positions of the nuclei and electrons due to their Coulomb attraction is correctly accounted for in our calculation.

B. Electron-nuclear dynamics of H2 in an intense laser field
In this subsection, we present numerical applications of the implementation of our method to H 2 subject to an intense laser pulse linearly polarized along the molecular axis in the z direction, both in the fixed-nuclei and quantum-nuclei case. We adopt VG and assume the vector potential of the laser field of the following form: where I 0 is the peak intensity, T the laser period, ω = 2π/T the angular frequency and N 0 the total number of  For time-dependent calculations, a grid width larger than that used in the ground-state calculation is employed so that the computational time falls in an affordable range. We use grid spacing ∆ξ = 0.2 a.u. for all the three directions of the electronic coordinates and   Table V. The ETDRK2 method described in the previous section is used to propagate EOMs with 10 000 time steps per optical cycle. We have confirmed the convergence of the results with respect to spatial and temporal resolutions. We first simulate a fixed-nuclei H 2 molecule subject to a near infrared laser field with a wave length of 800 nm, 3 × 10 14 W/cm 2 intensity, and total duration of 8 optical cycles. For such a process, the direct exact numerical solution of the six-dimensional TDSE is still unavailable to our knowledge, due to the unfavorable scaling of the problem size with the laser wavelength. So we choose the result of the largest configuration space considered in our time-dependent simulation, M = 7, as a reference.
The simulation results are shown in Figs. 3 and 4. We present here the time evolution of the n-electron ionization probability (Fig. 3) and the high harmonic generation (HHG) spectra (Fig. 4). The n-electron ionization probability, P n (n = 0, 1, 2), is conveniently defined as a probability to find n electrons located outside a cube with a side length of 10 a.u, which can be evaluated using the method introduced in Ref. [24]. The HHG spectra are obtained as the modulus squared of the Fourier transform of the dipole acceleration, which, in its turn, is evaluated by double numerical differentiation of the dipole moment with respect to time. Figure 3 plots the survival probability P 0 , the single-ionization probability P 1 , and double-ionization probability P 2 with different number of orbitals, respectively. While P 0 and P 1 are converged for M ≥ 3, the results with M = 1 and 2 show strong underestimations, which indicates the electron correlation plays a nonnegligible role during the interaction with the laser pulse. For the double ionization probability, in great contrast, we can see from Fig. 3(c) that it is not fully converged even with M = 7 and still varies after the end of the laser pulse (t > 8T 0 ). This observation indicates that the double ionization is much more sensitive to the electron-electron correlation than single ionization and HHG (see below) [66,67]. The high-harmonic spectra calculated with different numbers of orbitals are shown in Fig. 4. The TDHF (M = 1) result already provides the correct cutoff position but underestimates the harmonic intensity in the plateau region. As the number of orbitals is increased, the TD-MCSCF method provides improved description of electron-electron correlation, and thus the harmonic spectrum displays a fairly rapid convergence with respect to M . In particular, HHG spectra with M = 6 and 7 show remarkable agreement with each other on the scale of the figure.
Next, let us turn to the quantum-nuclei case. The laser parameters are the same as in the fixed-nuclei case. We first analyze the time-dependent nuclear probability density, defined as the modulus squared of the total wave function integrated over the electronic coordinates.   tion of electron-nuclear correlation, and thus molecular vibrations as well as the dissociative ionization can be well described by the TD-MCSCF method with a sufficient number of orbitals. Figure 6 shows n-electron ionization probabilities of a quantum-nuclei H 2 molecule. For this calculation, we have confirmed that the simulation box is large enough to hold the nuclear density, thus the nuclear coordinate can be safely integrated out. The results converge to that with M = 7 (taken as a reference) as the number of orbitals increases. However, quantitatively, the ionization probability differs from that of the fixed-nuclei model (Fig. 3) under the same laser condition. In Fig. 7, we show the single ionization probability of fixed-nuclei H 2 , quantum-nuclei D 2 , and quantum-nuclei H 2 with M = 7. It is found that the single ionization probability of quantum-nuclei H 2 is higher than that of D 2 , and both of them are higher than that of fixed-nuclei H 2 . This result suggests the effect of enhanced ionization due to nuclear motion [68,69]. Figure 8 presents the TD-MCSCF HHG spectra from the quantum-nuclei H 2 molecule. The results show an improvement with increasing number of orbitals as in the fixed-nuclei case. However, the harmonic spectrum is not fully converged even with M = 7. A similar trend with quantum treatment of nuclei has previously been reported for reduced-dimensional H + 2 [44]. Nevertheless, the quantitative agreement is already sufficient to address the effects of nuclear motion during the HHG process. While the present TD-MCSCF method treats the nuclear dynamics in a full quantum way, Saenz et al. [70][71][72][73][74] have proposed to take nuclear vibration into account by the "frozen-nuclei approximation", where the nuclear wave function is frozen to the vibrational ground state χ(R) of the electronic Born-Oppenheimer ground-state potential during the time propagation. Within this approximation, the R-integrated ionization probability can be calculated as the sum of fixed-nuclei signals weighted by the vibrational nuclear density [74], where P FN (t; R) is the fixed-nuclei time-dependent ionization probability. We calculate P FN (t; R) using fixednuclei MCTDHF method with M = 7 for a range of internuclear distances where the vibrational wave function χ(R) is nonvanishing, i.e., R = 1.0 − 2.5 a.u. The averaged ionization probability is compared with the quantum-nuclei and fixed-nuclei results in Fig. 9. On one hand, the averaged single-ionization probability is slightly higher than the fixed-nuclei one at the equilibrium internuclear distance, reflecting ionization enhanced at larger values of R. On the other hand, quantitatively, there is a substantial discrepancy between the averaged results and those from the full quantum calculation. This observation indicates that the dynamical aspects of the electron-nuclear correlation play an important role during the interaction with the laser pulse. To further examine the importance of dynamical electron-nuclear correlations, let us calculate the coherently averaged HHG spectrum. We first obtain the timedependent dipole acceleration averaged over internuclear distance, where a FN (t; R) denotes the fixed-nuclei dipole acceleration, calculated with the TD-MCSCF method with M = 7. Then, the averaged harmonic spectrum is obtained as the modulus squared of the Fourier transform of a(t). The averaged spectrum is clearly different from the full quantum result and close to the fixed-nuclei result (Fig. 10), due to nuclear dynamics during the quiver motion of the electron. These results unambiguously demonstrate the importance of dynamical electron-nuclear correlations and the significance of our full quantum treatment. Experimentally, it has been observed that highharmonic emission is stronger in D 2 than that in H 2 due to the faster nuclear motion in the lighter H 2 molecule [75,76]. We will show that the present TD-MCSCF method, which allows full quantum treatment of nuclei, can well reproduce the experimental results. We 11  consider an 800-nm laser field with a peak intensity of 2 × 10 14 W/cm 2 . The laser parameters are similar to those in the experiments. In order to reveal the odd order harmonic peaks, we use a multi-cycle laser pulse with total duration of 30 periods. The HHG spectra of quantum-nuclei H 2 and D 2 are calculated with M = 7, as shown in Fig. 11(a). Up to the 29th harmonics can be well resolved in the spectra. In Fig. 11(b), the spectra of H 2 and D 2 in the plateau region are plotted on a linear scale. It can be clearly seen that the harmonic signal is stronger in D 2 than that in H 2 . Fig. 11(c) shows the intensity ratio between D 2 and H 2 as a function of harmonic order. The ratio is calculated by using the peak values of odd order harmonics. The ratio is larger than unity at all orders and increases with the harmonic order. Our calculation shows quantitative agreement with the experimental measurement in Ref. [76] under the same laser condition. The experimental result from Ref. [75] using a few-cycle laser pulse also exhibits a similar trend, albeit with relatively small values compared to our calculation possibly because of the difference in the laser parameters. The present calculation using the TD-MCSCF method nicely reproduces the experimentally observed isotope effects on HHG, emphasizing the importance of dynamical electron-nuclear correlation during HHG process.

V. SUMMARY
We have presented the numerical implementation of the recently formulated TD-MCSCF method [44] for diatomic molecules in the full-CI expansion. By full quantum treatment of the nuclei, the present implementation allows one to investigate coupled electron-nuclear dynamics of molecules subject to intense laser fields. As a first numerical test, we have applied this method to the ground state as well as the laser-driven dynamics of H 2 molecule. The ground-state properties, including ground-state energy, natural orbitals and natural occupation numbers show remarkable agreement with benchmark data [38] available in the literature, which indicates the correctness of our implementation. For laser-driven dynamics, we have shown that the method can systematically improve the accuracy for describing the highly nonlinear HHG phenomenon by increasing the number of electronic and nuclear orbital functions. Furthermore, the TD-MCSCF method with sufficiently large number of orbitals has been applied to study the isotopic effects of HHG to justify experimental observations. Our results indicate that HHG is sensitive to the laser-induced nuclear vibrational motion. The harmonic emission is more intense in heavier isotopes and the role of nuclear motion in suppressing the intensity of high harmonics is increased with increasing harmonic order.
While the present paper has been mainly focused on the electron-nuclear dynamics in HHG, we have also shown that the TD-MCSCF method is capable of de-scribing photoionization and dissociation [37,77,78] of molecules interacting with a laser field. In order to properly explain experimental results such as electron-ion coincidence measurements [79,80] during molecular dissociation, extending the functionality of the present implementation to calculate momentum distribution of ionized electrons and the joint electron-nuclear-energy spectrum will be the further aspects.
where laser-electron interactionV ext within the dipole approximation either in the LG or in the VG is given bŷ V VG ext (r, t) = −iA(t) · ∇.
In the TD-MCSCF method, the N -electron wave function is expressed by a linear combination of Slater determinants as where the Slater determinant Φ I (t) is built from M oc- whereQ = 1 − µ |ψ µ ψ µ | is the projector against the occupied orbital space, D and P the one-and twoelectron RDMs,Ŵ µ ν the mean-field potential with the same definition as in Eq. (24), andX the constraint operator defined in Eq. (28). Here we also use the naturalorbital constraintX as derived in Appendix A.