Nominally brittle cracks in inhomogeneous solids: from microstructural disorder to continuum-level scale

We analyze the intermittent dynamics of cracks in heterogeneous brittle materials and the roughness of the resulting fracture surfaces by investigating theoretically and numerically crack propagation in an elastic solid of spatially-distributed toughness. The crack motion split up into discrete jumps, avalanches, displaying scale-free statistical features characterized by universal exponents. Conversely, the ranges of scales are non-universal and the mean avalanche size and duration depend on the loading microstructure and specimen parameters according to scaling laws which are uncovered. The crack surfaces are found to be logarithmically rough. Their selection by the fracture parameters is formulated in term of scaling laws on the structure functions measured on one-dimensional roughness profiles taken parallel and perpendicular to the direction of crack growth.


INTRODUCTION
Understanding how solids break continues to pose significant fundamental challenges. For brittle solids broken under tension, Linear Elastic Fracture Mechanics (LEFM) tackles the difficulty by reducing the problem to the destabilization and subsequent growth of a dominant pre-existing crack (see e.g., [1] for an introduction to LEFM). The theory is based on the idea that, in an elastic material, all dissipative and damaging processes are localized in a small zone around the crack tip, so-called fracture process zone, FPZ. Crack destabilization and further motion are then governed by the balance between the flux of mechanical energy released in the FPZ from the surrounding material and the dissipation rate into this zone. The former is computable within linear elasticity theory and connects to the stress intensity factor, which characterizes the near-tip stress field and depends on the external loading and specimen geometry only. The dissipation rate is quantified by the fracture energy required to expose a new unit area of cracked surfaces. Equivalently, it can be quantified by the fracture toughness, i.e., the onset value for the stress intensity factor above which crack starts to grow. Within LEFM theory, both the fracture energy and the fracture toughness are assumed to be material constants, to be measured experimentally.
LEFM framework provides a coherent framework to describe fracture in homogeneous solids. In contrast, heterogeneous solids remain unclear. Stress concentration at the crack tip makes the behavior observed at the continuum-level scale extremely sensitive to material's small scale inhomogeneities. Consequences include crackling dynamics for fracture with discrete pulses (avalanches) of a variety of sizes and erratic crack paths (see [2] for review). A generic observation in this field is the existence of scale free statistics: (i) Direct imaging of the crack motion in paper sheets [3,4] or along heterogeneous interfaces [5,6] has revealed power-law distribution for the size of the crack jumps; (ii) the acoustic [7][8][9] or seismic [10,11] events going along with fracture are characterized by power-law distribution for the energy; and iii) fracture surfaces are found to exhibit self-affine morphological features [12][13][14][15][16] or logarithmic [17] roughness.
These observations, by essence, cannot be addressed by conventional LEFM continuum approaches. In this context, over the past 25 years, promising alternative approaches have emerged from statistical and non-linear physics: Statistical lattice models like fiber bundle models (see [57] for review) or random fuse models (see [18] for review), for example, were found to reproduce, with a minimal set of ingredients, crackling dynamics in failure [19][20][21] and self-affine fracture roughness [22][23][24], in qualitative agreement with some of the experimental observations. However, these approaches remain too minimal to provide quantitative predictions for situations of engineering interests.
A different approach consists in considering the crack propagation in a solid with spatially-distributed toughness. The availability of asymptotic formulas [25][26][27][28] for the variations of stress intensity factors along a slightly distorted crack can explicitly take into account the microstructural disorder in a continuum-level scale LEFM-like description [29][30][31][32]. Within this framework, referred thereafter to as the Random-Toughness Continuum-Mechanics (RT-CM) approach, the fracture onset can be mapped to a critical depinning transition [2,[33][34][35]. Depending on the situation, this approach yields logarithmically rough [32] or self-affine [36] fracture surfaces. It can also explain the crackling dynamics sometimes observed as a consequence of a selfadjustment of the driving force around the depinning value [37]. The main advantage here is that the external parameters involved www.frontiersin.org November 2014 | Volume 2 | Article 70 | 1 PHYSICS in the depinning model can be mapped to conventional LEFM parameters. This has been used to relate effective toughness and microstructural disorder [38,39], or to unravel the specific conditions required to observed crackling in brittle fracture [40].
In the present work, depinning models unravel how the loading rate, specimen geometry, and microstructural disorder quantitatively select the statistics of continuum-level scale avalanches (as measured in conventional experimental fracture mechanics) and of the crack roughness (as recorded in conventional fractography). Section 2 presents the derivation of the depinning model from LEFM. This section details the successive assumptions and their implications. Section 3 assesses the statistics of the continuum-level scale avalanches, i.e., bursts evidenced in the time evolution of the spatially-averaged crack velocity. Statistics on the avalanche's size and duration demonstrate a power-law characterized by universal exponents. Conversely, the associated cutoffs do depend on the loading, microstructure, and specimen parameters according to scaling laws which are uncovered. Two regimes can be distinguished: A regime of pseudo-isolated avalanches not so different from what is predicted in the quasi-static limit (vanishing loading rate); and a regime where avalanches coalesce with each other. Section 4 analyzes the statistics of fracture roughness. The structure functions measured along profiles parallel and perpendicular to the direction of crack growth exhibit logarithmic scaling with prefactors and characteristic length-scales depending on the Poisson's ratio and microstructure parameters according to scaling which are uncovered.

MATERIALS AND METHODS
In LEFM theory, the crack velocity v is governed by the balance between the mechanical energy G and the (material constant) fracture energy (Griffith criterion). For slow cracks, v/μ = G − where μ = c R / relates the mobility μ to the Rayleigh wave speed c R . The principle of local symmetry (PLS) governs the growth direction [41]. This imposes a zero stress intensity factor condition for mode II (K II = 0) all along the fracture. As consequences, an initially straight crack in an ideal homogeneous material gently loaded in mode I would continuously grow within its plane, without any jerky motion or roughness. But inhomogeneities at the microstructure scale yield distortions of the crack front, which, in turn, induce perturbations δK I , δK II , and δK III in the local loading of the front.
In the RT-CM approach, the effect of material inhomogeneities is taken into account by introducing a random spatiallydistributed component for the fracture energy and for the shear part of the loading. Then, Griffith criterion and PLS combined with an asymptotic estimation of the perturbations in the loading induced by the front distortion describes the crack growth. This approach was pioneered by Gao and Rice [29] and subsequently developed [30-32, 37, 40, 42-44]. Here, the model is reviewed with an emphasis on the relation between the model parameters and experimentally measurable quantities. We also pay special attention to list the various assumptions and discuss the implied limits.

FROM FRONT DISTORTIONS TO THE PERTURBATIONS IN THE CRACK-TIP LOADING
Let one consider a crack embedded in an isotropic elastic solid of size L × H × W under tension. In the following, we adopt the usual convention of fracture mechanics and the axes e x , e y and e z align with the direction of crack propagation (L direction), tensile loading (H direction), and mean crack front (W direction), respectively. For a straight crack, the mode I stress intensity factor K 0 fully characterizes the near-tip stress field loading. Let one now consider the situation depicted in Figure 1: left where the presence of inhomogeneities yields small in-plane and out-of-plane crack distortions, f (z, t) and h(x = f (z, t), z) respectively. The first step in the RT-CM approach is to estimate how these distortions perturb locally the loading. To make the problem tractable, simplifying hypothesis are needed: Hyp. 1 The Young modulus E and Poisson ratio ν of the solid are considered as homogeneous. Hyp. 2 The spatial distribution of (x, y, z) and the PLS noise are narrow enough such that a first-order perturbation analysis can characterize the coupling of the local perturbations δK I (z), δK II (z) and δK III (z) to the front distortions {f , h}.
Then, Rice [25] and Movchan et al's [26] formula relate δK i (z) to the front distortions: where denotes a principal value integral. Note that δK I depends only on the in-plane distortions {f }, and δK II , δK III depend on the out-of-plane distortions {h}. Here, several additional assumptions were required: Hyp. 3 In the expressions of δK II and δK III , we have omitted the terms brought by the non-singular stresses near the reference straight crack front (the so-called T-stress and the third term A of the asymptotic expansion of neartip stress field). This assumption remains valid for stable crack paths (T ≤ 0 [45]) when the roughness wavelength is small with respect to the system size and the characteristic distances defined by the external loading [31,36].

GRIFFITH CRITERION AND EQUATION OF MOTION
The above expression for the local perturbations δK i leads to a first order expression of the energy release rate G(z) as a function of the front distortions: where G 0 (t, f ) = K 0 (t, f ) 2 /E is the reference energy release rate which, as K 0 (t, f ), would result from the same loading with a straight crack front within the mean plane y = 0 at the mean position x = f (t). Henceforth, the operator a indicates averaging of the variable a(z) over the spatial coordinate z. The two quantities G 0 and K 0 coincide with the ones that would have been defined in the conventional LEFM approach when ignoring front distortions. They only depend on the specimen geometry and imposed load (of which both evolve with t) and can be determined using continuum mechanics (e.g., finite element analysis). Slow cracking (as considered herein) implies that the solid is loaded by imposing external displacements rather than external stresses (the opposite would yield dynamic fracture with crack at speed on the order of the sound speed). This induces additional constraints: G 0 should decrease with f (t) (specimen compliance always increases with crack length) and increase with t (to drive a crack, external displacement can only increase with time). Without loss of generality, we set t = 0 at a time when the crack has just stopped and x = 0 is the crack tip's position at this time. G 0 (t = 0, f = 0) = where = (x, y, z) is the mean value of fracture energy. Considering the subsequent variation f (z, t) to be small with respect to the crack length at t = 0, we can write: whereĠ = ∂G 0 /∂t (driving rate) and G = −∂G 0 /∂f (unloading factor) are positive constants characterizing the loading rate and the specimen geometry, respectively. Inserting Equations 1 and 3 in Equation 2, and the resulting expression for G(z) into Griffith criterion, the equation of motion for the in-plane front displacement is: where γ (x, y, z) = (x, y, z) − is the fluctuating part of the fracture energy. The solution of this equation provides the spacetime dynamics of the in-plane projection of the crack front. Subsequently, it gives the time evolution of the fracture velocity at the continuum-level scale, v(t) = df /dt. This is the relevant observable characteristics of the crack dynamics in standard experiments of fracture mechanics.

PRINCIPLE OF LOCAL SYMMETRY AND EQUATION OF PATH
We now derive an equation of path by making use of the principle of local symmetry [41]. To do so, we take inspiration from the work of Katzav et al. [44,47] to model crack path in a model 2D situation and extend it for three-dimensional solids. The idea is to decompose the front propagation along x-axis into infinitesimal straight kinks of length , identified with the microstructural length-scale characterizing the spatial distribution of (x, y, z) (see Figure 1:right). Consider now the kink occurring at the front location z between time t and t + dt. Leblond et al's formulas [48] relate the stress intensity factors after kinking K i (z) (i = I, II, III) to the stress intensity factors before kinking K i (z) and to the kink where F i,j (δθ) are universal functions that were computed in Leblond [27] for three-dimensional solids. The principle of local symmetry implies that K II must be uniformly zero, which leads to: Now, to first order in δθ, F II,I (δθ) and F II,II (δθ) are given by Leblond [27]: F II,I (δθ) = δθ + O(δθ 2 ) and F II,II (δθ) = 1 + O(δθ 2 ). The kink angle is then related to K 0 , δK I (z) and δK II (z) via: δθ ≈ −δK II (z)/(K 0 + δK I (z)). Finally, by keeping only the first order loading perturbations , it writes: www.frontiersin.org November 2014 | Volume 2 | Article 70 | 3 To introduce the effect of microstructural disorder, a spatially distributed noise term K 0 ϑ(x, y, z) is added to the right-handed side of the equation. In the situation of inter-granular crack growth in a material made of sintered grains (e.g. sandstone), for example, this additional term would translate the difference between the kink angle predicted by the principle of local symmetry and that truly selected so that the crack propagates between the cemented grains. Finally, combining the resulting equation with the asymptotic formula for stress intensity factors (Equation 1), one gets: The solution of this equation provides the topography of the postmortem fracture surfaces, h(x, z). This is the quantity of interest in fractography science. Note that this equation differs from that given in Larralde [31] and Bonamy et al. [36] as it includes an additional curvature term ∂ 2 h/∂x 2 .

RELEVANT PARAMETERS AND NUMERICAL ASPECTS
Equations 4 and 8 predict deterministically the fracture dynamics and the morphology of fracture surfaces using the following inputs: Input 1 The loading rateĠ; Input 2 The specimen geometry (unloading factor G and specimen width W) Input 3 The material constants (fracture energy , Poisson ratio ν, and mobility μ), Input 4 The microstructure disorder (microstructure lengthscale and the spatial distribution of the two noise terms γ (x, y, z) and ϑ(x, y, z)).
Statistically, the two noise terms are characterized by the probability density functions p(γ ) and p(ϑ) and the spatial correlation functions C γ ( r) = γ ( r 0 + r)γ ( r 0 ) r 0 and C ϑ ( r) = γ ( r 0 + r)γ ( r 0 ) r 0 . Additional assumptions simplify the problem: Hyp. 5 The two noise terms are not correlated; At first glance, such an assumption may appear odd since both terms originate from the material heterogeneities. But due to the tensorial nature of elasticity, these heterogeneities will affect the equation of motion and the equation of path independently. As an illustrative example, let us consider again the situation of intergranular crack growth in a solid made of sintered grains. Two distinct space-dependent noise terms are required to describe the microstructural texture: The first one quantifies the local variations of adhesion between grains and mainly affects the equation of motion; the second one describes the local variations of joint orientation and affects the equation of path. Hyp. 6 Both probability density functions are Gaussian of zero mean and standard deviationγ andθ, respectively; Hyp. 7 In both cases, the correlation function C decreases linearly with |r| over the microstructure distance and is zero above.
Note that, in elastic interfaces equations as that proposed in Equations 4 and 8, the scaling properties remain unaffected upon changes in the probability functions (p(γ ) and p(ϑ)) and in the correlation functions (C γ (|r|) and C ϑ (|r|)) [43]. Thanks to this series of assumptions, the two spatially distributed noise terms are fully characterized with 3 parameters:γ ,θ and . At this point, the Equations 4 and 8 are coupled via the two 3D noise terms γ (x, y, z) and ϑ(x, y, z), since they both depend on the in-plane and out-of-plane positions of the crack front. A last assumption permits the decoupling of the two equations: Hyp. 8 The 3D spatially distributed terms γ (x, y, z) and ϑ(x, y, z) reduce to their in-plane projection γ (x, z) and ϑ(x, z).
As shown in Ramanathan et al. [32], Equation 8 with a 2D projection for the noise term yields logarithmically smooth fracture surfaces. Hence, zooming out on the fracture surfaces makes them appear flatter and flatter. When the zooming out is sufficient, fluctuations in h become negligible when compared to in-plane length scales. The two noise terms reduce to γ ( . This validates Hyp. 8. At this point, the selection of the fracture behavior brings into play 8 parameters: μ, ,Ġ, G ,γ ,θ, and W. The introduction of dimensionless time t → t/ /(μ ) and length {x, y, z, f , h} → {x/ , y/ , z/ , f / , h/ }) reduces this number. Equations 4 and 8 become: where the dimensionless noise term η is characterized by a standard deviationη =γ / and a spatial correlation length equal to unity. The final two decoupled dimensionless equations require only 6 parameters: • The dimensionless driving rate c =Ġ /μ 2 , • the dimensionless unloading factor k = G / , • the parametersη andθ characterizing the disorder strength, • the continuum-level vs. microstructure scale ratio N = W/ .
The following sections subsequently address the continuum-level scale statistics of crack dynamics v(t) and the morphology of the post-mortem fracture surfaces h(x, z). They also reveal their

CRACKLING DYNAMICS AND AVALANCHE STATISTICS
We first look at the crack dynamics. Equation 9a describes the motion of a long-range elastic chain driven in a frozen random pinning potential with a driving force F = ct − kf (t) selfadjusting around the depinning value [37]. For {c, k} → 0, the chain propagates while remaining at the critical depinning point and the crack moves through irregular jumps, or avalanches. The size S and duration D of the avalanche 1 follow a universal powerlaw distribution [37,49]: P(S) ∝ S −τ and P(D) ∝ D −α where the exponents τ and α can be estimated using renormalization [49][50][51] or numerical [52,53] methods. The most precise estimations yields τ = 1.280 ± 0.01 and α = 1.500 ± 0.010 [2]. These scale free statistics extend to upper cutoffs S 0 and D 0 set by the system 1 In problems of depinning of interfaces, the definition of the avalanche size differ depending on the publication. Sometimes, the size is defined as the area A swept by the front between two successive depinning configurations. Here, the size S is defined as the integral of the continuum-level scale velocity v(t) between the two successive depinning configurations. The two definitions are related by S = A/N when N is the system size.
size N: S 0 ∝ N ζ and D 0 ∝ N κ where ζ and κ are the roughness exponent and dynamic exponent, respectively. The most precise estimations of their value yield [52,53]: ζ = 0.385 ± 005 and κ = 0.770 ± 0.005. Here, both c and k are finite. The effects of a finite unloading factor k keeping c → 0 is now fairly well documented [2,54]: It modifies the upper cutoffs S 0 and D 0 to the preceding scale-free distributions: Here, the form of the functions f (u) and g(u) is expected to be universal. The exponents 1/σ and 1/ are also predicted to be universal: 1/σ = 0.69 ± 0.010 and 1/ = 0.385 ± 0.010 [2,54]. Conversely, the effect of a finite c is not uncovered yet. By yielding some overlap between the avalanches [55], it can significantly alter the dynamics. In particular, a recent work [40] has evidenced a transition line between a crackling-like dynamics made of irregular power-law distributed jumps and a continuum-like dynamics ruled by the conventional LEFM theory. Figure 2 presents typical times profiles of the continuum-level scale velocity v(t) for different values of c and k in the crackling regime. Here and thereafter, the system size N and the disorder strengthη are set constant, N = 1024 andη = 1, returning at the end of this section to a brief discussion on their effect. Note the irregular jumps characteristics of the underlying avalanching dynamics. Note also the qualitative changes in the signal appearance as k and c are modulated: Pulses become shorter as k increases, and the pulse density increases as c increases. Note finally that, due to the finite value of c, v(t) does not vanish

www.frontiersin.org
November 2014 | Volume 2 | Article 70 | 5 between the pulse, but becomes equal to a small value proportional to c (prefactor dependent of the Runge Kunta scheme). The avalanches are then identified with the bursts where v(t) is above a prescribed reference level v c = 10 −3 . Their duration D is defined as the interval between the two successive intersections of v(t) with v c , and their size S is defined as the integral of v(t) − v c between the same points. Figure 3 reports the probability density function of avalanche size P(S|c, k) and duration P (D|c, k) for different values of c and k. Power-laws are observed. The exponents τ and α associated with the power-law decrease are independent of c and k. Moreover, they compare well with the values τ = 1.280 ± 0.01 and α = 1.500 ± 0.010 expected for {c, k} → 0. On the other hand, the valid region of the power-law is observed depends on c and k. Both the lower cutoffs S min and D min are roughly independent of c and k. However, the precise values S min and D min were observed to depend on the reference level v c and, thus, are to be associated with the procedure to extract the avalanches. We will consider in the following only the part of the distributions above these lower cutoffs.
The power-law distributions observed in Figure 3 also exhibit upper cutoffs S 0 and D 0 . These cutoffs decrease with k in all cases. The effects of c is of two types: The upper-cutoff is found to decrease with k and to increase with c (resp. to be independent of c) when c is large enough (resp. when c is small enough).

Frontiers in Physics | Interdisciplinary Physics
November 2014 | Volume 2 | Article 70 | 6 ( Figure 4A) and D = ∞ D min D × P(D|c, k)dD (Figure 4B). At large enough k, both S and D are independent of c. This large k regime is attributed to a regime of pseudo-isolated (pi) avalanches. The distributions are then expected to take forms similar to that of Equation 10. As a result, S pi is expected [56] to take the form S pi ≈ S τ −1 These two scaling are compatible with the observations at large k (dash line in Figures 4A,B). As a synthesis, the mean avalanche size and duration are found to take the following form at large k: Let us try now to characterize the effects of the avalanche overlap when k becomes small or c becomes large. Previous work [40] evidenced a transition between the crackling dynamics studied here and a continuum-like dynamics when c becomes large enough or k small enough. This transition is believed [40] to coincide with the point where the avalanche overlap percolates throughout the entire system. At constantη and N, this transition was shown [40] to be fully driven by the ratio c/k. We hence plotted, in Figure 5, S / S pi and D / D pi (first order estimation of the number of individual avalanches having merged together to form the bursts detected from the signal v(t)), as a function of the control parameter c/k. A coarse collapse is observed. As expected, the master curves diverge at the transition value between crackling and continuum-like dynamics (materialized by the vertical dash lines in the main panels of Figures 5A,B). Note that significant deviations to the collapse are observed. They are believed to stem from a qualitative change in the distribution shape as k decreases/c increases and the avalanche overlap increases. Actually, for small k/large c a bump develops at large scales, for both the distribution in size and in duration. Due to this change in shape, several length-scales (resp. time-scales) intervenes in the distributions of S (resp. D). Thus, the recording of the mean values alone is not sufficient to capture their evolutions with c and k. Ongoing work aims at characterizing these effects more accurately.
To make the analysis complete, we looked at the effects of the system size N and disorder strengthη. Figure 6 plots the probability density function of avalanche size P(S|N,η) and avalanche duration P (D|N,η) as measured for different values of N andη at fixed values of c and k. For both size and duration, the decrease of N yields an increase of the lower-cutoffs S min and D min (main panel of Figures 6A,A ). The two dependencies are well fitted by power-laws: S min ∝ N −a SN and D min ∝ N −a DN with a SN ≈ 1.7 and a DN ≈ 0.6 (inset of Figures 6A,A ).η does not seems to affect S min . Conversely, D min decreases withη as D min ∝η −a Dη with a Dη ≈ 1.2 (inset of Figure 6B ). The effects of N andη of S and D are analyzed in Figure 7, by plotting the curves S vs c and D vs c for different N at constant k andη (insets in Figures 7A,A ), and for differentη at constant k and N (insets in Figures 7B,B ). Increasing N yields a decrease of the low-c plateau and a leftward shift of the divergence location for both S and D . Increasingη yields an increase of the low-c plateau for S , a decrease of the low-c plateau for D , and a rightward shift of the divergence location for both S and D . All curves can then be superimposed by making Figure 7B . By combining this with Equation 11 and the collapse obtained in Figure 5, we anticipate the following form for mean size and duration:  graphs (A,B). In the main panels of both graphs, the vertical dash lines indicates the transition value c/k between crackling and continuum-like dynamics as determined in Nukala et al. [40] for N = 1024 andη = 1.
where the two functions F(u) and G(u) exhibit a plateau at small u, and both diverge at the same value u c . The value of the exponents τ , α, 1/σ and 1/ are well known [2]: They can e.g., be related to the so-called roughness exponent ζ and dynamic exponent κ classically defined in the realm of critical depinning transition: Conversely, the precise origin of the exponents b SN , b Sη , b DN , b Dη , c N , and cη and their link with ζ and κ remain to be uncovered. Equation 12 quantitatively relate the material parameters to quantities that are accessible in conventional experimental mechanics, namely the mean size and duration of the avalanches. In this context, it is of interest to rewrites the equation with the original variables, before the non-dimensionalization procedure: where the mean avalanche size S and duration D are expressed with real length and time units, andĠ, G , W, , μ, andγ are recalled to be the loading rate, the unloading factor, the specimen width, the fracture energy, the mobility, the microstructure length-scale, and the contrast in local fracture energy. The value of the exponents are recalled to be τ = 1.280 ± 0.01 (predicted), α = 1.500 ± 0.01 (predicted), 1/σ = 0.69 ± 0.01 (predicted),

ROUGHNESS OF FRACTURE SURFACES
We turn now to the topography h(x, z) of the post-mortem fracture surfaces as predicted by equation 9b. Figure 8 reports typical topographies for different values of the two external parameters A andθ. When A is close to 1, the surface seems to be statistically isotropic while as A gets smaller, the surface appears more elongated in the direction of z. Conversely, the parameterθ only affects the range swept by the roughness. Note that in almost all elastic solids, Poisson ratio ν lies between 0 and 0.5, which impose a finite interval for A = (2 − 3ν)/(2 − ν): 1/3 ≤ A ≤ 1. Herein, only A within this interval are considered.
To characterize quantitatively the spatial distribution of fracture roughness, we adopted the classical procedure [2] and computed the structure function S( r) = (h( r + r) − h( r)) 2 . Here, the operator denotes averaging over all positions r = (x, z). First, we computed the structure function S z ( z) along 1D profiles taken parallel to z (mean direction of the crack front). The procedure is the following: (i) an initially straight front was first propagated over a distance equal to 10N to obtain a statistically stationary regime; (ii) 10000 subsequent profiles h(x i , z) separated by a distance x i+1 − x i = 1 were recorded; (iii) the structure function S i z ( z) was computed for each of these profiles; and (iv) finally, these 10000 individual structure functions were averaged to get S z ( z). Figure 9 depicts typical examples of S z ( z) curves for different values of N, θ and A. S z goes as S z = p z log ( z/λ z ) up to an upper cutoff set by the system size N. This logarithmic scaling is anticipated to extend over the whole range of length-scales as N → ∞. Note that logarithmically rough crack surfaces were also predicted in earlier theoretical works [31,32] analyzing crack propagation through a three-dimensional heterogeneous solid. As a plus, the present model allows relating the prefactor p z and characteristic length-scale λ z with the fracture parameters: p z is found to scale asθ 2 /A (Figure 10A), while λ z is independent of both A andθ ( Figure 10B). Finally, the structure function along z is: log ( z/λ z ), with C = 0.32 ± 0.01 and λ z = 0.24 ± 0.03 (14) We now look at the structure function S x ( x) along x. A direct computation of S x following the standard procedure proposed for S z was found to give a large scattering even for a constant set of parameters {N, A,θ}. Hence, the computation procedure was modified as follows: (i) an initially straight front was first propagated over a distance equal to 10N to obtain a statistically stationary regime; (ii) the x evolution ofh(x, z i ) = h(x, z i ) − h(x) was subsequently recorded at 100 locations z i uniformly distributed along the specimen width N (h(x) denotes averaging over the specimen width N); (iii) the structure function 2 was computed for each of these profiles; (iv) these 100 individual structure functions were averaged to get the S x ( x) for a single specimen; and (v) the so-obtained structure functions were further averaged over 100 specimens. This procedure produces accurate and reproducible curves S x vs. x. Figure 11 depicts typical examples of S x ( x) curves for different values of N, θ and A. The behavior resembles that observed for S z , with a logarithmic scaling S x = p x log ( x/λ x ) up to an upper cutoff set by the system size N. Note that S x saturates above the cutoff, and does not decrease down to zero as was observed for S z due to periodic boundary conditions. As for S z , the prefactor (now referred to as p x ) goes asθ 2  ( Figure 12A). The characteristic length-scale (λ x ) is independent ofθ (Figure 12B:inset). But contrary to what is observed for S z , this characteristic length λ x depends on A: λ x ∝ 1/A. This dependency is responsible for the apparent stretching along s of the images in Figure 10 observed as A decreases. Finally, the structure function along x is:  λ z associated with the curve S z vs. z as a function of A at constantθ = 1 (main) and as a function ofθ at constant A = 1 (inset). In both graphs, the red lines correspond to fits λ z = 0.24 ± 0.03. Here, ± indicates a 95% confident interval.
where S x , S z , x and z are expressed with real length units, and ν, andθ are recalled to be the Poisson ratio, the microstructure scale, and disorder contrast.

CONCLUDING DISCUSSION
Stress enhancement at crack tips makes the macroscale failure behavior observed extremely sensitive to the presence of disorder at the microstructure scale. This translates into crackling dynamics and rough fracture surfaces, which, by essence, cannot be addressed within the conventional LEFM framework. In this paper, we have used the RT-CM approach to obtain quantitative relations between some statistical observables characteristic of these two aspects and the fracture parameters: Loading rate (time derivative of the energy release rate), specimen geometry (specimen thickness and unloading factor), conventional mechanical constants (fracture energy, Poisson ratio), and microstructural disorder (microstructure scale and disorder strength).
Over a certain range of the fracture parameters, this RT-CM approach predicts crackling dynamics [40]: The crack growth splits up into discrete jumps, which are power-law distributed in size and duration. The characteristic exponents associated to these power-laws are universal. Conversely, the scales covered by these scale-free features are non-universal and, in particular, the mean size and duration of the crack jumps are found to depend on the fracture parameters according to scaling laws that are uncovered. These scaling laws can be understood over a certain range of the fracture parameters, in the regime of pseudo-isolated avalanches addressable via standard functional renormalization theory [49,[51][52][53]. Conversely, the effect of the avalanche overlapping is not understood. On-going work aims at analyzing the distribution of the local avalanches as detected in the space-time diagrams of the front dynamics, in order to understand the coalescence process.
Also, this RT-CM approach predicts rough fracture surfaces. The fracture roughness can be characterized by computing the structure function, which exhibits logarithmic scaling. The associated prefactor and characteristic length-scale are found to depend on the Poisson ratio, microstructure www.frontiersin.org November 2014 | Volume 2 | Article 70 | 11

FIGURE 12 | (A)
Slope p x associated with the curve S x vs. x as a function of A at constantθ = 1 (main) andθ at constant A = 1 (inset). In the inset, the axes are logarithmic. In both graphs, the red lines are fits p 1 = C x /A (main) and p 1 = Cθ 2 where the fitted parameter is found to be C = 0.31 ± 0.02 (95% confident interval). (B) Characteristic length-scale λ x associated with the curve S x vs. x as a function of A at constantθ = 1 (main) andθ at constant A = 1 (inset). In both graphs, the red lines are fits λ x = 0.21 ± 0.02/A (main) and λ x = 0.21 ± 0.02 (inset). Here, ± indicates a 95% confident interval.
length-scale, and disorder strength according to laws that were uncovered. This may have interesting applications: It allows one to infer the microstructure parameters (the access of which could be made difficult otherwise, due to the smallness of the length-scales involved) from the analysis of postmortem fracture surfaces at larger scale. Work in progress aims at testing the scaling predicted here for the structure functions against fractography experiments achieved in oxide glasses. Note finally that the RT-CM model studied here call upon a variety of assumptions (see Section 2). An interesting perspective would be to measure to which extend these assumptions can be released. Work in this direction is currently under progress. The model is also limited to nominally brittle fracture, with a single macroscopic crack propagating in an otherwise intact material. Promising alternative approaches have emerged from statistical and non-linear physics [18,57] and may permit to address quasi-brittle fracture, with many microcracking events interacting with each-others. A central challenge in the field would be to bridge the gap between these approaches and engineering damage mechanics.