Statistical analysis of fracture of silica glass using molecular dynamics

Disordered solids exhibit intermittent avalanches when slowly driven by an external load. These avalanches are associated with plastic rearrangements of the atoms at the nanoscale that manifest as stress and energy drops in the loading curve. The complexity arising from their interactions through long‐range elastic fields and the disorder makes statistical approaches suitable for studying their behavior by considering yielding and fracture as a phase transition. To investigate the avalanche statistics, we perform nanoscale simulations using the athermal quasistatic deformation protocol to induce fracture to silica glass samples, which were simulated using the melting‐quenching technique. We identify and measure the avalanches and voids and investigate their evolution in the process of fracture. Our results confirm that the avalanches exhibit scale‐free power law statistics, indicative of a critical phenomenon.


INTRODUCTION
Silica glass is among the most used materials in a wide range of industries, such as electronics, optics, telecommunications, and manufacturing [1].It possesses a wide range of properties that make it highly versatile, such as transparency, high thermal resistance, chemical inertness, and excellent electrical insulation capabilities.However, its brittleness and susceptibility to fracture greatly limit its uses, as it tends to fail in a catastrophic manner.Understanding the fracture process in disordered materials such as silica glass is challenging due to its inherent randomness [2,3].However, statistical mechanics provides a valuable framework to deal with the complexity by drawing an analogy between fracture and phase transition [4].
Studies have shown that, when a disordered solid is driven by an external field, it undergoes a phase transition from the intact to the broken phase, where the quenched disorder plays the role of temperature [4].During the deformation course, the potential energy landscape is continuously deformed while the system remains in an energy minimum.However, with increasing stress, the minimum where the system resides vanishes, and the system is driven to a lower energy minimum [5].This is known as a plastic rearrangement event, where the atoms are rearranged to a new configuration with lower energy.Moreover, the occurrence of a single event can trigger a cascade of subsequent events, akin to snow avalanches.Consequently, the sudden burst of plastic events is commonly referred to as an avalanche.In the athermal case, the avalanches are expected to display spinodal criticality due to the long range elastic fields [6], and, therefore, follow power law statistics [4].Moreover, they are believed to originate from soft spots in the system [7] and could be classified as angle changing or bond breaking [8].The latter leads to the creation of voids and could be related to irreversible processes in bulk silica, resembling macroscopic damage [9].Therefore, exploring the voids and their growth in the system, is crucial to understand the fracture behavior in silica glass.
In this paper, we use molecular mechanics simulations to study the fracture of silica glass at the nano-scale under tensile stress.We employ the athermal quasi-static protocol (AQS) on pre-notched systems for four different sizes.A novel method is introduced for recording and measuring stress and energy avalanches in the system.The statistical analysis reveals that these avalanches follow power law statistics, indicating critical behavior.The critical exponents are determined and the effect of the initial notch is investigated.Additionally, we capture the crack growth process and conduct a statistical analysis of its corresponding avalanches.

Samples preparation
We employ classical molecular dynamics in LAMMPS to simulate our systems using quenching-melting techniques.We simulate four different sample sizes to study the scaling behavior., where   is the inter-atomic distance,   refers to the effective charge while   and   are parameters taken from Jakse et al. [11].The three terms in the equation represent the Coulombic, repulsive, and Van der Waals interactions respectively.The cut-off distance was set to 10.17 Å.
Using NVT and periodic boundary conditions each system is heated to 8000 K far above the silica glass melting temperature of 1470 K. Then the system is kept at constant temperature for 2 ns to ensure that it has no memory of the initial configuration.Subsequently, the samples are generated by quenching at equidistant time spans of 250 ps at cooling rate of 10 K/ps until reaching a temperature close to zero  ≈ 0.001 K.After that, the sample is equilibrated for additional 10 ps to minimize the pressure fluctuations, followed by setting the velocities to zero and minimization of the potential energy.At this point, the samples are validated by comparing the measured densities and the radial distribution functions (RDFs) to experimental data [12].A notch is created in every sample by removing atoms from a cylindrical volume with an elliptical cross-section of 2.5 × 5 Å in the -plane which spans across the entire sample in the -direction.It is worth noting that the notch was kept at minimum volume, while still constituting a critical initial crack that leads to material failure.The same process is repeated for four different system sizes with 48,000, 24,624, 12,150, and 6100 atoms.Fifty samples are generated for each system size.Each sample is loaded in the -direction using the AQS protocol: at each step, the system is deformed by imposing a strain of Δ = 1 × 10 −4 in the -direction, followed by minimizing the potential energy of the system using the conjugate gradient method until reaching a strain of  = 0.25.
Finally, for analyzing avalanche statistics, we employ two measures: stress and energy dissipated during avalanche events.However, simply taking the difference in stress or energy before and after the event is not sufficient, as it doesn't account for the elastic energy introduced into the system by straining the sample.To compute the associated part, one can utilize the elastic constant of the material, denoted as Δ.However, in our case, the elastic constant continuously changes throughout the deformation process.Therefore, we employ a new method to record the avalanches.First, the avalanches are identified as local maxima in the stress difference curve Δ  =  +1 −   .Then, we construct a running average curve of Δ  during times without events and interpolate it at the time of the events.Finally, the size of each avalanche is computed as the difference between Δ  and the running average curve.A sketch of this approach is shown in Figure 1.
Avalanches in the course of mechanical deformation; the black line is the stress difference between two consecutive steps Δ  =  +1 −   plotted versus the strain.The avalanche sizes and frequencies increase exponentially as one approaches failure; (B) A stress drop caused by an avalanche; The red arrow represents the difference in the stress after and before the event.The blue arrow represents the correction to the event that should be added according to the stiffness of the material at the time of the event.The green arrows are what should be the total size of the avalanche; (C) An illustrative plot showing the method used in this work for recording and measuring the avalanches; The red arrow shows the avalanche measurements if one considers avalanches for only Δ  > 0 while green arrows show the avalanche measurements if one records avalanches as fluctuations of Δ  .

Void detection algorithm
The microstructure of vitreous silica is highly complex.Thus, conventional methods for calculating the crack depth do not provide the sufficient accuracy we need in order to evaluate the distribution of the discrete jumps occurring during the crack propagation.Therefore, we have chosen an algorithm based on the research work done by Sastry et al. by using a library provided by Voloshin et al. [13,14].The algorithm employed in this study operates by surrounding the atoms by exclusion spheres then breaking down the empty space between them into individual cavities.It then proceeds to calculate the volumes and surface areas of these cavities in three-dimensional packing of spheres.In our case, the system consists of SiO4 tetrahedra [15] with approximately 1.69 Å bond length between silica and O atoms.An exclusion sphere radius has to be accordingly chosen, so the space between the bonded atoms is not considered free volume.Thus, a radius of 2 Å has been chosen.

Statistical models
For this study, we utilized the scaling relations provided by two simple statistical models commonly used to investigate fracture phenomena: the fiber bundle model (FBM) and the mean-field depinning model.In the FBM, the susceptibility, defined as , where  is the order parameter, is expected to scale asymptotically with the average avalanche size ⟨⟩.Moreover, for a wide range of distributions of the strength of individual fibers, the avalanches are expected to follow the scaling relation (, ) ∼  −  (  −)  [16], where  represents the avalanche size,  denotes the applied force,   represents the critical force leading to bundle failure, as well as  = 3∕2 and  = 1 are critical exponents.However, when recording the avalanches throughout the entire loading regime until the critical force, the probability density function (PDF) takes the following form ()  ∼  − ′ (, ), where  ′ = 5∕2 and  denotes the tail correction function.Additionally, as the bundle size increases, the system becomes more brittle and more susceptible to experiencing catastrophic events.If the simulations are repeated sufficiently, these catastrophic events are expected to follow a Gaussian distribution [17].
Applying the mean-field depinning model, the avalanche PDF predicts the same form as the FBM but with a different exponent value of  = 2 if it is solved analytically for brittle systems, [18].Consequently, when recording the avalanches over the entire loading regime, the avalanche exponent becomes  ′ = 2.However, these predictions are expected to deviate in the vicinity of the critical force, as system-spanning events start to emerge [18].

RESULTS
The stress-and energy-strain curves exhibit elastic branches interspersed with sudden stress drops for individual samples.Moreover, when examining the averaged curves across different samples, it becomes evident that the brittleness of the system increases with larger system sizes as can be seen in Figures 2A,B.Additionally, the critical strain, indicating sample failure, decreases as the system size increases, following a power law relationship:   ∼  −0.183 .Furthermore, the volume growth of the initial crack shows similar trends as depicted in Figure 2C.Each drop in energy and stress corresponds to a plastic rearrangement event, where one or more clusters of atoms undergo non-affine displacements.Notably, our simulations reveal that these events predominantly exhibit a dipolar shape, in contrast to sheared glasses where quadrupolar-shaped events are more prevalent [5].The shape of these events influences the waiting time between avalanches, underlining the importance of their characterization [19].
Moreover, we observe that, in the elastic regime, these events are distributed throughout the system, indicating the competition between inherent disorder and the localizing effect of the initial notch, as illustrated in Figure 3A, where the accumulated non-affine displacement field during the simulation is plotted for different samples.However, beyond a certain threshold, an increasing number of events concentrate around the notch, leading to its growth and eventual nucleation, ultimately resulting in system failure.This behavior is depicted in Figure 3B,C   an analysis of the average decomposition in the vicinity of the upper tip of the initial crack (for a cutoff radius  = 8 Å) using the RDF for Si-O atom pairs, as illustrated in Figure 4. Our observations reveal that, as the material approaches failure, there is a slight broadening of the distribution near the first peak, indicative of bond stretching.Furthermore, the second peak in the distribution, corresponding to the second coordination shell, exhibits a flattening effect, suggesting an increased level of mid-range disorder.Subsequently, following crack propagation, the atoms partially revert towards their original state, albeit with some remaining deformations.

Avalanches
A total number of 23,025 stress and energy avalanches were analyzed in our study.The average waiting time between the events exhibits a power law relationship with the system size, ⟨Δ⟩ ∼  −0.42 .The partial distribution function (PDF) of energy drops is shown in the inset of Figure 5A.It is evident that the observed avalanches in our systems follow scalefree power laws.However, an exponential cutoff function is observed, which depends on the system size.By assuming a function of the form (ΔΠ) =   (ΔΠ∕ ∕ ), we computed the critical exponents  ′ ΔΠ = 1.75 and  ΔΠ = 2.57 through finite-size scaling (FSS) [20].These exponents were then used to scale the data in the main plot of Figure 5A, where it can be seen that the data collapses onto a single curve.Similarly, the stress drop avalanches for different system sizes follow power law statistics, as depicted in the inset of Figure 5B.We computed the exponents to be  ′ Δ = 1.52 and  Δ = 2.52, which were used to scale the data in the main plot, resulting in a satisfactory data collapse.Furthermore, it is known that F I G U R E 5 Avalanche statistics for the energy drops in (A) and the stress drops in (B).Where the avalanches have been recorded until the critical strain   of each system size.The insets show the unscaled PDFs while the main figures have been scaled with the exponents  ′ ΔΠ = 1.75 and  ΔΠ = 2.57,  ′ Δ = 1.52 and  Δ = 2.52 which were obtained by FSS.Both figures show satisfactory collapse of the data indicating the validity of the computed exponents.PDFs, probability density functions; FSS, finite-size scaling.

(A) (B)
F I G U R E 6 (A) Crack growth avalanches the has been recorded during the whole loading protocol until the failure of system.The dashed lines represent the exponents that has been computed through direct fitting of the lines for the biggest system.(B) The average number of cavities across various samples is plotted against the applied strain.The inset shows that the number of cavities for different system sizes collapses into a single curve when normalized by the total number of atoms in the system.
the fractal dimension is related to the avalanche exponent as   ∕ = ∕, where  = 3 is the spatial dimension of the system.Therefore, we computed the fractal dimension from the energy and stress avalanches to be approximately 4.4 and 4.97, respectively.The obtained results of the avalanches exhibit a significant deviation from the mean-field exponents so far, where one would expect an exponent of  ′ = 2.5 when integrating over the entire strain regime, as in the FBM, and an exponent of  ′ = 2 for the mean-field depinning model.Moreover, the computed values of the fractal dimension are much higher than the values typically reported in fracture models, which usually lie between 1.0 and 1.25 [17].However, it is important to note that our simulations have been conducted on pre-notched systems.Therefore, we performed additional tensile simulations on un-notched specimens with 48,000 atoms to study the effect of the initial notch.We computed the stress and energy exponents for the un-notched specimen to be approximately  ′ ΔΠ ≈ 2 and  ′ Δ ≈ 2, indicating that the initial notch does not only deviate from the avalanche exponent  ′ from the mean-field value but also creates a discrepancy between the energy and stress avalanches.
Regarding crack growth, we recorded all events that resulted in a volume increase larger than a threshold of 1 × 10 −3 nm 3 .A total number of 11,253 events were recorded and analyzed.The PDF is plotted in Figure 6.As shown in this plot, we find that events smaller than 0.3 nm 3 follow a scale-free power law for all system sizes while larger events exhibit different behavior for each system size.When considering only the largest system, the PDF resembles a double power law with a cutoff for the large events.By direct fitting we computed the exponents  Δ1 = 0.9 and  Δ2 = 2 for the small regime and big regime accordingly.Moreover, it is warranted to investigate cavities during the deformation process, as they represent regions with clusters of bond breakage.In the context of first-order transition, we expect the number of clusters to scale as   =   (∕) [4].Figure 6B illustrates that the scaling relationship for the number of clusters is obeyed.Additionally, it demonstrates that the average number of clusters peaks and then rapidly decreases before the material fails, indicating the coalescence of cavities with the main crack, followed by a phase of pure crack propagation.

CONCLUSION
In conclusion, our study sheds light on the statistics of avalanches in fracture of silica on the nanoscale.We have observed that these avalanches exhibit power-law scaling, indicating the presence of critical phenomenon.We have computed the critical exponents for stress and energy avalanches,  ′ Δ = 1.52 and  ′ Δ = 1.75.Furthermore, the crack growth avalanches are recorded showing a peculiar scaling behavior, where only the small avalanches follow a power law.
Our study demonstrates that our systems exhibit criticality but does not conform to mean field calculations.However, we show that this is (at least in part) due to the presence of an initial notch that favors a localizing effect of the plastic events around it.This effect leads to deviation from the mean-field exponents, especially, in the cutoff exponent  and creates a discrepancy between the stress and the energy avalanche statistics.However, it is shown that the avalanche statistics can be qualitatively captured by the fiber bundle and the mean-field depinning model and that the latter yields reasonable predictions for the exponents of the unnotched systems.However, it is important to note that the power-law scaling in crack volume growth avalanches is confined to a limited regime, prompting further research to fully understand the extent of this phenomenon.Finally, the number of damage clusters was analyzed and shown to conform to first order transition.However, a more comprehensive analysis is warranted in future research, that includes average cluster size near the propagating crack tip in addition to studying unnotched systems for confirming the universality of the system.
Looking ahead, our study opens up interesting possibilities for further research in this field such as exploring the relationship between crack growth avalanche statistics and crack surface morphology.This could provide valuable insights into the underlying mechanisms.Furthermore, our study prompts further inquiry into the type of critical phenomenon observed during crack growth.Understanding the critical phenomenon at play has the potential to enhance our understanding of crack dynamics and bridge the gap in our knowledge between the nano and the macroscale.

A C K N O W L E D G M E N T S
Open access funding enabled and organized by Projekt DEAL.

F I G U R E 2
Loading curves of the different sample sizes.The thin lines are from the individual samples of each system and the thick lines represent the average across the different samples of one system size.(A) Stress-strain curves.(B) The potential energy-strain curves where the energy of each system size has been shifted to zero for facilitating a comparison between the different system sizes.(C) Crack growth during loading, where  is the volume of the initial notch.
. In addition, we conducted F I G U R E 3 Cumulative non-affine displacement field mapped to the atoms in the reference configuration in the -plane, at strain  = 0.107 in (A),  = 0.1318 in (B) and  = 0.1492 in (C).The plot is repeated in both directions to show the periodic boundary condition, where the black box represent the original system.

F I G U R E 4
The radial distribution function for Si-O atom pairs within a cutoff radius of 8 Å from the upper tip of the initial crack in the largest system.
[11]the largest system, a rectangular parallelepiped simulation box with dimensions of 135 × 85 × 54 Å is used, which is then divided into 16,000 cells.A total of 48,000 atoms (16,000 Si and 32,000 O) are randomly distributed in the box, with each cell containing one Si atom and two O atoms.The dimensions and number of atoms are chosen to approximate the experimental density of silica glass, which is  ≈ 2.2 [g/cm 3 ].The other system sizes are created in the same manner, maintaining the box proportions and the same density.We used a potential based on the Born-Mayer-Huggins potential, first proposed by Matsui[10], and re-parameterized by Jakse et al.[11].The potential between two atoms of type  (Si)  (O) is written as   (  ) =