Molecular Dynamics Simulation of Silicon Carbide Indentation Behavior

1. Introduction Silicon carbide and its various polymorphs are known for their exceptional chemical, physical-mechanical and electrical properties, e.g. chemical inertness, hardness, heat resistance, low thermal expansion and semiconductor properties due to different band gaps [1] [2]. Because of these remarkable properties silicon carbide (SiC) is used in a variety of applications and industrial usage. It shows a high sublimation temperature and a large bandgap, making it useful in high temperature electronics. As semiconductor device SiC plays a prominent role in electronic vehicles and is of growing interest in industrial fields where high voltages and high temperatures have to be controlled. From a physical point of view SiC has a high hardness and is chemically inert, ensuring its durability even in harsh environments. Silicon carbides are essential in the synthesis of ceramics due to specific mechanic properties (e.g. hardness) or (thermo-)dynamic features, e.g. thermal expansion and thermal conductivity. Heat-shock resistance is of particular interest for various applications of SiC-ceramics, e.g. in turbines. Another characteristic feature of SiC is its polymorphic occurrence, which comprises Bravais lattice types, i.e. rhombohedric, hexagonal and cubic structures. This structural diversity gives rise to a broad range of applications caused by variability in properties. The different structures are classified by Ramsdell and can be described in detail using numbers of layers and an indication of the Bravais lattice type [3]. Differences in stacking sequences (e.g. ABCB or ABCACB, for 4H- or 6H-SiC, resp.) lead to differences in properties of polytypes. In this work we will focus on 3C-SiC and both 4H- and 6H-SiC, which are the most common polytypes. 3C-SiC has a cubic structure, while 4H- and 6H-SiC belong to the hexagonal polymorphs. The number refers to the repeating layers in the unit cell. 3C-SiC shows the smallest bandgap among the three and the highest thermal conductivity, which enables use in the area of air conditioning and LED lighting system [4]. 4H-SiC has the greatest bandgap while 6H-SiC has a lower hardness compared to the others [5]. Among the hexagonal forms, 6H-SiC is of particular interest because of its pronounced anisotropy and complex deformation mechanisms. Its structural response to mechanical loading, including elastic-plastic transition and amorphization, has been reported to differ from that of 3C- and 4H-SiC [6] [7]. In a previous study, we investigated the energetic bonding situation between different SiC structures and graphenes and their mutual influence on their properties [8]. In this work, we model the structural, mechanical and dynamic properties of silicon carbides as well as the influence of mechanical stresses or disturbances by hemispherical indenters using molecular dynamics. Fracture behaviour, hardness and mechanic stability of SiC are experimentally examined properties, and various studies for the different polymorph (e.g. 3C-SiC [9], 4H-SiC [10], 6H-SiC [11]) exist. The aim of our simulations is to illustrate the internal processes that occur during fractures or damage to the material at the atomic level. Therefore, we simulate the changes in properties caused by the application of indenters to the aforementioned SiC crystals. Molecular dynamics simulations are used to investigate indentations in the three different polytypes. For this, we employ a method similar to the Brinell test, in which spheres are pressed into the material under test. In our simulations, the size of the indenter and the number of time steps during which the indentation and load act on the structures are varied. The responses of the different SiC polymorphs are compared with respect to energy changes and the pair correlation function. The atomic displacements are analyzed using the root mean squared displacements (RMSD) values. We gain further insights into the dynamics of the loaded structures by determining the (atomic) vibration spectra using Fourier transformation of the velocity autocorrelation function (VACF). In the next section, we will discuss in detail the potential used to describe the interatomic interaction, the simulation parameters for the indenter, the boundary conditions used, and the observable quantities. The most important results will be presented and discussed in the following section. Finally, we will summarize the most important results and observations. 2. Methods The simulations are performed using classical molecular dynamics, which has been state of the art for over three decades [12]-[15]. However, it should be noted that in nanosystems, the chemical and physical processes are determined and dominated by quantum effects, which are not considered in this type of simulation. Using numerical integration, the Newtonian equations of motion are solved, and only atomic/vibronic behavior is observed. However, since the system sizes under consideration exceed 1500 atoms, these simulations go beyond the typical scale of quantum mechanical calculations. We investigated the 3C-SiC system with N = 1728 atoms, the 4H-SiC system with N = 2400 atoms, and the 6H-SiC system with a total of N = 2400 atoms. The MD integration time step is δt=1 fs. The systems are heated from 8 K to 300 K within 50,000 fs and then equilibrated at this temperature for a further 50,000 time steps. In the simulations, the temperature is averaged over 20 fs and then adjusted to the specified temperature by scaling the atomic velocities. The indentation tests are then performed using these equilibrated structures. 2.1. Potential For the simulations of the silicon-carbon structures, we use a potential which was initially developed to describe graphite or diamond structures as well as hydrocarbon molecules. This potential was proposed and adapted by D.W. Brenner [16] and was extended to describe SiC systems [17]. This empirical potential is well suited to reproduce the structural properties and dynamic characteristics of a variety of Si-C systems. The binding energy E is given by a potential consisting of two-body terms (considering short-range interactions of attractive and repulsive nature). E= ∑ i ∑ j>i [ V R ( r ij )− B ¯ ij V A ( r ij ) ] (1) In the equation is V R ( r ij ) the repulsive part and V A ( r ij ) the attractive part of the potential. The bonding situation for the atoms i and j is averaged and described with the last factor B ¯ ij . The detailed functional formulas of the equations are as follows: V R ( r ij )=f( r ij ) c ij s ij −1 e − 2 s ij β ij ( r ij − r ij,0 ) V A ( r ij )=f( r ij ) c ij s ij −1 s ij e − 2 s ij β ij ( r ij − r ij,0 ) B ¯ ij = ( B ij + B ji )/2 + F( N i t , N j t , N ij conj )/2 (2) The potential function with parameters c ij , s ij , β ij and r ij has the same type as a Morse potential. These parameters were adjusted to represent type-specific properties (energies and length scales). The function f( r ij ) has only a short interaction range and is defined by: f( r )={ 1 forr< r ij,1 1 2 ( 1+cos( π r− r ij,1 r ij,2 − r ij,1 ) ) for r ij,1 ≤r< r ij,2 0 forr≥ r ij,2 (3) This function ensures that the interatomic interaction within the interaction range steadily decreases to zero. The parameters of the extended Brenner potential are specified in [17]. As a reactive empirical bond order potential (REBO potential), the Brenner potential can particularly model and describe bond breaking and changes in hybridization; therefore, it is very well suited to represent the interatomic interactions in SiC crystals under external influences, e.g. high-stress phenomena. 2.2. Conditions and “Experimental Setup” By applying periodic boundary conditions, we want to reduce the surface effects generated due to the size of the system. For this purpose, the initial volume is kept constant during the simulation and the virial of the configuration is monitored to consider the internal pressure [18]. To generate the input structure for the penetration test, the crystalline structures are heated to 300 K within 50,000 fs and equilibrated at this temperature for another 50,000 fs. The process of indentation is described by placing a rigid hemisperical indenter on top of the surface of the SiC-crystals. The indenter size is varied between 5 and 30% of the overall sample size. At the beginning of the simulation, the indenter merely touches the surface. The indenter is then slowly shifted downwards into the crystals and the interaction with the atoms of the silicon carbides is modelled via a harmonic potential which repells and shifts the atoms lying in the zone touched by the indenter [19]. The energy of the indenter acts only on atoms within the indentation zone. This energy is described by the equation E( r )=0.5⋅k ( r−R ) 2 , where R is the current penetration depth of the hemisphere and r is the distance of the atoms from the center of the sphere. The force constant is set to 7 eV/Å2, which is similar to diamond. Atoms outside the zone are only indirectly affected by the intruder through interatomic interactions with the displaced particles. To simulate indentation into the test sample, the boundary condition in the respective direction is replaced by reflective boundary conditions. For the analysis of the indentation of silicon carbide, varying parameters are used. In addition to the indenter size mentioned above, the dependence of the sample behavior on the direction from which the indenter impacts the sample is also investigated. Furthermore, the simulation time required for the indenter to penetrate the SiC crystal is varied. Three different polymorphs of silicon carbide are investigated, namely the cubic 3C-SiC and the hexagonal 4H- and 6H-SiC. For 3C-SiC, all tests are performed in the z-direction. To consider the anisotropic properties of the 4H- and the 6H-SiC, tests are performed from all directions. The duration of the computer experiments ranges from 15,000 to 200,000 fs. The system sizes of the structures under consideration range from 1728 atoms (3C-SiC) to 2400 atoms (4H-SiC and 6H-SiC). For all tests, the temperature is kept at 300 K. 2.3. Properties To evaluate the material response of different SiC polytypes (3C-, 4H-, and 6H-SiC) under indentation, we analyzed several physical properties during the molecular dynamics simulations. These include energy evolution, atomic displacements, pressure and temperature development, as well as structural correlation measures and vibrational spectra. In molecular dynamics simulations, the total energy of a system is typically expressed as the sum of the kinetic and potential energy contributions. The kinetic term reflects the atomic motions and is directly related to the system temperature, while the potential energy arises from interatomic interactions described by the chosen force field [16] [20]. For covalently bonded materials such as silicon carbide (SiC), these potentials are parametrized to reproduce the directional C-C, Si-C, and Si-Si bonds, ensuring that the computed energetics reliably represent the response of the crystal under external loading [20]. During nanoindentation, the evolution of the potential energy provides direct insight into the deformation mechanism. In the early stage of loading, the increase in potential energy reflects the elastic regime, where atoms are displaced from their equilibrium positions but the overall lattice remains intact [7]. This phase is reversible, and upon unloading the energy would return close to its initial value. Once the indentation depth exceeds a critical value, however, a deviation from the elastic behavior is observed: the energy increase is accompanied by local bond breaking, dislocation nucleation, and the formation of amorphous regions, marking the onset of plastic deformation [7]. In this case, part of the energy remains stored in the system even after unloading, consistent with irreversible processes. The radial distribution function, often denoted as g( r ) , is a statistical measure that describes the probability of finding an atom at a distance r from a reference atom, relative to the probability expected for a completely random distribution at the same density [20]. The determination of the partial pair-correlation functions is made by monitoring: g( r )=〈 n( r ) 4π r 2 ρΔr 〉 (4) With n( r ) being the number of atoms at a distance r in a sphere with thickness Δr . This sphere is located around a reference atom, ρ is the number density of the atoms in the systems and 〈 〉 stands for the averaging over configurations which are calculated throughout the MD simulations. In crystalline materials, sharp peaks appear at well-defined positions corresponding to the nearest-neighbor distances and lattice periodicities (e.g., C-C or Si-C bonds). The height and sharpness of these peaks provide information about the degree of short- and medium-range order. Well-defined peaks indicate a crystalline lattice, while broadened or flattened peaks point to increasing disorder or amorphization [7]. In nanoindentation, the analysis of the radial distribution function (RDF) can reveal the structural changes caused by elastic and plastic deformation. A reduction in the height of the first coordination peak typically signals bond distortion or breaking, while the emergence of broader or additional peaks indicates the formation of amorphous regions [6]. In this study, we also determine the partial pair correlation functions for those atoms in the uppermost layers that are likely to be most strongly influenced by the indenter. This allows for a better understanding of the structural changes in the affected region. The root mean squared displacements (RMSD) is a statistical measure that quantifies the average displacement of atoms from a reference configuration during molecular dynamics simulations. It is calculated as the square root of the sum of the squared differences between the instantaneous and reference atomic positions: ΔR= ∑ i ( R i ( t )− R i ( 0 ) ) 2 (5) In crystalline materials, the RMSD therefore reflects the degree of deviation from the ideal lattice structure over time. A gradual increase in RMSD corresponds to elastic distortion, while a sharp rise marks the onset of plasticity and irreversible structural change [6] [7]. Vibrational spectra derived from molecular dynamics simulations provide insight into the dynamic behavior of atoms and the stability of crystal lattices under external perturbations. The most common approach is the calculation of the vibrational density of states, which describes the distribution of vibrational modes over frequency. The vibrational density of states can be obtained from the Fourier transform of the velocity autocorrelation function, linking the atomic trajectories generated in MD to the collective vibrational dynamics of the material [20]. f v ( t )= 〈 ∑ i=1 N v i ( 0 ) v i ( t ) 〉 〈 ∑ i=1 N v i ( 0 ) v i ( 0 ) 〉 (6) By using a cosine-transformation on the velocity-autocorrelation function, we determine the vibration density of states (DOS): Z( ν )= 1 Z 0 ∫ 0 t obs cos( 2πνt )exp( −λ t 2 ) f v ( t )dt (7) λ is a damping factor which broadens the δ peaks in the frequency spectrum [21] [22]. t obs is the observation time for the velocity-autocorrelation function and Z 0 a normalisation constant [18]. In crystalline systems, the vibrational spectrum exhibits well-defined peaks corresponding to optical and acoustic phonon modes. A sharp and well-structured spectrum is characteristic of ordered states and a mechanically stable lattice. Broadening or shifting of these peaks, on the other hand, indicates bond weakening, anharmonicity, and the development of disordered or amorphous regions [7] [20]. 3. Results and Discussion In this section we first present and discuss the results for the simulations of the cubic 3C-SiC structure for which the indentations are only performed in z-direction (due to the isotropy of the system). As second system we discuss the observations in the 4H-SiC crystal. Last but not least our findings for the 6H-SiC configuration are given and discussed. The simulations yield the following results. 3.1. 3C-SiC For the cubic form of SiC we simulate the indentation process from only one direction (i.e. from z-direction) due to the isotropy of the structure. The test procedure is varied with respect to size of the indenter (a hemisphere) and the duration of time during which the indenter is pressed into the crystal. For 3C-SiC several indentation hardness tests are simulated. Table 1 shows the parameters used in the simulations of the penetration of a hemisphere into 3C-SiC (size of the indenter relative to the crystal height, duration of the simulation and indentation speed). Table 1. Simulation parameters for the indentation tests of 3C-SiC. Indenter size [%] Number of time steps Indentation speed [m/s] 5 15,000 8.72 10 20,000 13.08 15 25,000 15.70 20 50,000 10.46 25 62,500 10.46 30 75,000 10.46 First let us consider the time dependence of the energy. In Figure 1, we show the energy development depending from the indenter size. At the beginning of the simulation, the energies behave similarly for all intruder sizes. The initial energy drop is due to the change in boundary conditions associated with the placement of the hemisphere and the start of the indenter test. For small intruders (up to 15%), no significant structural changes occur, and the energy remains almost constant. For larger intruders, the crystal reacts with volume reductions, followed by deformations or partial amorphization of the upper layers near the edge. The larger the intruder, the greater the volume reduction and the more pronounced the structural deformation. For the largest intruders (25% and 30%), the energy increases at the end of the simulation [23] [24]. As already mentioned, at the beginning of the simulations (up to approximately 15,000 fs), no significant differences between the different indenter sizes on the energy profile can be observed. The energy decreases from 25,000 fs onwards for the 30% indenter. As we can clearly see when the energy profile is plotted against the simulation time (cf. Figure 1). For the indentation with 25%, this happens after about 35,000 fs and for 20% it starts after 40,000 fs. Therefore, a larger indenter leads to an earlier response in the elastic and inelastic deformation. For the indentation with smaller indenter sizes, only minor changes in energy appear [25]. Figure 1. The energy profiles of 3C-SiC at room temperature are plotted against the simulation time for different indentation depths in the z-direction. Atomic shifts are taken into account by determining the RMSD (root mean squared displacement). The RMSD values at different indenter sizes are shown in Figure 2. The RMSD plotted over time show that differences also start to occur from about 25,000 fs onwards. This is also visible in the graphic representations generated with VMD [26]. Here, the atoms move closer together in the simulations with an indenter size of at least 20%. A displacement from the original position of the atoms would lead to a higher RMSD. The largest atomic displacements upto 90Å are found for the indentation simulations with indenter of 30% size, whereas for indenter sizes less than 15% displacements less than 20 Å occur. In Figure 3, we show the last configuration of 3C-SiC after the application of an indenter with size 30% and a simulation time of 75,000 fs. Despite the disturbance by the indenter only a roughening of the surface can be identified. The main part of the structure remains unchanged. Figure 2. The root mean square displacements—RMSD in Å of 3C-SiC at room temperature are plotted against the simulation time for different indentation depths. Figure 3. The last configurations of 3C-SiC after the application of indenter with sizes of 10%, 20%, and 30% and simulation times of 20,000 fs, 25,000 fs and 75,000 fs, resp. (from left to right). Silicon atoms are represented by yellow spheres, while carbon atoms are represented by blue spheres. The pair correlation functions are considered for the carbon-carbon pairs, the silicon-silicon pairs and the carbon-silicon pairs. The pair correlations functions obtained from the MD runs with different indenter sizes are compared the pair corrleations funcions of the undisturbed crystalline structure, as shown in Figure 4. In the case of C-C pairs and Si-Si pairs, some shorter distances occur before the respective first main peak, this can be observed particularly well in the pair correlation functions of those atoms that lie in the layers most strongly influenced by the indenter (see right side of this figure). The C-C distances are particularly affected in the disturbed region, and significantly shorter distances occur between these atoms. The Si-Si distances are only significantly reduced at high penetration depths of 20% or more. However, this effect affects only a few atoms compared to the entire solid, so the positions of the main peaks remain unaffected. Therefore, none of the pair correlation functions show any deviations in the distances between the atoms, i.e. the peak positions remain unchanged. Only the pair correlation function for the indentation of 30% is slightly different in its course (cf. Figure 4). This could indicate to the fact that the indenter sizes were too small to lead to changes in the atomic order for the 3C-SiC. Figure 4. On the left side the pair distribution functions g( r ) of C-Si pairs (top), C-C pairs (middle) and Si-Si pairs (bottom) are shown for different indentation depths in 3C-SiC. For comparison, the distribution from a simulation without indentation is shown. The right panel depicts the corresponding pair distribution functions for the subset of atoms most strongly influenced by the indenter, for better illustration, the curves have been enlarged. In Figure 5, the density of states (DOS) is plotted versus the frequencies (in THz). The DOS is calculated from the Fourier transformation of the velocity-autocorrelation functions. For better interpretation, the density of states from the simulations of the different indenter sizes are presented together with the “undisturbed” density of states. In comparison to the density of states of the reference structure, a significant broadening of the high-frequency peak can be observed, which is amplified with increasing indenter size. The less pronounced peaks of the reference density of states from the mid and higher frequency range are also influenced by the simulations with indenters, and a corresponding broadening of these peak structures occurs in the simulations with different indenter sizes. Figure 5. The vibrational densities of states Z( ν ) of silicon (left) and carbon (right) are shown for different indentation depths in 3C-SiC. For comparison, Z( ν ) from a simulation without indentation is shown. 3.2. 4H-SiC In the 4C-SiC simulation several indentation hardness tests are simulated. Table 2 lists the simulation parameters of the hemisphere penetration tests in the 4H-SiC crystal; the left column shows the size of the indenter in relation to the crystal height, in the middle column the respective simulation durations, and in the left column the indentation speed. Table 2. Simulation parameters for the indentation tests of 4H-SiC. Indenter size [%] Number of time steps Indentation speed [m/s] 5 15,000 9.73 10 25,000 11.67 15 60,000 7.30 20 100,000 5.84 25 150,000 4.87 30 200,000 4.34 The molecular dynamics simulations conducted on the 4C-SiC lattice reveal a clear transition from elastic to plastic deformation with increasing load. Figure 6 shows the final configurations of 4H-SiC after the application of indenters with sizes of 10%, 20%, and 30%, and simulation times of 15,000 fs, 100,000 fs, and 200,000 fs, respectively. In the results presented here, the indenters were applied in the x-direction. At lower penetration depths (indentation of 5% - 10%), the atomic arrangement remains largely crystalline, showing only minor distortions beneath the indenter. However, beginning at indentation of 20%, a fundamental structural transformation becomes evident. The atoms in the immediate contact region lose their long-range order, and a distinctly disordered zone emerges. This amorphous region expands further both laterally and in depth at higher loads (25% - 30%). Figure 6. The final configurations of 4H-SiC after the application of an indenter with a size of 10%, 20% and 30% and a simulation time of 15,000 fs, 100,000 fs and 200,000 fs respectively. The indenters were applied in the x-direction. Silicon atoms are represented in yellow spheres and carbon atoms in blue spheres. The energy profiles clearly show that low indentation depths cause only very minor changes. In contrast, larger indenters result in significant energy changes in the 4H-SiC system; these changes are particularly large when the sample is loaded with indenters of 25% from the x-side and indenters of 30% from the y- or z-side (see Figure 7). These observations are consistent with literature reports demonstrating that the deformation of SiC polytypes under nanoindentation is not governed solely by dislocation activity but also by the formation of amorphous zones once a critical load or indentation depth is exceeded. Zhu et al. (2019) [27] showed that nanoindentation of 4H-SiC induces amorphization along with phase transformation and dislocation generation. Similarly, Zhao et al. (2020) [28] reported that plastic deformation in 4H-SiC is dominated by the development of an amorphous region beneath the indenter, which persists after unloading. Wang et al. (2022) [29] further confirmed that 3C-, 4H-, and 6H-SiC all exhibit amorphization under indentation, alongside stacking faults and partial dislocations. From a mechanistic perspective, this transition occurs because the extremely high local stresses beneath the indenter surpass the energy barrier for conventional dislocation glide and stacking fault formation. Beyond this threshold, structural collapse into a metastable amorphous configuration becomes energetically more favorable than continued accumulation of crystalline defects. The progressive growth of the amorphous region with increasing indentation depth observed here is in excellent agreement with experimental nanoindentation studies on SiC, where transmission electron microscopy revealed extended amorphous/plastic zones in the contact region. Figure 7. The energy profiles of 4H-SiC at room temperature are plotted against the simulation time for different indentation depths each in x-, y- and z-direction from top to bottom, resp. As shown in Figure 8, the root mean squared displacement (RMSD) is a very good method for capturing the increasing displacements of atoms (due to dislocations) and the sliding of layers within structures. This behavior occurs particularly at greater penetration depths in the x and y directions. Figure 8. The root mean square displacements—RMSD in Å of 4H-SiC at room temperature are plotted against the simulation time for different indentation depths each in x- y- and z-direction from top to bottom, resp. To illustrate the structural changes caused by the indenter in 4H-SiC, we determine the partial pair correlation functions g( r ) for the C-Si, C-C, and Si-Si pairs both for the total structure (left part of the figure) and for the subset of atoms which are directly affected by the influence of the indenter (right side of the figure) (see Figure 9 (C-Si), Figure 10 (C-C), and Figure 11 (Si-Si), resp.). It is readily apparent that low indentation depths of 10% and 20% have virtually no effect on the pair correlation functions, as the observed interatomic distances correspond to those of the unstretched reference sample, with the first coordination peak largely preserved in position and intensity. This can be interpreted as indicating that the local bonding network of the sublattices is only slightly disturbed under moderate strain. At a penetration depth of 30%, a significant change in the C-C pair distribution functions is observed for all penetration directions (see Figure 10). Particularly striking is the splitting of the first coordination peak into two separate peaks, which are not present in the reference system. This splitting can be explained by a broadening of the nearest neighbor distance distribution: some atoms are compressed closer together by the applied compression, while others are displaced outwards. This behavior suggests the formation of new metastable bonding configurations, which can even lead to partial amorphization of the structure, consistent with local/partially local structural changes. Figure 11 shows the influence of the indenter on the Si-Si distances in the disturbed region; for all penetration depths, very short Si-Si distances are evident. The altered interatomic distances thus indicate plastic deformations and the subsequent formation of a (local) amorphous structure. The frequency spectra obtained from the simulations with different indenter sizes are shown in Figure 12 in comparison to the reference system. A significant broadening of the prominent high-frequency peak is evident, which, for indenter sizes above 20%, is accompanied by a slight shift towards higher frequencies. The splitting of the distances between the (nearest) C-C atoms (cf. Figure 10) and, to a somewhat lesser extent, also between the (nearest) Si-Si atoms (cf. Figure 11) leads, on the one hand, to a higher frequency, but on the other hand, also to a broadening due to the overall increased distances. Overall, these simulations emphasize that the deformation behavior of SiC cannot be explained by classical dislocation mechanisms alone. Instead, amorphization, defect accumulation, and local phase transformations play a central role in governing the response of SiC under nanoindentation. This highlights the complex, multi-mechanism nature of plasticity in SiC and aligns closely with recent molecular dynamics and experimental findings. 3.3. 6H-SiC The third crystal structure of SiC that is investigated in this work is 6H-SiC. As in the simulations shown above, various indentation tests were also carried out with the 6H-SiC structure. Table 3 specifies the parameters for the size of the indenter, the simulation duration during which the load is applied to the sample of a 6H-SiC crystal and the resulting indentation speed. Table 3. Simulation parameters for the indentation tests of 6H-SiC. Indenter size [%] Number of time steps Indentation speed [m/s] 5 15,000 9.73 10 25,000 11.67 15 50,000 8.76 20 50,000 11.67 25 75,000 9.73 30 100,000 8.76 Figure 9. The left panel shows the pair distribution functions g( r ) of C-Si pairs for 10%, 20% and 30% indentation depths in 4H-SiC, from top to bottom. For comparison, the distribution from a simulation without indentation is shown. On the right side the corresponding pair distribution functions for the subset of atoms most strongly influenced by the indenter are depicted, for better illustration, the curves have been enlarged by factors 20, 10 and 5, for 10%, 20% and 30% indentation, respectively. Figure 10. The left side shows the pair distribution functions g( r ) of C-C pairs for 10%, 20% and 30% indentation depths in 4H-SiC, from top to bottom, resp. For comparison, the distribution from a simulation without indentation is shown. The right panel depicts the corresponding pair distribution functions for the subset of atoms most strongly influenced by the indenter, for better illustration, the curves have been enlarged by factors 20, 10 and 5, for 10%, 20% and 30% indentation, respectively. First, we discuss the energy profiles exhibited by the system in various indenter tests, which are shown in Figure 13. In the y-direction, the equilibrium energy is slightly less negative compared to the x- and z-directions, indicating a lower initial stability of the lattice in this orientation. Across all directions, the indentation process leads first to an increase in the energy magnitude with simulation time and indenter size, most pronounced up to approximately 60,000 fs at 30% indenter size. This trend reflects the elastic displacement of atoms into new local configurations with higher potential energy. Beyond this point, a gradual decrease is observed, suggesting that atoms relax into energetically more favorable arrangements as the system approaches a metastable state. A special feature can be observed in the z-direction, where a sudden drop in energy occurs even at low penetration depths (see Figure 13). This discontinuity can be interpreted as the onset of plastic deformation, associated with dislocation nucleation or local amorphization. The subsequent stabilization of the curve close to the initial equilibrium energy indicates that the lattice along this axis accommodates deformation more effectively, retaining a configuration closer to its undisturbed state. The transition from elastic to plastic deformation in 6H-SiC has been reported to occur at characteristic indentation size, which coincide with abrupt changes in the load-displacement and energy-time curves. For example, Wu et al. (2023) observed that at a depth of approximately 3.8 nm, the system enters the elastic–plastic transition phase, accompanied by the formation of edge dislocations and partial amorphization [7]. With further indentation, the energy rises more steeply as dislocation rings expand and fracture occurs. This behavior is consistent with experimental observations of indentation-induced damage, where both brittle fracture and localized plasticity were detected depending on the applied load [30]. Figure 11. The left panel shows the pair distribution functions g( r ) of Si-Si pairs for 10%, 20% and 30% indentation depths in 4H-SiC, from top to bottom, resp. For comparison, the distribution from a simulation without indentation is shown. The right panel depicts the corresponding pair distribution functions for the subset of atoms most strongly influenced by the indenter, for better illustration, the curves have been enlarged by factors 20, 10 and 5, for 10%, 20% and 30% indentation, respectively. Figure 12. The vibrational densities of states Z( ν ) of silicon (left) and carbon (right) are shown for different indentation depths in 4H-SiC each in x-, y- and z-direction from top to bottom, resp. For comparison, Z( ν ) from a simulation without indentation is shown. Figure 14 shows the final configurations of 6H-SiC after the application of indenters with sizes of 10%, 20%, and 30%, and simulation times of 25,000 fs, 50,000 fs, and 100,000 fs, respectively. In the results presented here, the indenters were applied in the y-direction. This visualization of the structures with VMD shows a significant increase in bond breaks in the crystalline order and an order-disorder transition with growing amorphous regions, depending on the size of the indenter. The RMSD profiles exhibit a clear dependence on both indentation size and crystallographic orientation. For small penetration size (5% - 10%), the RMSD values remain relatively low in all directions, consistent with predominantly elastic deformation. With increasing size, the RMSD rises substantially, reflecting the onset of irreversible structural changes such as bond breaking, dislocation activity, and local amorphization. This trend is well documented in previous MD studies, where the elastic-plastic transition of 6H-SiC is associated with a steep increase in atomic displacements from their reference positions [6] [7]. A distinct feature is observed in the x-direction: at large indentation sizes, the RMSD first increases to a maximum (around 650 Å) and then decreases again, even though the indenter continues to penetrate deeper, see Figure 15. This non-monotonic behavior can be explained by atomic relaxation and structural reorganization. Part of the structural reorganization is due to the sliding of lattice planes caused by the applied pressure. Relaxations at the atomic level take place after the initial strong displacement of atoms; subsets of atoms rearrange into new, metastable configurations with lower average deviations from their local equilibrium positions. In other words, the system partially stabilizes by forming energetically favorable local structures despite ongoing indentation. Similar relaxation phenomena have been reported in MD simulations of SiC, where dislocation nucleation and subsequent reordering reduce the apparent RMSD after the onset of plasticity [7] [20]. By contrast, in the y- and z-directions the RMSD increases more steadily with the indenter size and shows less pronounced recovery. This indicates that the structural rearrangements in these orientations do not lead to equally effective re-stabilization, consistent with reports that anisotropy in 6H-SiC results in direction-dependent deformation pathways. Figure 13. The energy profiles of 6H-SiC at room temperature are plotted against the simulation time for different indentation depths each in x-, y- and z-direction from top to bottom, resp. Figure 14. The final configurations of 6H-SiC after the application of an indenter with a size of 10%, 20% and 30% and a simulation time of 25,000 fs, 50,000 fs and 100,000 fs respectively. The indenters were applied in the y-direction. Si atoms are represented as yellow spheres and C atoms as blue spheres. Figure 15. The root mean square displacements—RMSD in Å of 6H-SiC at room temperature are plotted against the simulation time for different indentation depths each in x- y- and z-direction from top to bottom, resp. To elucidate structural changes caused by the indentation of the test specimen in 6H-SiC, we determine the partial pair distribution functions g( r ) for the pairs C-Si, C-C and Si-Si, see Figure 16, Figure 17 (C-C) and Figure 18 (Si-Si), for both the total crystalline structure (left panel) and the subset of atoms from the mostly affected layers (right panel), resp. At shallow indenter sizes of 10% and 20%, the pair distribution functions remain close to the unindented reference, with the first coordination peak largely preserved in position and intensity. This indicates that the local bonding network of the sublattices is only marginally perturbed under moderate loading. At 30% indentation a significant modification of the C-C pair distribution functions is observed for all indentation directions, see Figure 17 (right part). Most notably, the first coordination peak splits into two distinct features, which are absent in the reference system. This splitting reflects a broadening of the distribution of nearest-neighbor distances. While some atoms are forced into closer proximity under compression, others are displaced outward. Such behavior suggests the formation of new metastable bonding configurations, consistent with localized structural rearrangements. The effect is particularly pronounced in the z-direction, where the pair distribution functions exhibit increased noise and additional minor peaks beyond the first coordination shell. These features indicate a more irregular arrangement of atoms, implying stronger displacement and a higher degree of amorphization. The emergence of such disorder highlights the loss of crystalline order at high indentation size and supports the conclusion that plastic deformation and local amorphization dominate the structural response in this regime. In addition, for the y- and z-directions a secondary peak emerges at a shorter distance, preceding the main coordination peak (cf. the graphs on the right side of Figure 17). This feature does not appear in the x-direction for 10% indenter size and indicates that, under loading along y and z, subsets of atoms adopt a new, energetically stable separation distance. The presence of this additional peak points to anisotropic deformation pathways and highlights that the carbon sublattice accommodates indentation differently depending on the crystallographic orientation. At 30% indentation size, all pair distribution functions (C-C, Si-C, and Si-Si) exhibit a pronounced anisotropy in the intensity of the first coordination peak. In the x-direction, the peak heights are consistently higher than in y and z, indicating that the nearest-neighbor distances remain more stable along this orientation. This suggests that the lattice in x-direction accommodates indentation with less disruption of the short-range order compared to the other crystallographic directions. A further distinction arises in the splitting behavior of the first peak. Clear splitting is observed in the C-C and Si-Si pair distribution functions, reflecting the emergence of two characteristic nearest-neighbor distances. This behavior points to bond bifurcation within the homonuclear sublattices, where some bonds are compressed while others are elongated under load. In contrast, the Si-C pair distribution functions show no comparable splitting but instead retains a broadened, unimodal first peak. This indicates that the heteronuclear Si-C bonds maintain a relatively more stable configuration, even with high indenter size. As the size of the indenter increases, the peaks broaden and secondary features appear, reflecting the nucleation of dislocations, the rearrangement of bonds and the progressive amorphization of the lattice [7]. Tian et al. (2019) further demonstrated that pair distribution function analysis can distinguish between different sublattices in SiC. Their simulations showed that the Si-C network remains comparatively stable under load, while homonuclear distances (C-C and Si-Si) are more susceptible to distortion. This highlights the role of the heteronuclear Si-C bonds as the structural backbone of hexagonal SiC, resisting bifurcation even under high stress [6]. The influence of the indentation on the dynamics of the system is illustrated by determining the corresponding frequency spectra in comparison to the reference system, see Figure 19. A clear broadening of the prominent high frequency peak is observed, which is accompanied with a slight shift to higher frequencies for indenter sizes larger than 20%. This finding is consistent with the experimental results. A small shift of the high-frequent peak due to nanoindentation is found experimentally by Rohbeck [31]. The indenter exerts the greatest influence on the dynamics of the carbon atoms. In the spectra of the C atom dynamics, particularly in the y and z directions, strong deviations from the reference spectrum can be observed for an indenter with a size of 30%. This manifests itself in the loss of the peak structure in the low and mid-frequency range and a broadening and shift in the high-frequency range. The dynamics of the C-C bond thus reflects the result we found for the corresponding pair distributions for 6H-SiC and points to the special role played by the C atoms in structural transformations, plastic deformation and amorphization in a heavily stressed 6H-SiC system. 4. Summary Molecular dynamics simulations were performed to investigate the indentation behavior of different SiC polymorphs (3C-, 4H-, and 6H-SiC). The simulations were carried out with varying indenter sizes (5% - 30% of the sample height) and simulation times of up to 200,000 fs at a constant temperature of 300 K. For 3C-SiC, loading was applied along the z-direction, whereas for 4H- and 6H-SiC the anisotropic response was additionally considered. In the simulations, energies and energy changes were monitored to gain information about the stability of the systems. Furthermore, partial pair correlation functions were evaluated to reveal changes in the structures relative to the reference structure. Finally, RMSD and frequency spectra were determined in the simulations to provide insights into the dynamic processes. For 3C-SiC, the pair distribution functions show no significant deviations up to an indenter size of 25%, indicating a predominantly elastic response. Only at 30% do slight changes appear, accompanied by rising RMSD values and a noticeable decrease in energy after approximately 25,000 fs. This marks the transition from elastic to plastic deformation, with larger indenters inducing an earlier response. Figure 16. The left panel shows the pair distribution functions g( r ) of C-Si pairs for 10%, 20% and 30% indentation depths in 6H-SiC, from top to bottom. For comparison, the distribution from a simulation without indentation is shown. On the right side the corresponding pair distribution functions for the subset of atoms most strongly influenced by the indenter are depicted, for better illustration, the curves have been enlarged by factors 20, 10 and 5, for 10%, 20% and 30% indentation, respectively. In 4H-SiC, the transition is more pronounced: at small indentation depths (5% - 10%) the crystal lattice remains largely intact, while at around 20% indenter size a locally amorphous zone forms beneath the indenter. This zone expands laterally and in depth with increasing load. These observations are consistent with literature reports that, in addition to dislocations, the formation of amorphous regions is the dominant plastic mechanism in hexagonal SiC. Figure 17. The left side shows the pair distribution functions g( r ) of C-C pairs for 10%, 20% and 30% indentation depths in 6H-SiC, from top to bottom, resp. For comparison, the distribution from a simulation without indentation is shown. The right panel depicts the corresponding pair distribution functions for the subset of atoms most strongly influenced by the indenter, for better illustration, the curves have been enlarged by factors 20, 10 and 5, for 10%, 20% and 30% indentation, respectively. For 6H-SiC, the energy profiles and structural measures reveal a strong directional dependence. Along the x-direction, partial relaxation and re-stabilization are observed after initial displacement, whereas in the y- and z-directions RMSD values and amorphization increases steadily. The pair distribution functions further show that homonuclear bonds (CC, SiSi) are more strongly distorted and may split, while heteronuclear SiC bonds remain comparatively stable. This indicates that the plastic response of 6H-SiC is governed by anisotropic deformation pathways and the stabilizing role of SiC bonds. Figure 18. The left panel shows the pair distribution functions g( r ) of Si-Si pairs for 10%, 20% and 30% indentation depths in 6H-SiC, from top to bottom, resp. For comparison, the distribution from a simulation without indentation is shown. The right panel depicts the corresponding pair distribution functions for the subset of atoms most strongly influenced by the indenter, for better illustration, the curves have been enlarged by factors 20, 10 and 5, for 10%, 20% and 30% indentation, respectively. In summary, the simulations have shown that the SiC configurations in the 3C, 4H and 6H SiC structure prove to be very hard, very stable and elastic, since only relatively large intruders cause plastic deformations, which in extreme cases lead to the formation of amorphous regions. As a consequence for practical applications, it must be taken into account that the anisotropic deformations observed in the simulations and the onset of amorphization can cause changes in the conductivity of the hexagonal polytypes and influence or modify the transport mechanisms in SiC-based electronic devices or mechanical components.