Parametrization in Models of Subcritical Glass Fracture: Activation Offset and Concerted Activation

There are two established but fundamentally different empirical approaches to parametrize the rate of subcritical fracture in brittle materials. While both are relying on a thermally activated reaction of bond rupture, the difference lies in the way as to how the externally applied stresses affect the local energy landscape. In the consideration of inorganic glasses, the strain energy is typically taken as an off-set on the activation barrier. As an alternative interpretation, the system’s volumetric strain-energy is added to its thermal energy. Such an interpretation is consistent with the democratic fiber bundle model. Here, we test this approach of concerted activation against macroscopic data of bond cleavage activation energy, and also against ab initio quantum chemical simulation of the energy barrier for cracking in silica. The fact that both models are able to reproduce experimental observation to a remarkable degree highlights the importance of a holistic consideration towards non-empirical understanding.


INTRODUCTION
The chemistry of glass fracture has been under scientific debate for a long time (Fairbarn and Tate, 1859;Griffith, 1921). However, a quantitative description of the underlying mechanisms of crack propagation, critical to the further development of crack-resistant glasses, is presently unavailable (Wondraczek et al., 2011;Freiman, 2012). This is largely due to the fundamental nature of cracking as it occurs in the broad group of brittle oxides, and especially in glassy silicates: crack propagation reflects a complex and, as of yet, not fully understood interplay of chemical and physical interactions (Ciccotti, 2009). In particular, below a critical value of applied stresses (the subcritical regime), chemical reactions and the effects of mechanical loading are so intertwined that the separation of the governing parameters becomes an extremely challenging task. The subcritical crack growth and its dependence on environmental conditions have profound consequences for the glass industry and the suitability of glass as a structural material. Here, parametric knowledge is needed not only for the focused development of new material compositions and toughening procedures but also for enabling better time-to-failure estimations (Freiman et al., 2009). Beyond materials development, this need encompasses other related areas such as the understanding of geophysical fracture processes (Atkinson, 1982;Brantut et al., 2013;Zhang and Zhao, 2014).
Generally speaking, the kinetics of crack propagation below Mode I critical stress concentration K Ic can be divided into three regimes: in region I, even though the applied stresses are much smaller than the critical value, the crack propagates through a directional corrosion reaction on stressed bonds. For the case of molecular water as the reaction partner, a theoretical framework has been FIGURE 1 | Reactive crack propagation in silica. (A) Atomistic view of a crack tip in silica glass; (B) subcritical crack growth mechanism (Michalske and Freiman, 1982). In the vicinity of the crack tip, one or more water molecules adsorb at a stressed Si-O-Si linkage (I), the Si-O bond is cleaved via a transition structure (II) and finally converted to surface silanol groups (III).
developed (Michalske and Freiman, 1982;Bunker, 1984, 1987) (Figure 1). In this region, the crack propagation velocity depends strongly on the applied stresses, temperature, and environment. In region II, the crack growth velocity becomes increasingly limited by the transport of the reactive molecules to the crack tip, resulting in a constant crack velocity with increasing stress intensity. Finally, in region III, the applied stresses are high enough to break the stressed atomic bonds without environmental help, so the crack velocity is again strongly dependent on the applied stresses and temperature, but not on the environment.

ACTIVATION OFFSET
Generally, the rate of crack growth v is modeled as ν = λ·k (Michalske and Bunker, 1984;Schoeck, 1990), where λ is a geometrical parameter depending on the matrix through which the crack is growing, and k is the rate of success for the crack to advance by λ units. The reaction rate is then framed within the transition state theory, and written in the form of an Arrhenian expression, k = k 0 exp(−ΔF/RT), where k 0 is the attempt frequency based on the oscillations of the system around an equilibrium position and ΔF is the molar energy difference between the initial and final states, in this case representing the energy barrier to break the bond. While working on the delayed fracture of glasses, Wiederhorn and coworkers (Wiederhorn, 1967;Wiederhorn et al., 1974a,b) pioneered a semiempirical equation to describe the phenomenon of subcritical crack growth in the presence of water: . (1) Here, v 0 is the pre-exponential term, a(H 2 O) is the water activity, n is the number of water molecules which interact with the Si-O-Si bond during its cleavage, E* is the activation energy, K I is the stress concentration factor at the crack tip, and b is a function of the activation volume and the crack tip radius (Wiederhorn et al., 1974a,b). The mathematical formulation for the activation energy in Eq. 1 is analogous to the enthalpy of gasses under constant pressure, only the pressure term is substituted by the applied stresses (Wiederhorn, 1967;Wiederhorn and Bolz, 1970). This expression has found widespread use in the literature, influencing a great part of the developments on the understanding of subcritical crack growth since its inception.
For the moment, we note that the primary feature of Eq. 1 is the reduction of E* by the strain energy U strain = bK I . This can be understood as an offset contribution in the potential well of the considered bond, increasing the probability for breakage (case I in Figure 2).

CONCERTED ACTIVATION
The derivation of Eq. 1 assumes that the applied mechanical field offsets the potential energy at the considered bond, thus reducing the activation energy of the breakage event and increasing breakage probability. However, this appears as only one of two extremes: the second is an increase in the overall excitation of the system (case II in Figure 2) without directly affecting the activation barrier. This interpretation is a natural result of the democratic fiber bundle model (DFBM), in which the material is represented as a series of springs connecting two parallel bars and oriented perpendicular to the direction of crack propagation. These springs are considered to be ideally elastic, ideally brittle with random failure threshold (modeling inherent structural disorder), and the load on each spring is also subjected to thermal noise (Kun et al., 2000;Roux, 2000;Ciliberto et al., 2001;Scorretti et al., 2001).
Here, we consider first the energetic equilibrium for cracks growing at rates below the limit of dynamic stability (Fineberg and Marder, 1999), FIGURE 2 | Schematic representation of how temperature and applied stresses interact with the potential energy well of the chemical bond at the crack tip. <U> is the average potential energy of the bond as a function of distance, U thermal is the thermal energy, U strain is the strain energy, and E*, E* ,1 , and E* ,2 are the activation energies for bond cleavage for non-strained bonds and strained bonds, case I and II, respectively.
where G is the strain release rate, μ is Poisson's ratio, E is the elastic modulus, A I is a material and crack speed dependent variable, and γ is the surface energy. Equation 2 can be reasonably approximated by Freund (1990), with c R being the Rayleigh wave speed. However, Eq. 4 is only valid for K I > K Ic , so to extend the range for which the equation is applicable, it has to be rewritten. Here, we consider the term as equivalent to the first two terms of a Taylor series expansion, leading to this original equation: .
Equation 5 fulfills the requirement of convergence to Eq. 4 for K I > K Ic , and additionally it produces a non-zero crack growth velocity at K I = K Ic and v = 0 for K I = 0. From the transition state theory, we can write that where the energy difference between the initial and final states ΔF is approximated by the ratio between the energy needed to create two surfaces and the number of chemical bonds that make up this surface, so N A is Avogadro's constant and surface bond density (ρ A ). This approach is consistent with the approximation of ideally brittle fracture and the lattice trapping understanding of fracture (Hsieh and Thomson, 1973;Lawn et al., 1980;Michalske and Bunker, 1984;Schultz et al., 1994;Zhu et al., 2004). Putting Eqs 5 and 6 together yields: or alternatively, where v ′ 0 is the pre-exponential term, E* = 2γN A /ρ A and C = (1 − μ 2 )N A /Eρ A . As noted above, Eq. 8 shares the Arrhenian scaling with the approach of activation offset, and by including both thermal as well as structural noise into its variance, there is a clear parallel with the average time-to-failure predictions derived through DFBM (Kun et al., 2000;Roux, 2000;Ciliberto et al., 2001;Scorretti et al., 2001). The main difference to Eq. 1 is that instead of decreasing the magnitude of the energy barrier to bond failure, the strain energy U strain = CK I 2 acts concertedly with the thermal energy RT, increasing the excitation level of the bond in question. This collaborative action can be intuitively thought as stemming from the fact that elastic deformation and thermal expansion have its origin on how the acting stress and temperature interact with the interatomic bond potential, leading to the same outcome: a shift on the interatomic equilibrium distance.

APPLICATION TO EXPERIMENTAL DATA
While the empirical usefulness of Eq. 1 is widely accepted in glass science, here, we focus on testing the parallel validity of Eq. 8. This is challenging because literature provides a relatively broad scatter in the values of activation energy for Si-O bond rupture, ranging from approximately 50-150 kJ mol −1 (Lindsay et al., 1994;Walsh et al., 2000;Del Bene et al., 2003;Zhu et al., 2005). Such a variation may result in prediction errors of crack velocity of several orders of magnitude. Therefore, we used a highly accurate multilevel CCSD(T):MP2//DFT approach for generating data on the homogeneous dissociation of the Si-O bond. Specifically, we evaluated the energy of 624 kJ mol −1 for the central Si-O bond of a dimeric hydridosilsesquioxane molecule. This model is used to reflect the energetics of straining and breakage of a single silica bond at the crack tip because it is small enough to be computationally accessible, but provides a level of accuracy of <2 kJ mol −1 , unmatched by alternative methods on larger systems (Feller et al., 2011). Since the energy was used to create two "surfaces, " we will be using the value of 312 kJ mol −1 for comparison to experimental data on region III cracking. Calculations were also performed for Si-O bond cleavage in the presence of one, two, and three H 2 O molecules (mimicking region I cracking), yielding values of 164.2, 166.6, and 158.3 kJ mol −1 , respectively. No simulation was performed for more than three water molecules interacting with the DHS molecule; experimental data suggest that the hydration number of the Si-O-Si bond in amorphous silica lies between 2.4 and 1.8 (Fournier and Marshall, 1983). These values will be used in the following for reference.
Applying Eq. 8 to published experimental data on subcritical crack growth of silica (Suratwala and Steele, 2003) and a soda-lime silicate glass (Wiederhorn, 1967;Wiederhorn and Bolz, 1970), we first calculate the constant C. Theoretical calculations yield values of 1.13 × 10 −6 and 1.05 × 10 −6 m 4 (N mol) −1 , respectively, for the two different glasses. In parallel, rearranging Eq. 8 allows for the experimental determination of C, i.e., using the slope of a plot of T versus K I 2 , In Figure 3, the relation above is applied to subcritical crack growth data of soda-lime silicate glass measured between 2 and 90 • C (Wiederhorn and Bolz, 1970). The value of C does not remain constant through the complete range of temperature-stress intensity data, and even as it stabilizes there is a disparity of 4 orders of magnitude between the experimental fit and the theoretical calculation. This effectively means that the actual "conversion factor" of strain to thermal energy is much smaller than the theoretical one, which might be due to many different effects which are not covered by our approach: first, for lower crack growth velocities, the applied stress intensity factor is much closer to the threshold value K th [the static fatigue limit, below which the crack ceases to grow (Fett et al., 2005)]. This means that the stresses around the crack tip become increasingly shielded from the applied stresses due to diffusive exchange between mobile ions in the glass and the hydronium ions in the environment, leading to an increased slope of the T × K I 2 graph. Second, in parallel with the DFBM approach , the parameter C can be interpreted as a measure of structural disorder, related to how the applied stresses are preferentially concentrated on the bonds located closer to the crack tip. However, the stress concentration is assumed to follow a Gaussian distribution, and that might not reflect the actual distribution, leading to the disagreement between the experimental and theoretical values. It is also worth noting that Eq. 8 is derived under the assumption of static crack growth with constant K I , while the experimental tests were conducted under constant load. Therefore, following the observed disparity, we used the experimentally fitted values to construct the Arrhenius plots shown in Figure 4. FIGURE 3 | Determination of molar strain energy in subcritical crack growth. The temperature is plotted as a function of the square of the stress intensity factor at the crack tip for constant crack velocity. A linear fit (straight lines) yields the C-parameter. The inset shows the apparent dependence of C on crack velocity. Here, the dotted line is just a guide to the eye.
FIGURE 4 | Rate of subcritical crack growth as a function of total stored energy based on experimental data for soda-lime silicate glass on region I (right) (Wiederhorn and Bolz, 1970) and on region III (left) (Wiederhorn, 1967) and for silica glass on region I (Suratwala and Steele, 2003)  The authors are unaware of published crack growth rates of silica glass under vacuum; therefore, the comparison has to be made using crack growth in water vapor. Silica glass exhibits accelerated water diffusion under stress (Tomozawa et al., 1991), which causes the crack surfaces to expand, creating compressive stresses and shielding the crack tip from the applied stresses (Fett et al., 2005;Wiederhorn et al., 2011). This explains the apparently anomalous results reported for silica (Suratwala and Steele, 2003), and silica aerogels (Despetis et al., 2004), where the measured crack propagation speed decreases with increasing temperature. Wiederhorn et al. (2013b) have developed a formalism that enables the direct calculation of the shielding stresses as a function of both applied stresses and temperature. The estimated activation energies for region I are 114 ± 8 kJ mol −1 for the soda-lime silicate glass and 156 ± 5 kJ mol −1 for silica glass; for region III the activation energy for the soda-lime silicate is 314 ± 9 kJ mol −1 .
As commented before, on a first approximation of ideal brittleness (Célarié et al., 2003a,b;Wiederhorn, 2004, 2006;Wiederhorn et al., 2013a,b), the only contribution to the strain release rate G during fracture is the creation of the fracture surface. Therefore, with proper knowledge of the Si-O-Si surface bond density, the energy required to break each Si-O bond can be estimated by the surface energy. As its value depends only on the type of bond that is being broken, the calculated values for glassy and crystalline silica should be similar. The surface energy of silica glass under N 2 atmosphere and at 300 K is reported (Wiederhorn, 1969) to be approximately 4.4 J m −2 ; while MD simulations (Leed and Pantano, 2003) yield an area density of Si-O-Si bonds of 7.2 bonds nm −2 , resulting on an activation energy of 360 kJ mol −1 . This is remarkably consistent with monocrystalline α-quartz, with values of 310 and 340 kJ mol −1 [for cracks propagating on the a < 1120 > plane, on directions normal to the planes z < 0111 > and r < 1011 >, respectively] calculated using fracture (Atkinson, 1979;Timms et al., 2010) and bond density data (Bloss and Gibbs, 1963) under inert atmosphere.
The activation energy values calculated from cleavage in monocrystals and the simulation data obtained by CCSD(T):MP2//DFT agree to a remarkable extent with the results from Eq. 8, which leads us to believe that our proposed equation provides accurate predictions of crack propagation rates.
As noted initially, the principal difference between the two approaches rests on how the applied mechanical strain affects the chemical bonds directly at the crack tip. The offset approach [as seen in Eq. 1, but also revised elsewhere (Henderson et al., 1970;Bartenev, 1973)] considers that the applied stresses decrease the energy barrier for bond breakage. Some simulation studies (Gagnon et al., 2001) support this interpretation, but they also give a very different mechanistic description of crack propagation, following the coalescence of voids forming ahead of the crack tip; a model that seems to be at odds with the most recent experimental investigation (Wiederhorn et al., 2013a). A possible route for resolving this apparent discrepancy would be to combine these considerations with the aspect of concerted activation, in which the applied stresses play a role analogous to temperature, increasing the overall excitation level of the system without altering the energy profile of the chemical bonds (Stillinger, 1995;Sciortino, 2005). Thereby, it is important to acknowledge the empirical nature of both Eqs 1 and 8, if applied individually. We hypothesize that in reality, the physical reactions taking place may be best understood on the basis of simultaneous reactions of activation offset and concerted activation.

CONCLUSION
We considered parametrization of models for crack propagation in brittle glasses. The action of strain energy can be modeled in two empirical ways, i.e., acting either as an offset to the activation barrier for bond rupture, or concertedly with the thermal energy of the system. Both approaches provide accurate reproductions of experimental data. However, non-empirical treatment is expected to require simultaneous implementation.