Time-dependent optimized coupled-cluster method for multielectron dynamics
Takeshi Sato Himadri Pathak Yuki Orimo Kenichi L. Ishikawa
physics.chem-ph




















[Electronic mail:]sato@atto.t.u-tokyo.ac.jp

Photon Science Center, School of Engineering, 
The University of Tokyo, 7-3-1 Hongo, Bunkyo-ku, Tokyo 113-8656, Japan


Department of Nuclear Engineering and Management, School of Engineering,
The University of Tokyo, 7-3-1 Hongo, Bunkyo-ku, Tokyo 113-8656, Japan


[Electronic mail:]pathak@atto.t.u-tokyo.ac.jp

Department of Nuclear Engineering and Management, School of Engineering,
The University of Tokyo, 7-3-1 Hongo, Bunkyo-ku, Tokyo 113-8656, Japan


[Electronic mail:]ykormhk@atto.t.u-tokyo.ac.jp

Department of Nuclear Engineering and Management, School of Engineering,
The University of Tokyo, 7-3-1 Hongo, Bunkyo-ku, Tokyo 113-8656, Japan


[Electronic mail:]ishiken@n.t.u-tokyo.ac.jp

Department of Nuclear Engineering and Management, School of Engineering,
The University of Tokyo, 7-3-1 Hongo, Bunkyo-ku, Tokyo 113-8656, Japan


Photon Science Center, School of Engineering, 
The University of Tokyo, 7-3-1 Hongo, Bunkyo-ku, Tokyo 113-8656, Japan










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 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 to intense-laser driven correlated electron dynamics in Ar atom is reported.




Time-dependent optimized coupled-cluster method for multielectron dynamics


    Kenichi L. Ishikawa
    December 30, 2023
============================================================================








§ INTRODUCTION

The strong-field physics and ultrafast science are rapidly progressing

towards a world-changing goal to directly measure and control the electron motion in
matters.<cit.> 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,<cit.> 
and more generally the time-dependent multiconfiguration self-consistent-field (TD-MCSCF) method.<cit.> 
In TD-MCSCF, the electronic wavefunction is given by the
configuration-interaction (CI) expansion, Ψ_ CI(t)=∑_ IC_ I(t)Φ_ I(t), 
and both CI coefficiets {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. 
for TDCI methods using fixed orbitals, and Ref. 
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<cit.> significantly reduces the cost by 
flexibly classifying occupied orbitals into frozen-core,

dynamical-core, 

and active orbitals as defined in Fig. <ref>. 
More approximate and thus less demanding models have been actively investigated,
<cit.> which rely on the truncaed CI expansion
within the active orbitals given, e.g, as

    Ψ_ CI
    
    
    = (C_0+C^a_iÊ^a_i+C^ab_ijÊ^ab_ij+C^abc_ijkÊ^abc_ijk+⋯)Φ,

where (Einstein's summation convention is implied throughout for orbital indices.) 
Ê^abc⋯_ijk⋯=ĉ^†_aĉ^†_bĉ^†_c⋯ĉ_kĉ_jĉ_i
[c^†_μ (c_μ) is a creation (anihilation) operator for
ψ_μ] 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), achieving a polynomial cost scaling.
Even though state-of-the-art real-space implementations<cit.>
have proved a great utility of these methods,

all the methods based on a truncated-CI expansion share a common drawback of
not being size extensive.<cit.>



Therefore one naturally seeks an alternative to Eq. (<ref>), where the truncated CI is
replaced with the size-extensize coupled-cluster (CC) 
expansion<cit.>

    Ψ_ CC = exp(τ^a_iÊ^a_i+τ^ab_ijÊ^ab_ij+τ^abc_ijkÊ^abc_ijk+⋯)Φ,

where both excitation amplitudes {τ^abc⋯_ijk⋯(t)} and orbitals are propagated in time, 
just as the stationary orbital-optimized CC (OCC) method<cit.>
optimizes both excitation amplitudes and orbitals to minimize the CC energy functional.

This idea was recently put forward by Kvaal<cit.>, who, based on Arponen's seminal work,<cit.> 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.<cit.>

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

    Ĥ(t)=h^μ_ν(t)Ê^μ_ν+1/2g^μγ_νδÊ^μγ_νδ,
     
                 h^μ_ν(t)=∫ dxψ^*_μ [h_0+V_ ext(t)]ψ_ν, 
          g^μγ_νδ=∫∫ dx_1 dx_2ψ^*_μ(x_1)ψ^*_γ(x_2) r^-1_12ψ_ν(x_1)ψ_δ(x_2),

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)·r in the length gauge or
V_ ext=-iA(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. (<ref>) not only a given number N_o of occupied orbitals
{ψ_p} 
but also N_v virtual orbitals {ψ_α} 
which are orthogonal to {ψ_p}. The number N=N_o+N_v of all
orbitals {ψ_μ} is equal to the number of basis functions to
expand orbitals (Fig <ref>), and in the basis-set limit (N→∞), Ĥ of Eq. (<ref>) is equivalent to the first-quantized
Hamiltonian. The index convention given in Fig. <ref> is used for orbitals.





Following the previous work<cit.>, we begin with the
coupled-cluster Lagrangian L,

    L    =   ⟨Ψ_ L|(Ĥ-i∂_t)|Ψ_ R⟩
       =   ⟨Φ|(1+Λ̂)e^-T̂
    (Ĥ-i∂_t)
    e^T̂|Φ⟩, 
    
    
    T̂   =   τ_i^a Ê_i^a +
    τ_ij^abÊ_ij^ab +
    τ_ijk^abcÊ_ijk^abc + ⋯ = τ_I^A Ê_I^A 
       =   T̂_1+T̂_2+T̂_3+⋯, 
    Λ̂   =   λ^i_a Ê^i_a +
    λ^ij_abÊ^ij_ab +
    λ^ijk_abcÊ^ijk_abc + ⋯  = λ^I_A Ê^I_A 
       =   Λ̂_1+Λ̂_2+Λ̂_3+⋯,

where
|Φ⟩=∏_i^''i^'iĉ^†_i^''ĉ^†_i^'ĉ^†_i|⟩ 
is a reference determinant.
For the notational brevity, we use composite indices I=ijk⋯ and A=abc⋯ to denote general-rank excitation
(deexcitation) operators and amplitudes as Ê^A_I (Ê^I_A)
and τ^A_I (λ^I_A). Note that we use the orthonormal orbitals {ψ_μ} and
⟨Φ|=|Φ⟩^†.

As the guiding principle,

we adopt the manifestly real action functional, 

    S    =   ∫_t_0^t_1 Ldt = 1/2∫_t_0^t_1( L + L^*) dt,

and require the action to be stationary

    δ S    =   ⟨δΨ_ L|Ĥ|Ψ_ R⟩
    +⟨Ψ_ L|Ĥ|δΨ_ R⟩
       -   i⟨δΨ_ L|Ψ̇_ R⟩
    +i⟨Ψ̇_ L|δΨ_ R⟩ + c.c = 0,


with respect to the variation of left- and right-state wavefunctins.
This real-action formulation using orthonormal orbitals 
(also adopted in Ref.  for a gauge-invariant response theory)
distinguishes our method from the
OATDCC method,<cit.> which is based on the bivariational principle with the complex-analytic action functional S=∫ dt L using biorthonormal orbitals.




The time derivative and variation of the left- and right-state wavefunctions can be written as<cit.>
 

    |Ψ̇_ R⟩ = 
    (T̂ +  X̂)e^T̂|Φ⟩,   
       |δΨ_ R⟩ = 
    (T̂ + Δ̂)e^T̂|Φ⟩,   
       ⟨Ψ̇_ L| =⟨Φ|[-(1+Λ)e^-T̂(T̂+X̂)+Λ̂e^-T̂],   
       ⟨δΨ_ L| =⟨Φ|[-(1+Λ)e^-T̂(T̂+Δ̂)+Λ̂e^-T̂],


where 
T̂=τ̇^A_IÊ^A_I,
T̂=δτ^A_IÊ^A_I,
Λ̂=λ̇^I_AÊ^I_A,
Λ̂=δλ^I_AÊ^I_A,
X̂=X^μ_νÊ^μ_ν,
Δ̂=Δ^μ_νÊ^μ_ν, 
X^μ_ν = ⟨ψ_μ|ψ̇_ν⟩, and
Δ^μ_ν = ⟨ψ_μ|δψ_ν⟩.





The matrices X and Δ
are anti-Hermitian to enforce the orthonormality of orbitals
⟨ψ_μ(t)|ψ_ν(t)⟩=δ^μ_ν at any time.







We insert Eqs. (<ref>) into Eq. (<ref>) and
group terms for each kind of variation to obtain


 

    δ S 
    
    
    
    
    
       =   δτ^A_I{⟨Φ|(1+Λ)e^-T̂
    [Ĥ-iX̂,Ê^A_I]e^T̂|Φ⟩
    +iλ̇^I_A}
    +δλ^I_A{⟨Φ^A_I|e^-T̂(Ĥ-iX̂)e^T̂|Φ⟩
    -iτ̇^A_I} + c.c 
    
    
    
    
    
    
    
    
    
       +   Δ^μ_ν{⟨Φ|(1+Λ)e^-T̂[Ĥ-i(T̂)-iX̂,Ê^μ_ν]e^T̂|Φ⟩
    +i⟨Φ|(Λ̂)e^-T̂Ê^μ_νe^T̂|Φ⟩. 
          .-
    ⟨Φ|(1+Λ)e^-T̂[Ĥ-i(T̂)-iX̂,Ê^ν_μ]e^T̂|Φ⟩^*
    -i⟨Φ|(Λ̂)e^-T̂Ê^ν_μe^T̂|Φ⟩^*
    },



where the anti-Hermiticity of Δ is
used for deriving Eq. (<ref>). 

The EOMs for the amplitudes
τ^A_I and λ^I_A are readily derived from the
conditions δ S/δλ^I_A=0 and δ S/δτ^A_I=0, 
respectively, as

    iτ̇^A_I    =   ⟨Φ_I^A| e^-T̂ (Ĥ-iX̂)e^T̂|Φ⟩, 
    
    -iλ̇^I_A    =   ⟨Φ|(1+Λ̂) e^-T̂ [Ĥ-iX̂,Ê^A_I] e^T̂|Φ⟩.


The EOMs for orbitals are to be obtained from δ S/δΔ^μ_ν=0.




After straightforward yet tedious algebraic
manipulations<cit.>, we found that (i) orbital pairs within the same
orbital space (with no further subclassification in Fig. <ref>)
are redundant as well known for the stationary single-reference CC methods, i.e,
X^i^''_j^'',X^i^'_j^',X^i_j,X^a_b can be arbitrary anti-Hermition
matrix elements, and (ii) the expressions of all the other X^μ_ν terms except
the hole-particle contributions X^a_i (and X^i_a=-X^a*_i) are formally identical to those for the TD-CASSCF method.<cit.>







Thus we follow the discussion in Ref. , and write
the final expression for the orbital EOM, with the hole-particle
terms {X^a_i} left unspecified until next section, 

    i|ψ̇_̇ṗ⟩ = 
    (1-P̂) F̂_p|ψ_p⟩ + i|ψ_q⟩ X^q_p,   
       F̂_p|ψ_p⟩ = ĥ |ψ_p⟩ + Ŵ^r_s|ψ_q⟩ P^qs_or (D^-1)^o_p,   
    
    
       D^p_q=1/2(ρ^p_q+ρ^q*_p), 
    P^pr_qs=1/2(ρ^pr_qs+ρ^qs*_pr),   
       ρ^p_q = ⟨Φ|(1+Λ)e^-T̂Ê^q_pe^T̂|Φ⟩,   
       ρ^pr_qs = ⟨Φ|(1+Λ)e^-T̂Ê^qs_pre^T̂|Φ⟩,

where P̂=∑_p|ψ_p⟩⟨ψ_p|, 
Ŵ^r_s is the mean-field operator<cit.>,
and {X^p_q} except {X^a_i} are given in Ref. . 
(Note that the Hermitian matrix R≡ iX was
used as the working variable in Ref..
See also Ref.  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 consequece 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
{X^a_i}, 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. (<ref>) and (<ref>). We, therefore, define the TD-OCC method by the ansatz

    T̂ = T̂_2+T̂_3+⋯, Λ̂ = Λ̂_2+Λ̂_3+⋯,

that is, with single amplitudes T̂_1 and Λ̂_1 omitted.
The same form has been adopted in the ground-state OCC methods<cit.> and OATDCC method. 

We have also tried to retain single amplitudes (T̂=T̂_1+T̂_2+⋯, 
Λ̂=Λ̂_1+Λ̂_2+⋯) and apply TDVP of the previous section. 
However, the resultant EOMs were stable neither in real time nor imaginary time propagation.<cit.>
A similar difficulty has been reported for the stationary OCC
method,<cit.> which is in essence due to the similar physical
roles played by single excitations and hole-particle rotations. 

For the ansatz of Eq. (<ref>), we derive the system of equations to be solved for X^a_i
[by coupling Eqs. (<ref>), (<ref>) and
δ S/δΔ^i_a = 0],<cit.> as

    i{(δ^a_bD^j_i-D^a_bδ^j_i)X^b_j+1/2A^aj_ib  X^b*_j} = 
    B^a_i,   
       A^aj_ib = ⟨Φ|(1+Λ̂) e^-T̂ [Ê^A_I,Ê^j_b] e^T̂|Φ⟩⟨Φ^A_I|e^-T̂Ê^i_ae^T̂|Φ⟩   
       -
    ⟨Φ|(1+Λ)e^-T̂[Ê^A_I,Ê^i_a]e^T̂|Φ⟩⟨Φ_I^A| e^-T̂Ê^j_be^T̂|Φ⟩,   
       B^a_i = 
    F^a_jD^j_i - D^a_bF^i*_b   
         +
    1/2⟨Φ|(1+Λ̂) e^-T̂ [Ê^A_I,Ĥ] e^T̂|Φ⟩
    ⟨Φ^A_I|e^-T̂Ê^i_ae^T̂|Φ⟩   
         -
    1/2⟨Φ|(1+Λ)e^-T̂[Ê^A_I,Ê^i_a]e^T̂|Φ⟩
    ⟨Φ_I^A| e^-T̂Ĥe^T̂|Φ⟩.

A special simplification arises when only double amplitudes are
included, in which case Eq. (<ref>) and the last two terms of
Eq. (<ref>) vanish, and X̂ needs not be included in
Eqs. (<ref>) and (<ref>). 

In summary, the TD-OCC EOMs consist of
Eq. (<ref>) for T̂_2,T̂_3,⋯, 
Eq. (<ref>) for Λ̂_2,Λ̂_3⋯, and
Eq. (<ref>) for orbitals, with {X^a_i} obtained by
solving Eq. (<ref>). 

The present TD-OCC method is, as the name suggests, a direct time-dependent extension of
the stationary OCC method<cit.>.
































§ 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<cit.>. These
codes share the same implementation for the basis-independent
procedures [Eqs. (<ref>), (<ref>),
(<ref>)-(<ref>), (<ref>)-(<ref>)].

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 <ref>, which compares the
total energies of the ground state of BH molecule with
double-ζ plus polarization (DZP) basis<cit.> obtained by imaginary
time propagation 
with results of corresponding stationary methods<cit.>.

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.  clearly demonstrates our correct
implementation and the feasibility of the imaginary-time method.







§ 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 3s3p3d4p4s.





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×1014
W/cm2, and a foot-to-foot pulse duration of three optical cycles with a
sin2-shaped envelope, within the dipole approximation in the velocity gauge.
The EOMs are propageted
using the fourth-order exponential Runge-Kutta method.<cit.> 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 absorbing boundary<cit.>) and l_ max=63.

Figure <ref> 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ẑ/dt^2⟩(t) which, in turn, 
is obtained with the modified Ehrenfest expression<cit.> using RDMs of Eq. (<ref>).
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,<cit.> which was
experimentally observed.<cit.> 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. <ref>. 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. <ref> (a)] and two  [Fig. <ref> (b)] electron(s) outside a sphere of radius R_0=20, 
which measure the single and double ionization probabilities, respectively.
These probabilites are much more sensitive to the description of correlated motions of electrons
than HHG,

which hinders a correct description with TDHF. 
The rapid converngence 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 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.<cit.>


This research was supported in part by a Grant-in-Aid for Scientific Research 
(Grants No. 25286064, No. 26390076, No. 26600111, No. 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).



