# Surface properties of alkali silicate glasses: Influence of the modifiers

Zhen Zhang, Simona Ispas, Walter Kob

► To cite this version:

Zhen Zhang, Simona Ispas, Walter Kob. Surface properties of alkali silicate glasses: Influence of the modifiers. The Journal of Chemical Physics, 2023, 158 (24), pp.244504. 10.1063/5.0155497. hal-04856797

HAL Id: hal-04856797

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

Submitted on 20 Mar 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>

Distributed under a Creative Commons Attribution 4.0 International License
# Fracture of silicate glasses: Micro-cavities and correlations between atomic-level properties

Zhen Zhang, Simona Ispas, and Walter Kob*

*Laboratoire Charles Coulomb (L2C),*

*University of Montpellier and CNRS, F-34095 Montpellier, France*

*(Dated: May 6, 2022)*

## Abstract

We use large-scale simulations to investigate the dynamic fracture of silica and sodium-silicate glasses under uniaxial tension. The stress-strain curves demonstrate that silica glass is brittle whereas the glasses rich in Na show pronounced ductility. A strong composition dependence is also seen in the crack velocity which is on the order of 1800 m/s for glasses with low Na concentration and decreases to 700 m/s if the concentration is high. We find that during the fracture of Na-rich glasses very irregular cavities as large as 3-4 nm form ahead of the crack front, indicating the presence of nanoductility in these glasses. Before fracture occurs, the local composition, structure, and mechanical properties are heterogeneous in space and show a strong dependence on the applied strain. Further analysis of the correlations between these local properties allows to obtain a better microscopic understanding of the deformation and fracture of glasses and how the local heating close to the crack tip, up to several hundred degrees, permits the structure to relax.

<footnote>* Corresponding author: walter.kob@umontpellier.fr</footnote>
# I. INTRODUCTION

Silicate glasses constitute a class of amorphous materials that is not only of high scientific interest but is also most important for many technical applications [1–3]. For many of these applications, the mechanical reliability of the material is crucial and hence this subject has been in the focus of interest of numerous studies in glass science and technology [4, 5]. Silicate glasses are known to be archetypal brittle materials which show little or no plastic deformation before fracture, although they do display a pronounced non-linear elastic behavior [6–10]. However, whether or not this well-known property is also observed on the micro-and meso-scale has been controversial for decades, as documented in several review articles [11–13]. One of the questions that is highly debated is whether or not the fracture of oxide glasses is accompanied by the formation, growth and coalescence of microscopic cavities in front of the crack tip [14–23], features that are usually associated with ductile fracture. In experiments glass fracture is frequently investigated by means of atomic force microscopy (AFM), often in a stress-corrosive environment, i.e., when water interacts with the glass structure [24]. This experimental setup allows to control the velocity for the crack propagation over a wide range, but usually smaller than $10^{-2}$ m/s (thus sub-critical fracture), and permits the direct observation of crack propagation in situ [25]. Cavities with sizes ranging from 5 nm to more than 100 nm have been reported, and this size depends on various factors, such as the measuring method, stress condition, and crack velocity, although the exact dependence on these factors is not really understood [15–17, 26]. Also molecular dynamics (MD) simulations have been used to study the fracture of silicate glasses, usually in a vacuum environment. These studies reported the existence of voids having sizes around 0.5 nm, weakly depending on glass composition [20–22, 27]. However, the comparison between such simulations and the experimental studies of stress-corrosive fracture has to be done with caution, since these two techniques differ strongly in the relevant parameters such as crack velocities, chemical conditions, strain rates, etc. Hence at the end our microscopic understanding of the fracture dynamics in oxide glasses is not at all in a satisfying state.

A further important question concerns the modification of the local properties of the glass during fracture. Since close to the crack tip the stresses are high, one can expect that in this region the structure and the mechanical properties are very different from the ones in the bulk. Unfortunately experimental techniques do not allow to probe these quantities with
high precision since the process zone near the crack tip in which energy dissipation occurs
has been reported to have a size of less than 10 nm [28, 29], i.e., is at the resolution limit
of typical experimental probes. Recent atomistic simulations have allowed to gain some
insight into atomic-scale properties, such as the dissipation of energy in the vicinity of the
crack tip using the atomistic J-integral approach [30–32], as well as the heterogeneities in
local structure and mechanical properties and the correlations between them [22, 27, 33–35].
However, at present there is still a lack of understanding on how various local properties
evolve and correlate with each other as a function of strain and hence how these correlations
influence the fracture behavior of oxide glasses.

The goal of the present study is to obtain a microscopic understanding of the fracture
behavior of silicate glasses at conditions encountered in real life, i.e., glass breaks when a
surface flaw meets tension. This situation corresponds to dynamic fracture, i.e., the speed
of the crack front is so fast that environment factors, such as humidity, do not influence
the properties of the crack tip and hence the fracture dynamics. For this case it should
therefore be realistic to model the fracture process of the sample in vacuum. We study this
dynamic fracture behavior by means of large-scale MD simulations of sodium silicate glasses
with different Na content. In particular we focus on the possible presence of nanoscale
cavities during fracture, a question which has been debated intensively in the past but for
which so far no conclusive answer has been obtained. Furthermore we determine, for the
case of a modifier-rich glass, the evolution of various atomic-scale properties focusing to
unravel the relationships between local composition, structure, and mechanical properties,
hence obtaining new insight on the fracture mechanisms of silicate glasses on the microscopic
scales.

The rest of the paper is organized as follows. In Sec. II we introduce the simulation
setups. In Sec. III we present and discuss the results regarding the cavitation process and
various atomic-scale properties during deformation. Finally, we summarize and conclude
this work in Sec. IV.

## II. METHODOLOGY

We performed molecular dynamics simulations for a series of sodium silicates samples
with the composition SiO₂ and Na₂O-xSiO₂ (denoted in the following by NSx), with x =
3, 4, 5, 7, 10, and 20, corresponding to the variation of Na₂O from 0 to 25 mole %. The interactions between the atoms were given by the pairwise effective potential proposed by Sundararaman et al. [36, 37]. The functional form of this potential consists of a Buckingham term for short-range interactions and a Coulombic term which is evaluated by the Wolf truncation method [38]. More details regarding the development of the potential and its ability to reproduce experimental data and *ab initio* calculations are given in Refs. [36, 37]. Recent studies have shown that this potential allows for a good quantitative description of the structure, mechanical, as well as surface properties of sodium silicate glasses [39–41], and therefore it can be expected that it is also able to give a realistic description of the fracture process.

For each sample, approximately 2,300,000 atoms were placed randomly in the simulation box which had a fixed volume determined by the experimental value of glass density at room temperature [42], corresponding to boxes with typical dimensions of 20 nm×30 nm×50 nm. Using periodic boundary conditions in three dimensions, these samples were first melted and equilibrated at 6000 K for 80 ps in the canonical ensemble (NVT) and then cooled and equilibrated at a lower temperature T₁ (still in the liquid state) for another 160 ps. The temperature T₁ ranges from 3000 K for SiO₂ to 2000 K for NS3 (25 mole % Na₂O), see Ref. [40, 43] for details. Subsequently we cut the sample orthogonal to the z-axis, and inserted an empty space, thus creating two free surfaces, i.e., the sample had the geometry of a slab. Periodic boundary conditions were applied in all three directions. In order to ensure that the two free surfaces do not interact with each other, the thickness of the vacuum layer varied from 6 nm for silica to 14 nm for NS3. These samples were then equilibrated at T₁ for 1.6 ns, a time span that is sufficiently long to allow the reconstruction of the surfaces and the equilibration of the interior of the samples. Subsequently the liquid samples were cooled via a two-stage quenching: A cooling rate of γ₁ = 0.125 K/ps was used to quench the samples from T₁ to a temperature T₂ and a faster cooling rate γ₂ = 0.375 K/ps to cool them from T₂ to 300 K. Finally, the samples were annealed at 300 K for 800 ps. The temperature T₂ at which the cooling rate changes was chosen to be at least 200 K below the simulation glass transition temperature T_g, i.e., depending on composition, 1500 K ≥ T₂ ≥ 800 K, see Refs. [40, 43] for details. At T₂, we also switched the simulation ensemble from NVT to NPT (at zero pressure) so that the generated glass samples were not under macroscopic stress at room temperature.
After the preparation of the glass samples we introduced on one of its free surfaces a linear “scratch” in the form of a triangular notch spanning the sample in the y-direction of width and depth of 3 nm and 2 nm, respectively. Subsequently we applied to the sample a strain in the x-direction, with a constant rate of 0.5 ns⁻¹, until it broke. Due to the presence of the notch, the place at which the fracture initiated could be changed at will. More details can be found in Ref. [43].

Temperature and pressure were controlled using a Nosé-Hoover thermostat and barostat [44–46]. All simulations were carried out using the Large-scale Atomic/Molecular Massively Parallel Simulator software (LAMMPS) [47] with a time step of 1.6 fs.

The results presented in this study correspond to one melt-quench sample for each composition. However, we emphasize that the system sizes considered are sufficiently large to make sample-to-sample fluctuations negligible. A systematic evaluation of the influence of various simulation parameters (system size, cooling rate, strain rate, etc.) on the mechanical behavior has been performed and documented elsewhere [43]. The results presented in this paper correspond to system sizes for which the fracture behavior have basically become independent of the system size.

### III. RESULTS AND DISCUSSION

#### A. Stress-strain behavior

We first demonstrate that the composition of the glass has a strong influence on the strength and fracture dynamics of the samples. In Fig. 1a we present the stress-strain curves for the different glasses and one sees a clear dependence of these curves on the concentration of sodium. These curves are qualitatively similar to the ones presented in Ref. [10] for the case of bulk systems, but we find that here the failure strain and strength (corresponding to the point of maximum stress) of the glasses are considerably smaller than the corresponding values for the bulk systems [39] or in experiments [48]. This difference can be attributed to the fact that in certain experiments and simulations of the bulk systems one probes the intrinsic values of the failure strain and strength, while here the fracture is strongly affected by the presence of the surface flaw.

The stress-strain curves demonstrate that silica glass is stiff (i.e., high elastic modulus)
<figure><img src="image_2.png" /><figcaption>Figure 1. (a) Stress-strain behavior of the NSx glasses. Note that because of the periodic boundary conditions in the tensile direction, the elastic waves generated by the fracture will bounce back and forth between the two surfaces. This explains the small oscillations in the stress-strain curves after fracture. (b) Crack velocity $v_c$ (filled symbols) and the ratio between $v_c$ and the estimated Rayleigh wave speed $v_R$ (open symbols, right scale).</figcaption></figure>

and brittle (i.e., abrupt stress drop after the maximum strength). The glass becomes softer and more ductile (i.e., can withstand more deformation before failure) upon increasing the concentration of Na₂O. Related to this behavior is the fact that the fracture process becomes more gentle with the increase of sodium concentration as can be quantified via the crack velocity $v_c$. This velocity can be estimated via the expression $v_c = L_c/(t_0 - t_m)$, where $L_c$ denotes the total crack length and $t_0 - t_m$ is the duration of crack propagation, i.e., stress drops from the maximum (at $t_m$) to zero (at time $t_0$). The dependence of $v_c$ on the NaO₂ concentration is shown in Fig. 1b and one sees that $v_c$ decreases from ≈1800 m/s for silica to ≈650 m/s for NS3. Note that the crack velocity depends on various factors such as composition and measurement method and hence it is difficult to make a quantitative comparison between experimental and simulation results. The crack velocities measured in our simulations are on the order of 10³ m/s which is what one typically finds in experiments probing the dynamic fracture of silicate glasses [29, 49–51]. Our data show that for low concentration of Na this dependence is rather mild but it starts to accelerate once the concentration exceeds 10 mol %. A previous study has found that at around this threshold
the glass samples have a significant change in their non-linear elastic properties [10] and hence it can be expected that these changes influence also the magnitude of the crack velocity.

It is instructive to compare the values of $v_c$ with the sound velocities of silicate glasses. For pure silica at room temperature the experimental values for the transverse and longitudinal sound velocities are, respectively, 3743 m/s and 5953 m/s [52] and the ones for sodium silicate glasses are a few percents lower [43, 53]. Hence one concludes that the crack velocities are significantly lower than the speed of sound, but of course still orders of magnitude faster than the sub-critical fracture processes studies in experiments [25]. Furthermore we mention that the crack velocity generally increases with strain rate, but this dependence is weak, although somewhat more pronounced for Na-rich compositions [43].

Within the theory of linear elastic fracture mechanics one finds that mode I (tensile) cracks cannot propagating faster than the Rayleigh wave speed, $v_R$, since this would correspond to a negative energy release rate [54]. To get an idea of how our measured crack velocities compare with $v_R$, we estimate the latter from the expression [49, 55–57]

$$v_R = \frac{0.87 + 1.12\nu}{1 + \nu} \sqrt{\frac{G}{\rho}}, \quad (1)$$

where G is the shear modulus, $\nu$ is Poisson’s ratio, and $\rho$ is density. (In practice we measured the Young modulus E from the slope of the stress-strain curve at small $\epsilon$ and $\nu$ from the geometry of the sample under strain. G was then obtained from $G = E/2(1+\nu)$.) In Fig. 1b we also include the ratio $v_c/v_R$ as a function of the Na₂O content. One sees that this ratio is high for glasses with low Na content and it decreases significantly with the addition of sodium. This dependence is thus consistent with the view that the fracture process of glasses with low Na concentration is dissipating less energy than it is the case for glasses with high Na content, in agreement with the finding that the former glasses are brittle while the latter are more ductile.

Also of interest is the observation that $v_c/v_R$ reaches a maximum at around 9.5% Na₂O (NS10), i.e., the composition at which a significant change in the elastic properties of the sodium silicate glasses is observed [10]. Thus this finding supports the view expressed above that the mechanical response before failure has consequences for the fracture process. We also mention that for the Na-poor compositions, i.e., Na₂O% < 13%, the estimated $v_c/v_R$ is compatible with the value of 0.5 ~ 0.6 that one usually finds in experimental measure-
ments [58], while for the Na-rich glasses, $v_c/v_R$ is considerably smaller. The reduction of $v_c/v_R$ can be related to the fact that the local deviation of the crack front (crack arresting) is more pronounced in these Na-rich glasses, which is also reflected in the change of fracture surface roughness [41]. Finally, we note that all the estimated $v_c/v_R$ are smaller than the experimental upper bound of 0.66 found for inorganic glasses [59, 60].

## B. Cavitation

<figure><img src="image_3.png" /></figure>

Figure 2. Number density (left scale) and number of the voids (right scale) in the unstrained samples versus the radius of the probing sphere, *R*. In the search of voids, the examined samples are slabs of dimensions 20 nm×30 nm×50 nm. The green dashed line represents $N_{void}=0.5$, i.e., with 50% probability to find only one void with size *R* in the unstrained sample. Error bars are standard error of the mean over the six fracture simulations.

To understand the deformation and fracture behavior of the glasses on the microscopic scale, we explored first whether or not glass fracture is accompanied by the formation and growth of cavities. To this aim we need an operational definition of what such a cavity really is since one has to be able to distinguish it from the random voids that are already present in the open network structure of the non-strained glass. To start, we considered all atoms to be
point-like objects, i.e., no volume is assigned to them. Then a 3D cubic grid with a mesh size
of 1.25 Å was superposed to the sample and a probing sphere of radius *R* was moved over all
grid points. This mesh size was several times smaller than the probing sphere (see below),
thus ensuring that the whole space is properly probed. Those grid points for which no atoms
were found inside the probing spherical regions were labeled as empty. If two neighboring
grid points were both classified as empty, the sum of their encompassed probing spherical
regions was defined to form one void. Figure 2 shows that for unstrained glasses the void
density ρ<sub>void</sub>, i.e., the number of voids per unit volume, decreases basically exponentially
with increasing *R*, and that this decrease is independent of the Na concentration. From this
data we recognize that the void density decreases by about a factor of 10 if one increases the
Na content from 0% (silica) to 25% (NS3). This trend can be rationalized to some extent by
the fact that the atomic number density ρ<sub>N</sub> of the glass increases with the addition of Na<sub>2</sub>O
(ρ<sub>N</sub> = 66.82 nm<sup>-3</sup> for silica and 73.83 nm<sup>-3</sup> for NS3). However, this modest difference of
10% is not able to explain the factor of 10 the density of the void size. Instead this large
factor shows that the addition of Na makes that the very open network structure of silica
is replaced by a more uniform structure that has significantly smaller holes, in agreement
with the observation that the addition of Na leads to a strong decrease of the first sharp
diffraction peak in the static structure factor [61].

From the graph one can conclude that, independent of the Na concentration, voids with
a size of *R* ≈ 4.0 Å become very rare, i.e., less than one void in the samples we consider, and
hence we can use this size as a threshold to distinguish between cavities that are already in
the unstrained sample from the ones that are due to the applied stress. We also mention that
in a previous simulation study in which the same definition of void was used, it was found
that the largest void in unstrained silica glass has a radius of ≈ 3.5 Å [18], i.e., significantly
less than the threshold we find here. However, in that work the samples were cubic boxes
with a size of ≈ 3.7 nm, and our data for silica shows that at *R* = 3.5 Å one can expect to
find about 10<sup>-2</sup> voids per nm³, i.e., 0.5 voids in a cube of size 3.7 nm, in good agreement
with the estimate of Ref. [18].

Having determined the upper size of the cavities in the unstrained samples, we can now
characterize the fracture process of the glasses in terms of voids. To start we present in
Figs. 3(a) and (b) snapshots of the system that show the evolution of the crack front for
the silica and Na-rich NS3 glasses, respectively (movies of this process can be found in
<figure><img src="image_4.png" /><figcaption>Figure 3. (a) and (b): Snapshots showing the empty space during the fracture of silica, (a), and NS3, (b), glasses (sample size 20 × 30 × 50 nm³). A notch is introduced on the top surface to initiate the fracture. Color coding is based on the coordinate in the *y*-direction. Movies of the fracture process can be found in Ref. [62]. (c) and (d): Stress-strain curves corresponding to the fracture of the two glasses. The red circles indicate the strains at which the snapshots in panels (a) and (b) are shown. (e) and (f): The total number (left scale) and volume of voids (right scale) during fracture of the two glasses. The red circles highlight the specific strain points shown in the snapshots of panels (a) and (b). The insets in panels (e) and (f) show the ratio between the volume (in nm³) and number of the voids.</figcaption></figure>
Ref. [62]). Panels (c) and (d) depict the stress-strain curves of the two glasses and the strains corresponding to the snapshots are marked by circles. The snapshots demonstrate that silica shows typical brittle behavior in that the fracture path is straight and no relevant cavities are found in front of the crack front during the fracture process. (We note that the white dots seen in the snapshots are transient voids which have a size of $R \approx 4$ Å, i.e., they should be considered as local fluctuations.). To quantify this last observation, we present in Fig. 3(e) the number of voids, $N_{void}$, as a function of the strain (left scale) and one recognizes that this number is small and fluctuates strongly. Also included in the graph is the total volume of the voids, $V_{void}$ (right scale), and one sees that this quantity tracks $N_{void}$ perfectly well. The fact, see inset, that the ratio between $V_{void}$ and $N_{void}$ is around $0.286 \text{ nm}^3$ during the entire fracture process (except for the two jumps at around 10% strain due to local structural fluctuations) indicates that all the voids have a size of $R \approx 4.1$ Å, i.e., they are isolated cavities that have a volume just above the cutoff threshold used to define non-trivial voids, in agreement with the observation from the snapshots. Hence one can conclude that silica glass breaks in a nearly perfect brittle manner by sequential bond breaking at the crack tip, i.e., no cavities are formed ahead of the crack tip. This finding agrees with the experimental results by Guin and Wiederhorn [16] and a recent simulation study using reactive force field [32].

The scenario is very different for glasses that are rich in Na, such as NS3 shown in panels (b), (d), and (f): Before the crack starts to propagate, i.e., $\varepsilon < 0.16$, no stable voids are detected. However, once the fracture process has started, one observes the formation of small cavities close to the crack tip and subsequently these voids grow and merge with the crack front (see panel for $\varepsilon = 0.172$). The crack tip also shows micro-branching (see panel for $\varepsilon = 0.192$) and this gives rise to a rougher fracture surface. Panel (f) confirms quantitatively the observations from the snapshots: For low strain the number of voids is small and the ratio $V_{void}/N_{void}$, see inset, is basically independent of $\varepsilon$. Once $\varepsilon$ exceeds 0.15, i.e., the crack starts to propagate, the number of voids quickly shoots up and so does the ratio $V_{void}/N_{void}$, indicating that the typical void grow in size. Hence, one concludes that the fracture of NS3 is accompanied by the formation of nanoscale cavities ahead of the crack front, which subsequently grow and coalescence with the fracture surface. The presence of these cavities is also one of the reasons why the resulting fracture surface has a roughness that is significantly higher than the one found for silica [41].
Although all silica-based glasses are usually considered to be prime examples for brittle materials, the observed differences between the silica and NS3 glasses indicate that the composition does have a crucial role for the fracture behavior on the microscopic scales, in accordance with previous simulations [21, 22]. Hence on small length scales silicate glasses with a significantly depolymerized structure can be considered as quasi-brittle materials, i.e., they are not brittle on all length scales [63]. (We note here that brittleness will depend to some extent on the quench rate with which the samples has been produced as well as on the strain rate used for the fracture. However, these dependencies are relatively mild and do not affect the present results [43].)

A further important question concerns the maximum size of the cavities and their shape. To address these points we have identified the largest void during the fracture of the NS3 glass, see Fig. 4. It is found that the biggest void has in the x-direction (plane parallel to the crack plane, i.e., perpendicular to the tension) an extent of ≈3.5 nm, while in the y-direction (parallel to the stress direction) the length is only ≈1.5 nm. This highly anisotropic shape indicates that once a small cavity has formed somewhat ahead of the crack front (since this is where the local stress is highest), it expands in the direction orthogonal to the stress, thus forming a (roughly) lens-like structure. This structure grows and subsequently merges with the approaching fracture front. Note that once the voids have been absorbed by the passing crack front, it can be expected that they relax elastically in the y-direction, at least partially, since the stress is released, i.e., their depth (as seen from the surface) will decrease [17]. This recoil might explain why no remnants of voids were detected in an experimental study which compared the two postmortem fracture surfaces of soda-lime-silicate glasses [16].

### C. Atomic-level properties

A better understanding of the deformation and fracture behavior of the glass can be obtained by exploring how various atomic-scale properties evolve as a function of strain and how they are related among each other. In this section, we therefore focus on the atomic-scale quantities concerning local compositions, structure, and mechanical properties. More specifically, the local properties considered and their definitions are the following ones: (i) The local number fraction of Na, $f_{\text{Na}} = N_{\text{Na}}/N_{\text{tot}}$, where $N_{\text{Na}}$ and $N_{\text{tot}}$ are, respectively, the number of Na atoms and the total number of atoms within a given probe volume (defined
<figure><img src="image_5.png" /><figcaption>Figure 4. Enlarged view of the largest void found at $\epsilon = 0.196$ for the fracture of a NS3 glass.</figcaption></figure>

below). We recall that in NS3 the stoichiometric value of $f_{Na}$ is equal to 0.167; (ii) The local network connectivity as represented by the coordination number of oxygen, $Z_O$, i.e., the number of Si in the nearest neighbor shell of a given O; (iii) The local variation of the inter-tetrahedral angle, $\Delta\theta_i(\epsilon) = \theta_i(\epsilon) - \theta_i(0)$, where $\theta_i(0)$ is the SiOSi angle of O atom $i$ in the initial configuration (0% strain), and $\theta_i(\epsilon)$ is the SiOSi angle of the same O atom in the sample at strain $\epsilon$. (Note that it is possible that a SiOSi bridge present at 0% strain does not exist anymore at strain $\epsilon$, since one of its Si-O bonds has broken. In this case, we simply assign $\Delta\theta_i(\epsilon) = 0$); (iv) The atomic shear strain as first introduced by Shimizu et al. [64] for measuring the local plastic deformation, i.e., a quantity that carries basically the same information as the non-affine squared displacement $D_2^{\min}$ introduced by Falk and Langer [65]; (v) The atomic stress defined via the per-atom stress tensor [66]. (vi) The local temperature as obtained by using the equipartition theorem to map the instantaneous kinetic energy of the atoms in the probing volume to a temperature.

In order to obtain three dimensional maps of these atom-based quantities we coarse-grained them by mean of a Gaussian weight function $\phi(r) = A(\omega) \exp[-r^2/\omega^2]$, where $r$ is the distance from the center of the probing sphere and $A(\omega)$ is a normalization factor. We have examined the influence of the coarse-graining length-scale $\omega$ on the maps and found that $\omega = 8\sqrt{2}$ Å gives stable results that are also consistent with previous simulation studies [27, 67]. In order to obtain a visual map that matches the glass matrix, it is useful to get rid of the empty spaces. To this aim we have used the mentioned coarse-graining procedure to get a map of the mass density and then applied a threshold of 1.5 g/cm$^3$ to
define the glass-vacuum interface. This procedure removes close to the surface a layer of
about 0.3 nm and the resulting maps match well the shape of the deformed glass sample
(see also Ref. [43]).

Figure 5 shows the maps of the various local properties at two representative strains during the fracture of the NS3 glass, i.e., ε = 0.172 (the fracture has just started) and ε = 0.192 (the fracture front has reached the middle of the sample). (Movies showing the evolution of these local properties during the entire deformation and fracture process are included in Ref. [62].) Firstly, one observes that the local compositional, structural, and mechanical quantities exhibit different degrees of heterogeneities, in that their extent and shape depends strongly on the observable considered. Secondly, a comparison of the snapshots demonstrates that certain properties are correlated to each other. For instance, the local Na concentration $f_{Na}$ and the local coordination number of O, $Z_O$, are clearly anti-correlated, and in the vicinity of the crack tip a large change of the SiOSi angle, $\Delta\theta$, seems to be correlated with the large tensile stress. Below we will discuss the details of these correlations as a function of the strain. The last column of Fig. 5 shows the map of the (kinetic) temperature $T$. One recognizes that close to the crack tip this temperature rises to about 700 K, i.e. several hundred degrees above the average temperature of the sample. This marked local heating, which extends over a distance of 3-7 nm from the fracture surface, is due to the release of the potential energy that was stored in the bonds before fracture and with increasing time this energy will diffuse into the sample. In the context of Fig. 6 we will discuss this process in more detail but already here we mention that this marked local heating can be expected to allow the created surface to relax at least to some extent and thus to become smoother.

To quantify the properties of these local observables, we investigate their distributions as a function of strain. Figure 6(a) shows that the distributions of local Na concentration and network connectivity are basically independent of the strain, implying that the topology of the local atomic configuration is not changing even if the sample is put under high stress. This observation is in agreement with the results from Ref. [10] where it was shown that up to high strains not only the local topology of the network does not change but that also the alkali atoms are displaced remarkably little (< 3 Å), i.e., they will not escape from their surrounding cage. In contrast to this, Fig. 6(b) shows that the distribution of the local shear strain $\eta^{\text{Mises}}$ depends strongly on the global strain. For small $\epsilon$ the distribution has a narrow
<figure><img src="image_6.png" /><figcaption>Figure 5. Maps of various local properties (see labelling of the panels) at (a) $\varepsilon = 0.172$ and (b) $\varepsilon = 0.192$ for the NS3 glass. Sample size is 20 × 30 × 50 nm³. See the text for the definition of each quantity. These maps are shown for the middle plane of the simulation box in the direction orthogonal to the crack front.</figcaption></figure>

peak at small values of $\eta^{\text{Mises}}$ (note the logarithmic scale) and it shifts to larger values if $\varepsilon$ is increased, indicating that the sample is increasingly more shear deformed. Once the crack starts to propagate, i.e., $\varepsilon \ge 0.16$, one finds that the distribution has a noticeable tail at large $\eta^{\text{Mises}}$ which becomes increasingly pronounced with the progression of the crack. This tail arises from the large strains near the fracture surfaces (see the maps in Fig. 5) and thus it can be expected that the weight of this tail is directly proportional to the surface area of the fractured surface. From the graph we also recognize that, once the crack has started to propagate, the main peak of the distribution at around $\eta^{\text{Mises}} = 0.12$ does not depend anymore strongly on $\varepsilon$. This indicates that there is a maximum amount of local strain that the sample can support before it fractures. Furthermore we notice that in the same range of $\varepsilon$, the distribution is non-zero only if $\eta^{\text{Mises}}$ is larger than 0.05, demonstrating that each part of the sample has significant non-affine strain.

Figure 6(c) shows the distribution of local tensile stress for different strains. For $\varepsilon = 0$ we find that this distribution is symmetric around zero (local) stress and has a full width at half maximum of approximately 5 GPa, indicating that the local stresses are surprisingly
<figure><img src="image_7.png" /><figcaption>Figure 6. Distribution of various local properties in NS3 at various values of the strain. (a) Atomic number fraction of Na (left/bottom axis) and coordination number of O (right/top axis). (b) Atomic shear strain $\eta^{\text{Mises}}$. (c) Atomic tensile stress. (d) Kinetic temperature $T$. Note that for this analysis the layers, ≈ 2.5 nm in thickness, near the free surfaces of the sample were not taken into account for calculating the distributions. In panels (b)-(d) “FP” denotes the failure point.</figcaption></figure>

high. If the strain is increased, the distribution shifts to the right and at the failure point has a maximum at around 6 GPa. Note that this value is significantly higher than the maximum global stress in the sample, which is, see Fig. 3(d), around 4.5 GPa. The location of the maximum in the distribution of the local stress does, however, coincide with the intrinsic failure stress of NS3 as obtained from previous simulation studies using periodic boundary conditions in all three directions, i.e., without a free surface [10, 39]. In other words, these observations show that the reduction of the failure stress measured here is due to the presence of a surface defect of the sample, i.e., the notch, while the intrinsic failure stress is significantly higher.

The graph also shows that the shift of the distribution to higher stresses with increasing
$\epsilon$ is accompanied with a widening of the distribution, i.e., the local stresses become more heterogeneous, thus reflecting the disorder of the sample on the microscopic scale. Once the fracture begins, one observes that the peak position at large stress remains unchanged whereas the peak intensity decreases and the distribution starts to show a weak tail toward small stresses. As the crack advances, this tail is transformed quickly into another peak at small stress and its height increases. Since the maximum of this new peak is located at the same stress as the distribution for $\epsilon = 0$, we conclude that it corresponds to the local stresses in those parts of the sample at which the crack has already passed, i.e., which are no longer under stress, in agreement with Fig. 5(b). For the highest values of the strain the distribution is therefore identical to the one at $\epsilon = 0$, i.e., a globally stress-free state, since the sample has been broken.

The distribution of the local kinetic temperature is shown in panel Fig. 6(d). For small strains one finds the expected Gaussian distribution peaking at $T = 300$ K, i.e., the temperature of the sample. The shape of the distribution changes only once the strain exceeds the fracture point in that it starts to show a pronounced tail to the right. From the temperature map in Fig. 5 we recognize that this tail is due to the atoms that are close to the crack front and whose kinetic energy increases rapidly because of the release of the potential energy that was stored in the highly stretched bonds. Note that at these high strains it is not possible to describe the total distribution as a sum of two (or few) Gaussian distributions since one probes here a system that is highly out of equilibrium and therefore the notion of temperature is not well-defined. So the temperature map we show in Fig. 5 should be interpreted with some caution. (However, the distribution of the kinetic energy is a well defined quantity.) Only at very large strains, i.e., once the sample is broken, the distribution starts to return back to its original Gaussian shape and it is possible to define again a temperature. The distribution in Fig. 6(d) also shows that for strains that are just above the failure strain, e.g. $\epsilon = 0.192$, the location of the main peak is at around 270 K, i.e., significantly below the temperature of the unperturbed sample, in agreement with the fact that in the temperature map in Fig. 5 the regions far away from the crack front do have a low $T$. This result might look odd since for distances far away from the crack front one expects that the local temperature is the mean temperature of the sample. The explanation for this somewhat unexpected behavior is the fact that in the simulation we use a thermostat which makes that the total kinetic energy of the particles is fluctuating around the fixed target
value. Thus the fact that a small portion of the sample becomes very hot, will automatically make that the rest of the sample will have an average kinetic energy, and thus temperature, that is below the one of the target value, thus explaining why the regions of the sample far away from the crack front are colder than expected. We emphasize, however, that this finite size effect does not influence in a significant manner the interpretation of the temperature map shown in Fig. 5.

<figure><img src="image_8.png" /><figcaption>Figure 7. Spearman's correlation coefficient, $r_s$, between various local properties of the NS3 glass, see labeling of curves for details. The vertical dashed lines indicate the failure point (i.e., at maximum stress). Error bars are standard error of the mean of 10 planes that are evenly spaced orthogonal to the direction of the crack front, i.e., orthogonal to the y-axis.</figcaption></figure>
Further insight into the mechanical properties of the glass can be obtained by quantifying
the strain dependence of the correlations between the local quantities. Since the propagation
of the crack depends strongly on the local structure of the material structure close to the
tip, we will focus on this region for characterizing the correlations between the various local
quantities. In practice we considered a plane of the sample orthogonal to the x-axis, i.e.,
orthogonal to the crack front, and probed a circular region around the crack tip having a
radius of 5.0 nm, i.e., a size that is comparable to the extent of the process zone as estimated
in previous dynamic fracture experiments of silicate glasses [29]. (We note that the size of
this zone, as determined from the temperature increase close to the crack tip, depends only
very weakly on glass composition and crack velocity, i.e., the chosen size is not a critical
parameter for our analysis [43].)

To quantify the correlation between various local properties we use the Spearman rank
correlation coefficient, r<sub>s</sub>(X, Y), which is a measure for the statistical dependence between
two variables X and Y [68]. By construction r<sub>s</sub> is in the range [-1, 1] and r<sub>s</sub> > 0 means
positive correlation, i.e., Y tends to increase when X increases, and r<sub>s</sub> < 0 indicates an
anti-correlation. The magnitude of r<sub>s</sub> gives the degree of correlation, with r<sub>s</sub> = 0 indicating
no correlation between the two variables.

To start, we show in Fig. 7(a) the correlation between local elemental concentrations.
One observes that Si and Na concentrations are strongly anti-correlated in that r<sub>s</sub>(Si, Na) ≈
-0.94. A strong anti-correlation is also found between O and Na with r<sub>s</sub>(O, Na) ≈ -0.90,
in agreement with previous studies on soda-silicate glasses [67]. From panel (a) one also
recognizes that f<sub>Na</sub> and Z<sub>O</sub> are strongly anti-correlated, i.e., Na-rich regions have a smaller
network connectivity in agreement with the expectation that the presence of an alkali atom
leads to the breaking up of the network. A slight reduction of the correlations can be seen
during the fracture process, which is likely due to the enhanced local disorder near the crack
tip, arising from the local stress-and heat-driven rearrangement of atomic positions. Overall,
we note that the curves in panel (a) are basically independent of the strain (both before and
during fracture), which indicates that the local compositions and network connectivity does
not change in a significant manner even if the glass is under strong deformation.

Figure 7(b) shows the correlations between local Na concentration and mechanical prop-
erties. Overall, we note that f<sub>Na</sub> is positively correlated with local shear strain, and anti-
correlated with local tensile stress. The former result can be rationalized by recalling that
Na-rich regions are more flexible due to the high mobility of Na and the reduced connec-
tivity of the network. Due to this flexibility these regions can deform relatively easily and
hence reduce stress, in agreement with the observed anticorrelation. We point out, however,
that the absolute value of $r_s$ that quantifies the correlation between $f_{Na}$ and $\eta^{\text{Mise}}$ is quite
small. This indicates that the local plastic deformation cannot be solely attributed to the
enrichment of Na since it depends also on other factors such as, e.g., the rotational motion
of a tetrahedron. In contrast to this, the correlation between $f_{Na}$ and local tensile stress
has a substantially stronger correlation, $r_s \approx -0.6$ at the early stage of the deformation and
$r_s \approx -0.5$ during fracture, in agreement with the maps shown in Fig. 5. Thus this indicates
that the local stress level is mainly related to the flexibility of the structure which in turn
is strongly influenced by the presence of Na.

The correlation between the local change of the SiOSi angle and the local tensile stress is
shown in Fig. 7(c). We observe that the correlation between these two quantities gradually
builds up (basically linearly), with $r_s$ increasing from close to zero at small strains to around
0.55 for $\epsilon$ beyond the failure strain. The strain dependence of this correlation can be un-
derstood by considering the coexistence of multiple deformation modes during tension: As
it was shown in Ref. [69] by probing the vibrational density of states, the softest modes are
the rotations of the tetrahedra and the change in the SiOSi angle. Hence it can be expected
that at low and intermediate strains these two modes become active and therefore one finds
an increasing correlation between stress and the change in the SiOSi angle. In contrast to
this the change in the Si-O bond length with increasing strain is small, see Ref. [10], and
hence one does not expect a significant correlation of this length with the value of the stress.
Thus, the SiOSi angle is likely the simplest local structural quantity that is sensitive to the
local stress.

Finally we present in Fig. 7(d) the correlation between local shear strain and tensile stress
and one finds that at small and intermediate strains there is basically no correlation between
these two quantities. This result can be rationalized by recalling that on the microscopic
scale the mechanical properties of a glass are quite heterogeneous [10, 22, 27, 67], i.e., there
are soft spots as well as regions that are relatively stiff, see Fig. 5. One can expect that the
application of an external stress will make that the soft regions will respond easily to an
applied stress, thus releasing it by rearranging the particles, i.e., the local strain is large. In
contrast to this, stiff regions will not change their geometries, i.e., the strain is small. Since,
however, at the interface of these two type of regions the local stresses need to be, more or less, balanced for the sake of mechanical stability, at the end the residual stress in the soft regions are not significantly smaller than the ones in the stiff regions, making that there is no substantial correlation between stress and strain, in agreement with the data in panel (d). This argument is no longer valid once the sample has reached its failure point, since in this case the presence of a propagating crack makes that there is no longer a local stress balance. As a result, one finds that there is a notable anti-correlation between local stress and shear strain with a correlation coefficient of around -0.35 during crack propagation.

## IV. SUMMARY AND CONCLUSIONS

To summarize, we have investigated the dynamic fracture of sodium silicate glasses under tension by using large scale molecular dynamics simulations. This study highlights the strong influence of the composition on the fracture behavior of glasses on the microscopic scales. The stress-strain curves demonstrate that while silica glass shows a nearly perfect brittle fracture behavior, the Na-rich glasses do exhibit a certain degree of ductility during the deformation and fracture process. A detailed analysis of void formation during the fracture process reveals that no growth and coalescence of voids are observed for silica glass, while the fracture of the Na-rich NS3 glass is accompanied by the nucleation of nanoscale voids, which subsequently grow and merge ahead of the crack tip. We have found that, for the NS3 glass, the linear dimension of the voids depend on the direction and that their largest extent is on the order of several nanometers. Further insight into the fracture of the glasses has been obtained by exploring the spatial and temporal evolution of various atomic-level properties and the correlations between them. By probing the evolution of local kinetic temperature during fracture, we find that close to the crack tip the kinetic temperature rises several hundred degrees above the average temperature of the sample, an increase that is observed in a region of size ≈ 10 nm around the fracture tip.

As a final remark we note that the fracture of glass is obviously a subject of high complexity since many factors, such as the length scale under investigation, the strain rate, and the chemical composition/microstructure, play an important role. This complexity might make it impossible to come up with universal conclusions regarding the fracture behavior of silicate glasses. Thus, although the present work provides a microscopic perspective on
the dynamic fracture behavior of sodo-silicate glasses, it remains to be explored how other
glass-formers used in practical applications, such as boro-silicates or alumino-silicates, will
fracture on the microscopic scale since the latter systems have several kinds of network-
formers which can be expected to influence the fracture dynamics. Further studies on more
complex chemical compositions and mechano-chemical conditions will also be valuable to
obtain a more comprehensive understanding of the failure of glasses, and for the future
design of glassy materials with improved mechanical performance.

V. ACKNOWLEDGMENTS

Z.Z. acknowledges financial support by China Scholarship Council (Grant No. 201606050112).
W.K. is member of the Institut Universitaire de France. This work was granted access to the
HPC resources of CINES under the allocations (Grant Nos. A0050907572 and A0070907572)
attributed by GENCI (Grand Equipement National de Calcul Intensif).

Z.Z. is currently at College of Materials Science and Engineering, Xi'an Jiaotong Univer-
sity.

[1] K. Binder and W. Kob. *Glassy materials and disordered solids: An introduction to their statistical mechanics*. World Scientific, 2011.

[2] J. E. Shelby. *Introduction to glass science and technology*. Royal Society of Chemistry, 2005.

[3] A. K. Varshneya. *Fundamentals of inorganic glasses*. Elsevier, 2013.

[4] L. Wondraczek, J. C. Mauro, J. Eckert, U. Kühn, J. Horbach, J. Deubener, and T. Rouxel. "Towards Ultrastrong Glasses." *Adv. Mater.*, **23**, 4578–4586, 2011.

[5] L. Wondraczek, E. Bouchbinder, A. Ehrlicher, J. C. Mauro, R. Sajzew, and M. M. Smedskjaer. "Advancing the mechanical performance of glasses: perspectives and challenges." *Adv. Mater.*, **34**, 2109029, 2022.

[6] F. Mallinder, B. Proctor, et al. "Elastic constants of fused silica as a function of large tensile strain." *Phys. Chem. Glas.*, **5**, 91–103, 1964.

[7] J. T. Krause, L. R. Testardi, and R. N. Thurston. "Deviations from linearity in the dependence of elongation on force for fibers of simple glass formers and of glass optical lightguides." *Phys.*
*Chem. Glas.*, **20**, 135–9, 1979.

[8] P. K. Gupta and C. R. Kurkjian. "Intrinsic failure and non-linear elastic behavior of glasses." *J. Non-Cryst. Solids*, **351**, 2324–2328, 2005.

[9] M. Guerette, C. R. Kurkjian, S. Semjonov, and L. Huang. "Nonlinear Elasticity of Silica Glass." *J. Am. Ceram. Soc.*, **99**, 841–848, 2016.

[10] Z. Zhang, S. Ispas, and W. Kob. "Origin of the non-linear elastic behavior of silicate glasses." *Acta Mater.*, **231**, 117855, 2022.

[11] M. Tomozawa. "Fracture of glasses." *Annu. Rev. Mater. Sci.*, **26**, 43–74, 1996.

[12] M. Ciccotti. "Stress-corrosion mechanisms in silicate glasses." *J. Phys. D*, **42**, 214006, 2009.

[13] S. M. Wiederhorn, T. Fett, J.-P. Guin, and M. Ciccotti. "Griffith Cracks at the Nanoscale." *Int. J. Appl. Glass Sci.*, **4**, 76–86, 2013.

[14] C. L. Rountree, R. K. Kalia, E. Lidorikis, A. Nakano, L. Van Brutzel, and P. Vashishta. "Atom-istic aspects of crack propagation in brittle materials: Multimillion atom molecular dynamics simulations." *Annu. Rev. Mater. Res.*, **32**, 377–400, 2002.

[15] F. Célarié, S. Prades, D. Bonamy, L. Ferrero, E. Bouchaud, C. Guillot, and C. Marlière. "Glass Breaks like Metal, but at the Nanometer Scale." *Phys. Rev. Lett.*, **90**, 075504, 2003.

[16] J.-P. Guin and S. Wiederhorn. "Fracture of Silicate Glasses: Ductile or Brittle?" *Phys. Rev. Lett.*, **92**, 215502, 2004.

[17] D. Bonamy, S. Prades, C. L. Rountree, L. Ponson, D. Dalmas, E. Bouchaud, K. Ravi-Chandar, and C. Guillot. "Nanoscale damage during fracture in silica glass." *Int. J. Fract.*, **140**, 3–14, 2006.

[18] K. Muralidharan, J. Simmons, P. Deymier, and K. Runge. "Molecular dynamics studies of brittle fracture in vitreous silica: Review and recent progress." *J. Non-Cryst. Solids*, **351**, 1532–1542, 2005.

[19] C. Rountree, S. Prades, D. Bonamy, E. Bouchaud, R. Kalia, and C. Guillot. "A unified study of crack propagation in amorphous silica: Using experiments and simulations." *J. Alloys Compd.*, **434**, 60–63, 2007.

[20] A. Pedone, G. Malavasi, M. C. Menziani, U. Segre, and A. N. Cormack. "Molecular Dynamics Studies of Stress-Strain Behavior of Silica Glass under a Tensile Load." *Chem. Mater.*, **20**, 4356–4366, 2008.

[21] A. Pedone, M. C. Menziani, and A. N. Cormack. "Dynamics of Fracture in Silica and Soda-
Silicate Glasses: From Bulk Materials to Nanowires." *J. Phys. Chem. C*, **119**, 25499-25507, 2015.

[22] B. Wang, Y. Yu, Y. J. Lee, and M. Bauchy. "Intrinsic Nano-Ductility of Glasses: The Critical Role of Composition." *Front. Mater.*, **2**, 11, 2015.

[23] L.-Q. Shen, J.-H. Yu, X.-C. Tang, B.-A. Sun, Y.-H. Liu, H.-Y. Bai, and W.-H. Wang. "Observation of cavitation governing fracture in glasses." *Sci. Adv.*, **7**, eabf7293, 2021.

[24] S. M. Wiederhorn, J.-P. Guin, and T. Fett. "The Use of Atomic Force Microscopy to Study Crack Tips in Glass." *Metall. Mater. Trans. A*, **42**, 267-278, 2011.

[25] M. Ciccotti and M. George. "In Situ AFM Investigations and Fracture Mechanics Modeling of Slow Fracture Propagation in Oxide and Polymer Glasses." In W. Andreoni and S. Yip, editors, *Handbook of Materials Modeling: Applications: Current and Emerging Materials*. Springer International Publishing, Cham, 2018.

[26] T. Waurischk, S. Reinsch, T. Rouxel, H. Behrens, J. Deubener, and R. Müller. "Vacuum crack growth in alkali silicate glasses." *J. Non-Cryst. Solids*, **572**, 121094, 2021.

[27] B. Wang, Y. Yu, M. Wang, J. C. Mauro, and M. Bauchy. "Nanoductility in silicate glasses is driven by topological heterogeneity." *Phys. Rev. B*, **93**, 064202, 2016.

[28] K. Han, M. Ciccotti, and S. Roux. "Measuring nanoscale stress intensity factors with an atomic force microscope." *EPL*, **89**, 66003, 2010.

[29] G. Pallares, C. L. Rountree, L. Douillard, F. Charra, and E. Bouchaud. "Fractoluminescence characterization of the energy dissipated during fast fracture of glass." *EPL*, **99**, 28003, 2012.

[30] J. M. Rimsza, R. E. Jones, and L. J. Criscenti. "Crack propagation in silica from reactive classical molecular dynamics simulations." *J. Am. Ceram. Soc.*, **101**, 1488-1499, 2018.

[31] J. M. Rimsza, R. E. Jones, and L. J. Criscenti. "Chemical effects on subcritical fracture in silica from molecular dynamics simulations." *J. Geophys. Res. Solid Earth*, **123**, 9341-9354, 2018.

[32] S. C. Chowdhury, E. A. Wise, R. Ganesh, and J. W. Gillespie. "Effects of surface crack on the mechanical properties of Silica: A molecular dynamics simulation study." *Eng. Fract. Mech.*, **207**, 99-108, 2019.

[33] T. Hao and Z. M. Hossain. "Atomistic mechanisms of crack nucleation and propagation in amorphous silica." *Phys. Rev. B*, **100**, 014204, 2019.

[34] E. J. Frankberg, J. Kalikka, F. G. Ferré, L. Joly-Pottuz, T. Salminen, J. Hintikka, M. Hokka,
S. Koneti, T. Douillard, B. Le Saint, et al. "Highly ductile amorphous oxide at room temperature and high strain rate." *Science*, **366**, 864–869, 2019.

[35] T. Du, H. Liu, L. Tang, S. S. Sørensen, M. Bauchy, and M. M. Smedskjaer. "Predicting Fracture Propensity in Amorphous Alumina from Its Static Structure Using Machine Learning." *ACS nano*, **15**, 17705–17716, 2021.

[36] S. Sundararaman, L. Huang, S. Ispas, and W. Kob. "New optimization scheme to obtain interaction potentials for oxide glasses." *J. Chem. Phys.*, **148**, 194504, 2018.

[37] S. Sundararaman, L. Huang, S. Ispas, and W. Kob. "New interaction potentials for alkali and alkaline-earth aluminosilicate glasses." *J. Chem. Phys.*, **150**, 154505, 2019.

[38] D. Wolf, P. Keblinski, S. R. Phillpot, and J. Eggebrecht. "Exact method for the simulation of Coulombic systems by spherically truncated, pairwise r-1 summation." *J. Chem. Phys.*, **110**, 8254–8282, 1999.

[39] Z. Zhang, S. Ispas, and W. Kob. "The critical role of the interaction potential and simulation protocol for the structural and mechanical properties of sodosilicate glasses." *J. Non-Cryst. Solids*, **632**, 119895, 2020.

[40] Z. Zhang, S. Ispas, and W. Kob. "Structure and vibrational properties of sodium silicate glass surfaces." *J. Chem. Phys.*, **153**, 124503, 2020.

[41] Z. Zhang, S. Ispas, and W. Kob. "Roughness and scaling properties of oxide glass surfaces at the nanoscale." *Phys. Rev. Lett.*, **126**, 066101, 2021.

[42] N. P. Bansal and R. H. Doremus. *Handbook of glass properties*. Orlando: Academic Press, 1986.

[43] Z. Zhang. *Fracture, surface, and structure of silicate glasses: Insights from atomistic computer simulations*. Doctoral dissertation, University of Montpellier, 2020.

[44] S. Nosé. "A unified formulation of the constant temperature molecular dynamics methods." *J. Chem. Phys.*, **81**, 511–519, 1984.

[45] W. G. Hoover. "Canonical dynamics: Equilibrium phase-space distributions." *Phys. Rev. A*, **31**, 1695–1697, 1985.

[46] W. G. Hoover. "Constant-pressure equations of motion." *Phys. Rev. A*, **34**, 2499–2500, 1986.

[47] S. Plimpton. "Fast parallel algorithms for short-range molecular dynamics." *J. Comput. Phys.*, **117**, 1–19, 1995.

[48] N. P. Lower, R. K. Brow, and C. R. Kurkjian. "Inert failure strain studies of sodium silicate
glass fibers.” *J. Non-Cryst. Solids*, **349**, 168–172, 2004.

[49] G. D. Quinn. "On terminal crack velocities in glasses." *Int. J. Appl. Glass Sci.*, **10**, 7–16, 2019.

[50] G. D. Quinn and J. J. Swab. "Fracture toughness of glasses as measured by the SCF and SEPB methods." *J. Eur. Ceram. Soc.*, **37**, 4243–4257, 2017.

[51] R. C. Bradt. "Applying fractography and fracture mechanics to the energy and mass of crack growth for glass in the mirror region." *J. Eur. Ceram. Soc.*, **34**, 3255–3262, 2014.

[52] A. Polian, D. Vo-Thanh, and P. Richet. "Elastic properties of a-SiO2 up to 2300 K from Brillouin scattering measurements." *EPL (Europhysics Letters)*, **57**, 375, 2002.

[53] Q. Zhao, M. Guerette, G. Scannell, and L. Huang. "In-situ high temperature Raman and Brillouin light scattering studies of sodium silicate glasses." *J. Non-Cryst. Solids*, **358**, 3418–3426, 2012.

[54] A. J. Rosakis. "Intersonic shear cracks and fault ruptures." *Adv. Phys.*, **51**, 1189–1257, 2002.

[55] L. Bergmann. *Ultrasonics and their scientific and technical applications*. Wiley, 1948.

[56] J. D. Achenbach. *Wave Propagation in elastic solids*. Elsevier, 1973.

[57] M. A. Meyers. *Dynamic behavior of materials*. John Wiley & Sons, 1994.

[58] G. D. Quinn. "Fractography of ceramics and glasses: A NIST recommended practice guide." *Special Publication 960-16*. Washington, DC: National Institute of Standards and Technology, 2007.

[59] H. Schardin. "Velocity effects in fracture." In *ICF0, Swampscott-MA (USA)*, 1959.

[60] S. R. Anthony, J. P. Chubb, and J. Congleton. "The crack-branching velocity." *Philos. Mag.*, **22**, 1201–1216, 1970.

[61] J. Horbach, W. Kob, and K. Binder. "Structural and dynamical properties of sodium silicate melts: an investigation by molecular dynamics computer simulation." *Chem. Geol.*, **174**, 87–101, 2001.

[62] https://zenodo.org/record/6499675#.YmuyT1xBw_0.

[63] M. Bauchy, B. Wang, M. Wang, Y. Yu, M. J. Abdolhosseini Qomi, M. M. Smedskjaer, C. Bichara, F.-J. Ulm, and R. Pellenq. "Fracture toughness anomalies: Viewpoint of topological constraint theory." *Acta Mater.*, **121**, 234–239, 2016.

[64] F. Shimizu, S. Ogata, and J. Li. "Theory of shear banding in metallic glasses and molecular dynamics calculations." *Mater. Trans.*, **48**, 2923–2927, 2007.

[65] M. L. Falk and J. S. Langer. "Dynamics of viscoplastic deformation in amorphous solids."
*Phys. Rev. E*, **57**, 7192, 1998.

[66] A. P. Thompson, S. J. Plimpton, and W. Mattson. "General formulation of pressure and stress tensor for arbitrary many-body interaction potentials under periodic boundary conditions." *J. Chem. Phys.*, **131**, 154107, 2009.

[67] G. Molnár, P. Ganster, J. Török, and A. Tanguy. "Sodium effect on static mechanical behavior of MD-modeled sodium silicate glasses." *J. Non-Cryst. Solids*, **440**, 12–25, 2016.

[68] C. Spearman. "The proof and measurement of association between two things." *The Am. J. Psychol.*, **15**, 30, 1901.

[69] S. N. Taraskin and S. R. Elliott. "Connection between the true vibrational density of states and that derived from inelastic neutron scattering." *Phys. Rev. B*, **55**, 117–123, 1997.