#  QED theory of electron beam-induced electronic excitation and its effect on
sputtering cross sections in 2D crystals

Anthony Yoshimura  Lawrence Livermore National Laboratory, Livermore, CA
94550, USA  Michael Lamparski  Department of Physics, Applied Physics, and
Astronomy, Rensselaer Polytechnic Institute, Troy, NY 12180, USA  Joel Giedt
Department of Physics, Applied Physics, and Astronomy, Rensselaer Polytechnic
Institute, Troy, NY 12180, USA
David Lingerfelt  Center for Nanophase Material Sciences, Oak Ridge National
Laboratory, Oak Ridge, TN 37831, USA  Jacek Jakowski  Center for Nanophase
Material Sciences, Oak Ridge National Laboratory, Oak Ridge, TN 37831, USA
Panchapakesan Ganesh  Center for Nanophase Material Sciences, Oak Ridge
National Laboratory, Oak Ridge, TN 37831, USA  Tao Yu  Department of
Chemistry, University of North Dakota, Grand Forks, ND 58202, USA
Bobby Sumpter  Center for Nanophase Material Sciences, Oak Ridge National
Laboratory, Oak Ridge, TN 37831, USA  Vincent Meunier  Department of Physics,
Applied Physics, and Astronomy, Rensselaer Polytechnic Institute, Troy, NY
12180, USA  Department of Materials Science and Engineering, Rensselaer
Polytechnic Institute, Troy, NY 12180, USA

###### Abstract

Many computational models have been developed to predict the rates of atomic
displacements in two-dimensional (2D) materials under electron beam
irradiation. However, these models often drastically underestimate the
displacement rates in 2D insulators, in which beam-induced electronic
excitations can reduce the binding energies of the irradiated atoms. This bond
softening leads to a qualitative disagreement between theory and experiment,
in that substantial sputtering is experimentally observed at beam energies
deemed far to small to drive atomic dislocation by many current models. To
address these theoretical shortcomings, this paper develops a first-principles
method to calculate the probability of beam-induced electronic excitations by
coupling quantum electrodynamics (QED) scattering amplitudes to density
functional theory (DFT) single-particle orbitals. The presented theory then
explicitly considers the effect of these electronic excitations on the
sputtering cross section. Applying this method to 2D hexagonal BN and MoS2
significantly increases their calculated sputtering cross sections and
correctly yields appreciable sputtering rates at beam energies previously
predicted to leave the crystals intact. The proposed QED-DFT approach can be
easily extended to describe a rich variety of beam-driven phenomena in any
crystalline material.

##  1 Introduction

The holy grail of materials engineering is atomic scale control of the
material structure. Towards this aim, electron irradiation by transmission
electron microscopy (TEM) can be an effective means of structural manipulation
with spatial control [1, 2, 3, 4]. Structural changes under electron
irradiation can arise via atomic displacement in which an incident electron
nudges a material atom from its initial site. We call this atom the primary
knock-on atom (PKA). Two-dimensional (2D) crystals provide an excellent
platform to measure the rates of these PKA displacements. When electron
irradiation is normal to a 2D crystal’s surface, a displacement likely propels
the PKA away from the crystal. As such, these displaced atoms are often
ejected from the crystal in a process called sputtering. Sputtering events
leave behind vacancies, which can then be counted using TEM. Counting the
number of vacancies for a given dosage and beam energy allows one to
experimentally determine the sputtering cross section of that crystal.

Sputtering occurs when the energy transferred to the PKA is greater than the
PKA’s displacement threshold $E_{d}$. This means that a displacement is
possible only if the kinetic energy of the beam electron exceeds some critical
energy $\epsilon_{c}$. Many computational models have been proposed to predict
both $E_{d}$ and $\epsilon_{c}$ to calculate electron beam-induced sputtering
rates in 2D crystals [5, 6, 7, 4]. However, the vast majority of current
methods focus solely on interactions between the beam electrons and material
nuclei, neglecting any coupling with the material’s electrons. Thus, while
present-day models give reasonable predictions for conductors [5], where
electronic relaxation is rapid, they often vastly underestimate the atomic
displacement rates in insulators. For example, the critical energy for
sputtering boron or nitrogen from hexagonal boron nitride (hBN) is predicted
to be 80 keV [8]. However, sputtering has been observed in hBN under 30 keV
irradiation [9]. Furthermore, selenium sputters from WSe2 and MoSe2 under
irradiation energies of 60 and 80 keV, respectively. These energies are almost
150 keV below their predicted critical energies [10, 11]. Lastly, while the
calculated critical energy for sulfur sputtering in MoS2 is about 90 keV,
sulfur has been shown to sputter under 20 keV beams [12]. Discrepancies like
these suggest that the displacement thresholds in insulating crystals are much
smaller than what is predicted by ground-state theory. Lehnert et al. have
proposed that the consideration of inelastic scattering, i.e., beam-induced
electronic excitation, can lead to such a reduction in the displacement
threshold [11]. This would increase the sputtering cross section for all beam
energies and enable sputtering for energies well-below the ground state
$\epsilon_{c}$.

To account for these effects, we combine quantum electrodynamics (QED) and
density functional theory (DFT) to derive the probability of beam-induced
electronic excitation in 2D insulating crystals. The basic idea is as follows:
DFT can provide effective single-particle states that can be decomposed into a
plane-wave basis [13, 14, 15], while QED is well-equipped to describe how each
plane-wave evolves in time through interactions with an electromagnetic field
[16, 17, 18]. Thus, a plane-wave decomposition of the Kohn-Sham orbitals can
allow for a component-by-component treatment of the interactions between the
beam and material electrons. This generalized QED-DFT approach enables, for
the first time, a first-principles description of any beam-matter interaction
process. The only limitations of this method are the order to which the time-
evolution operator is expanded and the sophistication of the theory used to
determine the material’s electronic structure. Additionally, while DFT is used
here, our method is compatible with any first-principles formalism that can
produce single-particle eigenstates and eigenvalues for a given material.

This paper is divided into three sections. First, we describe the key physical
processes that dictate the rate of beam-induced sputtering in a 2D insulating
crystal. Second, we derive the probability of beam-induced electronic
excitation in these materials as a function of beam energy. Finally, we show
how this excitation probability can be used to predict sputtering cross
sections in hBN and MoS2 that quantitatively agree with experiment.

##  2 Three interactions and three rates

The majority of present-day beam-damage models focus solely on the interaction
between the beam electron and target nucleus. These models are thus centered
on one process: energy transfer from the beam electron to the nucleus. For
this, one defines a differential cross section $d\sigma/dE(\epsilon_{b},E)$
providing the distribution of energy transfers $E$ for a given beam energy
$\epsilon_{b}$. The McKinley and Feshbach differential cross section can
adequately describe $d\sigma/dE$ for light nuclei (Z $<$ 20) [19, 20, 21, 22].
Setting $\hbar=c=1$,

$\frac{d\sigma}{dE}=\pi\left(\frac{Z\alpha}{|\mathbf{p}_{b}|\beta}\right)^{2}\left[\frac{E_{\text{max}}}{E^{2}}-\beta(\pi
Z\alpha+\beta)\frac{1}{E}+\pi
Z\alpha\beta\sqrt{\frac{E_{\text{max}}}{E^{3}}}\right],$ (1)

where $\mathbf{p}_{b}$ is the momentum of the beam electron, $\beta$ is its
velocity, $Z$ is the atomic number of the target nucleus, and $\alpha$ is the
fine structure constant.

One can then calculate the displacement cross section by integrating
$d\sigma/dE$ over all $E$ large enough to cause a displacement, so that

$\sigma_{0}(\epsilon_{b},{E}_{d})=\theta(E_{\text{max}}(\epsilon_{b})-E_{d})\int_{E_{d}}^{E_{\text{max}}(\epsilon_{b})}\frac{d\sigma}{dE}(\epsilon_{b},E)dE,$
(2)

where $E_{\text{max}}(\epsilon_{b})$ is the maximum possible energy transfer
for a given beam energy, i.e, the energy transfer resulting from a direct
collision, and the step function $\theta$ enforces that the cross section is
zero when $E_{\text{max}}<E_{d}$. We will leave the step function implicit
going forward for the sake of compactness. The critical energy $\epsilon_{c}$
is then the beam energy for which $E_{\text{max}}=E_{d}$. Therefore, the
observation of sputtering at energies well below $\epsilon_{c}$ is completely
at odds with equation (2). With this in mind, a fair amount of work has been
done to treat deviations from equation (2). Notably, several studies have
explored the effects of temperature on displacement cross sections [5, 6, 7].
This consideration involves calculating the degree to which the pre-collision
thermal motion of the nucleus increases the cross section. However, these
techniques essentially amount to smearing the beam energy dependence of the
cross section, so that the cross section only strays significantly from
equation (2) for beam energies very close to $\epsilon_{c}$. Thus,
temperature-induced increases in the cross section cannot account for the
disparities between the equation (2) model and experiment. This necessitates
the consideration of additional phenomena that can reduce $E_{d}$.

To address the limitations of equation (2), this work introduces a third
party: the material’s electrons. Doing so brings two new interactions into
play: one between the beam and material electrons and another between material
electrons and nuclei. This yields a total of three interactions between the
three pairs of particles (figure 1). Therefore, the rate of beam-induced
sputtering hinges on the rates of three processes mediated by these
interactions.

Figure 1:  Three interactions between three pairs of particles involved in
electron beam-induced sputtering. In insulators, the probability of material
electronic excitation ($P_{i}$) and the probability that those excitations are
substantially long-lived ($P_{f}$) can be large. Therefore, the effect of
exciting the material electrons on the sputtering rate should not be ignored
in these materials.

  1. 1.

Beam and material electrons: a beam electron can excite some number $n_{i}$
ground state electrons to the conduction band ($i$ denotes the initial
interaction with the beam). The probability of this event for a given beam
energy $\epsilon_{b}$ is $P_{i}(\epsilon_{b},n_{i})$.

  2. 2.

Material electrons and nuclei: Some number $n_{f}$ beam-induced excitations
can survive long enough for the target atom to leave its original site ($f$
denotes the final system at the completion of sputtering). This depends on the
nuclear kinetic energy $E$ and the excitation lifetime $\tau$. The probability
that $n_{f}$ of the $n_{i}$ excitations survive is
$P_{f}(E,\tau,n_{i},n_{f})$.

  3. 3.

Beam electrons and material nuclei: the energy transferred to a material
nucleus by the beam electron can exceed the PKA’s displacement threshold
$E_{d}(n_{f})$, which depends on the number of surviving excitations $n_{f}$.
We define $\\{{E}_{d}\\}$ as the set of all displacement thresholds for all
possible $n_{f}$. Sputtering occurs when $E>E_{d}(n_{f})$. The differential
cross section for an energy transfer $E$ from the beam electron to material
nucleus is $d\sigma/dE(\epsilon_{b},E)$.

The sputtering cross section can then be calculated by coupling $d\sigma/dE$
to $P_{i}$ and $P_{f}$ for all possible $n_{i}$ and $n_{f}$. With the terms
defined above, this excitation-sensitive sputtering cross section can be
written as

$\displaystyle\sigma(\epsilon_{b},\\{{E}_{d}\\},\tau)$
$\displaystyle=\sum_{n_{i}=0}^{\infty}P_{i}(\epsilon_{b},n_{i})\sum_{n_{f}=0}^{n_{i}}\int_{E_{d}(n_{f})}^{E_{\text{max}}(\epsilon_{b})}P_{f}(E,\tau,n_{i},n_{f})\frac{d\sigma}{dE}(\epsilon_{b},E)dE.$
(3)

If $P_{i}$ and $P_{f}$ are non-negligible when $n_{i}$ and $n_{f}$ are
nonzero, and $E_{d}$ depends strongly on $n_{f}$, then interactions with the
material electrons must be considered. We will later show that this makes
$\sigma$ in equation (3) larger than $\sigma_{0}$ in equation (2) for all beam
energies, most prominently when $\epsilon_{b}<\epsilon_{c}$.

The remainder of this paper focuses on the derivations of $P_{i}$, $P_{f}$,
and $\\{E_{d}\\}$. Section 3 describes how to combine QED with DFT to obtain
$P_{i}$. Section 4 then considers the evolution of the excited states during
the sputtering process to derive $\\{E_{d}\\}$ and $P_{f}$. It then
demonstrates how our formalism significantly improves the prediction of
sputtering rates in hBN and MoS2.

##  3 Probability of beam-induced excitation

For a crystal in its ground state, an occupied electron energy eigenstate has
zero overlap with any unoccupied state. However, the collision of a beam
electron can give an occupied state a momentum boost that breaks this
orthogonality. Thus, the boosted ground state has a nonzero probability of
being measured in an excited state. We can use this idea to derive
$P_{i}(\epsilon_{b},n_{i})$, the probability that a beam electron with kinetic
energy $\epsilon_{b}$ excites exactly $n_{i}$ material electrons. The
derivation can be broken down into four steps: (i) determine the amplitude for
a free electron to scatter from one momentum eigenstate into another after
collision with another free electron; (ii) generalize the formalism to obtain
the amplitude for scattering from one wave packet into another by summing over
the amplitudes for each momentum component of one wave packet to scatter into
each momentum component of the other; (iii) decompose a pair of occupied and
unoccupied crystal states into a momentum basis and plug them in as incoming
and outgoing wave packets respectively, then square the amplitude to obtain
the corresponding excitation probability for a particular transition; (iv)
Compute the sum of all transition probabilities and use combinatorics to
determine $P_{i}(\epsilon_{b},n_{i})$. The following subsections address each
step (i-iv) in detail.

###  3.1 Scattering of free electrons

Figure 2:  The lowest order electron-electron scattering perturbation includes
two Feynman diagrams called the (a) $t$-channel and (b) $u$-channel. The
incoming and outgoing electron states are represented by Dirac spinors
$u^{s}(p)$ and $\bar{u}^{s}(p)$ respectively, where $p$ and $s$ are the
electron’s 4-momentum and spin index respectively. Subscripts 1 and 2 label
components of the initial beam and material states respectively. The virtual
photon 4-momentum $q$ is the momentum transfer between the electrons.

We first derive the scattering amplitude for momentum transfer between two
free electrons via Møller scattering [23, 24, 25]. With this, the interaction
between actual material states can be written as a linear combination of these
free-particle interactions. Going forward, we label the 4-momenta of the
incoming electrons as $p_{1}$ and $p_{2}$, while the outgoing electrons have
momenta $p_{3}$ and $p_{4}$. We also choose to make $p_{1}$ and $p_{2}$
components of the initial beam and material states respectively. The
4-momentum of the $n$th electron can be written as
$p_{n}=(\epsilon_{n},p_{n}^{x},p_{n}^{y},p_{n}^{z})=(\epsilon_{n},\mathbf{p}_{n})$,
where $\epsilon_{n}$ is the particle’s energy and $\mathbf{p}_{n}$ is its
3-momentum. Dot products between 4-vectors are then taken over Minkowski
space, so that $p_{n}\cdot
p_{m}=g_{\mu\nu}p_{n}^{\mu}p_{m}^{\nu}=\epsilon_{n}\epsilon_{m}-\mathbf{p}_{n}\cdot\mathbf{p}_{m}$,
where
$\mathbf{p}_{n}\cdot\mathbf{p}_{m}=p_{n}^{x}p_{m}^{x}+p_{n}^{y}p_{m}^{y}+p_{n}^{z}p_{n}^{z}$.

To lowest order, the amplitude for free electron scattering can be represented
by two tree-level diagrams, which we call the $t$\- and $u$-channels (figure
2). Using Feynman’s rules [17, 18], we can write these diagrams in terms of
Dirac spinors, yielding the invariant matrix element

$\displaystyle\mathcal{M}(p_{4}p_{3}\leftarrow p_{2}p_{1})$
$\displaystyle=\frac{e^{2}}{2}\sum_{s_{1}}\sum_{s_{2}}\sum_{s_{3}}\sum_{s_{4}}$
(4)
$\displaystyle\bigg{[}\bar{u}^{s_{4}}\left(p_{4}\right)\gamma^{\mu}u^{s_{1}}\left(p_{1}\right)\left(\frac{1}{p_{3}-p_{2}}\right)^{2}\bar{u}^{s_{3}}\left(p_{3}\right)\gamma_{\mu}u^{s_{2}}\left(p_{2}\right)$
$\displaystyle+\bar{u}^{s_{3}}\left(p_{3}\right)\gamma^{\mu}u^{s_{1}}\left(p_{1}\right)\left(\frac{1}{p_{4}-p_{2}}\right)^{2}\bar{u}^{s_{4}}\left(p_{4}\right)\gamma_{\mu}u^{s_{2}}\left(p_{2}\right)\bigg{]},$

where $s_{n}=1\text{ or }2$ denotes the spin of the $n$th electron, $u^{s}(p)$
is a Dirac spinor, and $\bar{u}^{s}(p)$ is its conjugate (section
LABEL:app:M). The factor of 1/2 before the summation arises from the
assumption that the incoming states are spin unpolarized.

The first term in brackets is the $t$-channel describing momentum transfer
$p_{3}-p_{2}$ and the second is the $u$-channel describing momentum transfer
$p_{4}-p_{2}$. Because the DFT cutoff energy is much smaller than the beam
energy, it is always the case that $|\mathbf{p}_{2}|$ is much smaller than
$|\mathbf{p}_{1}|$. Furthermore, we need only consider outgoing momenta for
which the kinetic energy associated with either $|\mathbf{p}_{3}|$ or
$|\mathbf{p}_{4}|$ falls within the DFT cutoff energy. In these cases, the
magnitude of one outgoing momentum is similar to $|\mathbf{p}_{2}|$, while
that of the other is much greater. This means that one channel’s momentum
transfer is always much larger than the other’s. As the momentum transfers
reside in the denominators of either channel in equation (4), it follows that
one channel always contributes much more to $\mathcal{M}$ than the other.
Thus, when the $t$-channel is significant, the $u$-channel is negligible, and
vice versa. Additionally, when integrating over all possible outgoing momenta,
the contribution of the $t$-channel is equal to that of the $u$-channel.
Taking advantage of this along with the indistinguishably of the electrons, we
calculate only the $t$-channel and multiply the resulting amplitude by 2
instead of calculating both channels and adding them. We can then define the
4-momentum transfer between the electrons as that of the $t$-channel: $q\equiv
p_{3}-p_{2}$. Because the $t$-channel has a $q^{2}$ in the denominator, the
resulting scattering probability is proportional to $q^{-4}$. This makes large
momentum transfers statistically irrelevant, allowing us to only consider
momentum transfers inside the first Brillouin zone (BZ). The evaluation of the
$t$-channel in terms of the components of the electrons’ 4-momenta is
straightforward, though cumbersome, and is described in section LABEL:app:M.

We can use the resulting $\mathcal{M}$ to obtain the free electron scattering
amplitude

$\bra{p_{4}p_{3}}\hat{T}\ket{p_{2}p_{1}}=(2\pi)^{4}\delta(p_{1}+p_{2}-p_{3}-p_{4})\mathcal{M}(p_{4}p_{3}\leftarrow
p_{2}p_{1}),$ (5)

where $\hat{T}$ is the scattering operator [17, 18]. This gives the amplitude
for two free electrons with momenta $p_{1}$ and $p_{2}$ to scatter into
$p_{3}$ and $p_{4}$. Equation (5) is used to derive the scattering amplitude
between two arbitrary wave packets in the next subsection.

###  3.2 Scattering of wave packets

The free particle scattering amplitude in equation (5) can be used to
determine the amplitude for the scattering of two arbitrary electron states
$\phi_{1}$ and $\phi_{2}$ into $\phi_{3}$ and $\phi_{4}$. This is obtained by
sandwiching the scattering operator between the initial and final 2-particle
states, i.e.,

$\displaystyle\bra{\phi_{4}\phi_{3}}\hat{T}\ket{\phi_{2}\phi_{1}}=\int\frac{d^{3}p_{4}d^{3}p_{3}d^{3}p_{2}d^{3}p_{1}}{(2\pi)^{12}16\epsilon_{4}\epsilon_{3}\epsilon_{2}\epsilon_{1}}\bra{\phi_{4}\phi_{3}}\ket{p_{4}p_{3}}\bra{p_{4}p_{3}}\hat{T}\ket{p_{2}p_{1}}\bra{p_{2}p_{1}}\ket{\phi_{2}\phi_{1}}.$
(6)

On the right side, we have inserted two resolutions of the identity given in
equation (LABEL:eq:identity). Inserting equation (5) into the integrand, the
amplitude can be written in terms of the invariant matrix element
$\mathcal{M}$, becoming

$\displaystyle\int\frac{d^{3}p_{4}d^{3}p_{3}d^{3}p_{2}d^{3}p_{1}}{(2\pi)^{12}16\epsilon_{4}\epsilon_{3}\epsilon_{2}\epsilon_{1}}\bra{\phi_{4}}\ket{p_{4}}\bra{\phi_{3}}\ket{p_{3}}\bra{p_{2}}\ket{\phi_{2}}\bra{p_{1}}\ket{\phi_{1}}$
(7) $\displaystyle\qquad\times\mathcal{M}(p_{4}p_{3}\leftarrow
p_{2}p_{1})(2\pi)^{4}\delta(p_{1}+p_{2}-p_{3}-p_{4}).$

Using the delta function to integrate over $\mathbf{p}_{4}$ and $p_{3}^{z}$
yields

$\displaystyle\int\frac{d^{2}p_{3}^{\perp}d^{3}p_{2}d^{3}p_{1}}{(2\pi)^{8}16\epsilon_{2}\epsilon_{1}}\frac{\mathcal{M}(p_{4}p_{3}\leftarrow
p_{2}p_{1})}{\left|p_{3}^{z}\epsilon_{4}-p_{4}^{z}\epsilon_{3}\right|}\bra{\phi_{4}}\ket{p_{4}}\bra{\phi_{3}}\ket{p_{3}}\bra{p_{2}}\ket{\phi_{2}}\bra{p_{1}}\ket{\phi_{1}},$
(8)

where it is understood that $p_{4}$ and $p_{3}$ satisfy
$p_{1}+p_{2}=p_{3}+p_{4}$ (section LABEL:app:p3z). The normalization of
4-momentum states $\ket{p_{n}}$ in terms of 3-momentum states
$\ket{\mathbf{p}_{n}}$ as defined in equation (LABEL:eq:normalization) allows
us to rewrite the expression as

$\displaystyle\int\frac{d^{2}p_{3}^{\perp}d^{3}p_{2}d^{3}p_{1}}{(2\pi)^{8}}\sqrt{\frac{\epsilon_{4}\epsilon_{3}}{\epsilon_{2}\epsilon_{1}}}\frac{\mathcal{M}(p_{4}p_{3}\leftarrow
p_{2}p_{1})}{4\left|p_{3}^{z}\epsilon_{4}-p_{4}^{z}\epsilon_{3}\right|}\bra{\phi_{4}}\ket{\mathbf{p}_{4}}\bra{\phi_{3}}\ket{\mathbf{p}_{3}}\bra{\mathbf{p}_{2}}\ket{\phi_{2}}\bra{\mathbf{p}_{1}}\ket{\phi_{1}}.$
(9)

We can then discretize the momenta by replacing $d^{3}p_{i}/(2\pi)^{3}$ with
$V^{-1}$ and $d^{2}p_{i}^{\perp}/(2\pi)^{2}$ with $A^{-1}$, where $V$ and $A$
are the volume and cross sectional area of the simulated crystal, i.e., the
volume and cross sectional area of the unit cell times the number of k-points
used to sample the BZ. With this, the amplitude for electron states
$\ket{\phi_{1}}$ and $\ket{\phi_{2}}$ to scatter into $\ket{\phi_{3}}$ and
$\ket{\phi_{4}}$ takes the form

$\displaystyle\frac{1}{AV^{2}}\sum_{\mathbf{p}_{3}^{\perp}}\sum_{\mathbf{p}_{2}}\sum_{\mathbf{p}_{1}}\sqrt{\frac{\epsilon_{4}\epsilon_{3}}{\epsilon_{2}\epsilon_{1}}}\frac{\mathcal{M}(p_{4}p_{3}\leftarrow
p_{2}p_{1})}{4\left|p_{3}^{z}\epsilon_{4}-p_{4}^{z}\epsilon_{3}\right|}\bra{\phi_{4}}\ket{\mathbf{p}_{4}}\bra{\phi_{3}}\ket{\mathbf{p}_{3}}\bra{\mathbf{p}_{2}}\ket{\phi_{2}}\bra{\mathbf{p}_{1}}\ket{\phi_{1}}.$
(10)

In the next subsection, we replace $\ket{\phi_{1\dots 4}}$ with states
relevant to electron beam-induced excitation.

###  3.3 Probability of a crystal excitation

We now consider the specific case of beam-induced excitations to determine the
form of the four electron states in equation (10). We assign $\ket{\phi_{1}}$
and $\ket{\phi_{4}}$ to the initial and final beam states
$\ket{\mathbf{p}_{b}}$ and $\ket{\mathbf{p}_{b}^{\prime}}$ respectively.
States $\ket{\phi_{2}}$ and $\ket{\phi_{3}}$ are then the ground and excited
crystal states $\ket{n\mathbf{k}}$ and $\ket{n^{\prime}\mathbf{k^{\prime}}}$
respectively, where $n$ and $n^{\prime}$ are band indices and $\mathbf{k}$ and
$\mathbf{k^{\prime}}$ are k-points. Substituting these specific states into
expression (10), the amplitude for exciting $\ket{n\mathbf{k}}$ to
$\ket{n^{\prime}\mathbf{k^{\prime}}}$ becomes

$\displaystyle\bra{\mathbf{p}_{b}^{\prime},n^{\prime}\mathbf{k^{\prime}}}\hat{T}\ket{n\mathbf{k},\mathbf{p}_{b}}$
$\displaystyle=\frac{1}{AV^{2}}\sum_{\mathbf{p}_{3}^{\perp}}\sum_{\mathbf{p}_{2}}\sum_{\mathbf{p}_{1}}\sqrt{\frac{\epsilon_{4}\epsilon_{3}}{\epsilon_{2}\epsilon_{1}}}\frac{\mathcal{M}(p_{4}p_{3}\leftarrow
p_{2}p_{1})}{4\left|p_{3}^{z}\epsilon_{4}-p_{4}^{z}\epsilon_{3}\right|}$ (11)
$\displaystyle\times\bra{\mathbf{p^{\prime}}_{b}}\ket{\mathbf{p}_{4}}\bra{n^{\prime}\mathbf{k^{\prime}}}\ket{\mathbf{p}_{3}}\bra{\mathbf{p}_{2}}\ket{n\mathbf{k}}\bra{\mathbf{p}_{1}}\ket{\mathbf{p}_{b}}.$

The values of $\epsilon_{1\dots 4}$ need to be clarified before moving
forward. The zeroth components of the initial and final beam momenta obey the
free particle dispersion relations, so
$\epsilon_{1}=\sqrt{|\mathbf{p}_{1}|^{2}+m^{2}}$ and
$\epsilon_{4}=\sqrt{|\mathbf{p}_{4}|^{2}+m^{2}}$. The beam energy that appears
in equation (3) is then defined as the beam electron’s total energy minus its
rest mass: $\epsilon_{b}\equiv\epsilon_{1}-m$. Meanwhile, the momentum of the
crystal states can be treated nonrelativistically. Thus, the zeroth components
of the crystal state momenta are the energy eigenvalues of the crystal state
plus the electron rest mass, i.e., $\epsilon_{2}=\epsilon_{n\mathbf{k}}+m$ and
$\epsilon_{3}=\epsilon_{n^{\prime}\mathbf{k^{\prime}}}+m$. For the remainder
of this derivation, we continue to leave our expressions in terms of
$\epsilon_{1\dots 4}$ for compactness.

Sputtering from a 2D crystal often requires that the beam electron is
backscattered or nearly backscattered, in which case, its final trajectory
after collision with the nucleus is nearly antiparallel to its initial
trajectory and perpendicular to the crystal surface. Given that many 2D
materials (including hBN and MoS2) possess inversion and/or reflection
symmetry about the crystal plane, we assume that the likelihood of excitation
before and after the collision are about equal. In light of this, we calculate
the excitation probability during a sputtering event assuming the beam
electron’s trajectory is not altered by its collision with the nucleus. That
is, we impose that $\mathbf{p}_{1}=|\mathbf{p}_{1}|\hat{\mathbf{z}}$ until an
electronic excitation is induced.

We can now evaluate the bra-ket products in equation (11). The initial beam
state is highly localized on $\mathbf{p}_{b}$, meaning that

$|\bra{\mathbf{p}_{1}}\ket{\mathbf{p}_{b}}|^{2}=V\delta_{\mathbf{p}_{1},\mathbf{p}_{b}}\quad\Rightarrow\quad\bra{\mathbf{p}_{1}}\ket{\mathbf{p}_{b}}=V^{1/2}\delta_{\mathbf{p}_{1},\mathbf{p}_{b}}.$
(12)

Meanwhile, the ground and excited crystal states can be expanded into a plane-
wave basis, so that

$\displaystyle\ket{n\mathbf{k}}=V^{-1/2}\sum_{\mathbf{G}}C_{\mathbf{G+k}}^{n}\ket{\mathbf{G+k}}$
(13)
$\displaystyle\ket{n^{\prime}\mathbf{k^{\prime}}}=V^{-1/2}\sum_{\mathbf{G}}C_{\mathbf{G+k^{\prime}}}^{n^{\prime}}\ket{\mathbf{G+k^{\prime}}},$

where each $\mathbf{G}$ is a reciprocal lattice vector. By re-expressing
$\mathbf{p}_{2}$ and $\mathbf{p}_{3}$ as $\mathbf{G}_{2}+\mathbf{k}_{2}$ and
$\mathbf{G}_{3}+\mathbf{k}_{3}$ respectively, we find

$\displaystyle\bra{\mathbf{p}_{2}}\ket{n\mathbf{k}}=\bra{\mathbf{G}_{2}+\mathbf{k}_{2}}\left(V^{-1/2}\sum_{\mathbf{G}}C_{\mathbf{G+k}}^{n}\ket{\mathbf{G+k}}\right)=V^{1/2}C^{n}_{\mathbf{G}_{2}+\mathbf{k}}\delta_{\mathbf{k}_{2},\mathbf{k}}.$
(14)
$\displaystyle\bra{\mathbf{p}_{3}}\ket{n^{\prime}\mathbf{k^{\prime}}}=\bra{\mathbf{G}_{3}+\mathbf{k}_{3}}\left(V^{-1/2}\sum_{\mathbf{G}}C_{\mathbf{G+k^{\prime}}}^{n^{\prime}}\ket{\mathbf{G+k^{\prime}}}\right)=V^{1/2}C^{n^{\prime}}_{\mathbf{G}_{3}+\mathbf{k^{\prime}}}\delta_{\mathbf{k}_{3},\mathbf{k^{\prime}}}.$

Lastly, we do not care where the outgoing scattered electron ends up, so we
wish for $\ket{\mathbf{p}_{b}^{\prime}}$ to satisfy

$|\bra{\mathbf{p}_{4}}\ket{\mathbf{p}_{b}^{\prime}}|^{2}=V\quad\Rightarrow\quad\bra{\mathbf{p}_{4}}\ket{\mathbf{p}_{b}^{\prime}}=V^{1/2}.$
(15)

The excitation amplitude is then obtained by plugging in the bra-ket products
from equations (12), (14), and (15) into equation (11). This gives us the
excitation amplitude

$\displaystyle\bra{\mathbf{p}_{b}^{\prime},n^{\prime}\mathbf{k^{\prime}}}\hat{T}\ket{n\mathbf{k},\mathbf{p}_{b}}=\frac{1}{A}\sum_{\mathbf{G}_{3}^{\perp}}\sum_{\mathbf{G}_{2}}\sqrt{\frac{\epsilon_{4}\epsilon_{3}}{\epsilon_{1}\epsilon_{2}}}\frac{\mathcal{M}(p_{4}p_{3}\leftarrow
p_{2}p_{1})}{4\left|p_{3}^{z}\epsilon_{4}-p_{4}^{z}\epsilon_{3}\right|}C_{\mathbf{G}_{3}+\mathbf{k^{\prime}}}^{n^{\prime}\ast}C_{\mathbf{G}_{2}+\mathbf{k}}^{n},$
(16)

where it is understood that
$p_{2}=(m+\epsilon_{n\mathbf{k}},\mathbf{G}_{2}+\mathbf{k})$ and
$p_{3}=(m+\epsilon_{n^{\prime}\mathbf{k^{\prime}}},\mathbf{G}_{3}+\mathbf{k^{\prime}})$,
and $p_{4}$ and $p_{3}^{z}$ satisfy $p_{1}+p_{2}=p_{4}+p_{3}$. Squaring this
amplitude yields the probability of a single electronic excitation from the
valence band state $\ket{n\mathbf{k}}$ to the conduction band state
$\ket{n^{\prime}\mathbf{k^{\prime}}}$ for a given beam energy $\epsilon_{b}$,
that is,

$P(n^{\prime}\mathbf{k^{\prime}}\leftarrow
n\mathbf{k}|\epsilon_{b})=\left|\bra{\mathbf{p}_{b}^{\prime},n^{\prime}\mathbf{k^{\prime}}}\hat{T}\ket{n\mathbf{k},\mathbf{p}_{b}}\right|^{2}.$
(17)

###  3.4 Probability of $n_{i}$ excitations

We are finally ready to derive $P_{i}(\epsilon_{b},n_{i})$, the probability
that a beam electron excites a particular number of electrons $n_{i}$. First,
we define the sum of all transition probabilities

$S(\epsilon_{b})\equiv\sum_{\mathbf{k}}\sum_{\mathbf{k^{\prime}}}\sum_{n}\sum_{n^{\prime}}P(n^{\prime}\mathbf{k}^{\prime}\leftarrow
n\mathbf{k}|\epsilon_{b})\equiv\sum_{j}^{N}P_{j}(\epsilon_{b}),$ (18)

where $\mathbf{k}$ and $\mathbf{k^{\prime}}$ run over all k-points, $n$ over
the valence bands, and $n^{\prime}$ over the conduction bands (those
possessing states with energy between the Fermi level and the work function).
The index $j$ in the right-most expression labels the possible single-particle
excitations (e.g., $j=n^{\prime}\mathbf{k^{\prime}}\leftarrow n\mathbf{k}$).

Combinatorics tells us that the probability of exciting exactly one excitation
is

$\displaystyle P_{i}(n_{i}=1)$
$\displaystyle=\sum_{j_{1}}P_{j_{1}}\prod_{j_{2}\neq j_{1}}(1-P_{j_{2}})$ (19)
$\displaystyle=\sum_{j_{1}}P_{j_{1}}\left(1-\sum_{j_{2}\neq
j_{1}}P_{j_{2}}+\frac{1}{2}\sum_{j_{3}\neq j_{1},j_{2}}\sum_{j_{2}\neq
j_{1}}P_{j_{2}}P_{j_{3}}-\dots\right).$

In the large-crystal limit, the number of states, and thus the number of
transitions, is large so that the summations over $j_{i}$ in equation (19) are
approximately equal to one another. That is,

$\sum_{j_{2}\neq j_{1}}P_{j_{2}}\sim\sum_{j_{2}}P_{j_{2}}=S.$ (20)

In this limit, the probability of exactly one beam-induced excitation can be
written as

$\displaystyle
P_{i}(n_{i}=1)\sim\sum_{j}P_{j}\left(1-S+\frac{1}{2}S^{2}-\dots\right)=e^{-S}\sum_{j}P_{j}=Se^{-S}.$
(21)

In the same way, the probability of two excitations is

$\displaystyle P_{i}(n_{i}=2)$
$\displaystyle=\frac{1}{2}\sum_{j_{1}}\sum_{j_{2}\neq
j_{1}}P_{j_{1}}P_{j_{2}}\prod_{j_{3}\neq
j_{1},j_{2}}(1-P_{j_{3}})\sim\frac{S^{2}}{2}e^{-S}.$ (22)

In general, the probability of exactly $n_{i}$ beam-induced excitations is
approximately

$P_{i}(\epsilon_{b},n_{i})\sim\frac{S(\epsilon_{b})^{n_{i}}}{n_{i}!}e^{-S(\epsilon_{b})}.$
(23)

Thus, we see that the probability $P_{i}$ can be written purely in terms of
$S$.

Figure 3:  Probabilities of initially exciting a certain number of electrons
$n_{i}$ in (a) hBN and (b) MoS2 with a beam electron. The excitation
probabilities tend to decreases with increasing beam energy. The excitation
probabilities in MoS2 are larger than those in hBN because MoS2 has a
significantly smaller band gap.

We can now use formula (23) to calculate excitation probabilities in hBN and
MoS2 (figure 3). DFT is used to obtain the plane-wave coefficients
$C_{\mathbf{G}_{2}+\mathbf{k}}^{n}$ and
$C_{\mathbf{G}_{3}+\mathbf{k^{\prime}}}^{n^{\prime}}$ and eigenvalues
$\epsilon_{n\mathbf{k}}$ and $\epsilon_{n^{\prime}\mathbf{k^{\prime}}}$ for
the pristine unit cell of each material. These are plugged into equation (16)
to obtain the amplitude for each transition. We sum over the squares of all
resulting amplitudes to obtain $S$, which is then plugged into formula (23) to
obtain $P_{i}$ for both materials. We emphasize that the only DFT calculation
needed to determine $P_{i}$ is the electronic structure relaxation of a
pristine unit cell, a very inexpensive calculation.

The probabilities plotted in figure 3 reveal some notable trends. First, for
sufficiently large beam energies, $P_{i}(\epsilon_{b},n_{i})$ decreases with
increasing $\epsilon_{b}$ for all $n_{i}>0$. This is because a faster beam
electron has less time to interact with the material and cause an excitation.
In this regime, $P_{i}(\epsilon_{b})$ is proportional to $\epsilon_{b}^{-1}$,
a relationship originally predicted by Bethe [26, 27, 12]. This means that
multiple excitations are more likely at low beam energies. Furthermore, the
probability of remaining in the ground state $P(\epsilon_{b},n_{i}=0)$
vanishes as $\epsilon_{b}$ goes to zero. This implies that a stationary
electron in the vicinity of a material is guaranteed to interact with the
material’s electrons and affect its electronic structure. However, the
validity of formula (23) diminishes as $\epsilon_{b}$ approaches zero. In the
case of a slow beam electron, the interaction between the beam and material
electrons can no longer be approximated by the single virtual photon transfer
processes depicted in figure 2, as the amplitudes for higher-order processes
become more significant. The effects of processes beyond the tree-level should
be the subject of future work. Lastly, the excitation probability is inversely
proportional to the material’s band gap. This is because the zeroth component
of the momentum transfer $q$ depicted in figure 2 is the difference in energy
eigenvalues between the occupied and unoccupied states (section 3.3). Thus,
the smallest possible denominator of the $t$-channel in equation (4) is
proportional to the difference in eigenvalues squared. The experimentally
measured band gap of MoS2 about 1.9 eV [28] while that of hBN is about 6.1 eV
[29]. This means MoS2 hosts transitions with smaller eigenvalue differences,
making the summands in equation (16) larger.

With $P_{i}(\epsilon_{b},n_{i})$ derived in formula (23) and
$d\sigma/dE(\epsilon_{b},E)$ defined in equation (1), we now have two of the
three functions depicted in figure 1 needed to calculate the sputtering cross
section in equation (3). The final ingredients are
$P_{f}(E,\tau,n_{i},n_{f})$, the probability that $n_{f}$ excitations survive
the displacement event given $n_{i}$ initial beam-induced excitations, and
$\\{E_{d}\\}$, the set of all displacement thresholds for all $n_{f}$. The
derivations of these objects are described in the next section.

##  4 Sputtering cross section

We have demonstrated how the interaction between the beam and material
electrons can induce electronic excitations in the material. In this section,
we show how these excitations bring about much larger sputtering cross
sections than those predicted by a ground state theory. We start by showing
how electronic excitations can reduce the displacement threshold $E_{d}$. We
then show that longer excitation lifetimes increase $P_{f}$ for nonzero
$n_{f}$, giving the beam-induced excitations more opportunities to lessen
$E_{d}$. This motivates us to write the sputtering cross section in equation
(3) analytically in terms of the excitation lifetime $\tau$. Finally, the
resulting equation is used to predict the sputtering rates of boron and
nitrogen in hBN and sulfur in MoS2, which can be made to agree well with
experiment for appropriate values of $\tau$ for each material.

###  4.1 Effect of excitations on the displacement threshold

We begin by describing how beam-induced excitations can reduce the
displacement threshold $E_{d}$ in a process called bond softening. As $E_{d}$
is the lower bound of integration over $E$ in equation (3), and the
differential cross section in equation (1) behaves like $E^{-2}$ for small
$E$, reductions in $E_{d}$ can greatly increase the sputtering cross section.
Exactly how excitations change $E_{d}$ is an ambitious study on its own,
requiring a careful consideration of the excited electrons’ evolution and
various relaxation pathways [30]. Here instead, we make three simplifying
assumptions that allow us to calculate $E_{d}$ with only ground state DFT.

Before describing these assumptions, we first define some important terms.
Consider the moment immediately after the beam electron collides with a
material nucleus. The nucleus now has a velocity corresponding to the kinetic
energy $E$ transferred from the beam electron. The resulting nuclear motion
away from its equilibrium position causes the energy of the system to
increase. In this sense, the system climbs an energy surface from the bottom
of its equilibrium well. Far away from the well’s bottom, the energy surface
eventually plateaus. If the system reaches this plateau, the displaced PKA
moves freely away from its initial site without deceleration. At this point,
we consider the PKA to have sputtered. We call the energy at the well’s bottom
$E_{0}$ and that at the plateau $E_{s}$. If the energy surface is static
throughout the entire process, then the displacement threshold is simply
$E_{d}$ = $E_{s}$ \- $E_{0}$. When $E>E_{d}$, the system has enough energy to
climb out of the well, and the PKA sputters. Our task now is to determine how
beam-induced excitations in the material electrons and their subsequent
relaxation affect $E_{0}$ and $E_{s}$. To facilitate this, we make the
following assumptions.

Assumption 1: the excited electrons and holes occupy the band edges. To
justify this, we shift our focus to the material electrons immediately after
the collision. There are now $n_{i}$ electrons in the conduction band and
$n_{i}$ holes in the valence bands. Kretschmer et al. simulated the time-
evolution of these excitations using Ehrenfest dynamics [12, 31]. They found
that the excited electrons relax nonradiatively to the conduction band minimum
(CBM) in a few femtoseconds, while the holes in the valence band take a
similar amount of time to relax to the valence band maximum (VBM). In
contrast, the PKA takes several hundreds of femtoseconds to fully sputter [7].
Thus, the nonradiative relaxation of the electronic structure is essentially
instantaneous, and we can assume that all excited electrons and holes occupy
the CBM and VBM respectively before the PKA has been displaced by an
appreciable amount.

These findings greatly simplify the calculation of $E_{0}$. Given $n_{i}$
excitations, $E_{0}$ is just the energy of the pristine system with $n_{i}$
electrons and holes in the CBM and VBM respectively. In the large crystal
limit, this amounts to adding the pristine band gap $E_{g}$ to the system’s
ground state energy $n_{i}$ times. In other words,
$E_{0}(n_{i})=E_{0}(0)+n_{i}E_{g}$. This approximation of course ignores any
binding energy between the electron and hole. However, we believe that this
formalism allows for an efficient treatment of the lowest order effects of
excitation on $E_{0}$.

Assumption 2: the sputtered system is in its ground state. Finding the plateau
energy $E_{s}$ can be a bit more involved. To calculate it properly, one must
track how the excitations in the CBM and VBM evolve as the PKA moves away from
the crystal. This would require Ehrenfest dynamics of a supercell over
timescales of hundreds of femtoseconds, which is prohibitively expensive. Here
we seek a much less costly set of DFT calculations that can still provide a
reasonable approximation for $E_{s}$. To this end, we draw upon another
finding of Kretschmer et al. [12]. After excitation, the displacement of the
PKA causes the occupied and unoccupied CBM and VBM states to converge into the
band gap and localize on the resulting defect. These converging states are the
bonding and antibonding states that connect the PKA to the host crystal. Thus,
we should consider how the beam-induced excitations affect the electronic
structure on both the PKA and the remaining vacancy.

Figure 4:  Eigenvalues and occupations of the eigenstates of armchair hBN
throughout the sputtering process with respect to the distance between
sputtered (a) boron and (b) nitrogen and their pristine sites. The thickness
of each curve is proportional to the localization of the corresponding state
on the sputtered atom. In hBN, the electronegative nitrogen atoms borrow
charge from the neighboring boron atoms. Thus, the p-orbital states that would
be occupied on an isolated boron atom reside in the conduction band of hBN and
descend into the gap as boron sputters. On the other hand, most of the
p-orbital states of nitrogen lie in the valence band and rise into the gap
during nitrogen sputtering.

We start by considering how the ground state eigenvalues evolve as the PKA
moves away from the crystal. For hBN, the evolution of eigenvalues differ
depending on whether boron or nitrogen is sputtered (figure 4). Boron is
electropositive while nitrogen is electronegative. Thus, in an hBN crystal,
the nitrogen atoms borrow negative charge from the neighboring boron atoms.
This means that the p-orbital states that would be occupied on an isolated
boron atom are vacant, hovering in the conduction band of hBN. Conversely, the
p-orbital states of nitrogen lie occupied in the valence band. For both atoms,
separation from the crystal causes those states to converge to degenerate
p-orbitals on their respective atoms. When boron is sputtered, these states
come from the conduction band. This means that excited electrons residing in
the CBM tend to transfer negative charge to the sputtered boron. However, this
charge transfer is quite energetically unfavorable since boron is
electropositive. The energy cost of this is much larger than the band gap of
hBN, compelling any excess electrons on boron to relax to the host crystal. On
the other hand, the states that localize on sputtered nitrogen must rise from
the valence band. Thus, the holes at the VBM tend to transfer positive charge
to the sputtered nitrogen, which is again energetically unfavorable for the
electronegative atom. Thus, because charge transfer to the sputtered PKA has
such a large energy cost, we presume that in most cases, the PKA becomes
charge neutral by the time sputtering has occurred.

Furthermore, because the states that localize on the PKA converge to the same
degenerate p-orbital, the PKA assumes its ground state after it has sputtered,
regardless of how many electrons the beam excites. We call the energy of the
isolated ground state atom $E_{a}$. Meanwhile, the convergence of the bonding
and antibonding states localized on the remaining vacancy during the
displacement causes the excited electron and hole energy levels to cross
multiples times. This means that there are plenty of chances for an excited
electron to relax nonradiatively as the PKA separates from the material. Thus,
we assume that the host crystal also relaxes to its ground state for most
sputtering events.

Lastly, the sputtering threshold is always larger than the vacancy formation
energy, because some of the energy transferred to the PKA disperses to the
neighboring atoms. For this reason, we calculate the energy of the vacant
system without relaxing the structure, as it has been shown that the free
energy gained by leaving the structure unrelaxed roughly matches the energy
dissipated to the surrounding material [32]. We call the energy of the
unrelaxed ground state vacancy $E_{v}$. Thus, we calculate the plateau energy
using $E_{s}=E_{v}+E_{a}$, the sum of the ground state vacancy and isolated
atom free energies.

Assumption 3: the PKA travels at a constant velocity. Each number of
excitations creates its own energy surface. In this way, electronic excitation
and relaxation cause the system to hop from surface to surface. In our
problem, the initial beam-induced excitation perches the system on an elevated
surface. Once there, electronic relaxation via spontaneous emission (SE)
enables downward surface hopping at the cost of the emitted photon energy. In
this framework, $E_{d}$ is the energy gained along all energy surfaces that
the system traverses during sputtering. For example, suppose the system emits
one photon as it sputters. $E_{d}$ is then the energy gained along the portion
of the excited energy surface traversed before emission plus that gained along
the relaxed surface after emission. This is how SE can change $E_{d}$.
However, the effect of SE on $E_{d}$ must diminish as the PKA moves further
away from the crystal. This means there must be some distance $d$ beyond which
SE no longer affects the forces on the sputtered atom, leaving the energy
surface unchanged. We take $d=4.5$ Å in line with Kretschmer et al.’s study,
though we acknowledge that the exact meaning of this distance is different in
their work [12].

We wish to capture this idea while making the calculation of $P_{f}$ as
intuitive as possible. We do this by assuming that each energy surface is
constant except for a step when the PKA reaches a distance $d$ away from its
equilibrium site. This means that immediately after the collision that excites
$n_{i}$ electrons, the surface assumes a constant energy $E_{0}(n_{i})$. Each
relaxation via SE causes the system to drop to the surface beneath it,
decreasing its energy by $E_{g}$. We define $n_{f}$ as the number of
electronic excitations that survive the PKA’s traversal of distance $d$ to the
energy step. Thus, the energy of the system just before the PKA reaches the
step is $E_{0}(n_{f})$, Beyond the step, all energy surfaces have energy
$E_{s}$, and SE no longer changes the energy surface. With this formalism,
$E_{d}$ is simply the height of the step,

$E_{d}(n_{f})=E_{s}-E_{0}(n_{f})=E_{v}+E_{a}-E_{0}(0)-n_{f}E_{g}.$ (24)

Thus we have derived a simple relationship for how electronic excitations
affect $E_{d}$. However, before we show how this relationship affects the
sputtering cross section, we must first point out an important caveat.
Equation (24) suggests that $E_{d}$ can be negative for large $n_{f}$. In
these cases, the affected atom accelerates away from its pristine site even
without energy transfer from the beam electron. This would seem to suggest
that the sputtering cross section is infinite, since the PKA would sputter for
any energy transfer, no matter how small. This is corroborated by the fact
that the integral in equation (3) diverges as $E_{d}$ approaches zero, meaning
that the sputtering cross section $\sigma$ would approach infinity. However,
this is nonsensical, as it would imply an infinite sputtering rate for a
finite beam current. In reality, even though the beam electron does nudge all
of the material nuclei to some extent, one nucleus always receives a bigger
nudge than the rest. This is the atom onto which the beam-induced excitations
localize. Thus, for each beam electron, only one atom’s displacement threshold
is reduced by electronic excitations. For this reason, the displacement cross
section should never exceed the cross sectional area occupied by a single
target atom. This means that $E_{d}$, the lower integration bound in equation
(3), must have a lower bound itself, which we call $E_{\text{min}}$. When
$E_{d}$ falls below $E_{\text{min}}$, we replace it with $E_{\text{min}}$. We
can approximate $E_{\text{min}}$ by setting the displacement cross section in
equation (2) equal to the area occupied by the PKA and solving for the
displacement threshold (section LABEL:app:Emin). With this, we are finally
ready to use equation (24) to determine the set of displacement thresholds
$\\{E_{d}\\}$ to insert into equation (3).

All positive displacement thresholds computed for this work are listed in
table 1. Explanations for the calculations of $E_{v}$, $E_{a}$, $E_{0}(0)$,
and $E_{g}$ are given in subsctions 4.3 and 4.4. The three largest
displacement thresholds for MoS2 with $n_{f}=0$, 1, and 2 excitations are
similar to those calculated with DFT-based molecular dynamics simulations
[12]. This provides some assurance that our simplified approach to calculating
$E_{d}$ yields reasonable results. Excitation numbers $n_{f}$ greater than
those listed in table 1 make $E_{d}$ negative, in which case the exact value
of $E_{d}$ is not important since it is replaced with $E_{\text{min}}$ in the
calculation of $\sigma$. With that said, it is critical to consider large
$n_{f}$, even if the resulting $E_{d}$ is less than $E_{\text{min}}$. These
large $n_{f}$ can make appreciable contributions to the total cross section,
especially at small beam energies for which $P_{i}(\epsilon_{b},n_{i})$ is
significant for large $n_{i}$ (figure 3). As $n_{f}$ cannot exceed $n_{i}$,
this means that we must consider sufficiently large $n_{i}$ to acknowledge
these contributions. Eventually, these contributions diminish as we increment
$n_{i}$, because the $n_{i}!$ in the denominator of $P_{i}$ in formula (23)
eventually outgrows the $S^{n_{i}}$ in the numerator. We can therefore
truncate the summation over $n_{i}$ in equation (3) at some adequately large
$n_{i}^{\text{max}}$. This means that $n_{i}^{\text{max}}$ must be converged
for each material, as materials with greater $S$ require greater
$n_{i}^{\text{max}}$. In this work, we choose $n_{i}^{\text{max}}$ large
enough so that including more excitations increases the sputtering cross
section by less than 1% for the smallest experimental beam energy considered
for each matieral (figure LABEL:fig:nimax).

Thus, we have shown how $E_{d}$ depends only on the number of surviving
excitations $n_{f}$. We must now determine the likelihood that $n_{f}$
excitations survive given $n_{i}$ beam-induced excitations. This is done by
comparing the excitation lifetime $\tau$ to $t_{s}$, the time it takes for the
PKA to travel a distance $d$. We handle this task in the next subsection,
where we write the sputtering cross section in equation (3) in terms of
$\tau$.

| displacement threshold (eV)
---|---
$n_{f}$ | 0 | 1 | 2 | 3
B from hBN | 12.85 | 8.78 | 4.71 | 0.64
N from hBN | 12.71 | 8.64 | 4.57 | 0.50
S from MoS2 | 6.92 | 5.04 | 3.16 | 1.28
Table 1:  All computed positive displacement thresholds for sputtering from
the hBN armchair edge and MoS2 surface. The thresholds decrease as the number
of surviving excitations $n_{f}$ increases.

###  4.2 Sputtering cross section in terms of the excitation lifetime

In hBN and MoS2, occupation of the antibonding state localized on the PKA does
not affect the free energy of the sputtered system, since the antibonding and
bonding states to converge to the same degenerate p-orbital in the limit that
the PKA is isolated from the host material. Therefore, only the excitation
lifetime of the host material needs be considered in calculating
$P_{f}(n_{i},n_{f})$, the probability that $n_{f}$ of the $n_{i}$ beam-induced
excitations survive long enough to reduce the displacement threshold.

We define the ratio of surviving excitations as

$R(E,\tau)=e^{-t_{s}(E)/\tau},$ (25)

where $t_{s}$ is the time it takes for the sputtered atom to travel a distance
$d$, and the excitation lifetime $\tau$ is determined for each material by
fitting the cross section to experimental data. Given $n_{i}$ beam-induced
excitations, the probability that $n_{f}$ excitations survive is

$\displaystyle P_{f}(n_{i},n_{f})$ $\displaystyle=\begin{pmatrix}n_{i}\\\
n_{f}\end{pmatrix}R^{n_{f}}(1-R)^{n_{i}-n_{f}}=\begin{pmatrix}n_{i}\\\
n_{f}\end{pmatrix}\sum^{n_{i}-n_{f}}_{n=0}\begin{pmatrix}n_{i}-n_{f}\\\
n\end{pmatrix}(-1)^{n}R^{n_{f}+n}.$ (26)

Equation (26) allows us to rewrite the integral in equation (3) as

$\displaystyle\int_{E_{d}}^{E_{\text{max}}}P_{f}\frac{d\sigma}{dE}dE=\begin{pmatrix}n_{i}\\\
n_{f}\end{pmatrix}\sum^{n_{i}-n_{f}}_{n=0}\begin{pmatrix}n_{i}-n_{f}\\\
n\end{pmatrix}(-1)^{n}\int^{E_{\text{max}}}_{E_{d}}R^{n_{f}+n}\frac{d\sigma}{dE}dE.$
(27)

Assumption 3 from the previous subsection again aids us here. Because the PKA
travels at a constant velocity until it reaches the sputtering distance $d$,
we can write $t_{s}$ from equation (25) quite simply as $t_{s}=d\sqrt{M/2E}$.
This allows for the straightforward analytical integration of equation (27).
Using $d\sigma/dE$ defined in equation (1), the integral on the right can be
written as

$\displaystyle\int R^{n_{f}+n}\frac{d\sigma}{dE}dE$
$\displaystyle=\pi\left(\frac{Z\alpha}{|\mathbf{p}|\beta}\right)^{2}\bigg{[}\frac{4\tau^{2}E_{\text{max}}}{d^{2}M}\frac{\mu\xi+1}{\mu^{2}}e^{-\mu\xi}$
(28) $\displaystyle+2\beta(\pi Z\alpha+\beta)\text{Ei}(-\mu\xi)+\frac{\pi
Z\alpha\beta\tau}{\mu
d}\sqrt{\frac{8E_{\text{max}}}{M}}e^{-\mu\xi}\bigg{]}+\text{const}.,$

where we define

$\xi\equiv\frac{-d}{\tau}\sqrt{\frac{M}{2E}}\quad\text{and}\quad\mu\equiv
n_{f}+n,$ (29)

and the function Ei is the exponential integral,

$\text{Ei}(x)=\int_{-\infty}^{x}\frac{e^{t}}{t}dt.$ (30)

Looking back at equation (3), we now have everything we need to evaluate the
sputtering cross section. We derived $P_{i}(\epsilon_{b},n_{i})$ in formula
(23), and we performed the integration over $E$ analytically in equations (27)
and (28), setting $\\{E_{d}\\}$ in equation (24) as the lower integration
bounds. We also have a criterion to truncate the summation over $n_{i}$ at
$n_{i}^{\text{max}}$ for a given material, as described at the end of section
4.1. In the following subsections, we use these results to calculate the
sputtering cross sections of hBN and MoS2.

###  4.3 Boron and nitrogen sputtering from hexagonal boron nitride

Electron beam irradiation has been shown to bore nanopores in monolayer hBN at
beam energies far beneath the calculated ground state critical energy of
$\epsilon_{c}\sim 80$ keV [33, 34, 8, 9]. In a pristine hBN layer, these beam-
induced pores can be initialized from isolated boron and nitrogen vacancies.
The atoms surrounding these vacancies have reduced coordination numbers, and
thus, smaller displacement thresholds than those in the pristine material.
Therefore, the atoms lining the defect are more likely to sputter than their
surface counterparts for a given beam energy. As the edge atoms continue to
sputter away at a high rate, the nanopore grows, eventually extending up to a
few nm in diameter [34, 9, 35].

Cretu et al. measured the radial growth of these nanopores under electron
irradiation for several temperatures ranging from 673 to 1473 K [9]. By
dividing the radial growth rate by the beam current, they estimated the
sputtering cross section to be around 25 barn under both 30 and 60 keV beams
at temperatures of 1273 K and below (figure 5c). The cross sections were
fairly temperature independent at these relatively low temperatures. Cretu et
al. also found that the edges of these pores most often assume an armchair
structure. Therefore, in an attempt to reproduce these measurements, we
calculate the cross section of boron and nitrogen sputtering from an armchair
edge at 1273 K (figure 5). To evaluate the displacement thresholds defined in
equation (24) and listed in table 1, $E_{v}$ is the free energy of the
armchair edge supercell with a single boron or nitrogen vacancy, while $E_{a}$
is that of an isolated boron or nitrogen atom. $E_{0}(0)$ and $E_{g}$ are then
the calculated free energy and band gap of the pristine armchair edge
supercell. We plug the resulting set of displacement thresholds $\\{E_{d}\\}$
into equation (3) to calculate the sum of the boron and nitrogen sputtering
cross sections. The excitation probabilities of hBN plotted in figure 3a are
then used for $P_{i}$. Strictly speaking, the calculation of $P_{i}$ should
consider the effects of the edge state orbitals. However, in our formalism,
the beam electron is in a momentum eigenstate that is highly delocalized in
real space. We therefore assume that the radius of the beam is much larger
than that of the nanopore. This means that the majority of the beam-matter
interactions occur in regions of pristine material, validating the use of the
$P_{i}$ calculated for pristine hBN. Future work should consider the effects
of localized beam electron states to simulate beams with smaller focal points.
Lastly, temperature effects on the cross section are considered in the manner
described in our previous work [7].

Figure 5:  Sputtering from the armchair edge of hBN. The schematics in panels
(a) and (b) depict the sputtering of nitrogen from out-of-plane and in-plane
perspectives respectively. Panel (c) plots the sputtering cross sections of
both boron and nitrogen in green and blue respectively. The black curves are
the total sputtering cross section, i.e., the sum of their corresponding blue
and green curves. The solid curves account for beam-induced excitations
assuming an excitation lifetime of $\tau=240$ fs (labelled “exc.” for
excited). The dashed curves ignore the possibility of beam-induced excitation
(labelled “gro.” for ground). The black squares are experimentally observed
sputtering cross sections for hBN at 1273 K [9]. The consideration of beam-
induced excitations reduces the disagreement between theory and experiment
substantially.

Fitting our cross section curves to the data of Cretu et al. yields an
excellent agreement if the predicted excitation lifetime is set to $\tau\sim
240$ fs. This predicted lifetime is much shorter than the reported excitation
lifetime of $\sim$0.75 ns in pristine hBN [36], indicating that the sputtering
process can significantly reduce the excitation lifetimes of hBN. We suspect
that the atomic motion gives rise to non-radiative relaxation pathways that
are not explicitly accounted for here. This motivates a closer investigation
into the full electronic evolution of the hBN system post-collision. However,
such a study is beyond the scope of this work. With that said, the electronic
structure in the vicinity of a sputtering PKA differs significantly from that
of the pristine room-temperature systems in which excitation lifetimes are
experimentally measured. There is therefore no reason to expect the predicted
lifetimes in this work to match those obtained by experiment.

We also see that the sputtering cross section, and thus the growth rate of the
nanopore, is minimized for beam energies between 30 and 60 keV. Below these
energies, the sputtering cross section begins to grow as the beam energy
decreases. This is due to a non-negligible probability of final excitation
numbers $n_{f}>3$, for which $E_{d}(n_{f})$ falls below $E_{\text{min}}$
(table 1). In these cases, the sputtering cross section is $\Omega\sim 5\times
10^{8}$ barn, seven orders of magnitude larger than the measured cross section
at 30 keV. This suggests that one can expect the beam-induced nanopores in hBN
to grow under beam energies as low as 1 keV. However, one must again be
cautious of the predicted cross section at low beam energies in which the
tree-level theory starts to break down. Nevertheless, figure 5c demonstrates a
strong beam-energy dependence in the sputtering cross section across a wide
range of experimentally relevant beam-energies. These findings could
facilitate precise control of nanopore growth rates in hBN under electron
irradiation.

###  4.4 Sulfur sputtering from molybdenum disulfide

Kretschmer et al. measured the sputtering cross section of sulfur from MoS2
for several beam energies ranging from 20 to 80 keV [12]. They found a peak in
the cross section at 30 keV, much less than the predicted ground state
$\epsilon_{c}\sim 90$ keV [12]. To help explain this unexpected peak, we
calculate the cross section for sulfur sputtering from pristine MoS2 (figure
6). The vacant system free energy $E_{v}$ for the calculation of $\\{E_{d}\\}$
is that of a MoS2 supercell with a single sulfur vacancy. $E_{a}$ is then the
free energy of an isolated sulfur atom, and $E_{0}(0)$ is that of the pristine
MoS2 supercell. We then set $E_{g}$ equal to the experimental band gap of 1.88
eV [28]. Using the resulting $\\{E_{d}\\}$, we find that summing over the
contributions of final excitation numbers $n_{f}>2$ to the sputtering cross
section produces a peak just below 30 keV, matching Kretschmer et al.’s
findings remarkably well. In fitting to this peak, we predict an excitation
lifetime of $\tau\sim 81$ fs. This is again much shorter than the excitation
lifetime of pristine MoS2, which is on the order of a few picoseconds [37, 38,
39]. However, this is also much shorter than the fitted lifetime of hBN found
in the previous subsection, consistent with the fact that the excitation
lifetime of pristine hBN is much longer than that of MoS2.

The difference between the two materials’ lifetimes leads to markedly
dissimilar cross section behaviors at low beam energies. Below beam energies
of 30 keV, the cross section of MoS2 gradually drops to zero with decreasing
$\epsilon_{b}$. In contrast, hBN’s total cross section has a minimum at around
40 keV and begins to increase as $\epsilon_{b}$ decreases. Eventually, hBN’s
sputtering cross section peaks before dropping quickly to zero as the beam
energy goes to zero (figure LABEL:fig:edgePeaks). However, this peak occurs at
around 0.5 keV, far below the lower energy bound of figure 5c. These distinct
cross section behaviors can be explained by the amplified sensitivity of
$n_{f}$ to $\tau$ at low beam energies. Equations (25) and (26) tell us that
larger $\tau$ makes large $n_{f}$ more likely. At the same time, $n_{f}$
cannot exceed $n_{i}$. Thus, $n_{f}$ is more sensitive to $\tau$ at low beam
energies for which $n_{i}$ is large. Accordingly, because $\tau$ is greater in
hBN than in MoS2, the expected values of $n_{f}$ in hBN are much larger than
those of MoS2 at low beam energies. This explains why the effect of
considering excitations is much more pronounced in hBN than in MoS2.
Furthermore, the difference in the cross section behavior is exacerbated by
the fact that sulfur is heavier than both boron and nitrogen, so that its
post-collision velocity is always smaller for a given energy transfer $E$. It
follows that $t_{s}$ is larger, and thus, $P_{f}$ is smaller for a given $E$
in MoS2. This again increases the likelihood that hBN has more final
excitations than MoS2, meaning that the effects of beam-induced excitation on
hBN’s sputtering cross section are greater.

Figure 6:  Sputtering of sulfur from the outgoing MoS2 surface. The sputtering
of molybdenum from MoS2 can be ignored because its displacement threshold is
significantly larger than that of sulfur. The schematics in panels (a) and (b)
depict the sputtering process from out-of-plane and in-plane perspectives
respectively. Panel (c) plots the contributions of various numbers of final
excitations $n_{f}$ to the sputtering cross section assuming an excitation
lifetime of $\tau=81$ fs. The black dashed curve plots the predicted cross
section ignoring the possibility of beam-induced excitation. The black squares
are experimentally observed cross sections at 300 K [12]. The contribution
from the sum of all $n_{f}>2$ final excitations matches the experimental cross
section remarkably well.

We also plot the contributions of $n_{f}=0$, 1, and 2 excitations to MoS2’s
sputtering cross section separately. In doing so, we see that the
contributions of $n_{f}=1$ or 2 would conceal the peak at 30 keV. Thus, it
seems that the likelihoods of $n_{f}=1$ or 2 are somehow suppressed. This
suggests that the individual beam-induced excitations and subsequent
relaxation are in some cases correlated. That is, the excitation and/or
relaxation rates of a given electronic transition are affected by the
coinciding distribution of electronic excitations. Thus, a proper treatment of
these excitations and their effect on the sputtering cross section requires
that the excitation probabilities of every possible transition are calculated
for every possible excitation configuration. Such a nonlinear calculation is
beyond the scope of this work, but should certainly be pursued in a future
study. Nonetheless, the methods laid out here demonstrate that the
consideration of beam-induced excitation can provide a quantitative
justification for the sulfur sputtering rate to peak at a beam energy well-
below the expected ground state critical energy.

##  5 Conclusion

In this paper, we developed a first-principles method to more accurately
describe electron beam-induced sputtering cross sections in 2D insulating
crystals by accounting for beam-induced electronic excitations and their
subsequent relaxations. The method combines QED scattering theory with DFT
electronic structure calculations to determine the likelihood of beam-induced
excitation. The results show that the excitation probability is inversely
proportional to both the material’s band gap and beam electron’s kinetic
energy. We then show how these nonzero excitation probabilities increase the
predicted sputtering cross sections of both hBN in MoS2. These cross sections
can be made quantitatively similar to those obtained experimentally by
treating the excitation lifetime as a fitting parameter. The methods laid out
are computationally efficient, requiring only a few ground state electronic
optimization calculations for each cross section curve. Thus, the formalism
that we have developed can be easily applied to any 2D crystalline material to
simulate the rates of atomic displacement under electron irradiation.

With that said, several questions naturally arise from our study. For example,
why is the excitation lifetime reduced so drastically during sputtering? How
might excitation and relaxation rates of different transitions be correlated?
How might preexisting defects affect those rates? These questions urge follow
up work to address the full breadth of physical processes involved in beam-
matter interactions. Future studies should also consider additional electronic
relaxation pathways to determine their effect on $P_{f}$. Correspondingly,
other electronic responses such as ionization, core excitations, and second
order electronic excitation effects such as Auger scattering can be
incorporated into the calculation of $P_{i}$ [38, 40, 41, 42]. Moreover, spin
polarization effects can also be examined by making $\mathcal{M}$ spin-
dependent, i.e., not averaging over spins as is done in equation (4). The beam
electron path can also change significantly after collision with the nucleus.
Ensuing research should investigate how these altered trajectories generate
new excitation probabilities $P_{i}$, which must certainly play a role in 3D
bulk materials. Furthermore, the methods here can be modified to accommodate
more advanced DFT techniques. For example, the calculation of $P_{i}$ should
be made compatible with ultrasoft and projector augmented wave
pseudopotentials [43] in a manner similar to that employed in modern GW codes
[44, 45, 46]. In addition, the plane-wave coefficients of excited electronic
states can be calculated self-consistently using constrained DFT [47]. Perhaps
most importantly, progress in this field requires much more experimental data.
We therefore hope this paper encourages new experimental investigation into
beam-induced sputtering for beam energies below the predicted ground state
critical energy.

Clearly, the combined QED-DFT approach to modeling beam-induced excitations
and their effect on the sputtering cross section opens up a rich and diverse
field of physics for both theoretical and experimental exploration. We hope
that this work and the work it may stimulate can eventually enable the use of
electron beams for precise atomic-scale engineering of crystalline materials.

##  6 Methods

DFT [13, 14] was used to determine all electronic and ionic structures from
first-principles. Free energy calculations were carried out with the Vienna ab
initio Simulation Package (VASP) [48, 15] implementing the projector augmented
wave (PAW) method [43] along with the Perdew-Burke-Ernserhof (PBE) generalized
gradient approximation (GGA) to the exchange correlation functional [49]. Van
der Waals interactions were accounted for using the optB88-vdW density
functional methods [50, 51]. All parameters were converged so that any
increase in precision would change the total free energy by less than 1 meV
per atom. The cutoff energies for hBN and MoS2 were set to 800 and 550 eV
respectively, while the BZs of both materials’ pristine unit cells were
sampled with a $\Gamma$-centered $6\times 6\times 1$ Monkhorst-Pack meshes
[52], corresponding to 0.417 and 0.328 and k-points per inverse Å for hBN and
MoS2 respectively. To achieve the same k-point densities, surface vacancies in
MoS2 and edge vacancies in hBN were placed in respective $6\times 6\times 1$
and $4\times 1\times 1$ supercells whose BZs were sampled with a single
k-point on $\Gamma$. The heights of the hBN and MoS2 cells were 12 Å and 20 Å
respectively to provide sufficient separation from periodic images. Nanoribbon
structures were used to simulate isolated armchair edges in hBN. These ribbons
were more than 16 Å across and placed in cells 28 Å wide to avoid interactions
between opposing edges and periodic images. Lastly, relaxation iterations of
ionic positions and lattice constants persisted until the all Hellmann-Feynman
forces settled below 1 meV/Å.

Because equation (16) relies on the orthogonality of the Kohn-Sham orbitals,
the optimized norm-conserving Vanderbilt pseudopotentials [53, 54] implemented
in Quantum ESPRESSO [55] were used to determine all plane-wave coefficients
$C^{n}_{\mathbf{G+k}}$. The sum of all transition probabilities $S$ defined in
equation (18) was used to gage the convergence of all parameters, which were
deemed converged when any increase in precision changed $S$ by less than 5%
(figure LABEL:fig:convergences). The parameter values that satisfy this
criteria differed substantially from those needed for free energy convergence.
The cutoff energies for hBN and MoS2 were set to 299 and 286 eV respectively,
while their cell-heights were respectively set to 18 and 12 Å. For both
materials, the maximum virtual photon momentum $|\mathbf{q}_{\text{max}}|$
required for convergence fell well within their first BZs. We therefore chose
$|\mathbf{q}_{\text{max}}|$ to be the magnitude of high-symmetry point $K$ in
each respective BZ. The maximum number of initial excitations considered for
hBN and MoS2 were $n_{i}^{\text{max}}=5$ and 9 respectively (figure
LABEL:fig:nimax). Finally, convergence of $S$ requires extremely dense k-point
sampling of the BZ. This necessitates fitting a curve to $S$ calculated for
various k-point mesh densities and extrapolating to an infinitely fine mesh to
estimate the converged value of $S$ (section LABEL:app:fitting). The most
dense k-point meshes used to fit these curves had dimensions of $45\times
45\times 1$ and $36\times 36\times 1$ for hBN and MoS2 respectively.

##  7 Acknowledgements

Calculations were performed at the Center for Computational Innovations at
Rensselaer Polytechnic Institute, Livermore Computing Center at Lawrence
Livermore National Laboratory (LLNL), and Compute and Data Environment for
Science at Oak Ridge National Laboratory (ORNL), which is supported by the
Office of Science of the U.S. Department of Energy (DOE) under Contract No.
DE-AC05-00OR22725. This work was performed under the auspices of the U.S.
Department of Energy by LLNL under Contract DE-AC52-07NA27344. This material
is also based upon work supported by the U.S. DOE, Office of Science, Office
of Workforce Development for Teachers and Scientists, Office of Science
Graduate Student Research (SCGSR) program. The SCGSR program is administered
by the Oak Ridge Institute for Science and Education for the DOE under
contract number DE‐SC0014664. Funding was provided by the National Science
Foundation (Award 1608171). Work (BGS, PG, DL, JJ) was also supported by
ORNL’s Center for Nanophase Materials Sciences, a U.S. DOE Office of Science
User Facility.

## References

  * [1] F. Banhart, “Irradiation effects in carbon nanostructures,” Rep. Prog. Phys., vol. 62, pp. 1181–1221, jul 1999.
  * [2] R. F. Egerton, “Ultramicroscopy Control of radiation damage in the TEM,” Ultramicroscopy, vol. 127, pp. 100–108, apr 2013.
  * [3] X. Zhao, J. Kotakoski, J. C. Meyer, E. Sutter, P. Sutter, A. V. Krasheninnikov, U. Kaiser, and W. Zhou, “Engineering and modifying two-dimensional materials by electron beams,” MRS Bull., vol. 42, pp. 667–676, sep 2017\.
  * [4] T. Susi, J. C. Meyer, and J. Kotakoski, “Quantifying transmission electron microscopy irradiation effects using two-dimensional materials,” Nat. Rev. Phys., vol. 1, pp. 397–405, may 2019.
  * [5] J. C. Meyer, F. Eder, S. Kurasch, V. Skakalova, J. Kotakoski, H. J. Park, S. Roth, A. Chuvilin, S. Eyhusen, G. Benner, A. V. Krasheninnikov, and U. Kaiser, “Accurate measurement of electron beam induced displacement cross sections for single-layer graphene,” Phys. Rev. Lett., vol. 108, pp. 196102–196107, may 2012.
  * [6] T. Susi, C. Hofer, G. Argentero, G. T. Leuthner, T. J. Pennycook, C. Mangler, J. C. Meyer, and J. Kotakoski, “Isotope analysis in the transmission electron microscope,” Nat. Commun., vol. 7, p. 13040, 2016.
  * [7] A. Yoshimura, M. Lamparski, N. Kharche, and V. Meunier, “First-principles simulation of local response in transition metal dichalcogenides under electron irradiation,” Nanoscale, vol. 10, pp. 2388–2397, jan 2018.
  * [8] J. Kotakoski, C. H. Jin, O. Lehtinen, K. Suenaga, and A. V. Krasheninnikov, “Electron knock-on damage in hexagonal boron nitride monolayers,” Phys. Rev. B, vol. 82, p. 113404, sep 2010.
  * [9] O. Cretu, Y. C. Lin, and K. Suenaga, “Inelastic electron irradiation damage in hexagonal boron nitride,” Micron, vol. 72, pp. 21–27, may 2015.
  * [10] Y.-C. Lin, T. Björkman, H.-P. Komsa, P.-Y. Teng, C.-H. Yeh, F.-S. Huang, K.-H. Lin, J. Jadczak, Y.-S. Huang, P.-W. Chiu, A. V. Krasheninnikov, and K. Suenaga, “Three-fold rotational defects in two-dimensional transition metal dichalcogenides.,” Nat. Commun., vol. 6, pp. 6736–6741, oct 2015\.
  * [11] T. Lehnert, O. Lehtinen, G. A. Siller, U. Kaiser, G. Algara-Siller, and U. Kaiser, “Electron radiation damage mechanisms in 2D MoSe2,” Appl. Phys. Lett., vol. 110, p. 033106, jan 2017.
  * [12] S. Kretschmer, T. Lehnert, U. Kaiser, A. V. Krasheninnikov, and A. V. Krasheninnikov, “Formation of defects in two-dimensional MoS2 in transmission electron microscope at electron energies below the knock-on threshold : the role of electronic excitations,” Nano Lett., vol. 20, pp. 2865–2870, mar 2020.
  * [13] P. Hohenberg and W. Kohn, “Inhomogeneous electron gas,” Phys. Rev., vol. 136, pp. B864–871, nov 1964.
  * [14] W. Kohn and L. J. Sham, “Self-consistent equations including exchange and correlation effects,” Phys. Rev., vol. 140, pp. A1133–1138, nov 1965\.
  * [15] G. Kresse and J. Furthmüller, “Efficient iterative schemes for ab initio total-energy calculations using a plane-wave basis set,” Phys. Rev. B, vol. 54, pp. 11169–11186, oct 1996.
  * [16] H. Lourenço-Martins, A. Lubk, and M. Kociak, “Bridging nano-optics and condensed matter formalisms in a unified description of inelastic scattering of relativistic electron beams,” SciPost Phys., vol. 10, no. 2, p. 31, 2021.
  * [17] M. E. Peskin and D. V. Schroeder, An Introduction to Quantum Field Theory.  Reading, MA: Perseus Books Publishing, L.L.C., first ed., 1995.
  * [18] T. Lancaster and S. J. Blundell, Quantum Field Theory for the Gifted Amateur.  New York, NY, United States of America: Oxford University Press, first ed., 2014.
  * [19] N. F. Mott, “The Scattering of Fast Electrons by Atomic Nuclei,” Proc. R. Soc. A, vol. 124, pp. 425–442, jun 1929.
  * [20] W. A. McKinley and H. Feshbach, “The coulomb scattering of relativistic electrons by nuclei,” Phys. Rev., vol. 74, pp. 1759–1763, dec 1948.
  * [21] O. S. Oen, “Cross sections for atomic displacements in solids by fast electrons,” tech. rep., Oak Ridge National Laboratory, Oak Ridge, Tennessee, aug 1973.
  * [22] R. F. Egerton, R. McLeod, F. Wang, and M. Malac, “Basic questions related to electron-induced sputtering in the TEM,” Ultramicroscopy, vol. 110, pp. 991–997, jul 2010.
  * [23] C. Møller, “Zur Theorie des Durchgangs schneller Elektronen durch Materie,” Ann. Phys., vol. 406, no. 5, pp. 531–585, 1932.
  * [24] H. Kragh, “Relativistic collisions: The work of Christian Møller in the early 1930s,” Arch. Hist. Exact Sci., vol. 43, no. 4, pp. 299–328, 1992\.
  * [25] X. Roqué, “Møller scattering: a neglected application of early quantum electrodynamics,” Arch. Hist. Exact Sci., vol. 44, pp. 197–264, sep 1992\.
  * [26] H. Bethe, “Zur Theorie des Durchgangs schneller Korpuskularstrahlen durch Materie,” Ann. Phys., vol. 397, no. 3, pp. 325–400, 1930.
  * [27] M. Inokuti, Y.-K. Kim, and R. L. Platzman, “Total Cross Sections for Inelastic Scattering of Charged Particles by Atoms and Molecules. I. A Sum Rule for the Bethe Cross Sections and Its Application to the Helium Atom,” Phys. Rev., vol. 164, pp. 55–61, dec 1967.
  * [28] J. Gusakova, X. Wang, L. L. Shiau, A. Krivosheeva, V. Shaposhnikov, V. Borisenko, V. Gusakov, and B. K. Tay, “Electronic Properties of Bulk and Monolayer TMDs: Theoretical Study Within DFT Framework (GVJ-2e Method),” Phys. status solidi, vol. 214, no. 12, p. 1700218, 2017.
  * [29] C. Elias, P. Valvin, T. Pelini, A. Summerfield, C. J. Mellor, T. S. Cheng, L. Eaves, C. T. Foxon, P. H. Beton, S. V. Novikov, B. Gil, and G. Cassabois, “Direct band-gap crossover in epitaxial monolayer boron nitride,” Nat. Commun., vol. 10, no. 1, p. 2639, 2019.
  * [30] D. B. Lingerfelt, P. Ganesh, J. Jakowski, and B. G. Sumpter, “Electronically Nonadiabatic Structural Transformations Promoted by Electron Beams,” Adv. Funct. Mater., vol. 1901901, p. 1901901, 2019.
  * [31] X. Li, J. C. Tully, H. B. Schlegel, and M. J. Frisch, “Ab initio Ehrenfest dynamics,” J. Chem. Phys., vol. 123, no. 8, p. 84106, 2005.
  * [32] H. P. Komsa, J. Kotakoski, S. Kurasch, O. Lehtinen, U. Kaiser, and A. V. Krasheninnikov, “Two-dimensional transition metal dichalcogenides under electron irradiation: Defect production and doping,” Phys. Rev. Lett., vol. 109, pp. 035503–1–5, jul 2012.
  * [33] C. Jin, F. Lin, K. Suenaga, and S. Iijima, “Fabrication of a freestanding boron nitride single layer and Its defect assignments,” Phys. Rev. Lett., vol. 102, no. 19, p. 195505, 2009.
  * [34] J. C. Meyer, A. Chuvilin, G. Algara-Siller, J. Biskupek, and U. Kaiser, “Selective sputtering and atomic resolution imaging of atomically thin boron nitride membranes,” Nano Lett., vol. 9, no. 7, pp. 2683–2689, 2009\.
  * [35] M. Dogan, S. M. Gilbert, T. Pham, B. Shevitski, P. Ercius, S. Aloni, A. Zettl, and M. L. Cohen, “Electron beam-induced nanopores in Bernal-stacked hexagonal boron nitride,” Appl. Phys. Lett., vol. 117, no. 2, pp. 9–13, 2020.
  * [36] J. Li, X. K. Cao, T. B. Hoffman, J. H. Edgar, J. Y. Lin, and H. X. Jiang, “Nature of exciton transitions in hexagonal boron nitride,” Appl. Phys. Lett., vol. 108, pp. 122101–1–4, mar 2016.
  * [37] T. Korn, S. Heydrich, M. Hirmer, J. Schmutzler, C. Schüller, and C. Schller, “Low-temperature photocarrier dynamics in monolayer MoS2,” Appl. Phys. Lett., vol. 99, no. 10, p. 102109, 2011.
  * [38] D. Lagarde, L. Bouet, X. Marie, C. R. Zhu, B. L. Liu, T. Amand, P. H. Tan, and B. Urbaszek, “Carrier and Polarization Dynamics in Monolayer MoS2,” Phys. Rev. Lett., vol. 112, p. 047401, jan 2014.
  * [39] M. Palummo, M. Bernardi, and J. C. Grossman, “Exciton Radiative Lifetimes in Two-Dimensional Transition Metal Dichalcogenides,” Nano Lett., vol. 15, pp. 2794–2800, may 2015.
  * [40] D. Kozawa, R. Kumar, A. Carvalho, K. Kumar Amara, W. Zhao, S. Wang, M. Toh, R. M. Ribeiro, A. H. Castro Neto, K. Matsuda, and G. Eda, “Photocarrier relaxation pathway in two-dimensional semiconducting transition metal dichalcogenides,” Nat. Commun., vol. 5, no. 1, p. 4543, 2014.
  * [41] Z. Nie, R. Long, J. S. Teguh, C.-C. Huang, D. W. Hewak, E. K. L. Yeow, Z. Shen, O. V. Prezhdo, and Z.-H. Loh, “Ultrafast Electron and Hole Relaxation Pathways in Few-Layer MoS2,” J. Phys. Chem. C, vol. 119, pp. 20698–20708, aug 2015.
  * [42] H. Shi, R. Yan, S. Bertolazzi, J. Brivio, B. Gao, A. Kis, D. Jena, H. G. Xing, and L. Huang, “Exciton dynamics in suspended monolayer and few-layer MoS2 2D crystals,” ACS Nano, vol. 7, pp. 1072–1080, jan 2013.
  * [43] P. E. Blöchl, “Projector augmented-wave method,” Phys. Rev. B, vol. 50, pp. 17953–17979, dec 1994.
  * [44] M. Shishkin and G. Kresse, “Implementation and performance of the frequency-dependent GW method within the PAW framework,” Phys. Rev. B, vol. 74, p. 035101, jul 2006.
  * [45] M. Gajdoš, K. Hummer, G. Kresse, J. Furthmüller, and F. Bechstedt, “Linear optical properties in the projector-augmented wave methodology,” Phys. Rev. B, vol. 73, p. 045112, jan 2006.
  * [46] J. Paier, R. Hirschl, M. Marsman, and G. Kresse, “The Perdew–Burke–Ernzerhof exchange-correlation functional applied to the G2-1 test set using a plane-wave basis set,” J. Chem. Phys., vol. 122, no. 23, p. 234102, 2005.
  * [47] B. Kaduk, T. Kowalczyk, and T. Van Voorhis, “Constrained Density Functional Theory,” Chem. Rev., vol. 112, pp. 321–370, jan 2012.
  * [48] G. Kresse and J. Furthmüller, “Efficiency of ab-initio total energy calculations for metals and semiconductors using a plane-wave basis set,” Comput. Mater. Sci., vol. 6, pp. 15–50, jul 1996.
  * [49] J. P. Perdew, K. Burke, and M. Ernzerhof, “Generalized Gradient Approximation Made Simple,” Phys. Rev. Lett., vol. 77, pp. 3865–3868, oct 1996.
  * [50] J. Klimeš, D. Bowler, and A. Michaelides, “Chemical accuracy for the van der Waals density functional,” J. Phys. Condens. Matter, vol. 22, no. 2, pp. 022201–022205, 2009.
  * [51] J. Klimeš, D. R. Bowler, and A. Michaelides, “Van der Waals density functionals applied to solids,” Phys. Rev. B, vol. 83, pp. 195131–1–13, may 2011.
  * [52] H. Monkhorst and J. Pack, “Special points for Brillouin zone integrations,” Phys. Rev. B, vol. 13, pp. 5188–5192, jun 1976.
  * [53] D. R. Hamann, “Optimized norm-conserving Vanderbilt pseudopotentials,” Phys. Rev. B, vol. 88, p. 085117, aug 2013.
  * [54] M. Schlipf and F. Gygi, “Optimization algorithm for the generation of ONCV pseudopotentials,” Comput. Phys. Commun., vol. 196, pp. 36–44, 2015.
  * [55] P. Giannozzi, S. Baroni, N. Bonini, M. Calandra, R. Car, C. Cavazzoni, D. Ceresoli, G. L. Chiarotti, M. Cococcioni, I. Dabo, A. D. Corso, S. de Gironcoli, S. Fabris, G. Fratesi, R. Gebauer, U. Gerstmann, C. Gougoussis, A. Kokalj, M. Lazzeri, L. Martin-Samos, N. Marzari, F. Mauri, R. Mazzarello, S. Paolini, A. Pasquarello, L. Paulatto, C. Sbraccia, S. Scandolo, G. Sclauzero, A. P. Seitsonen, A. Smogunov, P. Umari, and R. M. Wentzcovitch, “QUANTUM ESPRESSO: a modular and open-source software project for quantum simulations of materials,” J. Phys. Condens. Matter, vol. 21, p. 395502, sep 2009.


postdoctoral researcher at CISPA under the supervision of Prof. Michael
Backes. He received his Ph.D. from KAIST in 2018, advised by Prof. Brent
Byunghoon Kang and his B.S. from The University of Texas at Austin. His main
research interests lie in retrofitting security in computing systems against
today’s advanced threats. His research interests include but are not limited
to Operating System Security, Trusted Execution Environments, Program
Analysis, Software Security, and Secure AI Computation in Cloud.
---|---