A modified and calibrated drift-diffusion-reaction model for time-domain analysis of charging phenomena in electron-beam irradiated insulators

This paper presents a modiﬁed self-consistent drift-diffusion-reaction model suitable for the analysis of electron-beam irradiated insulators at both short and long time scales. A novel boundary condition is employed that takes into account the reverse electron current and a fully dynamic trap-assisted generation-recombination mechanism is implemented. Sensitivity of the model with respect to material parameters is investigated and a calibration procedure is developed that reproduces experimental yield-energy curves for uncharged insulators. Long-time charging and yield variations are analyzed for stationary defocused and focused beams as well as moving beams dynamically scanning composite insulators. © 2018 Author(s). All article content, except where otherwise noted, is licensed under a Creative Commons (CC BY) license


I. INTRODUCTION
Charging phenomena in insulators have long been studied due to their importance in such areas as scanning electron microscopy (SEM), memory-based technologies, particle detectors, ceramic surfaces, industrial cables, and the safety of spacecraft. [1][2][3][4][5] Probably, the earliest systematic studies of electron-irradiation effects in solids and charging phenomena in insulators, as parts of research on electrets, were carried out by B. Gross who has had a great impact on this research field. In his seminal works on irradiation phenomena 6,7 Gross investigated the electron trapping and charge buildup in high-resistivity solid insulators bombarded with energetic electrons. Further studies by Gross and coworkers produced new experimental techniques and mathematical models. [8][9][10][11] These and more recent [12][13][14][15][16][17][18] studies have not yet been able to provide a complete and coherent account of all observed phenomena. This could be due to the prevailing emphasis on static (stationary) models [19][20][21] rather than time-domain analysis. Studying the dynamics of charging in time domain is especially important in the analysis of response times in particle detectors 2 and in designing novel scanning strategies for SEM. 16 The existing dynamic models are either one-dimensional 14,15 or do not include some of the essential physical processes, e.g., dynamic recombination, trapping, etc. 17,18 While the prevailing semi-classical Monte-Carlo (MC) method 22 makes very few assumptions about the complicated electron-sample interaction process, realizing its full theoretical potential is technically very challenging. First of all, MC simulations are slowed down by the need to continuously update the long-range electrostatic potential. Secondly, it is computationally difficult to keep track of all the trapped and de-trapped electrons. Finally, achieving acceptable variance not only in particle numbers, but also in the times of events (e.g. emission times), may require a prohibitive number of statistical realizations. a Electronic mail: B.Raftari@tudelft.nl Instead of sampling the probability space, the drift-diffusion-reaction (DDR) approach, mainly used to model low-energy transport in semiconductors, 23,24 directly describes the space-time evolution of a continuous probability density function. The pertaining partial differential equations are obtained from the semi-classical Boltzmann equation applying the method of moments and a few assumptions about the distribution of particles over the momentum space. From the mathematical point of view the DDR approach assumes that the symmetric part of the secondary electrons (SE) probability density function is isotropic about the origin of the momentum space and is well-described by a shifted Maxwellian distribution.
In our previous publication 25 we argued that the DDR approach can be applied to electron-beam irradiated insulators if the initial high-energy transport stage is approximated by an empirical source function. We showed that this pulsed source function allows modeling both the short-time processes immediately following the primary electron (PE) impact and the long-time charge evolution due to sustained bombardment. Importantly, we demonstrated that the sustained irradiation can also be modeled by a continuous current source, which gives practically the same secondary electron (SE) emission current as the time-averaged SE emission produced by many single-impact pulsed sources.
However, the original DDR model 25 had serious shortcomings as well. First of all, it was not calibrated against experimental data. Although we were able to reproduce any SE yield at a chosen PE energy by tuning a single parameter -the surface recombination velocity (SRV) at the sample-vacuum interface, it was not clear which yield should be taken as a reference, since yields tend to change over time and depend on beam currents. Secondly, using the same SRV for all PE energies resulted in curves not fully compatible with published SE yield data over the whole range of PE energies. And more seriously, the model produced nonphysical results in the case of prolonged irradiation. Namely, the surface potential at low PE energies could reach very large positive values, which is not possible, since positive potential attracts secondary electrons back to the sample leading to the neutralization of any potentials exceeding ∼ 10 V.
We have identified the main reasons behind the bad long-time behavior of the original DDR approach. 25 These were the employed steady-state generation-recombination model, which is not really suitable for the analysis of transient effects, and the neglected reverse electron current. The electrons that are being pulled back to the sample by a positive surface potential are called reverse electrons (RE's). Incorporating fully dynamic generation and recombination processes is relatively easy. Here we employ the so-called trap-assisted generation-recombination model, which also reduces the number of equations to be solved and charge species to be tracked.
In hybrid MC-DDR methods 17,18 reverse currents are estimated with direct MC simulations of particle trajectories. Here we propose an alternative approach that keeps intact the self-consistent nature of the DDR model. Namely, we introduce a novel boundary condition at the sample-vacuum interface that accounts not only for the total number of electrons returning to the sample, but also for the spatial distribution of this reverse current along the sample interface.
We have also developed and implemented a clear calibration procedure for our DDR model. It uses the fact that certain types of yield measurements -the so-called standard-yield measurementscorrespond to the situation where single PE impacts happen sufficiently far enough from each other across the sample surface for their mutual interaction to be neglected. As our code is able to simulate single impacts, its calibration can be performed in this single-impact mode.
In Section II we outline the mathematical details of the modified DDR model. Its physical applicability is further discussed in Section III, where we also investigate the sensitivity of the model, explain our parameter choices, develop a calibration procedure against published experimental data, and compare our results for defocused beams with an alternative one-dimensional approach. Section IV presents further quantitative analysis of more realistic scenarios with focused stationary and moving beams, including a dynamic line-scan of a laterally inhomogeneous target.

II. MODIFIED DDR MODEL
In this section we recall the main features of the DDR model 25 and describe several significant modifications that have been implemented since its introduction.

A. Basic equations
The continuum approximation of the equilibrium transport of charged particles in insulators consists of both partial (PDE) and ordinary (ODE) differential equations augmented with a semiempirical source function accounting for the initial ballistic transport stage. The PDE's are the Poisson equation for the potential and the transport equations for the free charge density: with the constitutive relations for the current densities given by where q is the elementary charge, V (x, t) is the electrostatic potential, n(x, t) is the density of free electrons, n T (x, t) is the density of trapped electrons, p(x, t) is the density of free holes. As explained below, N T /2 is the equilibrium value of the trapped electrons density. Thus, the local excess of trapped electrons n T > N T /2 causes additional negative charging, whereas, the local lack of trapped electrons n T < N T /2 causes additional positive charging. The constant ε 0 is the dielectric constant of vacuum, the function ε(x) is the (static) relative permittivity of the sample, µ n and µ p are the electron and hole mobilities, and D n and D p are the diffusion coefficients.

B. Generation, recombination, trapping, and de-trapping
In the present modification of the DDR approach the dynamic trap-assisted Shockley-Read-Hall (SRH) generation/recombination model is implemented. Here we explain it along the lines of the PhD study by Robert Entner conducted at TU Wien. 23 An attractive feature of this model is that there is no need to keep track of trapped holes as all the relevant physics is already contained in the single equation for the rate of electron trapping: This process is coupled to the basic equations (1)-(5) and can be divided into four subprocesses illustrated in Fig. 1.
(a) Electron capture: An electron from the conduction band gets trapped at the band-gap of the insulator and the surplus energy of E c E t is transmitted to the phonon emission. The average rate of this process is R n = σ n υ th n(N T − n T ).
In the above equations: σ n (x) and σ p (x) are the electron and hole mean trapping cross sections, N T (x) is the total density of traps, n T (x, t) is the density of trapped electrons, υ th (x) is the thermal velocity, and n i (x) is the intrinsic carrier density. The spatial variable x indicates the possibility of spatial inhomogeneity, i.e., the presence of different adjacent materials.
The initial conditions on n and p at t = 0 are set as the corresponding intrinsic carrier densities of the materials under consideration, whereas the initial condition for n T has been derived based on the assumption of the initial steady state for the density of trapped electrons prior to the start of irradiation (i.e. ∂n T /∂t = 0) and is set to (11)

C. Charge injection
As has been discussed in our earlier work, 25 there are two possibilities within the DDR approach to model charge injection by a low-to moderate-energy electron beam via the source terms S n and S p . The first fine-scale model captures the discrete nature of the electron beam. The rate of particle 015307-5 Raftari, Budko, and Vuik AIP Advances 8, 015307 (2018) production resolved at the level of pulses produced by individual PE impacts is given by: where L(t) is the logistic function, i is the number of the particular individual PE, t i is the i-th PE impact time, and t g is the generation time of the electron-hole pairs, whose choice is discussed in the next section. With single-impact events, due to a relatively small number of produced electron-hole pairs, the continuous results of the DDR model should be interpreted as probability densities rather than particle densities, especially at lower PE energies. The second model is designed for studying the sustained bombardment of the sample and is based on the temporal average of the above pulsed source function: where j 0 is the average electron beam current. Both source functions contain the semi-empirical distribution function of the electron-hole pairs at the end of the initial generation stage: where E 0 and E lan = E 0 + V s are the beam energy and the effective landing energy of PE's, V s is the surface potential at the point of PE impact, E i is the electron-hole pair creation energy, R is the maximum PE penetration depth (discussed further in detail), x 0 is the center of the Gaussian distribution with the distance of 0.3R from the sample-vacuum interface, and A is the constant corresponding to the backscattering rate. In the hole distribution function g p the constant B is zero, however, it is different from zero in the electron distribution function g n accounting for the remaining PE's. More details about the calculation of the normalization constants A and B can be found in Ref. 39 and in our earlier work. 25 In the continuous irradiation mode we consider two additional modifications of the source functions. One pertains to a defocused beam such that the computational domain is smaller than the beam radius. In this case we use the following distribution function derived from (14) by integrating over horizontal coordinates and enforcing the conservation of the amount of generated electron-hole pairs: where and δ is the radius of the irradiated area (computational domain) at the surface. Accordingly, the beam current can be calculated as where i 0 is the current density. The formula (18) adjusts the beam current to achieve results independent of δ. If, on the other hand, the radius of a partially focused beam is smaller than the radius of the computational domain we resort to the following distribution: Here δ denotes the beam radius rather than the radius of the computational domain.

D. Sample-vacuum interface and reverse current
The boundary conditions on V, n, p at the interfaces of the sample with its holder and at the walls of the vacuum chamber are standard: Dirichlet at ohmic contacts and Neumann to simulate isolation and prevent any currents from flowing through the corresponding interface.
The sample-vacuum interface, however, is not common in DDR-type simulations. Previously 25 we have used a Robin-type boundary condition J n · ν = n (n n i ) for n > n i at this interface, which sets the SE current density at the level proportional to the charge density at the boundary with the surface recombination velocity n ≤ th controlling the magnitude of the current (ν is the outward normal vector at the surface).
The boundary condition J n · ν = n (n n i ) for n > n i assumes that there is a nonzero probability for emission as long as the density of free electrons is above the equilibrium value. This means that in our model the density n(x, t) describes the relatively energetic electrons that are able to overcome the electron affinity gap of about 1 eV. While immediately after the impact all generated secondary electrons may, indeed, be considered sufficiently energetic, they will gradually loose their energy during the drift-diffusion process through the sample and some of the electrons that have managed to reach the sample boundary will no longer have enough energy to escape the sample. The exact treatment of this gradual energy loss would significantly complicate the model. However, since low-energy free electrons are being rapidly trapped, we may skip the whole intermediate energy-loss process and model it as direct trapping of energetic electrons. Therefore, while fitting to experimental data we may expect to arrive at somewhat higher trapping cross-sections and/or trapping site densities compared to tabulated values.
As mentioned in the Introduction the boundary condition J n · ν = n (n n i ) for n > n i does not account for the reverse current. This leads to nonphysical results -very strong positive charging of samples under prolonged irradiation with low-energy beams.
Experiments show 56 that the energy of secondary electrons, although greater than the electron affinity of the material, rarely exceeds 10 eV. Therefore, even a relatively weak positive potential at the surface will pull back some of the emitted secondary electrons. To account for this reverse electron current we propose the following modified version of the Robin-type boundary condition at the sample-vacuum interface: where ∂V ∂ν and where V g is the applied potential at the upper boundary, which in the present study is set to zero (20) represents the reverse electrons current density. The function α controls the total magnitude of this current and the factor − ∂V ∂ν − determines its spatial distribution. We assume that reverse electrons will re-enter the sample only through regions where the normal component of the electric field is negative. The stronger is the local attractive electric field, the higher is the density of reverse current at that location. The function α(t) is chosen here in such a way, see (23)- (26), that the magnitude of the total reverse current varies linearly from zero, when the maximum surface potential V + (t) is below a certain value V min , to the value of the total outward SE current, when V + (t) reaches V max . This means that the net current through the sample-vacuum interface will be zero if V + (t) ≥ V max as all SE's leaving the sample will re-enter the sample as reverse electrons. Typically this leads to the surface potential never rising above V max (or V max + V g ). This choice of α(t) is not unique and could be further refined to take the energy spectrum of the SE's into account. One should also mention that, from the computational point of view, the mesh along the sample-vacuum interface should be fine enough in order to capture the gradient of the potential at the surface.
Since in the present case of zero extraction potential the attractive surface potential does not typically exceed the value of 10 V, the electrons of the reverse current will not have enough energy to cause any further ionization in the studied materials (ionization energy is in the order of 28 eV for both alumina and silica). Thus, no tertiary electrons will be generated in the sample. In the situations where it is not the case, i.e., with positive surface potentials larger than the material ionization energy, our method would require further modification in the form of an additional source term -the creation rate of tertiary electron-hole pairs.

E. Numerical solution
The first step toward obtaining a numerical solution of an equation or a system of equations is to investigate the existence and uniqueness of the solution. With regards to the present model, the consistency analysis relies on previously published results. A detailed investigation concerning the existence and uniqueness of stationary drift-diffusion equations can be found in Ref. 24. In a study conducted by Jerome 26 a mathematical analysis of a system solution map for the weak form of the DDR model, which forms a basis for the numerical solution of the model, has been provided. Also, in a follow-up study by Busenberg et al. 27 the wellposedness of a DDR model similar to the present one (with different source/sink terms) has been demonstrated.
The multiscale nature of the problem calls for the same strategy as was used in Ref. 25 regarding the scaling of variables. We apply the finite element method (FEM) for the numerical solution of the model equations and implement it as a solver within the COMSOL Multiphysics package. To balance the accuracy and the computational costs a careful strategy is needed. Our investigations show that the best (i.e., most reliable) results are obtained when we use adaptive (or local) mesh refinement, Lagrange shape functions, the fully coupled approach with the Newton-Raphson solver, and an adaptive time-stepping algorithm. The use of the adaptive grid refinement, although costly, alleviates the need for more sophisticated approaches, such as the traditional exponential fitting applied in semiconductor studies. 28,29 To reduce the computational cost of the adaptive mesh refinement, a simple strategy has been followed by using a combination of both adaptive and local refinement methods. Namely, the adaptive refinement was applied only in the initial simulations to identify the regions where a fine mesh is needed and then the local refinement is used in the follow-up simulations.
In some cases, we reduce the original 3D problem to a 2D problem in the (r, z)-plane of the cylindrical coordinate system as the geometry, boundary conditions, and the source are all axially symmetric. In the simulations of beam scanning over laterally inhomogeneous samples a fully threedimensional implementation was used.

III. CALIBRATION AND COMPARISON
In this section we provide a detailed account of the calibration procedure for alumina and silica, which employs experimental data from the studies by Dawson 30  explaining the underlying approximations of the DDR approach. We also compare our time-domain simulations with the results of an alternative one-dimensional approach.

A. Reproducing the standard yield of insulators
There are two main kinds of yield measurements from insulators: dedicated measurements with homogeneous pure samples 30,31 and SEM scans of insulator-containing targets. 32 In the former case often a great care is taken to avoid the charging effects. Typically, a defocused beam, a weak beam current, and a pulse of short duration are used. We define the SE yield free of charging effects as the standard yield and calibrate our code to reproduce such data as close as possible.
Parameters of standard-yield experiments 30,31 (current, pulse duration, irradiation area) imply that the probability for the primary electrons to land anywhere close to each other on the surface is very small. In fact, with defocused beams and low, short-duration currents the expected distance between PE's is large enough to permit neglecting mutual interaction between any two impact zones. This is the main reason why standard-yield measurement are free of charging effects. The singleimpact source function (12) with i = 1 allows to compute directly the expected number of emitted SE's per single (isolated) PE impact and is, therefore, applicable for modeling the standard yield. The cross-section of the axially-symmetric configuration used for calibrating the code is shown in Fig. 2. With sufficiently large computational domain the boundary conditions at the sides of the sample have no influence on the SE yield from a single PE impact and were set to Neumann (zero current).
There are two classes of parameters that may be tuned within their physically admissible ranges: those that determine the shape of the source function approximating the initial pair generation and the short-time high-energy transport stage, and the material (bulk) parameters that determine the transport and trapping/de-trapping at much longer time scales. While these time scales may seem well-separated, in the DDR model material parameters, especially the SRV n , have some influence on the initial transport stage as well.
The electron-hole pair generation time t g , defined as the time when all electron-hole pairs have already been generated, determines the time width of the pulsed source functions S n,p (x, t) and of the resulting SE emission current pulse. According to theoretical and experimental investigations by D. I. Vaisburd et al, 59   away from the spectrum of the primary beam as the result of collisions. However, up to 10 14 s most of the generated electron-hole pairs still have energies above 20 eV. Since "true" SE's dominating the emission spectrum have energies below 20 eV, most of them must be emitted after 10 14 s. It has also been found that 10 11 s after impact all generated electron-hole pairs are already thermalized with their energy spectra tightly localized around the edges of conduction and valence bands and trapping becomes more pronounced. Hence, the SE emission current pulse following a single PE impact should start after 10 17 s and be almost finished by 10 11 s. Moreover, if one aims at modeling "true" SE's, then the relative contribution to the total emission between 10 17 and 10 14 s should be small, compared to the contribution between 10 14 and 10 11 s. The DDR method produces exactly this type of pulses for t g set between 10 16 and 10 14 s, see Fig. 3.
Notably, some parameters such as the SRV as well as the penetration depth influence only the magnitude, not the duration of the emission current pulse, see  R tuned is the tuned penetration depth by DDR model (explained below). As can be seen from Fig. 5, trap density and capture cross section influence not only the magnitude but also the duration of the emission pulse with the expected tendency for shorter pulses with the increase in the trapping rate. We emphasize that the curves of Fig's. 3-5 should be interpreted in the probabilistic sense. Namely, the integral of these curves between any two time instants t A ≤ t B is the number of particles expected to be emitted from the sample surface during the corresponding time interval. Thus, the expected yield at a given PE energy can be computed by numerically integrating the emission current between t = 0 and some sufficiently large t > 10 11 .
Among the material parameters the carrier mobilities µ n,p have been determined with the highest precision and are simply assumed here to have the same values as in Refs. 25, 54, 57, and 58. Strictly speaking, these are the so-called low-energy mobilities and a more rigorous approach would be to use femtosecond and picosecond mobilities to model the transport of particles during the corresponding time intervals after the impact. 59 However, mainly due to the absence of data about these high-energy mobilities, here we use the same low-energy mobility values at all times. Nevertheless, the extremely short duration [10 17 , 10 12 ] s of the high-energy regime allows us to expect the approximations made in the DDR approach concerning the mobility values during this stage to be appropriate at least in the numerical sense. Notice that the changes in the yield do not exceed 0.1 when we vary the generation time t g between 10 16 and 10 14 s in Fig. 3. Thus, to have any significant impact on the yield the mobility would have to vary dramatically during this interval of time.
Parameters σ n,p and N T related to trapping weakly influence the magnitude of the emission current pulse and have, generally, large uncertainties. For example, in a study set out to investigate electron trapping in alumina 33 a relatively large variation of 10 21 to 10 15 cm 2 was reported for the electron capture cross section σ n in polycrystalline alumina. The same study also revealed that polycrystalline metal oxide materials like sapphire (α-alumina) generally have trap site densities N T in the order of 10 18 cm 3 . Insulating solids are often grouped into three types: crystals, polycrystalline and amorphous. 34 The trap site density has been estimated to be around 10 16 cm 3 for an alumina crystal, from 10 17 to 10 20 cm 3 for polycrystals, and around 10 21 cm 3 for an amorphous sample.
Probably one of the most comprehensive and systematic studies on charge transport and trapping in silica has been done by DiMaria and co-workers, [35][36][37][38] where a strong link has been identified between the capture cross sections and the nature of traps and the capture cross sections have been estimated to range from 10 18 to 10 13 cm 2 . Confusingly, the values and ranges for these parameters are not limited to the above mentioned estimates. [60][61][62] Another parameter that strongly influences the magnitude of the emission current pulse is the maximum PE penetration depth R. It determines the spatial shape of the source functions S n,p (x, t) and, therefore, the expected number of particles in the neighborhood of the sample-vacuum interface -the main contributors to the emission current. Many semi-empirical expressions have been proposed for R with the following general form R(ρ, E 0 ) = CE Γ 0 , where the values for C and Γ vary from study to study. 39,40 The constant C depends on the material and the exponent Γ has been mostly assumed to have a certain material-independent value, although, in some studies Γ has also been considered material-dependent. 41 The exponential expression for R emanates from Bethe's theory for the stopping power of charged particles in matter. Bethe's formula involves the density, atomic number, and atomic weight of the material. However, with the exception of studies by Kanaya and Okayama 42 and by Feldman, 41 the density of the material is considered to be the only parameter influencing the electron penetration depth. As of now the estimation of R is far from being certain as can be seen from large discrepancies in the penetration depth estimates employed by different authors, see Fig. 6 (left) and Fig. 7 (left). Apparently, similar disagreement concerning the penetration depth exists for metals as well. 53 Having identified n , σ n,p , N T , and R as the most uncertain of the model parameters influencing the magnitude of the emission current pulse we have performed a series of numerical experiments to determine the sensitivity of the DDR model output (SE yield) with respect to changes in these parameters. During these simulations some of the factors would be held fixed while other were varied with the goal to achieve the best possible fit between the computed SE yields and the experimental data. Three key points emerged from this analysis: • The shape of the yield-energy curve is influenced by the capture cross section and the density of traps. Namely, the larger are the trap density and the capture cross section, the lower is the high-energy tail of the curve. • The SRV affects the height not the shape of the yield-energy curve.

015307-12
Raftari, Budko, and Vuik AIP Advances 8, 015307 (2018) • Following any one of the published penetration depth formulas together with adjusting the values of material parameters within their permitted ranges does not produce yield-energy curves fully compatible with the experimental data over the whole range of PE energies. In view of these facts and the aforementioned uncertainty about the energy dependence of the penetration depth, fine-tuning R for each PE energy against the available experimental data was deemed by us as not only admissible, but also necessary. While tuning R other parameters have been fixed at the best found fitting values within their reported ranges. In particular, the electron and hole capture cross-sections were set at the frequently used value of 10 15 cm 2 . The trap site density turned out to be slightly higher than the reported upper bound 10 19 cm 3 , namely, 3 × 10 19 cm 3 , leading to the initial (equilibrium) density of trapped electrons of 1.5 × 10 19 cm 3 , close to what was used by us previously. 25 For PE energies higher than 2 keV the tuned penetration depths for sapphire and silica presented in Fig. 6 (left) and Fig. 7 (left), perfectly match those of Lane and Zaffarano 43 and are well-described by the formula of Young: 44 However, according to Young 44  The value of 1.66 was assumed for sapphire in several other investigations as well. 47,48 As can be seen from the insets of Fig. 6 (left) and Fig. 7 (left) at energies below 2 keV the tuned penetration depths deviated from the formula (27) and did not follow any other published formulas, while remaining within their range. Least squares fitting of a separate exponential formula of the type (27) to the tuned penetration depths for alumina and silica did not provide a satisfactory fit. This suggests that below 2 keV the energy exponent Γ is indeed material dependent. Hence, for calibration purposes penetration depths below 2 keV must be determined by fitting to the corresponding standard yield data, whereas above this energy the depth may be safely deduced from the formula (27).
With the tuned penetration depths the DDR method provides practically exact yield-energy curves for the whole range of PE energies. As was mentioned previously, the height of the yield-energy curve is mainly controlled by the SRV at the vacuum-sample interface. In Fig. 6 (right) we compare the output of the calibrated DDR model with the standard-yield data 30 (reported also in the database of Joy 49 ), as well as Monte-Carlo simulations 50 and the empirical formula of Agarwal 51 for alumina samples. As far as the DDR model is concerned the only difference between the unpolished and polished alumina samples is the SRV at the sample-vacuum interface (1.35 × 10 5 cm/s and 2 × 10 5 cm/s, respectively), which sounds reasonable, since surface polishing should not affect the maximum penetration depth.
Comparison of the results by the calibrated DDR model with the experimental data, 31 Monte-Carlo simulations, 52 and the formula of Agarwal 51 for a silica sample is shown in Fig. 7 (right). The tuned SRV at the silica-vacuum interface (0.8 × 10 5 cm/s) is lower than the SRV at the alumina-vacuum interface, indicating that e depends on both the material and the surface properties.
DDR simulations indicate that the first and the second unit yields for sapphire occur around 50 eV and 10 keV, respectively. For silica, the unit yields are observed below 50 eV and again at 4.35 keV. The values of the calibrated material parameters used in the present study are listed in Table I.

B. Continuous irradiation with defocused beams
Sustained bombardment, even with defocused beams, increases the probability for an incoming PE to fall in a close proximity to a previous impact zone. This will introduce the interaction between the previously trapped charges and the newly generated electron-hole pairs, so that the yield will vary with time.
At the moment a standard experimental procedure for measuring yield variation during sustained bombardment does not exist. Therefore, here we compare predictions of the DDR model with the   15 FD model is a current-density based formalism incorporating a detailed recombination and trapping mechanism. For comparison purposes we have considered the same material (amorphous alumina), current density, and the penetration depth formula (energy exponent in (27) is set as Γ = 1.55). We switch now to the continuous (time-integrated) source function (13), (15)- (18) suitable for long-time modeling.
Since the sample is amorphous alumina rather than sapphire, we choose a higher trap density of 10 20 cm 3 pertaining to the so-called shallow traps. 15 We set the emission velocity to 1.4 × 10 5 cm/s, close to what we have obtained above for unpolished sapphire. We note that in time-domain investigations the quantity of interest is not the charge yield, but the instantaneous ratio of the net SE emission current to the incident beam current -SE emission rate.
Taking into account that our approach is fundamentally three-not one-dimensional, the results presented in Figs It appears that the distance to the closest Dirichlet boundary, where the electric potential is maintained at some fixed value, e.g., zero, strongly affects the value of the surface potential at the sample-vacuum interface. Apparently, the most important parameter controlling the magnitude of the potential is not the total charge density, as one would naively assume, but the proximity to an ohmic contact. Most likely this is due to the image-charge effect, which partially screens the charge accumulated in the sample.
Things are complicated by the fact that providing at least one Dirichlet boundary condition is essential for the numerical stability (possibly, existence and uniqueness of the solution as well) of the DDR equations. In fact, in the case of an isolated sample, the numerical solution of our nonlinear problem is only possible through the so-called fully coupled approach, since the Dirichlet condition is associated with the Poisson equation, which, therefore, must remain coupled to the rest of equations during the iterations. As the boundary conditions at the interfaces of the sample are not of Dirichlet type, the only available remote surface to impose this type of condition is Σ 1 . Numerically, the screening effect of the Dirichlet condition can be minimized by placing the ohmic contact Σ 1 as far as computationally possible from the sample surface Σ 2 . Thus, we have placed Σ 1 at various distances from Σ 2 and, as can be seen in Fig. 10, the surface potential does reach significant negative values when the Dirichlet boundary is far enough. However, the time of collapse of yield to unity becomes even shorter in these numerical experiments and remains at odds with the previous one-dimensional FD simulations.

IV. FOCUSED AND MOVING BEAMS
In this section we simulate sustained irradiation of sapphire, silica, and mixed targets by focused stationary and moving electron beams with beam currents typical in SEM. In the first series of numerical experiments we use the axially symmetric target of Fig. 2 illuminated in the middle by a focused stationary beam. The distance between Σ 1 and Σ 2 is set to 0.1 mm. The samples studied in Fig's. 11-13 are isolated in the sense that the only boundary penetrable for particles in the samplevacuum interface Σ 2 . We consider the worst case scenario -perfect focusing -where all PE's hit the same spot on the sample surface. It is easy to deduce that defocusing will affect low-energy PE's with their small impact zones much stronger than higher energy PE's with their extended impact zones. To anticipate the results for more realistic partially focused beams the reader is advised to compare plots of this Section IV with those presented in Section III B. Figure 11 pertains to an unpolished sapphire sample irradiated at 5 keV, where the standard yield is around 1.7 as can be deduced from Fig. 6 (right). The net SE emission rate -yield for short -starts at the standard yield value, but after a certain interval of time drops to unity for all beam currents. The stronger the current, the shorter is the standard yield interval preceding the drop. In fact, it is easy to calculate that the drop in the yield happens after a certain amount of charge has been injected into the sample by the beam, which confirms conclusions of many previous investigations. The point of fastest decline in the yield roughly corresponds to 3 × 10 18 C of injected primary charge, i.e., approximately 19 primary electrons.
The DDR model does not substantiate the usual intuitive explanation 15,54 concerning the reasons behind this seemingly inevitable convergence of yield to unity with time. Commonly it is argued that the charging of the sample leads to the change in the landing energy of PE's, so that the yield no longer corresponds to the standard yield of that energy, but rather to another point on the standard yield curve of Fig. 6 (right). If, for example, the standard yield is greater than one, then the sample accumulates positive charge. The landing energy increases and one should look to the right along the standard-yield curve to know what the new yield should be. If, on the other hand, the standard yield is less than one, then the accumulated negative charge reduces the landing energy of the PE's, thus, moving to the left along the standard-yield curve. Thus, it is argued, a yield larger than one would eventually lead to a positive potential high enough to shift the landing energy of primary electrons to the second unity-crossing point on the standard yield curve. This argument, while intuitively appealing, does not take into account the spatial distribution, the dynamics, and the screening of charges.
The unity-crossing argument for small PE energies and larger than one standard yields has been previously criticized in Ref. 63, where the significant role of the reverse current in the yield drop was pointed out. We also believe that the PE landing energy change alone cannot explain the yield collapse, since in all our simulations the accumulated potential was never strong enough for the landing energy to reach a unity-crossing point.
For example, Figure 11 (bottom) clearly shows that the value of the positive surface potential is insignificant with respect to the PE energy and cannot possibly change the landing energy by so much that it becomes 10 keV -the second point along the standard-yield curve where it crosses the 015307-18 Raftari, Budko, and Vuik AIP Advances 8, 015307 (2018) unity line. What the DDR model shows, though, is that the drop in the yield coincides with the rapid increase in the reverse current, caused by the relatively weak positive surface potential attracting low-energy SE's back to the sample. Figure 11 (top) compares the contribution of the positive part n (n n i ) of the emission current density (dashed lines) to the net SE emission rate (solid lines).
The onset of the reverse current can be deduced from the emergent discrepancy between the solid and dashed curves, which coincides with the positive surface potential reaching the value V min = 1 V in the bottom plot of Fig. 11. Moreover, reverse current remains significant even after the net yield reaches unity. Thus, the unity yield is the product of a neat dynamic balance between the PE injection, positive outward SE emission, and the reverse current. The result is a steady-state process and the conservation of total charge (on average): one PE in, one SE out, and a conserved 'circular' current at the sample-vacuum interface. Figure 12 (sapphire) and Figure 13 (silica) correspond to the beam current of 100 pA and show the time evolution of the yield and potential for various PE energies. Comparison with the defocused beam irradiation of Fig's. 8 and 9 reveals a larger discrepancy in convergence times of the yield to unity for different energies in the focused beam case. The yield drops much sooner at lower PE energies than it increases at higher PE energies.
Similarly, from the surface potential plots of Fig's. 12 and 13 we conclude that the rise of sub-unit yields (above 10 keV for sapphire and above 4 keV for silica) to unity cannot be explained by the change in the landing energy, as the associated potential is never negative enough for that. Minimizing the screening by removing the Dirichlet boundary Σ 1 farther away from the sample-vacuum interface Σ 2 we could bring the surface potential in silica down to 15 kV, which, however, was still not enough to decrease the landing energy of PE's from 30 keV down to the required 4.35 keV, where the standard yield of silica is equal to one.
The unity yield appears to be a stable equilibrium state for isolated samples. The precise mathematical nature of this state requires further theoretical analysis, beyond the scope of this paper. At the moment we can conclude that there are two processes -reverse current and trapping -that play significant role in the approach to equilibrium. Starting from the initial state corresponding to the standard yield above unity (low PE energies), the sample gradually acquires positive charge, which turns on the reverse current and reduces the yield towards unity. The sample ends up positively charged in equilibrium with the surface potential bounded above by V max .
Starting from the standard yield below unity (higher PE energies), the sample accumulates electrons, which are being transported by electrostatic drift and collision-induced statistical diffusion towards the sample boundary. While energetic electrons coming from higher depths will initially be completely lost to trapping (which accounts for both the energy loss and actual trapping in our model), as traps get filled and the trapping rate reduces more energetic electrons will survive the transport to the sample surface. This leads to the gradual increase in the emission through the sample-vacuum interface, thereby increasing the yield.
In addition, trapping plays the role of a damping factor in the swing of the yield towards the unity. If the trapping is relatively strong, then we have the so-called overdamped oscillation, where the unity yield is approached from below and remains at the unit value as soon this value is reached (30 keV in sapphire). The sample ends up negatively charged in equilibrium with the values of the surface potential reaching a few negative kV, depending on the proximity of the Dirichlet contact.
If the trapping is relatively weak, then we have a simple damped oscillation, where the yield overshoots the unity. The resulting positive charging switches on the reverse current and the yield returns to unity. The sample ends up weakly positively charged at equilibrium in this case.
We note that the relative strength of trapping (trapping rate) depends on the local balance between the free-electron density and the number of available trapping sites.
The yield does not always have to drop/increase to unity, though. If it was the case, all insulators would look exactly the same under SEM. One possible scenario, where the yield may not converge to unity, is a (partially) grounded sample. The condition on charge conservation that requires a unit yield in an isolated sample may be relaxed if the sample is grounded. It is, of course, an open question whether a contact between an insulator and, say, a metallic grounded holder can ever be made efficient enough to allow for an easy passage of charges. Assuming for simplicity a perfect ohmic contact, the charge conservation no longer requires the exact unit yield for the sample-vacuum interface as additional electrons may enter the sample via the ground channel. These additional particles do not have enough energy to be able to directly contribute to the emission current. However, they efficiently reduce the charging by neutralizing previously trapped particles, thereby indirectly influencing the yield.

015307-20
Raftari, Budko, and Vuik AIP Advances 8, 015307 (2018) This situation is illustrated in Fig. 14, where we have imposed an ohmic boundary condition on the side of sapphire sample. Although the yield in such a grounded sample does not stay at the level of the standard yield at that energy, after a few oscillations it stabilizes at a slightly lower value, well above unity. This effect is also observed in samples with a relatively poor ground contact described by a Robin-type boundary condition with a low SRV.
Although, the surface potential does take longer to build up in a sample with contact, Fig. 14  (bottom), the behavior of the surface potential at the injection point is not very revealing. It is, perhaps, more instructive to look at the distribution of the total charge at the surfaces of isolated and grounded samples under identical irradiation conditions. While the surface potential is weaker in the grounded case, Fig. 14 (bottom), the images of Fig. 15 explicitly show that the amount of accumulated positive charge at the surface of a grounded sample is higher. Also the spatial distributions of the surface charge are different. A large disk of positive charge surrounded by a ring of negative charge is seen in the grounded sample, whereas, in the isolated sample most of the positive surface charge is concentrated around the injection point followed by a weaker positive ring some distance away.
Another situation well-known to SEM practitioners where the yield does not drop/increase to unity is the rapid scanning of the sample by a moving focused beam. To simulate the scanning process the source function (14) has to be modified to account for the motion of the beam. This is achieved by setting x 0 (t) = x 0 + vt, where v is the velocity of beam displacement in the horizontal plane. Consider a 1 × 1 µm 2 sample surface imaged with a 1000 × 1000 pixels resolution at the rate of 30 frames per second. Then, the beam moves across the sample with the horizontal speed |v| ≈ 33 µm/s.
We consider an inhomogeneous sample consisting of adjacent blocks of sapphire and silica, see Fig. 16. Samples consisting of one insulator on top of another have been previously studied with a one-dimensional approach, 55 while vertical stacks of insulators, similar to the one considered here, have been recently investigated experimentally. 32 We simulate a single scan line through the middle of the sample perpendicular to the interface between the adjacent insulators. Across the vertical interface between the two different insulators the source function (14) exhibits a discontinuity due to the change in material density and the corresponding maximum PE penetration depth.  Since cylindrical symmetry is lost, the following DDR computations had been performed in the full three-dimensional mode. Figure 17 shows the yield as a function of the beam position along its trajectory for an isolated sample. These curves correspond to the intensity of pixels in a single-line FIG. 18. Surface potential at three beam locations (white circles) during sample scanning (5 keV, 10 pA). Top and bottom rows correspond to opposite scanning directions (indicated by arrows).

015307-22
Raftari, Budko, and Vuik AIP Advances 8, 015307 (2018) SEM image. The standard yields of both insulators at the considered PE energy are also shown as dotted lines. First of all we notice the difference between the left-to-right (from sapphire to silica) and the right-to-left (from silica to sapphire) scanning modes. This difference is easy to understand by looking at Fig. 18 where the images show the surface potential at the same beam locations during these two scans. Since the charging of sapphire is stronger than that of silica, the resulting residual charge strongly depends on the scan history.
Otherwise, the scans of Fig. 17 have several common features. One can notice higher yields in the neighborhood of the sample edges due to increased emission via the vertical interfaces. This is a well-known effect -the sample edges look brighter in SEM images compared to the rest of the sample surface. One can also see the drop of the yield to unity during left-to-right scanning due to continuous charging of the sapphire part. This charging also causes the yield in silica part to drop below its standard value. The relatively smaller charging during the right-to-left scanning does not allow the yield in silica to reach its standard below-unity value after the initial edge-related surge, and keeps the yield below the standard value when the beam crosses into the sapphire part. Additional simulations show that reducing the beam current (down to a few pA) while maintaining a high beam displacement velocity gives scans that truthfully reflect the standard yields of each part of the sample. Unfortunately, in practice this would, probably, result in a bad signal to noise ratio.
Our method as well as other simulation software would greatly benefit from publicly available high-quality time-domain data in addition to the already available standard yields. While we realize that direct time-domain sampling of detector currents may be difficult, it should be possible to collect and publish single line scans of ∼ 1 µm insulator targets for a range of dwell times (scan speeds).

V. CONCLUSIONS
The self-consistent DDR method proposed in Ref. 25 has been substantially modified in the present paper to include the dynamic trap-assisted generation-recombination model and a novel self-consistent boundary condition accounting for reverse current. The method has been calibrated against experimental data to deliver exact standard yields for alumina and silica samples over a large range of PE energies. For alumina and silica all calibrated parameters remain within or close to their reported uncertainty bounds, thereby further confirming the acceptability of model approximations. Time-domain simulations with defocused beams have been compared to the previously published results from a one-dimensional Flight-Drift model demonstrating similar long-time behavior. Our investigations so far show that the initial high-energy transport stage can, indeed, be approximated by a semi-empirical source function and low-energy material parameters, whereas, subsequent transport stages fall within the original domain of validity of the low-energy DDR method.
Simulations with stationary focused beams confirm that in electrically isolated samples the yield collapses to unity after a certain number of primary electrons, which depends on the PE energy, has been injected roughly at the same location on the sample surface. However, our simulations do not support the widespread intuitive explanation of this phenomenon in terms of the changing landing energy of PE's. The effect appears to have dynamic origins and is related to reverse currents and transient changes in the distribution of charges close to the sample surface.
The surface potential is strongly affected by the proximity of metallic grounded surfaces due to the associated charge screening. This may lead to misinterpretation of charging effects, if one relies solely on the surface potential measurements, but also may present an opportunity to alleviate SEM image distortions. Our simulations show that a good ground contact could also prevent the collapse of the yield to unity.
We have presented, probably, the first 3D simulations of a laterally inhomogeneous sample irradiated by a moving beam that take into account both the dynamic charge trapping/de-trapping and the reverse current. While to a certain extent the yield obtained during this realistic simulations could be interpreted on the basis of time-domain results with stationary beams, some effects are 015307-23 Raftari, Budko, and Vuik AIP Advances 8, 015307 (2018) unique to dynamic scanning. For example, the scan profile appears to depend on the direction of scanning. A recent review by Walker et al 16 mentions the lack of reliable simulations related to low-energy SEM studies. We hope to partly fill this gap with the present modified and calibrated version of the DDR method.