Nanoindentation of polysilicon and single crystal silicon: Molecular dynamics simulation and experimental validation

This paper presents novel advances in the deformation behaviour of polycrystalline and single crystal silicon using molecular dynamics (MD) simulation and validation of the same via nanoindentation experiments. In order to unravel the mechanism of deformation, four simulations were performed: indentation of a polycrystalline silicon substrate with a (i) Berkovich pyramidal and a (ii) spherical (arc) indenter, and (iii and iv) indentation of a single crystal silicon substrate with these two indenters. The simulation results reveal that high pressure phase transformation (HPPT) in silicon (Si-I to Si-II phase transformation) occurred in all cases; however, its extent and the manner in which it occurred differed significantly between polycrystalline silicon and single crystal silicon, and was the main driver of differences in the nanoindentation deformation behaviour between these two types of silicon. Interestingly, in polycrystalline silicon, the HPPT was observed to occur more preferentially along the grain boundaries than across the grain boundaries. An automated dislocation extraction algorithm (DXA) revealed no dislocations in the deformation zone, suggesting that HPPT is the primary mechanism in inducing plasticity in silicon.

(Some figures may appear in colour only in the online journal)

Introduction
An understanding of high pressure phase transformation (HPPT) in brittle materials is of high technological relevance as it induces ductility in brittle materials, thereby enabling their plastic deformation, akin to ductile metals [1]. Silicon is a classic example of a brittle material and it has dominated the consumer electronics market for much of the twentieth century. It has been used in many microelectronic devices, including: solar cells and conducting gate material for MOSFET, MOS and CMOS processing devices. Prior to the fabrication of such devices, crystals of bulk silicon are usually sliced, ground, and polished using cost intensive technologies with an intent to achieve an atomic level flatness [2]. However, at such a finite precision of flatness the surface force dominates, and this influences the response of the NEMS and MEMS devices. Particularly, the adhesion or 'stiction' exhibited by such fine surfaces affects the life of the components because it resists the relative motion between two parts. An understanding of the nanomechanical response of silicon is, therefore, vital to improve device fabrication methods. Accordingly, various techniques have been used [3] to investigate the nanomechanical response of silicon, which include: state-of-the-art veritable resolution using in-situ and ex-situ imaging, quasistatic nanoindentation [4], acoustic emission detection [5], scanning spreading resistance microscopy [6], high temperature contact loading [7], monitoring of electrical resistance [8], x-ray diffraction [9], Raman scattering [10], laser micro-Raman spectroscopy [11] and transmission electron microscopy [12,13]. In addition to the above experiments, simulations (for example molecular dynamics (MD) simulation [14,15], finite element method (FEM) [16][17][18] and multiscale simulation using quasi-continuum (QC) method) [19][20][21] have also provided significant insights into the deformation mechanisms of various materials, including silicon. However, each of these simulation methods has limitations. For example, Sanz-Navarro et al [22] have asserted that the QC method does not permit us to capture the atomistic aspects, such as defect formation, elastic-plastic transition and phase transformations in brittle materials while FEM simulation assumes the matter to be a continuum whereas phase transformation, adhesion, and cohesion mechanisms are of a discrete nature, thus making it difficult for them to be examined using classical continuum mechanics. Similarly, MD simulation suffers from limitations on the time and length scales, making it difficult to validate many parameters directly using the experimental approachable limits. Some of the current research has started to focus on hybrid approaches to MD simulation by validating MD simulation results with laboratory experiments [6,21]. Despite these limitations, MD simulation has gained wide popularity because it enables us to capture the atomistic tribology of simultaneously occurring processes, the foremost of which are structural transformations in the material and mechanochemistry of the process [23]. In the field of materials, although the current state of computational power still does not permit simulation of time scales similar to those in the laboratory, some of the fundamental material properties that are relatively insensitive to time scale can be compared between MD simulations and experiments. In this work, we compare the elastic modulus of silicon obtained from MD simulation and laboratory experiments, and find that the two results correlate well. The subsequent section reviews some of the important aspects of silicon reported in the literature.

Literature review on indentation response of silicon
Silicon is a relatively well-studied material and a wealth of publications have improved our overall understanding of the typical behaviour of single crystal silicon. Under ambient conditions, silicon (Si-I) is brittle because of the sp 3 bonding and diamond cubic crystalline structure, and it contains four nearest neighbours at an equal distance of 2.35 Å [24]. Upon loading (hydrostatic loading of 10-12 GPa or deviatoric loading of about 8.5 GPa), Si-I undergoes metallization, which results in the transformation of its original structure to a Si-II (ductile) lattice structure. Si-II (lattice parameter a = 4.684 Å and c = 2.585 Å) contains four nearest neighbours at a distance of 2.42 Å and two other near neighbours at 2.585 Å [25]. An increase in hydrostatic pressure by about 13 GPa results in the appearance of a newer phase of silicon (Si-XI or Imma silicon) whereas a further increase in pressure by about 16 GPa results in the formation of another phase of silicon which is an eightfold-coordinated simple hexagonal form of silicon [26]. The other phases of silicon recognized to date as a result of HPPT are Si-V, Si-VI, Si-VII and Si-X [27]. It may also be noted that non-hydrostatic pressure can directly transform Si-I to a stable bct-5 (five coordinated) phase of silicon. The bct5-Si crystalline structure contains one neighbour at a distance of 2.31 Å and four other neighbours at 2.44 Å [28]. While loading pressure enforces HPPT, unloading pressure, particularly, the rate of release of pressure results in the appearance of other forms or phases of silicon. For example, upon slow unloading, crystalline phases of Si-XII and Si-III may persist interspersed with an amorphous region. At a pressure of 2 GPa, Si-XII phase contains four nearest neighbours within a distance of 2.39 Å and another neighbour at a distance of 3.23 or 3.36 Å, while Si-III has four nearest neighbours within a distance of 2.37 Å and another unique atom at a distance of 3.41 Å. The main difference between these two phases is that while Si-XII is known to be a narrow band gap semiconductor that can be electrically doped with boron and phosphorus at room temperature to make electronic devices, Si-III is postulated to be a semi-metal [29]. Si-III has also been hypothesized to first transform to a six coordinated Si-XIII phase, which transforms further to either Si-IV or to amorphous-Si (a-Si) [6]. On the other hand, upon rapid release of pressure (i.e. upon rapid unloading) Si-II transforms to tetragonal Si-IX or tetragonal Si-VII phases of silicon. All these phases ultimately stabilize to form a-Si. For a visual understanding of this, several of the major crystal structures of silicon recognized to date are shown in figure 1. Apart from identification of different phases, the influence of indenters in driving these phases has also been explored. For example, a Berkovich indenter causes the formation of  crystalline R8 (Si-XII) and BC8 (Si-III) phases of silicon to appear at the bottom of the transformation zone. However, a spherical indenter will show these crystalline phases closer to the surface and in the middle of the transformation zone [5]. Using a finite element simulation, Zarudi et al [5] demonstrated that it is not the hydrostatic stress distribution but the octahedral shear stress distribution that changes with the shape of the indenter ( figure 2). This change may be responsible for the observation of different phases of silicon with a change in the indenter shape, as mentioned above.
Because silicon is an anisotropic material, HPPT in silicon is also highly direction sensitive; that is, the pressure required for transformation on the (1 1 1) orientation is lower than that required to induce the transformation on the (1 0 0) orientation [27]. Furthermore, most of the metastable phases of silicon are observed to form along the 1 1 0 direction, which is the direction of the slip plane in silicon. Aside from HPPT, nanotwinning has also been recognized to occur during nanometric cutting of silicon [30]. Another interesting observation comes from the work of Jang et al [31], in which the authors have observed a good correlation between the unloading curve and the phase of the material formed postindentation of silicon. Based on the micro-Raman spectra, they proposed that the unloading discontinuity, often known as the 'pop-out', corresponds to the formation of metastable Si-XII/Si-III crystalline phases, while the hysteresis, called the 'elbow' for one-cycle of loading and unloading, is associated with the formation of a-Si. They also reported that higher indentation loads (50 mN) lead to pop-out while lower loads (30 mN) produce the elbow. Clarke et al [32] proposed two arguments in support of their observation of a-Si during indentation of silicon with Knoop and Vickers indenter: (i) at rapid unloading rate, Si-II phase cannot transform to another crystalline phase due to the kinetic barriers in the nucleation and growth; and, (ii) direct amorphization of Si-I phase occurs, and this amorphous phase persists even upon unloading because of insufficient thermal energy to cause back transformation to Si-I. The first explanation was considered more plausible since direct amorphization will require a higher transformation pressure of 24 GPa, which is considerably higher than the hardness of Si [31].
Overall, a metabody of knowledge is available in the literature describing the response of single crystal silicon under various cases of contact loading. However, there seems to be no study in the literature that explores the nanomechanical response of polycrystalline silicon substrate, which is widely used in most real world applications (such as solar cells). Moreover, the size of the grains in most commonly used engineering materials (such as steel, aluminium, silicon and even reaction bonded silicon carbide) is in the range of 100 nm to 100 µm, which is around the feature size of micro-machined components. Therefore, the aspects of the microstructure of the workpiece material are important considerations for its machinability. For example, when the depth of cut and the cutting edge radius are of the order of the workpiece material grain size, the cutting tool will experience varying resistance from different grains and from the presence of grain boundaries, which could be critical when it is necessary to generate a smooth machined surface on a polycrystalline material. In this paper, we aim to contribute to the scientific knowledge by studying and comparing the responses of polycrystalline silicon and single crystal silicon, both with a spherical (arc) and a Berkovich pyramidal indenter.

MD simulation inputs
MD simulation comprises a numerical solution of Newton's equations of motion for a set of atoms [33].
The trajectory obtained from the simulation is analyzed for the thermodynamic changes that take place over a certain interval of time. In this work, the 'Large-scale atomic/molecular massively parallel simulator' (13 September 2013 version) [34] was used to perform a series of MD simulations. VMD [35] and OVITO [36] were used to visualize and analyze the atomistic simulation data while an automated 'dislocation extraction algorithm' (DXA) [37] was used for automated identification of crystal defects from the output of the MD data. The simulation model of the polycrystalline silicon substrate was generated by using an in-house developed Voronoi tessellation code [38,39]. The simulation uses an analytical bond order potential (ABOP) energy function proposed by Erhart and Albe [40], which is a trade off for the computational power and covalent bond interactions between silicon and carbon against the previously used Morse potential function in the past [41][42][43]. Further details of the parameters used to develop the  MD simulation model used in this work are shown in table 1, which can readily be used to replicate the simulation results. It may be noted that the MD simulation model developed in this work considers a vacuum surrounding, which might be a rare instance in practice. However, this consideration turned out to be fortuitous because discarding the role of ambient air permitted careful examination of the effect of the crystal structure alone. The MD algorithm and the necessary considerations required and used in this work are quite similar to other nanoindentation simulation studies done by the authors [1,44], and are not repeated here for brevity. Choi et al [45] have discussed some implications of the boundary conditions and accordingly the model developed in this work also assumes periodic boundary conditions along X and Z directions, where the bottom most layer in Y direction was kept fixed. The indenter material used was diamond, which was kept rigid throughout the simulation. A snapshot of the polycrystalline substrate for a spherical and a pyramidal indenter obtained from the simulation after the equilibration process is shown in figure 3.

MD simulation results
Gannepalli and Mallapragada [46] noted that a pop-in event in the P -h curve signifies plastic deformation in the brittle materials, in that the width of the peak during the pop-in event signifies the duration of the event while its height signifies the magnitude of the event. Also, each maxima in the P -h plot corresponds to the onset of the plastic activity, while each corresponding minima reflects the completion of the event. Unlike experiments, MD simulation provides a much better temporal and spatial resolution of each atomic scale event in the P -h plot with extreme accuracy. Figure 4 shows the P -h plots for the four different cases in the MD simulation. It is interesting to note that, although the peak depth (the maximum depth of indentation was kept fixed in all cases because it corresponds to displacement controlled nanoindentation) in all cases was the same, the peak maximum force was found to vary in each simulation run. In the simulation, the maximum force was found to be significantly different for the two indenters for both types of silicon (polycrystalline and single crystal). However, the maximum force was only marginally different for the same indenter. Another important difference in these plots is the residual depth of recovery (h f ). The residual depth of recovery governs the extent of plastic response in the substrate. By comparing the results shown in figure 4, it can be seen that the adhesion (negative indentation force) between the indenter and the substrate is larger for a polycrystalline substrate than for a single crystal silicon substrate.
Since the influence of the crystal structure is seen to influence the peak forces, it was necessary to explore this point further (this phenomenon is explored and discussed at length in the next section). A close examination of the MD simulation trajectory revealed that the grain distribution, grain size, shape and orientation of each grain in the polycrystalline material are arranged in a random manner. Indentation of such a polycrystalline material causes the indenter to encounter different crystallographic orientations and crystal directions of each grain. This causes the response (reaction force) of each grain to differ in magnitude. During indentation of a polycrystalline material, at infinitesimally smaller depths of indentation, only an individual grain may undergo deformation. Ironically, this situation is similar to indentation of a bulk single crystal material. However, at higher depths of indentation the combined response of several grains, in addition to the influence of grain boundaries, may cause a combined reaction force; thereby, making the response of a polycrystalline material different from that of a single crystal material. This results in changes in the indentation force as well as the extent of the sub-surface deformation (its extent and its properties). An instance for this phenomenon was captured during the simulation. The indentation of the single crystal material with a spherical indenter showed a symmetrical degree of amorphization on both sides of the indenter from the centre, whereas the degree of amorphization in the sub-surface of the polycrystalline material was found to be significantly different ( figure 5). The sub-surface of the material was seen to be non-symmetric in its properties from the centre of the indenter, as shown in figures 5(a) and (b). Figure 5 also shows the variation in the coordination number, which is now an established way of monitoring phase changes, although only up to a certain level of accuracy [15]. The bulk silicon has four nearest neighbours and, hence, the pristine silicon substrate has a coordination number of four, indicating the ideal diamond cubic lattice geometry of silicon. In figure 5(b), it can be seen that the lattice interspacing between the two kinds of atoms on the left side and on the right side of the indenter is different, which highlights differences in the two types of grains being stressed by the indenter simultaneously. The amorphization of silicon is an outcome of HPPT that occurs due to the inhomogeneous distortion of the tetrahedral bonding in the direction of the slip plane ( 1 1 0 in silicon). This explains why the microstructure of the mechanically processed polycrystalline silicon in the sub-surface and on the plastically deformed surface may exhibit different properties than those in the single crystal material.
To detect the existence of any dislocations induced by the indenter, we used an automated DXA [33], the results of which are shown in figure 6. Unlike metals where dislocation nucleation has been recognized as a primary source of plasticity, no dislocation nucleation was observed during nanoindentation of silicon. However, it may be noted that the event of dislocation is highly sensitive to the parameters of the indentation (indentation depth, size of the indenter, and crystal structure of the material such as defects, voids, as well as residual stresses in the substrate [47], environment etc.). Current limitations on MD software limit the scaling of indentation simulation to the experimental scale. This could possibly explain why no dislocations are evident in this work in contrast to findings in some microscale cutting experiments [48]. Table 2 details the important results obtained for all four of the test cases of the simulation. An interesting aspect to note in table 2 is that the ratio of residual indentation depth to the maximum depth of the indentation in all the cases was observed to be different. This ratio was higher for the single crystal material than for the polycrystalline material. A higher ratio of h f /h max in single crystal silicon shows that it deforms more plastically than the polycrystalline material.
Another interesting observation in table 2 is the measure of von Mises stress (the von Mises stress criterion is a very frequently used yield criterion to predict the yielding of a material due to the maximum deviatoric strain energy) to drive the phase transformation in Si-I silicon. In this particular simulation, the von Mises stress was found to be maximised for the polycrystalline material in comparison to that in the single crystal material when the spherical indenter was used. For the pyramidal indenter, this observation is in the opposite direction. Furthermore, the polycrystalline material showed different values of the peak von Mises stress in the indentation zone. When polysilicon was indented with a spherical indenter, the peak von Mises stress was about 13.76 GPa, whereas the peak von Mises stress during indentation of a polycrystalline material with a pyramidal indenter was only 7.16 GPa. To examine this effect more carefully, we focused on the granular structure (figure 7). Figure 7 highlights the indentation performed in a polycrystalline substrate (showing grain boundaries) with two different kinds of indenters.
It can be seen from figure 7 that the pyramidal indenter only entered a single grain, which was possibly aligned in a crystallographic orientation and a direction that is more amenable to deformation than the (0 0 1) crystallographic  orientation (used for single crystal silicon). Therefore, the maximum load and the von Mises stress in polysilicon were observed to be lesser than that in the single crystal material for a pyramidal indenter. It can also be seen that the presence of grain boundary resisted the process of structural transformation across the grain boundary. Therefore, it can be asserted that HPPT in polycrystalline silicon occurs more preferentially along the grain boundary than across the grain boundary. Consequently, indentation with the spherical indenter required extra energy to overcome the phase transformation across the grain boundaries, as well as to overcome the combined response of several grains having different properties than the single crystal material, thus requiring higher transformation pressure and exerting larger indentation forces on the indenter.  As early as 1975, Gilman asserted that nanoindentation hardness of silicon is insensitive to a temperature of about 773 K (∼0.45 times the melting point of silicon) [49]. The maximum temperature in the deformation zone (table 1) in all the cases was observed to be lower than 773 K. During nanoindenation, and even during nanometric cutting, the maximum temperature observed in the primary shear zone of silicon was up to about 759 K [50] (lower than 773 K). Combined together, these observations point to the fact that plastic deformation or yielding of silicon during contact loading may not be driven by temperature but by stress. This assertion lends further credence to the theory of Gilman [51], who proposed that twinning (an adiabatic diffusionless process occurring in monoatomic materials) is the mechanism underlying the occurrence of HPPT in brittle materials. Gilman argued that HPPT resembles a quantum tunnelling process in that both processes (HPPT and quantum tunnelling) are analogous to each other. An electric field (E) is responsible for tunnelling (qE bring the driving force; q is the electric charge) while the applied shear stress (τ ) is responsible for twinning, where the driving force is τ a (a = twin plane area per atom). Twinning differs from plastic glide because, unlike the slide, the shear displacement during twinning causes the lattice orientation to change. It may also be noted that high strain rates, low stacking-fault energy, and low temperature all facilitate twinning. Since twinning may occur at low temperatures, HPPT appears to be a stress-driven process rather than a thermally activated process. This seems particularly to be the case in a nanoindentation, which is a compressive stress-dominated process rather than a concentrated shear stress process, such as nanometric cutting, although both are contact loading processes. It is thus interesting to note that while nanoindentation and nanometric cutting are both classified as contact loading processes, they have certain major differences that have not previously been explicitly highlighted in the literature. To this end, we propose some differences (table 3) between nanoindentation and nanometric cutting.

Sensitivity of MD simulation results to the velocity of indenter
A drawback of the MD simulation is that the computational time and resource limitation do not permit scaling of the simulation parameters to the experimental scale. The length scales in this regard are restricted to up to few nanometres (or Table 3. Differences between nanoindentation and nanometric cutting.

S. No. Differences
Nanoindentation Nanometric cutting 1 Stress dominance Hydrostatic stress Deviatoric stress 2 Effect of the process Induced compressive strain results Induced shear strain results in pronounced in the pronounced change in the bond length. change in the bond angle. 3 Strain energy field Dimension of the indentation is Cognate to nano-indentation, depth of cut considered a criterion.
is normally chosen as a criterion. 4 Contact of the tool Point contact Line contact 5 Process outcome Pile-up formation (may be of irregular shape). Both pile up and cutting chips (continuous or discontinuous type). 6 Wear of the tool Trivial Non-trivial 7 Cutting action Role of the cutting edge radius is dominant Role of tool rake face is more dominant during a nano-indentation process.
during the nanometric cutting. 8 Direction of execution On the same cutting plane, nano-indentation is On the same cutting plane, cutting is done usually done in the direction normal to the plane.
along the direction of the cutting plane. 9 Duration of action Indentation depth is usually the parameter Track length is usually specified to used to specify the duration which is quantify the cutting length which is usually small in magnitude.
relatively quite large in magnitude. 10 Forces Unidirectional loading and unloading Two forces-tangential cutting forces and thrust forces are used to characterize the process.
forces-are used to characterize the process. 11 Toolpath The tool enters the substrate and retracts Usually the tool is retracted at a different along the same tool path.
tool path otherwise the finished surface may be destroyed by the tool. 12 Application To measure the hardness or elastic To measure specific cutting energy and/or modulus of the sample. study the chip flow process. few million of atoms). In addition, the velocity of indentation, or velocity of cutting, such as that used in the current simulation (50 m s −1 ) against an experimental speed of upto ∼1 m s −1 in nanometric cutting or 0.1-10 µm s −1 during nanoindentation, is a factor that might cause unexpected spurious effects. It may, however, be anticipated that a higher loading rate, such as the one used during the current work, leads to high strain rate deformation. It is, therefore, necessary to explore the effect of such strain rates on the MD simulation results. This section details the investigation of the effect of loading rate on Si-I in Si-II phase transformation. Several additional trials were carried out on single crystal silicon at different indentation and retraction speeds to explore such effects and the results of these trials are listed in table 4 and plotted in figure 8. Figure 8 shows that the loading curve is not affected due to the differences in the indentation speeds as much as the unloading curve, which leads to the appearance of different phases of silicon during unloading, as is well documented in the literature. To analyse these differences we have used the parameter Coefficient of Variation (CV), defined as the ratio of standard deviation to the arithmetic mean, which eliminates the effect of sample population since it is a normalized measure of dispersion of a probability distribution. It may be seen from table 4 that, even when the indentations were performed at varying speeds, the CV by contrast is much higher for the peak temperature (0.1891) than the peak von Mises stress (0.042). The variation in the stress and temperature is also plotted in figure 9, and it may be further seen that an increase in the indentation speed causes only a marginal change in the peak stress whereas the temperature rise is more significant. The above observation strengthens the earlier observation of Niihara [52], who asserted that the influence of large hydrostatic stresses (such as those existing during nanoindentation) could lead to the plastic deformation of almost any material (including super-hard substances like diamond), even at low temperatures. Especially in brittle materials (including in this work), HPPT is observed to induce ductility in silicon. Finally, to explore whether Si-I to Si-II transformation is a stress-driven or a temperaturedriven process, we plotted and projected the local stress and temperature during indentation on the pressure-temperature phase diagram of silicon, as shown in figure 10 with respect to  different speeds of nanoindentation. By extrapolating this line along a lower temperature it can be seen that at a sufficiently lowered indentation speed the formation of Si-II is achievable by the virtue of stress alone and is not due to the temperature and, hence, the observation of HPPT being a stress-driven process appears to be true.

Experimental validation of MD simulation
To validate the nanoindentation simulation results, we used the Oliver and Pharr method [54] to estimate the Young's modulus of the silicon (both single crystal and polycrystalline) in all the simulation test cases. The Oliver and Pharr method enables the estimation of the mechanical properties directly from the plot of P -h curve. As shown in table 2, MD simulation P -h plots revealed values of about 135 and 165 GPa as Young's modulus of single crystal and polycrystalline silicon, respectively, when a spherical indenter was used. In contrast, a much higher value of 220 GPa was obtained from the simulation indentation data of the pyramidal indenter. It may be seen that the MD simulation has provided key information concerning the role that the shape of the indenter and crystal structure of the substrate material plays in influencing the deformation mechanism of silicon. However, as ostensibly evident from table 1, the edge radius of the spherical indenter size was extremely small (∼10 nm), which is not comparable to the commercially available experimental scale. In an attempt to validate the simulation results, we performed a displacement controlled quasistatic nanoindentation experiment on a single crystal silicon specimen using a three sided pyramidal Berkovich indenter at extremely fine indentation depths of several nanometres. The nanoindentation tests were performed on a TI 900 Hysitron TriboIndenter at a room temperature of about 293 K. This setup takes advantage of an acoustic and thermal enclosure, which enables capture of precise and sensitive readings [55]. In addition, the patented capacitive transducer provides superior sensitivity and stability over other similar instruments. The specimen used was a single crystal silicon wafer with crystal orientation (0 0 1), with a diameter of 50 mm and thickness of 5 mm. It was polished on both sides. The tip of the indenter was noted to be blunt, having an edge radius of 300 nm as opposed to a newly procured tip radius of 150 nm. However, this blunt nature turned out to be a benefit rather than an experimental difficulty because the blunted geometry of the nanoindenter can often be approximated as spherical. It may be noted that the shape of the indenter may significantly influence the appearance of HPPT phases of silicon; that is, a sharp geometry as that of a Berkovich indenter (pyramidal indenter) promotes shear stress to cause plastic contact whereas a spherical indenter (blunt indenter) initiates elastic contact. This is due to the fact that at an identical penetration depth, the total deformation energy of silicon required by the spherical indenter is more than a conical indenter [28].
During the experiment, a displacement control feedback system was chosen over a load controlled feedback system so as to achieve a finite indentation depth [56]. The time allowed for reaching maximum displacement in all the cases was 10 s and the indenter was retracted immediately after attaining the peak indentation depth in the same duration of 10 s. Each indentation was conducted using a quick approach method so as to ensure the accuracy of depth measurements sensed by the indentation probe. Since the modulus of the material is a more fundamental property, it can be compared even at a relatively higher indentation depth. Therefore, we performed  an indentation at a depth of 15 nm so that the substrate size effects and sensitivity effects of the ambience on the results can be eliminated.
The indentation result obtained from the instrument is shown in figure 11. Only one result is shown here for brevity's sake, but in general the average value was observed to be quite close to this result. The reduced elastic modulus obtained from the instrument was E r = 132.8 GPa, which is in excellent agreement with the value obtained from the MD simulation during indentation of single crystal silicon from a spherical indenter using MD simulation (table 2). The small difference plausibly originated from several factors, such as sample roughness, air lubrication, sensitivity of the equipments, purity of the material and accuracy of the measurements etc. Moreover, the experimental value of Young's modulus of silicon reported in the literature is in the range of 120-170 GPa [57], depending on the orientation of the wafer, the parameters used for the indentation, and the sensitivity of the measurements. Thus, an excellent agreement in the values of elastic modulus of silicon obtained from MD simulation and the experiment provides confidence in the simulation results.

Conclusions
We have used a MD simulation to simulate the nanoindentation of a polycrystalline silicon substrate and a single crystal silicon substrate to understand: (a) the influence of the shape of the indenter (Berkovich pyramidal and spherical), and (b) the role of the crystal structure of silicon in influencing the state of the HPPT, which induces the well-known Herzfeld-Mott transition in silicon. HPPT in silicon (Si-I to Si-II phase transformation) was observed to be a common occurrence in all the cases, but its extent and the manner in which it occurs differed significantly. The varying sub-surface in a mechanically processed polysilicon arises due to this phenomenon. The most significant result obtained from the simulation was that the HPPT occurs more preferentially along the grain boundaries than across the grain boundaries. This single phenomenon results in an increase of the transformation pressure to a significantly higher extent and, consequently, even the indentation forces were observed to increase. The use of an automated dislocation extraction algorithm (DXA) revealed no dislocations, thus confirming HPPT to be the sole mechanism of plastic deformation of silicon in all the cases. Furthermore, HPPT appears to be a stress-dominated process rather than a thermally activated process and twinning appears to be the mechanism underlying the occurrence of HPPT in brittle materials (twinning resembles a quantum tunnelling process). The Oliver and Pharr method was extended to validate the simulation results and an excellent agreement was found between the MD simulation and the experiments.