An Investigation into the Densification-Affected Deformation and Fracture in Fused Silica under Contact Sliding

Subsurface damage of fused silica optics is one of the major factors restricting the performance of optical systems. The densification-affected deformation and fracture in fused silica under a sliding contact are investigated in this study, via three-dimensional finite element analysis (FEA). The finite element models of scratching with 70.3° conical and Berkovich indenters are established. A refined elliptical constitutive model is used to consider the influence of densification. The finite element models are experimentally verified by elastic recovery, and theoretically verified by hardness ratio. Results of densification and plastic deformation distributions indicate that the accuracy of existent sliding stress field models may be improved if the spherical/cylindrical yield region is replaced by an ellipsoid/cylindroid, and the embedding of the yield region is considered. The initiation sequence, and the locations and stages of radial, median, and lateral cracks are discussed by analyzing the predicted sliding stress fields. Median and radial cracks along the sliding direction tend to be the first cracks that emerge in the sliding and unloading stages, respectively. They coalesce to form a big median–radial crack that penetrates through the entire yield region. The fracture behavior of fused silica revealed in this study is essential in the low-damage machining of fused silica optics.


Introduction
Fused silica, or silica-rich glass optics, are widely used in laser nuclear fusion devices [1], large astronomical telescopes [2], semiconductor technology [3,4], and consumer electronics. Subsurface damage has plenty of negative effects on the performance of optics, e.g., increasing optical scatter, reducing mechanical strength, and increasing laser-induced damage (LID). For instance, the subsurface damage is one of the precursors resulting in LID. The LID of fused silica optics is one of the key factors restricting the output power and a key challenge for the long-term and stable operation of high-power laser facilities [5]. Therefore, an in-depth understanding of the material removal and damage formation mechanisms of fused silica subjected to machining is essential to fabricate damage-free optics.

Scratching Tests
As shown in Figure 1, fused silica samples (Corning UV 7980, Corning Corp., Corning, NY, USA) were scratched by an edge-forward Berkovich tip on a nanoindentation machine (TI-950 TriboIndenter, Hysitron Inc., Eden Prairie, MN, USA). The samples were carefully polished to achieve a surface roughness small than 2 nm. Scratching tests were performed under constant normal loads of 1 mN, 2 mN, 4 mN, 200 mN, 400 mN, 600 mN, 1 N, and 1.2 N. The sliding length was 250 µm, which is significantly greater than the scratching depth. The scratching process consists of the approaching stage A , the preliminary profiling stage B (to obtain the original surface profile), the indentation stage C , the scratching stage D , the unloading stage E , and the postmortem profiling stage F (to obtain the residual surface profile). The variations of normal load, normal displacement, lateral load, and lateral displacement with time were recorded during scratching.
Micromachines 2022, 13, x FOR PEER REVIEW 3 of 17 and scratching hardness is investigated, and the cracking behavior of fused silica under scratching explored.

Scratching Tests
As shown in Figure 1, fused silica samples (Corning UV 7980, Corning Corp., Corning, NY, USA) were scratched by an edge-forward Berkovich tip on a nanoindentation machine (TI-950 TriboIndenter, Hysitron Inc., Eden Prairie, MN, USA). The samples were carefully polished to achieve a surface roughness small than 2 nm. Scratching tests were performed under constant normal loads of 1 mN, 2 mN, 4 mN, 200 mN, 400 mN, 600 mN, 1 N, and 1.2 N. The sliding length was 250 μm, which is significantly greater than the scratching depth. The scratching process consists of the approaching stage Ⓐ , the preliminary profiling stage Ⓐ (to obtain the original surface profile), the indentation stage Ⓐ , the scratching stage Ⓐ , the unloading stage Ⓐ , and the postmortem profiling stage Ⓕ (to obtain the residual surface profile). The variations of normal load, normal displacement, lateral load, and lateral displacement with time were recorded during scratching. After scratching, the samples were measured by an atomic force microscope (AFM) (Innova, Veeco, Plainview, NY, USA) to obtain the three-dimensional topography of the impression. After etching by the buffered HF solution to open the surface cracks [31], the morphology of the cracks was characterized by a scanning electron microscope (SEM) (SU3500, Hitachi, Japan).

Finite Element Modeling
The finite element analysis of scratching with an edge-forward Berkovich indenter and a conical indenter was performed on a commercial finite element code ABAQUS. A modified elliptical constitutive model [28,30] was used to consider the influence of densification on the deformation in fused silica: After scratching, the samples were measured by an atomic force microscope (AFM) (Innova, Veeco, Plainview, NY, USA) to obtain the three-dimensional topography of the impression. After etching by the buffered HF solution to open the surface cracks [31], the morphology of the cracks was characterized by a scanning electron microscope (SEM) (SU3500, Hitachi, Japan).

Finite Element Modeling
The finite element analysis of scratching with an edge-forward Berkovich indenter and a conical indenter was performed on a commercial finite element code ABAQUS. A modified elliptical constitutive model [28,30] was used to consider the influence of densification on the deformation in fused silica: where q is equivalent shear stress; p is hydrostatic pressure; d is the von Mises yield stress under pure shear; and p b is the hydrostatic yield stress for pure compression. The relationship between hydrostatic pressure p and the densification ζ is modeled by: where ζ max (%) is the saturated densification under compression, and p 0 (GPa) is the hydrostatic pressure under which a densification of ζ max /2 is produced. The parameters of the modified elliptical model used in this study are taken from the ref. [28].
In the finite element model, an infinitely sharp Berkovich edge-forward indenter slides along the x-axis on the top surface of a deformable parallelepiped with a dimension of W × W × (l + 2W), as shown in Figure 2. The diamond indenter is assumed to be rigid because its Young's modulus [32] and hardness [33] are much higher than those of the fused silica samples [14]. An eight-node linear brick element with reduced integration and hourglass control is used. Refined FE mesh with an element size of l e is used in a parallelepiped with a dimension of a × a × (l + 2W), and graded FE mesh is used in the residual region. For conical scratching, the semi-included angle α of the conical indenter is set as 70.3 • , to ensure that the projected area-to-indentation depth function is the same as the commonly used Vickers and Berkovich indenters. A = 2.79 h max is the nominal contact radius for 70.3 • conical scratching.
where (%) is the saturated densification under compression, and p0 (GPa hydrostatic pressure under which a densification of /2 is produced. The par of the modified elliptical model used in this study are taken from the ref. [28]. In the finite element model, an infinitely sharp Berkovich edge-forward i slides along the x-axis on the top surface of a deformable parallelepiped with a dim of W × W × (l + 2W), as shown in Figure 2. The diamond indenter is assumed to because its Young's modulus [32] and hardness [33] are much higher than thos fused silica samples [14]. An eight-node linear brick element with reduced integrat hourglass control is used. Refined FE mesh with an element size of le is us parallelepiped with a dimension of a × a × (l + 2W), and graded FE mesh is used residual region. For conical scratching, the semi-included angle α of the conical i is set as 70.3°, to ensure that the projected area-to-indentation depth function is th as the commonly used Vickers and Berkovich indenters. A = 2.79 hmax is the n contact radius for 70.3° conical scratching. As shown in Figures 2 and 3, the sliding process is divided into three stages, indentation stage ①, the sliding stage ②, and the unloading stage ③. The scr depth, length, and speed are denoted as hmax, l, and v, respectively. The Coulomb model is used to model the adhesion friction behavior between the indenter sample. The coefficient of adhesion friction f was determined to be 0.04, by com FEA and scratching tests. As shown in Figures 2 and 3, the sliding process is divided into three stages, i.e., the indentation stage 1 , the sliding stage 2 , and the unloading stage 3 . The scratching depth, length, and speed are denoted as h max , l, and v, respectively. The Coulomb friction model is used to model the adhesion friction behavior between the indenter and the sample. The coefficient of adhesion friction f was determined to be 0.04, by comparing FEA and scratching tests. relationship between hydrostatic pressure p and the densification  is modeled by: where (%) is the saturated densification under compression, and p0 (GPa) is the hydrostatic pressure under which a densification of /2 is produced. The parameters of the modified elliptical model used in this study are taken from the ref. [28].
In the finite element model, an infinitely sharp Berkovich edge-forward indenter slides along the x-axis on the top surface of a deformable parallelepiped with a dimension of W × W × (l + 2W), as shown in Figure 2. The diamond indenter is assumed to be rigid because its Young's modulus [32] and hardness [33] are much higher than those of the fused silica samples [14]. An eight-node linear brick element with reduced integration and hourglass control is used. Refined FE mesh with an element size of le is used in a parallelepiped with a dimension of a × a × (l + 2W), and graded FE mesh is used in the residual region. For conical scratching, the semi-included angle α of the conical indenter is set as 70.3°, to ensure that the projected area-to-indentation depth function is the same as the commonly used Vickers and Berkovich indenters. A = 2.79 hmax is the nominal contact radius for 70.3° conical scratching. As shown in Figures 2 and 3, the sliding process is divided into three stages, i.e., the indentation stage ①, the sliding stage ②, and the unloading stage ③. The scratching depth, length, and speed are denoted as hmax, l, and v, respectively. The Coulomb friction model is used to model the adhesion friction behavior between the indenter and the sample. The coefficient of adhesion friction f was determined to be 0.04, by comparing FEA and scratching tests.  The single-variable method is used to optimize the cross-section dimension W, sliding length l, and element size l e . The appropriate parameters result in stable and convergent normal and tangential loads, and apparent coefficient of friction in the sliding stage. h max is assumed to be 1 µm. Results show that a cross-section dimension of 5a × 5a, a sliding length of 10 h max , and a mesh size of 1/8 h max are appropriate for the simulations.

Experimental Verification of Elastic Recovery
The elastic recovery ratio f e for scratching reflects the extent of elastic deformation relative to the whole deformation. In addition, f e can be conveniently measured by AFM. Therefore, f e is used to verify the finite element model in this study.
As shown in Figure 4, the leading end of the impression induced by scratching with an edge-forward Berkovich indenter is measured by AFM to obtain its three-dimensional topography. Pile-up is obvious on the lateral sides of the impression. Figure 4 also indicates that the residual depth h f (see Figure 5) slightly decreases with the distance d to the unloading position of indenter tip. The profile shown in Figure 5 is obtained by averaging five equally spaced cross-section profiles of the middle part of the impression. The residual scratch depth h f after elastic recovery, determined from Figure 5, is 668 nm. The single-variable method is used to optimize the cross-section dime sliding length l, and element size le. The appropriate parameters result in s convergent normal and tangential loads, and apparent coefficient of friction in t stage. hmax is assumed to be 1 μm. Results show that a cross-section dimension a sliding length of 10 hmax, and a mesh size of 1/8 hmax are appropriate for the sim

Experimental Verification of Elastic Recovery
The elastic recovery ratio f e for scratching reflects the extent of elastic de relative to the whole deformation. In addition, f e can be conveniently measured Therefore, f e is used to verify the finite element model in this study.
As shown in Figure 4, the leading end of the impression induced by scratc an edge-forward Berkovich indenter is measured by AFM to obtain its three-dim topography. Pile-up is obvious on the lateral sides of the impression. Figu indicates that the residual depth hf (see Figure 5) slightly decreases with the dis the unloading position of indenter tip. The profile shown in Figure 5 is ob averaging five equally spaced cross-section profiles of the middle part of the im The residual scratch depth hf after elastic recovery, determined from Figure 5, is  In order to determine the scratching depth directly from the displacement normal displacement of the indenter (see Figure 1c) during the scratching pro  The single-variable method is used to optimize the cross-section dimension W sliding length l, and element size le. The appropriate parameters result in stable and convergent normal and tangential loads, and apparent coefficient of friction in the sliding stage. hmax is assumed to be 1 μm. Results show that a cross-section dimension of 5a × 5a a sliding length of 10 hmax, and a mesh size of 1/8 hmax are appropriate for the simulations.

Experimental Verification of Elastic Recovery
The elastic recovery ratio f e for scratching reflects the extent of elastic deformation relative to the whole deformation. In addition, f e can be conveniently measured by AFM Therefore, f e is used to verify the finite element model in this study.
As shown in Figure 4, the leading end of the impression induced by scratching with an edge-forward Berkovich indenter is measured by AFM to obtain its three-dimensiona topography. Pile-up is obvious on the lateral sides of the impression. Figure 4 also indicates that the residual depth hf (see Figure 5) slightly decreases with the distance d to the unloading position of indenter tip. The profile shown in Figure 5 is obtained by averaging five equally spaced cross-section profiles of the middle part of the impression The residual scratch depth hf after elastic recovery, determined from Figure 5, is 668 nm.  In order to determine the scratching depth directly from the displacement curve, the normal displacement of the indenter (see Figure 1c) during the scratching process was In order to determine the scratching depth directly from the displacement curve, the normal displacement of the indenter (see Figure 1c) during the scratching process was corrected by the original profile of the sample surface. The corrected normal displacement in the scratching stage D was calculated by subtracting the uncorrected displacement from t 2 to t 1 in the stage B from the uncorrected displacement from t 3 to t 4 . Similarly, the corrected normal displacement in the postmortem profiling stage F was calculated by subtracting the uncorrected displacement from t 1 to t 2 from the uncorrected displacement from t 5 to t 6 . The evolution of the corrected normal displacement with time is shown in Figure 6. The maximum scratching depth and residual depth are 1063 nm and 479 nm, respectively. It is worth noting that this value of residual depth is smaller than that measure by AFM (i.e., h f = 668 nm). This is possibly because the indenter did not strictly follow the scratching path in stage F , due to the movement of the sample or the motion error of the indentation test in the lateral direction. By contrast, the AFM probe accurately detects the lowest positions of the residual scratching profiles for two reasons. First, the tip radius of the AFM probe is much smaller than that of the Berkovich indenter. Second, the AFM probe is scanning across the impression. Using the AFM-measured h f , the elastic recovery ratio is calculated to be f e = 1−h f /h amx = 37.2%. Ba analyzing the FEA-predicted profiles of the scratching impression at the fully-loaded and fully-unloaded states shown in Figure 7, the value of f e predicted by FEA is 37%, which is very close to the experimental value.
Micromachines 2022, 13, x FOR PEER REVIEW 6 of 17 corrected by the original profile of the sample surface. The corrected normal displacement in the scratching stage Ⓓ was calculated by subtracting the uncorrected displacement from t2 to t1 in the stage Ⓑ from the uncorrected displacement from t3 to t4. Similarly, the corrected normal displacement in the postmortem profiling stage Ⓕ was calculated by subtracting the uncorrected displacement from t1 to t2 from the uncorrected displacement from t5 to t6. The evolution of the corrected normal displacement with time is shown in Figure 6. The maximum scratching depth and residual depth are 1063 nm and 479 nm, respectively. It is worth noting that this value of residual depth is smaller than that measure by AFM (i.e., hf = 668 nm). This is possibly because the indenter did not strictly follow the scratching path in stage Ⓕ, due to the movement of the sample or the motion error of the indentation test in the lateral direction. By contrast, the AFM probe accurately detects the lowest positions of the residual scratching profiles for two reasons. First, the tip radius of the AFM probe is much smaller than that of the Berkovich indenter. Second, the AFM probe is scanning across the impression. Using the AFM-measured hf, the elastic recovery ratio is calculated to be fe = 1−hf/hamx = 37.2%. Ba analyzing the FEA-predicted profiles of the scratching impression at the fully-loaded and fully-unloaded states shown in Figure 7, the value of fe predicted by FEA is 37%, which is very close to the experimental value.   corrected by the original profile of the sample surface. The corrected normal displacement in the scratching stage Ⓓ was calculated by subtracting the uncorrected displacement from t2 to t1 in the stage Ⓑ from the uncorrected displacement from t3 to t4. Similarly, the corrected normal displacement in the postmortem profiling stage Ⓕ was calculated by subtracting the uncorrected displacement from t1 to t2 from the uncorrected displacement from t5 to t6. The evolution of the corrected normal displacement with time is shown in Figure 6. The maximum scratching depth and residual depth are 1063 nm and 479 nm, respectively. It is worth noting that this value of residual depth is smaller than that measure by AFM (i.e., hf = 668 nm). This is possibly because the indenter did not strictly follow the scratching path in stage Ⓕ, due to the movement of the sample or the motion error of the indentation test in the lateral direction. By contrast, the AFM probe accurately detects the lowest positions of the residual scratching profiles for two reasons. First, the tip radius of the AFM probe is much smaller than that of the Berkovich indenter. Second, the AFM probe is scanning across the impression. Using the AFM-measured hf, the elastic recovery ratio is calculated to be fe = 1−hf/hamx = 37.2%. Ba analyzing the FEA-predicted profiles of the scratching impression at the fully-loaded and fully-unloaded states shown in Figure 7, the value of fe predicted by FEA is 37%, which is very close to the experimental value.

Theoretical Verification of Hardness Ratio
As shown in Figure 8, FEA predicts that the hardness ratio k

Theoretical Verification of Hardness Ratio
As shown in Figure 8, FEA predicts that the hardness ratio p is slightly larger than utility, and nearly independent of the sharpness of the indenter, where p T H and p s H are the ploughing hardness along the sliding and vertical directions, respectively. This is consistent with theoretical analysis, as detailed below. For an infinitesimal contact area dA, the contact force is normal to dA. Therefore, the forces along the sliding and vertical directions, i.e., dFT and dFN, have the following relationship: where Apl and Apv are the laterally and vertically projected contact areas, respectively; p(h,β) is the contact pressure at the point (h,β) on the indenter surface; and β and h are the phase and the height measured from the indenter tip, respectively. According to the definition of hardness, the ploughing hardness can be expressed by: The contours of the contact stress induced by conical scratching and Berkovich scratching resemble concentric circles and triangles, respectively [30]. Therefore, p(h,β) is nearly independent of β. As the geometries of the conical and Berkovich indenters are selfsimilar, the vertically and laterally projected areas of the contact zone with a height from h to h + dh, i.e., , ∼ and , ∼ have the following relationship: Therefore, the hardness ratio: For an infinitesimal contact area dA, the contact force is normal to dA. Therefore, the forces along the sliding and vertical directions, i.e., dF T and dF N , have the following relationship: where A pl and A pv are the laterally and vertically projected contact areas, respectively; p(h,β) is the contact pressure at the point (h,β) on the indenter surface; and β and h are the phase and the height measured from the indenter tip, respectively. According to the definition of hardness, the ploughing hardness can be expressed by: and The contours of the contact stress induced by conical scratching and Berkovich scratching resemble concentric circles and triangles, respectively [30]. Therefore, p(h,β) is nearly independent of β. As the geometries of the conical and Berkovich indenters are self-similar, the vertically and laterally projected areas of the contact zone with a height from h to h + dh, i.e., dA pl,h∼h+∆h and dA pv,h∼h+∆h have the following relationship: Therefore, the hardness ratio: By combining Equations (3), (12) and (18) in ref. [30], the following expression can be obtained for a conical indenter: where k H = H T /H s is the ratio of tangential hardness and scratching hardness; and µ 0 is the friction coefficient induced by ploughing. As k p H ≈ 1, we can conclude from Equation (8) that k H > 1 when friction exists, i.e., H T > H s . This is consistent with the scratching tests [34]. It is worth noting that the above analysis considers the non-uniform distribution of the contact pressure. This is more accurate than the widely adopted assumption that the contact pressure is uniformly distributed [35].

Scratching Hardness
Scratching hardness is widely used to model the scratching load that is a key factor determining the fracture behavior. Although it is reported that friction only plays a small role in indentation hardness for blunt indenters [36], Figure 9 shows that the indentation hardness H i (the hardness at the end of stage 1 ) for edge-leading Berkovich scratching is slightly increased with the rise in friction. By contrast, the scratching hardness (the hardness in the right red box) is nearly independent of friction.  (3), (12) and (18) in ref. [30], the following expression obtained for a conical indenter: where kH = HT/Hs is the ratio of tangential hardness and scratching hardness; and 0  friction coefficient induced by ploughing. As ≈ 1, we can conclude from Equat that kH > 1 when friction exists, i.e., HT > Hs. This is consistent with the scratching tes It is worth noting that the above analysis considers the non-uniform distribution contact pressure. This is more accurate than the widely adopted assumption th contact pressure is uniformly distributed [35].

Scratching Hardness
Scratching hardness is widely used to model the scratching load that is a key determining the fracture behavior. Although it is reported that friction only plays a role in indentation hardness for blunt indenters [36], Figure 9 shows that the inden hardness Hi (the hardness at the end of stage ①) for edge-leading Berkovich scratc slightly increased with the rise in friction. By contrast, the scratching hardne hardness in the right red box) is nearly independent of friction. In order to understand the above phenomena, the evolutions of the contact ar normal force with time are plotted in Figure 10a,b, respectively. During the tra from indentation to sliding, i.e., in the stage ②-1, the contact area is considerably re because the support for the indenter rear is removed. If no elastic recovery occurs, t face of the indenter is entirely separated from the sample surface, and the conta decreases to two-thirds. However, Figure 10a shows that the contact area duri steady sliding stage (in the right red box) is significantly bigger than two-thirds during indentation (in the left red box). This is due to elastic recovery. In order to understand the above phenomena, the evolutions of the contact area and normal force with time are plotted in Figure 10a,b, respectively. During the transition from indentation to sliding, i.e., in the stage 2 -1, the contact area is considerably reduced because the support for the indenter rear is removed. If no elastic recovery occurs, the rear face of the indenter is entirely separated from the sample surface, and the contact area decreases to two-thirds. However, Figure 10a shows that the contact area during the steady sliding stage (in the right red box) is significantly bigger than two-thirds of that during indentation (in the left red box). This is due to elastic recovery.
The moving direction of the indenter relative to the sample determines the influence of friction on the contact area. The contact area during static indentation is slightly decreased by friction, because the friction-induced downward shear stress applied to the sample surface results in a decrease in the contact area, as shown in Figure 10a. By contrast, friction leads to an increase in the contact area during the steady sliding stage because the frictioninduced shear stress on the sample surface is upward. The indentation hardness increases with friction, while the contact area decreases with friction. Therefore, the normal force during indentation is nearly independent of friction, as shown in Figure 10b. During the steady sliding stage, the scratching hardness is nearly independent of friction because the friction increases both the contact area and the normal force. The moving direction of the indenter relative to the sample determines the influence of friction on the contact area. The contact area during static indentation is slightly decreased by friction, because the friction-induced downward shear stress applied to the sample surface results in a decrease in the contact area, as shown in Figure 10a. By contrast, friction leads to an increase in the contact area during the steady sliding stage because the friction-induced shear stress on the sample surface is upward. The indentation hardness increases with friction, while the contact area decreases with friction. Therefore, the normal force during indentation is nearly independent of friction, as shown in Figure 10b. During the steady sliding stage, the scratching hardness is nearly independent of friction because the friction increases both the contact area and the normal force.
Although the scratching hardness Hs for Berkovich indenter is independent of f, Hs for conical indenter is linearly decreased with f, as demonstrated in Figure 11. The scratching hardness induced by ploughing, i.e., , remains nearly unchanged when f increases from 0 to 0.2. Figure 11. The hardness during scratching as a function of adhesion friction coefficient f for scratching with a 70.3° conical indenter. The value pairs of (Hs, f) is fitted to obtain the dash line.  The moving direction of the indenter relative to the sample determines the infl of friction on the contact area. The contact area during static indentation is sl decreased by friction, because the friction-induced downward shear stress applied sample surface results in a decrease in the contact area, as shown in Figure 1 contrast, friction leads to an increase in the contact area during the steady sliding because the friction-induced shear stress on the sample surface is upward indentation hardness increases with friction, while the contact area decreases friction. Therefore, the normal force during indentation is nearly independent of fr as shown in Figure 10b. During the steady sliding stage, the scratching hardness is independent of friction because the friction increases both the contact area and the n force.
Although the scratching hardness Hs for Berkovich indenter is independent o for conical indenter is linearly decreased with f, as demonstrated in Figure 11 scratching hardness induced by ploughing, i.e., , remains nearly unchanged w increases from 0 to 0.2. Figure 11. The hardness during scratching as a function of adhesion friction coefficien scratching with a 70.3° conical indenter. The value pairs of (Hs, f) is fitted to obtain the dash Figure 11. The hardness during scratching as a function of adhesion friction coefficient f for scratching with a 70.3 • conical indenter. The value pairs of (H s , f ) is fitted to obtain the dash line.

Plastic Deformation
It is reported that the stress and deformation are relieved by densification [13,25]. The stress distribution is significantly influenced by the geometry and location of the elasticplastic boundary [14,37]. In the analytical models of sliding stress fields, e.g., the Ahn model [16] and the Wang model [17], the plastic region is assumed to be a sphere with the center on the sample surface.
As shown in Figure 12, the maximum densification locates in the region that the indenter tip passes. The maximum value of densification predicted by FEA, i.e., 22.6%, is close to the measured saturated densification, i.e., 21% [38]. By comparing the densification contours in the top surface and the xz-cross-section shown in Figure 12, it is found that the contours in the yz-cross-sections are flat ellipses. Figure 12 also indicates that the densification in fused silica caused by scratching is bigger than that caused by indentation. It is reported that the stress and deformation are relieved by densification [13,25]. The stress distribution is significantly influenced by the geometry and location of the elastic-plastic boundary [14,37]. In the analytical models of sliding stress fields, e.g., the Ahn model [16] and the Wang model [17], the plastic region is assumed to be a sphere with the center on the sample surface.
As shown in Figure 12, the maximum densification locates in the region that the indenter tip passes. The maximum value of densification predicted by FEA, i.e., 22.6%, is close to the measured saturated densification, i.e., 21% [38]. By comparing the densification contours in the top surface and the xz-cross-section shown in Figure 12, it is found that the contours in the yz-cross-sections are flat ellipses. Figure 12 also indicates that the densification in fused silica caused by scratching is bigger than that caused by indentation.  Figure 13a shows that the elastic-plastic boundary in fused silica induced by scratching deviates from a circle. By contrast, Figure 13b demonstrates that the boundary can be tightly fitted by an ellipse. Figures 13 and 14 indicate that the yield region is an ellipsoid in the front of the indenter (x/hmax > 10), and a cylindroid at the rear of the indenter (x/hmax ≤ 10).  Figure 13a shows that the elastic-plastic boundary in fused silica induced by scratching deviates from a circle. By contrast, Figure 13b demonstrates that the boundary can be tightly fitted by an ellipse. Figures 13 and 14 indicate that the yield region is an ellipsoid in the front of the indenter (x/h max > 10), and a cylindroid at the rear of the indenter (x/h max ≤ 10).   As shown in Figure 13b, the shape of the elastic-plastic boundary is characterized by the length m of the semi-major axis, the length n of the semi-minor axis, and the depth ξ of the elastic-plastic boundary center. The semi-major and semi-minor axes are along the lateral and vertical directions, respectively. These parameters can be determined by fitting the elastic-plastic boundary by the following formula: The fitted results are shown in Figure 15. The real contact radius ar predicted by FEA equals 2.24 hmax. Figure 15a shows that m is bigger than a, while n is smaller than a. After unloading, m remains nearly unchanged, but n in the cross-section close to point B increases, due to the significant elastic recovery. m is significantly larger than n in the As shown in Figure 13b, the shape of the elastic-plastic boundary is characterized by the length m of the semi-major axis, the length n of the semi-minor axis, and the depth ξ of the elastic-plastic boundary center. The semi-major and semi-minor axes are along the lateral and vertical directions, respectively. These parameters can be determined by fitting the elastic-plastic boundary by the following formula: The fitted results are shown in Figure 15. The real contact radius a r predicted by FEA equals 2.24 h max . Figure 15a shows that m is bigger than a, while n is smaller than a. After unloading, m remains nearly unchanged, but n in the cross-section close to point B increases, due to the significant elastic recovery. m is significantly larger than n in the steady sliding stage. Therefore, the prediction accuracy of the existent sliding stress field models may be greatly improved if the spherical/cylindrical yield region is replaced by an ellipsoid/cylindroid. Figure 15b shows that the depth of the elastic-plastic boundary center in the yz-cross-section behind the indenter decreases rapidly when the indenter moves far away from it. ξ for scratching is much bigger than that for indentation (close to ξ at l B = 9.5 h max ) at both the fully loaded and the fully unloaded states. This indicates that the Ahn and Wang models should be refined to allow for the embedding of the center of the plastic zone. models may be greatly improved if the spherical/cylindrical yield region is replaced by an ellipsoid/cylindroid. Figure 15b shows that the depth of the elastic-plastic boundary center in the yz-cross-section behind the indenter decreases rapidly when the indenter moves far away from it. ξ for scratching is much bigger than that for indentation (close to ξ at lB = 9.5 hmax) at both the fully loaded and the fully unloaded states. This indicates that the Ahn and Wang models should be refined to allow for the embedding of the center of the plastic zone.

In the Sliding Stage
The scratching-induced stress contours under a conical indenter at the end of the sliding stage ② are shown in Figure 16. In the front of the indenter tip, i.e., in the region x ≥ 0, the shape of the contour lines induced by scratching is similar to those induced by indentation [29]. By contrast, at the back of the indenter tip, the contours are flattened, due to the plastic deformation left by the sliding indenter.

In the Sliding Stage
The scratching-induced stress contours under a conical indenter at the end of the sliding stage 2 are shown in Figure 16. In the front of the indenter tip, i.e., in the region x ≥ 0, the shape of the contour lines induced by scratching is similar to those induced by indentation [29]. By contrast, at the back of the indenter tip, the contours are flattened, due to the plastic deformation left by the sliding indenter.

In the Sliding Stage
The scratching-induced stress contours under a conical indenter at the end of the sliding stage ② are shown in Figure 16. In the front of the indenter tip, i.e., in the region x ≥ 0, the shape of the contour lines induced by scratching is similar to those induced by indentation [29]. By contrast, at the back of the indenter tip, the contours are flattened, due to the plastic deformation left by the sliding indenter. The maximal values of σx and σy are identified in the regions below the indenter tip just outside the elastic-plastic boundary, i.e., the regions R1 and R2 shown in Figure 16, respectively. They are the driving forces of median cracks. As the maximal σy is higher than the maximal σx, the median crack along the sliding direction tends to initiate prior to that along the lateral direction. The maximal value of σz locates at the far rear of the indenter (region R3), which is the driving force of lateral cracks. In the sliding stage, the driving force of median cracks is higher than that of lateral cracks.

At the Fully-Unloaded State
After unloading, as shown in Figure 17a, the stress contours predict a high σy on the sample surface at the front of the indenter (region R4), which is the driving force of radial cracks along the sliding direction, i.e., the radial crack 1 in Figure 18a. σy at the bottom of the yield region remains nearly unchanged after unloading. Therefore, median cracks along the sliding direction remain open in the unloading stage if they initiate in the sliding stage. The maximal σz increases from 0.087H to 0.137H during the unloading process. This indicates that the lateral crack emerges more easily in the unloading stage compared with the sliding stage. The maximal values of σ x and σ y are identified in the regions below the indenter tip just outside the elastic-plastic boundary, i.e., the regions R 1 and R 2 shown in Figure 16, respectively. They are the driving forces of median cracks. As the maximal σ y is higher than the maximal σ x , the median crack along the sliding direction tends to initiate prior to that along the lateral direction. The maximal value of σ z locates at the far rear of the indenter (region R 3 ), which is the driving force of lateral cracks. In the sliding stage, the driving force of median cracks is higher than that of lateral cracks.

At the Fully-Unloaded State
After unloading, as shown in Figure 17a, the stress contours predict a high σ y on the sample surface at the front of the indenter (region R 4 ), which is the driving force of radial cracks along the sliding direction, i.e., the radial crack 1 in Figure 18a. σ y at the bottom of the yield region remains nearly unchanged after unloading. Therefore, median cracks along the sliding direction remain open in the unloading stage if they initiate in the sliding stage. The maximal σ z increases from 0.087H to 0.137H during the unloading process. This indicates that the lateral crack emerges more easily in the unloading stage compared with the sliding stage.

Maximum Principal Stress
The contours of σ1/H shown in Figure 19 indicate that the median crack along the sliding direction tends to be the first crack to appear during the sliding stage. It initiates below the yield region at the rear of the indenter, i.e., in the region R2. Once initiated in the xz-cross-section, the median crack propagates along the sliding direction during scratching. As the value of σ1 in region R3 is smaller than that in region R2, the initiation load of median crack under indentation is higher than that under scratching. This is consistent with experimental observations. Though median crack is absent during indentation tests under the normal load of 40 N [39], it is observed during scratching tests under the normal load of 600 mN, as shown in Figure 18b. During the unloading stage, radial cracks may initiate from the sample surface at the front of the indenter, i.e., region R4, as shown in Figure 18a. As σy in region R4 is small before unloading, and significantly increases during the unloading process, radial cracks tend to emerge in the unloading stage. If both radial and median cracks form, they coalesce to form a big median-radial crack that penetrates through the entire yield region, as verified by the experiments (see Figure 18).

Maximum Principal Stress
The contours of σ1/H shown in Figure 19 indicate that the median crack along the sliding direction tends to be the first crack to appear during the sliding stage. It initiates below the yield region at the rear of the indenter, i.e., in the region R2. Once initiated in the xz-cross-section, the median crack propagates along the sliding direction during scratching. As the value of σ1 in region R3 is smaller than that in region R2, the initiation load of median crack under indentation is higher than that under scratching. This is consistent with experimental observations. Though median crack is absent during indentation tests under the normal load of 40 N [39], it is observed during scratching tests under the normal load of 600 mN, as shown in Figure 18b. During the unloading stage radial cracks may initiate from the sample surface at the front of the indenter, i.e., region R4, as shown in Figure 18a. As σy in region R4 is small before unloading, and significantly increases during the unloading process, radial cracks tend to emerge in the unloading stage. If both radial and median cracks form, they coalesce to form a big median-radial crack that penetrates through the entire yield region, as verified by the experiments (see Figure 18).

Maximum Principal Stress
The contours of σ 1 /H shown in Figure 19 indicate that the median crack along the sliding direction tends to be the first crack to appear during the sliding stage. It initiates below the yield region at the rear of the indenter, i.e., in the region R 2 . Once initiated in the xz-cross-section, the median crack propagates along the sliding direction during scratching. As the value of σ 1 in region R 3 is smaller than that in region R 2 , the initiation load of median crack under indentation is higher than that under scratching. This is consistent with experimental observations. Though median crack is absent during indentation tests under the normal load of 40 N [39], it is observed during scratching tests under the normal load of 600 mN, as shown in Figure 18b. During the unloading stage, radial cracks may initiate from the sample surface at the front of the indenter, i.e., region R 4 , as shown in Figure 18a. As σ y in region R 4 is small before unloading, and significantly increases during the unloading process, radial cracks tend to emerge in the unloading stage. If both radial and median cracks form, they coalesce to form a big median-radial crack that penetrates through the entire yield region, as verified by the experiments (see Figure 18).

Conclusions
This study developed three-dimensional finite element models of scratching on fused silica with 70.3° conical and Berkovich indenters. A refined elliptical model was used to consider the influence of densification on elastic properties and the saturation of densification with hydrostatic pressure. By analyzing the predicted scratching hardness, plastic deformation, and stress fields, the following conclusions are obtained: (1) The tangential hardness is slightly larger than the scratching hardness, and their ratio linearly increases with the adhesion friction coefficient f. The scratching hardness for an edge-leading Berkovich indenter is nearly independent of f, because the friction increases both the contact area and the normal force. By contrast, the scratching hardness for a conical indenter linearly decreases with f. These findings are helpful to model the friction-affected forces induced by scratching; (2) The densification effect should not be ignored if one aims for a damage-free process in fabrication. It is found that the material's densification under scratching is bigger than that under indentation. This indicates that the prediction accuracy of the sliding stress analyses can be improved if the material's densification is properly integrated into the modelling, e.g., replacing the spherical/cylindrical plastic zone with an ellipsoid/cylindroid, and considering the embedding of the plastic zone center; (3) Median cracks along the sliding direction tend to be the first cracks that emerge in the sliding stage. Radial cracks may initiate under a smaller load in the sample surface at the front of the indenter during the unloading stage. Radial and median cracks coalesce to form a big median-radial crack that penetrates through the entire yield region.

Conclusions
This study developed three-dimensional finite element models of scratching on fused silica with 70.3° conical and Berkovich indenters. A refined elliptical model was used to consider the influence of densification on elastic properties and the saturation of densification with hydrostatic pressure. By analyzing the predicted scratching hardness, plastic deformation, and stress fields, the following conclusions are obtained: (1) The tangential hardness is slightly larger than the scratching hardness, and their ratio linearly increases with the adhesion friction coefficient f. The scratching hardness for an edge-leading Berkovich indenter is nearly independent of f, because the friction increases both the contact area and the normal force. By contrast, the scratching hardness for a conical indenter linearly decreases with f. These findings are helpful to model the friction-affected forces induced by scratching; (2) The densification effect should not be ignored if one aims for a damage-free process in fabrication. It is found that the material's densification under scratching is bigger than that under indentation. This indicates that the prediction accuracy of the sliding stress analyses can be improved if the material's densification is properly integrated into the modelling, e.g., replacing the spherical/cylindrical plastic zone with an ellipsoid/cylindroid, and considering the embedding of the plastic zone center; (3) Median cracks along the sliding direction tend to be the first cracks that emerge in the sliding stage. Radial cracks may initiate under a smaller load in the sample surface at the front of the indenter during the unloading stage. Radial and median cracks coalesce to form a big median-radial crack that penetrates through the entire yield region.
Author Contributions: Conceptualization, C.L. and L.Z.; methodology, C.L., L.S., and X.W.; Figure 19. The contours of σ 1 /H at the (a) fully loaded and (b) fully unloaded states in the xz-crosssection. The thick dash line is the boundary between σ 1 ≈ σ y and σ 1 ≈ σ*, where σ* is the in-plane principal stress.

Conclusions
This study developed three-dimensional finite element models of scratching on fused silica with 70.3 • conical and Berkovich indenters. A refined elliptical model was used to consider the influence of densification on elastic properties and the saturation of densification with hydrostatic pressure. By analyzing the predicted scratching hardness, plastic deformation, and stress fields, the following conclusions are obtained: (1) The tangential hardness is slightly larger than the scratching hardness, and their ratio linearly increases with the adhesion friction coefficient f. The scratching hardness for an edge-leading Berkovich indenter is nearly independent of f, because the friction increases both the contact area and the normal force. By contrast, the scratching hardness for a conical indenter linearly decreases with f. These findings are helpful to model the friction-affected forces induced by scratching; (2) The densification effect should not be ignored if one aims for a damage-free process in fabrication. It is found that the material's densification under scratching is bigger than that under indentation. This indicates that the prediction accuracy of the sliding stress analyses can be improved if the material's densification is properly integrated into the modelling, e.g., replacing the spherical/cylindrical plastic zone with an ellipsoid/cylindroid, and considering the embedding of the plastic zone center; (3) Median cracks along the sliding direction tend to be the first cracks that emerge in the sliding stage. Radial cracks may initiate under a smaller load in the sample surface at the front of the indenter during the unloading stage. Radial and median cracks coalesce to form a big median-radial crack that penetrates through the entire yield region. Data Availability Statement: Some or all data, models, or codes generated or used during the study are available from the corresponding author by request.

Conflicts of Interest:
The authors declare no conflict of interest.