The polymer consistent force field (PCFF), one of the consistent family force field developed for synthetic polymers was used to define the atomic interactions in the vitrinite substrate. This force field is widely used for simulating organic polymers and is sufficiently accurate to describe the coal molecular system, which consists of C, H, O, N, and S atoms (Meng et al., 2020). In the PCFF force field, the total energy is composed of valence, cross, and nonbonded energies. The 6/9 Lennard-Jones (LJ) potential was selected to describe the pair interactions between vitrinite atoms, which may be expressed as follows (Sun, 1994): E V D W = ε i j 0 [ 2 ( r i j 0 r i j ) 9 − 3 ( r i j 0 r i j ) 6 ] ( r < r c ) (1) r i j 0 = [ ( r i 0 ) 6 + ( r j 0 ) 6 2 ] 1 6 (2) ε i j 0 = 2 ε i 0 ⋅ ε j 0 [ ( r i 0 ) 3 ⋅ ( r j 0 ) 3 ( r i 0 ) 6 + ( r j 0 ) 6 ] (3) where E _VDW is the nonbonded van der Waals interactions between atoms i and j, r i j 0 is the equilibrium distance between atoms i and j, ε i j 0 is the coefficient of well-depth energy, r _c is the global cutoff distance for LJ interactions, which was set to 10 Å in this study, r _ij is the instantaneous distance between atoms i and j, r i 0 and r j 0 refer to the equilibrium distances between atoms i and atoms j, respectively, and ε i 0 and ε j 0 refer to the well-depth energies of atoms i and j, respectively.

A ghost spherical indenter with a radius of 15 Å was used in the simulations (Figure 2). The ghost indenter repels all atoms in the vitrinite substrate, and the repulsive force may be expressed as follows: F ( r ) = − K ( r − R ) 2 (4) where K is the specified force constant (which was set to 10 in the simulations), r is the distance between vitrinite atoms and the center of the ghost indenter, and R is the radius of the ghost indenter (= 15 Å).

The bottom layer, with a height of 10 Å, was fixed as a rigid body to avoid vertical translation of the vitrinite substrate (Figure 2). Atoms in the Newtonian layer were put into a microcanonical (NVE) ensemble, in which the particles are allowed to move freely as the ghost indenter is pushed into substrate. Atomic motion in the Newtonian layer strictly follows Newton’s laws of motion. The loading and unloading rate of the ghost indenter along the Z-direction was set to 200 m/s, and the whole indentation procedure lasted for 7,500 fs. The maximum indentation depth was set to 1.5 nm.

The plastic deformation of metallic materials mainly occurs through the dislocation of crystalline structures, while that of amorphous materials such as vitrinite occurs via shear transformation, which refers to the rearrangement of a pack of atoms overcoming the energy barrier (Greer et al., 2013). The rearranged atoms usually normally have very large displacements compared to the surrounding atoms, and the local domain containing rearranged atoms called the shear transformation zone (STZ) (Hua et al., 2020).

The atomic local shear strain, η i M i s e s , of each atom i was introduced to quantify the atomic shear deformation of the vitrinite substrate (Shimizu et al., 2007). Figure 3 shows the shear strain distributions of the vitrinite atoms at different indenter loading stages. The cross-sectional strain distributions reveal that, as the indentation depth increases, the shear strain of the substrate atoms under the ghost indenter rapidly increases, propagating from the contact region to the outside of the indenter. In addition, the overhead view of the vitrinite substrate indicates that the shear strain of atoms expands rapidly from the center to periphery as the indentation depth increases. The proliferation of the STZ is mainly concentrated directly below the ghost indenter, indicating that the shear strain spreads directly into the interior of the vitrinite substrate. Essentially, the proliferation of the STZ represents the spreading and rearrangement of vitrinite atoms, which allows the stress-strain energy produced from elastic loading to be released.

To determine the motion of atoms in the vitrinite substrate, the atomic displacement vectors in a vertical slice with a width of 10 Å were calculated using the OVITO software package (Version 3.4.4) (Figure 4). The green arrows in Figure 4 show the directions of vitrinite atom motion, and the length of the green line represents the displacement of the atom. The results show that as the loading depth of the ghost indenter increases, the displacement vectors of vitrinite atoms around the indenter gradually lengthen under the effect of the repulsive forces between the substrate and the indenter (Figures 4A–D). Further, the magnitudes and directions of the displacement vectors vary with location. At region B, atoms directly below the ghost indenter move vertically downwards during the indentation process, which indicates that the substrate atoms in this region are mainly subjected to the force along the Z-direction (f _Z ). On the other hand, the substrate atoms at regions A and C move downwards at an angle to the vertical, indicating that the atoms in these regions are subjected to the forces along both the X- and Z-directions (f _X and f _Z ). The displacements of atoms at region B are also much larger than those at A and C, which illustrates significant variation in the magnitudes of the forces acting on atoms at different locations.

Given the variations in displacement seen in Figure 4, the local substrate atoms around ghost indenter can be roughly subdivided into six distinct regions (Figure 5), with the atoms in each region moving in more or less the same direction, but in a different direction to those in the other regions. Molecular chains or atoms of neighboring regions slip at the region boundary, leading to the formation of shear bands. Figure 5 shows that the displacement of vitrinite atoms in regions I, II, and III is generally much larger than that in other regions (Lin et al., 2021). This is likely due to the fact that, in the initial indentation stage, atoms in these regions are the first to be subjected to repulsive forces and start moving. These atoms then squeeze the ones in regions IV, V, and IV, causing those peripheral atoms to move in the following indentation stage. However, the directions of motion of the squeezed atoms diverge due to the heterogenous molecular structures and inner stress differentiation within the substrate, resulting in slippage at the region boundaries, viz. the shear bands.

The response of coal molecular structures to tectonic stress-induced deformation, which is similar to the molecular variation during indenter loading process, has been well-researched (Cao et al., 2007; Pan et al., 2015; Song et al., 2018). In terms of the plastic deformation of coal induced by tectonic stress, the motion of molecular chains is the most significant deformation mechanism (Liu and Jiang, 2019), which also leads to the formation of STZs in TDCs. During the indentation process, the adjustment of molecular chains is another important deformation mechanism of the vitrinite substrate. At the initial loading stage, the emergence of the STZ in vitrinite substrate may be initiated in the “weak spots” (atoms with larger free volume) in the vitrinite substrate (Wang et al., 2010). Under the effect of shear stress, the STZ begins to propagate via the adjustment of vitrinite molecular chains. The motion of vitrinite atoms in distinct regions (Figure 5) indicated that the adjustment of molecular chains is due to factors such as sliding between aromatic layers, reorientation of aromatic lamella, and bending of aromatic chains (Figure 6), which is similar to the molecular variations seen for the plastic deformation of coal (Ju et al., 2014; Han et al., 2017).

Moreover, the development of structural defects, which mainly includes various types of point defects, is promoted in TDCs (Pan et al., 2017; Song et al., 2019). Similarly, for the indentation of vitrinite substrate, the STZ also appears be a kind of point defect, with a greater degree of structural disorder and a larger free volume (Greer et al., 2013). Apart from the development of defects, nanocrystallization also appears in the shear bands during nanoindentation (Kim et al., 2002), which is similar to the aromatization of ductile deformed coal (Hou et al., 2012). During the vitrinite indentation process, sliding between the aromatic layers can enhance the degree of stacking of the vitrinite molecular structures (Figure 6A); hence, it is evident that the atomic evolution of the vitrinite structure is consistent with that of TDCs.

A 10 Å slice extracted from the middle of the simulation box (as shown in Figure 4A) was used to identify the surface morphology of the vitrinite substrate (Figures 7A–C). Following penetration by the indenter, different substrates tend to show one of there will be two main modes of indentation behavior, either sink-in or pile-up modes (Hou et al., 2020) (as illustrated in Figures 7D,E). In this case, Figure 7 shows that the surface morphology of vitrinite atoms after penetration is inclined towards the sink-in mode, which is consistent with previous experiments (Hou et al., 2020). The surface morphology of substrate materials is actually determined by the density of substrate atoms (Hua et al., 2020), with a higher atomic density of substrate material likely to lead to greater atomic pile-up around the indenter. Vitrinite is of lower atomic density and has a loose molecular structure relative to metals and other materials. Consequently, there is enough free space for vitrinite atoms to move and for molecular chains to adjust, which leads to the sink-in type surface.

Figure 8 shows the mechanical properties of vitrinite obtained from the simulations. As the contact depth increases, the repulsive forces between the indenter and the vitrinite atoms are enhanced. In the loading process, the load-depth curve (Figure 8) is serrated, which may be due to the shear transformation of vitrinite substrate (Schuh and Nieh, 2003). Serrations on the load-depth curve correspond to sudden load drops of vitrinite indentation. For each serration, at the initial elastic loading stage, vitrinite atoms are not able to overcome the motion energy barrier; therefore, the loading force continues to increase. Once the generated strain energy crosses the energy barrier, a pile of atoms under the indenter is able to move freely, giving rise to the formation of the STZ. As the loading displacement is constant rate, the activation of shear transformation leads directly to a decreasing in the load drop, which corresponds to fluctuation in the load-depth curve. In the unloading process, a sudden drop of load appears at the initial unloading stage (Figure 8), which is caused by the high unloading velocity (Zhu and Fang, 2012). The repulsive force between the indenter and substrate atoms disappears once the distance between substrate atoms and the center of spherical indenter is larger than the indenter radius. Therefore, the unloading of indenter with high velocity directly leads to the separation between substrate atoms and indenter, which results in the sudden drop of load.

The hardness, H, which is a commonly-used parameter to characterize the resistant ability to allow hard objects press into (Zhu and Fang, 2012), is defined as the ratio of load force to the projected contact area when the indentation depth is maximum, as follows: H = F A c (5) where F is the load of ghost indenter, and A _c is the projected contact area of the ghost indenter embedded into the substrate material. For a spherical ghost indenter, A _c is calculated as follows (Qiu et al., 2014): A c = π ( 2 R − h ) h (6) where R is the radius of ghost indenter and h is the indentation depth. According to Eqs. 5, 6, the simulated hardness of coal vitrinite at indentation depth of 1.5 nm is calculated to be 1.31 GPa. Experiments have shown that the vitrinite hardness in coal with similar rank was in the range of 0.52–0.72 GPa (Hou et al., 2020), which is a little lower than the simulated value from the present study. For nano-scale simulations, the simulated hardness of vitrinite be influenced by factors such as the indentation size effects, temperature, indenter radius, and loading rate (Luu et al., 2021; Valencia et al., 2021). However, as there is basically no existing literature regarding the influence of the aforementioned factors on coal vitrinite hardness, it is hard to clarify exactly which of these factors led to the overestimation of the hardness in the current simulations. This issue will be examined in more detail in a future study.