Comparative Study of Interatomic Interaction Potentials for Describing Indentation into Si Using Molecular Dynamics Simulation

We compare the performance of three interatomic interaction potentials for describing the evolution of plasticity and phase transformations in Si: the well established Stillinger-Weber potential, a recent modification used in the description of Al/Si composites, and a modification of the well known Tersoff potential. We show that the generation of dislocations and the evolution of plasticity are well described by the Stillinger-Weber potential and its modification, while the phase transformation to the high-pressure bct5 modification and the subsequent amorphization are better included in the modified Tersoff potential.


I. Introduction
Nanomachining and in particular nanoindentation studies are performed increasingly often using molecular dynamics (MD) simulations [1]. These need only the interatomic interaction potentials as input into the simulation; thus the results obtained are essentially as good as the potentials used. The reliability of such potentials has been investigated in many cases, in particular for fcc [2], bcc [3,4], and hcp [5] metals, but also in Si [6,7].
More recently, indentation into composites, such as Al/Si metal/ceramic composites, has raised interest. Here the generation of suitable interatomic potentials and their assessment is more complex than for the elemental materials mentioned above.
For the Al/Si system, an interaction potential has been developed by Saidi et al. [8]. It is of the so-called angular embedded-atom method (AEAM) type and will for brevity be denoted as the AEAM potential. It combines an Al potential developed by Mendelev et al. [9] with a modified form of the renowned Stillinger-Weber (SW) potential for Si [10]; the modification was introduced in order to stabilize the cubic diamond (cd) structure over the wurtzite structure. The AEAM finally adds an interspecies term designed to reproduce the experimental phase diagram and the heat of mixing as a function of composition. It appears tempting to use this AEAM potential for simulating machining of the composite material.
However, before doing so, it is mandatory to study the performance of the potential for simulating the plastic behavior of the pure elements, Al and Si. Since the Al part of the potential is well tested [4,9], we concentrate here on its performance in pure Si. Previous simulation studies of Si nanoindentation [6,11,12] did not compare the performance of different potentials. Recently Goel et al. [13] reviewed MD simulations of the machining of pure Si, but did not discuss the AEAM potential.
In this paper, we explore the performance of the AEAM potential to describe nanoindentation into elemental Si. Our focus is on the evolution of plasticity and of pressure-induced phase transformations in this material. To this end, we compare a simulation using the AEAM potential with that for using the standard form of the SW potential [10]. In addition, we compare to a simulation using a variant of the renowned Tersoff potential [14]; since the latter describes plasticity only poorly, we use here the recently formulated modified Tersoff (MT) potential [15].

II. Method
We use molecular dynamics (MD) simulation to study the plastic deformation and phase transformation of Si upon nanoindentation. Si has a cd structure. We orient the Si crystal to have a free {100} top surface. We denote the direction of indentation, perpendicular to the top surface, by −z; the lateral coordinates x and y are aligned with the ⟨100⟩ crystal axes. The Si crystallite extends 21.7 nm in depth and 38.0 nm laterally. We prepare three Si crystallites, subjected to the AEAM potential [8], the SW potential [10] and the MT potentials [15], respectively.
Before starting the simulations, the crystals are relaxed such that all components of the stress tensor reach values < 12 MPa and the temperature drops below 1 mK. We use periodic boundary conditions in the lateral directions. The bottom layers with a width of 1 nm are fixed in order to prevent any rigid-body motion of the substrate during the indentation. The next 1-nm layers at the bottom as well as the outermost 1-nm layers at the sides of the substrate are thermostatted to 0 K using velocity scaling.
A spherical indenter with a radius of R = 10 nm is used to perform the indentation. According to the recipe of Kelchner et al. [16], the interaction between a Si atom and the indenter is purely repulsive and follows the potential where r is the distance between the atom and the center of the indenter; for r > R, the force vanishes. The force constant is set to k = 11.7 eV/ Å3 [7]. Upon loading, the indenter moves perpendicularly down into the surface to a final depth of d = 8 nm with a velocity of 20 m/s. The calculations are performed using the open-source code LAMMPS [17]. The Crystal Analysis Tool (CAT) is used to identify dislocations and stacking faults [18][19][20]. Phase transformations are analyzed by calculating the coordination numbers of atoms; these are determined using a cutoff distance of 3.0 Å.

III. Results
Fig. 1 displays the evolution of the force and the contact pressure with indentation depth. The results for the three potentials are quite close to each other. In the elastic regime, below 6-8 Å, the results coincide, as they must since all potentials have been fitted to the elastic constants of Si. In the plastic regime, the curves are characterized by a sequence of peaks and pop-in events which are due to both dislocation generation and phase transformations occurring under the indenter [7]. Note that the contact pressure does not yet saturate, since the material hardness of around 12-13 GPa [21] has just been reached at the end of the indentation. The amount of defect formation and phase transformation occurring in the materials can be read off Fig. 2a, which shows the fraction of remaining cd phase; here again all potentials agree with each other. This primarily demonstrates that the volume which is affected by the high pressure exerted by the indenter nicely coincides for the three potentials studied.
The amount of atoms in stacking-fault planes generated and the length of dislocations formed in the cd phase are shown in Figs. 2b and c. Here strong differences are seen. In particular the MT potential differs from the two other potentials in that the amount of defects in the cd phase is by a factor of 5-10 smaller than in the other potentials. This finding is in agreement with earlier results that the Tersoff potential predicts little plasticity [22]; also Kim and Oh [11] reported that they found no dislocation nucleation upon Si indentation. Pastewka et al. [23] actually state that the Tersoff potential has a tendency to describe Si as too ductile, since the narrow potential cut-off impedes bond breaking. On the other hand, the predictions of the SW and the AEAM potential are rather close to each other; in fact the AEAM even predicts more defects than the SW potential from which it was modified. This may be connected to the modification of the pair interaction term in the modified SW potential that exhibits a shallow minimum at the position of the third neighbor in the cd structure [8]. Godet et al. [24] compare the response of Si to shear of the {111} planes along the ⟨110⟩ direction for SW, Tersoff and ab-initio Si and find that SW better reproduces the ab-initio results. Finally, in their simulation study of the yielding of compressed Si nanospheres, Hale et al. [22] reported earlier that the SW potential favors dislocation plasticity over phase transformationin contrast to the Tersoff potential -, since on the hand the generalized stacking-fault energies are better represented in SW [24] while this potential on the other hand overestimates the energy barriers to high-pressure phases.
The portfolio of distinct types of dislocations generated in Si is shown in Fig. 2d. We observe mainly Burgers vectors ½ ⟨110⟩ (perfect dislocation) and 1/6 ⟨112⟩ (partial dislocation); they are produced in similar amounts in the SW and in the AEAM potential, but considerably less in the MT potential, as discussed above. The snapshots shown in Fig. 3 give a graphical view of the dislocations formed in the cd phase. Fig. 4 allows us to discuss the phase transformation behavior under indentation. The main high-pressure phase that forms is five-fold coordinated; it has a body-centered tetragonal structure and has been termed bct5 [25]. We see in all three potentials evidence that it nucleates first, see the snapshots at 25 Å indentation depth. Under higher pressure, amorphization sets in which is characterized in Fig. 4 by the occurrence of atoms with higher coordination. The bct5 phase remains stable and survives in small crystallites or platelets even at 75 Å indentation depth for the MT potential, while it is largely destroyed for the SW and AEAM potential; note the irregular boundaries of the 5-fold coordinated zones that indicate that these atoms do not form regular crystals. We conclude that the phase transformation to bct5 is not well described in the SW and AEAM potentials.

IV. Summary
The atomistic modeling of indentation into Si shows a clear dependence on the interatomic potential used. The force and contact pressure agree for the three potentials investigated, and so does the atomic fraction of surviving cd phase; this means that the volume affected by phase transformation is similar in all systems studied.
However, the three potentials show clear differences in the development of plasticity and phase transformation. The SW potential and its modification used in the AEAM potential feature a larger number of stacking faults and dislocations formed; the latter potential is particularly abundant in defect formation. The MT potential, on the other hand, generates almost a factor of 5 less dislocations and even less stacking faults. We note that the MT potential has been formulated in order to improve over the even worse performance of the original Tersoff potential [14] to describe plasticity generation [15]; thus the improvement of the MT appears to be only minor.

6
Physical Modeling for Virtual Manufacturing Systems and Processes All potentials show a phase transformation of the virgin cd structure to the high-pressure bct5 modification, followed by amorphization. The bct5 phase is most stable for the MT potential, where pockets of the new phase survive within amorphized material. This feature demonstrates that the well-known capability of the Tersoff potential to describe high-pressure modifications of Si carries well over to its modification, the MT potential.
We conclude that the Al/Si potential can be considered a promising potential to describe machining of Al/Si composites. The Al part of the potential has been tested elsewhere [4]. The interspecies part of the potential has been optimized to describe the enthalpy of mixing for different compositions, the ground state energies of the L12 and B1 structures as well as the eutectic point and other features of the Al/Si phase diagram [8]. It describes Si quite similarly as the SW potential from which it was derived and predicts ample dislocation formation, while improving on the description on the high-temperature cd and wurtzite phase. Only in the description of phase transformations, the generation of high-pressure phases may be underestimated.