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

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¹⁴ W/cm² 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.⁹ 

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) = ∑_IC_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)³² have been developed. These methods have proven to be of great utility³⁴ 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.³⁸ 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.³⁸ In an earlier development, Kvaal reported an orbital adaptive time-dependent coupled-cluster (OATDCC) method⁴² built upon the work of Arponen using biorthogonal orbitals.⁴³ 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).⁵⁶ This method is much more efficient than time-dependent optimized coupled-cluster with double (TD-OCCD) excitations, even though both scale as O(N⁶), 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⁵), 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.⁵⁸ The stationary variant of this method (CC2)⁵⁹ 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.⁵⁷ 

Thus, the TD-OCCD method and its lower-cost approximations such as TD-OCEPA0⁵⁶ and TD-OMP2⁵⁷ 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.³⁸ 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⁸) as in the stationary CCSDT method.⁶³ 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.³⁷ 

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⁷), 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

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 $ĉμ†$ (⁠$ĉμ$⁠) represents a creation (annihilation) operator in a complete, orthonormal set of 2n_bas spin orbitals {ψ_μ(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 μ, ν, γ, λ) 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 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²⁷ (see Fig. 1 of Ref. 38 for a pictorial illustration of the orbital subspacing).

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

where

where $τij…ab…$ and $λab…ij…$ are excitation and de-excitation amplitudes, respectively. The stationary conditions, δS = 0, with respect to the variation of amplitudes $δτij…ab…$⁠, $δλab…ij…$⁠, 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̂=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

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

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

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

The full TD-OCCDT Lagrangian is given by³⁸ 

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 L_CCDT by retaining those terms contributing up to the fourth order in the many-body perturbation theory,

where $f̄=f̂−iX̂$⁠, $f̂=(hqp+vqjpj){Êqp}$⁠, $v̂=vqspr{Êqspr}/4$⁠, and $vqspr=uqspr−usqpr$⁠.

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

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̂=∑μ|ψμ〉〈ψμ|$ 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 $Xpq$ describes orbital rotations among various orbital subspaces. Intra-space orbital rotations (⁠${Xj′i′}$⁠, ${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 equation³⁸ 

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²⁸ 

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

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:

for 0 ≤ t ≤ 3T, and E(t) = 0 otherwise, with a central wavelength of λ = 2π/ω₀ = 800 nm and a period of T = 2π/ω₀ ∼ 2.67 fs. We consider three different intensities 5 × 10¹⁴, 8 × 10¹⁴, and 1 × 10¹⁵ W/cm² for Ne and 1 × 10¹⁴, 2 × 10¹⁴, and 4 × 10¹⁴ W/cm² 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⁶⁷ 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.

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¹⁴ W/cm² 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(ω)|² of the Fourier transform of the expectation value of the dipole acceleration, which, in turn, is evaluated with a modified Ehrenfest expression.²⁸ 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.

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

The matrix B is given by

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