Physical mechanisms encoded in photoionization yield from IR+XUV setups

We theoretically examine how and to which extent physical processes can be retrieved from two-color pump-probe experiments of atomic and molecular gases driven by an attosecond XUV pulse train and an infrared (IR) pulse. The He atom, the N2 molecule and Na clusters are investigated with time-dependent density functional theory. Results are interpreted on the basis of a simple model system. We consider observables most commonly used in experiments: ionization yield, photo-electron spectra, and angular distributions. We find that the basic time-dependent signatures are dominated by the interplay of IR laser and continuum electrons. System information, contained in the signal, will in general require careful disentangling from the effects of photon-electron dynamics.


Introduction
Photo-induced ultrafast dynamics has become a major field of scientific research, following recent developments in laser technology [1,2]. Times down to the attosecond (1 as = 10 −18 s) scale can now be reached. Attosecond pulses realized in high harmonic generation processes [3,4] open the door to a fully time-resolved analysis of light matter interaction at the electronic level. Such ultrafast irradiation processes are involved in many mechanisms such as vision, photosynthesis or radiation damage of biomolecules [5]. Understanding them is essential, for example to control chemical reactivity down to the electronic, sub-fs scale, but also in astrochemistry or even solar cell developments [6]. Attosecond pulses have already been used to address phenomena such as electron tunneling [7], photoemission delay effect [8] and time-resolved valence electronic motion [9].
A standard two-color attosecond experimental setup consists of a non-ionizing (pump) infrared (IR) laser inducing electronic dynamics followed (with controllable delay) by an ionizing attosecond (probe) extreme ultraviolet (XUV) single pulse [10,11] or pulse train [12,13]. Attopulse trains consist of XUV pulses separated by half the IR period (τ I /2) [4]. The commonly recorded signals are total ionization and Photo-Electron Spectra (PES) [14][15][16][17] as a function of delay time. Until now, a significant number of two-color experiments have focused on the ionization yield of atoms [12,[18][19][20] and small molecules [21][22][23][24][25], as a function of delay. These experiments are driven by the goal to extract system-relevant information from the measured data. Experimental and theoretical studies performed in these simple systems [12,18,24] show a regular modulation of the ionization signal by half an IR period as the dominant response. The modulation coincides with the temporal separation of attopulses inside the XUV train and is rather robust for a wide range of systems. The main finding of our theoretical analysis is that the modulation is dominantly determined by the IR laser driven dynamics of the continuum electron. As a result, it needs to be removed in order to uncover system related properties (polarization dynamics, excitation energies) [26]. Thus, our work presents a starting point towards extraction of system relevant dynamics from the above measurements.
Our theoretical approach is two-pronged. We develop a simple model of two-color attosecond pump-probe experiments based on an analytical and numerical solution of the time-dependent Schrödinger equation (TDSE) [27,28] in the Single-Active-Electron (SAE) approximation [29]. The results are corroborated by Time-Dependent Density Functional Theory (TDDFT) [30][31][32] analysis of complex systems with many interacting electrons driven by two color fields [16,31,33].
We will first consider a typical scenario (IR (I) pulse + XUV (X) attosecond train) in two realistic examples (He atom, Na + 9 cluster), and then employ a simple (mostly analytical) model to work out dominant mechanisms in the measured ionization signal as a function of delay. On the way to figure out system specific information in signals from IR+XUV analysis, we consider the N 2 dimer for which measurements had been done [24] and Na + 21 as a more complex system, in both cases also inspecting more detailed observables such as average kinetic energy from PES and anisotropy in Photo-electron Angular Distribution (PAD). Finally, we come back to Na + 9 with replacing the IR pulse by a frequency close to the Mie plasmon resonance in order to force a signal from the system's response.

Analysis of IR+XUV photo-ionization
We consider a typical experimental setup [13,18,24] composed of an intense IR pulse and an XUV pulse train. Both pulses have wavelengths much larger than the system sizes considered here. Thus we take the external electromagnetic field in dipole approximation as (we use atomic units = m = e = 4π 0 = 1) where τ I = 2π/ω I and f a cosine-squared envelope. Here E I , E X T I , T X , ω I , ω X , are peak electric field, pulse duration, and frequency of the IR and XUV pulses, respectively. n X represents the number of pulses of the XUV train, and g m and t m denote weight and center of each individual train pulse. The IR and XUV pulses are both assumed linearly polarized along the z-direction. A crucial parameter is t d , the delay time taken here between the maxima of the XUV train and the IR pulse.
In such an experimental setup, a standard observable is the total ionization Y, which emerges after applying the pulse. It is, of course, a function of all laser parameters, but the most important one is usually t d . The typical outcome for Y(t d ) is an oscillatory function overlaid with some global trend. Its information content is most efficiently extracted by taking its Fourier transform and calculating a power spectrum P, or total power-of-yield (PoY): where Y lin (t d ) subtracts the global linear trend of the signal to reduce background in the Fourier spectrum. We will mostly focus in the first examples on Y(t d ) and P(ω). Later, we extend the analysis to observables from photoemission, the photo-emission spectrum PES σ PES (E kin ) and photo-emission angular distribution (PAD) σ PAD (ϑ). For their definition and computation see, e.g., [33]. However, PES and PAD are now considered additionally as functions of delay which become quickly an unoverseeable landscape. Therefore, we characterize each one by a single number. The PES is compressed to average kinetic energy E kin : and the PAD to the anisotropy parameter β [34,35]: Both observables are function of delay time t d and the associated power spectra is evaluated in the same manner as for ionization yield, see equation (2).

Two realistic examples
We begin with two realistic examples, (i) the He atom as a simple noble gas system and (ii) the small metal cluster Na + 9 with its strong dynamical polarization. It was observed in reference [36] that in metal clusters the ionization signal might be related to their dynamical polarizability, following some experimental results [24]. We use here TDDFT at the level of time-dependent local-density approximation [37] augmented by a self-interaction correction [38]; see reference [16] for details. Figure 1 shows, for He atom (right column) and Na + 9 (left column), the ionization yields Y(t d ) (bottom row) and the associated PoY (top row). To check the impact of attotrain pulse duration we compare in each case a long train (n X = 22) with one single XUV pulse (n X = 1). While signals in time domain (panels c and d) are not easy to interpret, the Fourier analysis (panels a and b) provides a clearer insight. The predominant features are peaks in the PoY spectrum only at integer multiples of the IR frequency ω I . The zero frequency component reflects the overall trend of the signal (IR pulse envelop) in time domain and is irrelevant here. In most cases, the spectra display peaks only at even multiples of ω I ; but, odd multiples are also observed occasionally for n X = 1. Appearance or absence of such odd multiples may be an indicator for system properties. The relative strength of the P(ω) peaks might also contain quantitative information on the system. The next question then is how to extract information on the system from these signals? This requires to identify how much from Y(t d ) can be attributed to the laser and how much to the system itself. We thus consider a simple analytical model delivering estimates of Y(t d ) and allowing to add or to remove system properties to test their impact.  n X =1 n X =22 Fig. 1. Analysis of ionization following IR+XUV pulses (1a)-(1d) in the He atom and in Na + 9 , and for two different lengths of the XUV train (nX = 1 and nX = 22). For Na + 9 : ωI = 0.4 eV, ωX = 6.8 eV, II = 10 11 W/cm 2 , IX = 10 11 W/cm 2 (nX = 22) and IX = 10 14 W/cm 2 (nX = 1). For He: ωI = 1.4 eV, ωX = 27.2 eV, II = 10 12 W/cm 2 , IX = 10 13 W/cm 2 (nX = 22) and IX = 10 16 W/cm 2 (nX = 1). The IR pulse length covers 30 IR cycles and the XUV pulse 1/10 of an IR cycle. Lower panels: ionization yield Y as a function of delay time t d of the XUV pulse. Upper panels: total Power-of-Yield P as a function of frequency (IR frequency ωI unit).

A simple model
Our schematic model relies on the Single Active Electron (SAE) approximation and utilizes the Strong Field (SF) approximation, in which the dynamics of the (free) continuum electrons is dominated by the laser field and the Coulomb potential can be neglected. The unperturbed system consists of 2 bound states |j and the plane waves of the electron continuum |p with HamiltonianĤ 0 . The dynamical electron wave function is expanded as Its time evolution is determined by the time-dependent Schrödinger equation with the coupling to the photon fieldĤ ext as given in equations (1a)-(1d) and the system's HamiltonianĤ 0 defined byĤ Starting point is |Ψ(t = −∞) = |1 and we are interested on b (p, t = ∞) from which we read off the ionization probability. The time-dependent Schrödinger equation is mapped into equations-of-motion for the expansion coefficients b (p) and c j as: where Ω is the Rabi frequency and the phases of the bound states |j were chosen such that the dipole matrix element 1|r|2 = 2|r|1 is real. Thereby we have assumed that feedback from the continuum electrons to the bound states is negligible (never come back approximation). Furthermore, following experimental conditions, we have assumed that the coupling between the bound states is dominated by the IR laser field and that their population change through XUV photoionization can be neglected. The initial condition reads in terms of the expansion coefficients An approximate solution of equations (6a) and (6b) can be found in the adiabatic limit of slow IR field, i.e., ω I ∆ε. To further simplify the analytical derivation, we assume weak coupling Ω ∆ε. A closed solution is possible without that assumption. But the emerging expression becomes too bulky to be illuminating and the simplified case contains already all features we need. First, we resolve equation (6b) by time integration where we have exploited ∂ t (Ωc 1 )/(∆ε) ≈ 0 in accordance with the assumption of slow IR field and weak coupling. This c 2 is inserted into equation (6a) which so produces an oscillator equation with frequency η∆ε and the result As argued above, the adiabatic approximation ω I ∆ε is fulfilled in the examples of Figure 1. In Na + 9 the HOMO-LUMO gap lies around 1.3 eV and the plasmon at around 2.7 eV far above the chosen IR frequency at 0.4 eV. In He, our ionization potential is 25 eV and the dominant dipole transition around 20 eV, again both well above the IR frequency. Setting the coupling strength d 21 = 0 delivers a minimal model with one single bound state. Except the ionization potential, it contains no information on the system and solely represents the dynamics of a free electron in the external field. By including d 21 , we can then examine the impact of the system's internal structure on the observed outcome.
Inserting amplitudes (8) To further simplify this expression, we exploit the above assumption that direct IR-driven transitions to the continuum are negligible (i.e., F I · D ≈ 0) and that the XUV pulse duration T X is short against the IR cycle. Furthermore, we consider a single XUV pulse instead of a train (n X = 1) which leaves wherefrom only the e iωX(t −t d ) contributes in the integral (9a) while the part with e −iωX(t −t d ) cancels due to the very fast oscillations which it produces together with the e iε1t factor (rotating wave approximation). The pre-factor confines the dt integration to a small region around t = t d due to smallness of T X . We extract a trivial phase factor such that the whole integrand in equation where we have already exploited that D(p, t ) varies very slowly and can be taken effectively at the time t d where the XUV pulse is active. The first line represents a trivial phase factor which becomes irrelevant in the probability |b (p, ∞)| 2 . The integrand in the dt integral is slowly varying within time T X . We set η(t ) ≈ η(t d ) and expand (p(t ) + A(t )) 2 up to first order in t − t d . This leaves us with an analytically solvable integral, leading to a total yield where the integration variable was shifted from p to P = p + A(t d ) and where all dependencies on delay time t d are explicit. It is interesting to note that |b (p, ∞)| 2 in the above expression is nothing else than a PES. The schematic model thus also provides information on PES and its average kinetic energy. One could formally accomplish the P integration. But we can discuss the outcome already at the present stage which avoids the bulky analytical expression. We thus disentangle the time dependencies directly in the expression (10). The key entry in E is and F X is independent of t d as observed previously. Similarly, we have A(t d ) = A I (t d ) + A X in the transition element D. The XUV contribution in Y is weak and thus can be expanded to first order in F X and A X . This yields where E I is defined as in equation (10) with the replacement . Backtracking the various t ddependent terms, we see that they all depend on powers of cos (ω I t d ). Thus, Y(t d ) is a function of cos (ω I t d ).
Consequently, we find the first qualitative result that the PoY P(ω) has peaks only at multiples of ω I , in accordance with Figure 1. Upon closer inspection, we realize that the majority of time-dependent entries (η 2 , E, and ηA) depend on cos 2 (ω I t d ). These terms thus contribute only to the peaks at even integer multiples of ω I in P(ω); they are related to laser field only and do not contain any system information. The terms linear in cos (ω I t) induce a qualitative change in the spectrum, as they contribute to peaks at odd integer multiples of ω I . Three sources produce linear terms in equation (11). First, the terms linear in η are associated with the coupling strength d 12 and thus, are sensitive to the system properties. Second, the term ∝ A · ∇ p d 1 is sensitive to the p-dependence of the transition element. Finally, we have the term ∝ F I (t d ) which already exists at the most basic level, without coupling to intrinsic excitations d 12 and has no reference to p-dependence of the d. There is thus no qualitative feature exclusively reserved to system properties. The latter are imprinted into the PoY, however, in subtle, quantitative manner regulating the relative height of the peaks at the prescribed frequencies.
We have also solved equations (6a)-(6f) numerically. Typical results are shown in Figure 2. The single pulse case (n X = 1, right panel) shows peaks at every integer multiple of ω I for the case with only one bound state, in line with the former analytic prognosis. But the impact of odd multiples of ω I shrinks when stepping up to the 2-level case (coupling strength d 12 = 0.5). The two terms in equation (11) Fig. 2. Total Power-of-yield P(ω) from the simple model with an ionization potential −ε1 = 1 also taken as unit for the energy scale, energy difference ∆ε = 0.2, frequencies ωI = 0.08, ωX = 1, field strengths EI = 0.02, EX = 0.005 and IR pulse width TX = 30τI, under varying conditions (1-level versus 2-level model, two nX) as indicated.
In the model, we can discriminate the two contributions by varying systematically IR strength E I and coupling strength d 12 . The latter is not possible in a realistic system where system information can only be extracted by varying E I . The step to an attosecond XUV train with n X = 22 (left panel) wipes out all signals at odd multiples of ω I . This is understandable as the train's repetition rate of half an IR cycle enhances the weight of even multiples of ω I . But this inhibits any simple (frequency) analysis of the signal in the pulse train limit. Only the relative weights of P(ω) peaks might contain specific system's information. To disentangle information on system properties from the trivial overhead of laser-electron dynamics this requires to perform detailed simulations of the process.

Variation of systems and conditions
The two introductory examples and the model calculations discussed above have shown that the laser-electron dynamics as such dominates the emerging structures of yield in time domain and power spectrum. To find possible signals from properties of the irradiated system, we now consider two other systems, the N 2 dimer with covalent binding and the Na + 21 cluster as a more complex molecule, and we drive the example Na + 9 into the resonant regime hoping to see signals from the Mie plasmon in the IR+XUV analysis. Moreover, we extend the analysis also to the more subtle observables from PES and PAD, namely average kinetic energy E kin and anisotropy β.

The N 2 dimer
So far, there exist only few experimental explorations of attosecond dynamics in complex molecules [25]. We take up here the example of N 2 which, together with a few other molecules of comparable size, was studied in [24] and for which we had performed TDDFT computations leading to good agreement with experimental findings. We have now extended these calculations over longer series of time delays and also analyzed more observables.  [24] with TDDFT calculations. The grey band corresponds to the frequency window used experimentally to recover ionization oscillations.  Figure 3 shows the experimental total PoY together with the theoretical result. Both agree in producing one isolated, strong peak at twice IR frequency which is, in fact, the frequency corresponding to the repetition rate of XUV pulses in the train (two pulses per IR cycle). This is the same result as found above for the He atom, the Na + 9 cluster, and the schematic model. This is thus a generic result dictated by the IR+XUV setup. Unfortunately, experimental delays displayed in [24] do not allow to extract the strength of higher harmonics. Theoretical TDDFT calculations allow to explore this and to include also other observables into the analysis. Figure 4 shows signal in time domain (left panel) and power spectrum thereof (right panel) for N 2 and the three observables ionization yield Y, average kinetic energy E kin , and anisotropy β. The window of delay times is taken in the second half of the pulse and somewhat after, similar as in the previous test cases, see Figure 1. And indeed, the ionization signal resembles these previous test cases. The three different observables show different global trends within the given time window. But the detailed pattern regularly repeated within intervals of one IR cycle are very much the same. This leads to the result that the PoY and power spectra (right panel) show the same pattern for all three observables. The peaks appear only at even multiples of the IR frequency and their strength decreases from one peak to the next. The amount of decrease and the signal to background ratio differs for the three observables. Interesting is, however, that all three observables agree again in predicting a particularly well developed peak at 8 ω I . The latter is not the case for other systems. We are seeing here probably a feature specific for N 2 , thus a glimpse of a system's property.
Note again that laser parameters used here in the N 2 case again match conditions of derivation of the analytical expressions, as were the cases of Na + 9 and He in Figure 1. Indeed, our IR frequency is ω I = 0.4 eV far below the lowest electronic excitation at about 6 eV and the IP at 15.6 eV.

The Na + 21 cluster
As a next example in the variation of systems, we step forward to Na + 21 as a larger system with high spectral density from which we know that it delivers an involved ionization signal [36]. To improve the signal over noise ratio, we consider, again, an XUV train rather than a single pulse. Figure 5 shows the three observables (yield, average kinetic energy, anisotropy) in time domain (left panel) and in frequency domain as power spectrum (right panel). The three observables deliver, up to details, again very similar information. Thus we see once more that more detailed information as delivered by the PES or PAD do not change the frequency content which remains heavily dominated by the laser frequency and its harmonics. Comparing the time signal (left panel) with previous cases, we see more structures than for He (Fig. 1) or N 2 (Fig. 4). This, however, is not necessarily a sign of complexity because the  Fig. 6. Total PoY from IR+XUV pulses (1a)-(1d) in Na + 9 with an "IR" pulse close to plasmon resonance and for two different lengths of the XUV train (nX = 1 and nX = 22). The pulse parameters were: ωI = 2.72 eV, ωX = 9.5 eV, II = 10 9 W/cm 2 , and IX = 10 11 W/cm 2 . The IR pulse length covers 30 IR cycles and the XUV pulse length is 1/10 of an IR cycle. The Fourier spectra are evaluated in two different windows of delay time t d . The left panel shows results from a window t d ∈ [165 fs , 210 fs] during the IR pulse and the right panel from a window t d ∈ [250 fs , 330 fs] after the IR pulse.
smaller Na + 9 cluster shows similarly structured pattern, see Figure 1.
The Fourier signals have, again, peaks only at even multiples of the IR frequency. The sequence of peak heights differs from previous cases. The second peak at 4 ω I is comparatively high which reflects the pronounced substructures in the time signals. The pattern of power spectra are much the same for all three observables although the global trends in the time signal looked rather different. This is achieved by subtracting a linear global trend in the Fourier transformation (2) which helps to concentrate on the fast oscillations of the signal.

Excitation near plasmon resonance
So far, we have considered IR pulses far off any system resonance and so it may be not surprising that we see always a dominance of the IR frequency, or multiples thereof, in the signals. This holds even for metal clusters in which the Mie surface plasmon is a prevailing feature in the dynamical response [16]. To see some impact of the plasmon resonance, we have to move the IR frequency closer to resonance. To this end, we consider now ω I = 2.72 eV which is very close to, but still different from, the plasmon frequency for our test case Na + 9 at ω plas = 2.77 eV. This is, in fact, not a true IR pulse any more. But we continue to use the notion IR to distinguish it from the XUV probe pulses. Note that we need ω I that close to ω plas to obtain sufficient plasmon excitation with the rather long (thus highly selective) IR pulse.
The resulting total PoY are shown in Figure 6. The near-resonant excitation leaves a strong dipole oscillation after the IR pulse terminated which, in turn, allows to see significant Y(t d ) amplitude for XUV delay times after the pulse, much unlike the off-resonant case where the signal fades away quickly with the IR pulse going down, see lower panels in Figure 1. We thus look at two qualitatively different windows of delay time t d , one taken during the IR pulse and another one after. The resulting total PoY are shown in the left and right panels respectively. They should be compared with the off-resonant excitation of Na + 9 in the left panel of Figure 1. There, all peaks in total PoY are centered at multiples of the IR ω I , only even integer multiples for n X = 22 and both even and odd for n X = 1. This changes in the resonant case shown in Figure 6. Already total PoY's collected during the IR pulse develop a broader shoulder at the right side, for higher harmonics even a slight upshift. This stems already from the plasmon resonance contributing more and more to the dipole response. The PoY's taken after the IR pulse show this even more pronounced. All peak positions now show an upshift, which corresponds exactly to the slightly higher Mie plasmon frequency ω plas , and multiples thereof. The effect is small but visible. The major difference (not a small one) between on and off resonance cases is the fact that the oscillatory signal does not fade out once the IR pulse is over. This is a standard behavior for typical fs pulses with frequencies approaching a resonance. It should be noted that we also checked the properties of the power spectra extracted from average kinetic energies and anisotropies. They do exhibit the same features as the total PoY and we thus do not show them. This is so far the only direct system property observed in our IR+XUV analysis. Although the plasmon frequency is already known from other measurements, the observed frequency shift towards the plasmon frequency also contains information on the population dynamics of the collective plasmon oscillation. This effect will be subject to future studies.

Summary and conclusions
We have discussed time-resolved analysis of atomic/molecular systems (He atom, N 2 molecule, Na clusters) at the attoscale by combined IR+XUV pulses, with tunable time separation. To analyze general effects of photonelectron dynamics from system properties, we developed a schematic model with one or two bound states, electron continuum, and photon field. We looked mainly at the ionization yield as a function of delay time and its Fourier analysis (total Power of Yield, total PoY). Additionally, we also explored two other observables deduced from the distributions of emitted electrons, their average kinetic energy and the anisotropy of their angular distribution.
We have shown that the universal modulation by the IR frequency strongly dominates the signals (ionization, kinetic energy, anisotropy). This raises the key question on how to access specific system's properties through such attosecond setups. All results (schematic model with or without internal structure, realistic systems) produce rather similar patterns, namely regular oscillations about a global trend, sometimes cosine-like, sometimes with extra beats. In all cases, the oscillations are aligned with the IR cycle which becomes apparent in frequency space in that the ionization signals contain only integer multiples of the IR frequency, in most cases even only at even multiples. All scenarios for signal structure (regular oscillations, beating, . . . ) and corresponding power spectra can be produced already in the simplest model consisting only of one electron in the photon field. We could not identify a direct, ideally qualitative, signal of system properties and time scales. The impact of the system is surely contained in the signals, e.g., in the relative relations of peak heights in the power spectra. However, it requires tedious, and model-dependent, analysis to extract the wanted information from the kinematic background of electron-photon dynamics. The situation is the same as often encountered in involved reactions namely that only model-dependent analysis paves the path to the underlying system. There is one interesting exception where we seem to see a system property directly. That happens if we drive the IR pulse close to a system's resonance which delivers for delay times toward the end of the IR pulse, and more so after, a possible indication of the resonance energy. One can even trace back how the system frequency gains in weight during the process which then is time resolved information at least at fs time scale level.
The results of our survey call for further developments in two directions. At the experimental side, one needs more detailed setups and observables to eliminate kinematic background, e.g., by coincidence [39] experiments. Such experiments are already underway at several places. At the theoretical side, one needs to work on reliable modeling of the reaction process to look closer into the system. More elaborate measurements become particularly challenging for theory as many-body correlations may then also play a role [40][41][42] and need further disentangling. Work along these directions is in progress.