Introduction

Yielding and fracture of amorphous materials are widely observed for hard1,2,3,4 and soft materials5,6. Material fatigue is a unique fracture mode in which structures fail when subjected to cyclic loading. It is known that fatigue failure occurs even when the applied strain (or stress) range is far below the yield strain (or stress) of the material. Because of this feature, we often use this type of fracture intentionally in our daily life when we intend to gently fracture a material like plastics by using relatively weak force. Fatigue is the most common type of mechanical structure failure and, thus, of significant importance in industrial applications of amorphous materials.

The fatigue fracture proceeds as follows. First, cyclic stress induces a microscopic fatigue crack, which gradually grows for each cycle, leading to the macroscopic failure of the material. Such fatigue-type failure is universally observed in various materials, covering metals, plastics, ceramics, and composites in both crystalline and amorphous forms. For crystalline materials, preexisting defects and grain boundaries are widely known to play a significant role in fatigue fracture initiation. On the other hand, the initiation mechanism of fatigue fracture has remained elusive for amorphous materials. The most critical question is why the fatigue fracture limit is apparently much smaller than the fracture strength, i.e., the yield or fracture stress for monotonic loading. This problem is critical in the damage design of various amorphous materials. For example, it has recently attracted particular attention in the field of bulk metallic glasses7,8,9,10,11.

Recently, the effect of periodic deformation has also attracted considerable attention in the physics community, based on numerical simulations12,13,14,15,16,17,18,19,20,21,22,23,24. Under periodic shear, amorphous solids transition from repetitive, predictable behaviour to irreversible behaviour13 of chaotic, irregular nature14,16,17 with an increase in the strain amplitude. The transition appeared sharp19. It was shown20,24,25 that the increase in the glass stability changes the nature of yielding under cyclic shear from ductile- to brittle-type. Furthermore, memory effects and interesting transitions in the particle trajectory were also revealed21,22,23,26,27,28,29.

Yielding can occur without accompanying catastrophic enhancement of density fluctuation, whereas fracture must accompany it. There are many pieces of evidence for significant density decrease and the resulting cavitation upon fracture of various amorphous materials (e.g., polymers30, rubbers31, oxides32 and metallic glasses33,34,35,36,37,38). Very recently, Sheng et al.39 succeeded in directly observing atomic-scale density alternation of deformed metallic glasses with transmission electron microscopy, revealing the autocatalytic generation of shear transformation zones. Molecular dynamics (MD) simulations can access microscopic information such as quadrupolar localized rearrangements, their coupling, microscopic plastic events, and shear band formation40,41,42,43,44. MD simulations also reported that density modulation and cavitation occur upon fracture for model glasses45,46, polymers47,48,49, silica50, and metallic glasses51,52,53,54,55,56. The occurrence of cavitation was also reported for fatigue fracture experiments in rubbers57 and metallic glasses58. Thus, considering the coupling of mechanical deformation to density fluctuation is essential for describing materials fracture. The shear transformation zone (STZ) theory is one of the most successful theories describing yielding behaviours59,60,61. It is a phenomenological theory based on microscopic T1 processes induced by deformation for amorphous metallic materials. For describing the fracture process after yielding, a condition for cavitation is needed to account for a crack opening along an operating shear band, but it requires some assumptions62,63,64. For example, Rycroft and Bouchbinder64 recently coupled a continuum STZ model with a condition for cavitation to describe the fracture of metallic glasses. They found that cavitation plays an essential role in the initiation of fracture. The STZ density should be related to the density, but the density is more suitable to describe cavitation since (1) it satisfies the mass conservation, and (2) the gas-liquid (glass) transition45,46 can be described by the well-established thermodynamic free energy65,66,67.

Here, we note that the fracture of materials must accompany the enhancement of density fluctuation and the resulting cavitation in a broad sense. Recently, molecular dynamics simulations have revealed many critical microscopic features, such as elastic and chemical heterogeneities, to impact fracture on a microscopic scale. This might give the impression that a phenomenological approach is not helpful in describing fracture. However, we emphasise that any prediction concerning the fracture of a material can be made only through its experimentally measurable macroscopic physical quantities. We can naturally describe cavitation in the framework of the coarse-grained phenomenological model, which consists of the mass conservation, momentum conservation, and material’s constitutive equations, and thermodynamic free energy describing the liquid (glass)-gas transition65,66,67. This approach cannot access the microscopic information such as particle trajectories, local orientational degrees of freedom, and particle-level avalanches (see, e.g., refs. 59,61). However, it provides a simple physical picture to understand the fundamental mechanical fracture mechanism and predicts the onset of instability from macroscopic measurable physical quantities. The power of the phenomenological approach was verified by the predicting ability of shear-induced cavitation (see ref. 66 and also below).

Here we consider fatigue fracture of amorphous materials under cyclic shear deformation on a phenomenological level. It is evident that extensional and compressional deformations accompany density change. On the other hand, simple shear (i.e., purely transverse) deformation does not involve any volume change directly, thus providing an impression of not causing density fluctuation. However, it was shown that this physical intuition breaks down even for simple liquids and glasses if the structural relaxation time and the elastic moduli depend on the material density66,67. Since supercooled liquids near the glass transition and amorphous materials are generally densely packed, these physical quantities should increase with increasing density. In other words, the elastic and relaxational properties have strongly asymmetric responses for the density increase and decrease. This situation is similar to what we experience on a packed train. Such dynamic and elastic asymmetries with respect to the density change lead to a nontrivial coupling between shear deformation and density fluctuation. For ordinary systems, higher density leads to higher elasticity and slower dynamics; however, for materials such as silica and water, it is the opposite (see, e.g.,25,68 for the difference). So, even simple shear deformation can induce material fracture by enhancing density fluctuation. As stressed above, density inhomogenisation is essential for fracture. This physical picture is reminiscent of the importance of the classical free volume concept in fatigue damage of bulk metallic glasses69,70, although the free volume is not necessarily a conserved quantity, unlike density.

For steady shear deformation, the following fracture criteria were obtained66,67: The liquid-type (ductile) fracture occurs when the shear rate \(\dot{\gamma }\) exceeds the critical value \({\dot{\gamma }}_{{{{{{{{\rm{c}}}}}}}}}={(\partial \eta /\partial p)}_{T}^{-1}\), where η is the shear viscosity, p is the pressure, and T is the temperature. Above \({\dot{\gamma }}_{{{{{{{{\rm{c}}}}}}}}}\), the enhancement of initially infinitesimal density fluctuation by shear-induced anisotropic compression and extension overwhelms its stabilisation by the isothermal compressibility. The \({\dot{\gamma }}_{{{{{{{{\rm{c}}}}}}}}}\) estimated from the experimentally measured (∂η/∂p)T qualitatively coincides with the critical shear rate for cavitation experimentally observed for lubricants66 and ductile fracture of metallic glasses67. This liquid-type fracture happens when tτ0, where τ0 is the structural relaxation time. In the opposite case of tτ0, the solid-type (brittle) fracture takes place when the strain γ exceeds the critical value \({\gamma }_{{{{{{{{\rm{c}}}}}}}}}={(\partial G/\partial {p}_{{{{{{{{\rm{L}}}}}}}}})}^{-1}\), where G is the shear modulus, and pL is the longitudinal pressure (see below on the precise definition). This criterion for brittle fracture, \({\gamma }_{{{{{{{{\rm{c}}}}}}}}}={(\partial G/\partial {p}_{{{{{{{{\rm{L}}}}}}}}})}^{-1}\), estimated from the experimentally measured P-dependence of G, is consistent with the critical strain for the brittle fracture of poly(methyl methacrylate)67. However, the predicted critical strain ( ~10%) was much larger than experimental ones (a few %) for metallic glasses67, which might be due to pre-existing defects or voids71,72,73 or nucleation-growth-type fracture, although it has to be examined carefully (see also the later discussion). We note that the intrinsic critical strain can be much larger than the experimental values obtained from macroscopic experiments, as supported by recent simulations and small-sample-size experiments on metallic glass fracture74,75. In the intermediate case, viscoelastic-type fracture takes place. Thus, the critical question of fatigue fracture is how materials fracture for oscillatory shear deformation; more precisely, the critical shear strain amplitude inducing fatigue fracture and the dependence of fracture behaviour on the cycle period.

This article focuses on fatigue fracture of amorphous materials under periodic shear deformation and aims to answer these fundamental questions by theoretical modelling and simulations. Our fundamental standpoint is to view material fracture as strain-induced enhancement of density fluctuation and the resulting void (or crack) formation and its development65,66,67. We have made a linear stability analysis of a coarse-grained fracture model and performed numerical simulations based on this model. We find that the fracture’s critical strain amplitude is about the same for monotonic and fatigue fractures, contrary to the common belief. We also reveal four types of fracture mode depending upon the cycle frequency, including a unique class of fatigue fracture (i.e., acoustic-type mode). Among them, the most important case is fatigue fracture of amorphous materials far below the glass transition temperature Tg, in which solid-type fracture takes place after many deformation cycles. In this case, the critical strain amplitude is given by γc, i.e., identical to the criterion for solid-type monotonic fracture in our definition. On the contrary, it has been widely believed that the critical strain amplitude γc is much smaller for fatigue fracture than monotonic fracture. This is because the practical fracture strain, γfrac, where the fracture occurs, has been regarded as the threshold strain where the macroscopic brittle fracture occurs. For yielding, the upper yield strain γyield is often used to characterise the yielding strength. We demonstrate that the critical strain γc for brittle fracture, the onset of irreversible deformation, is actually much smaller than γfrac. That is, even when the stress-strain relation apparently obeys the Hookean behaviour, the nonlinear effect already starts to cause irreversible deformation when γ exceeds γc (<γfrac). Therefore, the critical strain should be γc, not γfrac, for fracture under simple loading. These facts provide a clear physical explanation for why the critical strain of fatigue fracture looks much less than the apparent fracture strain of monotonic fracture. Our finding of the same critical strain amplitude for fatigue and monotonic fracture will provide crucial information for designing materials resistant to fatigue damage and risk management.

Results and discussion

Fundamental equations describing mechanical fracture

Here we consider fatigue fracture of amorphous materials based on our phenomenological theory of mechanical fracture66,67. As mentioned above, the critical point is that we view mechanical fracture as the development of spontaneous density fluctuations induced by mechanical deformation via the density dependence of rheological quantities, i.e., elastic moduli and relaxation times. In other words, we regard the mechanical fracture as a deformation-induced liquid (or solid)-to-gas transition, whose order parameter is density ρ. Local coupling between density and deformation is induced through the density dependence of the rheological properties. Here we consider simple shear deformation as a mechanical perturbation to a system, neglecting direct volume deformation. Then, the critical gross variables describing the fracture phenomena are density ρ, velocity v, and stress fields \(\overleftrightarrow{{{{{\boldsymbol{\sigma}}}}}}\), which obey the mass conservation, momentum conservation, and upper-convective Maxwell-type constitutive equations, respectively76:

$$\frac{\partial \rho }{\partial t}=-{\nabla }_{i}(\rho {v}_{i}),$$
(1)
$$\rho \left(\frac{\partial }{\partial t}+{v}_{j}{\nabla }_{j}\right){v}_{i}=-{\nabla }_{j}{{{\Pi }}}_{ij}+{\nabla }_{j}{\sigma }_{ij}+{\eta }_{0}{\nabla }_{j}^{2}{v}_{j},$$
(2)
$$\left(\frac{\partial }{\partial t}+{v}_{j}{\nabla }_{j}\right){\sigma }_{ij}-({\sigma }_{ik}{\nabla }_{k}{v}_{j}+{\sigma }_{kj}{\nabla }_{k}{v}_{i})=G(\rho )({\nabla }_{i}{v}_{j}+{\nabla }_{j}{v}_{i})-\frac{1}{\tau (\rho )}{\sigma }_{ij},$$
(3)

where \(\overleftrightarrow{{{{{\boldsymbol{\Pi }}}}}}\) is the osmotic stress. We assume that a system obeys the equation of the state for a van-der-Waals liquid (see Methods). In Eq. (2), \({\eta }_{0}{\nabla }_{j}^{2}{v}_{j}\) represents a small intrinsic viscous dissipation term due to short time- and small length-scale dynamics before the coarse-graining (see Methods for the detail). We note that the significant dissipation comes from viscoelastic relaxation represented by the upper-convective Maxwell equation Eq. (3). The nonlinear mechanical response, e.g., plasticity, emerges from the density dependences of the shear modulus G(ρ) and relaxation time τ(ρ) in our model (see Methods, Supplementary Note 1, and Supplementary Fig. 1 on their density dependences). Here we point out that the ρ-dependences of G(ρ) and τ(ρ) affect the values of the critical strain rate and strain but do not affect the fracture behaviours qualitatively. The decrease of the shear modulus and relaxation time with decreasing the density, i.e., dG(ρ)/dρ, dτ(ρ)/dρ > 0, is the origin of the self-catalytic enhancement of density fluctuation upon shear deformation (see below). Then, we consider fatigue fracture of an amorphous material (the shear relaxation time τ0) induced by cyclic shear deformation of \(\gamma \sin (\omega t)\) (t: duration time). We can show by linear stability analysis (see below) that there are four essential cases, depending on the relationship between the three characteristic time scales, t, 1/ω, and τ0 (see Table 1): (i) liquid instability for a case of ωτ0 1. (ii) solid instability for a case of ωτ0 1. (iii) viscoelastic instability for a case of ωτ0 1. (iv) acoustic instability for an elastic limit (τ0 → ).

Table 1 Criteria of fatigue fracture. Here τL = τ0(1 + 2G0KT), where G0 is the shear modulus and KT is the isothermal compressibility, and cL is the velocity of longitudinal sound wave (see below for the details).

We put small-amplitude Gaussian density fluctuation in the initial state77. We also neglect the effect of thermal and mechanical noises44,78 for simplicity. Thermal noise should play a critical role in the nucleation-and-growth-type cavitation78. However, the nature of the density fluctuation, the initial spatial heterogeneity of the material’s properties, and the thermal fluctuation have little influence on the results for spinodal-type cavitation considered here (see Supplementary Note 3, Supplementary Figs. 3 and 4). This is because even infinitesimal density fluctuation should be enhanced drastically above the stability limit and decay otherwise for the spinodal type. We may say that our prediction for the spinodal-type instability provides the upper limit of the material stability under oscillatory shear deformation.

Linear stability analysis

First, we make the linear stability analysis of Eqs. (1)-(3) under oscillatory shear deformation. Suppose that a system is under homogeneous oscillatory shear deformation (the flow direction x, the shear gradient direction y, and the angular frequency ω), whose deformation velocity vx(y) is given by

$${v}_{x}(y)=\gamma \omega \cos (\omega t)y.$$
(4)

Then, the stress field is given by

$$\langle {\sigma }_{xy}(t)\rangle ={G}_{0}\gamma \cdot \frac{\omega {\tau }_{0}}{1+{(\omega {\tau }_{0})}^{2}}\left[\cos \omega t+(\omega {\tau }_{0})\sin \omega t-{e}^{-t/{\tau }_{0}}\right],$$
(5)

where G0 = G0(ρ0) and τ0 = τ(ρ0) are the shear modulus and relaxation time of the initial homogeneous state of density ρ0, respectively.

Next, we consider a case of a system with a small-amplitude density fluctuation δρ. Hereafter we neglect the inertia and advection terms and assume γ 1. Then, we obtain the following time-evolution equations for density, velocity and stress fields in the linear order:

$${\rho }_{k}^{\prime}(t)+{\rho }_{0}{Z}_{k}=0$$
(6)
$$\frac{{k}^{2}}{{\rho }_{0}{K}_{T}}{\rho }_{k}+{W}_{k}=0$$
(7)
$$\left(\frac{\partial }{\partial t}+\frac{1}{{\tau }_{0}}\right){W}_{k}=-2{G}_{0}{k}^{2}{Z}_{k}-2{k}_{x}{k}_{y}{\rho }_{k}\left(G^{\prime} ({\rho }_{0})\gamma \omega \cos \omega t+\frac{\tau ^{\prime} ({\rho }_{0})}{{\tau }_{0}^{2}}\langle {\sigma }_{xy}(t)\rangle \right),$$
(8)

where KT is the isothermal compressibility of the system, ρk is the Fourier component of δρ with wavenumber k, Z = ivi, W = ijσij, and Zk and Wk are the Fourier component of Z and W, respectively.

Then, we obtain the time evolution of ρk as follows:

$${\rho }_{k}^{\prime}+\frac{1}{{\tau }_{0}(1+2{G}_{0}{K}_{T})}\left[1-\sin 2\theta \left({\tau }_{0}\left(\frac{\partial G}{\partial p}\right)\gamma \omega \cos \omega t+\frac{1}{{\tau }_{0}}\left(\frac{\partial \tau }{\partial p}\right)\langle {\sigma }_{xy}\rangle \right)\right]{\rho }_{k}=0,$$
(9)

where θ is the angle from the flow direction x. By considering the extensional direction (\(\sin 2\theta =1\)) and setting τL = (1 + 2G0KT)τ0, ρk is obtained as

$${\rho }_{k}(t) = \,{\rho }_{k}(0)\exp \left[-\frac{t}{{\tau }_{{{{{{{{\rm{L}}}}}}}}}}+\frac{\gamma \sin \omega t}{{\gamma }_{{{{{{{{\rm{c}}}}}}}}}}\right] \\ \cdot \exp \left[\frac{1}{{\tau }_{0}}\cdot \frac{\gamma }{1+{(\omega {\tau }_{0})}^{2}}\left(\frac{\partial \tau }{\partial {p}_{{{{{{{{\rm{L}}}}}}}}}}\right){G}_{0}\left[\sin \omega t-\omega {\tau }_{0}(\cos \omega t-{e}^{-t/{\tau }_{0}})\right]\right].$$
(10)

In the above, the longitudinal pressure in a quiescent state (i.e., \(\dot{\gamma }=0\)), pL, is given as (1 + 2G0KT)p.

Unlike the case of constant shearing, oscillatory shearing does not cause any instability in this linear regime. However, we can extract a necessary condition for fatigue fracture from the linear analysis. The key is whether the strain exceeds the stability limit during oscillation or not. Once the strain exceeds the critical value, then nonlinear, irreversible growth of density fluctuation occurs for a certain period of each cycle but then starts to decay in the subsequent recovery period. The linear response of density fluctuation to oscillatory shear is perfectly reproducible, which never contributes to the instability. However, the growth of density fluctuation induced by the nonlinearity can never be recovered entirely and must be accumulated. Thus, we can predict the critical threshold strain amplitude γc in the following linear stability analysis.

Using Eq. (10), we can discuss under what condition the mechanical perturbation enhances density fluctuation. There are three characteristic time scales in this problem: t, 1/ω, and τ0. Now we discuss four essential cases (see Table 1). Here we do not consider a case of t 1/ω and tτK (τK: the bulk (volume) relaxation time) since this case is unrealistic for fatigue fracture.

First, we consider a case of ωτ0 1. When tτ0, ρk(t) is given by

$${\rho }_{k}(t) \simeq {\rho }_{k}(0)\exp \left[-\frac{t}{{\tau }_{{{{{{{{\rm{L}}}}}}}}}}\right]\cdot \exp \left[\gamma \sin \omega t\left(\frac{1}{{\gamma }_{{{{{{{{\rm{c}}}}}}}}}}+\frac{{G}_{0}}{{\tau }_{0}}\cdot \left(\frac{\partial \tau }{\partial {p}_{{{{{{{{\rm{L}}}}}}}}}}\right)\right)\right]\\ ={\rho }_{k}(0)\exp \left[-\frac{t}{{\tau }_{{{{{{{{\rm{L}}}}}}}}}}+\frac{\gamma \sin \omega t}{{\tau }_{{{{{{{{\rm{L}}}}}}}}}{\dot{\gamma }}_{{{{{{{{\rm{c}}}}}}}}}}\right].$$
(11)

Furthermore, when t 1/ω (i.e., τLt 1/ω), ρk becomes

$${\rho }_{k}(t)={\rho }_{k}(0)\exp \left[\frac{t}{{\tau }_{{{{{{{{\rm{L}}}}}}}}}}\left(\frac{\gamma \omega }{{\dot{\gamma }}_{{{{{{{{\rm{c}}}}}}}}}}-1\right)\right].$$
(12)

Thus, the enhancement of density fluctuation happens if

$$\gamma \omega \, \gtrsim \, {\dot{\gamma }}_{{{{{{{{\rm{c}}}}}}}}}\cdot \left(1+\frac{{\tau }_{{{{{{{{\rm{L}}}}}}}}}}{t}\right) \sim {\dot{\gamma }}_{{{{{{{{\rm{c}}}}}}}}}.$$
(13)

This condition is the instability criterion for viscous fluids. We refer to this type of instability as liquid instability.

Secondly, we consider a case of ωτ0 1. We consider the case of t 1/ω. Then, ρk(t) is given by

$${\rho }_{k}(t)= \,{\rho }_{k}(0)\exp \left[-\frac{t}{{\tau }_{{{{{{{{\rm{L}}}}}}}}}}+\frac{\gamma \sin \omega t}{{\gamma }_{{{{{{{{\rm{c}}}}}}}}}}\right]\\ \cdot \exp \left[\frac{1}{{\tau }_{0}}\cdot \frac{\gamma }{\omega {\tau }_{0}}\left(\frac{\partial \tau }{\partial {p}_{{{{{{{{\rm{L}}}}}}}}}}\right){G}_{0}\left({e}^{-t/{\tau }_{0}}-\cos \omega t\right)\right].$$
(14)

The enhancement of density fluctuation should take place when the exponent of each exponential function is positive. The condition for the first exponent to be positive is that for tτL (i.e., 1/ωtτ0),

$$\gamma \;\gtrsim \;{\gamma }_{{{{{{{{\rm{c}}}}}}}}}.$$
(15)

The condition for the second exponent to be positive is then

$$\gamma \;\gtrsim \;\frac{\omega {\tau }_{0}}{\frac{{G}_{0}}{{\tau }_{0}}\left(\frac{\partial \tau }{\partial {p}_{{{{{{{{\rm{L}}}}}}}}}}\right)}={\gamma }_{{{{{{{{\rm{c}}}}}}}}}\cdot \frac{\omega {\tau }_{0}}{\frac{{\gamma }_{{{{{{{{\rm{c}}}}}}}}}}{{\dot{\gamma }}_{{{{{{{{\rm{c}}}}}}}}}{\tau }_{{{{{{{{\rm{L}}}}}}}}}}-1}.$$
(16)

The second condition (16) is satisfied more easily than (15) when

$$\frac{{\gamma }_{{{{{{{{\rm{c}}}}}}}}}}{{\dot{\gamma }}_{{{{{{{{\rm{c}}}}}}}}}{\tau }_{{{{{{{{\rm{L}}}}}}}}}}-1\;\gtrsim \;\omega {\tau }_{0}\gg 1.$$
(17)

This condition becomes more difficult to be satisfied with an increase in ω. The condition (15) is the instability criteria for elastic solids. We refer to this type of instability as solid instability.

Thirdly, we consider a case of ωτ0 1. In this case, ρk(t) is expressed as

$${\rho }_{k}(t) \simeq\, {\rho }_{k}(0)\exp [-\omega t]\cdot \exp \left[\frac{\gamma \sin \omega t}{{\gamma }_{{{{{{{{\rm{c}}}}}}}}}}\right]\\ \cdot \exp \left[\frac{1}{{\tau }_{0}}\cdot \frac{\gamma {G}_{0}}{2}\left(\frac{\partial \tau }{\partial {p}_{{{\mbox{L}}}}}\right)(\sin \omega t-\cos \omega t+{e}^{-\omega t})\right].$$
(18)

For ωt ~ 1, the instability condition is given by

$$\gamma \;\gtrsim \;\min ({\gamma }_{{{{{{{{\rm{c}}}}}}}}},{\dot{\gamma }}_{{{{{{{{\rm{c}}}}}}}}}{\tau }_{{{{{{{{\rm{L}}}}}}}}}).$$
(19)

Then, using ωτ0 1, it can be rewritten as

$$\gamma \omega \;\gtrsim \;\min ({\gamma }_{{{{{{{{\rm{c}}}}}}}}}/{\tau }_{{{{{{{{\rm{L}}}}}}}}},{\dot{\gamma }}_{{{{{{{{\rm{c}}}}}}}}}).$$
(20)

This instability condition is a consequence of the crossover between viscous and elastic cases; thus, it may be called viscoelastic instability.

Lastly, we consider a case of elastic limit (τ0 → ). For this elastic limit, the time evolution equation of ρk, including the inertia term, is given by

$${\rho }_{k}^{{\prime\prime} }+{k}^{2}{c}_{{{{{{{{\rm{L}}}}}}}}}^{2}\left[1-\sin 2\theta \cdot \frac{\gamma \sin \omega t}{{\gamma }_{{{{{{{{\rm{c}}}}}}}}}}\right]{\rho }_{k}=0,$$
(21)

where cL is the velocity of longitudinal sound wave. Equation (21) has the same form as the Mathieu equation describing parametric resonance79. This tells us that for γ ≠ 0, the system must become unstable for a wavenumber k satisfying kcL ~ ω. Practically, however, since there is a lower bound for k associated with the system size, this instability occurs when ω is rather high. For a system of size L, this condition for ω is given by

$$\omega \;\gtrsim \;\frac{2\pi {c}_{{{{{{{{\rm{L}}}}}}}}}}{L}.$$
(22)

We refer to this type of instability as acoustic instability.

Simulation results for fracture under simple shear deformation

Now we show the results of numerical simulations. We perform simulations by numerically solving the phenomenological equations, Eqs. (1)–(3), including nonlinear effects, for a two-dimensional system, under a periodic boundary condition. The initial condition and the parameter settings are described in Methods.

First, we study the mechanical response of a system to simple shear deformation with a constant shear rate, \(\dot{\gamma }\). Figure 1a shows fracture behaviours, or temporal development of density fluctuation, under a slow steady shear deformation of \(\dot{\gamma }=0.001\). We can see the growth of density fluctuation with an increase in strain, γ. This case corresponds to ductile fracture. We note that the density inhomogeneity is directly related to the inhomogeneities of the shear modulus G and relaxation time τ in our model since these material quantities only depend on the normalised density ϕ. The decreases of G and τ with a reduction of ϕ lead to the self-catalytic growth of density fluctuation under shear deformation: Lower density regions become softer and faster, i.e., easier to deform. Such behaviour has recently been directly observed in metallic glasses by transmission electron microscopy39.

Fig. 1: Fracture behaviour under simple shear.
figure 1

a Enhancement of density fluctuation as a function of γ for \(\dot{\gamma }=0.001\). The colour in the images reflects local density (see the colour bar). b Stress-strain curves for various \(\dot{\gamma }\). The purple dashed straight line is the stress-strain curve for an ideal elastic body. c γ-dependence of the variance of density fluctuation, 〈(δϕ)2〉. The vertical dashed lines in panels b and c indicate the critical strain γc. The green and red arrows indicate the upper yield strain γyield for \(\dot{\gamma }=0.00005\), 0.00010 and 0.00020 and the fracture strain γfrac for \(\dot{\gamma }=0.00020\), 0.00050, and 0.00100. For \(\dot{\gamma }=0.00200\), the fracture occurs at γ slightly beyond 0.5. Here, the characteristic strain, γchar, is defined as the strain at 〈(δϕ)2〉 = 0.06 in Supplementary Note 2 and Supplementary Fig. 2.

Next, we show in Fig. 1b the stress-strain curves for various \(\dot{\gamma }\). As we discussed above, the curve approaches the linear Hookian line σ = G0γ for higher \(\dot{\gamma }\), whereas exhibiting the yielding behaviours for lower \(\dot{\gamma }\). The fracture is brittle for high \(\dot{\gamma }\), whereas ductile for low \(\dot{\gamma }\)67. We also show how density fluctuation grows with an increase in γ for various \(\dot{\gamma }\) in Fig. 1c. We can see that lower \(\dot{\gamma }\) leads to faster growth of density fluctuation as a function of γ, which is a consequence of the fact that viscous-liquid-type instability becomes effective for low \(\dot{\gamma }\). In these figures, we also show the critical strain value, γc, by the vertical dashed lines. Moreover, we define the upper yield strain, γyield, and the fracture strain, γfrac, as the strains at the first local-maximum stress and at actual fracture, respectively. The upper yield and fracture strains are indicated by the arrows in Fig. 1b. We can see that the upper yield and fracture strains, γyield and γfrac, are larger than the critical strain, γc. We also confirm that the characteristic strain, γchar, defined as the strain at 〈(δϕ)2〉 = 0.06 is larger than γc (see Fig. S2). This is consistent with the results of previous molecular dynamics simulations that plastic events already occur much before the yielding or fracture limit (see, e.g., refs. 40,41,42,43,44). For solid-like fracture occurring at large \(\dot{\gamma }\), even above the critical strain γc for which density fluctuation is enhanced irreversibly, the stress-strain curves almost show the Hookean behaviour. This indicates that it is challenging to detect fracture onset (i.e., irreversible deformation) experimentally from the stress-strain curve. However, we might be able to detect it by measuring acoustic emission with high sensitivity if the local deformation rate is fast enough.

Simulation results for viscous- and viscoelastic-type fatigue fracture

First, we show the behaviour of fatigue fracture under oscillatory shear deformation \(\tilde{\gamma }(t)=\gamma \sin \omega t\) with ωτ0 = 0.10 in Fig. 2. In this case, density fluctuation grows, and cavitation eventually occurs in a correlated manner perpendicularly to the extensional direction after some incubation time, leading to fracture. This cavity chain formation can be regarded as a manifestation of liquid-type instability.

Fig. 2: Fatigue fracture of viscous type.
figure 2

This is observed for γ = 0.38 and ωτ0 = 0.10 (see also Supplementary Movie 1). The meaning of the colour is the same as Fig. 1a.

Next we show fatigue fracture observed for ωτ0 = 1.0 in Fig. 3. We can see that cracks are formed perpendicularly to the extensional direction (along both ~45 and 135) and gradually grow in size. There seems to exist a characteristic length scale for cracks. This fracture behaviour belongs to the viscoelastic type.

Fig. 3: Fatigue fracture of viscoelastic type.
figure 3

This is observed for γ = 0.14 and ωτ0 = 1.0 (see also Supplementary Movie 2). The meaning of the colour is the same as Fig. 1a.

Simulation results for elastic-type fatigue fracture

Next, we show the elastic type fatigue fracture under periodic shearing given by \(\gamma \sin \omega t\). Figure 4a shows the process of fatigue fracture for γ = 0.22 and ωτ0 = 31.1. There is no sign of fracture in the initial several cycle periods, but density fluctuation starts to grow after t 5T (T: the period of oscillation). Once the enhancement starts, fracture proceeds rather quickly, accompanying crack formation in both ~45 and 135 directions, eventually followed by fine crack formation, i.e., shear band formation along 0 or 90 from the shear direction.

Fig. 4: Fatigue fracture of elastic type.
figure 4

a Time evolution of density fluctuation in fatigue fracture of elastic type. This is observed for γ = 0.22 and ωτ0 = 31.1 (see also Supplementary Movie 3). The meaning of the colour is the same as Fig. 1a. Initially, density fluctuation and micro-cracks are induced along 45° and 135°, and eventually, their mechanical couplings lead to the shear band formation along  0° and 90° (see white and black arrows). b A crack opening process in fatigue fracture of elastic type observed at γ = 0.21 (see also Supplementary Movie 4). We can see the cooperative formation of micro-crack array (see the circled regions) and its development towards the final shear-band formation along 90°. c Growth of the variance of density fluctuation, 〈(δϕ)2〉, with time under oscillatory shear with ωτ0 = 31.1 for six different values of γ. d The dependence of the number of cycles at the fracture onset, Nonset, on Δγ = γ − γc. The solid line has a slope of −3.

This process of the shear band formation is similar to the one observed in MD simulations. For example, Dasgupta et al.43 showed that shear deformation activates Eshelby inclusions when their four-leaf patterns align around 45 or 135 from the shear direction and the coupling between them leads to the formation of the shear band along 0 or 90 degrees. Our simulations do not consider Eshelby inclusions since they are based on the coarse-grained model. However, we point out that there is a similarity: mechanical coupling between micro-cracks. Shear deformation creates low density regions along 45 or 135 from the shear direction under an oscillatory deformation (see Fig. 4a around t = 6T ~ 7T. When the fracture is liquid-like or viscoelastic fracture, the opened cracks align either 45 or 135 degrees (see Figs. 2 and 3). This is because the oscillatory deformation does not provide enough time for the cracks to be advected by shear flow, unlike monotonic shear deformation. However, when the fracture is solid-like, micro-cracks form due to their mechanical couplings and eventually align these low-density regions in the form of shear bands along either 0 or 90, as shown in the final panel of Fig. 4a. We can see the process of cooperative formation of micro-cracks and their development to shear band more clearly in Fig. 4b. We note that all fatigue fractures of amorphous solids belong to the solid-like case, i.e., ωτ0 1, and thus, shear bands are formed along 0 or 90, as observed experimentally (see, e.g., ref. 11).

We also show the time dependence of the variance of density fluctuation, 〈(δϕ)2〉, in Fig. 4c. The initial density fluctuation monotonically decays with time for a case without shearing (γ = 0). For high enough γ, we can see that a system becomes unstable rather quickly. For γ = 0.20, for example, density fluctuation starts to grow after several tens of oscillations. We also analyse the relationship between the excess strain, Δγ = γ − γc, and the number of cycles at the onset of instability, Nonset, in Fig. 4d. We find the following power-law relation:

$${N}_{{{{{{{{\rm{onset}}}}}}}}}\propto {{\Delta }}{\gamma }^{-\nu },$$
(23)

with ν = 3. This relation indicates the critical-like divergence of Nonset towards γ → γc. Because of the computational cost, we cannot approach γ close enough to γc, for which Nonset should be large; thus, our simulations are limited to low-cycle fatigue. We expect that this relation holds even for high-cycle fatigue, as observed experimentally. This point needs to be confirmed in the future.

Here we note that such stress (strain)-life relation has been universally observed in various amorphous materials such as metallic11 and polymeric glasses80. The phenomenological relation has been known as a fatigue-life model, Model I. It has been associated with the so-called two-parameter Weibull function based on a bimodal failure picture, i.e., the statistical origin of fatigue fracture (see, e.g., ref. 11). In our terminology, the two parameters are γc and ν. This model has been widely used to analyse and predict the fatigue life of metallic glasses.

On the other hand, our model provides an entirely different physics behind Eq. (23). In our linear stability analysis, we ignore the nonlinear and advection terms in the momentum conservation equation, which may play a critical role in the fatigue lifetime. In particular, there is always a relaxation path of density fluctuation with a diffusive origin (see the first exponential term in Eq. (10)). Thus, the competition between suppression and enhancement of density fluctuation under the influence of the advection and nonlinear terms present in our simulations may be the origin of the delayed growth of density fluctuation after the incubation period of many oscillation cycles. This point needs to be studied in detail in the future.

Simulation results for acoustic-type fatigue fracture

Finally, we show the case of very high-frequency shearing. Figure 5 shows fatigue fracture observed for γ = 0.14 and ωτ0 = 160.

Fig. 5: Fatigue fracture of acoustic type.
figure 5

This is observed for γ = 0.14 and ωτ0 = 160 (see also Supplementary Movie 5). The meaning of the colour is the same as Fig. 1a.

The fracture behaviour looks somewhat similar to that in Fig. 4, but this type of fracture can take place for a smaller amplitude of γ. The bands of micro-cracks are formed due to their mechanical couplings under the cyclic deformation. Although cracks are formed along y-direction in this particular case, we confirm by large-system-size simulations that cracks can be formed along other directions, including 45 and 135, a characteristic feature of acoustic fracture.

Instability threshold map for fatigue fracture

Here we show the instability threshold for fatigue fracture in the γ-ωτ0 plane (see Fig. 6). We plot the theoretically estimated threshold value of γ, above which fatigue fracture should occur against the angular frequency ω of oscillatory shear. When γ is larger than the threshold value, density fluctuation grows, leading to mechanical instability; otherwise, density fluctuation decays after a long time. The red dashed curve shows the threshold for viscous fracture, the green dashed line is for elastic fracture, and the black vertical arrow indicates the angular frequency above which the acoustic fracture occurs.

Fig. 6: Instability threshold map for fatigue fracture.
figure 6

The filled circles are the threshold γ estimated from numerical simulations (red-filled circles: liquid-type; blue-filled circles: viscoelastic-type; green-filled circles: solid-type; violet-filled circle: acoustic-type). The red dashed curve (\({{\dot{\gamma }}_{{{{{{{{\rm{c}}}}}}}}}/\omega}\)), blue dashed horizontal line (\({\dot{\gamma}}_{{{{{{{{\rm{c}}}}}}}}}{\tau }_{{{{{{{{\rm{L}}}}}}}}}\)), and green dashed horizontal lines (\({{\gamma}_{{{\rm{c}}}}}\)) are the theoretically predicted thresholds for viscous-type, viscoelastic-type, and elastic-type instability, respectively. The black vertical arrow indicates the threshold angular frequency above which acoustic-type instability should occur.

In Fig. 6, we also show the results of our numerical simulations. We can see that the theoretical predictions (dashed curves and horizontal lines) based on the linear stability analysis are consistent with those (circular symbols) determined from the simulations. For a low angular frequency (ωτ0 1), viscous-type fracture occurs, obeying γ 1/ω. For a high angular frequency (ωτ0 1), elastic-type fracture occurs, obeying γγc. For an intermediate angular frequency (ωτ0 1), fracture occurs at smaller γ, reflecting the competition between viscous-type and elastic-type fractures. Finally, we can see that fracture occurs at smaller γ again for a very high angular frequency. This frequency is close to the threshold angular frequency for acoustic-type instability: ω 2πcL/L (see the vertical black arrow).

Conclusion

To summarise, we have studied, based on the phenomenological density-based coarse-grained model, how fatigue fracture occurs theoretically and numerically. We have estimated the critical strain amplitude for fatigue fracture using the linear stability analysis of the fundamental equations describing the time evolution of density, velocity, and stress fields. The coupling between the density and the rheological properties leads to the self-catalytic enhancement of density fluctuation and the resulting mechanical fracture. We have also numerically solved the equations and estimated the threshold strain amplitude from the simulations. We have found reasonable agreements between the theoretical predictions based on the linear stability analysis and the values estimated from simulations with nonlinearity, as summarised in Fig. 6. We have discovered four types of fatigue fracture, i.e., liquid-, viscoelastic-, elastic-, and acoustic-type fracture modes. Our study provides a simple physical picture of how and why amorphous materials fracture under oscillatory shear deformation. Furthermore, our theory allows us to predict the critical strain amplitude for various loading conditions, i.e., ωτ0, from macroscopic physical quantities experimentally measurable.

Most fatigue fracture relevant to material science belongs to the elastic type. This is because amorphous materials are usually used at a temperature much lower than the glass-transition temperature Tg, where the structural relaxation time τ0 is enormously long. Thus, the condition ωτ0 1 is typically satisfied. If we compare the apparent fracture strain limit for monotonic deformation, γfrac, with the fatigue strain limit γc, the former is significantly larger than the latter (see Fig. 1b). Roughly speaking, for the oscillatory strain whose amplitude is less than γc, there is no fatigue fracture, and the deformation is reversible. Once γ exceeds γc, the fatigue fracture should take place, but the number of cycles required for fracture increases in a diverging manner when γ → γc, as shown in Fig. 4d. The equivalence of γc for fatigue and monotonic fracture indicates that we can infer the critical physical parameters of fatigue fracture from simple monotonic fracture experiments without long-term fatigue fracture experiments, which should have a significant practical impact.

Finally, we note some issues to be considered for quantitative comparison with actual experiments. In the present study, we consider fracture under homogeneous deformation (a homogeneous mechanism)11,81. However, the surface-specific mechanism may be critical in actual experiments since small surface cracks or defects (e.g., voids or preexisting cracks) generally exist in a sample. Thus, crack-like shear-band propagation from preexisting defects (a heterogeneous mechanism) also plays a critical role. For example, an experimental work by Lei et al.73 showed significant roles of preexisting cavities in fracture of metallic glasses. Furthermore, the presence of nanocavities was observed for metallic glasses via electron microscopy71, which may be responsible for more than two types of shear band nucleation sites that operate at different stress levels72. These observations and/or the nucleation-growth-type fracture mechanism (not considered here) may explain why the critical strain (~10%) predicted by our model67 is much larger than the experimental ones (a few %) for metallic glasses. Cheng and Ma also showed that the intrinsic critical strain can be much larger than the experimental values obtained from macroscopic experiments. Small-sample-size experiments of mechanical fracture of metallic glasses have strongly supported such a possibility (see ref. 75 for review). We have also confirmed that the introduction of the preexisting density fluctuation beyond the thermal-noise level, δϕ, substantially lowers the upper yield strain and the fracture strain for ductile and brittle fractures (see Figs. S3 and S4, respectively). Thus, avoiding such mechanical defects in glasses is crucial to estimating the material-specific fracture strength. Generally, the presence of such defects or cavities can be checked much more sensitively in transparent glasses, such as polymer and silicate glasses, using a light scattering technique. For example, optical fibres made of polymer or silicate glasses must be free from such defects and may be suitable for the fracture study free from defects.

Another factor to be considered is the spatial fluctuation of γc. Recent experiments82 and simulations83,84 showed that elastic moduli are spatially fluctuating in glasses. Soft regions85 (with anisotropy) trigger plastic events under shear86,87, which may also lead to the fluctuation of γc. In such a case, the minimum value of local γc should be the critical strain for fracture, but the effect may be minor. We also need to consider nucleation-growth-type fracture under thermal fluctuation. Finally, this article focuses on simple shear deformation, but various deformation modes accompanying volume deformation are employed in real experiments. The deformation mode inducing volume change should induce fracture much more easily than the simple shear mode since it directly enhances density fluctuations and void formation. We may need to consider elastic deformation’s effect on the apparent density change for such a case. For solids, the relevant free-volume-like variable m is not determined by density change δρ alone but affected by the elastic volume deformation as m = δρ/ρ0 + u65,88,89, where u is the deformation vector. This is because the density change is coupled to the elastic volume deformation in solids, unlike liquids. This may not be so important for shear deformation we consider here, but critical when the material volume is deformed. All these issues are challenging yet interesting topics for future study.

Our model predicts that the onset of irreversible deformation should be the same for fatigue and monotonic fractures. Currently, we do not have any experimental evidence for this prediction. However, (1) our theoretical estimation of the onset of irreversible deformation exists in the apparently linear regime much before actual yielding or fracture occurs, which has also been seen in molecular dynamics simulations (see, e.g., refs. 90,91,92), and (2) the fatigue fracture magnitude is much smaller than that for simple monotonic loading. These two facts are at least qualitatively consistent with our claim. The key issue is how to detect the onset of irreversible deformation experimentally. The acoustic emission is one possible way to detect it sensitively but might not be so helpful if the irreversible deformation occurs slowly. This poses an interesting question concerning the sensitive detection of the onset of irreversible density enhancement under deformation. Optically transparent amorphous solids such as optical fibres may be helpful for this purpose since the onset of irreversible deformation may be detected by light transmission or scattering. Once we can determine the onset of nonlinear deformation experimentally by simple loading, we can infer the critical strain of fatigue fracture without long-term oscillatory experiments. We hope our study will shed fresh light on the fundamental physical understanding of fatigue fracture and the design of materials resistive to fatigue fracture.

Methods

Simulation method

We use the Ginzburg-Landau-type free-energy functional for the liquid and glass states:

$${{{{{{{\mathcal{F}}}}}}}}=\int dV\left[f(\rho )+\frac{C}{2}| \nabla \rho {| }^{2}\right],$$
(24)

where the gradient term arises from the energy cost associated with the spatial density variation, and C is a positive constant. Then, the pressure tensor is given by Πij = (ρf/∂ρ − f − Cρ2/2 − Cρ2ρ)δij + Ciρjρ. To describe the gas-liquid phase transition properly, we use the van der Waals formula \(f(\rho )=\alpha \{\rho \ln [\rho /({\rho }_{{{{{{{{\rm{c}}}}}}}}}-\rho )]-\epsilon {\rho }^{2}\}\), where α and ϵ are positive constants, and ρc is the density at the closed-packed state.

We assume a simple density dependence of the shear elastic modulus: \(G(\rho )={G}_{0}\{\tanh [{{{\Lambda }}}_{{{{{{{{\rm{G}}}}}}}}}(\rho /{\rho }_{{{{{{{{\rm{G}}}}}}}}}-1)]+1\}+0.05{G}_{0}\), where G0, ΛG, and ρG are positive constants. The structural relaxation time is assumed to be given by \(\tau (\rho )={\tau }_{{{{{{{{\rm{B}}}}}}}}}\{\exp [\mu (\rho )\rho /({\rho }_{{{{{{{{\rm{c}}}}}}}}}-\rho )]-1\}+0.05{\tau }_{{{{{{{{\rm{B}}}}}}}}}\), where \(\mu (\rho )={{\Gamma }}\tanh [{{{\Lambda }}}_{\tau }({\rho }_{{{{{{{{\rm{c}}}}}}}}}/\rho -1)]\), τB, and Λτ are positive constants. In the above, 0.05G0 and 0.05τB are the shear modulus and relaxation time in the low-density limit, respectively. We note that this density dependence mimics the Vogel-Fulcher-Tammann law, observed in real systems93,94. We used G0 = 15, ρG/ρc = 0.835, ΛG = 28.4, τB = 200, Λτ = 10.5 and Γ = 0.48 for numerical simulations (see Supplementary Fig. 1 for the density dependences of G(ϕ) and τ(ϕ)).

Length and time are measured in units of \(\ell ={({\rho }_{{{{{{{{\rm{c}}}}}}}}}C/10\alpha )}^{1/2}\) and \(\lambda =(1/2\alpha ){({\rho }_{{{{{{{{\rm{c}}}}}}}}}C/10)}^{1/2}\), respectively. The stress tensors are scaled by ρc2/2λ2. We set ϵρc = 10.5. To avoid cumbersome expressions, we will use the same characters for the scaled variables in the text. We introduce the normalized density ϕ = ρ/ρc. In our numerical analysis, the density is bound as 0.02 ≤ ϕ ≤ 1.0, by introducing infinite potential walls at ϕ = 0.02 and ϕ = 1.0, to avoid numerical overflow. Time integration of Eqs. (1)-(3) is performed on 256 × 256 square lattices under periodic boundary conditions. The uniform shear flow \(\langle {{\mbox{v}}}{{\mbox{}}}\,({{\mbox{r}}}{{\mbox{}}},t)\rangle =\dot{\gamma }y\,\hat{{{\mbox{x}}}}\) (\(\hat{{{\rm{x}}}}\): a unit vector along x-direction) is applied after t = 0, using the remesh procedure95. At t = 0, v = 0, \(\overleftrightarrow{{{{{{{{\boldsymbol{\sigma }}}}}}}}}=0\), and Gaussian random values of ϕ are assigned at each lattice point77, with the mean value ϕ0 = 〈ϕ〉 = 0.85 and the average variance 0.04. This Gaussian density fluctuation, which represents thermal fluctuation determined by the isothermal compressibility77, is necessary for the density fluctuation to be amplified under deformation in our simulations.

In Eqs. (1)-(3), only the last term on the r.h.s. of Eq. (3) represents the dissipation process, which is due to large scale structural rearrangement with the relaxation time τ(ρ). However, strictly speaking, in Eq. (2) we should include a small intrinsic viscous dissipation term (η02v) due to short time- and small length-scale dynamics. In the present analysis, thus, we add 52v to Eq. (2), taking this effect of the intrinsic viscous dissipation into account. We note that the choice of this term does not affect the main results.