Quantum theory of electronic excitation and sputtering by transmission electron microscopy †

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 under-estimate 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 too small to drive atomic dislocation by many current models. To address these theoretical shortcomings, this paper develops a ﬁ rst-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 e ﬀ ect of these electronic excitations on the sputtering cross section. Applying this method to 2D hexagonal BN and MoS 2 signi ﬁ cantly 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.


Introduction
The holy grail of materials engineering is atomic scale control of the material structure.2][3][4] These structural changes 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).Twodimensional (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 PKA 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 dose 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 ε c .5][6][7] 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 pre-dicted to be about 80 keV. 8However, sputtering has been observed in hBN under 30 keV irradiation. 9Furthermore, selenium sputters from WSe 2 and MoSe 2 under irradiation energies of 60 and 80 keV, respectively.These energies are almost 150 keV below their predicted critical energies. 10,11Lastly, while the calculated critical energy for sulfur sputtering in MoS 2 is about 90 keV, sulfur has been shown to sputter under 20 keV beams. 12Discrepancies 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. 11This would increase the sputtering cross section for all beam energies and enable sputtering for energies well-below the ground state ε c .
To account for these effects, this paper combines quantum electrodynamics (QED) and density functional theory (DFT) to derive the probability of beam-induced electronic excitation in 2D insulating crystals.7][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 organized as follows.First, we derive the theoretical model to determine beam-induced sputtering cross sections in 2D insulating crystals.We then use this model to calculate sputtering cross sections in hBN and MoS 2 that quantitatively agree with experiment.

Theoretical model
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σ/dE(ε b ,E) providing the distribution of energy transfers E for a given beam energy ε b .0][21][22] Setting ħ = c = 1, where p b is the momentum of the beam electron, β is its velocity, Z is the atomic number of the target nucleus, and α is the fine structure constant.One can then calculate the displacement cross section by integrating dσ/dE over all E large enough to cause a displacement, so that where E max (ε b ) is the maximum possible energy transfer for a given beam energy, i.e., the energy transfer resulting from a direct collision.The step function θ (which we will leave implicit going forward for the sake of compactness) enforces that the cross section is zero when E max < E d .The critical energy ε c is then the beam energy for which Therefore, the observation of sputtering at energies well below ε c is completely at odds with eqn (2).With this in mind, a fair amount of work has been done to treat deviations from eqn (2).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 eqn (2) for beam energies very close to ε c .Thus, temperature-induced increases in the cross section cannot account for the disparities between eqn (2) and experiment.This necessitates the consideration of additional phenomena that can reduce E d .
To address the limitations of eqn (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 (Fig. 1).Therefore, the rate of beam-induced sputtering hinges on the rates of three processes mediated by these interactions.
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 ε b is P i (ε b ,n i ).
2. Material electrons and nuclei: some number n f beaminduced 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 τ.The probability that n f of the n i excitations survive is P f (E,τn i n f ).
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σ/dE(ε b ,E).
The sputtering cross section can then be calculated by coupling dσ/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 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 σ in eqn (3) larger than σ 0 in eqn (2) for all beam energies, most prominently when ε b < ε c .
The remainder of this section focuses on the derivations of P i , P f , and E d .Section 1 describes how to combine QED with DFT to obtain P i .Section 2 then considers the evolution of the excited states during the sputtering process to derive E d and P f .

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 (ε b ,n i ), the probability that a beam electron with kinetic energy ε 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 eigen-state 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 (ε b ,n i ).The following subsections address each step (i-iv) in detail.
2.1.1.Scattering of free electrons.4][25] Going forward, we label the 4-momenta of the incoming electrons as p 1 and p 2 , while the outgoing electrons have momenta 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 nth electron can be written as p n = (ε n , p n x , p n y , p n z ) = (ε n , p n ), where ε n is the particle's energy and p n is its 3-momentum.Dot products between 4-vectors are then taken over Minkowski space, so that 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 (Fig. 2).Using Feynman's rules, 17,18 we can write these diagrams in terms of Dirac spinors, yielding the invariant matrix element where s n = 1 or 2 denotes the spin of the nth electron, u s ( p) is a Dirac spinor, and u ˉs( p) is its conjugate (section S1 †).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 (hundreds of eV) is much smaller than the beam energy (tens of keV), it is always the case that |p 2 | is much smaller than |p 1 |.Furthermore, we need only consider outgoing momenta for which the kinetic energy associated with either |p 3 | or |p 4 | falls within the DFT cutoff energy.In these cases, the magnitude of one outgoing momentum is similar to |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 eqn (4), it follows that one channel always contributes much more to 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 Fig. 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.
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 ≡ 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 S1. † We can use the resulting M to obtain the free electron scattering amplitude where T ˆis the scattering operator. 17,18This gives the amplitude for two free electrons with momenta p 1 and p 2 to scatter into p 3 and p 4 .Eqn ( 5) is used to derive the scattering amplitude between two arbitrary wave packets in the next subsection.
2.1.2.Scattering of wave packets.The free particle scattering amplitude in eqn (5) can be used to determine the amplitude for the scattering of two arbitrary electron states φ 1 and φ 2 into φ 3 and φ 4 .This is obtained by sandwiching the scattering operator between the initial and final 2-particle states, i.e., On the right side, we have inserted two resolutions of the identity given in eqn (S8).† Inserting eqn (5) into the integrand, the amplitude can be written in terms of the invariant matrix element M, becoming Using the delta function to integrate over p 4 and p 3 z yields where it is understood that p 4 and p 3 satisfy p 1 + p 2 = p 3 + p 4 (section S3).† The normalization of 4-momentum states |p n 〉 in terms of 3-momentum states |p n 〉 as defined in eqn (S9) † allows us to rewrite the expression as ffiffiffiffiffiffiffiffi ffi We can then discretize the momenta by replacing d 3 p i /(2π) 3 with V −1 and d 2 p i ⊥ /(2π) 2 with A −1 , where V and A are the volume and cross sectional area of the simulated crystal respectively, 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 |φ 1 〉 and |φ 2 〉 to scatter into |φ 3 〉 and |φ 4 〉 takes the form In the next subsection, we replace |φ 1…4 〉 with states relevant to electron beam-induced excitation.
2.1.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 eqn (10).We assign |φ 1  specific states into expression (10), the amplitude for exciting | nk〉 to |n′k′〉 becomes The values of ε 1⋯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 . The beam energy that appears in eqn ( 3) is then defined as the beam electron's total energy minus its rest mass: ε b ≡ ε 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., ε 2 = ε nk + m and ε 3 = ε n′k′ + m.For the remainder of this derivation, we continue to leave our expressions in terms of ε 1⋯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 MoS 2 ) 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 p 1 = |p 1 |z ˆuntil an electronic excitation is induced.
We can now evaluate the bra-ket products in eqn (11).The initial beam state is highly localized on p b , meaning that Meanwhile, the ground and excited crystal states can be expanded into a plane-wave basis, so that where each G is a reciprocal lattice vector.By re-expressing p 2 and p 3 as G 2 + k 2 and G 3 + k 3 respectively, we find Lastly, we do not care where the outgoing scattered electron ends up, so we wish for |p′ b 〉 to satisfy The excitation amplitude is then obtained by plugging in the bra-ket products from eqn ( 12), (14), and (15) into eqn (11).This gives us the excitation amplitude where it is understood that p 2 = (m + ε nk ,G 2 + k) and p 3 = (m + ε n′k′ ,G 3 + k′), 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 |nk〉 to the conduction band state |n′k′〉 for a given beam energy ε b , that is, 2.1.4.Probability of n i excitations.We are finally ready to derive P i (ε 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 where k and k′ run over all k-points, n over the valence bands, and n′ 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′ k′ ← nk).
Combinatorics tells us that the probability of exciting exactly one excitation is 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 eqn (19) are approximately equal to one another.That is, In this limit, the probability of exactly one beam-induced excitation can be written as In the same way, the probability of two excitations is In general, the probability of exactly n i beam-induced excitations is approximately Thus, we see that the probability P i can be written purely in terms of S.
We can now use formula (23) to calculate excitation probabilities in hBN and (Fig. 3).DFT is used to obtain the planewave coefficients C G 2 +k n and C G 3 +k′ n′ eigenvalues ε nk and ε n′k′ for the pristine unit cell of each material.These areMoS 2 plugged into eqn ( 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 Fig. 3 reveal some notable trends.First, for sufficiently large beam energies, P i (ε b ,n i ) decreases with increasing ε 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 (ε b ) is proportional to ε b −1 , a relationship originally predicted by Bethe. 12,26,27This means that multiple excitations are more likely at low beam energies.Furthermore, the probability of remaining in the ground state P(ε b ,n i = 0) vanishes as ε 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 ε 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 Fig. 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 Fig. 2 is the difference in energy eigenvalues between the occupied and unoccupied states (section 3).Thus, the smallest possible denominator of the t-channel in eqn ( 4) is proportional to the difference in eigenvalues squared.The experimentally measured band gap of MoS 2 is about 1.9 eV (ref.28)  while that of hBN is about 6.1 eV. 29This means MoS 2 hosts transitions with smaller eigenvalue differences, making the summands in eqn (16) larger.With P i (ε b ,n i ) derived in formula (23) and dσ/dE(ε b ,E) defined in eqn (1), we now have two of the three functions depicted in Fig. 1 needed to calculate the sputtering cross section in eqn (3).The final ingredients are P f (E,τ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.

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 eqn (3) analytically in terms of the excitation lifetime τ.Finally, in the next section, the resulting equation is used to predict the sputtering rates of boron and nitrogen in hBN and sulfur in MoS 2 , which can be made to agree well with experiment for appropriate values of τ for each material.2.2.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 eqn (3), and the differential cross section in eqn (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. 30ere instead, we make three simplifying assumptions that allow us to calculate E d with only ground state DFT (see section S7 † for a comment on the validity of these assumptions).
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,31They 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. 7Thus, 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 fully 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 for each excitation considered, 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.
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 (Fig. 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.This means that 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 trans-ferred 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. 32We 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 velocity is initially constant.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 beaminduced 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. 12e 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, 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.Eqn (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 eqn (3) diverges as E d approaches zero, meaning that the sputtering cross section σ 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 eqn (3), must have a lower bound itself, which we call E min .When E d falls below E min , we replace it with E min .We can approximate E min by setting the displacement cross section in eqn (2) equal to the area occupied by the Fig. 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.
PKA and solving for the displacement threshold (section S8 †).With this, we are finally ready to use eqn (24) to determine the set of displacement thresholds E d to insert into eqn (3).
All positive displacement thresholds computed for this work are listed in Table 1.Explainations for the calculations of E v , E a , E 0 (0), and E g are given in subsections 1 and 2. The three largest displacement thresholds for MoS 2 with n f = 0, 1, and 2 excitations are similar to those calculated with DFT-based molecular dynamics simulations. 12This 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 min in the calculation of σ.With that said, it is critical to consider large n f , even if the resulting E d is less than E min .These large n f can make appreciable contributions to the total cross section, especially at small beam energies for which P i (ε b ,n i ) is significant for large n i (Fig. 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 in the numerator.We can therefore truncate the summation over n i in eqn (3) at some adequately large n i max .This means that n i max must be converged for each material, as materials with greater S require greater n i max .In this work, we choose n i 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 (Fig. S1 †).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 τ 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 eqn (3) in terms of τ.
2.2.2.Sputtering cross section in terms of the excitation lifetime.In hBN and MoS 2 , 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 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 to be considered in calculating P f (n i n f ), the probability that n f of the n i beaminduced excitations survive long enough to reduce the displacement threshold.
We define the ratio of surviving excitations as where t s is the time it takes for the sputtered atom to travel a distance d, and the excitation lifetime τ 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 Eqn (26) allows us to rewrite the integral in eqn (3) as Assumption 3 from the previous subsection again aids us here.Because we posit that the PKA travels at a constant velocity until it reaches the sputtering distance d, we can write t s from eqn (25) quite simply as . This allows for the straightforward analytical integration of eqn (27).Using dσ/dE defined in eqn (1), the integral on the right can be written as where we define and the functions I 1 , I 2 , and I 3 are and Lastly, the function E i in eqn (31) is the exponential integral, Looking back at eqn (3), we now have everything we need to evaluate the sputtering cross section.We derived P i (ε b ,n i ) in formula (23), and we performed the integration over E analytically in eqn (27) (28), setting E d in eqn (24) as the lower integration bounds.We also have a criterion to truncate the summation over n i at n i max for a given material, as described at the end of section 1.In the following section, we use these results to calculate the sputtering cross sections of hBN and MoS 2 .

Results and discussion
We now apply our theoretical model to calculate the sputtering cross sections of boron and nitrogen from hBN and sulfur from MoS 2 .In doing so, we show that the consideration of electronic excitations is necessary to make quantitative predictions about beam-induced sputtering rates in 2D insulators.

Boron and nitrogen sputtering from hexagonal BN
Electron beam irradiation has been shown to bore nanopores in monolayer hBN at beam energies far beneath the calculated ground state critical energy of ε c ∼ 80 keV. 8,9,33,34In 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. 9,34,35etu 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 (Fig. 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 (Fig. 5).To evaluate the displacement thresholds defined in eqn (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 eqn (3) to calculate the sum of the boron and nitrogen sputtering cross sections.The excitation probabilities of hBN plotted in Fig. 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 beammatter 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, Fig. 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 inplane 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 τ = 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.
temperature effects on the cross section are considered in the manner described in our previous work. 7itting our cross section curves to the data of Cretu et al. yields an excellent agreement if the predicted excitation lifetime is set to τ ∼ 240 fs.This predicted lifetime is much shorter than the reported excitation lifetime of ∼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 postcollision.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 experiments on systems not undergoing sputtering.
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 min (Table 1).In these cases, the sputtering cross section is Ω ∼ 5 × 10 8 barn, seven orders of magnitude larger than the measured cross section at 30 keV.This suggests that one can expect the beaminduced nanopores in hBN to grow under beam energies as low as 1 keV.However, one must again be cautious of the pre-dicted cross section at low beam energies in which the treelevel theory starts to break down.Nevertheless, Fig. 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.

Sulfur sputtering from MoS 2
Kretschmer et al. measured the sputtering cross section of sulfur from MoS 2 for several beam energies ranging from 20 to 80 keV. 12They found a peak in the cross section at 30 keV, much less than the predicted ground state ε c ∼ 90 keV. 12To help explain this unexpected peak, we calculate the cross section for sulfur sputtering from pristine MoS 2 (Fig. 6).The vacant system free energy E v for the calculation of E d is that of a MoS 2 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 MoS 2 supercell.We then set E g equal to the experimental band gap of 1.88 eV. 28Using 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 τ ∼ 81 fs.8][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 MoS 2 .
Fig. 6 Sputtering of sulfur from the outgoing MoS 2 surface.The sputtering of molybdenum from MoS 2 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 τ = 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.
The difference between the two materials' lifetimes leads to markedly dissimilar cross section behaviors low beam energies.Below beam energies of 30 keV, the cross section of MoS 2 gradually drops to zero with decreasing ε b .In contrast, hBN's total cross section has a minimum at around 40 keV and begins to increase as ε b decreases.Eventually, hBN's sputtering cross section peaks before dropping quickly to zero as the beam energy goes to zero (Fig. S2 †).However, this peak occurs at around 0.5 keV, far below the lower energy bound of Fig. 5c.These distinct cross section behaviors can be explained by the amplified sensitivity of n f to τ at low beam energies.Eqn (25)  and (26) tell us that larger τ makes large n f more likely.At the same time, n f cannot exceed n i .Thus, n f is more sensitive to τ at low beam energies for which n i is large.Accordingly, because τ is greater in hBN than in MoS 2 , the expected values of n f in hBN are much larger than those of MoS 2 at low beam energies.This explains why the effect of considering excitations is much more pronounced in hBN than in MoS 2 .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 MoS 2 .This again increases the likelihood that hBN has more final excitations than MoS 2 , meaning that the effects of beam-induced excitation on hBN's sputtering cross section are greater.
We also plot the contributions of n f = 0, 1, and 2 excitations to MoS 2 '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.

Conclusions
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 beaminduced 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 MoS 2 .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 .][42] Moreover, spin polarization effects can also be examined by making M spin-dependent, i.e., not averaging over spins as is done in eqn (4).The beam electron path can also change significantly after collision with the nucleus.4][45][46] Furthermore, the methods here can be modified to accommodate more advanced DFT techniques.9][50] In addition, the plane-wave coefficients of excited electronic states can be calculated self-consistently using constrained DFT. 51Perhaps most importantly, progress in this field requires much more experimental data.We therefore hope this paper encourages new experimental investigation into beaminduced 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.

Computational details
DFT 13,14 was used to calculate (i) the plane wave coefficients C G+k n and eigenvalues ε nk of each material's electronic orbitals to determine the excitation probabilities (section 2.1.3)and the free energies of the pristine and sputtered structures to determine the displacement thresholds for sputtering (section 2.2.1).Plane-wave coefficients and eigenvalues were calculated using Quantum ESPRESSO. 52Because eqn (16) relies on the orthogonality of the Kohn-Sham orbitals, the optimized normconserving Vanderbilt pseudopotentials 53,54 were employed.The sum of all transition probabilities S defined in eqn (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% (Fig. S5 †).The cutoff energies for hBN and MoS 2 were set to 490 and 408 eV respectively.While these cutoffs are lower than what is typically used for normconserving psuedopotentials, both cutoffs yield well-converged values of S (Fig. S5b and e †).Meanwhile, the hBN and MoS 2 unit cells were given heights of 18 and 12 Å respectively.For both materials, the maximum virtual photon momentum | q max| required for convergence fell well within their first BZs.We therefore chose |q max| to be the magnitude of high-symmetry point K in each respective BZ.The maximum number of initial excitations considered for hBN and MoS 2 were n i max = 5 and 9 respectively (Fig. S1 †).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 S6 †).The most dense k-point meshes used to fit these curves had dimensions of 45 × 45 × 1 and 36 × 36 × 1 for hBN and MoS 2 respectively, corresponding to very similar linear k-point spacings of 0.056 and 0.055 Å −1 respectively.Next, free energy calculations were carried out with the Vienna ab initio Simulation Package (VASP) 15,55 implementing the projector augmented wave (PAW) method 47 along with the Perdew-Burke-Ernserhof (PBE) generalized gradient approximation (GGA) to the exchange correlation functional. 56van der Waals interactions were accounted for using the optB88-vdW density functional methods. 57,58All parameters were converged so that any increase in precision would change the total free energy by less than 1 meV per atom.The parameter values that satisfy this criteria generally differed from those listed in the previous paragraph.The cutoff energies for hBN and MoS 2 were set to 800 and 550 eV respectively, while the BZs of both materials' pristine unit cells were sampled with Γ-centered 6 × 6 × 1 Monkhorst-Pack meshes, 59 corresponding to linear k-point spacings of 0.417 and 0.328 Å −1 for hBN and MoS 2 respectively.To achieve the same k-point densities, surface vacancies in MoS 2 and edge vacancies in hBN were placed in respective 6 × 6 × 1 and 4 × 1 × 1 supercells whose BZs were sampled with a single k-point on Γ.The heights of the hBN and MoS 2 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 Å −1 .

〉 and |φ 4 〉
to the initial and final beam states |p b 〉 and |p′ b 〉 respectively.States |φ 2 〉 and |φ 3 〉 are then the ground and excited crystal states |nk〉 and |n′k′〉 respectively, where n and n′ are band indices and k and k′ are k-points.Substituting these

Fig. 2
Fig.2The 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 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.

Fig. 3
Fig. 3 Probabilities of initially exciting a certain number of electrons n i in (a) hBN and (b) MoS 2 with a beam electron.The excitation probabilities tend to decreases with increasing beam energy.The excitation probabilities in MoS 2 are larger than those in hBN because MoS 2 has a significantly smaller band gap.

Table 1
All computed positive displacement thresholds for sputtering from the hBN armchair edge and MoS 2 surface.The thresholds decrease as the number of surviving excitations n f increases