isotropic and infinite medium (we can omit in that case the dependence on the
momentum vector and retain only that of $|{\bm{q}}|=q$), Eqs. (127) and (128)
can be integrated numerically or solved analytically as follows [61, 104,
140]:

$\varepsilon^{0}(q,\omega)\,=\,1+\frac{\chi^{2}}{z^{2}}[f_{1}(u,z)\,+\,i\,f_{2}(u,z)]\,,$
(129)

where $u=\omega/(qv_{\mathrm{F}})$, $z=q/(2q_{\mathrm{F}})$,
$\chi^{2}=e^{2}/(\pi\,\hbar\,v_{\mathrm{F}})$, $v_{\mathrm{F}}$ is the Fermi
velocity and $q_{\mathrm{F}}=mv_{\mathrm{F}}/\hbar$. The functions
$f_{1}(u,z)$ and $f_{2}(u,z)$ are given by

$f_{1}(u,z)\,=\,\frac{1}{2}\,+\,\frac{1}{8z}\,[g(z-u)\,+\,g(z+u)]\,,$ (130)

$f_{2}(u,z)=\left\\{\begin{array}[]{lll}&&\frac{\pi}{2}u\,,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,z+u<1\\\
&&\frac{\pi}{8z}[1-(z-u)^{2}]\,,\,\,\,\,\,\,|z-u|<1<z+u\\\
&&0\,,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,|z-u|>1\,,\end{array}\right.$
(131)

where

$g(x)\,=\,(1-x^{2})\,\ln\left|\frac{1+x}{1-x}\right|\,.$ (132)

####  5.5.2 The Drude-Lorentz and the Generalized Oscillator Strength methods

The Drude–Lorentz (DL) [44] and the Mermin Energy Loss Function Generalized
Oscillator Strength (MELF-GOS) [60, 61, 140, 142, 143] represent other
effective and accurate numerical approaches to calculate the ELF over the
entire momentum and energy surface. They include valence and inner-shell
electron excitations in the model ELF.

In both methods, the optical or long-wavelength (${\bm{q}}\rightarrow 0$)
limit of the ELF of the target material, typically obtained by optical
experiments, is fitted to an analytical form given for the valence (or outer)
electrons by a linear combination of DL functions

$\rm{Im}\left[\frac{1}{\varepsilon_{\rm{M}}(A_{i},\omega_{i},\gamma_{i};{\bm{q}}=0,\omega)}\right]=\frac{A_{i}\gamma_{i}\omega}{(\omega_{i}^{2}-\omega^{2})^{2}+(\gamma_{i}\omega)^{2}},$
(133)

as follows:

$\displaystyle\rm{Im}\left[\frac{1}{\bar{\varepsilon}({\bm{q}}=0,\omega)}\right]_{\rm
outer}=$
$\displaystyle\sum_{i}F(\omega-\omega_{\mathrm{th},i})\rm{Im}\left[\frac{1}{\varepsilon_{\rm
M}(A_{i},\omega_{i},\gamma_{i};{{\bm{q}}=0},\omega)}\right],$ (134)

where

$F(\omega-\omega_{\mathrm{th},i})=\frac{1}{1+e^{-\Delta_{i}(\omega-\omega_{\mathrm{th},i})}}.$
(135)

In Eq. (5.5.2) $\bar{\varepsilon}({{\bm{q}}=0},\omega)$ is the macroscopic
optical dielectric function (see discussion in section 5.4). $A_{i}$,
$\omega_{i}$, $\gamma_{i}$ and $\Delta_{i}$ are fitting parameters that are
determined by a best fit. They correspond to the relative weight, position and
width of the peaks observed in the optical ELF, while $\Delta_{i}$ is added to
more accurately reproduce the onset of electronic excitation from the inner
shells.

We note that the ELF obtained by using the Mermin dielectric function in Eq.
(126) reduces to the DL ELF of Eq. (133) in the optical limit, so that the
Mermin and DL oscillators provide the same values for the fitting parameters
in the optical limit [144].

Inner shell electrons are included in the representation of the ELF using
atomic generalised oscillator strengths (GOS) as follows [60]:

$\mathrm{Im}\left[\frac{1}{\bar{\varepsilon}({\bm{q}}=0,\omega)}\right]_{\mathrm{inner}}=\frac{2\pi^{2}{N}}{\omega}\sum_{j}\alpha_{j}\sum_{nl}\frac{df_{nl}^{j}({\bm{q}},\omega)}{d\omega}\Theta(\omega-\omega_{{\rm
th},nl}^{j}),$ (136)

where $\frac{df_{nl}^{j}({\bm{q}},W)}{dW}$ are the atomic GOS hydrogen wave
functions obtained using an effective nuclear charge for each inner orbital
$(n,l)$ of the $j$th atomic constituent of the target material with
stoichiometric weight $\alpha_{j}$ and ionisation energy $\omega_{{\rm
th},nl}^{j}$. In Eq. (136), ${N}$ is the atomic or molecular number density of
the target material.

Finally, a third way to determine the ELF is the Penn model [59, 145]. In this
approach, the ELF is represented by a convolution of the imaginary part of the
inverse Lindhard dielectric function with a spectral density function obtained
from experimental optical data.

The accuracy of the fit in both the DL and MELF-GOS methods can be verified by
checking the fulfilment (i) of the $f$-sum rule, according to which the
integral of the ELF multiplied by the energy loss gives an effective number of
electrons per atom (typically the number of valence or total electrons):

$Z_{\mathrm{eff}}(\omega_{\mathrm{max}})=\frac{2}{\pi\omega_{p}^{2}}\int_{0}^{\omega_{\mathrm{max}}}\omega\times\mathrm{Im}\left[\frac{1}{\bar{\varepsilon}(\omega,q\rightarrow
0)}\right]\,d\omega\,,$ (137)

(please note that $Z_{\mathrm{eff}}(\omega_{\mathrm{max}})\rightarrow Z$ is to
be understood as $\omega_{\mathrm{max}}\rightarrow\infty$, where $Z$ is the
target atomic number); ii) of the perfect-screening sum-rule, where the
integral:

$P_{\mathrm{eff}}(\omega_{\mathrm{max}})=\frac{2}{\pi}\int_{0}^{\omega_{\mathrm{max}}}\frac{1}{\omega}\times\mathrm{Im}\left[\frac{1}{\bar{\varepsilon}(\omega,q\rightarrow
0)}\right]\,d\omega,$ (138)

is an increasing monotonic function that asymptotically goes to 1 when
$\omega_{\mathrm{max}}\rightarrow\infty$.

Once the best fit parameters have been determined using the experimental data
for the optical ELF and carefully checked using the sum rules, the DL and
MELF-GOS methods differ in the way they deal with the extension to finite
momentum transfer necessary to account for the dispersion of the excitation
spectrum in the evaluation of energy losses (see Eq. (5.3.2)).

In the context of the MELF-GOS, it is assumed that the energy-momentum
dispersion law is given by the Mermin function of Eqs. (127) and (128).

Conversely, the extension of the ELF to the finite momentum transfer in the DL
approach is achieved by introducing a functional dependence of the energy on
the modulus of the electron momentum $\omega=\omega({\bm{q}})=\omega(q)$ in
Eq. (133) [60, 146], assuming, of course, that the solid is homogeneous and
isotropic. In these conditions, the most general dispersion law that can be
derived considering the correct limit forms for the momentum transfer
$q\rightarrow 0$ (Mermin or Drude formula) and $q\rightarrow\infty$ (single-
particle dispersion) is the following [146]:

$\hbar\omega_{n}(q)=\hbar\omega_{n}+\alpha_{\mathrm{RPA}}\frac{\hbar^{2}q^{2}}{2m},$
(139)

where $\omega_{n}$ is the $n$ peak plasmon frequency for $q=0$,
$\alpha_{\mathrm{RPA}}=6E_{\mathrm{F}}/(5\omega_{P})$ ($E_{\mathrm{F}}$ is the
Fermi energy of the material), $\omega_{P}=\sqrt{4\pi\rho e^{2}/m}$ is the
nominal bulk plasmon frequency, where $e$, $\rho$, $m$ are the electron
charge, density and mass, respectively. We note that for $q\rightarrow 0$ this
expression approaches the value $\omega_{n}$, which corresponds to the plasma
frequency in a metallic bulk system; for medium $q$ it approaches the Lindhard
dispersion relation of the electron gas; for high $q$ the energy-momentum
dispersion resembles the quadratic dispersion of the random-phase
approximation (RPA), which corresponds to the $n$-independent single-particle
excitation for a 3D electron gas. In the most general case, the damping
coefficients $\gamma_{i}$ in the DL relation (133) can also be varied over the
$q$-space.

For semi-infinite or finite systems, such as slabs, the occurrence of surface
plasmon excitation requires an additional term of the plasmon frequency, which
is linear with the momentum transfer $\hbar q$, as follows [147]:

$\hbar^{2}\omega^{2}=\hbar^{2}\omega_{P}^{2}+\hbar^{2}\alpha
q+\beta\frac{\hbar^{2}q^{2}}{m}+\frac{\hbar^{4}q^{4}}{4m^{2}}\,,$ (140)

where $\alpha=\sqrt{3E_{\mathrm{F}}/(5m)}\omega_{P}$. In this respect, the
equal inclusion of bulk and surface plasmons is important to obtain a correct
interpretation of the measured REEL spectra, which contain both contributions.
In Fig. 28 in the Results section we show the effects of the different
dispersion models with respect to the calculation of the IMFPs [50, 67, 146,
148]. Of course, the extension to finite momentum transfer using RPA is a
source of error, as it is based on the Fermi electron gas theory, in which it
is assumed that the particles are independent or interact via a mean field. In
particular, the short-range exchange and correlation effects, which are
expected to become relatively important at low energy, i.e. scattering at
large angles, which tends to be characterised by a finite $q$, are neglected
in RPA. Such many-body interactions, related to the Coulomb repulsion between
the electrons in the systems and between the incoming and atomic electrons and
the Pauli exclusion principle, can be easily accounted for in the DL and MELF-
GOS models in a phenomenological way. In particular, the exchange interaction
resulting from the indistinguishability of the incident, target and ejected
electrons is typically treated via the Born–Ochkur approximation [149, 150,
151]. In this model, the ELF is modified by rescaling the direct scattering in
the calculation of the differential inverse inelastic mean free path as
follows

$\frac{d\lambda_{\mathrm{inel}}^{-1}}{dW}\,=\,\frac{1}{\pi
a_{0}E}\,\int_{q_{-}}^{q_{+}}\frac{dq}{q}\,[1\,+\,f_{\mathrm{ex}}(q)]\,\mbox{Im}\left[\frac{1}{\bar{\varepsilon}(q,W)}\right]\,,$
(141)

where the exchange function is

$f_{\mathrm{ex}}(q)\,=\,\left(\frac{\hbar
q}{mv}\right)^{4}\,-\,\left(\frac{\hbar q}{mv}\right)^{2}\,,$ (142)

$v$ is the velocity of the electrons, $\hbar q$ is the momentum transfer and
$a_{0}$ is the Bohr radius. In Eq. (141) the integration limits can be
deducted by the energy-momentum conservation laws as

$\hbar\,q_{\pm}\,=\,\sqrt{2\,m\,E}\,\pm\,\sqrt{2\,m\,(E\,-\,W)}\,,$ (143)

where $E$ is the electron energy, $m$ is the electron mass and $W$ the energy
loss. This approximation may require corrections for incident energies below
200 eV, where exchange effects are important and can change the inelastic
scattering cross-section. This correction can be implemented by using an
interference phase [152] to account for the reduced scattering rates due to
the Pauli exclusion principle. We note that the inclusion of exchange effects
in the ELF is generally particularly important for insulators and
semiconductors beyond the optical limit, while they have less dramatic effects
for metals and semimetals due to more effective screening.

We emphasise that all these methods for calculating the ELF suffer from at
least three common shortcomings, namely: (i) the need for experimental optical
data, which are only available for a few materials; (ii) the implicit use of
Fermi electron-gas theory to treat the electron-electron interaction,
neglecting or crudely approximating the exchange-correlation effects; (iii)
the need for an analytical (see Eqs. (139) and (140)) or a built-in (see Eq.
(127)) relation for modelling the dispersion law at finite momentum.

The ab initio methods offer a practicable solution to these disadvantages.
Although these approaches are computationally more expensive than the DL and
MELF-GOS models, they can in principle provide the most accurate results with
regard to the ELF [22, 49]. These methods are discussed in the next section.

####  5.5.3 First principles simulations of the ELF

For the description of electronic motion, the ab initio methods rely only on
the laws of quantum mechanics, without resorting to experiments. In
particular, the calculation of the ELF can be performed at different levels of
theory, with increasing accuracy from the RPA via the Many-Body Perturbation
Theory (MBPT) [3, 4] and the Time-Dependent Density Functional Theory (TDDFT),
[49, 147, 153] over the entire momentum dispersion with the same computational
effort as the ${\bm{q}}\rightarrow 0$ calculations. These models are based on
a different description of the correlated motion of electrons. However, they
are all based on numerical or analytical models of the Coulomb repulsion
between electrons. In this regard, various exchange-correlation functionals
have been developed over the years to treat the Coulomb interaction using
known constraints in the low and high density regions, leading to different
approximating potentials [154]. While the search for universal and all-purpose
functionals remains a challenge, artificial intelligence-based approaches have
also been successfully applied in this regard recently [155, 156].

In general, the ab initio methods within the dielectric formalism provide a
more rigorous treatment of screening effects in the evaluation of the ELF than
the semi-empirical approaches implemented in DL and MELF-GOS [157, 158]. They
can treat the inclusion of excitonic effects, polarons and electron-phonon
interactions at the same level of theory. Their accuracy is higher than that
of DL and MELF-GOS, both in the long wavelength limit and in finite momentum
transfer. In fact, DL and MELF-GOS start with the fitting of experimental data
and thus inherit the finite resolution of the measurement apparatus, which can
hide many transitions [22]. Most importantly, such experimental results must
be available to elaborate these methods. A discussion of the effect of
different dielectric response models on the calculation of REEL spectra, IMFP
and stopping power of graphite and diamond can be found in Ref. [49]. There it
is shown that the use of ab initio calculated ELF leads to a better agreement
between calculated and experimental data, in the low energy range ($<$ 100 eV)
compared to the widely used DL and MELF-GOS models.

The most important ingredient for accessing the ELF is the macroscopic
dielectric function. The latter can be obtained by calculating the microscopic
polarisation function (or density-density response)
$\chi({\bm{r},\bm{r^{\prime}}},W)$, which is defined as follows:

$\rho^{\mathrm{ind}}({\bm{r}},W)=\int
d{\bm{r^{\prime}}}\chi({\bm{r},\bm{r^{\prime}}},W)V^{\rm{ext}}({\bm{r^{\prime}}},W)\,\mbox{,}$
(144)

where $\rho^{\mathrm{ind}}({\bm{r}},W)$ is the electron density induced by the
external Coulomb potential $V^{\mathrm{ext}}$ (see Eq. (113)) and $W$ is the
electron energy transfer. $\chi({\bm{r},\bm{r^{\prime}}},W)$ is associated
with the electronic band structure of the target material. Assuming periodic
boundary conditions, one can take advantage of the Bloch theorem for pristine
crystals, according to which one-particle states can be conveniently expanded
using a plane wave basis as follows:

$\phi_{{\bm{q}},n}({\bm{r}})=\frac{1}{\sqrt{V}}\sum_{\bm{G}}u_{{\bm{q}},n}({\bm{G}})\exp^{i({\bm{q}}+{\bm{G}})\cdot{\bm{r}}},$
(145)

where $V$ is the volume of the simulation cell, ${\bm{q}}$ is the electron
momentum transfer located in the 1BZ, ${\bm{G}}$ is a reciprocal lattice
vector, the quantum number $n$ characterises the band and the function
$u_{{\bm{q}},n}$ has the same periodicity as the crystal lattice. The
polarisability in Eq. (144) can also be conveniently expressed for periodic
crystals by using the Fourier transform from real to reciprocal space as a
matrix $\chi_{\bm{G},\bm{G^{\prime}}}({\bm{q}},W)$ in which the reciprocal
lattice points $\bm{G},\bm{G^{\prime}}$ identify the rows and columns for each
momentum and energy transfer, respectively. The full polarisability
$\chi_{\bm{G},\bm{G^{\prime}}}({\bm{q}},W)$ of a periodic solid in reciprocal
space can then be obtained by solving the following Dyson-like screening
equation [56, 159]:

$\displaystyle\chi_{\bm{G},\bm{G^{\prime}}}({\bm{q}},W)=\chi^{0}_{\bm{G},\bm{G^{\prime}}}({\bm{q}},W)+\sum_{\bm{G^{\prime\prime}},\bm{G^{\prime\prime\prime}}}\chi^{0}_{\bm{G},\bm{G^{\prime\prime}}}({\bm{q}},W)$
$\displaystyle\times\Big{[}v_{\bm{G^{\prime\prime}}}({\bm{q}})\delta_{\bm{G^{\prime\prime}},\bm{G^{\prime\prime\prime}}}+K^{xc}_{\bm{G^{\prime\prime}},\bm{G^{\prime\prime\prime}}}({\bm{q}},W)\Big{]}\chi_{\bm{G^{\prime\prime\prime}},\bm{G^{\prime}}}({\bm{q}},W),$
(146)

which relates the non-interacting density-density response function
$\chi^{0}_{\bm{G},\bm{G^{\prime}}}({\bm{q}},W)$ for independent particles to
the full response. The single-particle states appearing in
$\chi^{0}_{\bm{G},\bm{G^{\prime}}}({\bm{q}},W)$ are usually derived by solving
the mean-field Hartree-Fock equations or using the Kohn-Sham DFT (KS-DFT)
scheme. In the KS-DFT scheme, the polarisability of the independent particles
is, for example:

$\chi^{0}_{\bm{G},\bm{G^{\prime}}}({\bm{q}},W)=2\sum_{vc}\frac{\bra{v}e^{-i({\bm{q}+\bm{G}})\cdot\bm{r}}\ket{c}\bra{c}e^{i({\bm{q}+\bm{G^{\prime}}})\cdot\bm{r^{\prime}}}\ket{v}}{W-(\varepsilon_{c}-\varepsilon_{v})+i\eta},$
(147)

where $v,c$ and $\varepsilon_{v},\varepsilon_{c}$ are the KS eigenfunctions
and eigenvalues corresponding to the valence and conduction bands of the
crystal. Using the expression (147) instead of the Dyson Eq. (146), also known
as RPA, can lead to results that agree well with the experimental data (see
section 6.2), even though the single-particle KS polarisability is generally
not satisfactory for the description of electronic spectra.

In Eq. (146) the exchange-correlation kernel
$K^{xc}_{\bm{G^{\prime\prime}},\bm{G^{\prime\prime\prime}}}({\bm{q}},W)=\left\\{\frac{d}{d\rho}v_{\rm
xc}[\rho]\right\\}_{\rho=\rho(\bm{r},t)}$ encodes all non-trivial many-body
effects, while
$v_{\bm{G^{\prime}}}({\bm{q}})=4\pi/|{\bm{q}}+{\bm{G}^{\prime}}|^{2}$ is the
Fourier-transformed bare Coulomb interaction. The kernel
$K^{xc}_{\bm{G^{\prime\prime}},\bm{G^{\prime\prime\prime}}}({\bm{q}},W)$ can
be modelled using the adiabatic local density approximation (ALDA) [160], by
which the time dependence is local over the instantaneous DFT ground state
density

$v_{\rm xc}[\rho]({\bm{r}},t)=v_{\rm xc}^{\rm adia}[\rho]({\bm{r}},t)=v_{\rm
xc}^{\rm
GS}[\rho]({\bm{r}})|_{\rho({\bm{r}^{\prime}},t)=\rho(\bm{r}^{\prime})_{GS}},$
(148)

where $v_{\rm xc}^{\rm
GS}[\rho]({\bm{r}})|_{\rho({\bm{r}^{\prime}},t)=\rho(\bm{r}^{\prime})_{GS}}$
is the exchange-correlation potential for the ground-state density
$\rho(\bm{r}^{\prime})_{GS}$. $v_{\rm xc}^{\rm GS}[\rho]({\bm{r}})$ can be
calculated with any variant of the standard exchange correlation functionals,
such as LDA, GGA or hybrid functionals to account for exact exchange [160].
Recently, other functionals with memory (non-local in time) have also been
developed [161].

The full polarisability ${\chi_{\bm{G},\bm{G^{\prime}}}}({\bm{q}},W)$ is
related to the microscopic dielectric function as follows:

$\varepsilon_{\bm{G},\bm{G^{\prime}}}({\bm{q}},W)=1-v_{\bm{G}^{\prime}}({\bm{q}})\chi_{\bm{G},\bm{G^{\prime}}}({\bm{q}},W)\,\mbox{,}$
(149)

where $1=\delta_{\bm{G},\bm{G^{\prime}}}$ is the identity operator. Eq. (149)
shows that the microscopic $\varepsilon_{\bm{G},\bm{G^{\prime}}}$ in the
entire Bethe surface can ultimately be obtained from electronic structure
calculations, i.e. from the knowledge of the electronic wave functions (see
Eq. (100)). Finally, we recall that the absorption coefficient according to
Eq. (108) is proportional to ${\rm Im}(\bar{\varepsilon})$:

$\alpha\propto{\rm Im}(\bar{\varepsilon})=-v\times{\rm Im}(\bar{\chi})\,,$
(150)

while the electron energy loss from Eq. (5.3.2) is proportional to the
imaginary part of the reciprocal of the macroscopic dielectric function:

${\mathrm{ELF}}\propto{\rm
Im}\left(\frac{1}{\bar{\varepsilon}}\right)=-v\times{\rm
Im}\left(\frac{\bar{\chi}}{1-v\bar{\chi}}\right).$ (151)

where the macroscopic polarizability $\bar{\chi}$ is the function for which
$\bar{\varepsilon}=1-v\bar{\chi}_{\bm{G}=\bm{G^{\prime}}=0}({\bm{q}},W)$. We
note that macroscopic observables are typically obtained from experiments,
such as the absorption coefficient in Eq. (108) or the EEL spectra. In order
to be able to directly compare computer simulations and experiments, it is
therefore necessary to carry out a spatial average of the variables in Eq.
(149) (see the discussion on the average values in section 5.4).

###  5.6 Inelastic scattering: the Bethe-Bloch formula

The dielectric theory describes energy loss as a polarisation phenomenon that
can be modelled macroscopically by Maxwell’s electromagnetic field equations
and microscopically by the quantum response of a target to an external
perturbation, such as the collision of an electron beam. Previously, Bethe
[110, 162, 163] had proposed a different approach to calculating stopping
power. This approach was based on the idea that electromagnetic interactions
triggered by the impact of a charged projectile on a target system can be
described by the classical response of a collection of independent harmonic
oscillators with given energy spacing and dipole strength. In particular, the
excitation cross-section for such an electron-atom collision is [104]:

$\sigma_{j}=\frac{2\pi
e^{4}}{mv^{2}}\int\frac{dW}{W^{2}}|F_{j0}(\bm{q})|^{2},$ (152)

where $\bm{q}$ and $W=\hbar^{2}q^{2}/2m$ are the momentum and energy transfer
and

$F_{n0}(\bm{q})=\bra{n}\sum_{\nu=1}^{Z}e^{i\bm{q}\cdot{\bm{r}_{\nu}}}\ket{0}$
(153)

is the form factor, the square of which is the generalised atomic oscillator
strength, for a transition from the ground state $0$ to the excited state $n$
of the system of $Z$ charges in the positions $\bm{r}_{\nu}$. The presence of
$F_{n0}(\bm{q})$ in the following equations suggests the interpretation of the
variable ${\bm{q}}$ as momentum transfer. This cross-section can be derived
using semiclassical arguments, where bound electrons can be considered as
quantised harmonic oscillators [44]; however, we will discuss the rigorous
full quantum mechanical derivation, which can also be more easily extended to
relativistic regimes [104].

In this context, we assume that a homogeneous, isotropic and infinite medium
is perturbed by an external point charge, as in Eq. (109). From Eq. (119) the
stopping power results as follows:

$-\frac{dW}{dx}=-\frac{e}{v}{\bm{v}}\cdot{\bm{E}_{\mathrm{tot}}}({\bm{v}},t),$
(154)

as the stopping force is opposite to the velocity of the particles. By using
the scalar and vector potentials ($\Phi({\bm{r}},t),{\bm{A}}({\bm{r}},t))$ to
represent the electromagnetic fields and assuming Coulomb gauge
(${\bm{\nabla}}\cdot{\bm{A}}=0$), the field equations can be written in
Fourier space (see Eq. (114)):

$\displaystyle q^{2}\Phi({\bm{q}},\omega)=4\pi\rho_{\rm
tot}({\bm{q}},\omega),$ (155)
$\displaystyle\left(q^{2}-\frac{\omega^{2}}{c^{2}}\right){\bm{A}}({\bm{q}},\omega)-\frac{\omega}{c}{\bm{q}}\Phi({\bm{q}},\omega)=\frac{4\pi}{c}{\bm{j}}_{\rm
tot}({\bm{q}},\omega).$ (156)

In a polarisable medium, the total charge and current densities can be divided
into induced and external contributions as follows (see Eqs. (112), (113)):

$\displaystyle\rho_{\rm tot}=\rho_{\rm ind}+\rho_{\rm ext},$ (157)
$\displaystyle{\bm{j}}_{\rm tot}={\bm{j}}_{\rm ind}+{\bm{j}}_{\rm
ext}=\sigma{\bm{E}_{\rm tot}},$ (158)

where the induced components represent the (linear) response of the system to
the external perturbation of Eq. (109). Eq. (158) is the Ohm’s law (see also
Eq. (98)). Furthermore, the vector potential and the wave vector are
orthogonal in the Coulomb gauge, so that by projecting Eq. (156) into
transverse (orthogonal to ${\bm{q}}$) and longitudinal (parallel to
${\bm{q}}$) components we obtain

$\displaystyle\frac{\omega}{c}{\bm{q}}\left(1+\frac{4\pi i\sigma_{\rm
long}({q},\omega)}{\omega}\right)\Phi({\bm{q}},\omega)=\frac{4\pi}{c}{\bm{j}}_{\rm
ext,long}({\bm{q}},\omega),$ (159)
$\displaystyle\left(q^{2}-\frac{\omega^{2}}{c^{2}}\right)\left(1+\frac{4\pi
i\sigma_{\rm
trans}({q},\omega)}{\omega}\right){\bm{A}}({\bm{q}},\omega)=\frac{4\pi}{c}{\bm{j}}_{\rm
ext,trans}({\bm{q}},\omega),$ (160)

where $\sigma_{\rm trans},\sigma_{\rm long}$ and ${\bm{j}}_{\rm
ext,trans},{\bm{j}}_{\rm ext,long}$ are the conductivity and external current
density along the transverse and longitudinal directions, respectively.

Using the continuity equation:

${\bm{q}}\cdot{\bm{j}}_{\rm ext}({\bm{q}},\omega)=\omega\rho_{\rm
ext}({\bm{q}},\omega),$ (161)

one can rewrite Eq. (159) as:

$q^{2}\left(1+\frac{4\pi i\sigma_{\rm
long}({q},\omega)}{\omega}\right)\Phi({\bm{q}},\omega)=4\pi\rho_{\rm
ext}({\bm{q}},\omega)=k^{2}\varepsilon_{\rm long}\Phi({\bm{q}},\omega),$ (162)

and, analogously, Eq. (160) as:

$\left(q^{2}-\frac{\omega^{2}}{c^{2}}\varepsilon_{\rm
trans}\right){\bm{A}}({\bm{q}},\omega)=\frac{4\pi}{c}{\bm{j}}_{\rm
ext,trans}({\bm{q}},\omega),$ (163)

where we have defined the dielectric functions in the longitudinal and
transverse directions $\varepsilon_{\rm long}=\left(1+\frac{4\pi i\sigma_{\rm
long}({q},\omega)}{\omega}\right)$ and $\varepsilon_{\rm
trans}=\left(1+\frac{4\pi i\sigma_{\rm trans}({q},\omega)}{\omega}\right)$,
respectively. By inserting Eqs. (162), (163) and (5.3.2) in Eq. (154) one
obtains [104]:

$\displaystyle-\frac{dW}{dx}=\frac{ie^{2}}{\pi
v^{2}}\int_{0}^{\infty}\frac{dq}{q}\int_{-qv}^{qv}d\omega\omega\left(\frac{1}{\varepsilon_{\rm
long}({k},\omega)}-\frac{v^{2}}{c^{2}}\frac{q^{2}-\omega^{2}/v^{2}}{q^{2}-\varepsilon_{\rm
trans}({q},\omega)\omega^{2}/c^{2}}\right).$

In general, the second term in Eq. (5.6), which represents the contribution of
the transverse excitation to the stopping force, can be neglected in the non-
relativistic (up to the quasi-relativistic, typically for $\beta=v/c<0.5$)
regime, i.e. $dW/dx\approx dW/dx|_{\rm long}$.

To simplify the derivation of the Bethe-Bloch stopping formula, Eq. (152), we
assume that the target system is represented by a one-electron atom at
position ${\bm{R}}$ and the impinging projectile is a non-relativistic
electronic charge unit, so that only the first term in Eq. (5.6) can be
retained. For materials, this hypothesis naturally assumes that the electronic
clouds of the atomic components do not overlap significantly. The generalised
oscillator strength $|F_{j0}(\bm{q})|^{2}$ of Eq. (153) for this problem can
be determined using first-order perturbation theory. The interaction between
the perturbing electromagnetic field generated by the incoming charged
particle and the one-electron target atom is described by the potential energy

$V({\bm{r}},t)=-e\Phi({\bm{r}}-{\bm{R}},t)=\frac{-e^{2}}{|{\bm{r}}-{\bm{p}}(t)|},$
(165)

where $\Phi({\bm{r}}-{\bm{R}},t)$ is the potential generated by the atomic
electron at the position $({\bm{r}}-{\bm{R}})$ and
${\bm{c}}(t)={\bm{r}}_{0}+{\bm{v}}t-{\bm{R}}={\bm{p}}(t)-{\bm{R}}$ is the
trajectory of the impacting projectile characterised by the impact parameter
${\bm{r}}_{0}$. Within the framework of the first-order perturbative approach,
the wave function of the bound electron can be expanded in the eigenstates of
the unperturbed Hamiltonian $H_{0}\ket{n}=\varepsilon_{n}\ket{n}$ of the
target as follows:

$\Psi({\bm{r}},t)=\sum_{n}c_{n}(t)e^{-i\varepsilon_{n}t/\hbar|}\ket{n},$ (166)

where

$c_{n}(t)\simeq\delta_{n0}+c_{n}^{(1)}=\delta_{j0}+\frac{1}{i\hbar}\int_{-\infty}^{t}dt^{\prime}e^{i\omega_{n0}t^{\prime}}\bra{n}V({\bm{r}},t^{\prime})\ket{0}.$
(167)

In Eq. (167) $\delta_{j0}$ is the Kronecker symbol and
$\hbar\omega_{n0}=\varepsilon_{n}-\varepsilon_{0}$ is the energy jump between
the ground state and the $n$th excited state of the system. In Fourier space,
the matrix element in Eq. (167) is given using Eq. (165) by

$\displaystyle\bra{n}V({\bm{r}},t)\ket{0}$ $\displaystyle=$
$\displaystyle-e\int d{\bm{q}}\int d\omega\int
d{\bm{r}}\Phi({\bm{q}},\omega)\psi^{*}_{n}({\bm{r}}-{\bm{R}})\psi_{0}({\bm{r}}-{\bm{R}})\rangle
e^{-i\omega t}e^{i{\bm{q}}\cdot{\bm{r}}}$ (168) $\displaystyle=$
$\displaystyle-e\int d{\bm{q}}\int
d\omega\Phi({\bm{q}},\omega)e^{i{\bm{q}}\cdot{\bm{R}}}F_{n0}({\bm{q}})e^{-i\omega
t},$

where $F_{n0}({\bm{q}})$ is given in Eq. (153) and
$\psi_{n,0}({\bm{r}}-{\bm{R}})$ are the electron wavefunctions corresponding
to the $n,0$ states. We emphasise that the presence of $F_{n0}(\bm{q})$
implies the interpretation of $\hbar\bm{q}$ as the momentum transfer
associated with the electronic excitation.

At time $t=\infty$, which is asymptotically far away from the scattering
centre, where the potential has a negligible effect on the target, we thus
obtain from Eq. (167):

$c_{n}^{(1)}(\infty)=-\frac{e}{\hbar}\int d{\bm{q}}\int
d\omega\Phi({\bm{q}},\omega)e^{i{\bm{q}}\cdot{\bm{R}}}F_{n0}({\bm{q}})\frac{e^{i(\omega_{n0}-\omega)t}}{\omega_{n0}-\omega-i\Gamma},$
(169)

where the damping constant $\Gamma$ was added to force the convergence of the
integral in the complex plane near the singular points.

From Eqs. (155) and (162), we can also write:

$\rho_{\rm ind}({\bm{q}},\omega)=-\frac{q^{2}}{4\pi}[\varepsilon_{\rm
long}-1]\Phi({\bm{q}},\omega).$ (170)

For a medium with a number density of one-electron atoms $N$, the total
induced charge density can be assessed to first order in $\Phi$ by
perturbation theory (Eq. (169)) as:

$\displaystyle N\int d{\bm{R}}~{}\rho_{\rm ind}({\bm{r}},t)$ $\displaystyle=$
$\displaystyle-Ne\int
d{\bm{R}}~{}(\psi^{(0)*}({\bm{r}},t)\psi^{(1)}({\bm{r}},t)+\psi^{(0)}({\bm{r}},t)\psi^{(1)*}({\bm{r}},t))$
(171) $\displaystyle=$ $\displaystyle-\frac{Ne^{2}}{\hbar}\int
d{\bm{R}}~{}\Big{[}\sum_{n}\Big{(}\psi_{0}^{*}({\bm{r-R}},t)\psi_{n}({\bm{r-R}},t)\times$
$\displaystyle\times$ $\displaystyle\int d{\bm{q}}\int
d\omega\Phi({\bm{q}},\omega)e^{i{\bm{q}}\cdot{\bm{R}}}F_{n0}({\bm{q}})\frac{e^{-i\omega
t}}{\omega_{n0}-\omega-i\Gamma}\Big{)}{\rm+c.c.}\Big{]},$ $\displaystyle=$
$\displaystyle-\frac{Ne^{2}}{\hbar}\Phi({\bm{q}},\omega)\Bigg{(}\sum_{n}\frac{F_{0n}({-\bm{q}})F_{n0}({\bm{q}})}{\omega_{n0}-\omega-i\Gamma}+$
$\displaystyle+$
$\displaystyle\sum_{n}\frac{F_{0n}({\bm{q}})F_{n0}({-\bm{q}})}{\omega_{n0}+\omega+i\Gamma}\Bigg{)}.$

Finally, by inserting Eq. (170) into the left side of Eq. (171), we obtain the
longitudinal dielectric function for an isotropic medium as [104]:

$\displaystyle\varepsilon_{\rm long}(q,\omega)$ $\displaystyle=$
$\displaystyle 1+\frac{m\omega^{2}_{P}}{\hbar
q^{2}}\sum_{n}F_{n0}^{2}({\bm{q}})$ (172) $\displaystyle\times$
$\displaystyle\left(\frac{1}{\omega_{n0}-\omega-i\Gamma}+\frac{1}{\omega_{n0}+\omega+i\Gamma}\right),$

where $\omega^{2}_{P}=4\pi Ne^{2}/m$.

We note that for a non-relativistic Fermi gas of free electrons, the wave
function is represented by a plane wave, i.e:

$\ket{n}=\frac{1}{\sqrt{V}}e^{i{\bm{k}}_{n}\cdot{\bm{r}}},$ (173)

where $V$ is the volume occupied by the gas. Eq. (153) then results in
$F_{n0}({\bm{q}})=\delta_{{\bm{k}}_{n},{\bm{q}}}$, which used in Eq. (172)
yields:

$\varepsilon_{\rm
long}(q,\omega)=1+\frac{\omega^{2}_{p}}{\omega_{q}^{2}-(\omega+i\Gamma)^{2}},$
(174)

which is the DL function given in Eq. (133) to fit the experimental optical
ELF with $\omega_{q}=\hbar q^{2}/2m$. From this expression one can derive

$\frac{i}{\varepsilon_{\rm
long}(q,\omega)}=-\frac{i\omega^{2}_{p}}{\omega_{q}^{2}+\omega_{P}^{2}-(\omega+i\Gamma)^{2}},$
(175)

which can be inserted into the non-relativistic limit of Eq. (5.6) by
neglecting the second term to obtain :

$-\frac{dW}{dx}\approx-\frac{dW}{dx}{\Bigg{|}}_{\mathrm{long}}=\frac{e^{2}\omega_{P}^{2}}{v^{2}}\int\frac{dq}{q}=\frac{e^{2}\omega_{P}^{2}}{2v^{2}}\ln\left(\frac{\zeta+\sqrt{\zeta^{2}-1}}{\zeta-\sqrt{\zeta^{2}-1}}\right),$
(176)

where $\zeta=mv^{2}/\hbar\omega_{P}$. For $\zeta>>1$, i.e. for projectiles
with kinetic energies much greater than the mean excitation energy, the
asymptotic expansion of Eq. (176) can be written as follows:

$-\frac{dW}{dx}=\frac{4\pi
Ne^{4}}{mv^{2}}\Bigg{(}\ln\left(\frac{2mv^{2}}{\hbar\omega_{P}}\right)-\left(\frac{\hbar\omega_{P}}{2mv^{2}}\right)^{2}-\frac{3}{2}\left(\frac{\hbar\omega_{P}}{2mv^{2}}\right)^{4}+...\Bigg{)},$
(177)

which is Bethe’s asymptotic formula for a non-relativistic charge that hits a
target at high speed ($mv^{2}>>\hbar\omega_{P}$), which can be approximated as
a Fermi gas with free electrons. As an aside, it should be noted that if
Cauchy’s main partial equation

$\frac{1}{\omega-\omega_{n0}+i\Gamma}={\cal
P}\frac{1}{\omega-\omega_{n0}}-i\pi\delta({\omega-\omega_{n0}})$ (178)

is used in Eq. (175), we obtain (the stopping power is a real quantity):

${\rm Re}\Bigg{(}\frac{i}{\varepsilon_{\rm
long}(q,\omega)}\Bigg{)}=-\frac{\pi\omega^{2}_{p}}{2\alpha_{q}}\Big{(}\delta(\omega-\alpha_{q})-\delta(\omega+\alpha_{q})\Big{)},$
(179)

where $\alpha_{q}^{2}=\omega^{2}=\omega_{q}^{2}+\omega_{P}^{2}=(\hbar
q^{2}/2m)^{2}+\omega_{P}^{2}$, which reveals the quadratic dependence of the
energy-momentum transfer relation (see Eq. (139)). Due to the presence of the
Dirac functions, this expression of the inverse dielectric function, which is
proportional to the ELF, shows well-resolved resonances or plasmon poles and
is therefore referred to as the plasmon pole approximation.

The Bethe stopping formula can also be obtained in the limiting case of low
target density by using the non-relativistic limit of Eq. (5.6) and inserting
the expansion of Eq. (172) up to first order in $N$:

$\frac{1}{\varepsilon_{\rm long}}\approx 1-\frac{m\omega^{2}_{P}}{\hbar
q^{2}}\sum_{j}|F_{n0}({\bm{q}})|^{2}\times\left(\frac{1}{\omega_{n0}-\omega-i\Gamma}+\frac{1}{\omega_{n0}+\omega+i\Gamma}\right).$
(180)

Eq. (5.6) can be integrated over $\omega$ by using Eq. (178) to obtain the
Bethe-Bloch formula for the stopping power of a non-relativistic target with
low density:

$-\frac{dW}{dx}=\frac{2\pi
e^{4}N}{mv^{2}}\sum_{n}\hbar\omega_{n0}\int\frac{dW}{W^{2}}|F_{n0}({\bm{q}})|^{2},$
(181)

which is equivalent to Eq. (152) by interpreting the integrand as a
differential cross-section.

In the relativistic regime (typically for $\beta=v/c>0.9$), we are reminded
that (i) the transverse contribution to the stopping power (second term in Eq.
(5.6)) cannot be neglected; (ii) the wave functions are characterised by four-
component spinors as in Eq. (46). Assuming that we are dealing with a gas of
free relativistic electrons, the derivation of the longitudinal dielectric
function for relativistic particles follows the same reasoning that we have
presented to obtain Eq. (172), which in this case leads to the following
expression [104]:

$\displaystyle\varepsilon_{\rm long}(q,\omega)$ $\displaystyle=$
$\displaystyle 1+\frac{m\omega^{2}_{P}}{\hbar
q^{2}}B_{q}^{2}\Bigg{(}\frac{1}{\omega_{q}^{+}-\omega-i\Gamma}+\frac{1}{\omega_{q}^{+}+\omega+i\Gamma}+$
$\displaystyle+$
$\displaystyle\frac{(b_{q}q)^{2}}{\omega_{q}^{-}-\omega-i\Gamma}+\frac{(b_{q}q)^{2}}{\omega_{q}^{-}+\omega+i\Gamma}\Bigg{)}=$
$\displaystyle=$ $\displaystyle
1+\frac{\omega^{2}_{P}}{\omega_{q}^{2}-(1+\hbar\omega_{q}/mc^{2})(\omega+i\Gamma)^{2}+\hbar^{2}(\omega+i\Gamma)^{4}/(2mc^{2})^{2}},$

where $\hbar\omega^{\pm}_{q}=\pm E_{q}-mc^{2}$ represents the excitation
energy for states with positive (+) or negative (-) energy,
$B_{q}^{2}=(E_{q}+mc^{2})/(2E_{q})$ and $b_{q}=\hbar c/(E_{q}+mc^{2})$. We
note that Eq. (5.6) is far more complex than Eq. (172) valid for the non-
relativistic case; however, it is reduced to the first term of the expansion
in Eq. (177) under the assumption of $\hbar\omega_{P}/(2mc^{2})<<1$.

The contribution of the transverse excitation to the stopping power can be
rewritten as follows

$\displaystyle-\frac{dW}{dx}{\Bigg{|}}_{{\rm trans}}$ $\displaystyle=$
$\displaystyle\frac{2\pi
e^{4}N}{mv^{2}}\sum_{n}\hbar\omega_{n0}\int\frac{dW}{W^{2}}\times|F_{n0}({\bm{q}})|^{2}\frac{2mv^{2}W/(\hbar\omega_{n0})^{2}-1}{[2mc^{2}W/(\hbar\omega_{n0})^{2}-1]^{2}}$
$\displaystyle=$ $\displaystyle\frac{2\pi
e^{4}NZ}{mv^{2}}\sum_{n}\int\frac{dW}{W}\times|f_{n0}({W})|^{2}\frac{2mv^{2}W/(\hbar\omega_{n0})^{2}-1}{[2mc^{2}W/(\hbar\omega_{n0})^{2}-1]^{2}},$

for $W>\hbar\omega_{n0}$, where

$f_{n0}(W)=\frac{1}{Z}\frac{\hbar\omega_{n0}}{W}|F_{n0}({\bm{q}})|^{2}$ (184)

are the generalised oscillator strengths that can be approximated by a Poisson
distribution
$f_{n0}(W)\propto(W/\hbar\omega_{0})^{(n-1)}e^{-W/(\hbar\omega_{0})}$ in a
semi-classical theory in which the bound electron is modelled as a harmonic
oscillator.

Given the velocity-dependent weight factor in the integrand, one can
conveniently express Eq. (LABEL:rel_5) in the limits of low (essentially
within the validity of the dipole approximation, a regime typical of weak or
distant interactions in real space), medium (where the non-relativistic
description of particle dynamics is still valid) and high (relativistic,
typical of strong or close collisions in real space) momentum transfers [104].
Summing up all these contributions gives the relativistic Bethe formula for a
charged projectile:

$-\frac{dW}{dx}=\frac{4\pi
e^{4}NZ}{mv^{2}}\sum_{n}f_{n0}\Bigg{[}\ln\left(\frac{2mv^{2}}{\hbar\omega_{n0}}\right)-\ln\left(1-\frac{v^{2}}{c^{2}}\right)-\frac{v^{2}}{c^{2}}\Bigg{]}$
(185)

A comparison with Eq. (177) shows that the second and third terms in Eq. (185)
come from relativistic corrections. We also note that for $v<<c$, i.e. far
away from the relativistic limit, the expansion in Eq. (185) consists only of
the first term, which can be rewritten as

$-\frac{dW}{dx}=\frac{4\pi e^{4}NZ}{mv^{2}}\ln\frac{2mv^{2}}{I},$ (186)

where $\ln I=\sum_{n}f_{n0}\ln(\varepsilon_{n}-\varepsilon_{0})$ is commonly
referred to as the mean logarithmic excitation energy. 333Please note that if
we take into account the indistinguishability of the scattered electrons from
the ejected ones, Eq. (186) becomes  $-\frac{dW}{dx}=\frac{4\pi
e^{4}NZ}{mv^{2}}\ln\sqrt{\frac{\rm{e}}{2}}\,\frac{mv^{2}}{2I}=\frac{2\pi
e^{4}NZ}{E}\ln\frac{1.166E}{I}$ where $E=mv^{2}/2$ is the electron kinetic
energy and $\rm{e}$ is Euler’s number [44, 57, 58, 104].

Further details on the derivation of the relativistic Bethe relation with an
in-depth analysis of the different momentum transfer regimes can be found in
Ref. [104]. The Bethe-Bloch stopping power formula can be used in conjunction
with the continuous-slowing-down approximation to model charge transport in
solids as part of an MC approach.

###  5.7 Inelastic scattering: the Auger decay

Depending on the kinetic energy of the impinging projectile, inelastic
electron scattering can lead to excitation and ionisation of the atomic
constituents of a solid. Especially if the excitation or ionisation affects
inner shells, Auger states can be generated. These are quasi-bound states that
can decay through the emission of a secondary electron, the Auger electron, or
through radiation. The decaying system is in a singly or doubly ionised state,
depending on whether the initial excitation is neutral (autoionising state) or
ionising (Auger state), i.e. whether the projectile has enough energy to
promote an electron into an empty bound orbital or into the continuum. Often
the primary excitation process can lead to simultaneous ionisation (shake-off)
or excitation (shake-up) of a valence electron, resulting in satellite lines
in the spectrum (typically these processes account for 10% of the total
spectral intensity). In this section, dedicated to Auger decay, we use atomic
units to simplify the equations.

In general, we consider ionisation and Auger emissions as resonance processes
in which the excited system decays to the higher charge state. Auger decay in
particular is typically treated as a two-step process, where the first
ionisation is followed by the non-radiative decay. In order to determine the
contribution of these processes to the line shape (see Fig. 2), the
corresponding autoionisation or Auger cross-sections must be evaluated. This
task can be fulfilled within the framework of Fano’s theory of decay processes
[164] in conjunction with the projected potential approach developed in
section 5.2 for elastic scattering. In the framework of Fano, the scattering
wave function, which is an eigenstate of the full Hamiltonian of the system,
can be represented by the following linear combination of configurations

$|\Psi^{-}_{\alpha,{\bm{k}}}\rangle=\sum_{j}^{m_{d}}a_{j\alpha}({\bm{k}})|\Phi_{j}\rangle+\sum_{\beta=1}^{m_{c}}\int\limits_{0}^{\infty}\frac{d{\bm{p}}}{(2\pi)^{3}}\ket{\chi^{-}_{\beta,{\bm{p}}}}b_{\beta,\alpha}({\bm{p}},{\bm{k}}),$
(187)

where $|\Phi_{j}\rangle$ is a discrete bound resonance state generated by the
primary collision and degenerate with several continua
$\\{|\chi_{\beta,{\bm{p}}}^{-}\rangle\\}$ that asymptotically represent the
decay channels of the ionised target [165] identified by the quantum number
$\beta$ (e.g. the angular momentum and its projection along an axis in
spherically symmetric systems) and the momentum $\bm{p}$ of the emitted
electron. $|\Psi^{-}_{\alpha{\bm{k}}}\rangle$ is a stationary scattering state
at energy $E=E_{\alpha}+{k^{2}\over 2}$, which is asymptotically defined by
the momentum ${\bm{k}}$ of the emitted electron and the state $|\alpha\rangle$
of the singly ionised target (for an autoionisation process) or doubly ionised
target (for an Auger process). Eq. (187) is of course written under the
assumption that the Born–Oppenheimer approximation is valid, i.e. that the
typical time scale of electronic excitation/de-excitation is much faster than
that of the nuclei. Therefore, the discrete intermediate states can represent,
for example, electronic bound states with the corresponding roto-vibrational
levels of the target molecule. The continuum wave functions
$\\{|\chi_{\alpha,{\bm{k}}}^{-}\rangle\\}$ are characterised by the incoming
wave boundary conditions to represent a free electron asymptotically, which is
ejected with the kinetic energy ${\epsilon=k^{2}/2}$ in the direction
${\bm{k}}$ and a system in the ionised state $\alpha$ with the energy
$E_{\alpha}$, so that $E=E_{\alpha}+\epsilon$. The coefficients of Eq. (187)
are obtained by solving the following equations

$\langle\Phi_{l}|H-E|\Psi^{-}_{\alpha,{\bm{k}}}\rangle=\langle\chi^{-}_{\gamma,{\bm{\tau}}}|H-E|\Psi^{-}_{\alpha,{\bm{k}}}\rangle=0~{}~{}~{}~{}\forall~{}l,\gamma,{\bm{\tau}},$
(188)

which explicitly read:

$(E_{l}-E)a_{l\alpha}({\bm{k}})+\sum_{\beta=1}^{m_{c}}\int\limits_{0}^{\infty}\frac{d{\bm{p}}}{(2\pi)^{3}}M_{l\beta}(E,{\bm{p}})b_{\beta\alpha}({\bm{p}},{\bm{k}})=0\,,$
(189) $\sum_{j}^{m_{d}}a_{j\alpha}({\bm{k}})M_{\gamma
j}^{-}(\tau,E)^{*}+(\epsilon+E_{\gamma}-E)b_{\beta\alpha}({\bm{\tau}},{\bm{k}})=0\,,$
(190)

where the elements of the Hamiltonian matrix constructed using discrete and
continuum functions are the following

$\displaystyle\langle\Phi_{l}|H-E|\chi_{\beta,{\bm{p}}}^{-}\rangle=M_{l\beta}(E,{\bm{p}}),$
(191) $\displaystyle\bra{\Phi_{l}}H-E|\Phi_{j}\rangle=(E_{l}-E)\delta_{lj},$
(192)
$\displaystyle\bra{\chi^{-}_{\gamma,{\bm{\tau}}}}H-E\ket{\chi^{-}_{\beta,{\bm{p}}}}=(E_{\gamma}+\frac{\tau^{2}}{2}-E)\delta_{\gamma\beta}(2\pi)^{3}\delta({\bm{\tau}}-{\bm{p}}).$
(193)

The matrix element $M_{l\beta}$ in Eq. (191) couples the discrete electronic
state with the continuum state and from Eq. (192) ${E_{l}}$ is the energy of
the discrete level corresponding to the quasi-bound states
${|\Phi_{l}\rangle}$. The states $\ket{\chi^{-}_{\beta,\tau}}$ in Eq. (187)
are continuum states obtained by diagonalising the full Hamiltonian of the
system

$\displaystyle\langle\chi_{\beta,\epsilon}|H-E|\chi_{\gamma,\epsilon^{\prime}}\rangle=(E_{\beta}+\epsilon-E)\delta(E_{\beta}+\epsilon-
E_{\gamma}-\epsilon^{\prime})+V_{\beta\gamma}(\epsilon,\epsilon^{\prime};E),$

where $\\{|\chi_{\beta,\epsilon}\rangle\\}$ are interacting continuum states.

To avoid singularities in the coefficients $\\{a_{\alpha},b_{\beta\alpha}\\}$,
Eqs. (189) and (190) are shifted into the complex plane, where we obtain the
following expression for the scattering wave function

$\displaystyle|\Psi^{-}_{\alpha,{\bm{k}}}\rangle=|\chi^{-}_{\alpha,{\bm{k}}}\rangle+\sum_{jl}{M_{l\alpha}(E,{\bm{k}})\over(E-E_{l})\delta_{lj}-\Delta_{lj}-i{\Gamma_{lj}\over
2}}\times$ $\displaystyle\left[|\Phi_{j}\rangle+\lim_{\nu\rightarrow
0}\sum_{\beta}^{m_{c}}\int\limits_{0}^{\infty}\frac{d\bm{p}}{2\pi^{3}}{|\chi^{-}_{\beta,{\bm{p}}}\rangle
M^{*}_{\beta j}({\bm{p}},E)\over E-E_{\beta}-\frac{p^{2}}{2}-i\nu}\right],$
(195)

where

$\Gamma_{lj}=2\pi\sum_{\beta}^{m_{c}}\int\limits_{0}^{\infty}\frac{d\bm{p}}{2\pi^{3}}M_{l\beta}(E,{\bm{p}})M^{*}_{j\beta}(E,{\bm{p}})\delta(E-E_{\beta}-\frac{p^{2}}{2}),$
(196) $\Delta_{lj}=\sum_{\beta}^{m_{c}}{\cal
P}\int\limits_{0}^{\infty}\frac{d\bm{p}}{2\pi^{3}}\frac{M_{l\beta}(E,{\bm{p}})M^{*}_{j\beta}(E,{\bm{p}})}{E-E_{\beta}-\frac{p^{2}}{2}}.$
(197)

represent the width and energy shift of the transition.

Finally, if we know for each quantum number $\alpha$ and each momentum
${\bm{k}}$, which asymptotically define the remaining ionised system and the
emitted electron, the stationary states
$\\{|\Psi^{-}_{\alpha,{\bm{k}}}\rangle\\}$, we can derive the cross-section
for the autoionisation (or Auger) process. In this process, thus, the initial
collision perturbs the system by promoting an electron from an inner-shell
occupied orbital of the ground state $\ket{0}$ to an empty resonance state
$|\Phi_{j}\rangle$, which decays non-radiatively into the channel
$|\chi^{-}_{\alpha,{\bm{k}}}\rangle$ representing a single (double) ionised
atom and an electron in the continuum. In particular, the autoionisation
cross-section is proportional to the square element of the following on-shell
$T$ matrix ($H_{1}$ is the transition operator from Eq. (122)):

${\bm{T}}^{+}_{0\rightarrow\alpha}({E})=\langle\chi^{-}_{\alpha,{\bm{k}}}|H_{1}|0\rangle+{\langle\chi^{-}_{\alpha,{\bm{k}}}|H-E|\Phi_{j}\rangle\langle\Phi_{j}^{-}|H_{1}|0\rangle\over{E}-{E_{l}}\delta_{lj}-\Delta_{lj}-i{\Gamma_{lj}\over
2}},$ (198)

which can be obtained using Eqs. (191) and (197) with the definition

$\ket{\Phi_{j}^{-}}=\left[|\Phi_{j}\rangle+\lim_{\nu\rightarrow
0}\sum_{\beta}^{m_{c}}\int\limits_{0}^{\infty}\frac{d\bm{p}}{2\pi^{3}}{|\chi^{-}_{\beta,{\bm{p}}}\rangle
M^{*}_{\beta j}({\bm{p}},E)\over E-E_{\beta}-\frac{p^{2}}{2}-i\nu}\right],$
(199)

which represents the state $\ket{\Phi_{j}}$ modified by a mixture of continuum
states. In Eq. (198), the first term represents the direct ionisation of the
atom in the solid (it represents a direct double ionisation in the Auger
process). In the second term, the transition operator connects the initial
ground state $|0\rangle$ with the final continuum state
$|\chi_{\alpha,\bm{k}}^{-}\rangle$ via the resonance state $\ket{\Phi_{j}}$.
Finally, the cross-section of the resonant autoionisation process is

$\displaystyle{\partial\sigma_{0\rightarrow\alpha}\over\partial{\bm{k}}}({\bm{k}};\omega)=\left({4\pi^{2}\omega\over
c}\right)|<0|H_{1}|\chi^{-}_{\alpha\vec{k}}>$ $\displaystyle+$
$\displaystyle\sum^{m_{c}}_{\beta=1}\lim_{\nu\rightarrow\infty}\int{d{\bm{p}}\over(2\pi)^{3}}{<0|H_{1}|\chi^{-}_{\beta{\bm{p}}}>{\overline{\Gamma}}_{\beta\alpha}({\bm{p}},{\bm{k}})\over
E-E_{\beta}-{p^{2}\over 2}-i\nu}$ $\displaystyle+$
$\displaystyle\sum^{m_{d}}_{jl=1}<0|H_{1}|\Phi_{j}>\Lambda_{jl}^{-1}<\Phi_{l}|H-E|\chi^{-}_{\alpha{\bm{k}}}>|^{2}\delta[(E_{0}+\omega)-(E_{\alpha}+{k^{2}\over
2})],$

where

$\Lambda_{lj}=(E-E_{l})\delta_{lj}-\Delta_{lj}-i{\Gamma_{lj}\over 2},$ (201)

and

$\displaystyle{\overline{\Gamma}}_{\beta\alpha}(\vec{k},\vec{p})$
$\displaystyle=$ $\displaystyle\sum^{m_{d}}_{jl}M^{+}_{\beta
j}(\vec{p},E)\Lambda_{jl}^{-1}M_{l\alpha}(E,{\bm{k}})=$ (202) $\displaystyle=$
$\displaystyle\sum^{m_{d}}_{jl}<\chi^{-}_{\beta{\bm{p}}}|H-E|\Phi_{j}>\Lambda_{jl}^{-1}<\Phi_{l}|H-E|\chi^{-}_{\alpha{\bm{k}}}>.$

The first two terms on the right-hand side of Eq. (LABEL:TF11e) describe the
direct ionisation process from the ground state $|0\rangle$ to the final decay
state $|\chi^{-}_{\alpha\vec{k}}\rangle$. The third term, on the other hand,
provides the resonance contribution, i.e. the contribution that can be
attributed to the decay process via bound intermediate states
$\\{|\Phi_{j}\rangle\\}$. In the case of an Auger process, the direct process
is a double ionisation process and therefore the first two terms on the right-
hand side of Eq. (LABEL:TF11e) contain matrix elements connecting a bound
state (the ground state) to the continuum states with two unbound electrons.
The third term, on the other hand, is characterised by matrix elements between
discrete and continuum states with only one unbound electron. The interference
between these different terms leads to asymmetric Fano profiles in the
electron spectrum, which are specific to the transition. Depending on the
kinetic energy of the projectile, the first two terms are usually smaller than
the third term, which can usually be neglected. Within this approximation, Eq.
(LABEL:TF11e) shows a Lorentzian profile (rather than a Fano profile)
characterised by the decay rate
$\Gamma_{\alpha}=\sum_{\beta}\Gamma_{\beta\alpha}$ into the channel
$|\alpha\rangle$.

An expression similar to that in Eq. (LABEL:TF11e) with
$\bra{0}H_{1}\ket{\chi^{-}_{\alpha\vec{k}}}$ replaced by
$\bra{0}H_{1}\ket{\chi^{-}_{\alpha{\bm{k}}}\eta_{\bm{k}_{1}}}$, and
furthermore, $\bra{0}H_{1}\ket{\Phi_{j}}$ substituted by
$\bra{0}H_{1}\ket{\Phi_{j}\eta_{\bm{k}_{1}}}$ and
$\delta[(E_{0}+\omega)-(E_{\alpha}+{k^{2}\over 2})]$ by
$\delta[(E_{0}+\omega)-(E_{\alpha}+{k^{2}+k_{1}^{2}\over 2m})]$ gives the
detailed expression of the cross-section for an Auger process. In an Auger
process, the global $N$-electron stationary state with the primary electron at
the energy $\epsilon_{1}={{k_{1}}^{2}\over 2}$ and momentum $\bm{k}_{1}$, can
be represented by the following antisymmetric product:

$|\Theta_{\alpha}({\bm{k}},{\bm{k}_{1}})\rangle=\hat{\cal
A}\left[|\Psi^{-}_{\alpha{\bm{k}}}\rangle|\eta_{\bm{k}_{1}}\rangle\right]$
(203)

where the wave function $|\eta_{\bm{k}_{1}}\rangle$ describes the primary
electron, which only interacts very weakly with the secondary Auger electron
when the kinetic energy of the projectile is much higher than the ionisation
energy of the inner-shell electrons.

In the same framework, the derivation of the cross-section can be generalised
to a molecular system in which an isolated discrete electronic state
supporting a range of vibrational levels is coupled to several electronic
continua (open channels) with their vibrational levels. In fact, in order to
be compared with the experimental data, the theoretical cross-section, Eq.
(LABEL:TF11e), must take into account the intrinsic properties of the target,
the finite resolution of the electron spectrometer and the broadening of the
autoionisation (or Auger) lines due to the energy losses of the emitted
electrons on their way out of the solid [39] before they leave the material
surface. In the latter case, the ab initio spectrum can be used as input for
an MC scheme [71, 166]. An application of this idea is developed in the
Results section. The cross-section in Eq. (LABEL:TF11e) can be determined by
calculating the resonance state and the continuum wave functions. To this end,
we use the static exchange approximation, which amounts to calculating the
bound intermediate states using the methods of quantum chemistry (e.g. HF, CI,
MBPT) [39] or condensed matter physics (e.g. DFT) [167]) and the scattering
wave function using the projected potential approach separately.

The interacting decay channels are represented by the antisymmetric normalised
product

$\chi_{\alpha\vec{k}}(1,..N-1,N)=\hat{\cal
A}[\Theta_{\alpha}(1,...,N-1)\psi_{\alpha\vec{k}}(N)]\,,$ (204)

where $\psi_{\alpha\vec{k}}(N)$ is the spin orbital of the emitted electron
and $\Theta_{\alpha}(1,...,N-1)$ is the Slater determinant representing the
bound state of the remaining $(N-1)$ system ($(N-2)$ for Auger processes),
which is the solution of the following secular problem

$\langle\Theta_{\alpha}|\hat{H}^{N-1}|\Theta_{\beta}\rangle=E_{\alpha}\delta_{\alpha\beta}\,,$
(205)

in a functional space spanned by the one-particle orbitals
$|\theta_{j}\rangle,~{}j=1,2,...,M$ with $M>>N-1$, which are expanded in a GBF
(see Eq. (91)).

To construct the continuum orbital in Eq. (204), which represents a particle
emitted in the effective field generated by the bound electrons of the system,
we use the method of projected potentials developed in section 5.2.
Specifically, we define a model Hamiltonian for an $N$-electron system as
follows:

$\hat{H}(1,..,N)=\sum^{N}_{i=1}[\hat{T}(i)+\hat{V}_{en}^{\pi}(i)]+{1\over
2}\sum^{N}_{i\neq j}\hat{v}^{\pi}(i,j),$ (206)

where

$\hat{T}(i)=-{1\over
2}\nabla^{2}_{i}~{},~{}~{}~{}\hat{V}_{en}^{\pi}(i)=\hat{\pi}(i)\hat{V}_{en}(i)\hat{\pi}(i),$
(207)
$\hat{v}^{\pi}(i,j)=\hat{\pi}(i)\hat{\pi}(j)\hat{v}(i,j)\hat{\pi}(i)\hat{\pi}(j),$
(208)

and $\hat{\pi}(i)=\sum_{l=1}^{m}~{}|g_{l}(i)\rangle\langle g_{l}(i)|$ is a
projector into an $m$-dimensional space spanned by the Gaussian functions and
$\hat{V}^{en}(i)=-\sum_{\mu}{Z_{\mu}\over|\vec{r}_{i}-\vec{R}_{\mu}|}$ and
$\hat{v}(i,j)={1\over|\vec{r}_{i}-\vec{r}_{j}|}$ are the electro-nuclear and
electron-electron Coulomb interaction potentials.

Instead of directly solving the secular problem of the Hamiltonian (206), the
continuum eigenfunctions can be found (as we have demonstrated in the case of
the elastic scattering in section 5.2) by using the LS with the appropriate
boundary conditions (see Eq. (92), valid for the elastic collisions):

$\psi^{-}_{\alpha{\bm{k}}}(\bm{r})=e^{i\bm{k}\cdot\bm{r}}+\hat{G}_{0}^{-}(\epsilon)\hat{T}_{\alpha}(\epsilon)e^{i\bm{k}\cdot\bm{r}},$
(209)

where
$G_{0}^{-}(\epsilon=\frac{k^{2}}{2})=\frac{1}{\epsilon-\frac{\hat{k}^{2}}{2}-i\varepsilon}$
is the Green’s function of the free particle and $\hat{T}_{\alpha}$ is the
transition operator given by the equation:

$\hat{T}_{\alpha}=\hat{V}^{\pi}_{\alpha}+\hat{V}^{\pi}_{\alpha}\hat{G}_{0}^{-}(\epsilon)\hat{T}_{\alpha}.$
(210)

The interacting channels defined in Eq. (204) are coupled by the following
interchannel interaction

$\displaystyle\langle\chi_{\alpha{\bm{k}}}|\hat{H}^{N}|\chi_{\beta{\bm{p}}}\rangle=(2\pi)^{3}\delta({\bm{k}}-{\bm{p}})\delta_{\alpha\beta}\left({k^{2}\over
2}+E_{\alpha}\right)+\langle\psi_{\alpha{\bm{k}}}|\hat{V}^{en}_{\pi}\delta_{\alpha\beta}+\hat{W}_{\pi}^{\alpha\beta}|\psi_{\beta{\bm{p}}}\rangle,$
(211)

where

$\hat{W}_{\pi}^{\alpha\beta}(1)=\sum_{j=2}^{N}\langle\Theta_{\alpha}(2,.j.,N)|\hat{v}_{\pi}(1,j)(\hat{1}-\hat{\cal
P}_{1,j})|\Theta_{\beta}(2,.j.,N)\rangle,$ (212)

and $\hat{\cal P}_{1,j}$ is the operator that determines the exchange of the
$(1,j)$ variables. In Eq. (212) $\Theta_{\alpha}$ are the solutions of the
bound-state problem of Eq. (205). The matrix elements defined in Eq. (211),
through which the channels are mixed via the effective interchannel
interaction, can be viewed as a representation of an effective one-particle
Hamiltonian over a set of basis vectors
$\\{|\varphi_{\alpha\vec{k}}\rangle\\}$ with internal degrees of freedom,
which moves in an effective potential depending on the internal states
$\\{|\alpha\rangle\\}$ of the particle. The diagonalisation of Eq. (211),
projected onto the GBFs, yields the wave functions for the different channels
that are open to decay and can be used in Eq. (LABEL:TF11e) to determine the
autoionization or Auger cross-sections. Finally, the spectrum of
autoionisation (Auger) can be plotted as a function of the kinetic energy of
the emitted electrons.

##  6 Results and Discussion

In this section, we present a number of applications of the previously
discussed methods for modelling elastic and inelastic collisions as well as MC
calculations of energy losses and secondary spectra.

###  6.1 Differential and total elastic scattering cross-sections

To test the accuracy of Mott’s cross-section (see section 5.1) in reproducing
the experimental results, we present the DESCS obtained using the Eq. (73) for
an electron beam with a kinetic energy of 1000 eV scattered by bulk aluminum
(Fig. 13) and gold (Fig. 14) as a function of the scattering angle measured
from the surface normal (see Fig. 1 for the experimental setup). In Fig. 15 we
show the same calculation for mercury at lower kinetic energy (300 eV). The
agreement with experimental data and other calculations shown by dots and
stars in the figures is excellent at all energies and scattering angles.

Figure 13: Differential elastic scattering cross-section calculated with
Mott’s formula from Eq. (73) (solid line) for an electron beam with a kinetic
energy of 1000 eV scattered by bulk aluminum as a function of the scattering
angle measured from the surface normal (see Fig. 1 for the experimental
setup). The dots represent the calculations from Ref. [103]. Figure 14:
Differential elastic scattering cross-section calculated with Mott’s formula
from Eq. (73) (solid line) for an electron beam with a kinetic energy of 1000
eV scattered by bulk gold as a function of the scattering angle measured from
the surface normal (see Fig. 1 for the setup of the experiment). The dots
represent the calculations from Ref. [103]. Figure 15: Differential elastic
scattering cross-section, calculated with Mott’s formula from Eq. (73) (solid
line) for an electron beam with a kinetic energy of 300 eV scattered by bulk
mercury as a function of the scattering angle measured from the surface normal
(see Fig. 1 for the setup of the experiment). The dots [168] and stars [169]
represent the experimental data. Figure 16: Differential elastic scattering
cross-section calculated with Mott’s formula from Eq. (81) (solid line) for an
electron beam with a kinetic energy of 1000 eV scattered by a water molecule
as a function of the scattering angle measured from the surface normal (see
Fig. 1 for a plot of the experimental setup). The dots represent experimental
data from Ref. [170]. Figure 17: Total elastic scattering cross-section
calculated using Mott’s theory and the formula of Eq. (74) (solid line) of
electrons scattered by a water molecule as a function of the electron energy.
The dots [170] and stars [171] represent experimental data.

Furthermore, with the help of Eq. (81), which generalises Mott’s atomic
formula to molecular systems, we can also plot in Fig. 16 the DESCS of a
single water molecule. Finally, in Fig. 17 the Mott’s ESCS (solid line) in a
large kinetic energy range, which is obtained by integrating the DESCS over
the elevation angle $\theta$ [from 0 to $\pi$] (scattering angle) and the
azimuth angle $\phi$ [from 0 to $2\pi$] (see Eq. (74) for a centre-symmetric
potential) at each energy is compared with the experimental data (dots and
stars). The ESCS is still calculated for a single water molecule.

Nevertheless, it may be important to increase the accuracy of the ESCS
calculation by including solid-state effects, i.e. the effects of the
environment surrounding the atom, e.g. in a solid, or the molecule, e.g. in
the gas phase or, more importantly, in the liquid phase. This task can be
accomplished by using the first-principles relativistic approach developed in
section 5.2. In this context, we plot in Figs. 18, 19 and 20 (see Ref. [51])
the calculation of the ESCS of water for a single molecule, a cluster of three
molecules and the Zundel cation (protonated water).

Figure 18: Theoretical elastic scattering cross-section (continuous
multicoloured line) of electrons hitting a single water molecule with
different kinetic energies compared to mixed experimental/theoretical results
from Ref. [172] (full blue squares with error bars), with the recommended
elastic (full purple circles) and total (green triangles) cross-sections from
Ref. [173]. The recommended elastic cross-section from Ref. [173] is based on
the results obtained in Ref. [174] up to 6 eV and by interpolation of the
cross-sections obtained with different methods [171]. The recommended total
cross-section for energies between 0.1 and 7 eV is determined by R-matrix
calculations [175, 176], for 7–50 eV it is based on experiments from Refs.
[177, 178, 179] and above 50 eV on experimental measurements from Ref. [180].
The cyan-coloured circles show the results obtained with the ELSEPA code
[101]. In the inset, we give the statistical error bar of the cross-section
compared to the total cross-section from Ref. [173]. The data in the $x$\- and
$y$-axes are given in $\log$-scale, $a$ is in bohr-2. Source: Reprinted from
[51]. Figure 19: ESCS of electrons hitting clusters of three water molecules
in different molecular configurations (3 [H2O]1,2, multicoloured dots) with
different kinetic energies compared to the results for the single water
molecule (1 H2O, multicoloured line). The cyan coloured dots show the values
calculated with the ELSEPA code [101]. The total cross-sections are normalised
to the number of molecules in the clusters. The data on the $x$\- and $y$-axis
are given in $\log$-scale, $a$ is in bohr-2. Source: Reprinted from [51].
Figure 20: ESCS of electrons colliding with the Zundel cation (multicoloured
empty squares) at different kinetic energies compared to water clusters (3
[H2O]1,2, multicoloured dots). The cyan coloured dots show the values
calculated with the ELSEPA code [101] for the Zundel molecule, where the
charge of the shared proton is compensated by an additional electron to
achieve charge neutrality. The total cross-section of the water clusters was
divided by three to normalise with the number of molecules in the clusters and
by two for the Zundel cation. The data in the $x$ and $y$ axes are in $\log$
scale, $a$ is in bohr-2. Source: Reprinted from [51].

We emphasise that these calculations are performed by freezing the roto-
vibrational motion. We also stress that a single set of GBFs for projecting
the interaction potential and wave functions unfortunately cannot accurately
describe the elastic scattering cross-section in the whole spectrum, i.e. one
has to include several GBFs of type $s$ with an energy-dependent width to
accurately describe the oscillations of the continuum wave function in the
whole energy range $[0,1000]$ eV. In Fig. 18, each line colour actually
represents a different Gaussian width (see the figure legend, where seven
widths were needed to reproduce the experimental data). The width of the
Gaussian functions is set to capture the oscillations of the scattering wave
functions at different energies. Especially at low energy, the GBFs must
capture the behaviour of the long wavelengths; therefore, the width of the
GBFs must be somehow small. At higher energies, larger values for the width
are required to capture the narrow oscillations of the continuum wave
function. Above 4 eV, the theoretical spectrum is in excellent agreement with
the experimental data [172], showing that the use of a theoretical method
based on the self-consistent calculation of the molecular orbitals, taking
into account polarisation and the permanent dipole, is effective. At higher
energies ($>$ 10 eV), the ab initio simulations agree with the data indicating
the total cross-section (green triangles in Fig. 18), which also contain
inelastic contributions, although they only take into account the elastic
channel. This can be explained by the fact that the total scattering
probability does not change by opening several channels, but is merely
redistributed between them (see Refs. [39, 40]). In Fig. 18 ab initio
calculations are compared with the Mott’s cross-section using a semi-empirical
treatment of the dipole moment [101] (cyan coloured dots), which clearly shows
the effect of the different treatment of the interaction potential between
these two approaches.

Furthermore, the projected potential method can also be successfully applied
to study polycentric systems in their interaction with an electron beam. In
this context, we show in Fig. 19 the ESCS from ab initio simulations (upper
multicoloured dots) for two clusters of three water molecules in different
configurations with molecular dipole moments pointing in different directions
compared to the calculations for a single water molecule (lower multicoloured
line), ELSEPA (cyan dots) and experimental data recorded for water vapour
under different conditions of relative humidity (symbols). This result shows
the effect of the environment, as the ESCS of the water clusters are
normalised to the number of molecules in the simulation cell. In particular,
we notice the similarity between the ESCS of the two different water clusters.
In the medium energy range ($1<E<30$ eV), the ESCS of the water clusters
differ by almost 100% with respect to the single molecule, showing a clear
influence of multiple scattering and the need of a more accurate
representation of the inter- and intramolecular potential provided by first-
principles simulations.

Finally, we show in fig. 20 the ESCS of the Zundel cation (upper multicoloured
empty squares) obtained with the projected potential relativistic ab initio
method compared to the ESCS of the water clusters. The influence of the
environment is also clearly visible in this plot, which once again emphasises
the importance of an accurate treatment of the interatomic potential.

###  6.2 Inelastic scattering: ELF and IMFP

The inelastic scattering collisions can be treated with different accuracy
within the dielectric theory by using either the DL/Mermin model (see section
5.5.2) or first-principles simulations (see section 5.5.3) to determine the
ELF. We recall that the dielectric theory is based on the evaluation of the
ELF in the entire Bethe surface.

Figure 21: Long-wavelength limit (${\bm{q}}\rightarrow 0$) of the ab initio ELF of Au (red curve) compared to experimental optical data from Ref. [181] (black curve). The LSDA for the exchange correlation functional was used to include the SO coupling in addition to the ALDA approximation for the time-dependent exchange correlation kernel. Source: Reprinted from [22]. peak label |  $\omega^{a}$ (eV) |  $\omega^{b}$ (eV) | This work | Exp. | Origin
---|---|---|---|---|---
1’ | 2.2 |  | 2.66 | 2.5 | ME
1 | 5.1 | 5.5 | 5.94 | 5.9 |  $5d$ plasmon
2’ | 10.2 | 11.5 | 11.49 | 11.9 | Mainly-$6s$ plasmon
2 | 15.5 | 18 | 15.67 | 15.8 | IT
3 | 23.8 | 25 | 23.74 | 23.6 | ME
4 | 30.8 | 35 | 31.35 | 31.5 | ME
5 | 36.9 | 40.5 | 37.5 | 39.5 | IT
6 | 43.5. | 43.5 | 44 | 44 | IT
7 |  |  | 59.2 |  | IT
8 |  |  | 65.62 |  | IT
Table 3: Optical ELF (${\bm{q}}\rightarrow 0$) of bulk Au. First column: Peak
designations according to Fig. 21. The second ($\omega^{a}$) and third
($\omega^{b}$) columns report the peak position in the ELF according to Refs.
[182] and [137]. In the fourth column we give the values obtained in this work
with the number of digits representative of the numerical accuracy, while in
the fifth column (Exp.) we give the experimental EELS data from [181]. In the
last column (Origin), the theoretical interpretation of each energy loss peak
is outlined. IT stands for an interband transition, while ME stands for a
mixed excitation. Source: Reprinted from [22].

Figure 22: Upper panel: ELF of Au for finite momentum transfers along the
(1,1,1) direction. The calculated $q$-vectors are indicated in the legend in
lattice coordinates (e.g. 0.0625 means the $q$-vector (0.0625, 0.0625,
0.0625)). Bottom panel: ELF of Au for finite momentum transfers along the
(1,1,0) direction. The calculated $q$-vectors are indicated in the legend in
lattice coordinates (e.g. 0.0625 means the $q$-vector (0.0625, 0.0625, 0)).
The LSDA approximation for the exchange correlation functional was used to
account for the SO coupling in addition to the ALDA approximation for the
time-dependent exchange correlation kernel. Source: Reprinted from [22].

In Fig. 21 we show the optical ELF of bulk gold (red line, see [22])
calculated with the ab initio approach described in section 5.5.3 compared to
the experimental optical data [181, 183] (black line). This shows excellent
agreement, although the surface plasmon excitations were neglected with
respect to the experimental EELS data. We note that another approach to
calculate the ELF is to use a Liouville-Lanczos method [184, 185], which
avoids the computationally intensive summation over empty bands. The latter
approach was described in Ref. [182] and achieved comparable results.

An interpretation of the spectral features of the ELF lineshape can be found
in Table 3 and a detailed discussion can be found in Ref. [22]. We only
emphasise here that ab initio methods based on DFT are able to reproduce the
well-resolved double-peak structure due to the localised $d-$bands that are a
feature of other transition metals. In Au, the presence of these localised
orbitals acts as a dielectric background and is responsible for both the
lowering of the $5d$ plasmon energy from the free electron plasma frequency
($\simeq$9 eV) to 5.94 eV (peak 1 in Fig. 21) and the well-resolved $5d$ and
$6s$ valence electron peaks.

$i$ |  $A_{i}$ (eV2) |  $\gamma_{i}$ (eV) |  $\omega_{i}$ (eV)
---|---|---|---
1 | 1.969 | 3.463 | 5.382
2 | 17.508 | 6.121 | 11.644
3 | 26.815 | 4.948 | 15.909
4 | 106.586 | 5.784 | 24.449
5 | 43.843 | 3.086 | 30.993
6 | 98.402 | 5.881 | 33.355
7 | 90.344 | 5.670 | 37.938
8 | 177.047 | 7.852 | 44.489
9 | 551.321 | 18.326 | 59.799
10 | 29.773 | 3.667 | 68.164
11 | 128.517 | 7.449 | 74.020
12 | 65.4861 | 6.292 | 87.151
13 | 46.896 | 7.192 | 95.522
14 | 2500.0 | 350.0 | 330.0
15 | 2000.0 | 2000.0 | 2500.0
16 | 200.0 | 5000.0 | 14000.0
Table 4: Fitting parameters of Eq. (133) obtained for the optical ELF
(${\bm{q}}\rightarrow 0$) in the bulk. Source: Reprinted from [22].

Fig. 22 shows the peak dispersion for finite momentum transfers along two
symmetry directions (1,1,1) (upper panel of Fig. 22) and (1,1,0) (lower panel
of Fig. 22). The dispersionless characteristic of the double peak indicates
that it is formed by interband transitions and not by plasmon excitation.
Other peaks at 2.66 eV (peak 1’ in Fig. 21) and at 15.67 eV (peak 2 in Fig.
21) show the same dispersionless behaviour with finite momentum transfer,
while the peaks at 5.94 and 11.49 eV (peaks 1 and 2’ in Fig. 21), which are
5$d$ and 6$s$ plasmon signals, show a clear momentum dispersion. The peaks at
higher energies (5, 6, 7 and 8 in Fig. 21) are also weakly dispersed,
indicating their nature as interband transitions.

Figure 23: Ab initio ELF extended to high excitation energies [186] (green
line) together with the DL fit (purple line). The matching point between ab
initio and experimental data was set to the threshold energy of the semi-core
transitions $E_{5p1/2}=74.2\ $ eV [187]. The fit parameters are shown in Table
4. Source: Reprinted from [22].

We also emphasise that (see section 5.5.3) the first-principles simulations
are very accurate, but at the price of a high computational cost that makes it
impossible to extend the energy range and take into account the energy-
momentum dispersion beyond a few tens of eV. However, in order to perform MC
simulations for charged particle beams with high kinetic energy, it is
necessary to determine the inelastic scattering cross-section far beyond these
limits [32, 48, 49, 50, 67].

Figure 24: IMFP (solid line [22]) as a function of incident electron energy in
Au. Dashed line from Ref. [188].

In this regard, the extension of the optical ELF along the excitation energy
axis can be performed by matching the ab initio results with the experimental
data, which typically cover a larger energy. In particular, in the NIST
database [186] one can find X-ray form factors, attenuations and scattering
tables from which the ELF up to 30 keV can be derived. In Fig. 23 we show such
an extension for Au, where the matching point between ab initio and
experimental data was set to the threshold energy of the semi-core transitions
$E_{5p1/2}=74.2\ $ eV [187]. In Table 4 we also show the DL parameters of Eq.
(133), which were determined by a best-fit procedure.

Finally, to determine the IMFP the ab initio data must also be extended beyond
the optical limit for finite momentum transfers, e.g. by choosing one of the
dispersion laws described in section 5.5.2. The IMFP of Au can then be
determined using the formula (141) and is shown in Fig. 24 for a quadratic
dispersion (blue solid line) as by Eq. (139).

Figure 25: Elf of bulk a) Cu, b) Ag and c) Au in the optical limit
($q\longrightarrow 0$), obtained by best fitting the experimental data of zero
momentum transfer with the DL model of Eq. (133). These data can be found in
[189] for Ag, in [190] and [191] for Cu and Au respectively. Source: Reprinted
from [67]. Figure 26: $Z_{\mathrm{eff}}$ (left panel) and $P_{\mathrm{eff}}$
(right panel) versus the transferred energy ($\Delta E_{\mathrm{max}}$) for
Cu, Ag and Au. Source: Reprinted from [67].

As a further application of the DL model, we show in Fig. 25 the bulk ELFs in
the optical limit of Cu, Ag and Au obtained by best fitting the experimental
data for zero momentum transfer from Ref. [191] for Au, Ref. [190] for Cu and
Ref. [189] for Ag. In Table 5 we also give the best fit parameters for fitting
the experimental optical ELF by DL functions.

|  | Cu |  |  |  | Ag |  |  |  | Au |
---|---|---|---|---|---|---|---|---|---|---|---
i |  $\omega_{i}$ (eV) |  $\gamma_{i}$ (eV) |  $A_{i}$ (eV2) | i |  $\omega_{i}$ (eV) |  $\gamma_{i}$ (eV) |  $A_{i}$ (eV2) | i |  $\omega_{i}$ (eV) |  $\gamma_{i}$ (eV) |  $A_{i}$ (eV2)
1 | 4.08 | 1.09 | 0.33 | 1 | 7.89 | 3.37 | 12.80 | 1 | 9.52 | 14.97 | 18.49
2 | 10.07 | 5.99 | 22.10 | 2 | 38.20 | 42.93 | 1109.46 | 2 | 15.92 | 6.26 | 25.85
3 | 19.05 | 8.16 | 88.91 | 3 | 59.58 | 29.93 | 480.38 | 3 | 25.58 | 2.18 | 11.12
4 | 27.21 | 8.16 | 112.54 | 4 | 73.81 | 20.12 | 300.6 | 4 | 38.09 | 26.67 | 973.52
5 | 78.91 | 152.38 | 2216.74 | 5 | 85.70 | 27.70 | 226.83 | 5 | 64.49 | 30.48 | 507.39
|  |  |  |  |  |  |  | 6 | 99.32 | 19.05 | 88.88
|  |  |  |  |  |  |  | 7 | 402.71 | 612.23 | 337.32
Table 5: Fitting parameters of Eq. (133) for the optical ELF (${\bm{q}}\rightarrow 0$) in the bulk. The best fitting parameters for Cu were provided by C.C. Montanari et al. [190], while for Au by C.D. Denton et al. [191]. In the case of Ag, the parameters were determined by best-fitting the optical measurements of Smith et al. [189]. Source: Reprinted from [67]. |  Cu |  |  Ag |  |  Au |
---|---|---|---|---|---|---
$E$ (eV) | $Z_{\mathrm{eff}}$ | $P_{\mathrm{eff}}$ | $Z_{\mathrm{eff}}$ | $P_{\mathrm{eff}}$ | $Z_{\mathrm{eff}}$ | $P_{\mathrm{eff}}$
1 | 5.52e-05 | 0.017 | 9.24e-05 | 0.020 | 1.29e-04 | 0.028
10 | 0.09 | 0.275 | 0.17 | 0.336 | 0.1 | 0.27
108.5 | 7.28 | 0.950 | 18.63 | 1.056 | 16.65 | 1.021
1012 | 19.68 | 0.987 | 36.36 | 1.067 | 44.17 | 1.046
9886.4 | 27.07 | 0.987 | 45.46 | 1.067 | 70.74 | 1.047
29779 | 28.23 | 0.987 | 46.51 | 1.067 | 78.94 | 1.047
Table 6: $Z_{\mathrm{eff}}$, $P_{\mathrm{eff}}$ limit values of Cu, Ag and Au
represented in Fig. (26). Source: Reprinted from [67].

To test the accuracy of the DL model, we compared this data with the $f$-sum
rule and perfect-screening (ps-sum or $P_{\mathrm{eff}}$) sum rule (see Eqs.
(137) and (138) respectively). In Fig. 26 we show the $Z_{\mathrm{eff}}$ (left
panel) for the ELF of Cu, Ag and Au, while in the right panel we show the
$P_{\mathrm{eff}}$. The limit values $Z_{\mathrm{eff}}$, $P_{\mathrm{eff}}$
are also reported in the Table 6.

(a) IMFP of Cu.

(b) IMFP of Ag.

(a) IMFP of Au. Figure 28: IMFP of Cu, Ag and Au with (dashed curves) and
without (solid curves) the Born–Ochkur electron exchange corrections (see Eq.
(142)) compared to simulations by Tanuma et al. [188] (blue curve). The IMFPs
are calculated using the bulk ELFs at high $q$ limit (black curves), the full
3D bulk ELF (red curves, RPA, see Eq. (139)) and the effective ELF with the
bulk surface extension to finite momenta (green curves, BS, see Eq. (140)),
respectively. Abscissa and ordinate are given in logarithmic scale to
emphasise the difference between the different treatments of the
$q$-dispersion. Source: Reprinted from [67].

To show how the use of different dispersion laws (see section 5.5.2 for the
mathematical details) and the addition of the electron exchange correction
(see Eqs. (141) and (142)) affect the calculated inelastic scattering
quantities, we plot the IMFPs of the three metals Cu, Ag and Au in Fig. 28
[67]. We find that the discrepancy between the dispersion laws is significant
below 100 eV, while the inclusion of exchange in the calculation of the IMFP
becomes significant below 10 eV.

Figure 29: IMFP as a function of incident electron energy in PMMA. Solid line
from Ref. [192], dots from Ref. [57] and stars from Ref. [188]. Figure 30:
IMFP as a function of incident electron energy in liquid water [32]. Solid
line: TDDFT calculations. Dashed line: MELF-GOS. Figure 31: IMFP as a function
of the energy of the incident electrons in CeO2. Solid line: TDDFT simulations
[48]. Dots: experimental data using the Ni standard [193]. Stars: experimental
data using the Au standard [193].

Further examples of IMFPs of organic molecular materials, such as polymethyl
methacrylate (PMMA) [194], water [32] and inorganic solids, such as cerium
oxide [48], using the DL, MELF-GOS and TDDFT methods are shown in the Figs.
29, 30, and 31, respectively. A comprehensive collection of experimental and
calculated data on IMFPs for different combinations of projectiles and
materials can be found in the recently published Ref. [195].

###  6.3 Auger electron spectra of nanostructures

The ab initio method developed in section 5.7 can be used to calculate Auger
spectra in molecules and solids. The spectrum of the non-radiative O $K-LL$
decay of silicon dioxide is shown in Fig. 32 [39, 71, 166]. We find that the
theoretical spectrum (orange line) accurately reflects the experimental
transition energies (blue line), which are mainly distributed in three
regions, while the intensities are lower than the measured values. In
particular, the most intense peaks are found in the $[494-510]$ eV region,
corresponding to $2p-2p$ double holes in the final states of the oxygen atom
($K-L_{23}L_{23}$ transitions), with the maximum intensity at $502.76$ eV,
which is due to a singlet transition, in good agreement with the experimental
value of $502.58$ eV; between $[475-492]$ eV a group of features corresponding
to Auger emissions with final double-hole configurations $2s-2p$
($K-L_{1}L_{23}$), with the main peak corresponding to the singlet transition
at $481.1$ eV, in agreement with the experimental peak at $481.38$ eV; finally
between $[460-470]$ eV, with the final double-hole configuration in the
$2s-2s$ orbitals of oxygen ($K-L_{1}L_{1}$) and the strongest peak at 463.86
eV, which corresponds to the experimental value of 465.25 eV. In general, the
ab initio method based on the theory of projected potentials to calculate the
continuum wave functions shows that it is a valid tool to reproduce the
transition energy of Auger decay.

Nevertheless, Auger electrons emitted from atomic centres in solids find an
environment that can change both their spectral intensity and their
distribution, e.g. due to inelastic scattering events. In this respect, the MC
method can be used to account for the energy losses that Auger electrons
suffer while travelling in the solid and before leaving the surface.

In this context, one can use the ab initio electron distribution (orange line)
as input for an MC routine to account for such changes in the spectrum caused
by energy losses of the electrons on their way out of the solid. Using this
approach and assuming that Auger generation occurs at a depth of 40 Å [196],
the Auger line profile that includes inelastic events (green line in Fig. 32)
can be determined. We find that such inclusion increases and broadens the
Auger decay probability so that the measured experimental results are
satisfactorily reproduced over the entire energy range.

Figure 32: O $K-LL$ Auger spectrum in SiO2. Comparison between the ab initio
(orange line) and MC (green line) calculations with experimental data (blue
line) [39, 71, 166].

In particular, there is a strong broadening of the $K-L_{1}L_{23}$ peak, which
is due to the SiO2 plasmon excitation at $\approx$23 eV. Finally, the Auger
spectrum can also be determined in molecular systems by using the same
approach of projected potentials [39], e.g. the Auger spectrum of ozone was
recently calculated [72].

###  6.4 Spectra of backscattered electrons, secondary electrons and yields

Figure 33: Monte Carlo simulation of the energy loss spectrum of electrons
from an Au sample for a kinetic energy of 1000 eV of the primary beam. Figure
34: Monte Carlo simulation of the REEL energy region of an Au sample for a
kinetic energy of the primary beam of 1000 eV, showing the plasmon and part of
the elastic peak.

Knowledge of the elastic and inelastic scattering cross-sections, which can be
derived from the energy- and momentum-dependent ELF of a particular sample,
enables the MC simulation in the entire energy spectrum, taking into account
all processes such as single-electron and collective excitation, ionisation
and elastic collisions. Each energy loss mechanism and each generation of
secondary electrons leads to a well-resolved peak. In Fig. 33 we show such a
MC-simulated spectrum of both backscattered and secondary electrons emerging
from Au. From right to left of the spectrum, one finds: the elastic – or zero-
loss – peak, characterised by a width at half maximum of a few decimal eV to
account for the intrinsic width of the beam (typically 0.8 eV); the bulk
plasmon peaks, zoomed in the REEL spectrum of Fig. 34, at $\approx$ 5.77 and
11.6 eV from the elastic peak, slightly shifted with respect to the maxima of
the ELF (5.94 and 11.49 eV, see also Table 3); finally, a series of clearly
resolved spectral features (see also Fig. 34 for an enlarged view) in the
range 15-50 eV from the elastic peak at $\approx$ 15.7, 24.7, 31.7 and 44 eV,
which are due to $5d,4f$ $\rightarrow$ 6s interband electronic transitions
mixed with plasmon-like features (Re$(\bar{\varepsilon})\simeq 0$,
Im$(\bar{\varepsilon})$ small but $\neq 0$) at 24.7 and 31.7 eV. To generate
these spectra, a number of electron trajectories of more than 109 was used.

Figure 35: Comparison between the theoretical REEL spectra of bulk CeO2 (solid
line) and Ce2O3 (dashed line), simulated with the Monte Carlo method [48], and
the measured values (stars) [197]. The spectra are normalised with respect to
the area. The primary beam has a kinetic energy of 1600 eV and an angle of
incidence of 60∘ with respect to the surface normal. Figure 36: REEL spectrum
for 1000 eV electrons impinging on Al as a function of energy loss (the zero-
loss peak is shown on the far left of the spectrum). The primary beam strikes
the material with an angle of incidence of 30∘ with respect to the surface
normal. The spectra were normalised to a common area of the zero-loss peak.
Plasmon peaks up to the 5th order of scattering can be observed in both the
experimental and simulated spectra. Please note the presence of surface
plasmons and bulk plasmons in the spectra. Solid line: Monte Carlo simulations
[198]. Dashed line: experimental data from [199]. Figure 37: REEL spectrum for
1000 eV electrons impinging on Si as a function of energy loss (the zero-loss
peak is shown on the far left of the spectrum). The primary beam strikes the
material with an angle of incidence of 30∘ with respect to the surface normal.
The spectra are normalised to a common area of the zero-loss peak. Plasmon
peaks up to the 4th order of scattering can be observed in both the
experimental and simulated spectra. Please note the presence of surface
plasmons and bulk plasmons in the spectra. Solid line: Monte Carlo simulations
[44]. Dashed line: experimental data [199].

To show the general applicability of the TMC approach to materials with
different electronic properties, further examples of calculated REEL spectra
for different kinetic energies of the incident electron beam are shown in Fig.
35, 36 and 37 for cerium oxides [48], which are insulating materials, metallic
aluminium [44] and silicon semiconductors [44] in comparison to experimental
results. The spectral behaviour reflects the presence of single and plasmon
excitations for all samples. To accurately reproduce the experimental results
and reduce statistical noise, a number of 109 trajectories were used. In
particular, we note the presence of surface and bulk plasmon peaks in both the
theoretical and experimental spectra of Al and Si. In addition, several
equally spaced plasmon peaks can be seen in Figs. 36 and 37, the intensity of
which decreases with increasing energy loss, in the case of Al up to the 5th
energy level and in the case of Si up to the 4th energy level.

Figure 38: Monte Carlo simulation of the spectral range of the secondary
electron emission of an Au sample for a kinetic energy of the primary beam of
1000 eV, which shows a wiggling behaviour.

The energy straggling strategy can also be used in TMC for quantitative
modelling of the secondary electron spectra and yields (SEY) of samples, such
as Al [200], Si [200, 201, 202], Cu [148, 202], Ag [148, 200], Au [148, 22,
200], PMMA [194], P3HT [118], Al2O3 [65] and SiO2 [203]. Essentially, this
amounts to integrating the spectra of the secondary energy distribution (see
e.g. Fig. 33) between 0 and 50 eV by changing the incident kinetic energies of
the electron beam. Recall that these so-called secondary electrons are, by
definition, the particles produced by inelastic scattering events
characterised by energy losses above the atomic ionisation threshold and
ejected with an energy of less than 50 eV. Of course, they must also overcome
the energy barrier at the solid-vacuum interface, i.e. the work function (see
Eq. (38)). In this regard, there are also primary backscattered electrons that
leave the target with a kinetic energy of less than 50 eV after suffering
strong energy losses; however, their number is typically negligible compared
to the number of secondary electrons that are emitted with an energy of less
than 50 eV.

In Fig. 38 we show the MC simulation of the secondary electrons emitted by an
Au sample for a primary beam with a kinetic energy of 1000 eV in the energy
range where most of the secondary electrons are emitted. The wiggling
behaviour of the lineshape reflects the ELF peaks. We also recall that the
conditions for the emission of electrons at the sample surface are met if the
transmission coefficient, which depends on the work function of the material,
is greater than a random number (see Eq. (39)). The work function of Au
depends strongly on the type of impurities present on the surface, the
duration of exposure of the sample to air and the sputtering processes used to
clean it [204, 205]: here we have assumed a value of $\chi=4.7$ eV, which can
be obtained after $\mathrm{Ar}^{+}$ sputter-cleaned Au.

Figure 39: Comparison between MC calculations of the SEY of Ar+-sputtered Au
samples for different kinetic energies of the primary beam and the
experimental data [206]. The calculations were performed for normal incidence.
Figure 40: Angular distribution of the SE emerging from the Au surface for a
1000 eV electron beam.

In Fig. 39 we show a comparison between the SEY of Au simulated by the MC
(blue line) [22] and the experimental data (orange dashed line) [206] as a
function of the kinetic energy of the primary beam at perpendicular incidence
with respect to the target surface. The agreement is remarkable. We observe
that the SEY increases with increasing primary energy up to a maximum and then
decreases again. This behaviour can be explained by the fact that at a very
low primary energy only a few secondary electrons can be generated and thus
emitted from the surface; this number naturally increases with the primary
energy together with the average depth at which the secondary electrons are
generated. Finally, with increasing energy, the point at which the secondary
electrons are generated is so deep inside the material that hardly any can
escape from the solid surface.

The angular emission pattern of the secondary electrons is also shown in Fig.
40 for the same conditions (1 keV electron beam hitting an Au sample) [22]. We
note that the angular distribution closely follows a cosine function that is
almost inversion symmetric around the maximum at $\approx$ 45°.

Figure 41: Secondary electron yield of PMMA. The calculations were performed
for normal incidence. Solid line: MC simulations [194]. Dots: experimental
data [207]. Figure 42: Total electron yield (= secondary electron yield +
backscattering coefficient) (TEY) of PMMA. Calculations were performed for
normal incidence. Continuous line: MC simulations [194]. Dots: Boubaya and
Blaise experimental data [208]. Stars: Rau et al. experimental data [209].
Figure 43: Ionisation cross-sections of carbon ions in water. The solid and
dashed lines correspond to the calculations performed with the LR-TDDFT (see
section 5.5.3) and MELF-GOS (see section 5.5.2) approaches. The symbols show
the experimental data for water vapour [210]. (a) Energy distributions of
electrons emitted during collisions of carbon ions with water molecules as a
function of their kinetic energy $W$ for different carbon incident energies
$T$. (b) Angular distributions of the electrons emitted in the collision of 6
MeV/u C ions as a function of the ejection angle $\theta_{2}$ for different
values of the emission energy $W$. (c) Angular distribution of electrons
emitted with a kinetic energy of 100 eV that were created in liquid water by
the collision with C ions for different energies $T$. Source: Reprinted from
[32, 33].  Figure 44: Energy density radially deposited by secondary electrons
from a 2 MeV/u carbon ion beam in H2O. Continuous line from Ref. [33]. Dashed
line from Ref. [211]. Dots: Geant4-DNA [212]. Symbols: $\times$ from Ref.
[213], stars from Ref. [214]. Figure 45: Average cluster size in a sensitive
volume of liquid water as a function of the impact parameter from the ionic
pathway is shown for two kinetic energies [32, 33]: (a) $T$ = 0.2 MeV/u and
(b) $T$= 1 GeV. The symbols correspond to the clusters produced by different
types of events: Ionisations (red triangles), excitations (blue circles),
dissociative excitations (blue hollow circles) and dissociative electron
attachments (DEA) (green diamond). The total average cluster size, which
includes all events, is represented by black squares, while the size that
includes only the damaging events is represented by grey squares. The letters
in the inset of panel (a) represent the sensitive cylindrical volume and refer
to excitation (e), ionisation (i), and DEA (d). Source: Reprinted from [32,
33]. Figure 46: Relative contribution of ionisation, dissociative excitation
and DEA to the average cluster size for different kinetic energies of the
carbon ions as a function of the radial distance of the test volume,
represented by a cylinder with a dimension of two DNA turns. Source: Reprinted
from [32, 33].

Further examples of the MC calculations of SEY and TEY for the organic
material PMMA are shown in Figs. 41 and 42 respectively. The behaviour is
similar to that previously observed for Au, i.e. the yield of secondary
electron emission (which for TEY also includes the backscattered electrons of
the primary beam that escape with a kinetic energy ¡ 50 eV) increases until a
maximum is reached with increasing primary energy; thereafter, the yield
decreases with increasing primary energy.

###  6.5 The radial dose in liquid water

By accessing the elastic and inelastic scattering cross-sections, event-by-
event MC simulations can also be performed to model secondary electron energy
deposition [32, 33, 36, 44]. Especially for the development of novel and more
effective radiotherapies to cure cancer, it is of utmost importance to
investigate how charged particle beams such as protons and carbon ions
interact with the biological material and deposit their energy. This study is
also important for radiation protection and to prevent cosmic rays, which
consist mainly of protons and nuclei, from damaging electronic equipment on
board and harming personnel on manned space missions. In the specific case of
carbon ion beams, their successful use in hadron therapy is related to the low
angular deflection they suffer on their path through the biological medium, as
well as the well-defined penetration region characterised by a sharp peak of
the depth-dose curve (known as Bragg peak), which means that a large part of
their energy is deposited towards the end of their trajectory [215, 216, 217].
In addition, carbon ions have an increased relative biological effectiveness
(RBE) [218], i.e. the generation of numerous low-energy secondary electrons
and chemically reactive species increases their ability to kill cells more
effectively than photons, for example, for the same delivered dose [26, 215].

This property, which depends on the stopping power, enables precise focussing
on the tumour area, which can be treated without severely damaging the
surrounding healthy tissues. We note here that a relevant part of the damage
inflicted on cancer cells is mainly due to the secondary electrons generated
by the carbon ions on their way through the biological medium, which in
computational biology is mimicked by liquid water, as it represents 60% of the
human body weight. The secondary electrons generated by the collisions between
ions and atomic targets undergo elastic and inelastic processes that are
actively involved in the generation of further secondary electrons, excitation
of phonons, dissociative electron attachment (DEA) [38, 219] and trapping. The
relative importance of these phenomena in terms of damage to the cells depends
on the energy of the impacting electron. In order to evaluate the dose
deposited radially in the material along the ion track, the first and most
important mechanism to be simulated, which triggers the chemical-physical
mechanisms behind the energy transfer into the material, is thus the
generation of secondary electrons.

In this context, MC calculations for the specific kinetic energies of carbon
ions [32, 34] must first be performed. This calculation will be carried out a
number of statistically significant times with different random seeds to
obtain the structure of the trajectory, typically more than 1000. The length
of the ion track depends on the energy. It must be long enough for virtually
all secondary electrons generated along the carbon ion track to reach a
sensitive volume (with the dimensions of a DNA-like target) located at a
certain radial distance from the trajectory. Using dielectric theory, each ion
shot provides precise ionisation sites, the emission energy of the secondary
electrons and their ejection angle.

To achieve an acceptable compromise between computational cost and low signal-
to-noise ratio, the carbon ions should undergo a number of collisions with the
target atoms that can generate on average $10^{5}$–$10^{6}$ electrons, which
in turn produce an average 100 more electrons each. The avalanche of secondary
electrons can easily lead to the simulation of billion trajectories per ion
energy. We emphasise that these values depend crucially on the charge state of
the carbon ion, which can change dynamically during the motion to a fully
ionised state. The charge state also has a strong influence on the stopping
power, especially at energies around the maximum. The Bragg peak of the carbon
ion is also influenced by fragmentation into nuclear fragments. In Fig. 43 we
show the results of MC simulations to determine the ionisation cross-section
(panel (a)) of carbon ions in water as well as the energy and angular
distributions of the electrons emitted after the collision (panels (b) and
(c)). Once the energy and angular distributions of the secondary electrons
generated by the impact of carbon ions on liquid water and the relevant
elastic (see section 5.1 or 5.2) and inelastic (see section 5.3) scattering
cross-sections for electron transport in water are known, the simulation of
the dose deposited by carbon ions and the relevant mechanisms triggered by
their passage through the biological medium can be accurately performed.

In particular, we have access to: (i) the radial dose resulting from energy
deposition around the ion’s path and (ii) the clustering of damaging events
generated in DNA-sized volumes located at different radial distances from the
ion’s path. While the former is a macroscopic quantity that can be useful in
semi-empirical radiobiological models [30], the latter provides fundamental
information for the assessment of radiobiological effectiveness (RBE) for ion-
induced radiation damage [24, 27].

In general, the absorbed dose represents the average energy deposited per unit
mass by ionising radiation. The radial dose is an approximation of this. It is
estimated microscopically using concentric rings centred around the ion
trajectory with increasing radius. Using a MC algorithm [36] with the
previously calculated elastic and inelastic scattering cross-sections, one can
calculate the amount of energy deposited in cylinders of different radii by
secondary electrons divided by the mass of such an element. In Fig. 44 we show
the radial dose around carbon ions impinging on liquid water with a kinetic
energy of 2 MeV/u [33].

The average cluster size of damaging events on the nanoscale takes into
account processes such as ionisation, electronic excitation and DEA. To
account for DEA, the cross-section for liquid water was taken from
experimental data [171]. The average cluster size is commonly used to
represent the distribution of damage events generated in nanovolumes. These
are modelled by a cylinder of 2.3 nm (diameter) and 6.8 nm (height) positioned
with the axis perpendicular to the ion path at different impact parameters to
mimic sensitive biological targets that resemble the dimensions of two turns
of DNA molecules (a strand of 20 base pairs in length). In Fig. 45(a) and (b)
we show the average cluster size for impact parameters from 0 to 100 nm for
0.2 MeV/u and 1 GeV carbon ions, respectively.

First, we find that at shorter radial distances the total damage is larger for
0.2 MeV/u carbon ions, a tendency that weakens at larger distances. Second, we
find that ionisation processes (triangles in Fig. 45) dominate for all impact
parameters, with electronic excitation (full circles) having a significant
impact on causing severe damage, although only $\approx$40% of electronic
excitations are effective in dissociating water. Conversely, the DEA plays a
subordinate role as it is almost two orders of magnitude smaller than the
ionisation clusters.

In Fig. 46 finally, we plot the relative contribution to the average size of
damage clusters produced by ionisation events, dissociative electronic
excitations and DEA, respectively, for different carbon ion energies as a
function of the impact parameter. We find that, in general, ionisation
processes are responsible for $\approx$80% of the cluster size, followed by
dissociative excitations ($\approx$15-17%) and DEA $\approx$5%. Thus,
ionisation makes the largest contribution to the clusters of damaging events
triggered by carbon ions in liquid water, while DEA, which is normally
considered a very relevant biological damage mechanism [220], surprisingly
plays a minor role.

##  7 Summary, conclusions and outlook

This article examines the mechanisms of electron scattering, including elastic
and inelastic interactions with target atoms, molecules and solids, by
reviewing in detail the theoretical and computational methods that can be used
to interpret the experimental results. In particular, the application of
quantum mechanics techniques to the analysis of the interaction processes of
charged particles with various forms of matter is thoroughly explored. The
text also introduces the Monte Carlo method and provides a large number of
examples and extensive comparisons between simulated results and experimental
data from the literature. Modern derivations are discussed and ab initio
methods are explained in detail, providing a rigorous approach to the
calculation of elastic and inelastic scattering cross sections. Aspects such
as the energy and angular distributions of secondary electrons, the
backscattering coefficient, the yield of secondary electrons and their role in
hadron therapy are discussed.

In particular, Monte Carlo results on the radial distribution of energy
deposited in solids by secondary electrons generated by high-energy ion beams
are presented. The macroscopic nature of the radial dose resulting from the
energy deposition around the ion’s trajectory is emphasised and its usefulness
in semi-empirical radiobiological models is demonstrated. Furthermore, the
presented Monte Carlo results on the clustering of damage events within DNA-
sized volumes at different radial distances from the ion’s trajectory provide
fundamental information for the evaluation of radiobiological effectiveness in
the assessment of ion-induced radiation damage.

In this context, we point out that despite the maturity of the field of energy
loss spectroscopy and backscattered electron analysis and the thorough
understanding of electron transport and emission mechanisms, there are still
approximations, assumptions and challenges to overcome.

We believe, in particular, that future experimental research in this area
should focus on key aspects of fabrication, e.g. developing instruments such
as detectors that can collect more secondary electrons. This will enable high-
resolution imaging and improve the signal-to-noise ratio for better
visualisation of surface morphology and fine topographical details.

In addition, advances in the miniaturisation of devices at the nanoscale and
the ability to non-destructively analyse materials to extract energy and
momentum information from low-energy electrons in SEM and helium ion scanning
microscopes (HIM) are essential for improved characterisation of
nanomaterials. In particular, efforts should be made to develop non-
destructive microscopy techniques for imaging and mapping quantum electronic
properties and band structure information in nanostructures. Indeed, obtaining
electronic band structure information from SEM devices is still in its infancy
compared to other electron spectroscopy techniques such as angle-resolved
photoemission electron spectroscopy (ARPES).

On the theoretical side, there is a clear need to improve the integration of
ab initio methods into Monte Carlo simulations for a spectrally resolved study
of the generation and emission of secondary electrons in response to
excitation by an electron/ion primary beam, as can be observed in SEM or HIM.
Specifically, the aim is to improve the calculations of differential and total
elastic scattering cross sections for solids using first principles methods in
order to replace Mott’s theory. In particular, current methods are inadequate
to account for solid-state effects such as multiple scattering from
neighbouring atomic centres, especially in the low energy range (¡ 10 eV). In
addition, a relativistic dielectric theory that goes beyond the limitations of
the Ritchie approach is needed to determine the inelastic mean free path, the
differential inverse inelastic mean free path and the cross sections to
describe the interaction between electrons and ions with matter.

Theoretical models should also aim to establish a quantitative and robust link
between the band structure, the density of states (DOS) and the information
obtained from the energy spectra of secondary electrons in solids, which is
still lacking. Furthermore, extending the energy range of secondary electrons
to a few tens of electron volts could reveal additional features related to
plasmon and core-level energy loss and deliver information on conduction bands
that may help to interpret experimental results. These advances would
contribute enormously to the development of a comprehensive theory for the
generation of secondary electrons produced by the interaction between electron
and ion beams in solids.

##  8 Declaration of competing interest

The authors declare that they have no known competing financial interests or
personal relationships that could have appeared to influence the work reported
in this paper.

##  9 Acknowledgments

This action has received funding from the European Union under grant agreement
No 101046651. The authors gratefully acknowledge discussions with I. Abril, P.
de Vera, G. Garberoglio, R. Garcia-Molina, T. Morresi, S. Simonucci and P.
Trevisanutto.

## References

  * Werner [2023] W. Werner,  Electron beams near surfaces: The concept of partial intensities for surface analysis and perspective on the low energy regime,  Front. Mater. 10 (2023) 1202456. doi:https://doi.org/10.3389/fmats.2023.1202456.
  * Li et al. [2024] Z. Li, J. Gong, Y. Harada, B. Da, R. Zeng, D. ZJ,  Determination of the energy loss function of tungsten from reflection electron energy loss spectroscopy spectra,  Results Phys. 56 (2024) 107247. doi:https://doi.org/10.1016/j.rinp.2023.107247.
  * Taioli et al. [2009] S. Taioli, P. Umari, M. De Souza,  Electronic properties of extended graphene nanomaterials from gw calculations,  Phys. Status Solidi (B) 246 (2009) 2572–2576. doi:https://doi.org/10.1002/pssb.200982339.
  * Umari et al. [2012] P. Umari, O. Petrenko, S. Taioli, M. De Souza,  Communication: Electronic band gaps of semiconducting zig-zag carbon nanotubes from many-body perturbation theory calculations,  J. Chem. Phys. 136 (2012). doi:https://doi.org/10.1063/1.4716178.
  * Pedrielli et al. [2022] A. Pedrielli, N. M. Pugno, M. Dapor, S. Taioli,  In search of the ground state crystal structure of $\rm{Ta_{2}O_{5}}$ from ab-initio and $\rm{M}$onte $\rm{C}$arlo simulations,  Comput. Mater. Sci. 216 (2022) 111828. doi:https://doi.org/10.1016/j.commatsci.2022.111828.
  * Ding and Shimizu [1996] Z. Ding, R. Shimizu,  A $\rm{M}$onte $\rm{C}$arlo modeling of electron interaction with solids including cascade secondary electron production,  Scanning 18 (1996) 92–113. doi:https://doi.org/10.1002/sca.1996.4950180204.
  * Joy and Joy [1996] D. Joy, C. Joy,  Low voltage scanning electron microscopy,  Micron 27 (1996) 247–263. doi:https://doi.org/10.1016/0968-4328(96)00023-6.
  * Ding et al. [2001] Z. Ding, X. Tang, R. Shimizu,  $\rm{M}$onte $\rm{C}$arlo study of secondary electron emission,  J. Appl. Phys. 89 (2001) 718–726. doi:https://doi.org/10.1063/1.1331645.
  * Vos [2001] M. Vos,  Observing atom motion by electron-atom $\rm{C}$ompton scattering,  Phys. Rev. A 65 (2001) 012703. doi:https://doi.org/10.1103/PhysRevA.65.012703.
  * Orosz et al. [2004] G. Orosz, G. Gergely, M. Menyhard, J. Tóth, B. Varga D, Lesiak, J. A,  Hydrogen and surface excitation in electron spectra of polyethylene,  Surf. Sci. 566–568 (2004) 544–548. doi:https://doi.org/10.1016/j.susc.2004.05.106.
  * Yubero et al. [2008] F. Yubero, N. Pauly, A. Dubus, S. Tougaard,  Test of validity of the $\rm{V}$-type approach for electron trajectories in reflection electron energy loss spectroscopy,  Phys. Rev. B 77 (2008) 245405. doi:https://doi.org/10.1103/PhysRevB.77.245405.
  * Mao et al. [2008] S. Mao, Y. Li, R. Zeng, Z. Ding,  Electron inelastic scattering and secondary electron emission calculated without the single pole approximation,  J. Appl. Phys. 104 (2008) 114907. doi:https://doi.org/10.1063/1.3033564.
  * Larciprete et al. [2013] R. Larciprete, D. Grosso, M. Commisso, R. Flammini, R. Cimino,  Secondary electron yield of $\rm{C}$u technical surfaces: dependence on electron irradiation,  Phys. Rev. ST Accel. Beams. 16 (2013) 011002. doi:https://doi.org/10.1103/PhysRevSTAB.16.011002.
  * Salvat‐Pujol and Werner [2013] F. Salvat‐Pujol, W. Werner,  Surface excitations in electron spectroscopy. $\rm{P}$art $\rm{I}$: dielectric formalism and $\rm{M}$onte $\rm{C}$arlo algorithm,  Surf. Interface Anal. 45 (2013) 873–894. doi:https://doi.org/10.1002/sia.5175.
  * Bellissimo et al. [2020] A. Bellissimo, G. Pierantozzi, A. Ruocco, G. Stefani, O. Ridzel, A. Vytautas, W. Werner, M. Taborelli,  Secondary electron generation mechanisms in carbon allotropes at low impact electron energies,  J. Electron Spectrosc. Relat. Phenom. 241 (2020) 146883. doi:https://doi.org/10.1016/j.elspec.2019.07.004.
  * Tougaard et al. [2022] S. Tougaard, N. Pauly, F. Yubero,  $\rm{QUEELS}$: Software to calculate the energy loss processes in $\rm{TEELS}$, $\rm{REELS}$, $\rm{XPS}$ and $\rm{AES}$ including effects of the core hole,  Surf. Interface Anal. 54 (2022) 820–833. doi:https://doi.org/10.1002/sia.7095.
  * Polak and Morgan [2021] M. Polak, D. Morgan,  $\rm{MAST-SEY}$: Material simulation toolkit for secondary electron yield. $\rm{A}$ $\rm{M}$onte $\rm{C}$arlo approach to secondary electron emission based on complex dielectric functions,  Comput. Mater. Sci. 193 (2021) 110281. doi:https://doi.org/10.48550/arXiv.2108.11582.
  * Khan et al. [2023] M. Khan, S. Mao, Y. Zou, D. Lu, B. Da, Y. Li, Z. Ding,  An extensive theoretical quantification of secondary electron emission from silicon,  Vacuum 215 (2023) 112257. doi:10.1016/j.vacuum.2023.112257.
  * Pauly et al. [2023] N. Pauly, F. Yubero, S. Tougaard,  Determination of the primary excitation spectra in $\rm{XPS}$ and $\rm{AES}$,  Nanomaterials 13 (2023) 339. doi:https://doi.org/10.3390/nano13020339.
  * Dapor et al. [2023] M. Dapor, R. Garcia-Molina, P. E. Trevisanutto, Z. Ding, S. Taioli,  Editorial: Methodological and computational developments for modeling the transport of particles within materials,  Front. Mater. 10 (2023) 1302000. doi:https://doi.org/10.3389/fmats.2023.1302000.
  * Durante and Cucinotta [2011] M. Durante, F. Cucinotta,  Physical basis of radiation protection in space travel,  Rev. Mod. Phys. 83 (2011) 1245–1281. doi:https://doi.org/10.1103/RevModPhys.83.1245.
  * Taioli et al. [2023] S. Taioli, M. Dapor, F. Dimiccoli, M. Fabi, V. Ferroni, C. Grimani, M. Villani, W. J. Weber,  The role of low-energy electrons in the charging process of $\rm{LISA}$ test masses,  Class. Quantum Gravity 40 (2023) 075001. doi:https://doi.org/110.1088/1361-6382/acbadd.
  * Solov’yov et al. [2009] A. Solov’yov, E. Surdutovich, E. Scifoni, I. Mishustin, W. Greiner,  Physics of ion beam cancer therapy: A multiscale approach,  Phys. Rev. E 79 (2009) 011909. doi:https://doi.org/10.1103/PhysRevE.79.011909.
  * Surdutovich and Solov’yov [2014] E. Surdutovich, A. V. Solov’yov,  Multiscale approach to the physics of radiation damage with ions,  Eur. Phys. J. D 68 (2014) 353. doi::https://doi.org/110.1063/1.4802412.
  * Surdutovich and Solov’yov [2015] E. Surdutovich, A. Solov’yov,  Transport of secondary electrons and reactive species in ion tracks,  Eur. Phys. J. D 69 (2015) 193. doi:https://doi.org/10.1140/epjd/e2015-60242-1.
  * Nikjoo et al. [2016] H. Nikjoo, D. Emfietzoglou, T. Liamsuwan, D. Taleei, R andt, S. Uehara,  Radiation track, $\rm{DNA}$ damage and response: $\rm{A}$ review,  Rep. Prog. Phys. 79 (2016) 116601. doi:https://doi.org/10.1088/0034-4885/79/11/116601.
  * Solov’yov [2017] A. V. Solov’yov (Ed.), Nanoscale Insights into Ion-Beam Cancer Therapy, Springer International Publishing Switzerland, 2017. URL: https://link.springer.com/book/10.1007/978-3-319-43030-0.
  * Conte et al. [2017] V. Conte, A. Selva, P. Colautti, G. Hilgers, H. Rabus,  Track structure characterization and its link to radiobiology,  Radiat. Meas. 83 (2017) 506–511. doi:https://doi.org/10.1016/j.radmeas.2017.06.010.
  * Conte et al. [2018] V. Conte, A. Selva, P. Colautti, G. Hilgers, H. Rabus, A. Bantsar, M. Pietrzak, S. Pszona,  Nanodosimetry: Towards a new concept of radiation quality,  Radiat. Prot. Dosim. 180 (2018) 150–156. doi:https://doi.org/10.1093/rpd/ncx175.
  * Friedrich et al. [2018] H. Friedrich, K. Ilicic, C. Greubel, S. Girst, J. Reindl, M. Sammer, B. Schwarz, C. Siebenwirth, D. W. M. Walsh, T. E. Schmid, M. Scholz, G. Dollinger,  DNA damage interactions on both nanometer and micrometer scale determine overall cellular damage,  Sci. Rep. 8 (2018) 16063. doi:https://doi.org/10.1038/s41598-018-34323-9.
  * Surdutovich and Solov’yov [2019] E. Surdutovich, A. Solov’yov,  Multiscale modeling for cancer radiotherapies,  Cancer Nanotechnol. 10 (2019) 6. doi:https://doi.org/10.1186/s12645-019-0051-2.
  * Taioli et al. [2020] S. Taioli, P. E. Trevisanutto, P. de Vera, S. Simonucci, I. Abril, R. Garcia-Molina, M. Dapor,  Relative role of physical mechanisms on complex biodamage induced by carbon irradiation,  J. Phys. Chem. Lett. 12 (2020) 487–493. doi:https://doi.org/10.1021/acs.jpclett.0c03250.
  * de Vera et al. [2022a] P. de Vera, S. Taioli, P. E. Trevisanutto, M. Dapor, I. Abril, S. Simonucci, R. Garcia-Molina,  Energy deposition around swift carbon-ion tracks in liquid water,  Int. J. Mol. Sci. 23 (2022a) 6121. doi:https://doi.org/10.3390/ijms23116121.
  * de Vera et al. [2022b] P. de Vera, S. Simonucci, P. E. Trevisanutto, I. Abril, M. Dapor, S. Taioli, R. Garcia-Molina,  Simulating the nanometric track-structure of carbon ion beams in liquid water at energies relevant for hadrontherapy,  J. Phys. Conf. Ser. 2326 (2022b) 012017. doi:https://doi.org/10.1088/1742-6596/2326/1/012017.
  * Dapor et al. [2015] M. Dapor, I. Abril, P. de Vera, R. Garcia-Molina,  Simulation of the secondary electrons energy deposition produced by proton beams in pmma: influence of the target electronic excitation description,  Eur. Phys. J. D 69 (2015) 165. doi:https://doi.org/10.1140/epjd/e2015-60123-7.
  * Dapor et al. [2017] M. Dapor, I. Abril, P. de Vera, R. Garcia-Molina,  Energy deposition around swift proton tracks in polymethylmethacrylate: How much and how far,  Phys. Rev. B 96 (2017) 064113. doi:https://doi.org/10.1103/PhysRevB.96.064113.
  * Gorfinkiel et al. [2005] J. Gorfinkiel, A. Faure, S. Taioli, C. Piccarreta, G. Halmova, J. Tennyson,  Electron-molecule collisions at low and intermediate energies using the r-matrix method,  Eur. Phys. J. D 35 (2005) 231–237. doi:https://doi.org/10.1140/epjd/e2005-00179-4.
  * Taioli and Tennyson [2006] S. Taioli, J. Tennyson,  Waterwaves: wave particles dynamics on a complex triatomic potential,  Comput. Phys. Commun. 175 (2006) 41–51. doi:https://doi.org/10.1016/j.cpc.2006.01.008.
  * Taioli et al. [2010] S. Taioli, S. Simonucci, L. Calliari, M. Dapor,  Electron spectroscopies and inelastic processes in nanoclusters and solids: Theory and experiment,  Phys. Rep. 493 (2010) 237–319. doi:https://doi.org/10.1016/j.physrep.2010.04.003.
  * Taioli and Simonucci [2015] S. Taioli, S. Simonucci,  A computational perspective on multichannel scattering theory with applications to physical and nuclear chemistry,  in: Annual Reports in Computational Chemistry, volume 11, Elsevier, 2015, pp. 191–310. doi:https://doi.org/10.1016/bs.arcc.2015.09.005.
  * Shimizu and Ding [1992] R. Shimizu, Z.-J. Ding,  Monte $\rm{C}$arlo modelling of electron-solid interactions,  Rep. Prog. Phys. 55 (1992) 487–531. doi:https://doi.org/10.1088/0034-4885/55/4/002.
  * Joy [1995] D. C. Joy, Monte Carlo modeling for electron microscopy and microanalysis, volume 9, Oxford University Press, 1995. URL: https://academic.oup.com/book/54583.
  * Dapor [2003] M. Dapor, Electron-Beam Interactions with Solids. Application of the Monte Carlo Method to Electron Scattering Problems, volume 186 of Springer Tracts in Modern Physics, Springer-Verlag Berlin Heidelberg, 2003. doi:https://doi.org/10.1007/3-540-36507-9.
  * Dapor [2023] M. Dapor, Transport of Energetic Electrons in Solids. Computer Simulation with Applications to Materials Analysis and Characterization, volume 290 of Springer Tracts in Modern Physics, 4th ed., Springer Nature Switzerland AG, 2023. doi:https://doi.org/10.1007/978-3-031-37242-1.
  * Liljequist [2008] D. Liljequist,  A study of errors in trajectory simulation with relevance for 0.2–50 e$\rm{V}$ electrons in liquid water,  Radiat. Phys. Chem. 77 (2008) 835–853. doi:10.1016/j.radphyschem.2008.03.004.
  * Liljequist [2013] D. Liljequist,  Contribution from inelastic scattering to the validity of trajectory methods,  J. Electron Spectrosc. Relat. Phenom. 189 (2013) 5–11. doi:https://doi.org/10.1016/j.elspec.2013.04.011.
  * Ritchie [1957] R. H. Ritchie,  Plasma losses by fast electrons in thin films,  Phys. Rev. 106 (1957) 874–881. doi:https://doi.org/10.1103/PhysRev.106.874.
  * Pedrielli et al. [2021] A. Pedrielli, P. de Vera, P. E. Trevisanutto, N. M. Pugno, R. Garcia-Molina, I. Abril, S. Taioli, M. Dapor,  Electronic excitation spectra of cerium oxides: from ab initio dielectric response functions to $\rm{M}$onte $\rm{C}$arlo electron transport simulations,  Phys. Chem. Chem. Phys. 23 (2021) 19173–19187. doi:https://doi.org/10.1039/D1CP01810H.
  * Azzolini et al. [2017] M. Azzolini, T. Morresi, G. Garberoglio, L. Calliari, N. M. Pugno, S. Taioli, M. Dapor,  Monte $\rm{C}$arlo simulations of measured electron energy-loss spectra of diamond and graphite: role of dielectric-response models,  Carbon 118 (2017) 299–309. doi:http://dx.doi.org/10.1016/j.carbon.2017.03.041.
  * Azzolini et al. [2018] M. Azzolini, T. Morresi, K. Abrams, R. Masters, N. Stehling, C. Rodenburg, N. M. Pugno, S. Taioli, M. Dapor,  Anisotropic approach for simulating electron transport in layered materials: computational and experimental study of highly oriented pyrolitic graphite,  J. Phys. Chem. C 122 (2018) 10159–10166. doi:https://doi.org/10.1021/acs.jpcc.8b02256.
  * Triggiani et al. [2023] F. Triggiani, T. Morresi, S. Taioli, S. Simonucci,  Elastic scattering of electrons by water: an ab initio study,  Front. Mater. 10 (2023) 1145261. doi:https://doi.org/10.3389/fmats.2023.1145261.
  * Newton [1982] R. Newton, Scattering Theory of Waves and Particles, Springer Berlin, Heidelberg, 1982\. doi:https://doi.org/10.1007/978-3-642-88128-2.
  * Tennyson [2010] J. Tennyson,  Electron–molecule collision calculations using the R-matrix method,  Phys. Rep. 491 (2010) 29–76. doi:https://doi.org/10.1016/j.physrep.2010.02.001.
  * Mott [1929] N. F. Mott,  The scattering of fast electrons by atomic nuclei,  Proc. R. soc. Lond. Ser. A-Contain. Pap. Math. Phys. Character 124 (1929) 425–442. doi:https://doi.org/10.1098/rspa.1929.0127.
  * Salvat et al. [1987] F. Salvat, J. D. Martinez, R. Mayol, J. Parellada,  $\rm{A}$nalytical $\rm{D}$irac-$\rm{H}$artree-$\rm{F}$ock-$\rm{S}$later screening function for atoms ($\rm{Z}$=1–92),  Phys. Rev. A 36 (1987) 467–474. doi:https://doi.org/10.1103/PhysRevA.36.467.
  * Onida et al. [2002] G. Onida, L. Reining, A. Rubio,  Electronic excitations: density-functional versus many-body green’s-function approaches,  Rev. Mod. Phys. 74 (2002) 601–659. doi:https://doi.org/10.1103/RevModPhys.74.601.
  * Ashley [1988] J. C. Ashley,  Interaction of low-energy electrons with condensed matter: stopping powers and inelastic mean free paths from optical data,  J. Electron Spectrosc. Relat. Phenom. 46 (1988) 199–214. doi:https://doi.org/10.1016/0368-2048(88)80019-7.
  * Ashley [1990] J. C. Ashley,  Energy loss rate and inelastic mean free path of low-energy electrons and positrons in condensed matter,  J. Electron Spectrosc. Relat. Phenom. 50 (1990) 323–334. doi:https://doi.org/10.1016/0368-2048(90)87075-Y.
  * Tanuma et al. [1993] S. Tanuma, C. J. Powell, D. R. Penn,  Calculation of electron inelastic mean free paths,  Surf. Interface Anal. 21 (1993) 165–176. doi:https://doi.org/110.1002/sia.740210302.
  * Garcia-Molina et al. [2012] R. Garcia-Molina, I. Abril, I. Kyriakou, D. Emfietzoglou,  Energy loss of swift protons in liquid water: Role of optical data input and extension algorithms,  in: Radiation Damage in Biomolecular Systems, 2012, pp. 239–261. doi:https://doi.org/10.1007/978-94-007-2564-515.
  * Abril et al. [1998] I. Abril, R. Garcia-Molina, C. D. Denton, F. J. Pérez-Pérez, N. R. Arista,  Dielectric description of wakes and stopping powers in solids,  Phys. Rev. A 58 (1998) 357–366. doi:https://doi.org/10.1103/PhysRevA.58.357.
  * Egerton [2011] R. F. Egerton, Electron Energy-Loss Spectroscopy in the Electron Microscope, Springer, New York, Dordrecht, Heidelberg, London, 2011. doi:https://doi.org/10.1007/978-1-4419-9583-4.
  * Fröhlich [1954] H. Fröhlich,  Electrons in lattice fields,  Adv. Phys. 3 (1954) 325–361. doi:https://doi.org/10.1080/00018735400101213.
  * Zhou et al. [2021] J.-J. Zhou, J. Park, I.-T. Lu, I. Maliyov, X. Tong, M. Bernardi,  Perturbo: A software package for ab initio electron–phonon interactions, charge transport and ultrafast dynamics,  Comput. Phys. Commun. 264 (2021) 107970. doi:https://doi.org/10.1016/j.cpc.2021.107970.
  * Ganachaud and Mokrani [1995] J. Ganachaud, A. Mokrani,  Theoretical study of the secondary electron emission of insulating targets,  Surface Science 334 (1995) 329–341. doi:https://doi.org/10.1016/0039-6028(95)00474-2.
  * Franchini et al. [2021] C. Franchini, M. Reticcioli, M. Setvin, U. Diebold,  Polarons in materials,  Nat. Rev. Mater. 6 (2021) 560–586. doi:10.1038/s41578-021-00289-w.
  * Azzolini et al. [2020] M. Azzolini, O. Y. Ridzel, P. S. Kaplya, V. Afanas’ev, N. M. Pugno, S. Taioli, M. Dapor,  A comparison between $\rm{M}$onte $\rm{C}$arlo method and the numerical solution of the $\rm{A}$mbartsumian-$\rm{C}$handrasekhar equations to unravel the dielectric response of metals,  Comput. Mater. Sci. 173 (2020) 109420. doi:https://doi.org/10.1016/j.commatsci.2019.109420.
  * Ibach [1977] H. Ibach, Electron Spectroscopy for Surface Analysis, Springer, Berlin, Heidelberg, 1977\. URL: https://link.springer.com/book/10.1007/978-3-642-81099-2. doi:https://doi.org/10.1007/978-3-642-81099-2.
  * Gergely [2002] G. Gergely,  Elastic backscattering of electrons: determination of physical parameters of electron transport processes by elastic peak electron spectroscopy,  Progr. Surf. Sci. 71 (2002) 31–88. doi:https://doi.org/10.1016/S0079-6816(02)00019-9.
  * Jablonski [2003] A. Jablonski,  Analytical applications of elastic electron backscattering from surfaces,  Progr. Surf. Sci. 74 (2003) 357–374. doi:https://doi.org/10.1016/j.progsurf.2003.08.028.
  * Taioli et al. [2009] S. Taioli, S. Simonucci, M. Dapor,  Surprises: when ab initio meets statistics in extended systems,  Comput. Sci. Discov. 2 (2009) 015002. doi:https://doi.org/10.1088/1749-4699/2/1/015002.
  * Taioli and Simonucci [2021] S. Taioli, S. Simonucci,  The resonant and normal auger spectra of ozone,  Symmetry 13 (2021) 516. doi:https://doi.org/10.3390/sym13030516.
  * Colle et al. [2004a] R. Colle, D. Embriaco, M. Massini, S. Simonucci, S. Taioli,  Ab initio calculation of the normal $\rm{A}$uger spectrum of $\rm{C}_{2}\rm{H}_{2}$,  J. Phys. B 37 (2004a) 1237. doi:https://doi.org/10.1088/0953-4075/37/6/008.
  * Colle et al. [2004b] R. Colle, D. Embriaco, M. Massini, S. Simonucci, S. Taioli,  Auger-electron angular distributions calculated without the two-step approximation: Calculation of angle-resolved resonant $\rm{A}$uger spectra of $\rm{C}_{2}\rm{H}_{2}$,  Phys. Rev. A 70 (2004b) 042708. doi:https://doi.org/10.1103/PhysRevA.70.042708.
  * Colle et al. [2004c] R. Colle, D. Embriaco, M. Massini, S. Simonucci, S. Taioli,  Ab initio calculation of the c1s photoelectron spectrum of $\rm{C}_{2}\rm{H}_{2}$,  Nucl. Instrum. Methods Phys. Res. B 213 (2004c) 65–70. doi:https://doi.org/10.1016/S0168-583X(03)01535-0.
  * Dapor et al. [2010] M. Dapor, M. Ciappa, W. Fichtner,  Monte carlo modeling in the low-energy domain of the secondary electron emission of polymethylmethacrylate for critical-dimension scanning electron microscopy,  Journal of Micro/Nanolithography, MEMS, and MOEMS 9 (2010) 023001. doi:https://doi.org/10.1117/1.3373517.
  * Garcia-Molina et al. [2006] R. Garcia-Molina, I. Abril, C. D. Denton, S. Heredia-Avalos,  Allotropic effects on the energy loss of swift $\rm{H}$\+ and $\rm{H}$e+ ion beams through thin foils,  Nucl. Instrum. Methods Phys. Res. B 249 (2006) 6–12. doi:https://doi.org/10.1016/j.nimb.2006.03.011.
  * Dapor et al. [2011] M. Dapor, N. Bazzanella, L. Toniutti, A. Miotello, S. Gialanella,  Backscattered electrons from surface films deposited on bulk targets: A comparison between computational and experimental results,  Nucl. Instrum. Methods Phys. Res. B 269 (2011) 1672–1674. doi:https://doi.org/10.1016/j.nimb.2010.11.016.
  * Everhart [1960] T. E. Everhart,  Simple Theory Concerning the Reflection of Electrons from Solids,  J. Appl. Phys. 31 (1960) 1483–1490. doi:https://doi.org/10.1063/1.1735868.
  * Thomson [1906] J. J. Thomson, Conduction of Electricity Through Gases, 2nd ed., Cambridge University Press, Cambridge, England, 1906\.
  * Whiddington [1914] R. Whiddington,  The transmission of cathode rays through matter,  Proc. R. Soc. Lond. 89 (1914) 554–560. doi:https://doi.org/10.1098/rspa.1912.0028.
  * Niedrig [1982] H. Niedrig,  Analytical models in electron backscattering,  Scanning Electron Microscop y 1982 (1982) 51–68. URL: https://digitalcommons.usu.edu/electron/vol1982/iss1/5.
  * Archard [2004] G. D. Archard,  Back Scattering of Electrons,  J. Appl. Phys. 32 (2004) 1505–1509. doi:https://doi.org/10.1063/1.1728385.
  * Neubert and Rogaschewski [1980] G. Neubert, S. Rogaschewski,  Backscattering coefficient measurements of 15 to 60 kev electrons for solids at various angles of incidence,  Phys. Status Solidi A 59 (1980) 35–41. doi:https://doi.org/10.1515/9783112501726-005.
  * Vicanek and Urbassek [1991] M. Vicanek, H. M. Urbassek,  Reflection coefficient of low-energy light ions,  Phys. Rev. B 44 (1991) 7234–7242. doi:https://doi.org/10.1103/PhysRevB.44.7234.
  * Coleman et al. [1992] P. G. Coleman, L. Albrecht, K. O. Jensen, A. B. Walker,  Positron backscattering from elemental solids,  J. Phys. Condens. Matter 4 (1992) 10311–10322. doi:https://doi.org/10.1088/0953-8984/4/50/018.
  * Dapor [1996] M. Dapor,  Elastic scattering calculations for electrons and positrons in solid targets,  J. Appl. Phys. 79 (1996) 8406–8411. doi:https://doi.org/10.1063/1.362514.
  * Wolff [1954a] P. A. Wolff,  Theory of secondary electron cascade in metals,  Phys. Rev. 95 (1954a) 56–66. doi:https://doi.org/10.1103/PhysRev.95.56.
  * Wolff [1954b] P. A. Wolff,  Theory of electron multiplication in silicon and germanium,  Phys. Rev. 95 (1954b) 1415–1420. doi:https://doi.org/10.1103/PhysRev.95.1415.
  * Amelio [1970] G. F. Amelio,  Theory for the Energy Distribution of Secondary Electrons,  J. Vac. Sci. Technol. 7 (1970) 593–604. doi:https://doi.org/10.1116/1.1315884.
  * Streitwolf [1959] H. Streitwolf,  On the theory of secondary electron emission from metals the excitation process,  Ann. Phys. 458 (1959) 183–196. doi:https://doi.org/10.1002/andp.19594580308.
  * Chung and Everhart [1974] M. S. Chung, T. E. Everhart,  Simple calculation of energy distribution of low‐energy secondary electrons emitted from metals under electron bombardment,  J. Appl. Phys. 45 (1974) 707–709. doi:https://doi.org/10.1063/1.1663306.
  * Ceperley and Alder [1986] D. Ceperley, B. Alder,  Quantum $\rm{M}$onte $\rm{C}$arlo,  Science 231 (1986) 555–560. doi:https://doi.org/10.1126/science.231.4738.555.
  * Becca and Sorella [2017] F. Becca, S. Sorella, Quantum Monte Carlo approaches for correlated systems, Cambridge University Press, 2017. doi:https://doi.org/10.1017/9781316417041.
  * Yubero and Tougaard [2 II] F. Yubero, S. Tougaard,  Model for quantitative analysis of reflection-electron-energy-loss spectra,  Phys. Rev. B 46 (1992-II) 2489–2497. doi:https://doi.org/10.1103/PhysRevB.46.2486.
  * Giannozzi et al. [2009] P. Giannozzi, et al.,  $\rm{QUANTUM}$ $\rm{ESPRESSO}$: a modular and open-source software project for quantum simulations of materials,  J. Phys. Condens. Matter 21 (2009) 395502. doi:https://doi.org/10.1088/0953-8984/21/39/395502.
  * The Elk Code [2024] The Elk Code, 2024. http://elk.sourceforge.net/.
  * Giannozzi and Baroni [2005] P. Giannozzi, S. Baroni, Density-Functional Perturbation Theory, Springer Science+Business Media B.V., 2005, p. 195, 2005, pp. 195–214. doi:https://doi.org/110.1007/978-1-4020-3286-811.