# Relativistic EELS scattering cross-sections for microanalysis based on Dirac
solutions

Zezhong Zhang<EMAIL_ADDRESS>Ivan Lobato  Hamish Brown  Dirk
Lamoen  Daen Jannis  Johan Verbeeck  Sandra Van Aert
<EMAIL_ADDRESS>Peter D. Nellist
<EMAIL_ADDRESS>Electron Microscopy for Materials Research
(EMAT), University of Antwerp, Groenenborgerlaan 171, 2020 Antwerp, Belgium
NANOlab Center of Excellence, University of Antwerp, Groenenborgerlaan 171,
2020 Antwerp, Belgium  Department of Materials, University of Oxford, 16 Parks
Road, Oxford OX1 3PH, United Kingdom  Rosalind Franklin Institute, Harwell
Research Campus, Oxfordshire OX11 0FA, United Kingdom  Ian Holmes Imaging
Centre, Bio21 Molecular Science and Biotechnology Institute, University of
Melbourne, Victoria, Australia

###### Abstract

The rich information of electron energy-loss spectroscopy (EELS) comes from
the complex inelastic scattering process whereby fast electrons transfer
energy and momentum to atoms, exciting bound electrons from their ground
states to higher unoccupied states. To quantify EELS, the common practice is
to compare the cross-sections integrated within an energy window or fit the
observed spectrum with theoretical differential cross-sections calculated from
a generalized oscillator strength (GOS) database with experimental parameters.
The previous Hartree-Fock-based and DFT-based GOS are calculated from
Schrödinger’s solution of atomic orbitals, which does not include the full
relativistic effects. Here, we attempt to go beyond the limitations of the
Schrödinger solution in the GOS tabulation by including the full relativistic
effects using the Dirac equation within the local density approximation, which
is particularly important for core-shell electrons of heavy elements with
strong spin-orbit coupling. This has been done for all elements in the
periodic table (up to Z = 118) for all possible excitation edges using modern
computing capabilities and parallelization algorithms. The relativistic
effects of fast incoming electrons were included to calculate cross-sections
that are specific to the acceleration voltage. We make these tabulated GOS
available under an open-source license to the benefit of both academic users
as well as allowing integration into commercial solutions.

###### keywords:

Electron energy loss spectroscopy (EELS) , Generalized oscillator strength
(GOS) , Ionization , Inelastic electron scattering

††journal: Ultramicroscopy

##  1 Introduction

Electron energy-loss spectroscopy (EELS) involves analyzing the energy loss
distribution of electrons after their interaction with a specimen. With the
underlying inelastic electron scattering event, the fast-moving incident
electron transfers part of its energy and momentum by exciting the atomic
system from its initial quantum ground state to various excited final states
[1]. Since these transitions are characteristic of the type of elements and
their local electronic structure, EELS can fingerprint the composition and
bonding of the specimen. Thanks to advances in instrumentation such as
aberration correction [2], monochromation [3, 4] and more recently direct
electron cameras [5, 6], EELS can now be performed with high energy resolution
(up to meV) and spatial resolution (up to Å). Recognized as a powerful
analytical tool for site-specific materials characterization [7, 8], there is
a general interest in the accurate and precise quantification of EELS and this
therefore necessitates a reliable theoretical model to calculate the inelastic
scattering cross-sections, which describes the probability of inelastic
scattering events to take place.

To quantify EELS, one has to model the spectra as a superposition of the
background signal and multiple excitation edges with characteristic shapes.
The traditional method involves background removal using a power-law
extrapolation from the pre-edge region and integration of the edge signal
within a certain energy window [9]. For the analysis of fine structures,
deconvolution with the low-loss spectrum is performed to remove plural
scattering (known as the Fourier-ratio method)[10]. However, this methodology
can result in various artifacts during deconvolution and can be very sensitive
to the choice of suitable windows for background removal and signal
integration [11]. A better approach is to perform a model-based fitting of the
observed spectrum with theoretical scattering cross-sections. This has clear
benefits: (a) effective separation of overlapping components, (b)
consideration of multiple scattering through convolution of the core-loss
model with the low-loss spectra instead of deconvolution, and thereby its
artifacts, (c) most significantly, the model-based parameter estimation theory
grants statistical superiority, resulting in improved accuracy and free from
bias [12, 11]. The model-based method is implemented in EELSModel [11] (and
its python version pyEELSModel [13]), HyperSpy [14] (the EDX/EELS analysis is
recently renamed exspy) and Gatan Microscopy Suite SoftwareTM, which are
widely used in the EELS community. Continuous efforts have been devoted to
improving the model-based EELS, including developing a better background model
[15, 16], addressing the correlated noise effects from the camera [17] and
incorporating a phenomenological model for the fine structure during signal
processing [18]. However, the day-to-day scattering cross-section computation
has long relied on the Hartree-Fock (HF) based generalized oscillator strength
(GOS) database, which involves some limitations: (1) the HF solutions are not
fully relativistic, which can be critical for heavy elements with strong spin-
orbit coupling (SOC); (2) not all (though a good coverage) elements and edges
are included in the database; (3) the energy-momentum space is limited and
sampling is not fine enough for accurate cross-section calculations; (4) the
database is under a commercial license limiting its use in community-driven
development.

Theoretical core-shell cross-sections are computed based on a set of
experimental parameters (such as acceleration voltage, convergence angle, EELS
collection angle, etc) and a GOS database, which couples the atomic physics
with electron scattering. The name "generalized oscillator strength" comes
from the extension of the "oscillator strength", which describes the
probability of absorption or emission of electromagnetic radiation caused by
the electron orbital transitions for an atom or molecule [19]. Indeed the
precise analysis of atomic light spectra, tracing back to the starting point
of quantum mechanics, lays the foundation of atomic physics. The transition
probability can be calculated using Fermi’s golden rule, which states that the
transition probability is proportional to the strength of the coupling between
the initial and final states under weak perturbation (first Born
approximation). Interestingly, the idea was first formulated by Dirac in 1927
[20], just a year before the publication of the Dirac equation [21]. On its
importance for quantum transitions, Fermi referred to it as the "golden rule
No. 2" in a university nuclear physics class [22] and thereafter was known as
the "Fermi’s golden rule". In the context of core-shell electrons, the initial
state represents the occupied bound state of a given orbital, while the final
states can include unoccupied bound orbitals above the Fermi level or
continuum states above the vacuum level after ionization. The local density of
unoccupied states highly depends on the valence electron redistribution in the
crystal field (known as bonding), leading to complicated energy loss near-edge
fine structures (ELNES). In addition, the reflection of outgoing excited
electron waves by surrounding atoms will generate the extended energy loss
fine structures (EXELFS). In contrast, at a moderate higher energy loss
(typically tens of eV above the edge to be free from ELNES and EXELFS), the
perturbation caused by the local environment is much decayed, and the
excitation to the continuum states defines the general spectra shape for each
edge [10]. Thus, the GOS can be computed for a single isolated atom with final
states limited to the continuum state and applied to any material system for
the quantification of such an element.

The computation of the theoretical EELS cross-sections dates back to 1932 when
Bethe first used the Born approximation to consider the angular-dependent
probability of a collision of an electron with matter at relativistic speed
using hydrogenic solutions[23]. In the 1950-60s, ionization cross-section
studies still often relied on the hydrogenic solutions [24], which describe
the intermediate and large energy loss fairly well for K-shells but cannot
capture the near-edge behaviors for higher shells. Since the 1970s,
improvements were made by McGuire [25], Manson [26] and Scofield [27], who
used HF solutions of atomic orbitals to calculate K and L-edges cross-sections
of light elements. These HF-GOS efforts were systematically extended to a much
broader range of elements and shells by Leapman, Rez and co-workers [28, 29,
30], later becoming the basis for the Gatan GOS database. Many studies in the
atomic physics community focused on more accurate computational scattering
methods going beyond the plane wave Born approximation [31, 32]. For instance,
as the name suggests, the plane wave Born approximation assumes the incoming
and outgoing electron to be in the form of a plane wave, which is reliable for
the fast electron with kinetic energy well above the ionization threshold, but
becomes a problem when the electron is slow and noticeably distorted by the
presence of the target atom. In this case, the distorted wave Born
approximation method is developed to take the finite influence of atomic
potential into account. Meanwhile, much of the attention in the electron
microscopy community was devoted to dynamical scattering to explain the
complex phenomenon in energy-filtered imaging[33, 34], diffraction [35], EDS
[36] and EELS [37, 38, 39, 40, 41, 42]. In particular, Yoshioka [43]
demonstrated that the inelastic scattering and the associated absorption can
be treated as an additional transition potential with an imaginary part,
forming the basis for the implementations of dynamic inelastic scattering in
various algorithms, including Bloch wave [44, 45, 46], multislice [33, 47, 48,
49, 50] and more recently the PRISM algorithm [51]. The fast electron used in
mid/high voltage electron microscopes has considerable relativistic effects,
which boosts the cross-sections in an anisotropic way across the energy-
momentum space. Such relativistic electrodynamics behaviors of the incoming
electron beam on the differential scattering cross-sections was first studied
by Fano in 1956 [52] and have been later considered in several studies [50,
53, 54, 55, 56] but with orbital wavefunctions based on the Schrödinger
solution. In addition, significant interest has emerged in the density
functional theory (DFT) community to calculate the fine structure with the
density of states explicitly included for EELS [57] and X-ray absorption
spectroscopy (XAS) [58], all based on the Schrödinger picture. Given our focus
is on the absolute cross-sections for EELS microanalysis, the dynamic
scattering effects, and fine structures are beyond the scope of this paper. We
noticed a clear pattern in the literature that electron microscopy research
typically adheres to the Schrödinger framework, whereas the atomic physics
community has long utilized the Dirac framework. Therefore, this paper also
aims to bridge the conventions, recognizing that both fields are engaged in
the same endeavor.

To go beyond the limitation of the current GOS databases for accurate EELS
quantification, here we report the development of an open-source GOS database
using the Dirac equation. We want to note that this work at its inception is
in parallel with the recent Schrödinger GOS [59, 60], which solved the radial
wavefunction with DFT generalized norm-conserving pseudopotential potential
[61] (hereafter referred to as the DFT-based GOS). Once we learned about this
parallel development, we agreed to adopt the same file standard as GOSH [62]
for the convenience of the user. Different from existing Schrödinger-based GOS
databases (i.e. HF and DFT), we chose the Dirac equation for its intrinsic
incorporation of relativity. This is particularly crucial when the electron
approaches the speed of light, especially in heavy elements. In this paper, we
will introduce the theoretical background of atomic orbital calculations and
inelastic electron scattering. First, we will take the reader through the
inelastic scattering from the Schrödinger equation in Sec. 2, followed by a
discussion of the relativistic Dirac equation in Sec. 3 to highlight the
connections and differences. In the relativistic inelastic scattering theory,
we will consider the relativistic effects on the cross-sections due to fast-
moving incident electrons. In Sec. 4, we will explain the inclusion of
experimental parameters of EELS geometry (i.e. convergence angle and
collection angle) in the cross-section calculations. Then, we will present the
computational details for GOS database generation in Sec. 5 followed by the
results section in Sec. 6. Specifically, Sec. 6.1 will give an overview of the
Dirac-based GOS in terms of the energy-momentum sampling compared to existing
databases. Then we will explicitly explore the effects of SOC and the
contribution of small components which is unique in the Dirac solution. In
Sec. 6.2, we will demonstrate the relativistic electrodynamics effects of the
incoming fast electron on the cross-sections for typical acceleration voltages
and EELS collection angles. In Sec. 6.3, we will compare the Dirac-based GOS
with the HF-based GOS and the DFT-based GOS for different elemental edges.
Sec. 6.4 will briefly demonstrate how to use the Dirac-based GOS database for
EELS quantification. Finally, we will conclude the paper in Sec. 7 with a
summary and outlook for future development.

##  2 Theory: non-relativistic atomic orbitals and inelastic electron
scattering

In this section, we provide an overview of the Schrödinger equation for atomic
orbitals and inelastic electron scattering. For a comprehensive treatment of
this topic, readers are referred to Cowan’s book on atomic spectroscopy [63],
here we only briefly outline the key equations for calculations.

###  2.1 Schrödinger equation for atomic orbitals

In non-relativistic quantum mechanics, the Schrödinger Hamiltonian
$\hat{H}_{S}$ consists of the kinetic energy operator ($\hat{T}$) and the
potential energy operator ($\hat{V}$) for an electron in the presence of a
radially symmetric central potential of an atom $V(r)$ in spherical
coordinates:

$\begin{split}\hat{H}_{S}&=\hat{T}+\hat{V}=\frac{\hat{p}^{2}}{2m_{e}}-eV(r)\\\
&=-\frac{\hbar^{2}}{2m_{e}}\frac{1}{r}\frac{\partial}{\partial
r^{2}}r+\frac{\hat{L}^{2}}{2m_{e}r^{2}}-eV(r)\end{split}$ (1)

In this equation, $\hbar$ is the reduced Planck’s constant, $m_{e}$ is the
mass of the electron, linear momentum operator is given by
$\hat{p}=-i\hbar\vec{\nabla}$, where $\vec{\nabla}$ is the gradient operator
for the wavefunction in space. For a central potential in spherical
coordinates, the linear momentum can be further decomposed into a radial part
and an angular part. The angular momentum operator $\hat{L}$ is defined as:

$\hat{L^{2}}=-\hbar^{2}\left[\frac{1}{\sin\theta}\frac{\partial}{\partial\theta}\left(\sin\theta\frac{\partial}{\partial\theta}\right)+\frac{1}{\sin^{2}\theta}\frac{\partial^{2}}{\partial\phi^{2}}\right],$
(2)

with the $\theta$ and $\phi$ being the polar and azimuthal angles,
respectively. The projected angular momentum operator $\hat{L}_{z}$ along the
z-axis is given by:

$\hat{L_{z}}=-i\hbar\frac{\partial}{\partial\phi}.$ (3)

The angular momentum operator $\hat{L}=\vec{r}\times\hat{p}$ commutes with the
Hamiltonian, with the eigenvalues of $L^{2}=l(l+1)\hbar^{2}$ and
$L_{z}=m\hbar$ respectively. Here we follow the convention that $n$ is the
principal quantum number, $l$ is the orbital angular momentum quantum number,
and $m$ is the magnetic quantum number. The eigenfunctions of the Hamiltonian
$\psi_{nlm}(r,\theta,\phi)$ can also be separated into a radial part
$P_{nl}(r)$ and an angular part $Y_{l}^{m}(\theta,\phi)$:

$\psi(r,\theta,\phi)=\frac{1}{r}P_{nl}(r)Y_{l}^{m}(\theta,\phi),$ (4)

The angular part $Y_{l}^{m}(\theta,\phi)$ are the well-known spherical
harmonics that have analytical expressions. In contrast, the radial
wavefunction $P_{nl}(r)$(except for hydrogen) has to be solved numerically in
general with the following ordinary differential equation (ODE):

$\frac{d^{2}P(r)}{dr^{2}}=[\frac{l(l+1)}{r^{2}}-\frac{2m_{e}}{\hbar^{2}}(E+eV(r))]P(r).$
(5)

The central potential $V(r)$ includes the Coulomb potential from the nucleus
and the potential from the other electrons and is spherically averaged:

$\begin{split}V(r_{i})&=\frac{Ne}{r_{i}}-\sum_{j\neq
i}^{N}\frac{e}{|\vec{r_{i}}-\vec{r_{j}}|}\\\ &=\frac{Ne}{r_{i}}-\sum_{j\neq
i}^{N}\int\frac{e|\psi_{j}(\vec{r_{j}})|^{2}}{|\vec{r_{i}}-\vec{r_{j}}|}d\vec{r_{j}}\\\
&=\frac{Ne}{r_{i}}-\int\frac{e\rho(\vec{r^{\prime}})}{|\vec{r_{i}}-\vec{r^{\prime}}|}d\vec{r^{\prime}}\end{split}$
(6)

where $N$ is the positive charge of the nucleus and also the number of total
electrons for a neutral atom, $\vec{r_{i}}$ is the current position of $i$th
electron that we want to solve the wavefunction for and
$\psi_{j}(\vec{r_{j}})$ is the wavefunction of the $j$th electron with the
position $\vec{r_{j}}$. To calculate the potential for $i$th electron, we need
to sum the Coulomb interaction with all other electrons $j\neq i$.
Furthermore, this summation is equivalent to the integral of the charge
density $\rho(\vec{r^{\prime}})$ as a mean field shown in Eq. 6. One may find
the equations are nested in a quite interesting manner: solving the
wavefunctions of the current $i$th electron requires the potential in Eq. 5,
which requires solving the wavefunctions and corresponding charge density of
all other electrons in Eq. 6. As a consequence, the solution involves an
iterative process known as the self-consistent field (SCF) method which is the
key step in the Hartree-Fock-Slater approach or modern density functional
theory. For bound states ($E<0$) the radial function is normalized as follows:

$\int\psi^{*}(r,\theta,\phi)\psi(r,\theta,\phi)dr=\int_{0}^{\infty}P_{nl}^{2}(r)dr=1.$
(7)

For continuum states ($E>0$), the radial function is normalized by matching
the asymptotic form as a plane wave at a large distance away from the nucleus:

$P_{El^{\prime}}(r)\rightarrow\frac{1}{\sqrt{\pi
k}}\sin(kr-\frac{\pi}{2}l^{\prime}-\eta ln(2kr)+\delta_{l}),$ (8)

where $k=\sqrt{2m_{e}E}/\hbar$ is the wave vector and $\eta$ is the Sommerfeld
parameter written as:

$\eta=\frac{Z_{\infty}e^{2}}{\hbar v},$ (9)

with $Z_{\infty}=\lim_{r\rightarrow\infty}V(r)re^{-2}$ which approaches 0 for
finite-range potential, and $v$ is the non-realtivistic speed far away from
the atom. The phase shift $\delta_{l}$ can be calculated by solving the ODE
with the boundary condition $P_{El}(0)=0$ using the Numerov method and
matching the asymptotic form at large $r$.

###  2.2 Non-relativistic inelastic electron scattering

Following the derivations given in [19], starting from Fermi’s golden rule,
the transition rate from initial state $\ket{i}$ to final state $\ket{f}$ can
be written as:

$\tau_{i\rightarrow
f}=\frac{{2\pi}}{{\hbar}}|\bra{f}H^{\prime}\ket{i}|^{2}\rho(E_{i},E_{f}),$
(10)

where $\tau_{i\rightarrow f}$ represents the transition rate (transition
probability per unit time). $\bra{f}H^{\prime}\ket{i}$ denotes the matrix
element (in bra–ket notation) of the interaction Hamiltonian $H^{\prime}$
between the initial states $\ket{i}$ and final state $\ket{f}$ during the
transition. $\rho(E_{i},E_{f})$ represents the density of states at the
initial energy $E_{f}$ and the final energy. In this study,
$\rho(E_{i},E_{f})$ is taken as a delta function
$\rho(E_{i},E_{f})=\delta(E_{i}-E_{f}+\Delta E)$ with $\Delta E=E_{f}-E_{i}$
as the energy difference between initial and final states. Since one can
distinguish the probe electron and atomic electron (consequently the exchange
effects between them are negligible), the initial or final state can be
factorized as the wavefunction of the initial and final orbitals $\varphi_{i}$
and $\varphi_{f}$ (with N electrons in their atomic reference system
$r_{1}\ldots r_{N}$) coupled with the incoming or outgoing plane waves
$\Psi_{i}$ and $\Psi_{f}$ (in the laboratory reference system $R$), with
associated wave vector $\vec{k_{i}}$ and $\vec{k_{f}}$:

$\begin{split}\ket{i}=\ket{\varphi_{i}}\ket{\Psi_{i}}=\ket{\varphi_{i}(\vec{r_{1}}\ldots\vec{r_{N}})}\ket{e^{i\vec{k_{i}}\cdot\vec{R}}},\\\
\ket{f}=\ket{\varphi_{f}}\ket{\Psi_{f}}=\ket{\varphi_{f}(\vec{r_{1}}\ldots\vec{r_{N}})}\ket{e^{i\vec{k_{f}}\cdot\vec{R}}}.\end{split}$
(11)

To satisfy the anti-symmetry requirement (Pauli principle) upon exchange of
two electrons in the many-electron system, one typically uses the Slater
determinant to represent the atomic electron wavefunctions $\varphi_{i}$ and
$\varphi_{f}$. The Slater determinant is written as:

$\ket{\varphi}=\frac{1}{\sqrt{N!}}\left|\begin{array}[]{cccc}\psi_{1}\left(\vec{r}_{1}\right)&\psi_{2}\left(\vec{r}_{1}\right)&\cdots&\psi_{N}\left(\vec{r}_{1}\right)\\\
\psi_{1}\left(\vec{r}_{2}\right)&\psi_{2}\left(\vec{r}_{2}\right)&\cdots&\psi_{N}\left(\vec{r}_{2}\right)\\\
\vdots&\vdots&\ddots&\vdots\\\
\psi_{1}\left(\vec{r}_{N}\right)&\psi_{2}\left(\vec{r}_{N}\right)&\cdots&\psi_{N}\left(\vec{r}_{N}\right)\end{array}\right|,$
(12)

The interaction Hamiltonian is taken as the Coulombic potential between the
probe electron and the target atom similar to Eq. 6, we can evaluate the
transition matrix element by integrating in real space. The real space
integral can be simplified using the Fourier transform of the Coulombic
potential to the momentum space:

$\int
d^{3}\vec{R}\frac{e^{2}}{|\vec{R}-\vec{r_{j}}|}e^{i\vec{q}\cdot\vec{R}}=\frac{4\pi
e^{2}}{q^{2}}e^{i\vec{q}\cdot\vec{r_{j}}},$ (13)

with $\vec{R}-\vec{r_{j}}$ defined as the vector between the probe electron
and the $j$-th atomic electron. After the integration of the laboratory real
space $\vec{R}$, the transition matrix element in the momentum space can be
evaluated by integration of the atomic reference system for each electron
$(\vec{r_{1}}\ldots\vec{r_{N}})$ as:

$\begin{split}\bra{f}V\ket{i}&=\frac{e^{2}}{(2\pi)^{2}q^{2}}(\bra{\varphi_{f}}\sum_{j}^{N}e^{i\vec{q}\vec{r_{j}}}\ket{\varphi_{i}}-N\bra{\varphi_{f}}\ket{\varphi_{i}})\\\
&=\frac{e^{2}}{(2\pi)^{2}q^{2}}\bra{\psi^{j}_{f}}e^{i\vec{q}\vec{r_{j}}}\ket{\psi^{j}_{i}}.\end{split}$
(14)

Several observations can be made: (1) the plane wave with wave vector of the
momentum transfer $\vec{q}=\vec{k_{i}}-\vec{k_{f}}$ emerges as a result of the
substitution of the explicit form of the initial and final states defined in
Eq. 11 into the matrix element in Eq. 10 and Fourier transform of the
Coulombic potential in Eq. 13; (2) the contribution to the transition rate
from the nuclear electrostatic potential is zero when the initial and final
state are orthogonal to each other; (3) only the excited orbital (for
instance, the $j$th electron in the equation above) contributes to the
transition matrix element due to orthonormality, if we assume all the other
orbitals remain unchanged under the frozen core approximation.

The transition rate is dependent on the current density of the quantum system,
a more useful definition is the so-called double differential cross-section,
which is defined as the probability of the transition per unit solid angle and
energy collected. The double differential scattering cross-section (DDSCS) can
be written as:

$\begin{split}\frac{\partial^{2}\sigma}{\partial
E\partial\Omega}&=(\frac{2\pi}{\hbar})^{4}(\gamma
m_{e})^{2}\sum_{i,f}\frac{k_{f}}{k_{i}}|\bra{f}V\ket{i}|^{2}\delta(E_{i}-E_{f}+\Delta
E)\\\ &=\frac{4\gamma^{2}}{a_{0}^{2}q^{4}}\frac{k_{f}}{k_{i}}|S(q,E)|^{2},\\\
&=\frac{4\gamma^{2}}{q^{2}}\frac{R}{\Delta E}\frac{k_{f}}{k_{i}}f(q,E)\\\
\end{split}$ (15)

where $\gamma$ is the Lorentz factor for relativistic scalar correction,
$a_{0}$ is the Bohr radius, $R$ is the Rydberg energy, $k_{i}$ and $k_{f}$ are
the wave number of the incoming and outgoing plane waves, respectively. To
separate the part that is independent of the incoming electron energy, we
define the so-called GOS as $f(q,E)$:

$f(q,E)=\frac{\Delta E}{R}\frac{\left|S(q,E)\right|^{2}}{(qa_{0})^{2}}.$ (16)

Here $|S(q,E)|^{2}$ is known as the dynamic form factor, which is the squared
magnitude of the transition matrix element summed over the final states as:

$S(q,E)=\sum_{\psi_{i},\psi_{f}}\bra{\psi_{f}}e^{i\vec{q}\cdot\vec{r}}\ket{\psi_{i}}$
(17)

Note that we take the transition matrix only contributed by the excited
orbital and the j symbol in Eq. 14 is dropped thereafter for simplicity. The
initial and final states of the atomic orbitals are decomposed in terms of
radial and angular parts, which has been explained in Sec. 2.1. For the
initial states, the summation is applied for the magnetic quantum number $m$
since we already chose the orbital. For the final states, the summation
extends over both the angular quantum number $l^{\prime}$ and the magnetic
quantum number $m^{\prime}$ for each specified final state energy $E$.

To derive the dynamic structure factor, one can expand the plane wave using
spherical harmonics and integrate the radial and angular parts separately. The
plane wave expansion can be written as [64]:

$e^{iqr}=\sum_{\lambda=0}^{\infty}\sum_{m_{\lambda}=-\lambda}^{\lambda}4\pi
i^{\lambda}(2\lambda+1)j_{\lambda}(qr)Y_{\lambda}^{m_{\lambda}}(\theta,\phi)$
(18)

where $j_{\lambda}(qr)$ is the spherical Bessel function of the order
$\lambda$. First, we can collect all the radial components from the initial
states, the final states and the plane wave expansion for the radial integral
$R^{\lambda}_{nl,El^{\prime}}$:

$R^{\lambda}_{nl,El^{\prime}}(q,E)=\int_{0}^{\infty}\left[P_{nl}(r)P_{El^{\prime}}(r)\right]j_{\lambda}(qr)dr.$
(19)

Then we can find the angular integral of the products of three spherical
harmonics leads to the Wigner 3j symbol $(^{...}_{...})$, also known as Gaunt
coefficients:

$\int
Y_{\lambda}^{m_{\lambda}}(\theta,\phi)Y_{l^{\prime}}^{m_{l^{\prime}}}(\theta,\phi)Y_{l}^{m_{l}}(\theta,\phi)d\Omega=\sqrt{\frac{(2l+1)(2l^{\prime}+1)(2\lambda+1)}{4\pi}}\left(\begin{array}[]{ccc}l^{\prime}&\lambda&l\\\
0&0&0\end{array}\right)\left(\begin{array}[]{ccc}l^{\prime}&\lambda&l\\\
m^{\prime}&m_{\lambda}&m\end{array}\right).$ (20)

When summing the magnetic quantum numbers $m$ and $m^{\prime}$ of the initial
and final states in Eq. 20, the expression can be simplified by applying the
orthogonal relation of the Wigner 3j symbols, which results in a delta
function and a triangle condition:

$\begin{split}\sum_{m,m^{\prime}}\left(\begin{array}[]{ccc}l^{\prime}&\lambda&l\\\
m^{\prime}&m_{\lambda}&m\end{array}\right)\left(\begin{array}[]{ccc}l^{\prime}&\lambda^{\prime}&l\\\
m^{\prime}&m_{\lambda}^{\prime}&m\end{array}\right)&=\frac{1}{2\lambda+1}\delta_{\lambda^{\prime}\lambda}\delta_{m_{\lambda}m_{\lambda}^{\prime}}\Delta(l,l^{\prime},\lambda),\\\
\text{with triangle condition
}\Delta\left(l,l^{\prime},\lambda\right)&=\begin{cases}1&\text{ if
}\left|l-l^{\prime}\right|\leq\lambda\leq l+l^{\prime}\\\ 0&\text{ otherwise
}\end{cases}\end{split}$ (21)

With the radial and angular integral simplified, the dynamic form factor can
be written as:

$\left|S(q,E)\right|^{2}=\sum_{l^{\prime}}\sum_{\lambda}[l,l^{\prime},\lambda]{R^{\lambda}_{nl,El^{\prime}}}^{2}\left(\begin{array}[]{ccc}l^{\prime}&\lambda&l\\\
0&0&0\end{array}\right)^{2},$ (22)

where the abbreviated notion of
$[l,l^{\prime},\lambda]=(2l+1)(2l^{\prime}+1)(2\lambda+1)$ is used. The
dynamic form factor can be physically interpreted as the space-time Fourier
transform of the density autocorrelation during the transition, as derived in
[19].

The calculations of theoretical spectra, incorporating experimental parameters
and Generalized Oscillator Strengths (GOS), will be detailed in Section 4
since the EELS geometric considerations are independent of the choice between
Schrödinger or Dirac solutions.

##  3 Theory: relativistic atomic orbitals and inelastic electron scattering

The Dirac equation combines special relativity and quantum mechanics – the
cornerstone of modern physics – for relativistic electrons [21]. In this
section, we provide an overview of the Dirac equation for atomic orbitals and
inelastic scattering. For a comprehensive treatment of this topic, readers are
referred to books on relativistic quantum mechanics [65] and recent papers on
inelastic electron scattering [50, 66] as here we only outline the essential
equations.

###  3.1 Dirac equations for atomic orbitals

The Dirac equation is a first-order differential equation that describes the
behavior of relativistic electrons in the presence of an electromagnetic field
(minimum coupling), which is written as:

$(c\gamma^{\mu}(i\hbar\partial_{\mu}+eA_{\mu})-m_{e}c^{2})\psi=0,$ (23)

where the summation is implicitly applied over the values of the index
$\mu=0,1,2,3$ in the Einstein summation notation. For instance, the
multiplication of two arbitrary 4-vectors $\mathbf{A}$ and $\mathbf{B}$ is
written as:

$\mathbf{A}\cdot\mathbf{B}=A_{\mu}B^{\mu}=A^{\nu}B_{\nu}=A_{\mu}\eta^{\mu\nu}B_{\nu}=a^{0}b^{0}-\vec{a}\cdot\vec{b},$
(24)

where $A_{\mu}$ are the covariant 4-components and $A^{\mu}$ are the
contravariant 4-components of the vector $\mathbf{A}$ in the tensor index
notation, which relates to each other via the Minkowski metric $\eta^{\mu v}$.
The multiplication is written in such a manner for the invariance under the
Lorentz transformation of the space-time. Back to the Dirac equation,
$\partial_{\mu}=(\frac{1}{c}\frac{\partial}{\partial t},\vec{\nabla})$ are the
4-gradient for space-time vectors $x^{\mu}=(ct,\vec{r})$, which gives the
4-momentum $p_{\mu}=(\mathbb{E}/c,-\hat{p})$ where $\mathbb{E}$ is the
relativistic energy of the electron. $A_{\mu}=(V/c,-\vec{A})$ is the
4-potential consists of electrostatic potential $V$ and magnetic vector field
potential $\vec{A}$. The Dirac gamma matrices $\gamma^{\mu}$ are defined as:

$\gamma^{0}=\begin{pmatrix}I_{2}&0\\\
0&-I_{2}\end{pmatrix},\gamma^{i}=\begin{pmatrix}0&\sigma_{i}\\\
-\sigma_{i}&0\end{pmatrix},$ (25)

where $I_{2}$ is the 2 $\times$ 2 identity matrix and $\sigma_{i}$ represents
the Pauli matrices:

$\sigma_{1}=\begin{pmatrix}0&1\\\
1&0\end{pmatrix},\sigma_{2}=\begin{pmatrix}0&-i\\\
i&0\end{pmatrix},\sigma_{3}=\begin{pmatrix}1&0\\\ 0&-1\end{pmatrix}.$ (26)

We can express the Dirac equation in the 2-component form as follows:

$\begin{split}\begin{pmatrix}\mathbb{E}+eV(r)-m_{e}c^{2}&-c{\sigma}\cdot(\hat{p}+e\vec{A})\\\
-c{\sigma}\cdot(\hat{p}+e\vec{A})&\mathbb{E}+eV(r)+m_{e}c^{2}\end{pmatrix}\begin{pmatrix}\psi_{1}\\\
\psi_{2}\end{pmatrix}&=\begin{pmatrix}0\\\ 0\end{pmatrix},\end{split}$ (27)

where the upper component $\psi_{1}$ and lower component $\psi_{2}$ are
coupled in Eq. 27. We may define the energy in the excess of the electron rest
energy as $E=\mathbb{E}-m_{e}c^{2}$ for ease of comparison with the energy in
the Schrödinger equation. For an atomic electron in the presence of a radially
symmetric electrostatic potential, we may only consider the electrostatic
potential $V(r)$ and ignore the magnetic vector field potential $A_{\mu}$. In
the low-energy limit $E+eV(r)<<mc^{2}$, we can see that the lower component
$\psi_{2}$ is much smaller than $\psi_{1}$ (as
$\psi_{2}\approx\frac{\sigma\cdot\hat{p}}{2m_{e}c}\psi_{1}$), which is at the
order of expectation value of the speed relative to the light $\frac{\langle
v\rangle}{c}$. For the H atom, this magnitude is at the fine-structure
constant ($\sim$ 1/137) but gradually increases with atomic number as the
momentum increases. Hence, the upper and lower components are often referred
to as the large and small components, respectively.

The solution can be factorized into radial part and angular parts:

$\psi=\frac{1}{r}\begin{pmatrix}P_{n\kappa}(r)\mathscr{Y}_{jlm}(\theta,\phi)\\\
iQ_{n\kappa}(r)\mathscr{Y}_{jl^{\prime}m}(\theta,\phi)\end{pmatrix},$ (28)

As $l$ is no longer a good quantum number in the Dirac equation. Here, we
introduce the total angular momentum operator defined as
$\hat{J}=\hat{L}+\hat{S}$ with its eigenvalue $J^{2}=(j+1)j\hbar^{2}$ and
$J_{z}=m_{j}\hbar$, and $j$ is bounded by $|l-s|\leq j\leq|l+s|$ in magnitude.
$\hat{S}$ is the spin operator with eigenvalues of $S^{2}=s(s+1)\hbar^{2}$ and
$S_{z}=m_{s}\hbar$, where $s=1/2$ and $m_{s}=-1/2$ or $1/2$ for electron.
$P_{n\kappa}(r)$ and $Q_{n\kappa}(r)$ are the radial parts of the large and
small components, respectively. The relativistic quantum number $\kappa$ is
defined as:

$\begin{split}\kappa&=(l-j)(2j+1)=\begin{cases}-l-1,&\text{ if
}j=l+\frac{1}{2}\text{ (spin up)}\\\ l,&\text{ if }j=l-\frac{1}{2}\text{ (spin
down)}\end{cases}\\\ \end{split}.$ (29)

For the angular part, $\mathscr{Y}_{jlm}(\theta,\phi)$ represents the spinor
spherical harmonics (note in the small component $l^{\prime}=2j-l$ in Eq. 28),
which can be written as:

$\mathscr{Y}_{jlm}(\theta,\phi)=\sum_{m_{l},m_{s}}\langle
l,1/2,m_{l},m_{s}\mid j,m_{j}\rangle
Y_{l}^{m_{l}}(\theta,\phi)\chi_{\frac{1}{2},m_{s}},$ (30)

where $\langle l,1/2,m_{l},m_{s}\mid j,m_{j}\rangle$ is the Clebsch-Gordan
coefficient, $Y_{l}^{m_{l}}(\theta,\phi)$ is the usual spherical harmonics
with $m_{l}=m_{j}-m_{s}$ and $\chi_{\frac{1}{2},m_{s}}$ is the spin
eigenfunction of a spin 1/2 particle:

$\chi_{\frac{1}{2},+\frac{1}{2}}=\left(\begin{array}[]{c}1\\\
0\end{array}\right)\quad\text{ and
}\quad\chi_{\frac{1}{2},-\frac{1}{2}}=\left(\begin{array}[]{c}0\\\
1\end{array}\right)$ (31)

For the radial part, the Dirac equation can be written as coupled ODE
equations:

$\begin{split}\frac{\mathrm{d}P_{n\kappa}}{\mathrm{d}r}&=-\frac{\kappa}{r}P_{n\kappa}+\frac{E+2\mathrm{~{}m}_{\mathrm{e}}c^{2}+eV}{c\hbar}Q_{n\kappa},\\\
\frac{\mathrm{d}Q_{n\kappa}}{\mathrm{d}r}&=\frac{-E-eV}{c\hbar}P_{n\kappa}+\frac{\kappa}{r}Q_{n\kappa}.\end{split}$
(32)

For the bound state ($E<0$), the radial part can be normalized as:

$\int_{0}^{\infty}\left(P_{n\kappa}^{2}(r)+Q_{n\kappa}^{2}(r)\right)\mathrm{d}r=1.$
(33)

For the free state ($E>0$), the radial part can be normalized by matching the
asymptotic behavior:

$P_{E\kappa^{\prime}}\rightarrow\frac{1}{\sqrt{\pi
k}}\sin\left(kr-\frac{\pi}{2}\kappa^{\prime}-\eta
ln(2kr)+\delta_{\kappa^{\prime}}\right),$ (34)

where $k=\sqrt{E(E+2m_{e}c^{2})}/\hbar c$ is the relativistic wave number, and
$\delta_{\kappa^{\prime}}$ is the phase shift. $\eta$ is the relativistic
Sommerfeld parameter with the corresponding relativistic speed in Eq. 9. In
the low-energy limit, $E+eV<<m_{e}c^{2}$, the radial part can be written as
[66]:

$Q_{n\kappa}=\frac{c\hbar}{2\mathrm{~{}m}_{\mathrm{e}}c^{2}}\left(\frac{\kappa}{r}P+\frac{\mathrm{d}P}{\mathrm{~{}d}r}\right)$
(35)

and

$\frac{d^{2}P_{n\kappa}(r)}{dr^{2}}=[\frac{\kappa(\kappa+1)}{r^{2}}-\frac{2m_{e}}{\hbar^{2}}(E+eV(r))]P_{n\kappa}(r)$
(36)

Note that $\kappa(\kappa+1)=l(l+1)$, reducing the upper component
$P_{n\kappa}(r)$ radial equation to the Schrödinger radial equation as in Eq.
5.

###  3.2 Relativistic inelastic electron scattering

From the perspective of electrodynamics, the motion of a charged particle
induces electric and magnetic fields and hence perturbs the interaction
Hamiltonian. Instead of taking the rigorous quantum field theory, here we
follow perturbation theory due to its simplicity, as well outlined in [31, 54,
55, 56]. The DDSCS for inelastic electron scattering is given by:

$\frac{\partial^{2}\sigma}{\partial
E\partial\Omega}={\left(\frac{2\gamma}{a_{0}}\right)^{2}\frac{1}{\left(q^{2}-(\Delta
E/\hbar
c)^{2}\right)^{2}}\frac{k_{f}}{k_{i}}\sum_{\psi_{i},\psi_{f}}\mid\langle\psi_{f}|e^{i\vec{q}{\vec{r}}}\left(1-\frac{\hat{p}\vec{v}_{0}}{m_{e}c^{2}}\right)}\left.|\psi_{i}\rangle\right|^{2}\delta\left(E_{i}-E_{f}+\Delta
E\right),$ (37)

where $v_{0}$ is the velocity of the incident probe electron with energy
$E_{0}$. By comparing Eq.37 with the non-relativistic case in Eq.
LABEL:eq:ddscs, we can note that $q^{2}$ is replaced by $q^{2}-(\Delta E/\hbar
c)^{2}$ and there is an extra term $\frac{\hat{p}\vec{v}_{0}}{m_{e}c^{2}}$,
which originates from the perturbation of the fast moving charge under the
Coulomb gauge leading to the contraction of the scattering vector in the
direction of the incident beam (retardation effect) [55]. This is equivalent
to the z-component of the 4-transition current under the Lorentz gauge [56].
The relativistic effect of incident fast electron plays a critical role in the
experimental observation of the magic angle in EELS [55] and the detailed
intensity distribution in energy-filtered electron diffaraction [56], thus it
is important to include it correctly for quantitative work. By projecting the
contribution of $\frac{\hat{p}\vec{v}_{0}}{m_{e}c^{2}}$ to directions that are
parallel and normal to the incident electron, we can further simplify the
DDSCS as proposed in [31, 54]:

$\frac{\partial\sigma}{\partial
E\partial\Omega}=\left(\frac{2\gamma}{a_{0}}\right)^{2}\frac{k_{f}}{k_{i}}\left[\frac{1}{q^{4}}+\frac{\beta_{t}^{2}(\Delta
E/\hbar c)^{2}}{(q^{2}-(\Delta E/\hbar
c)^{2})^{2}}\right]\sum_{\psi_{i},\psi_{f}}|\langle\psi_{f}|e^{i\vec{q}\vec{r}}|\psi_{i}\rangle|^{2},$
(38)

where $\beta_{t}$ is the transverse component of the incident electron
velocity (in units of the speed of light) $\beta=v_{0}/c$ with respect to the
scattering vector $\vec{q}$:

$\beta^{2}_{t}=\beta^{2}-\frac{\Delta E^{2}}{(cq)^{2}}(1+\frac{(cq)^{2}-\Delta
E^{2}}{2\Delta E(E_{0}+m_{e}c^{2})})$ (39)

For the orbital relativistic effects, the transition matrix element must
account for the Dirac solution, which includes both the large and small
components [31]. As a consequence, the radial overlap integral is changed from
Eq. 19 to:

$R^{\lambda}_{n\kappa,E\kappa^{\prime}}(q,E)=\int_{0}^{\infty}[P_{n\kappa}(r)P_{E\kappa^{\prime}}(r)+Q_{n\kappa}(r)Q_{E\kappa^{\prime}}(r)]j_{\lambda}(qr)dr.$
(40)

For the angular part, the spherical harmonics in the Schrödinger solution are
now replaced by the spinor spherical harmonics in the Dirac solution. This
addition introduces a Wigner 6j $\\{_{...}^{...}\\}$ term when integrating
over the angle and summing over the states. Overall, the dynamic form factor
in the Dirac solution is given by:

$\left|S(q,E)\right|^{2}=\sum_{l^{\prime}}\sum_{\lambda}\sum_{\kappa^{\prime}}[l,l^{\prime},j^{\prime},\lambda]{R^{\lambda}_{n\kappa,E\kappa^{\prime}}}^{2}\left(\begin{array}[]{ccc}l^{\prime}&\lambda&l\\\
0&0&0\end{array}\right)^{2}\left\\{\begin{array}[]{ccc}j&\lambda&j^{\prime}\\\
l^{\prime}&\frac{1}{2}&l\end{array}\right\\}^{2},$ (41)

The GOS still follows the same definition as in the non-relativistic case in
Eq. 16.

To check we have the correct implementation of the relativistic correction
defined in Eq. 38, we can use the analytical expression for the relativistic
correction ratio under dipole approximation, which is defined as the ratio of
the scattering cross-sections with and without the relativistic
electrodynamics effects. This ratio is given by [50, 55, 56]:

$\sigma_{\text{ratio }}(\Delta E)=\frac{\sigma_{\mathrm{rel}}(\Delta
E)}{\sigma_{\text{conv }}(\Delta
E)}=\left[\ln\left(1+\frac{x_{\theta}}{1-\beta^{2}}\right)-\frac{\beta^{2}x_{\theta}}{1-\beta^{2}+x_{\theta}}\right]\frac{1}{\ln\left(1+x_{\theta}\right)},$
(42)

where $x_{\theta}=\theta_{0}^{2}/\theta^{2}_{E}$ is the squared ratio of the
collection angle $\theta_{0}$ to the characteristic angle $\theta_{E}$, which
is written as:

$\theta_{E}=\frac{\Delta E}{\gamma mv_{0}^{2}}$ (43)

In Section 6.2, we will examine the relativistic correction ratio for cross-
sections under different experimental conditions.

##  4 EELS geometry parameters in momentum space integration

In this section, we will explain how to consider the experimental EELS
geometry (i.e. EELS collection angle and STEM convergence angle) for cross-
section calculations. From the dynamic form factor in Eq. 22 and Eq. 41 (or
the associated GOS in Eq. 16), one can integrate over the EELS collection
angle to obtain the differential cross section as a function of energy loss
[19]. The differential cross section is given by:

$\frac{d\sigma}{dE}=\frac{4\pi\gamma^{2}}{k_{i}^{2}}\int^{Q_{\text{max}}}_{Q_{\text{min}}}\frac{|S(q,E)|^{2}}{Q}d(ln(Q)),$
(44)

where we take $Q=(a_{0}q)^{2}$ for ease of integration. From the scattering
geometry, we can set the upper and lower limits of the integration as:

$\begin{split}q_{\text{min}}=k_{i}-k_{f},\\\
q_{\text{max}}=\sqrt{{k_{i}}^{2}+{k_{f}}^{2}-2{k_{i}}{k_{f}}\cos(\beta_{0})},\\\
\end{split}$ (45)

as the minimum scattering angle is zero and the maximum scattering angle is
bound by the EELS collection aperture $\beta_{0}$. The differential cross
section as a function of energy loss is used to fit and quantify the
experimental EELS spectra.

For plane wave illumination in TEM, The above Eq. 44 is sufficient. For STEM
convergent beam illumination, we need to take into account the finite angle of
both the convergence angle and collection angle, which can be performed by the
geometric correction procedure proposed by Kohl [67]. In this approach, we
consider the STEM illumination as a collection of plane waves with their
corresponding wave vectors bounded by the convergence angle $\alpha_{0}$. The
EELS collection angle $\beta_{0}$ defines the maximum angle relative to the
incident direction to be recorded by the detector. For a given scattering
angle $\theta$, there could be different combinations of incident angle and
exit angles allowed by the geometry. As a consequence, the DDSCS should be
multiplied by a correction factor corresponding to those combinations. Based
on the geometric interpretation, this problem is equivalent to the
intersection area of two circles – one with a radius of $\alpha_{0}$ and
another with a radius of $\beta_{0}$ and they are seperated by a distance of
$\theta$. The resulting correction factor is a cross-correlation function,
which is written as [67]:

$\displaystyle F_{\mathrm{BF}}\left(\alpha_{0},\beta_{0},\theta\right)$ (46)
$\displaystyle\quad=\left\\{\begin{array}[]{l}\theta_{<}^{2}/\alpha_{0}^{2},\quad\text{
for }0\leqslant\theta\leqslant\left|\alpha_{0}-\beta_{0}\right|,\\\
\pi^{-1}\left[\arccos(x)+\left(\beta_{0}^{2}/\alpha_{0}^{2}\right)\arccos(y)\right.\\\
\left.-\left(1/2\alpha_{0}^{2}\right)\sqrt{4\alpha_{0}^{2}\beta_{0}^{2}-\left(\alpha_{0}^{2}+\beta_{0}^{2}-\theta^{2}\right)^{2}}\right],\text{
for }\left|\alpha_{0}-\beta_{0}\right|<\theta<\alpha_{0}+\beta_{0},\\\
0,\text{ otherwise. }\end{array}\right.$

where

$\begin{split}\theta_{<}&=\min\left(\alpha_{0},\beta_{0}\right),\\\
x&=\frac{\alpha_{0}^{2}+\theta^{2}-\beta_{0}^{2}}{2\alpha_{0}\theta},\\\
y&=\frac{\beta_{0}^{2}+\theta^{2}-\alpha_{0}^{2}}{2\beta_{0}\theta}.\end{split}$
(47)

The EELS correction factor $F_{BF}$ (closely related to the contrast transfer
function of bright field incoherent imaging) is used to correct the DDSCS for
a given set of scattering angles in the STEM-EELS geometry. The effective
partial scattering cross-section after integration in the momentum space is
written as [67]:

$\frac{d\sigma_{\text{eff}}}{dE}=\int_{0}^{\alpha_{0}+\beta_{0}}F_{BF}(\alpha_{0},\beta_{0},\theta)\frac{d^{2}\sigma}{dEd\Omega}2\pi\theta
d\theta,$ (48)

or in terms of the Q space integral as in Eq. 44

$\begin{split}\frac{d\sigma_{\text{eff}}}{dE}&=\frac{4\pi\gamma^{2}}{k_{i}^{2}}\int^{Q_{\text{max}}}_{Q_{\text{min}}}\frac{|S(q,E)|^{2}}{Q}F_{BF}(\alpha_{0},\beta_{0},\theta)d(ln(Q)),\\\
\text{ with
}\theta&=\arccos\left(\frac{{k_{i}}^{2}+{k_{f}}^{2}-{q}^{2}}{2{k_{i}}{k_{f}}}\right),\\\
q_{\text{min}}&={k_{i}}-{k_{f}},\\\
q_{\text{max}}&=\sqrt{{k_{i}}^{2}+{k_{f}}^{2}-2{k_{i}}{k_{f}}\cos(\alpha_{0}+\beta_{0})},\end{split}$
(49)

##  5 Computational details

In this study, we only perform GOS calculations for single atoms, thereby
disregarding the solid state effects in EELS. The atomic orbital calculations
are performed using the Flexible Atomic Code (FAC) [68] package. To solve the
Dirac radial wavefunction as described in Eq. 32, FAC considers the local
atomic potential as a combination of the electron-nucleus and electron-
electron Coulomb interactions. The solution is computed through a self-
consistent Dirac–Fock–Slater method within a modified local density
approximation for the correct asymptotic behavior of the exchange energy [68].
We take the bound state of the target orbital as the initial state and the
continuum state as the excited state. Continuum states are treated in the
distorted-wave approximation [68]. Fermi’s golden rule employed a delta
function to represent the density of states, thus ignoring the fine
structures. In addition, the orthogonality of the initial and final states is
necessary to exclude the contribution of the nucleus in Eq. 14, which
otherwise causes a divergence of the cross-sections when momentum transfer
approaches zero. To ensure orthogonality, we take the frozen core
approximation that assumes the potential remains unchanged for the initial and
final states, which means only the exited orbital contributes to the
transition matrix in Eq. 14. The calculated large and small components are
used in the Eq. 40 to compute the transition matrix element. The convergence
is checked by including final states (summation of $\kappa^{\prime}$ in Eq.
41) until the contribution of the last final state falls below 0.1%.

The GOS calculation is computationally intensive since we need to sum over
different final states for each energy loss. To accelerate the process, we
developed a highly efficient parallelization scheme to take advantage of
modern computer with multiple CPUs. Specifically, the calculation of the
transition matrix is evaluated using the multi-processing of different energy
losses. For a given energy loss, multi-threading is used for all the final
states $\kappa^{\prime}$ until convergence. We implemented a dynamic
predicting, caching, and dumping strategy for the number of final state
wavefunctions needed to manage memory usage during parallelization. After
optimization, the computation of the entire database takes about 2 days on a
desktop (intel i9-7900X 10 CPUs @4.5GHz), and less than 4 hours on a single
node on a high-performance computing cluster (2x 64-core AMD Epyc 7H12 CPU
@2.6 GHz).

We constructed the GOS database for all elements with atomic number Z ranging
from 1 to 118 and all available shells. This resulted in a total of 2143
entries tabulated as a function of energy loss $\Delta E$ and momentum
transfer $q$ with fine sampling. More specifically, we computed the GOS in an
energy range of 0.01 - 4000 eV above the ionization energy with 128 sampling
points and an adaptive momentum range of ($q_{\text{min}}$, $q_{\text{max}}$)
with 256 sampling points. Both energy and momentum sampling follows a log
sampling scheme, as also done in [26, 28, 60]. The minimum momentum transfer
$q_{\text{min}}$ is in the direction of the incident beam, defined as:

$q_{\text{min}}=\frac{\Delta E_{\text{min}}}{\hbar c},$ (50)

where $\Delta E_{min}$ corresponds to the minimum energy loss which is the
ionization energy. The maximum momentum sampling $q_{\text{max}}=2q_{r}$,
where $q_{\text{r}}$ is the momentum transfer at the Bethe ridge for a given
energy loss [10]:

$(a_{0}q_{\text{r}})^{2}=\frac{\Delta E}{R}+\frac{\Delta
E^{2}}{2m_{e}c^{2}R}.$ (51)

Adaptive momentum sampling is developed in such a manner to maximize the
physical information for a given finite number of sampling points. For
instance, an ionization event is always associated with a scattering vector
larger than $q_{\text{min}}$, thus sampling below the lower limit would be
pointless. Also, the transition probability beyond $q_{\text{max}}$ is well
decayed to negligible. Note that the ($q_{\text{min}}$, $q_{\text{max}}$) is
elemental edge dependent, so the sampling has to be adaptive instead of fixed
for all the edges. The database was saved in an HDF5 file following the GOSH
file format [62] so that users can switch easily between different variations
of GOS databases.

We also used the GPAW software to calculate the all-electron Schrödinger DFT
using the projector augmented wave method [69] with generalized gradient
approximation of the potential [70]. The DFT calculations are performed to
compare the orbital energies with Dirac results in this work but not for GOS
computation.

##  6 Results and discussions

In this section, we will first give an overview of the Dirac-based GOS
database compared to other Schrödinger solutions in Sec. 6.1, highlighting the
effects of SOC on the ionization energy, the wavefunctions, and the GOS. In
Sec. 6.2, we will demonstrate the necessity of considering the relativistic
effects of fast electrons in quantitative EELS. In Sec. 6.3, we will compare
the computed scattering cross-sections from the Dirac-based GOS with other
Schrödinger-based GOS databases, showing the general agreement and fine
differences. In Sec. 6.4, we will provide a step-by-step guide on how to use
the Dirac-based GOS database for EELS quantification.

###  6.1 Overview of the Dirac-based GOS database

Figure 1: Comparison of the GOS databases for the C K-edge calculated from (a)
the Hartree-Fock solution, (b) the DFT solution, and (c) our Dirac solution as
a function of momentum transfer and energy loss above the ionization. The
sampling grid is indicated by the white lines. The log sampling is dense at
low energy loss and low momentum transfer region, where most experiments focus
on.

A distinct improvement of our GOS database is the large energy-momentum range
with a fine sampling. Such sampling requires significant computational
efforts, which once limited the HF GOS, are now more manageable owing to the
rapid development of computation hardware and highly optimized parallelization
algorithms. For comparison, Fig. 1 shows the GOS database for the C-K edge
constructed previously based on (a) the HF solution (used in Gatan Digital
Micrograph), (b) the DFT solution, and (c) the Dirac solution. For this
particular case, HF GOS reaches a maximum energy of 200 eV above the onset
with 12 sampling points and a maximum momentum of 12.7 $\textup{~{}\AA}^{-1}$
with 20 sampling points. For conventional cross-section integration, one may
argue that 200 eV should be sufficient for an integration energy window but
not enough for model-based fitting. In contrast, the Dirac GOS extends the
sampling to 4000 eV above the ionization and 67 $\textup{~{}\AA}^{-1}$ with
128 and 256 sampling points, respectively. The sampling scheme enables the
clear visualization of the Bethe ridge with locally high transition
probabilities in a parabolic region ($\Delta E\sim q^{2}$) in the energy-
momentum space. In practice, the extended energy loss range of GOS can be
useful for the simultaneous model-based fitting of multiple elements. For
example, the presence of carbon support is common in TEM sample preparation.
In the case of a limited GOS energy range, the background fitting may only use
the pre-edge region to extrapolate a power law function instead of taking the
exact decaying tail of the C K-edge, which introduces additional bias for
model-based fitting of multiple elements. Fine sampling is also helpful for
the accurate generation of the EELS cross-sections and consequently spectrum
quantification. Note that the large sampling range and fine sampling step is a
common practice in modern computing, which is also available at recent DFT-
based GOS [60] shown in Fig. 1 (b) but with a momentum sampling fixed around
63 $\textup{~{}\AA}^{-1}$ and a smaller energy range (up to $\sim$ 1900 eV)
for all elements.

Table 1: The onset energies of different edges, with the experimental measurement from EELS Atlas, XPS, and the theoretical predictions from Dirac and all-electron Schrödinger solutions. Element | Orbital | Edge | Experimental Measurement (eV) | Theoretical Calculation (eV)
---|---|---|---|---
|  |  | EELS Atlas[71] | XPS[72] | Dirac | Schrödinger
Ti |  2p_1/2 |  L2 | 462 | 460 | 454 | 444
2p_3/2 |  L3 | 456 | 454 | 448
Ag |  3d_3/2 |  M4 | 373 | 374 | 370 | 358
3d_5/2 |  M5 | 367 | 368 | 364
Au |  4f_5/2 |  N6 | 86 | 88 | 84 | 81
4f_7/2 |  N7 | 83 | 84 | 80

The key benefit of employing the Dirac equation lies in its intrinsic
inclusion of relativistic effects, rather than applying corrections to the
Schrödinger equation. In the context of EELS, a notable advantage is the
precise prediction of spin-orbit splittings of different elements. Table 1
presents Ti L2/L3, Ag M4/M5, and Au N6/N7 edges with experimental measurements
from the EELS Atlas [71], X-ray Photoelectron Spectroscopy (XPS) [72] and
theoretical ionization energies calculated using the Dirac and Schrödinger
equations for each orbital. For instance, the Ti L2 edge is observed to be 6
eV higher than the L3 edge in both EELS and XPS, with the splitting known due
to spin-orbit coupling (SOC) [63]. From the perspective of applying
relativistic corrections to the Schrödinger equation, a spin-down electron has
lower energy than that of a spin-up electron with a SOC contribution that
scales with $\sim\langle{\hat{L}\cdot\hat{S}}\rangle$. Such measurements align
well with predictions from the Dirac equation (454 / 448 eV). In contrast, the
Schrödinger equation predicts a single onset energy for the 2p orbital
regardless of spin variations. However, it’s important to note that the Dirac
values don’t exactly match experimental measurements. This discrepancy arises
because we only performed the calculation for a single neutral atom, thereby
ignoring chemical shifts (as observed in ions with different oxidation states,
for example, the Ti edges were measured from TiO2 in the EELS Atlas) and solid
state effects (resulting from many-body perturbations by surrounding atomic
potentials). Nevertheless, the Dirac equation accurately predicts the energy
difference from splitting as observed experimentally. This trend is similarly
observed in the M4/M5 edge of Ag and N6/N7 edge of Au. Note that we did not
yet consider the effect of Zeeman splitting of the strong magnetic field in
the objective lens of a modern microscope, which could be interesting for
further study. To further verify the match between theory and experimental
measurements systematically, Fig. 2(a) plots the Dirac predicted ionization
energy against the edge onset energy all available in EELS Atlas, containing
1105 pairs of data in total and covering a wide range of atomic numbers (1-96)
and energies (11-34561 eV). As discrepancies between atomic calculations and
experimental measurement are expected due to the aforementioned chemical
shifts and solid-state effects, small deviations from the exact match (as
indicated by the dashed line) are observed in Fig. 2(b). A histogram of the
relative residual between the Dirac predictions and EELS Atlas values is
plotted inside of Fig. 2(b), showing the deviations (mean: 0.2%, standard
deviation: 7%) are well bounded. Further investigation shows that those cases
with large deviations are mainly from the O edges of heavy elements (i.e.
transition elements in the La and Ac families), which are also visible in Fig.
2(a). The general trend of residual can be captured fairly well with linear
regression as indicated by the red solid line. We can also make a reasonable
estimation for other edges not recorded in the EELS Atlas by applying the same
slope ($\sim$ 1.01) from the linear regression to the Dirac ionization energy.

Figure 2:  (a) Comparison of the Dirac predicted ionization energy with the
EELS Atlas for all available elements (Z=1-96). (b) A plot of the predicted
ionization energies against the edge onset energies from EELS Atlas to
demonstrate the linear correlation. The variation of atomic number can be
visualized by the color bar. The black dashed line indicates the exact match
between theory and experiment. The red solid line represents the linear
regression of the data. The probability histogram of the relative residual
between the Dirac calculation and EELS Atlas is plotted inside for the exact
match and after linear regression.

In addition to the precise onset energy, the Dirac solution is expected to
capture the accurate charge distribution due to relativity. The Dirac solution
has large and small components, where the large component is much more
significant and can be approximated to the Schrödinger solution at the low-
speed limit as shown in Eq. 5 and Eq. 36. However, this does not mean the
small component is not important, especially when the core-shell electron of
heavy elements approaches the speed of light. The expectation value of the
speed of an orbital electron $\langle v\rangle$ is defined as:

$\langle v\rangle=\frac{\langle p\rangle c^{2}}{\sqrt{m_{e}^{2}c^{4}+\langle
p^{2}\rangle c^{2}}}.$ (52)

Fig. 3(a) shows that orbital electrons (up to 5g orbital) $\langle v\rangle$
increase with the atomic number, reaching a considerable fraction of the speed
of light for heavy elements. For instance, the 1s electron of Au can reach
$\sim$ 57% of the speed of light. Ranking the speed among different orbitals,
one may observe that the 1s_1/2 electron is always associated with the highest
speed as it directly faces the nucleus (with the lowest expectation value of
the distance from nucleus $\langle r\rangle$ and therefore least screening
effects caused by other electrons). This is then followed by the 2s_1/2 and
2p_1/2 spin-down electrons with similar velocities. The 2p_1/2 electron has a
noticeably lower speed compared to that of 2p_3/2. There are additional
orbital pairs with similar speeds for 3s_1/2-3p_1/2, 3p_3/2-3d_3/2, etc. This
is not surprising as such pairs have the same total angular momentum $j$ and
relativistic quantum number term $(\kappa+1)\kappa$. Consequently, they have
the same effective potential and kinetic energy following Eq. 36 in the low
energy limit. For the Dirac solution, the fast core-shell electrons are
accompanied by the increasing contribution of the small component shown in
3(b). Here we define the contribution of the small component $\zeta$ for a
given orbital $n\kappa$ as:

$\zeta=\frac{\int_{0}^{\infty}\left(Q_{n\kappa}^{2}(r)\right)\mathrm{d}r}{\int_{0}^{\infty}\left(P_{n\kappa}^{2}(r)+Q_{n\kappa}^{2}(r)\right)\mathrm{d}r}=\int_{0}^{\infty}\left(Q_{n\kappa}^{2}(r)\right)\mathrm{d}r,$
(53)

since the wavefunctions are already normalized in Eq. 33. Fig. 3(c-e) shows
the radial charge density for three typical elements ranging from light to
heavy (Si, Ag, and Au). As expected, the large component always dominates the
charge distribution. In contrast, the small component is negligible for Si,
but becomes noticeable for Ag and significant for Au, mainly for core-shell
electrons close to the nucleus as indicated by the orbital decomposition
shaded in different colors. The direct consequence of the Dirac solution for
GOS is reflected in the radial integral of Eq. 40, where the small component
contributes to the transition matrix element. Further analysis of the L2/L3
edges shows that the contribution of small components to the GOS is negligible
for Si, approximately 1% for Ag, and 3-5% for Au (depending on the spin),
which is proportional to the contribution of the small components to the
charge distribution (Si: 0.03%, Ag 0.6%, Au 2%).

Figure 3: (a) The expectation value of relativistic speed $\langle v\rangle$
of the orbital electrons (up to 5g) for elements (Z=1-118) in the periodic
table. The spin-down electron is indicated by the dotted line. (b)
Contribution of the small component to the charge density for each orbital of
all elements. (c-e) The radial charge density for three typical elements from
light to heavy (Si, Ag, and Au). The large component is shown in blue and the
small component is shown in red. The orbital decomposition is indicated by the
shaded area in different colors. Figure 4: (a) Orbital energy level for Au as
a function of angular quantum number with spin variations. The red/blue shaded
line indicates spin up/down Dirac solution, while the green solid line
indicates the Schrödinger all-electron DFT solution. The radial wavefunctions
for (b) the 2p orbital, (c) the 3d orbital, and (d) continuum state d orbital
at 10 eV with spin variations.

The spin variation also yields different wavefunctions and hence different
transition matrixes and GOS, in contrast to a single matrix from the
Schrödinger solution. In Fig. 4(a), the electron ionization energy is plotted
as a function of the orbital and spin variations, wherein Fig. 4(b-d) plots
the radial wavefunctions for the bound and continuum states. The significance
of energy splitting and deviation from the Schrödinger solution is a direct
indication of the strength of SOC, which is particularly pronounced for p
orbitals and less obvious for higher orbitals (like d orbitals) in general due
to electron screening (note energies plotted on a log scale). For radial wave
functions of bound states, the 2p_1/2 orbital is contracted inward compared to
the 2p_3/2 orbital for both the large and small components as shown in Fig.
4(b). In contrast, the large components for 3d orbitals are similar regardless
of the spin, while its small component is slightly contracted inward for
3d_3/2 electron compared to 3d_5/2 as shown in Fig. 4(c). To calculate the
matrix elements, the bound state wavefunctions will be integrated with
different final continuum states for each energy, one of which is shown in
Fig. 4(d) for a final state angular quantum number $l^{\prime}=2$ with spin-
up/down at 10 eV above the ionization energy. The highly oscillating nature of
the continuum state (with frequency increases with increasing energy loss) and
spherical Bessel function (with frequency increases with increasing momentum
transfer) suggests that a slight shift in the radial wavefunction of the
highly localized bound state can result in a remarkable change in the integral
wavefunction overlaps for bound states with different spin quantum numbers.

To investigate how spin variation affects the GOS, Fig. 5 demonstrates the
contributions of final states with different angular quantum numbers and spin
variations, plotted as a function of energy loss and momentum transfer. Fig.
5(a-b) shows the decomposition of the contributions from final states for the
Au L2/L3 case at 10 eV and 4000 eV above the edge onset. At low energy loss,
the dipole transition $\Delta l^{\prime}=|l-l^{\prime}|=1$ dominates the GOS
as expected. However, for high energy loss, there are notable contributions
from final states with other $l^{\prime}$ leading to the emergence of a peak
at high momentum transfer known as the Bethe ridge [10]. In the case of Au,
the difference of spin reaches 20% for L2/L3 edges, but it can be higher for a
heavier element (like actinium can be over 40%). This difference in GOS is not
surprising, as the spin-down electron is closer to the nucleus with a higher
ionization energy due to SOC, its excitation is expected to have a lower
transition probability compared to the spin-up electron. In contrast, the
ratio of spin-down/up is close to 0.95 for M4/M5 edges which is a result of
negligible spin variations in their wavefunctions shown in Fig. 4(c).

Figure 5: The decomposition of the contributions from different final states
for the Au L2 and L3 case at (a) 10 eV and (b) 4000 eV above the ionization
energy and similarly for the Au M4 and M5 case in (d-e). The spin-down/up GOS
ratio against the energy loss and momentum transfer for (c) L2/L3 and (f)
M4/M5 edges.

###  6.2 Relativistic electrodynamics for EELS cross-section calculations

Figure 6: (a) Relativistic and non-relativistic double differential cross
sections for silicon K-ionization at 10 eV above the ionization energy with
the incident electron of 100 - 300 keV. (b) The ratio of relativistic to non-
relativistic differential double differential cross sections for silicon
K-ionization at 300 keV incident beam energy as a function of momentum
transfer and energy loss. The white dashed line indicates the minimum momentum
transfer.

The relativistic effects in electrodynamics significantly enhance the
scattering cross-section, particularly noticeable at low angles and high
energies. Fig. 6(a) shows the relativistic and non-relativistic double
differential cross sections for silicon K-edge at different acceleration
voltages, showing a small angle peak, which becomes more pronounced with
increasing acceleration voltage, negligible at 100 keV but not for higher
voltages. As explained in previous studies [54, 55, 56], this relativistic
correction is obtained when applying Lorentz gauge for a moving charge in Eq.
37-38. It also leads to the interference of the longitudinal and transverse
components of the excitation, resulting in anisotropic differential cross-
sections which is critical when measuring near-edge fine structures [55]. Fig.
6(b) shows the ratio of relativistic to non-relativistic differential cross
sections for silicon K-ionization at 300 keV incident beam energy as a
function of momentum transfer and energy loss, in contrast to the single
energy loss in Fig. 6(a), which agrees with previous results [54, 56]. The
correction is significant at the small scattering angles and varies slowly
against increasing energy loss. A white dashed line in Fig. 6(b) indicates the
minimum momentum transfer. To appreciate this effect on the spectra under
conventional experimental conditions, we can integrate DDSCS the scattering
within the EELS collection aperture via Eq. 44 to differential scattering
cross-section. To check we have the correct implementation of the relativistic
correction, Fig. 7 plots the relativistic correction ratio of the Si K-edge at
different acceleration voltages with (a) 10 mrad and (b) 100 mrad collection
angles using numerical calculation in Eq. 38 and analytical expression with
dipole approximation in Eq. 42. The numerical calculation of the relativistic
correction ratio is consistent with the analytical dipole solution in Eq. 42
for small angles of 10 mrad, yielding up to a 15% increase in the cross-
sections at 300 keV. However, with increasing collection angle, the
relativistic double differential cross-sections are very close to the non-
relativistic ones (i.e. only about 5% for 300 keV and 100 mrad) as shown in
Fig. 6. It is important to note that the dipole approximation is expected to
fail at high angles and one has to use numerical relativistic corrections.
Overall, we want to emphasize that the relativistic correction is important
for all elements (in contrast to the relativistic atomic orbital effects,
which are only obvious for core-shell electrons of heavy elements). Since this
relativistic effect varies with energy loss and momentum transfer, the
elemental quantification will be affected differently for different edges
under varying experimental conditions. These corrections should be routinely
performed when quantifying EELS with typical TEM acceleration voltages at
200-300 keV, particularly when using a small collection aperture.

Figure 7: The relativistic differential scattering cross-section (DSCS)
correction ratio of the Si K edge against the energy loss at different
acceleration voltages (100-300 keV) with (a) 10 mrad and (b) 100 mrad
collection angles using the full numerical solution and the analytical
solution with dipole approximation.

###  6.3 Comparison of the Dirac-based GOS database with others

Figure 8: Comparison of the GOS databases for (a) O K, (b) Ti L2, (c) Sr L3,
(d) La M4, (e) Bi M5 and (f) U O5 edges using different GOS databases with 300
keV incident beam, parallel illumination, and 50 mrad collection angle.

We compared the Dirac-based GOS with the existing GOS database used in the
EELS community, including the HF solutions, and DFT solutions based on the
Schrödinger equation. For a fair comparison evaluation of the GOS databases
themselves, the relativistic electrodynamics effect in Sec. 6.2 is not
included in this comparison. In Fig. 8, the calculated spectra for O K, Ti L2,
Sr L3, La M4, Bi M5, and U O5 edges using different GOS databases with 300 keV
incident beam energy, parallel illumination, and a 50 mrad collection angle.
Overall, these spectra exhibit similar shapes and absolute amplitude except
for a few details. Specifically, the HF spectra are more fragmented due to the
limited sampling in both energy and momentum, as shown in Fig. 1. The
Schrödinger cross-sections are in good agreement with the Dirac ones for light
elements where SOC is expected to be weak. However, notable differences are
observed in the heavy elements for the Bi M5 and U O5 edges. Interestingly for
the U O5 edge, HF and Dirac curves are in agreement while the DFT curve
deviates away at the near edge region. Also, we could see a sharp peak at the
start of the DFT Ti L2 edge. Since those GOS are computed using different
internal numerical routines, it is difficult at this stage to pinpoint whether
the discrepancy is attributed to numerical issues or the underlying physics.
It could be useful to calculate the Schrödinger and Dirac solutions and
construct the GOS under the same protocols to trace the causes. We noticed
that there are packages (e.g. RIDIAL[66] and dftatom [73]) that can compute
both the Schrödinger and Dirac atomic radial wavefunctions, which will enable
the direct comparisons.

We expect that the inner core-shell edges (i.e. K and L-edge) of heavy
elements will offer a better opportunity to observe the relativistic effects,
although this corresponds to extremely high energy loss that is typically not
explored often in EELS. With the development of direct electron cameras and
advanced spectrometers, experimental verification of the Dirac cross-sections
can be expected in the future.

###  6.4 Use of the Dirac-based GOS database

The Dirac-based GOS database is made available open-source on Zenodo under a
CC-BY license [74]. Similar to using all GOS databases, quantifying the EELS
cross-sections involves a three-step process. Step 1, we need to compute the
DDSCS from the GOS database given in Eq. LABEL:eq:ddscs for the non-
relativistic scattering but with atomic orbitals obtained from the Dirac
equation. Then we use Eq. 38 to account for the relativistic nature of the
fast incoming electron. Step 2, we need to integrate the DDSCS within the EELS
collection aperture to obtain the differential cross-sections in Eq. 44 for
parallel illumination. If convergent beam illumination is used in STEM-EELS,
we have to correct it with a geometric cross-correlation function in Eq. 49
[67]. The momentum space integration for both parallel and convergent beam
illumination is implemented in pyEELSModel. Step 3, the computed differential
cross-sections for all edges are fitted together to the experimental spectrum
for quantification. During the model-based fitting process, the core-loss
spectrum can be convoluted with the low-loss spectrum to consider the source
energy dispersion and plural plasmon scattering. The fitting procedure can be
performed using the existing packages such as EELSModel, HyperSpy, or Gatan
Digital Micrograph, which now all support the Dirac-based GOS. We also provide
a Jupyter notebook [75] to demonstrate the fitting process of the Dirac-based
GOS for the Ti L2/L3 and O K edges in SrTiO3 using pyEELSModel [13]. We
encourage inquiries from the community and industry to integrate the Dirac-
based GOS into their software packages.

##  7 Conclusions and outlook

EELS is a powerful technique for the characterization of materials composition
and electronic structures. The accuracy of the EELS quantification relies on
the precise calculation of the cross-sections facilitated by the GOS database,
which computes the transition probability as a function of energy loss and
momentum transfer. In this study, we developed an open-source Dirac-based GOS
database for all elements and shells for large energy and momentum space with
fine sampling. Unlike previous GOS databases relying on the Schrödinger
equation, the Dirac-based GOS database better captures the relativistic
effects. We demonstrated that the ionization energies calculated from the
Dirac equation align well with the experimental measurements with a relative
error of $0.2\%\pm 7\%$, including the energy splitting caused by SOC.
Notably, the large component of the Dirac radial wave is shown to be
equivalent to the non-relativistic Schrödinger solution at low energies, while
the small component of the Dirac solution becomes significant for the core-
shell electrons of heavy elements. We also showed that SOC leads to different
wavefunctions and hence different transition matrix elements and GOS for spin
variations. We compared the Dirac-based GOS database with existing GOS
databases commonly used in the EELS community. Despite fine differences
usually at the near-edge region, these spectra are generally similar for the K
and L edges of light elements or the higher edges of heavy elements. But for
the core-shell electrons of heavy elements, theoretical results suggest
considerable difference for spin up/down (up to 20% for Au and 40% for Ac),
which is rarely measured experimentally by EELS. Further experimental work is
needed to verify the Dirac-based GOS for absolute cross-sections of the heavy
elements at high energy loss, which requires an advanced EELS spectrometer and
detector capable of this precise measurement. Overall, the Dirac-based GOS
database is expected to be useful for microanalysis with relativistic effects
considered across the periodic table. In addition, the differential cross-
sections could be also useful for the Monte Carlo simulations of electron
scattering.

The calculations in this work are dedicated to the GOS database for EELS
microanalysis. However, this computational framework can be easily extended to
generate the oscillator strength database for XAS. As the momentum transfer is
negligible for photons, the optical dipole oscillator strength is not
dependent on the scattering angle. This can be obtained by substitution of the
Hamiltonian in the transition matrix of Fermi’s golden rule, which is the
Coulombic interaction for an electron, with the position operator and the
polarization vector of the X-ray beam. Indeed, EELS and XAS are often used in
tandem to provide complementary information on materials. A pair of EELS GOS
and XAS oscillator strength databases should be of interest to material
scientists who use both techniques.

In analytical electron microscopy for microanalysis, the most common
techniques are energy-dispersive X-ray spectroscopy (EDS) and EELS. The EDS
collects the characteristic X-ray generated during the de-excitation of the
core-shell electron, which is a physically coupled process and therefore its
theoretical cross-sections are deeply related to the EELS calculations
presented here. Experimentally, EDS elemental mapping is much simpler than
EELS and favored by material scientists. However, accurate and precise EDS
quantification is not easy, particularly at the atomic scale, due to various
theoretical and experimental challenges. We have generated the Dirac-based
ionization potential database for atomic EDS calculations in the Bloch wave or
multislice algorithms, which includes the relativistic effects from both the
orbital and fast incident electron. The Dirac-based EDS database will be made
available to the community in the future.

We notice that the ionization energy depends on the chemical environment,
which is also expected to change the cross-sections as well. Thus, the
calculation of the GOS for different oxidation states is currently under
development. Moreover, our current analysis, limited to continuum states of
individual atoms, fails to encapsulate the intricate unoccupied DOS inherent
in complex, multi-atom systems, thereby neglecting the fine structures in the
near-edge region. Given that the near-edge region often contains most of the
signal counts against background and noise, the exclusion of physics-based
fine structure from spectrum fitting represents a significant loss of
information. Hence, we are strongly motivated to incorporate the unoccupied
bounded states in our future generation of theoretical EELS cross-sections.

##  8 Data availability

The database is available on Zenodo under a CC-BY license [74].

## Acknowledgement

The authors acknowledge financial support from the Research Foundation
Flanders (FWO, Belgium) through Project No.G.0502.18N. This project has
received funding from the European Research Council (ERC) under the European
Union’s Horizon 2020 research and innovation programme (Grant Agreement No.
770887 PICOMETRICS and No. 823717 ESTEEM3). The authors acknowledge the
computational resources from Vlaams Supercomputer Centrum (VSC) through Tier-1
project No.2023.023. ZZ acknowledges the travel funding from FWO (Grant
No.V463823N) and the consultation with Dirk Van Dyck, Peter Schattschneider,
Christian Dwyer, and Christian Rossouw. ZZ thanks Linfeng Zhang and Chu-ping
Yu for reading the manuscript and offering helpful comments.

## References

  * [1] R. F. Egerton, Electron energy-loss spectroscopy in the electron microscope, Springer Science & Business Media, 2011.
  * [2] H. H. Rose, Historical aspects of aberration correction, Journal of electron microscopy 58 (3) (2009) 77–85.
  * [3] O. L. Krivanek, T. C. Lovejoy, M. F. Murfitt, G. Skone, P. E. Batson, N. Dellby, Towards sub-10 mev energy resolution stem-eels, in: Journal of Physics: Conference Series, Vol. 522, IOP Publishing, 2014, p. 012023.
  * [4] O. L. Krivanek, T. C. Lovejoy, N. Dellby, T. Aoki, R. Carpenter, P. Rez, E. Soignard, J. Zhu, P. E. Batson, M. J. Lagos, et al., Vibrational spectroscopy in the electron microscope, Nature 514 (7521) (2014) 209–212.
  * [5] J. L. Hart, A. C. Lang, A. C. Leff, P. Longo, C. Trevor, R. D. Twesten, M. L. Taheri, Direct Detection Electron Energy-Loss Spectroscopy: A Method to Push the Limits of Resolution and Sensitivity, Scientific Reports 7 (1) (2017) 8243.  doi:10.1038/s41598-017-07709-4.
URL https://www.nature.com/articles/s41598-017-07709-4

  * [6] B. Plotkin-Swing, G. J. Corbin, S. De Carlo, N. Dellby, C. Hoermann, M. V. Hoffman, T. C. Lovejoy, C. E. Meyer, A. Mittelberger, R. Pantelic, et al., Hybrid pixel direct detector for electron energy loss spectroscopy, Ultramicroscopy 217 (2020) 113067.
  * [7] D. A. Muller, L. F. Kourkoutis, M. Murfitt, J. H. Song, H. Y. Hwang, J. Silcox, N. Dellby, O. L. Krivanek, Atomic-scale chemical imaging of composition and bonding by aberration-corrected microscopy, Microscopy and Microanalysis 14 (SUPPL. 2) (2008) 132–133.  doi:10.1017/S1431927608087448.
  * [8] J. A. Hachtel, J. Huang, I. Popovs, S. Jansone-Popova, J. K. Keum, J. Jakowski, T. C. Lovejoy, N. Dellby, O. L. Krivanek, J. C. Idrobo, Identification of site-specific isotopic labels by vibrational spectroscopy in the electron microscope, Science 363 (6426) (2019) 525–528.
  * [9] R. Egerton, M. Malac, Improved background-fitting algorithms for ionization edges in electron energy-loss spectra, Ultramicroscopy 92 (2) (2002) 47–56.
  * [10] Egerton, EELS textbook, 2011, publication Title: Electron spectroscopy for surface analysis.
URL
http://link.springer.com/chapter/10.1007/978-3-642-81099-2_6%0Ahttp://link.springer.com/chapter/10.1007/978-3-642-81099-2%7B_%7D6

  * [11] J. Verbeeck, S. Van Aert, Model based quantification of eels spectra, Ultramicroscopy 101 (2-4) (2004) 207–224.
  * [12] R. Leapman, C. Swyt, Separation of overlapping core edges in electron energy loss spectra by multiple-least-squares fitting, Ultramicroscopy 26 (4) (1988) 393–403.
  * [13] D. Jannis, joverbee/pyeelsmodel: v0.1.1 (Apr. 2024).  doi:10.5281/zenodo.10992986.
URL https://doi.org/10.5281/zenodo.10992986

  * [14] F. de la Peña, T. Ostasevicius, V. T. Fauske, P. Burdet, P. Jokubauskas, M. Nord, M. Sarahan, E. Prestat, D. N. Johnstone, J. Taillon, et al., Electron microscopy (big and small) data analysis with the open source software package hyperspy, Microscopy and Microanalysis 23 (S1) (2017) 214–215.
  * [15] P. Cueva, R. Hovden, J. A. Mundy, H. L. Xin, D. A. Muller, Data processing for atomic resolution electron energy loss spectroscopy, Microscopy and Microanalysis 18 (4) (2012) 667–675.
  * [16] W. Van den Broek, D. Jannis, J. Verbeeck, Convexity constraints on linear background models for electron energy-loss spectra, Ultramicroscopy 254 (2023) 113830.
  * [17] J. Verbeeck, G. Bertoni, Model-based quantification of EELS spectra: Treating the effect of correlated noise, Ultramicroscopy 108 (2) (2008) 74–83.  doi:10.1016/j.ultramic.2007.03.004.
  * [18] J. Verbeeck, S. Van Aert, G. Bertoni, Model-based quantification of EELS spectra: Including the fine structure, Ultramicroscopy 106 (11-12 SPEC. ISS.) (2006) 976–980.  doi:10.1016/j.ultramic.2006.05.006.
  * [19] P. Schattschneider, Fundamentals of inelastic electron scattering, Springer Science & Business Media, 1986.
  * [20] P. A. M. Dirac, The quantum theory of the emission and absorption of radiation, Proceedings of the Royal Society of London. Series A, Containing Papers of a Mathematical and Physical Character 114 (767) (1927) 243–265.
  * [21] P. A. M. Dirac, The quantum theory of the electron, Proceedings of the Royal Society of London. Series A, Containing Papers of a Mathematical and Physical Character 117 (778) (1928) 610–624.
  * [22] E. Fermi, Nuclear physics: a course given by Enrico Fermi at the University of Chicago, University of Chicago press, 1950.
  * [23] H. Bethe, Bremsformel für elektronen relativistischer geschwindigkeit, Zeitschrift für Physik 76 (5-6) (1932) 293–299.
  * [24] B. L. Moiseiwitsch, S. J. Smith, Electron impact excitation of atoms, Reviews of Modern Physics 40 (2) (1968) 238.
  * [25] E. J. McGuire, Inelastic scattering of electrons and protons by the elements he to na, Physical Review A 3 (1) (1971) 267.
  * [26] S. T. Manson, Inelastic Collisions of Fast Charged Particles with Atoms: Ionization of the Aluminum L Shell, Physical Review A - Atomic, Molecular, and Optical Physics 6 (3) (1972) 1013.
  * [27] J. H. Scofield, K-and l-shell ionization of atoms by relativistic electrons, Physical Review A 18 (3) (1978) 963.
  * [28] R. D. Leapman, P. Rez, D. F. Mayers, K, L, and M shell generalized oscillator strengths and ionization cross sections for fast electron collisions, The Journal of Chemical Physics 72 (2) (1980) 1232–1243.  doi:10.1063/1.439184.
  * [29] D. K. Saldin, P. Rez, The theory of the excitation of atomic inner-shells in crystals by fast electrons, Philosophical Magazine B: Physics of Condensed Matter; Statistical Mechanics, Electronic, Optical and Magnetic Properties 55 (4) (1987) 481–489.  doi:10.1080/13642818708217957.
  * [30] P. Rez, Inner-shell spectroscopy: An atomic view, Ultramicroscopy 28 (1-4) (1989) 16–23.  doi:10.1016/0304-3991(89)90263-5.
  * [31] S. Segui, M. Dingfelder, J. M. Fernández-Varea, F. Salvat, The structure of the Bethe ridge. Relativistic Born and impulse approximations, Journal of Physics B: Atomic, Molecular and Optical Physics 35 (1) (2002) 33–53.  doi:10.1088/0953-4075/35/1/303.
  * [32] D. Bote, F. Salvat, Calculations of inner-shell ionization by electron impact with the distorted-wave and plane-wave Born approximations, Physical Review A 77 (4) (2008) 042701.  doi:10.1103/PhysRevA.77.042701.
URL https://link.aps.org/doi/10.1103/PhysRevA.77.042701

  * [33] H. Kohl, H. Rose, Theory of image formation by inelastically scattered electrons in the electron microscope, in: Advances in Electronics and Electron Physics, Vol. 65, Elsevier, 1985, pp. 173–227.
  * [34] P. Stallknecht, H. Kohl, Computation and interpretation of contrast in crystal lattice images formed by inelastically scattered electrons in a transmission electron microscope, Ultramicroscopy 66 (3-4) (1996) 261–275.  doi:10.1016/S0304-3991(97)00007-7.
  * [35] Z. L. Wang, Dynamical inelastic scattering in high-energy electron diffraction and imaging: A new theoretical approach, Physical Review B 41 (18) (1990) 12819.  doi:10.1016/0003-4916(58)90036-8.
  * [36] C. J. Rossouw, V. M. Maslen, Implications of (e, 2e) scattering for inelastic electron diffraction in crystals II. Application of the theory, Philosophical Magazine A: Physics of Condensed Matter, Structure, Defects and Mechanical Properties 49 (6) (1984) 743–757.  doi:10.1080/01418618408236561.
  * [37] D. A. Muller, J. Silcox, Delocalization in inelastic scattering, Ultramicroscopy 59 (1-4) (1995) 195–213.  doi:10.1016/0304-3991(95)00029-Z.
  * [38] E. C. Cosgriff, M. P. Oxley, L. J. Allen, S. J. Pennycook, The spatial resolution of imaging using core-loss spectroscopy in the scanning transmission electron microscope, Ultramicroscopy 102 (4) (2005) 317–326.  doi:10.1016/j.ultramic.2004.11.001.
  * [39] L. J. Allen, A. J. D’Alfonso, S. D. Findlay, M. P. Oxley, M. Bosman, V. J. Keast, E. C. Cosgriff, G. Behan, P. D. Nellist, A. I. Kirkland, Theoretical interpretation of electron energy-loss spectroscopic images, AIP Conference Proceedings 999 (2008) (2008) 32–46, iSBN: 9780735405196.  doi:10.1063/1.2918115.
  * [40] M. P. Oxley, S. J. Pennycook, Image simulation for electron energy loss spectroscopy, Micron 39 (6) (2008) 676–684.  doi:10.1016/j.micron.2007.10.020.
  * [41] A. J. D’Alfonso, S. D. Findlay, M. P. Oxley, L. J. Allen, Volcano structure in atomic resolution core-loss images, Ultramicroscopy 108 (7) (2008) 677–687.  doi:10.1016/j.ultramic.2007.10.011.
  * [42] J. Rusz, J. Spiegelberg, S. Muto, T. Thersleff, M. Ohtsuka, K. Leifer, P. M. Oppeneer, Localization of magnetic circular dichroic spectra in transmission electron microscopy experiments with atomic plane resolution, Physical Review B 95 (17) (2017) 1–12.  doi:10.1103/PhysRevB.95.174412.
  * [43] H. Yoshioka, Effect of Inelastic Waves on Electron Diffraction, Journal of the Physical Society of Japan 12 (6) (1957) 618–628.  doi:10.1143/JPSJ.12.618.
  * [44] V. W. Maslen, C. J. Rossouw, The inelastic scattering matrix element and its application to electron energy loss spectroscopy, Philosophical Magazine A: Physics of Condensed Matter, Structure, Defects and Mechanical Properties 47 (1) (1983) 119–130.  doi:10.1080/01418618308243112.
  * [45] V. W. Maslkn, C. J. Rossouw, Implications of (e, 2e) scattering for inelastic electron diffraction in crystals I. Theoretical, Philosophical Magazine A: Physics of Condensed Matter, Structure, Defects and Mechanical Properties 49 (6) (1984) 735–742.  doi:10.1080/01418618408236560.
  * [46] L. J. Allen, S. D. Findlay, M. P. Oxley, C. J. Rossouw, Lattice-resolution contrast from a focused coherent electron probe. Part I, Ultramicroscopy 96 (1) (2003) 47–63.  doi:10.1016/S0304-3991(02)00380-7.
  * [47] C. Dwyer, Multislice theory of fast electron scattering incorporating atomic inner-shell ionization, Ultramicroscopy 104 (2) (2005) 141–151.  doi:10.1016/j.ultramic.2005.03.005.
  * [48] C. Dwyer, S. D. Findlay, L. J. Allen, Multiple elastic scattering of core-loss electrons in atomic resolution imaging, Physical Review B - Condensed Matter and Materials Physics 77 (18) (2008) 1–8.  doi:10.1103/PhysRevB.77.184107.
  * [49] I. Lobato, S. van Aert, J. Verbeeck, Progress and new advances in simulating electron microscopy datasets using MULTEM, Ultramicroscopy 168 (2016) (2016) 17–27, publisher: Elsevier.  doi:10.1016/j.ultramic.2016.06.003.
URL http://dx.doi.org/10.1016/j.ultramic.2016.06.003

  * [50] S. Majert, H. Kohl, Simulation of atomically resolved elemental maps with a multislice algorithm for relativistic electrons, in: Advances in Imaging and Electron Physics, Vol. 211, Elsevier, 2019, pp. 1–120.  doi:10.1016/bs.aiep.2019.04.001.
URL https://linkinghub.elsevier.com/retrieve/pii/S1076567019300266

  * [51] H. G. Brown, J. Ciston, C. Ophus, Linear-scaling algorithm for rapid computation of inelastic transitions in the presence of multiple electron scattering, Physical Review Research 1 (3) (2019) 33186, arXiv: 1908.07628 Publisher: American Physical Society.  doi:10.1103/physrevresearch.1.033186.
URL https://doi.org/10.1103/PhysRevResearch.1.033186

  * [52] U. Fano, Differential Inelastic Scattering of Relativistic Charged Particles, Physical Review 102 (2) (1956) 385–387.  doi:10.1103/PhysRev.102.385.
URL https://link.aps.org/doi/10.1103/PhysRev.102.385

  * [53] H. Kurata, P. Wahlbring, S. Isoda, T. Kobayashi, Importance of relativistic effect on inelastic scattering cross-sections for quantitative microanalysis, Micron 28 (5) (1997) 381–388.  doi:10.1016/S0968-4328(97)00044-9.
URL https://linkinghub.elsevier.com/retrieve/pii/S0968432897000449

  * [54] R. Knippelmeyer, P. Wahlbring, H. Kohl, Relativistic ionisation cross sections for use in microanalysis, Ultramicroscopy 68 (1) (1997) 25–41.
  * [55] P. Schattschneider, C. Hébert, H. Franco, B. Jouffrey, Anisotropic relativistic cross sections for inelastic electron scattering, and the magic angle, Physical Review B - Condensed Matter and Materials Physics 72 (4) (2005).  doi:10.1103/PhysRevB.72.045142.
  * [56] C. Dwyer, J. S. Barnard, Relativistic effects in core-loss electron diffraction, Physical Review B - Condensed Matter and Materials Physics 74 (6) (2006) 28–32.  doi:10.1103/PhysRevB.74.064106.
  * [57] K. Jorissen, J. Rehr, J. Verbeeck, Multiple scattering calculations of relativistic electron energy loss spectra, Physical Review B 81 (15) (2010) 155108.
  * [58] J. J. Rehr, J. J. Kas, M. P. Prange, A. P. Sorini, Y. Takimoto, F. Vila, Ab initio theory and calculations of x-ray spectra, Comptes Rendus Physique 10 (6) (2009) 548–559.
  * [59] L. Segger, Berechnung generalisierter Oszillatorenstärken für die Quantifizierung von EEL-Spektren, Bachelor’s thesis, WWU-Münster (2019).
  * [60] L. Segger, G. Guzzinati, H. Kohl, Generalised oscillator strengths for the simulation of eels spectra, with a broader coverage of high energy and minor edges (2023).  doi:10.5281/zenodo.7645765.
URL https://doi.org/10.5281/zenodo.7645765

  * [61] D. R. Hamann, Generalized norm-conserving pseudopotentials, Physical Review B 40 (5) (1989) 2980.
  * [62] G. Guzzina, GOSH: Generalised Oscillator Strengths in HDF5, gitLab repository (2023).
URL https://gitlab.com/gguzzina/gosh

  * [63] R. D. Cowan, The theory of atomic structure and spectra, University of California Press, 1981, tex.date-added: 2021-06-01 21:26:57 +0800 tex.date-modified: 2021-06-01 21:26:57 +0800.
  * [64] F. W. Olver, NIST handbook of mathematical functions, Cambridge university press, 2010.
URL "https://dlmf.nist.gov/"

  * [65] W. Greiner, Relativistic quantum mechanics, Springer Berlin, Heidelberg, 2000.
  * [66] F. Salvat, J. M. Fernández-Varea, RADIAL: A Fortran subroutine package for the solution of the radial Schrödinger and Dirac wave equations, Computer Physics Communications 240 (2019) 165–177, publisher: Elsevier B.V.  doi:10.1016/j.cpc.2019.02.011.
URL https://doi.org/10.1016/j.cpc.2019.02.011

  * [67] H. Kohl, A simple procedure for evaluating effective scattering cross-sections in stem, Ultramicroscopy 16 (2) (1985) 265–268.
  * [68] M. F. Gu, The flexible atomic code, Canadian Journal of Physics 86 (5) (2008) 675–689.  doi:10.1139/P07-197.
  * [69] P. E. Blöchl, C. J. Först, J. Schimpl, Projector augmented wave method: ab initio molecular dynamics with full wave functions, Bulletin of Materials Science 26 (2003) 33–41.
  * [70] J. P. Perdew, K. Burke, M. Ernzerhof, Generalized Gradient Approximation Made Simple, Physical Review Letters 77 (18) (1996) 3865–3868, arXiv: 10.1016/0927-0256(96)00008 ISBN: 9780596529321.  doi:10.1103/PhysRevLett.77.3865.
URL https://link.aps.org/doi/10.1103/PhysRevLett.77.3865

  * [71] C. C. Ahn, Transmission electron energy loss spectrometry in materials science and the EELS atlas, John Wiley & Sons, 2006.
  * [72] A. Y. Lee, D. M. Blakeslee, C. J. Powell, J. R. Rumble Jr, Development of the web-based nist x-ray photoelectron spectroscopy (xps) database, Data Science Journal 1 (2002) 1–12.
  * [73] O. Čertík, J. E. Pask, J. Vackář, dftatom: A robust and general schrödinger and dirac solver for atomic structure calculations, Computer Physics Communications 184 (7) (2013) 1777–1791.
  * [74] Z. Zhang, I. Lobato, D. Jannis, J. Verbeeck, S. Van Aert, P. Nellist, Generalised oscillator strength for core-shell electron excitation by fast electrons based on dirac solutions (2024).  doi:10.5281/zenodo.11199911.
URL https://doi.org/10.5281/zenodo.11199911

  * [75] Z. Zhang, Gos demo: Tutorial on how to use dirac gos for eels quantification (2024).
URL https://github.com/zezhong-zhang/gos_demo