Monte Carlo modeling of low-energy electron-induced secondary electron emission yields in micro-architected boron nitride surfaces

TL;DR

This study uses Monte Carlo simulations to analyze how micro-porous boron nitride surfaces can significantly reduce secondary electron emission yields, improving the durability of Hall-effect thruster walls.

Contribution

It introduces a novel stochastic modeling approach combining multi-scattering Monte Carlo simulations with ray-tracing to evaluate SEE yields on micro-architected BN surfaces.

Findings

Micro-porous BN surfaces reduce SEE yields by over 50%.

The model accurately predicts SEE behavior across relevant energy ranges.

Micro-architected surfaces are promising for mitigating SEE in electric propulsion.

Abstract

Surface erosion and secondary electron emission (SEE) have been identified as the most critical life-limiting factors in channel walls of Hall-effect thrusters for space propulsion. Recent wall concepts based on micro-architected surfaces have been proposed to mitigate surface erosion and SEE. The idea behind these designs is to take advantage of very-high surface-to-volume ratios to reduce SEE and ion erosion by internal trapping and redeposition. This has resulted in renewed interest to study electron-electron processes in relevant thruster wall materials. In this work, we present calculations of SEE yields in micro-porous hexagonal BN surfaces using stochastic simulations of electron-material interactions in discretized surface geometries. Our model consists of two complementary parts. First we study SEE as a function of primary electron energy and incidence angle in flat surfaces using Monte Carlo simulations of electron multi-scattering processes. The results are then used to represent the response function of discrete surface elements to individual electron rays generated using a ray-tracing Monte Carlo model. We find that micro-porous surfaces result in SEE yield reductions of over 50% in the energy range experienced in Hall thrusters. This points to the suitability of these micro-architected surface concepts to mitigate SEE-related issues in compact electric propulsion devices.

Tables

Table 1. Table 1 : Parameters used to model the energy loss function of h ℎ h -BN.

$h$-BN	$\hslash {\omega }_{0i}$ $\left[\mathrm{eV}\right]$	$A_{0i}$ $\left[{\mathrm{eV}}^2\right]$	${\gamma }_{0i}$ $\left[\mathrm{eV}\right]$
($\zeta$ = 0.05)
1	8.65	6.6	0.5
2	19.0	18.0	5.0
3	25.2	270.0	9.5
4	35.8	140.0	10.0
5	51.0	80.0	20.0
6	65.0	20.0	20.0

Table 2. Table 2 : Parameters for h ℎ h -BN in present calculation. Among these, α c subscript 𝛼 𝑐 \alpha_{c} , W p ​ h subscript 𝑊 𝑝 ℎ W_{ph} , C 𝐶 C and η 𝜂 \eta are free parameters.

	Parameters
$h$-BN	${\alpha }_c$	$E_g$ [eV]	$\epsilon \left(\infty \right)$	$\epsilon \left(0\right)$	$W_{ph}$ [eV]	$C\left[{\mathrm{nm}}^{-1}\right]$	$\eta \left[{\mathrm{eV}}^{-1}\right]$	$\chi$ [eV]
	0.5	5.2	4.10	5.09	0.1	1.0	0.1	4.5

Equations

$$\frac{d{\sigma }_{el}}{d\Omega }=|f(\theta ){|}^2+|g(\theta ){|}^2$$

$$f(\theta )$$

$$+l[\exp (2i{\delta }_l^{-})-1]\}{P}_l(\cos \theta )$$

$$g(\theta )=\frac{1}{2iK}\sum _{l=1}^{\infty }[\exp (2i{\delta }_l^{-})-\exp (2i{\delta }_l^{+})]P_l^1(\cos \theta ).$$

$$P_l^1(x)={\left(1-x^2\right)}^{\frac{1}{2}}\frac{d{P}_l(x)}{dx}.$$

$$\tan {\delta }_l^{\pm }=\frac{K{j}_{l+1}(Kr)-j_l(Kr)[\xi \tan {\varphi }_l^{\pm }+(1+l+k^{\pm })/r]}{K{n}_{l+1}(Kr)-n_l(Kr)[\xi \tan {\varphi }_l^{\pm }+(1+l+k^{\pm })/r]}$$

$$\xi =\frac{E+m{c}^2}{\hslash c}.$$

$${\varphi }_l^{\pm }=\underset{r\to \infty }{\lim }{\varphi }_l^{\pm }(r)$$

$$\frac{d{\varphi }_l^{\pm }(r)}{dr}=\frac{k^{\pm }}{r}\sin [2{\varphi }_l^{\pm }(r)]-\frac{m{c}^2}{\hslash c}\cos [2{\varphi }_l^{\pm }(r)]+\frac{E-V(r)}{\hslash c}$$

$$\bigg{(}\frac{d\sigma_{el}(E,\vartheta)}{d\Omega}\bigg{)}_{\rm BN}=\bigg{(}\frac{d\sigma_{el}(E,\vartheta)}{d\Omega}\bigg{)}_{\rm B}+\bigg{(}\frac{d\sigma_{el}(E,\vartheta)}{d\Omega}\bigg{)}_{\rm N}.$$

$${\sigma }_{el}(E)=\frac{{\lambda }_{el}^{-1}}{N}=2\pi {\int }_0^{\pi }\frac{d{\sigma }_{el}(E,\vartheta )}{d\Omega }\sin \vartheta d\vartheta$$

$$R_c(E)=\tan [{\alpha }_c(E/E_g{)}^2]$$

$$R=P_{el}(E,\theta )=\frac{1}{{\sigma }_{el}}{\int }_0^{\theta }\frac{d{\sigma }_{el}(E,\vartheta )}{d\Omega }\sin \vartheta d\vartheta .$$

$$\operatorname{Im}\bigg{[}-\frac{1}{\varepsilon(\textbf{q},\Delta E)}\bigg{]}=\sum_{i=1}^{n}\frac{A_{i}\gamma_{i}\Delta E}{((\hslash\omega_{0i\textbf{q}})^{2}-\Delta E^{2})^{2}+\gamma_{i}^{2}\Delta E^{2}}\times\theta(\Delta E-E_{g}),$$

$$\hslash {\omega }_{0i\text{q}}=\hslash {\omega }_{0i}+\zeta \frac{{\hslash }^2{\text{q}}^2}{2m}$$

$$\frac{d\sigma_{inel}(E,\Delta E)}{d\Delta E}=\frac{me^{2}}{2\pi\hslash^{2}NE}\operatorname{Im}\bigg{[}-\frac{1}{\varepsilon(0,\Delta E)}\bigg{]}S\Big{(}\frac{\Delta E}{E}\Big{)}$$

$$S(x)=(1-x)\ln \frac{4}{x}-\frac{7}{4}x+x^{3/2}-\frac{33}{32}{x}^2$$

$${\sigma }_{inel}(E)$$

$$\displaystyle=\frac{\lambda_{inel}^{-1}}{N}=\frac{me^{2}}{2\pi\hslash^{2}NE}\int_{W_{min}}^{W_{max}}\operatorname{Im}\bigg{[}-\frac{1}{\varepsilon(0,\Delta E)}\bigg{]}S\Big{(}\frac{\Delta E}{E}\Big{)}d\Delta E$$

$$\lambda_{inel}^{-1}(E)=\frac{me^{2}}{2\pi\hslash^{2}E}\int_{W_{min}}^{W_{max}}\operatorname{Im}\bigg{[}-\frac{1}{\varepsilon(0,\Delta E)}\bigg{]}S\Big{(}\frac{\Delta E}{E}\Big{)}d\Delta E$$

$$-\frac{dE}{ds}=\frac{me^{2}}{\pi\hslash^{2}E}\int_{0}^{W_{max}}\operatorname{Im}\bigg{[}-\frac{1}{\varepsilon(0,\Delta E)}\bigg{]}G\Big{(}\frac{\Delta E}{E}\Big{)}\Delta Ed\Delta E$$

$$G(x)=\ln \frac{1.166}{x}-\frac{3}{4}x-\frac{x}{4}\ln \frac{4}{x}+\frac{1}{2}{x}^{3/2}-\frac{x^2}{16}\ln \frac{4}{x}-\frac{31}{48}{x}^2.$$

$$R=P_{inel}(W,E)=\frac{1}{{\sigma }_{inel}}{\int }_0^W\frac{d{\sigma }_{inel}}{d\Delta E}d\Delta E$$

$$\lambda_{ph}^{-1}=\frac{1}{a_{0}}\bigg{[}\frac{n(T)+1}{2}\bigg{]}\bigg{[}\frac{\varepsilon(0)-\varepsilon(\infty)}{\varepsilon(0)\varepsilon(\infty)})\bigg{]}\frac{\hslash\omega}{E}\ln\bigg{\{}\frac{[1+\sqrt{1-\hslash\omega/E}]}{[1-\sqrt{1-\hslash\omega/E}]}\bigg{\}}$$

$$\displaystyle\cos\theta=\bigg{(}\frac{E+E^{\prime}}{2\sqrt{EE^{\prime}}}\bigg{)}(1-B^{R})+B^{R}$$

$$B=\frac{E+E^{\prime }+2\sqrt{E{E}^{\prime }}}{E+E^{\prime }-2\sqrt{E{E}^{\prime }}}$$

$${\lambda }_{pol}^{-1}(E)=C{e}^{-\eta E}$$

$$\Delta s=-{\lambda }_T\ln R,$$

$$\frac{1}{{\lambda }_T}=\frac{1}{{\lambda }_{el}}+\frac{1}{{\lambda }_{inel}}+\frac{1}{{\lambda }_{ph}}+\frac{1}{{\lambda }_{pol}}.$$

$$0<R\le \frac{1/{\lambda }_{el}}{1/{\lambda }_T}⟹\text{elastic scattering}$$

$$\frac{1/{\lambda }_{el}}{1/{\lambda }_T}<R\le \frac{1/{\lambda }_{el}+1/{\lambda }_{inel}}{1/{\lambda }_T}⟹\text{inelastic scattering}$$

$$\frac{1/{\lambda }_{el}+1/{\lambda }_{inel}}{1/{\lambda }_T}<R\le \frac{1/{\lambda }_{el}+1/{\lambda }_{inel}+1/{\lambda }_{ph}}{1/{\lambda }_T}⟹\text{phonon excitation}$$

$$\displaystyle\frac{1/\lambda_{el}+1/\lambda_{inel}+1/\lambda_{ph}}{1/\lambda_{T}}<R\leq 1\Big{(}\equiv\frac{1/\lambda_{el}+1/\lambda_{inel}+1/\lambda_{ph}+1/\lambda_{pol}}{1/\lambda_{T}}\Big{)}$$

Full text

1 Introduction

Advances in electrode, chamber, and structural material technology will enable breakthroughs in future generations of electric propulsion and pulsed power (EP $\&$ PP) concepts [1, 2]. Although significant advances have been achieved during the past few decades, much of the progress has relied on empirical development of materials through experimentation and trial-and-error approaches [3, 4]. Materials and channel wall designs are sought by optimizing performance against weight, power density, and cost. Under extreme operating environments, the discharge channels and cathode insulators used in Hall thrusters require a very demanding set of properties. Typically, the range of materials of interest in EP $\&$ PP include refractory metals, such as W, Mo, and their alloys, ceramic composites, such as BN and $\rm Al_{2}O_{3}$, high-strength copper alloys, and carbon-carbon composites [5]. These classes of materials possess great mechanical strength, thermal shock resistance, refractoriness and machinability. Surface erosion and secondary electron emission (SEE) from the channel walls have been identified as among the most critical life-limiting factors for Hall thrusters [6]. To mitigate surface erosion and SEE, a promising route is to modify the surface roughness. Recently, demonstration designs have been developed, including surface architectures based on metal micro-spears, micro-nodules and micro-velvets [7, 8, 9, 10, 11]. See reviews on the topic for further information [12, 13, 14, 15, 16, 5].

The main objective of our work is to develop reliable physical models of SEE, to be applied to the calculation of effective SEE yields in micro-architected surfaces for space propulsion thrusters. The methodology relies on two distinct but complementary elements. The first is an experimentally validated theoretical model of electron scattering in solids. This model is built as a transport Monte Carlo simulator of individual electron trajectories in a solid, capturing the pertinent scattering mechanisms in terms of an interaction differential cross section that is integrated across the relevant energy and angular ranges. These cross sections reflect different elastic and inelastic (e.g. core electrons, valence electrons, polarons, plasmons, etc) scattering mechanisms in each material, and is formulated according to the best available physics. The second element of the methodology is a discretization procedure to represent arbitrary surface geometries in terms of discrete boundary elements. Both modules (physics and geometry) are coupled by way of a raytracing algorithm that captures the intersection of primary electron rays with different boundary elements. The material considered here is BN (used in current thrusters thanks to its low mass density, low thermal expansion, and highly dielectric properties). The procedure is described in two papers from our group focused on W [17, 18].

The structure of the paper is as follows. First, the theoretical models employed to study electron interaction with insulators is described. Subsequently, these models are validated experimentally for SiO2 and hexagonal BN. The resulting SEE yield and emitted energy distributions for BN are then used as response functions in ray-tracing Monte Carlo simulations of SEE in discretized complex foams. We finalize with a summary of the main findings and the conclusions and acknowledgments.

2 Theory and Methods

2.1 Electron-Insulator Interaction Model

The energy loss mechanisms for internal secondary electrons differ between metals and insulators. In metals, the internal secondary electrons lose energy through interactions with conduction electrons, lattice vibrations –known as polarons–, and defects. In order to escape, the kinetic energy of an internal secondary electron must be above the energy barrier –taken to be equal to the Fermi level $E_{F}$ plus work function $\Phi$ (typically $>$10 eV)–when it reaches the surface. This large threshold escape energy, as well as the high collision probability due to the large number of conduction electrons, result in the low SEE yields (usually $<$1) found in metals. Conversely, due to the low number of conduction electrons in insulators, internal secondary electrons lose energy primarily through the excitation of valence electrons into the conduction band. This prevents secondary electrons with kinetic energies below the bandgap energy from participating in such electron-electron collisions, significantly increasing their mean escape depths compared to that of metals [19]. The mean escape depth for insulators ranges from 10-50 nm, compared to 0.5-1.5 nm for conductors [20]. Consequently, SEE yields in insulators are typically substantially higher than in metals. In some reported cases, the SEE yields of certain insulators can exceed those of standard conductors by as much as a factor of 20 [20].

When the energy of the primary electron beam is considerably higher than the bandgap energy $E_{g}$ of an insulator material, the elementary scattering processes are essentially those encountered for metals. However, at low energies, several new aspects of the electron-material interaction processes become important. For example, at a low-energy electron within an insulator can locally distort the lattice and greatly reduce its mobility by forming a polaron. As well, defects and impurities can also act as traps. In any case, a distinctive feature of secondary electron emission in insulators is the buildup of charge, such that subsequent SEE must be considered in the context of the existence of an internal electric field. In the present work, our calculations consider the following electron-material interaction processes: (i) Mott’s theory for electron-atom interactions, (ii) Ritchie’s theory for electron-electron interactions, (iii) Fröhlich’s perturbation theory for electron-phonon interactions, and (iv) Ganachaud and Mokrani’s semi-empirical model for electron-polaron interaction.

2.2 Elastic Scattering

Due to the large mass difference of electrons and atomic nuclei, electron-atom collisions can be approximated as being perfectly elastic. A commonly used elastic scattering cross-section is the screened Rutherford cross-section, which has a convenient analytical form and is straightforward to implement in a Monte Carlo calculation. However, the screened Rutherford cross-section can be applied only to high-energy electrons and solids with a low atomic number. An alternative to the screened Rutherford cross-section is the relativistic partial wave expansion method of the Mott scattering cross-section. In Mott’s theory [21], the differential elastic scattering cross section (DESCS) with respect to solid angle $\Omega$ can be calculated as

[TABLE]

where $|f(\theta)|$ and $|g(\theta)|$ are the direct and indirect scattering amplitudes, respectively, given by

[TABLE]

[TABLE]

In these equations, $K$ is the momentum of the electron, $E$ the total energy, $m$ the electron mass, $c$ the speed of light, $P_{l}(\cos\theta)$ the Legendre polynomials, and $P_{l}^{1}(\cos\theta)$ the first-order associated Legendre polynomials:

[TABLE]

The phase shifts $\delta_{l}^{\pm}$ can be computed by using the equation

[TABLE]

where

[TABLE]

$k^{+}=-l-1$, $k^{-}=l$, $j_{l}$ are the regular spherical Bessel functions, $n_{l}$ the irregular spherical Bessel functions and

[TABLE]

where $\phi_{l}^{\pm}$ is the solution of the Dirac’s equation which can be reduced, as shown by Lin, Sherman, and Percus [22] and by Bunyan and Schonfelder [23], to the first-order differential equation

[TABLE]

with $V(r)$ being the electron-atom potential.

In the present calculation, the atomic differential cross-sections for elastic scattering are extracted from the NIST Electron Elastic-Scattering Cross-Section Database [24] ranging from 0*∘* to 180*∘* for 24 incident energies between 50 eV and 1 keV (in increments of 10 eV from 50 to 100 eV, and in increments of 50 eV from 100 eV to 1 keV) and put into a data file in tabulated form. The elastic differential cross-sections for energies and angles other than those in the table can be calculated by linear interpolation accurate to two decimal places.

For compounds, the DESCS can be approximated by using the additivity rule, as the sum of the atomic DESCS of all atoms in the molecule. With the DESCSs of the single elements B and N, the electron-molecule DESCS is obtained for hexagonal boron nitride ($h$-BN) through

[TABLE]

Figure 1 shows the DESCS of 50, 100, 500 and 1000-eV electrons scattered by $h$-BN as a function of the scattering angle.

The total elastic scattering cross section for electron-molecule interaction can then be calculated as:

[TABLE]

where $N$ is the number of molecules per unit volume in the target, $\lambda_{el}$ is the elastic mean free path and $\vartheta$ is the polar scattering angle.

At low energies, the elastic mean free path calculated using the phase-shift method becomes very small –on the order of 0.1 nm–, which is much smaller than the interatomic distances and therefore unphysical. This is probably a consequence of using a rigid static potential, as well as neglecting dynamic effects such as that associated with the polarization of the electron cloud. In any case, electrons undergoing elastic collisions become confined in a very small region of space, effectively becoming trapped. Consequently, these processes do not contribute appreciably to secondary electron emission, but may slow down calculations significantly. To address this issue, in this work the elastic cross-sections calculated by the phase-shift method are multiplied by a cut-off function whose role is to gradually decrease the importance of the elastic effect at very low energies. This function must tend towards unity as the energy increases so that the behavior of the static potential used in the phase-shift method is recovered. The choice of function is of course not unique but an expression of the form

[TABLE]

is seen to behave well for this purpose ($\alpha_{c}$ is a dimensionless parameter and $E_{g}$ is bandgap).

The angle of electron scattering at a certain step can be expressed by a random number

[TABLE]

2.3 Inelastic Scattering

Inelastic scattering is characterized by the energy loss function (ELF) \operatorname{Im}\big{[}-\frac{1}{\varepsilon(\textbf{q},\Delta E)}\big{]}, where q is the momentum transfer and $\Delta E$ is the energy loss. The dielectric function $\varepsilon(\textbf{q},\Delta E)$ in the ELF reflects the response of a solid to an external electromagnetic perturbation. Due to the difficulties of determining the energy loss function experimentally, Ritchie and Howie suggested an approximate function from the optical dielectric constants by fitting the measured optical data into a finite sum of Drude-Lindhard model functions in the optical limit ($\textbf{q}=0$) [25]. In this fashion one can extend the explicit formula to the required \operatorname{Im}\big{[}-\frac{1}{\varepsilon(0,\Delta E)}\big{]} for finite q-values.

The ELF of the material is parametrized in terms of an expansion of Drude-Lindhard-type oscillators at the optical limit with $N$-term analytic form, which are directly obtained from the features observed in the reflection electron energy loss spectroscopy (REELS) spectrum:

[TABLE]

where

[TABLE]

and $A_{i}$, $\gamma_{i}$, and $\hslash\omega_{0i\textbf{q}}$ are the oscillator strength, the damping coefficient, and the excitation energy of the $i$th oscillator, respectively. The step function $\theta(\Delta E-E_{g})$ simulates a bandgap for the case of semiconductors and insulators so that $\theta(\Delta E-E_{g})=0$ if $\Delta E<E_{g}$ and $\theta(\Delta E-E_{g})=1$ if $\Delta E>E_{g}$. Although the dependence of $\hslash\omega_{0i\textbf{q}}$ on q is generally unknown, eq. (14) is generally accepted using $\zeta$ as an adjustable parameter. The value of $\zeta$ is related to the effective mass of the electrons, so that for free electrons $\zeta=1$ and for insulators with flat bands $\zeta=0$. The calculated fitting parameters of energy loss function for some pure elements and oxides are stored in an open online database [26, 27]. The parameters used to model the energy loss function of $h$-BN are listed in Table 1.

To obtain ELF(q,$\Delta E$) from ELF(0,$\Delta E$) we use Ashley’s model [28] by which the electron differential inelastic scattering cross section can be defined as

[TABLE]

Here $m$ is the electron mass, $e$ the electron charge, $N$ the number of molecules per unit volume in the target, $E$ the electron energy, and $\Delta E$ the energy transfer. The function $S(x)$ takes the form:

[TABLE]

so that the inelastic scattering cross section $\sigma_{inel}(E)$ for the electron-electron interactions can be written as

[TABLE]

$W_{min}$ is set to zero for conductors and to the bandgap energy for semiconductors and insulating materials; $W_{max}=E/2$ is the maximum energy transfer.

[TABLE]

Figure 2(a) is the ELF of $h$-BN in the optical limit; Figure 2(b) is the inelastic mean free path of $h$-BN as a function of primary electron energy.

For the stopping power $SP=-dE/ds$ , the following expression is used [29]:

[TABLE]

where

[TABLE]

In order to find the energy loss, $W$, of an inelastic collision of an incident electron with kinetic energy $E$, it is necessary to calculate the function $P_{inel}(W,E)$ providing the fraction of electrons losing energy less than or equal to $W$.

[TABLE]

where $R$ is a random number uniformly distributed in the range (0,1].

2.4 Phonon Excitation

At low energies, when $E$ does not exceed two or three times the value of the bandgap $E_{g}$, an electron has a high likelihood of interacting with the lattice vibrations. The interaction of a quasi-free electron with the longitudinal optical (LO) phonons in a polar medium can be treated by Fröhlich’s perturbation theory [30]. The interaction with the lattice is accompanied by the creation or by the absorption of a phonon. For the optical branch, it is reasonable to ignore the dispersion relation of the longitudinal phonon and to characterize it by the unique frequency $\omega_{LO}$. Then, an electron with energy $E$ has a probability per unit of path length to create a phonon of frequency $\omega$ (thus losing an energy $\Delta E=\hslash\omega$) given by

[TABLE]

where $a_{0}$ is the Bohr radius, $k_{B}$ is the Boltzmann constant, $\hslash\omega$ is the electron energy loss (on the order of 0.1 eV), $\varepsilon(0)$ is the static dielectric constant, $\varepsilon(\infty)$ is the high frequency dielectric constant and $n(T)=\frac{1}{e^{\hslash\omega/k_{B}T-1}}$ is the occupation number for the phonon level at temperature $T$, taken here equal to 300 K. For the present calculation, we assume $\varepsilon(\infty)=4.5$ and $\varepsilon(0)=7.1$ and only one LO phonon mode has been considered (with energy $\Delta E=\hslash\omega_{LO}=0.1$ eV). Since the phonon generation probability is higher than the absorption probability by a factor of about 10, the annihilation of the LO phonons along the electron path is neglected.

The polar scattering angle is given according to Llacer et al. [31] by

[TABLE]

where $E$ and $E^{\prime}$ are the electron energy before and after electron-phonon scattering, respectively. $R$ is also a random number uniformly distributed in the range [0,1].

2.5 Polaronic Effects

A low-energy electron moving in an insulating material induces a polarization field that has a stabilizing effect on the moving electron. This phenomenon can be described as the generation of a quasi-particle called polaron. The polaron has a relevant effective mass and mainly consists of an electron (or a hole created in the valence band) with its polarization cloud around it. The polaronic effect is important in the description of low-energy electron transport in insulators as it allows the description of the electron trapping and de-trapping necessary for the investigation of electric current inside insulators and charging-up phenomena. However, cross sections of these channels of electron interaction with matter are scarce in the literature. Here, a semi-empirical formula proposed by Ganachaud and Mokrani [32, 33] is adopted. They assume that the inverse inelastic mean free path that rules the phenomenon –and which is proportional to the probability for a low-energy electron to be trapped in the ionic lattice– is given by

[TABLE]

where $C$ and $\eta$ are constants depending on the dielectric material. [32, 33]

3 Monte Carlo Calculation

The probability for any given scattering to occur is in proportion to its cross section. Thus, specifying the cross section for a given reaction is a proxy for stating the probability that a given scattering process will occur. Figure 3 shows various scattering cross sections as a function of primary electron energy.

The stochastic process for multiple scattering is assumed to follow Poisson statistics. If $R$ is a random number uniformly distributed in the interval $(0,1]$, the step length $\Delta s$ is given by

[TABLE]

where $\lambda_{T}$ is the electron mean free path, given by

[TABLE]

The procedure to sample different scattering events follows the sequence:

[TABLE]

In this fashion, events are sequentially sampled and executed, conforming effective trajectories that are tracked until an electron reaches the surface with an energy larger than the workfunction or until the electron is thermalized inside the material. Electrons that escape the surface are tallied and the net yield s computed as the ration of the number of escaped electrons relative to the total number of primary trajectories generated.

Our model accounts for the main physical processes on which the secondary electron emission of metal oxides depends. However, several parameters appear in the empirical laws we have proposed. The physical meaning of these parameters is quite clear but their values can be found only by comparing the simulation results to the experimental measurements. Preliminary calculations have allowed us to estimate what can be considered as a set of reference parameters. These values will be further varied to check the influence of their choice, particularly on the secondary electron emission yields. The parameters for $h$-BN in present calculation can be found in Table 2.

4 Results

4.1 Flat Surfaces

The total secondary electron yield for ideally-flat $h$-BN surfaces is calculated for incident angles of 0*∘, 15∘, 30∘, 45∘, 60∘, 75∘* and 89*∘* measured off the surface normal, and incident energies in the range 50-1000 eV. In this work, the typical number of primary particles simulated ranges between $10^{4}$ and $10^{5}$, which generally results in statistical errors around 3$\%$.

Experimental data on $h$-BN secondary electron emission yield is scarce and subjected to high uncertainty, so first we test our models on a material system such as SiO2 for which such data exist. This is a first step aimed at validating our codes before the study of BN. The results for SiO2 are given in the Appendix, which convincingly demonstrates the validity of our models. On this basis, the simulations of 50 eV-1 keV electron irradiation on flat $h$-BN surfaces are performed.

At low energy regime, the simulation results are found to agree reasonably well with experimental data, which corresponds exactly to normal operating conditions ($<$100eV) of Hall thrusters. The agreement is slightly worse in the intermediate energy regime, which we rationalize in terms of the charging effect of insulators under electron irradiation. Again, at high energy regime, the simulation results converge to the experimental data. Note that under steady state, the total SEE yield at high temperatures should tend to unity. This can be attributed to the onset of a charge gradient due to the existence of holes created by the departure of secondary electrons from the lattice. This charge gradient creates a ‘shielding’ electric field that captures further SEE until charge neutrality is achieved again. Once the material is neutral, the process starts again, leading to an oscillatory steady state that keeps SEE balanced [34]. This picture can be altered by factors such as radiation-induced conductivity changes, sheath potential modifications, or slowly evolving internal charge distributions. In addition, the roughness of ’real’ experimental surfaces compared to the ideally-smooth ones in the model surely plays a role in the comparison. SEE yields as a function of $E$ for all angles of incidence considered here are given in Figure 4(b).

Surface plots of both the SEE energy distributions and the yields are given in Figures 7(a) and 7(b). As mentioned earlier, these data will be used in ray-tracing Monte Carlo simulations of SEE in arbitrary surface geometries. The details of the function fitting process can be found in our prior publications [17, 18]. The final expressions for the total SEE yield and energy distributions are generated using machine learning software [35]:

[TABLE]

[TABLE]

4.2 Micro-Architectured Foam Structures

Next we calculate SEE from microfoam structures with various porosities. The details of these structures are given in our past studies [17, 18]. A finite element reconstruction of the material is used to extract surface elements that may be intersected by electron trajectories. Special algorithms are then used to identify intersections between primary and secondary rays. Daughter rays are generated from parent rates using correlations (28) and (29).

Figures 6(a) and 6(b) show the secondary electron emission yield for solid volume fractions, $V_{f}$, of 4, 6, 8, and 10% in the 50-to-1000-eV energy range for normal and random incidence. The inset to figure 6(a) shows the dependence of the yield with the material volume fraction at energies of 50, 100, 200, 300, 400, 500, and 600 eV for normal incidence. In the high porosity range explored here the dependence of the SEE yield on $V_{f}$ is clearly linear in the high porosity range explored here. For the sake of comparison, the maximum SEE yield for $V_{f}$ = 4% (which occurs for $E$ = 700 eV) is approximately 1.3, compared with a value of 2.3 for the flat surface (from Fig. 4(b)). This decrease in SEE yield by about a factor of two is indicative of the potential performance gains that micro-architected surfaces might offer relative to fully dense surfaces.

5 Discussion and conclusions

The main objective of our work is to develop reliable physical models of secondary electron emission (SEE), to be applied to the calculation of effective SEE yields in micro-architected surfaces for space propulsion thrusters. The methodology relies on two distinct but complementary elements. The first is an experimentally validated theoretical model of electron scattering in solids. This model is built as a transport Monte Carlo simulator of individual electron trajectories in a solid, capturing the pertinent scattering mechanisms in terms of an interaction differential cross section that is integrated across the relevant energy and angular ranges. These cross sections reflect different elastic and inelastic (e.g. core electrons, valence electrons, polarons, plasmons, etc) scattering mechanisms in each material, which are formulated according to the best available physics. The second element of the methodology is a discretization procedure to represent arbitrary surface geometries in terms of discrete boundary elements. Both modules (physics and geometry) are coupled by way of a raytracing algorithm that captures the intersection of primary electron rays with different boundary elements. The material considered here is BN (used in current thrusters thanks to its low mass density, low thermal expansion, and highly dielectric properties).

Our main findings are that, in the primary electron energy range of interest ($<$100 eV) micro-foam architected surfaces are seen to decrease the SEE yield by about a factor of two. This is already a significant advantage over flat surfaces (or with as-fabricated surface roughness). As well, these micro-foams are seen to suppress the SEE yield peak typically observed at 400 eV of primary energy. Our results have been validated at two different levels. As mentioned above, the electron scattering model has been compared to experiments in smooth surfaces for both SiO2 and BN. As this is the ‘physics’ building block of the methodology, great care has been placed on ensuring that the parameters of the model are consistent with available measurements in each case. Second, we have compared results for the W micro-foam to measurements in ‘fuzz’ W surfaces (smooth W surfaces pre-exposed to a He plasma), which have a similar degree of porosity and morphology as foams.

Our approach has several advantages. First, it allows us to study any arbitrary surface morphology thanks to the finite element discretization scheme and the raytracing Monte Carlo method to generate electron trajectories. This endows the methodology with an extraordinary versatility, as complex surface geometries such as foams, pillars, fuzz, cells, etc., of any size can be treated in a straightforward manner. As well, our physical scattering model being capable of describing both metals and ceramics, we can now study materials of high relevance for electric propulsion. We believe that this approach will enable a rapid parsing of thruster lining material surface concepts prior to costly development and characterization, to narrow down the parametric space and accelerate materials development.

To summarize, our main findings are:

In the primary electron energy range of interest ($<$100 eV) BN microfoam architected surfaces are seen to decrease the SEE yield by about a factor compared to ideally flat surfaces. 2. 2.

These micro-foams are seen to suppress the SEE yield peak typically observed at 400 eV of primary energy. 3. 3.

Our results have been validated by comparison to experiments in smooth surfaces for both SiO2 and BN.

Acknowledgement

The authors acknowledge support from the Air Force Office of Scientific Research (AFOSR), through award number FA9550-11-1-0282 with UCLA.

Appendix A List of Symbols

[TABLE]