Crater formation by nanoparticle impact: contributions of gas, melt and plastic flow

The processes underlying crater formation by energetic nanoparticle impact are investigated using molecular dynamics simulations. Both metallic and van-der-Waals-bonded targets are studied. We find a transition from crater formation by melt flow at small impact energies to an evaporation (gas flow) mechanism at higher energies. The transition occurs gradually at impact energies per atom of a few tens of the cohesive energy of the target. van-der-Waals-bonded solids do not exhibit the melt flow cratering regime, in agreement with the narrow liquid zone in their phase diagram. We find that the size of the target region heated above the critical temperature roughly corresponds to the crater volume. The transition shows up most clearly in the increase of the volume of ejected material relative to the crater volume. Finally, we demonstrate the punching of dislocations below the crater for high-velocity impact into ductile targets, leading to a contribution of plastic flow to the crater volume.


Introduction
The impacts of energetic nanoparticles on surfaces may lead to erosion of the target and the formation of a crater. This process is important in the space environment, where dust particles impinge on the surfaces of atmosphereless planets, moons, asteroids or comets [1]. A special case of interest is metallic (spacecraft) surfaces, while relevant projectiles are dust grains (composed of rock and/or ice), but also metal and other debris in the low-Earth-orbit environment [2][3][4][5]. Dust grains may measure between 1 nm and 1 µm [6], while debris may have considerably larger sizes. Typical velocities are in the range of 3-100 km s −1 ; this is conventionally called the hypervelocity impact regime. Further interest in such impacts can be found in applications of cluster-surface modification [7][8][9][10][11].
The volume of the excavated crater scales in good approximation linearly with the total kinetic energy of the projectile. Such a linear dependence is found in experiments of µmand mm-sized metal dust particles or bullets on various materials [12][13][14]. Hydrocode simulations typically do not include an intrinsic length scale and also display this linear behavior. Molecular dynamics simulations of cluster impacts [15,16] have also found this linear scaling. In a recent publication [17], we could show that this linear scaling extends over a wide range of impact energies and cluster sizes. However, the question of which mechanisms lie beneath crater formation has not been satisfactorily answered until now. Undoubtedly, there is not a single mechanism responsible, but-depending on the projectile size, its energy and the nature of the target-several mechanisms may be involved.
Molecular dynamics simulations may prove useful in studying nanoparticle impacts since they allow us to extract detailed (atomistic and thermodynamic) information on the crater formation mechanisms for these small projectiles. Here, we use these simulations to demonstrate that for a metallic target both melt flow, at lower energies, and gas flow, at higher energies, contribute to crater formation. For van-der-Waals-bonded materials, which exhibit only a narrow melting regime, the contribution of melt flow could not be identified. For ductile materials, such as the Cu target studied here, plastic flow by dislocation motion can also contribute to crater growth.

The method
The molecular dynamics method employed is standard [18,19]. Results for two materials, Ar and Cu, are presented. The bombarding clusters always consist of the same material as the irradiated target, i.e. only the case of self-bombardment is studied. The clusters of a spherical shape are chosen and consist of N = 1000 atoms. Their diameters amount to 2.8 and 4.1 nm for Cu and Ar, respectively. For the Ar system, a Lennard-Jones potential [20,21], and for the Cu system a many-body potential of the embedded-atom type [22] has been chosen [23]. The latter potential describes equilibrium and non-equilibrium structures in good agreement with ab initio calculations [23]. In both cases, the potentials have been splined to an appropriate high-energy potential [24,25] in order to accurately model close collisions. The bonding properties of the material are represented by the cohesive energy U , which amounts to U = 3.54 eV in the case of Cu and 0.0815 eV for Ar. The total kinetic energy of the impacting cluster will be denoted by E. A scaled energy,  3 . We verified that these target sizes were sufficiently large by performing additional simulations with larger target sizes. In order to study the cluster-induced crater sizes and shapes, simulations have to be performed for sufficiently large times so that stable results are obtained; we employed simulation times up to t end = 120 ps. At the lateral and bottom sides of the simulation target, we employ damped boundary conditions in order to mimic energy dissipation to the surrounding target material. The Cu target consists of an fcc crystal with a (100) surface. In the case of Ar, we employ an amorphous target; test simulations showed no strong systematic differences between results for amorphous and crystalline targets.
At the end of the simulation time, we monitor the size of the bombardment-induced crater. We define the crater volume V (in dimensionless units) as the number of atoms missing below the original surface plane [15]. The dimensional volume can be recovered from it with the help of the atomic volume of fcc Cu, = 11.8 Å 3 (and analogously 36.1 Å 3 for Ar).
Relevant thermodynamic data of the two materials, Cu and Ar, are gathered in table 1. We also plot in figure 1 the phase diagrams of the two materials. The Ar results have been taken from the generic Lennard-Jones phase diagram as described in [26][27][28][29][30], while the Cu data are due to [31]. By normalizing temperature and density to their values at the critical point, the main difference between a van-der-Waals-bonded material and a metal becomes evident: the metal has a relatively small triple point temperature and hence a broad range of temperatures and densities where a liquid phase exists. We may discuss the width of the liquid region with the help of the ratio of the critical to the melting (or triple-point) temperature. This ratio is 2.1 for the van der Waals material and 9.0 for Cu. Thus, the widths of the liquid zone differ by a factor of 4.2 (on a temperature scale). Table 1. Thermodynamic data of Ar and Cu. U , cohesive energy; T m and T v , melting and vaporization temperatures at 10 5 Pa; and T c , p c and n c , critical temperature, pressure and density.

Melt and gas flow
In figure 2, we display the volume (measured as the number of missing atoms) of the crater excavated by a 1000-atom nanoprojectile of the same species as the target as a function of the scaled impact energy . It is seen that-beyond a threshold energy th -over the wide range of energies investigated the crater volume increases linearly with the impact energy, The coefficient of proportionality, a, is termed the cratering efficiency. From the simulation results it can be read off that a = 0.56 (0.52) for Cu (Ar). The fact that the cratering efficiency assumes such similar values for different materials is indicative of a common crater formation mechanism. The thresholds are th = 5600 and 4700 for Cu and Ar, respectively. We note that the concept of a 'threshold' to crater formation was introduced and discussed in a previous publication [19]. Dependence of the crater volume V (measured as the number of missing atoms) on the scaled impact energy per projectile atom, E/NU . Selfbombardment of clusters containing N = 1000 atoms on Ar and Cu has been simulated. Some of the data has been taken from [35]. The error bars denote the standard deviation of the simulation results for the cases when several impacts have been simulated.
Linear energy dependences above a threshold energy, similar to equation (2), have been observed earlier in the context of crater formation, plastic surface deformation and target atom displacements induced by cluster impact on metal and Si surfaces [15,17,19,36,37]. Such a linear dependence shows that a fixed fraction of the impact energy is used for crater excavation; it is expected to break down if part of the impact energy is deposited deep inside the target-this might be the case for small clusters impinging at high energy. It is a salient feature of this and previous investigations [17] that the energy linearity of the crater volume holds over a wide range of impact energies and cluster sizes.
At and below the threshold, the physics changes [19]. Evidently, for very small impact energies, the cluster will land on the surface without damaging it, V = 0. A rough estimate of the threshold is therefore given by E th = NU , where it is assumed that each cluster atom delivers the energy necessary for bond breaking, U , to the surface. Our simulation results for the N = 1000 clusters gave E th /(N U ) = 5.6 and 4.7 for Cu and Ar, respectively. At the energies between E th /(N U ) = 1 and ∼5, the surface becomes damaged, but not all the impact energy can be efficiently used for cratering and there is only a shallow crater. Only above E th /(NU ) ∼ 5 do our simulations show a linear scaling with energy.
We display in figure 3 snapshots of the early phase of crater formation in Ar. We apply a particular color code for displaying temperature in these data. We differentiate material below the melting temperature T m (solid-blue) and above the critical temperature (gaseous-red). The liquid regime is subdivided by introducing the boiling temperature T v (at 10 5 Pa): green-molten and yellow-'boiling'; this subdivision is physically irrelevant and only serves to indicate the steepness of temperature gradients within the liquid.
The crater formation process can be described as gasification: the material has been heated above the critical temperature of Ar; thus it is gaseous and free to flow out into the vacuum. We note that the same features are seen at smaller impact energies; however, melt flow is difficult to detect in the rather soft Ar target. The absence of a pronounced melt flow is understandable from the phase diagram of Ar, figure 1; the molten phase occupies only a tiny region in the phase diagram. Note also that even at t = 60 ps, where the crater size has stabilized, evaporation is still continuing. This is characteristic of a Lennard-Jones material with its high vapor pressure. This late evaporation stage should, however, not be confused with the gas-flow stage, where the supercritical material flows out of the forming crater. While the importance of gas flow in a weakly bound solid such as Ar appears plausible, it is interesting to inquire whether this phenomenon is also seen in Cu. We display the results of our simulations for two impact energies, /N = 6.1 (figure 4) and /N = 70 (figure 5). The  snapshots for the higher-energy impact were selected at (by a factor of 1/3) earlier times, in agreement with roughly three times higher impact velocity. At the lower impact energy (figure 4), we clearly recognize that melt flow is the decisive mechanism for crater excavation. Supercritical Cu exists in the crater only during a small time window (<1.5 ps); the emission of material during this early time is better described as sputtering than as gas flow in the sense in which it was used for Ar. Note, in particular, that the crater rims were always in the molten phase, but not supercritical. Quite in contrast, the snapshots for Cu at the higher impact energy ( figure 5) show signs of a strong gas flow. The region heated to supercritical temperatures is largest at times around t = 0.5 ps; at this time and also later, the gas flow of the gasified material can be directly observed. Here, part of the developing rim-i.e. that part of the material pushed up by the huge compressive pressures developing under the projectile-has been heated above T c and flows out into the vacuum. Note that the position of the final crater rim is farther away than the position of the early rim; while it can be assumed that part of the early rim will have survived at the surface and only be pushed to the sides, part of it will undoubtedly have been vaporized and ejected. Note that the crater volume in figure 5 appears to correspond roughly to the maximum size of the supercritical volume (at time t = 0.5 ps). Remarkably, this also appears to be true for Ar (figure 3) at around t = 4 ps and even in the 'melt-flow' case of Cu (figure 4) at t = 0.75 ps. This motivates us to calculate the maximum possible target volume, V crit , that can be heated to supercritical temperatures by the projectile energy. With table 1 we see that kT c and U are roughly proportional to each other, kT c = αU , where α ∼ = 0.15 for a Lennard-Jones material such as Ar and α = 0.16 for Cu. Taking 3kT c as the energy necessary per atom for heating to T c , we estimate and hence we estimate the cratering efficiency as a ∼ = 2. This is in good order-of-magnitude agreement with the simulated values of a ∼ = 0.5 ( figure 2). We attribute the difference between the simple estimate, equation (3), and the simulation result to three factors: (i) the target material is on average heated to temperatures above T c ; (ii) the lower part of the supercritical material has obtained momentum toward the target rather than toward the vacuum; (iii) part of the energized volume will recondense on the walls of the forming crater, while heat conduction cools it, rather than be emitted. We note that a previous gas-flow model of ion bombardment of rare-gas solids was able to describe sputter yields and sputtered particle energy distributions satisfactorily [38]. Tubes: dislocation lines enclosing stacking fault planes, detected using the DXA algorithm [41]. Tube color changes from blue to red with increasing dislocation length. Visualization has been prepared using Paraview [42].
A similar calculation of the volume that can be heated up to melting temperature gives a considerably larger value, a = 3.8 for Ar and a = 10.1 for Cu, which is evidently far away from the simulation result.
For gas-flow sputtering, models are available [38][39][40] and might as a first approximation also be used to estimate the crater volume. For the contribution of melt flow to crater excavation, no simple model appears to be available. We therefore present such a model in the appendix.

Plasticity
Besides melt flow, plastic flow can also occur in the bombarded target. We analyzed the dislocation structures generated by cluster impact using the Dislocation eXtraction Algorithm (DXA) of Stukowski and Albe [41]. Generally, we find that for lower velocities, below ∼5 km s −1 , clusters smaller than ∼10 000 atoms do not produce any dislocations. There is a production of point defects, and there might be some transient defective structures. For instance, for the particular event of 21.7 keV Cu 1000 impact (8.1 km s −1 ), the resulting pressure pulse is not strong and long enough to lead to true dislocation punching, since it falls below the homogeneous nucleation threshold for Cu of ∼30 GPa [43] after only 1 ps, and the material below the crater recovers to a nearly perfect crystal after ∼5 ps.
The situation is quite different at higher impact speeds. The particular case of 250 keV Cu 1000 impact (27.5 km s −1 ) is analyzed in figure 6. The immense dislocation density below the crater, ∼3 × 10 16 m −2 in figure 6(a), together with fast dislocation velocities driven by large pressure gradients, leads to a complex network of dislocations that will survive after the crater region has relaxed. Thus we can talk about a threshold for plasticity that is not yet reached at 8.1 km s −1 , but is surpassed at 27.5 km s −1 . The emerging plasticity can be associated with a hardness increase in the near-crater region, as compared to the pristine material before bombardment. Hardness values can be estimated from dislocation densities [17], using a Taylor hardening model, and lead to reasonable agreement with experimental values for macroscopic impactors [14].
We present in figure 6(b) a close-up view of the emerging plasticity showing developing shear loops and the enclosed stacking faults. Such a dislocation structure, with shear loops emitted near the crater surface in all available {111} planes, has many similarities to the structure produced by nanoindentation [44,45], and to void growth/collapse induced by tensile or compressive strain [46,47], including the presence of 'multi-planar' loops, as shown in figure 6(b).
Plastic flow can be quantified using Orowan's equation [43,48], relating plastic strain s p to dislocation density ρ and velocity v: where b is the magnitude of the dislocations' Burgers vector. Assuming dislocation densities such as those shown in figure 6(a) and typical dislocation velocities of around a tenth of the sound velocity, together with the crater formation times discussed above, one obtains plastic strains that imply a significant contribution of plastic flow to the total crater volume [17].

Correlation of cratering and sputtering
Finally, we investigate how the crater volume V correlates with the ejection (or sputter) yield Y . This yield is defined as the number of atoms which have been emitted from the target and have lost contact with any other target atoms. This correlation helps us to understand how much of the crater material has been ejected by the crater excavation into the vacuum; the remainder has been either deposited on top of the surface (as crater rim), preferably by melt flow, or pushed into the target (by plastic flow). The data of the ratio Y/V in figure 7 show that above the threshold energy of E/NU = 5.6, Ar exhibits a rather stable ratio of Y/V ∼ = 0.13. We take this as an indication that gas flow contributes an approximately constant fraction to the excavation of the crater. For Cu, however, immediately above threshold (E/N U = 5) a minimum in Y /V with values of only ∼ 0.06 is seen; this characterizes the melt-flow regime of figure 4. The values increase steadily until, at around E/N U = 20, values similar to Ar-and even larger-are obtained, indicating the onset of the gas-flow regime. These data indicate that for Cu the transition from melt-flow to gas-flow cratering occurs at around E/N U = 20.
In a recent publication, Delcorte and Garrison [49] reported on simulation results for the sputtering and crater formation of a polyethylene target bombarded by C 60 projectiles. They estimated the crater volume by adding the sputtered material to the material deposited on the crater rim above the surface. For this system they found that the sputter yield increases roughly linearly with energy (above a threshold at around 2 keV). The ratio Y/V shows a maximum of 0.4 at threshold (2 keV) and then slowly decreases to values of 0.25 at the highest energy simulated, 10 keV. It thus appears that the crater volumes obtained in the polymer target are relatively smaller than in our cases. However, the true craters in these simulations may actually be larger than estimated-polymers will be pushed out of the energized zone into the crater walls, thus widening and deepening the crater-and we hence consider the difference in results between the polymer target and ours as not too striking. We note that high-energy small-cluster impacts may lead to an even more pronounced gasflow contribution to crater formation. As an example, 16 keV impact of Au n clusters (4 n 12) on an Au target, i.e. E/N U = 350-1050, leads to values of Y /V = 0.52-0.61 [50]. These values are indicative of a strong gas-flow contribution to sputtering, and indeed, Colla et al [50] present further evidence in the form of temperature, pressure and flow distributions. We note that for even smaller cluster impacts (n < 4), Y/V was found to decrease when energy is increased. In these cases, the projectile buries part of its energy deep inside the target; while still a crater is excavated, only a little material is sputtered. We also note that for single-ion bombardment, the ratio is also energy dependent, and it is ∼0.1 for both keV [18] and MeV ions [51], with a significant amount of material flowing to the surface.

Conclusions
Due to the short time scales and to the extreme states of matter involved, performing experimental and theoretical investigations on the mechanisms underlying the crater formation caused by nanoparticle impact on surfaces is difficult. Using molecular dynamics simulations, we explore these processes. We find that at high impact energies, the early stages of crater formation are characterized by the flow of supercritical gas out of the impact zone, thus excavating the crater. For ductile targets, plastic flow contributes with dislocation nucleation and motion into the target. At smaller impact energies, melt flow around the impacting projectile toward the surface is the dominant mechanism. For van-der-Waals-bonded targets, we could not identify the melt-flow regime, in agreement with the narrow liquid zone in their phase diagram. The transition occurs gradually at scaled impact energies (per projectile atom) of E/NU ∼ = 20. The transition from melt to gas-flow cratering shows up in particular in an increase of ejected (sputtered) material.

Acknowledgments
EMB and CJR acknowledge support through a grant from the SeCTyP U N Cuyo and grant no. PICT2009-0092 from the Argentinean Agencia Nacional de Promoción Científica y Tecnológica.

Appendix. A piston model for melt flow
A model to estimate the contribution of melt flow to crater excavation is presented here. It uses the idea that melt is created below the cluster impact point; the downward pushing cluster acts as a piston to drive out the liquid below it; see schematics shown in figure A.1. An alternative model could consider Poiseuille flow, driven by the large pressure just below the crater surface. We note that this model can actually be used only for clusters larger than those used in the molecular dynamics simulations of this paper; as demonstrated in the main part, our clusters are too quickly destroyed to act as a piston. This will change for larger cluster sizes.
To simplify our calculation we can assume a cylindrical geometry. Using a spherical geometry changes numbers only slightly. Here we also assume that the crater diameter does not change during the liquid flow. This is justified by our molecular dynamics simulations showing a rapid evolution of crater diameter, D cr , with a slower evolution of crater depth h. A liquid layer is rapidly formed below the projectile due to the immense amount of kinetic energy transfer. We assume that the liquid is incompressible and therefore the liquid mass displaced is given by the volume of piston submerged in the liquid, V . This is a reasonable assumption for a molten metal under these conditions. The submerged volume is We can take as the initial condition a crater depth h ∼ (D cl /2) + h 1 . The depth h 1 of the molten region, from our large-scale molecular dynamics simulations, is ∼3 nm. For the case of a projectile with D cl = 30 nm, and a Cu target to specify a 0 , this gives h ∼ 18 nm, and 1/(α 2 − 1) ∼ 20. We can assume that the piston pushes until reaching the bottom of the liquid layer, i.e. h 1 = h, resulting in h 1 > h, and therefore melt flow will be leaving the crater cavity.
The total volume of liquid pushed out (V melt ) can be calculated as a function of h 1 and D cl , and it is given in figure A.2. Very small projectiles do not lead to melt flow. A relatively small projectile, for instance with D cl = 10 nm, will only lead to a thin liquid layer (∼1 nm) and the contribution of melt flow to crater volume will be negligible. For the values above and D cl = 30 nm, the additional volume of the crater produced by such melt flow amounts to ∼15% of the resulting crater volume.
Assuming the projectile has already slowed down significantly, to a roughly constant velocity of 1/10 of its initial velocity v during melt flow, the time to push the liquid out is For v = 5 km s −1 , this gives a time of ∼15 ps for melt flow.