{'fraction_non_alphanumeric': 0.038064400071955384, 'fraction_numerical': 0.013563590573844217, 'mean_word_length': 4.538155010958358, 'pattern_counts': {'":': 1, '<': 1, '<?xml version=': 0, '>': 0, 'https://': 0, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 2, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
Calculation of photoelectron spectra within the time-dependent configuration interaction singles scheme

We present the extension of the time-dependent configuration interaction singles (TDCIS) method to the computation of the electron kinetic-energy spectrum in photoionization processes. Especially for strong and long ionizing light pulses the detection of the photoelectron poses a computational challenge because propagating the outgoing photoelectron wavepacket requires large grid sizes. Two different methods which allow for the extraction of the asymptotic photoelectron momentum are compared regarding their methodological and computational performance. The first method follows the scheme of Tong et al. \cite{tong} where the photoelectron wavefunction is absorbed by a real splitting function. The second method after Tao and Scrinzi \cite{scrinzi} measures the flux of the electron wavepacket through a surface at a fixed radius. With both methods the full angle- and energy-resolved photoelectron spectrum is obtained. Combined with the TDCIS scheme it is possible to analyze the dynamics of the outgoing electron in a channel-resolved way and, additionally, to study the dynamics of the bound electrons in the parent ion. As an application, one-photon and above-threshold ionization (ATI) of argon following strong XUV irradiation are studied via energy- and angle-resolved photoelectron spectra.


I. INTRODUCTION
With the development of new light sources such as free-electron lasers (FELs) and attosecond laser sources the interest in strong-field physics and multiphoton processes has grown, because they provide the experimental means to control and image atomic and molecular systems and to test theoretical predictions of nonlinear processes [3][4][5][6][7]. The photon energies of FELs extend from the UV to the X-ray range, and the intensities are such that they permit the investigation and control of inner-shell processes, Auger decay or above-threshold ionization (ATI) [8,9]. Above-threshold ionization, first observed in the 1970s [10], is a highly nonlinear phenomenon in which more photons are absorbed than are needed for ionization. The pulse durations of FELs can be as short as a few femtoseconds [11]. With these pulse properties typical atomic timescales which extend from a few attoseconds to tens of femtoseconds can be accessed in order to study electronic dynamics in atoms, molecules and clusters [12][13][14]. In the strong-field regime, multiphoton processes play a significant role, especially if the photon energies lie in the UV to X-ray range [15][16][17]. In general, in this frequency range a diversity of processes must be faced. The removal of a deep inner-shell electron is followed by various processes depending on the atomic states and the photon energy [18]. If the laser pulse is strong enough, also multiphoton inner-shell ionization [19] as well as ATI processes can occur.
Experimentally, photoelectron spectroscopy is a powerful tool to analyze and quantify the processes that happen due to the irradiation of complex systems [20,21]. For instance, in early experiments with intense light sources in the 1980s the angular distribution in ATI of xenon was measured [22] in order to understand the ATI phenomenon. Synchrotron radiation was used to obtain high-quality angular distributions of electrons in photoionization of atoms [23,24]. Also in recent experiments photoelectron spectroscopy has been used to reveal decay mechanisms and multiphoton excitations in deep shells of atoms [25] and to understand the origin of the low-energy structure in strong-field ionization [26,27].
The process of photoionization has been studied extensively [28][29][30][31][32], e.g., in argon or xenon [33,34]. In the weak-field limit, where the light-matter interaction can be treated perturbatively, the photoelectron spectrum has been calculated with methods that also include correlation effects. Most prominent examples are post-Hartree-Fock methods that use reference states, e.g. correlation methods like the configuration interaction [29,35], the coupled clusters method [36] and the random-phase approximation [37][38][39]. Furthermore, approaches constructing continuum wavefunctions, like R-matrix theory [40], have been applied to calculate photoionization cross sections [41] and photoelectron angular distributions [42]; taking into account the interaction of the liberated electron with other atomic orbitals has led to the explanation of the giant dipole resonance of the 4d subshell in xenon [29].
In the strong-field regime the description of the ionized wavepacket is challenging due to the nonperturbative interaction between the electrons and the light pulse. Therefore, the calculation of photoelectron spectra is numerically more demanding than in the weakfield limit. Furthermore, many-body processes are often neglected in the strong-field regime and single-active electron (SAE) approaches have become a standard tool [1,[43][44][45] where correlation effects are omitted. Nevertheless, recently, extensions to many-body dynamics have been presented, e.g., R-matrix theory [46][47][48][49][50], two-active electron [51][52][53] and timedependent restricted-active-space configuration interaction theory [54,55].
Generally, the calculation of the photoelectron spectrum can be done after the pulse is over by projecting the photoelectron wavepacket onto the eigenstates of the field-free continuum.
However, this approach requires large numerical grids and its application is very limited even in the SAE cases. For this reason, new methods were developed to calculate the spectrum using wavepacket information in a fixed spatial volume much smaller than the volume that would be needed to fully encapsulate the wavepacket at the end of the strong-field pulse.
There exist several approaches to overcome the obstacle of large grids, e.g., by measuring the electronic flux through a sphere at a fixed radius [56] or splitting the wavefunction into an internal and an asymptotic part [57,58] where the latter is then analyzed to yield the spectrum. The first implementation of the flux method in the strong-field case is the time-dependent surface flux ("tsurff") method introduced by Tao and Scrinzi [2]. It has recently been extended to the description of dissociation in molecules [59]. Tong et al. [1] applied the splitting approach to strong-field scenarios. With both methods doubledifferential photoelectron spectra can be calculated.
Our method for treating the electron dynamics within atoms is based on the timedependent configuration interaction singles (TDCIS) scheme [60]. The Schrödinger equation is solved exactly by wavepacket propagation in the configuration interaction singles (CIS) basis. The TDCIS approach [61] includes interchannel coupling and allows investigating the wavepacket dynamics and, in particular, the impact of correlation effects between the photoelectron wavepacket and the remaining ion as discussed e.g. in Refs. [62][63][64]. It is versatile with respect to the electric field properties (also multiple pulses can be chosen) and it has proven especially successful for strong-field studies [63][64][65][66].
In Sec. II, we present the theoretical details of how the photoelectron spectrum is obtained within the TDCIS scheme: The wavefunction splitting method [1] is described in Sec. II B and the time-dependent surface flux method [2] in Sec. II C. In Sec. III the two methods are analyzed with respect to their efficiency within TDCIS and are compared briefly. As an application we calculate and study the angle-and energy-resolved photoelectron spectrum of argon irradiated by strong XUV radiation. A summary and short outlook in Sec. IV conclude the article. Atomic units are used throughout except otherwise indicated.

A. Theoretical background
The time dependent Schrödinger equation of an N-electron system is given by Considering linearly polarized light, the Hamiltonian takes the form where A(t) is the vector potential 1 . Here,Ĥ(t) is the full N-electron Hamiltonian,Ĥ 0 = T +V nuc +V MF − E HF contains the kinetic energyT , the nuclear potentialV nuc , the potential at the mean-field levelV MF and the Hartree-Fock energy E HF ,Ĥ 1 = 1 |r 12 | −V MF describes the Coulomb interactions beyond the mean-field level, andˆ p · A(t) is the light-matter interaction within the velocity form in the dipole approximation.
Within the CIS approach only one-particle-one-hole excitations |Φ a i with respect to the Hartree-Fock ground state |Φ 0 are considered. Therefore, the wavefunction (now omitting the superscript N) is expanded in the CIS basis as where the index i symbolizes an initially occupied orbital and a denotes an unoccupied (virtual) orbital to which the particle can be excited: pσ andĉ pσ create and annihilate electrons, respectively, in the spin orbitals |ϕ pσ . The total spin is not altered in the considered processes (S = 0), so that only spin singlets occur. Therefore, and for the sake of readability, we drop the spin index and treat the spatial part of the orbitals |ϕ p . Inserting the wavefunction expansion (3) into the Schrödinger equation (1) and projecting onto the states |Φ 0 and |Φ a i yields the following equations of motion for the expansion coefficients α a i (t): where ε p denotes the energy of the orbital |ϕ p (Ĥ 0 |ϕ p = ε p |ϕ p ). As introduced in Ref. [67], for each ionization channel all single excitations from the occupied orbital |ϕ i may be collected in one "channel wavefunction": These channel wavefunctions may now be used to calculate all quantities in a channelresolved manner. In this way, effectively one-particle wavefunctions are obtained, which will be used in the following to derive the formulae for the photoelectron spectra. A detailed description of the TDCIS method can be found in Refs. [60,67]. The TDCIS method provides the coefficients of the wavefunction in the CIS basis, which are propagated in time.
During the propagation, quantities that are needed for the calculation of the photoelectron spectrum are prepared using the channel wavefunction coefficients. After the propagation, these quantities are then used in the subsequent analysis step to determine the spectral components of the channel wavefunctions. At the end, an incoherent summation over all ionization channels is performed to obtain the photoelectron spectrum. The two analysis methods are described in the following.

B. Wavefunction splitting method
We describe the concrete implementation of the splitting method introduced by Tong et al. in Ref. [1] within our time-dependent propagation scheme. A real radial splitting function of the formŜ is used to smoothly split the channel wavefunction (5). The parameter r c denotes the radius where the splitting function is centered, and ∆ is a "smoothing" parameter controlling the slope of the function. At the first splitting time step t 0 the channel wavefunction is split into two parts (for each channel i): |χ i,in (t) is the wavefunction in the inner region 0 < r r c and |χ i,out (t) is the wavefunction in the outer region r c r ≤ r max . Then, the following procedure is performed at t 0 : The outer part of the wavefunction |χ i,out (t 0 ) is analytically propagated to a long time T after the laser pulse is over using the Volkov HamiltonianĤ V (τ ) with the time propagator under the assumption that far from the atom the electron experiences only the laser field and not the Coulomb field of the parent ion. It is also assumed that, at the splitting radius, the electron is sufficiently far away to not return to the ion.
The inner part of the wavefunction |χ i,in (t 0 ) is propagated on a numerical grid using the full CIS Hamiltonian [see Eqs. (4a) and (4b)]. For the splitting function the ratio r c /∆ ≫ 1 must be chosen such that the ground state |Φ 0 is not affected by the splitting:Ŝ|Φ 0 = 0.
At the next splitting time t 1 the inner part of the wavefunction which was propagated from t 0 to t 1 is split again. Thus, the following prescription is obtained: This is now repeated for every splitting time t j , until all parts of the electron wavepacket that are of interest have reached the outer region. Each |χ i,out (t j+1 ) is again propagated Computationally, |χ i,out (t j ) is initially expressed in the CIS basis. For this purpose, we define new expansion coefficients for the outer wavefunction and express the wavefunction in the outer region as |χ i,out (t j ) = a β a i (t j )|ϕ a . During the propagation, at every splitting time step t j , which can be -and for computational efficiency should be-a multiple of the actual propagation time step, the splitting functionŜ is applied and the expansion coefficients (10) are calculated and stored. Since the outer wavefunction is split from the inner part and treated analytically, the grid size needed for the description of the wavefunction is automatically reduced. Later, when the spectrum is calculated, the coefficients β a i are inserted and used for the analysis. The Volkov states |Ψ V p ≡ | p V are eigenstates of the Volkov Hamiltonian and form a basis set in which the channel wavepacket at time T can be expanded: In the velocity form the Volkov states are nothing but plane waves Ψ V p ( r) = (2π) −3/2 e i p· r . The photoelectron spectrum is obtained by calculating the spectral components of the outer wavefunction. For this purpose, the following coefficients are evaluated: First, we calculate the c i ( p, t j ) for each splitting time t j where the orbital is now explicitly given in the spatial representation by r|ĉ † a |0 = r|ϕ a = ϕ a ( r) = u na,la (r) r Y la,ma (Ω r ) and, thus, possesses a radial and an angular part. We use the multipole expansion for the exponential function where j l (pr) denotes the spherical Bessel function of order l. The orthonormality relations of the spherical harmonics reduce the three-dimensional integrals in Eq. (13) to one-dimensional radial integrals. Finally, propagating to a long time T after the pulse, we obtain the coeffi- These coefficients can be used to calculate the angle and energy distribution of the ejected electron because at time T the canonical momentum equals the kinetic momentum. One can choose now a homogeneous momentum grid and calculate these coefficients for each splitting time step. In order to obtain the full electron wavepacket at time T all contributions from splitting times t j must be summed up coherently to obtain the coefficientsC i ( p ) in Eq. (11) for each ionization channel i. Then, incoherent summation over all possible ionization channels yields the photoelectron spectrum: The extra factor of p results from the conversion from the momentum to the energy differential. As long as the time T is chosen to be after the pulse the result is T independent.
Of course, one needs to choose a sufficiently large T such that the parts of the electron wavefunction that one wants to record have entered the outer region and can be analyzed.
C. Time-dependent surface flux method ("tsurff ") The second method for the calculation of photoelectron spectra is based on the approach presented by Tao and Scrinzi in Ref. [2] where it was used to calculate strong-field infrared photoionization spectra in combination with infinite-range exterior complex scaling [68].
In this approach the electron wavefunction is analyzed during its evolution when crossing the surface of a sphere of a given radius r c . Again, it is assumed that the wavefunction can be split into two parts: One part is bound to the atom and is a solution to the full Hamiltonian, the other part can be viewed as free from the parent ion and is a solution to the Volkov Hamiltonian. Therefore, the method also relies conceptually on a splitting procedure. Nevertheless, and in contrast to the splitting method, the wavefunction is not altered in this process. As above, the key idea is to obtain the spectral components of the wavefunction by projecting onto plane waves.
The surface radius r c is chosen such that the electron can be considered to be free, which vanishes for r ≤ r c . Thus, the photoelectron spectrum is the sum, over all channels, Here, the Heaviside step function θ enters (we adopted the notation by Tao and Scrinzi [2]). In order to avoid the need for a representation of χ i,out ( r, T ) at large r (because T is large, a fast electron moves far out during this time), this 3D-integral is converted into a time integral involving the wavefunction only at r = r c . For that, we must know the time evolution of the asymptotic part of the wavefunction after it has passed the surface.
Inserting the Schrödinger equation where necessary and using the Volkov solutions in the velocity form outside the sphere with radius r c we obtain [2] The commutator, which vanishes everywhere except at r = r c , is easily evaluated in polar coordinates (assuming linear polarization) and we obtain More details can be found in Ref. [2] as well as in Ref. [69]. We shuffle the derivative in the first operator term to the left, via integration by parts, and obtain the operator where " ← − ∂ " means that the derivative acts to the left on the Volkov state and " ∂" means that the derivative acts to the right on the channel wavefunction. In order to implement this operator acting on the channel wavefunctions we have to calculate also the first derivative of the wavefunctions with respect to r at the radius r c . After the propagation, during which the coefficients of the channel wavefunctions χ i (r c , t) as well as of their first derivatives [∂ r χ i (r, t)| r=rc ] have been calculated, the expression (19) can be computed. Since we introduce the multipole expansion [see Eq. (14)] for the Volkov states we have to calculate also derivatives of the spherical Bessel functions at the radius r c . This calculation is performed during the analysis step for each angular momentum l. In the last term we express the cosine as a spherical harmonic and use the identity for the integral over three spherical where the Clebsch-Gordan coefficients are given by C l 3 m 3 l 1 m 1 ,l 2 m 2 = l 1 m 1 , l 2 m 2 |l 3 m 3 . Thus, we obtain the spectral components in their final form: C l,ma la,ma;1,0 C l,0 la,0;1,0 Y l,ma (Ω p ) , where j ′ la (pr c ) = ∂ z j la (z) z=prc and u ′ na,la (r c ) = ∂ r u na,la (r) r=rc . Although this expression may seem fairly complicated, it involves only quantities evaluated at one single radius r = r c .
The photoelectron spectrum is then obtained as where, as in the splitting method, the incoherent sum over all ionization channels i is performed. In the present implementation, a complex absorbing potential (CAP) absorbs the wavefunction near the end of the numerical grid [70][71][72]. We use a CAP of the form W (r) = θ(r − r CAP )(r − r CAP ) 2 , where θ is again the Heaviside step function and r CAP is the radius where the CAP starts absorbing. It is added to the Hamiltonian in Eq. (2) in the form −i ηŴ , where η is the CAP strength. As will be discussed in Sec. III B the absorption via a CAP has to be optimized carefully, because reflections from the end of the numerical grid as well as from the CAP itself have to be minimized in order to obtain an accurate photoelectron spectrum.

III. APPLICATION: ARGON UNDER STRONG XUV RADIATION
With two methods for the calculation of photoelectron spectra implemented in TDCIS, we investigate one-photon and above-threshold ionization processes of argon in the XUV regime. Motivated by a recent experiment carried out at the free-electron laser facility FLASH in Hamburg [73] we assume a photon energy of 105 eV, which is far above the threshold for the ionization out of the 3p and 3s subshells. In the following we examine the functionality of the splitting and the surface flux methods by means of the specific example of ionization of argon in the XUV.

A. Wavefunction splitting method
In the splitting method, three parameters have to be adjusted: the splitting radius, the smoothness of the splitting function and the rate at which the absorption is applied. A first criterion for verifying that the absorption through the masking function is performed correctly is the comparison of the total ground state population obtained via splitting with the population obtained with the CAP. We use a Gaussian pulse with 9 × 10 13 Wcm In the studied parameter cases the agreement between the splitting results and that CAP result is better than 3 × 10 −3 relative difference (choosing, e.g., r max = 150 a.u., r c = 80 a.u., ∆ = 10 a.u., and varying the splitting time step between 0.2 a.u. and 10 a.u.). With more frequent absorption the agreement gets slightly better.
Analyzing the splitting method, we find that the splitting radius r c and the smoothing parameter ∆ can be varied rather freely without changing the (physical) spectrum. Although the total radial grid size can be chosen as small as 100 a.u. (cf. Fig. 1) we choose also a larger radial grid extension with r max = 250 a.u. and vary the splitting radius in the wide range from 80 a.u. to 230 a.u. Exemplarily, results in the direction θ = 0 (along the XUV polarization axis) when the radial grid size and splitting radius are varied are shown in Fig. 1(a). The spectrum shows the one-photon absorption peaks at the energy corresponding to the difference between photon energy and binding energy of the corresponding orbital (3s and 3p, respectively). The second part of the spectrum, in Fig. 1(b), is separated from the first part by the photon energy and is, therefore, attributed to above-threshold Since the splitting radius is not very crucial for the spectrum we proceed to the variation of the other parameters. Spectra for various smoothing parameters ∆ are shown in Fig. 2.
Here, the radial grid size is kept fixed at r max = 150 a.u., the splitting radius is 80 a.u., and absorption is performed every 10 a.u. of time. The physical peaks are reproduced correctly for all ∆, the noise amplitude, however, is changing. A value of ∆ = 10 seems to be the optimum, for ∆ = 15 the amplitude of unphysical peaks is higher, while the steeper slope corresponding to ∆ = 5 a.u. produces higher oscillations near the physical peaks, which should be avoided. The method is particularly sensitive to the splitting rate, i.e., how often the splitting is applied. The more frequently the splitting function is applied the less noise is obtained. This is shown in Fig. 3, where only the splitting time step is varied, while the radial grid size is kept constant at 150 a.u., the splitting radius is set to 80 a.u. and From the derivation of the splitting method in subsection II B it can be seen that the electron spectrum is normalized to the total ionization probability (because only normalized wavefunctions are used). Therefore, the integrated spectrum represents a good measure of the quality of the spectrum; the fully integrated spectrum must agree with the total ionization probability. This can be verified for different parameter specifications. The relative difference to the CAP result is found to be smaller than 2% in all studied parameter cases.
In the following we apply a strong XUV pulse centered at 105 eV with 0.7 eV bandwidth (FWHM), which corresponds to a Fourier-transform limited pulse with 108 a.u. (2.6 fs) duration. The peak intensity of the pulse is 1.0 × 10 15 Wcm −2 . In the upper left panel of Fig. 4 the full angle-and energy-resolved photoelectron spectrum of argon after one-photon absorption is shown. The angle denotes the direction with respect to the polarization axis.
The peaks arise from ionization out of the 3s and 3p shells, respectively. The lower left panel shows the corresponding ATI spectrum. As expected, the angular distributions feature the corresponding contributions from the different channels, which can be seen on the right in the four cuts along the fixed peak energies: The one-photon peak from the 3s shows a pwave character, the 3p peak has both an s-and a d-wave contribution. Analogously, the two-photon peak of 3s exhibits an s-and d-wave character and the 3p peak a p-and f-wave character. In a nutshell, the splitting method is a well-working tool for the calculation of photoelectron spectra although the requirement to optimize three parameters (r c , ∆, dt spl ) can render calculations time-consuming. fs pulse duration. The grid size is r max = 100 a.u., r c = 20 a.u., ∆ = 5 a.u., and dt spl = 10 a.u.
The angle denotes the direction with respect to the polarization axis of the pulse. The angular distribution reflects the change in angular momentum by multiphoton absorption.
dependent [70] the tsurff spectrum can be optimized only for a limited energy range. In Fig. 6 the energy spectrum for θ = 0 is shown for different CAP strengths η. It is clear that reflections from the CAP as well as from the end of the radial grid leave a trace in the spectrum. A weak CAP cannot fully absorb a fast electron before the end of the radial grid.
On the other hand, a strong CAP will reflect the electron. For the kinetic energies of the electrons considered here the optimized CAP parameter lies at a value of about 10 −3 . Of course, the other parameter that must be optimized is the CAP radius r CAP . We find that the optimum is an absorption range of r max −r CAP = 30 a.u. For tsurff also the distance of r c to r CAP plays a role. In Fig. 7 the spectrum is shown for different r CAP − r c values. For However, for both methods the spectrum has a very low noise level, up to ten orders of magnitude smaller than the physical signal. Summarizing, the tsurff-method is in principle applicable with a CAP, although it requires a good quality absorption over a broad energy range. Qualitatively, the tsurff method reproduces exactly the same results as obtained with the splitting method. r CAP -r c =0 r CAP -r c =10 a.u. r CAP -r c =20 a.u. tsurff: r CAP =220, r c =200, η=1×10 -3 , t max =1000 splitting: r c =200, ∆=10, dt spl =0.2, t max =400 FIG. 8: The argon photoelectron spectra obtained with the splitting and the tsurff methods for a pulse with 105 eV photon energy, 9 × 10 13 Wcm −2 intensity and 1.2 fs duration along the polarization direction are compared. The radial grid size is r max = 250 a.u. for both methods.

IV. CONCLUSION
We have implemented two computational methods for the calculation of photoelectron spectra within the TDCIS scheme. Both methods can produce reasonable and quantitative energy-and angle-resolved spectra within our model. Subshell ionization can be quantified. We have applied and compared these methods for the high-intensity XUV regime.
Advantages of the splitting method are the good absorption characteristics through the splitting function and the short propagation time that is needed. A disadvantage is the long evaluation time of the radial integrals. The tsurff method needs a longer propagation time. However, the calculation of the photoelectron spectrum in the analysis step is much faster than with the splitting method due to the evaluation at one point. The comparison of the two methods shows that, in principle, the same spectra can be obtained after the appropriate optimization of the computational parameters. Although our application in the present work focuses on the XUV range, it is of course possible to study also processes in the strong-field regime in the infrared range by analyzing the photoelectron spectrum.
Interesting applications arise from the fact that information about the coherence and the entanglement of the ionic state and the photoelectron can be extracted from the properties of the outgoing wavepacket.

V. ACKNOWLEDGMENTS
This work has been supported by the Deutsche Forschungsgemeinschaft under Grant No.