Assessing trajectory-dependent electronic energy loss of keV ions by a binary collision approximation code

The inelastic energy deposition of energetic ions is a decisive quantity for numerous industrial-scale applications, such as sputtering and ion implantation, yet the underlying physics being governed by dynamic many-particle processes is commonly only qualitatively understood. Recently, transmission experiments on single-crystalline targets (Phys. Rev. Lett. 124, 096601 & Phys. Rev. A 102, 062803) revealed a complex energy scaling of the inelastic energy loss of low-energy ions heavier than protons along diﬀerent trajectories. We use a Monte Carlo like binary collision approximation code equipped with an impact-parameter-dependent modeling of the inelastic energy losses to assess the role of local contributions to electronic excitations in these cases. We compare angular intensity distributions of calculated trajectories with experimental results for 50-keV 4 He and 100-keV 29 Si ions transmitted in a time-of-ﬂight setup through single-crystalline silicon (001) foils with nominal thicknesses of 200 and 50 nm, respectively. In these calculations, we employ diﬀerent models of electronic energy loss, i.e., local and nonlocal forms for light and heavy projectiles. We ﬁnd that the vast number of projectiles are eventually channeled along their trajectories, regardless of the alignment of the crystal with respect to the incident beam. It is, however, only when local electronic energy loss is considered that the simulated two-dimensional maps and energy distributions show excellent agreement with the experimental results, where channeling leads to signiﬁcantly reduced stopping, especially for heavier projectiles. We demonstrate the relevance of these eﬀects for ion implantations by assessing the nonlinear and nonmonotonic scaling of the ion range with the thickness of a random surface layer.


I. INTRODUCTION
When energetic ions propagate within crystals, they can penetrate much deeper if they are directed along specific crystal directions [1].The structure of a crystal inherently features orientations in which the matrix atoms line up into columns and sheets with long-range order; these are referred to as channels and planes, respectively.Projectiles entering these are steered by the symmetric potential and develop oscillatory trajectories along the channel axis [2].The channeled projectiles are only subject to small-angle scattering at large impact parameters and interaction distances with respect to atomic nuclei.Elastic collisions with host matrix atoms and the associated energy dissipation, called nuclear stopping, are to a large extent minimized, and as a consequence, the energy loss the projectile suffers is dominated by the ion-electron dynamics.From an applied perspective, ion channeling, which affects range distributions of implanted species, is a well-established technique for detecting lattice locations of impurities and assessing damage [3], and it has been listed as one of 75 major breakthroughs at American National Laboratories by the Department of Energy [4].Consequently, ever since the discovery of the channeling effect in the 1960s, it has been described by analytical theories and atomistic computer simulations and has been the subject of many experimental investigations [5,6].To date, ion implantation remains a key technique in conventional semiconductor processing [7] as well as in the emerging field of quantum information science, with very high demands on single-atom placement in crystalline materials like Si, SiC, and diamond [8][9][10].
The relative simplicity of channeling trajectories enables control over the interaction distance, creating an ideal scenario for the study of the trajectory dependence of electronic excitations.At high ion energies, channeling was shown to result in significantly reduced electronic excitations due to the suppression of core-electron excitations in close collisions [11,12].With decreasing projectile energy, this difference is expected to slowly diminish due the bound core electrons of the target.The remaining target constituents, which can be excited in binary Coulomb collisions, are only the valence and conduction electrons, which are accessible along every trajectory [13,14].However, a strong trajectory dependence of the inelastic energy loss is a prevailing effect observed for a number of projectiles in Si and SiC single crystals over a wide energy range, even below the Bohr velocity [15][16][17][18][19].An intuitive explanation for the observed difference in energy loss along different trajectories would be geometrical effects of a longer path length along nonchanneling trajectories, together with an increasingly more significant contribution of nuclear stopping, especially in the keV energy regime.Nevertheless, experiments conducted in transmission geometry allow for the selection of only relatively straight trajectories to effectively minimize path-length effects [20].Even under such conditions, for helium ions, the ratio between energy loss of channeled and nonchanneled ions was shown to undergo a maximum and again decrease with the projectile velocity.Increasing effects were observed for heavier primary ions, additionally showing nonmonotonic scaling with the projectile's atomic number.We reason that the observed energy scaling of the energy-loss ratios originates from increased local energy loss due to electron-hole pair excitations in a denser electron gas [21] as well as in discrete, locally induced charge-exchange processes enabled by the increased interaction times for the dynamic rearrangement of the electronic structures of the colliding partners, i.e., the formation of molecular orbitals [22][23][24].Projectiles traveling in nonchanneled trajectories will experience several close impact collisions, resulting in direct energy losses and an effectively higher average projectile charge state, consequently leading to a stronger Coulombic interaction with the surrounding electrons, and thus, to higher energy loss [25].In particular due to the second mechanism, i.e., the rearrangement of localized electron orbitals, when going towards lower ion energies, the pronounced impact-parameter dependence of electronic excitation emerges again but is prompted by a different mechanism than that for fast projectiles [15,16].Categorizing and predicting these strong effects, with a complex dependence on ion species and energy, can be considered highly advantageous for semiconductor processing.One way of supporting the observed coupling of elastic and inelastic losses at low energies is by means of computer simulations.
To date, the high number and complexity of discrete electronic-excitation events in ion-solid interactions make it impossible for them to be included in any large-scale simulations of ion transport.The most accurate description of ion-electron dynamics is provided by density-functional theory (DFT) and time-dependent DFT calculations [26][27][28][29].These approaches, however, commonly consider only simple straight trajectories and do not include any atomiccollision-induced charge-exchange processes.Notably, both of these limitations strongly resemble the previously described channeling situations, where, even for complex ions, excellent agreement between theoretical calculations and experiments are found [16].For the general situation, considering all trajectories, the best candidates for calculating large-scale transport are Monte Carlo (MC) binary collision approximation (BCA) algorithms.Ion transport in a binary collision framework was shown to be accurate for many ion-target combinations at energies way below the Bohr limit, typically all the way down to a few 100 eV/amu [30,31], yielding almost identical results to today's more advanced molecular-dynamic codes, while retaining its computational effectiveness [32].Along with the well-known SRIM code [33], many different codes for crystalline targets exist, for example, MARLOWE [34], Crystal-TRIM [35], KING-IV [36], and IMSIL [37], but the core algorithm remains almost the same.Most of the current MC simulation codes based on the BCA stem from the algorithms described by Robinson et al. in the mid-1970s [34,38].Modeling of the electronic stopping force has commonly not been changed since.The electronic energy loss, E el , subtracted from the ion energy during the discretized ion-atom interaction, can be treated as a nonlocal stopping force dependent only on the atomic density, N, and path length, r; as a local energy loss depending on the impact parameter, p, in the collision; or a combination thereof: The expression for local stopping has the form suggested by Oen and Robinson [38], where S E is a stopping cross section, f 1 is a weighing factor set to a value between 0 and 1, and s is a constant treated as a fitting parameter.A serves as a normalizing constant, which ensures that the average stopping cross section for an ion with a random trajectory is equal to known S E independent of the values of f 1 and s: where p max is the maximum impact parameter probed and α u is the screening length in the scattering potential.Due to the absence of an accurate analytical model for low-velocity electronic stopping, the expression of local stopping, with its exponential dependence on the interaction distance, was originally introduced only in an "arbitrary" form, extrapolating the known impact-parameter dependence at high energies, to follow the electron density around the target atoms approximately.The computationally challenging complexity of electronic excitations discussed above is simply embedded in the stopping crosssection value, S E , derived from experimental or semiempirical data as an averaged effect on the projectile traveling along the most-representative nonchanneling trajectory.Even if the individual energy-loss processes are expected to be discrete, the large number of collisions along a typical trajectory renders such a parametrization a feasible approximation.
Here, we assess the applicability of this model by comparing scattering distributions and energy-loss contrast maps acquired in ion-transmission experiments through single-crystalline silicon membranes and the same kind of data generated by the MCBCA simulation code known as simulation of ion implantation (SIIMPL) [39].In the first step, by comparing the intensity distributions, we assess the validity of the trajectories calculated by the BCA code.Second, we apply the local and nonlocal energy-loss model to study the induced energy loss along these trajectories.By decoupling the nuclear and geometrical energy losses from purely inealstic ones, enabled in BCA by changing the stopping model, we gain valuable information on the trajectory dependence of electronic excitations, which could help to disentangle different contributions to energy loss.Following this direct matching and verification of the necessity of the local treatment of the energy-loss processes, we provide a detailed look at ion trajectories developed inside the crystal during ion transmission and perform simulations that can be used for better predictions of implantation profiles or, e.g., studying the effects of surface oxidation.

II. EXPERIMENT AND SIMULATION MODEL
The transmission experiments providing intensity and energy loss maps were performed with the time-of-flight medium-energy-ion scattering (TOFMEIS) setup at Uppsala University [40].Based on a Danfysik implanter, a chopped beam of 50-keV singly charged 4 He and 100-keV singly charged 29 Si ions was transmitted through a single-crystalline silicon (001) foil with a nominal thickness of 200 and 50 nm, respectively.The angular distribution and exit energy of the projectiles are recorded by a large position-sensitive detector located 290 mm behind the sample.The detector covers a solid angle of 130 millisteradians, corresponding to a scattering angle of ±11.5°.A set of two perpendicular delay lines registers the projectile position after the impact on the microchannel plates, which provide the stop signal for the flight-time detection.The three-dimensional (3D) data are further processed to allow for the visualization of the contrast maps [41].By binning data in three dimensions, i.e., in energy and two spatial dimensions, an energy-distribution histogram of particles arriving at the same detector position is formed in every spatial bin.These histograms permit extraction of statistical values, such as the mean energy, that can be, in return, assigned to their respective positions and plotted as a two-dimensional map.
The silicon samples, purchased from Norcada [42], did not undergo any specific surface treatment prior to the experiment and were thus expected to exhibit a layer of natural oxide along with other common organic contaminants.Two independent studies using various ion beam analysis (IBA) techniques were performed to quantify the sample thickness and the level of contamination [16,43].As a result, a 20-Å-thick amorphous layer of silicon was included on top of the crystal in the simulations to mimic the presence of an oxide.
The SIIMPL code used in this work was originally developed by Janson as an alternative to the listed codes, with the intention of creating user-friendly software applicable to any projectile-target combination [39].SIIMPL has been successfully applied in multiple implantation studies [44], and the code itself is now published as open-source material [45].The latest modification to the software allows for the extraction of the final energy of the projectile, along with the angular coordinates of the transmitted projectiles in sufficiently thin targets [46].
The MCBCA algorithm operates on an event-driven basis, where the ion trajectory is determined through consecutive binary collisions.This approach stands in contrast to molecular-dynamics simulations, where the ion's path is discretized in time [32,47].The core of any MCBCA code is the algorithm that searches for the next collision partner.For the code to run efficiently, it is necessary to minimize the length of the list containing potential target atoms scanned during each collision event.However, this must be done without compromising the accuracy and reliability of identifying the correct target atom.The volume of the crystal that is being searched for collision partners is defined by the parameter p max .In the silicon crystal with the lattice parameter a 0 = 5.43 Å, a value of p max = a 0 /2 considers a sufficient number of colliding partners [48].In SIIMPL, the scattering integral in these binary collisions can be solved by two methods: by solving the integral numerically using the Gauss-Mehler quadrature method [49], or by using the computation-effective so-called magic formula developed by Biersack and Haggmark [50].The latter method, together with a screening potential reported by Ziegler, Biersack and Littmark [30], was applied.Thermal vibrations of the target atoms are a major factor for the dechanneling of channeled ions and should therefore be included in MCBCA simulations of crystalline targets.Compared to MARLOWE, SIIMPL uses an alternative algorithm for thermal vibrations, where the displacement is added to the lattice positions first, prior to the search for collision partners.The amplitude of the vibrations at room temperature was set to 0.077 Å, according to the work of Flensburg and Stewart [51].
The electronic stopping, S E , is expressed as S E = kE b , where the constants k and b are determined from experimental data for random or nonchanneled experiments.By demanding a linear velocity dependence of electronic stopping, i.e., setting b = 0.5, the values of k resulting from the conducted transmission experiments were 4.4 eV/(10 15 cm −2 ) for 50-keV 4 He and 6 eV/(10 15 cm −2 ) for 100-keV 29 Si [15,16,41].The trajectory dependence of electronic stopping is modeled by Eq. ( 1).For the local, i.e., impact-parameter-dependent energy loss, we have adopted parameters from the work of Hobler and Murthy, who proposed a generalized stopping model for the silicon target [37].The fitting parameter, s, from Eq. ( 1) is further developed into a function of Z 1 : where 0.3 is an original value introduced by Oen and Robinson and 1.37 is the fitted value extracted by simultaneously fitting 100 (or 111 ) and 110 channeling profiles by Hobler and Murthy.Above that, the model was shown to perform best when the total energy loss comprised of both local and nonlocal contributions.The weighting factor, f 1 , was therefore set to 0.65, corresponding to a 65% contribution of impact-parameter-dependent energy losses.In the following, we refer to this case as local.When considering purely nonlocal stopping, the contribution of local losses is set to 0%, i.e., f 1 = 0.The resulting impact-parameter dependence of the averaged effect of inelastic energy losses of He and Si in silicon, modeled with s from Eq. ( 3), is shown in Fig. 1.
In the case of nonlocal stopping, a simultaneous collision algorithm [34,52] was shown to lead to more accurate angular and energy distributions when calibrated against molecular-dynamics simulations [48].A simultaneous collision interval of ξ m = 0.25 Å, below which all collisions are treated as simultaneous, was applied.

III. RESULTS
Experimental scattering maps shown in Fig. 2 illustrate how the trajectories of helium ions transmitted through a 200-nm silicon crystal develop when decreasing the initial energy of the projectile.All experiments were conducted in a similar geometry, where the beam was intentionally not aligned with any major crystal axis or plane.As the projectiles traverse through the crystal, they lose energy and scattering by target atoms leads to a broadening of their angular distribution.At high energies, forward scattering into certain directions is prevented by the blocking effect, which leads to an intensity minimum in the image plane [53].The combination of all blocking directions reveals a real-space image of the crystalline structure of the silicon sample, nicely visible in Fig. 2(b).However, with decreasing energy, due to increasing scattering cross sections and increasing critical angles, the very same atomic rows that lead to blocking will start to guide the projectiles into the channels and planes they form [54,55].In other words, for a given channel, we observe a nonlinear phenomenon, i.e., the transition from blocking to channeling in this specific direction.In the intensity maps in Fig. 2, this effect is visible as, first, reduced average intensity turning to an above-average intensity for specific crystal directions.
Starting from the channels and planes in the proximity of the incident beam direction this effect will, with decreasing energies, propagate into crystalline directions at a larger angular distance from the initial beam direction.At sufficiently low energy in Fig. 2(e), a large number of projectiles no longer exit in the azimuthal symmetry around the beam axis but develop channeling trajectories along the low-index channels and planes.
To fully examine the capabilities of the MCBCA code in calculating ion transport on a large scale, we choose the low-energy scenario as the most complex but also most relevant for ion implantation.
Figure 3 shows both experimental and simulated results for 50-keV He transmitted through 200-nm Si.The maps are presented as a combination of the angular distribution of detected ions along with either intensity (top row) or energy (bottom row), specifically the mean energy loss along those trajectories.The initial direction of the beam is pointing towards the center of the detector and is indicated by the yellow cross.The angles between the beam and the [100] crystal axis are θ x = 7.5°and θ y = 12°, corresponding to a rotation around the {110} and the {100} planes, respectively.In comparison with the wafer's geometry consensus, this orientation would correspond to a wafer tilt of about 14°.Along with experimental data, two MCBCA simulations are presented: one where local electronic stopping is considered, and one where the electronic energy loss is exclusively nonlocal.Three million projectiles were directed onto the sample, matching the absolute number of projectiles detected in the experiment.The observed triangular shape in a real-space image of the crystalline structure is formed by the crossing of the {100}, {110}, and {131} planes.The crossing of the {100} and {131} planes at the top of the triangle corresponds to the 103 channel, while the 100 channel, normal to the sample surface, would be found in the bottom-left corner just outside of the region covered by the detector.Overall, good agreement is found between simulated and experimental maps for both intensity and energy, which show very similar features.The agreement in angular distribution indicates a high accuracy of the simulated trajectories, effectively resulting in an impactparameter distribution that is very comparable to the experimental situation.Beam divergence (not included in the simulation), in combination with geometrical straggling in the disordered surface layers, permits a fraction of projectiles to enter channels at a comparably large angular distance from the unperturbed primary beam axis.These effects could explain the somewhat lower apparent resolution of features in the experiment e.g., smearing out of the donut shape around the 103 channel and of the ridge along the {110} plane.A quantitative analysis of the effect of a disordered surface layer on particle energies is made in the discussion section.Little or no difference is observed when comparing the two simulated intensity maps, regardless of the different stopping models.However, modifying the choice of the energy-loss model has a distinct effect on the final energy distribution of transmitted projectiles.The contrast observed in the experimental energy map is well reproduced by simulations with the impact-parameter-dependent energy loss but is weaker in the nonlocal treatment.
To further quantify the effect of local contributions to electronic stopping, we compare the energy distributions along specific directions.Figures 4(a)-4(c) feature a compilation of energy histograms extracted from the three different-colored circles marked in the energy map in Fig. 4(d).The first notable observation is the large underestimation of the width of the simulated energy distributions when compared to the experiment.The experimental energy distributions contain contributions from the finite energy resolution of the detection system, as well as energy-loss straggling due to the statistical nature of the inelastic energy-loss processes, neither of which are explicitly considered in the simulations [56].The effect of local electronic stopping becomes apparent when comparing the respective positions and shapes of the simulated distributions along different trajectories.The growing separation of the three distributions using local stopping is a direct consequence of the impact-parameter dependence of the electronic energy losses responsible for the increased contrast in the energy map in Fig. 3.The weaker crystal-lattice contrast observed when applying only nonlocal stopping originates from a combined effect of the trajectory dependence of nuclear stopping and the reduced path length of the channeled projectiles compared to random trajectories.A good qualitative match in the centroid positions of the energy distributions is achieved for two of the plotted directions using local stopping: center and the 105 channel.The energy loss along the 103 channel is somewhat overestimated by the simulation.In contrast, purely nonlocal treatment does not provide a match for any of the directions and fails to provide even the correct order of the distributions.Figure 4(e) shows a comparison of the energy distributions of all projectiles acquired over the whole detector.The local treatment of electronic stopping results in the formation of tails toward both high and low energy loss, which indicates an increase in contributions from efficiently channeled projectiles and those experiencing high energy losses in close collisions, mimicking the shape of the experimental spectra.The nonlocal distribution shows a highly asymmetric peak with a tail extending only towards higher energy losses characteristic for nuclear stopping.
Figure 5 presents the corresponding results for 100-keV 29 Si transmitted through 53-nm Si.The angles between the beam and the 100 crystal axis are θ x = −7.5°andθ y = 12°.The geometry is symmetrical to the one in the Experimental results are presented side by side with maps simulated in SIIMPL using different energy-loss models.Upper row shows a comparison of intensity distributions and bottom row shows energy maps.Good agreement is found for intensity distributions.In the energy maps, both stopping models reveal qualitatively matching features to the experiment, but higher contrast is achieved with the local stopping model.
previous system but mirrored across the {110} plane.The characteristic triangular shape can again be identified.A substantial number of the total 4.8 million projectiles were again scattered to larger angles and experienced channeling.Greatly enhanced lattice contrast emerges in the mean energy map.The simulated angular distributions agree very well with the experiment, while the energy-distribution map with nonlocal stopping shows reduced contrast.The qualitative agreement with the experiment is however not as clear as in the helium case.The effect and importance of the local treatment become readily apparent when focusing on certain trajectories, as provided in Fig. 6.Excellent agreement between simulation and experiment is found using local stopping.Most consequentially, the complex structure featuring two maxima found along the 115 channel, as a result of contributions from channeled and nonchanneled projectiles, could be successfully reproduced by the simulation.Nonlocal treatment did not allow for any projectiles to exit the crystal with an energy loss lower than 15 keV and failed to reproduce the low-energy-loss peak.The widths of the energy-loss distributions are again underestimated by the simulation, except for the lowenergy-loss peak due to significantly reduced straggling along well-channeled trajectories [56].The shape of the integral energy distribution from the entire area of the detector shows good agreement between simulated and experimental data when considering local stopping.Note that the additional peak at an energy loss of 10 keV observed in the energy distributions in Figs.6(a) and 6(b) is not very pronounced in the overall energy distributions.This can be understood when looking at the intensity distribution in the same area in Fig. 5.It is only a few particular trajectories that feature drastically reduced stopping along this direction.

IV. DISCUSSION
The simulations of ion trajectories by the MCBCA code revealed that part of the reason for reduced stopping along channels was indeed, even for minimal final deflection from the incidence beam direction, due to a shorter path length and reduced energy transfer in elastic collisions with target atoms at large impact parameters.However, this effect itself, though prominent, especially for heavy projectiles, fails to explain the experimentally observed large differences in energy loss along certain trajectories.The observation that simulations and experiments agree only when applying a local contribution to inelastic energy loss is well aligned with and in support of our current understanding of the underlying interaction dynamics, namely, the strong coupling between elastic collisions and inelastic energy losses.Despite the simplicity of the energy-loss model from Eq. ( 1), it qualitatively captures the essence of the coupling problem.The parametrization of the local model developed by Hobler and Murthy [37]  provides a very good fit, in particular for the silicon projectiles used in this work.The discrepancy found for helium can be partially excused by the model being calibrated for 3 < Z 1 < 33.However, this finding can also serve as a reminder that the local energy-loss model is still a simplification of, in reality, complex nonstochastic processes of electronic excitations that will never have a universal character.As an example, while the Z 1 oscillations observed in the inelastic energy loss [57] can be modeled within such approximations, no general validity for a wide energy range nor all experimental geometries can be claimed [37].
The nonlocal and local parts, here treated as being energy independent, are expected to feature an actual energy dependence due to the increasing contribution of dynamic charge-exchange processes at lower velocities [15].Nevertheless, the transmission experiments through thin targets benefit from relatively small changes in the energy of the projectile as it traverses the entire thickness of the target.The hypothetical energy dependence of the terms in brackets in Eq. ( 1) can be, on one hand, safely omitted when matching results from a single transmission experiment, but, on the other hand, it can be assessed by studying its scaling under various conditions.
While the experiments from Figs. 2, 3, and 5 were conducted on a thin silicon membrane, a similar qualitative development of trajectories can be expected inside a crystal thick enough to stop the projectiles entirely.As a consequence of the unavoidable reoccurrence of channeling for slow projectiles, regardless of their initial alignment, all implantations in crystalline targets will always lead to increased range compared to an amorphous system, with range distributions deviating from predictions by, e.g., SRIM.The magnitude of reduced stopping along the channels and planes will govern the degree of range prolongation.The established and common method of testing and calibrating an MCBCA code is a comparison of simulated depth-concentration profiles with secondaryion mass spectrometry (SIMS) conducted on implanted samples [58].This process is, however, indirect while relying on a combination of two techniques, ion implantation and SIMS, each hampered by a set of drawbacks, e.g., alignment imprecision, damage accumulation, and matrix effects.Moreover, the MCBCA simulation needs to further rely on stopping-power data for the given projectile-target combination, which are often scarce, especially along different crystal orientations.Nonlinearities in the electronic energy loss scaling with ion energy have been commonly found at low energies [59,60] and are expected to lead to ambiguities in the obtained description when using the accessible integral parameter, i.e., the range from SIMS measurements.In contrast, the present experiments provide access to energy-loss values for all trajectories covered by the detector at specific ion energies.Random S e values from Eq. ( 1) can thus be extracted directly from the experiment.A combination of both approaches, i.e., assessing particle ranges and assessments of the impact-parameter dependence at specific energies, can provide a more holistic model for the combined impactparameter and energy dependence of electronic excitations.
The results from the combination of transmission experiments and MCBCA simulation can be further applied to offer an alternative visual perspective on, e.g., the role of a surface oxide on projectile trajectories and resulting range of implanted species.The effect of a surface top layer consisting of natural or intentionally added oxides on altering ion-implantation profiles, including scaling with implantation energy, oxide thickness, and ion species, was previously studied by a combination of BCA codes and SIMS [55,61].Using a screening oxide is a common practice in the development of silicon-based semiconductor devices [62] but more recently has also been applied to improve the placement accuracy of nitrogen-vacancy centers in diamond [63].Employing improved modeling for such predictions could be highly impactful in the design of experiments and the assessment of the feasibility of qubit architectures, e.g., with nearest-neighbor coupled qubits [64].For the system concerned in this work, we illustrate the effect of random silicon layers of different thicknesses added to mimic an unordered surface oxide on the intensity and energy distributions of transmitted silicon projectiles simulated in SIIMPL (see Fig. 7).From the intensity map, we observe the maybe intuitively unexpected effect that the increasing thickness of the random layer directs more ions into the planes and channels instead of randomizing their trajectories.The corresponding energy maps exhibit higher mean energy values along these trajectories.The resulting effect on the overall energy distribution is illustrated in Fig. 7(g).Notably, the channeled projectiles are not becoming more energetic in the presence of a surface oxide; instead, the mean energy value increases due to changes in the population of projectiles subjected to channeling, as shown in the inset of Fig. 7(g).A population of projectiles experiencing reduced energy loss ranging from 9 to 16 keV becomes abundant and would, in the implantation profile, form a more prominent concentration tail at increased depths.The presence of the screening layer therefore reduces the placement accuracy in a crystal oriented as random.

V. SUMMARY
We have shown how the combination of 3D MEIS transmission experiments with MCBCA simulations can provide valuable insights into the trajectory dependence of ion-solid interactions, specifically, into electronic energy loss.The necessity of modeling the inelastic energy loss as a function of impact parameter to match the experimental results, in a velocity regime where the core-shell electrons are no longer accessible, is in good agreement with the current understanding of local charge-exchange processes contributing to energy loss.Even though the stopping model is oversimplified, the BCA code proved to be successful at providing accurate information on the nontrivial trajectories developed in the crystal and their correlation with energy loss.The narrow energy window probed by the projectiles during transmission through the target, due to relatively small energy loss, makes the experiment perfectly suitable for testing and further improving the trajectory-dependent modeling of inelastic energy loss.Using the calibrated simulations for 100-keV 29 Si in silicon, we have further explored the role of amorphous surface layers on the projectile trajectories during ion implantation, verifying that an amorphous overlayer leads to an increased chance of channeling, and consequently, to an increased range, despite initially higher specific losses.

FIG. 1 .
FIG.1.Scaling of the averaged inelastic energy-loss function with the impact parameter, p, for 50-keV helium and 100-keV Si projectiles in a silicon target.Magnitude of the exponential decay with distance is modeled by parameter s.Electronic stopping cross section using nonlocal treatment becomes independent of the impact parameter and equal to experimental values orSRIM.

FIG. 2 .
FIG.2.Series of experimental scattering intensity maps of helium projectiles, illustrating the development of trajectories inside a 200-nm self-supporting silicon crystal, oriented oblique with respect to any major crystal axes, when decreasing the initial energy of the helium projectiles.Alongside the general broadening of the angular distribution at lower energies, the increasing fraction of channeled projectiles can be observed.

FIG. 3 .
FIG.3.Contrast maps of 50-keV He transmitted through a 200-nm silicon silicon crystal.Experimental results are presented side by side with maps simulated in SIIMPL using different energy models.Upper row shows a comparison of intensity distributions and bottom row shows energy maps.Good agreement for both intensity and energy is found for local stopping using parameter s = 0.25.

FIG. 4 .
FIG. 4. (a)-(c) Energy distributions of 50-keV He transmitted through a 200-nm silicon crystal along different trajectories marked in the energy map (d).Separation between peaks in experimental spectra is well reproduced by the simulation using local stopping.(d) Experimental energy map of 50-keV He transmitted through 200-nm Si crystal.(e) Comparison of the energy distribution of all detected projectiles.Full width of the energy distribution in the experiment is increased by energy-loss straggling and detector resolution, which are not included in simulations.

FIG. 5 .
FIG.5.Contrast maps of 100-keV Si transmitted through a 53-nm silicon crystal in a random direction.Experimental results are presented side by side with maps simulated in SIIMPL using different energy-loss models.Upper row shows a comparison of intensity distributions and bottom row shows energy maps.Good agreement is found for intensity distributions.In the energy maps, both stopping models reveal qualitatively matching features to the experiment, but higher contrast is achieved with the local stopping model.

FIG. 6 .
FIG. 6. (a)-(c) Compilation of energy distributions of 100-keV Si transmitted through a 53-nm silicon crystal along different trajectories marked in the energy map (d).Peak positions and the occurrence of low-energy-loss peaks are well reproduced by the simulation with local stopping.(d) Experimental energy map of 100-keV Si transmitted through a 53-nm Si crystal.(e) Comparison of the energy distribution of all detected projectiles, showing a close qualitative match between the simulation and experiment on both sides of the energy distribution.

FIG. 7 .
FIG. 7. Simulated intensity (a)-(c) and energy (d)-(f) contrast maps of 100-keV Si transmitted through 53-nm silicon with various thicknesses of the surface random layer.(g) Energy-loss spectrum visualizing the increase in the number of channeled projectiles experiencing lower energy loss with a growing thickness of the layer.Inset shows the effect along trajectories near the 115 channel marked with a circle in (c).