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

TL;DR

This paper introduces an extension of the tSURFF method for multielectron systems within the TD-MCSCF framework, enabling efficient and accurate computation of photoelectron spectra under intense laser fields.

Contribution

The authors develop a numerical implementation of tSURFF for TD-MCSCF, allowing direct extraction of photoelectron spectra from multielectron wave functions with reduced computational cost.

Findings

Successfully applied to Ne photoionization with attosecond XUV pulses.

Accurately computed above-threshold ionization of Ar with mid-infrared lasers.

Demonstrated improved efficiency and accuracy over traditional projection methods.

Abstract

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.

Figures

Equations

$$\Psi (t)=\sum _I{\Phi }_I(t)C_I(t),$$

$$H(t)=H_1(t)+H_2,$$

$$H_1(t)$$

$$=\sum _{i=1}^N\left(-\frac{{\nabla }_i^2}{2}-\frac{Z}{|{\mathbf{r}}_i|}-i\mathbb{A}(t)\cdot {\nabla }_i\right),$$

$$H_2$$

$$\Psi (t)=\hat{A}\left[{\Phi }_{\mathrm{fc}}{\Phi }_{\mathrm{dc}}(t)\sum _I{\Phi }_I(t)C_I(t)\right],$$

$$i\frac{d}{dt}{C}_I(t)=\sum _J⟨{\Phi }_I⟩\hat{H}-\hat{R}{\Phi }_J$$

$$i\frac{d}{dt}|{\psi }_p⟩=\hat{h}|{\psi }_p⟩+\hat{Q}\hat{F}|{\psi }_p⟩+\sum _q|{\psi }_q⟩R_p^q,$$

$$\hat{F}|{\psi }_p⟩=\sum _{oqsr}(D^{-1}{)}_p^o{P}_{or}^{qs}{\hat{W}}_s^r|{\psi }_q⟩,$$

$$W_s^r\left(\mathbf{r}\right)=\int d{\mathbf{r}}^{\prime }\frac{{\psi }_r^{\ast }({\mathbf{r}}^{\prime }){\psi }_s({\mathbf{r}}^{\prime })}{|\mathbf{r}-{\mathbf{r}}^{\prime }|}.$$

$$R_p^q=i⟨{\psi }_q⟩\dot{{\psi }_p}-h_p^q,$$

$$i\frac{\partial }{\partial t}\Psi (\mathbf{r},t)=h_1(t)\Psi (\mathbf{r},t),$$

$$h_1(t)=-\frac{1}{2}{\nabla }^2+V(\mathbf{r})-i\mathbf{A}(t)\cdot \nabla ,$$

$$\rho (E,t)=\int d\Omega |a(\mathbf{k},t){|}^2|\mathbf{k}{|}^2(E=|\mathbf{k}{|}^2/2),$$

$$a(\mathbf{k},t)$$

$$\equiv \int {\chi }_{\mathbf{k}}^{\ast }(\mathbf{r},t)\theta (|\mathbf{r}|-R_s)\Psi (\mathbf{r},t)d^3\mathbf{r},$$

$$={\int }_{r>R_s}{\chi }_{\mathbf{k}}^{\ast }(\mathbf{r},t)\Psi (\mathbf{r},t)d^3\mathbf{r}.$$

$$h_{\mathrm{s}}=-\frac{1}{2}{\nabla }^2-i\mathbf{A}(t)\cdot \nabla .$$

$$-i\frac{\partial }{\partial t}a(\mathbf{k},t)$$

$$P(\mathbf{r},{\mathbf{r}}^{\prime })=\sum _{pq}⟨\mathbf{r}|\theta (R_s)|{\psi }_p⟩D_q^p⟨{\psi }_q|\theta (R_s)|{\mathbf{r}}^{\prime }⟩.$$

$$\begin{gathered} \tilde{P}(\mathbf{k},{\mathbf{k}}^{\prime })=\int d\mathbf{r}d{\mathbf{r}}^{\prime }⟨{\chi }_{\mathbf{k}}(t)⟩{\mathbf{r}}^{\prime }P(\mathbf{r},{\mathbf{r}}^{\prime })⟨\mathbf{r}⟩{\chi }_{\mathbf{k}}(t) \\ =\sum _{pq}⟨{\chi }_{\mathbf{k}}(t)|\theta (R_s)|{\psi }_p⟩D_q^p⟨{\psi }_q|\theta (R_s)|{\chi }_{{\mathbf{k}}^{\prime }}(t)⟩. \end{gathered}$$

$$\rho (E)=\int d\Omega \tilde{P}(\mathbf{k},\mathbf{k})|\mathbf{k}{|}^2.$$

$$a_p(\mathbf{k},t)=⟨{\chi }_{\mathbf{k}}(t)⟩\theta (R_{\mathrm{s}}){\psi }_p,$$

$$\frac{d}{dt}|{\psi }_p⟩\simeq -i\left[{\hat{h}}_{\mathrm{s}}|{\psi }_p⟩-\sum _q|{\psi }_q⟩\left\{⟨{\psi }_q|\hat{F}|{\psi }_p⟩-R_p^q\right\}\right],$$

$$-i\frac{\partial }{\partial t}{a}_p(\mathbf{k},t)=$$

$$+$$

$${\psi }_p(r,\theta ,\varphi )=\sum _{kl}{c}_{kl}^p{f}_k(r)Y_{l{m}_p}(\theta ,\varphi ),$$

$${\chi }_{\mathbf{k}}(\mathbf{r},t)$$

$$\times \sum _{lm}4\pi i^l{Y}_{lm}^{\ast }({\theta }_k,{\psi }_k)j_l(kr)Y_{lm}(\theta ,\psi ),$$

$$\Lambda (\mathbf{k},t)={\int }_0^t\frac{1}{2}[\mathbf{k}-\mathbf{A}(\tau ){]}^{2}d\tau .$$

$$[h_{\mathrm{s}},\theta (R_{\mathrm{s}})]$$

$$-\frac{\delta (r-R_{\mathrm{s}})}{r}+i{A}_z(t)\cos (\theta )\delta (r-R_{\mathrm{s}}),$$

$$\braket{\chi_{\bm{k}}(t)}{[h_{\rm s},\theta(R_{\rm s})]}{\psi_{p}(t)}=\frac{\exp{i\Lambda(\bm{k},t)}}{(2\pi)^{3/2}}4\pi\sum_{l}(-i)^{l}\Big{[}Y_{lm_{p}}(\theta_{k},\psi_{k})\frac{R_{\rm s}^{2}}{2}\{j^{*}_{l}{}^{\prime}(kR_{\rm s})g^{p}_{l}(R_{\rm s})-j^{*}_{l}(kR_{\rm s})g^{p}_{l}{}^{\prime}(R_{\rm s})\}\\ +iA_{z}(t)\sum_{l^{\prime}}(-i)^{l^{\prime}}Y_{l^{\prime}m_{p}}(\theta_{k},\psi_{k}))R_{\rm s}^{2}j^{*}_{l}(kR_{\rm s})g^{p}_{l}(R_{\rm s})\alpha_{l^{\prime}lm}\Big{]},$$

$$\begin{gathered} \mathbf{a}(\mathbf{k},t+\Delta t)=\mathbf{a}(\mathbf{k},t)\exp [i\Delta tX(t)] \\ +\mathbf{S}(t)\frac{\exp [i\Delta tX(t)]-1}{X(t)}, \end{gathered}$$

$$X_{qp}(t)$$

Full text

††preprint: Preprint

I introduction

The rapid progress in experimental techniques for high-intensity, ultrashort optical pulses, has lead to the advent and development of strong-field physics and attosecond science Agostini and DiMauro (2004); Krausz and Ivanov (2009), 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 time-dependent multiconfiguration self-consistent field (TD-MCSCF) methods have been developed Ishikawa and Sato (2015); Zanghellini et al. (2003); Kato and Kono (2004); Caillat et al. (2005); Nguyen-Dang et al. (2007); Miyagi and Madsen (2013); Sato and Ishikawa (2013); Miyagi and Madsen (2014); Haxton and McCurdy (2015); Sato and Ishikawa (2015) which express the total wave function as a superposition,

[TABLE]

of the Slater determinants (or electronic configurations) $\Phi_{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 Zanghellini et al. (2003); Kato and Kono (2004); Caillat et al. (2005), the computational cost increases factorially with the number of electrons. To overcome this difficulty, we have recently formulated and numerically implemented the time-dependent complete-active-space self-consistent-field (TD-CASSCF) method Sato and Ishikawa (2013), and even more general and less demanding time-dependent occupation-restricted multiple-active-space (TD-ORMAS) method Sato and Ishikawa (2015). 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 time-dependent 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 Tao and Scrinzi (2012), 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 Tao and Scrinzi (2012) and then applied to multielectron simulations with, e.g., the hybrid antisymmetrized coupled channels method Majety et al. (2015), the time-dependent configuration interaction singles method Karamatskou et al. (2014), and the time-dependent density functional theoryWopperer et al. (2017).

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) Scrinzi (2010); Orimo et al. (2018) 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_{\uparrow}$ up-spin and $N_{\downarrow}$ down-spin electrons ($N=N_{\uparrow}+N_{\downarrow}$) in a laser field ${\bf E}(t)$. Its dynamics is described within the dipole approximation by the Hamiltonian,

[TABLE]

with,

[TABLE]

where $Z$ denotes the atomic number and $\bm{A}(t)=-\int_{-\infty}^{t}\bm{E}(t^{\prime})dt^{\prime}$ the vector potential. We use the velocity gauge, since ECS works only with it, and not with the length gauge McCurdy et al. (1991).

In the TD-CASSCF and TD-ORMAS methods, the total wave function is given by,

[TABLE]

where $\hat{A}$ denotes the antisymmetrization operator, $\Phi_{\rm fc}$ and $\Phi_{\rm dc}$ the closed-shell determinants formed with $n_{\rm fc}$ frozen-core and $n_{\rm dc}$ dynamical-core orbitals, respectively, and $\{\Phi_{I}\}$ the determinants constructed from $n_{\rm a}$ active orbitals. Whereas the active electrons are fully correlated among the active orbitals in the TD-CASSCF method, the TD-ORMAS method Sato and Ishikawa (2015) 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 $\{\psi_{p}\}$ are derived on the basis of the time-dependent variational principle (TDVP) Moccia (1973); Sato and Ishikawa (2015) and read,

[TABLE]

where $\hat{Q}=1-\sum_{q}\Ket{\psi_{q}}\Bra{\psi_{q}}$ the projector onto the orthogonal complement of the occupied orbital space. $\hat{F}$ is a non-local operator describing the contribution from the interelectronic Coulomb interaction, defined as

[TABLE]

where $D$ and $P$ are the one- and two-electron reduced density matrices, and $\hat{W}^{r}_{s}$ is given, in the coordinate space, by

[TABLE]

The matrix element $R^{q}_{p}$ is given by,

[TABLE]

with $h^{q}_{p}=\braket{\psi_{q}}{\hat{h}}{\psi_{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. Sato and Ishikawa (2015), 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 Tao and Scrinzi (2012) for single-electron systems governed by TDSE,

[TABLE]

where $V(\bm{r})$ denotes the nuclear potential. The photoelectron momentum amplitude (PMA) $a(\bm{k},t)$ for momentum $\bm{k}$ and PES

[TABLE]

can in principle be approximately calculated by projecting the outgoing wave packet residing outside a given radius $R_{s}$ onto the plane wave $(2\pi)^{-3/2}\exp(i\bm{k}\cdot\bm{r})$ at a time $t$ sufficiently after the pulse. Then, PMA is given by,

[TABLE]

where $\chi_{\bm{k}}(\bm{r},t)$ denotes the Volkov wave function and $\theta(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 $|\bm{r}|>R_{\rm s}$ are evolved by the same nuclear-potential-free Hamiltonian

[TABLE]

By taking time derivative of Eq. (14), we obtain the EOM of the momentum amplitude,

[TABLE]

As shown in Ref. Tao and Scrinzi (2012), since all the terms appearing in the commutator $[h_{\rm s},\theta(R_{\rm s})]$ contain delta functions $\delta(r-R_{\rm s})$ [see Eq. (IV.3) below], wave functions and their spatial derivative only on the surface $|\bm{r}|=R_{\rm 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.

IV Application of the tSURFF method to the TD-ORMAS simulations

IV.1 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 $|{\bf r}|>R_{s}$ , the single particle reduced density matrix of a photoelectron in the coordinate space can be defined as,

[TABLE]

In the momentum space, its elements are given by,

[TABLE]

The diagonal part $\tilde{P}(\bm{k},\bm{k})$ is interpreted as photoelectron momentum distribution, and, then, PES is given by,

[TABLE]

A similar definition of the PRDM as Eq. (20) and the direct projection onto scattering states were used to compute photoelectron spectrum in Ref. Omiste et al. (2017).

IV.2 EOMs of momentum amplitudes of orbitals

In this subsection, we derive the EOM of the momentum amplitude of orbital $p$,

[TABLE]

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_{\rm s}$. Then, the orbital EOM Eq. (7) can be approximated as,

[TABLE]

so that EOMs of $a_{p}(\bm{k})$ is obtained as

[TABLE]

It should be noted that the second term in Eq. (IV.2) 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.

IV.3 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 Sato et al. (2016) 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 McCurdy et al. (2004),

[TABLE]

where $f_{k}(r)$ and $Y_{lm_{p}}(\theta,\phi)$ denote the FEDVR function and the spherical harmonics, respectively. The magnetic quantum number $m_{p}$ for orbital $\psi_{p}$ does not change due to the $z$-polarization of the laser pulse Sato et al. (2016). As an absorbing boundary, infinite-range exterior complex scaling is implemented Scrinzi (2010); Orimo et al. (2018), 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 $\bm{k}=(k,\theta_{k},\psi_{k})$ is given by,

[TABLE]

with $j_{l}(kr)$ being the spherical Bessel function of the first kind and $\Lambda(\bm{k},t)$ the Volkov phase given by,

[TABLE]

The commutator in Eq. (18) can be rewritten as,

[TABLE]

with $z$ component $A_{z}(t)$ of the vector potential $\bm{A}(t)$. Using Eqs. (IV.3) and (27) and introducing $g^{p}_{l}(r)=\sum_{k}c^{p}_{kl}f_{k}(r)$, we can decompose the first term of Eq. (IV.2) into,

[TABLE]

where $j^{*}_{l}{}^{\prime}(r)$ and $g^{p}_{l}{}^{\prime}(r)$ denote the radial derivative of $j^{*}_{l}(r)$ and $g^{p}_{l}(r)$, respectively. For the second term in Eq. (IV.2), since $\{\bra{\psi_{q}(t)}\hat{F}\ket{\psi_{p}(t)}-R^{q}_{p}\}$ is evaluated during the orbital propagation [Eq. (7)], we can reuse it.

We integrate Eq. (IV.2) with the first order exponential integrator Cox and Matthews (2002) by treating Eq. (IV.2) as simultaneous inhomogeneous linear differential equations. The evolution of a vector $\bm{a}(\bm{k},t)\equiv\{a_{p}(\bm{k},t)\}$ from the time $t$ to $t+\Delta t$ is described as

[TABLE]

where the matrix $X(t)=\{X_{qp}(t)\}$ and vector $\bm{S}(t)=\{S_{p}(t)\}$ are defined as,

[TABLE]

We compute $\exp[i\Delta tX(t)]$ by directly diagonalizing $X(t)$ at every time step, which is not demanding since the size of $X(t)$ is $[N_{\rm orb}\times N_{\rm orb}]$ and the number of orbitals $N_{\rm 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.1):

[TABLE]

and for Ar atom (Sec. V.2):

[TABLE]

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

V.1 Ne atom

We first calculate the PES of a neon atom subject to an attosecond pulse with a peak intensity of $2.5\times 10^{12}\rm{W/m}^{2}$, a wavelength of 12.398 nm corresponding to 100 eV photon energy, and $N_{\rm 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\,{\rm 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_{\text{fc}},n_{\text{dc}},n_{\text{a}})=(0,0,5)$, $(0,0,9)$, and $(1,0,8)$, where $n_{\text{fc}}$, $n_{\text{dc}}$, and $n_{\text{a}}$ are the number of frozen-core, dynamical-core, and active orbitals, respectively. Note that the first and second correspond to the time-dependent Hartree-Fock (TDHF) Dirac (1930) 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 NIST ASD Team , 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.

V.2 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\times 10^{14}\,\mathrm{W/cm^{2}}$, a wavelength of 532 nm, and a pulse width of $N_{\rm 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_{\text{a}}$ active orbitals into 4 and $(n_{\text{a}}-4)$ orbitals, as schematically illustrated in Fig. 3 for orbital decomposition $(n_{\text{c}},n_{\text{d}},n_{\text{a}})=(5,0,13)$. By setting the maximum number of electrons in the second subgroup [$(n_{\text{a}}-4)$ active orbitals] to 2 or 3, the configurations 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,

[TABLE]

where $I_{\rm p}$ is the ionization potential. The difference in peak position observed in Fig. 4 can be attributed to that in $I_{\rm 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 Kulander et al. (1992). 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_{\rm a}=25$. Especially, the TDHF method significantly underestimates the yield in the direction ($90^{\circ}$) perpendicular to the laser polarization. This indicates that electronic correlation is non-negligible in detailed discussion of ATI ARPES.

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 Scrinzi (2010); Orimo et al. (2018).

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 Sato et al. (2018) and TD-MCSCF methods including nuclear dynamics Anzaki et al. (2017). Such applications would enable us to compute ARPES from even more complicated systems and processes.

Acknowledgements.

This research was supported in part by a Grant-in-Aid for Scientific Research (Grants No. 16H03881, No. 17K05070, and No. 18H03891) from the Ministry of Education, Culture, Sports, Science and Technology (MEXT) of Japan and also by the Photon Frontier Network Program of MEXT. This research was also partially supported by the Center of Innovation Program from the Japan Science and Technology Agency (JST), by CREST (Grant No. JPMJCR15N1), JST, by Quantum Leap Flagship Program of MEXT, and by JSPS and HAS under the Japan-Hungary Research Cooperative Program. Y. O. gratefully acknowledges support from the Graduate School of Engineering, The University of Tokyo, Doctoral Student Special Incentives Program (SEUT Fellowship).