Implementation of gauge-invariant time-dependent configuration interaction singles method for three-dimensional atoms
IImplementation of gauge-invariant time-dependent configuration interaction singlesmethod for three-dimensional atoms
Takuma Teramura , ∗ Takeshi Sato , , , † and Kenichi L. Ishikawa , , ‡ Department of Nuclear Engineering and Management,Graduate School of Engineering, The University of Tokyo,7-3-1 Hongo, Bunkyo-ku, Tokyo 113-8656, Japan Photon Science Center, Graduate School of Engineering,The University of Tokyo, 7-3-1 Hongo, Bunkyo-ku, Tokyo 113-8656, Japan and Research Institute for Photon Science and Laser Technology,The University of Tokyo, 7-3-1 Hongo, Bunkyo-ku, Tokyo 113-0033, Japan (Dated: August 29, 2019)We present a numerical implementation of the gauge-invariant time-dependent configuration in-teraction singles (TDCIS) method [Appl. Sci. , 433 (2018)] for three-dimensional atoms. In ourimplementation, orbital-like quantity called channel orbital [Phys. Rev. A , 043420 (2006)] ispropagated instead of configuration-interaction (CI) coefficients, which removes a computationalbottleneck of explicitly calculating and storing numerous virtual orbitals. Furthermore, besides itsphysical consistency, the gauge-invariant formulation allows to take advantages of the velocity gaugetreatment of the laser-electron interaction over the length gauge one in the simulation of high-fieldphenomena. We apply the present implementation to high-harmonic generation from helium andneon atoms, and numerically confirms the gauge invariance and demonstrates the effectiveness ofthe rotated velocity gauge treatment. I. INTRODUCTION
Recent laser technology capable of generating stronglaser pulses with an intensity (cid:38) W / cm has en-abled us to explore electron dynamics in nonpertur-bative regime, e.g., high-harmonic generation (HHG),above threshold ionization, nonsequential double ioniza-tion, and attosecond pulse generation [1–3]. While laser-driven electron dynamics is rigorously described by thetime-dependent Schr¨odinger equation (TDSE), its directnumerical solution is practically unfeasible for systemswith more than two electrons. For theoretical investiga-tion of multielectron dynamics in intense laser field, var-ious tractable ab-initio methods have been developed,e.g., time-dependent multiconfiguration self-consistentfield (TD-MCSCF) methods [4–16], time-dependent cou-pled cluster method [17, 18], time-dependent R − matrixapproach [19–23], and time-dependent reduced two bodydensity matrix approach[24, 25].Among them, the time-dependent configuration inter-action singles (TDCIS) method is one of the promissingmethods [26–36]. This method has been successfully ap-plied to various electron dynamics such as giant enhance-ment in HHG in Xe [35] and decoherence in attosecondphotoionization [27]. In the TDCIS method, the totalelectronic wavefunction is approximated by a superposi-tion of time-independent Slater determinants | Ψ( t ) (cid:105) = | Φ (cid:105) C ( t ) + occ (cid:88) i vir (cid:88) a | Φ ai (cid:105) C ai ( t ) , (1) ∗ [email protected] † [email protected] ‡ [email protected] where | Φ (cid:105) is the Hartree-Fock (HF) ground state and | Φ ai (cid:105) is a singly-excited configuration replacing an occu-pied orbital φ i with a virtual orbital φ a unoccupied inthe ground state. The orbital functions are fixed andpropagation of configuration-interaction (CI) coefficients( C and { C ai } ) describes the system dynamics. Althoughapplications of the TDCIS method are limited to the sin-gle excitation or ionization due to the truncation of CIspace, its low computational cost and ease of analysis areattractive.The conventional TDCIS method with CI coefficientshas two major issues; the explicit calculation and storageof virtual orbitals { φ a } and a violation of gauge invari-ance. Virtual orbitals { φ a } should include both boundand continuum orbitals, whose number is infinite in prin-ciple. In a practical simulation with real-space grids,one has to prepare virtual orbitals in advance by nu-merically obtaining all eigenstates of the discretized HFequation. The number of the virtual orbitals increaseswith the number of the grid points. Thus, the calculationand storage of virtual orbitals become unacceptably de-manding for systems with a large number of grid pointslike molecules. To solve this problem, Rohringer et.al. have proposed an alternative but equivalent formulationof the TDCIS method in which time-dependent orbital-like quantity called channel orbital is propagated insteadof CI coefficient [26]. The channel orbital is defined byusing { φ a } and C ai as | χ i ( t ) (cid:105) ≡ vir (cid:88) a | φ a (cid:105) C ai ( t ) . (2)The equations of motions (EOMs) of CI coefficient areconverted to those of C and { χ i } . This reformulationremoves computational bottleneck of handling numerous a r X i v : . [ phy s i c s . a t o m - ph ] A ug virtual orbitals, while in princple including all the virtualorbitals within a grid space, and has enhanced utility ofthe TDCIS method. However, to the best of our knowl-edge, the applications of the channel orbital-based TD-CIS method have been limited to one-dimensional (1D)model in Refs. [26, 37, 38] and nobles gas atoms withHartree-Slater potential in Ref. [27].The TDCIS method, either CI coefficient-based orchannel-orbital-based, suffers from a violation of gaugeinvariance, as a general consequence of the truncation ofCI space. Although it is known that the velocity gauge(VG) offers efficient simulations of high-field phenomena,it was impossible to enjoy this advantage within the TD-CIS method. To overcome this difficulty, we have recentlyreported a gauge-invariant reformulation of the TDCISmethod [39]. In our reformulation, a rotated velocitygauge (rVG) transformed from the length gauge (LG) bya unitary operator has been introduced. This unitarytransformation ensures the gauge invariance between theLG and rVG, and Ref. [39] numerically confirmed theequivalence of these gauges for a model 1D Hamiltonian.In this paper, we report a three dimensional numericalimplementation of the gauge-invariant TDCIS methodfor atoms subject to a linearly polarized laser pulse. Weemploy a spherical harmonics expansion of orbital func-tions with the radial coordinate discretized by a finite-element discrete variable representation (FEDVR)[40–43]. We apply the present implementation to HHG fromhelium and neon atoms and asses the advantage of therVG over the LG and VG.This paper proceeds as follows. In Sec. II, we brieflyreview the TDCIS methods. The numerical implemen-tation of the gauge-invariant TDCIS method to three-dimensional atoms is given in Sec. III. We describe nu-merical applications in Sec. IV and conclude this work inSec. V. We use Hartree atomic units (a.u.) throughoutthe paper unless otherwise stated. II. THEORYA. The System Hamiltonian and Gaugetransformation
We consider an atom with N electrons with a nucleuslocated at the origin. The time evolution of the N -electron wavefunction | Ψ( t ) (cid:105) is governed by the TDSE, i ∂∂t | Ψ( t ) (cid:105) = ˆ H ( t ) | Ψ( t ) (cid:105) , (3)where ˆ H ( t ) is the time-dependent non-relativistic Hamil-tonian ˆ H ( t ) = ˆ H + ˆ H ext ( t ) , (4)decomposed into the field-free part,ˆ H = N (cid:88) i =1 ˆ h ( r i , p i ) + N (cid:88) i>j | r i − r j | (5) and the laser-electron interaction partˆ H ext ( t ) = N (cid:88) i =1 ˆ h ext ( r i , p i , t ) , (6)In these expressions, r i and p i = − i ∇ i are the positionand the canonical momentum of the electron i , respec-tively. ˆ h is given by,ˆ h ( r , p ) = p − Z | r | , (7)where Z is the atomic number. Within the electric dipoleapproximation, ˆ h ext for the LG and VG are given byˆ h LGext ( r , p , t ) = E ( t ) · r (8a)ˆ h VGext ( r , p , t ) = A ( t ) · p , (8b)where E ( t ) and A ( t ) = − (cid:82) t −∞ dt (cid:48) E ( t (cid:48) ) are the electricfield and the vector potential of the external laser field,respectively.The two gauges are physically equivalent, and anyphysical observable takes the same value, independent ofthe choice of the gauge. The LG wavefunction (cid:12)(cid:12) Ψ LG (cid:11) andVG wavefunction (cid:12)(cid:12) Ψ VG ( t ) (cid:11) are mutually transformed bya gauge transformation as, (cid:12)(cid:12) Ψ VG ( t ) (cid:11) = ˆ U ( t ) (cid:12)(cid:12) Ψ LG ( t ) (cid:11) (9)ˆ U ( t ) ≡ exp (cid:34) − i N (cid:88) i =1 (cid:18) A ( t ) · r i − (cid:90) t −∞ dt (cid:48) | A ( t (cid:48) ) | (cid:19)(cid:35) (10) B. The CI coefficient-based TDCIS method in thelength gauge
In the conventional TDCIS method based on CI coeffi-cients, orbitals satisfy the canonical, restricted HF equa-tion ˆ f | φ p (cid:105) ≡ ˆ h | φ p (cid:105) + occ (cid:88) i (2 ˆ W φ i φ i | φ p (cid:105) − ˆ W φ i φ p | φ i (cid:105) )= (cid:15) p | φ p (cid:105) , (11)where ˆ f is the Fock operator and ˆ W φφ (cid:48) is the potential dueto the product of two given orbitals φ and φ (cid:48) , defined inthe real space asˆ W φφ (cid:48) ( r ) ≡ (cid:90) d r φ ∗ ( r ) φ (cid:48) ( r ) | r − r | . (12) (cid:15) p is the orbital energy of orbital φ p . In the TDCIS wave-function in Eq. (1), | Φ (cid:105) is the HF ground state formedwith the occupied orbitals as | Φ (cid:105) = occ (cid:89) i ˆ c † i ↑ ˆ c † i ↓ |(cid:105) , (13)where ˆ c † pσ and ˆ c pσ are the creation and annihilation op-erators, respectively, of spin-orbital φ p ⊗ σ , and |(cid:105) is thevacuum. σ ∈ {↑ , ↓} denotes the spin function. | Φ ai (cid:105) is asingly-excited configuration replacing an occupied orbital φ i with a virtual orbital φ a | Φ ai (cid:105) = 1 √ c † a ↑ ˆ c i ↑ + ˆ c † a ↓ ˆ c i ↓ ) | Φ (cid:105) . (14)The EOMs of CI coefficients is derived through theDirac-Frenkel time-dependent variational principle [44],requiring Lagrangian L ( t ) L ( t ) = (cid:104) Ψ | ˆ H ( t ) − i ∂∂t − E | Ψ (cid:105) (15)to be stationary with respect to the variation of C ∗ and { C a ∗ i } . E = (cid:104) Φ | ˆ H | Φ (cid:105) denotes the HF energy. Thisconstant shift, introduced for the simple notation of theEOMs, does not affect physical results. In the LG case,in which the wavefuntion is written as, (cid:12)(cid:12) Ψ LG ( t ) (cid:11) = | Φ (cid:105) C + occ (cid:88) i vir (cid:88) a | Φ ai (cid:105) C ai . (16)the EOMs of CI coefficients are obtained as i ˙ C = √ E · occ (cid:88) j vir (cid:88) b (cid:104) φ j | r | φ b (cid:105) C bj (17a) i ˙ C ai = (cid:104) φ a | { occ (cid:88) j vir (cid:88) b ˆ F ij | φ b (cid:105) C bj + vir (cid:88) b E · r | φ b (cid:105) C bi + √ E · r | φ i (cid:105) C } − E · occ (cid:88) j (cid:104) φ j | r | φ i (cid:105) C aj , (17b)whereˆ F ij | φ b (cid:105) = δ ij ( ˆ f − (cid:15) i ) | φ b (cid:105) + 2 ˆ W φ j φ b | φ i (cid:105) − ˆ W φ j φ i | φ b (cid:105) . (18) C. The channel orbital-based TDCIS method inthe length gauge
The EOMs of CI coefficients [Eq. (17)] can be rewrit-ten, by substituting channel orbital Eq. (2) into Eq. (17),as, i ˙ C = √ E · occ (cid:88) j (cid:104) φ j | r | χ j (cid:105) (19a) i | ˙ χ i (cid:105) = ˆ P (cid:110) ( ˆ F + E · r ) | χ i (cid:105) + √ E · r | φ i (cid:105) C (cid:111) − E · occ (cid:88) j | χ j (cid:105) (cid:104) φ j | r | φ i (cid:105) , (19b)whereˆ F | χ i (cid:105) = ( ˆ f − (cid:15) i ) | χ i (cid:105) + occ (cid:88) j (2 ˆ W φ j χ j | φ i (cid:105) − ˆ W φ j φ i | χ j (cid:105) ) , (20) and ˆ P is the projection operator onto the space spannedby virtual orbitalsˆ P = vir (cid:88) a | φ a (cid:105) (cid:104) φ a | = ˆ1 − occ (cid:88) j | φ j (cid:105) (cid:104) φ j | , (21)with ˆ1 being the identity operator. D. Velocity gauge and rotated velocity gauge
One can, in principle, derive the EOMs for the VGcase in the same way as for the LG. Let us write thetotal wavefunction and channel orbital as, (cid:12)(cid:12) Ψ VG ( t ) (cid:11) = | Φ (cid:105) D + occ (cid:88) i vir (cid:88) a | Φ ai (cid:105) D ai (22) | η ( t ) (cid:105) = vir (cid:88) a | φ a (cid:105) D ai . (23) | Φ (cid:105) and | Φ ai (cid:105) are the same configurations as those usedin the LG case. Their EOMs are obtained through thesame procedures as in the LG case as, i ˙ D = √ A · occ (cid:88) j (cid:104) φ j | p | η j (cid:105) (24a) i | ˙ η i (cid:105) = ˆ P (cid:110) ( ˆ F + A · p ) | χ i (cid:105) + √ A · p | φ i (cid:105) C (cid:111) − A · occ (cid:88) j | η j (cid:105) (cid:104) φ j | p | φ i (cid:105) . (24b)It is known that TDCIS, which uses time-independentorbitals, is not gauge-invariant [39, 45, 46]. Instead ofthe conventional VG as described above, we have recentlyproposed the rVG [39], where we define the rVG wave-function by the gauge transformation from the LG wave-function as, (cid:12)(cid:12) Ψ rVG ( t ) (cid:11) = ˆ U ( t ) (cid:12)(cid:12) Ψ LG ( t ) (cid:11) (25)The rVG orbitals are related to the LG ones by, (cid:12)(cid:12) φ (cid:48) p ( t ) (cid:11) = ˆ u ( t ) | φ p (cid:105) (26) | χ (cid:48) i ( t ) (cid:105) = ˆ u ( t ) | χ i (cid:105) = vir (cid:88) a | φ (cid:48) a (cid:105) C ai , (27)whereˆ u ( t ) = exp (cid:26) − i (cid:18) A ( t ) · r − (cid:90) t −∞ dt (cid:48) | A ( t (cid:48) ) | (cid:19)(cid:27) . (28)They satisfy the following EOMs [39]: i ˙ C = √ E · occ (cid:88) j (cid:10) φ (cid:48) j (cid:12)(cid:12) r (cid:12)(cid:12) χ (cid:48) j (cid:11) (29a) i | ˙ χ (cid:48) i (cid:105) = ˆ P (cid:48) (cid:110) ( ˆ F (cid:48) + A · p ) | χ (cid:48) i (cid:105) + √ E · r | φ (cid:48) i (cid:105) C (cid:111) − E · occ (cid:88) j (cid:0)(cid:12)(cid:12) χ (cid:48) j (cid:11) (cid:10) φ (cid:48) j (cid:12)(cid:12) r | φ (cid:48) i (cid:105) + (cid:12)(cid:12) φ (cid:48) j (cid:11) (cid:10) φ (cid:48) j (cid:12)(cid:12) r | χ (cid:48) i (cid:105) (cid:1) , (29b)where ˆ P (cid:48) and ˆ F (cid:48) i are given by Eqs. (21) and (20), respec-tively, with { φ j } replaced with (cid:8) φ (cid:48) j (cid:9) . Although Eq. (29)contains the dipole operator E · r , this does not preventenjoying the advantages of the VG treatment, since itacts only on the localized occupied orbitals { φ (cid:48) i } . III. IMPLEMENTATION TOTHREE-DIMENSIONAL ATOMSA. Spherical-FEDVR basis
The present implementation is based on our TD-MCSCF code [15], which uses spherical-FEDVR basisfunctions ψ klm ( r ) = 1 r α k ( r ) Y lm (Ω) , (30)where r and Ω are the radial and angular coordinateof r , respectively, Y lm are spherical harmonics, and α k are radial FEDVR basis functions [40, 41]. The ra-dial coordinate of the simulation box [0 , R max ] is dividedinto K FE finite elements. Each finite element supports K DVR local DVR functions, and neighboring elementsare connected by a bridge function. In total, there are K rad = K FE K DVR − ( K FE −
1) radial grid points { r k } ,on which α k ( r k (cid:48) ) = δ kk (cid:48) / √ w k , with { w k } being the inte-gration weights.We expand the channel orbital χ i in the spherical-FEDVR basis as χ i ( r ; t ) = K rad − (cid:88) k =2 L max (cid:88) l =0 ψ klm i ( r ) g kli ( t ) , (31)where L max is the maximum angular momentum in-cluded. The FEDVR basis functions corresponding to r = 0 and r K rad = R max are removed to enforce thevanishing boundary condition for rχ i at both ends of thesimulation box.The electrostatic potentials for electron-electron inter-action, ˆ W φ j φ i ( r ) and W φ j χ i ( r , t ) required for the EOM ofchannel orbitals, are computed by solving Poisson’s equa-tion, ∇ ˆ W φ j φ i ( r ) = − πφ ∗ j ( r ) φ i ( r ) , (32a) ∇ ˆ W φ j χ i ( r , t ) = − πφ ∗ j ( r ) χ i ( r , t ) , (32b) using the method described in Ref. [15]. It should benoted that ˆ W φ j φ i ( r ) is time independent, and Eq. (32a)needs to be solved only once before the simulation. Onthe other hand, ˆ W φ j χ i ( r , t ) is time dependent, and shouldbe computed at every time step. However, since its source φ ∗ j ( r ) χ i ( r , t ) [See Eq. (32b).] and operand { φ j ( r ) } [SeeEq. (20)] are both localized around the atom due to thelocality of occupied orbitals, Eq. (32b) can be solvedwith less computational cost than the similar equationappearing, e.g, in time-dependent Hartree-Fock and TD-MCSCF method [15]. B. Time-propagation with exponential timedifferencing fourth-order Runge-Kutta scheme
For an efficient propagation of the EOM of channel-orbital based TDCIS, we use the exponential time differ-encing fourth-order Runge-Kutta scheme (ETDRK4) byKrogstad [47–49]. To this end, we arrange C and { χ i } into a unified vector χ = ( C , χ ) T and rewrite the EOMsof C and { χ i } by a matrix form i ∂∂t χ = hχ + W [ χ , t ] , (33)where h is a chosen stiff part of the right-hand side of theEOM (See below), and W [ χ , t ] is a nonstiff remainder.We choose the stiff part h to be either (i) the field-freeone-electron Hamiltonian ˆ h or (ii) the totality of theone-electron Hamiltonian ˆ h + ˆ h ext ( t ). For the first case(i) with time-independent h , the time propagation from χ n = χ ( t n ) to χ n +1 = χ ( t n + ∆ t ) is given by χ n +1 = ϕ ( − i h ∆ t ) χ n − i ∆ t [ f ( − i h ∆ t ) W n + f ( − i h ∆ t )( W a + W b ) + f ( − i h ∆ t ) W c ] , (34)where f , f , and f are defined as f ( z ) = ϕ ( z ) − ϕ ( z ) + 4 ϕ ( z ) (35a) f ( z ) = 2 ϕ ( z ) − ϕ ( z ) (35b) f ( z ) = − ϕ ( z ) + 4 ϕ ( z ) , (35c)where z = − i h ∆ t , ϕ ( z ) = e z , and ϕ k +1 ( z ) = 1 z (cid:18) ϕ l ( z ) − k ! (cid:19) ( k = 0 , , , · · · ) . (36) W n , W a , W b , and W c are given by W n = W [ χ n , t n ] (37a) W a = W [ a n , t n + ∆ t/
2] (37b) W b = W [ b n , t n + ∆ t/
2] (37c) W c = W [ c n , t n +1 ] , (37d)where a n , b n , and c n are intermediate vectors given as a n = ϕ ( z/ χ n − i ∆ tϕ ( z/ W n / b n = ϕ ( z/ χ n − i ∆ tϕ ( z/ W n / − i ∆ tϕ ( z/ W a − W n ) (38b) c n = ϕ ( z ) χ n − i ∆ tϕ ( z ) W n − i ∆ tϕ ( z )( W b − W n )(38c)The operator exponential ϕ ( z ) and ϕ ( z/
2) inthe spherical-FEDVR basis are approximated by the Pad´e (3/3) approximation, and higher-order ϕ k functionsare obtained by successively applying Eq. (36). The de-nominator of the Pad´e approximation is factorized andoperated by the matrix iteration method [15]. We fol-low Ref. [48] for the modification required for a time-dependent stiff part h . In the absence of linear part for C , time-propagation of C reduces to well-known fourth-order Runge-Kutta scheme. C. Expectation value
The expectation value of one-body operator (cid:104) O (cid:105) = (cid:104) Ψ | O | Ψ (cid:105) can be evaluated in the LG case as [26, 39] (cid:104) O (cid:105) = occ (cid:88) i { (cid:104) φ i | O | φ i (cid:105) + (cid:104) χ i | O | χ i (cid:105)} + 2 √ (cid:34) C (cid:88) i (cid:104) χ i | O | φ i (cid:105) (cid:35) − occ (cid:88) ij (cid:104) χ i | χ j (cid:105) (cid:104) φ j | O | φ i (cid:105) . (39)The VG expression is obtained by simply replacing { C , χ j } with { D , η j } , and the rVG one by replacing { φ j , χ j } with (cid:8) φ (cid:48) j , χ (cid:48) j (cid:9) . The Ehrenfest theorem ddt (cid:104) O (cid:105) = − i (cid:104) Ψ | [ O, ˆ H ] | Ψ (cid:105) does not hold for TDCIS. Instead, we evaluate thetime derivative of (cid:104) O (cid:105) as [39],˙ (cid:104) O (cid:105) ≡ d (cid:104) O (cid:105) dt = 2 Re (cid:34) occ (cid:88) i (cid:110) (cid:104) χ i | O | ˙ χ i (cid:105) + √ C (cid:104) χ i | O | φ i (cid:105) + √ C (cid:104) ˙ χ i | O | φ i (cid:105) (cid:111)(cid:35) (40)in the LG case. { C , χ j } is to be replaced with { D , η j } for VG. The rVG case needs extra terms [39]:˙ (cid:104) O (cid:105) =2 Re occ (cid:88) i (cid:110) (cid:104) χ (cid:48) i | O | ˙ χ (cid:48) i (cid:105) + √ C (cid:104) χ (cid:48) i | O | φ (cid:48) i (cid:105) + √ C (cid:104) ˙ χ (cid:48) i | O | φ (cid:48) i (cid:105) (cid:111) − occ (cid:88) ij (cid:10) χ (cid:48) i (cid:12)(cid:12) ˙ χ (cid:48) j (cid:11) (cid:10) φ (cid:48) j (cid:12)(cid:12) O | φ (cid:48) i (cid:105) − √ (cid:34) C (cid:88) i (cid:110) E · (cid:104) χ (cid:48) i | O ˆ r | φ (cid:48) i (cid:105) + | A | (cid:104) χ (cid:48) i | O | φ (cid:48) i (cid:105) (cid:111)(cid:35) − i E · occ (cid:88) ij (cid:0) δ ij − (cid:10) χ (cid:48) i (cid:12)(cid:12) χ (cid:48) j (cid:11)(cid:1) (cid:10) φ (cid:48) j (cid:12)(cid:12) [ ˆ r , O ] | φ (cid:48) i (cid:105) . (41)Equations (39), (40), and (41) are valid not only for atoms but also any multielectron system including molecules. D. Ionization probability
To conveniently analyze how ionization proceeds us-ing the TD-MCSCF wavefunctions with time-varying or-bitals, we have previously introduced [11] a domain-based n -fold ionization probability P n , defined as a probabilityto find n electrons in the outer region | r | > R ion andthe other N − n electrons in the inner region | r | < R ion with a given distance R ion from the origin. This quantityis gauge-invariant even during the pulse, unlike the pop-ulation of the (field-free) continuum levels [45, 46]. Toapply this approach to the TDCIS method with chan-nel orbitals, it is reasonable to assume that the occupiedorbitals { φ i } are localized inside the inner region, i.e, φ i ( r ) = 0 for | r | > R ion . Then, the yield of the neutral species, or the survival probability P is computed as, P ( t ) = | C ( t ) | + occ (cid:88) i (cid:104) χ i | χ i (cid:105) < , (42)where (cid:104) χ i | χ j (cid:105) < is the overlap integral in the inner region (cid:104) χ i | χ j (cid:105) < ≡ (cid:90) | r | We present numerical applications of the implemen-tation of the reformulated TDCIS method described inthe previous section and assess efficiency of the rVG. Inall simulations reported below, we assume a laser fieldlinearly polarized along the z axis of the following form: E ( t ) = (cid:112) I sin( ωt ) sin (cid:18) π tN opt T (cid:19) (0 ≤ t ≤ N opt T ) , (44)where I is the peak intensity, ω is the central frequency, T = 2 π/ω is the period, and N opt is the total number ofoptical cycles. A. Helium First, we consider helium atom exposed to a laser pulsewith I = 4 . × W/cm , λ = 400 nm, and N opt = 12.In this condition, exact numerical solution of the TDSEis available [50, 51], from which the expectation value ofdipole velocity and dipole acceleration can be calculatedby using the Ehrenfest theorem. For the TDCIS method,we apply O = ˆ z in Eqs. (39) and (40) to evaluate theexpectation value of dipole moment and velocity, respec-tively. Dipole acceleration is computed by numericallydifferentiating dipole velocity.Time evolution of the calculated dipole moment, dipolevelocity, and dipole acceleration are shown in Fig. 1, andHHG spectra obtained as the modulus squared of theFourier transform of the dipole acceleration is presentedin Fig 2. In these figures, one can see the perfect agree-ment between the LG and rVG results, which numericallyconfirms the gauge invariance between the two gauges.In contrast, the results of conventional VG with fixed or-bitals strongly deviate from them. It should be notedthat, from the comparison between LG (and rVG) andVG results alone, we cannot a priori tell which is moreaccurate. The comparison with the TDSE results now re-veals that the former reproduces the TDSE results muchbetter than the latter, which convinces us of an empirical preference of the LG and rVG treatments.We show the temporal evolution of the survival prob-ability P with R ion = 20 a.u. in Fig. 3 (See Appendixfor the R ion dependence of P ). The conventional VGtreatment strongly overestimates ionization. Recallingthat tunneling ionization is the first process of the three-step model [52, 53], this explains the substantial over-estimation of the HHG yield in Fig. 2. The LG andrVG results, on the other hand, underestimate tunnelingionization. Correspondingly, indeed, we notice that theharmonic intensity is slightly underestimated in Fig. 2. B. Neon We next consider a neon atom subject to a laser fieldwith λ = 800 nm, I = 8 . × W/cm , and N opt = 3and discuss convergence with respect to the maximumangular momentum L max . We show the HHG spectracalculated with various values of L max in the LG and rVG d i po l e m o m en t ( a . u . ) (a) TDSELGrVGVG d i po l e v e l o c i t y ( a . u . ) (b) TDSELGrVGVG0 2 4 6 8 10 12optical cycle0.0080.0060.0040.0020.0000.0020.0040.0060.008 d i po l e a cc e l e r a t i on ( a . u . ) (c) TDSELGrVGVG FIG. 1. Time evolution of (a) the dipole moment, (b) thedipole velocity, and (c) the dipole acceleration of He subjectto a laser pulse with λ = 400 nm , I = 4 × W/cm , and N opt = 12, obtained with the exact TDSE (courtesy of J.Burgd¨orfer) and the TDCIS method with length gauge (LG),conventional velocity gauge (VG), and rotated velocity gauge(rVG). in Fig. 4. Figure 4. (a) shows the equivalence between theLG and rVG for sufficiently large L max (= 100). As can beseen in Fig. 4. (b), which shows LG results, L max = 60 isnot sufficient to obtain a converged result. On the otherhand, the rVG requires far less L max ; even L max = 40 wellreproduces the result with L max = 100, and the spectrumis converged with L max = 60 [Fig. 4. (c)]. This observa-tion indicates that the rVG TDCIS is simultaneously as i n t en s i t y ( a r b . un i t s ) TDSELGrVGVG FIG. 2. HHG spectra from He exposed to a laser pulse withthe same conditions as Fig. 1, computed from the dipoleacceleration shown in Fig. 1. (c). Comparison of the TDSEresult (courtesy of J. Burgd¨orfer) and TDCIS ones with LG,rVG, and VG. p r obab ili t y TDSELGrVGVG FIG. 3. Time evolution of the survival probability P . Com-parison of the TDSE result (courtesy of J. Burgd¨orfer) andTDCIS ones with LG, rVG, and VG with R ion = 20 a.u.. accurate as the LG and as efficient as the VG. V. CONCLUSIONS We have presented a 3D numerical implementation ofthe recently formulated gauge-invariant TDCIS method[39] for atoms subject to a linearly polarized intense laserfield. Compared to the conventional TDCIS method thatuses CI coefficients as working variables, the present im-plementation introduces channel orbitals [26], avoidingcalculation and storage of numerous virtual orbitals. Wehave applied this to He and Ne atoms, and calculatedsurvival probabilities and HHG spectra for intense laserpulses. The perfect agreement of the LG and rVG re-sults obtained with a sufficiently large number of partialwaves numerically demonstrates the gauge invariance ofthe method. The comparison with the numerically ex- i n t en s i t y ( a r b . un i t s ) (a) L max = 100LGrVG0 20 40 60 80 100 12010 i n t en s i t y ( a r b . un i t s ) (b) LG L max = 10060400 20 40 60 80 100 120harmonic order10 i n t en s i t y ( a r b . un i t s ) (c) rVG L max = 1006040 FIG. 4. HHG spectra of Ne subject to an IR laser pulsewith a wavelength of 800 nm and an intensity of 1 . × W/cm (a) Results of the LG and rVG with L max = 100.(b) Results of the LG with various L max . (c) Results of therVG with various L max act TDSE results for He shows the rVG and LG’s su-periority to the conventional VG in terms of accuracy.The VG largely overestimates tunneling ionization and,then, harmonic intensity. The analysis with neon revealsthat the rVG has advantage in computational efficiencyover the LG in terms of the number of spherical harmon-ics required to obtain converged HHG spectrum. Thus,our gauge-invariant reformulation will make TDCIS apromising approach for multielectron dynamics not onlyin atoms but also in molecules driven by high-intensitylaser fields. ACKNOWLEDGMENTS We are grateful to Prof. J. Burgd¨orfer of Vienna Uni-versity of Technology for providing us with He TDSEbenchmark results. This research was supported inpart by a Grant-in-Aid for Scientific Research (GrantsNo. 16H03881, No. 17K05070, No. 18H03891, andNo. 19H00869) from the Ministry of Education, Cul-ture, Sports, Science and Technology (MEXT) of Japan.This research was also partially supported by JST COI(Grant No. JPMJCE1313), JST CREST (Grant No. JP-MJCR15N1), and by Quantum Leap Flagship Programof MEXT. T. T. gratefully acknowledges support fromthe Graduate School of Engineering, The Universityof Tokyo, Doctoral Student Special Incentives Program(SEUT Fellowship). Appendix: R ion dependence on P We show P from the LG TDCIS results with variousvalues of R ion in fig. 5. We can see stepwise ejection ofelectron wavepackets that propagate outwards; the larger R ion is, the later P is depleted. In the cases of R ion =10, 20, and 30 a.u., P nearly reaches approximately thesame final value in the end of the simlation (after twelveoptical cycles). p r obab ili t y R ion = 10 a.u.20 a.u.30 a.u.40 a.u.50 a.u. FIG. 5. Temporal evolution of P for various values of R ion ,extracted from LG TDCIS results.[1] F. Krausz and M. Ivanov, Rev. Mod. Phys. , 163(2009).[2] M. Nisoli, P. Decleva, F. Calegari, A. Palacios, andF. Mart´ın, Chem. Rev. , 10760 (2017).[3] L. Gallmann, C. Cirelli, and U. Keller, Annu. Rev. Phys.Chem. , 447 (2013).[4] J. Zanghellini, M. Kitzler, C. Fabian, T. Brabec, andA. Scrinzi, Laser Phys. , 1064 (2003).[5] T. Kato and H. Kono, Chem. Phys. Lett. , 533 (2004).[6] D. J. Haxton and C. W. McCurdy, Phys. Rev. A ,012509 (2015).[7] H. Miyagi and L. B. Madsen, Phys. Rev. A , 063416(2014).[8] J. Caillat, J. Zanghellini, M. Kitzler, O. Koch,W. Kreuzer, and A. Scrinzi, Phys. Rev. A , 012712(2005).[9] H. Miyagi and L. B. Madsen, Phys. Rev. A , 062511(2013).[10] H. Miyagi and L. Madsen, Phys. Rev. A , 062511(2013).[11] T. Sato and K. L. Ishikawa, Phys. Rev. A , 023402(2013).[12] H. Miyagi and L. B. Madsen, Phys. Rev. A , 063416(2014).[13] D. J. Haxton and C. W. McCurdy, Phys. Rev. A ,012509 (2015).[14] T. Sato and K. L. Ishikawa, Phys. Rev. A , 023417(2015). [15] T. Sato, K. L. Ishikawa, I. Bˇrezinov´a, F. Lackner,S. Nagele, and J. Burgd¨orfer, Phys. Rev. A , 023405(2016).[16] R. Anzaki, T. Sato, and K. L. Ishikawa, Phys. Chem.Chem. Phys. , 22008 (2017).[17] S. Kvaal, J. Chem. Phys. , 194109 (2012).[18] T. Sato, H. Pathak, Y. Orimo, and K. L. Ishikawa, J.Chem. Phys. , 051101 (2018).[19] M. A. Lysaght, H. W. van der Hart, and P. G. Burke,Phys. Rev. A , 053411 (2009).[20] M. A. Lysaght, P. G. Burke, and H. W. van der Hart,Phys. Rev. Lett. , 253001 (2008).[21] P. G. Burke and V. M. Burke, J. Phys. B: At., Mol. Opt.Phys. , L383 (1997).[22] L. R. Moore, M. A. Lysaght, L. A. Nikolopoulos, J. S.Parker, H. W. Van Der Hart, and K. T. Taylor, J. Mod.Opt. , 1132 (2011).[23] D. D. Clarke, G. S. Armstrong, A. C. Brown, and H. W.Van Der Hart, Phys. Rev. A , 1 (2018).[24] F. Lackner, I. Bˇrezinov´a, T. Sato, K. L. Ishikawa, andJ. Burgd¨orfer, Phys. Rev. A , 023412 (2015).[25] F. Lackner, I. Bˇrezinov´a, T. Sato, K. L. Ishikawa, andJ. Burgd¨orfer, Phys. Rev. A , 033414 (2017).[26] N. Rohringer, A. Gordon, and R. Santra, Phys. Rev. A , 043420 (2006).[27] A. Gordon, F. X. K¨artner, N. Rohringer, and R. Santra,Phys. Rev. Lett. , 223902 (2006).[28] A. Karamatskou and R. Santra, Phys. Rev. A , 013415(2017). [29] S. Pabst, L. Greenman, P. J. Ho, D. A. Mazziotti, andR. Santra, Phys. Rev. Lett. , 053003 (2011).[30] L. Greenman, P. J. Ho, S. Pabst, E. Kamarchik, D. A.Mazziotti, and R. Santra, Phys. Rev. A , 023406(2010).[31] M. Grosser, J. M. Slowik, and R. Santra, Phys. Rev. A , 062107 (2017).[32] S. Pabst and R. Santra, J. Phys. B: At., Mol. Opt. Phys. , 124026 (2014).[33] A. Sytcheva, S. Pabst, S.-K. Son, and R. Santra, Phys.Rev. A , 023414 (2012).[34] S. Pabst, A. Sytcheva, A. Moulet, A. Wirth, E. Gouliel-makis, and R. Santra, Phys. Rev. A , 063411 (2012).[35] S. Pabst and R. Santra, Phys. Rev. Lett. , 233005(2013).[36] E. Heinrich-Josties, S. Pabst, and R. Santra, Phys. Rev.A , 043415 (2014).[37] J. A. You, J. M. Dahlstr¨om, and N. Rohringer, Phys.Rev. A , 023409 (2017).[38] J. A. You, N. Rohringer, and J. M. Dahlstr¨om, Phys.Rev. A , 033413 (2016).[39] T. Sato, T. Teramura, and K. L. Ishikawa, Appl. Sci ,433 (2018).[40] T. N. Rescigno and C. W. McCurdy, Phys. Rev. A ,032706 (2000).[41] C. W. McCurdy, M. Baertschy, and T. N. Rescigno, J.Phys. B: At., Mol. Opt. Phys. , R137 (2004).[42] B. I. Schneider, L. A. Collins, and S. X. Hu, Phys. Rev.E , 036708 (2006). [43] B. I. Schneider, F. Johannes, S. Nagele, R. Pazourek,S. Hu, L. A. Collins, and J. Burgd¨orfer, Quantum Dy-namic Imaging , edited by A. D. Bandrauk and M. Ivanov(Springer New York, New York, NY, 2011) p. 149.[44] P. Kramer and M. Saraceno,