Time-dependent optimized coupled-cluster method for multielectron dynamics III: A second-order many-body perturbation approximation

We report successful implementation of the time-dependent second-order many-body perturbation theory using optimized orthonormal orbital functions called time-dependent optimized second-order many-body perturbation theory [TD-OMP2] to reach out to relatively larger chemical systems for the study of intense-laser-driven multielectron dynamics. We apply this method to strong-field ionization and high-order harmonic generation (HHG) of Ar. The calculation results are benchmarked against the time-dependent complete-active-space self-consistent field (TD-CASSCF), time-dependent optimized coupled-cluster double (TD-OCCD), and the time-dependent Hartree-Fock (TDHF) methods to explore the role of electron correlations.


I. INTRODUCTION
Laser-driven multielectron dynamics has become an active area of research, thanks to the remarkable advance in laser technologies, which has made it possible to measure and control the electronic motion [1][2][3][4][5][6][7][8][9][10] . Atoms and molecules interacting with laser pulses of intensity 10 14 W/cm 2 in the visible to mid-infrared spectral range, exhibit highly, even nonperturbatively nonlinear response such as above-threshold ionization (ATI), tunneling ionization, nonsequential double ionization (NSDI), and high-order harmonic generation (HHG). The HHG process is one of the key elements in the study of light-matter interaction in the attosecond time-scale, delivering ultrashort coherent light pulses in the extremeultraviolet (XUV) to the soft x-ray regions, which carry the information on the electronic structure and dynamics of the generating medium. 11 The HHG spectra are characterized by a plateau where the intensity of the emitted radiation remains nearly constant up to many orders, followed by an abrupt cutoff.
In principle, the multielectron dynamics and electron correlation [12][13][14][15] are exactly described by the timedependent Schrödinger equation (TDSE). However, direct numerical integration of TDSE is not feasible for systems with more than two electrons [16][17][18][19][20][21][22][23][24][25] . As a result, single-active-electron (SAE) approximations has been widely used 26,27 , in which only the outermost electron is explicitly treated under the effect of the other electrons modeled by an effective potential. Whereas SAE has been useful in numerically exploring different strongfield phenomena, the electron correlation is missing in this approximation. a) Electronic mail: pathak@atto.t.u-tokyo.ac.jp b) Electronic mail: sato@atto.t.u-tokyo.ac.jp c) Electronic mail: ishiken@n.t.u-tokyo.ac.jp Therefore, various tractable ab initio methods have actively been developed for theoretical description of correlated multielectron dynamics in intense laser fields.
Among the most reliable approaches to serve the purpose are the multiconfiguration timedependent Hartree-Fock (MCTDHF) method [28][29][30][31][32] and time-dependent complete-active-space self-consistentfield (TD-CASSCF) method 33,34 . In MCTDHF, the electronic wavefunction is expressed in terms of full configuration interaction (FCI) expansion, Ψ(t) = I C I (t)Φ I (t), where both CI coefficients {C I (t)} and occupied spin-orbitals {ψ p (t)} constituting Slater determinants {Φ I (t)} are propagated in time. The TD-CASSCF method flexibly classifies occupied orbital space into frozen-core (doubly occupied and fixed in time), dynamical-core (doubly occupied but propagated in time) and active (fully correlated and propagated in time) subspaces. Though accurate and powerful, the computational cost of the MCTDHF and TD-CASSCF methods scales factorially with the number of correlated electrons.
More approximate and thus less demanding methods such as the time-dependent restricted-activespace self-consistent field (TD-RASSCF), [35][36][37] and timedependent occupation-restricted multiple-active-space (TD-ORMAS) 38 have been developed by further flexibly classifying active orbital sub-space to target larger chemical systems by limiting CI expansion of the wavefunction up to a manageable level. The TD-RASSCF and TD-ORMAS methods achieve a polynomial, instead of factorial, cost scaling, and state-of-the-art realspace implementations have turned out to be of great utility. 34 To regain the size-extensivity at the same level of truncation, we have recently derived and numerically implemented a time-dependent optimized coupled-cluster (TD-OCC) method 43 using optimized orthonormal orbitals where both orbitals and amplitudes are timedependent. Our TD-OCC method is the time-dependent formulation of the stationary orbital optimized coupledcluster method. 44,45 We have implemented the TD-OCC method with up to triple excitation amplitudes (TD-OCCD and TD-OCCDT), and applied it to multielectron dynamics in Ar induced by a strong laser pulse to obtain a good agreement with the fully-correlated TD-CASSCF method within the same active orbital space. In earlier work, Kvaal reported 46 an orbital adaptive coupledcluster (OATDCC) method built upon the work of Arponen using biorthogonal orbitals; 47 though, applications to laser-driven dynamics have not been addressed. The polynomial scaling TD-OCC 43 or the OATDCC 46 method can reduce the computational cost to a large extent in comparison to the general factorially scaling MCTDHF methods.
As a cost-effective approximation of the TD-OCCD method, we have recently developed a method called TD-OCEPA0 48 based on the simplest version of the coupled-electron pair approximation, [49][50][51][52][53][54] or equivalently the linearized CCD approximation, 55 popular in quantum chemistry. The computational cost of the TD-OCEPA0 method scales as N 6 , with N being the number of active orbitals. Although this is the same scaling as that of the parent TD-OCCD method, TD-OCEPA0 is much more efficient than TD-OCCD due to the linearity of amplitude equations and avoidance of a separate solution for the de-excitation amplitudes, resulting from the Hermitian structure of the underlying Lagrangian. 48 To enhance the applicability to even larger chemical systems, we are looking for further approximate versions with a lower computational scaling in the TD-OCC framework. The coupled-cluster method is intricately connected with the many-body perturbation theory. It allows one to obtain finite-order perturbation theory energies and wavefunction from the coupled-cluster equations. 56 The computation of the second-order energy requires the first-order wavefunction, and only doubly excited determinants have contributions in the first-order correction to the reference wavefunction. Thus, secondorder many-body perturbation theory (MP2) can be seen as an approximation to the coupled-cluster double (CCD) method, having a lower N 5 scaling.
The MP2 method with optimized orbitals has been developed by Bozkaya et al. 57 for the stationary electronic structure calculations by minimization of the socalled MP2-Λ functional. In earlier work, the optimization of the MP2 energy was based on the minimization of the Hylleraas functional with respect to the orbital rotation. 58,59 While both of these techniques provide identical energy at the stationary point, the Λ functionalbased derivation has an advantage that it can be eas-ily extended to higher-order many-body perturbation theory. 60 In this article, we propose time-dependent orbital optimized second-order many-body perturbation theory (TD-OMP2) as an approximation to the TD-OCCD method, based on the time-dependent (quasi) variational principle. The TD-OMP2 method inherits important features of size extensivity and gauge invariance from TD-OCC, with a reduced computational cost of O(N 5 ). As a first test case, we apply the TD-OMP2 method to electron dynamics in Ar atom irradiated by a strong laser field, and compare the results with those computed at TDHF, TD-OCCD, and TD-CASSCF levels to understand where really TD-OMP2 stands in describing the effects of correlation in the laser-driven multielectron dynamics. It should be emphasized that the present TD-OMP2 method, albeit named after a perturbation theory, can be applied to laser-induced, nonperturbative electron dynamics, since the laser-electron interaction is fully (nonperturbatively) included in the zeroth-order description with time-dependent orbitals.
This paper is organized as follows. Our formulation of the TD-OMP2 method is presented in Sec. II. Numerical applications are described and discussed in Sec. III. The concluding remarks are given in Sec. IV. Hartree atomic units are used unless otherwise stated, and Einstein convention is implied throughout for summation over orbital indices.

II. TD-OMP2 METHOD
Let us consider the electronic Hamiltonian H of the following form, where N e is the number of electrons, r i and p i are the position and canonical momentum, respectively, of electron i, h 0 is the field-free one-electron Hamiltonian, and V ext is the laser-electron interaction. The Hamiltonian in the second-quantization notation can be written as, is the creation (annihilation) operator for the set of orthonormal 2n bas spin-orbitals {ψ µ }, with n bas being the number of basis functions (or grid points) to represent the spatial part of ψ µ . The operator matrix elements h µ ν , u µγ νλ , and v µγ νλ are defined as and The full set of spin-orbitals (labeled with µ, ν, λ, γ) are divided into n occ occupied (o, p, q, r, s) and 2n bas − n occ virtual spin-orbitals, respectively. The occupied orbitals are further divided into n core core, and n act = n occ − n core active spin-orbitals (t, u, v, w). The active electrons are correlated in the active orbital subspace, which is further split into those in the hole space (i, j, k, l) and the particle space (a, b, c, d), which are occupied and unoccupied, respectively, with respect to the reference determinant. The core spin-orbitals can further be split into frozen-core space (i , j ) and the dynamicalcore space (i , j ). 33

A. Review of TD-OCCD
The stationary MP2 method can be considered as an approximation of the CCD method, where the full CCD energy functional, is approximated by retaining the terms giving the secondorder correction to the reference energy E 0 = Φ|Ĥ|Φ .
Here,T 2 = τ ab ijÊ ab ij andΛ 2 = λ ij abÊ ij ab , with τ ab ij and λ ij ab being the excitation and de-excitation amplitudes, respectively. Therefore, we start with the time-dependent CCD Lagrangian of the following form, which is a natural time-dependent extension of the energy functional, Eq. (7). Following Ref. 43, we consider the real-valued action functional, and make it stationary, δS = 0, with respect to the variation of amplitudes δτ ab ij , δλ ij ab and variations of orthonormality-conserving orbitals δψ ν . The equations of motion (EOMs) for the amplitudes are obtained as whereX = X µ νĉ † µĉν with X µ ν = ψ µ |ψ ν , and those for orbitals as, where ρ p q and ρ pr qs are the one-and two-body reduced density matrices, respectively, defined as ρ qs pr = Φ|(1 +Λ 2 )e −T2Êpr qs eT 2 |Φ .
B. Derivation of TD-OMP2 as an approximation to TD-OCCD Now we derive the TD-OMP2 method as an approximation to the TD-OCCD method, based on the partitioning of the electronic Hamiltonian, into the zeroth-order partĤ (0) =f = f p qÊ p q and the perturbationĤ (1) where (h 0 ) µ ν and (V exe ) µ ν are the matrix elements ofĥ 0 andV ext , respectively.
Based on this partitioning, we apply the Baker-Campbell-Hausdorff expansion to TD-OCCD Lagrangian of Eq. (8), and retain those terms up to quadratic in v, τ , and λ (thus, contributing to first-and second-order corrections to Lagrangian) to obtain Inserting this TD-OMP2 Lagrangian into Eq. (9) and making it stationary with respect to amplitude variations derives TD-OMP2 amplitude equations, where p(µν) is the cyclic permutation operator. Importantly, Eqs. (20) and (21) reveals that the EOM for λ ij ab is just complex conjugate of that for τ ab ij , resulting in λ ij ab = τ ab * ij . Therefore, we do not need a separate solution for the λ 2 amplitudes.
We also make the action stationary with respect to the orthonormality-conserving orbital variation to derive formally the same orbital EOMs as Eqs. (12) and (13), with one-particle reduced density matrices (1RDM) and two-particle reduced density matrices (2RDM) given explicitly as ρ qs pr = (ρ 0 ) qs pr + Γ qs pr , where (ρ 0 ) q p = δ q j δ j p and (ρ 0 ) qs pr = γ q p δ s j δ j r + γ s r δ q j δ j p − γ q r δ s j δ j p − γ s p δ q j δ j r + δ q j δ j p δ s k δ k r − δ s j δ j p δ q k δ k r are the reference contributions, and non-zero elements of γ q p and Γ qs pr are given by In summary, the TD-OMP2 method is defined by the EOMs of double excitation amplitudes τ ab ij [Eq. (20)], the relation λ ij ab = τ ab * ij , and the EOMs of orbitals [Eqs. (12) and (13)] with 1RDM and 2RDM elements given by Eqs. (23)- (26). The orbital time-derivative terms −iX can be dropped from Eq. (22), if one makes an arbitrary choice of ψ i |ψ j = ψ a |ψ b = 0 for the redundant orbital rotations. It should be noted that (i) our partitioning scheme [Eqs. (18)] is consistent with the standard Møller-Plesset perturbation theory in the absence of the external field V ext , and (ii) in case with V ext , the zerothorder Hamiltonian is time dependent, both through the change of orbitals and an explicit time dependence of V ext (t), the latter implying nonperturbative inclusion of the laser-electron interaction.

C. Alternative derivation
Another, more perturbation theoretic derivation begins with the following Lagrangian, where the subscript c implies restriction to connected terms, andf N andv N are the normal-ordered part off andv, respectively. Then we follow the procedure of time-dependent variational principle based on the action of Eq. (9), to obtain identical method as derived in the previous section. The Lagrangian of Eq. (27) can be viewed as the time-dependent extension of the Hylleraas energy functional used in conventional OMP2 method. 58,59 Two expressions of the TD-OMP2 Lagrangian, Eqs. (19) and (27), take the same numerical value as a function of time when evaluated with the solution of TD-OMP2, {τ ab ij (t), ψ p (t)}. The latter expression, Eq. (27), makes it clearer that the laser-electron interaction, included inf N , is treated nonperturbatively.

A. Ground-state energy of BH
To assess the performance of the method described in the previous section, we do a series of the groundstate energy calculations taking BH molecule as an example with double-ζ plus polarization (DZP) basis. 63 The imaginary time relaxation method 43 is used to obtain the ground state. We have started our calculations with the bond length of 1.8 a.u. and gradually increased to 7.0 a.u., beyond which we could not achieve convergence. The values for both the MP2 and OMP2 are reported in Table I. The required matrix elements are obtained from the Gaussian09 program 62 package and interfaced with our numerical code. All the values are compared with the values from the PSI4 program package. 61 We obtained identical results for the reported values except in five cases, even for which the difference appears only at the eighth or ninth digit after the decimal point. In these calculations, the number of orbitals N are taken to be the same as the number of basis n bas to make a comparison with PSI4, 61 and all the orbitals are treated as active. However, our implementation allows to optimize orbitals with N ≤ n bas in general, with a flexible classification of the occupied orbital space into frozen core, dynamical core, and active.

B. Ar in a strong laser field
We present numerical applications of the present method to Ar subject to an intense laser pulse linearly polarized in the z direction. Calculations for intense longwavelength laser fields are computationally demanding, thus serving as a good test for the newly implemented methods. The laser-electron interaction is introduced to the one-body part of the electronic Hamiltonian within the dipole approximation in the velocity gauge, where Z is the atomic number, the vector potential, with E(t) being the laser electric field. While our method is gauge-invariant, we obtain faster convergence with the velocity gauge for the case of intense long-wavelength laser pulses 34,40 . The laser electric field is of the following form for 0 ≤ t ≤ 3T , and E(T ) = 0 otherwise, with the central wavelength λ = 2π/ω 0 = 1200 nm, the period T = 2π/ω 0 ∼ 4.00 fs, and the peak intensity I 0 = E 2 0 . We have considered three different intensities of 1, 2, and 4 ×10 14 W/cm 2 for Ar. The spherical finite-element discrete variable representation (FEDVR) basis 34,40 is used in our implementation. The convergence with respect to the maximum angular momentum l max is checked at the TDHF level, and l max is set to 72. The radial coordinate r is discretized by FEDVR consisting of 78 finite elements with 23 DVR functions each, to support 0 < r < R max = 300. A cos 1 4 mask function is switched on at 240 to avoid any reflection from the boundary. Fourth-order exponential Runge-Kutta integrator 64 is used to propagate equations of motions with 20000 time steps per optical cycle. The simulations are continued after the end of pulse for further 6000 time steps. The 1s2s2p core is kept frozen, and the remaining eight electrons are distributed among thirteen active orbitals.
In Fig. 1, we plot the dipole moment as a function of time, evaluated as a trace ψ p |ẑ|ψ q ρ q p . We compare  the TD-OMP2 results with the TDHF, TD-CASSCF, and TD-OCCD ones. With the same active space, TD-CASSCF produces highly accurate results, useful for performance analysis of the TD-OMP2 method. Whereas the initial dipole moment is zero due to the spherical symmetry of the atom, it starts to oscillate with the application of the laser field. At 10 14 W/cm 2 all the methods produce similar results. The TD-OCCD produces virtually the identical result with the TD-CASSCF, whereas TD-OMP2 slightly overestimates, and TDHF underestimates, considering TD-CASSCF as the benchmark. With increase in intensity, the difference among the four methods become more prominent. While all the methods give similar results in the early stage, TD-OMP2 and TDHF start to overestimate and underestimate, respectively, with the progress of tunneling ionization (Fig. 2). TD-OCCD slightly underestimates. In general, the performance of the TD-OMP2 method is better than the TDHF due to consideration at least a part of the electron correlation.
The general trends in Fig. 1 are also found in the single ionization probability (Fig. 2), evaluated as the probability of finding an electron outside a sphere of 20 a.u. radius. Again, we see overestimation by TD-OMP2 and underestimation by TDHF in comparison to the TD-CASSCF result; the performance of TD-OMP2 is better than that of TDHF. The deviation of TD-OCCD from TD-CASSCF is relatively small but increases with an increase in intensity, since higher excitations are important for the description of substantial ionization.
In Fig. 3 we compare HHG spectra, calculated as the modulus squared of the Fourier transform of the expectation value of the dipole acceleration, which, in turn, is obtained with a modified Ehrenfest expression 34 . All the four methods reproduce the HHG spectra relatively well, including an experimentally observed characteristic dip around the 52nd order (∼ 54 eV) at 2× and 4 × 10 14 W/cm 2 related to the Cooper minimum of the photoionization spectrum 65 at the same energy. However, TDHF systematically underestimates and fails to reproduce fine details. The agreement with the TD-CASSCF results is much better for the TD-OCCD method, which contains nonlinear terms in the amplitude equations, then followed by TD-OMP2 with slight overestimation. This trend is consistent with the capability of each method to treat the electron correlation. It is worth noting that, for an intensity as high as 4 × 10 14 W/cm 2 , the TD-OMP2 simulation completed stably. This should be the direct consequence of full inclusion of the laser-electron interaction in the zerothorder Hamiltonian and the orbital optimization (propagation) according to the time-dependent variational principle based on the total (up to second-order) Lagrangian, which keeps the instantaneous perturbationĤ (1) =Ĥ −f relatively small in the present simulation.
Finally, the computational cost of different parts of TD-OMP2 and TD-OCCD methods are compared in Table II for the same computational condition as in Fig. 3 (c). The evaluation of T 2 equation, Λ 2 equation, and 2RDM all scale as N 6 for the TD-OCCD method. On the other hand, the evaluation of T 2 equation scales as N 5 , and we do not need a separate solution for Λ 2 , as it is just the complex conjugate of T 2 . The greatest time saving for the TD-OMP2 method comes from the evaluation of 2RDM since it scales as N 4 , and does not involve any operator products. Overall, TD-OMP2 achieves a significant cost reduction compared to TD-OCCD, making it an attractive choice for simulations involving larger chemical systems.

IV. CONCLUSION
We have successfully implemented the TD-OMP2 method for the real-time simulations of laser-induced dynamics in relatively large chemical systems. The TD-OMP2 method retains the size-extensivity and gaugeinvariance of TD-OCC, and is computationally much more efficient than the full TD-OCCD method. As a first numerical test, we have applied the method to the ground state of BH and the laser-driven dynamics of Ar. The imaginary time relaxation for BH obtains the identical ground-state energies with those by the stationary theory, which indicates the correctness of the implementation. The performance of the present method is tested against highly correlated TD-CASSCF and TD-OCCD methods for the case of laser-driven Ar. The results suggest a decent performance with consistent overestimation of the electron correlation in such highly nonlinear phenomena. Remarkably, the TD-OMP2 method is stable and does not breakdown even in the presence of strong laser-electron interaction, thanks to the nonperturbative inclusion of external fields and time-dependent orbital optimization. Further assessment of the TD-OMP2 method for different systems and severer simulation conditions (e.g, with higher-intensity and/or longerwavelength lasers) will be reported elsewhere.