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 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 of Ar. The calculation results are benchmarked against ab initio time-dependent complete-active-space self-consistent field, time-dependent optimized coupled-cluster double, and time-dependent Hartree–Fock methods, as well as a single active electron model to explore the role of electron correlation.

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–10 Atoms and molecules interacting with laser pulses of intensity ≳1014 W/cm2 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 extreme-ultraviolet (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 correlation12–15 are exactly described by the time-dependent Schrödinger equation (TDSE). However, direct numerical integration of TDSE is not feasible for systems with more than two electrons.16–25 As a result, single-active-electron (SAE) approximations have been widely used26,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 strong-field phenomena, the electron correlation is missing in this approximation.

Therefore, various tractable ab initio methods have been actively developed for the theoretical description of correlated multielectron dynamics in intense laser fields. Among the most reliable approaches to serve the purpose are the multiconfiguration time-dependent Hartree–Fock (MCTDHF) method28–32 and time-dependent complete-active-space self-consistent-field (TD-CASSCF) method.33,34 In MCTDHF, the electronic wavefunction is expressed in terms of full configuration interaction (FCI) expansion, Ψ(t) = ∑ICI(tI(t), where both CI coefficients {CI(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. Although 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-active-space self-consistent field (TD-RASSCF)35–37 and time-dependent 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 real-space implementations have turned out to be of great utility34,37,39,40 (see Ref. 12 for a review of ab initio wavefunction-based methods for multielectron dynamics, Ref. 41 and references therein for extension to correlated electron-nuclear dynamics, and Ref. 42 for a perspective on multiconfiguration approaches for indistinguishable particles). Nevertheless, truncated-CI-based methods, even with time-dependent orbitals, have a general drawback of not being size extensive.

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) method43 using optimized orthonormal orbitals where both orbitals and amplitudes are time-dependent. Our TD-OCC method is the time-dependent formulation of the stationary orbital optimized coupled-cluster 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 good agreement with the fully correlated TD-CASSCF method within the same active orbital space. In earlier work, Kvaal46 reported an orbital adaptive coupled-cluster (OATDCC) method built upon the work of Arponen using biorthogonal orbitals,47 although applications to laser-driven dynamics have not been addressed. The polynomial scaling TD-OCC43 or the OATDCC46 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-OCEPA048 based on the simplest version of the coupled-electron pair approximation,49–54 or equivalently the linearized CCD approximation,55 popular in quantum chemistry. The computational cost of the TD-OCEPA0 method scales as N6, 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 the 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, second-order many-body perturbation theory (MP2) can be seen as an approximation to the coupled-cluster double (CCD) method, having a lower N5 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 so-called 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 Λ functional-based derivation has an advantage that it can be easily 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(N5). As a first test case, we apply the TD-OMP2 method to electron dynamics in an Ar atom irradiated by a strong laser field and compare the results with those computed at SAE, 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.

Let us consider the electronic Hamiltonian H of the following form:

H=iNeh(ri,pi)+i<jNe1|rirj|,
(1)
h(r,p)=h0(r,p)+Vext(r,p,t),
(2)

where Ne is the number of electrons, ri and pi are the position and canonical momentum, respectively, of electron i, h0 is the field-free one-electron Hamiltonian, and Vext is the laser–electron interaction. The Hamiltonian in the second-quantization notation can be written as

Ĥ=ĥ+v^,
(3)
ĥ=hνμÊνμ,v^=12uνλμγÊνλμγ=14vνλμγÊνλμγ,
(4)

where Êνμ=ĉμĉν, Êνλμγ=ĉμĉγĉλĉν, and ĉμ(ĉμ) is the creation (annihilation) operator for the set of orthonormal 2nbas spin-orbitals {ψμ}, with nbas 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

hνμ=dx1ψμ*(x1)h(r1,p1)ψν(x1),
(5)
uνλμγ=dx1dx2ψμ*(x1)ψγ*(x2)ψν(x1)ψλ(x2)|r1r2|,
(6)

and vνλμγ=uνλμγuλνμγ, where xi = (ri, σi) is a spatial-spin coordinate. The complete set of 2nbas spin-orbitals (labeled with μ, ν, γ, λ) is divided into noccoccupied (o, p, q, r, s) and 2nbasnoccvirtual spin-orbitals having nonzero and vanishing occupations, respectively, in the total wavefunction. The occupied spin-orbitals are classified into ncorecore spin-orbitals, which are occupied in the reference Φ and kept uncorrelated, and nact = noccncoreactive spin-orbitals (t, u, v, w) among which the Nact = Nencore active electrons are correlated. The active spin-orbitals are further split into those in the hole space (i, j, k, l) and the particle space (a, b, c, d), which are defined as those occupied and unoccupied, respectively, in the reference Φ. The core spin-orbitals can be further split into the frozen-core space (i″, j″) and the dynamical-core space (i′, j′). The frozen-core orbitals are fixed in time, whereas the dynamical core orbitals are propagated in time33 (see Fig. 1 of Ref. 43 for a pictorial illustration of the orbital subspacing).

The stationary MP2 method can be considered as an approximation of the CCD method, where the full CCD energy functional

E=Φ|(1+Λ^2)eT^2ĤeT^2|Φ
(7)

is approximated by retaining the terms giving the second-order correction to the reference energy E0 = ⟨Φ|Ĥ|Φ⟩. Here, T^2=τijabÊijab and Λ^2=λabijÊabij, with τijab and λabij being the excitation and de-excitation amplitudes, respectively. Therefore, we start with the time-dependent CCD Lagrangian of the following form:

L(t)=Φ|(1+Λ^2)eT^2(Ĥit)eT^2|Φ,
(8)

which is a natural time-dependent extension of the energy functional [Eq. (7)]. Following Ref. 43, we consider the real-valued action functional

S=Rt0t1L(t)dt=12t0t1L(t)+L*(t)dt
(9)

and make it stationary, δS = 0, with respect to the variation of amplitudes δτijab, δλabij and variations of orthonormality-conserving orbitals δψν. The equations of motion (EOMs) for the amplitudes are obtained as

iτ̇ijab=Φijab|eT^2(ĤiX^)eT^2|Φ,
(10)
iλ̇abij=Φ|(1+Λ^2)eT^2(ĤiX^)eT^2|Φijab,
(11)

where X^=Xνμĉμĉν with Xνμ=ψμ|ψ̇ν and those for orbitals as

i|ψṗ=(1P^)F^p|ψp+|ψqXpq,
(12)
i(δbaρijρbaδij)Xjb=FjaρijρbaFbi*,
(13)

where P^=p|ψpψp|, Fqp=ϕp|F^q|ϕq, with

F^p|ψp=ĥ|ψp+Ŵsr|ψqρorqs(ρ1)po,
(14)
Wsr(x1)=dx2ψr*(x2)ψs(x2)|r1r2|,
(15)

where ρqp and ρqspr are the one- and two-body reduced density matrices, respectively, defined as

ρpq=Φ|(1+Λ^2)eT^2ÊqpeT^2|Φ,
(16)
ρprqs=Φ|(1+Λ^2)eT^2ÊqspreT^2|Φ.
(17)

Now, we derive the TD-OMP2 method as an approximation to the TD-OCCD method based on the partitioning of the electronic Hamiltonian

Ĥ=Ĥ(0)+Ĥ(1)
(18a)

into the zeroth-order part Ĥ(0)=f^=fqpÊqp and the perturbation Ĥ(1) = ĤĤ(0), with

fqp=hqp+vqjpj=(h0)qp+vqjpj+(Vext)qp,
(18b)

where (h0)νμ and (Vexe)νμ are the matrix elements of ĥ0 and V^ext, respectively.

Based on this partitioning, we apply the Baker–Campbell–Hausdorff expansion to the 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

L=L0iλabijτ̇ijab+Φ|Λ^2(ĤiX^)|Φ+Φ|[ĤiX^,T^2]|Φ+Φ|Λ^2[f^iX^,T^2]|Φ,
(19)

where L0=Φ|(ĤiX^)|Φ is the reference contribution. Inserting this TD-OMP2 Lagrangian into Eq. (9) and making it stationary with respect to amplitude variations derives TD-OMP2 amplitude equations,

iτ̇ijab=vijabp(ij)f¯jkτikab+p(ab)f¯caτijcb,
(20)
iλ̇abij=vabijp(ij)f¯kiλabkj+p(ab)f¯acλcbij,
(21)
f¯qp=fqpiXqp,
(22)

where p(μν) is the cyclic permutation operator. Importantly, Eqs. (20) and (21) reveal that the EOM for λabij is just the complex conjugate of that for τijab, resulting in λabij=τijab*. 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

ρpq=(ρ0)pq+γpq,
(23)
ρprqs=(ρ0)prqs+Γprqs,
(24)

where (ρ0)pq=δjqδpj and (ρ0)prqs=γpqδjsδrj+γrsδjqδpjγrqδjsδpjγpsδjqδrj+δjqδpjδksδrkδjsδpjδkqδrk are the reference contributions, and non-zero elements of γpq and Γprqs are given by

γji=12λcbkiτkjcb,γab=12λcaklτklcb,
(25)
Γijab=τijab,Γabij=λabij.
(26)

In summary, the TD-OMP2 method is defined by the EOMs of double excitation amplitudes τijab [Eq. (20)], the relation λabij=τijab*, 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|ψḃ=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 Vext, and (ii) in case with Vext, the zeroth-order Hamiltonian is time dependent through both the change of orbitals and an explicit time dependence of Vext(t), the latter implying nonperturbative inclusion of the laser–electron interaction.

Another more perturbation theoretic derivation begins with the following Lagrangian:

L=L0iλabijτ̇ijab+Φ(1+T^2)(f^NiX^)(1+T^2)+T^2v^N+v^NT^2Φc,
(27)

where the subscript c implies restriction to connected terms and f^N and v^N are the normal-ordered part of f^ and v^, respectively. Then, we follow the procedure of the time-dependent variational principle based on the action of Eq. (9) to obtain the identical method, as derived in Sec. II B. 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, {τijab(t),ψp(t)}. The Lagrangian of Eq. (27) can be viewed as the time-dependent extension of the Hylleraas energy functional used in the conventional, time-independent OMP2 method.58,59

To assess the performance of the method described in Sec. II, we do a series of ground-state energy calculations taking the BH molecule as an example with double-ζ plus polarization (DZP) basis.63 The imaginary time relaxation method43 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 MP2 and OMP2 are reported in Table I. The required matrix elements are obtained from the Gaussian09 program62 and interfaced with our numerical code. All the values are compared with the values from the PSI4 program.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 is taken to be the same as the number of basis nbas to make a comparison with PSI4,61 and all the orbitals are treated as active. However, our implementation allows us to optimize orbitals with Nnbas, in general, with a flexible classification of the occupied orbital space into frozen core, dynamical core, and active.

TABLE I.

Comparison of ground state energy (in hartree) of the BH molecule in DZP basis with the PSI461 program. The Gaussian09 program62 has been used to generate the required one-electron, two-electron, and overlap integrals required for the imaginary time propagation of EOMs in the orthonormalized Gaussian basis. A convergence cutoff of 10−15 hartree of energy difference has been chosen in subsequent time steps.

MP2OMP2
Bond length (bohr)This workPSI461 This workPSI461 
1.8 25.149 288 192 25.149 288 192 25.149 565 428 25.149 565 428 
2.0 25.183 068 430 25.183 068 429 25.183 367 319 25.183 367 319 
2.2 25.196 855 310 25.196 855 310 25.197 186 611 25.197 186 611 
2.4 (re25.198 570 797 25.198 570 797 25.198 947 432 25.198 947 432 
2.8 25.183 573 786 25.183 573 786 25.184 093 435 25.184 093 435 
3.2 25.159 043 339 25.159 043 339 25.159 809 546 25.159 809 546 
3.6 25.133 018 128 25.133 018 128 25.134 185 728 25.134 185 728 
4.0 25.108 605 403 25.108 605 402 25.110 381 925 25.110 381 924 
5.0 25.059 981 388 25.059 981 388 25.064 548 176 25.064 548 176 
6.0 25.029 750 598 25.029 750 598 25.039 814 328 25.039 814 327 
7.0 25.016 553 779 25.016 553 780 25.037 103 689 25.037 103 689 
MP2OMP2
Bond length (bohr)This workPSI461 This workPSI461 
1.8 25.149 288 192 25.149 288 192 25.149 565 428 25.149 565 428 
2.0 25.183 068 430 25.183 068 429 25.183 367 319 25.183 367 319 
2.2 25.196 855 310 25.196 855 310 25.197 186 611 25.197 186 611 
2.4 (re25.198 570 797 25.198 570 797 25.198 947 432 25.198 947 432 
2.8 25.183 573 786 25.183 573 786 25.184 093 435 25.184 093 435 
3.2 25.159 043 339 25.159 043 339 25.159 809 546 25.159 809 546 
3.6 25.133 018 128 25.133 018 128 25.134 185 728 25.134 185 728 
4.0 25.108 605 403 25.108 605 402 25.110 381 925 25.110 381 924 
5.0 25.059 981 388 25.059 981 388 25.064 548 176 25.064 548 176 
6.0 25.029 750 598 25.029 750 598 25.039 814 328 25.039 814 327 
7.0 25.016 553 779 25.016 553 780 25.037 103 689 25.037 103 689 

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 long-wavelength 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,

h(r,p)=12|p|2Z|r|+A(t)pz,
(28)

where Z(=18) is the atomic number and A(t) = −∫tE(t′)dt′ is 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:

E(t)=E0sin(ω0t)sin2πt3T
(29)

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 I0=E02. We have considered three different intensities of 1 × 1014 W/cm2, 2 × 1014 W/cm2, and 4 × 1014 W/cm2 for Ar.

The spherical finite-element discrete variable representation (FEDVR) basis34,40 is used in our implementation. The convergence with respect to the maximum angular momentum lmax is checked at the TDHF level, and lmax 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 < Rmax = 300. A cos14 mask function is switched on at 240 to avoid reflection from the boundary. The fourth-order exponential Runge–Kutta integrator64 is used to propagate equations of motions with 20 000 time steps per optical cycle. The simulations are continued after the end of pulse for further 6000 time steps.

For ab initio TDHF, TD-OMP2, TD-OCCD, and TD-CASSCF methods, the 1s2s2p core is kept frozen, and the dynamics of remaining 8 electrons are actively taken into account with 4 (TDHF) or 13 (TD-OMP2, TD-OCCD, and TD-CASSCF) active orbitals. The SAE method first diagonalizes the following effective Hamiltonian:65,66

heff=12|p|2+Veff(|r|)
(30)

on the FEDVR basis, where the effective potential Veff(r) is taken to be65,66

Veff(r)=1r1+Aer+(Z1A)eBr,
(31)

with A = 5.4 and B = 3.682 for Ar.66 This potential correctly supports 1s, 2s, 2p, 3s, and 3p bound orbitals, with the 3p orbital energy of −15.82 eV with the present FEDVR basis. After obtaining the ground state, we solve the effective one-electron Schrödinger equation

iddt|χ=Q^heff+A(t)pz|χ,
(32)

starting from the 3p0 orbital. The projector Q^=1j|ϕjϕj|, with ϕj running over 1s, 2s, 2p, 3s, and 3p± orbitals [multiplied by the gauge factor eiA(t)z], keeps χ(t) orthonormal to the inner shell orbitals.

In Fig. 1, we plot the time evolution of the sign-flipped dipole moment ⟨z⟩ evaluated as a trace z=ψp|z^|ψqρpq for ab initio methods and ⟨z⟩ = ⟨χ|z|χ⟩ for SAE. We compare the TD-OMP2 results with SAE, TDHF, TD-OCCD, and TD-CASSCF ones. Within the same active space, TD-CASSCF produces highly accurate results, useful for performance analysis of the TD-OMP2 method.

FIG. 1.

Time evolution of the sign-flipped dipole moment ⟨z⟩ of Ar irradiated by a laser pulse with a wavelength of 1200 nm and a peak intensity of 1 × 1014 W/cm2 (a), 2 × 1014 W/cm2 (b), and 4 × 1014 W/cm2 (c) calculated with TDHF, TD-OMP2, TD-OCCD, TD-CASSCF, and SAE methods.

FIG. 1.

Time evolution of the sign-flipped dipole moment ⟨z⟩ of Ar irradiated by a laser pulse with a wavelength of 1200 nm and a peak intensity of 1 × 1014 W/cm2 (a), 2 × 1014 W/cm2 (b), and 4 × 1014 W/cm2 (c) calculated with TDHF, TD-OMP2, TD-OCCD, TD-CASSCF, and SAE methods.

Close modal

The lowest (1.0 × 1014 W/cm2) and highest (4.0 × 1014 W/cm2) intensity cases characterize the dynamics with small and substantial amount of ionization, respectively (Fig. 2). The SAE approximation is known26,27,65 to work better for the latter case, where the dynamics is dominated by tunneling the ionization of the single, most weakly bound electron (one of the 3p0 electrons in the present case). It is not well suited for describing the former case dominated by collective bound dynamics. On the other hand, TDHF serves as the reference multielectron method without the (Coulomb) correlation by definition and provides a qualitatively correct description of the bound dynamics. It, however, fails to accurately describe cases with sizable tunneling ionization (see, e.g., Ref. 67 and references therein). These trends of SAE and TDHF methods are well confirmed in the performance comparison to the reference TD-CASSCF method in Figs. 1 and 2.

FIG. 2.

Time evolution of single ionization probability of Ar irradiated by a laser pulse with a wavelength of 1200 nm and a peak intensity of 1 × 1014 W/cm2 (a), 2 × 1014 W/cm2 (b), and 4 × 1014 W/cm2 (c) calculated with TDHF, TD-OMP2, TD-OCCD, TD-CASSCF, and SAE methods.

FIG. 2.

Time evolution of single ionization probability of Ar irradiated by a laser pulse with a wavelength of 1200 nm and a peak intensity of 1 × 1014 W/cm2 (a), 2 × 1014 W/cm2 (b), and 4 × 1014 W/cm2 (c) calculated with TDHF, TD-OMP2, TD-OCCD, TD-CASSCF, and SAE methods.

Close modal

As shown in Fig. 1(a), at the lowest intensity of 1014 W/cm2, all the methods except for SAE 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 an increase in intensity, the difference among the methods becomes more prominent. While all the methods except for SAE 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). In general, the performance of the TD-OMP2 method is better than TDHF and SAE due to consideration at least a part of the electron correlation. It is noticed that the TD-OMP2 dipole moment agrees better with the TD-CASSCF one for the highest intensity [Fig. 1(c)] than for the intermediate one [Fig. 1(b)], which might indicate that the latter case, with both sizable ionization and nontrivial correlation effects coexisting, is theoretically more challenging.

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 a systematic 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 and SAE, except for the highest intensity case, where the SAE result is as accurate as the TD-OCCD one.

In Figs. 3–5, 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 ab initio methods reproduce the HHG spectra relatively well, including an experimentally observed characteristic dip around the 52nd order (∼54 eV) at 2 × 1014 W/cm2 and 4 × 1014 W/cm2 related to the Cooper minimum of the photoionization spectrum68 at the same energy. However, TDHF systematically underestimates and fails to reproduce fine details. The SAE method severely underestimates the HHG yields, although the overall shape of the HHG spectrum is well reproduced, and an intensity-dependent scaling brings the spectral intensity at the higher plateau close to that of TD-CASSCF, especially for the highest intensity. The agreement with the TD-CASSCF results is much better for the TD-OCCD method, which contains nonlinear terms in the amplitude equations, which is then followed by TD-OMP2 with slight overestimation. This trend is consistent with the capability of each method to treat the electron correlation.

FIG. 3.

The HHG spectra from Ar irradiated by a laser pulse with a wavelength of 1200 nm and a peak intensity of 1 × 1014 W/cm2. Comparison of the results of SAE (a), TDHF (b), TD-OMP2 (c), and TD-OCCD (d) methods with that of TD-CASSCF.

FIG. 3.

The HHG spectra from Ar irradiated by a laser pulse with a wavelength of 1200 nm and a peak intensity of 1 × 1014 W/cm2. Comparison of the results of SAE (a), TDHF (b), TD-OMP2 (c), and TD-OCCD (d) methods with that of TD-CASSCF.

Close modal
FIG. 4.

The HHG spectra from Ar irradiated by a laser pulse with a wavelength of 1200 nm and a peak intensity of 2 × 1014 W/cm2. Comparison of the results of SAE (a), TDHF (b), TD-OMP2 (c), and TD-OCCD (d) methods with that of TD-CASSCF.

FIG. 4.

The HHG spectra from Ar irradiated by a laser pulse with a wavelength of 1200 nm and a peak intensity of 2 × 1014 W/cm2. Comparison of the results of SAE (a), TDHF (b), TD-OMP2 (c), and TD-OCCD (d) methods with that of TD-CASSCF.

Close modal
FIG. 5.

The HHG spectra from Ar irradiated by a laser pulse with a wavelength of 1200 nm and a peak intensity of 4 × 1014 W/cm2. Comparison of the results of SAE (a), TDHF (b), TD-OMP2 (c), and TD-OCCD (d) methods with that of TD-CASSCF.

FIG. 5.

The HHG spectra from Ar irradiated by a laser pulse with a wavelength of 1200 nm and a peak intensity of 4 × 1014 W/cm2. Comparison of the results of SAE (a), TDHF (b), TD-OMP2 (c), and TD-OCCD (d) methods with that of TD-CASSCF.

Close modal

In order to investigate the performance of the TD-OMP2 method for shorter wavelengths, where electrons not only in the highest-occupied but also in the inner-shell orbitals are driven by the laser, we consider a shorter wavelength of 200 nm with an intensity of 4 × 1014 W/cm2. All the other simulation parameters are identical to those given above. The obtained time evolution of the dipole moment (Fig. 6) shows a slight overestimation of the oscillation amplitude by TD-OMP2 and underestimation by TDHF compared to the TD-CASSCF result as in the case for the longer wavelength. It is encouraging, however, that the TD-OMP2 result is much closer to the TD-CASSCF one for the shorter wavelength, which is more sensitive to the treatment of correlation effects.

FIG. 6.

Time evolution of the sign-flipped dipole moment ⟨z⟩ of Ar irradiated by a laser pulse with a wavelength of 200 nm and a peak intensity of 4 × 1014 W/cm2, calculated with TDHF, TD-OMP2, TD-OCCD, and TD-CASSCF methods.

FIG. 6.

Time evolution of the sign-flipped dipole moment ⟨z⟩ of Ar irradiated by a laser pulse with a wavelength of 200 nm and a peak intensity of 4 × 1014 W/cm2, calculated with TDHF, TD-OMP2, TD-OCCD, and TD-CASSCF methods.

Close modal

It is worth noting that, for an intensity as high as 4 × 1014 W/cm2, the TD-OMP2 simulation completed stably. This should be the direct consequence of full inclusion of the laser–electron interaction in the zeroth-order 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 is compared in Table II for the same computational condition as in Fig. 3(c). The evaluation of the T2 equation, Λ2 equation, and 2RDM all scales as N6 for the TD-OCCD method. On the other hand, for the TD-OMP2 method, the evaluation of the T2 equation scales as N5, and we do not need a separate solution for Λ2, as it is just the complex conjugate of T2. The greatest time saving for the TD-OMP2 method comes from the evaluation of 2RDM since it scales as N4 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.

TABLE II.

Comparison of the CPU timea (in seconds) spent for the evaluation of the T2 equation, Λ2 equation, and 2RDM for TD-OCCD and TD-OMP2 methods.

TD-OCCDTD-OMP2
T2Λ22RDMT2Λ22RDM
40.8 55.5 109.5 1.11 … 0.25 
TD-OCCDTD-OMP2
T2Λ22RDMT2Λ22RDM
40.8 55.5 109.5 1.11 … 0.25 
a

CPU time spent for the simulation of the Ar atom for 1000 time steps (0 ≤ t ≤ 0.05T) of a real-time simulation (I0 = 4 × 1014 W/cm2 and λ = 1200 nm) using an Intel(R) Xeon(R) CPU with 12 processors having a clock speed of 3.33 GHz.

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 gauge-invariance 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 numerically investigated in comparison to SAE, TDHF, TD-OCCD, and TD-CASSCF methods for the case of laser-driven Ar. The results suggest a decent performance with a consistent overestimation of the correlation effect 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 longer-wavelength lasers) will be reported elsewhere.

The data that support the findings of this study are available from the corresponding author upon reasonable request.

This research was supported, in part, by a Grant-in-Aid for Scientific Research (Grant Nos. 17K05070, 18H03891, and 19H00869) from the Ministry of Education, Culture, Sports, Science and Technology (MEXT) of Japan. This research was also partially supported by the JST COI (Grant No. JPMJCE1313), the JST CREST (Grant No. JPMJCR15N1), and the MEXT Quantum Leap Flagship Program (MEXT Q-LEAP; Grant No. JPMXS0118067246).

1.
P. B.
Corkum
and
F.
Krausz
,
Nat. Phys.
3
,
381
(
2007
).
2.
F.
Krausz
and
M.
Ivanov
,
Rev. Mod. Phys.
81
,
163
(
2009
).
3.
J.
Itatani
et al.,
Nature
432
,
867
(
2004
).
4.
E.
Goulielmakis
et al.,
Nature
466
,
739
(
2010
).
5.
G.
Sansone
et al.,
Nature
465
,
763
(
2010
).
6.
M.
Schultze
et al.,
Science
328
,
1658
(
2010
).
7.
K.
Klünder
 et al,
Phys. Rev. Lett.
106
,
143002
(
2011
).
8.
L.
Belshaw
et al.,
J. Phys. Chem. Lett.
3
,
3751
(
2012
).
9.
F.
Calegari
et al.,
Science
346
,
336
(
2014
).
10.
O.
Smirnova
et al.,
Nature
460
,
972
(
2009
).
11.
P.
Antoine
,
A.
L’huillier
, and
M.
Lewenstein
,
Phys. Rev. Lett.
77
,
1234
(
1996
).
12.
K. L.
Ishikawa
and
T.
Sato
,
IEEE J. Sel. Top. Quantum Electron.
21
,
1
(
2015
).
13.
I.
Tikhomirov
,
T.
Sato
, and
K. L.
Ishikawa
,
Phys. Rev. Lett.
118
,
203202
(
2017
).
14.
Y.
Li
,
T.
Sato
, and
K. L.
Ishikawa
,
Phys. Rev. A
99
,
043401
(
2019
).
15.
S.
Pabst
and
R.
Santra
,
Phys. Rev. Lett.
111
,
233005
(
2013
).
16.
J. S.
Parker
,
E. S.
Smyth
, and
K. T.
Taylor
,
J. Phys. B: At., Mol. Opt. Phys.
31
,
L571
(
1998
).
17.
J. S.
Parker
et al.,
J. Phys. B: At., Mol. Opt. Phys.
33
,
L239
(
2000
).
18.
M.
Pindzola
and
F.
Robicheaux
,
Phys. Rev. A
57
,
318
(
1998
).
19.
S.
Laulan
and
H.
Bachau
,
Phys. Rev. A
68
,
013409
(
2003
).
20.
K. L.
Ishikawa
and
K.
Midorikawa
,
Phys. Rev. A
72
,
013407
(
2005
).
21.
J.
Feist
 et al,
Phys. Rev. Lett.
103
,
063002
(
2009
).
22.
K. L.
Ishikawa
and
K.
Ueda
,
Phys. Rev. Lett.
108
,
033003
(
2012
).
23.
S.
Sukiasyan
,
K. L.
Ishikawa
, and
M.
Ivanov
,
Phys. Rev. A
86
,
033423
(
2012
).
24.
W.
Vanroose
,
D. A.
Horner
,
F.
Martin
,
T. N.
Rescigno
, and
C. W.
McCurdy
,
Phys. Rev. A
74
,
052702
(
2006
).
25.
D. A.
Horner
 et al,
Phys. Rev. Lett.
101
,
183002
(
2008
).
27.
K. C.
Kulander
,
Phys. Rev. A
36
,
2726
(
1987
).
28.
J.
Caillat
 et al,
Phys. Rev. A
71
,
012712
(
2005
).
29.
T.
Kato
and
H.
Kono
,
Chem. Phys. Lett.
392
,
533
(
2004
).
30.
M.
Nest
,
T.
Klamroth
, and
P.
Saalfrank
,
J. Chem. Phys.
122
,
124102
(
2005
).
31.
D. J.
Haxton
,
K. V.
Lawler
, and
C. W.
McCurdy
,
Phys. Rev. A
83
,
063416
(
2011
).
32.
D.
Hochstuhl
and
M.
Bonitz
,
J. Chem. Phys.
134
,
084106
(
2011
).
33.
T.
Sato
and
K. L.
Ishikawa
,
Phys. Rev. A
88
,
023402
(
2013
).
34.
T.
Sato
 et al,
Phys. Rev. A
94
,
023405
(
2016
).
35.
H.
Miyagi
and
L. B.
Madsen
,
Phys. Rev. A
87
,
062511
(
2013
).
36.
H.
Miyagi
and
L. B.
Madsen
,
Phys. Rev. A
89
,
063416
(
2014
).
37.
D. J.
Haxton
and
C. W.
McCurdy
,
Phys. Rev. A
91
,
012509
(
2015
).
38.
T.
Sato
and
K. L.
Ishikawa
,
Phys. Rev. A
91
,
023417
(
2015
).
39.
J. J.
Omiste
,
W.
Li
, and
L. B.
Madsen
,
Phys. Rev. A
95
,
053422
(
2017
).
40.
Y.
Orimo
,
T.
Sato
,
A.
Scrinzi
, and
K. L.
Ishikawa
,
Phys. Rev. A
97
,
023423
(
2018
).
41.
R.
Anzaki
,
T.
Sato
, and
K. L.
Ishikawa
,
Phys. Chem. Chem. Phys.
19
,
22008
(
2017
).
42.
A. U. J.
Lode
,
C.
Lévêque
,
L. B.
Madsen
,
A. I.
Streltsov
, and
O. E.
Alon
,
Rev. Mod. Phys.
92
,
011001
(
2020
).
43.
T.
Sato
,
H.
Pathak
,
Y.
Orimo
, and
K. L.
Ishikawa
,
J. Chem. Phys.
148
,
051101
(
2018
).
44.
C. D.
Sherrill
,
A. I.
Krylov
,
E. F.
Byrd
, and
M.
Head-Gordon
,
J. Chem. Phys.
109
,
4171
(
1998
).
45.
A. I.
Krylov
,
C. D.
Sherrill
,
E. F.
Byrd
, and
M.
Head-Gordon
,
J. Chem. Phys.
109
,
10669
(
1998
).
46.
S.
Kvaal
,
J. Chem. Phys.
136
,
194109
(
2012
).
48.
H.
Pathak
,
T.
Sato
, and
K. L.
Ishikawa
,
J. Chem. Phys.
152
,
124115
(
2020
); arXiv:2001.02206.
49.
W.
Meyer
,
Int. J. Quantum Chem.
5
,
341
(
1971
).
50.
R.
Ahlrichs
,
P.
Scharf
, and
C.
Ehrhardt
,
J. Chem. Phys.
82
,
890
(
1985
).
51.
F.
Wennmohs
and
F.
Neese
,
Chem. Phys.
343
,
217
(
2008
).
52.
F.
Neese
,
F.
Wennmohs
, and
A.
Hansen
,
J. Chem. Phys.
130
,
114108
(
2009
).
53.
C.
Kollmar
and
F.
Neese
,
Mol. Phys.
108
,
2449
(
2010
).
54.
J.-P.
Malrieu
,
H.
Zhang
, and
J.
Ma
,
Chem. Phys. Lett.
493
,
179
(
2010
).
55.
J.
Čížek
,
J. Chem. Phys.
45
,
4256
(
1966
).
56.
R. J.
Bartlett
,
Annu. Rev. Phys. Chem.
32
,
359
(
1981
).
57.
U.
Bozkaya
,
J. M.
Turney
,
Y.
Yamaguchi
,
H. F.
Schaefer
 III
, and
C. D.
Sherrill
,
J. Chem. Phys.
135
,
104103
(
2011
).
58.
L.
Adamowicz
and
R. J.
Bartlett
,
J. Chem. Phys.
86
,
6314
(
1987
).
59.
F.
Neese
,
T.
Schwabe
,
S.
Kossmann
,
B.
Schirmer
, and
S.
Grimme
,
J. Chem. Theory Comput.
5
,
3060
(
2009
).
60.
U.
Bozkaya
,
J. Chem. Phys.
135
,
224103
(
2011
).
61.
J. M.
Turney
 et al,
Wiley Interdiscip. Rev.: Comput. Mol. Sci.
2
,
556
(
2012
).
62.
M.
Frisch
 et al, Gaussian 09, revision d. 01,
2009
.
63.
R.
Harrison
and
N.
Handy
,
Chem. Phys. Lett.
95
,
386
(
1983
).
64.
M.
Hochbruck
and
A.
Ostermann
,
Acta Numer.
19
,
209
(
2010
).
65.
H. G.
Muller
and
F. C.
Kooiman
,
Phys. Rev. Lett.
81
,
1207
(
1998
).
66.
K.
Schiessl
,
E.
Persson
,
A.
Scrinzi
, and
J.
Burgdörfer
,
Phys. Rev. A
74
,
053412
(
2006
).
67.
T.
Sato
and
K. L.
Ishikawa
,
J. Phys. B: At., Mol. Opt. Phys.
47
,
204031
(
2014
).
68.
H. J.
Wörner
,
H.
Niikura
,
J. B.
Bertrand
,
P. B.
Corkum
, and
D. M.
Villeneuve
,
Phys. Rev. Lett.
102
,
103901
(
2009
).