Binary Collision Approximation Simulations Research Articles

Low energy ion-solid interactions: a quantitative experimental verification of binary collision approximation simulations

Ultra-low energy ion implantation has become an attractive method for doping of two-dimensional materials and ultra-thin films. The new dynamic Monte Carlo program IMINTDYN based on the binary collision approximation allows a reliable prediction of low energy implantation profiles and target compositional changes, as well as efficient simulation of high energy light ion scattering. To demonstrate the quality of these predictions and simulations, we present a model case experiment where we implanted W ions into tetrahedral amorphous carbon with low (10 keV) and ultra-low (20 eV) ion energies and analyzed the W implantation profiles with high resolution Rutherford backscattering spectrometry (HR-RBS). This experiment is compared with a complete simulation of all aspects of ion-solid-interactions of the experiment using the new IMINTDYN program. A unique novel simulation option, also relevant for implantation into 2D materials, is the inclusion of the vacancy as target species with dynamic vacancy generation and annihilation. Whereas simulations neglecting vacancy formation cannot reproduce the measured implantation profiles, we find excellent agreement between simulated and measured HR-RBS spectra. We also demonstrate the important role of simultaneous weak collisions in the binary collision approximation at low projectile energies.

High accuracy interatomic potential model for binary collision approximation and its application into sputtering yield estimation for amorphous carbon

Although binary collision approximation (BCA) and molecular dynamics (MD) are well used for plasma-material interaction simulation, the atomic collision in the energy above 10 eV is often out of scope in general potential model for MD. For BCA, the Ziegler–Biersack–Littmark (ZBL) potential has been often employed for high energy collision. In the present work, as one of modernizations of BCA, more accurate potential model for high energy collision, which is named ReGenerated ZBL (ReGZ) potential, was proposed. The function from of the ReGZ potential was analytically derived from the spherical electron density of an independent atom. To evaluate the potential, sputtering yield and reflection coefficient were compared between the BCA simulations using the ReGZ potential and the ZBL potential. In addition, by also improving the surface binding energy by density functional theory, the sputtering yield using the ReGZ potential becomes consistent with existing results.

Modeling experimental low energy ion scattering multi-angle maps with molecular dynamics FAN

Multi-angle maps obtained with low energy ion scattering (LEIS) are not only a convenient way to visualize the atomic-level surface structure of single crystals, but also provide a comprehensive data set that is useful for quantitative fitting with models. However, modeling these maps is often difficult due to the complexity of the scattering process, as well as the number of simulations needed to accurately reproduce a single experimental map. For example, traditional molecular dynamics (MD) and binary collision approximation (BCA) simulations typically require on the order of one billion simulated ion trajectories to obtain sufficient statistics. This problem is exacerbated when modeling impact-collision ion scattering spectroscopy (ICISS) maps, as backscattering cross sections are much smaller than forward-scattering cross sections. A workaround to this problem was given by Neihus and Spitzl in the form of FAN (Niehus and Spitzl, 1991), a simulation that approximates the angular dependence of the backscattering signal by only simulating ion trajectories that reach and emanate from lattice atoms, using a BCA code. This greatly reduces the number of trajectories that must be calculated, though it can oversimplify some of the more complex collision sequences. In this work, we develop a fully three-dimensional FAN code that uses MD to simulate trajectories rather than the BCA. Our approach exploits the improved speed of FAN, while retaining the more accurate ion trajectory calculations of MD. We apply the MD-FAN simulation to model LEIS measurements on two surfaces. First, we benchmark our MD-FAN simulation in a comparison to a BCA simulation for 3 keV He+ incident on W(111) that has been previously reported in Wong et al. (2020). The MD-FAN map had stronger agreement with the experimental ICISS measurements than did the BCA map, despite containing nearly 100-fold fewer ion trajectories. Second, to obtain experimental data sets for a more complex surface, we perform experimental ICISS measurements using a 3 keV He+ beam on a magnetite single crystal, Fe3O4(001) to generate a multi-angle map. The simulated MD-FAN map for the Fe3O4(001) surface once again agreed well with the experimental map, allowing for quantitative comparisons. The MD-FAN simulation also provided insight into the scattering processes, revealing that some of the observed pattern arises from 3 keV He+ that have backscattered into the detector from Fe atoms more than 2 nm below the surface.

Binary collision approximation simulations of ion solid interaction without the concept of surface binding energies

The sputtering yield of atoms removed from a solid upon ion irradiation can be calculated using Monte Carlo simulations based on the binary collision approximation. These simulations rely on the assumption of a planar surface potential barrier similar to the work-function for electron emission. Although this model yields good agreement with experimental data, the surface binding energy often has to be arbitrarily adjusted. We propose the bulk binding energy model to calculate sputter yields using sublimation enthalpies and a spherically symmetric potential barrier at the position where the sputtered atoms are generated in the bulk. A planar potential barrier at the surface no longer exists. We also consider the energy density deposited in nuclear collisions, which may lead to a lower sublimation enthalpy compared to standard room temperature conditions. The new bulk binding energy model is able to quantitatively explain many experimental sputtering data.

Quantitative investigation of surface structure and interatomic potential with impact-collision ion scattering spectroscopy

Helium ion beam interactions with materials have important implications for magnetic confinement fusion, material modification, and helium ion microscopy. These interactions depend on the precise physics of how helium ions channel into the materials, which can vary greatly based on the local crystalline orientation. In this work, we performed a dedicated experiment to investigate helium ion channeling in a well-characterized tungsten single crystal. Time-of-flight impact-collision ion scattering spectroscopy was used to obtain multi-angle maps of the backscattering intensity for 3 keV He+ → W(111). We found that the backscattering intensity profile arising from helium ion channeling could be well described by a shadow cone analysis. This analysis revealed that subsurface W atoms as deep as the ninth monolayer contributed to the backscattering intensity profile. Binary collision approximation simulations were performed with MARLOWE to model the experimental maps with sufficient accuracy to allow for quantitative comparisons using reliability factors. These quantitative comparisons were applied to investigate how the W lattice structure and He–W interatomic potential affect the multi-angle maps.

Ion-induced nanopatterning of silicon: Toward a predictive model

We review recent progress toward the development of predictive models of ion-induced pattern formation on room-temperature silicon, with a particular emphasis on efforts to eliminate fit parameters in the linear regime by means of experimental measurements or atomistic simulations. Analytical approaches considered include “mechanistic” models of the impact-induced collision cascade, the Crater Function Framework, and continuum treatments of ion-induced stress and viscous flow. Parameter evaluation methods include molecular dynamics and binary collision approximation simulations, as well as wafer curvature measurements and grazing incidence small-angle x-ray scattering. Mathematical detail is provided in the context of key results from pattern formation theory, which are also briefly summarized.

Atomic insight into concurrent He, D, and T sputtering and near-surface implantation of 3C-SiC crystallographic surfaces

Abstract We present results from atomistic simulations of sputtering and near-surface implantation of concurrent He, D, and T bombardment of cubic silicon carbide (3C-SiC). This is achieved by first establishing a many-body interatomic potential parameter set to treat interactions of He and hydrogenic species in 3C-SiC informed by ab-initio calculations. To obtain sputtering yields we perform both classical molecular dynamics and binary collision approximation simulations for normal incident particles having energies ranging from 25 to 800 eV. We find that due to differences in species surface binding energy of various crystallographic surfaces in 3C-SiC, the sputtering yield of Si is significantly less than that of C, but sputtering yields show limited sensitivity to crystallographic surface orientation. An exception to this occurs when the terminating crystallographic surface plane is more rich in Si rather than C, resulting in comparable sputtering yields of Si and C. The influence of temperature on sputtering is explored and shows limited effect. Finally, the nature of implanted He, D, and T within 3C-SiC surfaces is investigated to understand implantation profiles and stability of defects.

Molecular dynamics studies of ion beam implantation and patterning of silicon: Effect of noble gas cluster formation

The use of energetic ion beams to induce nanopattern formation at surfaces has been well studied both experimentally and theoretically. However, the influence on morphological evolution of the implanted species themselves remains little understood, particularly in the case when the incident ion species does not interact chemically with the target material. In this work, MD simulation results are presented for cumulative ion bombardment of Si to a fluence of $3\ifmmode\times\else\texttimes\fi{}1{0}^{15}\phantom{\rule{0.28em}{0ex}}\mathrm{c}{\mathrm{m}}^{\ensuremath{-}2}$ or more for a range of incident ion energies (20--1000 eV), angles ($0\ensuremath{-}{85}^{\ensuremath{\circ}}$), and species (Ne, Ar, Kr, Xe). For most cases, the implanted ions are observed to form gas clusters or bubbles beneath the surface as the fluence increases. The implantation and cluster formation decrease in magnitude with increasing ion incidence angle, and remain fairly similar for the heavier-than-Si species (Ar, Kr, and Xe). However, the implantation and cluster formation are much more prominent for Ne irradiation. As the fluence continues to increase beyond $\ensuremath{\sim}1{0}^{15}\phantom{\rule{0.28em}{0ex}}\mathrm{c}{\mathrm{m}}^{\ensuremath{-}2}$, the gas clusters begin to become exposed to the vacuum as the Si layers trapping the gas atoms are eroded by the incident ions. The exposed gas clusters then degas very rapidly, leading to disruption at the surface and to viscous material flow of Si into the void left behind. Comparison to dynamic binary collision approximation (BCA) simulations indicates that cluster formation and degassing contributes to a wide distribution of single-impact emission yields of implanted ions, contrary to intuitive expectations based on BCA simulations. Notably, the increased size and frequency of many-atom implanted ion emission events contributes to a much lower concentration of the implanted species than is otherwise expected from BCA simulations. Additionally, this cluster degassing phenomenon is conjectured to provide a potential ``antipatterning'' mechanism by disrupting or destroying nanopattern ``seeds'' at the surface. This could provide an additional mechanism to improve model predictions of critical angles for patterning transitions, and may also provide at least a partial explanation for the difficulty of obtaining patterns on Ne-bombarded Si surfaces.

Si amorphization by focused ion beam milling: Point defect model with dynamic BCA simulation and experimental validation

A Ga focused ion beam (FIB) is often used in transmission electron microscopy (TEM) analysis sample preparation. In case of a crystalline Si sample, an amorphous near-surface layer is formed by the FIB process. In order to optimize the FIB recipe by minimizing the amorphization, it is important to predict the amorphous layer thickness from simulation. Molecular Dynamics (MD) simulation has been used to describe the amorphization, however, it is limited by computational power for a realistic FIB process simulation. On the other hand, Binary Collision Approximation (BCA) simulation is able and has been used to simulate ion-solid interaction process at a realistic scale. In this study, a Point Defect Density approach is introduced to a dynamic BCA simulation, considering dynamic ion-solid interactions. We used this method to predict the c-Si amorphization caused by FIB milling on Si. To validate the method, dedicated TEM studies are performed. It shows that the amorphous layer thickness predicted by the numerical simulation is consistent with the experimental data. In summary, the thickness of the near-surface Si amorphization layer caused by FIB milling can be well predicted using the Point Defect Density approach within the dynamic BCA model.

Study of tritium desorption by energetic ion bombardment from tungsten material

To solve the problem of tritium retention in a nuclear fusion reactor, the removal of hydrogen isotopes retained in a tungsten material by exposure to the other hydrogen isotope plasma was investigated. We proposed the balance equation between the absorption and desorption of the hydrogen isotope atoms in the tungsten material. Absorption ratio and desorption yield were calculated by binary collision approximation simulation. Simulation results showed that desorption yield was proportional to the concentration of the retained hydrogen isotope. Therefore, the decrease in the retained hydrogen isotope atoms by hydrogen plasma irradiation can be regarded as approximately exponential decay. The decay rate in the case that retained deuterium atoms are removed by light hydrogen irradiation was two thirds of the decay rate in the case that retained tritium atoms are removed by deuterium irradiation.

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 1keV 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Å for amorphous Si, whereas at small displacements a difference between BCA and MD arises due to a material flow 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 difference due to displacements in amorphous pockets.

Atomistic modelling of the interface of Si nanocrystal structures in a-SiO2 before and after ion irradiation

As the basis of the conventional microelectronics industry, silicon is the most appropriate material for integration of electronic and optical functionalities in the same chips. Recent investigations in the field have lead to a drive to overcome the problem of poor light amplification properties of silicon. It has been shown that Si nanocrystals (Si-nc) embedded into a dielectric matrix emit light with a wavelength in the visible range. Ion implantation is widely used to form nanocrystals, and a subsequent step of inverse Oswald ripening can be used to control their size uniformity and flat space distribution. By means of molecular dynamics, atomistic models of Si-nc inserted into amorphous silica (a-SiO2) have been constructed. The defect structure of the interface region is discussed. The influence of ion implantation on the Si-nc atom distribution is analyzed and compared to binary collision approximation simulation results. The amorphization of Si-nc’s due to high-dose ion implantation is also discussed.

MD simulation of sputtering on surface index and surface roughness dependence

Low-energy sputtering is studied through molecular dynamics (MD) and binary-collision approximation (BCA) simulations. First, MD simulations as a preliminary study are performed for 5 keV Ar impacts on a Cu target. When random numbers gave the direction of the target crystal axis for each impact, the total sputtering yield almost agreed with the experimental yield of a polycrystalline Cu target. Second, we performed MD and BCA simulations for self-sputtering of W with an incident energy of 100 eV and incident angles of 0°–85°. We found that the many-body collision results in high-yield sputtering at grazing incidence. The incident angle dependence of sputtering yields fade out for an MD result with a rough surface.

Vacancy-related defect distributions in 11B-, 14N-, and 27Al-implanted 4H–SiC: Role of channeling

The defect distributions in 11B-, 14N-, and 27Al-implanted epitaxial 4H–SiC are studied using monoenergetic positron beams. At least three types of defects are needed to account for the Doppler broadening annihilation spectra and two of the defects are tentatively identified as VSi, and VSiVC. By comparing the defect profiles extracted from the annihilation spectra to the chemical profiles determined by secondary ion mass spectrometry, and to the primary defect profiles obtained from binary collision approximation simulations, it is concluded that the defects found at depths considerably deeper than the projected range of the implanted ions mainly originate from deeply channeled ions.

Atomic displacement cascade distributions in iron

Full molecular dynamics (MD) and its binary collision approximation (BCA) are used in a complementary way in order to study displacement cascade distributions in iron. Frenkel pair distributions are particularly narrow and symmetrical. They are thus well described by their first moment. Therefore, quantitative estimates by MD are possible. The comparison between the dependence of Frenkel pair production on primary energy predicted by both computational techniques suggests a post-cascade recombination model. Its physical grounds are discussed. The variability between spatial distributions of individual cascades is particularly large as a consequence of instability, which takes place in the early stage of the cascade development. The subsequent loss of correlation with initial conditions is statistically demonstrated on the basis of BCA simulations of 5000–15 000 cascades. Sufficient statistics can be reached by MD in order to characterise spatial distributions within cascades. It comes out of systematic comparison between MD and its BCA that, after the ballistic phase, the spatial extent of both vacancies and interstitials tends to increase. This phenomenon correlates with atomic mixing in the cascade core. This mixing is not predicted in the BCA. It is suggested to be responsible for the fragmentation of vacancy clusters formed during the ballistic phase of the cascades.

Prediction of the morphology of the as-implanted damage in silicon using a novel combination of BCA and MD simulations

In order to predict type and amount of defects created by keV ions under realistic implantation conditions, a combination of computer simulations based on the binary collision approximation (BCA) with classical molecular dynamics (MD) calculations is proposed. Time-ordered BCA simulations are applied to ballistic processes with characteristic energies above several 10 eV. Athermal, rapid thermal, and thermally activated processes with lower characteristic energies are treated by MD simulations. They yield the as-implanted defect state formed several 10 ps after ion impact. The MD calculations are performed in cells which are much smaller than the entire volume of the collision cascade of an incident ion but much larger than the distance between nearest-neighbor atoms in the lattice. The as-implanted damage produced by a single ion in a certain cell is found to be completely determined by the nuclear energy deposition of the ion into the cell. Therefore, the MD calculations need to be performed only in one cell for different values of nuclear energy deposition, and statistical considerations based on BCA simulations can be employed to obtain the depth profile and the total number of different defect species (vacancies, interstitials, disordered atoms, etc.) created on average per incident ion. The novel simulation method is applied to investigate the damage morphology produced by 15 keV B +, 30 keV P +, and 15 keV As + implants.

Improving the Understanding of Ion-Beam-Induced Defect Formation and Evolution by Atomistic Computer Simulations

AbstractThe morphology of the as-implanted damage in silicon is investigated using a recently developed combination of time-ordered computer simulations based on the binary collision approximation (BCA) with classical molecular dynamics (MD) calculations. The method is applied to determine the type and the amount of defects formed within the first nanosecond after ion impact. The depth profile and the total number of different defect species (vacancies, interstitials, disordered atoms, etc.) produced on average per incident ion are calculated for B+ (15 keV), P+ (5, 10, 20, 30 keV), and As+ (15 keV) implantations. It is shown that the as- implanted defect structure depends not only on the nuclear energy deposition per ion but also explicitly on the ion mass. Therefore for each ion species the damage morphology exhibits characteristic features. For heavy ions the percentage of extended defects is higher than for light ions. In all cases investigated the number of free or isolated interstitials exceeds the amount of free vacancies. The results obtained allow a microscopic interpretation of the phenomenological model for the as-implanted damage employed in conventional BCA simulations in order to describe the dose dependence of the shape of ion range profiles. They can be also applied to get more realistic initial conditions for the simulation of the defect kinetics during post-implantation annealing.

The correlation of defect distribution in collisional phase with measured cascade collapse probability

The spatial distributions of atomic displacement at the end of the collisional phase of cascade damage processes were calculated using the computer simulation code MARLOWE, which is based on the binary collision approximation (BCA). The densities of the atomic displacement were evaluated in high dense regions (HDRs) of cascades in several pure metals (Fe, Ni, Cu, Ag, Au, Mo and W). They were compared with the measured cascade collapse probabilities reported in the literature where TEM observations were carried out using thin metal foils irradiated by low-dose ions at room temperature. We found that there exists the minimum or “critical” values of the atomic displacement densities for the HDR to collapse into TEM-visible vacancy clusters. The critical densities are generally independent of the cascade energy in the same metal. Furthermore, the material dependence of the critical densities can be explained by the difference in the vacancy mobility at the melting temperature of target materials. This critical density calibration, which is extracted from the ion-irradiation experiments and the BCA simulations, is applied to estimation of cascade collapse probabilities in the metals irradiated by fusion neutrons.

The correlation of defect distribution in collisional phase with measured cascade collapse probability

The spatial distributions of atomic displacement at the end of the collisional phase of cascade damage processes were calculated using the computer simulation code MARLOWE, which is based on the binary collision approximation (BCA). The densities of the atomic displacement were evaluated in high dense regions (HDRs) of cascades in several pure metals (Fe, Ni, Cu, Ag, Au, Mo and W). They were compared with the measured cascade collapse probabilities reported in the literature where TEM observations were carried out using thin metal foils irradiated by low-dose ions at room temperature. We found that there exists the minimum or “critical” values of the atomic displacement densities for the HDR to collapse into TEM-visible vacancy clusters. The critical densities are generally independent of the cascade energy in the same metal. Furthermore, the material dependence of the critical densities can be explained by the difference in the vacancy mobility at the melting temperature of target materials. This critical density calibration, which is extracted from the ion-irradiation experiments and the BCA simulations, is applied to estimation of cascade collapse probabilities in the metals irradiated by fusion neutrons.

Slowing-down time of energetic atoms in solids

Computer simulation models for the slowing down of energetic atoms in matter are often based on the binary collision approximation (BCA). Such models typically ignore the temporal aspects of the problem. A method is described whereby calculations of the times at which energetic particles reach their collision points may be included in such BCA simulations without altering the ``event-driven'' structure of the model. The enhanced model is illustrated by calculations of the time required by 10-eV to 10-keV Cu atoms to slow down in both crystalline and amorphous Cu targets. In addition to the mean slowing-down time, distributions of slowing-down times and correlations between slowing-down times and particle ranges are presented. The utility of range-time correlations in studying penetration problems is demonstrated.