Rare event molecular dynamics simulations of plasma induced surface ablation

Sigmund's theory of high energy ion bombardment of surfaces¹ has been successfully used to understand surface penetration of ions and sputtering in many systems. The theory is based on Boltzmann transport and the binary collision approximation. The principle physics in this theory is that after the collision of an ion with the surface, surface atoms undergo several random collisions. Eventually, one of the surface atoms gets sufficient energy to break its bonds, and it ablates. This picture can be accurate at high energy, but the theory generally fails near threshold energies for ablation. Another approach is based on the generalized Langevin equation (GLE),^2,3 originally introduced for surface interactions by Doll and Adelman⁴ and later applied to three-dimensional systems by Tully.^5,6 Within the GLE approach, interactions between primary atoms in the region of interest are treated explicitly and the secondary atoms in the environment are governed by friction and stochastic forces. Here, we adopt a more traditional molecular dynamics (MD) approach, which is frequently used for surface sputtering calculations.^7–15 MD simulations can provide detailed atomic mechanisms of ablation for a specific material, dependence of ablation yield on energy, and the angle of incident ions that give rise to ablation events. This information is very useful in plasma etching for semiconductor technologies, where the interest is in removing surface atoms by colliding atoms or ions to the surface. The sputtering rate in plasma etching is relatively high and theoretical studies focus on high energy mono-energetic ion beams; temperature effects in plasma-surface interactions have not been as well-studied.

In this work, we focus on ablation rates from surfaces exposed to lower energy plasmas produced from a capillary arc discharge. Plasma from these devices have a density on the order of 10²⁶/m³, are highly collisional and rapidly thermalize with kinetic energies on the order of a few electron volts.¹⁶ Capillary discharge is used in space propulsion applications such as in a pulsed plasma thruster. A problem in such devices is surface erosion due to collisions of the plasma particles with the material making up the walls of the thruster. Accordingly, there is an interest in developing a fundamental understanding of plasma-surface interactions in these low-temperature plasmas. At low temperature, however, low ablation rates present a computational challenge. Direct MD simulation of collision events for a 1 eV thermal plasma, for example, results in ablation yields below 10⁻⁶ so that a more than a million trajectories are required to eject just a single atom from the surface. In this paper, we have implemented a computational method which allows for sampling these rare events and calculating low ablation rates.

Specifically, we consider a plasma composed of neutral Ar atoms. Although a plasma consists of many ionized particles, the degree of ionization is expected to be less than 1%.¹⁶ We assume that our plasma is in local equilibrium, i.e., has a Maxwell-Boltzmann distribution of particle velocities. As model materials, we consider Si and W as representative covalent and metallic materials. We distinguish between the surface and plasma temperatures because the plasma has orders of magnitude higher temperature than the surface. We use two approaches for modeling ablation. The first method uses direct MD simulations to simulation ablation events for plasma temperatures above 5 eV. At lower plasma temperatures, where ablation events are rare, we use a second method in which we sample ablation events as a function of plasma particle energy and then integrate over a Boltzmann distribution to give the thermal ablation yield. In this latter approach, the mechanism and threshold energy of ablation are shown to give predictive ablation yields from material surfaces down to low plasma energies, where ablation is an extremely rare event.

All calculations were performed using classical MD simulations using the LAMMPS code.¹⁹ In each simulation, an Ar atom was placed above a material surface and given an appropriate thermal velocity for a plasma particle in thermal equilibrium. Si was chosen as a model covalent material and W as a model metal. Two Si surfaces (100) and (111) were chosen to investigate the importance of surface structure for the probability of ablation. The inter-atomic potential between the atoms in the material were taken to be of the Stillinger-Weber (SW) form for Si²⁰ and the modified embedded atom method (MEAM) form for W.²¹ Ar was assumed to interact with surface atoms according to a Lennard-Jones potential with bond length σ, well depth ε, and cut-off values r_cut given in Table I. While the details of the Ar-surface interaction may not be quantitatively accurate, collisions at high energy (tens to hundreds of eV) are largely governed by the transfer of momentum which are insensitive to the details of the attractive portion of the potential. Numerical evidence of this is provided later.

We used periodic boundary conditions in the plane of the material surface and non-periodic in the direction perpendicular to the surface. The bottom layer of the slab was frozen. The top five to seven layers were simulated using constant energy dynamics. The middle layers, which constitute most of the simulation cell, were modeled with constant temperature dynamics using the Nosé-Hoover thermostat.^22,23 Collisions of Ar atoms heat the system, so the thermostat is essential for modeling dissipation of heat to the bulk of the material. Limiting the thermostat to the bulk prevents artificial interactions between the incident and surface atoms. The middle layer is similar to secondary zone atoms in the GLE approach, but different in that we explicitly calculate the motion of these atoms with constant-temperature MD.^5,6 In all simulations, the temperature of the thermostat was set to 1000 K. The time-step used to integrate the equations of motion was 0.2 fs. This small time-step was needed to ensure the conservation of energy for the high energy collisions.

Plasma interactions with three different surfaces were studied. The first, Si(100) surface model contained 1280 atoms in 20 layers of Si in the (100) direction, giving a total slab thickness of 25.85 Å; the lateral cell dimensions were a = b = 34.74 Å. After optimization of the atomic positions, the Si atoms dimerize at the surface, as shown in Fig. 1, as found in experiment. The dimers, however, do not have the buckling that is experimentally observed.²⁴ This is one of the shortcomings of the SW potential. The surface energy was calculated as 1.44 J/m². The second surface is Si(111). We do not consider the complex 7×7 reconstruction;²⁵ instead only the simpler 1×1 surface was studied.²⁶ The surface structure of (111) is simpler than (100); there are no dimers and the surface is nearly flat (Fig. 1(b)). The lateral cell dimensions are a = 26.90 Å and b = 30.74 Å. The slab is composed of 26 layers with 1664 atoms total. The third surface is W(100) (Fig. 1(c)). The lateral cell dimensions are a = b = 18.99 Å; the 720 atom slab is 35.59 Å thick.

To simulate surface collisions, an Ar atom was positioned randomly in a plane 5–6 Å above the surface. A velocity vector v = (v_x, v_y, v_z) was drawn from a thermal distribution at the plasma temperature to initialize a MD trajectory. In equilibrium, plasma particle velocities follow the Maxwell Boltzmann (MB) distribution

The incident flux of particles to the surface is

where ρ is the number density of the plasma. The factor of 1/2 accounts for the fact that the particles must be moving towards the surface to collide with it. These velocities towards the surface (with v_z >0) follow the Maxwell flux (MF) distribution,

The flux at which atoms are ablated from the surface can be calculated from a MD simulation as

where N_A is number of ablated atoms in time Δt from surface area A. The ablation flux can also be expressed as an integral over incident particle velocities

where Y_v, the velocity-resolved ablation yield, is the probability that a surface atom is ejected given an incident plasma particle velocity v.

The thermal ablation yield Y_T is the ensemble average of Y_v over incident particle velocities at the plasma temperature T,

Thus, the ablation flux is simply

the incident flux of atoms F_I times the thermal ablation yield Y_T.

As mentioned in the Introduction, we perform our calculations of Y_T using two different approaches. In the direct approach, for a given temperature, we chose an initial flux-weighted velocity of the Ar atom towards the surface. We simulated MD for 2 ps, during which either a surface atom was ejected or the energy of the impact dissipated into the bulk. Ten to 15 000 such trajectories were repeated over a range of temperature to calculate the fraction of trajectories that resulted in ablation events. This direct approach works well at high temperature but becomes ineffective at low plasma temperature where the ablation yield approaches zero. For low temperature plasma conditions, we use the integration method to target ablation events near the threshold energy.

The method described here is not limited to the thermal plasma conditions. If additional effects such as a plasma drift velocity or sheath effects at the material surface can be quantified, a non-thermal distribution of incident particle velocities can be sampled using this same approach.

In the direct approach for modeling ablation, one can consider two flux regimes. In the high flux limit, the time between of collisions is shorter than the time scale for local surface relaxation, and surface damage from sequential collisions will accumulate. In the low flux limit, the time between collisions is sufficient for the surface to relax back to a crystalline structure. An estimate of the surface relaxation time and the flux value which separates these two limits can be made with long MD simulations. Figure 2 shows the potential energy of the surface during the first 100 ps in a 1 ns trajectory following the collision of a 100 eV Ar atom into the Si(100) surface. The system was minimized at 10 ps intervals giving potential energy values that are a direct measure of the surface damage. A surface area of roughly 10 Ų was damaged in the collision. Over 100 ps the surface anneals back to a defect-free structure. Based upon these conditions we can express the high-flux condition as

where A_D is surface area damaged in a single collision and τ is surface response time. This critical flux will depend upon the surface temperature and structure, but it provides an order-of-magnitude estimate for when the collisions between plasma particles and the surface can be considered independent. If an experiment is performed at atmospheric pressure for a 1 eV thermal plasma, the flux of Ar atoms is 2×10²⁶/m² s. Under these conditions, the low flux limit is the appropriate regime.

In our second approach for modeling low temperature plasmas and low ablation yields, we study collisions due to incident atoms with a specified kinetic energy. By characterizing the ablation yield as a function of energy, Y_E, we can integrate with the proper Boltzmann weight factor to obtain the average thermal ablation yield, Y_T, at the plasma temperature of interest. In probability theory, this method is known as law of total expectation.²⁷ According to this law, the average or expected value of Y_T for the distribution function f_MF(v) is Y_E averaged over the appropriate distribution function for E, where Y_E is the expected value of the ablation rate conditional that the energy of the incident atom is E. To do this, we first need to know the velocity for a given incident energy, f_B. Following Eq. (6), the constant energy average ablation yield in polar coordinates is

where f_B is the (yet unknown) velocity distribution and Y_B is the ablation yield from a beam of atoms with kinetic energy E and velocity v with angles θ$and$ϕ to the surface. Integrating Y_E with the appropriate flux distribution in polar coordinates gives

where

Equation (11) gives the velocity distribution for incident kinetic energy E. The distribution function f_B does not depend on the polar angle ϕ, which means that v_x and v_y should be chosen isotropically. The cos θ term in the integrand shows that f_B depends linearly on v_z, since v_z = vcos θ. This defines the distribution function for v_z as

where v_z ∈ [0, |v|]. Note here that v² = 2E/m is fixed. Summarizing this method: for a given energy E we choose v_z from the linear distribution in Eq. (12), and then the parallel components of the velocity according to

where ϕ ∈ [0, 2π] is a uniform random number.

Collisions at each energy are repeated 10 to 15 000 times to give the ablation yield Y_E. As discussed previously, the low flux approximation is more appropriate and we simulate every collision on a damage-free surface. We then use Eq. (10) to compute the thermal ablation yield Y_T for the range of temperatures of interest. Figure 3 shows how the Boltzmann distribution overlaps with Y_E at k_BT = 7 eV; the overlap region of the two functions determines Y_T.

The direct simulation is particularly relevant for high temperature plasmas for which the ablation yield is high enough to quantified in the thousands of MD simulations that can be run at each temperature. Particularly at high flux, one could worry about accumulated disorder of the surface increasing the ablation probability. Figure 4(a) shows that the ablation probability does indeed change as a function of the collision number, when simulated at 2 ps intervals without allowing additional time for the surface to anneal. At high thermal energies of 20 eV, there is actually an increased probability of ablation for the pristine surface. After some oscillation, steady state is reached after 1000 collisions, with a surface structure as shown in Fig. 4(b). At lower thermal energies of 5 eV, the ablation rate increases with damage and steady state is reached in a similar number of collisions. At 3 eV, no ablation events are observed, which places an upper bound of 10⁻⁴ on the ablation yield; we consider this below the threshold energy for ablation. Above 20 eV the surface of Si becomes disordered, as in amorphous Si; similar structures were reported in Ref. 28.

Importantly, we repeated the ablation calculations in the low flux limit where the surface was reset to pristine between collision events. Figure 5 shows the difference in ablation yield between the low and high flux cases. For k_BT = 20 eV the difference is only 12%. Similarly, the ablation yield is weakly sensitive to changes in the surface temperature, with threshold energies changing by less than 1 eV for a surface at 0 K as compared to 1000 K. This means that surface disorder does not play a critical role for the plasmas of interest in this study with thermal energies k_BT < 20 eV.

At low plasma energies, the direct method of simulating ablation events becomes increasingly difficult as the thermal yield Y_T approaches zero. The integration method focusses instead upon sampling the ablation yield as a function of energy Y_E. A subsequent reweighing with the Boltzmann distribution gives Y_T at arbitrarily low temperatures.

Figure 6 shows Y_E sampled for Ar on the Si(100) surface. At high energies, above 100 eV for Si(100), there is linear dependence of Y_E on energy. For lower energies, below 20 eV, there is we see also linear dependence with a smaller slope, as shown in Fig. 6(inset). At 17.3 eV the ablation yield becomes zero. This is the threshold energy, E₀. By numerically integrating Eq. (10) we obtain the average ablation thermal yields in Fig. 7. The ablation yield at low temperature is linear so that the dependence on plasma temperature is of the Arrhenius form,

Note that threshold energy E₀ = 17.3 eV is almost the same as the activation energy, E_A = 17.2 eV, as calculated from the slope of the Si(100) plot in Fig. 7. The reason the ablation kinetics follow an Arrhenius law can be understood from the linear dependence of Y_E on the incident energy near threshold. From Fig. 6(inset) we see that energies in the neighborhood of 17 eV, the ablation yield is

where α ≈ 10⁻⁴/eV. Using Eq. (10), we obtain the thermal ablation yield

For low temperatures, with k_BT on the order of an eV, the first term in the prefactor dominates,

This shows that E₀ ≈ E_A (from Eq. (14)) and the prefactor ln Y₀ = ln (αE₀) = −6.3 is close to what we obtain from the direct fit in Fig. 7. Interestingly, the Si(111) surface has a much higher activation energy of 40.6 eV, yet the prefactor is the same as for Si(100). The activation energy and prefactors for W are given in Table II.

Ablation yields for all three surfaces are shown in Fig. 7. Si(111) is the most stable surface and has an activation energy 2.3 times greater than Si(100). Even though the cohesive energy of Si is 4.3 eV and W is 8.8 eV (for our chosen potentials) Si(111) is more stable than W(100). Thus, the cohesive energy alone cannot explain the ablation trends from our MD simulations. The binding energy of a surface atom (the surface vacancy formation energy) is better correlated with the stability of the surface. The binding energies for Si(100), Si(111), and W(100) are 6.1, 10.6, and 9.5 eV, respectively. While the relative order of these surface atom binding energies are the same as the activation energy for ablation, they are not strongly correlated.

A second important consideration is that activation energies for ablation are several times larger than the surface bond energies and cohesive energies of the materials. If there were a disparity between the mass of the plasma particle and surface atom masses, energy transfer efficiency could play an important factor. In a direct elastic collation between Ar and Si, however, 97% of the incident Ar kinetic energy is transferred to a Si surface atom, calculated as 4m_Sim_Ar/(m_Si + m_Ar)². Therefore, the mass ratio between Ar and Si cannot explain this difference between the bond energy and activation energy.

Individual MD trajectories near the threshold energy can give more clues about the energy loss. Figure 8 shows the path of an Ar collision on a Si(100) surface for a 20.7 eV trajectory. An Ar atom collides with one Si atom (A) and gives it horizontal momentum. The Si atom collides with its neighbor Si atom (B) from the same dimer. The second Si atom is ejected, but before ejection it collides with a third Si atom (C) from the neighboring dimer. Over ten trajectories near threshold were examined and found to have the same ablation mechanism, which gives us confidence that this is the most likely ablation mechanism. Even for incident energies of 400 eV, where we see ablation of two Si atoms, the mechanism is the same. In such cases, the Ar atom pushes two Si atoms to the side, then these Si atoms collide with their neighbors. During each of these collisions, the incident kinetic energy is dissipated to the surface. Figure 9 shows the energy transfer between atoms in the lowest energy ablation process of Fig. 8. In the first collision 67% of the energy is transferred to the primary Si atom (A), then (A) transfers about half of its energy to the neighboring Si atom (B).

To understand the dissipation process, consider a diatomic molecule. If initially the first atom has the kinetic energy E and the second atom has energy zero, after the collision, on average half of the energy goes to vibrations and half of the energy goes into translation. Vibrations are responsible for breaking bonds. Hence, we have to give to the first atom twice the bond energy to break the bond. The actual energy transfer in our case is more complicated, but the mechanism of energy dissipation is the same. Overall, the kinetic energy of the Ar atom has to be 3-4 times the surface bond energy to eject an atom from the Si(100) surface.

For Si(111) the situation is similar. One of the surface Si atoms gets momentum from Ar and collides with its neighbor; that second atom is ejected. On Si(111), however, there are no dimers and the surface is less corrugated than Si(100). In order for ablation to occur, a Si atom must obtain upward momentum, and for that, the Ar atom must dive deeper into the surface. These considerations make the Si(111) surface more resistant to ablation.

The ablation mechanism for the W(100) surface is qualitatively different. Near threshold, an Ar atom penetrates the surface, and bounces off a subsurface W atom. On its way back to the surface it collides with a surface W atom, which is ejected without colliding with any other atoms. This mechanism was previously reported in Ref. 13 for other metal surfaces, and it is possible that this ablation mechanism is general for most metals. For semiconductors and insulators complex surface reconstructions can make the ablation mechanism more complicated.

The angular distribution of incident Ar atoms that cause ablation provides additional information about the ablation mechanism. At incident kinetic energies near threshold, Figs. 10(a) and 10(b), show ablation only for Si atoms in Si(100) when the incident Ar atom approaches near 45°. The 45° incidence angle is required to transfer as much horizontal momentum as possible to surface Si atom. At higher energies, Figs. 10(c) and 10(d) we see a wider distribution and more ablation at horizontal and vertical angles, since the concentration of energy along the ablation reaction coordinate is less essential.

Earlier experimental²⁹ and theoretical¹³ work have also shown that the activation (threshold) energy is several times larger than the binding energy. Reference 29 suggested that activation energy is 4H, where H is the sublimation heat of sublimation for materials with a high energy transfer factor (4m_pm_s/(m_p + m_s)²), where m_p and m_s are the masses of the plasma and surface atoms. However, our simulations show that the correlation between the binding energy and the activation energy is not so simple, and can be better understood from the mechanism of ablation.

We have implemented a sampling method that makes it possible to simulate low temperature plasma surface interaction. Although our focus is plasma in thermal equilibrium, the approach can be easily modified to include non-thermal velocity distributions. Drift velocities, for example, can be important for simulating plasma-surface interactions in fusion devices. Most importantly, we observed that ablation rates follow an Arrhenius law at low temperature. The activation energy is found to correspond to the threshold energy of ablation and the prefactor is largely material-independent. The main shortcoming of our integration method is that it assumes a steady-state surface structure, which we took to be ideal. As atoms ablate from the surface, the surface will become rougher or some complex patterns will form.^14,30,31 While understanding effects of surface roughness goes beyond the scope of the paper, it will be interesting to understand how the mechanism and energy of ablation depends upon the surface structure.

This work was supported by the (U.S.) Air Force Office of Scientific Research (USAFOSR) (FA9550-11-1-0062), the (U.S.) Department of Energy (DOE) (Program No. DE-FOA-0000611), and the Welch Foundation (F-1841). We gratefully acknowledge Laxminarayan Raja, Francis Stefani, Roger Bengtson, and Sam Chill for many helpful discussions. Computational resources were provided by the Texas Advanced Computing Center and the National Energy Research Scientific Computing Center.