# Direct simulation of ion beam induced stressing and amorphization of silicon

## I Introduction

Mechanical response of materials to ion irradiation has implications for many materials applications. Ion processing of silicon is of widespread fundamental and technological importance due to its central role in the semiconductor industry. The high dose, large atomic number wafer implants used today lead to a significant amount of amorphization, which must be removed by subsequent annealing. Ion implant induced stresses can lead to such problems as substrate bending, delamination and cracking, and anomalous diffusion of dopants. The creation and modification of stress in silicon by ion irradiation has been examined experimentally by several research groups. Experimental measurements have shown that Si expands by approximately 1.8% upon ion beam induced amorphization. Little atomic scale modeling of the amorphization process has been conducted, although many attempts have been made to simulate the structure of amorphous Si (a-Si). These simulations usually involve a rapid quench of liquid Si and suffer from drawbacks such as limited system size and short simulation times. The vast majority of these simulations result in a structure with similar properties to experimental a-Si, except for the fact that the density is typically around 4% greater than crystalline Si (c-Si). Here, we study the formation and structure of a-Si by direct MD simulation of the ion beam amorphization process, rather than by a simulation of quenching.

## II Molecular Dynamics Simulation Scheme

We employ a minimal atomistic model of the radiation damage process and investigate stressing and deformation of the substrate. A full MD simulation incorporating a large enough section of the target material to completely contain the paths of many implanted ions, and to run for a realistic time is neither feasible, nor particularly desirable. What we wish to study is the response of a section of the material to the damage induced by the irradiation. In order to achieve this, we simulate a model system, that incorporates all the necessary features of the complete system.

At the ion energies used, the material response is a bulk, rather than a surface effect, so it is not necessary to explicitly include the surface in simulations. Hence, periodic boundary conditions are applied in all directions to minimize finite-size effects. Energy dissipation and equilibration of pressure are controlled by coupling the system to an external bath. This involves two parameters, $`\tau _T`$ and $`\tau _P`$, that are time constants controlling the rate at which temperature and pressure tend towards set values. This non-physical coupling gives us control over quenching and expansion of the system, but does limit the inferences we can draw about dynamic processes. We have analyzed the dependence of simulation results on the strength of these couplings in order to ensure that we do not introduce any systematic errors into our simulations. We assume that the ion dose is sufficiently low that the presence of ions in the material has a negligible effect and may be ignored. The effect of the irradiation is modeled by choosing primary knock-on atoms (PKAs) within the simulated region and assigning velocities to them from a distribution that can be dependent upon ion type, ion energy, and the depth in the substrate. PKAs are produced by collisions with the implanted ion, but the majority are due to collisions of knock-on atoms in sub-cascades. Hence our PKA velocity distribution must take account of all knock-on atoms that are created by ion implantation. In order to account for the free surface, the material is allowed to expand or contract in the direction normal to the surface; expansion is not allowed in the in-plane directions. The passage of high energy ions through the simulated region are rare events, between which the material will relax as energy is dissipated, then enter an essentially constant state for a long period of time. As this constant state contributes nothing to the response of the material, we omit it and start another PKA at this time. By taking this approach, we can reduce the system size to the extent that a representative set of simulations can be conducted on a few workstations in a matter of days.

## III Simulation Details

MD based ion implant simulations were used to to determine realistic PKA velocity distributions for various ion types and energies, at differing substrate depths; the resulting distributions were similar. The magnitude of PKA velocities could be approximately described by a $`\chi ^2`$ distribution. For an incident ion beam perpendicular to the surface, the in-plane distribution of PKA velocity vectors was uniform. The angle between the PKA velocity vector and the ion beam could be approximated by a Gaussian distribution with mean and standard deviation of $`85^{}`$ and $`10^{}`$ respectively.

PKA energies were chosen such that the size of displacement cascades was less than the system dimensions, to ensure that parts of the same damage cascade could not interact through the periodic boundaries. The selection of coupling constants for pressure and temperature control is critical; too long coupling times lead to excessive cpu requirements, whilst too short coupling times lead to unphysical system behavior, e.g., rapid quenching of hot disordered material will produce an excessive number of defects. Therefore, preliminary simulations were conducted with various combinations of PKA energy, system size and temperature and pressure coupling constants to investigate and minimize the sensitivity of the simulation on these parameters.

The final simulations were conducted using PKA energies with a mean of 1 keV. All silicon atoms were coupled to a 300 K heatbath with a coupling constant, $`\tau _T`$, of 1.0 ps; this was determined to be large enough that increasing it had no effect on simulation results. A new PKA was initialized as soon as the system temperature was quenched to below 320 K. Stress in the direction normal to the surface was relaxed to zero with a time constant, $`\tau _P`$, of 0.1, 1.0, or 10.0 ps. No single value of this parameter could be used as it effectively represents the ability of the material to locally deform, both elastically and through plastic flow. As such it is affected by the size of the ion path, ion mass and energy, the ion dose rate and the total ion dose. We therefore ran three sets of simulations spanning two orders of magnitude range in this parameter. This range is expected to contain most experimentally accessible conditions. In order to relate stress to change in atomic positions, a Young’s modulus of 113.0 GPa was used for both crystalline and damaged Si; the value of this parameter is expected to vary during amorphization, but its exact value is not critical as variations will only affect the value of the the non-critical time constant, $`\tau _P`$. The initial system was crystalline Si, with dimensions of approximately 50.0 Å $`\times `$ 50.0 Å $`\times `$ 50.0 Å, containing 6240 atoms, at zero pressure. In order to be able to simulate a sufficiently large sample, we use a many-body empirical potential for silicon. To our knowledge, the derivation of the pressure virial for this potential has not been previously published. A method of obtaining the virial, which may be applied to any many-body potential, is given in Appendix A.

## IV Amorphization of silicon during ion irradiation

Three sets of simulations are presented, the only difference between them being the magnitude of the pressure control coupling constant. In each case 9 separate simulations were conducted for each parameter set and the final results averaged. Data was recorded at the end of each quench, i.e., immediately prior to initializing a PKA, to generate radial distribution functions (RDFs), and distributions of bond angles, bond lengths, and atomic coordination. The simulation parameters chosen resulted in a mean time between PKA initialization of approximately 6.25 ps.

During irradiation, damage accumulates within the sample, resulting in the development of stress. As the simulation progresses, the material expands, or contracts to relieve the normal stress, as shown in Figs. 1 to 3. Anisotropic deformation occurs by flow of material to relieve the in-plane stress. The normal stress is relaxed to zero in the case of the two shorter $`\tau _P`$ values and remains negative for the longer $`\tau _P`$. The stress would be expected to reach an equilibrium state for any value of $`\tau _P`$, given a sufficiently long simulation time, but cpu constraints prevented us from achieving this for the largest value of $`\tau _P`$. The in-plane stress is intially positive (compression) due to damage accumulation. Once a sufficient density of damage is present the sample is able to relax via flow of material and the stress becomes negative (tension). The in-plane negative stress is due to the outflow of material during the ‘thermal spike’ phase of radiation damage. As the sample cools during expansion a tensile state is reached, but material is not able to flow back to relieve the stress. A similar response to ion irradiation has been previously observed in wafer curvature experiments. At long times, the damage appears to saturate and the material maintains a dynamic equilibrium, with a constant in-plane stress of around $`0.3`$ GPa. The percentage expansion of the system is dependent upon the value of $`\tau _P`$; smaller values (faster relaxation) result in expansion of the material, whilst the largest value gives an initial expansion followed by contraction. Due to cpu time constraints only the intermediate case was followed until the system reached an equilibrium condition. In this case the material reaches a density around 1.9% less than c-Si, approximately the same as the experimental value of 1.8%.

The distributions shown in Figs. 4 to 7 were calculated for the final state of each amorphization simulation, i.e., they correspond to the final data points in Figs. 1 to 3. The distribution of bond angles is shown in Fig. 4. To calculate the distribution, angles were weighted by the product of the Tersoff cut-off functions of the two bonds involved. The small peak at $`60^{}`$ corresponds to three-fold rings and five-fold coordinated atoms, the shoulder at around $`80^{}`$ is due to four fold rings, and the main peak is due to five and six fold rings. The distribution due to c-Si is also shown to illustrate the effect of damage and amorphization; the frequency is scaled by $`1/3`$ in order to allow the sharp peak to appear on the same plot as the a-Si results. The mean bond angle in the amorphized samples is around $`2^{}`$ less than the tetrahedral angle, with a standard deviation of approximately $`19^{}`$.

The distribution of bond lengths is shown in Fig. 5. When interpreting this it is useful to recall that the Tersoff switching function acts between 2.7 and 3.0 Å. The small peak between 2.9 and 3.0 Å is therefore due to non-bonding (repulsive) interactions. Bond lengths become about 0.09 Å greater than those in c-Si. The angle and length distributions due to the different relaxation rates are so similar that they cannot be easily distinguished. In order to count the coordination of atoms, we have to specify a criterion to decide if any pair of atoms are bonded. Based on the bond-length distribution we count all atoms with interaction distances within the mid point of the switching function (2.85 Å) to be bonded. The resulting coordination distribution is shown in Fig. 6. Very few ($``$3%) undercoordinated atoms are produced during amorphization, but approximately 26% become five-fold coordinated, and occasional ($``$2%) six-fold coordination occurs. Even though the average coordination of atoms increases to approximately 4.27 during amorphization, expansion is possible as bond-lengths also increase.

Table I contains structural information on the samples produced by these simulations. The mean and standard deviation of the structural parameters are given, and the percentage expansion of the irradiated samples and potential energy per atom are given relative to c-Si at 300 K. The Tersoff potential gives a cohesive energy for Si of 4.63 eV/atom at zero K and 4.59 eV/atom at 300 K (this is in agreement with equi-partition of kinetic and potential energy, as 300 K corresponds to a kinetic energy of 0.04 eV/atom). In the final damaged structures, the potential energy per atom is approximately 0.41 eV higher than that of c-Si at 300 K. This stored energy is approximately 8.9% of c-Si cohesive energy and is primarily due to bond angle distortions of around $`19^{}`$ as shown in Fig. 4 and Table I.

The resulting radial distribution functions, shown in Fig. 7 are indistinguishable between the various simulation sets, but differ from the experimental data obtained by neutron diffraction experiments. The broader peaks in the simulated RDFs indicate that the structures contain a higher degree of disorder than the experimental samples. There is also a feature at around 3 Å, where the artifact in the shape of the first minimum reflects the abrupt cutoff in the empirical potential. The first two peaks of the a-Si RDF correspond to those for c-Si showing the existence of short-range order and a tetrahedral like environment for each atom. There is little correlation between the positions of the other peaks, except for the third peak in the experimental data of Kugler et al.; this may be an indication of the presence of micro-crystallites in their sample.

## V Annealing

The damage simulations were carried out at room temperature with extremely limited times for structural relaxation. Therefore, the amorphized samples are expected to contain a large number of high energy defects; this is demonstrated by the difference between the experimental and simulated RDFs. In order to directly compare experimental and simulated samples, we have therefore carried out annealing simulations. The final structures from each set of damage simulations were taken, and subjected to four separate anneals. Anneals were done at room temperature, between room temperature and the a-Si to c-Si transition temperature, between the transition temperature and the c-Si melting point, and at the melting point. The Tersoff potential over-predicts the c-Si melting temperature, giving a value of approximately 3000 K. Therefore we applied the following scaling between actual and simulated temperatures:

$`T_{Tersoff}=\{\begin{array}{cc}T_{Experiment}\hfill & T_{Experiment}300\mathrm{K}\hfill \\ & \\ 1.947T_{Experiment}283.994\mathrm{K}\hfill & T_{Experiment}>300\mathrm{K}\text{.}\hfill \end{array}`$ (1)   
This resulted in annealing temperatures of 300 K, 1121 K, 2194 K, and 3000 K, which approximately correspond to experimental temperatures of 300 K, 722 K, 1273 K, and 1687 K respectively. Anneals were continued until no further change in any of the measured structural parameters was observed. All simulation parameters were kept the same as those used for the damage simulations apart from the anneal temperature. After annealing, the samples were quenched to 300 K over a time of 40 ps before the final structural properties were measured.

The structural changes during annealing are shown in the following plots for the $`\tau _P`$=1 ps damage simulation. The high temperature anneals resulted in large density changes as shown in Fig. 8. After cooling to room temperature to remove the effect of thermal expansion, the final densities of the annealed samples corresponded to 2.1% and 1.9% expansion, and 0.4% and 1.0% contraction relative to c-Si respectively for the 300 K, 1121 K, 2194 K, and 3000 K anneals. The room temperature anneal gives very little structural modification in the time simulated; all distributions are indistinguishable from the unannealed sample.

The 1121 K anneal results in modifications to bond angle, bond length, and coordination distributions and to the RDF, whilst the density is almost unchanged. This indicates that local defects are being removed without any global structural reorganization. This is in agreement with experimental studies which indicate that defects in a-Si are similar to defects in c-Si, i.e., they behave as interstitials and vacancies within the amorphous network. Annealing therefore occurs through diffusion and annihilation of defects, rather than by a reorganization of the network. Fig. 9 shows the final distribution of bond angles; the mean increases by approximately $`1^{}`$, to around $`1^{}`$ less than the c-Si angle, while the standard deviation of bond angles is reduced to $`16^{}`$. Bond lengths are reduced by an average of 0.03 Å, to 0.06 Å longer than those in c-Si, as shown in Fig. 10. The anneal results in a reduction in 3-fold and 4-fold rings as evidenced by the reduction in the frequency of bond angles of $`60^{}`$ and $`80^{}`$, and the reduction in 5-fold coordinated atoms. After annealing and quenching, approximately 1% of atoms are 3-fold coordinated, 83% of atoms are 4-fold coordinated, and 16% are 5-fold coordinated, as shown in Fig. 11. The calculated radial distribution function, shown in Fig. 12, is in very good agreement with the experimental data, with the only discrepancy due to the short range of the potential cutoff. Annealing leads to an energy gain of approximately 0.13 eV per atom, to give an energy in the final relaxed structure 0.28 eV higher than that of c-Si. This stored energy is approximately 6.1% of c-Si cohesive energy and is primarily due to bond angle distortions of around $`16^{}`$ as shown in Fig. 9.

The 2194 K anneal causes partial melting of the amorphous material and subsequent reconstruction into a higher density structure. The final quenched structure contains more defects than the 1121 K annealed sample, but this may be in part due to the rapid quench from a semi-liquid state. We observe no indication of any recrystallization in the time simulated. The 3000 K anneal clearly involves complete melting of the sample, as the density increases to above that of crystalline Si. Again, the annealed sample contains more defects than the 1121 K annealed sample, but less than the un-annealed sample.

## VI Conclusions

We have investigated the structural response of crystalline silicon to ion-irradiation. We are able to simulate stress generation and structural relaxation with a relatively simple model. Simulation of continual radiation followed by annealing generates amorphous silicon with a low level of defects. Amorphous silicon prepared by various experimental methods usually has a density between 2% and 10% less than that of crystalline material, whereas atomistic simulations usually produce samples with a density around 4% higher than c-Si. The fact that a-Si prepared by simulations of quenching liquid Si has a higher density than c-Si is not surprising. Amorphous silicon has a structure similar to that of the liquid, which is 5% denser than the crystal form at the melting point . Therefore a defect free continuous random network (CRN) structure for a-Si would be expected to result in a density greater than c-Si.

We conclude that by conducting a direct simulation of radiation induced amorphization, we have produced the experimentally observed structure of ion beam amorphized Si. The structure differs from the proposed defect free CRN models, and differs from a-Si structures formed by simulations of quenching. Although this metastable structure is at a higher energy than CRN a-Si, it cannot be annealed to produce that structure, as transformation to c-Si, or melting will occur at a lower temperature. This leads us to question what is meant by ‘amorphous’ Si, as an experimental sample termed amorphous may contain micro-crystals, CRN structures, and defects such as vacancies, interstitials, c-a boundaries, etc. In some sense the never experimentally observed high density structure can be regarded as true amorphous material, whereas other less dense materials must be regarded as containing defects at the very least, with a detailed structure that is preparation dependent.

## VII Acknowledgments

We thank M. Nastasi for illuminating discussions during initial parts of this work. This work was performed under the auspices of the United States Department of Energy.

## A Derivation of the Virial for Many-Body Potentials

### 1 Introduction

While the pressure virial is easy to calculate for a thermodynamic system described by pair potentials, the extension to systems modeled by many body potentials is not obvious. Here we derive the virial for Tersoff type potentials, but note that the method is readily modified to other many body, or three body potentials so long as the configuration energy can be written down as a function of atomic coordinates.

Following Smith we write the pressure with an explicit dependence on volume. This is achieved by relating absolute atomic coordinates, $`𝐫`$, to scaled coordinates, $`𝝆`$ through the volume, $`V`$, so that under isotropic expansion, or contraction of the system, the scaled coordinates of atoms are unchanged:

$$𝐫=V^{1/3}𝝆.$$
(A1)   
The pressure, $`P`$, is then given by:

$$P=NkT/V\mathrm{\Phi }/3V,$$
(A2)   
where the virial, $`\mathrm{\Phi }`$, is:

$$\mathrm{\Phi }=3VU([V^{1/3}𝝆]^N)/V.$$
(A3)   
### 2 Application to the Tersoff Potential

The configuration energy, $`U`$, is given by:

$$U=\underset{i=1}{\overset{N}{}}\underset{j<i}{}u_{ij}(𝐫^N),$$
(A4)   
where $`u_{ij}`$ is the energy of the bond between atoms $`i`$ and $`j`$, as a function of $`𝐫^N`$, the set of atomic coordinates. The functional form is:

$$u_{ij}=f(r_{ij})[V_R(r_{ij})b_{ij}V_A(r_{ij})],$$
(A5)   
where $`f(r_{ij})`$ is a cutoff function that makes the interactions short ranged, $`V_R(r_{ij})`$ and $`V_A(r_{ij})`$ are pair terms, $`r_{ij}`$ is the bond length, and $`b_{ij}`$ is a many body term. As $`b_{ij}`$ is a function of the positions of neighboring atoms of $`i`$ and $`j`$, these atoms experience a force due to the $`i`$-$`j`$ bond. Let the set of atoms involved in the calculation of $`u_{ij}`$, i.e. atoms $`i`$, $`j`$ and neighbors, be denoted by $`S_{ij}`$. Substituting Eq. (A4) into Eq. (A3) gives:

$$\mathrm{\Phi }=3V\underset{i=1}{\overset{N}{}}\underset{j<i}{}u_{ij}([V^{1/3}𝝆]^N)/V.$$
(A6)   
Since from Eq. (A1):

$$𝐫=\frac{1}{3}V^{2/3}𝝆V,$$
(A7)   
this can be rewritten as:

$$\mathrm{\Phi }=3V\underset{i=1}{\overset{N}{}}\underset{j<i}{}\frac{u_{ij}(𝐫^N)}{𝐫^N}\frac{V^{1/3}𝝆^N}{3V}.$$
(A8)   
Explicitly writing the dependence on interacting atoms gives:

$$\mathrm{\Phi }=3V\underset{i=1}{\overset{N}{}}\underset{j<i}{}\underset{kS_{ij}}{}\frac{u_{ij}(𝐫^N)}{𝐫_k}\frac{V^{1/3}𝝆_k}{3V},$$
(A9)   
or, in terms of atomic positions:

$$\mathrm{\Phi }=\underset{i=1}{\overset{N}{}}\underset{j<i}{}\underset{kS_{ij}}{}𝐟_k^{ij}𝐫_k,$$
(A10)   
where $`𝐟_k^{ij}`$ is the force on atom $`k`$ due to the bond between atoms $`i`$ and $`j`$. The atomic positions, $`𝐫_k`$, are local to the bond being considered, e.g., periodic boundary conditions are accounted for before the potential and force calculations. All local atomic coordinates can in fact be translated to be relative to the position of atom $`i`$ for each bond, so the force on atom $`i`$ can be neglected:

$$\mathrm{\Phi }=\underset{i=1}{\overset{N}{}}\underset{j<i}{}\underset{kS_{ij}^{}}{}𝐟_k^{ij}𝐫_{ik}.$$
(A11)   
If the interaction potential is pairwise, i.e., $`b_{ij}`$ depends only on atoms $`i`$ and $`j`$, Eq. (A11) reduces to the pair virial:

$$\mathrm{\Phi }=\underset{i=1}{\overset{N}{}}\underset{j<i}{}𝐟_j^{ij}𝐫_{ij}.$$
(A12)   
Individual elements of the pressure tensor, $`𝚽`$, can be obtained by replacing Eq. (A1) with:

$$𝐫=𝐋𝝆,$$
(A13)   
where $`𝐋`$ is a $`3\times 3`$ matrix with determinant $`V`$. For example, by using:

$$𝐋=\left(\begin{array}{ccc}L_{xx}& 0& 0\\ 0& 1& 0\\ 0& 0& 1\end{array}\right),$$
(A14)   
one obtains:

$$\mathrm{\Phi }_{xx}=3\underset{i=1}{\overset{N}{}}\underset{j<i}{}\underset{kS_{ij}^{}}{}f_{x_k}^{ij}r_{x_{ik}}.$$
(A15)