Molecular Dynamics Simulations of Nanoparticle-Surface Collisions in Crystalline Silicon

We present a microscopic description for the impacting process of silicon nanospheres onto a silicon substrate. In spite of the relatively low energy regime considered (up to 1 eV/atom), the impacting process exhibits a rich behavior: A rigid Hertzian model is valid for speeds below 500 m/s, while a quasi-ellipsoidal deformation regime emerges at larger speeds. Furthermore, for speeds up to 1000 m/s the particle undergoes a soft landing and creates a long-lived coherent surface phonon. Higher speeds lead to a rapid attenuation of the coherent phonon due to a partial diamond cubic to-tin phase transformation occurring in the particle.

of reducing the contact pressure. Recent nanoindentation data [12] confirmed that silicon nanoparticles can accommodate plastic deformations. Thus, following conventional wisdom, one conjunctures that plastic yield may occur upon impact, in agreement with macroscopic predictions [9]. However, one concern is that hertzian theory predicts that the duration of the impact decreases with the particle radius r 0 as r 0 6/5 [13]. Thus, the resulting strain rate experienced by the nanoparticle could be too high for inducing plastic deformations.
This paper is organized as follows: Section II gives a description of the employed MD scheme. Section III presents the obtained numerical results, which are summarized in Section IV.

II. Methodology
We conducted detailed atomistic simulations to describe the nonadiabatic process of the impact of Si nanoparticles on (001) Si surfaces. The method of investigation is classical MD, where the motion of the Si cores is described with the environment-dependent potential proposed by Tersoff [14], already incorporated in the computational package Trocadero [15]. The classical treatment of Si covalent bonding allows simulating rather substantial systems. Here we are dealing with spherical Si particles of radius r 0 =5 nm and containing N p =31,075 atoms, and a Si substrate with dimensions 16.3 nm x 16.3 nm x 5.4 nm containing 72,000 atoms. An O(N) linked-cell list strategy to search for the interacting atoms was implemented and the system was evolved in time with a 2 fs time step with a velocity Verlet algorithm.
The dynamical treatment of the finite-size substrate requires special care. We employed artificial absorbing boundaries [3,16] outside the region of interest. During the impact simulation the two bottom layers (3,600 atoms) of the substrate were kept fixed in time, the next ten layers (18,000 atoms) were evolved with a damped Langevin dynamics [3] with a the friction coefficient γ=πω D /6,  where ω D is the  Debye frequency for silicon, while the rest of the substrate is treated with standard MD dynamics. Periodic boundary conditions are applied to the horizontal x and y directions (z impact direction is vertical). In preparation for the simulations, one nanocluster has been spherically cut from bulk Si, annealed at T p =1000 K or T p =300 K and optimized with a conjugate gradient algorithm. A Si substrate has been cut from bulk along the (001) direction and a similar annealing/optimization procedure has been applied. The procedure leads to the expected (001) Si surface reconstruction by dimerization.
We attempted to choose simulation parameters (impact velocities, particle and substrate temperatures) and measure quantities that are relevant for the impact experiments of Rao et al. [6,7]. In a series of simulations, the speed of the 5 nm incident nanoparticle was taken in the range of 550-1640 m/s. To explore the role of temperature we carried several simulations on combinations of the cold (T p =300 K) and hot (T p =1000 K) particle impacting on a cold (T s =300 K) or hot (T s =800 K) substrate. For a better understanding of the process we found it necessary to conduct several simulations outside this main velocity-temperature range.
All simulations were conducted for 40 ps and were followed by structural optimizations performed with a conjugate gradient scheme. During the simulation we have monitored a standard set of indicators such as instantaneous particle temperature (proportional to instantaneous kinetic energy), the pair correlation function (essentially a density-density correlation function), the particle cohesive energy E p , and the particle-surface adhesion E a . When pieced together, these indicators should converge to provide a microscopic picture for the dynamics of the impact.

III. Simulation Results
In all our simulations the impacting particles remained attached to the surface. Fig. 1 depicts, for two incident speeds, three representative stages from the impact process: the initially undistorted structure, the instant of maximum penetration into the substrate, and the structure after 20 ps. The selected utmost impacting conditions vividly illustrate that even in this relatively low-energy impact regime there is a strong qualitative dependence on the impact speed: At 550 m/s the particle penetrates very little into the room-temperature substrate and overall suffers little structural damage following the impact. After the maximum penetration point, on the rebound stage, a small junction with a neck is formed at the particle-surface interface, as it can be observed in the last sequence. On the other hand, at an incident speed of 1640 m/s, the particle shows evident signs of structural changes even at the maximum penetration point. The spherical shape is partially lost in favor of a quasi-ellipsoidal shape, and a junction with a large interfacial contact surface is formed. In both cases, the substrate appears compliant and does not show signs of significant structural deformation or any degree of implantation with atoms from the particle.
The Contact Area and Adhesion Energy. An interesting aspect about the dynamics of the collision process transpires from Fig. 2(a), which displays the time evolution of the contact radius, a [17]. As the particle impinges onto the substrate, a increases up to the maximum penetration point instant marked by arrows. Next, for the lowest impinging speed case, a practically saturates at the value attained at the maximum penetration distance into substrate. The higher impact velocity leads to a less quasi-static behavior, as after the maximum penetration point the contact area slightly decreases on the rebound. Most importantly, this dynamical result suggests that the contact area is largely determined by the incident speed and the mechanical properties of both the particle and substrate, and to a lesser extent by the subsequent events that follow the maximum penetration.
A broader view about the role of the impacting velocity emerges from Fig.  2(b), which plots the final contact radius a as a function of the incident speed under room temperature conditions. To elucidate the speed dependence it was rewarding to perform additional low KE simulations, outside our main v p range of interest. After a careful analysis we concluded that for speeds below 500 m/s the impact largely follows the hertzian predictions: The measured contact  radii (empty circles) fits into the typical v p 2/5 dependence, which assumes a rather localized particle deformation at the maximum penetration point. (The free term correction reflects the atomistic aspect of the contact of nonconforming shapes.). Additionally, the dynamically measured a values are very close to the predictions (shown with empty triangles) made with the formula [9] a=π p m r 0 /2Y, where p m represents the mean contact stress at the maximum penetration and Y=65 GPa is the computed [18] lumped particle-substrate modulus. A qualitative scaling change can be seen in the a data above the 500 m/s, where the hertzian predictions become less valid. The better agreement with the stronger power v p 1/2 scaling suggests the impact enters into an extended deformation regime. As noticed in Fig. 1(b) the particle exhibits a quasi-ellipsoidal deformation at the maximum penetration point, similar with the experimental observation on the macroscopic impact of gel balls [19].
In discussing the energetics of the bound system E s-p it is instructive to distinguish between the particle-substrate adhesion energy E a and the individual cohesive energies. For this propose, following the concluding relaxation we have artificially detached the particle and separately computed the particle E p and substrate E s cohesive energies. Thus, E s-p = E a + E p + E s . The energetically favorable aspect of the deposition process is portrayed in the linear v p trend of the obtained adhesion energy E a , as presented in Fig. 3(a). A higher impact speed causes a larger contact surface, thus "eliminating" the less energetically favorable surface. We notice that the thermal conditions are secondary and E a is mainly determined by the incident kinetic energy.
Collision Modes. We next turn to a detailed analysis of the dynamics of the impact, aiming to identify characteristic collision modes. Particularly useful for this propose is the E p data presented in Fig. 3(b), where a threshold-like dependence on the particle speed can be noted. For collision speeds of up to ~1000 m/s, the E p variation with the nanoparticle speed exhibits a plateau. Beyond the ~1000 m/s value there is a rapid E p increase that clearly indicates that the nanoparticle suffered a significant structural change after landing. Although a clear-cut critical speed value was not identified, our data suggest that function of the incident speed there are two rather distinct regimes of impact, having very different consequences on the final structure of the nanoparticle. Again, the thermal conditions appear to influence the impact outcome very little.
To better understand these two distinct collision modes we have concentrated on the dynamical features exhibited by two extreme cases, namely a slow-cold impact at v p =550 m/s with T p =T s =300 K, and a fast-hot impact v p =1640 m/s with T p =1000 K and T s =800 K. To gain insight about the transfer of the incident kinetic energy we have plotted in Fig. 4 the instantaneous speed of the center of mass of the nanosphere. modeling the resulting oscillatory motion as an under-damped second-order system, i.e., v p (t) ~ exp(-ζ ω n t) cos(ω n (1-ζ 2 ) 1/2 t + φ).
Here ω n is the system's natural frequency, φ a phase factor, and ζ is the dephasing parameter. Fitting to the atomistic data we obtained ζ =0.04. Thus, at this speed most of the incident KE is transferred to the coherent transversal vibrations of the particle. Fig. 4b corresponding to an 1640 m/s incidence, indicates the creation of the same type of mechanical surface oscillations but the energy dissipation is much more efficient. This is reflected by the increase in the dephasing parameter, which now becomes ζ =0.37. We have further characterized the dynamics of the collision by monitoring the instantaneous particle temperature T p in Fig. 5 and the spatial T distribution across the system, in Fig. 6. Overall, in the low energy regime T p is simply enhanced to about 400 K. The spatial decomposition of Fig.  6a reveals that at the touchdown instant only small portion of the interface material reaches a temperature of about 550 K. Next, T p is seen to rise rather uniformly, meaning that part of the initial incident kinetic energy is converted into thermal energy that moderately heats up the system. However, a very different dependence is observed in the fast impact, Figs. 5 and 6b. After the touchdown instant there is a sharp T p rise to 1450 K. This is followed, during the next few ps, by an ample T p increase up to a maximum of 1600 K. For the remaining 30 ps of the simulation T p slowly drops, as the system is driven towards equilibrium by the Langevin dynamics of the substrate. Fig. 6b clearly indicates that the early global T p peak can be attributed to the localized increase of temperature at the interface, which can reach values up to 2400 K. A distinguished behavior is seen during the next few ps: Instead of a homogenization, the 2500 K front expands into the nanoparticle only, specifically in the conical region corresponding to where the amorphization is found at the end of the simulation as seen in Fig. 1(b). As in the slow impact case, a rather uniform T distribution can be seen towards the end of the simulation. Particularly relevant are the results for the time evolution for the contact stress p m , shown in Fig.  7. Again, there are obvious differences in the two impact cases: While p m exhibits large oscillations for the slow impact, in the fast collision the first pressure peak is quickly dampened and, shortly after the impact, the stress at the interface becomes zero. In both cases, very large pressures are reached, around 13 GPa for the slow-elastic regime, and around 18 GPa for the fast-damped one. Such stresses are amply within the range of pressures capable of determining a phase transition in silicon, i.e., 10-15 GPa for the cubic diamond β-tin phase transformation [20,21].

m/s
Microscopic Mechanism for Inelastic Collision. What is the microscopic mechanism causing the above-threshold behavior? To answer this question we have analyzed in detail the structural dynamics of our system. We found that the extreme pressures reached at the nanoparticle-substrate interface provide the explanation for why such a localized increase in temperature takes place only in the nanosphere, leaving the substrate almost unaffected. Specifically, the nanosphere responds to the high pressure-rates experienced during the initial compression stage of the impact with creation of a new β-tin phase in the most strained regions. This way, a large amount of the elastic energy accumulated during compression is not returned as elastic energy for the particle to rebounce.
Moreover, this relaxation is accompanied by a KE release, as noted in Fig. 5.
To better illustrate this mechanism we have performed an additional simulation with the particle impacting with 2200 m/s under very low temperature conditions (10K for both the particle and the substrate). The evidence for a dynamical structural transition is presented in Fig. 8, where the fraction of nanoparticle's atoms having 4, 5, 6, and 7 nearest neighbors of the nanoparticle is plotted for the first 20 ps of the impact. Almost 90% of the atoms are initially 4-fold coordinated, as expected since most of the particle has a cubic diamond structure; however, during the impact, the fraction of atoms having a higher coordination number (up to 7) increases rapidly, and at around 20 ps more than 40% of the particle is made up of atoms with a coordination number between 5 and 7, with less than 15% of the atoms having only 4 nearest neighbors. Note that the βtin phase is characterized by a coordination number of 6.  The precise evidence of an ephemeral β-tin phase is given in Fig. 9 showing the radial distribution function g(r) obtained from the nanoparticle structure at t = 3.9 ps after a structural optimization. In the insert we have also plotted the g(r) for the perfect infinite β-tin crystal computed based on the lattice parameters of Ref. [22]. Aside from the first peak caused by the residual original structure, located at roughly 2.35 Å (i.e., the nearest neighbor distance in the Si cubic diamond lattice), the next spikes well correspond to those obtained from the perfect infinite βtin crystal.
As the simulation progressed, we have found that the rebounce stage produces structural amorphization, leading to a further release of kinetic energy. In fact, hardly any trace of the β-tin phase is found at the end of the simulation. The observed amorphization during the rapid upward bouncing motion appears in agreement with the amorphization noted in bulk silicon during a rapid retraction of the nanoindentor [23,24].

IV. Summary
We carried out MD simulations using the Tersoff potential to study the deposition process of a 5 nm Si nanoparticle onto a (001) Si surface. Our goal was to gain insight into the particle adhesion and the mechanism of dissipation in the KE range of below 1 eV/atom. All simulations clearly indicated the deposition speed as the important parameter and the particle and/or substrate temperature conditions as secondary.
We found it rewarding to examine the same aspects as in a macroscopic impact [9,19]. In all our simulations the particle remained attached to the surface. We quantified the degree of adhesion through the microscopic contact radius a, and revealed some important regularities: The contact radius at the end of the simulation is equal with the contact radius initially attained at the maximum penetration point. For velocities below 500 m/s, the obtained scaling a ~ v p 2/5 is in agreement with hertzian predictions. Above this speed the scaling law changes to a ~ v p 1/2 as the incident particle underwent extensive ellipsoidal changes during the impact. Next the particle-surface adhesion energy showed a linear scaling E a ~ v p .
Regarding the aspect of energy dissipation we identified two distinct regimes: For speeds below ~1000 m/s the impact is mainly elastic, and the nanoparticle structure and shape are well preserved. The impacting particle produced a transversal coherent phonon mode with a small dephasing time. For a more violent impact, above ~1000 m/s, there is an efficient transfer of the incident KE into the heating of the particle. A distinguished nanoscale feature was revealed: in the inelastic collision the KE transfer is not mediated by a plastic yield as in the macroscopic case but by an ephemeral cubic diamond to β-tin phase transformation occurring in the nanoparticle during the initial stages of the impact. Fig. 9. The radial distribution function g(r) for a 5 nm particle at the 3.9 ps maximum penetration instant. The vertical lines correspond to the bulk β-tin phase.
Beside the fundamental aspect of describing at the microscopic level the impact in a new velocity -nanosize regime, our study contributes towards the understanding of nanoparticulate materials, especially by way of a hierarchical multiscale modeling anchored on the presented atomistic data.