Gauge-invariant formulation of time-dependent configuration interaction singles method

We propose a gauge-invariant formulation of the channel orbital-based time-dependent configuration interaction singles (TDCIS) method [Phys. Rev. A 74, 043420 (2006)], one of the powerful ab initio methods to investigate electron dynamics in atoms and molecules subject to an external laser field. In the present formulation, we derive the equations of motion (EOMs) in the velocity gauge using gauge-transformed orbitals, not fixed orbitals, that are equivalent to the conventional EOMs in the length gauge using fixed orbitals. The new velocity-gauge EOMs avoid the use of the length-gauge dipole operator, which diverges at large distance, and allows to exploit computational advantages of the velocity-gauge treatment over the length-gauge one, e.g, a faster convergence in simulations with intense and long-wavelength lasers, and the feasibility of exterior complex scaling as an absorbing boundary. The reformulated TDCIS method is applied to an exactly solvable model of one-dimensional helium atom in an intense laser field to numerically demonstrate the gauge invariance. We also discuss the consistent method for evaluating the time derivative of an observable, relevant e.g, in simulating high-harmonic generation.


I. INTRODUCTION
Time-dependent configuration interaction singles (TDCIS) method is one of the powerful ab initio methods to investigate laser-driven electron dynamics in atoms and molecule . In the TDCIS method, the time-dependent electronic wavefunction is given by the configuration interaction (CI) expansion, where Φ is the ground-state Hartree-Fock (HF) wavefunction, and Φ ia is a singly-excited configuration-state function (CSF), replacing an occupied HF orbital φ i in Φ with a virtual (unoccupied in Φ) orbital φ a , and the electron dynamics is described through the time evolution of the CI coefficients, C 0 and {C ia }. Compared to more involved ab initio wavefunctionbased approaches [25] such as time-dependent multiconfiguration self-consistent-field (TD-MCSCF) methods [26][27][28][29][30][31][32][33], time-dependent R-matrix based approaches [34][35][36], or timedependent reduced density-matrix approach [37,38], distinct advantages of the TDCIS method include a low computational cost and the conceptual simplicity to analyze simulation results. Furthermore, an equivalent, effective one-electron theory with coupled channels has been developed [2], which introduces the orbital-like quantity, called channel orbital, and rewrites EOMs for CI coefficients with those for channel orbitals {χ i (r, t)} with no reference to virtual orbitals. This reformulation removes the bottleneck of the CI coefficientbased TDCIS method to compute all (or, at least sufficiently many, including bound and continuum) virtual orbitals prior to the simulation, and thus particularly useful in grid-based simulations.
Despite this advantage, numerical applications of the channel orbital-based TDCIS method has been limited to Refs. [2,14,15] for a one-dimensional Hamiltonian and Ref. [1] for noble gas atoms with a Hartree-Slater potential, as far as we know, and the vast majority of applications to date have adopted the CI coefficient-based approach  , except for the use of {χ i } as intermediate quantities in evaluating photoelectron spectra [18]. The preference of CI coefficientbased approach might be partially due to the high symmetry of atomic systems, for which the stationaly Hartree-Fock operator decouples for different angular momenta [4], making it a relatively feasible task to obtain all virtual orbitals (within a given radial grids or radial basis functions) for the lowest few angular momenta. The channel orbital-based approach would be more suited, on the other hand, to simulations of electron dynamics with intense and/or long-wavelength laser fields, requiring much longer angular momentum expansion [39][40][41], and moreover to grid-based molecular applications, where obtaining a sufficient spectrum of virtual levels could be unacceptably expensive.
However, the TDCIS method, either in the CI coefficientbased or channel orbital-based formulation, suffers from the lack of gauge invariance, as a general consequence of relying on truncated CI expansion with fixed orbital functions. Previously, the length gauge (LG) has been employed e.g, in Ref. [2][3][4][5][6][7][8][9][10][11][12][13][14][15][16], and the velocity gauge (VG) in Ref. [17][18][19][20][21][22][23][24]. Although gauge dependence of the TDCIS method using fixed orbitals has been noted already in Ref. [2], comparative assessment of the LG and VG treatments (within the grid-based TDCIS) has not been reported to the best of our knowledge, except for being briefly mentioned in Ref. [42]. In particular, the channel orbital-based approach [2] has been applied only in the LG [1,2], and as shown below in this paper, the VG treatment with fixed orbitals is not very appropriate for applications to high-field phenomena. This is a serious drawback, since for an efficient simulation of molecules, it is highly ap-preciated to take advantage of the velocity-gauge treatment, e.g, the feasibility of exterior complex scaling [43,44] as an absorbing boundary, to reduce the computational cost related to the number of grid points.
In the present work, we propose a gauge-invariant reformulation of the channel orbital-based TDCIS method. To this end, instead of applying the fixed-orbital TDCIS ansatz to the velocity-gauge time-dependent Schrödinger equation (TDSE), we adopt the formulation using unitary-rotated orbital φ p (t) = U (t)φ p , where U (t) is the gauge transformation operator connecting the (exact) solution of TDSE in the LG and VG. The resulting EOMs in the reformulated VG is equivalent to the LG ones with fixed orbitals by construction, and at the same time allows to exploit advantages of the velocitygauge simulations as mentioned above.
This paper proceeds as follows. In Sec. II, after defining the target Hamiltonian and the gauge transformation in Sec. II A and reviewing the TDCIS method using fixed orbitals both in the CI coefficient-based [Sec. II B] and channel orbital-based [Sec. II C] approaches, we present the gaugeinvariant reformulation in Sec. II D, and a consistent method for evaluating the time derivative of one-electron observables in Sec. II E. Then in Sec. III we apply the channel orbitalbased TDCIS method, using LG with fixed orbitals, VG with fixed orbitals, and the reformulated VG, to the model onedimensional (1D) Hamiltonian to compare the results of various TDCIS approaches with numerically exact TDSE results, and demonstrate the importance of non-Ehrenfest method to compute dipole acceleration. Finally, concluding remarks are given in Sec. IV. The Hartree atomic units are used throughout unless otherwise noted.

A. System Hamiltonian and gauge transformation
Let us consider an atom or a molecule consisting of N electrons interacting with an external laser field. In this work, we restrict our treatment in the clamped-nuclei approximation and the electron-laser interaction within the electric dipole approximation. Then the exact description of the system dynamics is given by the solution Ψ L (t) of TDSE, with the system Hamiltonian H L (t) = H 0 + H ext L (t), where H 0 is the field-free electronic Hamiltonian where r k and p k = −i∇ k are the coordinate and canonical momentum of an electron, h(r, p) = 1 2 p 2 + v n (r), with v n being the electron-nucleus interaction. Here we are considering the LG treatment, where the electron-laser interaction H ext L is given by where E(t) is the laser electric field. As well known, the system dynamics is equivalently described in the VG, of which the wavefunction Ψ V is connected with the LG one through with a unitary transformation where is the vector potential, and we arbitrarily include the second term in the exponential, which is a c-number, to avoid appearance of terms proportional to |A| 2 in subsequent equations. Then we substitute with One should carefully note that the present proof of equivalence of the LG and VG treatments, Eqs. (3) and (8), with the transformation of Eq. (7), applies only to the exact solution of TDSE. See e.g, Ref. [45][46][47] for deeper discussions on the gauge transformation within TDSE, and Ref. [25] for the gauge invariance of TD-MCSCF methods. For a compact presentation of the many-electron theory, we rewrite the system Hamiltonian in the second quantization, where {ĉ † pσ } and {ĉ pσ } are the creation and annihilation operators, respectively, for the set of spin-orbitals given as a direct product {φ p } ⊗ {s ↑ , s ↓ } of orthonormal spatial orbitals {φ p } and up-spin (down-spin) functions s ↑ (s ↓ ). The operatorsĥ,r, andp are defined, respectively, asĥ = where h pq , r pq , and p pq are the matrix elements of h, r, p, respectively, in terms of {φ p }, and The TDSE of the LG, Eq. (3), and VG, Eq. (8), read with the transformation µσĉµσ is the number operator. In this work, we consider a closed-shell system with even number of electrons, and choose as {φ p } the timeindependent Hartree-Fock (HF) orbitals satisfying the canonical, restricted HF equation where p is the orbital energy, andŴ φ φ is the electrostatic potential of a product φ * (r)φ (r) of given orbitals, defined in the real space as As usual, we separate the full set of HF orbitals {φ p } into the occupied orbitals {φ i } which are occupied in the HF groundstate wavefunction (also referred to as the reference , and the virtual orbitals {φ a } which are unoccupied in |Φ .

B. Review of CI coefficient-based TDCIS with fixed orbitals
We write the second-quantized version of Eq. (1), for the LG case, as where The equations of motion for the CI coefficients have been derived [2] by inserting Eq. (19) into the LG TDSE, Eq. (14a), and closing from the left with the reference and singly-excited CSFs, Conceptually more proper derivation of Eqs. (20) is based on Dirac-Frenkel variational principle, which considers the Lagrangian and requires ∂L L /∂C * 0 = ∂L L /∂C * ia = 0. SubstitutingĤ L of Eq. (10a) into Eqs. (20), using the Slater-Condon rule for the Hamiltonian matrix elements, and noting the canonical condition f pq = p δ pq , the EOMs for the length gauge are derived as [2] where the action of the operatorF i on a given orbital φ is defined aŝ References [17][18][19][20][21][22][23][24] have used the same expansion in terms of fixed CSFs also in the VG case, and required Eqs. (20) to hold, withĤ L , C 0 , and C ia replaced withĤ V , D 0 , and D ia . This is equivalent to consider the following Lagrangian, C. Review of Channel orbital-based TDCIS with fixed orbitals An interesting reformulation of the above-described TD-CIS method, as mentioned in Sec. I, has been proposed in Ref. 2, which introduces the time-dependent channel orbitals |χ i that collects all the single excitations originating from an occupied orbital |φ i , and rewrites the EOMs in terms of C 0 and {|χ i } as whereP =1 − j |φ j φ j |. According to these EOMs and the initial conditions [C 0 (t → −∞) = 1, and {C ia (t → −∞) = 0} ⇐⇒ {χ i (t → −∞) ≡ 0}], the channel orbitals |χ i gets gradually populated along with the laser-electron interaction, measuring an excitation of an electron out of |φ i . See Ref. [2] for interesting properties of the channel orbitals. It is also possible to formulate the channel orbital-based scheme based on the velocity gauge TDCIS using fixed orbitals, although not previously considered. We, therefore, introduce the analogous quantity and rewrite Eqs. (26) as Hereafter, we refer to the method based on Eqs. (28), i.e, the channel orbital-based TDCIS in the length gauge with fixed orbitals, simply as LG method, and that based on Eqs. (30), i.e, the channel orbital-based TDCIS in the velocity gauge with fixed orbitals, as VG method, for notational brevity.

D. Channel orbital-based TDCIS in the velocity gauge with rotated orbitals
The gauge dependence of the LG and VG treatments, Eqs. (28) and (30), results from the fact that the ansatz of Eqs. (19) and (24), both using fixed orbitals, cannot be connected with the transformation, Eq. (16), as is generally the case for truncated CI expansion using fixed orbitals. For a method to be gauge invariant, the underlying Lagrangian in LG and VG cases should be numerically the same when evaluated with the solution of respective EOMs, which does not hold in the present case, L L (t) = L V (t), with Eqs. (21) and (25).
Thus we define the total wavefunction |Ψ V (t) , transformed from |Ψ L (t) to the velocity gauge, as with |Ψ L (t) constructed with the solution of CI coefficientbased EOMs in the LG, Eqs. (22). Here |Φ =Û (t)|Φ and |Φ ia =Û (t)|Φ ia = σĉ † aσĉ iσ |Φ / √ 2 are the reference and singly-excited CSF constructed with unitary rotated orbitals, i.e, |φ p =Û |φ p andĉ pσ =Û (t)ĉ pσÛ −1 (t). It should be noted that |Ψ V cannot be rewritten into the form of Eq. (24) in general. Associated with this wavefunction, we consider the following Lagrangian, The equivalence of this approach to the LG treatment is readily confirmed by seeing One may naively expect that L V of Eq. (32), which differs from L V of Eq. (25) only by the replacement of Ψ V with Ψ V , leads to the EOMs of Eqs. (26) This is not the case, however, due to the time dependence of the rotated CSFs, e.g, Φ |Φ ia = iE(t) · Φ |r|Φ ia , and after extracting these time dependence, Eq. (32) reads where ∂ c t time differentiates CI coefficients only. Now requiring ∂L V /∂C * 0 = ∂L V /∂C * ia = 0, or equivalently, substituting the back transformation |φ p =Û −1 |φ p into Eqs. (22) derives whereF i is given by Eq. (23) with {φ j } replaced with {φ j }.
Equations (35) are the CI coefficient-based TDCIS EOMs based on the Lagrangian of Eq. (32). Although this approach is guaranteed to be equivalent to the CI coefficient-based LG TDCIS, it brings no numerical gain over Eqs. (22), peculiarly including both E · r and A · p, and requiring extensive gauge transformation of all occupied and virtual orbitals. None the less, a useful method can be derived, if one switches to the channel orbital-based scheme by defining the rotated channel functions, Then we use dÛ /dt = i(E ·r +N |A| 2 /2)Û , and notê . Although several terms in the EOMs still involve the dipole operator, they all apply to the (rotated) occupied orbital which is localized around nuclei, thus posing no difficulty in enjoying the same advantages of VG propagations of orbitals [39][40][41].

E. Evaluation of the time derivative of an observable
Let us next consider how to compute expectation value of a one-electron operator Ô (t) = Ψ(t)|Ô|Ψ(t) , and its time derivative d Ô /dt. For exact solution of TDSE, |Ψ = −iĤ|Ψ , the time derivative is given by known as the Ehrenfest expression. For an approximate method, however, the Ehrenfest theorem, Eq. (38b), generally does not hold, and one should explicitly evaluate the time derivative as Eq. (38a). Important exceptions include those theories using time-dependent orbitals evolving to satisfy the time-dependent variational principle, such as time-dependent Hartree-Fock (TDHF), TD-MCSCF, and time-dependent density functional theory. See Ref. [41] for more details.
The TDCIS expectation value of a one-electron operatorÔ is given [2] by in the LG case. That for the VG is given by replacing C 0 with D 0 in the above equation, and for the rVG by replacing {φ j , χ j } with {φ j , χ j }. The expression for the time derivative, in the LG case, is derived by using Eqs. (28) in Eq. (38a) as The VG expression is also given by the above equation with C 0 replaced with D 0 , and that for the rVG is Although Eqs (40) and (41) look rather complicated, their evaluations are straightforward given the time derivatives of working variables C 0 , {χ i }, etc, which are necessary, in any case, to propagate the EOMs.

III. NUMERICAL EXAMPLES
In this section, we numerically apply the channel orbitalbased TDCIS method in the LG, VG, and rVG to the 1D model Helium atom, using the computational code developed by modifying an existing TDHF code used in our previous work [30,33,48]. The field-free electronic Hamiltonian is given by for two electronic coordinates z 1 and z 2 , and the laser-electron interaction E(t) · r and A(t) · p are replaced with E(t)z and A(t)p z = −iA(t)∂/∂ z , respectively, in Eqs. (28), (30) and (37). Orbitals are discretized on equidistant grid points with spacing ∆z = 0.4 within a simulation box −1000 ≤ z ≤ 1000, with an absorbing boundary implemented by a mask function of cos 1/4 shape at 10% side edges of the box. Each EOM is solved by the fourth-order Runge-Kutta method with a fixed time step size (1/10000 of an optical cycle). Spatial derivatives are evaluated by the eighth order finite difference method, and spatial integrations are performed by the trape- zoidal rule. We consider a laser electric field given by for 0 ≤ t ≤ τ , and E(t) = 0 otherwise, with a wavelength λ = 2π/ω 0 = 750 nm, a foot-to-foot pulse length τ of three optical cycles, and a peak intensity I 0 = E 2 0 for I 0 = 5 × 10 14 W/cm 2 and I 0 = 10 15 W/cm 2 . The 1D Hamiltonian, computational details, and the applied laser field are the same as used in Ref. [48] to facilitate comparison with TDSE results in Ref. [48].
First, we compare the time-dependent dipole moment z (t) obtained with TDCIS approaches with that of TDSE in Fig. 1, which immediately reveals a strong gauge dependence of fixed-orbital approaches, i.e, the large difference between LG and VG results. One should note that the comparison of LG and VG results alone can tell nothing about the preference of either approach; TDCIS method in both LG and VG are the first approximation in the hierarchy of CI expansions, which, at the full-CI limit, would be gauge invariant. The point here is that the LG scheme outperforms the VG scheme in comparison to the exact TDSE result as clearly seen in Fig. 1, which convinces one an empirical preference of the LG treatment. On the other hand, the results of LG and rVG agree perfectly within the graphical resolution, numerically demonstrating the theoretical gauge invariance.
Next, we consider the dipole acceleration a (t) defined as the time derivative of the kinematic momentum, whereπ =p z for the LG, andπ =p z +A(t) for the VG. In the exact TDSE case, applying Eqs. (38) forÔ =π (also taking into account the trivial, explicit time dependence of π(t) in the VG case) derives where ∂v nuc /∂z = −∂/∂ z 2(z 2 + 1) −1/2 = 2z(z 2 + 1) −3/2 for the 1D Hamiltonian. Numerically achieving the theoretical equivalence of Eq. (44) and (45) (45), with that of TDSE in Fig. 2, clearly showing a better agreement of the results of the former approach with that of TDSE. From this result, and also by the fact that being based on Eq. (44) guarantees that the HHG spectra obtained from the velocity π (t) and the acceleration a (t), at the convergence, properly relate to each other [45], we consider that Eq. (44), together with Eq. (40) or Eq. (41), should be adopted as a consistent method for evaluating the dipole acceleration.
Then we compare the time evolution of the dipole acceleration [ Fig. 3] and the HHG spectrum [ Fig. 4] obtained as the modulus squared of the Fourier transform of the dipole acceleration obtained with TDCIS method in LG, VG, and rVG [based on Eq. (44)] with those of TDSE. We observe that (1) the LG and rVG results are identical to within the scale of the figure, (2) they also show a good agreement with TDSE  . HHG spectrum of 1D-He exposed to a laser pulse with a wavelength of 750 nm and an intensity of (a) 5×10 14 W/cm 2 and (b) 1×10 15 W/cm 2 . Comparison of the results with TDCIS in the LG, VG, and rVG with that of TDSE. results, (3) and in contract, the VG results strongly deviate from all the other results. Especially, Fig. 4 shows a remarkable agreement of the TDCIS spectra in the LG and rVG and the TDSE one, suggesting that the TDCIS method would be a useful computational method for studying HHG process in more complex atoms and molecules, in particular, when the present rVG treatment is combined with advanced, velocity gauge-specific computational techniques.

IV. CONCLUSIONS
In this work, we propose a gauge-invariant formulation of the channel orbital-based TDCIS method for ab initio investigations of electron dynamics in atoms and molecules. Instead of using fixed orbitals both in length-gauge and velocitygauge simulations, we adopt, in the velocity-gauge case, the EOMs derived with unitary rotated orbitals |φ p (t) = U (t)|φ p using gauge-transforming operatorÛ (t), which replaces the length-gauge operator E ·r appearing in the lengthgauge EOMs with the velocity-gauge counterpart A · p, while keeping the equivalence to the length-gauge treatment. This would make it possible to take advantages of the velocitygauge simulation over the length-gauge one, e.g, the faster convergence of simulations of atoms interacting with an intense and/or long-wavelength laser field, with respect to the maximum angular momentum included to expand orbitals, and the native feasibility of advanced absorbing boundaries such as the exterior complex scaling. Applications to real atoms and molecules with the three-dimensional Hamiltonian will be presented elsewhere.