Electron population dynamics in resonant non-linear x-ray absorption in nickel at a free-electron laser

The modern understanding of complex materials relies on suitable approximations to the unabridged quantum mechanical description of the full, correlated many-body problem. To assess the predictive power of theoretical models and the selected approximations, detailed experimental studies of systems driven out of equilibrium are particularly insightful.

X-ray free-electron lasers (FELs) with their high-power (about 10¹⁰ photons/pulse after monochromatization) and ultrashort pulses (tens of fs) are uniquely suited for probing the non-equilibrium dynamics of matter in all states from solid state to plasma.^1,2

While schemes with separate pump and probe pulses can be used to observe dynamics slower than the pulse durations, transmission measurements of individual intense pulses can be used to infer the x-ray interaction and material dynamics within the pulse duration. With increasing pulse energy, the interaction with an exposed solid becomes highly non-linear, as the material state changes during the interaction. In non-linear absorption measurements, a single intense x-ray pulse acts as both a pump and a probe at the same time, and its total absorption depends on the changes it induces in the electronic system.^3–9 The highly excited solid evolves toward a state of warm dense matter (WDM) where individual electronic excitations reach up to hundreds of eV and excitation energies average out to many eV per atom, while the nuclear lattice still resembles that of the cold solid during the femtosecond pulse.^10–19 

In this paper, we present new fluence-dependent x-ray absorption spectra recorded with monochromatic x rays on metallic nickel thin films around the nickel 2  $p 3 / 2$ (⁠ $L 3$⁠) edge, revealing changes in the valence electron system around the Fermi level, driven by FEL excitation densities spanning several orders of magnitude from the linear regime up to 60 J/cm² (corresponding to 2  $× 10 15$ W/cm²).

The electronic processes that ensue after the absorption of core-resonant photons trigger a complex series of dynamical processes of photon-absorption, electron excitation, and subsequent electronic scattering. Such non-equilibrium processes with a large amount of correlated particles involved and spanning orders of magnitude in internal energy are challenging to treat in ab initio simulations.^20–24 Rather than striving for a direct time-dependent solution to the quantum-mechanical many-body problem, our analysis, therefore, explores how much of the observed effects can be explained purely by the evolution of electronic populations within the ground-state density of states (DOS) and scaling the known ground-state probabilities of absorption, decay, and scattering processes with the current populations of the participating states. Therefore, we develop a rate equation model to provide an intuitive understanding of the electronic processes. Our model successfully describes the largest part of the non-linear changes in the spectra, affirming the dominant role of electron redistribution. However, especially in close vicinity of the resonance, our measurements deviate from the predictions of the rate model and indicate the need for more sophisticated theories. Nevertheless, our straightforward picture of an intense core-resonant x-ray pulse interaction with the valence system of a 3d metal lays a solid knowledge-based foundation for the planning and interpretation of non-linear x-ray spectroscopy experiments at FELs.

This paper is structured as follows: Sec. II describes the experimental setup used to acquire non-linear XANES spectra of the L₃-edge of metallic nickel films. In Sec. III, we introduce the working principle of the rate model used to interpret the experimental data. In Sec. IV, we show the experimental data in direct comparison to the simulation results and continue to interpret our findings in Sec. V. We conclude in Sec. VI.

X-ray absorption spectra of the nickel 2 $p 3 / 2$ (⁠ $L 3$⁠) edge were recorded at the Spectroscopy and Coherent Scattering Instrument (SCS) of the European XFEL.²⁵ 

The XAS (x-ray absorption spectroscopy) spectra were measured by continuously scanning the SASE3 monochromator²⁶ (synchronized with the undulator gap) back and forth many times in the range 846–856 eV. The photon bandwidth was about 420 meV full-width half-maximum (FWHM) and the FEL pulse duration on the sample was about 30 fs FWHM. The polarization was linear horizontal. The overall beam intensity was controlled using a gas attenuator filled with nitrogen and monitored using an x-ray gas-monitor (XGM) downstream of the monochromator.^27,28 Figure 1 illustrates the experimental concept. Panels (a) and (b) display the DOS (integrated over spin-up and spin-down states) next to the resulting XAS spectrum (including non-resonant background); panel (a) shows the situation of a low-fluence measurement, panel (b) shows the situation where sufficient photons are used to alter the DOS occupation and hence the measured XAS spectrum. Panel (c) shows the geometry of zone-plate, sample, and detector.

For x-ray absorption measurements at FELs based on self-amplified spontaneous emission (SASE), beam-splitting schemes can deliver optimal normalization of SASE-fluctuations.^30–32 Here, we utilize such a scheme using a focusing and beam-splitting zone plate which also creates the required tight focusing to achieve extreme fluences. The zone plate combines an off-axis Fresnel structure for focusing and a line grating for beam-splitting in a single optical element.²⁹ It thus produces two μm-sized, identical foci in the sample plane, 1.9 mm apart, originating from the first-order diffraction of the zone plate, as well as the positive and negative first orders of the line grating. The sample has a square support of 25 mm size, containing square Si₃N₄ membrane windows (orange in Fig. 1) of (0.5 × 0.5) mm² size and 200 nm thickness with a distance of 2 mm between adjacent windows. Every second pair of rows (blue in Fig. 1) was additionally coated with a 20  $nm$ thick sample layer of polycrystalline metallic Ni by sputter deposition, on top of a 2  $nm$ bonding layer of Ta; a 2  $nm$ Pt capping layer was applied to prevent oxidation during sample-handling.

The sample frame was positioned such that one zone plate focus impinged on a nickel-coated membrane, while the other hit a bare silicon-nitride membrane. Thus, the difference in transmission of both beams can be dominantly attributed to absorption in the nickel film.

The detector was a fast readout-speed charge-coupled device (FastCCD) with a high dynamic range, enabling 10 Hz readout and increasing the fluence range available to the experiment.^33–35 Due to an unstable detector temperature, significant retroactive calibration of the detector was necessary (see Subsection 2 of  Appendix B). To prevent detector saturation during measurements with the unattenuated beam, an additional aluminum filter of about 13  $μ m$ thickness was used between sample and detector.

During these high-intensity measurements, sample and reference films were locally damaged by intense individual FEL shots. Thus, the FEL was operated in single-shot mode at 10 Hz repetition rate, and the sample was scanned through the beam continuously at 0.5  $mm s − 1$⁠, resulting in ten shots per membrane window.

The shot craters in the reference membranes were later analyzed with scanning electron microscopy (SEM) to determine the effective focal size at specific photon energies. The resulting spot sizes were used to calibrate ray-tracing calculations which delivered the photon-energy-dependent spot size, ranging from 0.4  $μ m 2$ to about 3  $μ m 2$ (see Subsection 4 of  Appendix B for details on the spot size determination).

Various approaches have been proposed to describe the interplay between photon absorption and the electronic structure evolution during the absorption of an FEL pulse. Ab initio methods such as Monte Carlo calculations, which explicitly calculate a large number of individual particles' interaction pathways,³⁶ and time-dependent density functional theory (TD-DFT), which sets out to solve the full quantum-mechanical many-body problem in terms of the electron density,³⁷ generally scale poorly with particle number. In contrast, rate models provide a simpler yet useful tool by describing the interplay between photon absorption and electronic system using non-quantized volume-average quantities and rates directly on a macroscopic scale.^5,38–41

Away from material resonances, rate models have been successfully used to describe fluence-dependent x-ray absorption in three-level systems, representing the ground, core-excited, and intermediate valence-excited states.^5,38 When probing the valence bands around material resonances, however, the evolution of the electronic system requires explicit modeling of the energy-resolved valence state populations.⁴² Tracking the full non-thermal population history proved crucial for accurately describing the non-linear absorption changes near and around the Fermi level. We assume that on the modeled femtosecond timescale, the DOS does not change significantly, which is motivated by the slower lattice reaction. This approach though cannot capture subtle changes in electron correlations.⁴³ Furthermore, since our measurements are done with linear polarization and are therefore not sensitive to the magnetization, we do not consider the exchange-split DOS but rather integrate over minority and majority electrons.

We describe the propagation of x-ray photons through the sample as well as the dynamics of electron populations within the sample using a set of ordinary differential equations. These are assembled from terms that each describe the rate of a specific physical process. The rate of each process is based on a tabulated or measured ground-state parameter, such as the Auger-lifetime or the absorption cross section, scaled with the appropriate fractional occupations at the simulated time.⁴⁴ The relevant process rates are compiled into differentials of electronic populations and photon density in space and time and implemented in a finite-element simulation to derive the electron population history and ultimately the x-ray transmission of a three-dimensional sample. Only the time constants for the valence band thermalization and the scattering cascades of free electrons are treated as free parameters and adapted to fit the experimental data.

The model considers an idealized three-dimensional sample traversed by an x-ray pulse with Gaussian shape in space and time. We make key approximations to reduce computational effort, such as neglecting any movement of electrons within the sample [consider the inelastic mean free path of relevant photoelectrons of about 1.3 nm (Ref. 45)] and describing photon propagation exclusively in the forward direction. The temporal evolution is solved using the fourth-order Runge–Kutta method with adaptive time-stepping. The propagation of photons in space is calculated as if it happened instantaneously in between the time-steps using the explicit fourth-order Runge–Kutta method.

In order to account for the two-dimensional Gaussian transversal intensity profile of the FEL spot, we first calculate the transmission of the sample for transversally uniform illumination for different fluences. Since we omit transversal coupling, the response to the Gaussian beam profile can then be reconstructed by appropriate radial integration over many values obtained for constant illumination. With these simplifications, the overall computational complexity is drastically reduced, as we simplify a problem with partial differentials in four dimensions into two separable one-dimensional initial value problems, one for photon propagation in space and one for the evolution of electronic populations in time.

The model describes the interaction between three types of population densities of electrons as well as incident photons via six distinct physical processes, listed in Tables I and II, respectively. Figure 2 schematically illustrates their relationships. The electron populations R_C and R_V describe the total number of electrons bound in the core and valence system, respectively, for an average single atom in the sample. Their values are limited by the number of available states, M_C and M_V. In the presented nickel L₃-edge spectra, the ground-state populations are R_C = 4, representing the 2 $p 3 / 2$-electrons and R_V = 10, representing electrons from the 3d and 4s states. We describe the electronic population of the valence system in an energy-resolved manner, splitting it up into a discrete number of densities ρ_j, where j represents the index along the valence energy axis. The number of available states m_j for each energy bin in the valence system is derived from the calculated ground-state DOS^46,47 up to 30 eV above the Fermi level E_F. Beyond this value up to 800 eV above the Fermi level, the DOS of a free electron gas is used.²² All electrons with even higher energies, such as photo-electrons created via non-resonant absorption and Auger-electrons from the decay of core-holes, are described in a separate pool of electrons R_free without energy resolution, although the total energy of electrons in this pool is tracked by the parameter E_free.

The processes of resonant absorption and stimulated emission are summarized in the resonant interaction process $P j res$⁠. This is possible because both processes describe a transfer of electrons between the 2 $p 3 / 2$ core level and the valence system following the same interaction cross section, just with inverted sign. The non-resonant absorption process $P i , j non − res$ describes photon absorption at all other electronic levels. The non-resonant absorption length of the ground state is derived from the measured pre-edge absorption. Because of the relatively small contribution of the intermediate 3p and 3s electrons, the model ascribes all non-resonant absorption events to valence electrons, transferring them to the pool of high-energy electrons R_free. Sequential two-photon absorption (TPA) processes are implicitly treated by the model as a resonant absorption event followed by a non-resonant absorption event. This description does not account for the coherence and resonant enhancement of the TPA process. However, the fluences used here do not exceed 1.5  $× 10 31 photons cm − 2 s − 1$⁠, where according to the scaling law proposed by Szlachetko et al.,⁴⁸ the relative contribution of the direct TPA process should be on the order of 1%. Thus, the inaccuracy caused by not treating this process explicitly should be very small. Apart from photoelectrons from the valence system, the Auger–Meitner decay process $P j Auger$ likewise contributes to the free electron pool. These free electrons proceed to scatter with the valence electron system. The total rate of scattering $P scatt$ is modeled by a simple lifetime parameter $τ scatt$⁠, which characterizes how quickly free electrons re-join the valence system. The redistribution process $P j red$ distributes the kinetic energy of each scattering free electron among the valence system. Avoiding the complexity of explicitly calculating the multiple collisions involved in these electron scattering cascades, the algorithm instead approximates the effect of such a cascade: The energy of the scattering electron is spent to elevate an evenly distributed fraction of all valence electrons from occupied states to locally available unoccupied states (as indicated in purple in Fig. 2). Note that this process is not independent but represents an immediate consequence of the free electron scattering process; the scattering time $τ scatt$⁠, thus, characterizes both $P scatt$ and $P j red$ together. Finally, electronic thermalization is modeled with a bulk timescale $τ th$ (essentially quantifying electron–electron scattering) that moves the non-thermal valence electron distribution toward a target Fermi–Dirac distribution that corresponds to the momentary internal energy and population of the valence system.

The full mathematical description of the process terms is given in Subsection 1 of  Appendix A, while the choice and derivation of input parameters is detailed in Subsection 3 of  Appendix A.

From these process terms we assemble the time-differentials of the populations of electrons and photons listed in Table I. The movement of electrons between states is represented by process terms of electronic transitions appearing symmetrically in these differentials with positive and negative signs, thus ensuring the conservation of particle number. For example, the term for Auger decay appears twice with a negative sign in the valence electron differential, and once each with a positive sign in the differential for core- and free electrons.

The number of photons is reduced or increased by resonant interaction and reduced by non-resonant absorption. The model allows for an arbitrary number of incident photon energies E_i, each of which must be resonant to a specific bin of the valence energy system E_j. The temporal intensity profile of all incident photons is modeled as Gaussian. For the presented calculations, only a single resonant photon energy was used, representing measurements with monochromatic x rays. Incident photons are the only source of energy flow into the system, and all energy eventually emerges as the thermal energy of the valence system.

The valence system interacts via all modeled processes. The resonant absorption rate $P j res$ changes the valence electron densities ρ_j at all incident photon energies E_i that are resonant. Via non-resonant absorption, incident photons can be absorbed by electrons from all valence states ρ_j. Auger processes depopulate the valence system. The thermalization drives electrons toward the Fermi–Dirac distribution based on the current internal energy and population of the valence band, without changing the total valence occupation. Electron scattering P^scatt causes electrons from the free electron pool to re-join the valence system in a random unoccupied state h_j, as well as redistributes electrons inside the valence system in an electron cascade triggered by the process, which is described via the $P j red$ term.

The population of core electrons is reduced (or increased, depending on the sign of $P i , j res$⁠) by resonant transitions of all incident photon energies E_i to states at all resonant energies ρ_j, and is increased by Auger decay from electrons of all energies j in the valence system. Note that the radiative emission channel is neglected in our model as it is designed for soft x-ray energies where Auger emission accounts for most core-hole decays (here specifically, 99.1% of the nickel L₃ core-hole decays^49,50). In another concession to the specific experiment simulated here, we further neglect fast electrons leaving the sample, since the electron mean free path of about 1.3 nm (Ref. 45) is much shorter than the sample thickness of 20 nm. While the model is generally suited to implement a loss process for free electrons, the total number of electrons in the system being strictly constant over time is a valuable indicator for the self-consistency of the calculation.

Free electrons are generated by non-resonant absorption from all incident photon energies E_i as well as Auger-decays from all energies in the valence band. The population is reduced by the free electron scattering rate $P scatt$⁠.

See Subsection 2 of  Appendix A for the differential equations.

Figure 3 shows the measured spectra for the nickel $L 3$-edge next to simulated spectra for increasing x-ray fluence over more than three orders of magnitude, from 0.03 to 60 J/cm². Each measured point represents an average of several FEL shots, sorted by x-ray fluence and photon energy. The varying statistical uncertainty is a result of the pulse intensity fluctuations of monochromatized SASE radiation⁵¹ in combination with photon energy-dependent spot sizes (see Subsection 3 of  Appendix B for details on the shot classification).

We observe four main fluence-dependent effects, which we quantify and compare to the simulated results in Fig. 4: (a) a redshift of the absorption edge of up to 0.9 ± 0.1 eV in the rising flank; (b) an increase in the pre-edge absorbance, as the rising edge of the absorption peak shifts and broadens; (c) a reduced peak absorbance; and (d) and (e), a reduced post-edge absorbance. The integration regions from which the effects (b), (d), and (e) are derived, are highlighted in Fig. 3 as (I), (II), and (III), respectively. Each is 0.6 eV wide. The shift of the absorption edge is quantified by the photon energy at which the absorbance reaches half of the peak value; its uncertainty is propagated from the statistical uncertainty of the absorption peak measurement.

As apparent from Fig. 4, the rise in absorption at the pre-edge (region I), the drop in post-edge absorption (region III), as well as the shift of the rising edge show good agreement within the measurement uncertainties between simulation and experimental data. The deviations observed in the absorption level of the resonance peak and just beyond will be discussed later.

We emphasize that this level of agreement with the experimental data are achieved across more than three orders of magnitude in fluence, based on a rather simplified description of well-known physical processes in combination with experimental or tabulated ground-state properties such as density, electronic configuration, and ground-state spectrum. Only the valence thermalization time $τ th$ and electron scattering time τ_scatt were varied to achieve the best match to the experimental results. We obtain a value of $τ th = 6$ fs, which compares well to recent estimates for excitations on this energy scale.^40,41,52,53 The scattering time constant $τ scatt = 1.5$ fs produces the best agreement with experimental data. This value appears reasonable as it summarizes a cascade of many individual electron scattering events, which we would expect to occur roughly every 100 as.⁴⁵ 

Before further interpreting the non-linear effects shown in Fig. 4, let us first consider the example of a local valence band population history as shown in Fig. 5. The example is drawn from the uppermost 4 Å thick voxel of the simulated sample, excited with a Gaussian pulse profile centered around t = 0 with 30 fs FWHM duration and 30 J/cm² fluence. As such, the example is selected from the upper range of extreme excitations in this simulation to showcase the effects clearly. While panel (a) of Fig. 4 shows the calculated DOS as used by the simulation and published in Refs. 46 and 47, the colormap in (b) shows the occupation of these states over time. It is apparent that the occupation function mostly resembles a Fermi–Dirac distribution evolving from cold to hot. However, the states at E_j = E_i (highlighted by the blue ellipse), show greater population as they are directly populated by the resonant absorption process. We also show the effective electron temperature T and chemical potential μ, which are calculated from the internal energy and population of the valence system at every time step. Panel (c) shows the number of electrons per atom in the valence band below and above the Fermi level (blue solid and dashed curves, respectively) as well as the average number of core holes and the number of free electrons over time. One general observation is that for the given 30 fs pulse duration, the number of simultaneously existent core holes remains very small, even for high fluences. This has two reasons: On the one hand, the natural lifetime of the core-holes of 1.4 fs is small compared to the pulse duration.⁵⁰ On the other hand, the monochromatic excitation near the material resonance implies that the photons couple the core-level to a narrow selection of localized valence states.⁵⁴ In this case, the number of resonant valence states is small in comparison to the number of core electrons. Since the core-level and resonant valence states operate like a two-level system in which absorption and stimulated emission compete, the resonant absorption process saturates due to occupied valence states long before the core level is significantly depleted. This bleaching of valence holes is amplified over the pulse duration by an increasingly heated Fermi–Dirac distribution, which also increases the occupation of states above the Fermi level. Since both core-holes and free electrons decay so quickly, a majority of the absorbed energy is quickly translated into a broadening of the valence electron distribution. By the end of the pulse in this example, more than half of the 3d valence electrons are excited to valence states above the Fermi level, while the highest instantaneous number of core holes was only about one per 100 atoms, as shown in Fig. 5(c).

With these general observations about the evolution of the electronic system within the pulse duration in mind, we can now proceed to interpret the mechanisms responsible for the non-linear features in the spectra. Above the absorption edge, the decrease in absorption with increased fluence [see Fig. 4(d)] can be understood as a depletion of valence states available to the resonant core-to-valence transition. Similarly, below the absorption edge, the increase in absorption [see Fig. 4(b)] can be attributed to valence holes below the Fermi-level becoming available due to the thermalization process, as soon as the tail of the Fermi–Dirac distribution reaches the probed energy. The shift of the absorption edge [see Fig. 4(a)] can be explained by a non-linear combination of the two effects above. Consider that below the absorption edge at the beginning of the pulse, the sample only interacts with the x rays via the comparatively weak process of non-resonant absorption. However, once the sample is sufficiently heated that valence holes become available, additional resonant absorption begins to occur and accelerates further electronic heating—and in turn additional pre-edge absorption. Since the onset of this exponential process occurs earlier near the absorption edge, it contributes significantly to the observed spectral redshift. Another cause of the observed edge shift is the shift of the chemical potential μ, which strongly depends on the exact shape of the DOS and is shown in Fig. 5(b) as a green line. Initially, μ increases with absorbed fluence, as thermally excited electrons from the 3d states must spread out in energy to the lower DOS above the Fermi level. With rising electronic temperature, the high DOS of the 3d states becomes less relevant and the chemical potential drops again as expected in regular metals. A similar evolution of the chemical potential and electronic temperature was predicted for optically excited nickel by previous experiments and calculations.^14,55–57 It is remarkable that the experimentally observed redshift of 0.9 ± 0.1 eV can be reproduced by the rate model based on this very simple mechanism. However, this mechanism applies specifically to non-linear absorption using monochromatic x rays. Qualitatively similar redshifts have been observed in nickel after excitation with optical lasers even at up to three orders of magnitude lower excitation fluence.^43,58,59 These redshifts have recently been linked to modifications of the band structure due to the interplay of electronic correlations and optically induced demagnetization.⁴³ While such subtle, spin dependent effects may also occur in our high-fluence study, they are evidently overshadowed by the electron population dynamics. This aspect of an initially non-thermal electron distribution evolving toward a Fermi–Dirac distribution was also observed as a critical aspect of the optically excited spectra.^43,60 The timescale of electron thermalization was estimated to just over 100 fs, which is about 20 times slower than our estimate of 6 fs. This apparent discrepancy results from a scaling of the thermalization time with excitation density. Such a scaling is supported both by theory⁵² and recent pump–probe studies at the nickel M_2,3 edge⁵³ where the electron thermalization time decreased from 34 to 13 fs with rising optical pump fluence from 8 to 62 mJ/cm². This implies that the value of 6 fs found in our study represents an average time constant for the excitation densities in our experiments.

A significant deviation between model and experiment can be observed at the resonance peak itself, where the simulated electron dynamics lead us to expect a much stronger saturation effect than observed experimentally [Fig. 4(c)]. This saturation is reduced as the peak position redshifts and more thermally vacated d-states become resonant, leading the model to predict a slight recovery of absorbance at the highest fluences. Overestimating the saturation effect may be related to a fluence-dependent decrease in the excited state lifetime due to stimulated emission as well as increased carrier mobility around the Fermi edge, both leading to an energetic broadening of the resonant core–valence transition. Such a broadening would increase the number of resonant valence states and thus delay saturation especially at the edge, but is not considered in our model. While it may be expected that a purely population-based model cannot fully represent resonance effects at the resonance peak itself, the lack of any significant saturation around 852 eV [Fig. 4(d)] is more surprising. Both disagreements point to additional physical effects and call for more sophisticated models.

We speculatively propose three mechanisms which could contribute to these discrepancies: First, the transition matrix elements could get modified at higher excitation densities, especially around the resonance, while we model the absorption only based on the ground-state spectrum. Second, an energy dependence of the electron–electron scattering cross section could allow for particularly fast scattering of electrons with certain energies, counteracting the saturation. Third, a collective, correlated response of the electronic system (local field effects) could modify the DOS or the transitions even on the fast timescale of the FEL pulse duration.^43,61

A more detailed discussion of the model can be found in the  Appendixes: In Subsection 4 of  Appendix A, we show how a variation or elimination of specific processes leads to different predictions for the spectra, and in Subsection 5 of  Appendix A, we discuss the limitations of the rate model and its suitability for future extension.

To summarize, we have analyzed fluence-dependent near-edge x-ray absorption spectra of the nickel 2 $p 3 / 2$ core level up to x-ray fluences of 60 J/cm². We have developed a rate-equation model based on differential equations that describes the excitation and decay processes connecting populations of core and valence electronic states. Process rates are quantified by scaling known ground-state properties with evolving electron populations.

The model enables an understanding of the electronic population history under strong x-ray fluences and characterization of the resulting non-linear absorption near a core resonance. It successfully predicts the observed increase in absorption before and its decrease beyond the resonance, as well as the fluence-dependent redshift of the absorption peak over three orders of magnitude. However, the bleaching of the absorption peak is overestimated by the population-based model and will require more sophisticated models to accurately quantify. Here, the population dynamics rate model also provides a valuable point of reference for more advanced theoretical frameworks.

Providing the fundamental fingerprints of how strong x-ray fluences alter the electronic system and thus the absorption spectra, our straightforward picture of intense core-resonant x-ray pulse interaction can inform the design and interpretation of future FEL experiments. On the one hand, our model can guide the decision up to which point to maximize fluence for good statistics while keeping the absorption process linear, and to recognize the principal spectral fingerprints emerging at the onset of non-linear absorption due to electron dynamics within the pulse. On the other hand, an understanding of the population dynamics within high-fluence pulses, and in particular, an awareness of the dominant influence of electronic scattering processes, is crucial for emerging techniques that aim to utilize x-ray wave-mixing processes, such as stimulated core hole emission, in solids.^40,62–71

We acknowledge European XFEL in Schenefeld, Germany, for provision of x-ray free-electron laser beamtime at the SCS instrument and would like to thank the staff for their assistance. Funding by the Deutsche Forschungsgemeinschaft (DFG, German Research Foundation)—Project-ID 278162697 - SFB 1242 and the Helmholtz Association (Grant No. VH-NG-1105) is gratefully acknowledged. Access to Synchrotron SOLEIL through proposal ID 20160880 for characterization of static properties of the Ni films is acknowledged. Parts of this work were funded by the Swiss National Science Foundation (Grant No. PZ00P2-179944).

The authors have no conflicts to disclose.

Robin Yoël Engel: Data curation (equal); Formal analysis (equal); Investigation (equal); Methodology (equal); Software (equal); Visualization (equal); Writing – original draft (equal); Writing – review & editing (equal). Oliver Alexander: Data curation (equal); Formal analysis (supporting); Investigation (equal). Kaan Atak: Investigation (equal). Uwe Bovensiepen: Writing – review & editing (equal). Jens Buck: Investigation (equal). Robert Carley: Investigation (equal). Michele Cascella: Investigation (equal); Methodology (equal). Valentin Chardonnet: Investigation (equal). Georghe Sorin Chiuzbaian: Investigation (equal). Christian David: Conceptualization (equal); Investigation (equal); Methodology (equal); Resources (equal); Writing – review & editing (supporting). Florian Döring: Conceptualization (equal); Investigation (equal); Methodology (equal); Resources (equal). Andrea Eschenlohr: Investigation (equal); Writing – review & editing (equal). Natalia Gerasimova: Investigation (equal); Methodology (equal); Writing – review & editing (supporting). Frank de Groot: Writing – review & editing (equal). Loïc Le Guyader: Conceptualization (equal); Data curation (equal); Investigation (equal); Methodology (equal); Software (equal); Validation (supporting); Writing – review & editing (equal). Oliver S. Humphries: Data curation (equal); Investigation (equal); Visualization (equal); Writing – review & editing (equal). Manuel Izquierdo: Investigation (equal); Resources (equal); Writing – review & editing (supporting). Emmanuelle Jal: Resources (equal); Writing – original draft (equal). Adam Kubec: Resources (equal). Tim Laarmann: Writing – review & editing (equal). Charles-Henri Lambert: Resources (equal). Jan Lüning: Conceptualization (equal). J. Marangos: Conceptualization (equal); Project administration (equal). Laurent Mercadier: Investigation (equal); Methodology (equal). Giuseppe Mercurio: Investigation (equal). Piter S. Miedema: Investigation (equal); Writing – review & editing (equal). Katharina Ollefs: Writing – review & editing (equal). Bastian Pfau: Writing – review & editing (equal). Benedikt Rösner: Conceptualization (equal); Data curation (equal); Investigation (equal); Methodology (equal); Resources (equal). Kai Rossnagel: Writing – review & editing (equal). Nico Rothenbach: Investigation (equal); Resources (equal). Andreas Scherz: Investigation (equal); Methodology (equal). Justine Schlappa: Investigation (equal). Markus Scholz: Writing – review & editing (equal). Jan Oliver Schunck: Writing – review & editing (equal). Kiana Setoodehnia: Methodology (equal). Christian Stamm: Methodology (equal). Simone Techert: Conceptualization (equal); Investigation (equal). Sam Vinko: Writing – review & editing (equal). Heiko Wende: Methodology (equal); Writing – review & editing (equal). Alexander A. Yaroslavtsev: Investigation (equal); Methodology (equal). Zhong Yin: Methodology (equal); Writing – review & editing (equal). Martin Beye: Conceptualization (equal); Funding acquisition (equal); Investigation (equal); Project administration (equal); Supervision (equal); Validation (equal); Writing – original draft (equal); Writing – review & editing (equal).

Raw data were generated at the European XFEL large scale facility and are available at http://doi.org/10.22003/XFEL.EU-DATA-002170-00. Derived data supporting the findings of this study are available from the corresponding author upon reasonable request.

The simulation code that supports the findings of this study is openly available in GitHub, Ref. 72.

This appendix provides an in-depth discussion of the rate model, including the formalism for each of the modeled processes (Subsection 1 of  Appendix A), differentials (Subsection 2 of  Appendix A) and the choice of input parameters (Subsection 3 of  Appendix A), as well as a study of the models' predictions under different conditions (Subsection 4 of  Appendix A). Finally, we present an extended discussion of the applicability and limitations of our rate model (Subsection 5 of  Appendix A) beyond the presented dataset.

For a better overview, we introduce the processes as individual terms and assemble them into differential equations in Subsection 2 of  Appendix A. In principle, each process is described using an absorption length or lifetime which is known from ground-state measurements and then scaled linearly with the changing electron populations with respect to the ground state. The normalization is such that the ground-state rate is reproduced for an undisturbed electron system and the rate vanishes when the corresponding transition cannot happen due to a lack of electrons or holes. We use the indices i and j to refer to specific energies, where the index i is used for photon energies of x rays and the index j for the energy of electronic states in the valence band.

The first terms (in brackets) represent the difference in the occupation of core states $R C / M C$ and resonant valence states $ρ j / m j$⁠. The dominance of absorption over stimulated emission or vice versa is determined solely by this difference, as they represent an optically driven two-level system in the incoherent limit. If the core level population is smaller than the valence population, the resonant interaction process becomes negative, representing the dominance of stimulated emission. The second term on the right is the number of irradiated photons divided by the penetration length. The Kronecker delta ensures that only photons and electrons in corresponding energy bins interact.

The interaction is normalized by the total valence band population in the ground state $R V 0$⁠, so that the sum of the first term over all j becomes unity if all $ρ j = ρ j 0$ (since $R V 0 ≡ ∑ j ρ j 0$⁠). The second term represents the non-resonant absorption in the ground state as can be experimentally determined sufficiently before the resonance in the spectrum.

This treatment does not explicitly differentiate the non-resonant absorption from core energy levels other than the one treated by R_C. In the given example with photons resonant to the nickel 2p-absorption, the 3s and 3p core electrons only contribute to a minority of the non-resonant absorption events. In this model, we choose for simplicity to scale this contribution together with the non-resonant scattering from the valence electrons.

An exemplary incidence profile and the resulting transmission over time are shown in Fig. 6.

The first factor (in brackets) is the number of unoccupied core states, i.e., core-holes. The second factor describes the relative population of electrons at the energy E_j and the third term is the relative population of the entire valence band to which the electron could transfer its energy. The latter two are normalized by the respective ground-state population. The last term is the decay rate in the ground state, where τ_C represents the ground-state lifetime of a single core-hole. Altogether, this describes Auger decays as interactions between two valence electrons, one emitted, and one filling the 2 $p 3 / 2$ core-hole. In reality, some fraction of Auger decay events will emit electrons from the 3s or 3p core levels instead, followed by further Auger processes which emit electrons with the remaining energy of the original core-hole. These are not treated separately in our description, since the indirect decay is, on the one hand, a minority contribution and on the other hand, ultimately results in the same energy transfer to the valence band, albeit with a slightly longer time delay due to the intermediate steps.

Inspired by earlier approaches to a simplified solution of the Boltzmann equation,⁷³ we approximate the scattering rates of electrons in terms of characteristic time constants τ_scatt and τ_th for the free electrons and valence electrons, respectively.

While the calculation of r_j from T and μ is straightforward, determining T and μ from U and R_V is an inverse problem. This is solved by iterative optimization using the Levenberg–Marquardt method.⁷⁴ 

Note that since, by definition, both electron densities—the momentary ρ_j and the thermalized goal distribution r_j—hold the same amount of electrons and internal energy, the sum over the valence band of all electron-distributing thermalization rates is always zero and the change in overall energy is also zero, i.e., $∑ j P j therm = 0$ and $∑ j E j P j therm = 0$⁠.

Figure 7 shows how the temperature (given as an energy $k B T$ in units of eV) and the chemical potential develop over time for a high fluence of 30  $J / cm 2$⁠. Although reaching higher temperatures, these results are in qualitative agreement with studies treating the nickel valence system heated with optical lasers.^55,56

Because photons can “kick out” electrons from arbitrary states in the valence system, it becomes necessary to track the total kinetic energy of the free electrons, even though the distribution of energy among these electrons is not tracked. The rate of energy transfer to the free electron bath is described as the sum of non-resonant absorption and Auger decay processes, each multiplied by their respective energies (first and second line, respectively). The energy of the secondary Auger electrons is calculated as the average energy of all electrons other than the primary Auger electrons, as those drop to the core level. Finally, each electron that leaves R_free reduces the energy of the bath by the average energy, which is $E free R free$⁠.

Explicitly tracking this energy also enables us to demonstrate the conservation of energy within the simulation; since there is no channel that allows energy to leave the sample, the energy held in the electronic sub-systems matches that of the absorbed photons at all times. Figure 8 shows the internal energy of the electronic subsystems over time, integrated over an area of 1 nm² and the full 20 nm thickness of the sample. The comparison of the energy in various sub-systems demonstrates how quickly the energy of the photon pulse is transferred to valence excitations. Furthermore, observing the energy conservation has proven to be an invaluable tool to select sufficiently fine binning in time, space and energy, as it is particularly sensitive to the accumulation of numerical errors. For example, the calculation for the homogeneous illumination with 30  $J / cm 2$ shows a cumulative error in energy of 0.35%. This is in contrast to the electron conservation, which is strictly kept to machine precision level due to the symmetric way the process terms are arranged to form the time differentials.

There are four categories of parameters: First, resolution parameters for the number of steps in time, space, and energy are chosen as a compromise between calculation time and numerical error; second, experimental parameters such as the pulse duration, peak fluence, and bandwidth of the interacting photons reflect experimental conditions; third, ground-state properties such as the atomic density and number of states and electrons, as well as Auger-decay limited core-hole lifetime are drawn from published literature, while the resonant and non-resonant absorption lengths are derived from the ground-state spectrum as described below. Fourth, the model-inherent phenomenological parameters are the valence thermalization time τ_th and electron scattering time τ_scatt, which are varied to achieve the best match to the experimental results. A list of the relevant parameters and the chosen values for the present calculations is shown in Table III.

The parameterization of several ground-state properties deserves further comment. First, the DOS (motivated at the start of this section) was subdivided into bins m_j of varying size, favoring a fine resolution for the bound states.

The size of the energy bin that is resonantly coupled to the core level by the incident photons is chosen such that it represents the interaction bandwidth of the photons. The interaction bandwidth is calculated as the convolution of the bandwidth of incident photons (i.e., the resolution of the experiment) and the natural linewidth of the core excitation. We further account for the final state broadening of excitations into less tightly bound states by enlarging the interaction bandwidth by 0.1 eV per eV above the Fermi level. Furthermore, the resonant and non-resonant absorption lengths are derived from the ground-state spectrum. We treat the non-resonant absorption length as constant, i.e., independent of photon energy, and derive it from the pre-edge absorption level. The transition matrix element of a core–valence transition exhibits a resonant enhancement close to the absorption edge, which translates into an energy dependence of the resonant absorption length. Above the Fermi level, where the DOS is unoccupied in the ground state, the resonant absorption length is encoded in the ground-state absorption spectrum. As the transition matrix element from the core level to states below the Fermi level is experimentally not straightforward to access, we use the approximation that the energy-dependence of the transition matrix element is symmetric around the Fermi level. To derive the resonant absorption length, the non-resonant absorption level is subtracted from the spectrum and line broadening is accounted for by deconvolution with a pseudo-Voigt-profile of 50% Gaussian and Lorentzian share and a width of 640 meV FWHM, representing 420 meV broadening from the experimental resolution and 480 meV from the core-hole lifetime.⁵⁰ The deconvolved resonant absorption spectrum above the Fermi level is then mirrored around the Fermi level and the discontinuity within 320 meV around it is reconstructed with cubic interpolation. This results in the mirrored resonant absorption spectrum shown in Fig. 9, which is used as the resonant absorption length parameter $λ i res$⁠. Note that the results of simulated spectra were finally re-convolved with the same pseudo-Voigt-profile to simulate the same experiment.

In Fig. 10, we show the measured non-linear x-ray absorption spectra labeled (I) together with sets of simulated spectra computed for different sets of parameters.

The first set of simulated spectra represents the best match with the experimental conditions with the parameters shown in Table III and is labeled (II), while the consecutive sets, labeled (III) to (VI), demonstrate how the results change when individual parameters are modified. We present this set of simulations to show how our model may be used to understand the relation between the non-linear changes and various parameters. To fit the experimental results, however, only the parameters τ_th and τ_scatt are treated as unknowns, while all other parameters are known experimental or ground-state parameters.

In the best matching simulation (II), the experimental observations of a redshifted rising edge, increased pre-edge absorption, as well as reduced absorption at and beyond the absorption peak, are reproduced. However, the saturation effect at the resonance is over-estimated and the lack of measured saturation around 852 eV cannot be explained by our model.

The next simulation (III) was performed without non-resonant absorption. While this eliminates the pre-edge absorption rise, some shift of the absorption onset is still predicted within the original peak, while the spectra above the resonance onset behave similarly to the best-matching simulation.

With simulations (IV) and (V), we demonstrate the effect of prolonging the thermalization time or shortening the pulse duration by a factor of ten, respectively. Both have the similar effect of strongly reducing the peak shift and a moderate decrease in saturation beyond the resonance.

Finally, we show non-linear spectra (VI) where the free electron scattering process was eliminated (⁠ $τ scatt = ∞$⁠). This prevents the majority of the excitation energy from entering the valence system and thus drastically reduces valence heating. This causes the pre-edge absorption rise and rising edge shift to vanish nearly entirely, and absorption-decrease due to saturation at and above the resonance is also strongly reduced. These differences underscore the importance of electronic scattering cascades to these phenomena.

While an extensive study of the correlations between the various model parameters and the observed effects is beyond the scope of this work, this brief parameter study allows an interpretation of how redshift and pre-edge absorption rise occur: Let us consider a case where the incident photon energy is slightly below the absorption edge. Initially, only non-resonant absorption transfers energy to the sample, specifically by creating free electrons in form of photo-electrons. This energy is then transferred to the valence band due to electronic scattering cascades, where it causes many small, non-thermal excitations. These are homogeneously distributed over all valence energies and, thus, only lead to a small increase in pre-edge absorption (see simulation IV with slow thermalization). However, this distribution of excitations develops toward the shape of a Fermi distribution over time, specifically with the thermalization time constant. If the pulses are long in comparison to the thermalization time, a significant number of states below the Fermi level becomes available for core–valence transitions during the pulse duration. This causes the observed increase in pre-edge absorption. Once the first empty valence states become available at the current photon energy, resonant absorption begins to occur in addition to the non-resonant absorption. The additional resonant absorption leads to more free electrons from Auger-decay, which in the same way as the photo-electrons further contribute to secondary electron scattering and thermalization. Because thermalization mostly creates free states just below the Fermi level, this process of self-enhancing rise of overall absorption manifests in a fluence-dependent shift of the rising edge to lower energies.

As demonstrated in the parameter study above, this kind of analysis can provide a straightforward interpretation of how non-linear changes to the absorption spectrum emerge, and what each change says about the development of the electronic populations. While the model already demonstrates good agreement with experimental results, below we discuss the limitations of our approach and possibilities for expanding or refining it.

One particular physical effect that has been observed in warm dense matter is ionization potential depression. Here, the absorption edge is lowered for those atoms which have already been ionized, and manifests typically in discrete peaks emerging before the edge, see, e.g., Ref. 75. This signature is not visible in the spectra, likely because the non-resonant absorption at the photon energies where this could have been observed was not strong enough to ionize a significant number of atoms. It was, therefore, also omitted in the simulations. It may, however, be relevant when simulating spectra measured with higher fluence, a broader bandwidth or two pulses with separate photon energies (e.g., below and above the edge). Including it in the simulation could be possible by adding additional resonant absorption channels with a shifted core level, scaling their rate proportional to the fraction of ionized atoms.

In the model as described here, the DOS is not only assumed constant but also does not differentiate between spin-up and spin-down states, although we are treating a magnetic material. Splitting up the DOS in spin-up and spin-down states would allow for the inclusion of angular momentum conservation and transfer in the various scattering rates.

While phonon-mediated energy transfer is usually negligible in the first 100 fs after excitations, spacial coupling due to ballistic electrons can become relevant if the electron mean free path becomes similar to the sample thickness or the absorption length.^76,77 In this case, the model assumption of strictly local electronic excitations would have to be carefully reevaluated. While an additional process for ballistic electrons leaving the sample or interacting across the interface with a substrate could be introduced, directly calculating spatial coupling of electronic excitations within the sample in three dimensions would dramatically increase the computational complexity.

Fluorescent decay should be accounted for when moving to harder x rays.⁴⁹ While introducing the decay channel itself would be simple, accurately accounting for reabsorption may be less straightforward, since the model only propagates light in one direction.

Furthermore, the thermalization time τ_th is used in this work as a global fitting parameter, although electron thermalization times have been suggested to depend on electronic temperature.^52,53 Since the electronic temperature and target distribution are calculated every time step, an arbitrary dependence could be easily introduced, albeit with the necessity of additional fitting parameters or predetermined knowledge about the scaling of electronic temperature.

The DOS is dominated by the crystal lattice, which is typically stable on the sub-100 fs timescale. However, recent (time-dependent) density functional theory (TD)DFT calculations show that electronic processes causing sub 100 fs demagnetization via spin–orbit coupling can also lead to modifications of the DOS due to the presence of local electronic correlations modeled by introducing an onsite Hubbard correlation U to the mean-field Hamiltonian.⁴³ Since the rate model approach is generally not suited to calculate the DOS, incorporating modifications to the DOS would require a close interplay of the rate model with (TD)DFT calculations, leading to an ultimately much more complex approach reminiscent of models developed for the study of radiation-induced damage mechanisms.^23,24

Furthermore, we derive the interaction bandwidth (the valence energy range to which core states can be resonantly coupled by incident photons) as a convolution of instrumental resolution and the lifetime of the core excitation, i.e., the Auger lifetime. The final state lifetime broadening of excitations into continuum states is described as a continuous broadening of 0.1 eV per eV above the Fermi level.⁷⁸ It is however reasonable to expect that the final state lifetime is further shortened at higher fluences, due to increased rates of both electronic scattering and stimulated emission; the latter would be particularly relevant at the resonance peak. Such further broadening would cause more valence states to be available for resonant interaction and reduce the observed saturation effect while increasing the number of core-holes that may be created by high fluences.⁷⁹ Taking into account for a fluence-dependent broadening of the interaction bandwidth in a refined rate model might help to remedy the overestimation of the saturation effect at the absorption peak in the presented calculations.

Another candidate for further refinement is an energy dependence of the electronic scattering rates. The presented model uses fixed rates for thermalization and scattering cascades, which both act on all valence states indiscriminately. This description would be especially inadequate when applied without modifications to a bandgap material. An advanced model could describe both the thermalization rate of the valence band and the energy of excitations from scattering cascades in an energy-resolved manner.⁸⁰ 

While such refinements may seem attractive, a core strength of the rate model approach is its relative simplicity and computational tractability, as well as the use of known ground-state parameters, which supports a straightforward physical interpretation. It is ultimately a mostly classical, phenomenological model which offers a complementary approach to ab initio calculations. Every added complexity should therefore be weighed against its relevance, as a simpler model facilitates a meaningful interpretation and avoids introducing ambiguity in the results due to correlations between redundant input parameters.

Since the model operates on widely applicable principles, we expect that it may be applied to a wide range of materials with some predictive power, while the limitations described above apply. The results will especially deviate from observations wherever multi-particle effects, such as electron correlations or quasiparticles become relevant.

This appendix provides details on the experiment, namely, about the procedure of the data acquisition (Subsection 1 of  Appendix B), the calibration of the main detector (Subsection 2 of  Appendix B), the classification of individual FEL shots (Subsection 3 of  Appendix B) and finally on the determination of the effective FEL spot size (Subsection 4 of  Appendix B).

The experimental x-ray absorption spectra presented here were collected in the scope of a community proposal as the first user-beamtime at the SCS instrument. The intensities of the two beams generated by the beam-splitting zone plate were recorded using a FastCCD detector³⁴ with 1920 × 960 pixels of 30 × 30  $μ m 2$⁠. The intensities of both beams were integrated over a region of interest (ROI) corresponding to 350 × 350 pixels each to retrieve the signal and reference intensities for each FEL shot. The high beam divergence due to the zone plate focusing distributed the signal on a 4 mm wide square on the detector 1 m downstream of the sample, thus greatly decreasing the fluence incident per detector area in order to avoid detector saturation. The gas-attenuator was filled with varying low pressure of nitrogen gas to regulate the transmission through the beamline to the required fluence. We refer to low-intensity spectra if the fluence was consistently below the sample damage threshold and the full measurement could be recorded on a single spot, without scanning the sample.

For measuring high-intensity spectra, the fluence often exceeded the material damage threshold, creating shot craters and sometimes causing larger fractures in the support membrane. For measurements at these fluences, the sample holder was scanned at a speed of 0.5  $mm s − 1$⁠. Therefore, only about 50% of all FEL shots were transmitted through the windows; in the other cases, one or both beams were blocked or clipped by the frame. The membranes were arranged on the frame in a periodic pattern of two rows of sample and two rows of reference membranes with a distance of 1 mm between rows. This ensured that the two FEL foci always impinged on one sample and one reference membrane; a third row was unused in between. Therefore, every time the currently scanned rows were incremented, the upper beam would switch from probing reference membranes to probing sample membranes or vice versa, while the opposite holds for the lower beam. To prevent detector saturation, an additional aluminum filter of about 13 μm thickness was installed in front of the detector during these measurements.

The temperature of the FastCCD rose consistently during operation and required cooldown periods in between measurements, leading to the temperature varying between −27 and −5 °C not only over time during the measurement, but also spatially over the detector area. The detector dark signal, as well as the gain coefficients for the three gain settings between which the detector pixels switch automatically, depend on the detector temperature. This made it necessary to reconstruct a temperature-dependent gain calibration. The three temperature-dependent background levels were drawn from dark images collected at various temperatures for each gain setting; the gain coefficients for each setting were drawn from a statistical analysis of the observed gain switching thresholds, such that the calibrated histogram of pixel intensities becomes continuous over all three gain levels. While the calibration accounts for the temperature measured using a temperature sensor on the detector, spatial variations over the detector area remain. The primary effect of this temperature variation was a varying background signal, following a spatial exponential distribution between the detector center and rim, with a higher baseline near the detector center. To account for this, an estimated background signal was derived from the measurements themselves: For a running average of 100 images, the illuminated area was cut out and interpolated using fits to the background level in the non-illuminated area. This additional background variation was then integrated into the gain calibration described above. Furthermore, a mask of hot and dark pixels with irregular behavior was generated from separate measurements and the respective pixels were excluded from the analysis. Despite these corrections, the detector inhomogeneities constitute a significant part of measurement uncertainty in the presented spectra. In particular, the uneven warming of the two detector halves on which the upper and lower beam impinged has the potential of introducing systematic uncertainties. Thus, the shot-sorting algorithm described below was applied separately for all rows where the sample was in the upper beam and the reference in the lower and then to all rows where the orientation was reversed. Fortunately, differences between the temperature of both detector hemispheres affect the two equally sized groups of data (sample up and sample down) with equal and opposite magnitude. Therefore, possible systematic deviations are eliminated in the average over both groups and instead contribute to the statistical uncertainty which is represented in the error shown in Fig. 3 of the paper.

It is useful to mention to potential future users of the beamline that the difficulties with hardware calibration and data pipeline described above were a feature of this first user experiment and have since been addressed. The evolved version of the here described setup is now in routine operation and has since been shown to operate with a sensitivity close to the photon shot-noise limit.³² 

Furthermore, whenever the sample or reference membrane was torn due to a particularly intense shot, subsequent shots sometimes impinged on the torn rim of the sample, possibly at an angle to the membrane surface, or did not hit any sample material at all. Shots affected in this way were not always trivial to identify from any single measurement parameter, which lead to the following procedure to identify and exclude faulty FEL shots: First, the detector image in the ROI around one beam was compared to an extended ROI around the other beam using a normalized two-dimensional cross correlation algorithm. For this, the images were first smoothed by convolution with a Gaussian kernel to remove the influence of the rough surface structure of the aluminum filter before applying the cross correlation function, both algorithms implemented in the scikit-image⁸¹ package. This procedure yields a correlation coefficient and a displacement vector. The correlation coefficient was used as an indicator that both beams were transmitted through a window without significant differences in the wavefront.

Since data acquisition of the motor encoder for the position transverse to the scanning direction was not working, the real path of the beam over the sample frame was reconstructed by combining knowledge of the beam position along the scanning direction and the manual notes in the laboratory book with the correlation coefficients between both beams on the detector. In the later analysis of the damaged samples, this allowed associating specific shot craters scrutinized with SEM to specific FEL shots.

From the FEL shots which hit sample windows according to this reconstruction, outliers were dropped if either the correlation coefficient dropped below 85% or the displacement vector deviated by a significant margin (manually calibrated for each measurement setting) from the expected beam splitting. Likewise, extreme outliers in the ratio of reference to sample intensity were also dropped at this stage. These criteria proved to be largely redundant as they mostly agree with each other on which events to exclude. Using all of these criteria, events where one or both beams are blocked or clipped can be excluded reliably. In Fig. 11, these excluded points are shown in gray. Each dot represents one FEL shot; its y-position is the logarithmic ratio between sample and reference intensity; the x-position is the photon energy setting of the monochromator, while the gray-scale encodes the Pearson correlation coefficient $C corr$ of the two regions of interest.

From the distribution of the remaining events, shown in color, it is obvious that further classification is needed. This is because shots onto a damaged or missing membrane can produce spots with good correlation and the expected displacement vector. These can only be distinguished by the fact that the apparent optical density deviates unreasonably from the expected value at the given photon energy. However, this transmission ratio is also the quantity of interest for the final spectrum. To disentangle events affected by prior sample damage from true measurements, an iterative approach using a Bayesian Gaussian mixture model was utilized: To start, an initial guess for a spectrum is computed using the events that passed filter conditions described above and shown as a green line in Fig. 11. Then, for each FEL shot, the deviation of the logarithmic ratio between sample and reference intensity from the initially guessed spectrum is computed. This allows for generating a histogram of these absorbance deviations from the initial guess. Assuming a good guess of the initial spectrum, one may expect that the “good” FEL shots are normally distributed around zero deviation, while shots affected by various sources of uncertainty are distributed with some other distribution, dependent on the type of uncertainty. Thus, the histogram is fitted with four Gaussian distributions which are used as prior probability distributions, the first of which corresponds to good FEL shots. Then, the posterior probability of belonging to the category of good shots is computed for each FEL shot. These posterior probabilities are used as statistical weights to compute an improved guess of the measured spectrum. This improved guess was convoluted with a Gaussian kernel to prevent over-fitting before using it as a new initial guess for the Gaussian mixture model. The spectrum is considered converged when the average change per iteration anywhere in the spectrum is less than 5  $μ OD$⁠. This procedure converges (if at all) within a few (between 3 and 20) iterations to a solution that is robust even against a strong variation of the initial guess. In Fig. 11, the shots analyzed in this way are shown in colors encoding the final posterior probability $P ok$⁠, indicating the estimate of the model about the validity of each shot. The final guess of the average spectrum is shown as an orange line.

This final procedure rejects outliers purely based on their deviation from the expected value, which requires further justification: In this case, the procedure should return valid results under the condition that all valid FEL shots at a certain photon energy measure a transmission value within a single continuous range. The width of the distribution is resulting from a combination of non-linear changes and measurement noise. This condition must be fulfilled in the present case as long as the fluence-dependent transmission curve of the sample is continuous and a continuous range of pulse energies is contained in the dataset (which is the case for SASE fluctuations). Apart from this logic, the rejected outliers do not appear systematic, further supporting the applicability of the algorithm.

To characterize the effective focal size, the reference membranes were analyzed using SEM. Figure 12 shows two SEM images: one overview image of an entire membrane and one example of a high-resolution image of a single imprint. Such high-resolution images were taken of 85 selected spots which were associated with the corresponding FEL photon diagnostics data by matching the reconstructed movement of the sample stage to the pattern of imprints on the sample. The reference membranes were chosen for the imprint analysis since their x-ray absorbance can be considered constant for the scanned photon energies.

Since the zone plate is a diffractive element, its properties such as focal length are energy-dependent. Furthermore, the fluence-dependent spectra presented in this paper are combined from measurements at two distinct object distances of 250.84 and 252.14 mm.

Thus, the SEM images were grouped into six groups (see Table IV) by focal length and photon energy of the associated shot, and the total damaged or ablated surface area was determined for each shot. For each group, the minimum FEL pulse energy at which damage is observed on the reference membranes was determined by a Liu's plot⁸⁵ as shown in Fig. 13: The ablated area was plotted over the logarithm of the shot energy and a linear fit was applied to the shots with less than 1  $μ m 2$ ablated area to determine the pulse energy damage threshold.

Following the concept outlined in Ref. 83, the pulse energy of all shots was then normalized using this damage threshold to derive the normalized fluence level f(S) as a function of the ablated area S. If the spot size were Gaussian, the function f(S) should be a simple negative exponential, and the effective area would correspond to the value of S where f(S) equals $1 / e$⁠, which is indicated in Fig. 13(b) as a horizontal dotted line. For the given non-Gaussian case, the function f(S) was fitted with a modified exponential $f ( S ) = e ( − a S b )$⁠, yielding for each group of shots fit parameters a and b and their uncertainties σ_a and σ_b. The effective area is then calculated as the integral of f(S). For the estimate of the uncertainty shown as error bars in Fig. 13(c), the integral was also evaluated for $a + σ a$ and $b + σ b$ as well as $a − σ a$ and $b − σ b$⁠. Groups (2) and (3) each contained only a single FEL shot with visible damage. While this does allow for an estimate of the damage threshold [compare Fig. 13(a)], the function f(S) shown in (b) is barely determined by fitting only to the normalization point and a single further point. Thus, no mathematical uncertainty estimate can be given and the resulting points (2) and (3) shown without error bars in Fig. 13(c) should be considered as tentative estimates.

The ray optics calculation tracks the beams in the vertical and horizontal planes. The focus area was calculated from the ray optics calculations assuming an elliptical beam shape. It revealed a slightly astigmatic focus on the sample, due to the beamline optics using separate horizontal and vertical focusing, thus illuminating the zone plate with slightly non-uniform beam divergence. The separate location of horizontal and vertical foci leads to the two minima visible in Fig. 13. Based on the measurements, a minimal beam waist radius of 80 nm was imposed on the ray tracing calculations, accounting for imperfections in beam quality and optics. Normalizing the beam damage thresholds by the thus determined effective area, we find on average that the reference membranes where damaged above fluences of about 10.5 J/cm².

The fluences for the high-fluence spectra shown in this paper are calculated based on the pulse energy measured by the XGM behind the monochromator, multiplied with the efficiency of the focusing optics (80%)⁸⁶ and zone plate (9%) and divided by the effective area derived with the presented ray optics calculations. Since the XGM was calibrated for high pulse energies, the fluences for the low-fluence spectra are based on the reference intensities measured on the CCD and calibrated to be consistent with the reference intensities measured for the high-fluence data, accounting for the calculated transmission of the aluminum filter.