# An Efficient Molecular Dynamics Scheme for the Calculation of Dopant Profiles due to Ion Implantation

## I Introduction

The principal reason for implanting ions<sup>*</sup><sup>*</sup>*Within this paper, *ion* is used to refer to the implanted species, and *atom* to refer to a particle of the target material; this has no implication to the charge state of either atom type. into silicon wafers is to dope regions within the substrate, and hence modify their electrical properties in order to create electronic devices. The quest for ever increasing processor performance demands smaller device sizes. The measurement and modeling of dopant profiles within these ultra shallow junction devices is challenging, as effects that were negligible at high implant energies become increasingly important as the implant energy is lowered. The experimental measurement of dopant profiles by secondary ion mass spectrometry (SIMS) becomes problematic for very low energy (less than 10 keV) implants. There is a limited depth resolution of measured profiles due to profile broadening, as the SIMS ion-beam produces ‘knock-on’s, and so leads to effects such as diffusion of dopants and mixing. The roughness and disorder of the sample surface can also convolute the profile, although this can be avoided to a large extent by careful sample preparation.

The use of computer simulation as a method for studying the effects of ion bombardment of solids is well established. Binary collision approximation (BCA), ‘event-driven’ codes have traditionally been used to calculate such properties as ranges of implanted species and the damage distributions resulting from the collision cascade. In this model, each ion trajectory is constructed as a series of repulsive two-body encounters with initially stationary target atoms, and with straight line motion between collisions. Hence the algorithm consists of finding the next collision partner, and then calculating the asymptotic motion of the ion after the collision. This allows for efficient simulation, but leads to failure of the method at low ion energies. The BCA approach breaks down when multiple collisions (where the ion has simultaneous interactions with more than one target atom) or collisions between moving atoms become significant, when the crystal binding energy is of the same order as the energy of the ion, or when the time spent within a collision is too long for the calculation of asymptotic trajectories to be valid. Such problems are clearly evident when one attempts to use the BCA to simulate channeling in semiconductors; here the interactions between the ion and the target are neither binary nor collisional in nature, rather they occur as many simultaneous soft interactions which steer the ion down the channel.

An alternative to the BCA is to use molecular dynamics (MD) simulation, which has long been applied to the investigation of ion bombardment of materials, to calculate the ion trajectories. The usefulness of this approach was once limited by its computational cost and the lack of realistic models to describe materials. With the increase in computational power, the development of efficient algorithms, and the production of accurate empirical potentials, it is now feasible to conduct realistic MD simulations. In the classical MD model, atoms are represented by point masses that interact via an empirical potential function that is typically a function of bond lengths and angles; in the case of Si a three-body or many-body potential, rather than a pair potential is required to model the stable diamond lattice and to account for the bulk crystal properties. The trajectories of atoms are obtained by numerical integration of Newton’s laws, where the forces are obtained from the analytical derivative of the potential function. Thus, MD provides a far more realistic description of the collision processes than the BCA, but at the expense of a greater computational requirement. Here we present a highly efficient MD scheme that is optimized to calculate the concentration profiles of ions implanted into crystalline silicon. The algorithms are incorporated into our implant modeling molecular dynamics code, REEDNamed for ‘Rare Event Enhanced Domain following’ molecular dynamics., which runs on many architectures either as a serial, or as a trivially parallel program.

## II Molecular Dynamics Model

The basis of the molecular dynamics model is a collection of empirical potential functions that describe interactions between atoms and give rise to forces between them. In addition to the classical interactions described by the potential functions, the interaction of the ion with the electrons within the target is required for ion implant simulations, as this is the principle way in which the ion loses energy. This is accomplished via a phenomenological electronic stopping-power model. Other ingredients necessary to the computation are a description of the target material structure and thermal vibration within the solid. It is also necessary to define a criterion to decide when the ion has come to rest in the substrate. We terminate a trajectory when the *total* energy of the ion falls below 5 eV. This was chosen to be well below the displacement threshold energy of Si (around 20 eV).

### A Empirical Potential Functions

Interactions between Si atoms are modeled by a many-body potential developed by Tersoff. This consists of Morse-like repulsive and attractive pair functions of interatomic separation, where the attractive component is modified by a many-body function that has the role of an effective Pauling bond order. The many-body term incorporates information about the local environment of a bond; due to this formalism the potential can describe features such as defects and surfaces, which are very different to the tetrahedral diamond structure.

ZBL ‘pair specific’ screened Coulomb potentials are used to model the ion-Si interactions for As, B, and P ions. Where no ‘pair specific’ potential was available, the ZBL ‘universal’ potential has been used. This is smoothly truncated with a cosine cutoff between 107% and 147% of the sum of the covalent radii of the atoms involved; the cutoff distances were chosen as they give a screening function that approximates the ‘pair specific’ potentials for the examples available to us.

The ZBL ‘universal’ potential is also used to describe the close-range repulsive part of the Tersoff Si-Si potential, as the standard form is not sufficiently strong for small atomic separations. The repulsive Morse term is splined to a shifted ZBL potential, by joining the two functions at the point where they are co-tangent. In the case of Si-Si interactions, the join is at an atomic separation of 0.69 Å, and requires the ZBL function to be shifted by 148.7 eV. The increase in the value of the short-range repulsive potential compensates for the attractive part of the Tersoff potential, which is present even at short-range.

### B Inelastic Energy Loss

The Firsov model is used to describe the loss of kinetic energy from the ion due to inelastic collisions with target atoms. We implement this using a velocity dependent pair potential, as derived by Kishinevskii. This gives the force between atoms $`i`$ and $`j`$ as:

$$𝐅_{ij}=\frac{2^{1/3}\mathrm{}}{2\pi a_\text{B}}(𝐯_j𝐯_i)\left(Z_1^2\times I\left(\frac{Z_1^{1/3}\alpha R}{a}\right)+Z_2^2\times I\left(\frac{Z_2^{1/3}(1\alpha )R}{a}\right)\right)$$
(1)   
where:

$`I(X)={\displaystyle _X^{\mathrm{}}}{\displaystyle \frac{\chi ^2(x)}{x}}𝑑x\text{, and }\alpha =\left(1+\left({\displaystyle \frac{Z_2}{Z_1}}\right)^{1/6}\right)^1`$ (2)   
and $`\chi (x)`$ is a screening function, $`Z`$ is atomic number ($`Z_1>Z_2`$), $`R`$ is atomic separation, and $`a=(9\pi ^2/128)^{1/3}a_\text{B}`$. For consistency with the ion-Si interactions, we use the ZBL ‘universal’ screening function within the integral; there are no fitted parameters in this model.

We have found that it is necessary to include energy loss due to inelastic collisions, and energy loss due to electronic stopping (described below) as two distinct mechanisms. It is not possible to assume that one, or other, of these processes is dominant and *fit* it to model all energy loss for varying energies and directions.

### C Electronic Stopping Model

A new model that involves both global and local contributions to the electronic stopping is used for the electronic energy loss. This modified Brandt-Kitagawa model was developed for semi-conductors and contains only one fitted parameter per ion species, for all energies and incident directions. We believe that by using a realistic stopping model, with the minimum of fitted parameters, we obtain a greater transferability to the modeling of implants outside the fitting set. This should be contrasted to many BCA models which require completely different models for different ion species or even for different implant angles for the same ion species, and that contain several fitted parameters per species.

Our model has been successfully used to describe the implant of As, B, P, and Al ions with energies in the sub MeV range into crystalline Si in the $``$100$``$, $``$110$``$, and non-channeling directions, and also into amorphous Si. While initially developed for use in BCA simulations, the only modification required to the model for its use in MD is to allow for the superposition of overlapping charge distributions, due to the fact that the ion is usually interacting with more than one atom at a time. The one fitting parameter is $`r_s^0`$, the ‘average’ one electron radius of the target material, which is adjusted to account for oscillations in the $`Z_1`$ dependence of the electronic stopping cross-section.

### D Structure of the Target Material

For the calculations presented here, the target is crystalline Si with a surface amorphous layer. The amorphous structure was obtained from a simulation of repeated radiation damage and annealing, of an initially crystalline section of material. Thermal vibrations of atoms are modeled by displacing atoms from their lattice sites using a Debye model. We use a Debye temperature of 519.0 K for Si obtained by recent electron channeling measurements. This gives an rms thermal vibrational amplitude in one dimension of 0.0790 Å at 300.0 K. Note, we do not use the Debye temperature as a fitting parameter in our model, as is often done in BCA models. The thermal velocity of the atoms is unimportant as it is so small compared to the ion velocity, and is set to zero.

At present there is no accumulation of damage within our simulations, as we wish to verify the fundamental model with the absolute minimum of parameters that can be fit. At a later date we will incorporate a statistical damage model into our simulations in a manner similar to that used in BCA codes. We also intend to include the capability of using amorphous, or polycrystalline targets in our simulations. 

## III Efficient Molecular Dynamics Algorithms

During the time that MD has been in use, many algorithms have been developed to enhance the efficiency of simulations. Here we apply a combination of methods to increase the efficiency of the type of simulation that we are interested in. We incorporate both widely used methods, which are briefly mentioned below, and new or lesser known algorithms for this specific type of simulation which we describe in greater detail.

### A Basic Algorithms

We employ neighbor lists to make the potential and force calculation O($`N`$), where $`N`$ is the number of particles. Coarse grained cells are used in the construction of the neighbor list; this is combined with a Verlet neighbor list algorithm to minimize the size of the list. Atoms within 125% of the largest interaction distance are stored in the neighbor list, which is updated only when the relative motion of atoms is sufficient for interacting neighbors to have changed.

### B Timestep selection

The paths of the atoms are integrated using Verlet’s algorithm, with a variable timestep that is dependent upon both kinetic and potential energy of atoms. For high energy simulations the potential energy as well as the velocity of atoms is important, as atoms may be moving slowly but have high, and rapidly changing, potential energies during impacts. The timestep is selected using:

$$\mathrm{\Delta }t_n=\frac{C_{DIS}}{\sqrt{\begin{array}{c}\mathrm{max}\\ 1iN\end{array}\left(\frac{2\times \left[KE_i+\mathrm{max}(0,PE_i)\right]}{M_i}\right)}}$$
(3)   
where $`KE_i`$, $`PE_i`$ and $`M_i`$ are the kinetic energy, potential energy and mass respectively of atom $`i`$, and $`C_{DIS}`$ is a constant with a value of 0.10 Å. Away from hard collisions, only the kinetic energy term is important, and the timestep is selected to give the fastest atom a movement of $`C_{DIS}`$ in a single timestep. When the timestep is increasing, it is limited by:

$$\mathrm{\Delta }t_n^{}=\mathrm{min}(1.05\times \mathrm{\Delta }t_{n1},\frac{3}{4}\mathrm{\Delta }t_{n1}+\frac{1}{4}\mathrm{\Delta }t_n)$$
(4)   
to prevent rapid oscillations in the size of the timestep, and the maximum timestep is limited to 2.0 fs.

The timestep selection scheme was checked to ensure that the total energy in a full (i.e., without the modifications described below) MD simulation was well conserved for any single ion implant with no electronic stopping; e.g. in the case of a non-channeling (10 tilt and 22 rotation) 5 keV As ion into a 21168 atom Si{100} target, the energy change was 3.6 eV (0.004%) during the 250 fs it took the ion to come to rest.

### C Domain following

Even with the computation resources available today it is infeasible to calculate dopant profiles by full MD simulation. Although the method is O($`N`$) in the number of atoms involved, the computational requirements scale extremely quickly with the ion energy. The cost of the simulation can be estimated as the number of atoms in the system multiplied by the number of timesteps required. Consider the case of an ion, subject to an energy loss proportional to its velocity, $`v(t)`$, which is then given by $`v(t)=u\mathrm{exp}(\alpha t)`$ where $`u`$ is its initial velocity and $`\alpha `$ is the loss coefficient. Each dimension of the system must scale approximately as the initial ion velocity, $`u`$, to fully contain an ion-path. If the timestep size is chosen so that the maximum distance moved by any particle in a single step is constant, the number of timesteps is approximately proportional to the ion distance.

Hence the method is roughly O($`u^4`$). Although it is possible to compute a few trajectories at ion energies of up to 100s of keV, the calculation of the thousands necessary to produce statistically reliable dopant profiles is out of the question. Therefore, we have concentrated on developing a restricted MD scheme which is capable of producing accurate dopant profiles with a much smaller computational overhead.

As we are only concerned with the path of the implanted ion, we only need to consider the region of silicon immediately surrounding the ion. We continually create and destroy silicon atoms, to follow the domain of the substrate that contains the ion. Material is built in slabs one unit cell thick to ensure that the ion is always surrounded by a given number of cells on each side. Material is destroyed if it is outside the domain defined by the ion position and the domain thickness. In this scenario, the ion sees the equivalent of a complete crystal, but primary knock-on atoms (PKAs) and material in the wake of the ion path behave unphysically, due to the small system dimensions. Hence we have reduced the cost of the algorithm to O($`u`$), at the expense of losing information on the final state of the Si substrate. This algorithm is similar to the ‘translation’ approach used in the MDRANGE computer code developed by Nordlund. The relationship between the full and restricted MD approaches is shown in Fig. 1.

Fig. 2 illustrates a single domain following trajectory. The ion is initially above a semi-infinite volume that is the silicon target. As the ion approaches the surface, atoms begin to be created in front of it, and destroyed in its wake. This process is continued until the ion comes to rest at some depth in the silicon substrate. Several thousand of such trajectories are combined to produce the depth profile of implanted ions.

### D Moving Atom Approximation

This was first introduced by Harrison to increase the efficiency of ion sputtering yield simulations. In this scheme atoms are divided into two sets; those that are ‘on’ have their positions integrated, and those that are ‘off’ are stationary. At the start of the simulation, only the ion is turned on, and is the only atom to have forces calculated and to be integrated. Some of the ‘off’ atoms will be used in the force calculations and will have forces assigned to them. If the resultant force exceeds a certain threshold, the atom is turned on and its motion is integrated. The simulation proceeds in this way with more and more atoms having their position integrated as energy becomes dispersed throughout the system.

We use two thresholds in our simulation; one for atoms interacting directly with the ion, and one for atom-atom interactions. We are, of course, mostly concerned with generating the correct motion for the ion, so the ion-atom interactions are of the most critical and require a lower threshold than the atom-atom interactions. In fact, for any reasonable threshold value, almost any ion-atom interaction will result in the atom being turned on, due to the large ion energy. Hence the ion-atom threshold is set to zero in these simulations, as adjusting the value gives no increase in efficiency.

In the case of the atom-atom threshold, we estimate a reasonable value by comparison to simulations without the moving atom approximation (MAA). Smith et al. found a force threshold of $`1.12\times 10^9`$ N for both atom-atom and ion-atom interactions gave the correct sputtering yield (when compared to simulations without the MAA) in the case of 1 keV Ar implant into Si. We have found a larger value ($`8.0\times 10^9`$ N) gives the correct dopant profile, when compared to a simulations without the approximation. Our ability to use a larger value is due to two reasons. The motion of atoms not directly interacting with the ion only has a secondary effect on its motion by influencing the position of directly interacting atoms, so small errors in the positions of these atoms has little consequence. Also, by dividing the interactions into two sets, we do not have to lower the threshold to give the correct ion-atom interactions.

### E Pair Potential Approximation and Recoil Interaction Approximation

While we use a many-body potential to describe a stable silicon lattice for low energy implants, this introduces a significant overhead to our simulations. For higher ion velocities, we do not need to use such a level of detail. A pair potential is sufficient to model the Si-Si interactions, as only the repulsive interaction is significant. Also, as the lattice is built at a metastable point with respect to a pair potential, with atoms initially frozen due to the MAA, and the section of material is only simulated for a short period of time, stability is not important. Hence, at a certain ion velocity we switch from the complete many-body potential to a pair potential approximation (PPA) for the Si-Si interactions. This is achieved in our code by setting the many-body parameter within the Tersoff potential to its value for undistorted tetrahedral Si, and results in a Morse potential splined to a screened coulomb potential.

We make a further approximation for still higher ion energies, where only the ion-Si interactions are significant in determining the ion path. For ion velocities above a set threshold we calculate only ion-Si interactions. This approximation, termed the recoil interaction approximation (RIA), brings the MD scheme close to many BCA implementations. The major difference that exists between the two approaches is that the ion path is obtained by integration, rather than by the calculation of asymptotes, and that multiple interactions are, by the nature of the method, handled in the correct manner.

We have determined thresholds of 90.0 eV/$`m_\text{u}`$ and 270.0 eV/$`m_\text{u}`$ for the PPA, and RIA, respectively are sufficiently high that both low and high energy calculated profiles are unaffected by their use. As the thresholds are based on the ion velocity, a single high energy ion simulation will switch between levels of approximation as the ion slows down and will produce the correct end of range behavior.

## IV Rare Event Algorithm

A typical dopant concentration profile in crystalline silicon, as illustrated in Fig. 3, has a characteristic shape consisting of a near-surface peak followed by an almost exponential decay over some distance into the material, with a distinct end of range distance. The concentration of dopant in the tail of the profile is several orders of magnitude less than that at the peak. Hence if we wish to calculate a statistically significant concentration at all depths of the profile we will have to run many ions that are stopped near the peak for every one ion that stops in the tail, and most of the computational effort will not enhance the accuracy of the profile we are generating.

In order to remove this redundancy from our calculations, we employ an ‘atom splitting’ scheme to increase the sampling in the deep component of the concentration profile. Every actual ion implanted is replaced by several virtual ions, each with an associated weighting. At certain *splitting depths* in the material, each ion is replaced by two ions, each with a weighting of half that prior to splitting. Each split ion trajectory is run separately, and the weighting of the ion is recorded along with its final depth. As the split ions see different environments (material is built in front of the ion, with random thermal displacements), the trajectories rapidly diverge from one another. Due to this scheme, we can maintain the same number of virtual ions at any depth, but their weights decrease with depth. Each ion could of course be split into more than two at each depth, with the inverse change in the weightings, but for simplicity and to keep the ion density as constant as possible we work with two.

To maximize the advantages of this scheme, we dynamically update the splitting depths. The correct distribution of splitting depths is obtained from an approximate profile for the dopant concentration. The initial profile is either read in (e.g. from SIMS data), or estimated from the ion type, energy and incident direction using a crude interpolation scheme based on known depths and concentrations for the peak and tail. Once the simulation is running, the profile and the splitting depths are re-evaluated at intervals. The algorithm to determine the splitting depths from a given profile is illustrated in Fig. 3. At the start of the simulation, we specify the number of orders of magnitude, $`M`$, of change in the concentration of moving ions over which we wish to reliably calculate the profile. We split ions at depths where the total number of ions (ignoring weighting) becomes half of the number of actual implanted ions. Hence we will use $`N`$ splitting depths, where $`N`$ is the largest integer $`M\times \mathrm{log}_210`$. The splitting depths, $`d_i`$ ($`1iN`$), are then chosen such that:

$$_0^{d_i}C(x)𝑑x=(1(\frac{1}{2})^i)\times _0^{\mathrm{}}C(x)𝑑x$$
(5)   
where $`C(x)`$ is the concentration of stopped ions (i.e., the dopant concentration) at depth $`x`$. Although we are using an approximate profile from few ions to generate the splitting depths, the integration is a smoothing operation and so gives good estimates of the splitting depths.

To minimize the storage requirements due to ion splitting, each real ion is run until it comes to rest, and the state of the domain is recorded at each splitting depth passed. The deepest split ion is then run, and further split ions are stored if it passes any splitting depths. This is repeated until all split ions have been run, then the next real ion is started. Hence the maximum we ever need to store is one domain per splitting depth (i.e., 16 domains when splitting over 5 orders of magnitude).

## V Simulation Details

All simulations were run with a Si{100} target at a temperature of 300 K. A surface amorphous layer of one, or three unit cells thickness was used. Dopant profiles were calculated for As, B, P, and Al ions; in each case it was assumed that only the most abundant isotope was present in the ion beam. The direction of the incident ion beam is specified by the angle of tilt, $`\theta ^{}`$, from normal and the azimuthal angle $`\varphi ^{}`$, as ($`\theta `$,$`\varphi `$). The incident direction of the ions was either (0,0), i.e. normal to the surface ($``$100$``$ channeling case), (7–10,0–30) (non-channeling), or (45,45) ($``$110$``$ channeling), and a beam divergence of 1.0 was always assumed. Simulations were run for 1,000 ions, with the splitting depths updated every 100 ions. A domain of 3$`\times `$3$`\times `$3 unit cells was used and the profile was calculated over either 3, or 5 orders of magnitude change in concentration. The simulations were run on pentium pro workstations running the Red Hat Linux operating system with the GNU g77 Fortran compiler, or SUN Ultra-sparc workstations with the SUN Solaris operating system and SUN Fortran compiler. The running code typically requires about 750K of memory.

## VI Results and Discussion

Two sets of results are presented; we first demonstrate the effectiveness and stability of the rare event enhancement scheme, and then give examples of data produced by the simulations and compare to SIMS data. Example timings from simulations are also given. A more extensive set of calculated profiles will be published separately .

### A Performance of the Rare Event Algorithm

An example of the evolution of splitting depths during a simulation is shown in Fig. 4, for the case of non-channeling 5 keV As implanted into Si{100}. The positions of the splitting depths near the peak stabilize quickly. Splitting depths near the tail take far longer to stabilize, as these depend on ions that channel the maximum distance into the material. Although atom splitting enhances the number of (virtual) ions that penetrate deep into the material, the occurrence of an ion that will split to yield ions at these depths is still a relatively rare event. The fact that all splitting depths do stabilize is also an indication that we have run enough ions to generate good statistics for the entire profile.

The paths of 5 keV As ions implanted at normal incidence into Si{100} are shown in Fig. 5, with the number of splittings shown by the line shading. The 1,000 implanted real ions were split to yield a total of 19,270 virtual ions. The paths taken by 27 split ions produced from the first real ion of this simulation, and the resulting distribution of the ion positions are shown in Fig. 6. The final ions are the result of between 3 and 6 splittings, depending upon the range of each trajectory. This is typical of the distribution of splittings for one real ion; the final depths of ions are not evenly distributed over the entire ion range, but are bunched around some point within this range. This reflects how the impact position of the ion and collisions during its passage through the amorphous layer affect its ability to drop into a channel once in the crystalline material. The weighting of the second 500 of the ions (after the splitting depths had stabilized) is plotted against final depth in Fig. 7 (note the log scale).

We have estimated the uncertainty in the calculated dopant profiles in order to judge the increase in efficiency obtained through the use of the rare event enhancement scheme. The uncertainty was estimated by dividing the final ion depths into 10 sets. A depth profile was calculated from each set using a histogram of 100 bins, with constant bin size. A reasonable measure of the uncertainty is the standard deviation of the distribution of the 10 concentrations for each bin. Fig. 8 shows calculated dopant profiles from 1,000 real ions for the case of 2 keV As at (7,0) into Si{100}, obtained with and without atom splitting over five orders of magnitude. The profiles are plotted with the uncertainty represented by the size of error bars; the length of each error bar corresponds to half the standard deviation of concentrations in that bin. The uncertainty is constant in the case of the profile obtained with the rare event scheme, whereas the profile obtained without the scheme is only reliable over one order of magnitude. 

Timings from these simulations, and a simulation with splitting to three orders of magnitude are given in Table I. From these timings, we can estimate the efficiency gain due to the rare event algorithm. We decrease the time required by a factor of 89 in the case of calculating a profile to three orders of magnitude, and by a factor of 886 when calculating a profile over 5 orders of magnitude, compared to the estimated time requirements without rare event enhancement. The gain in efficiency increases exponentially with the number of orders of magnitude in concentration over which we wish to calculate the profile.

### B Comparison of Profiles to SIMS Data

The remaining figures show the calculated concentration profile of B, As, and P ions for various incident energies and directions. Profiles were generated from a histogram of 100 bins, using adaptive bin sizes; the final ion depths were sorted and the same number of virtual ions assigned to each bin. No other processing, or smoothing of the profiles was done. Also shown are low dose ($``$ $`10^{13}`$ ions/cm<sup>2</sup>) SIMS data; for comparison, all profiles were scaled to an effective dose of $`10^{12}`$ ions/cm<sup>2</sup>. We have also examined Al ion implants, but were unable to match calculated profiles to the available SIMS data for a physically reasonable parameter value in our electronic stopping model. This may be due to one or more of the following reasons; Al is the only metal that we are implanting; the Al-Si interaction is the only interaction for which we do not have a pair specific ZBL potential; we only have a very limited set of SIMS data to compare to.

In the case of the low energy ($``$ 10 keV) implants, we compare to SIMS data obtained with a thin and well controlled surface layer; here we assume one unit cell thickness of surface disorder in our simulations. For the other cases considered here, the surface was less well characterized; we assume three unit cells of disorder at the surface, as this is typical of implanted Si. For the low energy implants, we have calculated profiles over a change of five orders of magnitude in concentration; for the higher energy implants we calculate profiles over 3 orders of magnitude. The results of the REED calculations show good agreement with the experimental data. In the case of the low energy implants, the SIMS profile is only resolved over two orders of magnitude in some cases, while we can calculate the profile over five orders of magnitude.

We give timing results from several simulations, as examples of the cpu requirements of our implementation of the model. Note, the results presented here are from a functional version of REED, but the code has yet to be fully optimized to take advantage of the small system sizes (around 200 atoms). Timing data are given in Table II, for profiles calculated over five orders of magnitude on a single pentium pro. Run times are dependent on the ion type and its incident direction, but are most strongly linked to the ion velocity. We estimate a runtime of approximately 30 hours per $`\sqrt{\text{keV}/m_\text{u}}`$, for this version of our code.

## VII Conclusions

In summary, we have developed a restricted molecular dynamics code to simulate the ion implant process and directly calculate ‘as implanted’ dopant profiles. This gives us the accuracy obtained by time integrating atom paths, whilst obtaining an efficiency far in excess of full MD simulation. There is very good agreement between the MD results and SIMS data for B, P, and As implants. We are unable to reproduce published SIMS data for Al implants with our current model. This discrepancy is currently being investigated; our findings will be published separately. We can calculate the dopant profile to concentrations one or two orders of magnitude below that measurable by SIMS for the channeling tail of low dose implants.

The scheme described here gives a viable alternative to the BCA approach. Although it is still more expensive computationally, it is sufficiently efficient to be used on modern desktop computer workstations. The method has two major advantages over the BCA approach: (i) Our MD model consists only of standard empirical potentials developed for bulk Si and for ion-solid interactions. The only fitting is in the electronic stopping model, and this involves *only one* parameter per ion species. This should be contrasted to the many parameters that have to be fit in BCA models. We believe that by using physically based models for all aspects of the program, with the minimum of fitting parameters, we obtain good transferability to the modeling of implants outside of our fitting set. (ii) The method does not break down at the low ion energies necessary for production of the next generation of computer technology; it gives the correct description of multiple, soft interactions that occur both in low energy implants, and high energy channeling.

We are currently working to fully optimize the code, in order to maximize its efficiency. The program is also being extended to include a model for ion induced damage, amorphous and polycrystalline targets, and to model cluster implants such as BF<sub>2</sub>. We also note that the scheme can be easily extended to include other ion species such as Ge, In and Sb, and substrates such as GaAs and SiC.

## VIII Acknowledgments

We gratefully acknowledge David Cai and Charles Snell for providing us with their insight during many discussions, and Al Tasch and co-workers for providing preprints of their work and SIMS data. This work was performed under the auspices of the United States Department of Energy.

## IX References

## X Figures

## XI Tables