Kai Nordlund
http://www.helsinki.fi/~knordlun
Related info in http://www.helsinki.fi/~knordlun/range.html
(Last updated Sep 30, 1996)
The energy loss of particles in matter is of considerable interest in several different branches of physics. Ion implantation is a widely used technique in for instance semiconductor fabrication and different branches of materials science.
At very high energies (larger than several hundred MeV) the slowing down of all charged particles is contributed to by bremsstrahlung, Cherenkov radiation and nuclear reactions [1]. Since we are primarily interested of lower energies, we will not discuss these processes in greater detail.
At lower energies, the slowing down of ions is traditionally separated into to distinct processes: electronic and nuclear slowing down or stopping power. The sum of these two processes is called the total stopping power
dE
S(E) = --
dr
which tells how much the ion loses energy on a certain distance it travels in a medium [2].
With electronic stopping one means slowing down due to the inelastic collisions between electrons in the medium and the ion moving through it. The term inelastic is used to signify that the collisions may result in excitations in the electron cloud of the ion; therefore the collision can not be treated as a classical scattering process between two charged particles.
Since the number of collisions an ion experiences with electrons is large, and since the charge state of the ion while traversing in the medium may change frequently, it is very difficult to describe all possible interactions for all possible ion charge states. Instead, the electronic stopping power is given as a simple function of energy Se (E) which is an average taken over all energy loss processes for different charge states. It can be theoretically determined to an accuracy of a few % in the energy range above several hundred keV from theoretical treatments, the best known being the Bethe-Bloch formula [1]. At energies lower than about 100 keV it becomes very difficult to determine the electronic stopping theoretically. Hence several semi- empirical stopping power formulas have been devised. The by far most used today is the model given by Ziegler, Biersack and Littmark (the ‡ZBLƒ stopping) [2].
Figure 1. Ratio between nuclear and electronic stopping power. The maximum of the nuclear stopping curve typically occurs at energies between 10-100 keV, of the electronic stopping power at MeV energies. For very light ions slowing down in heavy materials, the nuclear stopping is weaker than the electronic at all energies.
With nuclear stopping one means elastic collisions between the ion and atoms in the sample. If one knows the form of the repulsive potential V (r) between two atoms, it is possible to calculate the nuclear stopping power Sn (E). Nowadays, the repulsive potential between two atoms can be obtained to great accuracy from ab initio quantum mechanical calculations (cf. the next section). In fig. 1 the typical ratio between nuclear and electronic stopping powers is illustrated.
1 Z_1 Z_2
V(r) = -------------- ------- phi(r),
4 pi epsilon_0 r
where phi(r) -> 1 when r -> 0.
Here Z_1 and Z_2 are the charges of the interacting
nuclei, and r the distance between them.
A large number of different repulsive potentials and screening functions have been proposed over the years, some determined semi-empirically, others from theoretical calculations. A much used repulsive potential is the one given by Ziegler, Biersack and Littmark, the so called ZBL repulsive potential. It has been constructed by fitting a universal screening function to theoretically obtained potentials calculated for a large variety of atom pairs [5]. The ZBL screening function has the form
-3.2x -0.9423x -0.4029 -0.2016 x
phi(x) = 0.1818 e + 0.5099 e + 0.2802 e + 0.02817 e
where x = r / a_u and
0.8854 a_0
a_u = -------------------
Z_1^0.23 + Z_2^0.23
(a_0 is the Bohr atomic radius = 0.529 Å).
The standard deviation of the fit of the universal ZBL repulsive potential to the theoretically calculated potentials is 18 \% above 2 eV [5]. In cases where experimental range profiles can be determined with a good depth resolution, it is crucial to have a repulsive potential with a better accuracy than that given by the ZBL potential.
A more accurate repulsive potential can be obtained from self-consistent total energy calculations using density-functional theory and the local-density approximation (LDA) for electronic exchange and correlation [7,8]. In this approach the total energy, including both the electronic part and the inter-nuclear Coulomb repulsion, is obtained as a function of the distance between atoms in a dimer bond. In practice, we have used the DMol package well tested in calculations of energetics and structures of small molecules to obtain the interatomic potentials.
In the beginning of the slowing-down process at high energies, the ion slows down mainly by electronic stopping, and moves almost in a straight path (although, for very light ions like H and He the occasional nuclear collisions with heavy target atoms are so strong that the path can be better likened with a random walk). When the ion has slowed down sufficiently, the collisions with nuclei become more and more probable, finally dominating the slowing down. When atoms receive significant recoil energies, they will be removed from their lattice positions, and produce a cascade of further collisions in the lattice. These collision cascades are the main cause of damage production during ion implantation. The process is schematically illustrated in fig. 2.
Figure 2. Illustration of the slowing down of a single ion in a solid material. The inset shows a typical range distribution. The case shown here might for instance be the slowing down of a 1 MeV silicon ion in silicon. The mean range for a 1 MeV ion typically is about 1 - 2 µm.
From the total stopping power it is possible to calculate a value for the mean range by straightforward integration [2],
E
/ dE
R(E) = | ----
/ S(E)
0
Since the total stopping power S(E) is just the sum of the electronic and nuclear stopping
powers, a value for the range can be obtained when both the electronic and nuclear stopping
powers are known. However, it turns out that a mean range value obtained in this way can easily
be false by a factor of 2 or even more.
This is due to the fact that the simple stopping power treatment is very inaccurate in describing the
nuclear slowing down. In reality, the collisions an ion experiences may be strongly correlated,
especially in crystalline materials. This leads to the need to use computer simulations when
calculating the slowing down of an ion in a medium.
Computer simulation methods to calculate the motion of ions in a medium have been developed since the 1960's. The basic idea in them is to follow the movement of the ion in the medium by simulating the collisions with nuclei in the medium. The electronic stopping power is usually taken into account as a frictional force slowing down the ion.
Conventional methods used to calculate ion ranges are based on the binary collision approxi- mation (BCA) [3]. In these methods the movement of ions in the implanted sample is treated as a succession of individual collisions between the recoil ion and atoms in the sample. For each individual collision the classical scattering integral
infinity
/ p dr
|
Theta = pi - 2 | ---------------------------------------
| _______________________
/ 2 / V (r)
r r _ / 1 - ----------------
min \/ E_c - p^2 / r^2
is solved by numerical integration. Here
Theta is the scattering angle in center-of-mass coordinates, E_c the center
of mass energy, p the impact parameter of the collision and
V(r) is the interatomic repulsive potential.
The impact parameter p in the scattering integral is determined either from a stochastic distribution or in a way that takes into account the crystal structure of the sample. The former method is suitable only in simulations of implantation into amorphous materials.
The best known BCA simulation program is TRIM [4,5] (short of TRansport of Ions in Matter), which is based on the ZBL electronic stopping and interatomic potential. It has a very easy- to-use user interface, which has made it immensely popular. More info on TRIM is available in http://www.helsinki.fi/~knordlun/TRIM-answer in question-and-answer form. However, it doesn't take account of the crystal structure, which severely limits it's usefulness in many cases. Several BCA programs overcome this difficulty; some fairly well-known are MARLOWE [6], BCCRYS and crystal-TRIM.
Although the BCA methods have been successfully used in describing many physical processes, they have some obstacles for describing the slowing down process of energetic ions realistically. Due to the basic assumption that collisions are binary, severe problems arise when trying to take multiple interactions into account. Also, in simulating crystalline materials the selection process of the next colliding lattice atom and the impact parameter p always involve several unphysical parameters, which may affect the results 10-20 % even for quite reasonable-seeming choices of the parameter values. The best reliability in BCA is obtained by including multiple collisions in the calculations, which is not easy to do correctly. However, at least MARLOWE does this.
A fundamentally more straightforward way to model atomic collisions is provided by molecular dynamics (MD) simulations, in which the time evolution of a system of atoms is calculated by solving the equations of motion numerically. MD simulations are well suited for realistically calculating the time evolution of single full collision cascades at low (E < 10 keV) energies.
Unfortunately, full MD simulations are far too inefficient to calculate ion range distributions at keV or MeV energies. However, quite recently methods have been devised in which the number of interactions and atoms involved in MD simulations have been reduced in order to make them more efficient in calculating ion ranges. One such program is MDRANGE, which has been developed at our laboratory [7,8]. It overcomes several of the drawbacks of BCA methods, but is still efficient enough to obtain range distributions at keV energies in a reasonable amount (a few hours to a few days) of computer time.
In crystalline materials the ion may in some instances get "channeled", i.e. get focused into a channel between crystal planes where it doesn't experience almost any collisions with nuclei. Also, the electronic stopping power may be weaker in the channel. The crystalline BCA codes and the MDRANGE code automatically describe channeling effects. However, in order to take into account the weaker electronic stopping in simulations one must use a non-local electronic stopping power, i.e. one that depends on the strength of the collisions the ion experiences. These are still under development, but for instance MARLOWE has the possibility of using some model for it. The MDRANGE code does from V3.0 on have a non-local electronic stopping power. See the mdrange documentation for more info.
[2] J. Lindhard, M. Scharff, and H. E. Shi|tt. Range concepts and heavy ion ranges. Mat. Fys. Medd. Dan. Vid. Selsk., 33(14):1, 1963.
[3] M. T. Robinson and Ian M. Torrens. Computer simulation of atomic-displacement cascades in solids in the binary-collision approximation. Phys. Rev. B, 9(12):5008, 1974.
[4] J. P. Biersack and L. G. Haggmark. A Monte Carlo computer program for the transport of energetic ions in amorphous targets. Nucl. Instr. Meth., 174:257, 1980.
[5] J. F. Ziegler, J. P. Biersack, and U. Littmark. In The Stopping and Range of Ions in Matter, volume 1, New York, 1985. Pergamon.
[6] Mark T. Robinson. Computer simulation studies of high-energy collision cascades. Nucl. Instr. Meth. Phys. Res. B, 67:396, 1992.
[7] Kai Nordlund. Molecular dynamics simulation of ion ranges in the post-kev energy region. Comp. Mat. Sci., 3 (1995) 448-
[8] J. Keinonen, A. Kuronen, K. Nordlund, R. M. Nieminen and A. P. Seitsonen, First-Principles Simulation of Collision Cascades Si to Test Pair-Potential for Si-Si Interaction at 10 eV -- 5 keV, Nucl. Instr. Meth. Phys. Res. B 88 (1994) 382-.