We present a cost-effective treatment of the triple excitation amplitudes in the time-dependent optimized coupled-cluster (TD-OCC) framework called TD-OCCDT(4) for studying intense laser-driven multielectron dynamics. It considers triple excitation amplitudes correct up to the fourth-order in many-body perturbation theory and achieves a computational scaling of O(N7), with N being the number of active orbital functions. This method is applied to the electron dynamics in Ne and Ar atoms exposed to an intense near-infrared laser pulse with various intensities. We benchmark our results against the TD complete-active-space self-consistent field (TD-CASSCF), TD-OCC with double and triple excitations (TD-OCCDT), TD-OCC with double excitations (TD-OCCD), and TD Hartree–Fock (TDHF) methods to understand how this approximate scheme performs in describing nonperturbatively nonlinear phenomena, such as field-induced ionization and high-harmonic generation. We find that the TD-OCCDT(4) method performs equally well as the TD-OCCDT method, almost perfectly reproducing the results of the fully correlated TD-CASSCF with a more favorable computational scaling.

The advance in experimental techniques for high-intensity ultrashort optical pulses has contributed substantially to enrich strong-field physics and attosecond science.1–4 It also shows promise to measure and manipulate the motions of electrons in a quantum many-body system.5,6 One of the key processes in attosecond physics is high-order harmonic generation (HHG).7,8 It is a highly nonlinear frequency conversion process that delivers ultrashort coherent light pulses in the extreme-ultraviolet (XUV) to the soft x-ray regions interacting with laser pulses of intensity ≳1014 W/cm2 in the visible to the mid-infrared spectral range. The HHG spectra consist of a plateau where the intensity of the emitted radiation remains nearly constant up to many orders, followed by a sharp cutoff, beyond which harmonics are ceased to observe.9 

The time-dependent Schrödinger equation (TDSE) provides the most rigorous theoretical description of all these dynamical phenomena.10–17 However, it remains a major challenge to achieve the exact numerical solution of the TDSE for more than two-electron systems due to its high dimensionality. Thus, single-active electron (SAE) approximations have been widely used,18,19 in which only the outermost electron is explicitly treated. Though useful in obtaining a qualitative insight into different strong-field phenomena, an SAE is intrinsically unable to treat the electron correlation in the field-induced multielectron dynamics.

Therefore, the development of various time-dependent (TD), electron-correlation methods and the refinement of the existing ones is ongoing for a better understanding of intense laser-driven multielectron dynamics. See Ref. 20 for a review on various wavefunction-based methods for the study of laser-driven electron dynamics and Ref. 21 for multiconfiguration approaches for indistinguishable particles. The multiconfiguration time-dependent Hartree–Fock (MCTDHF) method22–26 and time-dependent complete-active-space self-consistent-field (TD-CASSCF) method27,28 are the best known for the purpose. Both these methods are based on the full configuration interaction (FCI) expansion of the wavefunction, Ψ(t) = ∑ICI(tI(t), where both CI coefficients {CI(t)} and occupied spin orbitals {ψp(t)} constituting Slater determinants {ΦI(t)} are time-dependent and propagated in time according to the time-dependent variational principle (TDVP). The TD-CASSCF method broadens the applicability by flexibly classifying occupied orbital space into frozen-core (do not participate in dynamics), dynamical-core (do participate in the dynamics remaining doubly occupied), and active (fully correlated description of the active electrons) subspaces. By such flexible classification, it helps get rid of a large number of chemically or physically inert electrons from the correlation calculation and, thereby, enhance the applicability. However, these methods scale factorially with the number of correlated electrons, restricting large scale applications. By limiting the configuration interaction (CI) expansion of the wavefunction up to an affordable level, various less demanding methods29–33 such as the time-dependent restricted-active-space self-consistent field (TD-RASSCF)29,30 and time-dependent occupation-restricted multiple-active-space (TD-ORMAS)32 have been developed. These methods have proven to be of great utility34 even though not size-extensive.

On the other hand, the coupled-cluster method is the most celebrated many-body method while dealing with the electron correlation. It is size-extensive at any given level of truncation, scales polynomially, and has proven its effectiveness on a wide range of problems in the stationary theory.35–37 Therefore, we have developed its time-dependent extension using optimized orthonormal orbitals, called the time-dependent optimized coupled-cluster (TD-OCC) method.38 The TD-OCC method is the time-dependent formulation of the stationary orbital optimized coupled-cluster method,39–41 and we have considered double and triple excitation amplitudes (TD-OCCDT) in our implementation. Here again, both orbitals and amplitudes are time-dependent and propagated in time according to the TDVP. We applied our method to simulate HHG spectra and strong-field ionization in Ar. The TD-OCCDT produced virtually identical results with the TD-CASSCF method in a polynomial cost-scaling within the same chosen active orbital subspace.38 In an earlier development, Kvaal reported an orbital adaptive time-dependent coupled-cluster (OATDCC) method42 built upon the work of Arponen using biorthogonal orbitals.43 However, this work does not address applications to laser-driven dynamics. We also take note of the very initial attempts of the time-dependent coupled cluster44–46 and its applications using time-independent orbitals.47–52 See Ref. 53 for the discussion on numerical stability of time-dependent coupled-cluster methods with time-independent orbitals and Refs. 54 and 55 for the symplectic integrator and interpretation of the time-dependent coupled-cluster method. The nonlinear, exponential parameterization combined with the time-dependent orbital optimization in TD-OCC has been found essential for the accurate description of high-field processes.38,56–58

To further reduce the computational scaling and, thereby, enhance the applicability to larger chemical systems, we have developed time-dependent optimized coupled electron pair approximation (TD-OCEPA0).56 This method is much more efficient than time-dependent optimized coupled-cluster with double (TD-OCCD) excitations, even though both scale as O(N6), with N being the number of active orbitals. The superior efficiency is due to the linear Hermitian structure of the Lagrangian, which helps avoid a separate solution for the de-excitation amplitudes. The above argument is also true for the time-dependent optimized second-order many-body perturbation (TD-OMP2) methods, which scale as O(N5), reported in Ref. 57. We have also reported the implementation and numerical assessments of the second-order approximation to the time-dependent coupled-cluster single double (TD-CCSD) method, called the TD-CC2 method.58 The stationary variant of this method (CC2)59 is one of the widely accepted lower-cost methods. However, the TD-CC2 method suffers from the lack of gauge invariance and instability in the real-time propagation under the high-intensity fields, both due to imperfect optimization of orbitals as opposed to the case in fully optimized TD-OMP2.57 

Thus, the TD-OCCD method and its lower-cost approximations such as TD-OCEPA056 and TD-OMP257 will be useful to investigate multielectron dynamics in large chemical systems driven by intense laser fields. Nevertheless, the inclusion of triple excitation amplitudes is certainly attractive to achieve decisive accuracy, as demonstrated by quantitative agreement of TD-OCCDT and TD-CASSCF results for strong-field ionization from Ar atom.38 Therefore, it is desirable to reduce the computational cost of a method including the effect of triple excitations. The full TD-OCCDT method scales as O(N8) as in the stationary CCSDT method.63 Numerous reports are available in the stationary coupled-cluster literature for reducing the computational scaling, retaining a part of the triple excitation amplitudes.60–64 The coupled-cluster theory and the many-body perturbation theory are intricately connected. Therefore, approximation of the coupled-cluster effective Hamiltonian based on many-body perturbation theory is an appealing one.37 

In this article, we report an efficient approximation to the TD-OCCDT method, treating triple excitation amplitudes within the flexibly chosen active space correct up to the fourth order in the many-body perturbation theory. We designate this method as TD-OCCDT(4). It is closely related to the stationary CCSDT1 method of Bartlett and co-workers.60,61 We do not include single excitation amplitudes following our previous work38,56–58 but optimize the orbitals according to the TDVP, which is ideally suited for studying strong-field dynamics involving substantial excitation and ionization. The TD-OCCDT(4) method inherits gauge invariance from the parent TD-OCCDT method due to the variational optimization of orbitals. This method scales as O(N7), which is an order of magnitude lower than the full TD-OCCDT method, rendering this method very much affordable while retaining important contributions of the triples.

We apply our newly developed TD-OCCDT(4) scheme to Ne and Ar atoms exposed to a strong near-infrared laser field. We compare TD-OCCDT(4) results with those of the TD-CASSCF, TD-OCCDT, TD-OCCD, and uncorrelated TDHF methods to assess the numerical performance in describing laser-driven multielectron dynamics. We find that the TD-OCCDT(4) method performs equally well as the TD-OCCDT and fully correlated TD-CASSCF methods while reducing the computational cost.

This paper is organized as follows: We describe the TD-OCCDT(4) method in Sec. II. Section III presents its numerical application to Ne and Ar. Conclusions are given in Sec. IV. Hartree atomic units are used unless otherwise stated, and Einstein convention is implied throughout for summation over orbital indices.

The time-dependent electronic Hamiltonian is given by

H=i=1Neh(ri,pi,t)+i=1Ne1j=i+1Ne1|rirj|,
(1)

where Ne denotes the number of electrons, ri and pi denote the position and canonical momentum, respectively, of an electron i, and h denotes the one-electron Hamiltonian operator including the laser–electron interaction.

In the second quantization representation, it reads

Ĥ=hνμĉμĉν+12uνλμγĉμĉγĉλĉν
(2)
=hνμÊνμ+12uνλμγÊνλμγ,
(3)

where ĉμ (ĉμ) represents a creation (annihilation) operator in a complete, orthonormal set of 2nbas spin orbitals {ψμ(t)}, which are, in general, time-dependent, and Eq. (3) defines one-electron (Êνμ) and two-electron (Êνλμγ) generators. nbas is the number of basis functions used for expanding the spatial part of ψμ, which, in the present real-space implementation, corresponds to the number of grid points (see Sec. III), and

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

where xi = (ri, σi) represents a composite spatial-spin coordinate.

The complete set of 2nbas spin orbitals (labeled μ, ν, γ, λ) is divided into noccoccupied (o, p, q, r, s) and 2nbasnoccvirtual spin orbitals. The coupled-cluster (or CI) wavefunction is constructed only with occupied spin orbitals, which are time-dependent in general, and virtual spin orbitals form the orthogonal complement of the occupied spin-orbital space. The occupied spin orbitals are classified into ncorecore spin orbitals, which are occupied in the reference Φ and kept uncorrelated, and N = noccncoreactive spin orbitals (t, u, v, w) among which the 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 further be split into frozen-core space (i″, j″) and dynamical-core space (i′, j′). The frozen-core orbitals are fixed in time, whereas dynamical-core orbitals are propagated in time27 (see Fig. 1 of Ref. 38 for a pictorial illustration of the orbital subspacing).

Following our previous work,38 we rely on the real action formulation of the time-dependent variational principle with orthonormal orbitals,

S=Ret0t1Ldt=12t0t1L+L*dt,
(6)
L=Φ|(1+Λ̂)eT̂ĤiteT̂|Φ,
(7)

where

T̂=T̂2+T̂3=τijabÊijab+τijkabcÊijkabc,
(8)
Λ̂=Λ̂2+Λ̂3=λabijÊabij+λabcijkÊabcijk,
(9)

where τijab and λabij are excitation and de-excitation amplitudes, respectively. The stationary conditions, δS = 0, with respect to the variation of amplitudes δτijab, δλabij, and orthonormality-conserving orbital variations δψμ, gives us the corresponding equations of motions (EOMs).

The coupled-cluster Lagrangian can be written down into two equivalent forms,

L=L0+Φ|(1+Λ̂)[{ĤiX̂}eT̂]c|Φiλabijτ̇ijab
(10a)
=(hqpiXqp)ρpq+12uqsprρprqsiλabijτ̇ijab,
(10b)

where X̂=XνμÊνμ, Xνμ=ψμ|ψ̇ν is anti-Hermitian, L0=Φ|(ĤiX̂)|Φ is the reference contribution, and [⋯]c indicates that only the contributions from connected coupled-cluster diagrams are retained. The one-electron and two-electron reduced density matrices (RDMs) ρpq and ρprqs are defined, respectively, by

ρpq=Φ|(1+Λ̂)eT̂ÊqpeT̂|Φ,
(11)
ρprqs=Φ|(1+Λ̂)eT̂ÊqspreT̂|Φ.
(12)

The one-electron and two-electron RDMs are separated into reference and correlation contributions,

ρpq=(ρ0)pq+γpq,
(13)
ρprqs=(ρ0)prqs+γprqs,
(14)

where the reference contributions (ρ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 (j, k running over core and hole spaces) are independent of the correlation treatment, and the correlation contributions are defined as

γpq=Φ|(1+Λ̂)[{Êqp}eT̂]c|Φ,
(15a)
γprqs=Φ|(1+Λ̂)[{Êqspr}eT̂]c|Φ,
(15b)

where the bracket {⋯} implies that the operator inside is normal ordered relative to the reference.

The full TD-OCCDT Lagrangian is given by38 

LCCDT=L0+Φ|(1+Λ̂2+Λ̂3)×{ĤiX̂}1+T̂2+12!T̂22+T̂3+T̂2T̂3c|Φiλabijτ̇ijabiλabcijkτ̇ijkabc,
(16)

where we truncate after Λ̂=Λ̂2+Λ̂3 and T̂=T̂2+T̂3, expand eT̂2+T̂3, and retain only the connected diagrams. For deriving the TD-OCCDT(4) method, we make a further approximation to LCCDT by retaining those terms contributing up to the fourth order in the many-body perturbation theory,

LCCDT(4)=L0+Φ|(1+Λ̂2)[(f̄+v̂)eT̂2]c|Φ+Φ|Λ̂2[(f̄+v̂)T̂3]c|Φ+Φ|Λ̂3(f̄T̂3)c|Φ+Φ|Λ̂3(v̂T̂2)c|Φiλabijτ̇ijabiλabcijkτ̇ijkabc,
(17)

where f̄=f̂iX̂, f̂=(hqp+vqjpj){Êqp}, v̂=vqspr{Êqspr}/4, and vqspr=uqsprusqpr.

Using LCCDT(4) of Eq. (17), the TD-OCCDT(4) amplitude equations are derived from the stationary conditions δS/δλabij(t)=0, δS/δλabcijk(t)=0, δS/δτijab(t)=0, and δS/δτijkabc(t)=0 as

iτ̇ijab=vijabp(ij)f̄jkτikab+p(ab)f̄caτijcb+12vcdabτijcd+12vijklτklab+p(ij)p(ab)vicakτkjcb12p(ij)τikabτjlcdvcdkl+12p(ab)τijbcτkladvcdkl+14τklabτijcdvcdkl+12p(ij)p(ab)τilbcτjkadvcdkl+f̄ckτijkabc+12p(ab)vcdkaτkijcdb12p(ij)vciklτkljcab,
(18)
iτ̇ijkabc=p(k/ij)p(a/bc)vdkbctijadp(i/jk)p(c/ab)vjklctilabp(k/ij)f̄klτijlabc+p(c/ab)f̄dcτijkabd,
(19)
iλ̇abij=vabijp(ij)f̄kiλabkj+p(ab)f̄acλcbij+12vabcdλcdij+12vklijλabkl+p(ij)p(ab)vkbcjλacik12p(ij)λcdikτklcdvabjl+12p(ab)λbcklτklcdvadij+14λabklτklcdvcdij+12p(ij)p(ab)λacjkτklcdvbdil12p(ij)λabikτklcdvcdjl+12p(ab)λbcijτklcdvadkl+14λcdijτklcdvabkl+12p(ab)vbkdcλadcijk12p(ij)vlkjcλabcilk,
(20)
iλ̇abcijk=p(k/ij)p(a/bc)vbcdkλadijp(c/ab)p(i/jk)vlcjkλabij+p(c/ab)f̄cdλabdijkp(k/ij)f̄lkλabcijl+p(i/jk)p(a/bc)f̄aiλbcjk,
(21)

where p(μν) and p(μ|νγ) are the permutation operators, p(μν)Aμν = AμνAνμ, and p(μ/νγ) = 1 − p(μν) − p(μγ).

Following the discussion in Ref. 28, the EOM for the orbitals can be written down in the following form:

i|ψṗ=(1̂P̂)F̂|ψp+i|ψqXpq,
(22)

where 1̂=μ|ψμψμ| is the identity operator within the space spanned by the given basis, P̂=q|ψqψq| is the projector onto the occupied spin-orbital space, and

F̂|ψp=ĥ|ψp+Ŵsr|ψqPorqs(D1)po,
(23)
Dqp=12(ρqp+ρpq*),
(24)
Pqspr=12(ρqspr+ρprqs*),
(25)
Wsr(x1)=dx2ψr*(x2)ψs(x2)|r1r2|.
(26)

The matrix element Xpq describes orbital rotations among various orbital subspaces. Intra-space orbital rotations ({Xji}, {Xji}, and {Xba}) are redundant and can be arbitrary anti-Hermitian matrix elements. On the other hand, inter-space rotations (Xia and Xai=Xia*) are non-redundant and determined by solving the linear equation38 

iδbaDijDbaδijXjb=Bia.
(27)

The algebraic expressions of RDMs and the matrix B are given in  Appendix A. The orbital rotations involving frozen-core orbitals within the electric dipole approximation are given by28 

iXμi=0(lengthgauge)E(t)ψi|r|ψμ(velocitygauge),
(28)

where E is the external electric field.

In this section, we apply the TD-OCCDT(4) method to the laser-driven electron dynamics in Ne and Ar. The field dependent one-electron Hamiltonian is given by

h(r,p)=p22Z|r|+A(t)pz,
(29)

within the dipole approximation in the velocity gauge, where Z is the atomic number and A(t) = −tE(t′)dt′ is the vector potential, with E(t) being the laser electric field linearly polarized along the z axis. Our methods are gauge invariant; both length-gauge and velocity-gauge simulations produce identical results for physical observables upon numerical convergence. The velocity gauge is known to be advantageous in simulating high-field phenomena.28,65,66

We assume a laser electric field of the following form:

E(t)=E0sin(ω0t)sin2πt3T,
(30)

for 0 ≤ t ≤ 3T, and E(t) = 0 otherwise, with a central wavelength of λ = 2π/ω0 = 800 nm and a period of T = 2π/ω0 ∼ 2.67 fs. We consider three different intensities 5 × 1014, 8 × 1014, and 1 × 1015 W/cm2 for Ne and 1 × 1014, 2 × 1014, and 4 × 1014 W/cm2 for Ar.

Our implementation uses a spherical finite-element discrete variable representation (FEDVR) basis28,65 for representing orbital functions,

χklm(r,θ,ψ)=1rfk(r)Ylm(θ,ϕ),
(31)

where Ylm and fk(r) are spherical harmonics and the normalized radial-FEDVR basis function, respectively. The expansion of the spherical harmonics continued up to the maximum angular momentum lmax, and the radial-FEDVR basis supports the range of radial coordinates 0 ≤ rRmax, with the cos1/4 mask function used as an absorbing boundary for avoiding unphysical reflection from the wall of the simulation box. We have used lmax = 47 for Ne and lmax = 63 for Ar, and the FEDVR basis supporting the radial coordinates 0 < r < 260 using 68 finite elements each containing 21 (for Ne) and 23 (for Ar) DVR functions. The absorbing boundary switched on at r = 180 in all our simulations. First, the ground state is obtained by the imaginary time relaxation. Starting from the ground state, we perform real-time simulations by switching the laser field. The fourth-order exponential Runge–Kutta method67 is used to propagate the EOMs with 10 000 time steps for each optical cycle. The simulations are continued for further 3000 time steps after the end of the pulse. The 1s orbital of Ne and 1s2s2p orbitals of Ar are kept frozen at the canonical Hartree–Fock orbitals, and for the correlated methods, the eight valence electrons are correlated among thirteen active orbitals.

In Figs. 1 and 2, we report the time evolution of the dipole moment, evaluated as a trace ψp|ẑ|ψqρpq using the 1RDM. In Figs. 1(a), 1(c), and 1(e), we observe that, whereas the TDHF results deviate substantially from the TD-CASSCF results with larger deviation for higher intensity, the TD-OCCDT produces almost identical results to those of TD-CASSCF irrespective of the laser intensity. The valence electrons are driven farther with increasing laser intensity, which renders the dynamical electron correlation (the effect beyond the TDHF description) more relevant.

FIG. 1.

Time evolution of the dipole moment of Ne irradiated by a laser pulse with a wavelength of 800 nm and a peak intensity of 5×1014 W/cm2 [(a) and (b)], 8×1014 W/cm2 [(c) and (d)], and 1×1015 W/cm2 [(e) and (f)], calculated with TDHF, TD-OCCD, TD-OCCDT(4), TD-OCCDT, and TD-CASSCF methods.

FIG. 1.

Time evolution of the dipole moment of Ne irradiated by a laser pulse with a wavelength of 800 nm and a peak intensity of 5×1014 W/cm2 [(a) and (b)], 8×1014 W/cm2 [(c) and (d)], and 1×1015 W/cm2 [(e) and (f)], calculated with TDHF, TD-OCCD, TD-OCCDT(4), TD-OCCDT, and TD-CASSCF methods.

Close modal
FIG. 2.

Time evolution of the dipole moment of Ar irradiated by a laser pulse with a wavelength of 800 nm and a peak intensity of 1×1014 W/cm2 [(a) and (b)], 2×1014 W/cm2 [(c) and (d)], and 4×1014 W/cm2 [(e) and (f)], calculated with TDHF, TD-OCCD, TD-OCCDT(4), TD-OCCDT, and TD-CASSCF methods.

FIG. 2.

Time evolution of the dipole moment of Ar irradiated by a laser pulse with a wavelength of 800 nm and a peak intensity of 1×1014 W/cm2 [(a) and (b)], 2×1014 W/cm2 [(c) and (d)], and 4×1014 W/cm2 [(e) and (f)], calculated with TDHF, TD-OCCD, TD-OCCDT(4), TD-OCCDT, and TD-CASSCF methods.

Close modal

Figures 1(b), 1(d), and 1(f) compare coupled-cluster results with three different approximation schemes, TD-OCCDT, TD-OCCDT(4), and TD-OCCD. We observe that the TD-OCCDT and TD-OCCDT(4) results are indistinguishable on the scale of Figs. 1(b), 1(d), and 1(f), while TD-OCCD slightly underestimates the oscillatory behavior. The TD-OCCDT(4) method takes into account the dominant part of the electron correlation, successfully approximating the triple excitation amplitudes.

Figure 2 shows similar results for Ar. The TDHF, TD-OCCDT, and TD-CASSCF methods give virtually the same results for 1014 W/cm2 intensity [Fig. 2(a)]. On the other hand, the TDHF result deviates from the others for the higher intensities [Figs. 2(c) and 2(e)], indicating that the electron correlation plays an increasingly important role as the electrons are driven farther from the nucleus. In Figs. 2(b), 2(d), and 2(f), we observe, analogous to the Ne case, that the TD-OCCDT(4) result excellently agrees with the TD-OCCDT one, while the TD-OCCD result slightly deviates from them.

Overall, for both Ne (Fig. 1) and Ar (Fig. 2), the TD-OCCDT(4) and TD-OCCDT methods with polynomial scalings deliver practically the same results with the TD-CASSCF method with a factorial scaling, irrespective of the employed intensities. The TD-OCCDT(4) scales by one order lower than the TD-OCCDT but works equally well with the latter. Thus, the TD-OCCDT(4) method is the most advantageous of the three.

In Figs. 3 and 4, we present the temporal evolution of single ionization of Ne and Ar, respectively, evaluated as the probability of finding an electron outside a radius of 20 a.u. For this purpose, we have used RDMs as defined in Refs. 68–70. Ionization probabilities are more sensitive to the electron correlation26,71 to assess the capabilities of newly implemented methods in comparison to the established ones. In Figs. 3(a), 3(c), 3(e), 4(a), 4(c), and 4(e), we compare the TD-CASSCF, TD-OCCDT, and TDHF methods, while we compare different time-dependent optimized coupled-cluster approaches in Figs. 3(b), 3(d), 3(f), 4(b), 4(d), and 4(f). For both Ne and Ar, the TDHF method leads to substantial underestimation, indicating an important role played by the electron correlation in ionization. On the other hand, the TD-OCCDT produces nearly identical results with the TD-CASSCF method. We see from Figs. 3(b), 3(d), 3(f), 4(b), 4(d), and 4(f) that the TD-OCCDT(4) and TD-OCCDT results are almost identical to each other, irrespective of the laser intensities, whereas the TD-OCCD results slightly underestimate. Especially for Ar [Figs. 4(b), 4(d), and 4(f)], the deviation grows with increasing laser intensity. This observation again suggests that the consideration of the triples becomes more important for higher intensity due to a greater role of the electron correlation.

FIG. 3.

Time evolution of single ionization probability of Ne irradiated by a laser pulse with a wavelength of 800 nm and a peak intensity of 5×1014 W/cm2 [(a) and (b)], 8×1014 W/cm2 [(c) and (d)], and 1×1015 W/cm2 [(e) and (f)], calculated with TDHF, TD-OCCD, TD-OCCDT(4), TD-OCCDT, and TD-CASSCF methods.

FIG. 3.

Time evolution of single ionization probability of Ne irradiated by a laser pulse with a wavelength of 800 nm and a peak intensity of 5×1014 W/cm2 [(a) and (b)], 8×1014 W/cm2 [(c) and (d)], and 1×1015 W/cm2 [(e) and (f)], calculated with TDHF, TD-OCCD, TD-OCCDT(4), TD-OCCDT, and TD-CASSCF methods.

Close modal
FIG. 4.

Time evolution of single ionization probability of Ar irradiated by a laser pulse with a wavelength of 800 nm and a peak intensity of 1×1014 W/cm2 [(a) and (b)], 2×1014 W/cm2 [(c) and (d)], and 4×1014 W/cm2 [(e) and (f)], calculated with TDHF, TD-OCCD, TD-OCCDT(4), TD-OCCDT and TD-CASSCF methods.

FIG. 4.

Time evolution of single ionization probability of Ar irradiated by a laser pulse with a wavelength of 800 nm and a peak intensity of 1×1014 W/cm2 [(a) and (b)], 2×1014 W/cm2 [(c) and (d)], and 4×1014 W/cm2 [(e) and (f)], calculated with TDHF, TD-OCCD, TD-OCCDT(4), TD-OCCDT and TD-CASSCF methods.

Close modal

Figures 5 and 6 display the temporal evolution of double ionization of Ne and Ar, respectively, evaluated as the probability of finding two electrons outside a radius of 20 a.u. Double ionization, which may involve recollision, is expected to be even more susceptible to electron correlation. For the case of Ne (Fig. 5), all the methods except for TDHF predict similar double ionization probability. However, in the case of Ar (Fig. 6) with more double ionization than that for Ne, the TD-OCCD method fails to reproduce the TD-OCCDT(4) and TD-OCCDT results, which, in turn, agrees well with the TD-CASSCF. The deviation is especially large at the highest laser intensity employed [Fig. 6(f)].

FIG. 5.

Time evolution of double ionization probability of Ne irradiated by a laser pulse with a wavelength of 800 nm and a peak intensity of 5×1014 W/cm2 [(a) and (b)], 8×1014 W/cm2 [(c) and (d)], and 1×1015 W/cm2 [(e) and (f)], calculated with TDHF, TD-OCCD, TD-OCCDT(4), TD-OCCDT, and TD-CASSCF methods.

FIG. 5.

Time evolution of double ionization probability of Ne irradiated by a laser pulse with a wavelength of 800 nm and a peak intensity of 5×1014 W/cm2 [(a) and (b)], 8×1014 W/cm2 [(c) and (d)], and 1×1015 W/cm2 [(e) and (f)], calculated with TDHF, TD-OCCD, TD-OCCDT(4), TD-OCCDT, and TD-CASSCF methods.

Close modal
FIG. 6.

Time evolution of double ionization probability of Ar irradiated by a laser pulse with a wavelength of 800 nm and a peak intensity of 1×1014 W/cm2 [(a) and (b)], 2×1014 W/cm2 [(c) and (d)], and 4×1014 W/cm2 [(e) and (f)], calculated with TDHF, TD-OCCD, TD-OCCDT(4), TD-OCCDT, and TD-CASSCF methods.

FIG. 6.

Time evolution of double ionization probability of Ar irradiated by a laser pulse with a wavelength of 800 nm and a peak intensity of 1×1014 W/cm2 [(a) and (b)], 2×1014 W/cm2 [(c) and (d)], and 4×1014 W/cm2 [(e) and (f)], calculated with TDHF, TD-OCCD, TD-OCCDT(4), TD-OCCDT, and TD-CASSCF methods.

Close modal

Let us now turn to high-harmonic generation. We present the HHG spectra calculated at three different laser intensities for Ne in Figs. 79 and for Ar in Figs. 1012. The HHG spectrum is obtained as the modulus squared I(ω) = |a(ω)|2 of the Fourier transform of the expectation value of the dipole acceleration, which, in turn, is evaluated with a modified Ehrenfest expression.28 We also plot the absolute relative deviation,

δ(ω)=a(ω)aTD-CASSCF(ω)aTD-CASSCF(ω),
(32)

of the spectral amplitude a(ω) from the TD-CASSCF value for each method in Figs. 7(d), 8(d), 9(d), 10(d), 11(d), and 12(d).

FIG. 7.

The HHG spectra [(a)–(c)] and the relative deviation (d) of the spectral amplitude from the TD-CASSCF spectrum from Ne irradiated by a laser pulse with a wavelength of 800 nm and a peak intensity of 5 × 1014 W/cm2 with various methods.

FIG. 7.

The HHG spectra [(a)–(c)] and the relative deviation (d) of the spectral amplitude from the TD-CASSCF spectrum from Ne irradiated by a laser pulse with a wavelength of 800 nm and a peak intensity of 5 × 1014 W/cm2 with various methods.

Close modal
FIG. 8.

The HHG spectra [(a)–(c)] and the relative deviation (d) of the spectral amplitude from the TD-CASSCF spectrum from Ne irradiated by a laser pulse with a wavelength of 800 nm and a peak intensity of 8 × 1014 W/cm2 with various methods.

FIG. 8.

The HHG spectra [(a)–(c)] and the relative deviation (d) of the spectral amplitude from the TD-CASSCF spectrum from Ne irradiated by a laser pulse with a wavelength of 800 nm and a peak intensity of 8 × 1014 W/cm2 with various methods.

Close modal
FIG. 9.

The HHG spectra [(a)–(c)] and the relative deviation (d) of the spectral amplitude from the TD-CASSCF spectrum from Ne irradiated by a laser pulse with a wavelength of 800 nm and a peak intensity of 1 × 1015 W/cm2 with various methods.

FIG. 9.

The HHG spectra [(a)–(c)] and the relative deviation (d) of the spectral amplitude from the TD-CASSCF spectrum from Ne irradiated by a laser pulse with a wavelength of 800 nm and a peak intensity of 1 × 1015 W/cm2 with various methods.

Close modal
FIG. 10.

The HHG spectra [(a)–(c)] and the relative deviation (d) of the spectral amplitude from the TD-CASSCF spectrum from Ar irradiated by a laser pulse with a wavelength of 800 nm and a peak intensity of 1 × 1014 W/cm2 with various methods.

FIG. 10.

The HHG spectra [(a)–(c)] and the relative deviation (d) of the spectral amplitude from the TD-CASSCF spectrum from Ar irradiated by a laser pulse with a wavelength of 800 nm and a peak intensity of 1 × 1014 W/cm2 with various methods.

Close modal
FIG. 11.

The HHG spectra [(a)–(c)] and the relative deviation (d) of the spectral amplitude from the TD-CASSCF spectrum from Ar irradiated by a laser pulse with a wavelength of 800 nm and a peak intensity of 2 × 1014 W/cm2 with various methods.

FIG. 11.

The HHG spectra [(a)–(c)] and the relative deviation (d) of the spectral amplitude from the TD-CASSCF spectrum from Ar irradiated by a laser pulse with a wavelength of 800 nm and a peak intensity of 2 × 1014 W/cm2 with various methods.

Close modal
FIG. 12.

The HHG spectra [(a)–(c)] and the relative deviation (d) of the spectral amplitude from the TD-CASSCF spectrum from Ar irradiated by a laser pulse with a wavelength of 800 nm and a peak intensity of 4 × 1014 W/cm2 with various methods.

FIG. 12.

The HHG spectra [(a)–(c)] and the relative deviation (d) of the spectral amplitude from the TD-CASSCF spectrum from Ar irradiated by a laser pulse with a wavelength of 800 nm and a peak intensity of 4 × 1014 W/cm2 with various methods.

Close modal

Again, TDHF underestimates the spectral intensity, presumably due to the underestimation of tunneling ionization (Figs. 3 and 4), i.e., the first step in the three-step model.72,73 The other four methods predict the virtually same spectra on the logarithmic scale of Figs. 3 and 4. The relative deviation of results for each method from the TD-CASSCF ones depends only weakly on the harmonic order and has a general trend of decreasing as TDHF > TD-OCCD > TD-OCCDT ≳ TD-OCCDT(4). This observation clearly demonstrates the usefulness of the TD-OCCDT(4) method.

Finally, we compare the computational cost of the different methods considered in this article. All our simulations have been run on an Intel(R) Xeon(R) Gold 6230 central processing unit (CPU) with ten processors with a clock speed of 2.10 GHz. Table I reports the total computation time spent for the simulation corresponding to Fig. 12. It also lists the relative cost reduction with respect to the TD-CASSCF method. The cost reduction for the TD-OCCDT(4) method (40%) is larger than that for the TD-OCCDT method (25%). As we have repeatedly seen above, the inclusion of triples is certainly important for highly accurate simulations of high-field phenomena, and the TD-OCCDT(4) method takes into account the essential part of the triples at an affordable reduced computational cost.

TABLE I.

Comparison of the total simulation timea (in min) spent for TD-CASSCF, TD-OCCDT, TDCCDT(4), and TD-OCCD methods.

MethodTime (min)Cost reduction (%)
TD-CASSCF  7638 ⋯ 
TD-OCCDT  5734  25 
TD-OCCDT(4)  4600  40 
TD-OCCD  3904  49 
MethodTime (min)Cost reduction (%)
TD-CASSCF  7638 ⋯ 
TD-OCCDT  5734  25 
TD-OCCDT(4)  4600  40 
TD-OCCD  3904  49 
a

Time spent for the simulation of Ar atom for 33 000 time steps (0 ≤ t ≤ 3.3T) of a real-time simulation (I0 = 4 × 1014 W/cm2 and λ = 800 nm), using an Intel(R) Xeon(R) Gold 6230 CPU with ten processors having a clock speed of 2.10 GHz.

We have developed the TD-OCCDT(4) method as a cost-effective and efficient approximation to the TD-OCCDT method. This method considers triple excitation amplitudes correct up to the fourth order in the many-body perturbation theory. Its computational cost scales as O(N7). As a numerical assessment, we have applied this method to Ne and Ar atoms irradiated by an intense near-infrared laser pulse and compared the calculated time-dependent dipole moment, single and double ionization probability, and HHG spectra with those obtained with other methods. We have found that the TD-OCCDT(4) method takes into account the major part of triple excitation amplitudes, which are important for obtaining consistently accurate results over a wide range of problems; the TD-OCCDT(4) method delivers the results virtually indistinguishable from those of its parent, full TD-OCCDT method. Furthermore, the TD-OCCDT(4) method achieves a 40% computational cost reduction in comparison to the TD-CASSCF method for the case of Ar calculation with eight active electrons with 13 active orbitals. The present TD-OCCDT(4) method will extend the applicability of highly accurate ab initio simulations of electron dynamics to larger chemical systems.

This research was supported, in part, by a Grant-in-Aid for Scientific Research (Grants Nos. JP18H03891 and JP19H00869) 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.

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

The algebraic expressions of non-vanishing correlation contributions to 1RDMs and 2RDMs are given by

γij=12λcdkjτkicd112λabckljτkliabc,
(A1a)
γab=12λcaklτklcb+112λcdaijkτijkcdb,
(A1b)
γia=14λbcjkτjkibca,
(A1c)
γabcd=12λabklτklcd,γijkl=12λcdklτijcd,
(A1d)
γciab=12λdcjkτjkidab,γjkia=12λbcliτljkbca,
(A1e)
γbjia=λcbkiτkjca,γbcia=12λdbcjkiτjkda,
(A1f)
γabij=λabij,γbjik12λcablikτljca,
(A1g)
γijab=τijab+12p(ij)p(ab)λcdklτkicaτjlbd12p(ij)λcdklτkicdτljab12p(ab)λcdklτklcaτijdb+14λcdklτijcdτklab.
(A1h)

The matrix B is given by

Bia=FpaDipDpaFpi*i8τ̇ijkabcλbcjki8τijkabcλ̇bcjk,
(A2)

where Fqp=ψp|F̂|ψq, with the operator F̂ defined by Eq. (23).

1.
J.
Itatani
 et al,
Nature
432
,
867
(
2004
).
2.
P. B.
Corkum
and
F.
Krausz
,
Nat. Phys.
3
,
381
(
2007
).
3.
F.
Krausz
and
M.
Ivanov
,
Rev. Mod. Phys.
81
,
163
(
2009
).
5.
E.
Goulielmakis
 et al,
Nature
466
,
739
(
2010
).
6.
G.
Sansone
 et al,
Nature
465
,
763
(
2010
).
7.
P.
Antoine
,
A.
L’Huillier
, and
M.
Lewenstein
,
Phys. Rev. Lett.
77
,
1234
(
1996
).
8.
F.
Calegari
 et al,
Science
346
,
336
(
2014
).
9.
K. L.
Ishikawa
, “
High-harmonic generation
,” in
Advances in Solid State Lasers Development and Applications
(
In Tech
,
2010
).
10.
J. S.
Parker
,
E. S.
Smyth
, and
K. T.
Taylor
,
J. Phys. B: At., Mol. Opt. Phys.
31
,
L571
(
1998
).
11.
J. S.
Parker
 et al,
J. Phys. B: At., Mol. Opt. Phys.
33
,
L239
(
2000
).
12.
K. L.
Ishikawa
and
K.
Midorikawa
,
Phys. Rev. A
72
,
013407
(
2005
).
13.
J.
Feist
 et al,
Phys. Rev. Lett.
103
,
063002
(
2009
).
14.
K. L.
Ishikawa
and
K.
Ueda
,
Phys. Rev. Lett.
108
,
033003
(
2012
).
15.
S.
Sukiasyan
,
K. L.
Ishikawa
, and
M.
Ivanov
,
Phys. Rev. A
86
,
033423
(
2012
).
16.
W.
Vanroose
,
D. A.
Horner
,
F.
Martin
,
T. N.
Rescigno
, and
C. W.
McCurdy
,
Phys. Rev. A
74
,
052702
(
2006
).
17.
D. A.
Horner
 et al,
Phys. Rev. Lett.
101
,
183002
(
2008
).
18.
19.
K. C.
Kulander
,
Phys. Rev. A
36
,
2726
(
1987
).
20.
K. L.
Ishikawa
and
T.
Sato
,
IEEE J. Sel. Top. Quantum Electron.
21
,
8700916
(
2015
).
21.
A. U. J.
Lode
,
C.
Lévêque
,
L. B.
Madsen
,
A. I.
Streltsov
, and
O. E.
Alon
,
Rev. Mod. Phys.
92
,
011001
(
2020
).
22.
J.
Caillat
 et al,
Phys. Rev. A.
71
,
012712
(
2005
).
23.
T.
Kato
and
H.
Kono
,
Chem. Phys. Lett.
392
,
533
(
2004
).
24.
M.
Nest
,
T.
Klamroth
, and
P.
Saalfrank
,
J. Chem. Phys.
122
,
124102
(
2005
).
25.
D. J.
Haxton
,
K. V.
Lawler
, and
C. W.
McCurdy
,
Phys. Rev. A
83
,
063416
(
2011
).
26.
D.
Hochstuhl
and
M.
Bonitz
,
J. Chem. Phys.
134
,
084106
(
2011
).
27.
T.
Sato
and
K. L.
Ishikawa
,
Phys. Rev. A
88
,
023402
(
2013
).
28.
T.
Sato
 et al,
Phys. Rev. A
94
,
023405
(
2016
).
29.
H.
Miyagi
and
L. B.
Madsen
,
Phys. Rev. A
87
,
062511
(
2013
).
30.
H.
Miyagi
and
L. B.
Madsen
,
Phys. Rev. A
89
,
063416
(
2014
).
31.
D. J.
Haxton
and
C. W.
McCurdy
,
Phys. Rev. A
91
,
012509
(
2015
).
32.
T.
Sato
and
K. L.
Ishikawa
,
Phys. Rev. A
91
,
023417
(
2015
).
33.
E.
Lötstedt
,
T.
Kato
, and
K.
Yamanouchi
,
Phys. Rev. A
97
,
013423
(
2018
).
34.
I. S.
Wahyutama
,
T.
Sato
, and
K. L.
Ishikawa
,
Phys. Rev. A
99
,
063420
(
2019
).
35.
T.
Helgaker
,
P.
Jorgensen
, and
J.
Olsen
,
Molecular Electronic-Structure Theory
(
John Wiley & Sons
,
2014
).
36.
H. G.
Kümmel
,
Int. J. Mod. Phys. B.
17
,
5311
(
2003
).
37.
I.
Shavitt
and
R. J.
Bartlett
,
Many-Body Methods in Chemistry and Physics: MBPT and Coupled-Cluster Theory
(
Cambridge University Press
,
2009
).
38.
T.
Sato
,
H.
Pathak
,
Y.
Orimo
, and
K. L.
Ishikawa
,
J. Chem. Phys.
148
,
051101
(
2018
).
39.
G. E.
Scuseria
and
H. F.
Schaefer
 III
,
Chem. Phys. Lett.
142
,
354
(
1987
).
40.
C. D.
Sherrill
,
A. I.
Krylov
,
E. F. C.
Byrd
, and
M.
Head-Gordon
,
J. Chem. Phys.
109
,
4171
(
1998
).
41.
A. I.
Krylov
,
C. D.
Sherrill
,
E. F. C.
Byrd
, and
M.
Head-Gordon
,
J. Chem. Phys.
109
,
10669
(
1998
).
42.
S.
Kvaal
,
J. Chem. Phys.
136
,
194109
(
2012
).
44.
K.
Schönhammer
,
Phys. Rev. B
18
,
6606
(
1978
).
45.
P.
Hoodbhoy
and
J. W.
Negele
,
Phys. Rev. C.
18
,
2380
(
1978
).
46.
P.
Hoodbhoy
and
J. W.
Negele
,
Phys. Rev. C.
19
,
1971
(
1979
).
47.
S.
Guha
and
D.
Mukherjee
,
Chem. Phys. Lett.
186
,
84
(
1991
).
48.
M. D.
Prasad
,
J. Chem. Phys.
88
,
7005
(
1988
).
49.
K. L.
Sebastian
,
Phys. Rev. B
31
,
6976
(
1985
).
50.
C.
Huber
and
T.
Klamroth
,
J. Chem. Phys.
134
,
054113
(
2011
).
51.
D. A.
Pigg
,
G.
Hagen
,
H.
Nam
, and
T.
Papenbrock
,
Phys. Rev. C.
86
,
014308
(
2012
).
52.
D. R.
Nascimento
and
A. E.
DePrince
 III
,
J. Chem. Theory Comput.
12
,
5834
(
2016
).
53.
H. E.
Kristiansen
,
Ø. S.
Schøyen
,
S.
Kvaal
, and
T. B.
Pedersen
,
J. Chem. Phys.
152
,
071102
(
2020
).
54.
T. B.
Pedersen
and
S.
Kvaal
,
J. Chem. Phys.
150
,
144106
(
2019
).
55.
T. B.
Pedersen
,
H. E.
Kristiansen
,
T.
Bodenstein
,
S.
Kvaal
, and
Ø. S.
Schøyen
,
J. Chem. Theory Comput.
17
,
388
(
2021
).
56.
H.
Pathak
,
T.
Sato
, and
K. L.
Ishikawa
,
J. Chem. Phys.
152
,
124115
(
2020
).
57.
H.
Pathak
,
T.
Sato
, and
K. L.
Ishikawa
,
J. Chem. Phys.
153
,
034110
(
2020
).
58.
H.
Pathak
,
T.
Sato
, and
K. L.
Ishikawa
,
Mol. Phys.
118
,
e1813910
(
2020
).
59.
O.
Christiansen
,
H.
Koch
, and
P.
Jørgensen
,
Chem. Phys. Lett.
243
,
409
(
1995
).
60.
Y. S.
Lee
,
S. A.
Kucharski
, and
R. J.
Bartlett
,
J. Chem. Phys.
81
,
5906
(
1984
).
61.
M.
Urban
,
J.
Noga
,
S. J.
Cole
, and
R. J.
Bartlett
,
J. Chem. Phys.
83
,
4041
(
1985
).
62.
J.
Noga
,
R. J.
Bartlett
, and
M.
Urban
,
Chem. Phys. Lett.
134
,
126
(
1987
).
63.
J.
Noga
and
R. J.
Bartlett
,
J. Chem. Phys.
86
,
7041
(
1987
).
64.
J. D.
Watts
,
J.
Gauss
, and
R. J.
Bartlett
,
J. Chem. Phys.
98
,
8718
(
1993
).
65.
Y.
Orimo
,
T.
Sato
,
A.
Scrinzi
, and
K. L.
Ishikawa
,
Phys. Rev. A
97
,
023423
(
2018
).
66.
Y.
Orimo
,
T.
Sato
, and
K. L.
Ishikawa
,
Phys. Rev. A
100
,
013419
(
2019
).
67.
M.
Hochbruck
and
A.
Ostermann
,
Acta Numer.
19
,
209
(
2010
).
68.
F.
Lackner
,
I.
Březinová
,
T.
Sato
,
K. L.
Ishikawa
, and
J.
Burgdörfer
,
Phys. Rev. A
91
,
023412
(
2015
).
69.
F.
Lackner
,
I.
Březinová
,
T.
Sato
,
K. L.
Ishikawa
, and
J.
Burgdörfer
,
Phys. Rev. A
95
,
033414
(
2017
).
70.
T.
Sato
,
Y.
Orimo
,
T.
Teramura
,
O.
Tugs
, and
K. L.
Ishikawa
, “
Time-dependent complete-active-space self-consistent-field method for ultrafast intense laser science
,” in
Progress in Ultrafast Intense Laser Science XIV
(
Springer
,
2018
), pp.
143
171
.
71.
T.
Sato
and
K. L.
Ishikawa
,
J. Phys. B: At., Mol. Opt. Phys.
47
,
204031
(
2014
).
72.
73.
K. C.
Kulander
,
K. J.
Schafer
, and
J. L.
Krause
, “
Dynamics of short-pulse excitation, ionization and harmonic conversion
,” in
Super-Intense Laser–Atom Physics
, NATO ASI Series B: Physics, Vol. 316, edited by
B.
Piraux
, et al
(
Springer
,
1993
), pp.
95
110
.