# Fracture of silicate glasses: Micro-cavities and correlations between
atomic-level properties

Zhen Zhang  Laboratoire Charles Coulomb (L2C), University of Montpellier and
CNRS, F-34095 Montpellier, France  Simona Ispas  Laboratoire Charles Coulomb
(L2C), University of Montpellier and CNRS, F-34095 Montpellier, France  Walter
Kob<EMAIL_ADDRESS>Laboratoire Charles Coulomb (L2C), University
of Montpellier and CNRS, F-34095 Montpellier, France

(August 27, 2024)

###### 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.

##  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 binder_glassy_2011 ; shelby_introduction_2005 ;
varshneya2013fundamentals . 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
wondraczek_towards_2011 ; wondraczek2022review . 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 mallinder1964elastic ; krause_deviations_1979 ; gupta_intrinsic_2005
; guerette_nonlinear_2016 ; zhang2022stiffness . 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
tomozawa1996fracture ; ciccotti_stress-corrosion_2009 ;
wiederhorn_griffith_2013 . 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
rountree2002atomistic ; celarie_glass_2003 ; guin_fracture_2004 ;
bonamy_nanoscale_2006 ; muralidharan_molecular_2005 ; rountree_unified_2007 ;
pedone_molecular_2008 ; pedone_dynamics_2015 ; wang_intrinsic_2015 ;
shen2021observation , 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 wiederhorn_use_2011 . 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
ciccotti_situ_2018 . 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 celarie_glass_2003 ;
guin_fracture_2004 ; bonamy_nanoscale_2006 ; waurischk2021vacuum . 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 wang_intrinsic_2015 ; wang_nanoductility_2016 ;
pedone_molecular_2008 ; pedone_dynamics_2015 . 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
han_measuring_2010 ; pallares2012fractoluminescence , 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 rimsza_crack_2018 ; rimsza2018chemical ;
chowdhury_effects_2019 , as well as the heterogeneities in local structure and
mechanical properties and the correlations between them wang_intrinsic_2015 ;
wang_nanoductility_2016 ; hao_atomistic_2019 ; frankberg2019highly ;
du2021predicting . 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 SiO2 and Na2O-$x$SiO2 (denoted in the following
by NS$x$), with $x=3,4,5,7,10,$ and 20, corresponding to the variation of Na2O
from 0 to 25 mole %. The interactions between the atoms were given by the
pairwise effective potential proposed by Sundararaman et al.
sundararaman_new_2018 ; sundararaman_new_2019 . 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
wolf_exact_1999 . More details regarding the development of the potential and
its ability to reproduce experimental data and ab initio calculations are
given in Refs. sundararaman_new_2018 ; sundararaman_new_2019 . 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 zhang_potential_2020 ; zhang_surf-vib_2020 ; zhang2021roughness , 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 bansal_handbook_1986 , corresponding to
boxes with typical dimensions of 20 nm$\times$30 nm$\times$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_{1}$ (still in the
liquid state) for another 160 ps. The temperature $T_{1}$ ranges from 3000 K
for SiO2 to 2000 K for NS3 (25 mole % Na2O), see Ref. zhang_surf-vib_2020 ;
zhang_thesis_2020  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_{1}$ 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 $\gamma_{1}=0.125$ K/ps was used to quench the
samples from $T_{1}$ to a temperature $T_{2}$ and a faster cooling rate
$\gamma_{2}=0.375$ K/ps to cool them from $T_{2}$ to 300 K. Finally, the
samples were annealed at 300 K for 800 ps. The temperature $T_{2}$ 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
$\geq T_{2}\geq 800$ K, see Refs. zhang_surf-vib_2020 ; zhang_thesis_2020  for
details. At $T_{2}$, 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-1, 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. zhang_thesis_2020 .

Temperature and pressure were controlled using a Nosé-Hoover thermostat and
barostat nose_unified_1984 ; hoover_canonical_1985 ; hoover_constant-
pressure_1986 . All simulations were carried out using the Large-scale
Atomic/Molecular Massively Parallel Simulator software (LAMMPS)
plimpton_fast_1995  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 zhang_thesis_2020 . 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

###  III.1 Stress-strain behavior

Figure 1: (a) Stress-strain behavior of the NS$x$ 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).

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. zhang2022stiffness  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
zhang_potential_2020  or in experiments lower_inert_2004 . 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) 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 Na2O. 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 NaO2 concentration is shown in Fig. 1b and one sees that $v_{c}$
decreases from $\approx$1800 m/s for silica to $\approx$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 103 m/s which is what one typically finds
in experiments probing the dynamic fracture of silicate glasses
quinn2019terminal ; quinn_fracture_2017 ; bradt_applying_2014 ;
pallares2012fractoluminescence . 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 zhang2022stiffness  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 polian2002elastic  and the ones for sodium silicate
glasses are a few percents lower zhao_-situ_2012 ; zhang_thesis_2020 . 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 ciccotti_situ_2018 . 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
zhang_thesis_2020 .

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
rosakis_intersonic_2002 . To get an idea of how our measured crack velocities
compare with $v_{R}$, we estimate the latter from the expression
bergmann_ultrasonics_1948 ; achenbach_wave_1973 ; meyers_dynamic_1994 ;
quinn2019terminal

$v_{R}=\frac{0.87+1.12\nu}{1+\nu}\sqrt{\frac{G}{\rho}},$ (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 $\varepsilon$ 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 Na2O 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% Na2O (NS10), i.e., the composition at which a significant change
in the elastic properties of the sodium silicate glasses is observed
zhang2022stiffness . 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.,
Na2O%$\leavevmode\nobreak\ <$ 13%, the estimated $v_{c}/v_{R}$ is compatible
with the value of $0.5\sim 0.6$ that one usually finds in experimental
measurements quinn_fractography_2007 , 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 zhang2021roughness .
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
schardin_velocity_1959 ; anthony_crack-branching_1970 .

###  III.2 Cavitation

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$\times$30 nm$\times$50 nm. The green dashed line represents $N_{\rm
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 $\rho_{\rm void}$, 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 $\rho_{N}$ of the glass increases with the addition of Na2O
($\rho_{N}=66.82$ nm-3 for silica and $73.83$ nm-3 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 horbach_structural_2001
.

From the graph one can conclude that, independent of the Na concentration,
voids with a size of $R\approx 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
$\approx 3.5$ Å muralidharan_molecular_2005 , i.e., significantly less than
the threshold we find here. However, in that work the samples were cubic boxes
with a size of $\approx$ 3.7 nm, and our data for silica shows that at $R=3.5$
Å one can expect to find about $10^{-2}$ voids per nm3, i.e., 0.5 voids in a
cube of size 3.7 nm, in good agreement with the estimate of Ref.
muralidharan_molecular_2005 .

.

Figure 3: (a) and (b): Snapshots showing the empty space during the fracture
of silica, (a), and NS3, (b), glasses (sample size $20\times 30\times 50$
nm3). 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. movies . (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 nm3) and number of the voids.

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 Ref. movies ). 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_{\rm 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_{\rm void}$ (right scale), and one sees that this quantity
tracks $N_{\rm void}$ perfectly well. The fact, see inset, that the ratio
between $V_{\rm void}$ and $N_{\rm void}$ is around 0.286 nm3 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 guin_fracture_2004
and a recent simulation study using reactive force field
chowdhury_effects_2019 .

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_{\rm
void}/N_{\rm 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_{\rm void}/N_{\rm
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 zhang2021roughness .

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 wang_intrinsic_2015 ; pedone_dynamics_2015 . 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 bauchy_fracture_2016 . (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
zhang_thesis_2020 .)

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 $\approx$3.5 nm, while in the
$y-$direction (parallel to the stress direction) the length is only
$\approx$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
bonamy_nanoscale_2006 . 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 guin_fracture_2004 .

Figure 4: Enlarged view of the largest void found at $\varepsilon=0.196$ for
the fracture of a NS3 glass.

###  III.3 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_{\rm Na}=N_{\rm Na}/N_{\rm tot}$, where
$N_{\rm Na}$ and $N_{\rm tot}$ are, respectively, the number of Na atoms and
the total number of atoms within a given probe volume (defined below). We
recall that in NS3 the stoichiometric value of $f_{\rm Na}$ is equal to 0.167;
(ii) The local network connectivity as represented by the coordination number
of oxygen, $Z_{\rm 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}(\varepsilon)=\theta_{i}(\varepsilon)-\theta_{i}(0)$, where
$\theta_{i}(0)$ is the SiOSi angle of O atom $i$ in the initial configuration
(0% strain), and $\theta_{i}(\varepsilon)$ is the SiOSi angle of the same O
atom in the sample at strain $\varepsilon$. (Note that it is possible that a
SiOSi bridge present at 0% strain does not exist anymore at strain
$\varepsilon$, since one of its Si-O bonds has broken. In this case, we simply
assign $\Delta\theta_{i}(\varepsilon)=0$.); (iv) The atomic shear strain as
first introduced by Shimizu et al. shimizu_theory_2007  for measuring the
local plastic deformation, i.e., a quantity that carries basically the same
information as the non-affine squared displacement $D_{2}^{\rm min}$
introduced by Falk and Langer falk_dynamics_1998 ; (v) The atomic stress
defined via the per-atom stress tensor thompson2009general . (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 molnar_sodium_2016 ;
wang_nanoductility_2016 . 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/cm3 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. zhang_thesis_2020 ).

Figure 5 shows the maps of the various local properties at two representative
strains during the fracture of the NS3 glass, i.e., $\varepsilon=0.172$ (the
fracture has just started) and $\varepsilon=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. movies .) 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_{\rm Na}$ and the local coordination number of O, $Z_{\rm
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.

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\times 30\times 50$ nm3. 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.

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. zhang2022stiffness
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^{\rm Mises}$ depends strongly on the global strain. For
small $\varepsilon$ the distribution has a narrow peak at small values of
$\eta^{\rm 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\geq
0.16$, one finds that the distribution has a noticeable tail at large
$\eta^{\rm 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^{\rm 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^{\rm Mises}$ is larger than 0.05,
demonstrating that each part of the sample has significant non-affine strain.

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^{\rm
Mises}$. (c) Atomic tensile stress. (d) Kinetic temperature $T$. Note that for
this analysis the layers, $\approx 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.

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 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
zhang_potential_2020 ; zhang2022stiffness . 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 $\varepsilon$ 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 $\varepsilon=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 $\varepsilon=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.
$\varepsilon=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 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.

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 pallares2012fractoluminescence . (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 zhang_thesis_2020 .)

To quantify the correlation between various local properties we use the
Spearman rank correlation coefficient, $r_{s}(X,Y)$, which is a measure for
the statistical dependence between two variables $X$ and $Y$
spearman_proof_1901 . By construction $r_{s}$ is in the range [-1, 1] and
$r_{s}>0$ means positive correlation, i.e., $Y$ tends to increase when $X$
increases, and $r_{s}<0$ indicates an anti-correlation. The magnitude of
$r_{s}$ gives the degree of correlation, with $r_{s}=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_{s}({\rm Si,Na})\approx-0.94$. A strong anti-correlation
is also found between O and Na with $r_{s}({\rm O,Na})\approx-0.90$, in
agreement with previous studies on soda-silicate glasses molnar_sodium_2016 .
From panel (a) one also recognizes that $f_{\rm Na}$ and $Z_{\rm O}$ 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 properties. Overall, we note that $f_{\rm Na}$ 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
connectivity 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_{\rm Na}$ and $\eta^{\rm 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_{\rm 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 observes 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
$\varepsilon$ beyond the failure strain. The strain dependence of this
correlation can be understood by considering the coexistence of multiple
deformation modes during tension: As it was shown in Ref.
taraskin_connection_1997  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.
zhang2022stiffness , 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 zhang2022stiffness ;
wang_intrinsic_2015 ; molnar_sodium_2016 ; wang_nanoductility_2016 , 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 $\approx 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 University.

## References

  * (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.  “Atomistic 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#.YmuyTlxBw_O.
  * (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.