Time-dependent optimized coupled-cluster method for multielectron dynamics IV: Approximate consideration of the triple excitation amplitudes

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 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 time-dependent complete-active-space self-consistent field (TD-CASSCF), time-dependent optimized coupled-cluster with double and triple excitations (TD-OCCDT), time-dependent optimized coupled-cluster with double excitations (TD-OCCD), and the time-dependent 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 fully-correlated TD-CASSCF with a more favorable computational scaling.


I. INTRODUCTION
The advance in experimental techniques for highintensity ultrashort optical pulses has contributed substantially to enrich strong-field physics and attosecond science [1][2][3][4] . It is also showing promise to measure and manipulate the motions of electrons in a quantum manybody 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 extremeultraviolet (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][11][12][13][14][15][16][17] . However, it remains a major challenge to achieve the exact numerical solution of 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 usea) Electronic mail: pathak@atto.t.u-tokyo.ac.jp b) Electronic mail: sato@atto.t.u-tokyo.ac.jp c) Electronic mail: ishiken@n.t.u-tokyo.ac.jp ful in obtaining a qualitative insight into different strongfield phenomena, SAE is intrinsically unable to treat the electron correlation in the field-induced multielectron dynamics.
Therefore, the development of various time-dependent, 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 wavefunctionbased 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][23][24][25][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 spinorbitals {ψ 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 to 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, re-stricting large scale applications. By limiting the configuration interaction (CI) expansion of the wavefunction up to an affordable level, various less demanding methods [29][30][31][32][33] like the time-dependent restricted-activespace self-consistent field (TD-RASSCF) 29,30 and timedependent occupation-restricted multiple-active-space (TD-ORMAS) 32 have been developed. These methods have proven to be of great utility 34 even though not sizeextensive.
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][36][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][40][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 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 (OAT-DCC) 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][45][46] , and its applications using timeindependent orbitals [47][48][49][50][51][52] . See Ref. 53 for the discussion on numerical stability of time-dependent coupled-cluster methods with time-independent orbitals, and Refs. 54,55 for symplectic integrator and interpretation of the timedependent coupled-cluster method. The TD-OCCDT method scales as O(N 8 ), with N being the number of active orbitals.
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 excitations (TD-OCCD), even though both scale as O(N 6 ). The superior efficiency is due to the linear Hermitian structure of the Lagrangian, which helps to 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 scales 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 a 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. Numerous reports are available in the stationary coupled-cluster literature for reducing the computational scaling retaining a part of the triple excitation amplitudes [60][61][62][63][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 coworkers 60,61 . We do not include single excitation amplitudes following our previous work 38, [56][57][58] but optimize the orbitals according to 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 asses the numerical performance in describing laser-driven multielectron dynamics. We find that the TD-OCCDT(4) method performs equally as well as the TD-OCCDT and fully correlated TD-CASSCF methods, while reducing the computational cost.
This paper is organized as follows. We derscribe 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.

A. Background
The time-dependent electronic Hamiltonian is given by, where N e denotes the number of electrons, r i and p i the position and canonical momentum, respectively, of an electron i, and h the one-electron Hamiltonian operator including the laser-electron interaction.
In the second quantization representation, it readŝ whereĉ † µ (ĉ µ ) represents a creation (annihilation) operator in a complete, orthonormal set of 2n bas spinorbitals {ψ µ (t)}, which are, in general, time-dependent, and Eq. (3) defines one-electron (Ê µ ν ) and two-electron (Ê µγ νλ ) 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 2n bas spin-orbitals (labeled with µ, ν, γ, λ) is divided into n occ occupied (o, p, q, r, s) and 2n bas −n occ virtual spin-orbitals. The coupled-cluster (or CI) wavefunction is constructed only with occupied spinorbitals 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 the dynamicalcore space (i , j ). The frozen-core orbitals are fixed in time, whereas dynamical core orbitals are propagated in time. 27  Following our previous work 38 , we rely on the real action formulation of the time-dependent variational principle with orthonormal orbitals, whereT where τ ab··· ij··· and λ ij··· ab··· are excitation and deexcitation amplitudes, respectively. The stationary conditions, δS = 0, with respect to the variation of amplitudes δτ ab··· ij··· , δλ ij··· ab··· 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, whereX = X µ νÊ µ ν , and X µ ν = ψ µ |ψ ν is anti-Hermitian, L 0 = Φ|(Ĥ − 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) ρ q p and ρ qs pr are defined, respectively, by ρ qs pr = Φ|(1 +Λ)e −TÊpr qs eT |Φ .
The one-electron and two-electron RDMs are separated into reference and correlation contributions, where the reference contributions (ρ 0 ) q p = δ q j δ j p and 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.
For deriving the TD-OCCDT(4) method, we make a further approximation to L CCDT by retaining those terms contributing up to fourth order in the many-body perturbation theory, L CCDT of Eq. (17), the TD-OCCDT(4) amplitude equations are derived from the stationary conditions δS/δλ ij ab (t) = 0, δS/δλ ijk abc (t) = 0, δS/δτ ab ij (t) = 0, and δS/δτ abc ijk (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, where1 = µ |ψ µ ψ µ | 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 The matrix element X q p describes orbital rotations among various orbital subspaces. Intra-space orbital rotations ({X i j }, {X i j }, and {X a b }) are redundant and can be arbitrary anti-Hermitian matrix elements. On the other hand, inter-space rotations (X a i and X i a = −X a * i ) are non-redundant, and determined by solving the linear 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 is given by 28 , where E E 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, 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 velocitygauge simulations produce identical results for physical observables upon numerical convergence. The velocity gauge is known to be advantageous in simulating highfield phenomena 28,65,66 .
We assume a laser electric field of the following form: 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×10 14 W/cm 2 , 8×10 14 W/cm 2 , and 1×10 15 W/cm 2 for Ne, and 1×10 14 W/cm 2 , 2×10 14 W/cm 2 , 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 coordinate 0 ≤ r ≤ R max , with 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 coordinate 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 10000 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 dipole moment, evaluated as a trace ψ p |ẑ|ψ q ρ q p using the 1RDM. In Fig. 1 (a)(c)(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 the identical results with 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. Figure 1 (b)(d)(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 the figures, while TD-OCCD slightly underestimates the oscillatory behavior. The TD-OCCDT(4) method takes account of 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 [ Fig. 2 (c)(e)], indicating that the electron correlation plays an increasingly important role as the electrons are driven farther from the nucleus. In Fig. 2 (b)(d)(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][69][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 panels (a)(c)(e) we compare the TD-CASSCF, TD-OCCDT, and TDHF methods, while we compare different time-dependent optimized coupledcluster approaches in panels (b)(d)(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 panels (b)(d)(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 [ Fig. 4(b)(d)(f)], the deviation grows with an 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 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 the panels (d) of Figs. 7-12. 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 the figures. The relative deviation of results for each method from 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 10 processors with a The HHG spectra (a, b, 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×10 14 W/cm 2 with various methods.
clock speed of 2.10 GHz. Table I reports the total computation time spent for the simulation corresponding to Fig. 12. This table 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 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 account of the essential part of the triples at an affordable reduced computational cost. The HHG spectra (a, b, 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×10 14 W/cm 2 with various methods.

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 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 account the major part of The HHG spectra (a, b, 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×10 15 W/cm 2 with various methods.
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 thirteen 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.  The HHG spectra (a, b, 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×10 14 W/cm 2 with various methods.

DATA AVAILABLITY
The data that support the findings of this study are available from the corresponding author upon reasonable request. The HHG spectra (a, b, 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×10 14 W/cm 2 with various methods.