Документ взят из кэша поисковой машины. Адрес оригинального документа : http://danp.sinp.msu.ru/LIIWM/ArXiv-0904.2151v1_Khodyrev2009.pdf Дата изменения: Sun May 23 02:38:24 2010 Дата индексирования: Mon Oct 1 22:32:52 2012 Кодировка: IBM-866

The shower approach in the simulation of ion scattering from solids
V.A. Khodyrev, R. Andrzejewski, A. Rivera, D.O. Boerma, and J.E. Prieto


Centro de Microan‡lisis de Materiales and Instituto Universitario "Nicol‡s Cabrera", a a Universidad Aut‡noma de Madrid, E-28049 Madrid, Spain o For basic studies of ion-solid interactions as well as for the retrieving of reliable structural information from ion scattering exp eriments, accurate simulations of swift ion-solid interactions are essential. This, however, faces a fundamental problem: due to the small cross sections for scattering, the probability of an incoming particle to reach the detector is always small. As a result, programs based on the "brute force" approach where ion tra jectories are calculated taking into account the large diversity of histories of collisions, require extremely long computation times to reach acceptable statistics. In this pap er, a new simulation approach is prop osed which combines an exact treatment of the interactions within the binary collision model with a high computational efficiency. The main idea is to use the strategy of imp ortance sampling in the simulation of atomic thermal displacements: the p ositions resulting in scattering events close to the direction of the detector are sampled with increased probability. As a result, the detector is illuminated by intensive showers of scattered ions, even if multiple scattering is imp ortant in the way out of the sample volume. The develop ed computer code provides the p ossibility to address wide classes of simulation problems. It is able to treat complex crystal and amorphous structures, to use different models of ion-atom interaction (p otential, inelastic energy loss, multiple scattering on electrons). The program provides information such as depth profiles of close collision yield, energy sp ectra and 2D angular distributions of scattered and recoiled ions. To illustrate these features we present the results of simulations for several, widely different conditions and compare them with exp erimental results. It is argued that, in fact for the first time, this method provides a p ossibility to reliably treat the rare events of plural scattering, the most intricate problem for the theoretical description as well as for computer simulations of ion-solid interactions.
PACS numbers: 02.70.-c, 61.85.+p, 34.35.+a

arXiv:0904.2151v1 [cond-mat.mtrl-sci] 14 Apr 2009

I.

INTRODUCTION

The classical description in the binary collision approximation is well applicable for the description of ion-solid interactions in the energy range above 1 keV; the multiple interaction with two or more atoms, if significant, may be treated in a simple first-order approximation. In the computer simulation of phenomena related to ionsolid interactions, this provides the possibility of an efficient reproduction of many experimental conditions playing with the parameters of the interaction model. Due to the importance of this sub ject great efforts have been performed towards the development of efficient simulation programs [1]. To illustrate the central role that computer simulations have played in the field, it will suffice to mention that one of the most prominent phenomena, the channeling of ions in crystals, was first observed in simulation results [2]. In this context, one should refer also to the paper of Barrett [3] which provides an important supplement to the theory of ion-crystal interaction allowing to address some aspects of the phenomenon which are difficult to treat theoretically. Currently, with ion beams being widely used as a precision tool for the



On leave from Institute of Nuclear Physics, Moscow State University, Moscow 119899, Russia. E-mail: khodyrev@anna19.sinp.msu.ru Corresp onding author. Email: joseemilio.prieto@uam.es

analysis and modification of materials, there has been an increase of the level of demands to simulation algorithms. One example is the use of low-energy ion beams for surface structure determinations. The classical picture of scattering together with the effects of blocking of scattering on atomic pairs form the basis for obtaining detailed information about the atom locations. The only way to extract this information is the comparison of intensities measured for different scattering geometries with the results of simulations for many trial structures. However, a fundamental problem occurs here due to the insufficient efficiency of existing algorithms. Due to the small scattering cross sections, the direct simulation by calculation of individual ion tra jectories (the program MARLOWE [4] is the most developed code of this type) is often not practicable. To understand this one has to keep in mind that the experimental procedure typically consists in the measurement of angular scans with a small-aperture energy detector. To acquire the necessary statistics in a measured spectrum, the required number of ions in the beam amounts to 109 for lowenergy (keV range) scattering of heavy ions like Ne or Ar and to 1013 for scattering of ions with energies in the MeV range. Specially in the latter case it is out of the question to calculate the histories of motion of such a number of particles directly. For this reason, the conclusion is reached [1] that even the power of supercomputers is by far not enough to perform such simulations. An overview of the existing approaches to this problem shows that two main ideas are used depending on


2 the ion energy. In the simulation of backscattering of high-energy channeled ions, one can rely on the singlescattering model assuming that the motion of the ion in the outgoing path can be described by a straight-line tra jectory. This allows to avoid a precise description of this segment of the ion tra jectory. As a result, each ion contributes to the statistics according to the probability of close collisions along the ingoing path. The inverse (blocking) condition can be treated analogously (see Ref. [5], for example). This model is well tractable and serves as the basis for a large number of algorithms which have been developed (see Refs. [6, 7, 8] for the most widely used). This algorithm does not work, however, if one is interested to simulate the channeling-blocking conditions or the scattering of low-energy and/or heavy ions. In these cases, significant influence of crystal structure and/or strong multiple and plural scattering result in the necessity to reproduce the whole history of motion of individual ions inside the sample. As a result, the singlescattering algorithm needs a serious modification: now, to avoid the calculation of outgoing tra jectories, we need to know the probability for an ion to end up in the detector for every position of the atom of the considered lattice site. A solution with emphasis for applications in the medium-energy ion scattering (MEIS) case ( 100 keV light ions) was proposed in Ref. [9] (the computer program VEGAS). According to the timereversibility rule [10], the mentioned probability can be represented by the flux of ions at the same point in the case that the beam enters, inversely, from the detector side. As a result, the simulation consists in the calculation of two fluxes, the direct flux and the inverse flux, followed by their convolution with the probability that the scattering atom occupies the required position. In this approach the electronic and nuclear energy losses are neglected, thus its applicability is guaranteed only at small depths of ion penetration and for relatively light ions. Due to purely technical difficulties (the main one is the need to store large sets of data for the fluxes), this method for medium energies is always used in a simplified form. Namely, the energy and angular variables are not considered during the convolution of fluxes. It is assumed, in fact, that there is a one-to-one correspondence between energy and depth and the scattering angle of the collision is taken to be the same as the angle between the direction of the incoming beam and the direction to the detector. Thus, such approach considers mainly the effect of the spatial distribution of fluxes and cannot claim to achieve a detailed description of energy spectra and multiple (plural) scattering effects. In the case of low-energy (E 1 keV) ion scattering (LEIS) where only a small number of atomic layers contribute to the detected yield, the full version of this "reversing approach" seems to be possible [11]. The only additional problem arising is that with the higher dimension of phase space the accumulated data arrays become very dilute. As a result, the convolution of the fluxes requires some special procedure. For this purpose, it was proposed to sum pairs of crossing tra jectories of incoming and (time-reversed) outgoing fluxes with weights determined by the scattering cross section. It is clear that some tolerances in coordinates and energy have to be assumed in order to identify crossings among the finite set of calculated tra jectories; the associated inaccuracies are equivalent to those present in the case of the numerical convolution of both fluxes. One might expect, however, that allowing for an admissible level of accuracy, the efficiency of the simulation can be considerably increased as compared to the direct approach. Indeed, a satisfactory description of experimental results for several LEIS conditions was found as a result of the application of this approach [12]. In the present paper a new simulation approach is proposed. The developed computer code has the name TRIC (Transport of Ions in Crystals). In this method, complete histories of ion collisions inside the sample are calculated as in the direct simulation approaches. To increase the efficiency we use a specific advance of the Monte Carlo method, the strategy of "importance sampling" [13]: the distribution of thermal displacements of crystal atoms is sampled giving preference to those positions which result in scattering in a direction close to that of the detector. The actual probabilities of such displacements are accounted for in the accumulated statistics by the use of appropriate weighing factors. This procedure results in an illumination of the detector by intensive "showers" of scattering events. Additionally, in order to achieve a smoother statistics of rare events due to plural scattering, the strategy of "stratified sampling" [14] is applied: the main fraction of close collisions which results in scattering to different directions are also sampled more often and these tra jectories are also allowed to send showers towards the detector. The same strategy is used for a preferable choice of the initial coordinates of the ions entering into the crystal among those resulting in higher efficiency in the accumulated statistics. The paper is organized as follows. We present in Sect. II a detailed description of the simulation procedure. In Sect. III the efficiency and quality of the procedure are illustrated by the description of several experimental results under different conditions, selected also based on general interest. Since the reversing approach, closely related with our proposed method, is in fact the only alternative claiming a similar efficiency we pay special attention in Sect. IV to a careful analysis of its basis; the main intention there is to recognize clearly the relative advantages of the two approaches. Some concluding remarks are made in Sect. V.

II.

THE SIMULATION PROCEDURE

The main idea that allows the increase of efficiency of the simulation is based on the fact that, due to ther-


3 In the accumulated statistics, events due to ions in the showers ("secondary" ions) must contribute with weights wi given by the probability for the atom to be displaced into the "hot" region. To calculate them we notice first that a probability should be adscribed also to the primary ion itself (Wi before the i-th collision) which decreases with depth as a result of emission of the secondary ions, (tra jectory "degradation"). Then the weights wi and Wi+1 after the collision are updated as wi = Pi Wi /n,
1

Wi+1 = (1 - Pi )Wi ,

FIG. 1: (color online) The principle of selection of the "hot" region of atomic displacements used in the shower generation. The density of the Gaussian distribution of atomic displacements is shown by the contrasted gray while the "hot" region is indicated as the tub e enclosing the region of impact parameters for scatterings into the chosen angular cone (left). Analogously, the picture at the right panel illustrates the case of showers of recoiled atoms

mal motion, each crystal atom can be found in different positions when the incoming ion crosses the considered lattice site. In the case that the atom is located close to the ion tra jectory, a significant change in the ion motion can take place. Events of scattering into directions close to the direction to the detector are obviously of special interest and it appears desirable to treat them selectively. This procedure should increase the number of relevant events at least in cases where scattering does not strongly perturb the further motion of the ion. To implement this simple idea we use the following procedure. First, the impact parameter bi which corresponds to scattering from the current direction of motion into the direction to the center of the detector (see the left-hand side of Fig. 1) is determined (i denotes the number of the current lattice site). Knowing this value, we can define the line (in the scattering plane) of atom positions necessary for such scattering. Then, assuming that the relevant scattering directions lie within a cone of width , we can define the corresponding region of impact parameters whose half-widths in the scattering plane, b , and in the perpendicular direction, b , can be estimated as b |d/db|-1 , b (b/ sin ) . (1)

In this way also a fraction of the distribution profile of atom displacements is separated, the "hot" region. Then, one or several (n) atom positions can be drawn from this region treating it as a part of the total (Gaussian) distribution of atomic displacements; as a result a shower of ion tra jectories is directed towards the detector. The remaining part of the distribution is used to continue the tra jectory of the "primary" ion. In such a way all possible displacements of the atom are sampled. Notice that the width of the shower is assumed to be significantly larger then the width of the detector .

where Pi is the integral fraction of the probability of the "hot" region. The sum of probabilities is conserved, Wi+1 + nwi = Wi , as necessary. It is worthwhile to note also that this approach provides absolute values of the scattering yield: the expectation value for a certain energy-angular range is equal to that obtained in a direct simulation with the same number of ions sent to the sample. One can easily recognize the similarity of this approach with the strategy of "importance sampling" used in the Monte Carlo numerical integration [13] where the values of the integration variable are sampled more densely in the regions that give the highest contribution to the integral. This similarity is not accidental because the scattering yield, as an average over the ion initial conditions and the thermal displacement of atoms, is, in fact, determined by an integral over a multidimensional space. Due to the inherent complexity of this integral, the importance-sampling strategy can be used only on an intuitive basis as described above. It is clear that this strategy should be effective at high energies where the single-scattering model describes well the process. At low energies the shower approach can also be efficient. In this case the showers emitted at large depths are significantly attenuated due to strong re-scattering on the way out of the crystal. This is a disadvantage but it is largely compensated by the fact that the primary ions themselves can leave with large probability the sample volume and the numerous accompanying showers from the top layers produce an intensive flux in the direction of the detector. It can be said that, in this case, the main part of the procedure is the direct simulation of ion tra jectories and the showers are used only as an improvement at the last stage. These arguments show that the most serious difficulties in the application of the shower approach can be met in the simulation for intermediate ion energies where the rare plural-scattering events result in violent fluctuations (see examples in the next section). As shown at the right-hand side of Fig. 1, a similar procedure can be used for the simulation of the yield of recoiled atoms. When an ion crosses a lattice site occupied by an atom of the considered species, a shower of recoils is emitted. The "hot" region selects now the atom displacements which result in emission of recoiled atoms within the chosen cone. One additional recoil atom is sent by sampling the remaining part of the Gaussian distribution of atomic positions. This recoil, in turn, is allowed


4 to produce new recoils, as in the case of the primary ion, and also scattering showers. The tra jectory of the impinging primary ion is followed further with the atom displacement drawn according to the total Gaussian distribution. To describe the whole cascade we consider also recoils produced by the ions in showers (not allowed to produce new showers). The result of these many possibilities is a strongly developed tree of the cascade. Since all particles in the cascades have similar histories of collisions, the treatment using a recursive algorithm turns out to be very efficient. Violent fluctuations in the accumulation of scattering events can have different causes. A first important point in this respect is that, when the showers are generated at every lattice site crossed by the primary ion, this ion itself cannot scatter into the detector direction. This excludes the possibility of detection of events with weights Wi which are ordinarily much larger than the weights of the secondary ions wi . Thus, the displayed fluctuations are entirely due to the dispersion of the values of wi . If the number of ions in one shower n is fixed, the value of the weight wi for a certain scattering angle depends on the position of the "hot" region with respect to the center of the Gaussian distribution of thermal displacements. Since the width of this distribution is rather small compared to interatomic distances, variations of several orders of magnitude in the values of wi are possible. To improve upon this situation, it is reasonable, instead of fixing n, to fix the value of the weight of one event w0 , treating accordingly the number of ions in a shower as variable, ni = wi /w0 (the fractional part is treated as the probability for sending one additional ion). When wi w0 the resulting effect is expressed as an increase of the number of detected events with a smaller weight w0 from intense showers (distributed to some extent over angles and energies) instead of one event with a large weight. An advantage in the inverse case, wi w0 , is that the load on the computer due to the simulation of non-significant events is avoided. In addition, the discrete counts of such "quanta of probability" closely mimics the aspect of experimental data. For orientation, a reasonable value of the quanta w0 can be estimated defining n in the input as the number of ions in the most intense shower (wi = max). In the case of crystalline structure of the sample one has a further possibility to increase the number of detected events. Clearly, tra jectories of ions entering into the crystal at different points of the surface can produce significantly different contributions to the useful statistics. Thus, in the sampling of the initial coordinates it seems reasonable to choose points with higher efficiency more often. In the program TRIC this is organized as a self-learning procedure: the distribution of initial coordinates, taken initially as uniform over the unit cell of the crystal surface, is deformed in the course of the simulation according to the accumulated information on the efficiency. As necessary, the initial value of the weight W0 is taken to be inversely proportional to the density of the simulated distribution. Effectively, such modification corresponds to a larger number of impinging ions and this results in faster (up to an order of magnitude) accumulation of statistics. Another source of fluctuations is the plural scattering. In order to develop a picture of how these fluctuations arise we should mention the following. An appropriate choice of the width of the shower is determined by the demand that the shower should be wide enough to encompass the whole profile of multiple scattering in the outgoing path. One may believe in this case that the transport inwards and the transport outwards the central region (of width ) are properly balanced in the shower. However, at low energies and/or for heavy ions, this condition is difficult to fulfill because the shower should be taken very wide and, as a result, only a small fraction of the ions in the shower reaches our small-aperture detector. But, if the above condition is not fulfilled, the reduction of ion flux in the shower cone due to diffusion outward is to be compensated by the plurally scattered primary ions. As this takes place, the main contribution comes from primary ions scattered accidentally to directions close to the cone mantle. Since the probability of further scattering into the cone is large for these ions (small |d/db| in (1)) they produce rare but very intense showers of secondary ions with almost equal energies. As a result the accumulated energy spectrum is disturbed by splashes in some channels. We solve this problem in a similar way as described above: the plural scattering events are sampled more frequently than they happen in reality. For this goal a second cone is considered, coaxial with the first and significantly wider than it. Wide showers are produced by the primary ion inside this cone in addition to the showers in the internal cone. Then, ions of the outer cone are allowed to send showers into the internal cone, similarly to the primary ion. As a result, our goal is achieved because the weights of ions in the outer cone are small compared with the weight of the primary ion Wi and the number of showers sent to the internal cone is now large to smooth sufficiently the accumulated energy spectrum. In principle, a whole hierarchy of such nested cones or even some smooth deformation of the density of sampling of atom displacements could be organized. It turns out, however, that in many cases the implementation of the above two-step approach solves practically all problems. This is illustrated in the next section (Fig. 7a). To estimate the efficiency of the shower approach compared to the direct simulation we should consider first the values of the weights of the accumulated events w0 . Since the weight of an event in the direct simulation is unity, the number of impinging ions in the latter should be larger approximately by a factor 1/w0 to achieve the same level of fluctuations in the statistics. Of course, this cannot serve as an accurate estimate of relative efficiency, the calculation of the history of motion with shower emission is more extensive while, on the other hand, the calculated exiting flux is concentrated near


5 the relevant region. Nevertheless, the direct comparison of two approaches, where possible, shows that the mentioned factor reflects well the situation. The value of this factor depends on ion energy and is typically 10-3 at keV energies while at hundreds of keV it is 10-7 . Therefore, the advantage of the shower approach in this respect is enormous. A typical time for simulation of one energy spectrum using an ordinary PC is several minutes, thus the proposed approach is capable to treat problems practically not solvable by the direct approach. Note that both approaches are based on the same, rather general, model of ion-solid interaction. It is worthwhile to note also that the sending of showers aiming to the detector mainly results in misses. This unavoidable drawback, when one is interesting only in the yield of a small-aperture detector, turns out to be an advantage if one needs to simulate the 2D angular distribution of the scattered ions. Such possibility is illustrated in the following section. In fact, we are capable to calculate with the present approach 3D-distributions including the energy scale. The implementation of this algorithm as a FORTRAN program includes all important elements of the binary collision model [1]. A wide possibility to choose the type of the interaction potential is provided. The scattering integrals for binary collisions as functions of impact parameter and energy are tabulated at a preliminary stage of the calculation and used afterwards by interpolating with splines. Additionally, the impact-parameter dependence of inelastic energy losses is considered. To account for the simultaneous interaction with two or more atoms we sum the deviations in ion motion due to the interaction with individual atoms. In principle, all this provides the possibility to perform simulations for energies in a wide range. At high energies the procedure of generation of showers could encounter difficulties due to restrictions in the accuracy of computer calculations. The reason is that the width of the "hot" region in the transverse plane becomes so small that an accurate transformation of impact parameter to scattering angle may be practically impossible. For this reason, if this width becomes smaller than 0.1 of the thermal vibration amplitude u1 , the procedure of shower generation is inverted: the scattering angle is sampled within the shower cone and then, in reverse order, the impact parameter is calculated. As a result, our approach can claim to represent an improvement of the single-scattering model as it additionally accounts for the multiple scattering on the outgoing path. Materials of any crystal structure, including compounds, with rectangular unit cells can be described in the input. On a lower level, however, the program works with a description in terms of the Wigner-Setz cell and, in principle, this provides the possibility to treat also wider classes of structures. An amorphous structure can be simulated by a random rotation of the crystal lattice after each collision. The surface of the sample is defined as an imaginary plane appropriately located with respect to the crystal lattice. If the detector position is defined to be at the back side of the sample, a simulation of transmission through a crystal slab of a given thickness is performed. As output data, the depth profiles of close collisions, the energy spectra and the 2D angular distributions of scattered ions or recoiled atoms of a certain species are calculated. Finally, the procedure is automated to calculate angular scans with step-wise rotation of the target, the detector or the beam direction around a certain axis. The FORTRAN code of the program TRIC is supplied with a graphical user interface allowing to comfortably supply the input data, to run the calculations and to inspect the simulation results.

III.

PROGRAM TESTS AND APPLICATIONS

In this section we present examples of applications of our program based on the shower approach to simulations under different experimental conditions. The first one concerns the scattering of 5 keV Ar+ ions from the (100) surface of a Fe4 N crystal as studied experimentally in Ref. [15]. The structure of this crystal can be considered as fcc Fe with an additional N atom located at the center of the unit cell. The topmost atomic layer at the surface contains nitrogen and, depending on the growth conditions, can exist in two types of surface reconstructions, c4 and 4pg. Figure 2 shows simulated angular scans for the c4 surface differing from the bulk-terminated one only by an outward displacement (0.23 ђ) of the nitrogen atoms. A In these scans, the crystal is rotated around the surface normal, which is coplanar with the beam and the detector directed respectively at 42 and 12 from the surface at opposite sides of the normal. The yield of scattered Ar as well as of Fe and N recoils were calculated. In Fig. 3 the energy spectrum of scattered Ar is shown for the orientation indicated in Fig. 2 by an arrow. Also shown is the angular distribution of scattered ions. The latter illustrates the blocking of scattering from the toplayer Fe atoms by re-scattering from nearby atoms; this largely explains the strong variation of the yield in the angular scan. These variations are strongest when either the beam or the detector are positioned close to the scattering "horizon" which is formed by the reflection from the surface as a whole (very small intensity at the bottom part of the angular distribution). In the case shown in Fig. 3 the detector is at the border of the shadow cone. The shape of the energy spectrum as an isolated peak is due to the blocking of the particles scattered from atoms of deep atomic layers by atoms of the top layer. This shape of the spectra favors the choice of the energy window for the performance of angular scans. The double-humped structure is explained by the effective competition of single- and double-scattering from Fe atoms of the top layer [16]. As seen in Fig. 2, the experimental results are well reproduced in the simulation both for scattered ions and for recoils of the two atomic species. Note that, in order to achieve this agreement, the screening radii in the


6

Fe4N(100)
unreconstructed surface Ar yield (arb. u.)

0

-40

-20

0

20

40

60

80

100

120

Fe yield (arb. u.) 0 N yield (arb. u.) -40

-20

0

20

40

60 exp. TRIC

80

100

120

0

-40

-20

0 20 40 80 60 Azimuthal angle (deg.)

100

120

FIG. 2: Comparison of simulated angular scans for the scattering of 5 keV Ar+ ions from the surface of a Fe4 N(100) crystal with the exp erimental results from Ref.[15] (top panel). The yield of recoiled Fe and N atoms is shown in the center and lower panels, resp ectively. See text for more details.

used "universal" potentials [17] were properly reduced, as commonly done for the description of scattering of heavy ions at low energies. With the same potential the data for the more complex reconstruction 4pg are also described with the same quality. In general, the present results demonstrate that the quality of the description of the experimental data is similar to that provided by the code MATCH [15]; in fact, the results of both simulation can hardly be distinguished when plotted together. The second example, shown in Fig. 4a, corresponds also to a LEIS experiment, now of 3 keV Ne+ scattering from a polycrystalline Cu sample. Under these conditions an intriguing peak is observed [18] in the experimental spectrum at the energy corresponding to a single Ne-Cu collision. Here the beam incidence was along the surface normal and the scattering angle was 129. For these strong-scattering conditions the code MARLOWE is capable to reproduce the general shape of the spectrum including the surface peak. Our simulations shown in the same picture also reproduce the spectrum shape; some differences with the MARLOWE simulation at low energies may be due to a difference in the used potentials of ion-atom interaction or in the treatment of inelastic energy losses. In general, the shape of the spec-

FIG. 3: Energy sp ectrum of Ar+ ions scattered from the (100) surface of a Fe4 N crystal for the orientation indicated in Fig. 2 by an arrow. The lower part of the figure shows the angular distribution of scattered ions. Results are shown using a gray scale with darker representing higher intensity.

trum under these conditions differs considerably from the spectrum predicted by the single-scattering model (also shown in the figure). The authors of Ref. [18] associated the appearance of this peak with the onset of plural and multiple scattering in the deeper layers. First, due to the the large variations of kinematic energy losses, the plural collisions in scattering from deeper layers yield a broad scattering spectrum (in particular, this explains the presence of the high-energy tail) while (single) scattering events from the top layers form the narrow peak. On the other hand, the latter is also favored by the relatively small inelastic energy losses: independently of the depth of scattering, all single-scattering events have similar energies at the exit. To demonstrate the decisive role of strong localization of energy losses at small impact parameters, we present in Fig. 4b the results of a simulation where both nuclear and electronic energy losses are replaced by equivalent continuous stopping forces (thin solid curve). As seen, with such description of energy losses, the spectrum changes drastically. First, the sur-


7
10 8 6 4 2 0 520 10 8 6 4 2 0

3 keV Ne -> Cu
Yield, arb.units

+

Exp. MARLOWE Single Scatt. TRIC

a

1020

1520
Full simulation Uniform stopping Single scatt.

b

Yield, arb.units

520

1020

1520

Energy, eV
FIG. 4: (a) Comparison of our simulation results for 3 keV Ne+ ions scattering from amorphous copp er with exp erimental data [18] and results of simulation with the MARLOWE code. The simulation results were folded with a gaussian distribution in order to simulate the exp erimental energy resolution of 40 eV. The dashed curve corresp onds to a calculation in the single-scattering approximation. Panel (b) shows simulations p erformed by applying additional approximations. The dashed curve corresp onds to a simulation using the single scattering approximation with impact-parameter-dep endent energy losses while the thin solid curve was calculated considering the energy losses as the result of a uniform stopping force. FIG. 5: Results of the simulation of 100 keV He+ ions scattering from a Si crystal. The scattering yield is accumulated from depth ranges of (a) 5-30 nm and (b) 30-55 nm. The displayed angular range lies around the <112> crystal direction.

face peak disappears since single scattering events from different depths have now different energies. And second, the variation of kinematic energy loss in plural scattering is now not effective and, as a result, the shape of the spectrum on both sides of the surface peak is also strongly modified. On the other hand, in the shower approach we can check also the role of plural and multiple scattering. For this purpose, the simulation was performed using a special procedure, in which in all collisions before and after the emission of the shower, the ion deflection was canceled (the single-collision model with impact-parameterdependent energy losses). These results are also shown in Fig. 4b as a dashed curve and demonstrate that the multiple and plural scattering events themselves are of minor importance for the formation of the surface peak. In summary we can conclude that the surface peak reflects mainly the correlation between scattering and energy loss: at this low energies ions scattered in deeper layers have more chances to exit from the sample if they do not re-scatter strongly on atoms of the upper lay-

ers and, consequently, the energy losses in the passage through these layers are also small. It is worthwhile to note that the program TRIC passes here a serious test: the simulation with shower generation produces results identical to those obtained when running the program in the direct simulation mode. Turning now to medium energy ion scattering (MEIS), we show in Fig. 5 the angular distribution of 100 keV He+ ions scattered on a Si(100) crystal. The geometry (the <112> axis is in the center of the position-sensitive detector) and the two depth ranges chosen for collection of the scattering events are equivalent to the experimental conditions used by Kobayashi [19]. The results of the simulations are also very similar to the experimental results: even at such relatively small depths the washingout of the blocking pattern is well seen, first of all for the


8

Si Y - - [100]/[111] a

Yield, arb.units

1.4 1.2 1 0.8 0.6 0.4 0.2 1.2 40 45

50

55

60

65

70

Yield, arb.units

-- [011]/[100]

b

1

0.8

0.6

70

75

80

85

90

95

100

Scattering angle, dg.
FIG. 6: Comparison of our simulations (solid curves) of angular scans for 100 keV p scattering from a monolayer-thick YSi2 film on Si(111) with the result of a VEGAS simulation (dashed curves) and exp erimental data [20] (dots). The direction of the incident b eam is parallel to (a) the <п00> and (b) 1 <0пп> directions of the Si substrate while the detector lies in 11 1 the plane (01п). The structure of surface layers is shown at the top. The yield is normalized to the Rutherford cross section. The dash-dotted curve in the lower graph shows results for a trial structure with Y atoms located at the top layer.

narrow channels. Although this effect of re-channeling is well predictable, its detailed demonstration in the referred experiment is rather interesting and our simulations support these results. Notice that, with 105 ions being sent to the crystal in the simulation, the total yield over the detector amounts here to 0.01 only (remind that, in the direct simulation, this would be an expected number of counts). In this sense, the results of our simulation are unique, it is hardly possible to obtain them using other known methods. In Fig. 6 we show MEIS angular scans for the yield of 100 keV protons scattered from Y atoms of a YSi2 monolayer epitaxially grown on Si(111) (a side-view of the structure is shown in the picture). The results of the

measurements and of the simulations with the program VEGAS are taken from Ref. [20]. The dips in the scans are due to blocking of the scattering from Y atoms by Si atoms of the upper layers. All parameters in the two simulations are taken to be identical. The achievement of a perfect agreement by optimization of parameters of crystal structure and lattice dynamics is, in fact, the actual goal of the MEIS analysis [20]. Here, we emphasize only that results of our simulation coincide well with the results of the VEGAS simulation. In principle, this is predictable for such thin layers (see Introduction). As an additional test of our simulation procedure we performed simulations also for the case when the top layer is terminated by Y atoms instead of Si atoms. As is seen from the figure, the scattering yield in this case simply reproduces the Rutherford cross section. The most difficult problem for simulations is the description of conditions where plural scattering is strongly effective. Besides low energies, this happens also in the medium and high-energy cases where heavy-ion beams and/or heavy-atom targets are considered. For the case of scattering from amorphous samples, Biersack et al. [21] proposed the fast version of the code TRIM [17] where the consumption of computer power is reduced at the cost of additional approximations. Figure 7a shows their results [22] for scattering of 100 keV protons from a 1000 ђA thick gold foil, together with the results of simulations done with the shower approach. As seen, the two approaches produce in fact identical results, also for similar computation efforts. The only visible difference is at the low-energy tail where yield is entirely due to the plural scattering. Therefore, these results can be regarded as a confirmation of the accuracy of the TRIM approach by comparison with an exact treatment of the model with a similar structure of the target. In Fig. 7b the simulation results are compared with results of an experiment [23] of 280 keV -particles scattering from a 1130 ђ-thick Pt foil. In general, the agreeA ment is satisfactory also here. In the two cases shown in Fig. 7 the shape of the spectra differs significantly from those predicted by the single-collision model: the yield is higher for the main part of the spectra, also the lowenergy tail which is entirely due to the plural scattering cannot be predicted at all by the single-collision model. To reproduce all the features above the fluctuation level, the showers must be generated in wide cones. In the considered cases the double-cone approach (see Sect. II) was used with a half-width of 40 for the internal cone and of 160 for the outer cone. In this way, at least the double-scattering events are treated with special efforts. The effect of the usage of the two cones is illustrated in Fig. 7a where, for comparison, the simulation results obtained in the single-cone approach are included. The level of fluctuations in this case suggests that much larger efforts are necessary to achieve the same level of precision of simulation as in the double-cone approach. To demonstrate the role of plural scattering in more detail, we show (Fig. 7b) separately the contribution of showers


9
10

4

Backscattering yield, arb.units

0.1 MeV p -> Au
8 6

Yield, arb.units

RBSTRIM Singl.sc. TRIC TRIC, 1 cone

a
3 2

1 MeV p -> <100>W

4

2

1
0.05 0.06 0.07 0.08 0.09 0.1

Channeling Blocking Blocking/cos

0 0.04

in

Energy, MeV
5

0

0

1

2

3

Backscatterinng yield, arb.units

4

0.28 MeV He -> Pt

3

Exp. Singl.sc. TRIC Partial (MS/PS)

b

Angle of misalignment, dg
FIG. 8: Channeling and blocking wells calculated for the case of irradiation with 1 MeV protons of a tungsten crystal along the <100> axial direction with detection at a 135 scattering angle and in the inverse geometry. The yield is calculated in the energy window corresp onding to scattering from a depth A range b etween 2500 and 3500 ђ. The blocking dip normalized to the thickness of the layer along the b eam direction is also shown.

2

MS PS

1

0 0.03

0.08

0.13

0.18

0.23

0.28

Energy, MeV
FIG. 7: Backscattering energy sp ectra of 0.1 MeV p scattered from a 1000 ђ Au film (a) and of 0.28 MeV He ion scattering A from a 1130 ђ Pt foil (b). The simulation results and the calA culations in the single-scattering approximation are compared with the RBSTRIM simulation [22] and with exp erimental data [23]. Results of simulation in the single-cone approximation are shown in the top panel. The lower graph shows the partial contributions of ions from the internal (MS) and outer (PS) cones.

emitted by ions moving in the outer cone. It is seen that these histories of collisions completely explain the nature of the low-energy tail and, in general, contain mainly the effect of plural scattering. On the other hand, the partial spectrum of ions in the showers produced directly by the primary ion is influenced by multiple scattering. It looks surprising at a first glance that, in contrast with the predicted increase of the yield compared to the singlescattering spectrum [24, 25], this partial spectrum shows the inverse ordering due to multiple scattering. This is, however, natural because in this partial spectrum we do not consider the compensation of the transport out of the internal cone by the transport coming from the external region. As a last example we reproduce in our simulations the experiment [26] performed in the early times of chan-

neling studies with the special aim to test the Lindhard time-reversibility rule [10]. Here a proton beam of high energy (1 MeV) is incident on a W crystal along the <100> direction (close to the surface normal) and the protons scattered by 135 into a random direction of the crystal lattice are detected. In the inverse situation the beam and detection directions were interchanged. The angular spread in the beam and the aperture of the detector were approximately equal, 0.1 , and the yield was measured in an energy window corresponding to a layer of 1000 ђ thickness at a depth of 3000 ђ. The yield was A A measured as a function of the beam misorientation in the first case and as a function of the detection angle relative to the <100> channel in the second. The simulated channeling and blocking wells are shown in Fig. 8. Their widths are similar and coincide well with the experimental results. To achieve also the coincidence of the absolute values we had to normalize the yield to the path length of the incoming ions inside the considered layer (by multiplying the yield by cos in in the second case, where in is the angle of beam incidence relative to the surface normal). The meaning of this procedure is explained in the next section. In general, the timereversibility is confirmed by the simulation similarly as in the referred experiment.

IV.

DISCUSSION

In this section we discuss the advantages of the proposed approach in comparison with those used currently


10 in simulations of swift ion-solid interaction. First, the shower approach is equivalent in accuracy with the direct simulation as performed by the code MARLOWE. The difference lies only in the much higher efficiency (lower level of statistical fluctuations) as illustrated in Sect. II and III. This is in fact of primary importance since, in many situations, extensive analysis of experimental data becomes feasible. Furthermore, for medium and higher ion energies a proper simulation of multiple and plural scattering is not possible at all by other methods. For medium-energy ion scattering from amorphous samples, the accelerated version of the code TRIM [21] seems competitive but this is achieved at the cost of additional approximations. Barrett's approach [3], followed by later developments like the programs FLUX [6] and UPIC [8], is based on the single-scattering model and has consequently its region of applicability restricted to the cases where multiple and plural scattering can be neglected. In principle, the calculation of the close-encounter probabilities in this approach bears some similarity with our procedure of shower generation. In this way, the picture of collisions with small impact parameters is well reproduced. However, only one of the two segments of the tra jectory is described realistically while the other is approximated by a straight line. The reversing procedure used in the programs VEGAS [9] and MATCH [11] is closely related with the shower approach and it is rather instructive to compare both strategies in more detail. First, to fill some gaps in the published descriptions, we present here a general analysis of the basis of this approach in a systematic way. The problem is formulated in terms of the fluxes of incoming and outgoing ions properly convoluted at certain lattice sites. Analogously as for the incoming flux, a certain distribution of initial conditions (see below) has to be assumed for the calculation of the time-reversed outgoing flux. Therefore, the problem is formulated from the outset differently than in the direct simulation: the simulation is aimed to obtain the integral yield in some range of exit directions and of energies of scattered ions E ? E + E (the detector acceptance). In order to treat the problem in the most general way, we will not assume in the following that or E are small. Further, as in any non-direct approach, special care must be taken in order to ensure that the statistics of thermal displacements of target atoms is treated correctly in order to guarantee that, in this respect, the approach is equivalent to the direct simulation. Basically, since the ion velocity is much larger than the thermal velocities of the target atoms, each ion passing through the sample volume sees a frozen pattern of atom displacements. The configurations of displacements must be chosen randomly according to the model assumed for the description of thermal vibrations. In the reversing approach, as also in the shower approach, we use the following circumstance [27]: when averaging the scattering yield over different configurations of atom displacements one may replace the contribution of each configuration by the average over the displacements of one of the atoms, arbitrarily chosen. Note that the probabilities of displacements of this atom can be dependent on the actual (fixed) positions of other atoms; in fact, accounting for correlations in the thermal vibrations is needed in any realistic description. The possibility of averaging can be used to smooth the statistics of scattering events in the simulation: the smoothing is achieved by multiple sampling of the displacements of an atom met along the ion path. For every resulting tra jectory the same procedure can be iterated at the next lattice site (the atom positions chosen in the previous collisions must be considered as a modification of the initial configuration of the surrounding). Note that the direct simulation is usually performed in the same manner, the only special feature here is that only one sample is taken at each lattice site. Multiple sampling of atom positions allows to obtain several scattering events for every impinging ion. The advantage of this is that a significant part of the calculations is common to all tra jectories. A dramatic increase of efficiency in simulation using the shower approach is achieved due to a special way of sampling which enhances the probability to finally reach the detector. To apply the strategy of importance sampling we appeal, in fact, to the single-scattering model assuming that it can be used at least as a first approximation. The reversing approach is based on a different idea. The random sampling of atom displacements has, in fact, the only purpose to find those positions of the atom which result in useful scattering events. Every such position defines simultaneously the corresponding outgoing tra jectory which starts from the vicinity of the considered lattice site and ends up in the detector. Therefore, the problem could be solved in the reversed order, starting with the identification of these tra jectories. A set of such tra jectories can be generated considering the time-reversed motion of ions starting from the detector. With this information the random sampling of atom displacements is not needed, the relevant positions of the atom are those which provide a connection of the tra jectory of the ingoing ion with one of the outgoing tra jectories. The tra jectories are to be connected as a result of the scattering on this atom. To recognize clearly the relation between the shower and reversing approaches we can formally write the considered contribution to the scattering yield Yi (i is the site number) as Yi =< d|T2 Si T1 |b >, (2)

where |b > denotes the flux of ions in the beam before entering into the sample volume, Si the operator of scattering on the atom and operators T1 and T2 describe the transformations of the flux before and after the collision with the i-th atom, respectively. In the shower approach the resulting transformation of |b > is pro jected on the detector acceptance profile |d > (|b > and |d > can be considered as elements of a real Hilbert state). The reversing approach consists in this picture in the action of


11 the two operators T1 and Si followed by the pro jection on the state T2-1 |d >. To specify the way this plan can be implemented, we need to represent the action of the operator Si explicitly. Denoting by 1 the phase variables of the incoming ion entering the lattice site volume, we can write the probability that further motion will end up in the detector as Pi (1 ) = dPsc d2 Pout . d2 (3) in such a way that dx2 = sin dzc . In addition, we take into account that, for a given scattering angle, zc is directly related to the coordinate z of the atom: dzc = dz . As a result we arrive at the familiar expression: d d3 Psc = D(R) dz d2 n2 d (6)

The integration is performed over the phase variables of the ion after the collision 2 , (dPsc /d2 )d2 is the probability of scattering into one of the states within d2 ; this probability is related to the probability that the atom is located at the appropriate position. Finally, Pout (2 ) is the probability that, in the course of its further motion, the ion in a state within d2 will leave the crystal with an energy and in a direction within the detector acceptance. The probability dPsc /d2 is explicitly determined in the case that both the scattering angle (b) and the energy loss E (b) in the binary ion-atom collision are uniquely determined by the impact parameter b. To derive the corresponding expression we describe the states 1 and 2 in a local coordinate system where the z axis is aligned with the ion velocity before the collision. As a result, dPsc 3R D(R), = (y2 - y1 ) (E2 - E1 + E (b)) d2 x2 2 n2 (4) where x and y are the coordinates in the scattering plane and in the perpendicular direction, respectively. The first -function satisfies the condition that the two tra jectories must intersect and the second takes into account the relation between the energies before and after the collision. In the case of potential scattering the state of ion motion after the collision (given by the coordinate in the scattering plane x2 and the velocity direction n2 ) is uniquely determined by the atom position R. The Jacobian of this relation appearing in (4) is expressed as 3R b 1 d db = , = 2 2n x2 2 sin d sin d

which expresses nothing but the concept of differential cross section d /d. The contribution of scattering from the considered atom to the yield at the detector is obtained as an integral of the probability (3) multiplied by the flux b (1 ) of ions of the beam reaching this lattice site: Yi = d1 b Pi . (7)

Further, the variables 1 and 2 are uniquely related to in and out respectively, the parameters of motion of the ion when it enters (exits) the sample volume. It is useful to refer directly to the latter parameters and, for this purpose, we replace the integrations over 1 and 2 in (3) and (7) by integrations over in and out , respectively. Then, the Jacobians J1 = |d1 /din | and J2 = |d2 /dout | must be added into the integrands. Note only that the relation 1 (in ) cannot be defined for all in (some initial conditions in result in tra jectories which never pass through the vicinity of the considered lattice site). The same is possible for the pair of variables 2 and out . To account for this, we assume the value of the respective Jacobian in such cases to be zero. Now, combining the above equations, we arrive at the expression: Yi = din Jin
b

d

out Jout Pout

§

§ (y2 - y1 ) (E2 - E1 + E (b))

1 d D(R). (8) sin d

(5)

where b is the impact parameter corresponding to the scattering angle and d /d is the scattering cross section. Finally, D(R) in (4) is the density of distribution of thermal displacements of the target atom. Note that R, the atom position which results in the considered collision, is a function of 1 and 2 . To confirm the validity of Eqs.(4) and (5), let us calculate the angular dependence of the probability of scattering on a single atom (the variation due to the uncertainty of the atom position R). First, we integrate both sides of Eq.(4) over y2 and E2 and this cancels the two functions. The coordinate x2 is related to the coordinate zc of the point of crossing of the scattering asymptotes,

As a function of out the probability Pout is easily evaluated: Pout = 1 if out is within the region of acceptance of the detector, otherwise Pout = 0. One can interpret this as the operation of pro jection of the outgoing flux on the exit states, the set determined by the detector acceptance; the "density of states" is here uniform. In general, all aspects of Eq.(2) are reproduced in the above equation. Particularly, the operators T1 and T2 are represented as the transformations of variables in (1 ) and out (2 ) with the corresponding Jacobians. Eqs. (2) and (8) show that the two segments of the ion tra jectory, before and after the collision with the considered atom, can be treated analogously. The flux T2-1 |d > can be represented as the result of the timereversed evolution of the "beam" with a uniform profile (uniform distribution of angles and energies within the detector acceptance and, also, uniform distribution of initial coordinates at the sample surface). This shows a way in which the integral (8) can be evaluated using the Monte-Carlo approach. The relations of variables


12 1 (in ) and 2 (out ) can be reconstructed by the calculation of a number of tra jectories of incoming ions with the corresponding distribution of initial conditions and a number of tra jectories of outgoing ions calculated in time-reversed mode (including negative energy losses). The densities of the two fluxes approaching the considered lattice site are represented by the factors Jin b and Jout Pout , respectively, in the integrand of (8). The integration itself, however, presents here a serious problem. First, it is difficult in practice to reconstruct in numerical form the distribution of the fluxes in the 5dimensional section of the phase space. And, second, the presence of -functions in (8) results in the need of some special procedure. For this reason, the standard numerical integration is practicable only in a simplified version of this approach, as performed by the computer code VEGAS, where the energy is not resolved and the scattering angle is assumed to be fixed. In this case the procedure is reduced to the convolution of two two-dimensional fluxes at a closed surface containing the considered atom. If the fluxes are represented by sets of tra jectories, the conditions expressed in (8) by the -functions (the "conditions of crossing" for a pair of ingoing and outgoing tra jectories) cannot be satisfied exactly. Therefore, to accumulate a representative statistics, one has to introduce in these conditions some tolerances, sacrificing thereby the accuracy of the simulation. In the numerical integration where the density in phase space is presented by histograms, the sizes of the bins of the histograms play the role of the tolerances. Another procedure for the integration is implemented in the program MATCH. Pairs of incoming and outgoing tra jectories are matched to identify that crossing at the considered lattice site. The tolerances for coincidence of coordinate y and for the energy relation, y and E , are specified beforehand and, in terms of Eq.(8), the whole procedure can be formulated as follows. The tolerances are introduced as a replacement of the -functions by normalized pulse functions (t) (t) = The integration over Monte-Carlo sum: d
out

1

1 if |t/ | < 1/2 0 if |t/ | > 1/2

(9)

out

in (8) is associated with the S0 E Nout
N
out

Jout Po

ut



,
i=1

(10)

where S0 is the area of the unit cell of the crystal surface and Nout is the number of calculated outgoing tra jectories. To obtain the yield of scattering for Nin ions in the beam, we replace the integral over in (together with the Jacobian Jin and the beam profile b ) simply by a sum over all impinging ions. The terms of the resulting double sum are determined by the rest of the integrand in (8): wi = S0 E 1 d D(R). Nout y E sin d (11)

Hence, each crossing of tra jectories contributes to the accumulated statistics with its own weight. The weight (11) differs from the weight actually applied in the program MATCH, where it is simply taken as the product of the cross section and the atomic density. The dependence on the angle in (11) includes, however, an additional factor 1/ sin . The necessity of such factor is immediately recognized when one tries to calculate the yield of scattering from one atom (6) by the simulation. The factor 1/ sin in (11) is compensated in this case by the number of crossings of tra jectories, which is proportional to sin ; the latter follows from a simple geometrical consideration. The first fractional factor in (11) plays the same role, the normalization on the number of crossings. As a result, the Monte Carlo sum, the sum of weights (11), gives the absolute value of the scattering yield Yi . Interestingly, the first fractional factor in the weight (11) is reminiscent of the quantum definition of the number of states, the volume of phase space divided by the Planck's constant. Here, analogously, the tolerances y and E represent the volume per state for states of ions exiting from the sample volume; the phase space attains thereby a coarse-grained structure. As stated above, the described procedure is equivalent to the simulation in the direct mode, the result represents the contribution of the chosen configuration of atoms displacements additionally averaged over displacements of the i-th atom. As stated above, such averaging does not bias the results of the simulation and it might be expected that the resulting smoothing of the accumulated statistics will improve significantly the efficiency of the simulation. However, one can meet serious problems in the practical implementation of this approach. The first problem is the choice of reasonable values for the tolerances y and E . In fact, the use of tolerances introduces some inaccuracy into the description of scattering and, consequently, the tolerances must be small to provide sufficient accuracy of the simulation results. On the other hand, small tolerances imply smaller probabilities of tra jectory crossings and, therefore, the efficiency of the calculation is lowered. Thus, the optimal values are the maximal values of y and E which still provide sufficient accuracy. The inaccuracies have the following origin. Replacement of the -functions in (8) by finite-width pulse functions is allowed only when the fluxes Jin and Jout are sufficiently smooth. In fact, we assume here that each calculated tra jectory of ion motion can be associated with a family of close tra jectories, the family limited by the tolerances, and they are treated simultaneously. The critical points in this assumption are the following. First, the simultaneous treatment of the family of tra jectories assumes that they all reach the detector. This means that the tolerances should correlate with the detector acceptance, the "resolution" in the connection of tra jectories, so to say, must be finer than that of the detector. In this connection, one should also take into account that, after many collisions with atoms on the outgoing path, the


13 ions of the family can exit from the sample in strongly different directions and, as a result, not all of them will be detected. It is clear that this difficulty increases when the ion energy is increased (scattering from large depths). Applying this approach one must decide also at which lattice site the desired matching procedure will be performed. Assuming, as in Eq.(8), that the site is the same for all incoming ions, one has to choose the site in the top layer, the first one met along the ion path. Otherwise, this atom should be considered as being at a fixed position and, therefore, direct scattering in the direction to the detector would be possible. However, in such implementation of the approach, the problem of strong statistical fluctuations is not entirely eliminated. The fluctuations manifest themselves as extremely large values of the weight (11) when two tra jectories intersect at a small angle . Such fluctuations rarely appear in the cases when a large-angle scattering (with small cross section) happens at some other lattice site. Here, such scattering events (happening at larger depths) are present in the outgoing tra jectories. In fact, this means that the largeangle scattering is treated here in a way equivalent as in the direct simulation and this is a clear disadvantage. The procedure can be modified in such a way that the connection of ingoing and outgoing tra jectories is performed only at the points of large-angle scattering. The modification consists in the following. First, the interesting site will be chosen as specific for each ingoing ion. The whole tra jectory of ion motion in the frozen lattice is inspected to find the site where the ion most closely approaches the equilibrium position of the atom. If the matching procedure is performed at this site, the range of impact parameters of possible collisions is restricted, it is limited by the amplitude of thermal vibrations of the atom. In the case of low-energy heavy ions (strong interaction) this means an effective restriction of the possibility of small-angle scattering. However, even at low energies, this does not ensure an appropriate low level of fluctuations of the weights wi . At higher energies this becomes impossible at all, the allowed values of the impact parameter correspond to large variations of . It is readily recognized that this is the same problem as that solved by the shower approach. Again, the main difficulties lie in the treatment of the plural scattering. To solve this problem we can use again the main idea of the shower approach, which is to emphasize the events of scattering in the direction to the detector. The only difference will be that the multiple sampling of atom displacements from the "hot" region is replaced by a corresponding selection of crossings of tra jectories. Namely, the position of the atom necessary for the connection of ingoing and outgoing tra jectories must be located within the "hot" region. Therefore, the reversing approach turns out to be nothing but a special version of the shower approach. In order to choose, for a given application, which of the two versions is more convenient, one should consider the following. First, the simulation with the reversing version cannot be considered as an exact treatment inside the binary collision model since the finite tolerances for connection of tra jectories represent an unavoidable source of inaccuracies. Practically, the only way to control the inaccuracies is the repetition of the simulation using different values of y and E . Obviously, this is a serious drawback of the approach. Another drawback is that the whole procedure of the reversing approach is too cumbersome for an extensive usage and for developments of the model of interaction. The arguments above show that, with some simplifications, this approach can be applied to the simulation of scattering of low-energy heavy ions (as it was used before) but is hardly practicable at higher energies. It is even difficult to estimate, for example, what kind of efforts would be necessary to reproduce the data presented in Fig. 7. For the program TRIC, this is also a difficult problem, the calculations with an ordinary PC take about 20 hours. One can conclude, however, that in general, this program is capable to address by simulations a wide class of problems of scattering and recoil emission. To finish this section we make some remarks intended to clarify possible consequences of the reversibility rule for the interpretation of experimental results. It is seen in Eq.(8) that, under the assumption of a purely potential motion of the ions (the Jacobians Jin = Jout = 1), the yield of scattering from one lattice site is symmetric under interchange of the beam and detector directions provided that the flux in the beam is also uniform and the respective phase-space volumes are equal. The latter condition is less relevant because the difference can be simply accounted for by a factor in the yield. Therefore, we can conclude that the yield of scattering from one atom for a given beam-detector configuration is determined if it is known for the inverse situation. However, in the measurements of the yield in a certain energy window, as ordinarily made, the effective number of contributing atoms can be different. This fact was taken into account in the transformation of the data shown in Fig. 8.

V.

CONCLUSION

The shower approach proposed in this paper solves effectively the main problem of simulation of ion scattering from solids within the binary collision model, which is elimination of the violent statistical fluctuations in the Monte-Carlo sum. This is achieved as a result of specific effective improvements of the direct simulation approach. The improvements consist in the use of the strategies of importance- and stratified sampling, special ingredients of the Monte-Carlo method. As a result, the computer power required for simulation is reduced by several orders of magnitude. This is, in fact, a decisive advantage allowing to address simulation problems which cannot be treated with other methods. As discussed in Sect. IV, our method avoids also many shortcomings inherent to alter-


14 native approaches. It is argued in particular that, in fact for the first time, this method allows a reliable treatment of the rare events of plural scattering. Such possibility is especially important because the plural scattering is also not amenable to theoretical treatment and the only proposed method of simulation, the modified version of the code TRIM, has significant restrictions. First of all, it is capable to treat only amorphous structures. We performed a detailed analysis of the basis of the reversing approach, as applied in the programs VEGAS and MATCH. The latter program offers the only alternative for an exact treatment of the binary collision model, including multiple and plural scattering. It is argued in particular that the corresponding procedure for matching must use a similar scheme of importance sampling as used in the shower approach. The main difference lies in the method of sampling. Basically, if the scattering from shallow depths is considered, the two methods are equivalent. The applicability of the reversing approach at larger depths needs, however, a special justification. Also, right from the beginning, this approach uses an approximation (the finite tolerances assumed for the connection of the two segments of the ion tra jectory). Note also that the procedure used in its practical application is much more cumbersome than that used in the shower approach. In fact, in simulations for medium and high energy ions this approach becomes impracticable at all. Simulations with the code TRIC can be performed for large classes of crystal structures of the samples and provide a detailed picture of scattering or recoil yields in the form of energy and angular distributions. All these qualities are promising for a wide use of the developed computer of materia efficiently for many structural code both in basic research and in the analysis ls. In particular, this provides the possibility to compare measurements with simulations made trial structures allowing in this way a precise analysis.

This paper shows several examples of the use of the program although, of course, it is not possible to cover all potential applications (e.g. simulations of the sputtering or total reflection yields, well possible in this approach, are also interesting applications of the code TRIC). The illustration examples in Sect. III are chosen to demonstrate the capabilities to solve specific problems and to show the accuracy of this method in comparison with others. In particular, we show the capability to simulate the interaction of different ions of low and medium energies with solid matter including complex structures, to calculate the yield of scattered ions and recoils and to reproduce their energy spectra and angular distributions. In addition, we address the interpretation of the time reversibility rule and provide additional insight into the origin of the surface peak.
Acknowledgments

V.A. Khodyrev gratefully acknowledges the hospitality provided by the CMAM during several research stays. Authors thank A. Guirao and J. Alvarez from CMAM for technical help with computer resources. This work was financed by pro jects MAT2005-03011 and FIS2008-01431; J.E. Prieto and A. Rivera also acknowledge support by the programs "Ram‡n y Ca jal" and "J. de la Cierva" of o the Spanish MEC, respectively.

[1] W. Ekstein, Computer Simulation of Ion-Solid Interactions. Springer-Verlag Berlin Heidelb erg, 1991. [2] M.T. Robinson and O.S. Oen, Phys. Rev. 132, 2385 (1963). [3] J.H. Barrett, Phys. Rev. B 3, 1527 (1971). [4] M.T. Robinson and I.M. Torrens, Phys.Rev. B 9, 5008 (1974). [5] S.A. Karamian, J.S. Forster, J.U. Andersen, W. Assmann, C. Broude, J. Chevallier, J.S. Geiger, F. Gruner, и V.A. Khodyrev, F. Malaguti, and A. Uguzzoni, Eur. Phys. J. A 17, 49 (2003). [6] P.J.M. Smulders and D.O. Boerma, Nucl.Instr.Meth. B 29, 471 (1987). [7] M. Posselt, Crystal-TRIM, May 2004, http://www.fz-rossendorf.de/ . [8] V.A. Khodyrev, V.Ya. Chumanov, K.K. Bourdelle and G.P. Pokhil, Nucl.Instr.Meth. B 94, 523 (1994). [9] R.M. Tromp and J.F. van der Veen, Surf.Sci. 133 159, (1983); J.W.M. Frenken, R.M. Tromp, and J.F. van der Veen, Nucl. Instr. Meth. B 17, 334 (1986). [10] J. Lindhard, Kgl. Danske Videnskab. Selskab, Mat.-Fys. Medd. 34, No.14 (1965). [11] G. Dorenb os, M. Breeman, and D.O. Boerma, Nucl. Instr. Meth. B 108, 173 (1996); M.H. Langelaar, M. Bree-

[12] [13] [14] [15]

[16] [17] [18] [19] [20] [21] [22]

man, A.V. Mijiritskii and D.O. Boerma, Nucl. Instr. Meth. B 132, 578 (1997). D.O. Boerma, Nucl. Instr. Meth. B 183, 73 (2001). J.P. Lepage, Journ. of Computational Physics 27, 192 (1978). H. Press and J.R. Farrar, Computers in Physics 4, 190 (1990). ‡ S.Y. Grachev, J.M. Gallego, D. Ecija, D.O. Boerma, R. Gonzalez-Arrabal and R. Miranda, Nucl. Instr. Meth. B 219, 593 (2004). E.S. Mashkova, V.A.Molchanov, E.S. Parilis, and N.Yu. Turaev, Phys.Lett. 18, 7 (1965); E.S.Mashkova and V.A.Molchanov, Rad.Effects 13, 183 (1972). J.F. Ziegler, J.P. Biersack and U. Littmark, The Stopping and Range of Ions in Solids (Pergamon, New York, 1985). M. Draxler, R. Beikler, E. Taglauer, K. Schmid, R. Grub er, S.N. Ermolov, and P. Bauer, Phys. Rev. A 68, 022901 (2003). T. Kobayashi, Nucl. Instr. Meth. B 249, 266 (2006). T.J. Wood, C. Bonet, T.C.Q. Noakes, P. Bailey, S.P. Tear, Surf.Sci. 598, 120 (2005). J.P. Biersack, E. Steinbauer, and P. Bauer, Nucl. Instr. Meth. B 61, 77 (1991). P. Bauer, E. Steinbauer, and J.P. Biersack, Nucl. Instr.


15
Meth. B 64, 711 (1992). [23] W.K. Chu, J.W. Mayer, and M.A. Nicolet, Backscattering Spectrometry (Academic, New York, 1978). [24] E.I. Sirotinin, A.F. Tulinov, V.A. Khodyrev, and V.N. Mizgulin, Nucl. Instr. Meth. B 4, 337 (1984). [25] Z. Smit, Phys. Rev. A 48, 2070 (1993). [26] E. Bњgh and J.L. Whitton, Phys. Rev. Lett. 19, 553 (1967). [27] This statement directly follows from the prop osition in probability theory known as the law of total exp ectation, see: P. Billingsley, Probability and Measure (John Wiley & Sons, New York, 1995).