Large fraction of crystal directions leads to ion channeling

It is well established that when energetic ions are moving in crystals, they may penetrate much deeper if they happen to be directed in some specific crystal directions. This ‘channeling’ effect is utilized for instance in certain ion beam analysis methods and has been described by analytical theories and atomistic computer simulations. However, there have been very few systematic studies of channeling in directions other than the principal low-index ones. We present here a molecular dynamics-based approach to calculate ion channeling systematically over all crystal directions, providing ion ‘channeling maps’ that easily show in which directions channeling is expected. The results show that channeling effects can be quite significant even at energies below 1 keV, and that in many cases, significant planar channeling occurs also in a wide range of crystal directions between the low-index principal ones. In all of the cases studied, a large fraction (∼20–60%) of all crystal directions show channeling. A practical implication of this is that modern experiments on randomly oriented nanostructures will have a large probability of channeling. It also means that when ion irradiations are carried out on polycrystalline samples, channeling effects on the results cannot a priori be assumed to be negligible. The maps allow for easy selection of good ‘nonchanneling’ directions in experiments or alternatively finding wide channels for beneficial uses of channeling. We implement channeling theory to also give the fraction of channeling directions in a manner directly comparable to the simulations. The comparison shows good qualitative agreement. In particular, channeling theory is very good at predicting which channels are active at a given energy. This is true down to sub-keV energies, provided the penetration depth is not too small.


I. INTRODUCTION
Ion implantation is widely used in the semiconductor industry for materials modification [1][2][3].Moreover, radiation damage in nuclear reactors is formed primarily by energetic atomic recoils induced by neutrons [4][5][6].It is well established that when energetic ions or atoms are moving in crystals, they may penetrate much deeper if they happen to be directed in some specific crystal directions [7][8][9][10][11].This 'channeling' effect is utilized for instance in ion beam analysis methods such as RBS/channeling [12] and has been studied theoretically by analytical theories [9,13,14] and atomistic simulations [7,11,15,16].
The issue is of increased current interest due to the use of increasingly low energies in industrial ion irradiation [17] and also the interest in examining ion modification of single nanostructures where the energies are small to maximize the irradiation effects on the nanostructure [18][19][20].Some of these studies have shown a large variability in the radiation response of seemingly identical nanorods or nanoparticles [18,19,21], an effect difficult to explain by other means than channeling.
There are numerous experimental studies of channeling (see, e.g., Refs.[8][9][10]22]).However, few of them attempted to examine systematically where channeling occurs outside the Published by the American Physical Society under the terms of the Creative Commons Attribution 4.0 International license.Further distribution of this work must maintain attribution to the author(s) and the published article's title, journal citation, and DOI.principal crystal directions.A few experimental channeling maps have been measured that show the amount of channeling as a function of incidence direction [8,10], however, these did not provide a quantitative scale of the degree of channeling, making them of limited usefulness.
Molecular dynamics (MD) methods are well suited to study ion penetration in materials at energies where also multiple simultaneous collisions may be significant [11,[23][24][25][26][27].In particular, molecular dynamics in the 'recoil interaction approximation' (RIA), where only the interactions of the energetic ion with the lattice atoms is taken into account, has been found to be an efficient yet accurate way for describing ion penetration also at quite low energies [11,23,26,27].The approach has also been shown to agree well with the binary collision approximation simulations widely used at higher energies [27,28].
In this paper, we present a molecular dynamics-based approach to calculate ion channeling systematically over all crystal directions, providing channeling maps in several different systems.We assess and discuss the possible systematic uncertainties of the method and point out from the results how significant axial and planar channeling can be even at very low ion energies.We also compare our results with the predictions of channeling theory as developed in Ref. [14] in order to examine the possibility of estimating the effects of channeling without expensive MD simulations.
The theory of channeling was pioneered by Lindhard [13] and was further developed by Morgan and Van Vliet [29].These early works took into account the interaction of the ion with single rows or planes of atoms.They were able to predict the critical angle of channeling for low-index directions or planes and not too low ion energies [30].They also could estimate the fraction of channeled ions when the incidence direction was close to one of these directions.
By considering the interaction of the ion with all nearby rows or planes of atoms and with some improvements on the channeling criteria, one of the present authors achieved a more accurate description of channeling at low ion energies and for high-index directions and planes [14,31].This made it possible to draw channeling maps which for any incidence direction provide the information whether channeling is possible.
Channeling theory has mostly been applied to light ions in the MeV energy range and to 'close encounter events' such as nuclear reactions, large-angle ion scattering, or the generation of high-energy recoils [32].One purpose of this paper is to test channeling theory in a wider range of conditions, in particular, to investigate to what extent channeling theory provides useful information on energy deposition in decananometer-sized crystalline targets and on projected ranges in bulk samples.

A. Channeling theory
The theory of channeling as described in Ref. [14] can be used to examine channeling for any ion energy E and for any ion-target combination.It predicts which axial and planar channels are active and what the critical angles are for the active channels.In addition, the fraction of channeled ions can be calculated for a given incidence direction.
The calculation of these quantities requires, in a first step, the determination of the critical approach distance r crit to the rows or planes of atoms under consideration.The critical approach distance defines the closest distance of approach of the ion to a target atom so that channeling is possible.According to Ref. [14], for axial channeling the critical approach distance r 0 crit in a static lattice is given by while for planar channeling it may be calculated from or, in the special case where atoms in adjacent planes exactly face each other (such as in (111) planes of the diamond lattice), from (3) Input to these equations are the interatomic potential via the continuum potential U R1 (r) (see below), the ion energy E, and various lengths ( dR , d R , d 2 , dP ) characterizing the arrangement of the atoms in the rows and planes [14].
The critical approach distance r crit is then calculated from r 0 crit by considering thermal vibrations via the mean square atomic displacement x rms .x rms is provided by the Debye model [33].In our calculations we use the room temperature thermal vibration amplitudes given in Ref. [34] except for silicon where a value of 0.083 Å is used, corresponding to a Debye temperature of 490 K [35].
For given Miller indices, channeling is possible for ions that are perfectly aligned with the crystallographic direction or plane, if positions exist within the crystal which are further away from all atomic rows or planes than the critical approach distance r crit [14].Defining the channel center as the point that is furthest away from all atomic rows or planes, and denoting the distance of the channel center from the rows or planes r ch , channeling is possible if r ch r crit .In particular, ions that enter the crystal parallel to the channel at a distance r inc from the nearest row or plane are channeled if r crit r inc ( r ch ).The fraction of channeled ions f chan of a homogeneous beam perfectly aligned with a channel is given by the ratio of the area in transverse space (the 2D or 1D space orthogonal to the channel) where r inc r crit to the total area exposed to the beam.
Since the critical approach distance according to Eqs. ( 1)-( 4) increases with decreasing ion energy, a lower energy limit E min to ion channeling is given by r crit = r ch .Thus, E min can be calculated by inserting r crit = r ch in Eq. ( 4) and using the resulting value of r 0 crit in Eqs. ( 1)-( 3).When the ion is incident at an angle ψ inc > 0 with respect to the channel axis or plane, its motion in the transverse space may be described by a continuum potential U [=U R ( x) and U P ( x) for axial and planar channeling, respectively] [13].U is given by the sum of the continuum potentials of the nearby atomic rows or planes.For the Ziegler-Biersack-Littmark (ZBL) interatomic potential [36], the single-row and singleplane continuum potentials read and respectively.These expressions are analogous to those for the Moliere potential [37].Here, r denotes the distance from the row or plane, Z 1 and Z 2 are the atomic numbers of the ion and the target atoms, respectively, d R the mean distance of the atoms in the row, N 2 = 1/d 2 2 the areal atomic density of the plane, and K 0 (x) the modified Bessel function of the second kind and order zero.a i , b i , and a ZBL are the coefficients and the screening length of the ZBL potential.
In the continuum approximation, the transverse energy is conserved [13].When the ion moves in the channel, its position x in the transverse space and the angle ψ between the ion's direction of motion and the row or plane, vary, however, in such a way that E ⊥ remains constant.The transverse energy is determined at the ion's entrance into the channel by its incidence position x inc and angle ψ inc .When the ion moves into the target, its angle will sometimes get close to ψ = 0, in which case it attains the maximum value of the continuum potential along its path.Since U (r) is a monotonically decreasing function (except around the channel center), the ion reaches its minimum distance r min from a row or plane at this point.Thus, Since r min > r crit must hold for channeling to take place, the condition for channeling may be written as U (r min ) < U(r crit ) or, using Eq. ( 8), The most favorable point of incidence for channeling is at or close to the channel center, where the continuum potential attains its minimum U min .The maximum angle at which channeling is possible in this most favorable case is called the critical angle and is given by which follows directly from Eq. ( 9).From Eq. ( 9) also follows that for a given incidence angle ψ inc < ψ crit channeling occurs if For a homogeneous ion beam with a given incidence direction, the probability for an ion to be channeled (the channeling fraction f chan ) is given by the ratio of the area in transverse space where Eq. ( 11) holds and the total area exposed to the beam.
We construct theoretical channeling maps by first determining which channels up to a maximum Miller index of 10 are active for the given implant energy E. This is done by comparing the implant energy E with the minimum energy for channeling E min of each channel.For the active channels we determine percentiles of the continuum potential by choosing 10 000 random points in the crystallographic unit cell and calculating the corresponding values of the continuum potential U .Then we iterate over the polar and azimuthal angles of the incidence direction analogous to the MD simulations.For each incidence direction the angles ψ with respect to all active channels are determined and checked for ψ < ψ crit .When this is the case, the channeling fraction is calculated as the fraction of U values that fulfill Eq. ( 11), using the precalculated percentiles of U .If an incidence direction falls within ψ crit of more than one channel, the maximum of the f chan values is taken.These maps will be shown in stereographic projection, where circles on the unit sphere are mapped into circles in the projection, although features at larger tilt angles θ appear larger than features at smaller θ .

B. Molecular dynamics simulation
We used the MDRANGE code [23] that has been widely used to simulate ion penetration in materials [38,39] and shown to be able to reproduce very well experimental ion penetration depth profiles even in channeling directions in Si.With the proper choice of electronic stopping model, the mean ranges are found to agree within the statistical uncertainties, and the depth of the channeling tails within ∼10% with experiments [11,38,[40][41][42].It reads in an arbitrary atom structure and hence is well suited to study ion channeling in any crystal structure and any direction in this structure.
We set up the MDRANGE calculations in two different ways to estimate channeling effects.The basic simulation setup was following the practice described in Ref. [23].In the simulations the time step is selected adaptively using [23] Here v max is the maximum velocity in the system and F max the maximum force that affects the ion from any other atom.The first term ensures atoms do not move further than k t during one time step, while the second is necessary to handle strong collisions accurately.The last term ensures the time step does not increase drastically from one MD time step to the next.In this work, we used the standard values of k t = 0.1 Å and E t = 300 eV, except for H and He ions for which it was necessary to use E t = 30 eV to ensure energy conservation in strong collisions.
One of the basic outcomes of channeling is that ions penetrate deeper in the material.Hence using the mean ion projected range R p is a natural way to estimate the magnitude of ion channeling.However, this approach is not very well suited for nanostructures, since in many cases practically all ions penetrate the structure, making the ion range a meaningless quantity.For such cases, one can instead sum up the total nuclear energy deposition F D n from the ion to primary knock-on atoms in the structure.This is the full energy available for further modification of the nanostructure and hence a natural quantity for quantifying channeling effects in nanostructures.Although it would be possible to consider any shape of a nanostructure, in the current work we study irradiation of a thin foil of a given thickness t foil in the decananometer range for the energy deposition calculations.
The MD simulations used the universal ZBL potential [36], consistent with the channeling theory.In Appendix A we present a comparison of the results with another, densityfunctional theory based interatomic potential.The results there show that the channel positions and widths are not affected by the choice of potential within the statistical uncertainties.
In all simulations presented in the main text, we included the nonlocal electronic stopping power from the ZBL96 parametrization [43].In Appendix B, we examine whether the choice of electronic stopping model can affect the results.The results presented there show that while the choice of electronic stopping can affect the mean range by ∼10-20%, (as expected from previous works [44,45]), the angular width depends much weaker ( 5%) on the electronic stopping.
We also included random thermal displacements for the atoms calculated from the Debye model of lattice vibrations [23,33].In the current work, this was implemented always for 300 K, but the model can be easily changed to deal with any temperature.
The simulation was set up with an MD simulation cell with a [001] surface normal and tilting (θ ) and twisting (ϕ) the incoming ion direction.The atom coordinates were perfect crystal, and we did generally not implement possible surface relaxation or reconstructions.However, in Appendix C we present a test of the effect of surface reconstruction on the results.The results detailed there show that including a surface reconstruction can slightly ( 2 nm) affect mean ion ranges in channels, but has no appreciable effect on the angular dependence of channeling.
The initial position of the ion was chosen 3 Å above the top atoms in the cell in the z direction, with a random position generated in the x and y directions over one or four crystallographic unit cells.The range R p was always calculated as the ion range projected on the initial incoming direction [5].In simulations of the 'foil' samples the ion was defined to have left the sample when R p > t foil .In this way, the simulations mimic irradiation of a foil or nanostructure of known thickness as a function of the crystallographic orientation, which often is not known in the experiments.
The simulations can be carried out in the full range of θ from 0 • to 89 • , and ϕ from 0 • to 360 • , scanning both angles at intervals of 1 • .For each (θ,ϕ) combination, 1000-5000 ions were simulated.Naturally, for cubic crystals the simulated ϕ range can be reduced to the (0 • ,90 • ) interval or for highsymmetry lattices like FCC or the diamond structure to the (0 • ,45 • ) interval.We chose to always simulate at least the (0 • ,90 • ) ϕ interval, however, since any visible asymmetry between the (0 • ,45 • ) and ( 45• ,90 • ) interval would be a good indication of insufficient statistics on the plotted color scale.In this simulation setup, for the largest θ values all ions will be reflected.As an aside, we note that counting one specific initial ion trajectory as one MD simulation, the total number of MD simulations within the RIA approximation carried out for this paper was roughly 300 million.
The results were gathered as the mean ion range for the bulk penetration cases and average nuclear energy deposition in the foil cases.For each case, the uncertainty of each data point was calculated as the standard error of the mean of the R p or F D n values.When plotting channeling maps, the pixel values were calculated from the data on the two-dimensional (θ,ϕ) grid by bilinear interpolation [46].The number of interpolation points was chosen such that it is higher than the number of pixels used to plot the results.To map the data defined on the unit sphere to the plane, we use equidistant polar projection, where distances in the polar direction are undistorted, while distances in the azimuthal direction appear amplified for larger polar angles θ .(Thus, circles appear elongated in the azimuthal direction.)The color scale was always chosen such that bluish colors correspond to low degree channeling, reddish colors to a high degree.
For the definition of a channeling fraction we also need the nuclear energy deposition or projected range in a random target.The random targets were made by generating atom positions randomly in a three-dimensional periodic cell, with the only constraint that the atoms should not be closer than 2.1 Å to each other (for a motivation of this value, see appendix D).The number of atoms was chosen to give the same density of the material as the corresponding crystal cell.In principle this approach can still lead to enhanced atom transmission in some directions due to the finite thickness of the cell, i.e., some directions can by pure coincidence lack atoms.Systematic testing of ranges in random cells of different size (see Appendix D) showed that a cubic random cell of size 4 nm was sufficiently large to completely prevent such 'finite-size channeling artifacts.'The nuclear energy deposition or projected range was then calculated in these 4 nm random cells with identical ion energy-material density combinations as those in the crystalline cells.

A. Comparison with experiments
The main aim of the current paper is to probe which crystal directions are 'channeling' ones.Hence we consider here a comparison of experiments and simulations with respect to the angular dependence of channeling.
Experimental studies on channeling have mostly been performed with the background of RBS/Channeling, and therefore most experimental data have been obtained in a limited range of conditions.We compare here our range calculations with one set of experiments where the dependence on channel width was measured systematically and quantitatively: the angular scans for 0.8-2.0MeV He near different channels in Si by Azevedo et al. [47,48].We note that the comparison with backscattering cannot be performed directly, since the measured quantity was the He stopping power in a backscattering experiment, which is computationally too heavy to simulate with the MDRANGE method.On the other hand, the MDRANGE simulations are carried out until the ions have stopped and therefore contain information about the stopping power at all energies between the initial energy and zero.Hence the comparisons made are all strictly on the channel positions and angular width, not the absolute values of the range or stopping data.
The 0.8-2.0MeV He channeling cases obtained for the systematic angular dependence data measured by Azevedo et al. [47,48] were simulated for all the energies and same principal directions as in the experiments.To simulate irradiation in the [011] orientation, the implantation angles were varied around θ = 45 • and ϕ = 0 • and for [111], around θ = 54.73 • and ϕ = 45 • .Note that in the 100 direction experiments, the results in the experiments were obtained by averaging over several different twist angles ϕ, [47] while for the 110 and 111 cases the θ angles were tilted towards the 'axes parallel to the {100} and {110} planes, respectively' [48].Hence in the simulations for the [001] direction, we averaged the results over several ϕ values, and for [011] and [111], we used a single fixed ϕ in the simulations.In the latter two cases, this corresponds to tilting towards the {100} and {110} planes as in the experiments.
The results of the comparison are shown in Fig. 1, which compares the channeling ratio data for two representative cases [(a) and (b)], as well as the full-width-half-maxima (FWHM) values of the angular lineshape for all the cases considered (c).The FWHM were obtained by fitting a Gaussian shape centered at zero to either the experimental α or simulated ion range data (for the experimental data, the fit was performed to −α to fit a maximum rather than a minimum).We emphasize that this FWHM calculation does not involve any rescaling of the data, and hence the comparison in Fig. 1(c) does not involve any adjustable parameters.The data shows very good agreement between the simulations and experiments.In all cases, the experimental and simulated channel angular widths agree within ∼10%, and considering the ∼5% uncertainty in the experimental data, the observed differences may be purely statistical.In the experiments, this ratio α is the fraction of channeled to random stopping power, while in the simulations it is the difference between the mean ion range in the channel to that in the offchanneling direction, rescaled in the y axis to match the experimental data.The y axis of the simulated data on mean ranges has been linearly rescaled to be comparable to the experimental stopping power ratio.The lines are fits of a Gaussian profile to the data.(b) Same for 1.2 MeV He in Si near a 110 crystal channel.(c) Full-width half maxima (FWHM) of the experimental and simulated channeling results compared to each other for all the seven cases studied.The

B. Energy deposition in thin Au foils
The model case used for testing is the nuclear energy deposition by 1.7 MeV Au ions into a 20 nm thick foil of Au in different crystal orientations.This particular case was chosen due to recent experiments of Au bombardment of Au nanowires of 20 nm thickness, which show a major variability in the results (the experimental results will be published elsewhere [49]).
The full data set for the case of 1.7 MeV Au ions on a 20 nm Au foil with (001) surface orientation is shown in Fig. 2 as a function of the tilt (θ ) and twist (ϕ) angles off the [001] crystal direction.Due to the cubic crystal symmetry, the results above ϕ = 90 • were identical to the ones below, and hence are not shown.In this case, the results are shown both in a Cartesian [Fig.2(a)] and a polar [Fig.2(b)] plot of the (θ,ϕ) axes.Note that in the Cartesian way of plotting, the solid angle is not represented equally in different crystal directions.This is particularly well visible for θ = 0 • , i.e., irradiation perpendicularly to the surface, for which all twist angle ϕ values are of course crystallographically equivalent.Hence in the remainder of the paper we use the equidistant polar plots.In all the cubic crystal systems, there should be a symmetry about ϕ = 45 • which is indeed observed, indicating that the statistic of the data points is sufficient.
The results show that the nuclear energy deposition can vary more than an order of magnitude depending on crystal direction, varying between a minimum of about 10 keV for irradiation straight into the [110] channel, to a maximum of around 300 keV in several nonchanneling directions.A very notable result is the large areas of planar channeling, i.e., channeling in crystal planes between the principal crystal directions.In particular, channeling along {111} planes (connecting the 110 and 211 directions), {100} planes (at ϕ = 0 • and ϕ = 90 • ), {110} planes (at ϕ = 45 • and connecting 110 and 111 directions), as well as {311} planes (the remaining, mostly green paths) is observed.Planar channeling is well known to be significant for light ion irradiation [10,38], but observing the high degree of planar channeling also for a heavy ion like Au is somewhat surprising.
We note that the observation of strong planar channeling in Au around the principal 110 channel towards the 111 directions and a similar, but much weaker effect in the 100 directions agrees qualitatively well with the channeling map of Chadderton for protons in Au [8].A quantitative comparison is unfortunately not possible since Ref. [8] did not provide any intensity scale.
To get a more comprehensive view of channeling effects in Au, we also simulated the cases of 80 keV Xe and 10 keV H irradiation of 20 nm Au foils, and 30 keV Ga irradiation of 10 nm Au foils, see Fig. nanostructures at these ion-energy combinations [18,19].The H case was chosen to add a view of low-energy light ion channeling.The results show that, in spite of the wide variation of ion mass and energy, in all cases there is significant axial and planar channeling.

C. Low-energy ion ranges in Au
As a test of what is the low-energy limit of channeling, we simulated ion ranges of 5 keV and 500 eV Au in Au.The results in Fig. 4 show that even for these very low energies, significant directional effects are going on.For 5 keV, channeling along the 110 , 100 , and possibly the 211 directions is observed.For 500 eV, channeling appears to occur at θ ≈ 34 • , ϕ = 0 • .However, this direction does not coincide with a low-index crystallographic direction, so this cannot be classical channeling.Inspection of a set of ion trajectories showed that for incoming angles around θ ≈ 34 • , ϕ = 0 • , a large fraction of the ions are 'steered' at the surface into the [101] channeling direction, where they move slightly deeper than in a random direction.This effect can be explained by ion shadow cones [50], which focus a range of incoming ions into a narrow spatial and angular region.If this happens to coincide with the center of a 110 channel, the ion moves within this channel for some distance with reduced stopping, even if the initial direction was different.
The mean ion range shows a variation of more than a factor of 3 in both cases, and the channel widths increase with decreasing energy, as expected from the classical theories [cf.Eq. ( 10)].

D. Ion ranges in bcc: W
As a test of channeling in a bcc metal, we simulated the cases of 150 keV W and 300 eV D irradiation of W. The 150 keV case was selected because experiments of irradiation of W foils at this energy have been carried out recently [51], while the 300 eV D irradiation of W is a typical condition in fusion reactors [52,53].The results in Fig. 5 show that the 111 and 100 channels are very wide, while the 110 channel is clearly narrower.At 150 keV the mean ion range was more than 1800 Å in the 111 and 100 channels, while the mean range in nonchanneling directions (blue areas in the figure) is only about 140 Å.This corresponds to a tenfold increase in range.The order of the penetration depths observed in the main channeling directions is R 111 > R 100 > R 110 > R 112 , which agrees with those reported experimentally in the first systematic ion channeling experiments carried out in W [22].In addition, strong {110} planar channeling is observed.
For 300 eV D, channeling is clearly visible for the 111 and 100 directions and to some extent for 110 .This is surprising in view of the 500 eV Au on Au results, where variations of the range could not be correlated to crystallographic directions.The reason probably is that the ranges of 300 eV D in W are much larger than of 500 eV Au in Au.Obviously, channeling needs some travel distance to develop.

E. Ion ranges in diamond structure: Si
Ion implantation is one of the key techniques used in the semiconductor industry for chip manufacturing, and generally there is a desire to avoid channeling effects [38].Hence knowing the channeling directions in Si is important.Therefore, we also simulated channeling maps for a few representative cases of 10 keV ion irradiation of Si.
The results are illustrated in Fig. 6.They show very strong channeling in the 110 channels, consistent with previous experiments and simulations [11,14,[40][41][42]54].It is noteworthy, however, that all the ions show significant channeling also in several other crystal directions and that the 211 directions have stronger channeling than the 100 or 111 directions, consistent with results on B in Si [14].There are also clear planar channeling effects for all the ions studied.Comparing Figs.6(a)-6(c), it is observed that the critical angles increase with ion mass.

F. Ion ranges in hcp: Zr
As a test of channeling in a hcp metal, we simulated the case of 10 keV Zr ions on Zr.Because the hcp metal has less symmetry than bcc, we did this for two cells, one with the hcp a [2 11 0] axis as the surface normal and another one with the c [0001] axis as the normal.Due to the hexagonal symmetry  The white areas close to θ = 90 • correspond to directions where all ions were reflected and hence ion ranges could not be obtained.

A. Ion reflection
As a test of both MDRANGE simulations and channeling theory, we have simulated the reflection of 10 keV H ions from an (unreconstructed) (100)-Si surface.Ion reflection from a surface may be considered as a subproblem of planar channeling, as planar channeling is the oscillation of the ion between two planes.In the MD simulations, the ions are started at a distance of a little more than 3 Å above the topmost atomic layer with glancing polar angles θ .The azimuthal angle was set to ϕ = 10 • in order to minimize effects caused by low-index atomic rows.Figure 8(a) shows the mean minimum approach distance r min from the surface (blue line) together with its 1σ interval as a function of the incidence angle θ .The prediction of r min by channeling theory (red line) is obtained by evaluating Eq. ( 8) with U ( x inc ) = 0 and U (r min ) = U P1 (r min ) from Eq. ( 6).The MD results perfectly agree with channeling theory around θ = 89 • .For more glancing angles (θ > 89 • ) the MD results slightly deviate from theory, which might be due to the use of a cutoff distance in the interatomic potential.For less glancing incidence (θ < 89 • ), MD results and theory also differ very slightly.
The critical angle of channeling is the incidence angle for which r min equals the critical approach distance r crit .According to Eqs. ( 2) and ( 4 ψ crit = 90 • − θ crit = 1.2 • is read.Strikingly, for all angles θ θ crit the penetration coefficient equals zero, i.e., all ions are reflected, while just below θ = θ crit the fraction of ions penetrating into the target becomes nonzero.However, the rise in the penetration coefficient is only gradual and reaches 10% at θ = 88.4 • .Figure 8(b) shows that the average exit angle equals the incidence angle for θ θ crit as predicted by the channeling theory.Just below θ crit , the exit angle starts to deviate significantly from the incidence angle, particularly the mean exit angle θ exit becomes smaller (less glancing) than the incidence angle θ inc , or ψ exit > ψ inc .This means that on average the transverse energy increases upon reflection.Transferring this result to channeling between two atomic planes, it means that the critical angle marks the onset of the violation of transverse energy conservation rather than the onset of massive penetration of the atomic plane first approached.With increased transverse energy, the ion will approach the opposing plane of atoms after half an oscillation closer than the first one, and will be scattered even more there, and so on.Thus, ions with angles just above ψ crit normally require several oscillations to be dechanneled.These observations are exactly consistent with the concepts and methods used in channeling theory [14].

B. High and medium energies
Because of symmetry, for the cubic lattice systems investigated, all information other than surface effects is contained in the triangle formed by the Note that the color scales in the two parts of the figures have been chosen such that one end of the scale corresponds to a channeling fraction of 0 and the other to 1.In particular, it has been assumed that energy deposition or projected range equal to that in a random target corresponds to a channeling fraction of zero.Further details are given in Sec.V B.
The two parts should be exact mirror reflections at the line connecting [001] and [111].This is indeed quite well the case for the two foil simulations [Figs.9(a) and 9(b)], except that the channels appear a bit wider in the MD maps.The wider channels seen by MD can be explained by scattering of initially nonchanneled ions into channels.This effect is not taken into account by the theory.
In the projected range calculations [Figs.9(c) and 10] the less prominent axial and the planar channels appear weaker but wider in the MD data.The latter is likely due to the same reason as in the foil simulations; the weaker channeling fraction may be due to the less prominent channels being narrower than the primary channel, thus the probability of dechanneling is higher and the mean projected range is shorter.However, in all cases shown in Figs. 9 and 10, channeling theory is good at predicting which channels are active.This means that the theoretical calculation of the minimum energy for channeling E min is reliable.

C. Low energies
For 500 eV Au in Au [Fig.11(a)], theory predicts no active channels at all, while the MD results seem to indicate channeling along [111] and a direction at θ = 34 • , ϕ = 0, which does not coincide with a low-index crystallographic direction.The key to understanding this is the extremely small ranges, which are less than 6 Å: Over a few atomic distances classical channeling, which is the smooth oscillation between rows or planes of atoms, is not possible.Instead, the 500 eV Au results can be understood by the shadow cone effects discussed in Sec.III C.
The situation is completely different for 300 eV D in W [Fig. 11(b)], where the maximum range is 168 Å: Here the channeling directions are correctly predicted by the theory, and the MD data are again smeared up compared to the theoretical

D. Isotope effect
Standard channeling theory as outlined in Sec.II A predicts that ψ crit is independent of the masses of the ion and the target atoms (M 1 and M 2 , respectively).Based on MD simulations of channeling in carbon nanotubes, Zheng et al. [55], however, have found ψ crit ∝ (M 2 /M 1 ) 1/2 , and this result has been confirmed by Takeuchi [56].In order to test our simulations on the isotope effect, we have repeated the 10 keV H simulation in Si with T ions.Since M T /M H = 3, according to Zheng et al. the critical angles in the T simulation should be a factor of √ 3 smaller than those of the H simulation.In Fig. 12 the mean ranges of our H and T simulations are shown as a function of polar angle for two azimuthal angles, ϕ = 0 • and ϕ = 12 • .Here, the H results have been scaled by a linear transformation such that the range observed in nonchanneling directions and the range in the [001] direction coincide with those of the corresponding T ranges.No indication of an ion mass effect on the widths of the channeling peaks is observed in any of the channeling directions.
Explanation of the disagreement with Zheng et al.'s and Takeuchi's work [55,56] is beyond the scope of the present work.While the results presented by these authors seem conclusive, they have not given all details of their simulations (e.g., the length of their nanotubes), and they have not given a sound physical explanation.Therefore, a more comprehensive study would be necessary to resolve the discrepancy.

A. Maximum energy deposition and minimum projected range
According to common perception, tilting and rotating the sample 'far enough' away from the major channeling directions mimics the conditions found in a random target.This, however, is only true if extended regions of constant energy deposition or projected range exist away from the channeling directions.From the channeling maps shown in Sec.III it is obvious that this definition of a random direction is problematic, since these regions do not exist in many cases.Rather, nuclear energy deposition and projected range seem to vary continuously.Comparison with results of MD simulations in random targets show that the maximum energy deposition F D n ,max or the minimum projected range R p,min across the map may deviate significantly from that in the random target.This is illustrated quantitatively for a few representative Au and Si cases in Fig. 13, which shows the distributions of energy deposition and projected ranges scaled to the values in the random targets.The graphs clearly show that there may be significant portions of the distribution at F D n /F D n ,ran > 1 and R p /R p,ran < 1, that there is no sharp separation between channeling and nonchanneling directions, and that a large fraction of all directions are channeled (F D n /F D n ,ran < 1 and R p /R p,ran > 1).In the remainder of this section, we develop a systematic analysis of this fraction.
The data for F D n ,max and R p,min are compiled for various ion-target combinations in the fourth column of Table I in random targets (fifth column) reveals quite remarkable differences exceeding 50% in one case (319 versus 202 keV for 1.7 MeV Au in Au).In most cases there is more interaction with the target than in the random case (F D n ,max > F D n ,ran , R p,min < R p,ran ).This is explained by increased ion-target interaction in the angular vicinity of planar channels, an effect referred to as 'channeling shoulder' in the literature [15,34,57].It has been ascribed to 'quasichanneled' ions [58], ions that travel in the direction of a planar channel within or very close to an atomic plane.In one case (300 eV D in W) the minimum projected range is larger than the random range.This may be due to a large probability for ions to be scattered into channels.

B. Fraction of channeled ions
In order to extract channeling fractions from the MD data, we assume that ions are either channeled or nonchanneled, with respective well-defined energy depositions and projected TABLE I. Nuclear energy deposition or projected range, average fraction of channeled ions, and fraction of channeling directions for various ion-target combinations: In the fourth and fifth column the maximum energy deposition or minimum projected range in the triangle formed by the [001], [111], and [101] direction are compared with the corresponding value in a random target.In the sixth and seventh and in the last two columns the average fraction of channeled ions and the fraction of channeling directions, respectively, are compared between MD and theory.A channeling direction is defined by a channeling fraction f chan > 0.1 in the MD results and by falling within the critical angle of a channel in the theory.ranges (F D n ,chan and R p,chan , F D n ,ran and R p,ran ).The fraction of channeled ions f chan may then be defined for a given incidence direction by and where F D n and R p denote the average energy deposition and projected range, respectively, for that incidence direction.While F D n ,ran and R p,ran may be determined by MD simulations of random targets, definition of F D n ,chan and R p,chan requires another assumption.We assume here that in the best channeling direction, where we find the minimum F D n ,min or the maximum R p,max , the channeling fraction is properly described by its theoretical value f chan,max .This yields and Equations ( 15) and ( 16) provide the formulas for the calculation of the channeling fraction f chan from the energy deposition F D n and the projected range R p , respectively.These equations are also used by the software plotting the channeling maps, if the color scale limits are specified as (F D n 1 ,F D n ,ran ) and (R p,ran ,R p2 ), respectively, with and as we did in Figs.9-11.We note that Eqs. ( 15) and ( 16) yield negative values when F D n > F D n ,ran and R p < R p,ran , respectively, i.e., in the regions of the channeling shoulders.This is a consequence of the simplifying assumptions which we have made, particularly the disregard of the effects leading to channeling shoulders.However, even if the interpretation as the fraction of channeled ions is not possible for f chan < 0, we believe f chan is still a useful parameter indicating the predominance of the channeling shoulder effect over channeling if f chan < 0.

C. Average fraction of channeled ions
The fraction of ions that are channeled when implanted into a large polycrystalline sample with many randomly oriented grains is given by the average of the channeling fraction f chan over all directions.To exclude surface effects in our MD simulations, which use externally started ions with incidence angles up to glancing, we restrict the domain of averaging to the triangle formed by the [001], [111], and [101] direction.Since we used a regular (θ,ϕ) grid for the calculation of the channeling maps, the averaging is done with a weight factor sin θ to account for the dependence of the solid angle on the polar angle.The average channeling fraction f chan thus obtained may be interpreted as the probability for an ion to be channeled if the orientation of the crystal is unknown.
The average fractions of channeled ions f chan calculated from the MD results for the various ion-target combinations are listed in the sixth column of Table I.The corresponding values obtained from theory as described in Sec.II A are listed in column 7. The agreement is reasonable, with the theoretical values usually underestimating the MD foil results and overestimating the MD results derived from projected range simulations.This is most apparent in the 150 keV W in W simulations which have been performed in both a foil and a bulk sample.Channeling theory largely underestimates the MD results in the low-energy cases (5 keV Au in Au and 300 eV D in W).This is probably due to substantial scattering into channels of initially nonchanneled ions.

D. Fraction of channeling directions
Another way of analyzing the channeling maps is to ask how many directions are channeling directions.In other words, what is the probability that the nuclear energy deposition is substantially less or the projected range is substantially larger than the ones in a random target, if the orientation of the crystal is unknown?We define here as a channeling direction a direction where f chan > 0.1.The results for the fraction of channeling directions f chandir extracted from the MD simulations are shown in column 8 of Table I.To obtain theoretical values of f chandir , we have determined for each (θ,ϕ) pair of the MD simulations whether the direction is within the critical angle of any channel.The results are shown in the last column of Table I.As for the fraction of channeled ions, the agreement is reasonable although not perfect.Again the MD values tend to be higher for foils and lower for bulk samples, compare, in particular, the two values for the 150 keV W in W simulations.For the low-energy cases (5 keV Au in Au and 300 eV D in W) the MD results are much higher than the theoretical ones.Overall, it may be concluded from the data that the fraction of channeling directions is on the order of 30%, and may be up to 60% in particular cases (300 eV D in W).

E. Effect on range profile in polycrystalline samples
The large fraction of channeling directions also has significant practical implications with respect to irradiation of polycrystalline materials.It is a common assumption that when a polycrystalline material is implanted, one can calculate the range profile as for an amorphous material.However, our finding of a large fraction of channeling directions implies that also the range profile in a polycrystalline material with random grain orientations may differ very significantly.This is especially significant for metals, which are typically polycrystalline (and in the case of elemental metals, never amorphous) with grain sizes of tens to hundreds of micrometers [59], and also for transistor technology which involves polycrystalline silicon parts [60].
To illustrate the possible significance of this effect, we show in Fig. 14 range profiles calculated for the same ionmaterial combinations for purely amorphous material and for a crystalline material in a completely random direction.The former were obtained with a binary collision approximation (BCA) method [61] recently modified to use exactly the same interatomic potential and electronic stopping as MDRANGE [28], and the latter with MDRANGE, selecting averaging over random crystal orientations between 0 • and 90 • in ϕ and 0 • and 89 • in θ .MDRANGE calculations for completely random cells gave essentially identical profiles as the BCA simulations (see Appendix D) except for the low energy 300 eV D ions.These MDRANGE calculations correspond to experiments on ion ranges or penetration in a polycrystalline material with random grain orientations at the surface and grain size much larger than the ion range.
The results show that there is a major difference between the amorphous-and polycrystalline-material range profiles.The latter have in all cases a long 'tail' caused by channeling, and this leads to maximal penetration depths more than a factor of 2 deeper than the maximal depth of the range profile in amorphous material.This can have major practical implications, e.g., if one wishes to implant a polycrystalline thin film layer with no ions penetrating below the layer.

VI. CONCLUSIONS
In conclusion, we have compared MD and BCA simulations, experiments, and theory on ion channeling in several different crystal systems.Comparisons of the MD range simulations with experiments showed that the simulated angular channel widths agree with experiments within 10%, and the observed small differences may be just statistical in origin.
Our results show that channeling theory is very good at predicting which channels are active at a given energy, not only at high energies, where channeling theory has historically been applied, but also down into the sub-keV regime, as long as the ion travels a sufficient distance.Quantitative comparison of the fraction of channeled ions and of the fraction of channeling directions between theory and MD is reasonable but not perfect.This should not surprise since channeling theory has to make several assumptions that are not strictly fulfilled: It has been assumed that the ions may be classified into channeled and nonchanneled, each with their respective, well defined nuclear energy loss and projected range.This neglects different dechanneling probabilities and thus different mean energy deposition and ranges of ions moving in different channels.It also neglects that nonchanneled ions may interact more with the crystal than they would in a random target ('channeling shoulder' effect), to a degree that depends on ion species, target, and energy.Moreover, channeling theory is based on a given ion energy, while the projected range is influenced by all energies between the implant energy and zero.Since the minimum energy for channeling is different for different channels, this leads to different projected ranges for each channel.Finally, channeling theory does not consider the scattering of initially nonchanneled ions into channels, which leads to a broadening of the apparent channeling areas in the channeling maps, an effect most prominent at low energies.
Thus, channeling theory corroborates our MD simulations but does not provide all the quantitative details.One practical application of channeling theory would be to find incidence directions that avoid channeling, since it correctly predicts which channels are active.For instance, in all maps studied, the direction (θ = 20 • ,ϕ = 20 • ) is nonchanneling.Alternatively, channeling theory may be used to find wide channels for beneficial uses of channeling.
Channeling theory is independent of the atom masses, while an isotope effect has been discussed in the literature [55,56].We did not find an indication of an isotope effect comparing 10 keV H and T bombardment of Si.The isotope effect certainly deserves a more comprehensive investigation.
Finally, we highlight the observation from both simulations and theory that under typical ion irradiation conditions, a huge fraction (20%-60%) of the incidence directions are channeling directions.This implies that even when ion irradiations are carried out on polycrystalline samples with random surface orientations, channeling effects on the results cannot a priori be assumed to be negligible.
To test the possible dependence of the results on the interatomic potential, two different interatomic potentials were used, the ZBL universal [36] and DMol [62] ones.We considered two different test cases on opposite ends of the energy scale: 1.7 MeV Au ions into a 20 nm thick foil of Au in different crystal orientations (also presented at length in the main text for the ZBL potential results), and 15 keV B in Si.The inputs for the latter case were identical to those used in Ref. [40], i.e., also a 1.5 nm amorphous Si layer was used at the surface.For the interatomic potential testing, the nonlocal ZBL electronic stopping was used.
The results for 1.7 MeV Au ions in Au show that the range results are the same within the statistical uncertainties, see Fig. 15(a).Importantly, all dips in the angular profile (i.e., the The twist angle ϕ was chosen randomly between 0 and 360 degrees.For both potentials, the nonlocal ZBL electronic stopping was used.The lines are fits of a Gaussian profile to the data.channeling directions) are reproduced in the same manner by both potentials.Also the full channeling maps (not shown) were essentially identical for both potentials.
The results for the case of 15 keV B in Si near the [001] channel are shown in Fig. 15(b).The plot shows that while the absolute value of the ion ranges differ by about 3% for the two interatomic potentials, comparison of fits of a Gaussian lineshape to the two potentials show that the angular channel widths are very similar (for easier comparison of the width, the figure also shows the ZBL potential results scaled along the y axis to match the DMol potential level).The fullwidth-half-maximum (FWHM) of the angular channel widths obtained from the (unscaled) Gaussian fits were 6.04 • ± 0.08 • for the DMol potential and 5.92 • ± 0.08 • for the ZBL potential, i.e., agree within the statistical uncertainty.Thus these two comparisons indicate that while choice of repulsive interatomic potential may in some cases affect the mean ion ranges by a few % (as expected from previous works [63,64]), the angular width is the same within the statistical uncertainty.

APPENDIX B: TEST OF ELECTRONIC STOPPING EFFECT ON CHANNELING
To get the absolute shape of range profiles and R p fully correct in channeling directions, it is sometimes necessary to consider the reduction of electronic stopping in channels [11,[40][41][42]54].Since the electronic slowing down does not change the ion movement directions, this is not expected to affect strongly the angular channel widths.It can, of course, somewhat affect the fraction of nuclear energy deposition in a thin foil or the relative difference in ranges.To test the possible effect of different electronic stoppings, we ran the test case of 15 keV B implantation of Si near the [001] crystal channel with both the standard nonlocal ZBL96 electronic stopping and compared it with the "Puska-Echenique-Nieminen-Ritchie" electronic stopping based on a 3D electron density of the Si crystal [40].The latter model has been shown to agree almost perfectly with experimental range profiles.The inputs were identical to those used in Ref. [40], i.e., also a 1.5 nm amorphous Si layer was used at the surface.For the tests of the electronic stopping, the DMol interatomic potential was used.To determine the angular channel width, we simulated the projected mean ion range as a function of tilt angle θ off the [001] direction.The results shown in Fig. 16 show that while the absolute value of the ion ranges differ for the two stopping models, comparison of the ZBL and scaled Puska Gaussian fits show that the angular channel widths are very similar.The fullwidth-half-maximum (FWHM) of the angular channel widths obtained from the (unscaled) Gaussian fits were 6.04 • ± 0.08 • for the ZBL stopping and 5.74 • ± 0.08 • for the Puska stopping.Hence this comparison shows that while choice of electronic stopping can affect the mean ion ranges by ∼20% (as expected from previous works [44,45]), the angular width depends much weaker on the electronic stopping, in the tested case ∼5%.FIG.16.Comparison of effect of different electronic stopping models on the results, for the case of 15 keV B implantation of Si off the [001] channel.The models compared are the ZBL96 nonlocal electronic stopping [36] and the electronic stopping based on a 3D electron density from Ref. [40].The interatomic potential was for both cases the DMol one.The twist angle ϕ was chosen randomly between 0 and 360 degrees.The lines are fits of a Gaussian profile to the data.surface oxide can have an effect of reducing ion channeling and hence the mean range [65].Here we address the, until now unaddressed, issue whether inherent surface reconstructions, in the absence of an oxide, can similarly affect ion ranges.
To simulate the effect of surface reconstruction, we used for the initial atom layer Si (001) surface cells with the well known (2 × 1) surface reconstruction from Ref. [66].Prior to the MDRANGE calculations, the 54.3 Å sized cells were relaxed at 300 K using the PARCAS MD code [24,66].After this, ion ranges were simulated with the MDRANGE multilayer simulation mode, using the PARCAS MD cell as the top 54.3 Å and perfect Si crystal as the rest of the cell.The two different cells were constructed to have matching lattice plane coordinates in the lateral directions.Since the PARCAS MD cell already had thermal displacements, additional thermal displacements were not added for it by MDRANGE.
The results on the mean range for Si ions of varying energy are shown in Fig. 17(a), comparing cells with and without the top layer surface reconstruction.For the perfect [001] channeling irradiation, the ranges with the surface reconstruction are slightly ( 2 nm) shorter than without it.This shows that the surface reconstruction can somewhat enhance the dechanneling probability at the surface, leading to shorter mean ranges in the channel.The results also show that in the nonchanneling condition (θ = 20 • , ϕ = 20 • ), the results are slightly larger with the surface reconstruction within the statistical uncertainty.A similar effect has been previously observed due to surface oxides [67] and explained to be due to scattering of ions into a channeling direction [68].
We also examined whether the surface reconstruction can affect the angular width of channels.The results in Fig. 17(b) show that although the ranges differ slightly near the center of the channel [consistent with Fig. 17 FWHM of the unscaled distribution were 7.37 • ± 0.07 • and 7.23 • ± 0.06 • , i.e., the same within the statistical uncertainty.From this analysis we conclude that surface reconstructions can affect the mean ion ranges slightly but have no appreciable effect on angular channel widths.generated by placing atom positions in random in 3D space at the experimental atomic density of Si, however with the constraint that no atoms should be closer than r min = 2.1 Å from each other.As usual in MDRANGE runs, these atoms were still given random thermal displacements corresponding to 300 K according to the Debye model, with 10 different sets of random displacements [69].However, the initial atom coordinates were always the same.
The resulting channeling maps (not shown) showed that for cell sizes of 30 Å and 40 Å, the mean ion ranges were completely evenly distributed for each (θ,ϕ) combination.The results in Fig. 18(a) selected for two values of ϕ show that the 30 Å and 40 Å cells appear almost identical, while the 20 Å cell shows statistically significant variations above and below the values for the larger cells.This shows that the size of the 20 Å cell is too small to represent a true amorphous material (that by definition has no long-range order).Careful inspection of the data for different ϕ values showed, however, that for the 30 Å cell the θ = 0 • ranges are systematically slightly (about 5%) larger than for θ = 2 • or θ = 4 • , indicating that even the 30 Å cell is not sufficiently large to fully avoid 'channeling' due to finite cell sizes.The 40 Å cell shows no such effect and hence clearly is sufficiently large.
The slight increase of the mean ranges with θ visible in all cases is simply due to the probability of ion reflection near the start of the penetration which increases for glancing angles and leads to a decreased probability of short ion ranges.On the other hand, the decrease in range above θ = 80 • is related to an increased probability of also ions near the end of range leaving the sample (the simulations showed that for θ > 80 • more than half of the ions are reflected), leading the mean range being dominated by a few ions that were scattered strongly into the material and hence had shorter ranges.
The value of r min = 2.1 Å was chosen because this is about the minimum interatomic separation between atoms in Si at room temperature (the nearest-neighbor separation in crystalline Si is 2.35 Å, but stress in amorphous materials reduces the minimum distance between atoms).We tested whether this choice of r min affects the results by using different values and simulating the mean range for 10 keV Si irradiation in Si, for 40 Å cells.The results were also compared with range profiles simulated with a realistic a-Si structure created previously by molecular dynamics [70].The results in the inset of Fig. 18(a) show that the mean ranges are identical within the statistical uncertainty for all r min values and also the MD cell.This confirms that the approach of generating random atom coordinates is appropriate for performing ion range calculations in amorphous materials.
As an alternative test whether the random atom coordinate simulation cell approach is suitable to mimic amorphous material, we also compared the MDRANGE range distribution for a 40 Å cell with those from a BCA range calculation.The latter generates each atom position separately in a Monte Carlo algorithm [28,61,71].The results in Fig. 18(b) show that both approaches give identical results within the statistical fluctuations.
From this model study, we conclude that to mimic a truly amorphous material with individual disordered simulation cells, one should use cells of at least 40 Å side length.For much higher ion energies or lighter ions, for which channeling is more pronounced, even larger cell sizes could be needed.

FIG. 1 .
FIG. 1. (a)Relative degree of channeling from experiments and simulations for 1.2 MeV He ions in single-crystalline Si near a 100 crystal channel.In the experiments, this ratio α is the fraction of channeled to random stopping power, while in the simulations it is the difference between the mean ion range in the channel to that in the offchanneling direction, rescaled in the y axis to match the experimental data.The y axis of the simulated data on mean ranges has been linearly rescaled to be comparable to the experimental stopping power ratio.The lines are fits of a Gaussian profile to the data.(b) Same for 1.2 MeV He in Si near a 110 crystal channel.(c) Full-width half maxima (FWHM) of the experimental and simulated channeling results compared to each other for all the seven cases studied.The FIG. 2. (Continued) error bars are 2σ errors of the FWHM valuesobtained by fitting a Gaussian profile to the channeling data.Note that the comparison in (c) does not involve any fitting or rescaling of the simulated data to the experimental ones.The experimental data is scanned in from Refs.[47,48].Following the experimental papers, the data in (a) and (b) is presented reflected around 0 degrees.

FIG. 3 .
FIG. 3. Channeling results for various ions in Au foils.The three different cases are (a) 80 keV Xe, (b) 30 keV Ga, and (c) 10 keV H irradiation of 20 nm, 10 nm, and 20 nm Au foils, respectively.The colors show the nuclear energy deposition by primary knock-on atoms to the Au foil as calculated by the MDRANGE range calculation software.The color scale shown to the right gives the energy deposition in units of eV.

FIG. 5 .
FIG. 5. Channeling results for (a) 150 keV W in W and (b) 300 eV D in W. The colors show the mean range of ions as calculated by the MDRANGE range calculation software.The color scale shown to the right gives the mean range in units of Å.The white areas close to θ = 90 • correspond to directions where all ions were reflected and hence ion ranges could not be obtained.

FIG. 6 .
FIG. 6. Channeling results for various ions and energies in Si.The colors show the mean range of ions as calculated by the MDRANGE range calculation software.The three different cases are (a) 10 keV Xe, (b) 10 keV Si, and (c) 10 keV H.The white areas close to θ = 90 • correspond to directions where all ions were reflected and hence ion ranges could not be obtained.
FIG. 8. (a)Minimum approach distance and penetration coefficient for 10 keV H incident on a (001)-Si surface as a function of incidence angle.The azimuthal incidence angle is 10 • with respect to the [010] direction.MD results of the minimum approach distance are compared to the predictions of channeling theory.The critical angle is obtained by equating the minimum approach distance with the critical approach distance.Immediately below the critical angle some ions may penetrate the surface atomic plane.(b) Average exit angle of the MD simulations as a function of incidence angle.Immediately below the critical angle the mean polar exit angle starts to deviate from the incidence angle, and the standard deviation increases sharply.
[001], [101], and [111] direction.This part of the theoretical channeling map is shown in the lower right parts of Figs.9-11.The MD energy deposition results are shown for the same implant conditions in the upper left part, in the triangle formed by [001], [111], and [011].

The 5 keV
Au in the Au case [Fig.10(c)] behaves slightly differently: [101], [001], and [112] channeling is correctly predicted by the theory as well as the absence of [111] channeling.The [101] channel, however, is extremely smeared in the MD data.

FIG. 9 .
FIG. 9. Comparison of MD simulation results with theory for high-energy ions: (a) 1.7 MeV Au in a 20 nm thick Au foil, (b) 150 keV W in a 10 nm thick W foil, (c) 150 keV W in bulk W. In (a) and (b) the nuclear energy deposition is shown, while in (c) the mean projected range is depicted.

FIG. 10 .
FIG. 10.Comparison of MD simulation results with theory for medium energy ions: (a) 10 keV H in Si, (b) 10 keV Si in Si, (c) 5 keV Au in Au.map.It may be concluded that channeling theory remains valid down to sub-keV energies as long as the ion ranges are sufficiently large.

FIG. 11 .
FIG. 11.Comparison of MD simulation results with theory for low-energy ions: (a) 500 eV Au in Au, (b) 300 eV D in W.

FIG. 12 .
FIG.12.Mean ranges of 10 keV H and T in Si as a function of polar angle θ for azimuthal angles of ϕ = 0 • and ϕ = 12 • , as calculated by MD simulations.

4 R
FIG. 13. (a)Statistics of the distribution of the obtained nuclear energy deposition F Dn over all crystal directions.The x axis is normalized relative to the energy deposition in an amorphous sample for the same ion-energy combinations.To allow for comparison of the different cases, the results are normalized.(b) Same for mean ion ranges R mean in Si.The x axis is normalized relative to the mean range in an amorphous sample for the same ion-energy combinations.To allow for comparison of the different cases, the maxima of the distributions are scaled to 1 on the y axis.

FIG. 14 .
FIG.14.Range profiles in Si, W, and Au calculated for amorphous material in a BCA approach and by MD in polycrystalline material of grain size much larger than the ion range.

FIG. 15 .
FIG. 15. (a)Results on nuclear energy deposition for 1.7 MeV Au ions on a 20 nm Au foil for ϕ = 0 • and ϕ = 24 • as a function of the polar angle θ, calculated with two different interatomic potentials, the Universal ZBL and DFT DMol potentials.The dips in the curve show the channeling directions, since in these much less energy is deposited into the foil than in nonchanneling ones.The ZBL curve is calculated with a statistic of 5000 ion trajectories per (θ,ϕ) combination, and the DMol one for 1000 ions.Since the error bars are 1σ errors of the mean, the results with the two potentials are identical within the statistical uncertainty.(b) Results for mean range of 15 keV B ion on Si off the [001] channel calculated with two different interatomic potentials, the Universal ZBL and DFT DMol interatomic potentials.The twist angle ϕ was chosen randomly between 0 and 360 degrees.For both potentials, the nonlocal ZBL electronic stopping was used.The lines are fits of a Gaussian profile to the data.

FIG. 17 .
FIG. 17. Test of effect of the (2x1) Si (001) surface reconstruction on mean ranges R. The error bars are 1σ errors of the mean.(a) The upper part shows the mean ion range R as a function of incoming ion energy, and the lower part the relative difference between the ranges for the unreconstructed and reconstructed surface cases.(b) Effect of surface reconstruction on the angular dependence of channeling for 10 keV Si ions.The twist angle ϕ was chosen randomly between 0 and 360 degrees.