

\chapter{Extraction of Electronic Structures from RIXS Measurements}






Understanding the electronic structure of a material is essential for predicting its physical and chemical properties, such as conductivity and reactivity, and in particular its interactions with electromagnetic radiation. In this context, as shown in this chapter, the electronic structure of a material is embedded in the cross-section of x-ray photons resonantly scattered by that material. Resonant inelastic x-ray scattering (RIXS) spectra, therefore, encode a wealth of information about the material’s density of states, albeit in a complex way. Extracting the electronic structures from RIXS measurements, which is connected to the inversion of a function operator, is feasible but at the same time challenging and, more importantly, unstable due to factors like experimental noise. This chapter presents a method that combines neural networks (NNs) and physics modelling, devised to perform this extraction. Additionally, this approach is capable of deconvolving the source spectra (and potentially the spectrometer's instrument function), thereby removing the major broadening effects and providing high-resolution measurements of the DOS. The method described here is generalisable to other physical processes, provided they satisfy certain conditions, which are discussed in detail. Since the RIXS cross-section also depends on thermodynamic parameters - particularly the sample temperature $T$ - it is possible to also extract $T$ from RIXS spectra, provided the sample's DOS is already known (e.g. through DFT calculations). The accuracy of the temperature estimate relies on various experimental parameters and material characteristics, including its electronic structure.

This chapter begins by describing the general theory of RIXS, focusing on the approximations relevant to our scientific cases and deriving the theoretical expressions used throughout this thesis to model our experiments. For completeness and to enable instructive comparisons, I then review alternative techniques for measuring electronic structures, including XAS (divided into EXAFS and XANES), XES, and x-ray Raman scattering (XRS), highlighting their respective advantages and disadvantages in relation to RIXS. The second part of this chapter provides an overview of the most common machine learning (ML) methods applied within physics, followed by a description of my approach for the extraction of physical quantities embedded in experimental measurements, with an emphasis on its capabilities and limitations. This method is then tailored to the RIXS case and validated on synthetic data. Application of the method to experimental data will be discussed in the final two chapters.




% do I cite the STEP paper




\section{\label{RIXS theory} RIXS Theory}







Resonant inelastic x-ray scattering is an experimental technique in which one scatters x-ray photons inelastically off matter (\textit{photon-in/photon-out} type of spectroscopy). During this process, the photons exchange energy, momentum, and polarisation with the probed physical system. By measuring how these quantities change upon scattering, one can extract information from the analysed matter (e.g. electronic energy levels). RIXS is a resonant technique in which the energy of the incident photons is chosen such that it coincides with, and hence resonates with, one of the atomic transitions of the sample. The resonance greatly enhance the inelastic scattering cross-section and therefore simplify the spectroscopic analysis of the system. Here, I start from a general treatment of RIXS and then specify the results to the specific cases of interest. 

The RIXS cross-section can be derived from standard quantum mechanics perturbation theory. The Hamiltonian of the full system (matter + radiation) is given by~\cite{landau}:

\begin{align}
& H = H_{0} + H_{int} \\
& H_{int}=\frac{e}{mc}\sum_j \textbf{p}_j \cdot \textbf{A}(\textbf{r}_j)+\frac{e^2}{2mc^2}\sum_j |\textbf{A}(\textbf{r}_j)|^2 = O(A) + O(A^2) \ ,
\label{RIXS theory 1}
\end{align}

\noindent where the index $j$ refers to the electrons of the sample, $H_{0}$ is the Hamiltonian of the full system neglecting the matter-radiation interactions, and $\textbf{A}$ is the quantised electromagnetic vector potential~\cite{landau}. In order to compute the differential cross-section of the scattered photons, one needs to calculate the following matrix element \cite{sturm}:

\begin{equation}
\langle i,\lambda_1|H_{int}|f,\lambda_2\rangle \ ,
\label{RIXS theory 2}
\end{equation} 

\noindent where $i$ and $f$ are the initial and final state of the probed system, and $\lambda_1$ and $\lambda_2$ are the initial and final state of the photon. The state $\lambda$ includes the photon wavevector $\textbf{k}$ (with associated energy $\hbar\omega = c|\boldsymbol{k}|$) and the polarisation vector $\epsilon$. This matrix element is associated with the scattering process $(i+\lambda_1 \longrightarrow f+\lambda_2)$ and can be computed with perturbation theory. The perturbation expansion is carried out in terms of $A$ and must be brought up to order $O(A^2)$ to retrieve the RIXS process\footnote{Actually, the expansion is done in an adimensional parameter of the form $(\kappa A)$, with $[\kappa]=[A]^{-1}$~\cite{berk}. This parameter is small for the typical photon energies and electrons momenta occurring in our systems.}. To do that, one must use quantum perturbation theory (or alternatively the Dyson series) to the first order for the $O(A^2)$ term in $H_{int}$ and up to second order for the $O(A)$ term. The first order of the $O(A)$ term describes the absorption/emission of a photon and therefore gives no contribution to the scattering process.

After performing the expansion, the total double differential cross-section of the scattering process is found to be (see~\cite{berk} for a step-by-step derivation):

\begin{equation}
\begin{split}
\frac{d^2\sigma}{d\Omega_{\textbf{k}_2} d\omega_2}= \hbar \frac{\omega_2}{\omega_1}\biggl(\frac{e^2}{mc^2}\biggr) & \sum_f \  \delta(E_f-E_i+ \hbar\omega_2-\hbar\omega_1) \times \\
\biggl|\langle f|\rho_{\textbf{k}_1-\textbf{k}_2}|i\rangle &(\boldsymbol{\epsilon}_1\cdot\boldsymbol{\epsilon}_2) + \\
 \frac{1}{m}\sum_n & \biggl(\frac{\langle f|\textbf{p}(\textbf{k}_2)\cdot\boldsymbol{\epsilon}_2|n\rangle\langle n|\textbf{p}(-\textbf{k}_1)\cdot\boldsymbol{\epsilon}_1|i\rangle}{E_n-E_i- \hbar\omega_1}+ \\
 & \frac{\langle f|\textbf{p}(\textbf{k}_1)\cdot\boldsymbol{\epsilon}_1|n\rangle\langle n|\textbf{p}(-\textbf{k}_2)\cdot\boldsymbol{\epsilon}_2|i\rangle}{E_n-E_i+ \hbar\omega_2}\biggr)\biggr|^2 \ ,
\end{split}
\label{RIXS theory 3}
\end{equation}

\noindent where $\ket{i}$, $\ket{n}$, and $\ket{f}$ indicate the initial, intermediate, and final electronic states of the probed system, with respective energies $E_i$, $E_n$ and $E_f$. The indexes `1' and `2' refer respectively to the incoming and outgoing photon (e.g. $d\Omega_{\textbf{k}_2}$ is the differential solid angle around the direction of the scattered photon). Here, the following notation was used:

\begin{align}
    & \rho_{\textbf{k}} = \sum_j \exp(-i\textbf{k}\cdot\textbf{r}_j) 
    & \textbf{p}(\textbf{k}) = \sum_j \textbf{p}_j\exp(-i\textbf{k}\cdot\textbf{r}_j) \ .
\end{align}

\noindent These two operators are often used in their multipole approximations~\cite{jackson2021classical}, expanding the exponentials in power series. This approximation is valid as long as the photon wavelength is much larger than the spatial extension of the electron wavefunctions (of order $a_B/Z$ for bound electrons, with $a_B$ the Bohr radius), which is usually the case. Eq. (\ref{RIXS theory 3}) is also known as the Kramers-Heisenberg dispersion formula, named after the two physicists who first derived it. The summation $\sum_n(\cdots)$ is the second order contribution of the $O(A)$ term in $H_{int}$, and it describes the absorption of a photon followed by the emission of another photon (first addend), and the reverse process (second addend). These two terms represent, respectively, RIXS and a non-resonant inelastic x-ray scattering (IXS) processes, once the incoming photon energy is properly tuned. The other term inside the modulus is instead the first order contribution of the $O(A^2)$ piece of $H_{int}$ and it describes both x-ray Thomson scattering~\cite{crowley2014quantum} and XRS, depending on the specific excitations probed (e.g. bound-bound transitions). In terms of Feynman diagrams, these three types of scattering are illustrated in Fig. \ref{Feyn_diag}.
\begin{figure} [t]
\begin{center}
  \includegraphics[width=\textwidth]{Figures/diagrams.pdf}  
  \caption[Feynman diagrams for different scattering processes]{\small{Feynman diagrams of (XRTS+XRS), RIXS, and non-resonant IXS processes. In RIXS the initial state $i$ pass through the intermediate state $n$ by absorbing a photon and then end up in the final state $f$ by emitting a photon. The subprocesses of emission and absorption are not individually energy-conserving, making the intermediate state a virtual state.}} 
   \label{Feyn_diag}
\end{center}
\end{figure}

In RIXS analysis, the incident photon energy is taken close to an electron excitation threshold, so that the first term in the sum over $n$ becomes dominant as the denominator vanishes. If the lifetime of the intermediate state is considered, the energy of the state is modified $E_n \to E_n+i\Gamma_{nf}$, with $\Gamma_{nf}$ inversely proportional to the lifetime, and the divergence is removed.
The RIXS double differential cross-section is then given by:

\begin{equation}
\begin{split}
\frac{d^2\sigma}{d\Omega_{\textbf{k}_2} d\omega_2}= \hbar \frac{\omega_2}{\omega_1}\biggl(\frac{e}{mc}\biggr)^4\sum_f \delta(E_f-E_i+\hbar\omega_2-\hbar\omega_1) \times \\
\biggl|\sum_n\frac{\langle f|\textbf{p}(\textbf{k}_2)\cdot\boldsymbol{\epsilon}_2|n\rangle\langle n|\textbf{p}(-\textbf{k}_1)\cdot\boldsymbol{\epsilon}_1|i\rangle}{E_n+i\Gamma_{nf}-E_i-\hbar\omega_1}\biggr|^2 \ .
\end{split}
\label{rixs cross-section}
\end{equation}

\noindent While the $\delta$ function ensures the energy conservation of the whole scattering process, the sum over $n$ describes the absorption of a photon $(\textbf{k}_1,\boldsymbol{\epsilon}_1)$, causing the system to transition from $\ket{i}$ to a generic intermediate state $\ket{n}$, followed by the emission of a photon $(\textbf{k}_2, \boldsymbol{\epsilon}_2)$, relaxing the system to the state $\ket{f}$. The transition matrix elements inside the sum give the probability for these two subprocesses, which are then summed over all the possible intermediate states. Note also that, even though these subprocesses are not required to conserve energy, the denominator penalises the non-energy conservation.

In this thesis, RIXS technique is applied to the study of pure Fe, iron oxides (FeO and Fe$_2$O$_3$) and gadolinium (Gd) under extreme thermodynamic conditions. The discussion now focuses on iron and iron oxides, with modifications for the gadolinium case discussed directly in Chapter 4. For RIXS to occur, as mentioned earlier, the incident photon energy must be set close to an absorption edge. For the experiments with iron compounds, this energy was tuned just below the Fe K-edge ($|\epsilon_K|=7112$ eV). Working below the edge allows one also to suppress direct absorption and the consequent radiative recombinations (e.g. $K_{\alpha}$ decay), leaving RIXS as the predominant process\footnote{If the incident photon energy is larger than $|\epsilon_K|$, that is the energy necessary to ionise a K-shell electron, the K-shell ($1s$ states) would be filled with holes, which would then be eliminated mainly by electrons decaying from $2p$ states via radiative recombination. Since these two processes combined have a larger cross-section than RIXS, one could risk weakening the RIXS signal.}. Under these circumstances, there are two dominant RIXS channels for these systems, which are schematically depicted in Fig. \ref{Channels}.

\begin{figure}[t]
\begin{center}
    \includegraphics[width=0.8\textwidth]{Figures/RIXS scheme.pdf}
    \caption[Scheme of the RIXS process in iron compounds]{\small{Scheme of the RIXS process in iron compounds. This is composed of an energy non-conserving absorption, which excites a K-shell electron (left), followed by an energy non-conserving emission produced by the decay of a L-shell electron into the K-shell core hole (right). The overall process conserves energy. The accessible states are the vacant states in the valence bands (first channel) or, for finite temperature systems, the thermally ionised states in the M-shell (second channel).}}
    \label{Channels}
\end{center}
\end{figure}

In the left-hand diagram, the initial state, which is assumed to be thermalised (i.e. the electrons' population follows a Fermi-Dirac distribution), absorbs a photon that excites, without conserving the energy, an electron from the K-shell ($1s$ states) into either a valence state (first channel) or an M-shell one (second channel). Then the hole in the K-shell is filled by an electron of the L-shell ($2p_{1/2}$, $2p_{3/2}$ states), as it is shown in the right-hand diagram. In order to conserve the energy over the whole process, the emitted photon has a lower energy than the K-L resonances (6392 and 6405 eV for iron compounds). Thus, in our case, the intermediate state $n$ is the initial state with a hole in the K-shell and an electron in an unoccupied valence or M-shell state, whereas the two possible final states $f$ are the initial state with an electron in the continuum or in the M-shell and a hole in one of the two L-states. Since one does not discern the final state of the system during the measurement, the total RIXS differential cross-section for these materials is given by Eq. (\ref{rixs cross-section}), with a sum over all the possible final states $f=1,2$. To activate the M-shell channel, the temperature of the system must be sufficiently high to have vacancies in the bound M-shell ($3p_{1/2}$, $3p_{3/2}$ states). This latter characteristic, as demonstrated shortly, can be exploited to estimate the temperature of the system.

A simpler form of Eq. (\ref{rixs cross-section}), specific to our cases of interest, will now be considered. The aim is to show that the electronic structure of a material can be efficiently extracted from its RIXS spectra. However, one is ultimately limited on the resolution of the extracted DOS by the instrument function of the setup, which is a few eV wide in our cases. One can therefore neglect the interactions between the core-hole potential and the electrons in the continuum (the intermediate state has a hole in the K-shell), as they introduce an uncertainty in the energy of the intermediate state given by the inverse of the core-hole lifetime~\cite{fetter2012quantum}, equivalent to a few eV in Fe~\cite{henke1982low}. These interactions are therefore responsible for modulations in the RIXS cross-section on an energy scale constrained by this lifetime\footnote{Potential future studies with higher resolution would require DFT calculations that incorporate a K-shell hole in the pseudopotential of the intermediate state.}. Starting from these considerations, Eq. (\ref{rixs cross-section}) is simplified through the following approximations and substitutions~\cite{humphries2020probing}:



\begin{itemize}
    \item Neglecting the core-hole potential, the transition matrix element of the $K_{\alpha}$ decay becomes independent of the intermediate state $n$, which has an additional electron in a continuum state, and therefore

    \begin{equation}
    |\bra{f} \textbf{p}(\textbf{k}_2)\cdot \boldsymbol{\epsilon}_2 \ket{n}|^2 = A_f,
    \label{approx1}
    \end{equation}
    
    \noindent where $A_f$ is the $K_{\alpha}$ line intensity. Note that in this approximation the emission is assumed to be isotropic;

    \item For analogous reasons, the lifetime of the intermediate state can be approximated to be dependent only on the final state $f$ and given by the core hole lifetime, $\Gamma_f$

    \begin{equation}
    \Gamma_{nf} \to \Gamma_f \ ;
    \label{approx2}
    \end{equation}

    \item One can substitute the summation over the intermediate states $\ket{n}$ with an integration over the accessible (via Pauli's principle) density of states $\rho(E_n)$

    \begin{equation}
    \sum_{n} \to \int_{\mathbb{R}} dE_n \ \rho(E_n)[1-f_{FD}(E_n;T)] \ ;
    \label{approx3}
    \end{equation}

    \item Assuming the absorption to be isotropic as well, and considering that the initial state of the electron is fixed (1s state), the absorption matrix element can be rewritten as

    \begin{equation}
    \bra{n} \textbf{p}(-\textbf{k}_1)\cdot \boldsymbol{\epsilon}_1 \ket{i}\vcentcolon= M(E_n) \ ;
    \label{approx4}
    \end{equation}

    \item 
    The energy difference between the initial and final states corresponds to the disparity between the energy of the ionised electron, $\epsilon_{\textbf{k}}$, and the binding energy of the L-shell electron $\epsilon_{L,f}$, that undergoes decay to occupy the core hole

    \begin{equation}
    E_f - E_i = \epsilon_{\textbf{k}} - \epsilon_{L,f} \ ;
    \label{approx5}
    \end{equation}

    \item 
    Likewise, for the intermediate state, the energy difference corresponds to the disparity between the energy of the ionised electron and the binding energy of the core hole state, $\epsilon_K$
    \begin{equation}
    E_n - E_i = \epsilon_{\textbf{k}} - \epsilon_K \ .
    \label{approx6}
    \end{equation}
    
\end{itemize}

\noindent Applying these substitutions to Eq. (\ref{rixs cross-section}) and identifying $E_n$ with $\epsilon_{\textbf{k}}$ (one-to-one correspondence), one obtains\footnote{One can bring the squared modulus inside the integral thanks to the $\delta$ function properties.}:

\begin{equation}
\begin{split}
\frac{d^2\sigma}{d\Omega_{\textbf{k}_2} d\omega_2}= & \hbar\frac{\omega_2}{\omega_1}\biggl(\frac{e}{mc}\biggr)^4\sum_f \int_0^{\infty} d\epsilon_{\textbf{k}} \ \delta(\epsilon_{\textbf{k}} - \epsilon_{L,f} + \hbar\omega_2 -\hbar\omega_1) \\
& \times [1-f_{FD}(\epsilon_{\textbf{k}};T)]\rho(\epsilon_{\textbf{k}})\biggr| \frac{\sqrt{A_f}M(\epsilon_{\textbf{k}})}{\epsilon_{\textbf{k}}-\epsilon_K - \hbar\omega_1+i\Gamma_f} \biggr|^2 \ .
\label{rixs cross-section 1}
\end{split}
\end{equation}

\noindent Finally, evaluating the integral and squaring the denominator, one finds:
\begin{equation}
\frac{d^2\sigma}{d\Omega_{\textbf{k}_2} d\omega_2}=\hbar\frac{\omega_2}{\omega_1}\biggl(\frac{e}{mc}\biggr)^4\sum_f\frac{[1-f_{FD}(\epsilon;T)]\rho(\epsilon)A_f|M(\epsilon)|^2}{[\hbar\omega_2-(\epsilon_{L,f}-\epsilon_K)]^2+\Gamma_f^2} \ ,
\label{rixs cross-section 2}
\end{equation}

\noindent where $\epsilon=\hbar\omega_1-\hbar\omega_2+\epsilon_{L,f}$. Under these approximations, the RIXS cross-section as a function of $\omega_2$ (RIXS spectrum) is approximately described as an absorption cross-section $[\propto \rho(\epsilon(\omega_2))|M(\epsilon(\omega_2))|^2]$ modulated by a Lorentzian centred in the K-L resonance ($\epsilon_{L,f}-\epsilon_K$) with a width given by $\Gamma_f$. The quantity $\rho|M|^2$ represents the density of states accessible from a 1s core state, and it is the part of the electronic structure, with the addition of a thermal modulation given by $(1-f_{FD})$, that RIXS is capable to investigate. This experimentally accessible quantity is denoted as the effective density of states, $\rho_{\rm eff}$, given by:

\begin{equation}
\rho_{\rm eff}(\epsilon) = [1-f_{FD}(\epsilon;T)]\rho(\epsilon)|M(\epsilon)|^2 \ .
\label{rho eff}
\end{equation}

\noindent Note that if either $T$ or $\rho|M|^2$ is known from other measurements or calculations, the other quantity can be estimated. For instance, when $T$ is already known, along with the matrix transition elements $M$, the actual density of states of the system can be determined. Since atomistic transition rules are still partially obeyed in solid systems~\cite{moore2020optical}, certain states of the DOS (e.g. 4s in this case) are suppressed in the density of accessible states (DAS).

Equation (\ref{rixs cross-section 2}) demonstrates that, with just one single shot, the RIXS analysis provides an entire spectrum, which contains information on the electron temperature via the shape of the Fermi distribution and on the DAS via the scattering intensity\footnote{In principle, relative binding energies between core shells—and thus the IPD—can also be investigated by examining global shifts in the RIXS spectrum.}. As discussed in the next section, this property of RIXS is a huge advantage with respect to the absorption/emission spectroscopy, where instead several shots with different XFEL photon energies across the edge are needed to reconstruct a spectrum. The wealth of information provided by RIXS in just a single shot, together with the short duration of XFEL pulse, makes this technique particularly well suited for time-resolved studies.

% This is due to the fact that one can simultaneously derive the material's physical properties and the corresponding thermodynamic conditions, which would be otherwise challenging to estimate in a dynamical process. 

Integrating Eq. (\ref{rixs cross-section 2}) over the solid angle detected for each outgoing energy $\Omega(\omega_2)$ and over energy bins of width $\Delta E$ given by the pixel size (see also Chapter 2), the intensity of scattered radiation measured by the detector at discrete energies $\hbar\omega_{2,i}=E_i$ is obtained:

\begin{equation}
I(E_i) = \int\limits_{E_i-\Delta E/2}^{E_i+\Delta E/2} \Omega(\omega_2) \overbrace{\frac{d\sigma}{d\omega_2}}^{S_{RIXS}} d(\hbar\omega_2) \ .
\label{finite_res}
\end{equation}

\noindent Here $\Omega(\omega_2)$ is a known function, calculable from the geometry of the experimental setup (see Chapter 2). Until this point, a monochromatic radiation incoming on the target has been considered. In a real experiment, the incident XFEL radiation possesses a finite spectral width. In this case, Eq. (\ref{finite_res}) is modified as:

\begin{equation}
I(E_i) = \int\limits_{E_i-\Delta E/2}^{E_i-\Delta E/2} d(\hbar\omega_2) \int\limits_{\mathbb{R}} d(\hbar\omega_1) \Omega(\omega_2) \Phi(\omega_1) S_{RIXS} (\omega_1,\omega_2) ,
\label{finite_FWHM}
\end{equation}

\noindent where $\Phi(\omega_1)$ is the XFEL spectrum. Recalling Eq. (\ref{rixs cross-section 2}), one observes that this integral has the characteristic form of a convolution of the scattering cross-section with the incident FEL pulse. It is intuitive that this convolution lowers the resolution with which one is able to extract $\rho(\epsilon)$ from $I(E)$ since different $\epsilon$ contribute to the same $I(E_i)$. Performing this deconvolution is a challenging task due to the complexity of Eq. (\ref{rixs cross-section 2}), combined with the typical stochastic and spiky profile of XFEL pulses~\cite{kayser2019core}. The approach to the solution of this problem is described in Section \ref{STEP for the RIXS Process}. To rigorously compare $I(E_i)$ with the experimental measurements, it ought to be convolved with the instrument function of the experimental setup.


% However, in the following, this broadening effect will be disregarded unless explicitly stated otherwise. 
% Because this effect was not focused on in our experiments
%The consideration of the instrument function in the RIXS model is the next logical step to improve our approach

It is finally worth noticing that in all the previous equations, the absorption in the target was not considered. This assumption is justified as long as the absorption does not change significantly with $\omega_1$ belonging to the probe bandwidth\footnote{This assumption holds true in our cases, where $\Phi(\omega_1)$ lies entirely below the absorption edge.}. In this case, there is just a multiplicative constant in front of the integral \ref{finite_FWHM}, which does not affect the results, as will be discussed later on.








\section{\label{Alternative Spectroscopic Techniques} Alternative spectroscopic techniques}






%In this section we will briefly outline other spectroscopic techniques that allow to investigate the electronic structures of materials, mainly focusing on x-ray absorption spectroscopy, which will be used for comparison with out results obtained trhough RIXS.






%\subsection{\label{Alternative spectroscopic techniques} X-ray absorption spectroscopy}




In absorption spectroscopy, the transmitted radiation intensity $I$ through a material is measured and compared to the incident intensity $I_0$. These two intensities are related by the Beer-Lambert law:

\begin{equation}
    I(\omega_1) = I_0(\omega_1) \exp{-n z \sigma(\omega_1)} \ ,
    \label{abs spec1}
\end{equation}

\noindent where $n$ is the number density of particles that absorb and scatter light, $z$ is the sample thickness, $\omega_1$ is the frequency of the photon source (assumed to be monochromatic here) and $\sigma$ is the total cross-section encompassing all scattering and absorption processes. Measuring the transmission function $\ln(I/I_0)(\omega_1)$ thus provides a measurement of $\sigma(\omega_1)$. By scanning the incoming photon energy $\hbar\omega_1$ across an absorption edge—a necessary condition to probe the electronic structure near the chemical potential—the photons cross-section becomes dominated by absorption ($\sigma(\omega_1)\approx\sigma_{{\rm abs}}(\omega_1)$). The x-ray absorption cross-section can also be calculated theoretically using methods similar to those for RIXS, as given in~\cite{brouder1990angular}:

\begin{equation}
     \sigma_{\rm abs}(\omega_1) = \frac{C}{\omega_1} \int_{\mathbb{R}} \ dE_f \delta(E_f - E_i -\hbar\omega_1)[1-f_{FD}(E_f)]\rho(E_f) |\langle f|\textbf{p}(-\textbf{k}_1)\cdot\boldsymbol{\epsilon}_1 | i \rangle |^2 \ ,
    \label{abs spec2}
\end{equation}

\noindent where $C$ is a physical constant. This cross-section, for an incident spectrum $\phi(\omega_1)$ with a finite bandwidth, becomes:

\begin{equation}
\begin{split}
     & \sigma_{\rm abs} = \frac{C}{\bar\omega_1} \int_{\mathbb{R}}\int_{\mathbb{R}} d(\hbar\omega_1) dE_f \ \delta(E_f - E_i -\hbar\omega_1)[1-f_{FD}(E_f)]\rho(E_f) |\langle f|\textbf{p}(-\textbf{k}_1)\cdot\boldsymbol{\epsilon}_1 | i \rangle |^2 \Phi(\omega_1) \\
     & = \frac{C}{\bar\omega_1} \int_{\mathbb{R}}d(\hbar\omega_1) \ [1-f_{FD}(E_i+\hbar\omega_1)]\rho(E_i+\hbar\omega_1) |M(E_i+\hbar\omega_1)|^2 \Phi(\omega_1) \\
     & = \frac{C}{\bar\omega_1} \int_{\mathbb{R}}d(\hbar\omega_1) \ \rho_{\rm eff}(E_i+\hbar\omega_1)\Phi(\omega_1) \approx \frac{C}{\bar\omega_1} \int_{\bar\omega_1-\Delta B/2}^{\bar\omega_1+\Delta B/2}d(\hbar\omega_1) \ \rho_{\rm eff}(E_i+\hbar\omega_1)\Phi(\omega_1) \ , 
\end{split}
    \label{abs spec3}
\end{equation}

\noindent where $\bar\omega_1$ is the central incoming photon frequency, pulled out from the integral by assuming a sufficiently narrow source bandwidth $\Delta B$ (as is typical in the analysed cases). Here, the notation introduced in the previous section was used. This equation shows that to obtain a spectrum with a form resembling that of Eq. (\ref{finite_FWHM}), which allows for the extraction of $\rho_{\rm eff}$ through deconvolution as with RIXS, multiple measurements are needed. Specifically, one has to shift the source spectrum $\Phi$ and compute $\sigma_{abs}$ for varying central energies $\hbar\bar\omega_1$. 

This requirement poses significant challenges in warm dense matter experiments, where thermodynamic conditions achieved through laser heating are difficult to stabilise and may vary widely between shots.  Consequently, there is a risk of assembling a spectrum from data obtained under varying conditions, compromising its reliability. Additionally, absorption spectroscopy requires crossing an edge. This not only introduces variability in thermodynamic conditions—as crossing the edge notably increases sample heating—but also causes fluctuations in signal intensity across shots, potentially leading to technical issues like detector saturation. In contrast, RIXS avoids these challenges by capturing an entire spectrum in a single shot while working below the absorption edge. However, because the cross-section for RIXS is smaller than for absorption, more shots are required to achieve the same SNR, which can reintroduce the issue of shot-to-shot variability, as seen for XAS. Nonetheless, rapid advancements in laser facility performance—especially in intensity and pulse control—are expected to mitigate these downsides of RIXS.

The expression for $\sigma_{abs}$ in the previous equation applies specifically to absorption in the states closest to the chemical potential, which are probed by what is known as x-ray absorption near edge structure (XANES). The structure of these states is highly sensitive to the local environment surrounding the atoms, allowing small variations in properties like coordination number and interatomic distance to be detected using XANES~\cite{newville2014fundamentals}. As the incoming photon energy increases, the kinetic energy of photoejected electrons rises, reducing their de Broglie wavelengths $k_{w}$ until they approach the scale of interatomic distances. At this point, the final state of the ejected electron $|f\rangle$ is altered by the interference of its wavefunction with the neighbouring atoms. This effect becomes relevant at about 40 eV above the edge. Above these energies, this interference modifies the absorption cross-section, introducing a modulation factor~\cite{newville2014fundamentals}:

\begin{equation}
    \Chi(k_{w}) = \sum_j \frac{N_je^{-2k_{w}^2\zeta_j^2}e^{-2R_j/\lambda(k_{w})}f_j(k_{w})}{R_j^2k_{w}}\sin(2k_{w}R_j+\delta_j(k_{w})) \ ,
    \label{EXAFS}
\end{equation}

\noindent where $j$ represents the individual coordination shell of $N_j$ identical atoms at approximately the same distance $R_j$ from the central atom, with scattering properties summarised by $f_j(k_{w})$ and $\delta_j(k_{w})$~\cite{newville2014fundamentals}). Here, $\zeta_j^2$ is the mean-square-displacement in the bond distance $R_j$ due to disorder sources like thermal agitation, and $\lambda(k_{w})$ is the mean-free-path of the ejected electron. This modulation factor, as mentioned earlier, is not valid for small $k_{w}$ values, but it is necessary for $\sigma_{\rm abs}$ to describe absorption when sufficiently far above the edge. This type of spectroscopy is known as extended x-ray absorption fine structure (EXAFS). Because of the $\lambda(k_{w})$ term and $1/R_j^2$ factor, the sum is limited only to the closest shells, making EXAFS a local probe of the material structure. In fact, since $f_j(k_{w})$, $\delta_j(k_{w})$ and $\lambda(k_{w})$ are known, one can fit this theoretical model of EXAFS (essentially given by $\sigma_{\rm EXAFS} = \sigma_{\rm abs}\Chi$) to the experimentally-measured absorption cross-section to determine values for $N_j$, $R_j$ and $\zeta_j$~\cite{newville2014fundamentals}. An example of an x-ray absorption spectrum for FeO is shown in Figure \ref{xas}. In developing the RIXS expressions, these distortions of the excited electron's final state $|f\rangle$ have been neglected. As for the simpler absorption process, this assumption holds as long as the interpretation of the extracted $\rho_{\rm eff}$ is limited to energies within approximately 50 eV above the chemical potential.

\begin{figure}[t]
\begin{center}
    \includegraphics[width=0.8\textwidth]{Figures/xas.pdf}
    \caption[Example of x-ray absorption spectrum for FeO]{X-ray absorption spectrum for FeO around the K-edge (blue curve), shown after baseline subtraction and edge jump normalisation. By subtracting the atomic absorption spectrum (red line), the resulting modulation function $\Chi(k_{w}(E))$ (shown in the inset) can be fitted to the theoretical model or Fourier transformed to extract the structural parameters of the material. The XAS spectrum is divided into two distinct regions (XANES and EXAFS), approximately separated 50 eV above the edge. Only the portion of the spectrum corresponding to the EXAFS regime is used to compute $\Chi$. Data is sourced from Ref.~\cite{newville2014fundamentals}.}
    \label{xas}
\end{center}
\end{figure}

In x-ray emission spectroscopy, core electrons of the probed system are excited by the incident photons, and then the photoemission, occurring during relaxation, is measured, revealing information about the system's energy levels~\cite{van2016x}. However, the large line broadening in the emission from continuum states—primarily due to collisional broadening—has limited the use of XES for investigating these states with detailed resolutions. Instead, XES has been primarily applied to analyse bound-bound transitions, enabling measurements of phenomena such as ionisation potential depression and pressure ionisation~\cite{ciricosta2012direct,hoarty2013observations, gawne2023investigating}. One advantage of XES is that it probes the filled states rather than vacant ones, therefore providing complementary information to those extracted with RIXS and XAS.

Lastly, I highlight x-ray Raman scattering (XRS) as an additional method for probing electronic structures~\cite{Voigt_PoP_20201}. This process is described by the first term in Eq. (\ref{RIXS theory 3}):

\begin{equation}
     \sigma_{\rm XRS}(\omega, \textbf{k}) \propto \sum_f |\langle f|\rho_{\textbf{k}}|i\rangle|^2 \delta(E_f-E_i-\hbar\omega) \ ,
    \label{xrs}
\end{equation}

\noindent where $\textbf{k}= \textbf{k}_2- \textbf{k}_1$ and $\omega = \omega_2 - \omega_1$ were defined. This expression encompasses various inelastic scattering processes, which are distinguished differently by the type of transitions and the energy transfers involved. For instance, XRTS is described by the same formula but pertains to smaller energy transfers than XRS and involves also collective excitations~\cite{kritcher2008ultrafast}.

The XRS process is always present and has the advantage of versatility, enabling the exploration of a wide range of energy and momentum transfers, and thus a variety of initial and final states. However, the XRS signal is extremely weak and, can often be obscured by other competing processes. By contrast, RIXS benefits from a larger cross-section due to resonance, making it more effective when the set of states under investigations can perform this type of resonant scattering.


% A figure with a XRS spectrum (fig.2 voigt paper), and an absorption spectrum (XANES+EXAFS) (fig2.6 (top) fundamental of exafs)







\section{\label{method} Method for RIXS Analysis}






Modern physics experiments increasingly depend on indirect or proxy measurements because direct access to quantities of interest is often limited. Instead, these experiments utilise known physical processes—dependent on the relevant quantities—to generate measurements that, in principle, allow for their extraction. However, these processes are frequently complex, making extraction a challenging inverse problem. The determination of $\rho_{\rm eff}$ from the RIXS spectrum, specifically through Eq.~(\ref{finite_FWHM}), is an example of such a problem. To address these challenges, a method based on neural networks designed to solve problems of this nature is now introduced. While addressing inverse problems with neural networks is not new, a variety of methods have evolved around it. Many of these methods involve fully inverting the forward process with end-to-end estimators~\cite{Jin2017DeepImaging}, which, in experimental analysis, essentially means re-learning known physics. This approach can improve speed compared to traditional numerical models, albeit often at the expense of accuracy.

A recent alternative, focused on accuracy, is emerging in differentiable modelling. Extending parameter fitting via gradient descent and backpropagation~\cite{Thuerey2022PhysicsLearning}, the use of neural networks embedded within differentiable physical models—hereafter referred to as surrogate training embedded in physics (STEP)—has gained momentum in fields such as robotics and quantum chemistry. Similar approaches based on differentiable modelling are also seen in fluid dynamics, where neural networks have been used to refine or accelerate solvers of partial differential equations (PDEs)~\cite{um2020solver}. Additionally, physics-informed neural networks (PINNs)~\cite{raissi2017physics} achieve similar results through differentiable programming, although PINNs rely on physics modelling (PDEs in their case) as additional constraints to those posed by experimental observations~\cite{karniadakis2021physics}, rather than directly backpropagating the experimental constraints through the physical model. Therefore, in PINNs physics modelling is used alongside experimental data in a parallel way, instead of a sequential way, as is the case with STEP. Likewise, most STEP applications in quantum chemistry use either extensive high-quality datasets~\cite{schoenholz2020jax}, specialised estimators~\cite{schoenholz2020jax, li2021kohn}, or both. A notable exception is the work by Kasim and Vinko~\cite{kasim2021learning}, which demonstrates that the exchange-correlation functional can be accurately learned with comparatively small datasets, a finding also supported by Li \textit{et al.}~\cite{li2021kohn}. Kasim and Vinko further show that a simple feed-forward neural network (FFNN) with three hidden layers of 32 nodes each is sufficient for this purpose. The simplicity of this neural network architecture and its ability to work with limited data are both significant for physics research, where data scarcity and the need for tailored ML architectures often pose challenges for end-to-end models to reach high accuracy.

Here, I first present a formal description of the STEP paradigm for inverse problems, demonstrating that STEP generalises well, offers interpretability, and resists overfitting—as long as the forward model is accurately defined. Then this method is applied to extract $\rho_{\rm eff}$ from RIXS spectra solving the deconvolution described by Eq. (\ref{finite_FWHM}). 



% In robotics, STEP is commonly applied to optimise reinforcement learning policies within simulations~\cite{toussaint2018differentiable, degrave2019differentiable, de2018end}, exposing neural networks to large datasets. 







\subsection{The STEP Paradigm}




Consider a known physical process $P$, which takes as input a set of $N$ functions $B_i$ with $i = 1,\hdots,N$, and returns another function $A$, i.e. $P$ is a functional which may be written as
\begin{equation}
    A(x) = P[B_1,\hdots,B_N](x).
    \label{ML 1}
\end{equation}

\noindent Now consider the case where one wants to determine one specific function $B_i$ given that one knows $A$, $P$, and $B_j$ for any $j\neq i$. Mathematically, this entails finding an inverse functional $P^{-1}$ w.r.t. $B_i$ such that
\begin{equation}
    B_i(y) = P^{-1}[A, B_1, \hdots, B_{i-1}, B_{i+1}, \hdots, B_N](y).
     \label{ML 2}
\end{equation}

\noindent To construct this inverse map, one can employ an estimator $\tilde{P}^{-1}$ to approximate $P^{-1}$, which gives us an estimate for $B_i$ as
\begin{equation}
    \tilde{B}_i(y) = \tilde{P}^{-1}[A, B_1, \hdots, B_{i-1}, B_{i+1}, \hdots, B_N](y).
    \label{ML 3}
\end{equation}

\noindent While this is a highly successful approach~\cite{stielow2021reconstruction}, particularly in image analysis~\cite{Jin2017DeepImaging}, it also relies on the approximate inversion of the potentially highly complex process $P$. Therefore, non-trivial design choices are usually made regarding the estimator, as it has to invert a potentially highly complex process. Furthermore, this method does not necessarily generalise beyond the domain of the data it was trained on, making predictions outside the range of training values for $A$ and $B_i$ non-trustworthy. Finally, many distinct data pairs $A$ and $B_i$ are also required for an estimator to learn the general inverse $P^{-1}$, data which may not be readily available.

To address these issues, the STEP method offers an alternative approach that emphasises accuracy, though at the expenses of longer computational times. Instead of approximating the complex object $P^{-1}$, STEP directly estimates the much simpler function $B_i$ as $\tilde{B}^{(STEP)}_i$, represented with a neural network. One then has

\begin{equation}
    \tilde{A}(x) = P[B_1, \hdots, \tilde{B}^{(STEP)}_i, \hdots, B_N](x),
\end{equation}

\noindent and may compute a loss $\mathcal{L}(A, \tilde{A})$ to measure the estimate's performance on the data. Provided that the process $P$ is known, mathematically differentiable, and implemented using automatic differentiation, backpropagation and gradient descent can be applied to iteratively refine $\tilde{B}^{(STEP)}_i$ toward an accurate approximation of the true $B_i$. In the physical applications, $P$ and $B_i$ can be defined as known and unknown processes, respectively, isolating $B_i$ as the physical quantity one aims to uncover. This distinction helps us focus on inferring unknown physical quantities rather than inverting well-understood processes $P$. Finally, STEP allows generalizing to arbitrary $B_i$ with confidence, and ensures that $\tilde{B}^{(STEP)}_i$ always has an intuitive and meaningful interpretation as the exact estimate of $B_i$, two important features in physical ML.

\begin{figure}[t]
     \centering
     \includegraphics[width=.9\textwidth]{Figures/step scheme.pdf}
     \caption[Scheme of the STEP paradigm]{An illustration of a supervised learning scheme under the STEP paradigm. Note the existence of the differentiable physical model, which allows us to represent the unknown learnable physical object $B_i$ as an estimator with much smaller architecture than the ones typically used for representing an inverse operator. Black arrows indicate forward passes, and red ones backpropagation.}
     \label{fig:diff_modelling}
 \end{figure}







\subsubsection{Limitations}







When using STEP, there are some important limitations to consider. The foremost requirement is the ability of the operator $P$ to accurately model the underlying physical process, and to be implemented in an automatically differentiable way. However, for measurement processes in physics, the former is in fact the underlying assumption in measurement routines. For automatic differentiation of $P$, one may leverage the plethora of tools for differentiable physical programming~\cite{schoenholz2020jax, Innes2019AComputing, Freeman2021BraxSimulation} developed in the context of supervised learning. These tools, embedded in growing frameworks like PyTorch~\cite{Paszke2019PyTorch:Library} and JAX~\cite{Bradbury2018JAX:Programs}, enable full differentiability in physical models by reducing intractable computational graphs and non-differentiable effects into feasibly differentiable forms using automatic differentiation (e.g. autograd~\cite{Paszke2019PyTorch:Library}). In practice, such frameworks allow gradient descent and backpropagation to be performed within physical processes, making it possible to apply STEP. However, care must be taken to manage issues like vanishing gradients in backpropagation (see below), which can affect optimisation.

Secondly, the use of STEP incurs a significant increase in computational time, as it relies on repeated fittings for every $B_i$, rather than the single evaluation of the approximate inverse map. This issue has become less significant with the continuing rise in available computing power, and importantly, it is most critical for real-time applications. Most physics experiments, however, record all data and rely primarily on post-experimental analysis, where accuracy supplants fitting time as the key performance metric. Additionally, the $B_i$ of interest are often constant for a given set of experiments (as in RIXS), or even for all physics (as in the case of the exchange correlation functional~\cite{kasim2021learning}), and consequently, different datasets can be treated together to extract a single $B_i$.







\subsection{\label{STEP for the RIXS Process} STEP for the RIXS Process}









% Normalisation solid angles : They can be included in the forward model




The expression for the RIXS spectra given in Eqs. (\ref{rixs cross-section 2}) and (\ref{finite_FWHM}) takes the form of a sum of convolutions with a dynamic kernel $\Phi(\omega_1)$.
The kernel is dynamic because it represents the XFEL pulse, which is commonly formed of a series of narrow spikes in photon energy that change stochastically from shot to shot. If a typical dataset is assumed to contain $N$ RIXS single shots indexed by $k$, with associated RIXS spectra $I_k(\omega_2)$, then each is given by Eq.~(\ref{finite_FWHM}) using the corresponding XFEL spectra $\Phi_k(\omega_1)$. The following primarily focuses on XFELs operating in SASE mode to demonstrate how this method enables the extraction of electronic structures with resolutions surpassing the probe's bandwidth. However, this approach is general and can be applied to XFEL spectra of any shape, such as those obtained in seeded mode~\cite{amann2012demonstration}. 

The relation that links the $k$ measured spectra $I_k$ to the effective density of states $\rho_{\rm eff}$ and the spectral FEL kernel $\Phi_k$, given by Eqs. (\ref{rixs cross-section 2} and (\ref{finite_FWHM}), can be formally written with the notation introduced in the previous section, defining an operator $P_{\rm RIXS}$ as:

\begin{equation}
    I_k(\omega_2) = P_{\rm RIXS}[\Phi_k, \rho_{\rm eff}](\omega_2),
\label{Eq:inductiveBias}
\end{equation}

\noindent which is an example of the general problem \ref{ML 1}. As already discussed, one possible solution to find the effective DOS is by inverting the RIXS operator $P_{\rm RIXS}$ with respect to $\rho_{\rm eff}$:
\begin{equation}
    \rho_{\rm eff}(\epsilon) = P^{-1}_{\rm RIXS}[I_k, \Phi_k](\epsilon),
\end{equation}

\noindent which is valid for each $k$.

\begin{figure*}
     \centering
     \includegraphics[width=0.8\textwidth]{Figures/RIXS STEP scheme.pdf}
     \caption[Scheme of the method for the extraction of electronic structures from RIXS spectra]{Representation of the scheme used to extract physical information on the electronic structure of a material from a SASE-based RIXS measurement. The estimate of $\rho_{\rm eff}$, represented by a feedforward neural network, is fed to the RIXS differentiable model alongside the measured XFEL spectrum. The model outputs the calculated RIXS spectrum, which is then compared with the experimental spectrum. This comparison generates a loss function used to adjust the neural network parameters, thereby refining the $\rho_{\rm eff}$ estimate. The resolution of the resulting density of states depends on the resolutions of the spectrometers used to measure the RIXS signal and the profile of the XFEL pulse, but is independent of the XFEL pulse structure or bandwidth.}
     \label{Fig:STEP}
 \end{figure*}

However, since the RIXS cross-section is relatively small, so that single-shot measurements typically have low SNR, the use of standard inversion techniques via deconvolution, such as the Richardson-Lucy method, is unsuitable for single-shot analysis. Integrating over many shots ($I_k \rightarrow \bar{I}$ and $\Phi_k \rightarrow \bar{\Phi}$) can mitigate such limitations in SNR, but, as discussed shortly, it also limits the resolution with which one can extract $\rho_{\rm eff}$. The lack of a robust approach to process low SNR data actually represents a considerable bottleneck for x-ray spectroscopy in high energy density physics applications~\cite{kayser2019core}. Alternatively, one could use the ML method introduced in the previous section and represent $P_{\rm RIXS}^{-1}$ with a neural network. Nevertheless, given the stochastic nature of the FEL pulse profile and the complexity of the RIXS operator, the construction of a NN estimator to approximate $P_{\rm RIXS}^{-1}$ from a large labelled dataset of known pairings $(\rho_{\rm eff}, I_k, \Phi_k)$, is not effective. This, in fact, would require a highly complex estimator and the collection and validation of a sufficiently large dataset, which is not a straightforward task.

Therefore, rather than searching for a general estimator to approximate the highly complex object $P_{\rm RIXS}^{-1}$, one can instead look for a suitable approximation $\tilde{\rho}_{\rm eff}$ to the much simpler function $\rho_{\rm eff}$, using the STEP method. A comparison of the performances, for the RIXS problem, between the STEP method and an end-to-end inversion estimator is presented in Ref.~\cite{alaa2024step}.
For the RIXS case, the scheme of the STEP method is illustrated in Fig.~\ref{Fig:STEP}. The known forward model $P_{\rm RIXS}$ is used to calculate the predicted spectral intensity given a known $\Phi_k$ and the estimate $\tilde{\rho}_{\rm eff}$:
\begin{equation}
    \tilde{I}_k(\omega_2) = P_{\rm RIXS}[\Phi_k, \tilde{\rho}_{\rm eff}](\omega_2).
\end{equation}

\noindent Considering that, for this cases with low photon rates, the uncertainties of the single shots are approximately constant over $\omega_2$, the spectral intensity $\tilde{I}_k$ is then used to compute an L$_2$ loss function 
\begin{equation}
    \mathcal{L}_k = || \tilde{I}_k(\omega_2) - I_k(\omega_2)  ||_2,
\end{equation}

\noindent which provides a measure of the quality of the approximate $\tilde{\rho}_{\rm eff}$ given the observed $I_k(\omega_2)$ and $\Phi_k(\omega_1)$. More generally, for uncertainties varying significantly across the spectrum, a weighted L$_2$ norm should be employed, where the weights for each point $\omega_{2,i}$ are taken inversely proportional to the uncertainty on $I(\omega_{2,i})$, which can be determined with the methods outlined in Chapter 2. Using the STEP methodology, $\tilde{\rho}_{\rm eff}$ is represented using a trainable feed-forward neural network acting as a universal approximator, and $P_{\rm RIXS}$ is implemented in an automatically differentiable form~\cite{schoenholz2020jax}. Improving the approximation for $\tilde{\rho}_{\rm eff}$ can now be viewed as a standard machine learning optimisation problem. Thanks to the differentiability of the forward model, the gradients of this loss function can be propagated back to the NN parameters, which are therefore adjusted through a gradient-descent optimisation, in order for the forward model predictions to match the experimental measurements. In this framework, the learning process is essentially just a fitting of the experimental RIXS spectra through the reconstruction of $\rho_{\rm eff}$.

Since convolution is an injective operation, there is a one-to-one correspondence between $\rho_{\rm eff}$ and the corresponding RIXS spectrum. This property ensures that, in absence of noise, the NN reconstructs the correct $\rho_{\rm eff}$ as long as the forward model describes accurately the physical process. In reality, due mainly to the signal suppression away from the $K_{\alpha}$ lines, caused by the Lorentzian factor in Eq. (\ref{rixs cross-section 2}), $P_{\rm RIXS}$ is an operator that maps very different $\rho_{\rm eff}$ into very similar RIXS spectra (contraction map). This fact, along with the noise in $I_k$, makes the solution to the inverse problem and the relative estimate of $\rho_{\rm eff}$ very challenging.

The energy resolution with which $\tilde{\rho}_{\rm eff}$ can be extracted is given by the SNR of the single RIXS spectra and the energy resolutions with which $I_k$ and $\Phi_k$ are measured. While the method can, in principle, be used to optimise single shot data, it is more effective to perform batching of the data to improve the SNR across multiple shots (see next section for more information). In contrast to standard averaging approaches, the resolution is not degraded by such a merging of multiple shots, since each pairing ($I_k$, $\Phi_k$) is still considered individually in the optimisation process.
















\subsubsection{Details of the STEP Method Implementation}










The feed-forward NN used to reconstruct $\rho_{\rm eff}$ has a single input ($\epsilon$), a single output ($\rho_{\rm eff}(\epsilon)$), 4 hidden layers with 40 nodes each, and the \textit{softplus} activation function~\cite{zheng2015improving}. These architectural choices were selected based on extensive testing and are used throughout this thesis. However, alternative architectures may be more suitable for specific problems. For example, achieving a higher resolution for $\rho_{\rm eff}$ might require increasing the number of hidden layers or nodes. The neural network is trained by means of the ADAM~\cite{Kingma2014Adam:Optimization} optimiser with the initial learning rate set to $10^{-3}$. This optimisation, together with the automatic differentiation of the forward model, is carried out using the library PyTorch~\cite{ketkar2021introduction}. To ensure a good performance of this scheme, it is necessary to address the vanishing gradient problem~\cite{hochreiter1998vanishing} by omitting exponentially small factors in the gradient backpropagation. The training is generally performed splitting the set of experimental spectra ($I_k$, $\Phi_k$ for $k=0,1,\dots,N$) into batches, whose size is increased dynamically during the training to shift from exploratory to exploitative training. The optimisation routine considers the loss function, and the respective gradients, on an entire batch to update model parameters at each step. Furthermore, the training is carried out over many epochs, i.e. rolling over all the N shots multiple times. Finally, during the training, the output of the neural network $\tilde{\rho}_{\rm eff}(\varepsilon)$ is constrained to be positive, a fundamental physical requirement. 




















\subsection{Validation on Synthetic Data}



% Studies for different NN architectures!!! (written above)




\begin{figure}[t]
    \centering
    \includegraphics[width=\textwidth]{Figures/Synthetic dataset.pdf}
    \caption[An illustration of the RIXS synthetics spectra]{An illustration of the RIXS synthetics spectra. The real XFEL pulses 1 and 2 (spectra shown in top row) are used as incoming radiation into a simulated RIXS processes in materials with DAS A and B (left column). The 4 panels in the bottom right show the simulated scattering spectra before (blue) and after (magenta) noise is applied.}
    \label{fig:artificial_data}
\end{figure}



To show the efficacy of the method, it was tested on synthetic data, generated using the following procedure:

\begin{enumerate}
    \setlength{\itemsep}{10pt}
    \item Create an artificial DAS ($\rho|M|^2$) designed to resemble a real one (Gaussian energy bands and an optional square-root continuum);
    \item Compute the relative $\rho_{\rm eff}$ choosing a temperature $T$ and an electron number density $n_e$. The chemical potential $\mu$ is then calculated from the following implicit equation

    \begin{equation}
        n_e = \int_0^{\infty} \rho(\epsilon)f_{FD}(\epsilon;T,\mu) \ .
        \label{chemical potential}
    \end{equation}
    
    \noindent For this data, I chose $T=5$ eV and $n_e = 16$ atoms/unit cell (which is the iron electron density);
    \item Evaluate $I_k(\omega_{2})$ using Eqs. (\ref{rixs cross-section 2}) and (\ref{finite_FWHM}) with real XFEL SASE pulses as $\Phi_k(\omega_1)$. For the physical parameters of the forward model $P_{\rm RIXS}$, appearing in Eq. (\ref{rixs cross-section 2}, I chose those specific of iron. These 
    parameters have been taken from~\cite{henke1982low} and are shown in table \ref{tab:params};
    \item Add Gaussian noise to the obtained spectra $I_k(\omega_{2})$. The noise is distributed with a standard deviation of $\sigma = \alpha\cdot \max\left(I(\omega_{2})\right)$.
\end{enumerate}



\begin{figure*}[t]
     \centering
     \includegraphics[width=\textwidth]{Figures/Synthetic.pdf}
     \caption[Reconstruction of a synthetic density of states from a RIXS spectrum at varying levels of noise]{Reconstruction of a synthetic density of states from a RIXS spectrum at varying levels of noise, using a realistic SASE XFEL profile. A total of 50 single shots were used. Examples of the synthetic RIXS spectra are shown in panels (a-c), with the corresponding reconstructions in panels (d-f). The technique is robust to noise, owning to the strong inductive bias given by the model, and the stochastic nature of the SASE pulse profile $\Phi$ that allows us to oversample the DOS. The zero of the DOS energies is centred at the resonance. Note that the original $\rho_{\rm eff}$ corresponds to the DAS B of Fig. \ref{fig:artificial_data} modulated by the Fermi-Dirac distribution.}
     \label{Fig:synthetic}
 \end{figure*}

\noindent I chose six spiky synthetic DAS, which are both challenging to tackle for traditional deconvolution algorithms, but also representative of narrow bound states and $d-$band features. Examples of resultant spectra for different combinations of DAS and XFEL pulses are shown in Fig. \ref{fig:artificial_data}. Note that with the noise applied, the features in the spectrum are barely recognisable. These synthetic RIXS spectra are then provided to the deconvolution scheme, alongside the corresponding SASE spectra, to extract the $\rho_{\rm eff}$ corresponding to a certain DAS (here $T$ is assumed to be unknown as a more realistic case). Fig \ref{Fig:synthetic} shows the results for DAS B of Fig. \ref{fig:artificial_data} under three different levels of noise analysed, considering that analogous behaviours are observed for all the other artificial DAS. 

% \bigskip

\begin{table}[h!]
    \renewcommand{\arraystretch}{1.2}
    \centering
    \begin{tabular}{|M{30mm}|M{20mm}|M{20mm}|}
        \hline
        Parameter[units] & $f$ & Value \\\hline\hline
        $\epsilon_K[eV]$ & 1, 2 & -7112\\\hline
        \multirow{2}{*}{$\epsilon_{L,f}[eV]$} & 1 & -723 \\
        & 2 & -708.5 \\\hline
        \multirow{2}{*}{$A_f$} & 1 & 50 \\
        & 2 & 100 \\\hline
        \multirow{2}{*}{$\Gamma_f[eV]$} & 1 & 3 \\
        & 2 & 2.55 \\
        \hline
        
    \end{tabular}
    \vspace{20pt}
    \caption{Physical parameters of iron, used for the synthetic spectra. Iron has two L-shell states, implying two possible final states $f=1,2$.}
    \label{tab:params}
\end{table}

This approach shows good convergence of the extracted DOS for all levels of noise up to $\alpha=30$\%. However, the accuracy of reconstruction decreases as the level of noise is raised, with the extracted $\rho_{\rm eff}$ starting to miss small or narrow features, especially for large detuning energies. This behaviour is expected, since the area under a DOS feature and its detuning, i.e., its distance from resonance, are the two main parameters that determine the intensity of the feature in the RIXS spectrum. Bright features closer in energy to the resonance are thus more robust to low SNR. Another common factor that deteriorates the accuracy of the reconstructions is the presence of spurious oscillations in the extracted $\rho_{\rm eff}$. These oscillations are a typical byproduct of deconvolution, and are present on energy-scales smaller than the kernel's bandwidth. Regardless, the successful extraction of features using this approach even from spectra dominated by noise, as shown in Fig. \ref{Fig:synthetic}c, remains significant.

Figure \ref{loss function} illustrates the loss function, defined as the L$_2$ distance between $\rho_{\rm eff}$ and $\tilde{\rho}_{\rm eff}$, during the neural network training for one of the analysed syntethic DAS without noise applied. The figure also shows the evolution of $\tilde{\rho}_{\rm eff}$ until its final comparison with the original $\rho_{\rm eff}$. Convergence is achieved after approximately 6,000 epochs, with an initial steep decrease in the loss function corresponding to the reconstruction of the coarse features of $\rho_{\rm eff}$. This is followed by a gradual decline, reflecting the refinement of finer details. Notably, the neural network encounters difficulty in increasing the number of inflection points, resulting in plateaus in the loss function. Similar trends are observed across all the synthetic DAS analysed without noise in the spectra. When noise is introduced, the loss function exhibits greater fluctuations and eventually converges to a higher asymptotic value.

 \begin{figure}[t]
     \centering
     \includegraphics[width=\textwidth]{Figures/loss evolution.pdf}
     \caption[Example of a loss function during the STEP training]{Loss function, defined as the L$_2$ distance between $\rho_{\rm eff}$ and $\tilde{\rho}_{\rm eff}$, during the neural network training for one of the analysed syntethic DAS (distinct from that of Fig. \ref{Fig:synthetic}) with no noise applied. The plot includes also the evolution of $\tilde{\rho}_{\rm eff}$ until its final comparison with the original $\rho_{\rm eff}$. Here, one epoch is completed when each $\Phi_k$ (50 in this synthetic dataset) is used for the calculation of the loss $||\tilde{I}_k-I_k||_2$.}
    \label{Fig:loss function}
 \end{figure}
% Maybe this paragraph can be removed
A note is in order here: in this thesis I denote as RIXS scattering processes with detuning energies up to 150 eV, following the notation adopted in~\cite{humphries2020probing}, however, other names (e.g. high energy resolution off-resonance spectroscopy) have also been used to describe such scattering processes in the literature~\cite{blachucki2014high,szlachetko2014communication}. 
Unlike in Ref.~\cite{humphries2020probing} where the incident photon energies are varied, in this work (both for experimental and synthetic data) the XFEL central energy is fixed and set at about 10s of eV below the material's absorption edge to optimise the extraction of the valence density of states. 

 \begin{figure}[t]
     \centering
     \includegraphics[width=0.8\textwidth]{Figures/numberofshots_studies_combined.pdf}
     \caption[Quality of the reconstructed $\rho_{\rm eff}$ as a function of the number of RIXS shots for different methods]{Quality of the reconstructed $\rho_{\rm eff}$ as a function of the number of shots N included in the inversion dataset. The DOS used is that shown in Fig.~\ref{Fig:synthetic}a. While the direct inversion works best for noiseless data, the quality rapidly deteriorates for real data with noise. In contrast, the STEP approach is fairly insensitive to noise. The loss function employed for this plot is a weighted L$_2$ distance between the reconstructed and the original DOS, with bigger weights for the energy regions where features are present. This is achieved defining the weighted L$_2$ norm as $1/F\sum_i|d\rho_{\rm eff}/d\epsilon|(\epsilon_i)(\tilde{\rho}_{\rm eff}(\epsilon_i)-\rho_{\rm eff}(\epsilon_i))^2$, where $\rho_{\rm eff}$ is the original DOS and $F$ is a normalisation factor.}
    \label{Fig:numberofshots_studies}
 \end{figure}



To demonstrate the efficacy of the STEP method, in Fig.~\ref{Fig:numberofshots_studies} the quality of the $\rho_{\rm eff}$ reconstruction is shown as a function of the number of shots in the dataset for four different approaches and two noise levels.
The blue and red curves employ the paradigm described in Fig.~\ref{Fig:STEP} to reconstruct $\rho_{\rm eff}$: the results using all individual experimental pairs ($I_k$,$\Phi_k$) are depicted with blue circles, while the results using average values ($\bar{I}$,$\bar{\Phi}$) are shown with red squares. Similarly to what has been done in~\cite{fuller2021resonant,kayser2019core}, Eq.~(\ref{Eq:inductiveBias}) is also directly inverted via standard numerical techniques for comparison (see next section for details), both for the individual pairs (diamond markers) and for the averaged spectra (triangle markers). It is observed that while the direct inversion method achieves the best results in the absence of noise, its performance degrades rapidly as noise is added. This is due to the characteristics of $P_{{\rm RIXS}}$, which maps very different $\rho_{\rm eff}$ into similar RIXS spectra, making its inversion an ill-conditioned problem. On the other hand, the STEP approach shows little variation on the quality of the reconstructed DOS with noise, and is thus more robust. In particular, it is observed that with the STEP method, the results obtained for a 15$\%$ noise level converge to those achieved without noise already for a relatively restrained N$\approx$50.
The comparison between single shot and averaged analysis indicates that averaged calculations suffice for broad features, but that shot-by-shot analysis is required to capture finer structure, and that the averaging process leads to a loss in resolution. This can be explained considering that deconvolution is more efficient with spiky XFEL spectra, characteristic of the single shot case, rather than with relatively broad Gaussian XFEL spectra, which occur in averaged calculations\footnote{To understand this, it is better to think about the deconvolution operation in the Fourier space, where it takes the form of a division by the Fourier transform of the kernel function. A spiky kernel function has a transform with a larger support, reducing the division instabilities at high frequencies.}. Specifically, in the absence of noise, where no gain on the SNR is obtained through averaging, the use of the averaged spectra deteriorates the quality of the reconstructions as N grows.













\subsubsection{\label{sec:numerical approach}Details of the Numerical Approach}





To extract $\rho_{\rm eff}$ from a set of RIXS and XFEL spectra without making use of machine learning techniques, the following numerical routine is employed:

\begin{enumerate}
    \item One first discretises the RIXS operator for each FEL spectrum $\Phi_k$, computing the associated matrices $M_k$ for $k=1,2,\dots,N$. Each $M_k$ is constructed taking as its i-th column $P_{\rm RIXS}[\Phi_k,\rho_i]$, where $(\rho_i)_j = \delta_{i,j}$, with $i,j=1,2,\dots,L$ and $L$ the length of $\rho_{\rm eff}$. Note that this discretisation is possible because $P_{\rm RIXS}$ is a linear operator in $\rho_{\rm eff}$;
    \item One then constructs the macro-linear system of equations in $\rho_{\rm eff}$ by stacking the matrices $M_k$ and the RIXS spectra $I_k$ as

    \begin{equation}
    \begin{bmatrix}
    M_1 \\
    M_2 \\
    \vdots \\
    M_N
    \end{bmatrix}
    \rho_{\rm eff} \vcentcolon = M \rho_{\rm eff} = 
    \begin{bmatrix}
    I_1 \\
    I_2 \\
    \vdots \\
    I_N
    \end{bmatrix}
    \vcentcolon = I
        \label{linear system}
    \end{equation}

    \noindent notice that the matrix $M$ is in general not square; 
    \item Next, one approximates the solution of this linear system, and hence $\rho_{\rm eff}$, using the conjugate gradient descent method~\cite{luenberger1984linear};
    \item Finally, one applies a smoothing filter to $\rho_{\rm eff}$ to mitigate the spurious oscillation. One can show that this procedure is equivalent to perform a Tikhonov regularisation~\cite{golub1999tikhonov} with a cost function penalizing too fast oscillations in the solution.
\end{enumerate}

\noindent This idea of discretizing the physical process operator and inverting the associated matrix has already been employed in the field of spectroscopic analysis~\cite{kayser2019core}. An alternative approach to addressing this inversion problem is the maximum entropy method~\cite{rietsch1976maximum}. When working with the averaged spectra, one just searchs for a solution of the linear system $\bar{M}\rho_{\rm eff}=\bar{I}$, where $\bar{M}$ is constructed using the averaged FEL spectrum $\bar{\Phi}$.






\subsection{Extracting the Temperature}





\begin{figure}[t]
    \centering
    \includegraphics[width=.9\textwidth]{Figures/temperature extraction.pdf}
    \caption[Extraction of temperatures from the RIXS spectra with the STEP method]{Temperatures extracted from RIXS measurements plotted against true temperatures, taken from Ref~\cite{alaa2024step}. Panels (a) through (d) correspond to noise levels of 0\%, 10\%, 20\%, and 30\% respectively. Note that all values are plotted with errorbars, most of which are too small to be seen.}
    \label{fig:temp_fits}
\end{figure}


As previously discussed, the STEP paradigm can also be used to extract temperature from RIXS spectra, provided the underlying DAS is known from other sources (e.g. DFT calculations)\footnote{While this assumption - especially at high pressures -is somewhat idealised (it actually represents the central challenge addressed in this thesis), reasonable approximations of the DAS are available and may be refined in the future.}. In this scenario, the extracted $B_i$ (refer to Eq.~\ref{ML 1}) becomes a scalar rather than a function, while the method itself remains unchanged. To test the extraction of 
$T$ from RIXS spectra using the STEP method, the same procedure as before is followed, but incorporating the known DAS into the forward model. Consequently, the backward propagation of gradients must go through additional operations, including the division of $\rho_{\rm eff}$ by the DAS and the extraction of $T$ from the Fermi-Dirac distribution. Furthermore, the chemical potential must be recalculated with Eq. \ref{chemical potential} after each update of $T$.

This method was applied to the previously introduced set of spectra, but with seven different values of temperatures ranging from 4 to 10 eV. The results, shown in Figure~\ref{fig:temp_fits}, demonstrate excellent temperature predictions that remain unaffected by noise. For each noise level, the extracted temperatures and their uncertainties were calculated as the mean and standard deviation of the values obtained across all the different DAS. On average, convergence was achieved after under 2,000 epochs, depending on the initial random value of temperature. These results highlight the successful functioning of the STEP method for scalars and its robustness against noise.








\section{Summary}




This chapter revisited the theory of RIXS, tailoring it to our specific cases of interest and demonstrating how it can, in principle, provide information about the density of states and the temperature of the material being probed. The RIXS technique was also compared with other spectroscopic methods, highlighting their relative advantages and disadvantages. I then introduced a method, called STEP, which combines neural networks with physical modeling to address inverse problems where a physical quantity must be extracted from experimental measurements in which it is embedded. This method was then specified for the RIXS problem to extract the electronic structures and temperatures from the RIXS spectra. Finally, the application of this method to synthetic data was shown, testing its performances and comparing them with numerical techniques. The key finding of this study was the STEP method's improved robustness to noise compared to numerical approaches. Potential enhancements of this approach and ideas for other physical problems where the STEP method could be employed are given in Chapter 6.




