Molecular dynamics of nanodroplet impact: The effect of particle resolution in the projectile model

magnetic properties and perpendicular magnetic anisotropy in the Heusler alloy Fe1.5CoGe ABSTRACT The energetic impact of projectiles with diameters between a few nanometers and microns can now be investigated with electrospays operating in the cone-jet mode, a particle source that produces beams of highly charged and monodisperse droplets with average diameters down to a few nanometers. The hypervelocity impact of these nanodroplets on ceramic targets cause sputtering, amorphization and cratering. This experimental phenomenology has been reproduced with molecular dynamics modeling the molecules of the projectile as large pseudo atoms. This model can be over simplistic, especially for liquids made of large molecules, and the goal of this article is to evaluate this uncertainty by comparing the impacts resulting from this coarse model with those of a full atomic model of the molecules. Impact simulations for projectiles of two liquids with dissimilar molecular complexity, formamide and 1-ethyl-3-methylimidazolium bis (trifluoro-methylsulfonyl) imide, show that sufficient resolution of the projectile is needed to reproduce the impact zone, which has a depth of the order of the diameter of the projectile.


I. INTRODUCTION
The energetic impact of electrosprayed nanodroplets modifies the surface of hard materials in different ways. It causes surface amorphization, produces craters comparable to the size of the projectile, and sputters a large number of atoms. 1,2 The energy of the nanoprojectile is transferred to the target within a depth of the order of the diameter of the projectile. Sputtering by atomic-level particles is well understood on the basis of Sigmund's collision cascade model, 3 while the impact by larger projectiles like molecular ions, gas cluster ions and nanodroplets has been studied more recently, often using Molecular Dynamics. Work in the projectile diameter range between a few nanometers and microns has lagged the upper and lower diameter ranges due to the lack of an appropriate particle source. 4 The impact of electrosprayed nanodroplets on silicon and other ceramic targets has been studied both experimentally, 2,5 and through Molecular Dynamics. 6,7 This computational technique integrates the equations of motion of all atoms within a simulation cell using prescribed interatomic potentials to reproduce the interaction forces. The computed atomic positions, velocities and forces are postprocessed to obtain thermodynamic properties as well as other parameters of interest. In the case of nanodroplet impact, the existing simulations model the molecules of droplets as single mass points, with a mass totaling that of the constituent atoms. 6,7 Although an useful approximation, this pseudo-atom model (PA) misses aspects of the impact associated with the actual and finer distribution of the projectile's energy and momentum among all constituent atoms. The main goals of this article are evaluating this shortfall and reproducing the physics of the impact more accurately, by employing a full atomic model (AM) for the molecules that accounts for all atoms and atomic bonds in the projectile. To this end we simulate the impact of a nanodroplet with a diameter of 5 nm on a (100) Si target at varying impact velocity, using two liquids of distinct molecular masses, formamide and 1-ethyl-3methylimidazolium bis (trifluoro-methylsulfonyl) imide (EMI-Im), and the AM and PA models. The phenomenologies of the impacts are compared, with focus on the transfer of energy to the target, the ARTICLE scitation.org/journal/adv depth of the amorphization region, and the sputtering yield. The breakup of atomic bonds and its effect on the transfer of energy is also investigated with the AM model.

II. SIMULATION METHODS
The impact of a nanodroplet on a silicon target is simulated with LAMMPS (Large-scale Atomic/Molecular Massively Parallel Simulator) 8,9 in a 27.15 nm x 27.15 nm x 46.16 nm simulation box. 6 The simulations are run for 100 ps with an integration timestep of 1 fs during the initial 10 ps, and 4 fs from this point forward. 10 The silicon target is a rectangular slab with 755,000 atoms, has a depth of 20.36 nm, and a 27.15 nm x 27.15 nm cross section. A Berendsen thermostat is applied to all sides except for the one receiving the impact, in order to prevent the reflection of shock waves. 11,12 The thickness of the thermostat is 2.17 nm in the lateral sides and 3.8 nm in the bottom one. 10 The thermostat is set at 293 K, and the slab relaxed for 200 ps. The resulting equilibrated configuration is used for all impacts. The atoms are arranged in the standard cubic diamond lattice (100) with a lattice constant of 5.431 Å. The Stillinger Weber potential is used for the interaction between silicon atoms. 13 The nanodroplets have a diameter of 5 nm. Two different liquids are simulated, formamide and EMI-Im. This choice is motivated by the very different molecular masses and complexities of both liquids, which aligns with the stated goal of studying the complexity of the molecular model on the impact; and available experimental data on nanodroplet impact for both liquids. 14 Only liquids like formamide and EMI-Im with very low vapor pressure and high electrical conductivity can be electrosprayed in vacuum to produce the charged nanodroplets required for hypervelocity impact experiments. The number of molecules per droplet is chosen to match the density of the liquid. Table I lists basic parameters of the nanodroplets.
The PA model treats a molecule as a point mass, 9 while the AM model takes into account all atoms and atomic bonds in the molecule. In the AM model the relative positions of the atoms are calculated with the Winmostar software for the formamide molecule, and taken from Ref. 15 for EMI-Im. Each molecule is then brought to equilibrium at 0 K. The molecules are evenly distributed inside a 5 nm diameter sphere using the PACKMOL software. 16 The Ziegler-Biersack-Littmark (ZBL) potential, 17,18 designed to reproduce the repulsive force in high energy atomic collisions, is used to describe the interaction between all atoms not sharing an atomic bond. In the AM model the potential of an atomic bond is constructed with a combination of the ZBL and Morse 19 potentials. The Morse potential, with its energy minimum, models the attractive force and the bond dissociation energy, while the ZBL potential is used to model the close-range repulsive force. Both potentials are smoothly connected in the close-range side of the energy minimum with an interpolating function of the form e B 0 +B 1 r+B 2 r 2 +B 3 r 3 +B 4 r 4 +B 5 r 5 , as shown in Fig. 1. r stands for the distance between the bonded atoms. The six parameters of the interpolating function are determined by requiring continuity of the function and of its first and second derivatives at the endpoints of the connecting interval. The positions of the endpoints is determined by minimizing the length of the interpolated gap. The Morse potential is given by: where r 0 is the bond equilibrium distance and De is the dissociation energy. The parameter a is related to the spring constant ke of the harmonic potential function, a = √ ke/2De. Tables II and  III collect Table IV shows impact parameters. The simulations are done for five projectile kinetic energies, with impact velocities that differ among the two liquids due to their different densities. In many of the plots describing the impact we will express the simulation time in units of the characteristic time tc = Dp/Vp to help collapsing the curves for different impact velocities; Dp and Vp stand for the diameter and the velocity of the projectile respectively.
The local temperature at the position of the j − th silicon atom is defined as the average of the kinetic energy of neighboring silicon atoms within a spherical control volume of radius 10.
N is the total number of neighboring atoms, kB is the Boltzmann constant, m is the mass of the silicon atom, vi is the velocity of the i − th neighbor atom, and ⟨vj⟩ is the average velocity of all atoms in the control volume. There are approximately 280 atoms inside the control volume, a number sufficiently large to obtain reliable statistics. Near the coordination number CN of the silicon atoms, computed with the radial distribution function, is used to determine the silicon phase. The crystalline phase has a coordination number of 4, 4.1 < CN < 4.3 corresponds to the amorphous phase, 4.3 < CN < 4.9 is indicative of a glassy phase, and the liquid phase is characterized by CN > 4.9. 9

ARTICLE
scitation.org/journal/adv a more elevated temperature field than the AM model. The EMI-Im projectile heats up the surrounding silicon atoms more than the formamide projectile when using the PA model, while the differences in the temperature fields for the AM model are much smaller. The peak temperature produced by the PA/EMI-Im projectile exceeds 3500 K, while the peak temperature for the AM model is near 1500 K.
Most of initial kinetic energy of the projectile is transferred to the target during the impact, increasing its thermal E Th , potential EP, and translational ECM energies. ECM is computed with the local average velocity field using the same control volume as in Eq. (2) where the summation is over all silicon atoms in the target. The atoms of the projectile eventually rebound out of the target with a total kinetic energy EK. Figure 5 plots the evolution of these energies for the four atomic-model/liquid combinations, and for an impact energy of 15.31 keV. All energies are given as the ratio over the initial projectile energy EK(0), and in the case of the target they are offset with the values before the impact. Figure 5(d) shows that at these high impact velocities most of the projectile's initial kinetic energy is transferred to the substrate. The transfer is almost 100% efficient for the case PA/EMI-Im, followed by PA/formamide, while the AM model yields a lower energy transfer with similar curves for both liquids (transfer efficiency of 85%). Note also that in all cases the transfer is completed within approximately 2 tc time units. Figure 5(b) shows that most of the energy transferred to the target is in the form of potential energy. Although in the initial stage of the impact, e.g. before t/tc = 2, this potential energy is associated with the compression of the target, the subsequently slow decay is indicative of an alternative storage mechanism such as a phase change (we will show that a significant amount of crystalline silicon melts, and later solidifies as an amorphous phase). Similarly to the EK curves, the conversion into potential energy is largest for PA/EMI-Im, followed by PA/formamide, while the AM model yields the lowest potential energy conversions and similar curves for both liquids. Figure 5(c) shows that, except for the initial stage of the impact in which the surface recedes while accommodating the momentum of the projectile, the translational energy of the substrate is small and becomes negligible after a few oscillations. The oscillations are indicative of the initial penetration of the projectile, and the push back of the projectile by the compressed substrate; subsequent oscillations of the target are rapidly damped. Figure 5(a) shows that the thermal energy of the target reaches values between 16% and 27% of the initial projectile's energy, and follows the same atomic-model/liquid trends as for the potential energy. The plots in Fig. 5 illustrate a clear connection between energy transfer processes and the particle resolution of the projectile: as the number of particles defining the projectile decreases, the energy transfer is more intense producing a hotter and larger number of silicon atoms with an excess potential energy. Furthermore the characteristics of the energy transfer asymptotes relatively fast with particle resolution, as the plots for AM/EMI-Im and AM/formamide are almost identical (the PA/EMI-Im, PA/formamide, AM/EMI-Im and AM/formamide configuration track 153, 989, 5202 and 5934 point masses respectively). Figure 6 shows how the average thermal energy per atom changes with time and distance from the surface of the crater, for an impact energy of 15.31 keV. To compute the average thermal energy we divide the target into consecutive layers 1.1 nm thick surrounding the interface, obtain the total thermal energy of each layer, and divide these values by the numbers of atoms in the layers. The projectile transfers its kinetic energy through collisions with atoms in the first layer, which therefore exhibits the highest thermal energy per atom. The PA/EMI-Im case, with a peak value of 12613 K, produces significantly higher temperatures than the other simulations, especially for t/tc < 2. The other three cases have similar temperature profiles. As in Fig. 5, there is an correspondence between the particle resolution of the projectile and the intensity of the thermal field: the PA/EMI-Im projectile, having its kinetic energy distributed among fewer and heavier point masses, produces a less uniform distribution of atomic collisions characterized by much more energetic knock-on collisions, a higher variance in the velocity of the silicon atoms and ultimately a more intense temperature field. Figure 7 shows cross sections of the impact for the PA/EMI-Im and AM/EMI-Im projectiles, at two different velocities (5 km/s and 7 km/s) and two different times (t/tc = 4 and t/tc = 8). At these impact velocities the atoms of the AM/EMI-Im projectile do not penetrate through the target, producing a sharp and symmetric crater surface. On the other hand the massive pseudo atoms of the PA/EMI-Im projectile do penetrate through the target, and as many as 95% of the pseudo atoms remain implanted in the target by the end of the simulation for the 7 km/s impact case. Figure 8 shows the shape of the crater at 100 ps from the start of the impact, for the four atomic-model/liquid combinations and all impacts velocities. All fields are stationary by this time and the profiles thus show the final shape of the crater. Although the PA/formamide, AM/formamide and AM/EMI-Im impacts always produce symmetric craters with volumes that increase with the impact velocity, there are differences among the profiles of the three projectiles: PA/formamide produces the largest craters, while AM/EMI-Im produces craters with the higher aspect ratio. The craters produced by PA/EMI-Im are much more irregular and shallow. This difference is caused by the implantation of the EMI-Im pseudo atoms which, being more significant at increasing impact velocity, raises the surface. Table V collects the atoms sputtered from the target for all impacts simulated, i.e. for all atomic-model/liquid and projectile kinetic energy combinations. PA/EMI-Im is the only projectile that sputters a large number of atoms, while none of the AM projectiles ejects any. On the other hand experiments show that electrosprayed nanodroplets of EMI-Im and formamide impacting on silicon at the velocities simulated in this article sputter large numbers of atoms. In particular, the experimental sputtering yields are between 0 and 0.1 atoms per molecule for formamide, and between 1 and 2.75 atoms per molecule for EMI-Im. 1,10,18 The inability of the AM projectiles to  sputter atoms is paradoxical because knock-on collisions is the main sputtering mechanism, 10 and therefore the realistic atomic resolution of the AM model should reproduce the sputtering conditions much more accurately than the PA model. The failure of the AM model to reproduce the experimental sputtering suggests that the conditions under which experimental sputtering occurs differ from those in the simulations. The following two differences, or a combination of both, are likely causing the disagreement: first, the use of the bulk Stillinger-Weber potential for the silicon atoms on the surface is inaccurate (these atoms have a lower coordination number, are easier to sputter, and participate in the most energetic knockon collisions with the atoms of the projectile), and using a specific potential for them to reconstruct the surface may be important; 23 second, although the AM simulations show that the energy transferred to silicon atoms at these impact velocities is insufficient to break atomic bonds in the bulk of a perfect silicon crystal, the impact does alter the crystalline structure creating defects and an amorphous phase, and a significant fraction of the projectile's kinetic energy remains stored near the surface of the crater in the form of potential energy. In sputtering experiments any spot on the target receives multiple impacts, and we expect that the increasing weakness of the region affected by the impact leads to conditions in which knock-on collisions are able to sputter target atoms from the modified target, at impact velocities that are not energetic enough to sputter atoms from the perfect crystal. For example, Ref. 1 shows that the surface of a silicon target develops a characteristic topography of large craters only after it receives a critical projectile dose, which supports the idea that the accumulated weakening of the crystal is an important factor in the impact phenomenology. Testing the accumulated weakening hypothesis is difficult and requires work beyond the scope of this article: first, although simulating a single impact on a virgin target is relatively simple, the associated experiment characterizing an isolated impact is far from trivial because of the need to isolate a nanodroplet and measure its diameter and velocity in flight to the target; second, the simulation of the impact on a target modified by previous impacts requires a potential for the silicon atoms other than the standard Stillinger-Weber, which reproduces well the crystalline phase but is not accurate for either the amorphous or the damaged crystalline phases created by the impacts. Finally, we note that pretreatment of the surface of the target (well polished single-crystal wafers) is not important in the outcome of experiments employing a significant projectile dose, 1 because the surface is rapidly altered by the projectiles. Thus the majority of the nanodroplets impact on a surface that is gradually receding due to the sputtering, and is characterized by indentations comparable to their sizes and an atomic arrangement different from that of the original single-crystal (amorphous phase and/or crystalline phase with defects). Figure 9 shows the coordination number of the silicon atoms at 100 ps from the time of impact, for the four atomic-model/liquid combinations and a projectile's kinetic energy of 15.31 keV. Dark blue atoms have a coordination number of 4 and are in the cubic diamond structure. The crater is surrounded by a large number of atoms with coordination numbers between 4.2 and 4.8, i.e. that are in either the amorphous or glassy phases. Both phases result from the fast quenching of molten silicon. The PA/EMI-Im projetile produces the largest quantity of glassy atoms due to the higher peak temperatures and faster cooling rates. The thickness of the noncrystalline region for this projectile is also the largest due to the deeper penetration of its atoms into the slab. The two AM projectiles and the PM/formamide projectile produce similar coordination number fields. Figure 10 shows the number of atoms in the amorphous and glassy phases as a function of the projectile's kinetic energy, for the four projectiles. The amount of atoms that are left in this higher energetic state is proportional to the energy of the projectile, and is similar for all atomic-model/liquid combinations. As previously found in Ref. 6, a threshold velocity of 3 km/s is required to amorphatize silicon. Thus a minimum projectile energy density is needed to increase the temperature of the target over the melting point, and increasing the energy above this value increases proportionally the volume of the target that melts and is amorphatized, regardless of the particle resolution used to model the projectile. In fact, the large majority of atoms in the amorphous/glassy phase, and which are previously in the liquid phase, reach this latter state through collisions with other silicon atoms rather than through collisions with the atoms of the projectile.  Tables VI and VII collect, for all the atomic bonds of the formamide and EMI-Im molecules, the fraction that is dissociated at each impact velocity. A bond is considered dissociated when the atoms are in the attractive side of the potential and the interaction energy is 1% of the dissociation energy De. Only a very small fraction (bellow 0.9%) of the bonds in the formamide projectile dissociate, the highest impact velocity is needed to produce these low dissociation levels, and there is no clear trend between dissociation fraction and De. The bonds of the EMI-Im molecule exhibit higher levels of dissociation, especially for the C-N and C-S bonds, which are the ones with lower energies. The highest dissociation fraction is 20.5% for the bond of lowest De at the highest impact velocity simulated. Note that bonds with higher De are generally less likely to dissociate in the case of EMI-Im. C-N bonds located in the cation ring are the ones with highest dissociation fractions, and are also the bonds with lowest dissociation energy. EMI-Im exhibits different dissociation fractions for some values of the bond dissociation energy. This corresponds to bonds that are located in multiple positions within the molecule. Thus bond dissociation not only depends on De, but also depends on the location of the bond and the neighboring atoms. For example at De = 2.385 eV and an impact velocity of 7 km/s,  1% of the C-N bonds located in the ethyl radical of the cation dissociate while 11.7% of those in the methyl radical break free. Furthermore, there are four bonds with almost identical dissociation energies near 3.903 eV, three of them are C-H located in the methyl and ethyl groups of the cation and one is the C-S bond located in the anion. The C-S bond exhibits a dissociation fraction significantly higher than the C-H bonds. The dissociation energies in Table III have been estimated with values for the same atomic bonds in similar chemical groups. To investigate the effect of this uncertainty in the transfer of energy we model the impact of an EMI-Im projectile with bonds forced to remain intact (dissociation energies are made tenfold those in Table III). Figure 11 plots the evolution of the kinetic and potential energies of the atoms of the projectile, and the thermal and potential energies of the target for both sets of dissociation energies and an impact velocity of 7 km/s. The differences are not significant: the excess potential energy of the projectile atoms modeled with our best-estimate dissociation energies (i.e. Table III) is 22.32% larger than for bonds with tenfold dissociation energies. This difference is just 1.58% of the initial kinetic energy of the projectile. Note also that the two potential energy curves separate in the early stage of the impact, between 1 and 2 tc time units, indicating that this is the time period when bonds dissociate. The energy released by the broken bonds increases very slightly both the kinetic energy of the EMI-Im atoms, and the potential energy of the silicon atoms.

IV. CONCLUSIONS
The particle resolution with which the nanoprojectile is modeled has a significant effect on the temperature field of the target, especially in the region near the surface where the kinetic energy of the projectile is directly transferred through knock-on and secondary atomic collisions. EMI-Im simulated with the pseudoatom model is the projectile with lowest particle resolution (153 point-mass particles), and yields the highest temperature fields. The same model applied to formamide has a particle resolution of 989, which is already sufficiently high to produce temperature fields similar to the full atomic model (5202 atoms for EMI-Im, and 5934 atoms for formamide). When the projectile is simulated with the full atomic model the basic impact phenomenology (energy transfer and temperature fields) of both liquids is almost identical. The results derived from the atomic model, being more realistic than the pseudo-atom model, must reproduce the actual impact more accurately.
The impacts simulated with the atomic model do not produce sputtering in the velocity range 3.00 km/s < Vp < 8.15 km/s, which is at odds with experimental results. On the other hand the pseudoatom model yield sputtering values similar to those in experiments. Since sputtering is caused by knock-on collisions, 10 and the unrealistically high momentum of the EMI-Im pseudo-atom is needed to reproduce experimental sputtering yields, we conclude that sputtering in experiments must be occurring from a target that significantly differs from that of the perfect crystal. The simulations show that the region of the target surrounding the impact is modified to produce a large density of defects and amorphous pockets which store, in the form of potential energy, a large fraction of the energy of the projectile. Sputtering in experiments is likely occurring when projectiles impact on this modified surface. This is consistent with the cratering phenomenology observed in silicon, which depends on the target receiving a threshold projectile dose. 1 In addition reconstruction of the surface with a specific potential (instead of using a bulk potential for all target atoms) will likely lead to better agreement between simulations and experiments.
The atomic model allows the study of atomic bond dissociation. The model is purely mechanical, and does not capture atomic recombination and chemistry. Bond dissociation is minimal in formamide, and reaches values of up to 20.5% for the C-N bonds in the imidazolium ring of EMI-Im, at a projectile velocity of 7 km/s. This is the bond with lowest dissociation energy in the molecule. Bonds with nearly identical dissociation energies for different atomic pairs can have very different dissociation ratios, as exemplified by the S-N, C-S and C-H bonds in EMI-Im: all three have a dissociation energy of 3.903 eV, but only the C-S bond exhibits significant dissociation. The mechanical effect of bond dissociation on the impact is minimal: when dissociation energies are artificially augmented to prevent dissociation, the amount of energy transferred to the target and the allocation in different energy modes remains virtually unchanged.
In conclusion, a full atomic description of the projectile is preferred over the pseudo-atom model. Sufficient particle resolution of the projectile is needed to reproduce realistically the primary knockon collisions in the projectile and in the target, as well as the subsequent collisions between secondary silicon atoms and the rest of the target. Therefore using a full atomic model is especially important to model large molecules like EMI-Im, and generally to replicate the dynamics of the impact in a region with a depth of the order of the diameter of the projectile.

ACKNOWLEDGMENTS
This work has been partially funded by the Balsells Fellowship program. The authors are grateful for the support of Pete J. Balsells.