Molecular dynamic simulations of the elastic and inelastic surface scattering of nanoparticles.

Molecular dynamics calculations have been undertaken to simulate the collision of a solid, rotating nanoparticle with a planar, two-dimensional surface at thermal velocities (linear and rotational) equivalent to 500 K. During the course of a collision, mechanisms have been introduced into the simulation process that allows for the dissipation of kinetic energy and for components of linear and angular velocity to couple. Although previous studies of particle-particle collisions have used a similar energy dissipation procedure, in these first calculations on particle-surface collisions, it is found that the mechanism actually facilitates the movement of particles across a surface. It is also shown that the direction of travel of particles on a surface is strongly influenced by their rotational motion.


Introduction
There are many examples in chemistry, physics, and biology where collisions between small particles and solid surfaces have important implications for the deposition and accumulation of nano-scale material. 1 Within this category are such diverse processes as the impact of inhaled aerosols on to the surface of the lung, 2 powder coating, 3 and the preparation of food. 4 Possibly one of the most significant developments in this general area is the emergence of methods that are capable of assembling surface layers of nanoparticles. At the centre of a number of these techniques is an electrostatic mechanism that either promotes particle formation or provides particle guidance during the self-assembly process. [5][6][7][8][9] However, what many of these processes lack is an understanding of the events that can take place when single particles collide with a surface, and how certain factors, such as size and impact velocity, might influence the trajectories particles take towards and across a surface to their final location. To date, most theoretical models depicting the collision of a nanoparticle or a large molecule, such as a peptide, with a solid surface have taken the form of finite collections of atoms held together with a suitable interatomic potential, colliding with a surface that it is frequently represented by an appropriate Lennard-Jones or van der Waals potential. [10][11][12][13][14][15][16][17][18][19][20][21] Molecular dynamics (MD) methods are then used to follow the trajectories of these particles as a function of time. Advantages of the atomic-scale approach are firstly, that it provides a realistic description of energy exchange between the centre of mass kinetic energy and the internal degrees of freedom of a particle during a collision with a surface.
Secondly, the individual atoms respond to any forces acting to deform the overall structure of a nanoparticle. The disadvantages are that only comparatively small particles (~ 10 nm) can be treated and the simulation timescale (~ 60 ps) is such that calculations are typically limited to a single particlesurface encounter. These limitations arise because of the time taken to integrate the equations of motion for large numbers of individual atoms, with the size of the integration time step being determined by the fastest motion, which in these cases are interatomic vibrations. However, even from such short simulations, atomic-scale calculations have characterized conditions under which particles can bounce, 17,18,20,21 and can become distorted and/or are deposited on a surface. [16][17][18]20,21 Experimental examples of particle deposition are quite common; less frequent are experimental observations of particles bouncing. 21,22 Reported here are the results of a series of molecular dynamics calculations where the trajectories of single solid particles impacting on a planar solid surface have been followed up to a point where the particles come to rest. There are two significant advantages of this approach; first, there are no restrictions on particle size and secondly, the timescales over which trajectories can be followed are very much longer than is possible with an atomistic MD simulation. In some of the simulations reported below, particle trajectories have been followed for up to 0.2 ms. The calculations have been designed to mimic the characteristics expected of a macroscopic, rotating object striking a solid surface, and it will be seen from the results that many of the behavior patterns exhibited by particles arise because of their rotational motion. A disadvantage of this "solid particle" treatment is that the trajectories do not have the same capability as collections of atoms for absorbing kinetic energy and/or fragmenting; however, many of the latter studies are undertaken at much higher particle kinetic energies than are considered here. 19 Allowance for energy loss by a solid particle comes in the form of a coefficient of restitution, which provides a damping mechanism that removes kinetic energy during the time that the particle is in close proximity with a surface. It is shown that the damping mechanism together with overall rotation facilitate the movement of particles across a surface. In a related study, Krinke et al. have used a stochastic approach to model the distribution patterns of charged particles that have been deposited on to a surface under the influence of an electric field. 5 For larger, micron-sized particles, several kinematic models have been presented in the literature to account for either the rebound of elastic spheres from hard surfaces, 23 or adapted to account for particle adhesion. 24

Molecular dynamic method
Interaction potential between a particle and a surface.
For the purposes of these calculations, the interaction between a particle and a surface is taken to be "soft" and is defined by a Lennard-Jones (9,3) potential of the form: 25,26 where z is the distance from the centre of the particle to the surface,  is the depth of the attractive well between the particle and the surface, and the distance from the surface at which the potential energy is a minimum is given by 3 1 6 = . Two further quantities are derived from Eq. (1), these are the force, FN, which is given by: and a force constant, = 2 2 , which is given by = 27/ 2 . The latter is used to estimate the duration of a collision tcoll = √ /(2 ) , [27][28][29] which in turn is used to provide a measure of frictional behaviour (see below). M is the mass of the particle; this is determined from the volume of the particle and the density of the material from which it is composed. If the interaction potential, ( ), corresponded to a harmonic oscillator, then tcoll would be equivalent to the period of a single oscillation. Table 1 summarises the parameters used in the calculations described below.

Particle rotation and quaternion parameters.
To develop a realistic molecular dynamics model of particlesurface collisions, the equations of motion have been formulated to include rotation of the particle. In order to avoid the problems associated with describing three-dimensional rotational motion in terms of Euler angles, quaternions have been adopted as a system of coordinates for following the time-dependent reorientation of a sphere that occurs as the result of a collision. [30][31][32][33][34][35] Quaternions provide equations for rotational motion that are free from singularities and their time-dependence, together with all necessary transformations between space-fixed and bodyfix axes, can be expressed in terms of matrix algebra. The four quaternion parameters, , , , and , can be written in terms of the three Euler angles as follows:  = cos(/2).cos((+)/2);  = sin(/2).cos(-)/2);  = sin(/2).sin((-)/2);  = cos(/2).sin((+)/2). With these parameters a rotation matrix can be defined as: Matrix A(q) is used to rotate the three space-fixed axes to be parallel to that of the body fixed frame. Since A(q) is orthogonal when q 2 = 1 and A -1 (q) = A T (q), rotation back to a spacefixed frame can be executed through multiplication by the transpose of A(q). For the matrix A(q) to achieve this condition, the quaternions are required to satisfy the equality q 2 =  2 + 2 + 2 + 2 = 1. It is customary for q to be continuously renormalized during the course of a molecular dynamics simulation. [31][32][33][34][35]

Equations of motion
To follow the time-dependent motion of a particle as it strikes a surface, a finite difference method in the form of the velocity Verlet method has been used to generate trajectories. 36 For the quaternions, their time derivatives can be expressed in terms of the principal angular velocity components, px, py and pz as follows: 35 An incremental change in q from t to t+t is then given by q(t+t) = q(t) + ̇(t)t + ̈(t)  , T is the torque being exerted on a particle during a collision and I is the inertia tensor. As the examples being treated are solid spheres, the scalar moment of inertia is given by I = 2/5MR 2 . The torque is given by Where F(t) is the force acting on a particle at time t, and is given by Eq. (2), and R is the radius of the particle. At t+t a new angular velocity is calculated from Where ̇( +  ) = ( + ) . Changes in the velocity, v, and position, r, of the centre of mass have been calculated as follows where a is an acceleration and is given by: A gravitational force was also included in the above equations, but for nanoscale particles this had no significant influence on their dynamics.

Selection of initial conditions
Initial values were assigned to each of the variables of position, linear velocity, angular velocity, and centre-of-mass orientation and these determined the behaviour of a particle during the course of each trajectory. The surface is located at z = 1500 nm and particles with a 200 nm diameter start at x = 0, y= 0 and z = 0. The Box-Muller method 37 has been used to generate pairs of random number, ux and uy, from a standard normal distribution and these were then used to assign the velocities vx = ( / ) 1/2 ux and vy = ( / ) 1/2 uy. To ensure each particle hits the surface, vz was always set equal to +(2 / ) 1

Results and Discussion
Coupling forces acting on a colliding particle For a freely rotating particle falling in the z direction on to a surface defined by x and y, there are opportunities for the exchange of momentum between rotational and translational motion. As a first step, three elastic forces that could serve to couple these different motions are identified and these are forces which can be closely associated with the patterns of behaviour expected of large macroscopic bodies, such as tennis or billiard balls. The term elastic is used here to denote the fact that the forces serve to exchanged energy and momentum between linear and angular motion, but these quantities remain conserved. The latter condition is an important constraint since it ensures that the mechanisms used to describe the exchange of rotational and translational motion have been correctly formulated. [39][40][41] Once these elastic forces have been characterized, a more realistic model is presented below whereby a dissipative frictional correction is added to each elastic force.
What follows has been developed from the very extensive formalism presented by Kondic to describe the motion of spherical particles on a two-dimensional surface. 28 In this first instance, simulations are limited to a single particle on a planar surface; subsequent publications will address the topics of corrugated surfaces and many-particle systems. 42 We begin by defining a force, FN = − ( )̂ where ̂ = z/z and z = z, which acts normal to the surface. For a particle falling from the direction z, FN couples the angular momentum components x and y with the velocity components vx and vy, respectively. The acceleration/deceleration in i = x and/or y due to this component takes the form: ai = where N is a parameter that reflects the degree of coupling that is promoted each time a particle hits a surface and R is the radius of the particle. tracked as a function of time during a collision between a particle and a surface. A rotating particle has been drop vertically onto the surface, but with no initial vx or vy velocity components. The two vertical lines in Fig. 1a represent the time window calculated for tcoll from the equation given above, and as can be seen, it provides a good estimate of the period over which there is a strong interaction between the particle and the surface. In selecting a value for N, the objective is for the coupling to be effective during the time a particle is in close proximity to a surface. The close match between tcoll and the duration over which changes in vx and x take place in Fig. 1a, would suggest that a value of N ~ 0.0005 is appropriate. Figure 2 shows the consequences of this type of collision in terms of changes in direction that a particle can undergo. As the collisions at this stage are elastic, conservation of momentum and energy mean that the particle must rebound from the surface. A macroscopic analogy to Fig. 2 would be the path taken by a spinning tennis ball that is dropped vertically onto a solid surface. 40 A second type of interaction between a particle and a surface is expressed in terms of a shear coupling, where particles approach the surface from a non-normal direction.
Components of such motion are tangential to the surface and can lead to an increase in angular momentum at the expense of linear momentum or vice versa if the incoming particle has, for example, considerable top spin. The torque generated by the collision is given by T =  (FNRS) and the corresponding transfer of linear momentum is again limited by the initial kinetic energy of a particle. S is a parameter that determines the magnitude of the coupling.
 denotes a unit vector that takes the sign of either the x or y component of linear velocity.
The corresponding change in linear momentum is calculated from the acceleration/deceleration: a = - (FNR 2 S /I). Figures 3a and 3b show this type of coupling in action where changes in vx and x are plotted as a function of time during a collision between a particle and a surface. In this example, the particle strikes the surface with an initial angular velocity of zero and S = 0.0005. As already noted, the coupling is elastic in the sense that momentum and total energy are conserved during the collisions and there is, at this stage, no frictional energy loss to the surface.
The final coupling to be considered as acting on a particle is concerned with rotational motion across a surface. Whereas the other two forces are limited by the total energy available, with rotational motion, there is an exchange of momentum between forward and angular velocity such that velocity of the centre of mass tangential to the surface has ultimately to be equal to the angular velocity of a particle at the point where it touches the surface. The acceleration/deceleration of a particle is taken as: ai = -RFNR 2 /I, where  = i -Ri and i is either x or y. The corresponding torque is then given by Ti = RFNR, where R, which in this case has the units of (velocity) -1 , determines the fraction of  that is partitioned between i and i in order for them to achieve a balance. Depending on the magnitudes of i and i,  can be either positive or negative. The duration of a single elastic collision of the type described above, is not sufficient for a particle to acquire rolling motion; each interaction with the surface transfers some fraction of momentum between translation and rotation and, depending on the direction of transfer, those single collisions produce dependences that look very similar to what has been plotted in Figs. 1 and 3. It is recognised that an accurate description of rolling motion is not easily achieved 28 and that there are many alternatives to the formalism given above. 43

Dissipative forces
In order to provide a more realistic description of the dynamics between an impacting particle and a surface, it is necessary to amend the above model by introducing a mechanism whereby particles could come to rest rather than undergoing elastic scattering. To achieve this objective a dissipative correction is introduced to moderate FN and this takes the form: [27][28][29] where = 2(1 − )/ is a coefficient of restitution, which provides a measure of kinetic energy loss by a particle to the surface during a collision, and ̂ = , , or ̂. When used in association with Eq. 1, has the form of a damped oscillator, which provides a mechanical connection between a particle and a surface, and where the damping force is proportional to the velocity of the particle. 44,45 In this instance, damping operates in conjunction with x, y, and z to moderate the motion of a particle across a surface. For the calculations presented here, varies between 0.95 and 0.99, and as will be shown, the ease with which particles come to rest on the surface is extremely sensitive to the assigned value. At present, this dissipative mechanism is only effective for the linear velocity of particles, and there is no direct equivalent for angular velocity; however, the latter can dissipate via the processes which couple individual components of  with . Figure 4a shows the consequences of introducing into the calculation of particle trajectories. In this example, the particle has no initial rotational momentum, the parameters N, S, and R have been set to zero, x and y have been set to +( / ) 1/2 with respect to the laboratory-frame and = 0.98. Since the interaction potential between the particle and the surface is soft, the particle is seen to bounce across the surface whilst gradually losing kinetic energy. The simulation, which lasted 200 s, finishes before the particle comes to a complete rest; however, what is most obvious is that the presence of the tangential velocity components x and y in association with the dissipative force, causes the particle to traverse across the surface. Following from the initial collision, the particle travels for at least a further 1000 nm whilst it loses kinetic energy. An alternative perspective on the energy dissipation processes in Fig.4a is shown in Fig. 4b, where changes in z have been plotted as a function of z, the distance the particle resides above the surface during the course of a trajectory. As can be seen, in the early stages of the dissipation process the particle undergoes large changes in both z and z; however, as the damping process proceeds, factors that contribute to diminish rapidly and the changes in z and z become smaller and smaller. In part, the declining efficiency of the dissipation process (for this set of initial conditions) contributes to the ability of the particle to move across the surface. Note also from Fig. 4b that the minimum distance between the particle and the planar surface located at 1500 nm is approximately 200 nm. For the parameter set given in Table 1 for the Lennard-Jones potential (Eq. 1), the repulsive wall begins to influence the path of a trajectory at around 1300 nm. Figure 5 shows the consequences of introducing particle rotation into a trajectory. For the purpose of reinforcing the outcome, each of the principal components of the angular velocity has been given a value pi = +(3 / ) 1/2 ; however, an angular velocity of this magnitude is not beyond the bounds of a thermal distribution at temperature T. Since a random procedure (see above) is used to assign the quaternion parameters, particles in the laboratory-frame could be rotating either clockwise or counter-clockwise. Comparing Fig. 4a and 5, (x and y are equal to +( / ) 1/2 in both cases) it can be seen that the effects of rotation (for this particular example) are to reduce particle travel in the x direction, but increase it in the y direction. In this example the x component of the laboratory-frame angular velocity, x, has a negative value, which means that the particle starts the trajectory rotating counter-clockwise to the direction of travel as determined byx. Through the dissipative coupling mechanism, this difference in sign has the effect of eventually causing the particle to change direction. Rotational excitation in particles has also been found to moderate surface behavior in atomic-scale simulations. 20 Figure 6 shows how x and y (Fig.6a) and x and y  Figs. 1 and 3. To achieve a change in direction requires the particle to experience a sequence of collisions with the surface, and this eventually brings about a change in the sign of one or more of the velocity components (x in Fig. 6b). Figure 7 shows an extreme example of the significance of rotation to particle trajectories; x has been given the equivalent of a high backspin which then serves to prevent the particle from bouncing forward even though x has an initial positive value. Similar behaviour is to be seen with rotating tennis and billiard balls. 40,41 The effects of varying the coefficient of restitution can be seen in Fig. 8, where simulations have been run with  in the range 0.95 -0.99. With a value of 0.95, the particle rapidly comes to a halt and moves just one diameter after the initial collision. In contrast, when   0.985 the dissipation of kinetic energy is far less effective and there is extensive movement of the particle across the surface. For   0.99 the collision becomes elastic and the particle rebounds away from the surface. Since the trajectory started at x = y = z = 0, the particle has bounced completely free of the surface by the end of the simulation. In figure 8a the red line represents a linear extrapolation of the initial trajectory of the particle, and as can be seen, the actual path taken shows the particle being drawn towards the surface as it makes contact. The outcome of the trajectory calculated for = 0.982 (Fig. 8b) is not too dissimilar to what is seen in Fig. 6 and where, as a consequence of rotational motion, the particle has changed direction as it travels across the surface. In this particular trajectory, it is y that has a negative value. Figure 9 shows an expand section of the final steps in the latter simulation where it can again be seen that the diminishing efficiency of the dissipation mechanism leads to the particle to travel ~ 600 nm in both the x and y directions before the trajectory stops. In total, the trajectories show that particles are capable of travelling distances of the order of a few m following their initial contact with a surface.
As far as the particles are concerned, the above calculations have been undertaken on what could be regarded as a thermal system at 500 K. In order to explore the effects of particle temperature (and, as a result, velocity) on sticking probability, a series of trajectories have been undertaken where the temperature has been varied between 200 K and 5000 K and for values for  of between 0.950 and 0.985. Fig. 10a shows examples of how trajectories run with  = 0.985 evolve as a function of time at two different temperatures. The trajectory at 200 K results in the particle adhering to the surface, whereas at 3000 K the particle escapes after a single collision. Note that the initial velocity of each particle is positive with respect to the surface and then becomes negative on the rebound. As each particle approaches the surface, the dissipative and/or coupling mechanisms lead to an initial reduction in velocity, but the latter then increase as the particle encounters the most attractive section of the energy surface. Similar trajectories for particles consisting of collections of atoms are to be seen in Ref. 18 and 20, and with reference to Eq. 3, atomicscale calculations have also characterized the pattern of behavior seen at 200 K as being that of an under-damped oscillator, 16 which is exactly what dissipative dynamics is designed to achieve. Fig. 10b summarises results for a series of trajectories run under the conditions identified above. In each case a coefficient of restitution that is an alternative to that defined for Eq. 3 has been calculated according to the  Fig. 10a. The results are summarized in Fig. 10b where ez for different values of  has been plotted as a function of the velocity a particle has at the very start of a trajectory (vint). For each value of , the velocity at which particles moves from being trapped at the first collision to rebounding and escaping has been identified and these results are bounded by grey shading. Regions between the four fixed values of  were determined by interpolation. As expected, a reduction in  has to be accompanied by an increase in velocity (higher temperature) if a particle is to rebound and escape from a surface. For   0.95, all particles approaching at a temperature of 5000 K or less can expect to become trapped; however, this conclusion neglects any thermal contribution a solid material may make to particles being released from the surface. At each value of , the gradual increase in -vfz/viz as a function of vint matches the behaviour seen for a non-or low-deformable solid particle; 16,18,21 what the current trajectories will not reproduce is the deformation seen in particles that strike surfaces at higher velocities than those discussed here. [16][17][18]21 An interesting feature of Fig. 10b, is that where a particle moves from being trapped to escaping, the crossing point on the shaded boundary corresponds to   -vfz/viz .

Conclusion
The molecular dynamics method has been used to follow the trajectories of single, rotating particles as they interact with a solid, planar surface. In order to model the possible accommodation of particles by the surface, a dissipative element that depends on the velocity of a particle has been introduced into the calculations. The results show that the latter relationship contributes to the ease with which particles move across a surface, but that exact pathways are strongly influenced by the presence of rotational motion. Certain patterns of behaviour identified here have also been observed in molecular dynamics simulations involving nanoparticles represented by large collections of atoms. The current approach would appear to be most suited to non-or low-deformable particles impacting on surfaces at low velocities, with the advantage that there are no restrictions on particle size or the requirement that particles adhere after a single collision.  used to illustrate operation of the mechanism which couples  and  through N.

Figure 2
As for Fig. 1, but illustrating the changes in direction that come about when a rotating particle undergoes an elastic collision with a plane surface.  Data taken from Fig. 4a showing how the velocity z and the distance of closest approach to the surface, z, change during the course of a dissipative trajectory.

Figure 5
Similar trajectory to that depicted in Fig. 4, but the particle now has rotational velocity and the coupling terms, N, S, and R all have a value of 0.0005.

Figure 6
Changes in linear and angular velocity that accompany the trajectory depicted in Fig. 5 and plotted as a function of time. (a) Changes in laboratory-frame angular velocity; (b) Changes in laboratory-frame linear velocity. Note how coupling to x causes x to become negative, which in turn makes the particle change direction as it moves across the surface.

Figure 7
Example of an elastic trajectory where the particle has been assigned a high angular velocity in the form of backspin (-x).  Expanded view of the trajectory shown in Fig. 8b showing how the diminishing efficiency of the dissipation process leads to the particle moving a significant distance across the surface. In this example, α = 0.982 and the remaining parameters have the values listed in Table 1.  Fig. 10b; (b) Plot of a coefficient of restitution, defined as ez =vfz/viz, for different values of  against the velocity a particle has at the start of a trajectory, vint. The shaded area corresponds to a region, inside of which, particles remain bound once they strike the surface.