Molecular Dynamics Simulation of the Crystal Orientation and Temperature Influences in the Hardness on Monocrystalline Silicon

A nanoindentation simulation using molecular dynamic (MD) method was carried out to investigate the hardness behavior of monocrystalline silicon with a spherical diamond indenter. In this study, Tersoff potential was used to model the interaction of silicon atoms in the specimen, and Morse potential was used to model the interaction between silicon atoms in the specimen and carbon atoms in the indenter. Simulation results indicate that the silicon in the indentation zone undergoes phase transformation from diamond cubic structure to body-centred tetragonal and amorphous structure upon loading of the diamond indenter. After the unloading of the indenter, the crystal lattice reconstructs, and the indented surface with a residual dimple forms due to unrecoverable plastic deformation. Comparison of the hardness of three different crystal surfaces of monocrystalline silicon shows that the (0 0 1) surface behaves the hardest, and the (1 1 1) surface behaves the softest. As for the influence of the indentation temperature, simulation results show that the silicon material softens and adhesiveness of silicon increases at higher indentation temperatures.


Introduction
Nanoindentation is one of most effective methods to measure mechanical characteristics of materials from microscale to nanoscale, such as Young's modulus, hardness, and creep performance [1][2][3][4].As the dimension approaches the nanoscale, the conventional continuum theory often fails in providing reasonable explanations for its mechanical responses.It is very difficult to experimentally monitor the deformation behavior and internal microstructure changes of the specimen's materials at atomic scale via conventional experimental devices.Since MD simulation method can trace the atomic structure behavior in the whole simulation process, it provides a powerful and effective approach to investigate atomic deformation behavior near the material's interfaces [5,6].
Nowadays, many studies of orientation-dependent deformation behaviors have also been carried out to investigate the anisotropy performance of hardness with different crystal orientations.Chen et al. analyzed the effect of crystal orientation on the extrusion formation of silicon surface [7].Results showed that the patterns on the surfaces of (0 0 1), (0 1 1), and (1 1 1) oriented c-Si appear to have four-, two-, and threefold symmetries, respectively.Kim and Oh simulated nanoindentation on the (1 0 0) and (1 1 1) silicon surfaces [8].It is found that the atoms of the amorphous structure take place during large lateral displacement when silicon was loaded on the (1 1 1) silicon surface.Lin et al. investigated the effect of crystal orientation on phase transformation [9].They have found that indentation on the (1 1 0) surface shows more significant internal slipping and spreading of phase transformation than that on the (0 0 1) surface.
Moreover, Goel et al. proposed that temperature gradient between the tool rake and flank face caused relatively higher flank wear [10].Fang et al. used MD to study the effects of temperature on atomic-scale nanoindentation process, and they found that both Young's modulus and hardness became smaller with temperature increase [11].Liu et al. conducted MD simulation to understand the atomistic mechanism of nanoindentation process under different loads, temperatures, and loading rates, and the simulation results indicated that both Young's modulus and hardness became lower as temperature increased [12].
The influence of crystal orientation and temperature is important to the hardness of monocrystalline silicon in the nanoindentation.Compared with the achievement in MD studies of nanoindentation simulation of monocrystalline material, a clear understanding of the influence of crystal orientation and temperature on the hardness of monocrystalline silicon is still lacking, particularly at the atomic scale.Therefore, the main objective of this study was to evaluate the influence of crystal orientation and temperature on the nanoindentation of monocrystalline silicon by using the molecular dynamics (MD).

MD Simulation
2.1.Simulation Model. Figure 1 shows the physical model of nanoindentation on monocrystalline silicon using a spherical diamond indenter.The spherical indenter consists of 7432 atoms with a radius of 21.4 Å, and it is treated as a rigid body due to the evidently higher hardness of diamond than the specimen of silicon.The size of the specimen is 108.6 Å × 108.6 Å × 81.45 Å along the , , and  directions, consisting of 49,420 silicon atoms.
The silicon atoms in the specimen are categorized into three kinds of atoms, namely, Newtonian atoms, Thermostat atoms, and Boundary atoms, respectively.Boundary atoms are fixed in space to reduce the edge effects and maintain the proper symmetry of the lattice.Thermostat atoms are used to ensure reasonable outward heat conduction away from the machined zone.The motion of Thermostat atoms and Newtonian atoms is governed by Newton's second law.In order to reduce the effect of simulation scale, periodic boundary condition is set along  and  direction in MD simulation model [13].The initial temperature is set at 296 K.At this temperature, the silicon atoms are organized in a perfect diamond cubic structure, and the lattice constant is 5.430 Å.Moreover, the indenter is loaded and unloaded at a constant velocity of 40 m/s to reduce the computational time and memory requirements [14].
The MD simulations are conducted with the Large-scale Atomic/Molecular Massively Parallel Simulator (LAMMPS) developed by Plimpton [15].

Selected Potential Energy Function.
The interaction potential function governs the accuracy of a MD simulation which in turn defines the reliability of simulation results.As the Tersoff potential takes the effects of bond angle and covalent bonds into account and has been shown to be particularly feasible in describing those with a diamond lattice structure such as carbon, silicon, and germanium, therefore, it is applied in the current simulation to describe the interaction between silicon atoms in the specimen [16][17][18].Tersoff potential can be expressed as follows: where  and   are, respectively, the total energy and the bond energy of all the atomic bonds; , ,  label the atoms of the system;   is the length of the  bond;   is the bond order term;   is the bond angle between the bonds  and ;   represents a repulsive pair potential;   represents an attractive pair potential;   merely represents a smooth cutoff function to limit the range of the potential;   counts  the number of other bonds to atom  beside the  bond.The parameters , , , , , , , , , , and ℎ are constants adopted from [19] (Table 1).
To avoid the influence of equivalent potential functions, the two-body Morse potential is adopted to simulate the  interaction between the carbon atoms in diamond indenter and silicon atoms in the specimen.Morse potential can be expressed as follows: where (  ) is a pair-potential function and , , and  0 correspond to the cohesion energy, the elastic modulus, and the atomic distance at equilibrium, respectively.The values of potential parameters adopted from reference are  = 0.435 eV,  = 46.487nm −1 , and  0 = 0.19475 nm [20,21].

Analysis of Phase Transformation in Deformed Region.
Figure 2 shows the phase transformation after indentation on the (0 0 1) silicon surface.The coordination numbers (CN) in the deformed zone just underneath the indenter are marked by different colors.It is observed that with indentation depth increasing, -silicon (CN = 6) structure mixed with amorphous phase appears.Thus, it can be known that the specimen's materials undergo plastic deformation and phase transformation from diamond cubic -silicon to -silicon and amorphous silicon.Upon the unloading of the indenter, the crystal lattice reconstructs and finally part of -silicon gradually disappears, transforming into original silicon or amorphous silicon.The results are consistent with the results from the previous paper [1].After the unloading of the indenter, the residual amorphous layers are left on the specimen's surface due to previous plastic deformation, which is unable to be restored.

Influence of Crystal Orientations on Hardness in Monocrystalline Silicon.
In the present work, three simulation cases are considered with different crystal orientations to investigate their influence on hardness.The evolution of phase transformation and Von Mises stress with different crystal orientations can be a useful criterion for identifying the appropriate crystal orientation for practical purposes.In our simulation, one region of the specimen underneath the indenter is selected to calculate the coordination number (CN) and hydrostatic pressure during loading and unloading of the indenter and the size is 50 Å × 50 Å × 25 Å.Figures 3  and 4 show the distributions of coordination number (CN) and variations in hydrostatic pressure with different crystal surfaces at different stages.From Figures 3 and 4, it can be observed that with the indentation depth increasing, when the hydrostatic pressure reaches 12 GPa, silicon atoms of the specimen underneath the indenter gradually undergo phase transformation from -silicon to -silicon.With releasing of the hydrostatic pressure, part of -silicon (CN = 6) structures gradually recover to the original -silicon and amorphous phase.Moreover, the number of atoms which take place during phase transformation on the (0 0 1) silicon surface is the largest, on the (1 1 0) silicon surface is smaller, and on the (1 1 1) silicon surface is the smallest.The hydrostatic pressure follows the same rule with the number of atoms which take place during phase transformation.This is mainly due to the fact that high hydrostatic pressure easily leads more atoms to take place during phase transformation.This reveals that the (1 1 1) silicon surface is the softest.
Figure 5 shows the load-displacement curves for silicon with three different orientation surfaces.From the unloading curve, the unloading curve of (1 1 1) silicon surface firstly approaches to zero.It is because the atoms on the (1 1 1) silicon surface mainly undergo plastic deformation and phase transformation, and elastic recovery is minimal.
Figure 6 shows the cross-section views of residual Von Mises stress distribution with three different orientation surfaces at different loading stages.From Figure 6, as the indentation depth increases, the Von Mises stress distribution area expands towards deep.The Von Mises stress distribution depth of the (0 0 1) silicon surface shows to be the deepest and the (1 1 1) silicon surface shows to be the shallowest, namely, ℎ 1 > ℎ 2 > ℎ 3 .It is because lower hardness requires smaller force at the same indentation depth, and as a result, shallower Von Mises stress distribution depth is produced.
Figure 7 shows the atomic distribution of three different crystal surfaces and orientations.The following is the analysis of three typical surface structures of monocrystalline silicon which explains why the (1 1 1) silicon surface behaves the softest.From Figure 7 and Table 2, the distribution and related parameters of three silicon surfaces are similar, but the atomic layer distance differs greatly.Different from AAA-AAA pattern of the distribution of the (0 0 1) and (1 1 0) silicon surfaces, the distribution of the (1 1 1) silicon surface presents ABAB-ABAB pattern.Meanwhile, the atomic layer distance A of the (1 1 1) silicon surface is obviously larger than that of the (0 0 1) and (1 1 0) silicon surfaces.When the indenter  loads along ⟨1 1 1⟩ crystal orientation, the atoms obtain more space to deform and the crystal lattice reconstructs in the deformation layer.Therefore, at maximum load, that is, the beginning of the unloading of the indenter, the resilience degree of the atomic deformation layer of the (1 1 1) silicon surface under the indenter is smaller than that of the other two silicon surfaces.As a result, the contact stiffness of the (1 1 1) silicon surface at maximum load is lower; thus the elasticity modulus and hardness is smaller.

Influence of Temperature on Hardness in Monocrystalline
Silicon.Load-displacement curves for nanoindentation on silicon with three different temperatures are shown in Figure 8.Because the temperature has little influence on the hardness of silicon substrate before 500 K, the temperature is set to 296 K, 800 K, and 1200 K in our simulation, respectively [22].From Figure 8, a residual dimple in the specimen is observed after indentation by the diamond indenter.The formation of this residual dimple is due to unrecoverable plastic deformation.The load fluctuates largely in 1200 K. Due to the atomic interaction forces, including repulsive and attractive forces, becoming weaker as the atomic distance increases with temperature increase, the crystal lattice is easy to break and reconstruct.In addition, the required indentation force in 800 K and 1200 K is smaller than that in 296 K at the same indentation depth.These are mainly caused by some defects, such as vacancies and plastic deformation on the surface region.It can be inferred that the specimen material softens at higher temperature.
In Figure 9, during the unloading stage, some materials of the specimen's surface adhere to the diamond indenter, especially in 1200 K.Moreover, the adhesive process leads to greater debris pileup on the indented surface region at higher temperature.It is observed that only the localized region beneath the tip presents high stress distribution due to irreversible plastic deformation in 296 K.However, as temperature increases, the high stress distribution depth of the localized region beneath the tip decreases, and other  regions of the specimen present high stress distribution.It verifies that the specimen material softens and adhesiveness increases at higher temperature.

Conclusion
Using the MD simulation approach, the influences of crystal orientations and temperature on the hardness of silicon were investigated.Based on the above discussions, we concluded the following.
(1) During the indentation process, behavior of local phase transformation is monitored by techniques of coordination number (CN) and RDF.It is observed that phase transformation from -silicon to -silicon and amorphous phase take place when the indentation loading is higher than the yield strength of material.Upon unloading of the indenter, some amorphous phase is usually reversible and the crystal lattice reconstructs.Finally, the residual amorphous layers lead to the formation of the indented surface with a dimple.
(2) Through analysis of three different crystal surfaces in monocrystalline silicon, the atoms on the (1 1 1) silicon surface can obtain more space to deform and the crystal lattice reconstructs in the deformation layer.As a result, the (1 1 1) silicon surface performs the softest.
(3) Under the influence of indentation temperature, the hardness of monocrystalline silicon becomes smaller, but adhesiveness increases as temperature increases.
It is also revealed that some materials of the specimen's surface adhere to the diamond indenter.

Figure 1 :
Figure 1: Molecular dynamics (MD) model of monocrystalline silicon for nanoindentation with a spherical indenter.

Figure 2 :
Figure 2: Front and upward views of the transformed region under nanoindentation at different stages.

Figure 3 :
Figure 3: Number of atoms with specified nearest number of neighbors during nanoindentation.

Figure 4 :
Figure 4: Variation in hydrostatic pressure with different crystal orientations during nanoindentation.

Figure 5 :
Figure 5: The load-depth curves with different crystal orientations during nanoindentation.

Figure 6 :
Figure 6: Cross-section views of residual Von Mises stress distribution with different crystal orientations during nanoindentation.

Figure 7 :Figure 8 :Figure 9 :
Figure 7: The atomic distribution of three different crystal surfaces and orientations.

Table 1 :
Parameters used in Tersoff potential for silicon.

Table 2 :
Crystal surface distribution and related parameter of monocrystalline silicon.