Dissipation and plastic deformation in collisions between metallic nanoparticles

The collisional dynamics of nano- and micron-scale particles is often described theoretically using the continuum Johnson-Kendall-Roberts (JKR) theory, which describes adhesion between elastic spheres[1], with the addition of dissipative processes. Various physical mechanisms have been proposed for dissipation, including elastic waves, viscoelastic dissipation, and dissipative mechanisms related to crack formation[2, 3, 4, 5, 6, 7]. While JKR has typically been compared favorably with experimental results[8, 9, 10], often the surface energy term is not well known and in some cases is fit to experimental results[11, 2]. In addition, no experimental results have clearly demonstrated which dissipative mechanism provides the most suitable description[11]. Finally the JKR model is not able to account for plastic deformation, which should become relevant especially at higher collision velocities.

Recently, atomic-scale simulation has been used to test some assumptions of the JKR model. In Nietiadi et al.[12], molecular dynamics (MD) simulations of collisions between amorphous silica nanoparticles demonstrated strong deviations with the JKR model. Specifically, the critical bouncing velocity was established by simulations to be a factor 3.4 greater than predicted. Furthermore, only sticking was observed for nanoparticles with radii less than 15 nm. This enhanced sticking behavior appeared to be connected to the presence of a larger contact area due to strong plastic deformation, including the generation of filaments, between the nanoparticles. In agreement with these results, it was shown in Quadery et al.[13] that using MD collisions of silica nanoparticles, the lack of adsorbed OH groups results in strong covalent bonds, with no clear bouncing threshold. However, in Quadery et al.[13], only small nanoparticles with radii 2 nm and below were simulated. Nevertheless, existing simulation results suggest strong deviations from JKR model predictions at least for very small nanoparticles.

The adhesive behavior of powders in turbulent flows is of industrial and theoretical interest as can be seen by the wide range of current research and review articles [14, 15, 16]. Most of these models utilize the JKR model of contact mechanics in some form. While Fe nanoparticles might not typically be a material considered for simulating collision dynamics, metallic nanoparticles, and specifically iron nanoparticles, are of interest due to their enhanced catalytic activity, which has been demonstrated for industrial processes [17, 18, 19, 20] and in the context of promoting chemical reactions of astrophysical interest[21]. Iron nanoparticles can even be used to assist in environmental cleanup efforts[22]. Metallic Fe nanoparticles represent a simple system and are a kind of extreme case where due to dangling bonds at the surface one might expect particularly strong dissipation and adhesion, and thus it represents a sort of “limiting case” in collision dynamics which would be interesting to understand. Velocities were chosen to span a range of regimes but with a special focus on lower velocity collisions from 10 to 500 m s^-1 which are relevant in the astrophysical context of protoplanetary dust cloud dynamics. For example, in astrophysics, additional dissipative mechanisms and enhanced adhesion could be helpful in surmounting the so-called millimeter bouncing barrier[23], whereas in a catalytic context the sintering of nanoparticles and consequent loss of surface area leads to a significant drop in effectiveness. Thus, we hope to offer suggestions towards the development of more physically relevant models of adhesion and dissipation for a range of interaction environments.

In the next section, the basic assumptions of the JKR model applied to nanoparticle collisions with various models of dissipation are briefly described. Next, the atomic-scale simulation methodology used in this paper is presented, followed by a section describing the results of extensive simulations. In the discussion section we attempt to determine scaling relations for relevant physical quantities that can potentially yield predictions for significantly larger length scales, including up to 1 $\mu$m particles. Finally, in the conclusion we tie these results back to potentially relevant applications, including planetary formation, other mineral systems, and previous findings in the field of tribology, including potential relationships to the modern understanding of Amonton’s laws of friction.

In the JKR theory[1], the energy associated with adhesion between two identical spheres of radius $R$ with reduced radius $R^{*}=\frac{R}{2}$ depends on the interfacial energy and the elastic strain energy. The interfacial energy $U_{S}$ depends on the surface energy $\gamma$ and contact radius $a$,

while the elastic strain energy $U_{E}$ is given by

For two identical spheres, the combined elastic modulus $E^{*}$ is defined by

where $E$ is the Young’s modulus and $\nu$ is the Poisson ratio of the material in bulk. The quantity $\delta$ is the length associated with compression of the two spheres in contact. Specifically, for two identical spherical objects with radius $R$ whose center-of-mass coordinates are separated by a distance $d$, the compression is given by

For a particular compression $\delta$, JKR theory predicts that the system will optimize the contact radius $a$ to minimize the total energy $U_{JKR}=U_{S}+U_{E}$, resulting in a relationship between the compression $\delta$ and contact radius $a$,

With this assumption, the JKR theory results in a force

The point where $F_{JKR}$ vanishes yields the equilibrium contact radius and length of compression,

As can be seen in the above expressions, the JKR theory applied to collisions predicts that equilibrium is attained with a contact area and compression (and thus an adhesive energy) that are independent of the collision velocity. These assumptions are not valid if significant plastic deformation occurs.

In order for adhesion to occur, complete dissipation of the incident kinetic energy is required. Dissipation always involves generation of internal thermal energy, while for large enough collision velocities, plastic deformation, generation of coordination defects, and melting can occur. Several attempts have been made to add dissipation to existing JKR-based adhesion models. For example, in Krijt at al.[2], models of viscoelastic growth of cracks, bulk viscoelastic dissipation, and plastic deformation were used to describe dissipation. The resulting model provides theoretical predictions for the coefficient of restitution and the sticking velocity. In Chokshi et al.[7], dissipation into elastic waves was used to describe dissipation, with the critical sticking velocity predicted from the requirement that the energy dissipation in elastic deformations be greater than the energy required to separate the nanoparticles. However, while many models for dissipation exist, no results have been reported which validate any particular dissipation model.

In the present article, the detailed atomic-scale mechanisms for dissipation are described using molecular dynamics (MD) simulation. This approach includes all length scales for dissipation, including vibrational modes with wavelengths comparable to the separation between the atoms. The advantage of this approach is that there is no requirement for a physical model with assumptions, but rather all atomic degrees of freedom are explicitly described.

The LAMMPS simulation code[24] with the embedded-atom method (EAM) potential for Fe from Mendelev et al. (potential 2)[25] was used. Visual renders of atomic structures were produced using the Visual Molecular Dynamics (VMD) software package[26]. The simulations were performed with an MD time step of 0.25 fs, small enough to ensure energy conservation, and much smaller than the timescale of the highest interatomic vibrational frequency in the system. Three sizes of amorphous nanoparticles were generated: a small nanoparticle with N = 1024 atoms and radius R = 1.4 nm, a medium nanoparticle with N = 50286 atoms and radius R = 5.2 nm, and a large nanoparticle with N = 470561 atoms and radius R = 11 nm. The Fe nanoparticles were melted in a constant temperature ensemble by increasing the temperature to T = 2200 K. This was followed by a slow anneal to T = 5 K. Melting was performed over 100 ps for the small nanoparticle, 120 ps for the medium nanoparticle, and 240 ps for the large nanoparticle. Annealing times were identical to melting times. In the top right of Figure 1, the atomic structure of the large amorphous Fe nanoparticle is depicted. After the annealing process was complete, the Fe-Fe radial distribution function (RDF) was calculated to ensure that the nanoparticles were amorphous. The RDF of the large nanoparticle is depicted in the bottom left of Figure 1. From the RDF data, the coordination number was determined by integration to the first minimum of the RDF at 0.33 nm. The calculated Fe-Fe coordination number of 13.2 is somewhat higher than for an FCC crystal due to the fact that the first minimum in the RDF is significantly greater than the Fe-Fe bond lengths $\sim$ 0.25 nm in either BCC or FCC iron.

The radii of the nanoparticles were determined first by computation of the $T=0$ K density of bulk amorphous Fe, and then assuming a spherical shape for the nanoparticles. For a nanoparticle with mass $m$ and mass density $\rho$, the radius is defined by

Using this expression, with the computed mass density $\rho=$ 7.82 g cm^-3 for amorphous Fe, the small nanoparticle with $N=$ 1024 corresponds to a radius $R=$ 1.4 nm. For the medium nanoparticle with $N=$ 50286, the radius is $R=$ 5.2 nm. Finally, for the large nanoparticle with $N=$ 470561 we obtained the radius $R=$ 11 nm. Given the radii of the nanoparticles, the surface energy $\gamma$ was determined using the computed energy of the nanoparticles and that of the bulk amorphous Fe solid. In each case, the surface energy was computed to be nearly $\gamma=$ 0.09 eV $\text{\AA }^{-2}$, indicating that there was no size-dependence to the surface energy.

Collisions between same-sized nanoparticles with relative velocities between 10 m s^-1 and 3000 m s^-1 were simulated. These velocities were chosen to span the entire range from moderately slow collisions up through collisions with enough kinetic energy to ensure a liquid final state. We cast a special focus on lower velocity collisions up to 500 m s^-1, as these velocities were found to result in minimal large-scale plastic deformation of the nanoparticle structure away from the inter-particle contact. For each relative velocity, multiple collisions were simulated for different random rotations of the two nanoparticles. For the small and medium nanoparticles, 30 simulations were performed at each velocity. For the large nanoparticle, 3 simulations were performed at each velocity. After relaxation, both nanoparticles were given a desired translational velocity which resulted in a head-on collision. Simulations were continued for at least 50 ps after the collision until equilibrium was achieved. The center-of-mass coordinates of the two nanoparticles were monitored to determine whether the nanoparticles adhered together. For each simulated nanoparticle size and incident velocity only sticking was observed, with no incidents of bouncing behavior.

To compute the work of adhesion $W_{adh}$ from the results of the irreversible collision simulations, the approach first described in Quadery et al.[13] was used. The basic assumption is that nanoparticle collisions result in a change in potential energy and thermal excitation. If the system remains in a solid state, it is reasonable to assume that the thermal energy can be described using the equipartition theorem applied to a system of harmonic oscillators. The work of adhesion $W_{adh}$ was therefore computed according to,

where $N_{tot}=2N$ is the total number of atoms in the simulation, $T_{f}$ and $T_{i}$ are the computed temperatures before and after the collision, and $K_{trans}$ is the translational kinetic energy before the collision. In determining $T_{f}$ and $T_{i}$, the kinetic energy in the center-of-mass reference frame of each nanoparticle was used. The physical meaning of $W_{adh}$ is that it approximately corresponds to the work required at $T=0$ K to adiabatically separate the nanoparticles by breaking bonds at the interface. While the interpretation of $W_{adh}$ is clear when adhesion occurs with minimal disruption of the surfaces, in instances with large plastic deformation or melting of the nanoparticles, significant amounts of energy can be stored in disruption of the atomic structure, and the concept of adiabatically separating the nanoparticles after the collision is somewhat ill-defined. Nevertheless, $W_{adh}$ is calculated using Eq. 10 for all collisions including those where the physical interpretation is more complicated.

Each simulated head-on collision resulted in sticking. Fragmentation was not observed at any of the simulated velocites or radii. For very large velocities, full melting was observed: for v_rel = 2750 m s^-1, T_f was around 1900K, and for v_rel = 3000 m s^-1, T_f was around 2200K, comparable to the solid to fully liquid transition temperature T_liq $\sim$ 1900 K of a single $R$ = 1.4 nm nanoparticle melted in an ancillary MD simulation. Consequently, there is no velocity range where bouncing of solid nanoparticles occurs, although at very large velocities either vaporization or splashing will occur. It may be that bouncing events are not impossible, but that they are at least extremely rare. This is in significant contrast to previous results for SiO₂ nanoparticles in Quadery et al.[13] where significant instances of bouncing were observed. Bouncing might also occur for very large nanoparticle sizes or when collisions are not exactly head-on.

In Figures 2, 3, and 4, we show the computed average values of $W_{adh}$ as a function of collision velocity $v_{rel}$ for lower velocities, with an inset showing the data for all velocities. The error bars represent the standard deviations obtained from all simulations performed at each value of $v_{rel}$. Along with the computed values, an image of a typical atomic structure for the adhered nanoparticles is shown for each collision velocity. Each structure depicted represents the final structure attained after equilibration. Predictions made by JKR for elastic, surface, and total energies are denoted by lines on Figures 2, 3, and 4,. Comparison to JKR theory will be presented in the next section.

From the data plotted in Figures 2, 3, and 4, several consistent trends emerge which are independent of either the model or the size of the nanoparticle. First, the computed value of $W_{adh}$ increases gradually with increasing collision velocity for values of $v_{rel}$ up to at least 1000 m s^-1. Corresponding to these cases, the structures in Figures 2, 3, and 4, show an apparent contact area which increases with $v_{rel}$. For larger values of $v_{rel}$ (i.e. significantly past 1000 m s^-1), the computed $W_{adh}$ begins to decrease and eventually becomes negative. In this regime, the atomic structure appears less like two adhered nanoparticles and progressively more like a single deformed nanoparticle, generally elliptical in shape, eventually becoming spherical at the highest simulated values of $v_{rel}$. Based on calculations of the self-diffusion coefficient, values of $v_{rel}$ above about 2000 m s^-1 were in a liquid state when fused into a spherical shape. For the largest nanoparticles with radius $R=$ 11 nm, only velocities $v_{rel}=$ 250 m s^-1 and below were simulated, and the increase in $W_{adh}$ with $v_{rel}$ is apparent but less dramatic. The simulations for the large nanoparticles with $R=$ 11 nm did not extend into the regime where large deformation and melting occurs due to the extensive computational resources required.

The general trends can be understood in a simple way. The changing nature of the contact area is indicative of significant plastic deformation. Specifically, at lower velocities (i.e. below $\sim$ 1000 m s^-1) the gradual increase in $W_{adh}$ with increasing $v_{rel}$ is due to plastic deformation which allows for greater contact between the two surfaces. Although this could occur by elastic deformation, it is important to note that elastic deformation would also involve significant elastic strain energy (see analysis in next section for further discussion). At higher velocities (above $\sim$ 1000 m s^-1), the decrease in $W_{adh}$ is due to the complete fusing of the two nanoparticles along with the generation of significant numbers of coordination defects, and possibly strain energy. The coordination defects correspond to stored energy, and thus result in negative contributions to $W_{adh}$. In these instances, separation of the two nanoparticles would result in nanoparticles with markedly different structures than before the collision. Nevertheless, $W_{adh}$ is a measure of the internal potential energy of the system with respect to the isolated nanoparticles before the collision. Finally, at the highest simulated values of $v_{rel}$, the change in internal energy is large enough to result in a phase transition to the liquid state along with very large negative values of $W_{adh}$.

If elastic strain energy and energy associated with coordination defects are neglected, it is possible to use the computed values of $W_{adh}$ to establish an effective contact area

As in Quadery et al.[13], we also compare this effective contact area to the cross-sectional area via

The unitless parameter $\eta$ captures both the relative area of the interface as well as the quality of the bonding. Specifically, values of $\eta\approx 1$ imply a perfectly-coordinated interface with both nanoparticles deformed such that the interfacial area corresponds to the entire cross-sectional area. By contrast, values of $\eta$ that approach 0 imply a weakly bonded interface, either with a small contact area or a significant number of defects. Generally, it is expected that $\eta$ will have a value intermediate to these extremes. From the data plotted in Figure 5, at $v_{rel}$ = 10 m s^-1, the value of $\eta=0.29$ for $R=$ 1.4 nm nanoparticles indicates a substantial contact area at the interface between the nanoparticles even at the lowest values of $v_{rel}$. For larger nanoparticles, the values of $\eta$ are significantly smaller, yet still are indicative of substantial contact and strong bonding. Additionally, the final compression length $\delta$ was directly determined from the MD simulations and Eq. 4. In Figure 6 $\delta$ is plotted as a function of $v_{rel}$ for each radius $R$. It is clear that $\delta$ increases strongly with increasing $v_{rel}$. This observation is consistent with plastic deformation.

To establish the dissipation mechanism, we will now focus on the low-velocity collisions. In all collisions, the center-of-mass coordinates of each nanoparticle were retained as a function of time, allowing for determination of the acceleration and hence net force during the collision. In Figure 7, the velocity and acceleration of both nanoparticles are plotted as a function of time for one simulated collision of two R = 1.4 nm nanoparticles with $v_{rel}$ = 100 m s^-1. The collision just after 20 ps is evident; the spike is an artifact of how the initial translational velocity was imparted. Due to the strong attraction, the nanoparticles initially accelerate towards each other once they enter the interaction range determined by the cutoff of the EAM potential. After significant compression, the acceleration changes sign and the velocities slow. However, before the velocities change sign for the rebound, the forces already have begun to decrease. If this were an elastic deformation, any increased compression of the two particles would lead to increased force as more elastic energy is stored. The point where the velocities change sign corresponds to the rebound phase. However, during the rebound phase, the magnitude of the acceleration is dramatically decreased from the compression phase, and consequently the rebound velocities are quite small. These observations clearly show that dissipation is correlated with strong plastic deformation, and in fact the plastic deformation itself represents the primary dissipation mechanism.

After the initial deformation and rebound phase, the subsequent behavior of the center-of-mass coordinates is consistent with that of an underdamped simple-harmonic oscillator. The results for the center-of-mass velocity (e.g. Figure 7) were used to establish the oscillation period $\tau$ for collisions with $v_{rel}$ = 10 m s^-1. As expected, the period increases with increasing radius $R$. In Figure 8, the kinetic energy in the oscillations of two nanoparticles of mass $m$ and center-of-mass velocity $v_{CM}$, $K_{CM}$ $=$ $m$ $v_{CM}^{2}$, is normalized by the energy to be dissipated and plotted as a function of time again for $v_{rel}$ = 10 m $s^{-1}$ collisions. It is evident from Figure 8 that the energy of the vibrational motion of the particles is very strongly damped. Exponential decay function fits were performed to determine the damping time $\tau_{d}$. In Table 1, the values of $\tau$ and $\tau_{d}$ are given for each radius $R$ for $v_{rel}$ = 10 m s^-1 collisions. For the particles to bounce, the requirement would be that $\tau_{d}>>\tau$, so that most of the incident energy is available during the rebound phase. In each case, $\tau$ is significantly greater than $\tau_{d}$, demonstrating very strong damping; however, it is evident that $\tau_{d}$ is increasing with $R$ faster than $\tau$, indicating the potential for bouncing at large enough $R$ values. In the next section, this condition will be explored to determine when bouncing might be expected to occur.

The stresses at the interface can be estimated from the computed forces as well as a reasonable estimate of the contact area. Using the values of $\eta$ shown in Figure 5 as an estimate of the contact area, the stresses during the collision, both tensile and compressive, are in the range of $3-10$ GPa in magnitude. This is significantly greater than the yield stress for bulk Fe of 80-100 MPa[27]. While yield stresses of nanoparticles can be somewhat larger than for bulk materials the differences are generally fairly small. For example, in Hawa et. al[28], simulations were used to compute the yield stresses for Ag nanoparticles, with yield stresses in the range 0.5-0.7 GPa for crystalline nanoparticles with $R\sim 5-10$ nm, and somewhat lower yield stresses for amorphous nanoparticles. Therefore, even taking into account the higher yield stresses typically exhibited by nanoparticles, the stresses exerted upon Fe nanoparticles in the simulation are large enough to result in plastic deformation. In Chokshi et al., the authors briefly discuss the sizes below which plastic deformation should be relevant for dissipation[7]. Following their arguments, for the currently presented simulations plastic deformation should only be important for sizes smaller than roughly R $\sim$ 1.3 nm, which is significantly smaller than the R = 11 nm particles reported here. However, it is also clear that more accurate calculations of stress would be worthwhile, including a calculation of the stress gradients in the vicinity of the contact region.

The results presented above demonstrate some important features. Specifically, at least in the case of particles of $R\leq$ 11 nm, plastic deformation and extremely strong dissipation occurs in a manner not consistent with the JKR theory. Though much of the interest in particle interaction dynamics lies in larger particles of at least $R\sim 1\mu$m and beyond, where direct MD simulations are not possible, some behaviors can be predicted based on how various quantities scale with $R$. In this section we first strengthen the understanding of how JKR fails at nanometer length scales, and then explore how the results scale with radius $R$ to determine how MD predictions can be used to develop theoretical understanding at larger length scales, including where JKR likely has more validity.

In the JKR model, elastic strain accommodates an increase in the contact area. In addition, it predicts that the final contact area does not depend on the collision velocity $v_{rel}$ and that above a certain velocity, the initial translational kinetic energy is too great to cause adhesion. In order to elucidate which contributions to $U_{JKR}=U_{S}+U_{E}$ were most important in determining the value of $W_{adh}$ calculated as per Eq. 10, we used the directly computed values of $\delta$ to numerically solve Eqn. 5 and to get an expression for contact radius $a$ as a function of $v_{rel}$. These contact radii were then used in Eqs. 1 and 2 to obtain values for the contributions to $U_{JKR}$. These are plotted on Figures 2, 3, and 4, for comparison to the MD simulation results for $W_{adh}$. Clearly, the surface energy contribution accounts very well for the values for $W_{adh}$ for velocities up to $v_{rel}$ = 500 m s^-1. Figure 6 shows that $\delta$ increases with $v_{rel}$, resulting in an increase in contact radius $a$ and consequently an increased magnitude for the surface energy contribution $U_{S}$. The JKR prediction for the elastic energy component $U_{E}$ is not consistent with the trend shown by $W_{adh}$. These observations present a clear picture of plastic deformation as the mechanism responsible for the increased contact area. Hence, both $\delta$ and $a$ are found to increase with increasing $v_{rel}$ in a manner which could not occur if the deformation included elastic strain energy. In addition, the computed values of $\delta$ plotted in Figure 6 are significantly greater than the predictions of JKR. Specifically, JKR theory predicts values of $\delta_{eq}=$ 0.23 nm, 0.35 nm, and 0.43 nm for the small, medium, and large nanoparticles, respectively, which in all cases underestimate the values shown in Figure 6. This demonstrates enhanced compression beyond the predictions of JKR, consistent with the observation of plastic deformation.

In Figure 9, we plot values of $\eta$ for three velocities as a function of nanoparticle radius $R$. The gradual increase in $\eta$ with $v_{rel}$ is consistent with the observation of greater plastic deformation. The decrease in $\eta$ with increasing radius can be understood to be in part simply a geometric effect. Because the interactions are finite range, when $R$ increases significantly beyond 0.53 nm, the cut off distance for interactions, simple geometric arguments suggest that $\eta$ should scale as $R^{-1}$. The actual data shows a trend with a somewhat different scaling exponent, likely due to the fact that the strongest plastic deformation happens for the smallest nanoparticles. The value of $\eta$ appears to scale with radius $R$ approximately as $R^{-0.83}$ for $v_{rel}=$ 10 m s^-1. For higher $v_{rel}$, the scaling changes somewhat. Specifically, for $v_{rel}=100$ m s^-1 , $\eta$ scales as $R^{-0.72}$, while for $v_{rel}=250$ms^-1 , $\eta$ scales as $R^{-0.57}$. The dependence on $v_{rel}$ is clearly due to the fact that $K_{trans}$ becomes dominant in comparison to surface interaction as $v_{rel}$ increases. For low enough velocities the contact area indicated by $\eta$ appears to increase approximately linearly with $R$. The relevance for collisions at very large scale, both in terms of the dissipation mechanism and the crossover towards bouncing behavior, will be addressed in the final section. However, we note that for $v_{rel}$ = 10 m $s^{-1}$ and R = 1 $\mu m$, the scaling behavior results in a prediction $\eta\approx 1.2\times 10^{-3}$, which indicates that the adhesion and dissipation occurs over a much smaller relative area than for $R=11$ nm and smaller particles.

For small particles $R=11$ nm and below, the results indicate that plastic deformation always occurs. Indeed, we observe that is not possible at these scales to lower $v_{rel}$ sufficiently to observe elastic behavior consistent with JKR theory. For small particles, surface attraction energy dominates $K_{trans}$ for the lower values of $v_{rel}$. The result of the strong attraction is that the nanoparticles accelerate significantly just before the collision, and the effective collision velocity $v_{rel,f}$ can often be substantially greater than the initial velocity $v_{rel}$. Assuming that the surface interaction is conservative, then the effective collision velocity $v_{rel,f}$ should depend on the initial velocity $v_{rel}$ and the size-dependent velocity $v_{c}$,

The value of $v_{c}$ is a parameter which depends on $R$ and was determined by examination of the maximum nanoparticle velocity during the collision. In Figure 10, $v_{rel,f}$ obtained from simulation is plotted as a function of $v_{rel}$ along with the fit curve from Eq. 13. All data was used to obtain fits, but only a subset is plotted for clarity. For the three radii $R=$ 1.4 nm, $R=$ 5.2 nm, and $R=$ 11 nm, the values for the fitted parameter $v_{c}$ are respectively $v_{c}=$ 213.9 $\pm$ 4.9 m s^-1, $v_{c}=$ 74.3 $\pm$ 1.7 m s^-1, and $v_{c}=$ 41.2 $\pm$ 1.5 m s^-1. This demonstrates that, for these small particle sizes, the velocity at collision $v_{rel,f}$ is substantially greater than the lowest value of $v_{rel}$ simulated. Consequently, a simulation with a very small $v_{rel}$ will result in $v_{rel,f}\approx v_{c}$ as a minimum effective collision velocity. Therefore, while it might be thought that a low enough value of $v_{rel}$ should exist where collisions are elastic, the present results demonstrate that for $R=11$nm and below, this is not the case. Specifically, for $R=11$nm, $v_{c}=$ 41.2 $\pm$ 1.5 m s^-1 is substantially greater than the lowest $v_{rel}$ = 10 m s^-1, and even much lower values of $v_{rel}$ would yield essentially the same collisions with plastic deformation. For smaller particles surface attraction is even more important, and the very large values of $v_{c}$ result in even more dramatic plastic deformation at all values of $v_{rel}$. In understanding scaling behavior, we note that $v_{c}\propto R^{-0.81}$. This indicates that as $R$ increases, surface attraction becomes a less significant factor. However, even for $R$ = 1 $\mu$m particles, the scaling of $v_{c}$ predicts $v_{c}\approx$ 1 m s^-1. Therefore, even if $v_{rel}$ is below 1 m s^-1, the effective collision velocity will be $v_{rel,f}\approx$ 1 m s^-1, which is still substantial and should involve some plastic deformation at the interface. This is further discussed in the last section of the paper.

It is possible to use scaling relations to predict the size $R$ where bouncing will occur. Specifically, one criterion for bouncing is that the kinetic energy during the first rebound should be larger than $W_{adh}$ in order for the particles to separate. We apply this criterion for $v_{rel}$ = 10 m s^-1 collisions. The results above for $\eta$ demonstrate that $W_{adh}\propto R^{1.1}$, whereas the kinetic energy during the first rebound scales $K\propto R^{3.4}$, hence eventually the rebound kinetic energy will be substantially greater than $W_{adh}$. For $v_{rel}$ = 10 m s^-1, this criterion results in the prediction that bouncing might occur for $R\approx$ 23 nm. However, the fact that substantial plastic deformation occurs casts some doubt on these results, since the neck formed by the deformed surfaces means the behavior is strongly irreversible. In fact, as the results plotted in Figures 2, 3, and 4, demonstrate, even negative values of $W_{adh}$ occur without bouncing, due to the very strong deformation of the particles.

Another approach to predict behavior at larger scales is to use the computed scaling of the oscillation period $\tau$ and the damping time $\tau_{d}$. As described above, bouncing would seem to require $\tau_{d}>>\tau$ in order for bouncing to occur, since strong dissipation of the kinetic energy of the two particle oscillations prevents separation of the particles. In Figure 11, the scaling of $\tau$ and $\tau_{d}$ with particle size $R$ is shown. From Figure 11, the fit indicates $\tau\propto R^{1.23}$, and $\tau_{d}\propto R^{2.50}$, which leads to a crossover at $R=17.3$nm. Beyond $R\approx$ 20 nm, $\tau_{d}$ eventually becomes substantially greater than $\tau$ and eventually bouncing should occur.

For metallic nanoparticles of radii $R=$ 11 nm and below, simulation results indicate adhesion occurs in every head-on collision. Simulation results show extremely large stress values and strong plastic deformation, indicating that JKR is not applicable in the present case. Because the observed deformation is plastic, it does not include large elastic strain energy contributions predicted by JKR, thereby resulting in a larger contact area that increases strongly with increasing $v_{rel}$. Moreover, the primary mechanism for the dissipation appears to be connected strongly to the atomic rearrangement that occurs at the interface associated with the plastic deformation. Scaling behavior suggest that for relatively low collision velocities, $v_{rel}\sim$ 10 m s^-1, nanoparticles would need to be at least $R\sim$ 20 nm in order for bouncing to occur.

For nanoparticles, surface attraction can often dominate incident kinetic energy, and plastic deformation is expected to always occur for any $v_{rel}$. The results show that even for $R\approx 1\mu m$, the relative velocity during the collision is at least $v_{rel,f}\approx$ 1 m s^-1. Using the scaling of $v_{c}\propto R^{-0.81}$ at $v_{rel}=$ 10 m s^-1, and the fact that the mass of a grain scales as $R^{3}$, the minimum kinetic energy in a collision scales as $R^{1.38}$. However, the scaling of $\eta$ at the same $v_{rel}$ shows that the kinetic energy associated with $v_{c}$ needs to be dissipated over an effective area which scales as $R^{1.1}$. This suggests that the minimum kinetic energy associated with $v_{c}$ increases more rapidly with $R$ than the area available to dissipate the energy. Therefore, it is quite possible that while dissipation becomes less effective as $R$ increases, and as expected bouncing becomes the dominant behavior, strong plastic deformation at the contact interface likely occurs. Hence, while $v_{c}$ decreases with $R$, the kinetic energy associated with $v_{c}$ actually increases, and because the relative area for the collision tends to decrease, it is reasonable to expect high stress and plastic deformation at the contact. However, this possibility remains a subject requiring more direct verification.

There is a direct relationship to the conventional explanation of Amonton’s laws of friction, wherein dissipative frictional forces are independent of apparent contact area and are solely dependent on the normal force at an interface. This has been contradicted based on the presence of nanoscale asperities[29] which result in an actual contact area which is generally much smaller than the apparent contact area. When normal forces exist, either due to surface attraction or some other applied force, the actual contact area grows often due to plastic deformation of the asperities, and hence the frictional force increases. When the actual contact area is much smaller than the apparent contact area, the stresses in the asperities can become quite large. This is very similar to the results found here. Specifically, since $\eta\propto R^{1.1}$ is less than $R^{2}$, the predictions here indicate that as $R$ increases, the stress at the contact point will tend to increase quite dramatically, thereby leading to plastic deformation. However, the results also show that any enhanced dissipation with increasing $R$ is not enough to prevent bouncing behavior at larger values of $R$.

In addition, the strong attraction which results in the larger effective collision velocities is the same reason plastic deformation occurs. In other words, even when the incident kinetic energy is relatively low, strong interaction tends to result in plastic deformation. As $R$ increases, surface attractive forces do not become weaker. Instead, the area where strong interactions occur increases less rapidly with $R$ than the overall mass of the particles does, and hence the value of $v_{c}$ decreases with increasing $R$. However, because the interactions are localized over a smaller relative area, they can still be very strong even when $v_{c}$ is small. In fact, we expect that not only stresses at the contact are very strong, but that very strong stress gradients likely are responsible for the observed plastic deformation. This point will be a focus of future efforts.

While the general picture of plastic deformation as the dominant mechanism for adhesion and dissipation is in stark contrast to JKR, some features remain consistent. Specifically, the simulations demonstrate the formation of a “neck” which tends to increase the adhesion and more strongly prevent rebound. The distinction is in the mechanism for the formation of the neck, which we find to be plastic deformation rather than elastic deformation. This view is also consistent with previous efforts in simulations of silica particles which demonstrated strong plastic deformation [12, 13]. It should also be noted that other works have explored viscoelastic dissipation mechanisms associated with plastic deformation[2], even though it appears that the dependence of the contact area on $v_{rel}$ and $R$ has not been previously considered. These insights could be of particular interest to the nuclear and pharmaceutical industries, where critical processes depend on tightly controlled powder flows. Additionally, the adhesion of nanoparticles could have a significant effect on their catalytic performance via a reduction in surface area.

Future work will be directed towards elucidating behavior of silicate and other oxide particles, using some of the same approaches here. It is expected that surface bonding occurs with more defects when cations and anions are present, and also that plastic deformation by the motion and generation of dislocations at the interface occurs less readily than with metals. It is also possible that the approach of determining the oscillation and damping times $\tau$ and $\tau_{d}$ can be extended to computation of the coefficient of restitution when bouncing occurs.

This research was made possible by support from the National Science Foundation, Division of Astronomical Sciences, under grant 1616511. The authors also acknowledge support and computational time provided by the Institute of Simulation and Training and the STOKES computer cluster at the University of Central Florida.

The raw and processed data required to reproduce these findings cannot be shared at this time as the data also forms part of an ongoing study.