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

TL;DR

This paper introduces a gauge-invariant, computationally efficient implementation of the time-dependent configuration interaction singles method for three-dimensional atoms, enabling improved simulations of high-field phenomena like high-harmonic generation.

Contribution

It presents a novel gauge-invariant TDCIS implementation using channel orbitals, enhancing computational efficiency and physical consistency in atomic high-field simulations.

Findings

Confirmed gauge invariance in simulations

Demonstrated effectiveness of velocity gauge treatment

Applied to high-harmonic generation in helium and neon

Abstract

We present a numerical implementation of the gauge-invariant time-dependent configuration interaction singles (TDCIS) method [Appl. Sci. 8, 433 (2018)] for three-dimensional atoms. In our implementation, orbital-like quantity called channel orbital [Phys. Rev. A 74, 043420 (2006)] is propagated instead of configuration-interaction (CI) coefficients, which removes a computational bottleneck of explicitly calculating and storing numerous virtual orbitals. Furthermore, besides its physical consistency, the gauge-invariant formulation allows to take advantages of the velocity gauge treatment of the laser-electron interaction over the length gauge one in the simulation of high-field phenomena. We apply the present implementation to high-harmonic generation from helium and neon atoms, and numerically confirms the gauge invariance and demonstrates the effectiveness of the rotated velocity gauge treatment.

Figures

Equations

$$|\Psi (t)⟩=|{\Phi }_0⟩C_0(t)+\sum _i^{\text{occ}}\sum _a^{\text{vir}}|{\Phi }_i^a⟩C_i^a(t),$$

$$|{\chi }_i(t)⟩\equiv \sum _a^{\text{vir}}|{\varphi }_a⟩C_i^a(t).$$

$$i\text{partialderivative}t|\Psi (t)⟩=\hat{H}(t)|\Psi (t)⟩,$$

$$\hat{H}(t)={\hat{H}}_0+{\hat{H}}_{\text{ext}}(t),$$

$${\hat{H}}_0=\sum _{i=1}^N{\hat{h}}_0(r_i,p_i)+\sum _{i>j}^N\frac{1}{\text{absolutevalue}{r}_i-r_j}$$

$${\hat{H}}_{\text{ext}}(t)=\sum _{i=1}^N{\hat{h}}_{\text{ext}}(r_i,p_i,t),$$

$${\hat{h}}_0(r,p)=\frac{p^2}{2}-\frac{Z}{\text{absolutevalue}r},$$

$${\hat{h}}_{\text{ext}}^{\text{LG}}(r,p,t)$$

$${\hat{h}}_{\text{ext}}^{\text{VG}}(r,p,t)$$

$$|{\Psi }^{\text{VG}}(t)⟩=\hat{U}(t)|{\Psi }^{\text{LG}}(t)⟩$$

$$\hat{U}(t)\equiv \exp \text{quantity}[-i\sum _{i=1}^N\text{quantity}(A(t)\cdot r_i-\frac{1}{2}{\int }_{-\infty }^{t}d{t}^{\prime }\text{absolutevalue}{A(t^{\prime })}^2)]$$

$$\hat{f}|{\varphi }_p⟩$$

$$={ϵ}_p|{\varphi }_p⟩,$$

$${\hat{W}}_{{\varphi }^{\prime }}^{\varphi }(r_1)\equiv \int d{r}_2\frac{{\varphi }^{\ast }(r_2){\varphi }^{\prime }(r_2)}{\text{absolutevalue}{r}_1-r_2}.$$

$$|{\Phi }_0⟩=\prod _i^{\text{occ}}{\hat{c}}_{i\uparrow }^{†}{\hat{c}}_{i\downarrow }^{†}|⟩,$$

$$|{\Phi }_i^a⟩=\frac{1}{\sqrt{2}}({\hat{c}}_{a\uparrow }^{†}{\hat{c}}_{i\uparrow }+{\hat{c}}_{a\downarrow }^{†}{\hat{c}}_{i\downarrow })|{\Phi }_0⟩.$$

$$L(t)=⟨\Psi |\hat{H}(t)-i\text{partialderivative}t-E_0|\Psi ⟩$$

$$|{\Psi }^{\text{LG}}(t)⟩=|{\Phi }_0⟩C_0+\sum _i^{\text{occ}}\sum _a^{\text{vir}}|{\Phi }_i^a⟩C_i^a.$$

$$i{\dot{C}}_0$$

$$i{\dot{C}}_i^a$$

$$+\sqrt{2}E\cdot r|{\varphi }_i⟩C_0\}-E\cdot \sum _j^{\text{occ}}⟨{\varphi }_j|r|{\varphi }_i⟩C_j^a,$$

$${\hat{F}}_{ij}|{\varphi }_b⟩={\delta }_j^i(\hat{f}-{ϵ}_i)|{\varphi }_b⟩+2{\hat{W}}_{{\varphi }_b}^{{\varphi }_j}|{\varphi }_i⟩-{\hat{W}}_{{\varphi }_i}^{{\varphi }_j}|{\varphi }_b⟩.$$

$$i{\dot{C}}_0$$

$$i|{\dot{\chi }}_i⟩$$

$$-E\cdot \sum _j^{\text{occ}}|{\chi }_j⟩⟨{\varphi }_j|r|{\varphi }_i⟩,$$

$$\hat{F}|{\chi }_i⟩=(\hat{f}-{ϵ}_i)|{\chi }_i⟩+\sum _j^{\text{occ}}(2{\hat{W}}_{{\chi }_j}^{{\varphi }_j}|{\varphi }_i⟩-{\hat{W}}_{{\varphi }_i}^{{\varphi }_j}|{\chi }_j⟩),$$

$$\hat{P}=\sum _a^{\text{vir}}|{\varphi }_a⟩⟨{\varphi }_a|=\hat{1}-\sum _j^{\text{occ}}|{\varphi }_j⟩⟨{\varphi }_j|,$$

$$|{\Psi }^{\text{VG}}(t)⟩$$

$$|\eta (t)⟩$$

$$i{\dot{D}}_0$$

$$i|{\dot{\eta }}_i⟩$$

$$-A\cdot \sum _j^{\text{occ}}|{\eta }_j⟩⟨{\varphi }_j|p|{\varphi }_i⟩.$$

$$|{\Psi }^{\text{rVG}}(t)⟩=\hat{U}(t)|{\Psi }^{\text{LG}}(t)⟩$$

$$|{\varphi }_p^{\prime }(t)⟩$$

$$|{\chi }_i^{\prime }(t)⟩$$

$$\hat{u}(t)=\exp -i\text{quantity}(A(t)\cdot r-\frac{1}{2}{\int }_{-\infty }^{t}d{t}^{\prime }\text{absolutevalue}{A(t^{\prime })}^2).$$

Full text

††preprint: APS/123-QED

I Introduction

Recent laser technology capable of generating strong laser pulses with an intensity $\gtrsim 10^{14}\,{\rm W/cm}^{2}$ has enabled us to explore electron dynamics in nonperturbative regime, e.g., high-harmonic generation (HHG), above threshold ionization, nonsequential double ionization, and attosecond pulse generation Krausz and Ivanov (2009); Nisoli et al. (2017); Gallmann et al. (2013). While laser-driven electron dynamics is rigorously described by the time-dependent Schrödinger equation (TDSE), its direct numerical solution is practically unfeasible for systems with more than two electrons. For theoretical investigation of multielectron dynamics in intense laser field, various tractable ab-initio methods have been developed, e.g., time-dependent multiconfiguration self-consistent field (TD-MCSCF) methods Zanghellini et al. (2003); Kato and Kono (2004); Haxton and McCurdy (2015a); Miyagi and Madsen (2014a); Caillat et al. (2005); Miyagi and Madsen (2013a, b); Sato and Ishikawa (2013); Miyagi and Madsen (2014b); Haxton and McCurdy (2015b); Sato and Ishikawa (2015); Sato et al. (2016); Anzaki et al. (2017), time-dependent coupled cluster method Kvaal (2012); Sato et al. (2018a), time-dependent $R-$matrix approach Lysaght et al. (2009, 2008); Burke and Burke (1997); Moore et al. (2011); Clarke et al. (2018), and time-dependent reduced two body density matrix approachLackner et al. (2015, 2017).

Among them, the time-dependent configuration interaction singles (TDCIS) method is one of the promissing methods Rohringer et al. (2006); Gordon et al. (2006); Karamatskou and Santra (2017); Pabst et al. (2011); Greenman et al. (2010); Grosser et al. (2017); Pabst and Santra (2014); Sytcheva et al. (2012); Pabst et al. (2012); Pabst and Santra (2013); Heinrich-Josties et al. (2014). This method has been successfully applied to various electron dynamics such as giant enhancement in HHG in Xe Pabst and Santra (2013) and decoherence in attosecond photoionization Gordon et al. (2006). In the TDCIS method, the total electronic wavefunction is approximated by a superposition of time-independent Slater determinants

[TABLE]

where $\ket{\Phi_{0}}$ is the Hartree-Fock (HF) ground state and $\ket{\Phi^{a}_{i}}$ is a singly-excited configuration replacing an occupied orbital $\phi_{i}$ with a virtual orbital $\phi_{a}$ unoccupied in the ground state. The orbital functions are fixed and propagation of configuration-interaction (CI) coefficients ($C_{0}$ and$\quantity{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 $\quantity{\phi_{a}}$ and a violation of gauge invariance. Virtual orbitals $\quantity{\phi_{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 orbital-like quantity called channel orbital is propagated instead of CI coefficient Rohringer et al. (2006). The channel orbital is defined by using $\quantity{\phi_{a}}$ and $C^{a}_{i}$ as

[TABLE]

The equations of motions (EOMs) of CI coefficient are converted to those of $C_{0}$ and $\quantity{\chi_{i}}$. This reformulation removes computational bottleneck of handling numerous 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 TDCIS method have been limited to one-dimensional (1D) model in Refs. Rohringer et al. (2006); You et al. (2017, 2016) and nobles gas atoms with Hartree-Slater potential in Ref. Gordon et al. (2006).

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 TDCIS method. To overcome this difficulty, we have recently reported a gauge-invariant reformulation of the TDCIS method Sato et al. (2018b). 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. Sato et al. (2018b) 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 finite-element discrete variable representation (FEDVR)Rescigno and McCurdy (2000); McCurdy et al. (2004); Schneider et al. (2006, 2011). 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 three-dimensional 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.

II Theory

II.1 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 $N$-electron wavefunction $\ket{\Psi(t)}$ is governed by the TDSE,

[TABLE]

where $\hat{H}(t)$ is the time-dependent non-relativistic Hamiltonian

[TABLE]

decomposed into the field-free part,

[TABLE]

and the laser-electron interaction part

[TABLE]

In these expressions, ${\bf\it r}_{i}$ and ${\bf\it p}_{i}=-i\nabla_{i}$ are the position and the canonical momentum of the electron $i$, respectively. $\hat{h}_{0}$ is given by,

[TABLE]

where $Z$ is the atomic number. Within the electric dipole approximation, $\hat{h}_{\text{ext}}$ for the LG and VG are given by

[TABLE]

where ${\bf\it E}(t)$ and ${\bf\it A}(t)=-\int^{t}_{-\infty}dt^{\prime}{\bf\it E}(t^{\prime})$ 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 $\ket{\Psi^{\text{LG}}}$ and VG wavefunction $\ket{\Psi^{\text{VG}}(t)}$ are mutually transformed by a gauge transformation as,

[TABLE]

II.2 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 equation

[TABLE]

where $\hat{f}$ is the Fock operator and $\hat{W}^{\phi}_{\phi^{\prime}}$ is the potential due to the product of two given orbitals $\phi$ and $\phi^{\prime}$, defined in the real space as

[TABLE]

$\epsilon_{p}$ is the orbital energy of orbital $\phi_{p}$. In the TDCIS wavefunction in Eq. (1), $\ket{\Phi_{0}}$ is the HF ground state formed with the occupied orbitals as

[TABLE]

where $\hat{c}^{\dagger}_{p\sigma}$ and $\hat{c}_{p\sigma}$ are the creation and annihilation operators, respectively, of spin-orbital $\phi_{p}\otimes\sigma$, and $\ket{}$ is the vacuum. $\sigma\in\quantity{\uparrow,\downarrow}$ denotes the spin function. $\ket{\Phi^{a}_{i}}$ is a singly-excited configuration replacing an occupied orbital $\phi_{i}$ with a virtual orbital $\phi_{a}$

[TABLE]

The EOMs of CI coefficients is derived through the Dirac-Frenkel time-dependent variational principle Kramer and Saraceno (1981), requiring Lagrangian $L(t)$

[TABLE]

to be stationary with respect to the variation of $C^{*}_{0}$ and $\quantity{C^{a*}_{i}}$. $E_{0}=\bra{\Phi_{0}}\hat{H}_{0}\ket{\Phi_{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,

[TABLE]

the EOMs of CI coefficients are obtained as

[TABLE]

where

[TABLE]

II.3 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,

[TABLE]

where

[TABLE]

and $\hat{P}$ is the projection operator onto the space spanned by virtual orbitals

[TABLE]

with $\hat{1}$ being the identity operator.

II.4 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,

[TABLE]

$\ket{\Phi_{0}}$ and $\ket{\Phi^{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,

[TABLE]

It is known that TDCIS, which uses time-independent orbitals, is not gauge-invariant Ishikawa and Sato (2015); Sato et al. (2018b, c). Instead of the conventional VG as described above, we have recently proposed the rVG Sato et al. (2018b), where we define the rVG wavefunction by the gauge transformation from the LG wavefunction as,

[TABLE]

The rVG orbitals are related to the LG ones by,

[TABLE]

where

[TABLE]

They satisfy the following EOMs Sato et al. (2018b):

[TABLE]

where $\hat{P}^{\prime}$ and $\hat{F}^{\prime}_{i}$ are given by Eqs.  (21) and (20), respectively, with $\quantity{\phi_{j}}$ replaced with $\quantity{\phi^{\prime}_{j}}$. Although Eq. (29) contains the dipole operator ${\bf\it E}\cdot{\bf\it r}$, this does not prevent enjoying the advantages of the VG treatment, since it acts only on the localized occupied orbitals $\{\phi^{\prime}_{i}\}$.

III Implementation to three-dimensional atoms

III.1 Spherical-FEDVR basis

The present implementation is based on our TD-MCSCF code Sato et al. (2016), which uses spherical-FEDVR basis functions

[TABLE]

where $r$ and $\Omega$ are the radial and angular coordinate of ${\bf\it r}$, respectively, $Y_{lm}$ are spherical harmonics, and $\alpha_{k}$ are radial FEDVR basis functions Rescigno and McCurdy (2000); McCurdy et al. (2004). The radial coordinate of the simulation box $[0,R_{\text{max}}]$ is divided into $K_{\text{FE}}$ finite elements. Each finite element supports $K_{\text{DVR}}$ local DVR functions, and neighboring elements are connected by a bridge function. In total, there are $K_{\text{rad}}=K_{\text{FE}}K_{\text{DVR}}-(K_{\text{FE}}-1)$ radial grid points $\{r_{k}\}$, on which $\alpha_{k}(r_{k^{\prime}})=\delta_{kk^{\prime}}/\sqrt{w_{k}},$ with $\{w_{k}\}$ being the integration weights.

We expand the channel orbital $\chi_{i}$ in the spherical-FEDVR basis as

[TABLE]

where $L_{\text{max}}$ is the maximum angular momentum included. The FEDVR basis functions corresponding to $r_{1}=0$ and $r_{K_{\text{rad}}}=R_{\text{max}}$ are removed to enforce the vanishing boundary condition for $r\chi_{i}$ at both ends of the simulation box.

The electrostatic potentials for electron-electron interaction, $\hat{W}^{\phi_{j}}_{\phi_{i}}({\bf\it r})$ and $W^{\phi_{j}}_{\chi_{i}}({\bf\it r},t)$ required for the EOM of channel orbitals, are computed by solving Poisson’s equation,

[TABLE]

using the method described in Ref. Sato et al. (2016). It should be noted that $\hat{W}^{\phi_{j}}_{\phi_{i}}({\bf\it r})$ is time independent, and Eq. (32a) needs to be solved only once before the simulation. On the other hand, $\hat{W}^{\phi_{j}}_{\chi_{i}}({\bf\it r},t)$ is time dependent, and should be computed at every time step. However, since its source $\phi^{*}_{j}({\bf\it r})\chi_{i}({\bf\it r},t)$ [See Eq. (32b).] and operand $\{\phi_{j}({\bf\it 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 Sato et al. (2016).

III.2 Time-propagation with exponential time differencing fourth-order Runge-Kutta scheme

For an efficient propagation of the EOM of channel-orbital based TDCIS, we use the exponential time differencing fourth-order Runge-Kutta scheme (ETDRK4) by Krogstad Krogstad (2005); Hochbruck and Ostermann (2010); Kidd et al. (2017). To this end, we arrange $C_{0}$ and $\quantity{\chi_{i}}$ into a unified vector ${\bf\it\chi}=(C_{0},\chi)^{T}$ and rewrite the EOMs of $C_{0}$ and $\quantity{\chi_{i}}$ by a matrix form

[TABLE]

where ${\bf\it h}$ is a chosen stiff part of the right-hand side of the EOM (See below), and ${\bf\it W}[{\bf\it\chi},t]$ is a nonstiff remainder. We choose the stiff part ${\bf\it h}$ to be either (i) the field-free one-electron Hamiltonian $\hat{h}_{\text{0}}$ or (ii) the totality of the one-electron Hamiltonian $\hat{h}_{\text{0}}+\hat{h}_{\text{ext}}(t)$. For the first case (i) with time-independent ${\bf\it h}$, the time propagation from ${\bf\it\chi}_{n}={\bf\it\chi}(t_{n})$ to ${\bf\it\chi}_{n+1}={\bf\it\chi}(t_{n}+\Delta t)$ is given by

[TABLE]

where $f_{1}$, $f_{2}$, and $f_{3}$ are defined as

[TABLE]

where $z=-i{\bf\it h}\Delta t$, $\varphi_{0}(z)=e^{z}$, and

[TABLE]

${\bf\it W}_{n}$, ${\bf\it W}_{a}$, ${\bf\it W}_{b}$, and ${\bf\it W}_{c}$ are given by

[TABLE]

where ${\bf\it a}_{n}$, ${\bf\it b}_{n}$, and ${\bf\it c}_{n}$ are intermediate vectors given as

[TABLE]

The operator exponential $\varphi_{0}(z)$ and $\varphi_{0}(z/2)$ in the spherical-FEDVR basis are approximated by the Padé (3/3) approximation, and higher-order $\varphi_{k}$ functions are obtained by successively applying Eq. (36). The denominator of the Padé approximation is factorized and operated by the matrix iteration method Sato et al. (2016). We follow Ref. Hochbruck and Ostermann (2010) for the modification required for a time-dependent stiff part ${\bf\it h}$. In the absence of linear part for $C_{0}$, time-propagation of $C_{0}$ reduces to well-known fourth-order Runge-Kutta scheme.

III.3 Expectation value

The expectation value of one-body operator $\expectationvalue{O}=\bra{\Psi}O\ket{\Psi}$ can be evaluated in the LG case as Rohringer et al. (2006); Sato et al. (2018b)

[TABLE]

The VG expression is obtained by simply replacing $\quantity{C_{0},\chi_{j}}$ with $\quantity{D_{0},\eta_{j}}$, and the rVG one by replacing $\quantity{\phi_{j},\chi_{j}}$ with $\quantity{\phi^{\prime}_{j},\chi^{\prime}_{j}}$. The Ehrenfest theorem $\frac{d}{dt}\expectationvalue{O}=-i\bra{\Psi}[O,\hat{H}]\ket{\Psi}$ does not hold for TDCIS. Instead, we evaluate the time derivative of $\expectationvalue{O}$ as Sato et al. (2018b),

[TABLE]

in the LG case. $\quantity{C_{0},\chi_{j}}$ is to be replaced with $\quantity{D_{0},\eta_{j}}$ for VG. The rVG case needs extra terms Sato et al. (2018b):

[TABLE]

Equations (39), (40), and (III.3) are valid not only for atoms but also any multielectron system including molecules.

III.4 Ionization probability

To conveniently analyze how ionization proceeds using the TD-MCSCF wavefunctions with time-varying orbitals, we have previously introduced Sato and Ishikawa (2013) a domain-based $n$-fold ionization probability $P_{n}$, defined as a probability to find $n$ electrons in the outer region $\absolutevalue{{\bf\it r}}>R_{\text{ion}}$ and the other $N-n$ electrons in the inner region $\absolutevalue{{\bf\it r}}<R_{\text{ion}}$ with a given distance $R_{\text{ion}}$ from the origin. This quantity is gauge-invariant even during the pulse, unlike the population of the (field-free) continuum levels Ishikawa and Sato (2015); Sato et al. (2018c). To apply this approach to the TDCIS method with channel orbitals, it is reasonable to assume that the occupied orbitals $\{\phi_{i}\}$ are localized inside the inner region, i.e, $\phi_{i}({\bf\it r})=0$ for $\absolutevalue{{\bf\it r}}>R_{\text{ion}}$. Then, the yield of the neutral species, or the survival probability $P_{0}$ is computed as,

[TABLE]

where $\innerproduct{\chi_{i}}{\chi_{j}}_{<}$ is the overlap integral in the inner region

[TABLE]

Noting that the atom described by a TDCIS wavefunction is at most singly ionized, we obtain the single-ionization 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:

[TABLE]

where $I_{0}$ is the peak intensity, $\omega$ is the central frequency, $T=2\pi/\omega$ is the period, and $N_{\text{opt}}$ is the total number of optical cycles.

IV.1 Helium

First, we consider helium atom exposed to a laser pulse with $I_{0}=4.0\times 10^{14}\ \text{W/cm}^{2}$, $\lambda=400\ \text{nm}$, and $N_{\text{opt}}=12$. In this condition, exact numerical solution of the TDSE is available Feist et al. (2009); Pazourek et al. (2011), 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=\hat{z}$ 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_{\text{ion}}=20$ a.u. in Fig. 3 (See Appendix for the $R_{\text{ion}}$ dependence of $P_{0}$). The conventional VG treatment strongly overestimates ionization. Recalling that tunneling ionization is the first process of the three-step model Corkum (1993); Piraux et al. (1993), 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.

IV.2 Neon

We next consider a neon atom subject to a laser field with $\lambda=$ 800 nm, $I_{0}=8.0\times 10^{14}$ W/cm2, and $N_{\text{opt}}=3$ and discuss convergence with respect to the maximum angular momentum $L_{\text{max}}$. We show the HHG spectra calculated with various values of $L_{\text{max}}$ in the LG and rVG in Fig. 4.

Figure 4. (a) shows the equivalence between the LG and rVG for sufficiently large $L_{\text{max}}(=100)$. As can be seen in Fig. 4. (b), which shows LG results, $L_{\text{max}}=60$ is not sufficient to obtain a converged result. On the other hand, the rVG requires far less $L_{\text{max}}$; even $L_{\text{max}}=40$ well reproduces the result with $L_{\text{max}}=100$, and the spectrum is converged with $L_{\text{max}}=60$ [Fig. 4. (c)]. This observation indicates that the rVG TDCIS is simultaneously as 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 Sato et al. (2018b) 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 Rohringer et al. (2006), 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 exact 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.

acknowledgments

We are grateful to Prof. J. Burgdörfer of Vienna University of Technology for providing us with He TDSE benchmark results. This research was supported in part by a Grant-in-Aid for Scientific Research (Grants No. 16H03881, No. 17K05070, No. 18H03891, and No. 19H00869) from the Ministry of Education, Culture, Sports, Science and Technology (MEXT) of Japan. This research was also partially supported by JST COI (Grant No. JPMJCE1313), JST CREST (Grant No. JPMJCR15N1), and by Quantum Leap Flagship Program of MEXT. T. T. gratefully acknowledges support from the Graduate School of Engineering, The University of Tokyo, Doctoral Student Special Incentives Program (SEUT Fellowship).

Appendix A $R_{\text{ion}}$ dependence on $P_{0}$

We show $P_{0}$ from the LG TDCIS results with various values of $R_{\text{ion}}$ in fig. 5. We can see stepwise ejection of electron wavepackets that propagate outwards; the larger $R_{\text{ion}}$ is, the later $P_{0}$ is depleted. In the cases of $R_{\text{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).