Accurate Simulation of Channeling Implants Using Monte Carlo Techniques


A number of different codes have been developed for modeling implantation in crystal structures. The MARLOWE [1], UT-MARLOWE [2], and CRYSTAL-TRIM [3] codes have gained popularity and are used by many groups for studying ion implantation, sputtering and scattering from crystal surfaces. In ATHENA the current Monte Carlo models are based on PEPPER. However new models are currently implemented based on CRYSTAL, a BCA code developed at the University of Surrey, [4], which was written to explicitly meet the requirements for modeling of channeling implantation plus damage accumulation in crystalline silicon. These new models will be available in ATHENA and later in 3D in ODIN.

New Model Description

To calculate the penetration and depth distribution of bombarding particles, the model uses the well-known Binary Collision Approximation (BCA), in which, the deflection of the trajectories of moving particles is calculated in a strict binary way - between the moving atom and the closest atom in the lattice. The code has been developed extensively in two aspects: i) to account for changes in the target structure due to the ion bombardment and, ii) to account for ion implantation through two and three dimensional topography at the surface. Since the principles of BCA simulation have been described elsewhere [5] (and references therein), we will only briefly note the details here.

The ion implantation is simulated by following the fate of a large number of sequentially generated pseudo-projectiles, each of which carries an equivalent dose corresponding to a fluence increment obtained by dividing the total fluence by the number of pseudo-projectiles and scaled for the topography of interest.

After the collision cascade has been completed (i.e. the energies of all moving particles have dropped below some predetermined cut-off values), the probability P for each grid point is calculated, the average value of which is proportional to the damage buildup. In the subsequent collision cascade the choice of partner for interaction depends on the value of P. P=1 represents undamaged, perfect crystal.

The atomic collisions are considered to be composed of a quasi elastic part and of an essentially separate electron excitation part. The barycentric scattering angle is evaluated with the MAGIC formula, [6]. The present program uses the screened Coulomb potential with the universal screening function described in [6]. In order to simulate the channeling, the model uses non-local and local inelastic energy models. The local, i.e. impact dependent inelastic energy loss model is either Firsov's 5th power. [7], or Oen-Robinson's exponential model, [8]. Usually, for metals, a small amount of non-local energy loss is needed in order to approximate the electronic density in the middle of the channels. At present, our model uses the same inelastic model for silicon, although the picture there might be quite different for light ions and/or high energies where valence electrons play the primary role in inelastic energy losses.

To be able to make direct comparisons to the experiments, the model takes into account damage accumulation with the ion dose. The amount of damage, at the end of each cascade, recurrently depends on the previous one and the current density of newly created disorder.

The model is highly embedded into Athena, allowing for modeling of ion implantation through two dimensional surface structures and the possibility of performing ion implantation in targets with any initial surface topography. Also, re-implantation due to scattering from different surface planes any number of inclusions with different structure and density are automatically taken into account.


Careful examination of all available experimental data shows big discrepancies in depth profiles among different research groups. The main differences are due to three main factors:

i) implantation geometry (beam divergence,beam alignment along crystallographic channels) and proper
control to maintain it,
ii) the quality of the target surface,
iii) the depth profile measurements.

Any kind of consolidation between experiment and simulation should take into account these factors (a good example of how the misalignment with the channel changes the profiles is given in [9]). Therefore we have selected published experiments with estimated misalignment, beam divergence and quality of the surface.

The figures are for simulation of ion implantation under channeling conditions. This is a challenging problem since channeled profiles are very sensitive to the electronic energy losses and the damage model, and any discrepancy between simulation and experiment indicates the inapplicability of the inelastic loss or damage models chosen. The choice of surface orientation, angle of bombardment and beam divergence depends on the setup in the available experimental data. Universal screening function, Oen-Robinson model local energy loss and dynamic accumulation of the damage were used throughout the simulations.


Figure 1. Phosphorus depth profiles for implantations at 100 KeV in <100> direction. The experimental profiles are collected from [10] (experiments were performed with VARIAN/EXTRON 220 implanter, 'well channeled' conditions were maintained across the wafer by comparing SIMS depth profiles from different spots). The simulation is with 0 tilt, 0.5 beam divergence and 16 Å surface oxide layer. Different curves represent different fluences
in the range of 10 cm to 1015cm.



Figure 2. Surface oxide thickness dependence of Phosphorus depth profiles.
Simulation conditions as in Figure 1, the ion dose for all profiles is 10 cm.
Currently the SVDP models do not include screen oxide dependence for Phosphorus. This model will be used to extend the SVDP tables in the absence of experimental data.


Figure 3. Low energy (5 Kev) Boron implant into crystalline silicon.
The calculated profile is compared with Sims Validated Dual Pearson model [11]
available in ATHENA. The simulation conditions are 0 tilt and 10Å screening oxide layer.




[1] M.T. Robinson and I.M. Torrens. Physical Review B, Vol.9, pp.5008-5024, 1974.

[2] S.J. Morris, B. Obradovic, S.-H. Yang, and A.F. Tasch. IEDM Technical Digest 1996

[3] M. Posselt. Radiation Effects and Defects in Solids, Vol.130-131, pp.87-119, 1994.

[4] I. Chakarov and R. Webb. Radiation Effects and Defects in Solids, Vol.130-131, pp.447-452, 1994.

[5] W. Moeller. Simulation of Stopping and Ranges , volume 155 of NATO ASI Series E: Applied Sciences. Kluwer Academic Publishers, Dordrecht, 1989.

[6] J.F. Ziegler, J.P. Biersack, and U. Littmark. The Stopping and Range of Ions in Solids, volume 1 of The Stopping and Ranges of Ions in Matter. Pergamon Press, Oxford, 1985.

[7] O. B. Firsov. Sov. Phys. - JETP, Vol.36, p.1076, 1959.

[8] O.S. Oen and M.T. Robinson. Nuclear Instruments and Methods , Vol.132, pp.647-653, 1976.

[9] E.H.A. Dekempeneer, P.C. Zalm, G. van Hoften, and J. Politiek. Nuclear Instruments and Methods}, Vol.B48, pp.224-230, 1990.

[10] R. J. Schreutelkamp, V. Raineri, F. W. Saris, R. E. Kaim, J. F. M. Westendorp, P. F. H. M. van der Meulen, and K. T. F. Janssen. Nuclear Instriments and Methods, Vol.B55, pp.615-619, 1991.

[11] Al F. Tasch, H. Shin, C. Park, J. Alvis and S. Novak, J. Electrochem. Soc., Vol. 136, pp 810-814, 1989.