First-principles simulations for attosecond photoelectron spectroscopy based on time-dependent density functional theory

We develop a first-principles simulation method for attosecond time-resolved photoelectron spectroscopy. This method enables us to directly simulate the whole experimental processes, including excitation, emission and detection on equal footing. To examine the performance of the method, we use it to compute the reconstruction of attosecond beating by interference of two-photon transitions (RABBITT) experiments of gas-phase Argon. The computed RABBITT photoionization delay is in very good agreement with recent experimental results from [Kl\"under et al, Phys. Rev. Lett. 106 143002 (2011)] and [Gu\'enot et al, Phys. Rev. A 85 053424 (2012)]. This indicates the significance of a fully-consistent theoretical treatment of the whole measurement process to properly describe experimental observables in attosecond photoelectron spectroscopy. The present framework opens the path to unravel the microscopic processes underlying RABBITT spectra in more complex materials and nanostructures.


I. INTRODUCTION
Recent progress of laser technologies has enabled the observation of ultrafast phenomena with attosecond resolution and offered novel opportunities to directly explore real-time electron dynamics in matter [1][2][3]. Broadly speaking one can assign the available attosecond timeresolved measurement techniques to two major groups: one is based on all-optical measurements such as the attosecond transient absorption spectroscopy [4][5][6]. The other is based on photoelectron spectroscopy such as the reconstruction of attosecond beating by interference of two-photon transitions (RABBITT) [7,8] and the attosecond streaking camera [9,10].
In the past decade, attosecond transient absorption spectroscopy has been applied to atomic and molecular systems, and ultrafast electron dynamics in these relatively small systems has been intensively investigated both experimentally [4,[11][12][13] and theoretically [14,15]. Recently, this technique has been extended to solid-state materials where nonequilibrium electron dynamics has been investigated towards future applications such as petahertz devices [16][17][18]. Although attosecond spectroscopy of solid-state materials provides in principle a wealth of information on novel aspects of ultrafast dynamics, experimental results are often hard to interpret directly, because of the strong nonlinearlity of lightmatter interactions combined with the complex electronic structure of solids. To extract microscopic insight from attosecond transient absorption spectra of solids, firstprinciples simulations based on the density functional theory (DFT) [19,20] and the time-dependent density functional theory (TDDFT) [21] have played a significant role [16,17,22].
Likewise, attosecond photoelectron spectroscopy has been applied to atomic systems [23][24][25][26] as well as recently to solid-state materials [27][28][29]. However, in spite of the intensive development of several ab-initio approaches for atomic and molecular systems [30][31][32], similar approaches for solids and surfaces have not yet been established. Therefore, in order to understand experiments on such complex systems, further development of first-principles approaches is required. In this regard one promising candidate is represented by real-time electron dynamics simulations based on TDDFT. Realtime TDDFT simulations of solids have been already applied to ultrafast as well as strong-field-induced phenomena such as attosecond transient absorption spectroscopy [17], high harmonic generation [33][34][35], laserinduced damage [36], and laser-induced magnetism [37].
In a recent work [38] the authors have introduced a method to compute the angle-resolved photoelectron spectrum of solid-state systems. This method is based on the calculation of the photoelectron flux trough a closed surface that can be simulated with TDDFT in a real-space real-time implementation. Here we illustrate how this approach can be employed to calculate attosecond photoelectron spectra of finite systems. This constitutes a fundamental benchmark towards the application to solid-state materials.
For this purpose, we perform the attosecond photoelectron spectroscopy simulation of an Argon atom and compare the theoretical results with recent experiments. In particular here we focus on RABBITT experiments, since the alternative approach, attosecond streaking, provides equivalent information [39,40] and we emphasize that the TDDFT attosecond photoelectron spectroscopy can be straightforwardly applied also to the attosecond streaking technique.
RABBITT has been originally introduced for the temporal characterization of attosecond pulses [7], and has then been employed to investigate the photoionization delay in atoms and molecules [24,26,39] as well as, in recent times, in solid-state surfaces [28,41]. The RAB-BITT technique employs two laser pulses in a pumpprobe fashion: an attosecond extreme ultraviolet (EUV) pulse train is used as a pump to ionize the system while a femtosecond infrared (IR) pulse is used as a probe. This configuration is designed for experiments where the pulse train is obtained by an high harmonic generation stage seeded by the IR pulse. As the attosecond pulse train consists of a frequency comb of odd IR frequency multiples it produces an energy comb of photoelectron spectra that are shifted by IR photons. Probing the system with the delayed IR brings two adjacent photoelectron peaks in contact and forms an interference pattern that oscillates as a function of the delay. This pattern encodes information on the emission delay with attosecond resolution. In this paper we will demonstrate how the entire process can be efficiently simulated with TDDFT.
The structure of this paper is as follows: In Sec. II, we describe the theoretical and numerical methods to compute electron dynamics and photoelectron spectra based on the TDDFT. In Sec. III, we demonstrate the firstprinciples simulation for attosecond photoelectron spectroscopy and compare the theoretical results with recent experiments. We further discuss the role of a many-body effects in the photoemission process. Finally, we summarize our findings and provide some perspective for future work in Sec. IV.
Hartree atomic units ( = e = m e = 4π 0 = 1) are employed throughout the paper unless otherwise specified.

II. METHOD
The fundamental concept of TDDFT is that all physical properties of a time-dependent system can be determined through their functional dependence on the timedependent interacting many-body density [21], n(r, t) and the initial many body state, which can be disregarded if we start from the ground state. The idea of both DFT and TDDFT, is to obtain this many-body density by mapping it to the density of a fictitious auxiliary system of non-interacting electrons: the Kohn-Sham (KS) system. The dynamics of the KS system can be obtained by propagating the one-particle equations for the orbitals ϕ i (r, t) of a single Slater determinant, according to the time-dependent KS (TDKS) equations To simplify notation we here only consider systems with an even number of electrons N , so that each spatial orbital ϕ i is doubly occupied with two electrons of opposite spin. The KS Hamiltonian governing the dynamics of the orbitals in (1) is defined as: where, due to the action of the KS potential v KS , the time-dependent density n(r, t) = 2 N/2 i=1 |ϕ i (r, t)| 2 corresponds both to the real and to the KS system. The KS potential is composed of three terms. The first term is the electron-ion potential provided by the nuclei, while the second term is the electrostatic potential generated by the electronic charge density v H [n](r, t) = dr n(r , t)/|r − r |. The last term v xc [n](r, t) is the so-called exchange and correlation (xc) potential that accounts for the many-body effects deriving from the electron-electron interaction; it is a functional of the density at all times n(r, t) and, since its explicit form is unknown, it must be approximated. In this work we employ the adiabatic local-density approximation (ALDA) [42,43] which is based on the xc potential of a homogeneous electron gas evaluated with the instantaneous density in time at every point in space. In order to compensate the self-interaction error [43] of the local approximation and obtain a correct ionization potential we employ the simplest scheme based on the averaged-density self-interaction correction (ADSIC) [44].
Given the energy range of the lasers employed in the simulations, it is well justified to invoke the dipole approximation for the light-matter interaction. Under this condition, the coupling with the laser field can be expressed in the velocity gauge which amounts to modifying the kinetic operator by adding the spatially homogeneous time-dependent vector potential A(t) of the classical laser field, as in Eq. (2)[? ]. The time profile of A(t) can accommodate any linear combination of laser fields and is therefore naturally suited to describe any kind of pump-probe configuration, including the one employed for RABBITT experiments.
To obtain the photoelectron spectrum from the timedependent KS orbitals, we use the fact that it can be expressed as a flux integral of the ionization current through a closed surface. This approach is based on the t-SURFF method, first introduced by Scrinzi [45] for one-electron systems and later extended to many electrons with TDDFT [46]. According to this formulation the momentum-resolved photoelectron probability P(p), i.e. the probability to measure an electron with a given momentum p, can be expressed as whereĵ is the single-particle gauge-invariant current density operator and χ p (r, t) = (2π) − 3 2 e i(p+A(t))·r e iΦ(p,t) . The phases Φ(p, t) = t 0 dτ 1 2 p + A(t) c 2 describe Volkov waves of momentum p that are the analytical solutions of the time-dependent Schrödinger equation for free particles in a field. The bracket notation in the equation is thus used as shorthand to indicate the evaluation of the current-density operator matrix element between KS orbitals and Volkov waves. The energy-resolved photoelectron spectrum, P(E), employed to build the RAB-BITT traces can be obtained by integrating the angular dependence of P(p) as follows This approach to photoemission is particularly suited to numerical implementations where the TDKS equations (1) are solved in real-space and propagated in realtime. In our implementation the spatial coordinates are discretized on a cartesian grid with spacing ∆ = 0.3 a.u. and the equations are solved on a spherical box of radius R = 30 a.u.. The TDKS equations are propagated under the influence of a time dependent field with a time step ∆t = 0.04 a.u. starting from the ground state configuration. The photoelectron probability is calculated with (4) by collecting the flux integral calculated on a spherical surface of radius R S = 20 a.u. while the KS orbitals are propagated over time. To prevent spurious reflections from the boundaries of the simulation box we employ a complex absorbing potential (CAP) acting on the region outside the surface S with parameters tuned in such a way to be maximally efficient in the energy region where we expect photoelectrons to be mostly emitted [47]. The geometry employed in the simulations is summarized in Fig. 1.
Scheme illustrating the geometry employed to calculated the photoelectron spectrum with t-SURFF and TDDFT.
Finally, since the inner shell electrons of Argon are not expected to take significant part in the ionization dynamics, we use the Hartwigsen-Goedecker-Hutter (HGH) pseudopotential [48] that effectively accounts for the core electrons and consider explicitly only the n = 3 electrons.
All the simulations presented are carried out with the Octopus code [49].

III. RESULT
In this section, we examine the performance of the TDDFT simulation for the attosecond photoelectron spectroscopy. For this purpose, we simulate the RAB-BITT measurement processes for an Argon atom. We first explain how to simulate the entire the RABBITT measurement. Then, we compare the theoretical results with the recent experimental data [24,25]. Finally, we investigate the role of many-body effects in the RABBITT photoemission delay.

A. RABBITT spectroscopy
Here, we revisit the RABBITT pump-probe technique from a computational point of view. The RABBITT measurement is a pump-probe experiment that employs an attosecond EUV pulse train as a pump and an IR femtosecond pulse as a probe. Importantly, the attosecond pulse train is generated by high-order harmonic generation of the same IR femtosecond pulse. Therefore, the attosecond EUV pulse train consists of odd harmonics of the IR field.
In this work, we employ the following form for the femtosecond IR pulse, in the domain −T IR /2 < t < T IR /2 and zero outside.
Here ω 0 is a mean frequency of the IR pulse, and T IR is the full duration of the pulse. The maximum field amplitude E IR is related to the peak laser intensity as I IR = cE 2 IR /8π. We set ω 0 to 1.55 eV, and T IR to 30 fs. The peak intensity I IR is set to 10 11 W/cm 2 . As a corresponding attosecond pulse train, we employ the following form: in the domain −T train /2 < t < T train /2 and zero outside.
Here ω EU V is the central frequency, and T train is the full duration of the pulse train. We set ω EU V to 25ω 0 , and T train to 10 fs. We also set the peak laser intensity I EU V = cE 2 EU V /8π to 5 × 10 10 W/cm 2 . Figure 2 shows the attosecond pulse train of Eq. (7) in the time-domain (a) and the frequency-domain (b). As seen from Fig. 2 (a), several attosecond pulses follow each other in a line with equal distance in time-domain. Each pulse has about 120 attoseconds full-width halfmaximum. This train forms a comb in the frequency domain, as seen in Fig. 2 (b). We note that the comb consists of the odd order harmonics of the IR probe pulse.
We then compute the photoelectron spectrum induced by the attosecond pulse train. Figure 3 shows the photoelectron spectrum as a function of kinetic energy of emitted electrons. Since the photoelectron spectrum is computed based on the time-dependent Kohn-Sham orbitals, one can naturally decompose the signal into each orbital contribution. In Fig. 3, the contribution from the Ar 3s shell is shown as a red-solid line, while that from the Ar 3p shell is shown as a green-dashed line. One sees that the two contributions are energetically well separated because of the large difference between the ionization potentials of the 3s and 3p shells. Each contribution shows the comb structure, reflecting the frequency comb feature of the attosecond pulse train in Fig. 2 (b).
In a RABBITT experiment, the photoelectron spectrum under both the attosecond pulse train and the femtosecond IR pulse is measured. In perfect analogy, we can compute in the theoretical simulation the photoelectron spectrum under both the attosecond pulse train and the IR femtosecond laser pulse. Figure 4 shows the photoelectron spectrum from the Ar 3p shell. Red-solid line shows the photoelectron spectrum created by both the attosecond pulse train and the IR femtosecond pulse, while blue-dashed line shows the signal solely due to the pulse train. One sees that the IR pulse results in additional peaks between those peaks that were created only by the pulse train. These additional peaks originate from a twophoton absorption process: one-photon from the attosecond pulse train and the other from the IR pulse. Adding the ionization potential of the Ar 3p shell to the photoelectron kinetic energy, the absorbed photon energy can be calculated. The calculated absorbed photon energy is shown as the secondary x-aixis of Fig. 4. Each additional peak due to the IR field consists of two excitation paths: one corresponds to the EUV photon energy plus the IR photon energy, while the other corresponds to the EUV minus IR photon energy. For example, as seen in the schematic picture of Fig. 4, the additional peak at the energy of 24ω 0 is created by the following two excitation paths: One is the 23rd harmonics plus the IR photon energy, while the other is the 25th harmonics minus the IR photon energy. As discussed above, this interference between the two excitation path is the central effect used in which RABBIT spectroscopy. We next perform the RABBITT pump-probe simulations by changing the time delay between the attosecond pulse train and the IR pulse. Figure 5 shows the calculated photoelectron spectrum as a function of the time delay. One sees that the even order side bands show an oscillating feature in time delay, reflecting the interference of the two different two-photon absorption paths, which is described in the schematic picture of Fig. 4. The frequency of the oscillation is twice that of the IR frequency ω 0 . Generally, each side band has its own timedelay with respect to the IR field. Because the difference of the delay of these side bands reflects the difference of the photoionization delay in each excitation channel, the time delay in the RABBITT trace has been used to investigate the photoionization processes [24][25][26]. Since TDDFT can directly simulate the whole RABBITT experimental process and provide the resulting RABBITT trace as seen in Fig. 5, it enables us to directly compare calculated results with experimental results. In the following subsection, we demonstrate the comparison of theory with experiment.  Fig. 2.

B. Comparison with experimental results
Here, we compare the computed time delays from TDDFT simulations with the experimental results [24,25]. For this purpose, we first numerically extract the delay from the RABBITT trace in Fig. 5. To extract the delay for each side band, we average the RABBITT trace around the central frequency of the side-band with width of ω 0 /2. For example, to extract the 26th side band in Fig. 5, we average the signal between 26ω 0 ± ω 0 /4. Figure 6 shows the extracted signal for the 26th side-band in Fig. 5. Each red-point shows the result from a single TDDFT simulation with the corresponding time delay. In order to extract the time delay, we further fit the numerical signal by an analytic function of the following form: (8) where A, C, σ, t 0 , and τ delay are fitting parameters. Here, τ delay is the time delay, which we aim to extract. In Fig.  6, the fitting function is shown as a blue line. One sees that the fitting function represents the signal very well over the entire time delay range. Even though the absolute time delay with respect to the IR field can be readily extracted from these theoretical simulations, it is in fact hard to extract this absolute delay from experimental results, since the absolute timezero cannot be determined in the experiments. Therefore, so far, experimental results only provide the relative time delay between two different excitation channels such as excitation from different atomic shells. Figure 7 shows the relative time delay for ionization of Ar from the 3s and 3p shells. Red-circles show the TDDFT results, while up and down-pointing triangles show recent experimental results [24,25]. One sees that the TDDFT results are in excellent agreement with the experimental results.
We note that while our work based on the TDDFT with ALDA and ADSIC shows very good agreement with the experiment, recent results by Magrakvelidze et al. [50], also based on TDDFT in the local density approximation, appears to disagree on the same experiment. Here, we discuss a possible origin for this apparent inconsistency and suggest that emerges from the separation of the photoemission delay into two consecutive steps -an approach shared by many RABBITT models. In many works that deal with modelling RABBITT [24,39], the time delay is first decomposed into two components: where τ W is so-called Wigner delay due to the EUV single-photon absorption process [51,52], and τ cc is the so-called continuum-continuum delay due to the additional contribution from the IR field [24,53]. These two delays are often treated and computed separately. In contrast, in the present work, we do not rely on such a decomposition of the RABBITT delay, but directly compute the total delay, by simulating the whole measurement processes, starting from two external laser fields and the system in its groundstate and by performing a time-propagation all the way to the detection of the emitted photo-electron. As a result our method treats the excitation, emission and detection process on the same footing and, as shown in the present work, succeeds to accurately reproduce the experiment. This fact indicates that a fully consistent treatment for all the delay components is significant to correctly understand experimental results, and thus, a direct simulation of the entire measurement processes is required. Furthermore, a separated treatment of the delay components is highly non-trivial for complex systems such as large molecules or solid-state surfaces. In the last case, an additional delay related to the electron transport towards the surface has to be taken into account [28].

C. Many-body effects
One of the strong points of the present TDDFT photoelectron spectroscopy is the capability to investigate the impact of many-body effects directly in the experimental observable. Therefore, the present TDDFT simulation of photoelectron spectroscopy offers novel opportunities to explore the role of many-body effects in the photoelectron emission process. To demonstrate this capability, we here investigate the role of the dynamical electron-electron interaction effects in the argon RABBITT spectroscopy. In our calculation we employ the local density approximation (LDA) for all exchange and correlation effects that are beyond the time-dependent Hartree approximation. While LDA is known to not be able to represent most exchange effects and only weak correlation, the dynamical Hartree potential we use in our calculation should be expected to have a large impact on the dynamics. Since, for the photo-emission process we require a functional with SIC-correction it is not possible to completely separate the effect of the xc-functional and the Hartree potential, but we posit that our results can be considered to be roughly equivalent to at least a random phase approximation level of theory.
To demonstrate the influence of some of the manybody effects as captured by the current approximation, we additionally perform TDDFT RABBITT simulations, where we neglect the time-dependence of the Hartree and the exchange-correlation potentials. That is to say that throughout the propagation of the KS equation, the Kohn-Sham potential is kept "frozen" to the ground state. This treatment corresponds to the independent particle (IP) approximation since all the electrons independently move in a common and fixed mean-field potential. Figure 8 (a) shows the relative delays τ 3s − τ 3p computed by the TDDFT and the IP calculations. One sees that, while the two calculations provide the similar relative delays in the lower photon energy region, they show a discrepancy in the higher energy region. In this high energy range around 42 eV, the Ar 3s photoionization cross section becomes very small due to many-electron effects [54]. This region is the so-called Cooper minimum, and the influence of the Cooper minima in the photoionization delay has been intensively discussed [39]. Previous TDDFT calculations with this ADSIC reported a photoionization cross section in good agreement with the experiment [55].
To obtain further insight of the impact of the manybody effects in the photoionization delay in atoms, we investigate the RABBITT delay for individual Ar 3s and 3p shells. Figure 8 (b) and (c) show the RABBITT delays for Ar 3s and 3p shells, respectively. As seen from Fig. 8 (b), many-body effects play different roles for the photoionization from Ar 3s shell in the lower and higher energy ranges: while the many-body interaction induces a positive delay in the lower energy range, it induces a negative delay in the higher energy range. This fact indicates a correlation among many-body effects, Cooper minima, and the photoionization delay. In contrast, as seen from Fig. 8 (c), the many-body effects uniformly increase the delay of the Ar 3p shell in all the investigated photon energy range. Importantly, one sees that it induces similar amount of positive delay for both Ar 3s and 3p shells in the low photon energy range. Therefore, the influence of many-body effects on the relative 3s-3p delay in the low photon energy range is cancelled out and appears to have no influence on the relative delay.

IV. SUMMARY
In this work, we developed an efficient first-principles attosecond photoelectron spectroscopy technique based on time-dependent density functional theory (TDDFT), focusing on the reconstruction of attosecond beating by interference of two-photon transition (RABBITT). We applied the TDDFT RABBITT simulation to investigate the photoemission from the 3s and 3p shells of Argon. We demonstrated that the TDDFT results nicely reproduce recent experimental results [24,25]. The good agreement of our TDDFT simulation with the experimental results is apparent inconsistency with previous work that also employs TDDFT [50]; Magrakvelidze, et al. reported that the results computed by TDDFT with the local density approximation disagrees with the measured relative Ar 3s-3p time delay. While the previous work computed only the Wigner delay with TDDFT but employed another theory to treat the continuum-continuum delay, the present work treats all the components of the delay at the same level in the TDDFT propagation. Therefore, the apparent inconsistency between the present and previous works may originate from the inconsistent treatment of individual delay contributions of the previous work. This fact indicates the significance of a consistent treatment for each delay contribution and the direct simulation of the whole measurement processes. Furthermore, once target systems become complex, such as large molecules and solid-state surfaces, this kind of step-wise approach to the complete delay becomes nontrivial or unfeasible. Therefore, the fully consistent simulations for the whole measurement processes naturally emerges as a significant tool to attain microscopic insight of such attosecond experiments. Furthermore, the presented TDDFT approach offers novel opportunities to investigate the role of microscopic many-body effects in the photoemission process. In this work we have shown how, by freezing the time-dependent Hartree and exchange-correlation potentials, the role of many-body interactions can be systematically investigated.
As a result, it turned out that many-body effects substantially affect the RABBITT photoionization delay. In particular, we found that the induced delay in Ar 3s photoionization changes its sign around the Cooper minimum. At the moment, accurate description of the exchange-correlation potential as well as electronion coupling is limited, and thus, many-body effects are not fully captured by our TDDFT simulation. However, once a better description for electron-electron and electron-ion interactions is developed, the TDDFT RAB-BITT simulation could be employed to investigate the role of decoherence due to electron-electron, electron-ion and electron-phonon scattering both in the photoemission as well as the transport processes.
While in this work we presented results on a simple system such as gas-phase Argon, the current technique can be readily employed to more complex targets. In particular, the current approach as well as our implementation can be already used to investigate attosecond photoelectron dynamics of solid surfaces. It therefore represents a very powerful and timely technique to guide state-of-theart experiments, and indeed, work along these lines with experiments is already underway.