Comparison of Molecular Dynamics and Binary Collision Approximation Simulations for Atom Displacement Analysis

Molecular dynamics (MD) and binary collision approximation (BCA) computer simulations are employed to study surface damage by single ion impacts. The predictions of BCA and MD simulations of displacement cascades in amorphous and crystalline silicon and BCC tungsten by 1 keV Ar + ion bombardment are compared. Single ion impacts are studied at angles of 50 ◦ , 60 ◦ and 80 ◦ from normal incidence. Four parameters for BCA simulations have been optimized to obtain the best agreement of the results with MD. For the conditions reported here, BCA agrees with MD simulation results at displacements larger than 5˚ A for amorphous Si, whereas at small displacements a diﬀerence between BCA and MD arises due to a material ﬂow component observed in MD simulations but absent from a regular BCA approach due to the algorithm limitations. MD and BCA simulation results for crystalline W are found to be in a good agreement even at small displacements, while in crystalline Si there is some diﬀerence due to displacements in amorphous pockets.


Introduction
Studying ion irradiation induced damage at surfaces is important for the construction of shielding materials for applications in which the the materials are subject to small particle bombardment [1,2].Ion irradiation at ∼keV energies causes surface and near-surface radiation damage due to atomic displacements, sputter erosion and stress generation, which may cause surface instability.Depending on incidence angle, irradiation by ions can roughen or smoothen the surface [3].Under certain conditions, ion bombardment is an effective approach for generation of self-organized, periodic nano-structures.[3,4,5,6,7] Both smoothening and nano-pattern formation by ion beam irradiation were recently concluded, by both theoretical [8] and experimental [9] lines of reasoning, to be a coherent effect of he impact-induced displacements of the atoms that are not sputtered away, but come to rest at new locations within or on the solid.The theoretical approach requires knowledge of the "crater function" -the average surface height-change profile induced by an ion impact [10,11].The most direct method of determining the crater function is molecular dynamics (MD) studies of individual ion impacts.
Craters induced by single ion impacts have been intensively studied over the last decades.[12,13,14,15] Previous studies have shown that small changes in the shape of craters can lead to significant changes in macroscopic pattern-forming behavior.[16] Hence well-converged, accurate simulations are an essential input to crater function theory for drawing conclusions about surface stability or large-scale pattern formation.
MD is limited in the size of systems that it can model.Because higher-energy impacts create larger collision cascades involving more atoms, a practical upper limit exists on the impact energy for which well-converged MD crater functions may currently be obtained.In our experience, with the materials we have studied, the limit is of order 1 keV -an energy at which the collision cascade strongly overlaps the surface.Yet it is important to understand whether the dominance of impact-induced mass redistribution over preferential sputtering persists to higher energies, where the collision cascade primarily interacts with atoms below the surface.This creates the motivation to employ faster atomistic simulation methods such as the binary collision approximation (BCA).The essential difference between these methods is how the atom motion is considered.In MD simulations the system of atoms evolves by solving numerically the equations of motion.This method simulates the full many-body dynamics in atomic system, with the accuracy limited only by the reliability of the interaction Hamiltonian employed [17,18,19].MD evolves a finitesized molecular configuration forward in time, in a step-by-step fashion.In conventional MD simulation methods [20,21,22,23,24] the movements of all atoms are calculated.Therefore, MD methods describe the interactions involved in ion implantation more realistically compared to other computer simulation methods, but require much larger computer capacity [23].Reasonable agreement with zero adjustable parameters has been found between MD-informed theory and experiment for large-scale pattern formation [25].
One of the widely accepted techniques employed for sputter erosion, which is important at higher energies, is the Binary Collision Approximation (BCA) and its modifications that include surface relaxation effects.[26].This computationally more efficient method approximates the full atomic dynamics of a material by a series of binary collisions, neglecting possible many body effects.For each collision, the classical scattering integral is solved [27] for a given impact parameter between the moving atom and a stationary atom in the target material.The impact parameter is chosen randomly within the radius of the circular area of interaction cross-section and is calculated based only on composition and atomic density of the target material, if the structure is amorphous or its order can be ignored, as in the SRIM code [28].Alternatively it is chosen by tracking the trajectory of the moving atom in the crystalline structure, as in the MARLOWE code [29,30].In both codes the scattering angles of colliding atoms and the energy transferred in a collision are deduced from the numerical solution of the scattering integral with the universal repulsive potential Ziegler-Biersack-Littmark (ZBL) [28].This operation requires far fewer calculation steps than solving the equations of motion for all the atoms at each time step.Moreover, in contrast to MD methods, the BCA simulations follow the motion of a single energetic atom at a time and its interaction with a single partner, which also makes this approach computationally more efficient.For ∼ 1 keV energies 50 000 BCA trials are about 10 6 times faster than a single ion trial in MD.
Unlike MD, the BCA approach becomes less accurate with decreasing kinetic energy in the collisions, where multi-body interactions can become significant [31].Because BCA is the most practical method at energies too high for widespread simulation with MD, it is important to explore the accuracy of BCA in energy regimes accessible to MD only with difficulty, in order to develop reliable ways of calibrating BCA for more widespread use.In this paper we compare the results of BCA and MD for single ion impacts in an energy regime accessible to both simulation methods: 1 keV Ar + on Si and W.
In both BCA and MD simulations, atom displacement statistics was collected for all the displacements that took place during a single ion impact event, and averaged over the total number of ions simulated.

Simulation methods
Single 1 keV Ar + ion irradiation impacts of amorphous and crystalline Si cells and a crystalline W cell were simulated with classical molecular dynamics code PARCAS [32,33] and compared to the Monte Carlo BCA code CASWIN [34], in which for every collision scattering integral is solved and distance to next colliding atom is chosen randomly taking into account the density corresponding to MD simulation cell.To describe interactions between atoms the universal repulsive ZBL potential is used and electronic stopping power is applied between collisions.
In CASWIN code the planar model of surface barrier was used [? ].In this model a particle was considered as sputtered if its energy in the direction perpendicular to the surface exceeded the surface binding energy (the surface barrier in this case) for the sputtering species.This surface barrier was deduced from the particle's kinetic energy in perpendicular to the surface direction, which was calculated as E cos 2 Θ, where Θ is the angle between the particle direction and the outward normal to the surface.If this component becomes less than zero while the total kinetic energy of the particle is still greater than the cut off energy assumed for the surface layers, the particle experienced the internal reflection (the direction of the particle was a mirror reflection with respect to the surface) back into the bulk.Otherwise, the particle sputtered after the refraction due to the energy loss on the surface barrier.If the total energy after the deduction was below the cut off energy in the surface layers the particle remained at the surface as an adatom.In the present simulations a zero surface barrier was assumed for the incoming ions.
The main principles and advantages of both simulation methods are discussed in section 1.
Statistics were defined by 50 ion impacts in MD simulations and 50,000 ion trials in BCA at the angles 50 • , 60 • and 80 • from normal incidence.In MD simulations, the initial position of ions was fixed at the distance 10 Å above the cell surface, while the x and y coordinates and azimuthal angles for impacts were chosen at random.The surface of all three targets was relaxed for 1 ps before the impact to minimize built-in stress in the structure.
A crystal silicon cell of the size 140 x 140 x 115 Å with 113569 atoms and an amorphous silicon cell of the size 400 x 400 x 100 Å with 877952 were used in the simulations.The interaction between Si atoms was described using the environment-dependent interatomic potential (EDIP) [35].The sizes of the cells were chosen to fully enclose the developing cascade within a single cell.The amorphous Si structure was optimized using the algorithm of Wooten, Winer, and Weaire (WWW) [36] and subsequently relaxed using the EDIP.This method gives reasonably optimized amorphous structure with realistic density, where there are almost no coordination defects in the initial structure [37].
For the purpose of comparing to a system that does not become amorphous during irradiation, we also used the BCC tungsten structure with a lattice unit a = 3.165 Å and comprising 54000 atoms.The size of simulation cell corresponded to 95 x 95 x 120 Å.The embedded atom model (EAM) potential in MD was applied to describe interactions between W atoms.This empirical model is widely used for describing metallic bonding.Pair-potentials were employed for atom -ion interactions.The universal repulsive ZBL potential was used at small interatomic distances for a realistic description of strong collisions.
In MD simulations an infinitely large surface was mimicked using an open surface in the incoming ion direction, fixed atoms within ∼ 20 Å from the bottom of the cell, and periodic boundary conditions in the two lateral directions.The duration of each MD simulation was 10 ps.No temperature control was used for simulation of an ion impact.After irradiation, a Berendsen thermostat was used to cool the system to 0 K over a duration of 1 ps.The cooling was applied to quench very small displacements caused by lattice thermal vibrations.The small change in displacement statistics that occurs during the cooling process is shown in Fig. 1.
The following BCA input parameters were adjusted to find the best agreement with MD results: displacement threshold energy, surface binding energy, bulk cut-off energy, and surface layer cut-off energy.In both BCA and classical MD simulations the interactions between atoms in the sample are described with an interatomic potential and electronic stopping is taken into account as a frictional force.For MD and BCA to be comparable, we used identical interatomic potential for high-energy (close distance) interactions and electronic stopping powers.

Results and Discussion
Here we present and compare the results obtained from the MD and the BCA simulations.A good agreement between MD and BCA for large distances is found by varying BCA simulation parameters as described in the following.
The displacement threshold energy is the energy that a target atom needs to leave its lattice site and Comparison between MD simulation results before and after system cooling to 0K.
form a stable interstitial [38,39].Its values given in the literature range from ∼9 to 35 eV for silicon [40,41,42,43,44,45,46] and for tungsten ∼40 -90 eV [47,48].A 13 eV displacement threshold energy was adopted in our simulations for amorphous Si, 20 eV for crystalline Si and 60 eV for BCC W as a result of this parameter optimization in our simulations.
Values of the binding energy in tungsten obtained from different calculations vary between 7.9 [49] and 10.09 eV/atom [50].The experimental cohesive energy is 8.9 eV/atom [51].In our BCA simulations, we assumed the surface binding energy to be 8.5 eV for W and 4.7 eV for Si.If the incident and scattered atoms have enough energy to overcome surface binding energy, these atoms are sputtered.
It is necessary to define a criterion when the ion has stopped in the target material.In BCA, an atom is counted as moving when its total kinetic energy is larger than the cut-off value [52].In our simulations we chose the bulk cut-off energy and the cut-off energy in the surface layer to be 3 and 1 eV, respectively, for all three target materials.Here surface layer is defined as the layer of atoms that have the highest z coordinate.
The input parameters used in BCA that were obtained from the optimization described above are summarized in Table 1.Bulk cut-off energy (eV) 3 3 3 The average penetration depth of Ar + ions was checked for both simulation methods.The fraction of reflected Ar + ions, and and average depth of not reflected ion and standard deviation of the mean value of depth are presented in following tables (Table 2, Table 3, Table 4).The error on the fraction of reflected ions is calculated for MD results, whereas BCA results are presented with accuracy 0.001.
The ion depth varied from ∼ 17 − 26 Å for 50 • and 60 • , while for the angle 80 • Ar + , the ion was reflected from the surface in most of the single ion impact simulations and the depth for not reflected ions was ∼ 14 Å for both methods in amorphous Si.The difference in penetration depth between MD and BCA in crystal targets arises due to the ordered structure of material.The average displacement profiles from single Ar + ion impact on amorphous Si for MD and BCA are shown in Fig. 2.
At large displacements -above about 5 Å -the MD and BCA displacement distributions are in a good agreement, whereas at small displacements -below about one bond length -the BCA and MD results differ by over an order of magnitude.This difference is attributed to the neglect in BCA of collective many-body effects, such as such as heat spikes and flow of matter.No combination of BCA parameters can reproduce this.
Displacement field analysis would be necessary to understand the nature of displacements of both simulation methods in detail, but this is not necessary for the purpose of the current paper.
The displacement profiles for the amorphous Si target irradiated at angles 50    Similarly to amorphous silicon, the crystal silicon and crystal tungsten BCC structure were studied with MD and BCA methods.The atom displacement profile for 1 keV Ar + impacts on crystalline Si is shown in Fig. 5, and on tungsten is shown in Fig. 6.In the case of W, the curves corresponding to the displacements in the MD cell and the BCA cell agree over the whole range, all the way down to the smallest displacements at which the two can be compared.This result is in marked contrast to that described above for amorphous Si.We interpret the different behavior to the small displacement field associated with material flow, which appears to be a key feature of the irradiation of amorphous materials.Unlike in materials that are crystalline, or remain crystalline during irradiation, the memory of the original lattice site is lost during flow, and a large number of small displacements occur that would have been exactly zero if the potential well at the lattice site could have been maintained during the collision cascade.
Displacements smaller than the nearest neighbor distance of 2.74 Å in W (in Fig. 6) correspond to self interstitial atom (SIA) formation in 111 crowdion configuration [53] and strained atoms around defects.
The results for crystalline Si are intermediate between those of W and amorphous Si.It is known from experiments and simulations that individual ion impacts partially amorphize the crystalline structure In our simulations the structure is reset to fully crystalline in between individual impact trials, so although there is no amorphization effect of cumulative ion impact in the simulations, the partial amorphization effect of individual ion impacts may be an important determining factor in the results.
To investigate whether this partial amorphization is indeed significant, we analyzed some of the postimpact structure by evaluating the angular structure-factor [57,58].In this analysis, a list of all angles between nearest-neighbour bonds is formed, and compared after sorting with angles in the ground state crystal structure [57].This approach has been shown to be able to distinguish very well atoms in a crystalline or disordered (liquid or amorphous) environment [58,59].A disordered region consisting of about 70 atoms for the 50 • angle, 80 atoms for the 60 • angle and 40 atoms for the 80 • angle were recognized in the crystalline Si simulation cells, which means that due to the bombardment amorphous pockets formed in the crystal.
Thus the development of local amorphous regions upon ion impact of c-Si allows for the small displacements seen in the MD but not the BCA in Fig. 5, and thereby explains the intermediate nature of the results in c-Si.
In light of these results, it appears possible that pattern formation on materials that remain crystalline during irradiation may behave in a fundamentally different manner than on materials that are amorphous or become amorphous during irradiation.In the former case, the effects of mass redistribution may be suppressed.It may turn out that BCA is adequate for determining crater functions in crystalline materials over a wider range of conditions than in amorphous or amorphizable materials.An important direction for future research is to study the actual crater functions, in order to determine the contributions of the very small displacements to the parameters governing surface stability or pattern formation.

Conclusions
The comparison between BCA and MD simulations of collision cascades in amorphous silicon by individual 1 keV Ar + ion impact at three different incidence angles reveals a significant shortcoming of the BCA approach when applied to materials that are amorphous or become amorphous during ion irradiation.
The BCA approach misses the very large number of very small displacements seen in MD and attributed to the collective phenomenon of amorphous material flow.This shortcoming could seriously compromise our ability to use BCA to evaluate crater functions of sufficient accuracy for use in crater-function theory of topographical pattern formation or surface stability.This material flow component is absent in the simulated impacts on tungsten, where the BCA and MD displacement distributions agree quite well.The discrepancy between BCA and MD results is reduced but still significant for impacts on crystalline Si, which becomes partly amorphous during irradiation.These findings may have significant implications for the contributions of crater functions to surface stability or topographic pattern formation.
Finally, we note that it is curious to find that Monte Carlo BCA simulations, which do not account for the crystal structure of a material in any way, describe the atom displacements clearly better for crystalline than amorphous materials.

Figure 1 :
Figure 1: Atom displacement statistics for 1 keV single Ar + ion impacts on amorphous Si at 60 • from normal incidence.

Figure 2 :
Figure 2: Atom displacement statistics for single 1 keV Ar + impact on amorphous Si at 50 • incidence angle -MD and BCA

Figure 5 :
Figure 5: Atom displacement statistics for single 1 keV Ar + impact on crystalline Si at 50 • incidence angle MD and BCA

Figure 6 :
Figure 6: Atom displacement statistics for single 1 keV Ar + impact on crystalline W at 50 • incidence angle.

Table 1 :
BCA input parameters after the optimization procedure.

Table 2 :
MD and BCA comparison.Reflected fraction and mean penetration depth of Ar + ion in amorphous Si target.Penetration depth of Ar + 13.54 ±0.07 14.4 ±1.94The sputtering yield obtained from MD and BCA are presented in Table5.The lowest sputtering yield was obtained at incidence angles closer to the grazing angles -80 • from normal incidence for both simulation methods due to Ar + ion reflection from the target surface.No sputtering occurred in crystal W [60,61] incidence angle in the MD case and only 1.36 atoms per incoming ion were sputtered in a-Si.Much larger sputtering yield was obtained at the angle 60 • and angle 50 • atoms per incoming ion for amorphous Si.The angles 50 • and 60 • are reported to be critical for ripple formation[60,61].

Table 3 :
MD and BCA comparison.Reflected fraction and mean penetration depth of Ar + ion in crystalline Si target.

Table 4 :
MD and BCA comparison.Reflected fraction and mean penetration depth of Ar + ion in crystalline W target.