Application of the time-dependent surface flux method to the time-dependent multiconfiguration self-consistent-field method

We present a numerical implementation of the time-dependent surface flux (tSURFF) method [New J. Phys. 14, 013021 (2012)], an efficient computational scheme to extract photoelectron energy spectra, to the time-dependent multiconfiguration self-consistent-field (TD-MCSCF) method. Extending the original tSURFF method developed for single particle systems, we formulate the equations of motion for the spectral amplitude of orbital functions constutiting the TD-MCSCF wave function, from which the angle-resolved photoelectron energy spectrum, and more generally, photoelectron reduced density matrices (RDMs) are readiliy obtained. The tSURFF method applied to the TD-MCSCF wave function, in combination with an efficient absorbing boundary offered by the infinite-range exterior complex scaling, enables accurate {\it ab initio} computations of photoelectron energy spectra from multielectron systems subject to an intense and ultrashort laser pulse with a computational cost significantly reduced compared to that required in projecting the total wave function onto scattering states. We apply the present implementation to the photoionization of Ne exposed to an attosecond extreme-ultraviolet (XUV) pulse and above-threshold ionization of Ar irradiated by an intense mid-infrared laser field, demonstrating both accuracy and efficiency of the present method.


I. INTRODUCTION
The rapid progress in experimental techniques for highintensity, ultrashort optical pulses, has lead to the advent and development of strong-field physics and attosecond science [1,2], with the ultimate goal to directly measure and control electron motion in atoms, molecules, and solids. Although the time-dependent Schrödinger equation (TDSE) provides a rigorous theoretical framework to investigate electron dynamics, solving it for multielectron systems poses a major challenge. To simulate multielectron dynamics in intense laser fields, the timedependent multiconfiguration self-consistent field (TD-MCSCF) methods have been developed [3][4][5][6][7][8][9][10][11][12] which express the total wave function as a superposition, of the Slater determinants (or electronic configurations) Φ I (t) constructed from single-electron orbitals. Both the configuration-interaction (CI) coefficients {C I } and orbitals are varied in time. If we consider all the possible configurations for a given number of orbitals as in the multiconfiguration time-dependent Hartree-Fock (MCTDHF) method [4][5][6], the computational cost increases factorially with the number of electrons. To overcome this difficulty, we have recently formulated and * ykormhk@atto.t.u-tokyo.ac.jp numerically implemented the time-dependent completeactive-space self-consistent-field (TD-CASSCF) method [9], and even more general and less demanding timedependent occupation-restricted multiple-active-space (TD-ORMAS) method [12]. The former introduces orbital classification into doubly occupied core orbitals and fully correlated active orbitals; the core orbitals are further divided into time-independent frozen core and timedependent dynamical core. The latter further divides active orbitals into an arbitrary number of subgroups, specifying the minimum and maximum numbers of electrons distributed in each subgroup. Flexible description of the wave function offered by the orbital classification and occupation restriction enables converged simulations of highly nonlinear, correlated multielectron dynamics in systems containing several tens of electrons.
Photoelectron energy spectra (PES) and their angular distribution, or angle-resolved photoelectron energy spectra (ARPES) are among important experimental probes. In principle, they could be calculated by projecting the departing photoelectron wave packet onto plane waves or Coulomb waves. This approach, however, requires retaining the complete wave function without being absorbed, leading to a huge simulation box and prohibitive computational cost. To circumvent this difficulty, Tao and Scrinzi have devised the time-dependent surface flux (tSURFF) method [13], which extracts PES by integrating the electron flux through a surface. Thus it allows one to use an absorbing boundary, bringing significant cost reduction. The tSURFF method was first developed for single electron systems [13] and then applied to multielec-tron simulations with, e.g., the hybrid antisymmetrized coupled channels method [14], the time-dependent configuration interaction singles method [15], and the timedependent density functional theory [16].
In this study, we combine the tSURFF method with the TD-CASSCF and TD-ORMAS methods in order to extract photoelectron energy spectra from ab initio simulations of multielectron dynamics in atoms subject to intense and/or ultrashort laser pulses. Under a physically reasonable assumption that the nuclear potential and interelectronic Coulomb interaction are negligible for photoelectron dynamics in the region distant from the nuclei, we derive the equations of motion for the momentum amplitudes of each orbital. They contain an additional term compared with the single-electron case. While we use the infinite-range exterior complex scaling (irECS) [17,18] as an efficient absorbing boundary, the present tSURFF implementation also supports other absorption boundaries such as mask functions and complex absorbing potentials. We achieve highly accurate calculations of angle-resolved PES with considerably reduced computational costs.
This paper is organized as follows. We briefly review the TD-CASSCF and TD-ORMAS methods in Sec. II and the tSURFF method for single-electron systems in Sec. III. In Sec. IV, we describe our theoretical application and numerical implementation of tSURFF to the TD-MCSCF method. Numerical examples are presented in Sec. V. Conclusions are given in Sec. VI. We use Hartree atomic units unless otherwise indicated.

II. THE TD-CASSCF AND TD-ORMAS METHODS
We consider an N -electron atom with N ↑ up-spin and N ↓ down-spin electrons (N = N ↑ + N ↓ ) in a laser field E(t). Its dynamics is described within the dipole approximation by the Hamiltonian, with, where Z denotes the atomic number and A(t) = − t −∞ E(t )dt the vector potential. We use the velocity gauge, since ECS works only with it, and not with the length gauge [19].
In the TD-CASSCF and TD-ORMAS methods, the total wave function is given by, whereÂ denotes the antisymmetrization operator, Φ fc and Φ dc the closed-shell determinants formed with n fc frozen-core and n dc dynamical-core orbitals, respectively, and {Φ I } the determinants constructed from n a active orbitals. Whereas the active electrons are fully correlated among the active orbitals in the TD-CASSCF method, the TD-ORMAS method [12] further subdivides the active orbitals into an arbitrary number of subgroups, specifying the minimum and maximum number of electrons accommodated in each subgroup. The equations of motion (EOMs) describing the temporal evolution of the CI coefficients {C I } and the orbitals {ψ p } are derived on the basis of the time-dependent variational principle (TDVP) [12,20] and read, whereQ = 1 − q |ψ q ψ q | the projector onto the orthogonal complement of the occupied orbital space.F is a non-local operator describing the contribution from the interelectronic Coulomb interaction, defined aŝ where D and P are the one-and two-electron reduced density matrices, andŴ r s is given, in the coordinate space, by The matrix element R q p is given by, with h q p = ψ q |ĥ|ψ p . R q p 's within one orbital subspace (frozen core, dynamical core and each subdivided active space) can be arbitrary Hermitian matrix elements, and in this paper, they are set to zero. On the other hand, the elements between different orbital subspaces are determined by the TDVP. Their concrete expressions are given in ref. [12], where iX q p = R p q + h q p is used for working variables.

III. THE tSURFF METHOD FOR SINGLE-ELECTRON SYSTEMS
In this section, we briefly review the tSURFF method [13] for single-electron systems governed by TDSE, where V (r) denotes the nuclear potential. The photoelectron momentum amplitude (PMA) a(k, t) for momentum k and PES can in principle be approximately calculated by projecting the outgoing wave packet residing outside a given radius R s onto the plane wave (2π) −3/2 exp(ik · r) at a time t sufficiently after the pulse. Then, PMA is given by, where χ k (r, t) denotes the Volkov wave function and θ(x) the Heaviside function.
On the other hand, the tSURFF method calculates PES by a time integration of the wave function surface flux, based on an assumption that the nuclear potential does not affect the time evolution of the distant photoelectron wave packet. Under this assumption, the Volkov wave functions and photoelectron wave packets in the region |r| > R s are evolved by the same nuclear-potentialfree Hamiltonian By taking time derivative of Eq. (14), we obtain the EOM of the momentum amplitude, As shown in Ref. [13], since all the terms appearing in the commutator [h s , θ(R s )] contain delta functions δ(r − R s ) [see Eq. (28) below], wave functions and their spatial derivative only on the surface |r| = R s are required to solve Eq. (18). Hence, it is no longer needed to keep the whole wave function and allowed to use an absorbing boundary, which leads to a significant computational cost reduction.

A. Photoelectron reduced density matrix
To obtain PES in multielectron systems described by the multiconfiguration expansion Eq. (1), we define the photoelectron reduced density matrix (PRDM). Since our definition of ionization is based on the spatial domain |r| > R s , the single particle reduced density matrix of a photoelectron in the coordinate space can be defined as, In the momentum space, its elements are given by, The diagonal partP (k, k) is interpreted as photoelectron momentum distribution, and, then, PES is given by, A similar definition of the PRDM as Eq. (20) and the direct projection onto scattering states were used to compute photoelectron spectrum in Ref. [21].

B. EOMs of momentum amplitudes of orbitals
In this subsection, we derive the EOM of the momentum amplitude of orbital p, which appears in Eq. (20). We assume that the nuclear potentials are negligible for photoelectrons as in the single-electron case and additionally that the interelectronic Coulomb interaction does not affect the dynamics of photoelectrons in the region beyond the radius R s . Then, the orbital EOM Eq. (7) can be approximated as, so that EOMs of a p (k) is obtained as It should be noted that the second term in Eq. (24) represents a significant difference from the single-electron case [Eq. (18)]. This term includes the effect of the interelectronic Coulomb interaction inside R s and is not negligible; for example, the phase variation due to the energy of the ionic core is reflected in photoelectron momentum spectra through this term. It may not be a priori obvious if the Coulomb interaction between electrons in the outer region is negligible. A numerical validation of this approximation will be given in Sec. V.

C. Implementation
In this subsection, we describe the implementation of the tSURFF method to TD-MCSCF simulations. The implementation is based on our TD-CASSCF code [22] for atoms subject to a laser pulse linearly polarized in the z direction. The orbitals are discretized with spherical finite-element discrete-variable-representation (FEDVR) basis functions [23], where f k (r) and Y lmp (θ, φ) denote the FEDVR function and the spherical harmonics, respectively. The magnetic quantum number m p for orbital ψ p does not change due to the z-polarization of the laser pulse [22]. As an ab-sorbing boundary, infinite-range exterior complex scaling is implemented [17,18], which shows high efficiency and accuracy. On the other hand, we discretize the orbital momentum amplitudes with grid points in the spherical coordinates, where the Volkov wave function for a momentum k = (k, θ k , ψ k ) is given by, with j l (kr) being the spherical Bessel function of the first kind and Λ(k, t) the Volkov phase given by, The commutator in Eq. (18) can be rewritten as, with z component A z (t) of the vector potential A(t).
Using Eqs. (26) and (27) and introducing g p l (r) = k c p kl f k (r), we can decompose the first term of Eq. (24) into, where j * l (r) and g p l (r) denote the radial derivative of j * l (r) and g p l (r), respectively. For the second term in Eq. (24), since { ψ q (t)|F |ψ p (t) −R q p } is evaluated during the orbital propagation [Eq. (7)], we can reuse it.
We integrate Eq. (24) with the first order exponential integrator [24] by treating Eq. (24) as simultaneous inhomogeneous linear differential equations. The evolution of a vector a(k, t) ≡ {a p (k, t)} from the time t to t + ∆t is described as where the matrix X(t) = {X qp (t)} and vector S(t) = {S p (t)} are defined as, We compute exp[i∆tX(t)] by directly diagonalizing X(t) at every time step, which is not demanding since the size of X(t) is [N orb × N orb ] and the number of orbitals N orb usually falls within the range from a several to several tens.

V. NUMERICAL RESULTS
In this section, we present numerical applications of the implementation of tSURFF to the TD-MCSCF method described in the previous section. The electric field of the laser pulse is assumed to have the following shape for simulations of Ne atom (Sec. V A): and for Ar atom (Sec. V B): where I 0 is a peak intensity, T is a period at the central frequency ω = 2π/T , and N opt is the total number of optical cycles.

A. Ne atom
We first calculate the PES of a neon atom subject to an attosecond pulse with a peak intensity of 2.5×10 12 W/m 2 , a wavelength of 12.398 nm corresponding to 100 eV photon energy, and N opt = 16. The results by tSURFF and direct projection on plane waves are compared. As an absorbing boundary we use irECS with a scaling radius R 0 of 40 a.u for tSURFF and 400 a.u for direct projection. The latter is large enough to hold the departing wave packet from two-photon ionization. R s = 40 a.u. for tSURFF, and the wave packet outside this radius is used for projection.
We do TD-CASSCF simulations with 3 kinds of orbital classifications (n fc , n dc , n a ) = (0, 0, 5), (0, 0, 9), and (1,0,8), where n fc , n dc , and n a are the number of frozencore, dynamical-core, and active orbitals, respectively. Note that the first and second correspond to the timedependent Hartree-Fock (TDHF) [25] and MCTDHF methods, respectively. The results are shown in Fig. 1. We see single photon ionization peaks (around 30 -90 eV) and two photon above threshold ionization (ATI) peaks (around 130 eV -190eV) from 2s and 2p orbitals. The agreement between the results by tSURFF and direct projection is excellent. This shows the validity of the neglect of the electron-electron and nucleus-electron Coulomb interaction beyond R s assumed in tSURFF. While the single photon ionization peaks from 2s and 2p orbitals are expected to be at 51.5 and 78.4 eV, respectively, based on the experimental values of the binding energies [26], the peaks in the calculated spectra are located at 47.5 and 76.7 eV in Fig. 1(a), 48.9 and 77.9 eV in (b), and 48.7 and 77.9 eV in (c), coming closer to the experimental positions with increasing number of orbitals. Otherwise, TDHF gives results similar to the MCTDHF and TD-CASSCF ones for this process.

B. Ar atom
Next, we simulate ATI of a Ar atom subject to an intense visible laser pulse using the TD-ORMAS method and discuss the effect of electronic correlation. We consider a pulse, which has a peak intensity of 2.0 × 10 14 W/cm 2 , a wavelength of 532 nm, and a pulse width of N opt = 14 optical cycles. The ponderomotive energy U p is 5.285 eV. The temporal shape and energy spectrum of the pulse are shown in Fig. 2.
Here, we subdivide n a active orbitals into 4 and (n a −4) orbitals, as schematically illustrated in Fig. 3  with up to double (SD) or triple (SDT) excitation are considered, respectively. Figure 4 shows PES calculated with different orbital classifications. We can recognize the direct cutoff at 2U p and rescattering cutoff at 10U p . In Fig. 4(a), the results with different numbers of active orbitals within single and double (SD) excitation are compared. The spectrum is nearly converged with 25 and 29 active orbitals. Figure 4(b) compares the results with SD and SDT excitation restriction and full CI (TD-CASSCF), with 13 active orbitals. The SDT and full CI results almost overlap each other, indicating that SDT is sufficient for numerical convergence. Thus, the result using 25 active orbitals with SDT excitation is expected to be numerically nearly exact. Then, in Fig. 4(c), we compare the PES calculated using 13, 20 and 25 active orbitals with SDT excitation restriction and also the TDHF result. The peak positions slightly depend on the number of orbitals, as we have also seen in Fig. 1. Moreover, in the TDHF case, the peaks are significantly broadened. The ATI peak position E n corresponding to n-photon absorption is given by, where I p is the ionization potential. The difference in peak position observed in Fig. 4 can be attributed to that in I p , which depends on the number of orbitals and excitation restriction. In addition, in mean-field approaches such as TDHF, the ionization potential effectively increases as ionization proceeds and the electron density near the nucleus decreases [27]. This results in the peak broadening.
Finally, in Fig. 5, we compare ARPES calculated with the TD-ORMAS method using 25 active orbitals with SDT excitation restriction and with the TDHF method using 4 dynamical-core orbitals. We see difference in detailed structure. In particular, the high-order (> 2U p ) rescattering contribution has much broader angular distribution in the TD-ORMAS result. The difference can also be clearly observed in Fig 6, which shows the photoelectron angular distribution (PAD) at 10 eV (< 2U p ) and 40 eV (> 2U p ), representatives of the lower and higher energy regions, respectively. At 10 eV photoelectron energy, where the main contribution is from direct ionization, all the results exhibit similar behavior. In contrast, at 40 eV photoelectron energy, for which rescattering from the parent ion is involved and thus strong electron correlation is expected, the calculated PAD varies with the number of orbitals till it approximately converges with n a = 25. Especially, the TDHF method significantly underestimates the yield in the direction (90 • ) perpendicular to the laser polarization. This indicates that electronic correlation is non-negligible in detailed discussion of ATI ARPES. The results with different number of orbitals (n fc , n dc , na) and excitation restrictions (SD, SDT, full-CI) are compared.

VI. SUMMARY
We have presented a successful numerical implementation of tSURFF to the TD-MCSCF (TD-CASSCF and TD-ORMAS) methods to extract angle-resolved photoelectron energy spectra from laser-driven multielectron atoms. We have derived the EOMs for photoelectron momentum amplitudes of each orbital. To obtain PES in systems described within the MCSCF framework, the photoelectron reduced density matrix has been introduced, whose diagonal elements in the momentum space correspond to PES. Since one of the biggest benefits of tSURFF is no need to hold the complete wave function within the simulation box, it allows combined use of an efficient absorbing boundary such as irECS [17,18].
We have applied the present implementation to Ne and Ar atoms subject to attosecond XUV pulses and intense visible laser pulses, respectively. We have demonstrated converged calculation of ATI spectra from Ar including electronic correlation, which would be prohibitive without tSURFF and irECS. It is expected that this achievement enables direct comparison with experiments and precise prediction of high-field and ultrafast phenomena.
While we have presented the application of tSURFF to TD-MCSCF methods in this study, it is straightforward to extend it to other multielectron ab initio methods using time-dependent orbitals such as the time-dependent optimized coupled-cluster method [28] and TD-MCSCF methods including nuclear dynamics [29]. Such applications would enable us to compute ARPES from even more complicated systems and processes.