Implementation of gauge-invariant time-dependent configuration interaction singles method for three-dimensional atoms

Takuma Teramura1,∗ Takeshi Sato1,2,3,† and Kenichi L. Ishikawa1,2,3‡ 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: May 23, 2019)

Among them, the time-dependent configuration interaction singles (TDCIS) method is one of the promissing methods [26][27][28][29][30][31][32][33][34][35][36]. This method has been successfully applied to various electron dynamics such as giant enhancement in HHG in Xe [35] and decoherence in attosecond photoionization [27]. In the TDCIS method, the total electronic wavefunction is approximated by a superposition of time-independent Slater determinants where |Φ 0 is the Hartree-Fock (HF) ground state and |Φ a i is a singly-excited configuration replacing an occupied orbital φ i with a virtual orbital φ a unoccupied in the ground state. The orbital functions are fixed and propagation of configuration-interaction (CI) coefficients (C 0 and{C a i }) describes the system dynamics. Although applications of the TDCIS method are limited to the single excitation or ionization due to the truncation of CI space, its low computational cost and ease of analysis are attractive.
The conventional TDCIS method with CI coefficients has two major issues; the explicit calculation and storage of virtual orbitals {φ a } and a violation of gauge invariance. Virtual orbitals {φ a } should include both bound and continuum orbitals, whose number is infinite in principle. In a practical simulation with real-space grids, one has to prepare virtual orbitals in advance by numerically obtaining all eigenstates of the discretized HF equation. The number of the virtual orbitals increases with the number of the grid points. Thus, the calculation and storage of virtual orbitals become unacceptably demanding for systems with a large number of grid points like molecules. To solve this problem, Rohringer et.al. have proposed an alternative but equivalent formulation of the TDCIS method in which time-dependent orbitallike quantity called channel orbital is propagated instead of CI coefficient [26]. The channel orbital is defined by using {φ a } and C a i as virtual orbitals, while in princple including all the virtual orbitals within a grid space, and has enhanced utility of the TDCIS method. However, to the best of our knowledge, 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 with Hartree-Slater potential in Ref. [27]. The TDCIS method, either CI coefficient-based or channel-orbital-based, suffers from a violation of gauge invariance, as a general consequence of the truncation of CI 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 recently reported a gauge-invariant reformulation of the TDCIS method [39]. In our reformulation, a rotated velocity gauge (rVG) transformed from the length gauge (LG) by a unitary operator has been introduced. This unitary transformation ensures the gauge invariance between the LG and rVG, and Ref. [39] numerically confirmed the equivalence of these gauges for a model 1D Hamiltonian.
In this paper, we report a three dimensional numerical implementation of the gauge-invariant TDCIS method for atoms subject to a linearly polarized laser pulse. We employ a spherical harmonics expansion of orbital functions with the radial coordinate discretized by a finiteelement discrete variable representation (FEDVR) [40][41][42][43]. We apply the present implementation to HHG from helium and neon atoms and asses the advantage of the rVG over the LG and VG. This paper proceeds as follows. In Sec. II, we briefly review the TDCIS methods. The numerical implementation of the gauge-invariant TDCIS method to threedimensional atoms is given in Sec. III. We describe numerical applications in Sec. IV and conclude this work in Sec. V. We use Hartree atomic units (a.u.) throughout the paper unless otherwise stated.

A. The System Hamiltonian and Gauge transformation
We consider an atom with N electrons with a nucleus located at the origin. The time evolution of the Nelectron wavefunction |Ψ(t) is governed by the TDSE, whereĤ(t) is the time-dependent non-relativistic Hamil-tonianĤ decomposed into the field-free part, and the laser-electron interaction part In these expressions, r i and p i = −i∇ i are the position and the canonical momentum of the electron i, respectively.ĥ 0 is given by, where Z is the atomic number. Within the electric dipole approximation,ĥ ext for the LG and VG are given bŷ where E(t) and A(t) = − t −∞ dt E(t ) are the electric field and the vector potential of the external laser field, respectively.
The two gauges are physically equivalent, and any physical observable takes the same value, independent of the choice of the gauge. The LG wavefunction Ψ LG and VG wavefunction Ψ VG (t) are mutually transformed by a gauge transformation as, B. The CI coefficient-based TDCIS method in the length gauge In the conventional TDCIS method based on CI coefficients, orbitals satisfy the canonical, restricted HF equationf wheref is the Fock operator andŴ φ φ is the potential due to the product of two given orbitals φ and φ , defined in the real space aŝ p is the orbital energy of orbital φ p . In the TDCIS wavefunction in Eq. (1), |Φ 0 is the HF ground state formed with the occupied orbitals as whereĉ † pσ andĉ pσ are the creation and annihilation operators, respectively, of spin-orbital φ p ⊗ σ, and | is the vacuum. σ ∈ {↑, ↓} denotes the spin function. |Φ a i is a singly-excited configuration replacing an occupied orbital φ i with a virtual orbital φ a The EOMs of CI coefficients is derived through the Dirac-Frenkel time-dependent variational principle [44], requiring Lagrangian L(t) to be stationary with respect to the variation of C * 0 and {C a * i }. E 0 = Φ 0 |Ĥ 0 |Φ 0 denotes the HF energy. This constant shift, introduced for the simple notation of the EOMs, does not affect physical results. In the LG case, in which the wavefuntion is written as, the EOMs of CI coefficients are obtained as C. The channel orbital-based TDCIS method in the length gauge The EOMs of CI coefficients [Eq. (17)] can be rewritten, by substituting channel orbital Eq. (2) into Eq. (17), as, wherê andP is the projection operator onto the space spanned by virtual orbitalŝ with1 being the identity operator.

D. Velocity gauge and rotated velocity gauge
One can, in principle, derive the EOMs for the VG case in the same way as for the LG. Let us write the total wavefunction and channel orbital as, |Φ 0 and |Φ a i are the same configurations as those used in the LG case. Their EOMs are obtained through the same procedures as in the LG case as, It is known that TDCIS, which uses time-independent orbitals, is not gauge-invariant [39,45,46]. Instead of the conventional VG as described above, we have recently proposed the rVG [39], where we define the rVG wavefunction by the gauge transformation from the LG wavefunction as, The rVG orbitals are related to the LG ones by, wherê They satisfy the following EOMs [39]: whereP andF i are given by Eqs. (21) and (20), respectively, with {φ j } replaced with φ j . Although Eq. (29) contains the dipole operator E · r, this does not prevent enjoying the advantages of the VG treatment, since it acts only on the localized occupied orbitals {φ i }.

A. Spherical-FEDVR basis
The present implementation is based on our TD-MCSCF code [15], which uses spherical-FEDVR basis functions where r and Ω are the radial and angular coordinate of r, respectively, Y lm are spherical harmonics, and α k are radial FEDVR basis functions [40,41]. The radial coordinate of the simulation box [0, R max ] is divided into K FE finite elements. Each finite element supports K DVR local DVR functions, and neighboring elements are 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 ) = δ kk / √ w k , with {w k } being the integration weights. We expand the channel orbital χ i in the spherical-FEDVR basis as where L max is the maximum angular momentum included. The FEDVR basis functions corresponding to r 1 = 0 and r K rad = R max are removed to enforce the vanishing boundary condition for rχ i at both ends of the simulation box. The electrostatic potentials for electron-electron interaction,Ŵ φj φi (r) and W φj χi (r, t) required for the EOM of channel orbitals, are computed by solving Poisson's equation, using the method described in Ref. [15]. It should be noted thatŴ φj φi (r) is time independent, and Eq. (32a) needs to be solved only once before the simulation. On the other hand,Ŵ φj χi (r, t) is time dependent, and should be computed at every time step. However, since its source φ * j (r)χ i (r, t) [See Eq. (32b).] and operand {φ j (r)} [See Eq. (20)] are both localized around the atom due to the locality of occupied orbitals, Eq. (32b) can be solved with less computational cost than the similar equation appearing, e.g, in time-dependent Hartree-Fock and TD-MCSCF method [15].

B. Time-propagation with exponential time differencing fourth-order Runge-Kutta scheme
For an efficient propagation of the EOM of channelorbital based TDCIS, we use the exponential time differencing fourth-order Runge-Kutta scheme (ETDRK4) by Krogstad [47][48][49]. To this end, we arrange C 0 and {χ i } into a unified vector χ = (C 0 , χ) T and rewrite the EOMs of C 0 and {χ i } by a matrix form where h is a chosen stiff part of the right-hand side of the EOM (See below), and W [χ, t] is a nonstiff remainder. We choose the stiff part h to be either (i) the field-free one-electron Hamiltonianĥ 0 or (ii) the totality of the one-electron Hamiltonianĥ 0 +ĥ 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 where f 1 , f 2 , and f 3 are defined as where z = −ih∆t, ϕ 0 (z) = e z , and 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) where a n , b n , and c n are intermediate vectors given as a n = ϕ 0 (z/2)χ n − i∆tϕ 1 (z/2)W n /2 (38a) b n = ϕ 0 (z/2)χ n − i∆tϕ 1 (z/2)W n /2 The operator exponential ϕ 0 (z) and ϕ 0 (z/2) in the spherical-FEDVR basis are approximated by the Padé (3/3) approximation, and higher-order ϕ k functions are obtained by successively applying Eq. (36). The denominator of the Padé approximation is factorized and operated by the matrix iteration method [15]. We follow Ref. [48] for the modification required for a timedependent stiff part h. In the absence of linear part for C 0 , time-propagation of C 0 reduces to well-known fourthorder Runge-Kutta scheme.

C. Expectation value
The expectation value of one-body operator O = Ψ| O |Ψ can be evaluated in the LG case as [26,39] The VG expression is obtained by simply replacing {C 0 , χ j } with {D 0 , η j }, and the rVG one by replacing {φ j , χ j } with φ j , χ j .
in the LG case. {C 0 , χ j } is to be replaced with {D 0 , η j } for VG. The rVG case needs extra terms [39]: 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 using the TD-MCSCF wavefunctions with time-varying orbitals, we have previously introduced [11] a domain-based n-fold ionization probability P n , defined as a probability to find n electrons in the outer region |r| > R ion and the other N − n electrons in the inner region |r| < R ion with a given distance R ion from the origin. This quantity is gauge-invariant even during the pulse, unlike the population of the (field-free) continuum levels [45,46]. To apply this approach to the TDCIS method with channel orbitals, it is reasonable to assume that the occupied orbitals {φ 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 0 is computed as, where χ i |χ j < is the overlap integral in the inner region Noting that the atom described by a TDCIS wavefunction is at most singly ionized, we obtain the singleionization probability as P 1 (t) = 1 − P 0 (t).

IV. NUMERICAL EXAMPLE
We present numerical applications of the implementation of the reformulated TDCIS method described in the previous section and assess efficiency of the rVG. In all simulations reported below, we assume a laser field linearly polarized along the z axis of the following form: where I 0 is the peak intensity, ω is the central frequency, T = 2π/ω is the period, and N opt is the total number of optical cycles.

A. Helium
First, we consider helium atom exposed to a laser pulse with I 0 = 4.0×10 14 W/cm 2 , λ = 400 nm, and N opt = 12. In this condition, exact numerical solution of the TDSE is available [50,51], from which the expectation value of dipole velocity and dipole acceleration can be calculated by using the Ehrenfest theorem. For the TDCIS method, we apply O =ẑ in Eqs. (39) and (40) to evaluate the expectation value of dipole moment and velocity, respectively. Dipole acceleration is computed by numerically differentiating dipole velocity.
Time evolution of the calculated dipole moment, dipole velocity, and dipole acceleration are shown in Fig. 1, and HHG spectra obtained as the modulus squared of the Fourier transform of the dipole acceleration is presented in Fig 2. In these figures, one can see the perfect agreement between the LG and rVG results, which numerically confirms the gauge invariance between the two gauges. In contrast, the results of conventional VG with fixed orbitals strongly deviate from them. It should be noted that, from the comparison between LG (and rVG) and VG results alone, we cannot a priori tell which is more accurate. The comparison with the TDSE results now reveals that the former reproduces the TDSE results much better than the latter, which convinces us of an empirical preference of the LG and rVG treatments.
We show the temporal evolution of the survival probability P 0 with R ion = 20 a.u. in Fig. 3 (See Appendix for the R ion dependence of P 0 ). The conventional VG treatment strongly overestimates ionization. Recalling that tunneling ionization is the first process of the threestep model [52,53], this explains the substantial overestimation of the HHG yield in Fig. 2. The LG and rVG results, on the other hand, underestimate tunneling ionization. Correspondingly, indeed, we notice that the harmonic intensity is slightly underestimated in Fig. 2.

B. Neon
We next consider a neon atom subject to a laser field with λ = 800 nm, I 0 = 8.0 × 10 14 W/cm 2 , and N opt = 3 and discuss convergence with respect to the maximum angular momentum L max . We show the HHG spectra calculated with various values of L max in the LG and rVG in Fig. 4. Figure 4. (a) shows the equivalence between the LG and rVG for sufficiently large L max (= 100). As can be seen in Fig. 4. (b), which shows LG results, L max = 60 is not sufficient to obtain a converged result. On the other hand, the rVG requires far less L max ; even L max = 40 well reproduces the result with L max = 100, and the spectrum is converged with L max = 60 [ Fig. 4. (c)]. This observation indicates that the rVG TDCIS is simultaneously as LG rVG VG FIG. 2. HHG spectra from He exposed to a laser pulse with the same conditions as Fig. 1, computed from the dipole acceleration shown in Fig. 1. (c). Comparison of the TDSE result (courtesy of J. Burgdörfer) and TDCIS ones with LG, rVG, and VG. accurate as the LG and as efficient as the VG.

V. CONCLUSIONS
We have presented a 3D numerical implementation of the recently formulated gauge-invariant TDCIS method [39] for atoms subject to a linearly polarized intense laser field. Compared to the conventional TDCIS method that uses CI coefficients as working variables, the present implementation introduces channel orbitals [26], avoiding calculation and storage of numerous virtual orbitals. We have applied this to He and Ne atoms, and calculated survival probabilities and HHG spectra for intense laser pulses. The perfect agreement of the LG and rVG results obtained with a sufficiently large number of partial waves numerically demonstrates the gauge invariance of the method. The comparison with the numerically ex- act TDSE results for He shows the rVG and LG's superiority to the conventional VG in terms of accuracy. The VG largely overestimates tunneling ionization and, then, harmonic intensity. The analysis with neon reveals that the rVG has advantage in computational efficiency over the LG in terms of the number of spherical harmonics required to obtain converged HHG spectrum. Thus, our gauge-invariant reformulation will make TDCIS a promising approach for multielectron dynamics not only in atoms but also in molecules driven by high-intensity laser fields. Appendix: R ion dependence on P0 We show P 0 from the LG TDCIS results with various values of R ion in fig. 5. We can see stepwise ejection of electron wavepackets that propagate outwards; the larger R ion is, the later P 0 is depleted. In the cases of R ion = 10, 20, and 30 a.u., P 0 nearly reaches approximately the same final value in the end of the simlation (after twelve optical cycles).