Time-dependent coupled-cluster method with time-varying orbital functions, called time-dependent optimized coupled-cluster (TD-OCC) method, is formulated for multielectron dynamics in an intense laser field. We have successfully derived the equations of motion for CC amplitudes and orthonormal orbital functions based on the real action functional, and implemented the method including double excitations (TD-OCCD) and double and triple excitations (TD-OCCDT) within the optimized active orbitals. The present method is size extensive and gauge invariant, a polynomial cost-scaling alternative to the time-dependent multiconfiguration self-consistent-field method. The first application of the TD-OCC method of intense-laser driven correlated electron dynamics in Ar atom is reported.

## INTRODUCTION

The strong-field physics and ultrafast science are rapidly progressing toward a world-changing goal to directly measure and control the electron motion in matters.^{1–5} Reliable theoretical and computational methods are indispensable to understand and predict strong-field and ultrafast phenomena, which frequently involve both bound excitations and ionizations, where dynamical electron correlation plays a key role.

One of the most advanced methods to describe such electron dynamics is the multiconfiguration time-dependent Hartree-Fock (MCTDHF) method^{6–13} and more generally the time-dependent multiconfiguration self-consistent-field (TD-MCSCF) method.^{14,15} In TD-MCSCF, the electronic wavefunction is given by the configuration-interaction (CI) expansion, Ψ_{CI}(*t*) = *∑*_{I}*C*_{I}(*t*)Φ_{I}(*t*), and both CI coefficients {*C*_{I}(*t*)} and *occupied* spin-orbital functions {*ψ*_{p}(*t*)} constituting the Slater determinants {Φ_{I}(*t*)} are propagated in time. (Hereafter spin-orbital functions are simply referred to as orbitals. See, e.g, Refs. 16–19 for TDCI methods using fixed orbitals and Ref. 15 for a broad review of *ab initio* wavefunction-based methods for multielectron dynamics.) Though powerful, the full CI-based MCTDHF method suffers from the factorial scaling of the computational cost with respect to the number of electrons, which restricts its applicability to the systems of moderate size.

The time-dependent complete-active-space self-consistent-field (TD-CASSCF) method^{20} significantly reduces the cost by flexibly classifying occupied orbitals into *frozen-core*, *dynamical-core*, and *active* orbitals as defined in Fig. 1. More approximate and thus less demanding models have been actively investigated,^{21–24} which rely on the truncated CI expansion within the active orbitals given, e.g, as

where $\xcaijk\cdots abc\cdots =\u0109a\u2020\u0109b\u2020\u0109c\u2020\cdots \u0109k\u0109j\u0109i$, with $c\mu \u2020$ and *c*_{μ} being creation and annihilation operators, respectively, for *ψ*_{μ}. (Einstein summation convention is implied throughout for orbital indices.) The operator $\xcaijk\cdots abc\cdots $ excites electron(s) from orbital(s) *i*, *j*, *k*, … occupied in the reference determinant Φ (hole) to those *a*, *b*, *c*, … not occupied in the reference (particle). These methods achieve a polynomial, instead of factorial, cost scaling, and state-of-the-art real-space implementations^{23,25–28} have proved their great utility; however, truncated-CI-based methods share a general drawback of not being size extensive.^{29}

Therefore, one naturally seeks an alternative to Eq. (1), where the truncated CI is replaced with the size-extensive coupled-cluster (CC) expansion^{29–32}

where both excitation amplitudes ${\tau ijk\cdots abc\cdots (t)}$ and orbitals are propagated in time, just as the stationary *orbital-optimized* CC (OCC) method^{33–36} optimizes both excitation amplitudes and orbitals to minimize the CC energy functional. This idea was recently put forward by Kvaal,^{37} who, based on Arponen’s seminal work,^{38} developed a time-dependent coupled-cluster method using time-varying biorthonormal orbitals, called orbital-adapted time-dependent coupled-cluster (OATDCC) method, and applied it to the collision problem in a one-dimensional model Hamiltonian. This method, however, has a difficulty in obtaining the ground state of the system, required as an initial state of simulations, and has not yet been applied to real three-dimensional or time-dependent Hamiltonian. We also notice recent investigations on the time-dependent coupled-cluster methods using fixed orbitals.^{39,40} Their primary focus, however, has been placed on retrieving spectral information rather than the time-domain electron dynamics itself.

In this communication, we report our successful formulation and implementation of the time-dependent coupled-cluster method using time-varying orthonormal orbitals, called time-dependent optimized coupled-cluster (TD-OCC) method, and present its first application to electron dynamics of Ar atom subject to an intense laser field. The Hartree atomic units are used throughout unless otherwise noted.

## TIME-DEPENDENT VARIATIONAL PRINCIPLE WITH REAL CC ACTION FUNCTIONAL

We are interested in the dynamics of a many-electron system interacting with external electromagnetic fields governed by the Hamiltonian written in the second quantization as

where [*x* = (** r**,

*σ*) labels spatial-spin coordinate.]

*h*

_{0}and

*V*

_{ext}are the field-free one-electron Hamiltonian and the interaction of an electron and external fields, given, e.g., in the electric dipole approximation as

*V*

_{ext}=

**(**

*E**t*) ·

**in the length gauge or**

*r**V*

_{ext}= −

*i*

**(**

*A**t*)·

**∇**

_{r}in the velocity gauge, with

**(**

*E**t*) and

**(**

*A**t*) = −

*∫*

^{t}

*dt*′

**(**

*E**t*′) being the electric field and vector potential, respectively.

Note that we include in Eq. (3) not only a given number *N*_{o} of occupied orbitals {*ψ*_{p}} but also $Nv$ *virtual* orbitals {*ψ*_{α}} which are orthogonal to {*ψ*_{p}}. The number *N* = *N*_{o} + $Nv$ of all orbitals {*ψ*_{μ}} is equal to the number of basis functions to expand orbitals (Fig. 1), and in the basis-set limit (*N* → *∞*), *Ĥ* of Eq. (3) is equivalent to the first-quantized Hamiltonian. The index convention given in Fig. 1 is used for orbitals.

Following the previous work,^{29,37,41} we begin with the coupled-cluster Lagrangian *L*,

where $\Phi =\u220fi\u2032\u2032i\u2032i\u0109i\u2032\u2032\u2020\u0109i\u2032\u2020\u0109i\u2020$ is a reference determinant. For the notational brevity, we use composite indices *I* = *ijk*… and *A* = *abc*… to denote general-rank excitation (de-excitation) operators and amplitudes as $\xcaIA$ ($\xcaAI$) and $\tau IA$ ($\lambda AI$). Note that we use the orthonormal orbitals {*ψ*_{μ}} and ⟨Φ| = |Φ⟩^{†}.

As the guiding principle, we adopt the manifestly real action functional,

and require the action to be stationary

with respect to the variation of left- and right-state wavefunctions. This real-action formulation using orthonormal orbitals (also adopted in Ref. 41 for a gauge-invariant response theory) distinguishes our method from the OATDCC method,^{37} which is based on the bivariational principle with the complex-analytic action functional *S* = ∫*dtL* using biorthonormal orbitals.

The time derivative and variation of the left- and right-state wavefunctions can be written as^{20,22,24,42}

where $\u2202t\u2032T^\u2009=\u2009\tau \u0307IA\xcaIA$, $\delta \u2032T^\u2009=\u2009\delta \tau IA\xcaIA$, $\u2202t\u2032\Lambda ^\u2009=\u2009\lambda \u0307AI\xcaAI$, $\delta \u2032\Lambda ^\u2009=\u2009\delta \lambda AI\xcaAI$, $X^=X\nu \mu \xca\nu \mu $, $\Delta ^=\Delta \nu \mu \xca\nu \mu $, $X\nu \mu =\u27e8\psi \mu |\psi \u0307\nu \u27e9$, and $\Delta \nu \mu =\u27e8\psi \mu |\delta \psi \nu \u27e9$. The matrices *X* and Δ are anti-Hermitian to enforce the orthonormality of orbitals $\u27e8\psi \mu (t)|\psi \nu (t)\u27e9=\delta \nu \mu $ at any time.

The equation of motion (EOM) for $\tau IA$ is derived from $\delta S/\delta \lambda AI(t)=0$ and that for $\lambda AI$ from $\delta S/\delta \tau IA(t)=0$ as

The EOMs for orbitals are to be obtained from $\delta S/\delta \Delta \nu \mu (t)=0$. After straightforward yet tedious algebraic manipulations, we found that (i) orbital pairs within the same orbital space (with no further subclassification in Fig. 1) are redundant as well known for the stationary single-reference CC methods, i.e., $Xj\u2032\u2032i\u2032\u2032,Xj\u2032i\u2032,Xji,Xba$ can be arbitrary anti-Hermitian matrix elements, and (ii) the expressions of all the other $X\nu \mu $ terms except the hole-particle contributions $Xia$ (and $Xai=\u2212Xia*$) are formally identical to those for the TD-CASSCF method.^{20} Thus we follow the discussion in Ref. 20 and write the final expression for the orbital EOM, with the hole-particle terms ${Xia}$ left unspecified until the next section,

where $P^=\u2211p\psi p\psi p$, $\u0174sr$ is the mean-field operator,^{20} and ${Xqp}$ except ${Xia}$ are given in Ref. 20. (Note that the Hermitian matrix *R* ≡ *iX* was used as the working variable in Ref. 20. See also Ref. 25 for the velocity gauge treatment of frozen-core orbitals.) We note the natural appearance of the Hermitialized reduced density matrices (RDMs) *D* and *P* is the direct consequence of relying on the real action functional with orthonormal orbitals.

## TD-OCC METHOD

In order to derive a working equation for the hole-particle rotations ${Xia}$, and thus complete the derivation of EOMs, one has to specify a particular ansatz of the CC Lagrangian, i.e., which terms to be included in Eqs. (7) and (8). We, therefore, define the TD-OCC method by the ansatz

that is, with single amplitudes $T^1$ and $\Lambda ^1$ omitted. The same form has been adopted in the ground-state OCC methods^{33–35} and OATDCC method. We have also tried to retain single amplitudes ($T^=T^1+T^2+\cdots \u2009$, $\Lambda ^=\Lambda ^1+\Lambda ^2+\cdots \u2009$) and apply the time-dependent variational principle. However, the resultant EOMs were stable neither in real time nor in imaginary time propagation. A similar difficulty has been reported for the stationary OCC method,^{33} which is in essence due to the similar physical roles played by single excitations and hole-particle rotations.

For the ansatz of Eq. (20), we derive the system of equations to be solved for $Xia$ [by coupling Eqs. (13) and (14) and $\delta S/\delta \Delta ai=0$], as

A special simplification arises when only double amplitudes are included where the cluster operator $eT^$ consists only of even-order (double, quadruple, …) excitations, in which case Eq. (22) and the last two terms of Eq. (23) vanish, and $X^$ needs not be included in Eqs. (13) and (14).

In summary, the TD-OCC EOMs consist of Eq. (13) for $T^2,T^3,\u2026$, Eq. (14) for $\Lambda ^2,\Lambda ^3,\u2026$, and Eq. (15) for orbitals, with ${Xia}$ obtained by solving Eq. (21).

The present TD-OCC method is, as the name suggests, a direct time-dependent extension of the stationary OCC method.^{33–35}

## FEASIBILITY OF IMAGINARY RELAXATION METHOD

We have implemented the TD-OCC method with double excitations (TD-OCCD) and double and triple excitations (TD-OCCDT) for atom-centered Gaussian basis functions and the atomic Hamiltonian with finite-element discrete variable representation (FEDVR) basis.^{25,28} These codes share the same implementation for the basis-independent procedures [Eqs. (13), (14), (17)–(19), and (21)–(23)].

It turns out, in contrast to the OATDCC method, that the imaginary-time relaxation is feasible for the TD-OCC method to obtain the ground state, which is quite convenient when using grids or locally supported basis like FEDVR to expand orbitals. This is confirmed in Table I, which compares the total energies of the ground state of BH molecule with double-*ζ* plus polarization (DZP) basis^{44} obtained by imaginary time propagation with results of corresponding stationary methods.^{35} The OCC and CASSCF methods correlate all six electrons among six active orbitals. The agreement of OCCD energies up to the reported (six) decimal places in Ref. 35 clearly demonstrates our correct implementation and the feasibility of the imaginary-time method.

. | This work^{a}
. | Reference 35 . |
---|---|---|

HF | −25.124 742 010 | −25.124 742 |

OCCD | −25.178 285 704 | −25.178 286 |

OCCDT | −25.178 312 565 | |

CASSCF | −25.178 334 889 | −25.178 335 |

. | This work^{a}
. | Reference 35 . |
---|---|---|

HF | −25.124 742 010 | −25.124 742 |

OCCD | −25.178 285 704 | −25.178 286 |

OCCDT | −25.178 312 565 | |

CASSCF | −25.178 334 889 | −25.178 335 |

^{a}

The overlap, one-electron, and two-electron repulsion integrals over Gaussian basis functions are generated using Gaussian09 program (Ref. 43), and used to propagate EOMs in imaginary time in the orthonormalized Gaussian basis, with a convergence threshold of 10^{−15} Hartree of energy difference in subsequent time steps.

## INITIAL APPLICATIONS

Finally, we report our initial application of the TD-OCC method to Ar atom in a strong laser field. We used the spherical FEDVR basis given as the product *Y*_{lm}(*θ*, *ϕ*)*f*_{k}(*r*) of spherical harmonics *Y*_{lm} with *l* ≤ *l*_{max} and the FEDVR basis *f*_{k}(*r*), which divides the range of radial coordinate 0 < *r* < *R*_{max} into finite elements of length 4 (with several finer elements near the nuclei), each supporting 23 local DVR functions. We first obtained the ground states of TDHF, TD-OCC, and TD-CASSCF methods by the imaginary time relaxation. For the latter two methods, the Ne core was kept frozen at HF orbitals and eight valence electrons were correlated among 13 optimized active orbitals, with their initial characters being 3*s*3*p*3*d*4*p*4*s*. Starting from the ground states, we then simulated the electron dynamics of Ar atom subject to a linearly polarized (in the *z* direction) laser pulse with a wavelength of 800 nm, a peak intensity of 6 × 10^{14} W/cm^{2}, and a foot-to-foot pulse duration of three optical cycles with a sin^{2}-shaped envelope, within the dipole approximation in the velocity gauge. The EOMs are propagated using the fourth-order exponential Runge-Kutta method.^{45} The Ne core was frozen in all the methods and 13 active orbitals were propagated in TD-OCC and TD-CASSCF methods. The basis parameters were calibrated to achieve convergence with *R*_{max} = 240 (with the mask absorbing boundary^{25}) and *l*_{max} = 63.

Figure 2 shows high-harmonic generation (HHG) spectra calculated as the modulus squared of the Fourier transform of the expectation value of the dipole acceleration ⟨*d*^{2}$z^$/*dt*^{2}⟩(*t*) which, in turn, is obtained with the modified Ehrenfest expression^{25} using RDMs of Eq. (17). All methods predicted qualitatively similar spectra, with a minimum at 53 eV (≈34th order, indicated with an arrow) close to the Cooper minimum of photoionization spectrum,^{46} which was experimentally observed.^{47} The TDHF method, however, failed to reproduce fine structures of the TD-CASSCF spectrum especially at higher plateau as seen in the inset of Fig. 2. The description is largely improved by TD-OCCD, and the TD-OCCDT spectrum well reproduces the TD-CASSCF one in virtually all details.

We also compare the probabilities of finding one [Fig. 3(a)] and two [Fig. 3(b)] electron(s) outside a sphere of radius *R*_{0} = 20, which measure the single and double ionization probabilities, respectively. These probabilities are much more sensitive to the description of correlated motions of electrons than HHG, which hinders a correct description with TDHF. The rapid convergence of TD-OCC results to the TD-CASSCF ones for such nonperturbatively nonlinear processes strongly promises the TD-OCC method to be a vital tool to investigate correlated high-field phenomena.

## SUMMARY

We have successfully formulated a new time-dependent coupled-cluster method called the TD-OCC method and implemented the first two, and the most important, series of approximations, TD-OCCD and TD-OCCDT. The present method is size extensive and gauge invariant, a polynomial cost-scaling alternative to the TD-MCSCF method. It would open new possibilities of high-accuracy first-principle investigations of multielectron dynamics in ever-unreachable large target systems. The rigorous derivation, details of the implementation, as well as other ansatz for TD-CC theories and comparison thereof, will be presented in a forthcoming article.

## ACKNOWLEDGMENTS

This research was supported in part by a Grant-in-Aid for Scientific Research (Grants Nos. 25286064, 26390076, 26600111, 16H03881, and 17K05070) from the Ministry of Education, Culture, Sports, Science and Technology (MEXT) of Japan and also by the Photon Frontier Network Program of MEXT. This research was also partially supported by the Center of Innovation Program from the Japan Science and Technology Agency, JST, and by CREST (Grant No. JPMJCR15N1), JST. Y.O. gratefully acknowledges support from the Graduate School of Engineering, the University of Tokyo, Doctoral Student Special Incentives Program (SEUT Fellowship).