Model sensitivity analysis of Monte-Carlo based SEM simulations

The sensitivity of simulated scanning electron microscopy (SEM) images to the various physical model ingredients is studied using an accurate, but slow simulator, to identify the most important ingredients to include in a reliable and fast SEM image simulator. The quantum mechanical transmission probability (QT) model and the electron-acoustical phonon scattering model are found to have the most significant effect on simulated 2D and 3D metrology results. The linewidth measurement error caused by not including these models in the simulation is less than 2 nm. Specifically, it was found from a comparison to experimental data that the QT model is essential in accurately predicting particular signal features in linescans such as “ shadowing ” . The simulator is compared with two other publicly available simulators, JMONSEL and CASINO, where the first one is also based on first-principle physics models and the latter one is using phenomenological models. CASINO is the fastest simulator on CPU, but Nebula on GPU is two orders of magnitude faster compared to a single threaded CPU simulation. Only up to 6% speed increase has been achieved by different model choices.


Introduction
As semiconductor devices shrink in size, it becomes crucial to inspect with optimized system parameters to minimize measurement errors [1,2]. For example, for high-end transistors accurate edge positioning [3][4][5] becomes essential when structure sizes approach a few nanometers. In addition, the dramatically changing height/width (aspect) ratio of lines [6] requires full 3D metrology [7][8][9]. Synthetic scanning electron microscope (SEM) images produced by Monte Carlo simulations can be a great aid in analyzing such complex metrology tasks [10][11][12]. But, the resulting images will depend on the physical models included in the simulator. Therefore, it is essential to know what their influence is on the results.
Typically, Monte Carlo simulators contain models for electron generation (inelastic scattering), electron transport (elastic and inelastic events) and boundary-crossing events. Simulators based on phenomenological models can provide fast solutions, but their accuracy is not guaranteed unless they are carefully calibrated for particular purposes. On the contrary, simulators based on first principle physics models should lead to accurate results without any calibration, but they can be impractically slow. The question is whether all model ingredients in such accurate, but slow, simulators have a noticeable influence on the results. Some of them may well be left out, or can be approximated, to speed up the simulations. This way, a fast simulator can be obtained which may still be quite accurate. The goal of this study is to study how sensitive the results are to the various model ingredients and identify those which are essential to keep in a fast and accurate simulator of SEM images for 2D and 3D metrology.
Previously, Villarrubia et al. have reported a detailed sensitivity analysis of the effects of inelastic scattering models and surface potential barrier models on 2D linewidth measurements [13]. They considered inelastic scattering models such as continuous slowing down approximation (CSDA), a phenomenological model, and the dielectric function theory (DFT), a first-principle physics model, and concluded that the edge positioning error due to the model differences is at most 3 nm.
Verduin et al. have done a model sensitivity analysis by artificially changing the elastic and inelastic scattering cross-sections in line edge roughness (LER) simulations [14]. They concluded that the crosssections do have an impact on SEM image contrast, but not on the determination of LER.
To obtain 3D information from top-down SEM images, Arat et al. developed a method to deduce step heights from signal contrast in SEM images [15]. They used simulations to show that part of the secondary electron signal, when scanning over a step edge, is determined by the step height rather than the geometry of the step edge (see Fig. 1). In the figure, f s (x) is the linescan signal where x is the beam position, f b is the baseline level determined by the SE yield at the particular beam energy and the brightness level setting of the SEM, a 1 and a 2 are the integration boundaries between which the signal is most sensitive to height as discussed in [15]. The integrated area is used to quantify the height. They also demonstrated this experimentally. To accurately predict step heights from simulations, it is important to identify those models that influence the signal contrast most (the green area in Fig. 1).
This study will focus on the simulation of linescans over a Si line of rectangular cross-section, and the influence of the various first principle physics models of the simulator on the linescan secondary electron (SE) signals. The primary energy dependence of the linescans will be studied, and visualized such that one can easily judge the importance of each model ingredient for the prediction of a certain property such as linewidth, or LER or step height. Furthermore, the calculation times will be taken into consideration, and the results will be compared with results from two well-known other simulators, JMONSEL [10] and CASINO [16]. The essential models for an accurate and fast simulator will be identified and the possibility to speed-up the simulator without changing the accuracy will be discussed. In the next section the various model ingredients will be introduced and briefly discussed.

Models
Primary electrons incident on a material sample interact with the material and will be scattered (Fig. 2). The scattering events in which the energy loss is negligible are termed elastic scattering events. In these events, electrons mostly experience trajectory changes. On the other hand, the events in which there is a considerable energy loss are termed inelastic scattering events. In those events, electrons transfer energy to other electrons and experience both energy loss and trajectory change, or they can be trapped. The event in which an electron is incident on an interface, e.g. between sample and vacuum, requires a boundarycrossing model to determine the probability for transmission or reflection.

Elastic scattering models
The sensitivity of simulated SE signals to the elastic scattering models has not been investigated in the literature as much as the sensitivity to inelastic scattering models.
Elastic scattering occurs when an electron is deflected by the nuclear potential of a specimen atom without any exchange of kinetic energy. The frequency of the events and the probability distribution of deflection angles can be modeled by the Mott cross-sections, which are obtained using a computer program ELSEPA based on relativistic partial-wave calculations of elastic scattering of electrons by atoms [17]. The program also allows to include solid-state effects, using a Muffin-Tin potential as a representation of the potential well in a densely packed crystal lattice, which influences the total cross-section.
At low energies the electron couples stronger with the lattice as its wavelength becomes comparable to the lattice distance. Therefore, modeling the electron-phonon interaction is essential for the completeness of the physical picture at low energies. The elastic scattering approach by Kieft [18] and Verduin [19] is followed, using the Mott cross-section model for energies above 200 eV, quasi-elastic acoustic (AC) phonon [20] scattering cross-sections below 100 eV and an interpolation between the two cross-sections for energies in between 100 and 200 eV.
In Fig. 3, the effect of different elastic scattering model choices on the mean free path (MFP) is shown for silicon (Si). The best model choice, from a physics point of view, is Mott cross-sections with Muffin -Tin potential, including AC phonons for low energies (shown by the red line with crosses). The default model in ELSEPA is Mott cross-sections for free atoms, i.e. not taking the Muffin-Tin potential into account. Combining this with AC phonon scattering cross sections results in the dashed blue line in Fig. 3. The MFP is roughly 1.5 times smaller at high energies. For example, the MFP of a 10 keV electron is 10 nm (free atom) instead of 15 nm (Muffin-Tin).
Mott cross-sections with Muffin-Tin potential, but excluding the AC phonon scattering, leads to the MFP depicted by the black solid line in Fig. 3, at very low energies a factor of 2 lower than the AC cross-sections.   It is worth noting that the energy scale adopted in the simulations refers to the kinetic energy of an electron in the material taken with respect to the bottom of the valence band. This means that the kinetic energy of an electron should at least be equal to the barrier energy (12.44 eV for Si) to escape from the material to the vacuum. Therefore, the simulator stops tracing electrons with energy below the barrier energy. These electrons are trapped, which has no consequences for conducting materials, but it may have for insulators leading to charging effects. The latter are not considered in this study.

Inelastic scattering models
In contrast to elastic events, an electron loses energy and momentum in an inelastic event. These events play a significant role in SE generation. Energy loss processes include interaction with outer-shell electrons and plasmons, inner-shell electrons, and longitudinal optical (LO) phonons. The energy loss of an electron can be modeled discretely for each event, or as an average energy loss per distance traveled, i.e. using the stopping power. The two conventional approaches to model inelastic scattering are dielectric function theory (a first-principle physics model) and continuous slowing down approximation (a phenomenological model).
Villarrubia et al. studied the sensitivity of CD measurements to eight different inelastic scattering models, such as variants of the dielectric function theory (based on the Penn [21] algorithm) and phenomenological models based on fitting experimental SE-yield-versus-energy curves [13]. In the third variant of the dielectric function model, DFT3 in their notation, they followed the approach of Mao et al. [22] and separated electron-electron and electron-plasmon types based on momentum transfer. In case the momentum is conserved, it is assumed that the event is electron-electron scattering and the initial momentum of the SE is randomly sampled from the Fermi Sphere. Otherwise, it is an electron-plasmon event and the SE energy is randomly sampled from the Fermi sea electrons based on the free electron density of states.
The inelastic model used here is similar to DFT3. We use the full Penn algorithm (FPA) as described by Shinotsuka [23]. However, we restrict the energy transfer such that the primary electron cannot end up with energy below the Fermi level. For the outer-shell electron excitation, we follow the approach described by Mao et al. [22] as in the DFT3 model. In the electron-plasmon event, the instantaneous (non-zero) momentum of the secondary electron is randomly chosen in an isotropic direction, as in DFT3. This has also been suggested by Dapor [24].
For the inner-shell electron excitation, the model is followed as described by Kieft et al. [18], where the energy transfer is larger than 100 eV. For energy transfers lower than the bandgap energy, again the approach by Kieft et al. [18] is followed, treating them as LO phonons [25]. These losses are about 100 meV [25].
In an elastic event, it is assumed that there is no energy loss. However, when an electron is deflected by the nucleus it does lose some energy due to the conservation of momentum [26]. Although these losses are very small (<1 eV) for light materials and low beam energies, they will still be modeled as "atomic recoil losses".

Boundary crossing models
When an electron is incident on an interface between two materials, e.g. the specimen-vacuum interface, it can either be reflected or transmitted. When it is transmitted the electron energy changes as a result of the different potential energy on both sides of the interface, and due to conservation of momentum its direction also changes, i.e. the electron is refracted. Quantum mechanically, part of the incident electron wave can be transmitted, and part of it can be reflected. This quantum mechanical transmission probability with a single barrier height model [19] has been implemented in the present simulator Nebula, which is available on Github [27].

Methods
The simulation code has been written such that individual parts of the physics models can be included or ignored; in other words, "knobs" were implemented, which can switch model parts on and off. To judge the model effects on topographic features, a silicon line of 32 nm × 32 nm cross-section was chosen as the specimen, shown schematically in Fig. 4a. An incident electron beam, of 0.1 nm (FW50) Gaussian beamwidth, is scanning an area of 64 nm × 8 nm with 500 electrons per 0.5 nm pixel. Each scanline is taken at a different energy, 16 energies in total. All electrons emitted from the specimen per pixel are collected by a detector above the specimen, and those with energy less than 50 eV (the SE's) are taken as the signal. The 16 scanlines are shown in Fig. 4b, stretched in width for visualization purposes, and stitched together for increasing energy ordered from top to bottom, basically a multi-energy SEM image. The detected signal is normalized to the number of primary electrons per pixel and shown on a color scale.
The physics models to be switched on and off in the simulator are summarized in Table 1. The first column lists all models which are assumed to describe the physics of the problem in the best possible way, the first-principle models. By default they are all switched on, i.e. they are included in the simulator. The second column describes the model status of the simulator when the first-principle model of the first column is switched off, i.e. not implemented.
The effects that these models are expected to have on the simulated images will be discussed first.
Disabling the Muffin-Tin model means that the simulator uses a freeatom model, resulting in a corresponding decrease of the elastic MFP. This causes high energy electrons to scatter and change direction more frequently. The effect of this is not expected to be significant for the SE emission since the cross-sections at low energies are not affected.
Switching from AC phonon cross-sections to Mott cross-sections decreases the mean free path at low energies by a factor of ~2. As a direct result, the total mean free path decreases and the first expectation is a smaller scattering volume. However, the scattering probabilities of different types of events (inelastic & trapping) change as well. It is therefore not straightforward to estimate the scattering volume without running the Monte Carlo process.
When acoustic phonon energy losses are switched off the scattering becomes purely elastic. This may lead to an increased emission. Because there are not many energy loss channels for low energy electrons deep in the material, their scattering lifetime will be long before they get trapped. Therefore, the simulation time is expected to be higher.
The effect of the LO phonon losses is similar to that of the AC phonon losses. However, the energy loss functions (ELFs) show that these events are less frequent than the AC phonon losses. Therefore, the effect should be small. Similarly, atomic recoil events [26] should have not an important effect for typical CD-SEM energies of 300-500 eV since the loss is a few orders of magnitude lower than the energy of the electron. For instance, the loss is about 100 meV for a 300 eV electron.
When the instantaneous momentum of the SE's is switched off, electrons are modeled as being stationary before the collision. Although this physically does not make sense, the electron will lose its initial momentum any way after several collisions. Therefore, it is expected to have no effect on the results. Switching off this model should speed-up the simulator.
In the default model momentum is conserved after collision with an inner-shell electron. However, the initial momentum of the inner-shell electron is not known. The model assumes it is zero. The other extreme would be that the primary electron does not change direction after inelastic scattering. Although this is not a valid assumption, it will provide insight into how these two extremes influence the SE image.
Similarly, switching off momentum conservation during boundarycrossing, i.e. refraction, is not physically correct. But its effect will show the sensitivity of the results to emission angles. It is expected that the SE signal will change, especially at the shadow regions next to the sidewalls of the line.
Finally, switching off the quantum mechanical transmission probability will lower the emission yield because electrons then have a higher probability to be reflected from the interface. It is expected that the effect is stronger in the shadow regions next to the line because there will be a lower detection probability for the electrons which escape from the surface and then scatter from the sidewall and finally hit the detector.
How to analyze the effect of switching off a model component is best illustrated with an example. Fig. 5a shows the multi-energy image of the silicon line when all models are switched on. This way of visualization provides a nice overview of what a top-down SEM image of a rectangular line will look like at various energies. Fig. 5b shows the linescan results when the quantum mechanical transmission model is switched off, and Fig. 5c shows the difference, i.e. Fig. 5b subtracted from Fig. 5a. The shape of the line appears in Fig. 5c especially at low energies, which indicates that the topographic features are sensitive to the disabled model.

Results
A sensitivity analysis was performed with the nine different models listed in Table 1 using the aforementioned method. The resulting difference images are shown in Fig. 6. It is seen that the line profiles are influenced most by the exclusion of the acoustic phonon model, refraction, and the quantum mechanical transmission of an electron through interfaces (Fig. 6b, h and i, resp.).

Consequences for 2D measurements
To quantify the impact of excluding these models on 2D metrology, linewidths were measured. The measurement method employed turns out to be crucial. The 50% threshold method is chosen here because of its accuracy and simplicity [28,29]. Although it is a relatively a simple method, the exact definition should still be clarified. Two different definitions of the threshold method were found in the literature. The first one assumes that the edge position is where the signal is halfway between the global maximum and the global minimum of the edge profile, as implemented by Shishido et al. [28]. The second one estimates the edge position where the signal is halfway between the global maximum of the edge profile and the signal baseline, as used by Villarrubia et al. [30]. The former is easy to apply to signals where the global minimum of the signal is easily identified. The latter technique can be more robust for signals where the shadowing effect is not that apparent, as in high energy or noisy experiments.
For this study, at higher energies (>2 keV), the uncertainty in the determination of the signal minimum increases (as large as 4 nm) and measurements using the first method show variations up to 10 nm. The second method predicts the linewidth with a much better overall performance (error <2 nm) at all energies. The measurement quality could be improved by fitting the signal to a more complex function such as a sigmoidal [31,32] or a double-Gaussian function [14,33].
The comparison here is focused on the energy range between 100 eV and 1 keV which is relevant for CD-SEMs. Fig. 7 shows that the CD deviation from its nominal value is least when all models are switched on and larger when quantum mechanical transmission probability or electron-acoustic phonons cross-sections are switched-off, although almost all measurements are within pixel uncertainty (±0.5 nm). The same comparison is given for the second method in Fig. 8. These measurements are slightly less accurate but qualitatively very similar.
From this comparison, it is concluded that the AC phonon model has a bigger effect on the CD than any other model and the measurement bias is less than 2 nm. The conclusion is that switching off the model ingredients shown to be relevant in the SE images does not have a significant effect on 2D measurements such as the critical dimension, in agreement with what has been observed before by Villarrubia [34] and Verduin [14]. Furthermore, at CD-SEM energies, the measurement method based on global maximum and minimum leads to slightly more precise results. However, this difference is within the uncertainty of discretization. Therefore, the effect is insignificant.

Consequences for 3D measurements
As an example of 3D measurements determination of step height is chosen. In Fig. 1 it was argued that step height can be estimated from a top-down SEM image using the scanline profile [15]. To evaluate the impact of model ingredients on 3D metrology applications the attention has to be on contrast changes in signals. For instance, the effect on a linescan profile at 300 eV of switching off acoustic phonon cross-sections is plotted in Fig. 9. The baseline signal has not changed in the intervals − 32 to − 16 nm and 16 to 32 nm. But the signal from the top of the line is significantly altered since the Mott MFP's (replaced model) are much smaller than the electron-acoustic phonon MFP's. In the low energy regime, electrons undergo more elastic scattering before they are trapped and spread over a larger area in the material. When the electron beam scans over the line, this leads to more SEs escaping from the sidewalls of the line. The magnitude of the signal has increased as well as the slope of the "edge blooming" part of the signal. Both are important changes that seriously influence model-based 3D measurements such as step heights [15].
Similarly as the AC phonon model, the quantum mechanical transmission (QT) model has an essential impact on the signal profile (Fig. 10). Switching off the QT model results in more electron emission from both the line and the substrate, be it with different ratios. In this case, electrons which would usually be reflected are now allowed to cross the interface. It is observed that the SE signal increases more on top of the line (from 0.7 to 0.85) than at the substrate level (from 0.6 to 0.65). The reason for this is the topography, allowing electrons to escape from both sidewalls of the line as well. For the same reason, the substrate signal increase is more, further away from the line. Although at each position next to the line the same number of electrons escapes from the substrate, the ones emitted close to the line end up in the line again and are not detected.
The conclusion is that the implementation of the QT model is important for the correct simulation of signal contrast.

Comparison with CASINO and JMONSEL
To put the results, obtained with Nebula and presented here, in perspective a comparison is done with two other widely used Monte Carlo simulators: JMONSEL [10] and CASINO (version 3.2.0.4). These packages include different scattering models or model assumptions. It is assumed that the default options provide the most accurate results to a basic user. Understanding and comparing the effects of model assumptions in these packages is beyond the context of this study. The present goal is to obtain results from these packages, as used with the default physics models they come with, and compare them to results obtained with Nebula such that we can understand the overall effect of the different models.
The CASINO simulations are used with the "MONSEL Defaults", a default user option in the program [16]. That means Mott cross-sections by Browning [35] for elastic scattering and inelastic scattering based on the stopping power as modeled by Joy and Luo [36], and modified by Lowney [37]. It is not stated in any of the publications on CASINO whether the quantum mechanical transmission probability has been included in the models. Therefore, it is assumed that it is not included.
The default version (the example script "LinesOnLayers.py") of JMONSEL makes use of the Mott cross-sections for elastic scattering, the full Penn algorithm for inelastic scattering and classical boundary transmission based on momentum conservation. In more detail, for elastic scattering the Browning empirical cross-sections were used for E < 50 eV, tabulated Mott cross-sections for 50 eV < E < 20 keV, and Screened Rutherford cross-sections for E > 20 keV, where E is the energy of the electron. The elastic scattering model is very similar to the one in CASINO. For inelastic scattering, there are two options: electron-electron scattering and electron-plasmon scattering. Electron-electron scattering has been modeled as described by Mao et al. [22]. In an electron-plasmon event, if the energy loss of the primary electron ΔE is larger than the binding energy (E b ) of an inner-shell, then the energy of the SE is ΔE -E b [38]. Otherwise, it is determined probabilistically based on electron densities between E f -ΔE and E f where E f is the Fermi energy [39]. The inelastic scattering model is similar to the model in Nebula.
The silicon line of Fig. 4a is used for the comparison. The signal profiles of a linescan at 300 eV are compared in Fig. 11. The signals are normalized to the number of primary electrons (80,000 for CASINO, 10,240 for JMONSEL, and 20,000 for Nebula). The electron beam, a Gaussian beam with 0.1 nm diameter (full width 50 -FW50), was scanned over 64 nm with a step size of 0.5 nm.
As seen from Fig. 11a the simulators give very different values for the SE yield at 300 eV. When the SE signals are normalized to the signal    Fig. 7, but the threshold (50%) method was applied between the global maximum and the baseline. value from Nebula in the middle of the line, the results from CASINO and JMONSEL agree surprisingly well (Fig. 11b), at the edge-blooming region and especially at the shadowing region, although very different inelastic scattering models were used. In contrast to the Nebula result, the signals obtained from the two other simulators, continue to rise between 20 nm and 32 nm although they should have reached their nominal emission value already. Fig. 12 now also shows the comparison to Nebula excluding the QT model, revealing the same signal rise as observed from CASINO and JMONSEL in Fig. 11.
The final results are very similar although the exact scattering crosssections are different. Apparently the different scattering models converge to similar results after several scattering events. But, based on this comparison, we conclude that including the QT-model makes a significant difference. So, perhaps there is no need to improve the accuracy of the cross-sections and phenomenological models work just as well for 2D and 3D metrology purposes. However, a QT-model should be included.

Comparison with experiment
For a final verification, the model effects were also compared with an experiment. Results were used from RM8820 [40], an industry-level calibration sample for metrology tools. The data were obtained from Fig. 2 in Postek et al. [41]. An integrated line scan of the middle polysilicon line is plotted in Fig. 13.     Many details about the inspection are given in the paper. For example, the sample geometry is described well and line cross-sections are also given. However, performing a quantitative comparison between experiment and simulation is an extremely challenging task. Many other parameters should also be known, such as beam spot size, or whether the region is imaged for the first time [42]. Besides, there are also details that are very hard to know such as detector collection efficiency and the surface condition of the sample. Last but not least, lack of knowledge of the detection settings such as brightness (voltage offset) and contrast (voltage gain) makes it very complicated to compare experimental data to simulated data quantitatively. Therefore, the best attempt to a qualitative comparison is performed here, to give an idea of the model effect on the signal trends.
The following parameters were used. The top and bottom CDs are modeled as 228 nm and 240 nm, respectively. The beam energy is 2 keV and 10,000 electrons are scanned over 430 nm with 1 nm pixel size. The spot size is 12 nm, which seems overestimated, but leads to a better match. Offset (brightness) and gain (contrast) have been arranged such that simulated signals match in the middle of the line and at the substrate baseline. As before, a comparison was done to a simulation with all models included and one with the QT model excluded. From Fig. 13 it is clearly seen that the quantum mechanical transmission probability has a significant impact on the signal. Only with the QT model included the experimental signal trend can be reproduced, especially in the shadowing area.

Model-affected simulation time
In literature, it is very hard to find explicit information on Monte Carlo simulation times. That is understandable because the computer system used, and the programming techniques, and such kind of nonmodel related issues, influence the results. To have a rough idea of how long a typical calculation can take with the packages used in this work, a timing benchmark was performed. For this comparison, it was deliberately avoided to do any user-level optimization, such as e.g. switching-off tracing secondary electrons which are far from the surface. Therefore, the results can be regarded as upper bounds to the calculation times.
A linescan of the same silicon line as above, of 128 pixels scan length, is simulated with 10,240 electrons/pixel for JMONSEL, and 80,000 electrons/pixel for CASINO and Nebula. It is worth mentioning that JMONSEL is a single-threaded program, whereas CASINO can utilize CPU parallelism, and Nebula can use both CPU and GPU parallelism. The CPU versions were run on a Windows 7 platform with Intel Xeon CPU E5-1620 v3 @3.5 GHz (8 threads), and the GPU version of Nebula was run on NVIDIA GeForce GTX 1080 with 8 GB memory. The CPU timings are normalized to the number of CPU cores for the comparison in Fig. 14. The Nebula GPU version timing was reported without any normalization to the number of threads because the number of active threads is changing during the runtime.
The comparison shows that CASINO is the fastest program running on a CPU. Assuming all the simulators were coded at a similar level of expertise, the phenomenological model leads to faster simulation times. However, Nebula has a very comparable speed and includes firstprinciple inelastic models. The GPU version of Nebula is almost two orders of magnitude faster than the CPU version and CASINO.
In Fig. 6, the effect that the model ingredients have on the linescans was shown. To speed-up the simulator, one could switch-off model ingredients that were shown to have no influence on the linescans. For instance, one could switch off acoustic phonon loss (Fig. 6c), optical phonon loss (Fig. 6d), atomic recoil losses (Fig. 6e), and instantaneous momentum of SEs (Fig. 6f). The GPU version of Nebula then becomes slower by 3% at 0.3 keV, 1 keV, and 20 keV, instead of becoming faster. In fact, switching off various loss mechanisms causes low energy electrons to scatter without losing much energy, such that they exist for long times at the expense of an increased computational load. Instead of completely switching off low energy electron loss mechanisms, approximating their net effect by a single energy loss value (50 meV in Fig. 15) in the AC phonon energy loss model, can speed up the simulations by 2% for 0.3 keV, 5% for 1 keV and 6% for 20 keV without significantly affecting the signal. The conclusion is that one has to be careful with switching off model components that do not seem relevant. It may even slow down the simulations, in which case it would be better to keep the models and optimize the simulation time in different ways.

Discussion
The results have shown that the scattering models have an impact on 2D and 3D metrology tasks. For instance, the linewidth measurement error due to the model variations presented here is less than 2 nm for a 32 nm line. This might still be acceptable but today, the semiconductor industry already targets sub-10 nm features [7,43]. In that size range a 2 nm error would not be acceptable as it would have a serious influence on the physical properties of devices, for instance on the thermal conductivity [44,45], or on the optical and electronic properties such as the dielectric function [7] or band-gap [46].
One may wonder whether phenomenological models can be used for dimensional metrology. So far, phenomenological inelastic scattering models were calibrated using experimental SE emission yields from semi-spherical surfaces [47]. But then still, a correct value of the SEyield does not necessarily mean that the electron scattering volume in the material is calculated properly. The scattering cloud of these models can be very misleading, as pointed out in literature [48]. Re-calibration, using more complex procedures, would then be required to make the simulation tool suitable for dimensional metrology tasks. In any case, as demonstrated in this work, the phenomenological simulation tools do need to include a quantum mechanical transmission probability model.  Dynamic effects were not investigated in this study. For example, if the material is not sufficiently conductive, the accumulated charge will lead to different signal contrast and critical dimension measurements. It that case, the impact of the model ingredients can be different. This should be investigated separately in a future study.

Conclusion
We studied the sensitivity of simulated SEM images to various firstprinciple physics models. For this purpose multi-energy SEM images were introduced, a very powerful visualization technique that allows a quick analysis over a broad energy range. The quantum mechanical transmission probability model (QT) was found to have an essential impact on images at low energies (up to 1 keV). It is also found that this model is essential in reproducing the "shadowing" features of the signal.
Including acoustic phonon (AC) cross-sections instead of Mott crosssections at low energies, yields significantly different results. A smaller (AC phonon) elastic mean free path leads to a smaller scattering volume for silicon. It affects the linewidth measurement error by less than 2 nm. We found that the difference between the two CD measurement methods based on 50% threshold is insignificant for beam energies up to 1 keV.
Both AC phonon and QT models are found to change the emission in a topography dependent way. This will affect 3D measurements based on signal intensities, such as height measurements from top-down SEM images.
Concerning simulation times, different model choices did not lead to a significant speed-up. For the effect of inelastic scattering model effect on CD, we compared two different simulators, JMONSEL and CASINO.
The comparison between these simulators shows that phenomenological inelastic scattering models are the fastest, but they often lead to different results. In a future study, it would be more accurate to implement these different models into the same simulator to analyze their effect on simulation times.
Comparing CASINO and JMONSEL to Nebula only led to similar results when the QT model was excluded. However, good agreement was found between experimental results and Nebula simulations with the QT model included. This demonstrates the relevance of the QT model and, assuming that no QT model has been included yet in CASINO or JMONSEL, if it would be included, both simulators would do much better in reproducing results from simulators based on first-principle physics models.

Declaration of Competing Interest
The authors declare that they have no known competing financial interests or personal relationships that could have appeared to influence the work reported in this paper.