# Indentation stress fields in brittle materials: A micro-photoelastic investigation in silicate glasses

Gustavo Alberto Rosales-Sosa, Etienne Barthel, Yoshinari Kato, Matthieu Bourguignon, Akihiro Yamada, Tomiki Inoue, Shingo Nakane, Hiroki Yamazaki

► To cite this version:

Gustavo Alberto Rosales-Sosa, Etienne Barthel, Yoshinari Kato, Matthieu Bourguignon, Akihiro Yamada, et al.. Indentation stress fields in brittle materials: A micro-photoelastic investigation in silicate glasses. Acta Materialia, 2025, 292, pp.120973. 10.1016/j.actamat.2025.120973. hal-05038745

HAL Id: hal-05038745

https://hal.science/hal-05038745v1

Submitted on 17 Apr 2025

HAL is a multi-disciplinary open access
archive for the deposit and dissemination of sci-
entific research documents, whether they are pub-
lished or not. The documents may come from
teaching and research institutions in France or
abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire HAL, est
destinée au dépôt et à la diffusion de documents
scientifiques de niveau recherche, publiés ou non,
émanant des établissements d'enseignement et de
recherche français ou étrangers, des laboratoires
publics ou privés.

<figure><img src="image_1.png" /></figure>
Indentation stress fields in brittle materials: a
micro-photoelastic investigation in silicate glasses

Gustavo Alberto Rosales-Sosa¹, Etienne Barthel², Yoshinari Kato¹, Matthieu Bourguignon², Akihiro Yamada³,⁴, Tomiki Inoue¹, Shingo Nakane¹, and Hiroki Yamazaki¹

¹Fundamental Technology Division, Nippon Electric Glass Co., Ltd., 2-7-1, Otsu, Shiga, 520-8639, Japan.

²Soft Matter Sciences and Engineering, ESPCI Paris, PSL University, CNRS, Sorbonne University, Paris, 75005, France.

³Department of Materials Chemistry, School of Engineering, The University of Shiga Prefecture, Hikone, Shiga, 522-8533, Japan

⁴The Center for Glass Science and Technology, The University of Shiga Prefecture, Shiga, Hikone, 522-8533, Japan

April 17, 2025

Abstract

Indentation experiments have helped understand non linear mechanical properties of brittle ma-
terials such as plasticity, damage and fracture. However, our understanding of stress fields under
indentation remains limited due to the general lack of direct measurements. This study intro-
duces a novel approach to characterizing indentation stress fields in silicate glasses by combining
high sensitivity birefringence measurements, photoelastic calculations and finite element analy-
sis (FEA). We extensively investigated the elastoplastic response of soda-lime-silicate (SLS) and
silica glasses under different indentation conditions, highlighting the effects of composition and
indenter geometry on the photoelastic distributions. To predict the photoelastic response under
indentation, we carefully calibrated an elastoplastic constitutive relation for silicate glasses using
a combination of high-pressure (up to 25 GPa) and nanoindentation experiments. Computed
indentation stress fields can then be validated through the comparison of full 3D photoelastic
calculations and measured birefringence patterns. One key finding is that the residual stresses
arising from the calibrated constitutive relations offered a far more realistic representation of
the indentation stress fields than a commonly used approximate analytical elastoplastic model.
With this method, stress fields can be investigated not only in oxide glasses but also in any other
transparent isotropic material. While the calculated stress fields were generally satisfactory for
both glass compositions, they also evidenced that improvements in the constitutive relation are
needed for amorphous silica, which is known to undergo significant densification.
# 1 Introduction

Since the 1970s, indentation has emerged as an invaluable tool for the investigation of the advanced mechanical properties of brittle or quasi-brittle materials [1–4]. Indeed, during the application of a sharp tip on a surface, all major non-linear mechanisms such as plasticity, damage and cracking can develop in a relatively controlled manner, due to confinement by the elastic half-space. They can even be selectively activated through the choice of the applied load or the indenter geometry. As a result, indentation is used intensively to derive insight into the non-linear mechanical properties of many brittle materials such as single crystal oxides and nitrides, semi-conductors, ceramics and glasses, and the effect of composition and processing [5]. However, during indentation, deformation proceeds in a very inhomogeneous manner which makes it difficult to assess indentation stress and strain fields. This major drawback makes the quantitative evaluation of plasticity, damage or fracture toughness through indentation very dependent on the availability of adequate models.

In brittle materials, unlike plastic materials, indentation involves a significant elastic contribution [6]. This concept has in fact been at the center of the development of most analytical models. For example, Lawn *et al.* explained the formation of some of the indentation crack systems in terms of stiffness, hardness and fracture energy through a Boussinesq-type stress distribution, i.e. a model essentially based on elastic considerations [5, 7]. Several models have included the plastically deformed zone more explicitly, using Hill’s expanding cavity solution [8, 9], or equivalently, using a center of dilation with various levels of correction to account for the free surface [10, 11]. Among these models, Yoffe’s [12] has enjoyed great popularity because of its simple analytical form and has been used consistently to discuss crack patterns as a function of material type [2, 4, 12].

Despite these successes, several questions have arisen which have not yet been solved. For example, due to the sensitivity of the predicted elastic stress fields to material parameter variations, there has been some controversy on the suitability of indentation to measure fracture toughness [3, 13]. For indentation, silicate glasses stand out as model systems because, unlike single or poly-crystals, they are perfectly homogeneous and isotropic, with a simple elastic description and no initial defect distribution. Interestingly, Chiang *et al.* noticed that predicting the experimentally observed cracking threshold in silicate glasses was beyond the capabilities of their model [10]. More recently Davis *et al.* [14] found that Yoffe’s model deviates considerably from FEM predictions for indentation in silicate glasses, to a point which could raise doubts about the analysis of observed crack patterns. In fact, it is clear that the wide range of fracture thresholds observable in oxide glasses [15] cannot be explained within the standard indentation cracking theories [7] since they involve only elastic modulus, hardness and toughness and all three exhibit very moderate variations within the usual glass families [16].

For more quantitative analyses of indentation in brittle materials, direct measurements of stress fields would be very valuable. Two-dimensional photoelastic methodologies have been effectively used for several materials, including metal, ceramics, and polymers, to measure surface stress distributions. These methodologies involve the application of coatings [17], slicing [18], or through thin slabs with relatively homogeneous fields along the normal direction [19]. In the present case of indentation, it is necessary to deal with a full 3D distribution which cannot be reduced to a simpler 2D field by any of these methods. One of the very few tools available for that purpose, micro-birefringence, has been pioneered by Yoshida *et al.* [20] in silicate glasses. They measured birefringence patterns around sphere contacts on silica and soda-lime-silicate glasses and they backed out the stress field from the photoelastic response through an iterative numerical analysis based on the axisymmetry of the stress field ("onion-peeling" method) [21]. However, measuring the elastoplastic indentation stress field by the same method is more difficult: due to cracking at higher loads, the birefringence pattern must be measured on low load (subthreshold) indents, so that the birefringence signal is
weak because the indents are small. As a result of large curvatures and low signal to noise ratios, extraction of the stress field becomes unreliable.

In this work, we have accurately measured the low-birefringence signal around shallow indents in two silicate glasses using a highly-sensitive birefringence imaging system. In addition, we have developed a robust method to analyze the data, based on the direct calculation of the birefringence pattern from computed stress fields. This new method, however, requires accurately calibrated constitutive relations, which we also provide for by a combination of nanoindentation and high pressure experiments. Comparison of predictions and experiments validates the method which thus provides precise inputs for refinement of the constitutive relations. In contrast, the predicted photoelastic response based on Yoffe's approximate analytical stress field is shown to differ significantly from the data, except in the far field, confirming the analysis by Davis *et al.* [14]. More generally, these results demonstrate that birefringence can efficiently be combined with micromechanics experiments for better quantification.
## 2 Experimental

The experiments were performed on two typical glass compositions: amorphous silica (SiO₂ - GE-124 - type I) and soda-lime-silicate (SLS - commercial glass) glass. Although these two glasses share similar Young's modulus they differ significantly in their flow stress, which is directly related to their hardness [16]. In addition, they also differ significantly in their ability to undergo permanent volumetric strain, which is expected to affect the residual stresses. The basic elasto-optical properties of both glasses are shown in table 1. The photoelastic constant *C* and refractive index *n*<sub>d</sub> of the glasses were taken as given by the providers while the elastic properties were measured.

### 2.1 Measurement of birefringence patterns from indentation experiments

The birefringence patterns were measured for two different types of loading. The first was a purely elastic sphere-plane contact measured in situ under load. The second was an elastoplastic indent for which the residual stress field was measured ex situ, after unloading. In this case, two types of indenters were used : cone and sphere with small diameter.

#### 2.1.1 In-situ elastic contact

The patterns were measured using a birefringence imaging system (WPA-micro, Photonics Lattice, Japan) with a wavelength of $\lambda = 543$ nm. The integrated photoelastic response can be measured with a maximum retardation range of 3500 nm with a repeatability better than 1 nm. Cubic blocks (10x10x10 mm) of SiO₂ and SLS glass were indented with a stainless steel (SUJ-2) sphere (radius $R = 2.75$ mm) as shown in figure 1 A. The images were taken using an objective lens with magnifications of 5X and 10X with a field of view of 1.04 mm x 0.80 mm and 0.52 mm x 0.40 mm, respectively. The resolution of the retardation and azimuth images was 384 x 288 pixels. The spatial x-y resolutions is estimated to be about half of the system wavelength (~ 260 nm). Photoelastic patterns were obtained at different applied loads, measured using a 20 N load cell (accuracy 0.01 N). The applied loads were 4.25 N, 6.30 N, 6.40 N, 8.58 N, and 10.03 N for silica glass and 6.57 N, 8.18 N, 8.29 N, 10.35 N, and 13.30 N for SLS glass.

#### 2.1.2 Residual stresses in elastoplastic contact

The photoelastic measurements were performed using a high sensitivity birefringence imaging system (PA-micro, Photonic Lattice, Japan) operating at a wavelength of $\lambda = 520$ nm with a maximum retardation range of 130 nm. Images were taken using an objective lens with a magnification of 100X and a field of view of 0.08 mm x 0.07 mm. The space resolution is estimated to be half the wavelength of the light source. Here, to limit spurious birefringence signals, the samples were redrawn silica [22] and redrawn SLS (commercial glass) glass fibers with 500 µm square sections (figure 1 B). The glass fibers were indented using a cone (diamond - semi-angle 70° and tip radius about 1 µm) at a load of 0.98 N (100 gf) and with a sphere (diamond - 47 µm radius) at a load of 5.88 N (600 gf). The indented fibers were immersed in an oil with a similar refractive index $n_d$ to the glass to avoid scattering of light at the interfaces. The photoelastic patterns were taken at a resolution of 616 x 514 pixels.

### 2.2 Prediction of birefringence patterns from constitutive relations

To calculate the birefringence patterns, our approach consists of the following steps: 1) calibrate existing constitutive relations for the two glasses by a combination of high-pressure (up to 25 GPa)
<figure><img src="image_2.png" /><figcaption>Figure 1: (A) Schematic diagram of the elastic contact using a glass cube (in-situ) and (B) elasto-plastic contact (residual) birefringence measurement using a squared-section glass fiber</figcaption></figure>

<table><caption>Table 1: Elastic and optical properties of used glasses: Young's modulus *E*, Poisson's ratio *ν*, stress optical coefficient *C*, and refractive index *n*<sub>d</sub>.</caption><thead><tr><th>Property</th><th>Silica glass</th><th>SLS glass</th><th>Stainless steel</th></tr></thead><tbody><tr><td><i>ene E</i> (GPa)</td><td>72</td><td>73</td><td>202</td></tr><tr><td><i>ene ν</i></td><td>0.17</td><td>0.21</td><td>0.3</td></tr><tr><td><i>C</i> (TPa<sup>-1</sup>)</td><td>3.50</td><td>2.70</td><td>-</td></tr><tr><td><i>n</i><sub>d</sub></td><td>1.46</td><td>1.52</td><td>-</td></tr><tr><td><i>ene</i></td><td></td><td></td><td></td></tr></tbody></table>

and nanoindentation experiments, 2) compute the indentation stress distributions through FEA, 3) compute the optical photoelastic patterns (retardation and azimuth angles) from the FEA stress distributions using the light-propagation equations proposed by Aben *et al.* [23, 24]

### 2.2.1 Generic constitutive relation for silicate glasses

The constitutive relation for the elastoplastic behavior of silicate glasses follows Kermouche *et al.* [25] with the addition of saturation for the densification, as detailed below. In this constitutive model, the yield surface depends upon pressure (figure 2). There is pressure hardening with volumetric strain while the flow stress remains constant with plastic deformation. The elliptic yield surface $f(\sigma_{ij})$ is described by

$$f(\sigma_{ij}) = \left(\frac{q}{q_c}\right)^2 + \left(\frac{p}{p_c(\epsilon_v^p)}\right)^2 - 1, \quad \text{if } p > 0 \qquad (1)$$

$$f(\sigma_{ij}) = q - q_c, \quad \text{if } p \le 0 \qquad (2)$$
<figure><img src="image_3.png" /><figcaption>Figure 2: Elliptical yield surface for silicate glasses with pressure-hardening due to densification ($\Delta\rho/\rho_0$). Under tension (negative pressures), the material follows a perfectly elastic-plastic von Mises model.</figcaption></figure>

where $\sigma_{ij}$ is the Cauchy stress tensor, $q$ is the von Mises stress, $q_c$ is the flow stress, $p$ is the hydrostatic pressure and $p_c$ is the onset pressure for densification under pure hydrostatic pressure. For the evolution of $p_c$ with volumetric plastic strain, Kermouche et al. [25] included only a simple hardening slope matching diamond anvil cell experiments and this form can also been used in Abaqus, either shoehorned into the cap of a Drücker-Prager model [26] or implemented as a UMAT by G. Molnar [27]. Here we used a modified version of Molnar’s UMAT and introduced the saturation of the densification $\xi$ at high pressures as a bi-linear process

$$p_c(\epsilon_v^p) = p_{c0} + \gamma\epsilon_v^p, \quad \text{if } \xi < 0.95 \xi_{sat} \qquad (3)$$

$$p_c(\epsilon_v^p) = p'_c + \gamma'(\epsilon_{v,sat}^p - \epsilon_v^p), \quad \text{if } 0.95 \xi_{sat} \le \xi \le \xi_{sat} \qquad (4)$$

$$\gamma' = \frac{p_{sat} - p'_c}{\epsilon_{v,sat}^p - \epsilon_v^{p'}} \qquad (5)$$

where $p_c$ is the current onset pressure for densification, $p_{c0}$ is the initial onset pressure, $\gamma$ is the hardening slope of the first regime, $p'_c$ is onset pressure for densification at $\xi = 0.95 \xi_{sat}$, $\gamma'$ is the hardening slope of the second regime, $\epsilon_{v,sat}^p$ is the volumetric plastic strain at saturation, and $\epsilon_v^{p'}$ is the volumetric plastic strain for $\xi = 0.95 \xi_{sat}$. As an approximation, we assumed that full saturation $\xi_{sat}$ is reached at a pressure $p_{sat}$ of 100 GPa. The transition between the first and second linear regimes was conveniently chosen to be 0.95 $\xi_{sat}$. The relationship between densification $\xi$ and volumetric plastic strain $\epsilon$ was assumed for finite strain formulation [28],

$$\epsilon_v^p = -\ln(1 + \xi) \qquad (6)$$

Associated plasticity was used to determine the plastic flow $\dot{\epsilon}^p$. The stresses and strains were obtained using a full-implicit scheme from the yield surface and the consistency condition [25, 29].

### 2.2.2 Calibration of the constitutive relation parameters

The calibration consists of four parameters: yield pressure $p_{c0}$, pressure-hardening slope $\gamma$, saturation volumetric strain $\epsilon_{v,sat}^p$ and flow stress $q_c$. The pressure-hardening behavior of silica glass
was taken from the literature [30, 31]. For SLS glass, the data available in the literature is much more limited, however. For this reason, high-pressure Diamond Anvil Cell (DAC) and Multi Anvil Cell (MAC) experiments were performed and combined with earlier results [15] to obtain a complete description of the pressure-hardening behavior. Finally, we performed nanoindentation experiments in silica and SLS glasses to determine $q_c$ [25].

### Evaluation of $p_{c0}$ through Diamond Anvil Cell experiments

The DAC experiments on SLS glass were performed at the University of Shiga prefecture (Japan). Glass fibers with a diameter of 100 µm were obtained by re-melting of the SLS glass sample. The fiber was then cut and sanded down to a 100 µm thickness. The glass sample was compressed statically within the symmetry-type DAC equipped with 800 µm culet single crystalline diamond anvils. Stainless steel (SUS304, 300 µm thickness) with a 250 µm hole (sample room), was used for the gasket. The mixture of methanol-ethanol (4:1 in volume) was put in the sample room as the pressure-transmitting medium together with a ruby sphere, which is used as the pressure indicator by the ruby fluorescence method [32].

The two-dimensional imaging technique using visible light was exploited to measure the linear strain of the sample under hydrostatic stress. The sample image was taken with a CMOS camera (AdvanCam-U3II, AdvanVision Co., Ltd.) equipped with an optical microscope (BX53M, Olympus). To make the boundary of the sample sharper, a green filter (43IF550-W45, Olympus) was put in the white light source by LED. An example taken at a pressure of 1.95 GPa is shown in Figure 3A.

Since glass is an isotropic material, it shrinks homogeneously under perfectly hydrostatic pressure. The volumetric strain ($V/V_0$), where $V$ and $V_0$ are the volume at high pressure and ambient condition, respectively, can be estimated from the linear strain determined from analysis of the sample image using,

$$ \left(\frac{l}{l_0}\right)^3 = \frac{V}{V_0} \qquad (7) $$

where $l$ and $l_0$ are the sample diameters at high pressure and ambient conditions, respectively. Finally, the isothermal (in situ) density of the glass $\rho_{st}$ can be obtained from the initial density $\rho_0$ by $\rho_{st}/\rho_0 = V_0/V$. The diameter of the sample was determined from the line profile shown in Figure 3A. The final sample diameter was the average of eight line profiles, and the error was the standard deviation. Since the sample shape is not perfectly circular, the difference in the length measured with the eight-line profile was checked at ambient pressure before compression. The scatter can be regarded as the difference from the perfect circle, which is subtracted from the error in the diameter at high pressure, $l$.

The measured density as a function of pressure for 4 independent runs is summarized in Figure 3B. The experimental data show a visible slope change at a pressure of about 4 GPa which is taken as the onset of densification. This onset pressure $p_{c0}$ seems smaller than that (7 GPa - 9 GPa) reported from Multi Anvil Cell experiments [30, 33].

### Evaluation of pressure-hardening slope through Multi Anvil Cell experiments

Multi Anvil Cell (MAC) experiments were performed on SLS samples compressed at pressures of 15 GPa, 20 GPa, and 25 GPa at room temperature at the Geodynamics Research Center at Ehime University in Japan. At first, SLS samples were processed to obtain cylinders with a diameter and height of about 3-4 mm. Samples were compressed using the 6-8-type double-stage compression method with a 3000-ton hydraulic press apparatus (Orange-3000, Sumitomo Heavy Industries, Ltd.,
<figure><img src="image_4.png" /><figcaption>Figure 3: (A) Example of the in-situ line-profile analysis for SLS sample diameter (l) determination. The disk sample image was taken at a pressure of 1.95 GPa. (B) Experimental in-situ DAC compression densities estimated from sample dimensions vs. applied pressure for SLS glasses for several runs. The experiments reveal that the onset pressure for densification is about $p_{c0} = 4.0$ GPa.</figcaption></figure>

Japan). The 6-8-anvil apparatus is a split sphere housing six hardened steel outer anvils which in turn enclose eight tungsten carbide inner anvils shaped as truncated cubes to form an octahedral sample cavity (4 mm edge length) filled with semi-sintered (Mg, Co)O pressure-transmitting medium [34]. Samples were compressed to the target pressure in less than 7 hours with a dwell time of 1 h. After that, the samples were released to ambient pressure gradually during a time ranging from 4 hours to 8 hours. The density of the samples compressed at 10 GPa was measured using the Archimedes method in toluene. The density of the samples compressed at 20 GPa and 25 GPa was measured using the sink-float method using properly mixed 1,1,2,2-tetrabromoethane and 1,1,2,2-tetrachloroethane.

From the residual density, it is possible to quantify the relative density change (densification $\xi$) as

$$ \xi = \frac{(\rho - \rho_o)}{\rho_o} \qquad (8) $$

The densities and estimated densification values of pristine and compressed glasses are shown in Table 2. As observed, the maximum densification is about 4% for an applied pressure of 25 GPa.

The experimental MAC results in this work were combined with previously reported results for SLS glass by Kato et al. [15] to obtain a complete pressure-hardening curve. This hardening curve (densification vs. applied pressure) for SLS glass is shown in figure 4 (right) along with the hardening curve for amorphous silica (left) taken from the literature [30, 31] and bilinear hardening approximations. The onset pressure obtained by DAC for SLS glass is also shown. The pressure
<table><caption>Table 2: Residual measured density of SLS glass samples as a function of applied pressure in MAC experiments.</caption><thead><tr><th>Pressure p (GPa)</th><th>Density ρ (g/cm³)</th><th>Densification ξ (%)</th></tr></thead><tbody><tr><td>0</td><td>2.491</td><td>0.0</td></tr><tr><td>15</td><td>2.521</td><td>1.2</td></tr><tr><td>20</td><td>2.587</td><td>3.9</td></tr><tr><td>25</td><td>2.585</td><td>3.8</td></tr><tr><td>ine</td><td></td><td></td></tr></tbody></table>

hardening parameters are shown in table 3 for both glasses. It is noted that the onset pressure and maximum densification of silica are significantly larger than those of SLS glass. In addition, the densification proceeds faster with pressure in silica compared to the case of SLS, as revealed by the hardening slope $\gamma$.

<figure><img src="image_5.png" /><figcaption>Figure 4: (left) Experimental MAC densification data for silica glass taken from Rouxel *et al.* [30] and Deschamps *et al.* [31] along with the proposed bilinear approximation. (right) Experimental MAC densification data for SLS glass from this work and the data points from Kato *et al.* [15] along with the bilinear approximation. The onset pressure $p_{c0}$ extracted from DAC for SLS glass is shown as a green triangle.</figcaption></figure>

### Evaluation of flow stress $q_c$ from nanoindentation experiments

Nanoindentation experiments were performed on silica and SLS glasses at the National Institute of Materials Science (NIMS) in Japan with a diamond Berkovich tip (Bruker-Hysitron TI950 Triboindenter, USA). The indentations were performed in air up to a maximum load of 500 mN (load-control) with a constant loading rate of 70 mN/s and no dwell time. Five indentations were performed at the same load to check the reproducibility of the test. The experimental and calculated load-displacement curves for silica and SLS glasses are shown in figure S4 of the supplementary materials.
<table><caption>Table 3: Pressure-hardening parameters $p_{c0}$, $\gamma$, and $\xi_{sat}$, and used for the constitutive model of silica and SLS glasses. The values of densification $\xi$ are expressed as fractions.</caption><thead><tr><th>Parameter</th><th>Silica glass</th><th>SLS glass</th></tr></thead><tbody><tr><td>ine $p_{c0}$ (GPa)</td><td>10.0</td><td>4.0</td></tr><tr><td>$\gamma$ (GPa)</td><td>40</td><td>330</td></tr><tr><td>$\xi_{sat}$</td><td>0.21</td><td>0.04</td></tr><tr><td>ine</td><td></td><td></td></tr></tbody></table>

The recorded indentation hardness $H_{IT}$ for silica and SLS glasses were 8.93 ± 0.10 GPa and 6.66 ± 0.05 GPa, respectively.

A solid axisymmetric FEA model was built (ABAQUS 2022) to simulate the nanoindentation experiments for both glass compositions. For this purpose, a perfectly sharp solid diamond cone indenter with a semi-angle of 70° was modeled along with a glass substrate. This axisymmetric geometry provides a response equivalent to the three-dimensional Berkovick indenter. The FEA solid model is shown in figure S5 of the supplementary materials. Large strain formulation was used and the contact was modeled using the penalty method under frictionless conditions. Four-node fully integrated elements (CAX4) were used for the indenter and the glass sample. A high mesh density was used near the contact area with an average element size of 0.25 µm.

Values of $q_c$ of 7.3 GPa and 4.8 GPa for silica and SLS glasses gave best fits to the experimental results. These values are slightly larger than reported by Kermouche *et al.* [25].

### 2.2.3 Calculation of the birefringence patterns from the FEA stress distributions

#### Indentation stresses from FEA calculations

All the stress fields were calculated by FEA under the assumption of axisymmetry. The model geometries and material properties for the elastic and residual (elastoplastic) calculations were taken from section 2.1. For the purely elastic case, the FEA solid model was about 1.5 mm x 2.4 mm in size with an average mesh size around the contact zone of 2.5 µm. The elastic properties were taken from table 1. The results were compared to the analytical solutions of the linear elastic frictionless (Hertz) contact problem. For the residual case (cone and small sphere), the properties and geometries of the models were the same as those shown in section 2.2.2.

#### Photoelastic calculations

The photoelastic patterns can be calculated by integration of Aben's equations [23, 24] for the propagation of light in a condensed medium. A similar method has been used by Yokoyama *et al.* for measurements of in-situ Hertzian stress fields in soft elastic materials [35]. This methodology is based on a linear stress-optical law; therefore, it cannot describe non-linearities and especially the optical effects of plastic deformation. For that purpose, it would be necessary to calibrate the effects of plasticity on the photoelastic response of the material, which is beyond the scope of the present paper. The conclusions drawn from this work are only based on the integrated optical response of elastically deformed media. The region of the image affected by plasticity, although it contains its most salient features, has been excluded from the comparison to the predictions. The region of the image affected by plasticity, although it contains some salient features, has been shown but excluded from the quantitative comparison between data and predictions.
Assuming an orthogonal $x_1, x_2, x_3$ coordinate system with light propagating in the $x_2$ direction, the components of the electric field perpendicular to the propagation direction obey $\bar{E} = (E_1, E_3)$ and $d\bar{E}/dx_2 = -i\pi\bar{P}\bar{E}(x_1, x_2, x_3)$ where $\bar{P}$ is the propagator

$$ \bar{P}(x_1, x_2, x_3) = \begin{pmatrix} \mathcal{P}_{11} & \mathcal{P}_{13} \\ \cdot & -\mathcal{P}_{11} \end{pmatrix} \qquad (9) $$

where

$$ \mathcal{P}_{11} = \frac{C(\sigma_{11} - \sigma_{33})}{2} \qquad (10) $$

$$ \mathcal{P}_{13} = C\sigma_{13} \qquad (11) $$

and $C$ is the stress-optical coefficient which is usually expressed in Pa$^{-1}$ or Brewsters ($C_{Pa^{-1}} = C_{Brewster} \times 10^{-12}$). The stress tensor at any position in the solid is easily obtained by rotation of the axisymmetric solution of the FEA. Denoting $\bar{E}_n = E_{n,r} + iE_{n,i}$ for $n = 1,3$ the electric field in the $x_1$ and $x_3$ directions, Aben's equations read,

$$ \frac{d}{dx_2} \begin{pmatrix} E_{1,r} \\ E_{1,i} \\ E_{3,r} \\ E_{3,i} \end{pmatrix} = \pi \begin{pmatrix} 0 & \mathcal{P}_{11} & 0 & \mathcal{P}_{13} \\ -\mathcal{P}_{11} & 0 & -\mathcal{P}_{13} & 0 \\ 0 & \mathcal{P}_{13} & 0 & -\mathcal{P}_{11} \\ -\mathcal{P}_{13} & 0 & \mathcal{P}_{11} & 0 \end{pmatrix} \begin{pmatrix} E_{1,r} \\ E_{1,i} \\ E_{3,r} \\ E_{3,i} \end{pmatrix} \qquad (12) $$

For each image point $(x_1, x_3)$, they are integrated numerically using a Runge Kutta scheme along the $x_2$ axis from $-\infty$ to $+\infty$ i.e. from far away behind the stressed area to far away ahead. As a result, the FEA model size has to be large enough for the solution to fulfill this condition. Finally the retardation $\Delta$ and azimuth $\phi$ are calculated through

$$ \Delta = 2 \arccos(E_{1,r}) \qquad (13) $$

$$ \phi = \arctan(-E_{3,i}/E_{1,i})/2 \qquad (14) $$

The photoelastic calculation and a Python implementation are detailed in the supplementary materials (S1 and S2). An example calculation is also given in the supplementary materials (S3). Note that the assumption of axisymmetry is not a strong requirement for the present method to work. When the stress tensor is known everywhere in space, integration of the Aben equations is as straightforward as in the axisymmetric case and will require only minor modifications of the Python code.
# 3 Results

We first present the results and predictions for the *in situ* elastic contact, which were obtained for validation purposes. Then we present the measured birefringence patterns for residual stresses in elastoplastic contacts and compare them first to the predictions based on Yoffe's analytical model, then to the predictions based on the FEA calculations.

## 3.1 In-situ elastic contact

The experimental and computed photoelastic patterns (retardation and azimuth) for silica glass indented at a load of 10.03 N with the large sphere indenter (purely elastic contact) are shown in figure 5. The distances have been normalized by the contact radius $a= 72.76$ µm. The azimuth angle at 0° corresponds to the direction of the slow optical axis. We find that the retardation forms a lobe with maximum intensity somewhat below the surface. The azimuth map is characterized by a roughly conical region of maximum extending from the edge of the contact at about 30° from the surface. The small white lines describe the orientation of the azimuth angle for the slow optical axis. The agreement between the experimental and computed patterns is remarkable regarding both distribution and intensity.

<figure><img src="image_6.png" /><figcaption>Figure 5: Experimental and computed *in situ* integrated photoelastic patterns for silica glass indented at a load of 10.03 N: Retardation $\Delta$ (top) and Azimuth $\phi$ (bottom). The images were obtained using a 5X objective lens. The plot field of view was normalized to each contact radius $a$ as obtained from FEA calculations. The small white lines represent the orientation of the azimuth angle; a horizontal line corresponds to zero azimuth angle.</figcaption></figure>
For a more detailed comparison between the experimental and the computed birefringence pat-
terns, profiles of the retardation Δ along the loading axis were plotted for both silica and SLS
glasses at several loads (figure 6). Without any adjustable parameters, the calculations almost per-
fectly reproduce the experimental results, revealing that the calculations scale correctly with applied
load. The maximum retardation is located at a position of about 0.5 times the contact radius as
expected from the Hertzian equivalent shear stress fields for sphere indentation [36]. For SLS glass,
the retardation intensity is always lower than silica glass at a similar applied load. This is expected
since the stress-optical coefficient for SLS glass is lower than for silica glass while they have nearly
identical elastic modulus. We also found experimentally that the maximum retardation scales with
applied load as $\Delta_{mxz} \propto P^{2/3}$, a law which can be established by assuming the retardation is ap-
proximately given by the product of typical hertzian contact stress, contact size and photoelastic
constant. Details can be found in the supplementary materials (figure S6).

<figure><img src="image_7.png" /><figcaption>Figure 6: Experimental *in-situ* retardation depth profiles along the normalized loading axis *z/a* for both silica (A) and SLS (B) glasses at several applied loads. Solid lines represent experiments, and dashed lines represent FEA calculations. The experimental data for some loads were taken using different objective lenses (5X and 10X). The inset figure shows the direction of the selected profile. The contact radii *a* were obtained from FEA calculations.</figcaption></figure>
## 3.2 Sphere and cone indentation - residual stresses

We now present the measured residual photoelastic patterns after indentation by cone and sphere indenters for both silica and SLS glass and compare them with the predictions based on first Yoffe's model then the FEA stress fields.

### 3.2.1 Experimental residual photoelastic patterns

The residual photoelastic patterns of indented silica and SLS glass are shown in figure 7 for cone indentation and in figure 8 for sphere indentation. Following the scaling principle mentioned above, in these plots, the retardation is normalized by the product of the photoelastic constants *C* and the contact radius *a* while distances are normalized by *a*. The normalized retardation is therefore a direct indication of the stress level. The white lines shown in figure 7 are examples of the orientation of the azimuth angle in the residual state.

Regarding the retardation distribution, for both glass compositions and both indenter types, the normalized retardation distribution is dominated by the well-known "bowl-shaped" region below the surface which roughly corresponds to the plastically deformed area. The overall stress level of SLS glass is more significant than that of silica glass by approximately a factor of two; a similar contrast in birefringence has been noticed by Arora et al. [37]. Also, we note that the stress levels are more significant for the cone indenter than for the sphere, although the latter was performed at a six times larger load. It is well documented that sharper indenters lead to more considerable residual stresses and easier cracking due to enhanced plastic flow [38, 39]. The residual stresses naturally extend outside the plastic zone on the sides and below where they exhibit a distinctive distribution. Two low-retardation spots or "black-holes" can be observed. For silica they are near the surface and just below the plastic zone ($z/a \simeq 1.0$). For SLS, they are located deeper, at the bottom of the plastic zone and at some distance below ($z/a \simeq 1.6$). Inside the plastic zone, the birefringence patterns are different for the two glasses. Silica shows a low-birefringence area at the center of the plastic zone, while for SLS glass, high-birefringence regions appear on the sides and top of the plastic zone.

The azimuth distributions are also globally similar for silica and SLS glasses and for the two types of indenter. They feature a distinctive "inverted $\Psi$" concave distribution which is associated with a vertical segment of 0°-180° azimuth orientation along the loading axis, similar to the elastic in situ cases. Below this segment a plateau joins two symmetrically inclined lines which are 0°-180° cuts. The "black-holes" mark the endpoints of the vertical segment. The four distributions are differentiated by some details, though: for silica glass, the "inverted $\Psi$" azimuth distribution is narrower than that of SLS glass and even displays a convex shape for sphere indentation.
<figure><img src="image_8.png" /><figcaption>Figure 7: Residual experimental photoelastic patterns: normalized retardation (top) and azimuth (bottom) after conical indentation ($\alpha = 70^{\circ}$) on the squared silica and SLS glass fibers at a load of 0.98 N. The image dimensions were normalized by the imprint radius $a$. The experimental residual imprint radii are $a_{SiO_2} = 5.4$ µm and $a_{SLS} = 6.4$ µm. Retardation was normalized by the product of stress-optical coefficient $C$ and the contact radius $a$. The small white lines represent the orientation of the azimuth angle; a horizontal line corresponds to zero azimuth angle.</figcaption></figure>
<figure><img src="image_9.png" /><figcaption>Figure 8: Residual experimental photoelastic patterns: normalized retardation (top) and azimuth (bottom) after spherical indentation ($R = 47\mu m$) at a load of 5.88 N on the squared silica and SLS glass fibers. The image dimensions were normalized by the imprint radius $a$. The experimental residual imprint radii are $a_{SiO_2} = 13.7\mu m$ and $a_{SLS} = 15.2\mu m$. Retardation was normalized by the product of stress-optical coefficient $C$ and the contact radius $a$.</figcaption></figure>
### 3.2.2 Photoelastic patterns computed from Yoffe's elastoplastic stress fields

Following Davis' observations [14] we first considered the residual stress distribution proposed by Yoffe [12]. This model describes the stress distribution around the indentation plastic zone as the sum of a Boussineq-type field, i. e. an elastic point load, for the external force and a term arising from the presence of the plastic core and the effect of the surface, hence the form and denomination of blister field. Yoffe's model in polar coordinates $(r, \theta, \phi)$ is expressed as

$$ \sigma_{rr}(r, \theta) = \frac{P(1 - 2\nu)}{2\pi r^2} \left( 1 - \frac{2(2 - \nu)}{(1 - 2\nu)} \cos \theta \right) + \frac{4B}{r^3} \left( (5 - \nu) \cos^2 \theta - (2 - \nu) \right) \quad (15) $$

$$ \sigma_{\theta\theta}(r, \theta) = \frac{P(1 - 2\nu)}{2\pi r^2} \left( \frac{\cos^2 \theta}{1 + \cos \theta} \right) - \frac{2B}{r^3} (1 - 2\nu) \cos^2 \theta \quad (16) $$

$$ \sigma_{\phi\phi}(r, \theta) = \frac{P(1 - 2\nu)}{2\pi r^2} \left( \cos \theta - \frac{1}{1 + \cos \theta} \right) + \frac{2B}{r^3} (1 - 2\nu)(2 - 3 \cos^2 \theta) \quad (17) $$

$$ \tau_{r\theta}(r, \theta) = \frac{P(1 - 2\nu)}{2\pi r^2} \left( \frac{\sin 2\theta}{2(1 + \cos \theta)} \right) + \frac{2B}{r^3} (1 + \nu) \sin 2\theta \quad (18) $$

where $\nu$ is the Poisson's ratio, $B$ the blister-field coefficient, and $P$ the load. The indices $r$, $\theta$ and $\phi$ correspond to the radial, tangential and azimuthal coordinates in the polar coordinate system, respectively. We can transform readily the stress in polar coordinates to a more convenient cylindrical coordinate $(r, z, \theta)$ system using the rotation matrix,

$$ \sigma_r = \sigma_{rr} \sin^2 \theta + \sigma_{\theta\theta} \cos^2 \theta + \tau_{r\theta} \sin 2\theta \quad (19) $$

$$ \sigma_z = \sigma_{rr} \cos^2 \theta + \sigma_{\theta\theta} \sin^2 \theta - \tau_{r\theta} \sin 2\theta \quad (20) $$

$$ \sigma_{\theta} = \sigma_{\phi\phi} \quad (21) $$

$$ \tau_{rz} = (\sigma_{rr} - \sigma_{\theta\theta}) \sin \theta \cos \theta + \tau_{r\theta} \cos 2\theta \quad (22) $$

The transformation allows for direct comparison with FEA calculations and for computing the photoelastic patterns using the proposed method. The blister field coefficient can be calculated through the following relation,

$$ B = cHa^3 \quad (23) $$

$H$ is the hardness, $a$ is the contact radius, and $c$ is a tunable dimensionless parameter defined as the blister field strength coefficient. Davis et al. [14] fitted the parameter $c$ in the case of silica glass under conical indentation to reproduce the maximum residual hoop stress observed by FEA. In this work, we used the same procedure, which led to a value of $c = 0.087$ for the case of SLS glass under conical indentation at a load of 0.98 N. The hardness used for the calculation was obtained experimentally from the nanoindentation experiments at a load of 500 mN. The residual photoelastic patterns (retardation and azimuth) computed with Yoffe's blister field are shown in figure 9

We observe that near the plastic zone, where cracks are expected to develop, usually during the unloading part of the indentation cycle [2], most of the features we have observed experimentally such as the "black holes" and "inverted $\Psi$," are not reproduced by Yoffe's model. The retardation has an almost perfect angular symmetry around the plastic zone, which differs markedly from the
experimental results. However, the patterns become more similar farther away from the plastic zone. Better agreement in the far field could be expected given that the blister term is a point expansion and the contact size is not included in the distribution apart from the prefactor. Also, when normalizing the obtained retardation intensity by the product of the photoelastic constant and the contact radius, we obtain stress values in the edge of the plastic zone about six times larger than the experimental results. This discrepancy could be corrected by a reduction of *c* but the overall distribution would remain unchanged because in the unloaded state, there remains only the blister term.

<figure><img src="image_10.png" /><figcaption>Figure 9: Residual retardation and azimuth distributions using Yoffe's blister field for the case of conical indentation in SLS glass at a load of 0.98 N. A value of *c* = 0.087 was set to reproduce the same maximum residual hoop stress at the surface as obtained from FEA calculations. Yoffe's model results are only valid for distances larger than *a*. The small white lines represent the orientation of the azimuth angle; a horizontal line corresponds to zero azimuth angle.</figcaption></figure>

### 3.2.3 Photoelastic patterns from the stress distributions calculated by FEA

We have also computed the photoelastic patterns using the stress fields from the FEA calculations with the calibrated constitutive relations of the two glasses. The results for conical indentation of SLS glass are shown in figure 10 and for silica glass in figure 11. The results for spherical indentation in SLS glass and silica glass are shown in the supplementary materials (figures S7 and S8).

#### Soda-lime-silicate glass

Figure 10 summarizes the comparison between birefringence experiments and calculation for the case of SLS glass indented with a cone. Radial plots (figure 10A and 10B) were used to evaluate the discrepancies between experiments and calculations in different directions starting from the contact
point. The profiles show satisfactory agreement between experiments and calculations. Large deviations in retardation are observed within the plastic zone and near the surface. It is noted that the agreement with the azimuth angles is remarkable, even inside the plastic zone. We do not expect the present approach to work in the plastic zone because the optical response of the plastically deformed material is unknown so that we note the qualitative agreement in this region but will not comment it further. Outside of the plastic zone, the face-to-face maps (figure 10D) show that most of the retardation and azimuth features are reproduced by the calculations. The agreement with the azimuth profile is better than the retardation profiles (stress orientations are reproduced better than their intensities). In addition, the calculation results on the assumption of different saturation densification $\xi_{max}$ (0.01 and 0.08) are also shown in figure 10A. It is observed that the maximum densification ability of the glass $\xi_{max}$ has little to no impact on the retardation distribution in agreement with our previous work [16]. For the sphere indenter case (figure S7 in supplementary materials), the agreement of retardation and azimuth seems satisfactory and similar to the case of the cone. The birefringence response under sphere indenter is relatively identical compared to the case of cone indentation.
<figure><img src="image_11.png" /><figcaption>Figure 10: (A) & (B) Radial plots of the experimental (blue solid lines) and calculated (red solid lines) photoelastic retardation $\Delta$ and azimuth $\phi$ along several directions $\psi$ (C) for the case of SLS glass indented with a cone at a load of 0.98 N. Results calculated with different maximum densification abilities ($\xi_{sat} = 0.01, 0.04, 0.08$)) are also shown to highlight the effect of densification on the birefringence calculation. The distance $D$ was normalized by the contact radius $a$ for each dataset; (D) face-to-face comparison of calculated and experimental retardation and azimuth mappings with dimensions normalized by their respective contact radius.</figcaption></figure>

### Silica glass

Figures 11 and figure S8 (supplementary materials) shows the comparison between birefringence experiments and calculations for the cases of cone and sphere indentation on silica glass. In both cases, the general features of birefringence could be reproduced. However, although the distribution of the retardation outside the plastic zone is modeled correctly, we find that in this case the predicted amplitude is too large. The calculations overall also reproduce the azimuth distributions well. The lower experimental values for retardation in the case of silica suggests that an additional stress relaxation process occurs which is not correctly accounted for in our constitutive relation. Silica glass is well known to densify significantly under indentation; however, the distribution of the densification under an indent is still a matter of debate [40].

In brief, the birefringence experiments demonstrate that FEA based on calibrated constitutive relations provide reasonably accurate predictions of residual stress fields for SLS, and also a correct
distribution for silica glass, although the overall magnitude is overestimated. In contrast, Yoffe's residual field misses most of the experimentally observed features.

<figure><img src="image_12.png" /><figcaption>Figure 11: (A) & (B) Radial plots of the experimental (blue solid lines) and calculated (red solid lines) photoelastic retardation $\Delta$ and azimuth $\phi$ along several directions $\psi$ for the case of silica glass indented with a cone at a load of 0.98 N. The line-profiles distance $D$ were normalized by the contact radius $a$ for each dataset; (C) face-to-face comparison of calculated and experimental retardation and azimuth mappings with dimensions normalized by their respective contact radius.</figcaption></figure>
# 4 Discussion

For the elastic contact (sec. 3.1), the signal-to-noise ratio under in-situ conditions is good and the sample required for the experiments can be relatively large (cubes with a side of 10 mm) and no immersion liquid is required. Our birefringence computations showed remarkable agreement with the experimental results for both silica and SLS glasses (figure 5). Moreover, the results scaled correctly with load and photoelastic constant (supplementary materials figure S6). These results validate the method. The remaining slight differences in intensity and position may have been caused by several factors, including friction (and maybe off-axis tilt) or geometrical inaccuracies (for example in the determination of the surface position).

The measurement of the residual photoelastic patterns around indentation imprints (sec. 3.2) were more challenging. The maximum size for crack-free imprints being of the order of 10 µm, the signal to noise ratio is much less favorable, with retardations in the range of a few tens of nm. Specific experimental strategies have to be applied with small (500 µm) square-section glass fibers immersed in index matching liquids to minimize spurious birefringence and optimize image collection. For glasses that are susceptible to crystallization such as some aluminosilicates, fiber redrawing is impossible. Therefore, our current method is limited to glasses with relatively good glass formation ability. Note also that we have collected top views of the birefringence pattern, imaging along the indentation axis, in addition to the side views shown in this paper. However, the data provided by the top views seemed less informative, compared to the side views: the distribution of retardation and azimuth displayed fewer characteristic features and less contrast between the two glass compositions. Although still valuable, the top view is expected to provide less information than the lateral view, although the experiments are easier to carry out.

As to the calculations, they rely on an elasto-optical relation so that they are expected to apply only to the region outside the plastic zone. Of course, the elastic stress distribution outside the plastic zone depends on the deformation resulting from plastic flow during indentation and so reflects the elastoplastic properties of the glasses. It would be interesting to model the birefringence pattern inside the plastic zone as well. This could be done by measuring the birefringence properties of the plastically deformed material, for example by external calibration on larger samples [Bechgaard17] (sec. 2.2.2) or simply by assuming dependencies of the optical properties on the irreversible deformations. This problem is beyond the scope of the present paper but remains a promising perspective: we already note some similarities between the measured birefringence patterns and the predictions of the standard elasto-optical law based on the elastoplastic stress field through regions involving the plastic zone.

In the regions outside the plastic zone, we find that the overall comparison between the measured and the calculated birefringence patterns (retardation and azimuth) is good for the FEA method. The geometrical features of calculated and measured distributions are similar, both for retardation and azimuth, for the two glasses and the two indenters. This agreement validates the overall approach, including the constitutive relations.

To further assess the potential of the method, we have also calculated the predicted photoelastic patterns during loading and unloading for the case of cone indentation of SLS glass (videos in the supplementary materials). The results reveal that the elastic field due to the applied force (in fact, quite similar to the elastic contact case) dominates the loading and most of the unloading stages, while the residual stress field only becomes evident in the very last stage of the unloading process where the birefringence pattern changes drastically. Clearly, the unloaded indent is not only easier to investigate experimentally: it is also the configuration which is most sensitive to the details of the plastic deformation.

Despite the general similarity in the patterns, we also found differences in the measured photoe-
elastic responses of the two glass compositions (outside the plastic zone), which reflects their different elastoplastic behaviors. The results for SLS were the most satisfactory both in the far field and near the boundary of the plastic zone. The calculations reproduced semi-quantitatively most of the features of the measured distributions. We conclude that the constitutive relation, once adequately calibrated, allowed for a suitable description of the plastic deformation process in SLS glass. The computation for silica glass showed less agreement with experiments, mostly in terms of retardation amplitude. This means the distribution is correct but the level of residual stress found both in cone and sphere experiments was lower than predicted. In fact, for silica, the retardation intensities came closer to the measurement limits set by the birefringence imaging system and background signal in the far field. This lower stress level in silica glass is certainly a manifestation of densification, which reduces plastic flow and limits the build-up of confining stresses. Our current constitutive relation predicts a relatively smooth decrease of density from the center to the boundary of the elastoplastic zone [25]. More recently, several works have challenged this result and evidenced experimentally that the densification distribution in silica glass seems much more homogeneous. Guin et al. [40], using a chemical etching method, showed that silica displays an almost fully saturated plastic zone. Berthelot et al. [41] built up density maps in indented silica glass using Brillouin spectroscopy: their results also suggest that the density distribution is mostly homogeneous in the plastic zone along the loading axis of the indenter, but changes gradually outside in the perpendicular direction. Both works suggest that, in silica glass, the densification distribution has a sharper gradient at its periphery than predicted by our constitutive relation. The lower birefringence level observed in the present work very likely confirms these results and adds to the growing evidence that our current constitutive relation underevaluates the extent of densification under indentation in silica glass.

Improving on the present constitutive relation is not trivial though, because the apparently excessive spread of the densification gradient directly results from the moderate slope in the well established hardening curve under hydrostatic pressure (figure 4). This apparent contradiction can probably be solved if in addition to density, another internal variable such as plastic shear strain, is included and coupled to the densification process. Setting up such a model requires extensive calibration and is beyond the scope of this paper. However, this discrepancy further demonstrates how birefringence data can provide valuable input for more precise quantification of advanced mechanical properties in brittle materials.

Finally, the birefringence patterns predicted based on Yoffe's analytical model [12] proved very different from the measured ones, except in the far field where the blister term indeed seems to produce the quadrupolar-like motive evidenced by the azimuth (figures 7 and 9 bottom). Otherwise, Yoffe's model, as expected by its author ("at a little distance" - sec. 3.2 in ref. [12]), fails to describe the region closer than 2-3 times the contact radius, as already pointed out by Davis et al. [14] in their direct comparison to FEA calculations. Particularly conspicuous in our data, the two "black holes" (in the retardation) and interconnecting vertical segment (in the azimuth) are completely absent in Yoffe's approach. These shortcomings should not come as a surprise because the model is one of the simplest possible, using a single point expansion in so-called force doublets, from which its tractability (and thus popularity) actually derives. In general, Yoffe's model is able to properly catch a number of features under load: the tensile stresses developing around the plastic zone explain the formation of typical ring and median cracks, which arise from the dominating Boussinesq term. In contrast, the representation of the effects of the plastic strain through the blister field remains comparatively poor. Higher order features such as proposed by other analytical models [10, 11] should improve the accuracy of the predictions. A full evaluation of these models is an involved process and beyond the scope of the present paper. In contrast, our results demonstrate that FEA calculations, as common now in mechanics laboratories as they were exotic in the early 80s, provide robust, though still imperfect, descriptions of the effect of plastic flow [25, 26, 29, 42].
# 5 Conclusion

We have demonstrated a novel approach to the measurement of indentation stress fields in transparent brittle materials. This method combines high sensitivity birefringence measurements with 3D linear stress-optical calculations based on FEA predicted stress fields. When successful, the method provides rich sets of data which help quantitatively evaluate the non linear mechanical properties such as plasticity, damage and rupture in brittle materials. After validation on a simple case, namely elastic contacts under load, we have successfully used it to investigate the residual indentation stress fields in two silicate glasses. We have demonstrated that the predicted stress field in SLS glass matches the measured birefringence pattern very well, which suggest that the constitutive relation is suitable and accurately calibrated. This calibration was obtained through independent high pressure micromechanics experiments. However, for amorphous silica, the archetypal anomalous glass, the predicted pattern was consistent with the data overall but the stress level was found to be overestimated, which points to necessary improvements for the standard form of the constitutive relation. This observation ties in with recent works reporting that silica undergoes maximum residual volumetric strain nearly throughout the plastic zone. Finally, Yoffe's model for indentation was found to provide poor predictions for the residual stress field except in the far field. In particular, we expect that this new method will provide valuable information on stresses and damage in silicate glasses and especially in the vicinity of the boundaries of the plastic zone where crack initiation often takes place.
## CRediT authorship contribution statement

**Gustavo Alberto Rosales-Sosa**: Conceptualization, Data curation, Formal Analysis, Investigation, Methodology, Software, Validation, Visualization, Writing - original draft. **Etienne Barthel**: Supervision, Funding Acquisition, Conceptualization, Investigation, Methodology, Software, Validation, Writing - review & editing. **Yoshinari Kato**: Supervision, Conceptualization, Investigation, Methodology, Writing - review & editing. **Matthieu Bourguignon**: Investigation, Software. **Tomiki Inoue**: Investigation. **Akihiro Yamada**: Investigation. **Shingo Nakane**: Conceptualization, Supervision, Funding Acquisition. **Hiroki Yamazaki**: Conceptualization, Supervision, Funding Acquisition, Writing - review & editing.

## Declaration of competing interest

The authors declare that they have no known competing financial interests or personal relationships that could have appeared to influence the work reported in this paper.

## Acknowledgements

This work was funded by Nippon Electric Glass and Agence Nationale de la Recherche, project GaLAaD (project ANR-20-CE08-0002). The present study was partially supported by Grant-in-Aid Scientific Research C (25400517) to Prof. Akihiro Yamada. Physical property measurements under hydrostatic stress conditions related to this work were performed at SPring-8 (Japan) under proposals 2019B1364 and 2021B1368. The authors would like to thank Prof. Takahito Ohmura and Dr. Paul Viola from the National Institute of Materials Science (NIMS) in Japan for the support with the nanoindentation experiments. The authors would like to thank B.Sc. Severin Barthel for the helpful discussion during the Python code implementation. We are indebted to Prof. Satoshi Yoshida for numerous discussions on indentation in glasses and birefringence measurements. The authors would also like to express their deepest gratitude to Prof. Chuck Kurkjian† from Rutgers University for his invaluable comments and inspiration to pursue this research. His legacy will continue to inspire researchers in the field of materials science.

## Supplementary materials

The supplementary materials associated with this article can be found in the online version

## References

[1] B. R. Lawn, A. G. Evans, and D. B. Marshall, Elastic/Plastic Indentation Damage in Ceramics: The Median/Radial Crack System, Journal of the American Ceramic Society 63 (9-10) (Sept. 1980), 574-581, ISSN: 1551-2916, 10.1111/j.1151-2916.1980.tb10768.x.

[2] Robert F Cook and George M Pharr, Direct observation and analysis of indentation cracking in glasses and ceramics, Journal of the American Ceramic Society 73 (4) (1990), 787-817.
[3] D.B. Marshall, R.F. Cook, N. P. Padture, M.L. Oyen, A. Pajares, J.E. Bradby, I.E. Reimanis, R. Tandon, T.F. Page, G.M. Pharr, et al., The compelling case for indentation as a functional exploratory and characterization tool, J. Am. Ceram. Soc. 98 (9) (2015), 2671-2680, https://doi.org/10.1111/jace.13729.

[4] T. Rouxel, J.I. Jang, and U. Ramamurty, Indentation of glasses, Prog. Mater. Sci 121 (2021), 100834, https://doi.org/10.1016/j.pmatsci.2021.100834.

[5] M. Lu B. R. Lawn H. Huang, O. Borrero-López, and Y. Zhang, Threshold damage mechanisms in brittle solids and their impact on advanced technologies, Acta Mater. 232 (2022), 117921, https://doi.org/10.1016/j.actamat.2022.117921.

[6] D.M. Marsh, Plastic flow in glass, P. Roy. Soc. Lond. A Mat 279 (1378) (1964), 420-435, https://doi.org/10.1097/00006534-199010000-00053.

[7] B.R. Lawn and M.V. Swain, Microfracture beneath point indentations in brittle solids, J. Mater. Sci. 10 (1975), 113-122, https://doi.org/10.1007/BF00541038.

[8] K.L. Johnson, The correlation of indentation experiments, J. Mech. Phys. Solids 18 (2) (1970), 115-126, https://doi.org/10.1016/0022-5096(70)90029-3.

[9] K.L. Johnson and K. Langstreth, Contact mechanics, Cambridge University Press, 1987, https://doi.org/10.1017/cbo9781139171731.

[10] S.S. Chiang, D.B. Marshall, and A.G. Evans, The response of solids to elastic/plastic indentation. I. Stresses and residual stresses, J. Appl. Phys. 53 (1) (1982), 298-311, https://doi.org/10.1063/1.329930.

[11] G. Feng, S. Qu, Y. Huang, and W.D. Nix, An analytical expression for the stress field around an elastoplastic indentation/contact, Acta Mater. 55 (9) (2007), 2929-2938, https://doi.org/10.1016/j.actamat.2006.12.030.

[12] E.H. Yoffe, Elastic stress fields caused by indenting brittle materials, Philos. Mag. A 46 (4) (1982), 617-628, https://doi.org/10.1080/01418618208236917.

[13] G.D. Quinn and R.C. Bradt, On the Vickers indentation fracture toughness test, J. Am. Ceram. Soc. 90 (3) (2007), 673-680, https://doi.org/10.1111/j.1551-2916.2006.01482.x.

[14] B.C. Davis, G.Scott Glaesemann, and I. Reimanis, Sharp indentation stress fields in fused silica: Finite element analysis and Yoffe analytic model, J. Am. Ceram. Soc. 103 (12) (2020), 7135-7146, https://doi.org/10.1111/jace.17399.

[15] Y. Kato, H. Yamazaki, S. Yoshida, Satoshi, J. Matsuoka, and M. Kanzaki, Measurements of density distribution around Vickers indentation on commercial aluminoborosilicate and soda-lime silicate glasses by using micro Raman spectroscopy, J. Non-Cryst. Solids 358 (24) (2012), 3473-3480, https://doi.org/10.1016/j.jnoncrysol.2012.04.035.

[16] E. Barthel, V. Keryvin, G. Rosales-Sosa, and G. Kermouche, Indentatioban cracking in silicate glasses is directed by shear flow, not by densification, Acta Mater. (2020), https://doi.org/10.1016/j.actamat.2020.05.011.

[17] W.W. Gerberich, Plastic strains and energy density in cracked plates: Part I, Exp. Mech. 4 (11) (1964), 335-344, https://doi.org/10.1007/BF02323544.

[18] G. Wang and G. Zhao J. Guan, A photo-plastic experimental study on deformation of rotary forging a ring workpiece, J. Mater. Process. Technol. 169 (1) (2005), 108-114, ISSN: 0924-0136, https://doi.org/10.1016/j.jmatprotec.2005.03.003.
[19] C. E. Athanasiou, C. D. Fincher, C. Gilgenbach, H. Gao, W. C. Carter, Y.M. Chiang, and B.W. Sheldon, Operando measurements of dendrite-induced stresses in ceramic electrolytes using photoelasticity, Matter 7 (1) (2024), 95–106, ISSN: 2590-2385, https://doi.org/10.1016/j.matt.2023.10.014.

[20] S. Yoshida, S. Iwata, T. Sugawara, Y. Miura, J. Matsuoka, A. Errapart, and C.R. Kurkjian, Elastic and residual stresses around ball indentations on glasses using a micro-photoelastic technique, J. Non-Cryst. Solids 358 (24) (2012), 3465–3472, https://doi.org/10.1016/j.jnoncrysol.2012.01.069.

[21] J. Anton, A. Errapart, H.K. Aben, and L. Ainola, A discrete algorithm of integrated photoelasticity for axisymmetric problems, Exp. Mech. 48 (5) (2008), 613–620, https://doi.org/10.1299/jsmeatem.2003.2._os01w0035.

[22] C.R. Kurkjian, S.M. Allameh, and S. Yoshida, Radial cracking at Vickers indentations in pristine flat/square silica fibers, Int. J. Appl. Glass Sci. 13 (1) (2022), 11–20, https://doi.org/10.1111/ijag.16382.

[23] H.K. Aben, Optical phenomena in photoelastic models by the rotation of principal axes, Exp. Mech. 6 (1966), 13–22, https://doi.org/10.1007/bf02327109.

[24] H.K. Aben and C. Guillemet, Photoelasticity of glass, Springer Science & Business Media, 2012, https://doi.org/10.1007/978-3-642-50071-8_3.

[25] G. Kermouche, E. Barthel, D. Vandembroucq, and Ph. Dubujet, Mechanical modelling of indentation-induced densification in amorphous silica, Acta Mater. 56 (13) (2008), 3222–3228, https://doi.org/10.1016/j.actamat.2008.03.010.

[26] S. Bruns, K. E. Johanns, H. Rehman, G. Pharr, and K. Durst, Constitutive modeling of indentation cracking in fused silica, J. Am. Ceram. Soc. 100 (5) (2017), 1928–1940, https://doi.org/10.1111/jace.14734.

[27] G. Molnar, Linear Elastic, density hardening plasticity (2015), https://www.molnar-research.com/tutorials_Glass.html.

[28] V. Keryvin, J.-X. Meng, S. Gicquel, J.-P. Guin, L. Charleux, J.-C. Sanglebœuf, P. Pilvin, T. Rouxel, and G. Le Quilliec, Constitutive modeling of the densification process in silica glass under hydrostatic compression, Acta Mater. 62 (2014), 250–257, https://doi.org/10.1016/j.actamat.2013.07.067.

[29] E. Barthel G. Molnar G. Kermouche, Plastic response of amorphous silicates, from atomistic simulations to experiments–A general constitutive relation, Mech. Mater. 114 (2017), 1–8, https://doi.org/10.1016/j.mechmat.2017.07.002.

[30] T. Rouxel, H. Ji, T. Hammouda, and A. A. Moréac, Poisson’s ratio and the densification of glass under high pressure, Phys. Rev. Lett. 100 (22) (2008), 225501, https://doi.org/10.1103/physrevlett.100.225501.

[31] T. Deschamps, A. Kassir-Bodon, C. Sonneville, J. Margueritat, C. Martinet, D. De Ligny, A. Mermet, and B. Champagnon, Permanent densification of compressed silica glass: a Raman-density calibration curve, J. Phys-Condens. Mat. 25 (2) (2012), 025402, https://doi.org/10.1088/0953-8984/25/2/025402.

[32] H.K. Mao, P.M. Bell, J.W. Shaner, and D.J. Steinberg, Specific volume measurements of Cu, Mo, Pd, and Ag and calibration of the ruby R1 fluorescence pressure gauge from 0.06 to 1 Mbar, J. Appl. Phys. 49 (6) (1978), 3276–3283, https://doi.org/10.1063/1.325277.
[33] T. Deschamps, C. Martinet, J.L. Bruneel, and B. Champagnon, Soda-lime silicate glass under hydrostatic pressure and indentation: a micro-Raman study, J. Phys-Condens. Mat. 23 (3) (2011), 035402, https://doi.org/10.1088/0953-8984/23/3/035402.

[34] N. Kawai and S. Endo, The Generation of Ultrahigh Hydrostatic Pressures by a Split Sphere Apparatus, Rev. Sci. Instrum. 41 (8) (Aug. 1970), 1178-1181, https://doi.org/10.1063/1.1684753.

[35] Y. Yokoyama, B.R. Mitchell, A. Nassiri, B.L. Kinsey, Y.P. Korkolis, and Y. Tagawa, Integrated photoelasticity in a soft material: phase retardation, azimuthal angle, and stress-optic coefficient, Opt. Laser Eng. 161 (2023), 107335, https://doi.org/10.2139/ssrn.4171604.

[36] F.C and Anthony C., Introduction to contact mechanics, 101, Springer, 2007, https://doi.org/10.1007/978-0-387-68188-7.

[37] A. Arora, D.B. Marshall, B.R. Lawn, and M.V. Swain, Indentation deformation/fracture of normal and anomalous glasses, J. Non-Cryst. Solids 31 (3) (1979), 415-428, https://doi.org/10.1016/0022-3093(79)90154-6.

[38] T.M. Gross, Deformation and cracking behavior of glasses indented with diamond tips of various sharpness, J. Non-cryst. solids 358 (24) (2012), 3445-3452, https://doi.org/10.1016/j.jnoncrysol.2012.01.052.

[39] J.F.S. Christensen, N.P.A. Krishnan, M. Bauchy, and M.M. Smedskjaer, Indenting glasses with indenters of varying stiffness and sharpness, J. Non-Cryst. Solids 603 (2023), 122111, https://doi.org/10.1016/j.jnoncrysol.2022.122111.

[40] J.-P. Guin, K. Han, L. Charleux, J.-C. Sangleboeuf, M. Ferry, and V. Keryvin, A new nanometer resolution method for probing densification ratio at nanoindentation sites in glass: Unraveling discrepancies in the literature, Acta Mater. (2024), 120005, https://doi.org/10.1016/j.actamat.2024.120005.

[41] A. Berthelot, R. Elodie, X. Dagany, T. Deschamps, E. Herry, R. Kineider, Y. De Leon, B. Emmeline, G. Kermouche, E. Barthel, and C. Martinet, Strength of Brillouin spectroscopy to identify spatial densification model in indented silica, J. Non-Cryst. Solids 639 (2024), 123058, https://doi.org/10.2139/ssrn.4704946.

[42] K. Hashiguchi, H. Yamazki, S. Nakane, Y. Kato, G. Rosales-Sosa, and M. Ueno, Subloading-elastoplastic constitutive equation of glass, J. Mater. Sci. Technol. 185 (2024), 221-232, https://doi.org/10.1007/978-3-030-93138-4_8.