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*(*N*^{7}), 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.

## I. INTRODUCTION

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 ≳10^{14} W/cm^{2} 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) method^{22–26} and time-dependent complete-active-space self-consistent-field (TD-CASSCF) method^{27,28} are the best known for the purpose. Both these methods are based on the full configuration interaction (FCI) expansion of the wavefunction, Ψ(*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 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 methods^{29–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 utility^{34} 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) method^{42} 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 cluster^{44–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*(*N*^{6}), 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*(*N*^{5}), 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-OCEPA0^{56} and TD-OMP2^{57} 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*(*N*^{8}) 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 work^{38,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*(*N*^{7}), 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.

## II. METHOD

### A. Background

The time-dependent electronic Hamiltonian is given by

where *N*_{e} denotes the number of electrons, *r*_{i} and *p*_{i} 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

where $c\u0302\mu \u2020$ ($c\u0302\mu $) represents a creation (annihilation) operator in a complete, orthonormal set of 2*n*_{bas} spin orbitals {*ψ*_{μ}(*t*)}, which are, in general, time-dependent, and Eq. (3) defines one-electron ($E\u0302\nu \mu $) and two-electron ($E\u0302\nu \lambda \mu \gamma $) generators. *n*_{bas} 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

where *x*_{i} = (*r*_{i}, *σ*_{i}) represents a composite spatial-spin coordinate.

The complete set of 2*n*_{bas} spin orbitals (labeled *μ*, *ν*, *γ*, *λ*) is divided into *n*_{occ} *occupied* (*o*, *p*, *q*, *r*, *s*) and 2*n*_{bas} − *n*_{occ} *virtual* 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 *n*_{core} *core* spin orbitals, which are occupied in the reference Φ and kept uncorrelated, and *N* = *n*_{occ} − *n*_{core} *active* 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 time^{27} (see Fig. 1 of Ref. 38 for a pictorial illustration of the orbital subspacing).

### B. Review of the TD-OCC method

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

where

where $\tau ij\u2026ab\u2026$ and $\lambda ab\u2026ij\u2026$ are excitation and de-excitation amplitudes, respectively. The stationary conditions, *δS* = 0, with respect to the variation of amplitudes $\delta \tau ij\u2026ab\u2026$, $\delta \lambda ab\u2026ij\u2026$, 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,

where $X\u0302=X\nu \mu E\u0302\nu \mu $, $X\nu \mu =\u27e8\psi \mu |\psi \u0307\nu \u27e9$ is anti-Hermitian, $L0=\u27e8\Phi |(H\u0302\u2212iX\u0302)|\Phi \u27e9$ 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) $\rho pq$ and $\rho prqs$ are defined, respectively, by

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

where the reference contributions $(\rho 0)pq=\delta jq\delta pj$ and $(\rho 0)prqs=\gamma pq\delta js\delta rj+\gamma rs\delta jq\delta pj\u2212\gamma rq\delta js\delta pj\u2212\gamma ps\delta jq\delta rj+\delta jq\delta pj\delta ks\delta rk\u2212\delta js\delta pj\delta kq\delta rk$ (*j*, *k* running over core and hole spaces) are independent of the correlation treatment, and the correlation contributions are defined as

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

### C. TD-OCCDT(4) method

The full TD-OCCDT Lagrangian is given by^{38}

where we truncate after $\Lambda \u0302=\Lambda \u03022+\Lambda \u03023$ and $T\u0302=T\u03022+T\u03023$, expand $eT\u03022+T\u03023$, and retain only the connected diagrams. For deriving the TD-OCCDT(4) method, we make a further approximation to *L*_{CCDT} by retaining those terms contributing up to the fourth order in the many-body perturbation theory,

where $f\u0304=f\u0302\u2212iX\u0302$, $f\u0302=(hqp+vqjpj){E\u0302qp}$, $v\u0302=vqspr{E\u0302qspr}/4$, and $vqspr=uqspr\u2212usqpr$.

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

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:

where $1\u0302=\u2211\mu |\psi \mu \u3009\u3008\psi \mu |$ is the identity operator within the space spanned by the given basis, $P\u0302=\u2211q|\psi q\u3009\u3008\psi q|$ is the projector onto the occupied spin-orbital space, and

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

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 by^{28}

where ** E** is the external electric field.

## III. NUMERICAL RESULTS AND DISCUSSION

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

within the dipole approximation in the velocity gauge, where *Z* is the atomic number and *A*(*t*) = −*∫*^{t}*E*(*t*′)d*t*′ 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:

for 0 ≤ *t* ≤ 3*T*, 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 × 10^{14}, 8 × 10^{14}, and 1 × 10^{15} W/cm^{2} for Ne and 1 × 10^{14}, 2 × 10^{14}, and 4 × 10^{14} W/cm^{2} for Ar.

Our implementation uses a spherical finite-element discrete variable representation (FEDVR) basis^{28,65} for representing orbital functions,

where *Y*_{lm} and *f*_{k}(*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 *l*_{max}, and the radial-FEDVR basis supports the range of radial coordinates 0 ≤ *r* ≤ *R*_{max}, with the cos^{1/4} mask function used as an absorbing boundary for avoiding unphysical reflection from the wall of the simulation box. We have used *l*_{max} = 47 for Ne and *l*_{max} = 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 method^{67} 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 1*s* orbital of Ne and 1*s*2*s*2*p* 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 $\u27e8\psi p|z\u0302|\psi q\u27e9\rho 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.

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 10^{14} W/cm^{2} 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 correlation^{26,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.

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)].

Let us now turn to high-harmonic generation. We present the HHG spectra calculated at three different laser intensities for Ne in Figs. 7–9 and for Ar in Figs. 10–12. 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,

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).

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.

Method . | Time (min) . | Cost reduction (%) . |
---|---|---|

TD-CASSCF | 7638 | ⋯ |

TD-OCCDT | 5734 | 25 |

TD-OCCDT(4) | 4600 | 40 |

TD-OCCD | 3904 | 49 |

Method . | Time (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.3*T*) of a real-time simulation (*I*_{0} = 4 × 10^{14} W/cm^{2} and *λ* = 800 nm), using an Intel(R) Xeon(R) Gold 6230 CPU with ten processors having a clock speed of 2.10 GHz.

## IV. CONCLUDING REMARKS

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*(*N*^{7}). 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.

## ACKNOWLEDGMENTS

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.

## DATA AVAILABLITY

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

### APPENDIX: ALGEBRAIC DETAILS OF THE REDUCED DENSITY MATRICES AND THE MATRIX *B* FOR THE TD-OCCDT(4) METHOD

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

The matrix *B* is given by

where $Fqp=\u27e8\psi p|F\u0302|\psi q\u27e9$, with the operator $F\u0302$ defined by Eq. (23).