Monte Carlo Ion Beam Simulations using the Program FLUX.

A talk presented at the RENIBEL workshop on ion beam techniques (Lisbon, 10 March 2004)

Historical.

Channeling and computer simulation have gone hand in hand from the very start.  In fact the discovery of the channeling effect was due to a computer simulation, by Oen an Robinson in the 60's, when an unexplained tail appeared in implantation depth profiles simulated  with their computer code MARLOWE.

Different simulation programs were developed in the 70's. Most of these were leaning heavily on the "continuum string" approximation introduced by Lindhard, which allows an analytical approach. Of course this analytical approach saves time, but it is not always clear what is the effect of the approximations involved.

The first succesful Monte Carlo program specifically aimed at the simulation of channeling phenomena was Barrett's program LAROSE. This program was considered the standard for a long time. In fact program FLUX might have never been developed if John Barrett had not been so reserved in giving away his code.

FLUX versions.

The initiative for developing a new Monte Carlo program was taken by Dik Boerma in the early 80's. In the first version, FLUX1, a strictly Monte Carlo treatment was limited to one row of atoms. The effect of surrounding atom strings was taken into account in the continuum string approximation. This was also the approach as used by Barrett. In subsequent versions, FLUX3 and later, more atom rows were given the Monte Carlo treatment with discrete binary collisions, and only atoms more then about 1 to 2 Å away from the beam ion were still treated as continuum strings. This greatly improved the accuracy. Thus the continuum string approximation only remained as a more or less negligible correction to the strict Monte Carlo treatment.

intermezzo

That this correction was neglible was proven a few years ago. A very serious bug was discovered  in the calculation of the correction, that had gone unnoticed for over 10 years. After fixing the bug, no noticeable differences were found when several old simulations were repeated.

Some people have developed special versions of FLUX. I should mention here the work of Luis Rebouta, who extended the use to multiple atomic species such as lithium niobate. His program is called FLUX4, I am curious if it has been used more often after the original LiNbO3 work. Ulrich Wahl has worked out the application to the channeling of alpha particles emitted by a radioactive nuclide. This approach used time reversal: the final coordinates of the emitted alpha particle are fed as the initial coordinates to FLUX and collisions now result in an energy gain rather then an energy loss.
Note added July 2004: all 3 of these versions are now integrated in flux7.8

Time and money.

 It should also be noted that the spread of Monte Carlo simulations has only been possible due to the rapid development in computer performance. Before the 80’s things like that were only feasible using a main-frame (of which there was typically one per university, and you had to compete with chemists and astronomers). You handed in your deck of cards at the computer center, and came back the next day to see if the program had been run already. Around 1980 ‘powerful’ computers were also introduced in the laboratory. The first version of FLUX was developed on the Perkin Elmer 3220, called a minicomputer, but filling up half a room.  After that the development went fast. For a typical job the times involved are (also shown a rough guess of the price of the computer):

1985	Perkin Elmer 3220	days	€100 000
1990	HP720 Unix workstation	10 hours	€ 10 000
1997	HP C180 Unix workstation	2.5 hours	€ 10 000
2004	3 GHz PC Pentium4, Linux/gcc3.3	20 minutes	€   1000

The physical model.

binary collisions.

The effect of one collision is a change in the linear-momentum vector. This change is assumed to take place at one spot in space. In between collisions the trajectory is assumed to be a straight line. The locations where the collisions are thus assumed to take place are in planes spaced at regular intervals along a main axis of the crystal.

the potential..

The potential from which the force is derived is a screened Coulomb potential. The program has 3 different options for the screening function. The Molière potential, the “universal potential” of Ziegler Biersack and Littmark, further called the ZBL potential, and a potential based on Hartree-Fock charge distributions, further called HF potential.

The general features of the potential for ions in the MeV range may be seen from an example illustrated in the figure. The bulk of the interaction takes place in a very small region a few tenth’s of Ångstroms around the atom.

At 2 Å the scattering angle is of the order 10^-5 degrees, and this may explain why the above-mentioned bug, or the field of surrounding atom rows in general, has little influence.

Energy loss.

The energy loss due to interactions with electrons is split into two parts: localized, impact parameter dependent, interaction with inner shell electrons,  and,  interactions with outer shell electrons, assumed to be uniformly distributed. Both treatments in FLUX are maybe old-fashioned,  and deserve improvement, especially if someone ever would use it to analyse an experiment where the electronic energy loss is of critcal importance.

Nuclear energy loss

is, of course, automatically done by the kinematics of the binary collisions.

Thermal vibrations

are accounted for by random sampling the position of the atom from a gaussian distribution at each collision event. This implies that correlations in the vibrations of adjacent atoms are neglected. Barrett and Jackson have developed a method to incorporate correlated vibrations. This method has been incorporated in FLUX1, the version for monoatomic crystals. However this method has never been worked out in detail for multi-atomic crystals, and is not available anymore in the current version of FLUX.

The philosophy of the program.

The quantities of interest, i.e. the position and velocity coordinates are updated at regular intervals, with a uniform spacing along the z-axis, i.e. a main axis of the crystal. All quantities that the user specified for output are also sampled at these planes. One of these quantities, essential for the application to lattice site location, is the flux of ions, integrated over the depth. This is where the name of the program stems from.

The ion is allowed to undergo collisions with any atom present in the plane perpendicular to the z-axis. Than the energy and the new direction vector are calculated, and than on to the next plane.  For this purpose the choice of crystal lattice and of the z-axis have to be specified. This is done in file COMBIN.FIG. Following is a detail of this file. This specimen, and some of the following pictures ar for the case ot the <0001> axis of 4H-SiC. I just show this as an illustration, to get the idea that it is quite a precise job to set up such a configuration.

 /   <0001> STRING of the 4H-SiC polytype

 Z=<0001> AXIS, Z-X PLANE: (1 0 -1 0),       Y-Z PLANE: (1 -2 1 0)

 X-unit=a*sqrt(3)/4

 Y-unit=3a/4

 Z-unit=c/16

 In array Iatom:

 'A' planes have index 10(1*6 + 4), 'B' 24(3*6+6), 'C' 17 (2*6+5)

 sequence : A B A C

 make sure to count 'A' rows twice as neighbouring atom rows!

****4HSIC0001 hexagonal

  1,  0,  0,        1,  -2,  0

  6,

  2,  2,  16,

5.333333 5.333333

12, 12,

60001

-1,-0.3333333, -1,0.3333333,  -1,1, -1,1

 2, 1.3333333,  2,0.6666667,   2,0, 2,0

 1,1.6666667   0,-0.6666667

-1,1.6666667   2,-0.6666667

-1,-0.3333333, -1,0.3333333,  -1,1, -1,1

 2, 1.3333333,  2,0.6666667,   2,0, 2,0

 1,1.6666667   0,-0.6666667

-1,1.6666667   2,-0.6666667

0.433012701 0.0

0.75, 0.0

0.0,  0.0625

16 16

10  0  0  0 24  0  0  0 10  0  0  0 17  0  0  0

 0  0  0 10  0  0  0 24  0  0  0 10  0  0  0 17

0 0

1 0.333333333

0 0.666666667

1 1

0, 1.3333333

1, -0.333333333

What these numbers do is define a unit cell, the location of the atom rows with respect to this cell, and the location of the atoms along the rows. They also define a set of “surrounding strings”. Fortunately we have a little program, plotfig, to visualize this structure. For the case above this gives the following picture (to the left).

The picture to the right was borrowed from the thesis of Sang-Kwon Lee, http://www.lib.kth.se/Sammanfattningar/lee020508.pdf

The unit cell.

The unit cell is indicated by the small rectangle. In the calculation of the ion trajectory we restrict ourselves to one such cell. If the ion leaves the unit cell, it is put back at an equivalent position in the original cell by a suitable symmetry operation, i.e. a translation of the coordinates of the ion and/or a reflection of its velocity components. This implies, that  the ion flux distribution in the XY plane, collected during the simulation, is actually an average over cells of different orientations!! However, as long as these orientations are equivalent for our purpose, the location of impurity atom sites, this does not matter. Meanwhile, the true coordinates and velocity components are also kept track of, and are used, e.g. for the output of exit coordinates and velocity distributions. The initial coordinates are sampled over a true primitive cell  that comprises all possible orientations of the unit cell. In the current example 4H-SiC <0001> this means 4 ‘unit cells’.

The user may also optionally specify the size of the beamspot instead. This might be a larger  area, as in the example that follows,  but also even half a unit cell!

Note that this unit cell concept does not imply that the binary collisions are limited to atoms within the unit cell. In this particular case we selected 6 atom rows (painted red), 2 of which outside the unit cell, for the binary collisions. The remaining rows are the surrounding strings whose contribution, as argued before, is probably close to negligible.

Example: transmission channeling

One of the options of FLUX is to output the position and velocity coordinates after passing through a layer of the crystal. The next figure shows the result of this, for the SiC case, with a beam of 1.2 MeV alpha’s entering along the 0001 axis. The “grey scale” is such that high beam intensity is shown dark. We see, that the beam is steered very nicely into the empty regions between the atom rows. The atom rows appear in this picture as white regions. Note that not all white regions have the same area. The reason for this is that the rows labelled “A” in the picture above, have twice the density of atoms compared to the B and C rows.

A giant focusing peak

This an example of an RBS application. One time, when doing RBS/channeling on Si we just hit by chance on this peculiar phenomenon. At a beam angle about 1 degree away from a main axial channeling direction a huge peak appeared in the spectrum. Although at first we thought this was some instrumental effect, it was fully confirmed by Monte Carlo simulations. The following 2 pictures show the yield (close encounter probability) as a function of depth from a FLUX simulation, and an RBS spectrum constructed from the FLUX output quantities yield, energy loss, and energy straggling.

The explanation of this phenomenon becomes evident if you look at the path of the ions inside the lattice. The following was constructed using FLUX option XYZ  , and processing the result with program track .

  This is the XY projection of a few hundred ion histories followed for the first 300 Å. All histories start in a rectangle in the lower left corner, with the same transverse velocity, directed in the (211) plane.

It is seen that the deflection due to atom rows happens to focus the ions just onto a subsequent row.

Some details: 1 MeV He+ ions in (211) plane of Si 1 deg away from <110>,

P.J.M.Smulders, A.Dygo, and D.O.Boerma, Nuc Instr and Meth B67 (1992) 185

The righthand picture  is a detail of the paper,  P.J.M.Smulders, A.Dygo, and D.O.Boerma, Nuc Instr and Meth B67 (1992) 185, and shows that the correspondance with experiment is remarkable

It turned out this effect was very sensitive to the potential.

The solid line is calculated with the Hartree Fock potential, the dotted line with the ZBL potential.

Yes, the line is a simulated RBS spectrum, not just  a guide to the eye!

Lattice site location.

After giving you a small taste of things you can do with FLUX, now comes the main course.

Lattice site location of impurity atoms in a single crystal using RBS/channeling is done by measuring a set of angular scans, both for the host atoms and for the impurity atoms. Using incident beam directions in and around several different main axial directions one hopes to get a discrimination between different possible impurity sites. The analysis involves comparing the measured angular scans with predicted ones. This means that we do the FLUX Monte Carlo simulations as much as possible for the same situation as the experiment,  same energy, beam direction with respect to the lattice, and so on. The fact that a beam direction is not defined by the angle with the main axis alone is often forgotten in an experiment. Also an azimuthal angle is required, and may be of utmost importance, as we see for instance from the focusing effects described above. So we assume the experimenter knows what he is doing and also determines the azimuthal angle as accurately as possible. A different analysis procedure is carried out for the host atoms and the impurity atoms.

Host atoms.

The yield, or rather the nuclear encounter probability, as a function of depth, is a standard quantity determined in the FLUX simulations. We may integrate this over some depth interval (program PRAAK) and subsequently compare this to the experimental yield over the same depth interval. This assumes that the depth-energy relation for the RBS spectrum, is constant for a scan over a small angular region. But, channeling and non-channeling directions have different energy loss! A better method is to simulate the RBS spectrum for each of the points of the scan, and then compare the yield in the same energy region of simulations and experiment.

Impurity atoms.

Here a different approach is followed. The simulation program produces a weighted flux. By this I mean a convolution of the flux density ρ(x,y,z) with the depth distribution w(z), which is stored in a 2-dimensional array f(x,y). This quantity is subsequently used in a second program named YIMP to calculate the yield for the impurity atoms. This involves convoluting the distribution of assumed lattice sites taking into account their thermal vibrations) with the weighted flux. The advantage of this approach is, that we have to do the time-consuming Monte Carlo simulations only once. For each of the different possible sets of impurity sites we want to compare we only have to run YIMP, which takes a lot less time. since no Monte Carlo simulation is involved here, but only analytical expressions.

The whole procedure is schematically represented in the following slide:

As an example, we will look at the thesis work of Wouter Segeth (Groningen 1987). The work concerned the trapping of oxygen in silver by 5sp impurities e.g. Sb or In. The Sb or In atoms are first built into the Ag lattice by ion implantation followed by annealing. Subsequently an oxidation process is carried out under well-defined circumstances. We want to study the microscopic structure of the impurity-oxygen complexes thus formed inside the silver. Ion beam analysis may give information about the depth profile as well as the lattice location of the trapped oxygen. Mossbauer spectroscopy and Perturbed Angular Correlation measurements were also carried out to study the electromagnetic fields at the impurity sites. 

Depth profiling, as well as channeling, of the O were both done by measuring alpha particles from the ¹⁸O(p,α) ¹⁵N reaction.

For Ag crystals, implanted with Sb, and oxidized at 550K, it was found from the depth profiles that precisely 2  oxygen atoms are present for every antimony atom. From the channeling measurements it was found that the oxygen  is localized at octahedral interstitial sites in the silver.

For Ag crystals, implanted with indium atoms, a more complex situation was found.  Two different complexes are formed, depending on the oxidation treatment. At low temperature, 550K, a complex is formed that resembles closely the result for Sb, where the O atoms were interstitial. However now the oxygen atoms are slightly displaced from the insterstitial site, resulting in a lower fluxpeak. At higher temperatures, 700K, the O atoms shift to substitutional lattice sites

The next figure shows the model for the isolated SbO2 atoms in Ag:

Channeling can not distinguish between the two possible configurations, but the stretched configuration is energetically favoured. In channeling the O atoms will be obscured by the Ag atoms for the <100> and <111> incident directions, but not for the <110> direction. The following 3 figures show the projected lattice for these 3 beam directions:

The circles are the lattice rows of the host Ag, the stars are the rows on which the O atoms may be found. Different planes perpendicular to the projections are indicated.

Some measured angular scans, together with fitted scans, calculated from Monte Carlo simulations, follow, first for the Sb in Ag, oxidized at 550 K.

The curves in the left half are for the Ag host atoms, measured by RBS

The curves in the right half are for O, measured using the ¹⁸O(p,α) ¹⁵N  reaction.

Incident main axes <110>, <100>, and <111>, respectively, as indicated.

The fluxpeaking effect in the <110> direction is even stronger then the Monte Carlo prediction (notice the lonely data point in the top of the picture!)

Below are results for InO2 in Ag, on the left side after an oxidation treatment at 550K, on the right after heating to 700K.

At 550K the results are similar to the previous case, although the flux peak is less intense, showing that the In have shifted from the interstitial sites. From the analysis it was found that this shift was 0.5 Å.

Flux in space.

Many of the successes of FLUX have been due to the cooperation with Mark Breese and his co-workers.

Many options in the program, such as multilayer samples and bent crystals are the result of this cooperation.

Some subjects of investigation over the years:

Transmission ion channeling images of crystal defects

interesting behaviour of stacking faults in Si: blocking to channeling transition

Study of Si(1-x)Ge(x) layers grown on Si substrate: direct determination of lattice strain

Ion beam bending using strained bilayers

Planar channeling oscillations: phase space analysis

I will now go into the last subject, phase space analysis of planar channeled ions. Let me first remark that for FLUX planar channeling does not exist as such. The binary collisons of each individual ion with lattice atoms are treated exactly the same whether or not the ion is incident in a planar direction. The collective effects of atoms arranged in a plane thus appear automatically, and are in no way specified in the input to the program. Yet, as we shall see, the collective results that appear may be very well interpreted in terms of a planar potential, and, in terms of a phase space where the coordinates are the distance from the plane, and the velocity component perpendicular to the plane.

It is a well known fact that the planar potential resembles a harmonic potential. If there is no energy loss then the orbit in phase space of a particle in such a potential is an ellipse. The particle moves through this ellipse with a fixed period, independent of the initial values of coordinate and linear momentum. This is the same as saying that the trajectory of the ion exhibits planar oscillations with a period (or equivalently, wavelength) independent of amplitude.

Consider for instance 3 MeV protons channeling in the {110} plane of silicon. We may look at the process by plotting a few 100 ions tracks in the zy plane:

The left picture shows the tracks for incidence in the planar direction, the right picture for a tilt of 0.06 °.

The striking difference is that the planar oscillations seems to be preserved much better at the larger tilt angle.  This is only apparently so. In the left picture we see contributions of all amplitudes superimposed, in the righthand picture there is only a narrow band of amplitudes. If we plot the situation at z=0 in phase space, we get something like the following picture:

For an ideal harmonic oscillator potential we would now, as a function of time, see this line rotate about the origin. At a depth of about 600 Å we see the first “waist” appear. We are here at λ/4, so the line is now in its vertical position. The spread in position is minimal here, and the spread in angle to maximal. Monte Carlo simulation of the phase space plot at this point gives qualitatively the expected result:

From these results, it is now simple to predict how the corresponding plots for a tilt of 0.06° would look:

The main difference with the previous case of tilt 0 is now that the plot at depth 0 is lifted vertically, while at a depth of   λ/4 it is shifted to the right.

 Mark Breese has sent me some “fluxmovies” that he made with Mukhtar Rana, one of his students, that illustrate the above.MOVIES\index.html

Planar channeling in buried nano-layers

One of the striking features of planar channeling, as we have seen above, is that the coherence of the ion trajectories persists for a large depth. This property was used in the thesis work of Louis Selen (Eindhoven , Leo van IJzendoorn’s group, 2001) in the study of buried strained films of only a few nanometers thick, of Si_1-xGe_x in Si. In the top of the picture we see schematically the sample, the nanolayer, and a capping layer.

The rotation of the {011} lattice planes at the top and bottom interfaces of the epi-layer has as a net effect a translation of these lattice planes. (The thickness of the epi-layer is only a few nm so the rotation itself is hardly felt by the ion beam). This shift causes a step in the channeled RBS spectrum. The step height depends on the magnitude of the translation, and, on the phase of the channeled ions. Thus, in general, we expect the step height as a function of incident angle of the ion beam to be asymetric. The shape of this curve depends not only on the magnitude of the translation, but also on the thickness of the capping layer. Furthermore it requires a knowledge of the wave length of the channeling oscillations, which is quite sensitive to the ion-atom potential. An example is shown in the next figure:

We see a very big difference between the behaviour for different potentials,  with the HF potential again being favoured.

The maximum asymmetry occurs for a tilt angle of 0.1°. From Monte Carlo simulated phase space plots, below, we may conclude that the main effect of taking different potentials is the wave length and thus the phase at the location of the buried layer.

Altogether we have a complicated set of correlated parameters: properties of the potential used, as well as properties of the sample. Therefore it is desirable to have independent determinations of some of these properties by other methods. One such method is axial channeling. In axial channeling we also have a dechanneling effect due to the lattice translation, but here the collective effects of the beam trajectories are absent at such a large depth, and the dependence on the choice of potential is much smaller. The observed step size in RBS spectra of axial channeling is now much more symmetric with repect to the incident beam angle.

In other words, although the step size in the axial channeling spectra is much smaller, this is compensated by the fact that the results are less model-dependent.

Final remarks.

I have shown a few possible applications of the FLUX package. Among its capabilities are.

RBS channeling spectra

Lattice site location of impurity atoms

Transmitted ion angular distributions and energy loss distributions

mixed crystals

multilayer samples: lattice rotations and translations

curved crystals

phase space analysis

emission channeling (Ulrich Wahl)

Early versions of FLUX were limited to very small angles between the ion trajectory and a main crystal axis, and still depended to some extent on the continuum string model. These limitations have been removed. Also, energy dependent energy loss of the beam ion was introduced in a recent version.

Some limitations that remain are:

The possibilies to treat lattice damage are limited. Basically FLUX assumes a perfect crystal. Lattice damage may to some degree be simulated by an increased thermal vibration amplitude. It is also possible to introduce voids.

The treatment of the electronic energy loss processes might be improved.

The strategy of program FLUX is to look for the next collsion partner in planes at fixed intervals. An alternative is to search for the next collision partner in 3-dimensional space. Such an approach has the advantage that it allows more versatility in de structure of the sample and the beam direction. A search algorithm of this type is not present in FLUX.

In the calculation of RBS channeling spectra the Monte Carlo approach is used only for the incoming trajectories, while the outgoing ions are assumed to experience no further channeling or focusing, and have ‘random’ energy loss. The increase of computer performance in the last decennia makes it feasible to treat both in- and outgoing trajectories in the binary collision approach. Dik Boerma   (D.O.Boerma, private communication) has told me that Vassily Khodyrev and he have just completed  such a program. It uses clever strategies for finding the next collision partner, and  ‘smart sampling’ of collision events, and follows ingoing and outcoming trajectories without approximations except for the models for the potential, thermal vibrations, and such. Maybe this will be the next generation of ion beam simulation programs.

Note added 2011: A description of the program mentioned above, called TRIC, may be found at http://arxiv.org/PS_cache/arxiv/pdf/0904/0904.2151v1.pdf
For information, mail to:
Khodyrev (main author)
Roch (user)

It is now 20 years ago since Dik Boerma and I started developing the program FLUX. We had no idea then that the result would still be going strong in the 21nd century!