Elastic Scattering

The elastic scattering between the incident electrons and the strong nuclear Coulomb field represents the main cause of the angular deflection in the electron trajectory. Owing to the large mass of the nucleus compared to that of the electrons, the energy transfer of electrons during elastic scattering is negligible. In this study, the electron elastic total cross sections and differential cross sections were calculated based on the Mott theory^[43] using the ELSEPA code developed by Salvat et al.^{[44, 45]}. The muffin-tin radius (1.44 Å) of gold was also considered in the calculation for solid gold. In the latest version of ELSEPA^[45], this low-energy limit was extended down to 5 eV. The modification of elastic scattering was introduced by the following electron-acoustic phonon scattering.

At very low energies, electrons tend to behave like Bloch waves rather than as particles^[47]. Thus, they can interact with the collective movement of atoms in a solid, that is, phonons. We applied the theory for electron-acoustic electron scattering described by Fitting et al.^[48] and modified by Verduin^[49] and Theulings^[50]. The electron-acoustic phonon scattering rate is given as follows: $\begin{array}{c}{P}_{\text{ac}}\left(E\right)=\{\begin{array}{ll}\frac{\pi {ϵ}_{\text{ac}}^2{k}_{B}T}{\hslash c_s^2{\rho }_{\text{m}}}\frac{A}{A+E}D\left(E\right), \hfill & E<\frac{E_{BZ}}{4}\hfill \\ \frac{4\pi \left(2N_{\text{BZ}}+1\right)m_{\text{D}}{ϵ}_{\text{ac}}^2}{\hslash {\omega }_{\text{BZ}}\hslash {\rho }_{\text{m}}}\frac{A^2}{E}\left[\ln \frac{A+E}{A}-\frac{A}{A+E}\right],\hfill & E\ge E_{\text{BZ}}\hfill \end{array}\end{array},$ (1) where ${ϵ}_{\text{ac}}$ is the acoustic deformation potential, $k_{\text{B}}$ is the Boltzmann constant, $T$ is the temperature, $c_{\text{s}}$ is the sound velocity, ${\rho }_{\text{m}}$ is the mass density, and $A=5E_{BZ}$ is the screening factor adopted by Bradford and Woolf^[51]. The variable $E_{\text{BZ}}$ is the Brillouin zone energy given by $E_{\text{BZ}}=\frac{{\left(\hslash k_{\text{BZ}}\right)}^2}{2m_{\text{e}}}$, where the boundary wave vector is found at $k_{\text{BZ}}=2\pi /a$ if the first Brillouin zone is assumed to be spherical. On the other hand, $a$ is the lattice constant, and $D\left(E\right)$ is the density of states (DOS), which is typically assumed to have a parabolic dispersion relation. In this case, $D\left(E\right)$ is expressed as follows: $\begin{array}{c}D\left(E\right)=\frac{\sqrt{2m_{\text{D}}^3\left(E-E_{\text{CB}}\right)}}{{\pi }^2{\hslash }^3},\end{array}$ (2) where $m_{\text{D}}$ is the effective DOS mass and $E_{\text{CB}}$ is the energy at the bottom of the conduction band and is equal to 0 for metals.

$N_{\text{BZ}}$ is the acoustic phonon population $N\left(k,T\right)$ at the Brillouin zone boundary ($k=k_{\text{BZ}}$). The acoustic phonon population $N\left(k,T\right)$ follows the Bose-Einstein distribution. Thus, $\begin{array}{c}N\left(k,T\right)=\frac{1}{\text{exp}\left(\frac{\hslash {\omega }_{\text{ac}}\left(k\right)}{k_{\text{B}}T}\right)-1} ,\end{array}$ (3) where ${\omega }_{\text{BZ}}$ is the acoustic phonon energy ${\omega }_{\text{ac}}\left(k\right)$ at the Brillouin-zone boundary. The acoustic phonon energy ${\omega }_{\text{ac}}\left(k\right)$ follows a dispersion relation that was modified to fit the experimental data^[46] as follows: $\begin{array}{c}{\omega }_{\text{ac}}\left(k\right)=c_{\text{s}}k-\alpha k^2-\beta k^3,\end{array}$ (4) where $\alpha$ and $\beta$ are the fitting parameters. There are three acoustic phonon modes: one longitudinal and two transverse. For the longitudinal mode, the chosen parameters are $c_{\text{s}}=\text{3240m}/\text{s}$, $\alpha =\text{3}\text{.5}\times {\text{10}}^{\text{–8}}{\text{m}}^{\text{2}}/\text{s}$, and $\beta =\text{3}\text{.5}\times {\text{10}}^{\text{–18}}{\text{m}}^{\text{3}}/\text{s}$. On the other hand, the chosen parameters for the transverse modes are $c_{\text{s}}=\text{1200 m}/\text{s}$, $\alpha =\text{–7}\text{.6}\times {\text{10}}^{\text{–8}}{\text{ m}}^{\text{2}}/\text{s}$, and $\beta =\text{5}\text{.3}\times {\text{10}}^{\text{–18}}{\text{ m}}^{\text{3}}/\text{s}$. The fitted phonon dispersion relations are demonstrated in Fig. 1(b).

Subsequently, the inverse mean free path (MFP) or macroscopic cross section can be derived by dividing the scattering rate by the electron velocity, which can be expressed as $\begin{array}{c}{\lambda }_{\text{ac}}^{-1}\left(E\right)=\{\begin{array}{ll}{\lambda }_0^{-1}\frac{A}{A+E}, \hfill & E<\frac{E_{\text{BZ}}}{4}\hfill \\ {\lambda }_0^{-1}\left(N_{\text{BZ}}+\frac{1}{2}\right)\frac{8m_{\text{D}}{c}_{\text{s}}^{2}A}{\hslash {\omega }_{\text{BZ}}{k}_{\text{B}}T}\left[\frac{A}{E}\ln \frac{A+E}{A}-\frac{A}{A+E}\right],\hfill & E\ge E_{\text{BZ}}\hfill \end{array}\end{array},$ (5) where ${\lambda }_0^{-1}=\frac{\sqrt{m_{\text{e}}^{*}{m}_{\text{D}}^3}{ϵ}_{\text{ac}}^2{k}_{\text{B}}T}{\pi {\hslash }^4{c}_s^2{\rho }_m}$.

In region [$E_{\text{BZ}}/4,E_{\text{BZ}}$], the inverse MFP is interpolated between both values at the upper limit of the top formula and the lower limit of the bottom formula. According to the angular differential scattering rate, the angular differential inverse MFP is given by $\begin{array}{c}\frac{\text{d}{\lambda }_{\text{ac}}^{-1}\left(E\right)}{\text{dΩ}}=\{\begin{array}{ll}\frac{{\lambda }_0^{-1}}{4\pi }\frac{1}{{\left(1+\frac{1-\cos \theta }{2}\frac{E}{A}\right)}^2},\hfill & E<\frac{E_{\text{BZ}}}{4}\hfill \\ \frac{{\lambda }_0^{-1}}{4\pi }\left(N_{\text{BZ}}+\frac{1}{2}\right)\frac{8m_{\text{D}}{c}_s^{2}A}{\hslash {\omega }_{\text{BZ}}{k}_{\text{B}}T}\frac{\frac{1-\cos \theta }{2}\frac{E}{A}}{{\left(1+\frac{1-\cos \theta }{2}\frac{E}{A}\right)}^2},\hfill & E\ge E_{\text{BZ}}\hfill \end{array}\end{array}.$ (6)

The angular distribution of electron-acoustic scattering can also be obtained from the equation above.

Both acoustic phonon emission and absorption exist in electron-acoustic scattering; however, the probability of the former is larger than that of the latter. Therefore, the weighted average energy loss is primarily considered. Although the energy loss is only a few meV, it cannot be neglected in low-energy regions. Thus, the average energy loss $E_{\text{ac}}$ during electron-acoustic phonon scattering can be approximated by the following expression: $\begin{array}{c}{E}_{\text{ac}}=\frac{{{\displaystyle \int }}_0^{k_{\text{BZ}}}\hslash {\omega }_{\text{ac}}\left(k\right)k^2\text{d}k}{{{\displaystyle \int }}_0^{k_{\text{BZ}}}\left[2N\left(k,T\right)+1\right]k^2\text{d}k}.\end{array}$ (7)

Finally, the acoustic deformation potential of the metal was derived using the resistivity based on the model of Kieft and Bosch^[52]. The acoustic deformation potential can then be derived as follows: $\begin{array}{c}{ϵ}_{\text{ac}}=\sqrt{\frac{2\hslash e^2{c}_{\text{s}}^2}{3\pi m_{\text{e}}{k}_{\text{B}}T}{\rho }_{\text{R}}{\rho }_{\text{m}}{E}_{\text{F}}},\end{array}$ (8) where ${\rho }_{\text{R}}$ is the resistivity. Using the parameters indicated above, the acoustic deformation potential of gold was found to be 4.86 eV and 1.80 eV for the longitudinal and transverse modes, respectively.

As suggested by Fitting et al.^{[47, 52]}, electrons should behave like Bloch waves and interact with phonons at low energies (< 100 eV); however, Mott cross sections are no longer valid because the unrealistically high scattering rate leads to a MFP that is shorter than the interatomic distance. Therefore, we modified the elastic scattering as follows: Mott cross sections were applied for energies above 100 eV, whereas electron-acoustic phonon cross sections were used for energies below 50 eV. The elastic cross sections between 50 eV and 100 eV were derived through interpolation between the two cross sections. As a result, a new elastic MFP is demonstrated in Fig. 1(c). Although the elastic scattering models may not be sufficient for describing the structure of small AuNPs, a more precise atomic potential for the Mott theory will be obtained in the future using more advanced techniques, such as ab initio and DFT.

Inelastic Scattering

This study mainly considers electron-electron inelastic scattering. Inelastic cross sections were computed using the dielectric theory^[53]. The core of this theory is the complex dielectric function $\epsilon \left(q,\omega \right)$, where $q$ is the momentum loss and $\omega$ is the energy loss. The imaginary part of the inverse dielectric function is called the energy loss function (ELF) $\text{Im}\left[\frac{-1}{\epsilon \left(q,\omega \right)}\right]$, which provides detailed energy loss information for the material electrons. From the optical data of experiments^[54] or ab initio and DFT calculations^[55], the dielectric function and ELF can only be obtained at the optical limit ($q\to 0$), which is called the optical dielectric function and optical ELF.

To extend the ELF from the optical limit to $q>0$, we used the extended Drude model^[56] and Mermin model^[57]. In these two models, the sum of the Drude-type ELFs in (4) is used to fit the optical ELF because both the extended Drude model and Mermin model coincide at the optical limit^[58], that is, the Drude-type ELFs: $\begin{array}{c}\text{Im}\left[\frac{-1}{{\epsilon }_{\text{D}}\left(0,\omega \right)}\right]=\text{Im}\left[\frac{-1}{{\epsilon }_{\text{M}}\left(0,\omega \right)}\right]\\ ={\displaystyle \sum _{i=1}^n{A}_i}\frac{{\gamma }_i\omega }{{\left({\omega }^2-{\omega }_{p,i}^2\right)}^2+{\gamma }_i^2{\omega }^2},\end{array}$ (9) where $A_i=a_i{\omega }_{p,i}^2$ is the oscillator intensity, ${\gamma }_i$ is the damping constant, and ${\omega }_{p,i}$ is the plasma frequency.

The fitting model was built based on several previous studies^{[39, 40, 59, 60]}. In the development of Geant4 MircoElec, Valentin^[60] and Gibaru^[40] et al. used a Heaviside step function to cut each oscillator in the fitting. Pablo et al.^[61] used an exponential smooth switching function instead of a simple step function for the outer-shell electrons. Based on these prior methodologies, the fitting model is expressed as follows: $\begin{array}{c}\text{Im}\left[\frac{-1}{\epsilon \left(0,\omega \right)}\right]={\displaystyle \sum _{i=1}^{n}F\left(\omega \right)A_i}\frac{{\gamma }_i\omega }{{\left({\omega }^2-{\omega }_{p,i}^2\right)}^2+{\gamma }_i^2{\omega }^2}.\end{array}$ (10)

For the outer-shell electrons, $F\left(\omega \right)$ represents a smooth exponential switching function for the electronic excitations^[61]: $\begin{array}{c}F\left(\omega \right)=\frac{1}{1+e^{-{\text{Δ}}_i\left(\omega -{\omega }_{th.i}\right)}},\end{array}$ (11) where ${\text{Δ}}_i$ and ${\omega }_{\text{th}.i}$ are the additional fitting parameters. For the inner shells, particularly the K and L shells, which have sharp edges in the optical ELF, $F\left(\omega \right)$ is chosen as a simple step function $\text{Θ}\left(\omega -{\omega }_{th,i}\right)$. For plasmons and conduction band electrons, $F\left(\omega \right)$ can be simply set to 1.

The experimental optical ELF, which covers the electron energies from 0.005 eV to 1 MeV, is derived from the databases of Palik^[62], Windt^[63], and the evaluated photo data library 1997 version (EPDL97^[64]), in which the high-energy optical ELF is calculated using the photo-ionization cross section^[65] as follows: $\begin{array}{c}\text{Im}\left[\frac{-1}{\epsilon \left(0,\omega \right)}\right]=\frac{n_{\text{m}}c{\sigma }_{\text{phot}}}{\omega },\end{array}$ (12) where $n_{\text{m}}$ is the atomic density, $c$ is the speed of light, and ${\sigma }_{\text{phot}}$ is the photoionization cross section.

The fitting results were validated by using two sum rules (f-sum and ps-sum) and the mean ionization potential. A partial f-sum rule for the fitted oscillators was also used to distinguish the contribution of the oscillators to the electronic shells. Initial values of the fitting (energy loss from 0 to 100 eV) were partially obtained by referring to the study of Pablo et al.^[61]. We managed to use 28 oscillators to describe the optical ELF. Figure 2 demonstrates that the fitted ELF data sufficiently agrees with the experimental ELF.

Fig. 2Full-screen View

(Color online) Comparison of the fitted ELF (solid) and experimental ELF (circles): (a) contribution of each electronic shell obtained from the combined oscillators at the energy loss range of 0–100 eV; (b) contribution to all electronic shells of the combined oscillators.

During fitting, the plasmon peak was chosen as 25.42 eV, which is near 24.8 eV in Egerton's book^[66]. We verified the f-sum rule of the fitted ELF, that is, $Z_{\text{eff}}=78.05$ (ideal value: 79.0), while the ps-sum rule resulted in 1.0908 (ideal value: 1.0). The mean ionization potential based on the fitted ELF was 736.45 eV. These values obtained from the fitted ELF coincide with those from the experimental ELF, with only a small deviation in the high-energy region (L and K shells). This is likely due to the restriction of the Drude functions. If a more accurate experimental ELF is available, then better fitting results can be obtained. We also checked the partial f-sum rule to determine the shell to which the fitted oscillators belong, and the results are also displayed in Fig. 2. Therefore, these categorized oscillators enabled the calculation of the inelastic mean free path (IMFP) of the different shells.

The calculation of IMFP is achieved by extrapolating the momentum loss q from the optical limit to $q>0$. For the extended Drude model, the extrapolation of q was achieved by using the dispersion relation. The extrapolation was performed using the following equation: $\begin{array}{c}\text{Im}\left[\frac{-1}{{\epsilon }_{\text{D}}\left(q,\omega \right)}\right]={\displaystyle \sum _{i=1}^n{A}_i}\frac{{\gamma }_i\left(q\right)\omega }{{\left({\omega }^2-{\omega }_{p,i}^2\left(q\right)\right)}^2+{\gamma }_i^2\left(q\right){\omega }^2}.\end{array}$ (13)

In this method, ${\omega }_{p,i}\left(q\right)$ describes the plasmon dispersion and ${\gamma }_i\left(q\right)$ is the q-dependent damping constant, which represents the broadening effect. The variable ${\gamma }_i$ is usually considered a constant; thus, one can simply choose (1) ${\omega }_{p,i}\left(q\right)={\omega }_{p,i}+\frac{q^2}{2}$ or (2) ${\omega }_{p,i}\left(q\right)=\sqrt{{\omega }_{p,i}^2+\frac{6E_{\text{F}}}{5}{q}^2+\frac{q^2}{4}}$ ^{[56, 67]}. If ${\gamma }_i$ is considered a function of $q$, another dispersion relation is given by^{[68, 69]} (3) ${\omega }_{p,i}\left(q\right)=\sqrt{{\omega }_{p,i}^2+\frac{2E_{\text{F}}}{3}{q}^2+\frac{q^2}{4}}$ and ${\gamma }_i\left(q\right)=\sqrt{{\gamma }_i^2+\frac{q^2}{4}}$. This was already presented in the study of Ritchie et al.^[56] and proved^[68] to be more comprehensive than the extrapolation considering only the ${\omega }_{p,i}\left(q\right)$. Therefore, we chose the third extrapolation method for the conduction band and outer shells and the first method for the inner shells as the first inelastic model, which is labeled as Ritchie in the following text.

Regarding the Mermin model, the extrapolation procedure was already achieved using the Mermin dielectric function^[70] as follows: $\begin{array}{c}{\epsilon }_{\text{M}}\left(q,\omega \right)=1+\frac{\left(1+\frac{i\gamma }{\omega }\right)\left[{\epsilon }_{\text{L}}\left(q,\omega +i\gamma \right)-1\right]}{1+\frac{i\gamma }{\omega }\left[{\epsilon }_{\text{L}}\left(q,\omega +i\gamma \right)-1\right]/\left({\epsilon }_{\text{L}}\left(q,0\right)-1\right)},\end{array}$ (14) where ${\epsilon }_{\text{L}}$ is the Lindhard dielectric function^[71] and ${\epsilon }_{\text{L}}\left(q,0\right)$ is the Lindhard dielectric function at the static limit. This is the second inelastic model, which is also labeled as Mermin.

A significant advantage of the extended Drude model and Mermin model is that fitted oscillators can approximate the contribution to particular electronic shells or the plasmon in the calculation. Therefore, combined with the divided shells in the ELF, we can calculate the IMFP of the different electronic shells.

The general calculation of IMFP considering the relativistic high energy correction is performed by $\begin{array}{c}{\text{λ}}_{\text{inel}}^{-1}\left(E\right)=\frac{{\left(1+\frac{E}{c^2}\right)}^2}{1+\frac{E}{2c^2}}\frac{1}{\pi E}\underset{W_{-}}{\overset{W_{+}}{{\displaystyle \int }}}\text{d}\omega \underset{q_{-}}{\overset{q_{+}}{{\displaystyle \int }}}{f}_{\text{ex}}\left(q,E\right)\frac{\text{d}q}{q}F\left(\omega \right)\text{Im}\left[\frac{-1}{\epsilon \left(q,\omega \right)}\right],\end{array}$ (15) where E is the kinetic energy of the incident electrons, c is the speed of light, $\frac{{\left(1+\frac{E}{c^2}\right)}^2}{1+\frac{E}{2c^2}}$ is the relativistic correction, $W_{+}$ is the maximum energy loss, $W_{-}$ is the minimum energy loss, and $q_{+}$ and $q_{-}$ are the maximum and minimum momentum losses, respectively. The variables $q_{+}$ and $q_{-}$ are given by $\begin{array}{c}{q}_{\pm }=\sqrt{E\left(2+\frac{E}{c^2}\right)}\pm \sqrt{\left(E-\omega \right)\left(2+\frac{E-\omega }{c^2}\right)},\end{array}$ (16) where $f_{\text{ex}}\left(q,E\right)=1-\frac{q^2}{2E}+{\left(\frac{q^2}{2E}\right)}^2$ is the Ochkur exchange correction^[72]. For a collective excitation, $f_{\text{ex}}=1$. As indicated above, $F\left(\omega \right)$ is an exponential switching function or a simple step function.

The inelastic cross sections of each shell were individually calculated. For metals such as Au, $W_{-}$ equals zero. As proposed by Pablo et al.^[61], the maximum energy loss $W_{+}$ depends on the different electronic shells. Plasmon and 5d electrons are considered collective excitations; thus, $W_{+}=E-E_{\text{F}}$. For the outer-shell and inner-shell electrons, $W_{+}$ equals $\text{min}\left[E-E_{\text{F}},\left(E+B_i\right)/2\right]$, where $B_i$ is the binding energy of the i-th electronic shell. The binding energies of the outer- and inner-shell electrons were derived from the FFAST database^[73]. Meanwhile, a classical Coulomb-field Born correction that accounts for the potential energy gained by the incident electrons is added for the calculation of the outer-shell and inner-shell electrons. This is achieved by replacing the kinetic energy $E$ with $E^{\prime }=E+2B_i$ in the calculation; however, the energy loss should remain unchanged.

Based on the aforementioned discussion, we compared calculated IMFP with previously reported experimental^[74-76] and theoretical data^[77-86], which are shown in Fig. 3.

Fig. 3Full-screen View

(Color online) (a) Calculated inelastic mean free path and comparison with other published data^[74-76]. (b) Calculated inner shell cross sections (K, L, M shells); solid, dashed, and dotted lines indicate the Ritchie model, Mermin model, and Gryzinski theory, respectively^{[87, 88]}.

Figure 3 (a) presents the calculated IMFP, which sufficiently agrees with the experimental and previous theoretical data. A relativistic effect is clearly observed at high energies. At low energies below 100 eV, the modified Ritchie model can derive a good IMFP similar to that of the Mermin model due to its improved estimation of the plasmon lifetime by considering $\gamma$ as a function of $q$. We also calculated the inner-shell ionization cross sections and compared them with the relativistic Gryzinski's formula^{[87, 88]} for validation. As shown in Fig. 3 (b), the calculation using the Ritchie model can yield better inner-shell ionization data compared to the Mermin model, which indicates that the Ritchie model is recommended for high-energy electron simulations. This may be due to the fact that for the inner-shells, the plasmonic consideration in the Mermin model is no longer valid, whereas in the Ritchie model, the simple extrapolation of $q$ without considering the plasmon is more applicable. Although the inner shell cross sections using both the Ritchie and Mermin models remain lower than those calculated by the Gryzinski theory, which is considered more accurate, this may be due to the fact that the fitted ELF is insufficiently accurate for the inner shells. Therefore, this lowers the accuracy of the calculated inner-shell cross sections.

Following an inelastic event, in normal considerations, a secondary electron will be generated with an energy of $E_{\text{S}}=W+E_{\text{F}}$ for Fermi electron excitation (conduction band and plasmon excitation as considered in this study) and $E_{\text{S}}=W-B_i$ for the inner- and outer-shell ionization, where W is the energy loss. In this case, we applied the improvement made by Ding et al.^[89], which assumes that the initial energy $E^{\prime }$ of the Fermi electron is proportional to the joint density of states (JDOS) of the free-electron gas, that is, $\text{P}\left(E^{\prime },W\right)\propto \sqrt{E^{\prime }\left(E^{\prime }+W\right)}$. In the expression, $E^{\prime }$ must satisfy $E^{\prime }<E_{\text{F}}$ and $E^{\prime }+\text{W}>E_{\text{F}}$. Therefore, the modified energy of the Fermi electron after inelastic scattering is $E_{\text{S}}=\text{W}+E^{\prime }$.

Surface Potential

For an electron to escape from the surface of metals, it must have a kinetic energy larger than the Fermi level E_F and the work function $W_{\text{F}}$, that is, the potential barrier between the metal surface and vacuum^[89]. Therefore, the cut-off energy in the simulation is set to $E_{\text{F}}+W_{\text{F}}$, and the electrons are locally terminated below this limit, which can also save memory.

To escape from the surface, once an electron reaches the gold surface, its kinetic energy must satisfy the following condition: $\begin{array}{c}E{\text{cos}}^2\theta >U,\end{array}$ (17) where $E$ is the energy that considers the Fermi level as the zero point, $\theta$ is the angle between the electron motion direction and is normal to the surface before its emission, and $U=W_{\text{F}}$ is the potential barrier above the Fermi level. Subsequently, according to the quantum theory, the emission will depend on the transmission coefficient T^[89-91], which is given by the following: $\begin{array}{c}T=\frac{4\sqrt{1-U/\left(E{\text{cos}}^2\theta \right)}}{{\left[1+\sqrt{1-U/\left(E{\text{cos}}^2\theta \right)}\right]}^2}.\end{array}$ (18)

A random number $R$ is generated for comparison with $T$ to determine the emission as follows: $\{\begin{array}{ll}\text{Emission},\hfill & \text{if} R<T\hfill \\ \text{Reflection},\hfill & \text{else}\hfill \end{array}.$ (19)

The emitted electrons have a kinetic energy of $E-U$, and electrons entering the solid will also gain potential energy $W_{\text{F}}$ above the Fermi level. In the SEE simulation, the work function between the gold surface and vacuum was set to 5.1 eV. Considering practical cases of water-gold interfaces for future development, we will use the difference between the work functions of these two materials and assume that the calculation of the transmission coefficient is similar to that considered in the surface-vacuum interface.

Secondary electron emission

3.1.1

Secondary electron yield

To validate the code, we simulated SEE experiments and calculated the secondary electron yield (SEY) and backscattering coefficient (BSC) of gold. This simulation was simply set up as multiple groups of mono-energetic electron beams vertically bombarding an infinite gold plane, and the information of the electrons emitted from the surface was collected. The SEY is defined as the ratio of the number of emitted secondary electrons (< 50 eV) to the number of incident primary electrons. The BSC is defined as the ratio of the number of backscattering electrons (> 50 eV) to the number of incident primary electrons. To obtain the SEY and BSC, 1×10⁶ particles were simulated. The uncertainty of this code in the SEE simulation is less than 1%. As presented in Fig. 4, we demonstrate the calculated SEY and BSC and compare them with the experimental data collected by Joy et al.^[92]. Regarding the experimental data, monoenergetic electron beams were used in most of the experiments to bombard the material under vacuum conditions, which was similar to the simulation setup.

Fig. 4Full-screen View

(Color online) (a) SEY and BSC derived from the SEE simulations and compared to the experimental data^[92]; M indicates the Mermin model and R indicates the Ritchie model. (b) Comparison of the SEY values calculated by two ways (solid lines: Mott-eph, dashed lines: Mott).

As shown in Fig. 4 (a), the simulated SEY values mostly coincide with the experimental values^[92], which is essential for the calculation of the emission of LEEs from AuNPs. The Ritchie model derived a larger SEY compared to the Mermin model, whereas the BSC values were similar in both models. In addition, if the modification for low-energy elastic scattering was not added, the SEY values calculated in both ways would mainly differ below 1 keV, as shown in Fig. 4 (b). The figure shows that the calculated SEY using only the Mott cross sections is larger than the experimental values for the primary electrons with a lower energy (< 400 eV). The modification of the elastic cross sections improved the description of the SEE. The high SEY also suggests that the ¹²⁵I emitted electrons can enhance the LEE emission from gold nanoparticles.

3.1.2

Secondary electron spectra

The secondary electron spectra are critical for the analysis of LEEs emitted from AuNPs because the cross sections of the LEEs damaging the DNA, such as DSBs and SSBs, are highly energy-dependent^[9]. The 1×10⁷ electrons were simulated to obtain stable secondary electron spectra, which are presented in Fig. 5. The figure demonstrates that most of the emitted secondary electrons were below 30 eV; a zero-loss elastic peak is also observed herein.

Fig. 5Full-screen View

(Color online) Emitted electron spectra from SEE simulations: (a) 200 eV primary electrons, (b) 500 eV primary electrons, (c) 1 keV primary electrons, and (d) 2 keV primary electrons.

Moreover, the peak of LEEs was found to be located at 2.25–3.25 eV, which is similar to the 2.2 eV peak value reported in the ion-bombardment SEE experiment^[93]. This large amount of LEE emission also suggests the feasibility of using sub-keV electrons combined with AuNPs as a high-flux LEE source.

Dosimetry test

Dosimetry testing is critical for the evaluation of further microdosimetry. The stopping power, $\text{S}\left(E\right)$, for electrons of energy E can be calculated using a method similar to the IMFP. The expression for $S\left(E\right)$ is shown below: $\begin{array}{c}S\left(E\right)=\frac{{\left(1+\frac{E}{c^2}\right)}^2}{1+\frac{E}{2c^2}}\frac{1}{\pi E}\underset{W_{-}}{\overset{W_{+}}{{\displaystyle \int }}}\omega \text{d}\omega \underset{q_{-}}{\overset{q_{+}}{{\displaystyle \int }}}{f}_{\text{ex}}\left(q,E\right)\frac{\text{d}q}{q}F\left(\omega \right)\text{Im}\left[\frac{-1}{\epsilon \left(q,\omega \right)}\right].\end{array}$ (20)

As shown in Fig. 6, we compared the calculated stopping power with those of other theoretical calculations^{[81, 85, 94-97]} and experiments^[98]. The effect of the relativistic corrections introduced in the calculation is demonstrated by the stopping powers at high electron energies, which is in good agreement with the ESTAR^[85] and the evaluated electron data library (EEDL) databases^[81]. With the modifications, the Ritchie model can produce an accurate stopping power comparable to that of the Mermin model at low energies.

Fig. 6Full-screen View

(Color online) Calculated stopping power and comparison with other published data^{[81, 85, 94-97]}.

The energy deposition was mainly tested using the semi-infinite gold model. Calculation of the energy deposition depth profiles can provide information regarding the energy deposited in gold and its uniformity. This verification also guarantees the correctness of dosimetry for future studies. Electron energy deposition as a function of the depth of semi-infinite gold was calculated for comparison using the developed code and Geant4 standard electromagnetic physics (option 4). In the calculation, 200 eV, 500 eV, 1 keV, 2 keV, and 10 keV electrons were simulated; 5e6 electrons were adopted in each test. The calculated results are shown in Fig. 7.

Fig. 7Full-screen View

(Color online) (a) Calculated total energy deposition. (b) Calculated energy deposition as a function of depth (1 – 1000 Å with 1 Å as the thickness of each layer) using the developed code (solid and dashed lines indicate the Mermin and Ritchie models, respectively) and Geant4 (dotted line).

The outcomes of both inelastic models in the code present a similar tendency. Figure 7 (a) demonstrates that the ratio of the total energy deposition results of this code to that of Geant4 was approximately 89%–93%. This slight decrease was due to the larger secondary electron emission from the surface. In contrast, the difference in the energy deposition depth profiles is shown in Fig. 7(b). At low energies, our inelastic cross sections differ from those of the Livermore model (EEDL^[81]). Gibaru et al.^[40] suggested that the Geant4 standard physics underestimates the range of electrons, which is also observed upon comparing our results. At high energies, the energy deposition depth profiles of both this code and Geant4 tended to coincide. Note that the more detailed energy-straggling strategy in this code is significantly slower than the multiple scattering strategy in Geant4 standard physics.

A function for calculating the energy deposition in the 3D voxels was also added to enable future microdosimetry studies. The method for calculating the energy deposition in the 3D voxels is demonstrated, and the calculation is displayed in Fig. 8. The length of the 3D voxels was set to 1 Å, and the maximum length of each axis was 1000 Å. The 1×10⁷ electrons were simulated here.

Fig. 8Full-screen View

(Color online) (a) Diagram of the calculation of the energy deposition in 3d voxels; projection of the calculated energy deposition (Ritchie) by the 3D voxels. (b) 1 keV X-Z plane, (c) 2 keV X-Z plane, (d) 10 keV X-Z plane, (e) 1 keV X-Y plane, (f) 2 keV X-Y plane, and (g) 10 keV X-Y plane.

In Fig. 8, the spatial distribution of the energy deposition demonstrates good symmetry, which verifies the particle transport in the Monte Carlo procedure. The intensity of the energy deposition in the X-Z plane is more comprehensive than that of the depth profiles. At low energies, the projection of energy deposition in the X-Z plane is an ellipse around the entry point of the electron beams. As the incident electron energy increases, the energy deposition distribution tends to become more cylindrical. The distribution of the energy deposition also demonstrates the nanoscale precision of this code. This function allowed the further visualization and analysis of the electron energy deposition inside the AuNPs when irradiated with electrons.