Enhanced phenomenological models of ion channeling contribution to doping profiles in crystals

New phenomenological models are proposed to describe the effect of an ordered lattice structure of crystalline targets on the as-implanted doping profiles of low-energy heavy ions. The models account for the channeling kinetics and clarify the effect of bi-directional transitions of ions between random-like and channeled modes of motion on the target depth dependencies of dopant concentration. They also incorporate a simple model of the target radiation damaging effect on doping profiles. The presented results of model validation against the experimental and Monte Carlo computer simulation data and the comparative analysis of the capabilities of the proposed and the existing models show that the application of a more physically grounded approach allows us to improve the quality of doping profile description. The theoretical models developed are useful for obtaining physical parameters of low-energy ion channeling kinetics from the experimental data.


Introduction
Ion implantation into solids is a mainstream technology of the modern semiconductor industry and a state-of-the-art method to modify mechanical, electrical, optical and other properties of solid surfaces and interfaces by a controlled introduction of definite concentrations of appropriate dopants (mainly heavy ions) into proper spatial regions of irradiated targets [1,2].The express prediction of spatial distributions of implanted ions (and, in particular, the calculation of the dopant concentration depth dependencies known as "doping profiles") is a key starting point in the development and optimization of these technologies, and a valuable constituent of the physically based semiconductor technology computer-aided design (TCAD) systems [3,4].
Except an implanter beam characterization issues, various problems relevant to solid state physics arise when predicting doping profiles theoretically.A specific problem that this paper deals with is the effect of the ordered atomic structure of irradiated targets on the penetration and stopping of heavy ions of low, keV-to-MeV, energies inherent to industrial implantation processes.
This problem has got a straightforward zero-order solution in the approximation of structureless "truly amorphous" targets having a random arrangement of atoms.Then it is resolved by numerical solution of the Boltzmann type kinetic equation of ion transport using the method of moments.In the most topical cases, doping profiles in an amorphous solid can be acceptably characterized by first four moments [5].So, one can parameterize the numerically calculated profiles by analytical distribution functions of Pearson type IV [6,7] well studied in mathematical statistics.
However, ion implantation into long-range ordered single or polycrystalline targets is much more complicated due to the decisive role of directional effects [8] commonly attributed to ion channeling [9,10].As compared to those in highly amorphized solids, doping profiles in crystals manifest long-range phenomena such as the so-called "channeling tails" [11,12].They are very essential for modern applications that frequently require a nanoscale spatial resolution and doping precision.
No general master kinetic equation formalism is currently developed to describe ion transport for this case of highly correlated atomic arrangement.Thus, two alternative methods are used to calculate as-implanted doping profiles in crystals.They are the atomistic computer simulation, and the application of empirical models.
The Monte Carlo binary collisions or molecular dynamics (MD) computer simulation is the most powerful method of calculation of 3D spatial distributions of implanted ions that adequately accounts for the lattice structure of targets.Numerous Monte Carlo codes [13][14][15][16][17] have been successfully validated against experimental Secondary Ion Mass Spectrometry (SIMS) data, and are incorporated into commercial TCAD simulators of industrial processes.An atomistic computer simulation is a perfect predictive tool.However, it is rather resource-intensive, and so it is hard to be applied routinely for the analysis of large arrays of experimental data.For these purposes, analytical models are generally favorable in providing the theoretical paradigm of compact description and intentional interpretation of experimental and simulation results.
Empirical models are designed to fit the implant concentration profiles C(z) by multivariate analytical functions of depth z in a target.The implantation conditions dependency of model parameters is obtained from the nonlinear least squares (NLSq) fit procedure.Of course, the information on the lattice driven directional effects is already contained in the underlying experimental or computer simulation data.Therefore, it is described appropriately if the model is flexible enough.Advanced empirical models are intended to provide interpolation and extrapolation of parameters in order to predict profiles at irradiation conditions different from those used for primary data fit.But one should emphasize that empirical models are not derived from first principles.Their quality and physical adequacy can only be evaluated by quantitative estimation of NLSq best-fit residuals and by numerical stability of inteprolation/extrapolation procedures.
Except pure mathematical models that fit the profiles by the properly chosen set of orthogonal (e.g.Legendre [18]) polynomials, the heuristic empirical models [19][20][21] are widely used for practical calculations.To some degree, they reflect the physics of ion transport and implantation into crystals.The basic assumption of the majority of these models is the separation of a profile onto "random-like" (w r ) and "channeled" (w ch ) fractions weighted by an adjustable parameter q: where φ is an ion fluence, p r,ch are vectors of fitting parameters.The models differ by a choice of w-functions.For example, within the scope of the most accepted "dual Pearson" model [20] both w r and w ch are Pearson-type distributions [22] having independent sets of parameters p r and p ch .The existing empirical models manifest a reasonable capability of fitting the doping profiles for a limited set of implantation setups (ion beam energy and direction, implantation dose).However, these models have certain drawbacks due to a rough treatment of directional effects.They presume (see (1)) that "random-like" and "channeled" fractions of a profile are completely independent of one another.Physically, it means that random and channeled fractions of an implanted beam are non-interacting, i.e. no particle transitions occur between them.As a result, non-physical dependencies of model parameters on irradiation conditions can frequently appear.It makes their independent validation difficult, impedes interpolation and extrapolation beyond the calibration domain, and considerably decreases the accuracy and predictive ability of conventional models.
One can expect that the quality of multivariate models can be improved if the kinetics of ion transport in crystals is a priory assumed as the basis of their formulation.It permits to transform empirical models into phenomenological ones with model parameters having an immediate physical meaning and, at least in principle, permitting independent estimations.Such an advanced physically based approach is developed in this paper with due account taken of the physics of directional effects relating to the kinetics of directed and random-like motion of ions in ordered lattices.A set of phenomenological models is developed (including a simple model of high-dose effects on implantation profiles) that result in the analytical formulae suitable for fitting of experimental doping profiles for a wide range of implantation conditions.It is shown that the fit accuracy is better (or at least not worse) than that of the state-of-the-art "dual Pearson" model at a lesser number of fitting parameters.It is also shown that the dependencies of these parameters on implantation conditions are physically reasonable and fairly consistent with the well-established concepts of the theory of directional effects in crystals.

Problem setup
The initial separation of ions onto random and channeled fractions at the entrance to a crystal does not uniquely determine the future of each particle.Channeling is not completely stable owing to various stochastic factors (thermal vibrations of atoms, lattice defects, etc.) that cause dechanneling [8].On the other hand, there exists a possibility of initially non-channeled ions to be captured into the channeling mode in the bulk of a crystal (the volume capture effect), and even the finite probability of dechanneled particle to return into a channel (the rechanneling effect).It can be followed by subsequent dechanneling, and so on.As a result, as-implanted doping profiles in crystalline targets are affected by the kinetics of (generically multiple) transitions of the moving implants between different modes of "aligned" and "random-like" motion.
To identify the subject matter of our problem, we have performed the computer simulation of 15 keV Boron implantation into the [100] axial channel aligned Silicon crystal.The MD code MICKSER [23] has been applied to the calculation of the as-implanted doping profile and to its separation into the specific contributions of ions implanted after a different number of transitions from/into the [100] axial channeling mode of motion.At 3D MD simulation, transition events were detected by calculation of transversal energy E ⊥ of ions and by its comparison with the critical transversal energy E ⊥c of [100]Si axial channeling (see reference [23] for algorithmic details).
The simulation results shown in figure 1 argue evidently that the doping profile definitely consists of three (but not of two, like in (1)) fractions.The "random" fraction dominates at small depths z < 30 nm near the crystal surface and is formed by ions that had never been captured into the [100] channel.The "channeled" fraction forms a long-range edge of a profile (z > 250 nm) due to the ions captured into the stable channeling at the crystal surface.For these fractions, the shapes of partial doping profiles are qualitatively similar to those of the asymmetric Pearson type IV distributions that are generally accepted in the physics of implantation and are used as basis functions w r,ch of the standard "dual Pearson" model [20].
However, in a wide intermediate range of implantation depths (50 ÷ 220 nm) the profile predominantly consists of the implants that have experienced at least one transition from one mode of motion to another (a dechanneling from, or a bulk capture into the [100] channel).This fraction is conceptually ignored by the existing empirical models.The significance of the dechanneling related fraction of a doping profile has been confirmed experimentally [24] for the case of highenergy (0.25÷1 MeV) implantation of Phosphorus ions into crystalline Silicon at high implantation doses when the dechanneling caused by an ion beam induced lattice defects dominates.However, our simulation argues that this is also significant in an undamaged crystal, and deserves special consideration.
The goal of the current study is to develop and to validate the phenomenological model that would take into account the existence of such a fraction and describe explicitly the implantation of ions that have experienced transitions between random-like and channeled modes of motion in crystals.

General model
It is well-known that the rate of the ion slowing-down is essentially different in random-like and channeling modes of motion in crystals [8,25].In reference [26], a phenomenological model of ion channeling has been built taking into account kinetics of transitions between these modes.For the sake of simplicity, this model assumes that ions can come to a stop only within the nonchanneled mode of motion having non-reduced rate of energy losses.Therefore, it is not well fit for a description of the case of a highly collimated implanter beam directed parallel to the channel axis when a considerable fraction of ions is captured into the long-range stable channeling mode and one cannot neglect their energy losses before the dechanneling event.The model proposed in reference [26] is reformulated herein below by taking into account the possibility of particles' coming to a stop in both modes of motion.
The occurrence of particle transitions between the modes of motion implies that up to the moment of a stop at a given depth, z ions can move in each of these modes over certain parts of their total path length.The trajectory of each ion is characterized by a specific chain of transitions that hereinafter is referred to as "transition history".
To build an as-implanted doping profile, one has to specify the ion's probability of coming to a stop in a given mode of motion provided that it was moving a definite part of the path in another mode.For this purpose, the "Numerical Range Scaling" (NRS) method is applicable.It was proposed in reference [27] for a description of doping profiles in layered structures.However, one can draw an analogy of each mode of motion with the layer of a material having a specified distribution of the particle probability to come to a stop inside the layer.
To shorten the notation of the formulae, let us indicate the quantities relevant to different modes of motion by indices "u" and "v".Both of them can possess the values "r" for random-like and "ch" for channeled mode (obviously, u = v).
According to the NRS approach, the conditional probability distribution function (PDF) p u (y u |y v ) of a trajectory stop at the path length y u in the u th mode of motion under the condition of passing the path length y v in the v th mode can be written in the form: where R (u,v) p and W u,v (y; p u,v ) are the projective ranges and the distribution functions of the particle probabilities to come to a stop in the u th and v th modes (the stop functions), respectively.The stop functions W u,v (y; p u,v ) are normalized as follows: ; p u,v )dy = 1.The normalization of p u (y u |y v ) for u th mode is determined by the non-stop probability in v th mode and the scaled range p ) that, according to NRS method assumption, depends on the total path length y v spent in the alternative mode of motion.
The partial doping profile C u (z) of the ions that started up in the u th state of motion at depth z = 0 can be obtained by integrating the conditional PDF (2) over all possible z-ended paths taking into account the statistical weight of a specific transition history each path is associated with.
To specify these weights, let us introduce new phenomenological parameters of the model, the characteristic lengths of transitions [26].Let us assume that the probability (per unit of path length) of transition from u th to v th mode of motion is described by the constant length R u = Σ −1 u→v that is inversely proportional to the macroscopic cross-sections Σ u→v of the correspondent transition.
One should note that the cross-section of the first transition is substantially dependent on the initial conditions of irradiation while the subsequent cross-sections are mainly governed by the physical mechanisms of dechanneling and rechanneling.Thus, for each u = (r, ch) one has to discriminate the primary and secondary characteristic lengths, R Then, the partial profile C u (z) can be represented as a series expansion in terms of the number n of particle transitions between the modes of motion: where The structure of expressions (3-5) follows from the straightforward probabilistic considerations.The statistical weight of a transition history having n transitions at depths

and those of an ion to keep the specific mode of motion at depth intervals ∆
Each of the latter ones is of exponential form exp (−∆ k /R k ) since R k are assumed to be depth independent in this model.Subject to the number of preceding transitions, R k can possess the values v for odd k, and R (s) u for even k.By multiple integration over all x k followed by summation over all n ⊂ [0, ∞) one takes into account all possible transition histories that correspond to different random values of total path lengths y u and y v that an ion has traveled in each mode of motion.By definition (2), they are distributed according to the conditional PDF p u (y u |y v ).The formulae (4) and ( 5) differ due to the fact that the mode of motion at depth z coincides with the initial mode of motion for an even number of transitions, contrary to the odd n.
Since in this model the transition probabilities are completely determined by three lengths v , the integrands of expressions (4-5) effectively depend only on three variables, namely the depth of the first transition and the subsequent total path lengths spent in each mode of motion after the first transition.Then, changing the variables and the order of integration one can reduce the infinite series (3) for C u (z) to the closed form that contains only a double quadrature: where the functions B 0,1 (a, b) are expressed in terms of modified Bessel functions I 0,1 : Let us denote the probability of ion capture into the "channeled" mode of motion at the entrance to a crystal by q.Then, the doping profile C(z) for a given ion fluence φ has the following form: where C ch (z) and C r (z) are described by equation ( 6) with substitution of formula (2) and with indices (u, v) replaced by (ch, r) and (r, ch), respectively.Evidently, the formula (9) transforms into the basic equation ( 1) of the conventional "dual Pearson" model if there are no transitions (all characteristic lengths tend to infinity), and W ch (z) and W r (z) are chosen as the Pearson-type distribution functions.

A simplified model with dechanneling
As against the conventional "dual Pearson" model, the full-featured model proposed in section 3 introduces four additional fitting parameters R (p,s) ch,r , the characteristic lengths of the correspondent transitions (including volume capture and rechanneling).However, in the first approximation one can neglect the volume capture of ions into the channeling mode, thus restricting ourselves to the consideration of dechanneling of initially channeled ions.It can be authorized by a substantial suppression of the capture of particles into the channeling mode of motion due to blocking [28] of scattering of initially unchanneled (or dechanneled) ions into the directions of open channels of crystal lattice.
The correspondingly simplified model of the doping profile C(z) can be readily derived from the general model of section 3. Let's note that in this approximation the only finite characteristic length of transitions is the primary dechanneling length R (p) ch ≡ R ch since the neglect of volume capture and rechanneling corresponds to the limiting case of infinite R The proceeding to these limits in equations ( 3)-( 5) results in C (n) u (z) → 0 for all n 2.Moreover, the infinite series of equation ( 3) reduces in these limits to its first two terms for initially channeled ions (u = ch), and to the single first term for initially non-channeled ions (u = r).By substitution of these remaining terms into equation ( 9) one obtains the following expression for the total doping profile that is now dependent on dechanneling only: where the NRS conditional PDF p r is described by formula (2) with indices (u, v) → (r, ch).Obviously, equation ( 10) also follows from the correspondent limiting cases of partial profiles C ch , C r described by equation ( 6) of the general model, at their substitution into formula (9).
Three terms of (10) directly correspond to three fractions of a doping profile (see figure 1).The first term describes the implantation of initially random fraction of a beam.When volume capture is neglected, it coincides with the stop function W r z − R (r) p ; p r scaled with the initial probability 1 − q of the random mode of motion.The second term corresponds to the implantation of perfectly channeled ions that come to a stop inside a channel.The model reduces to the dechanneling independent two-fractional "dual Pearson" form (1) if the phenomenological dechanneling length ch → ∞.However, at finite R ch the third term of formula (10) describes the "mixed" fraction of a profile that appears primarily due to single dechanneling of ions.
In our models, dechanneling is considered to be independent of the ions coming to a stop described explicitly by PDFs p r,ch (2).The exponential function qR −1 ch exp(−x/R ch ) in the integrand of the third term of ( 10) is the per unit path length dechanneling rate at depth x.This rate has been calculated assuming that ions cannot stop in a channeled mode of motion and using the approximation of constant R ch derived from the depth-averaged macroscopic cross-section of dechanneling.Without going beyond the single-dechanneling model, one can get rid of the latter approximation by substituting the exponential dechanneling rate by a generalized function |dP ch (x)/dx| where P ch (z) is the dechanneling function (the probability of finding a channeled particle at depth z) calculated without regard to the possibility of coming to a stop in a channel, and P ch (0) ≡ q.For a depth independent R ch , P ch (z) = q • exp(−z/R ch ) which immediately leads to (10).For a given P ch (z), one can rewrite (10) as follows: To fit the doping profiles, this formula makes it possible to apply different representations of P ch (z) including those obtained within the framework of microscopic models of channeling kinetics [29,23].

A model of high-dose doping profiles
At sufficiently high implantation doses (i.e. at ion fluence φ > 10 13 ÷ 10 14 cm −2 ) the dechanneling caused by the irradiation produced crystal lattice defects becomes considerable and has to be taken into account in order to describe the doping profile C(z).Due to spatial dependence of the concentration of defects n d (z, φ) this is just the case when the depth dependence of the damage-induced dechanneling should be taken into account.The averaged dechanneling function P ch (z) at given fluence φ of implanted particles can be calculated as follows: where σ d is the microscopic cross-section of the damage-induced dechanneling [30].Generally, the defect concentration profile n d correlates to the damage-dependent distribution of implanted ions.It results in the nonlinearity of the problem of calculation of high-dose dopant distributions.However, in a rough linear approximation, one can assume that the damage profile n d (z, φ) follows the doping profile C 0 (z) (see (10)) in an undamaged crystal, i.e. n d (z, φ) ∝ C 0 (z).
Within this approximation, one can introduce an effective cross-section Σ d of dechanneling that takes into account both the microscopic dechanneling cross-section σ d and the cascade function of the defect production by implanted ions.This new phenomenological fitting parameter permits to replace the product n d (x, φ) • σ d in (12) by the product containing C 0 (z): Consequently, to calculate the doping profile C(z) in a damaged crystal, one can apply the expression (11) with the following representation of P ch (z) that takes into account the dechanneling caused by defects:

Range distributions for random-like and channeled modes
Different representations having different number of parameters p can be used to specifically define the range distribution functions W r (y; p r ) and W ch (y; p ch ) of random-like and channeled modes of motion of implanting ions.We shall operate with W r,ch derived from the Pearson type IV PDF.Its applicability is confirmed qualitatively in figure 1, and is also fairly consistent with the standard "dual Pearson" model.
The Pearson type IV distribution is unambiguously determined by its first four moments, namely range R, straggling σ, asymmetry γ, and excess β.It has the following form [31]: where To fit the physically valid unimodal profiles, the third and fourth moments of this distribution must satisfy the conditions γ 2 < 32 and β > β 0 (γ) where The fourth moment β mainly determines the shape of a profile far from its maximum [7] where experimental profiles are typically known with large data errors.Its value only weakly affects the profile reconstruction accuracy.Thus, one can express the excess β by a function of the third moment γ using, e.g., the relation β = β 0 (γ) + γ 2 that was specifically proposed in reference [7] with regard to the profiles in amorphous targets.Such a parametrization effectively reduces the number of the Pearson type range distribution parameters to three, and is knowingly used in all calculations presented below in this paper.However, it is inappropriate for fitting the channeling profiles with the standard 9-parameter "dual Pearson" model (1) due to the lack of its flexibility in describing the dechanneling related part of a profile at intermediate depths (see figure 1).
Within the scope of our models, the channeling kinetics being explicitly taken into account lowers the requirements put forward to the flexibility of stop functions W r,ch (y; p r,ch ) and permits to apply, except for the above-mentioned reduced Pearson distribution, simpler PDFs such as the 3-parameter joined half-gaussian asymmetric distribution [19], or even the simplest 2-parameter normal PDF.On the other hand, the models introduce additional kinetic parameters.One can summarize that in the case of reduced 3-parameter Pearson type stop functions, the general model of section 3 has 11 parameters, the simplified model of section 4 depends on 8 parameters, and the high-dose model of section 5 is 9-parameter.Thus, the total number of fitting parameters in the computationally efficient models to be validated quantitatively in the next section does not exceed that of the "dual Pearson" model.

Model validation
To validate our models, a large array of experimental and simulation data available in publications [32][33][34][35][36][37][38][39][40] has been used.Their model description was carried out by means of the Levenberg-Marquardt method of NLSq fitting.The fit quality was estimated according to the reached value of the reduced chi-square (χ 2 ).The selected results of model validation are presented and discussed below with an emphasis made on the physical interpretation of the obtained model parameters.

Energy dependence of ion ranges
The results of NLSq fitting of experimental SIMS profiles of implantation of Boron ions into [100] axial channel of Silicon single crystals are presented in figure 2 for a wide range of ion energies spreading from 20 keV up to 2.8 MeV.Experimental SIMS data picked up from various sources [32][33][34][35][36][37][38] were fitted using both the 8-parameter model developed in section 4, and the standard 9-parameter "dual Pearson" model.It is clear that both models are qualitatively applicable to fit all the data.However, the dechanneling being taken into account in our model permits to achieve 2 ÷ 10 times smaller values of reduced χ 2 as compared with those of the standard model, and therefore works better at a lesser number of parameters.
The energy E dependencies of the ranges R (r,ch) p (E) at random-like and channeled modes of ion motion are depicted in figure 2(b).They have been obtained as parameters of NLSq best fits of figure 2(a) profiles using both kinds of models.
To evaluate these phenomenological adjustable parameters physically, they are compared in figure 2(b) with the results of a priori calculations of ranges.The values of "random" ranges were immediately taken from the state-of-the-art database of the SRIM2006 [41] Monte-Carlo code.The estimation of "channeled" ranges R (ch) p (E) were carried out in a continuous slowing-down approximation (CSDA) according to the formula where S ch (E) denotes the effective specific energy losses (the stopping power |dE/dz|) at channeling of ions.Only electronic S ch (E) (also borrowed from the SRIM2006 database) were incorporated into these calculations since for the stable channeling of sufficiently light Boron ions they are known to significantly exceed the nuclear ones due to a strong flux-peaking effect [8,29].One can see that the dependencies R (r,ch) p (E) obtained within the scope of the proposed model agree much better with the results of physically grounded estimations than the range predictions that have arisen from the standard "dual Pearson" model.The latter one systematically overestimates the ranges relevant to random-like fraction of an implanted beam.This is quite explicable since the disregard of dechanneling impels the standard model to shift the partial "random" profile w r much deeper as compared with that predicted by the ion stopping calculations for amorphous materials.
It is also significant that at low (E < 40 keV) energies our model predicts the deviation of the "random" range R (r) p (E) from the SRIM2006 code calculated range in a structureless material without whatever lattice ordering.This is due to the well known directional effect of strong blocking collisions with atoms of aligned atomic chains, and is frequently observed in experiments with crystalline targets as a shift of the near-surface peak of a profile toward smaller depths as compared with the peak in amorphous targets located at a depth close to the correspondent projective range R p .This effect cannot be reproduced by "dual Pearson" model that conceptually ignores the profile variations specific to dechanneling and blocking.
A reasonable agreement of the results of both models obtained for a "channeled" range at high energies is caused by the fact that R (ch) p is determined by the most stably channeled (the so-called hyperchanneled [42]) ions having maximal ranges and coming to a stop in a channeling mode.Thus, this range is practically dechanneling-independent.

Global fitting of profile families with parameters sharing
An important aspect of the analysis of experimentally measured doping profiles consists in their dependence on the implantation conditions such as the implanter beam orientation and angular spread, the thickness of amorphous oxide layer that covers the surface of target wafer, and so on.
From the physical point of view, it is quite probable that each of these experimental factors mainly affects certain parameters of the phenomenological model while its effect on other parameters is not so essential.These a priory considerations can be taken into account by the procedure of global NLSq fitting of all profiles that belong to a specific family and differ by the value of a certain characteristic variable of implantation process.At this procedure, certain parameters of a fitting model can be treated as the shared ones among all profiles of a family while other parameters remain to be the specific profile dependent.
The model can be validated by means of intentional physical interpretation of the dependency of the profile-dependent adjustable parameters on the implantation process characteristics.
The examples of applying such an approach are illustrated in figures 3 and 4. Six profiles of figure 3(a) were simulated in reference [39], using the binary collision method, by the Monte Carlo module of the ATHENA TCAD package [4] for the cases of Si target coverage with amorphous SiO 2 films that differed by a variable thickness t ranging from that of native oxide (≈ 1 nm) up to 80 nm.
The surface coverage with amorphous layer results in extra scattering of ions which hinders their surface capture into the channeling mode and spreads their angular and transversal energy distributions.This factor strongly affects the directionally dependent model parameters, the surface capture probability q and the dechanneling length R ch .As regards the parameters of range distributions, one should note that, by definition, they are mainly determined by the bulk stopping of ions having a definite mode of motion (random-like or channeled), and only weakly depend on the surface-localized target properties.Even at channeling, its kinetic theory argues that at sufficiently large depths the specific energy losses of channeled particles become independent of their startup conditions at a crystal surface due to the establishment of quasi-stationary statistical equilibrium transversal energy distribution of ions.
Therefore, the range-relating parameters depend mainly on the ion energy that only weakly changes within the thin layers of surface contaminations.Consequently, one can consider them to be equal for all profiles of figure 3(a), and they can be set as shared parameters of the global NLSq fit procedure.
At such a procedure, we have found the following best-fit values of these shared parameters: R (r) p = 133.13nm, σ r = 51.52 nm, γ r = 0.17 for the stop function W r of the random-like motion, and R (ch) p = 538.71nm, σ ch = 165.65 nm, γ ch = −4.03for W ch .As it is clear from figure 3(a), they provide a very reasonable fit of all profiles at any oxide film thickness t available.
The global fit derived dependencies q(t) and R ch (t) are depicted in figure 3(b).They are smooth monotonic functions of the oxide film thickness and are quite interpretable physically within the framework of the directional effects theory.Really, both the surface capture probability and the dechanneling length have to decrease as t increases.Both effects are due to angular spreading of a primary beam at scattering inside the oxide layer.The kinetic theory of channeling states that the dechanneling rate of initially well-channeled particles increases as their angular and transversal energy distributions broaden.However, at near-barrier states of transversal motion that are predominately populated at large t, the rate and the effective length of dechanneling become prac-tically independent of the details of initial angular distribution of ions, and therefore independent of the oxide film thickness.
Smooth behaviour of functions q(t) and R ch (t) permits to suppose profiles to have a continuous dependence on the fitting parameters.Then, the parameters can be interpolated for the case of intermediate values of t to predict the shape of the correspondent doping profile.
Similarly, the data of figure 4 illustrate the results of the global fit of directional dependencies of doping profiles at off-axis implantation of 12 keV B ions into Si single crystal target.At experiment [43], the angular scan of implanter beam by tilt angle θ with respect to (001) crystal facet was carried out in a plane that contains, at tilt angle θ = 35.26• , the axial channel [112].
The global fit of six profiles of figure 4(a) performed within the same assumptions as for figure 3 also has demonstrated a good agreement of the model with experiment.The following values of shared parameters of stop functions W r and W ch have been obtained: R

Fitting of high-dose implantation profiles
The validation of the 9-parameter high-dose implantation model of section 5 has been performed for two families of experimental doping profiles [32,40] that cover substantially different energies, ions and targets (see figures 5 and 6).It has been assumed that for each profile of this family all model parameters except for the experimentally declared value of ion fluence φ are equal.at such a strong restriction of the model flexibility the results of global fit are in quite reasonable agreement with experimental data, especially for low energy ions (see figure 5).The agreement is somewhat worse at high (∼MeV) energy (see figure 6).But this figure clearly demonstrates the capability of the model of section 5 of describing complicated profiles that consist of three peaks formed by (i) initially random beam fraction, (ii) ions dechanneled by defects at inter- mediate depths, and (iii) well channeled ions.It is absolutely out of the capabilities of the standard "dual Pearson" model.
The deviations of fitting curves from experimental data that one can observe in figure 6 we attribute to a rather approximate assumption of the section 5 model concerning the simple proportionality of a damage profile at a given dose to an implantation profiles in an undamaged crystal.Obviously, this assumption can be rated as too rough for high-energy implantation, and is subject to further refinement within the framework of the phenomenological approach developed.

Conclusions
The proposed physically based phenomenological parameterizations of as-implanted doping profiles in crystalline targets appropriately accounts for the kinetics of ion channeling that leads to bi-directional transitions between channeled and random-like modes of ion motion in crystals.These transitions considerably affect the target depth dependencies of implant concentrations.
The results of phenomenological models validation against the experimental data and their comparative analysis against the fitting capabilities of the standard "dual Pearson" model show that the application of a more physically grounded model approach permits to improve the quality of doping profiles description as compared to that of conventional models.
Except the enhanced profile fitting capabilities, the benefit of the proposed models is the straightforward physical interpretation of their adjustable phenomenological parameters (such as a dechanneling length etc.).Due to the well defined physical sense of these parameters, the application of these models to the description of ion implantation experiments can provide valuable data on fundamental characteristics of the ion channeling kinetics in crystals.Thus, the developed model approach can also be used for the investigation of directional effects at low energy heavy ion channeling where conventional experimental techniques (Rutherford backscattering, nuclear reactions, transmission energy and angular spectroscopy [8]) are generally inapplicable.

Figure 1 .
Figure 1.The MD simulated doping profile of 15 keV B ions implantation into the [100] axial channel of Si single crystal (closed markers) and partial profiles (open markers) of the fractions of ions moved exclusively in a random-like mode (only "random"), exclusively in a channeled mode (only "channeled"), and experienced transitions from one mode to another ("mixed").

Figure 2 .
Figure 2. (a) -the results of the NLSq fit of B→Si[001] doping profiles (markers) using the "dual Pearson" model (dashed curves) and model (10) (solid curves).(b) -the projective ranges R (r,ch) p (E) predicted by our (closed markers) and "dual Pearson" (open markers) models in comparison with estimates from the SRIM2006 databases (dashed and dotted curves).

Figure 3 .
Figure 3.The results (a, curves) of the global NLSq fitting of 100 keV P→Si[001] doping profiles[39] (a, markers) at φ = 10 13 cm −2 and at different thicknesses t of the surface oxide layer.The t-dependencies of the surface capture probability q (b, open markers) and the dechanneling length R ch (b, closed markers) derived from the global NLSq fitting.
.67 nm, σ r = 25.94 nm, γ r = 0.376, R (ch) p = 169.35nm, σ ch = 62.58 nm, γ ch = −3.335.The fit has made it possible to uncover how the parameters q and R ch are directionally dependent on the target tilt angle ψ (see figure 4(b)).The monotonous increase of their values correlates to the beam axis direction approaching the direction of the [112] axial channel.Moreover, the noticeable inflection at a tilt angle θ ≈ 31 • = 35.26• − 4.26 • agrees well with the value of a critical angle ψ c = 4.4 • of stable channeling in this channel of Si crystal lattice (see figure 4(b)).