Molecular Dynamics Simulation of Nanoindentation-induced Mechanical Deformation and Phase Transformation in Monocrystalline Silicon

This work presents the molecular dynamics approach toward mechanical deformation and phase transformation mechanisms of monocrystalline Si(100) subjected to nanoindentation. We demonstrate phase distributions during loading and unloading stages of both spherical and Berkovich nanoindentations. By searching the presence of the fifth neighboring atom within a non-bonding length, Si-III and Si-XII have been successfully distinguished from Si-I. Crystallinity of this mixed-phase was further identified by radial distribution functions.


Introduction
Silicon plays an important role in applications such as semiconductor devices, sensors, mechanical elements, and electronics. Its electronic characteristics have therefore been intensively investigated. Mechanical properties of Si, however, became a research focus only in the past few years owing to the development of the silicon on insulator (SOI) technology and microelectromechanical systems (MEMS), in which Si serves as a substrate. For these applications, deformation mechanisms of Si under nanocontact are essential.
It is well-known that diamond cubic Si (Si-I) undergoes pressure-induced phase transformations during mechanical loading using diamond anvil cell (DAC) or nanoindentation [1][2][3][4][5][6]. The Si-I transforms to the metallic b-Sn (Si-II) phase under a load of up to 11 GPa [1]. Upon pressure release, Si-II undergoes a further phase transformation to a mixed-phase of Si-III (bc8, body-centered-cubic structure) and Si-XII (r8, rhombohedral structure) at a low unloading rate while it transforms to the a-Si phase at a fast unloading rate [3,7,8]. Jang et al. [9] reported the extrusion and phase change mechanism using a sharp or blunt indenter with various indentation loads and rates. Phase transformations corresponding to repeated indentations were also studied by Zarudi et al. [10,11].
Comprehensive understanding of phase transformations in Si requires the use of experimental techniques such as cross-sectional transmission electron microscopy (XTEM), scanning electron microscopy (SEM), and Raman microspectroscopy [5,10]. On the other hand, molecular dynamics (MD) simulations have also been employed to identify the phase transformation mechanism. Among related MD studies, Cheong and Zhang [12] identified different phases through their coordination numbers and also performed the radial distribution function (RDF) analysis. A stress criterion for the onset of the transformation to Si-II was also proposed [13,14].
This study presents the MD approach toward mechanical deformation and phase transformation mechanisms of monocrystalline Si(100) subjected to nanoindentation. The MD simulations were performed to identify load-displacement characteristics of the nanoindentation process and nanoindentation-induced phase transformations during loading and unloading. Both spherical and Berkovich indenters were considered.

Molecular Dynamics Simulation
The interatomic potential function proposed by Tersoff [15][16][17][18] that considers the effect of bond angle and covalent bonds has been shown to be particularly feasible in dealing with IV elements and those with a diamond lattice structure such as carbon, silicon, and germanium. The Tersoff function was therefore adopted in this study to analyze the dynamic correlations in carbon-carbon and silicon-silicon atoms. In regard to the mutual interaction between carbon and silicon under the equivalent potential, we made use of the two-body Morse potential [12], which has been well described for carbon-silicon atoms. Although a two-body potential leads to less precise solutions than a many-body potential does, its parameters can be accurately calibrated by spectrum data, and hence is extensively employed in MD simulations. In addition to the periodic boundary conditions, a modified five-step methodology was used to incorporate Newton's equations of motion so that the position and velocity of a particle can be effectively evaluated. Moreover, the mixed neighbor list was applied to enhance computational efficiency.
Physical models for spherical and Berkovich indenters contained 46,665 and 29,935 carbon atoms, respectively, with covalent bonds. The 250 Å 9 250 Å 9 175 Å modeling region of the (001)-oriented Si substrate contained 518,400 silicon atoms with covalent bonds. We simulated the nanoindentation process by applying perpendicular loading along the (001) direction. Detailed MD modeling and calculation techniques of nanoindentation on monocrystalline Si(100) are referred to Lin et al. [19]. The maximum penetration depth in the present MD simulations was set at 3.5 nm.

Results and Discussion
Since the formation of metastable Si-III and Si-XII phases is strongly stress-dependent, different stress distributions induced by spherical and Berkovich indenters would result in different Si-III and Si-XII distributions within the nanoindentation-induced deformed region. Boyer et al. [20] have observed and discussed the presence of Si-I, Si-II, Si-III, Si-XII, and bct5-Si phases during nanoindentation. Among the several possible mechanisms of phase transformations in Si, it is generally acceptable that Si-I transforms to the metallic Si-II during the loading stage. The Si-I crystalline structure contains four nearest neighbors at a distance of 2.35 Å at ambient pressure. When the stress increases up to 10.3 GPa, Si-I transforms to Si-II, whose crystalline structure contains four nearest neighbors at a distance of 2.42 Å along with two others at 2.57 Å . Moreover, the bct5-Si crystalline structure contains one neighbor at a distance of 2.31 Å and four others at 2.44 Å [21]. The Si-III is constructed by four nearest neighbors within a distance of 2.37 Å and a unique one at 3.41 Å at 2 GPa. The Si-XII is with the four nearest neighbors within a distance of 2.39 Å and also a unique one at a distance of 3.23 or 3.36 Å at 2 GPa [22,23]. Upon pressure release, part of the highly pressured Si-II phase would transform to a mixed-phase of metastable Si-III and Si-XII. Although distinguishing of Si-III and Si-XII from Si-I apparently has been a difficulty in previous MD studies because the coordination numbers of these phases are identical at four, the two metastable phases can be readily identified from Si-I by searching the presence of the fifth neighboring atom within a non-bonding length.
Previous MD simulations showed that under nanoindentation, the bond angle along the (001)-oriented surface direction of monocrystalline Si could be gradually compressed from 90°to 70°, whereas the relative slip among atoms along the compression direction would slowly form Si-II [24]. A pop-in event encountered during the loading stage is an indicator of the occurrence of plastic deformation that leads to phase transformation from Si-I to Si-II in the severely compressed region [19]. Most of the previous studies that explored phase transformations of Si applied a spherical indenter capable of triggering large-scale phase transformations. In the present MD simulations, a spherical indenter was first adopted to interpret phase transformation features in monocrystalline Si. We then adopted a Berkovich indenter in the simulations to compare the difference of phases induced by the two indenters. Figure 1 shows the load-displacement curves led by spherical and Berkovich indenters. At an identical penetration depth, the total deformation energy of the spherical indenter is larger than that of the Berkovich indenter. An apparent pop-out event is also present for the spherical indenter during the unloading stage. However, the pop-out event is unapparent for the Berkovich indenter perhaps because the maximum penetration depth is not large enough in the MD simulations to trigger the event. Figure 2a shows phase distributions on the cross-sectional (011) plane under an indentation load induced by the spherical indenter along (001) at the moment when the maximum penetration depth is reached. Clearly, the highly pressured zone (in red) is surrounded from below by the Si-II phase (in yellow) while the Si-II phase is surrounded by the bct5-Si phase (in cyan). The tilted distributions of these phases follow the {110} slip planes of monocrystalline Si. It is particularly interesting to note that a ring representing a mixed-phase of bct5-Si and Si-I (blank) is present close to the boundary of Si-II. The presence of this mixed-phase implies that energy transfer during nanoindentation is non-continuous, indicating that the continuum assumption is no longer feasible under such a circumstance. Figure 2b shows phase distributions on the crosssectional (011) plane after the spherical indenter is completely withdrawn. Residual phases consist of a mixture of Si-III and Si-XII (in green), Si-II, and the amorphous phase. The presence of Si-III and Si-XII as well as the amorphous phase corresponds to the pop-out event occurred during the unloading stage. Furthermore, recrystallization upon unloading is the most active along the slip planes.
Phase distributions on the cross-sectional (011) plane induced by a Berkovich indenter, as shown in Fig. 3, are in general similar to the ones induced by a spherical indenter, while the phase transformation region of the former is smaller than the latter. A ring surrounding Si-II of a mixedphase of bct5-Si and Si-I is also present.
Crystallinity of Si-III and Si-XII for monocrystalline Si(100) subjected to spherical or Berkovich indentation along the (001) direction was identified by RDF, as shown in Fig. 4. For both indentations, there are obvious peaks at bond lengths of 2.3-2.4 Å , 3 Å , and 3.2-3.45 Å . The first peak corresponds to the fact that the mixed-phase of Si-III and Si-XII is concentrated at 2.37-2.39 Å while the third peak refers to the presence of the fifth neighboring atom of Si-III or Si-XII within a non-bonding length at 3.23-3.41 Å . The second peak at 3 Å should come from the   amorphous phase [25] whose atoms are separated at the critical bond length set in our MD simulations (3 Å ) as a result of atomic interactions between the indenter and Si. We need to emphasize that this particular peak would correspond to a slightly different bond length when a different potential function is followed. Moreover, minor peaks at bond lengths greater than 3 Å can be referred to thermal vibrations of Si atoms [25].

Conclusion
Nanoindentation-induced deformation and phase transformations in monocrystalline Si(100) were investigated through MD simulations. The Si-III and Si-XII were distinguished from Si-I by searching the presence of the fifth neighboring atom within a non-bonding length. Crystallinity of the mixed Si-III and Si-XII phase was further identified by RDF. The MD results also indicate that phase distributions induced by a Berkovich indenter are in general similar to the ones induced by a spherical indenter, while the phase transformation region of the former is smaller than the latter.