Mathias Lebihain
email: mathias.lebihain@enpc.fr

Manish Vasoya

Véronique Lazarus

Size effects in the toughening of brittle materials by heterogeneities: a non-linear analysis of front deformations

Keywords: Brittle failure, three-dimensional fracture, stress intensity factor, coplanar perturbation, circular crack, effective toughness

published or not. The documents may come

Introduction

Despite the rise of powerful numerical methods, like e.g. cohesive zone models (CZMs) [START_REF] Barenblatt | The mathematical theory of equilibrium cracks in brittle fracture[END_REF][START_REF] Xu | Numerical simulations of fast crack growth in brittle solids[END_REF][START_REF] Camacho | Computational modelling of impact damage in brittle materials[END_REF] or phase-field models (PFMs) [START_REF] Francfort | Revisiting brittle fracture as an energy minimization problem[END_REF][START_REF] Hakim | Crack path prediction in anisotropic brittle materials[END_REF][START_REF] Bourdin | The variational approach to fracture[END_REF] that can compute spontaneous three-dimensional crack initiation and propagation in complex structures and multiphysics environments, more "traditional" simulation techniques based on the Linear Elastic Fracture Mechanics (LEFM) like e.g. the eXtended Finite Element Method (XFEM) [START_REF] Moës | A finite element method for crack growth without remeshing[END_REF], have still some value in modeling crack propagation in elaborate settings [START_REF] Paul | 3d coupled HM-XFEM modeling with cohesive zone model and applications to non planar hydraulic fracture propagation and multiple hydraulic fractures interference[END_REF][START_REF] Lebihain | Effective toughness of disordered brittle solids: A homogenization framework[END_REF]. This is particularly true when dealing with sharp discontinuities of material properties, as the characteristic length scale associated with the failure process in CZMs and PFMs may interact with that of the fluctuations. In that case, numerical simulations may fail in reproducing the predictions of LEFM, even in simple situations of crack deflection at a bi-material interface [START_REF] Henry | Limitations of the modelling of crack propagating through heterogeneous material using a phase field approach[END_REF], and ah-hoc numerical compensation methods are required [START_REF] Hansen-Dörr | Phase-field modeling of crack branching and deflection in heterogeneous media[END_REF]. Meanwhile, LEFM-based methods, which build on the small-scale yielding assumption, do not suffer these limitations due to the innate localized nature of the dissipation [START_REF] Lebihain | Effective toughness of periodic heterogeneous materials: the effect of out-of-plane excursions of cracks[END_REF]. Yet, in these models, crack extension is not spontaneous, but must be predicted from the combination of ad-hoc propagation [START_REF] Griffith | The phenomena of rupture and flow in solids[END_REF][START_REF] Irwin | Fracture[END_REF] and direction [START_REF] Erdogan | On the crack extension in plates under plane loading and transverse shear[END_REF][START_REF] Hussain | Strain energy release rate for a crack under combined mode i and mode II[END_REF][START_REF] Gol'dstein | Brittle fracture of solids with arbitrary cracks[END_REF] criteria. They often require the a priori knowledge of the stress intensity factors (SIFs) along the whole crack front. Stress intensity factors are known analytically for a limited number of simple crack geometries and loading conditions [START_REF] Tada | The Stress Analysis of Cracks Handbook[END_REF], and they can be computed numerically from e.g. the J-integral [START_REF] Rice | A path independent integral and the approximate analysis of strain concentration by notches and cracks[END_REF], or the virtual crack extension (also called G-theta) method [START_REF] Destuynder | Quelques remarques sur la mécanique de la rupture élastique[END_REF][START_REF] Delorenzi | On the energy release rate and the j-integral for 3-d crack configurations[END_REF]. However, when a crack propagates in a 3D composite material, its front is locally pinned by material heterogeneities, and adopts a tortuous shape that departs from "'standard" geometries [START_REF] Gao | A first-order perturbation analysis of crack trapping by arrays of obstacles[END_REF]. As a result, the SIF distribution along the front cannot be inferred from exact analytical solutions, and the computational cost of SIF evaluation often proves too high to model 3D propagation from the scale of the smallest heterogeneity to that of the structure.

A good compromise can be found by resorting to perturbative approaches of LEFM [START_REF] Lazarus | Perturbation approaches of a planar crack in linear elastic fracture mechanics: A review[END_REF], in which the stress intensity factors along a front that is slightly distorted within its plane from a reference configuration can be estimated from the knowledge of that in the unperturbed configuration. Building on Bueckner-Rice's weight function theory [START_REF] Bueckner | Weight functions and fundamental fields for the penny-shaped and the half-plane crack in three-space[END_REF][START_REF] Rice | Weight function theory for three-dimensional elastic crack analysis[END_REF], [START_REF] Rice | First-order variation in elastic fields due to variation in location of a planar crack front[END_REF] derived the mode I stress intensity factor distribution along a semi-infinite planar crack whose front is slightly perturbed from its reference straight configuration, up to first order in the perturbation amplitude. The approach was then extended to (i) mixed mode conditions by [START_REF] Gao | Shear stress intensity factors for a planar crack with slightly curved front[END_REF], (ii) other crack geometries, like e.g. the penny-shaped crack (Gao and Rice, 1987b;[START_REF] Borodachev | Variational formula for a body with a crack and its applications (review)[END_REF], the circular connection (Gao and Rice, 1987a), tunnel crack [START_REF] Leblond | The tensile tunnel-crack with a slightly wavy front[END_REF]Lazarus and Leblond, 2002a,b), interaction of two tunnel cracks [START_REF] Pindra | In-plane perturbation of a system of two coplanar slit-cracks -i: Case of arbitrarily spaced crack fronts[END_REF][START_REF] Legrand | In-plane perturbation of a system of two coplanar slit-cracks -II: Case of close inner crack fronts or distant outer ones[END_REF], the crack lying between two plates [START_REF] Legrand | Coplanar perturbation of a crack lying on the midplane of a plate[END_REF], (iii) out-of-plane perturbations [START_REF] Moës | A finite element method for crack growth without remeshing[END_REF][START_REF] Leblond | Theoretical analysis of crack front instability in mode i+III[END_REF], (iv) dynamic ruptures [START_REF] Perrin | Disordering of a dynamic planar crack front in a model elastic medium of randomly variable toughness[END_REF][START_REF] Ramanathan | Dynamics and instabilities of planar tensile cracks in heterogeneous media[END_REF], and (v) cohesive materials [START_REF] Lebihain | Quasi-static crack front deformations in cohesive materials[END_REF]. The unparalleled computational performances offered by the perturbative methods, through the use of the Fast Fourier Transform (FFT), allowed to investigate the propagation of a crack front in heterogeneous media [START_REF] Gao | A first-order perturbation analysis of crack trapping by arrays of obstacles[END_REF], its roughening under the action of material disorder [START_REF] Schmittbuhl | Interfacial crack pinning: Effect of nonlocal interactions[END_REF][START_REF] Rosso | Roughness at the depinning threshold for a long-range elastic string[END_REF], and the intermittent statistics that emerge from its interaction with material heterogeneities [START_REF] Ponson | Crack propagation in brittle heterogeneous solids: Material disorder and crack dynamics[END_REF][START_REF] Laurson | Avalanches and clusters in planar crack front propagation[END_REF]. It has given quantitative means to rationalize the crack-inclusion interactions observed during failure experiments in patterned [START_REF] Bower | A three-dimensional analysis of crack trapping and bridging by tough particles[END_REF][START_REF] Dalmas | Crack front pinning by design in planar heterogeneous interfaces[END_REF][START_REF] Chopin | Crack front dynamics across a single heterogeneity[END_REF]Patinet et al., 2013a) and disordered materials [START_REF] Delaplace | High resolution description of a crack front in a heterogeneous plexiglas block[END_REF][START_REF] Barés | Aftershock sequences and seismic-like organization of acoustic events produced by a single propagating crack[END_REF]. It has also provided basic ingredients to formulate a homogenization framework for coplanar [START_REF] Roux | Effective toughness of heterogeneous brittle materials[END_REF]Patinet et al., 2013b;Démery et al., 2014a,b) and non-coplanar [START_REF] Lebihain | Effective toughness of disordered brittle solids: A homogenization framework[END_REF] brittle failure. Note that perturbative approaches have also been used in modeling contact and adhesion along homogeneous and heterogeneous surfaces [START_REF] Adda-Bedia | Crack-front instability in a confined elastic film[END_REF][START_REF] Xia | Toughening and asymmetry in peeling of heterogeneous adhesives[END_REF][START_REF] Budzik | Perturbation analysis of crack front in simple cantilever plate peeling experiment[END_REF][START_REF] Xia | Adhesion of heterogeneous thin films -ii: Adhesive heterogeneity[END_REF][START_REF] Argatov | Controlling the adhesive pull-off force via the change of contact geometry[END_REF][START_REF] Sanner | Crack-front model for adhesion of soft elastic spheres with chemical heterogeneity[END_REF].

However, the experimental study of Vasoya et al. (2016b) suggests that higher-order theories may help to reconcile theoretical predictions with experimental observations of crack front deformation. Moreover, material toughening by heterogeneities has been shown to be of second-order in the amplitude of the toughness fluctuations (Patinet et al., 2013b;Démery et al., 2014b), which is of the same order as the terms omitted by these authors in their asymptotic expansion of the stress intensity factor. As a result, their models only draw conclusions on how material toughening scales with the amplitude of the fluctuations, but they may fail in predicting it quantitatively. Higher-order theories for the crack front perturbations are thus crucial to improve our understanding of the interactions between cracks and heterogeneities. Yet, up to this day, second-order perturbative approaches have only been derived for the semi-infinite crack (Adda-Bedia et al., 2006;[START_REF] Leblond | Second-order coplanar perturbation of a semi-infinite crack in an infinite body[END_REF][START_REF] Vasoya | A geometrically nonlinear analysis of coplanar crack propagation in some heterogeneous medium[END_REF]Vasoya et al., , 2016b) ) loaded in tensile mode I. Their extension to other geometries (like e.g. the penny-shaped crack or the circular connection) may bring valuable insights on a wide variety of physical problems, like e.g. the onset of the fingering instability occurring during crack propagation (Vasoya et al., 2016a) or during the contact between an indenter and soft elastic films [START_REF] Yu | Adhesion-induced instability regulates contact mechanics of soft thin elastic films[END_REF], the adhesion hysteresis observed along chemically heterogeneous interfaces [START_REF] Sanner | Crack-front model for adhesion of soft elastic spheres with chemical heterogeneity[END_REF], or even the shape of fluid-driven shear ruptures along frictional interfaces governed by Coulomb's friction [START_REF] Sáez | Three-dimensional fluid-driven stable frictional ruptures[END_REF].

In this work, we extend the linear model of Gao and Rice (1987b) for the penny-shaped crack loaded in tensile mode I to second-order, following the reasoning of [START_REF] Leblond | Second-order coplanar perturbation of a semi-infinite crack in an infinite body[END_REF]. We build next on our non-linear theory to investigate material reinforcement arising from the presence of randomly arranged heterogeneities of toughness, and focus in particular on the size effects emerging from the finiteness of the penny-shaped crack geometry.

The paper is organized as follows: in Section 2, we derive a non-linear theory that predicts the mode I stress intensity factor variations arising from coplanar perturbations of the crack front from its reference circular configuration, up to second-order in the perturbation amplitude. This is performed by deriving first the linear expansion of the fun-damental kernel associated to the penny-shaped crack geometry. In Section 3, we validate our theory by comparing its predictions with analytical solutions obtained for elliptical cracks loaded by a uniform tensile stress. The validity range of the second-order model is then assessed, building on numerical SIF estimates along wavy crack fronts obtained by the method of [START_REF] Lazarus | Brittle fracture and fatigue propagation paths of 3d plane cracks under uniform remote tensile loading[END_REF] and [START_REF] David | On the key role of crack surface area on the lifetime of arbitrarily shaped flat cracks[END_REF] that computes the SIF variations at all orders. Section 4 is devoted to the analysis of quasi-static crack propagation in random fields of toughness, with a particular focus on the toughening arising from these fluctuations. Namely, we discuss how (i) accounting for non-linear terms in the SIF expansion and (ii) the finiteness of the penny-shaped crack geometry may change our results with respect to those available in the literature for the semi-infinite crack (Patinet et al., 2013b;Démery et al., 2014a).

Derivation of the second-order theory

Generalities

We consider a planar penny-shaped crack embedded in an infinite medium made of some isotropically linearly elastic material of Young's modulus E and Poisson's ratio ν. The crack front F 0 describes a circle of radius a 0 centered on a point O. The crack is contained within the xOz plane. The pair (r, θ) denotes the coordinates of a point of observation M in the polar coordinate system centered in O, and (e r , e θ ) the polar coordinate basis vectors at this point. The crack is loaded in pure mode I through some system of forces, giving rise to a distribution K 0 I (a 0 , θ) of the unperturbed stress intensity factor along the reference circular crack front F 0 .

The crack front is now perturbed within its plane in the direction e r normal to its reference configuration by a small amount:

δa(θ) = Aϕ(θ), (1)
where A is an infinitesimally small parameter quantifying the non-circularity of the perturbed front F , and ϕ is a shape function (see Fig. 1). δK I denotes the variations of the mode I SIF arising from the perturbation δa. Evaluated at a point M 1 of F indexed by θ 1 , δK I reads [START_REF] Panasyuk | Propagation of a crack whose shape is close to circular[END_REF]Gao and Rice, 1987b):

δK I (θ 1 ) = ∂K 0 I ∂a (a 0 , θ 1 )δa(θ 1 ) + PV 2π 0 Z 0 (a 0 ; θ 1 , θ)K 0 I (a 0 , θ) [δa(θ) -δa(θ 1 )] a 0 dθ + O(A 2 ), (2)
at first order in the perturbation δa. Here, the symbol PV denotes a Cauchy principal value, and Z 0 (a 0 ; θ 1 , θ 2 ) is the kernel associated with the reference circular crack of radius a 0 , evaluated at any points M 0 1 and M 0 2 of F 0 , indexed by θ 1 and θ 2 respectively (see Fig. 1). It takes the general form [START_REF] Rice | Weight function theory for three-dimensional elastic crack analysis[END_REF]:

Z 0 (a 0 ; θ 1 , θ 2 ) = 1 2π W 0 (a 0 ; θ 1 , θ 2 ) D(M 0 1 , M 0 2 ) 2 , (3)
where D(M 0 1 , M 0 2 ) is the Euclidean distance between the points M 0 1 and M 0 2 , and W 0 is the fundamental kernel1 associated with the reference crack front geometry. For a circular crack front, W 0 (a 0 ; θ 1 , θ 2 ) = 1, so that Z 0 reads:

Z 0 (a 0 ; θ 1 , θ 2 ) = 1 8πa 2 0 sin 2 [(θ 2 -θ 1 )/2] . ( 4
)
Our goal here is to go beyond the linear expansion (2) of Gao and Rice (1987b), and derive the expression of the perturbed SIF K I at second order in the perturbation δa. We follow here the approach of [START_REF] Leblond | Second-order coplanar perturbation of a semi-infinite crack in an infinite body[END_REF], and calculate first the variations of the kernels W and Z at first order in δa.

x

a 0 θ 2 θ 1 a ( θ ) O y M 1 M 2 θ M 0 2 M 0 1 Cracked Unbroken Reference crack front ℱ 0 Perturbed crack front ℱ z M 0 M Figure 1:
The reference penny-shaped crack with a circular front F 0 of radius a 0 (in dashed black line) is perturbed within its plane (xOz) by a small amount δa(θ) in the direction e r , giving rise to a perturbed crack front F parametrized by a(θ) = a 0 + δa(θ) (in solid red line). The points M 1 , M 2 , and M, indexed by θ 1 , θ 2 and θ respectively, are located on the perturbed front F . M 0 1 , M 0 2 and M 0 are their orthogonal projections on F 0 .

First-order expansion of the fundamental kernel

For tensile cracks, the first-order expansion of the fundamental kernel W has been derived by [START_REF] Rice | Weight function theory for three-dimensional elastic crack analysis[END_REF]. Provided that the perturbation δa vanishes at M 1 and M 2 indexed by θ 1 and θ 2 , the variations δW of the fundamental kernel write as:

δW(θ 1 , θ 2 ) D(M 1 , M 2 ) 2 = 1 2π PV 2π 0 W 0 (a 0 ; θ 1 , θ) D(M 1 , M 0 ) 2 W 0 (a 0 ; θ 2 , θ) D(M 2 , M 0 ) 2 δa(θ) a 0 dθ + O(A 2 ), (5)
where M 0 is the point indexed by θ on F 0 .

In the general case where δa does not vanish at M 1 and M 2 , one must find a perturbation δa * (i) for which the variations of the fundamental kernel δW * are known a priori, and (ii) that takes the same values as δa at two points M 1 and M 2 (see Fig. 2):

δa * (θ 1 ) = δa(θ 1 ), and δa * (θ 2 ) = δa(θ 2 ). ( 6)

One then has (Favier et al., 2006a):

δW(θ 1 , θ 2 ) D(M 0 1 , M 0 2 ) 2 = δW * (θ 1 , θ 2 ) D(M 0 1 , M 0 2 ) 2 + 1 2π PV 2π 0 W 0 (a 0 ; θ 1 , θ) D(M 0 1 , M 0 ) 2 W 0 (a 0 ; θ 2 , θ) D(M 0 2 , M 0 ) 2 [δa(θ) -δa * (θ)] a 0 dθ + O(A 2 ), (7)
where M 0 1 and M 0 2 are the locations defined by the orthogonal projection on F 0 along the direction e r of the points M 1 and M 2 , i.e. their position before the application of the transformation δa (see Fig. 1). The variation of the kernel Z then reads:

δZ(θ 1 , θ 2 ) = W 0 (a 0 ; θ 1 , θ 2 ) + δW(θ 1 , θ 2 ) 2πD(M 1 , M 2 ) 2 - W 0 (a 0 ; θ 1 , θ 2 ) 2πD(M 0 1 , M 0 2 ) 2 = Z 0 (a 0 ; θ 1 , θ 2 )              D(M 1 , M 2 ) D(M 0 1 , M 0 2 )       -2 -1        + δW * (θ 1 , θ 2 ) 2πD(M 0 1 , M 0 2 ) 2 + PV 2π 0 Z 0 (a 0 ; θ 1 , θ)Z 0 (a 0 ; θ 2 , θ) [δa(θ) -δa * (θ)] a 0 dθ + O(A 2 ), (8)
at first order in A. In the case of the semi-infinite crack, the first term of δZ is of second order in A, and reduces thus to zero in the first-order calculations of [START_REF] Leblond | Second-order coplanar perturbation of a semi-infinite crack in an infinite body[END_REF]. This is not the case for a reference circular crack front,

z x M 1 M 2 θ 2 θ 1 Cracked O* θ a 0 + Δ a * δa * (θ) Reference crack front ℱ 0 Perturbed crack front ℱ Shifted crack front ℱ* y Unbroken Figure 2:
The first-order variations of the kernel Z 0 (a 0 ; θ 1 , θ 2 ), associated with the reference circular crack front (in dashed black line), are calculated from a shifted circular crack front F * (in dashed blue line) of radius a 0 + ∆a * and centered in O * , that passes by M 1 and M 2 located on the perturbed crack front (in red solid line). The transformation from the reference crack front F 0 and shifted crack front F * is parametrized by the perturbation δa * (θ) defined in Eq. ( 10).

for which it writes as :

Z 0 (a 0 ; θ 1 , θ 2 )              D(M 1 , M 2 ) D(M 0 1 , M 0 2 )       -2 -1        = - δa(θ 1 ) + δa(θ 2 ) 8πa 3 0 sin 2 [(θ 2 -θ 1 )/2] + O(A 2 ). ( 9
)
The second and third terms of δZ in Eq. ( 8) depend on the actual choice of the perturbation δa * . To compute the second term analytically, the fundamental kernel of the crack front parametrized by a * (θ) = a 0 + δa * (θ) must be known explicitly. For unbounded solids, [START_REF] Rice | Weight function theory for three-dimensional elastic crack analysis[END_REF] proposed to express δa * as the combination of translations, rotations or scaling, leaving the fundamental kernel unchanged, i.e. δW * = 0. Such a transformation has been used by [START_REF] Leblond | Second-order coplanar perturbation of a semi-infinite crack in an infinite body[END_REF] to derive the first-order variation of the fundamental kernel W for the semi-infinite coplanar crack loaded in tensile mode I. However, one can show that, in the case investigated here, this transformation contains a rotation so that the perturbation δa * cannot be directly expressed as a normal extension in the direction e r .

Here, we follow another route and take δa * as the orthogonal projection in the direction e r , on the reference crack front F 0 , of another circular crack front F * that goes by M 1 and M 2 (see Fig. 2). F * corresponds to F 0 dilated by a factor (1 + ∆a * /a 0 ), and translated by an amount ∆r * in the direction π/2 + θ * . The perturbed crack front F * is thus a circle of radius a 0 + ∆a * , centered on the point O * of Cartesian coordinates (-∆r * sin θ * , ∆r * cos θ * ). The associated perturbation δa * can be expressed as:

δa * (θ) = ∆r * sin(θ -θ * ) + (a 0 + ∆a * ) 2 -∆r * 2 cos 2 (θ -θ * ) -a 0 . ( 10
)
Since the crack front F * is circular, its fundamental kernel W * is also equal to 1. As a result, δW * = 0, and the second term of δZ in Eq. ( 8) vanishes. It leaves us with the calculation of the principal value integral in Eq. ( 8). There is an infinite number of circles that go by the points M 1 and M 2 . Consequently, δa * is a 3-parameters transformation (∆a * , ∆r * , θ * ) that must satisfy the sole 2 conditions of Eq. ( 6). One may then carefully choose an appropriate θ * so that (i) δa * of Eq. ( 10) is well-defined, and (ii) one can compute the third term of δZ analytically. Here, we take:

θ * = θ 1 + θ 2 2 . ( 11
)
Note that, in equation ( 8), one only needs an expression of δa * at first order in A to compute the variations of Z. The conditions (6) then yield:

∆r * = δa(θ 2 ) -δa(θ 1 ) 2 sin [(θ 2 -θ 1 )/2] + O(A 2 ), and ∆a * = δa(θ 1 ) + δa(θ 2 ) 2 + O(A 2 ). ( 12
)
The perturbation δa * can finally be expressed as:

δa * (θ) = δa(θ 1 ) + δa(θ 2 ) 2 + δa(θ 2 ) -δa(θ 1 ) 2 • sin [θ -(θ 1 + θ 2 )/2] sin [(θ 2 -θ 1 )/2] + O(A 2 ) = δa(θ 1 ) + δa(θ 2 ) -δa(θ 1 ) 2 cos [(θ 2 -θ 1 )/2] sin [(θ 2 -θ 1 )/2] sin(θ -θ 1 ) + 2 sin 2 [(θ -θ 1 )/2] + O(A 2 ) = δa(θ 2 ) + δa(θ 2 ) -δa(θ 1 ) 2 cos [(θ 2 -θ 1 )/2] sin [(θ 2 -θ 1 )/2] sin(θ -θ 2 ) -2 sin 2 [(θ -θ 2 )/2] + O(A 2 ). ( 13
)
The second and third expressions of δa * (θ) in Eq. ( 13) only involves (θ -θ 1 ) and (θ -θ 2 ) respectively, thanks to the careful choice of θ * . They are analogous to [START_REF] Leblond | Second-order coplanar perturbation of a semi-infinite crack in an infinite body[END_REF]'s Eq. ( 7), and prove key in the evaluation of the principal value integral of our Eq. ( 8) (see Appendix A for details).

Next, we use another trick, analogous to the partial fraction decomposition of [START_REF] Leblond | Second-order coplanar perturbation of a semi-infinite crack in an infinite body[END_REF]'s Eq. ( 9), and propose the following decomposition:

1 sin 2 [(θ -θ 1 )/2] sin 2 [(θ -θ 2 )/2] = 1 sin 2 [(θ 2 -θ 1 )/2] cos 2 [(θ -θ 1 )/2] sin 2 [(θ -θ 1 )/2] + cos 2 [(θ -θ 2 )/2] sin 2 [(θ -θ 2 )/2] + 2 cos[(θ 2 -θ 1 )/2] sin 3 [(θ 2 -θ 1 )/2] cos 3 [(θ -θ 1 )/2] sin[(θ -θ 1 )/2] - cos 3 [(θ -θ 2 )/2] sin[(θ -θ 2 )/2] + 2 sin 2 [(θ 2 -θ 1 )/2] cos 2 [(θ -θ 1 )/2] + cos 2 [(θ -θ 2 )/2] (14)
Combining Eqs. ( 8), (13-14), and after some simplifications, one gets:

δZ (θ 1 , θ 2 ) = - 5 32πa 3 0 δa(θ 1 ) + δa(θ 2 ) sin 2 [(θ 2 -θ 1 )/2] + 1 (8π) 2 a 3 0 sin 2 [(θ 2 -θ 1 )/2] 2π 0 2 cos 2 θ -θ 1 2 + cos 2 θ -θ 2 2 δa(θ)dθ + 1 (8π) 2 a 3 0 sin 2 [(θ 2 -θ 1 )/2] × PV 2π 0 cos 2 [(θ -θ 1 )/2] sin 2 [(θ -θ 1 )/2] [δa(θ) -δa(θ 1 )] dθ + PV 2π 0 cos 2 [(θ -θ 2 )/2] sin 2 [(θ -θ 2 )/2] [δa(θ) -δa(θ 2 )] dθ + 2 cos[(θ 2 -θ 1 )/2] sin 3 [(θ 2 -θ 1 )/2] • PV 2π 0 cos 3 [(θ -θ 1 )/2] sin[(θ -θ 1 )/2] - cos 3 [(θ -θ 2 )/2] sin[(θ -θ 2 )/2] δa(θ)dθ + O(A 2 ) (15)
Details on the derivation of Eq. ( 15) are given in Appendix A. Note that the two first terms of Eq. ( 15) typically emerge from the finiteness of the crack geometry, and go to zero in the limit a 0 → +∞. One can show that [START_REF] Leblond | Second-order coplanar perturbation of a semi-infinite crack in an infinite body[END_REF]'s Eq. ( 10) for the semi-infinite crack is retrieved in the same limit.

We derived here the first-order variation of the kernel Z for a penny-shaped crack loaded in mode I. However, the determination of δa * provided here enables extending our results to mixed mode I+II+III (Favier et al., 2006a), to other geometries like e.g. the circular connection (Gao and Rice, 1987a), or to circular cracks propagating in cohesive materials [START_REF] Lebihain | Quasi-static crack front deformations in cohesive materials[END_REF].

Second-order expansion of the SIF

We consider again a crack front F that has been perturbed from a circular configuration of radius a 0 by an amount δa(θ) = Aϕ(θ), where A is now a small albeit not infinitesimal parameter. The stress intensity factor K I and the kernel Z along the perturbed crack front F can be expanded at several orders in A. Their expansion writes:

           K I (a 0 , Aϕ ; θ 1 ) = K 0 I (a 0 , θ 1 ) + A K 1 I (a 0 , [K 0 I ], [ϕ]; θ 1 ) + A 2 K 2 I (a 0 , [K 0 I ], [ϕ]; θ 1 ) + O(A 3 ) Z(a 0 , Aϕ ; θ 1 , θ 2 ) = Z 0 (a 0 ; θ 1 , θ 2 ) + AZ 1 (a 0 , ϕ ; θ 1 , θ 2 ) + O(A 2 ), (16)
where K 1 I (a 0 , [K 0 I ], [ϕ]; θ 1 ) is the first-order variation of the stress intensity factor induced at θ 1 by the perturbation ϕ applied to a reference circular configuration of radius a 0 and SIF K 0 I , and K2 I (a 0 , [K 0 I ], [ϕ]; θ 1 ) its second-order variation. Similarly, Z 1 (a 0 , ϕ ; θ 1 , θ 2 ) is the first-order variation of the kernel Z 0 induced at (θ 1 , θ 2 ) by the perturbation ϕ from the reference configuration of radius a 0 . K 1 I and Z 0 have been derived by Gao and Rice (1987b), and are recalled in Eqs. (2-3). Z 1 has been derived here in Eq. (15). Our goal here is to provide an analytical expression for K 2 I . To this end, we follow the steps of [START_REF] Leblond | Second-order coplanar perturbation of a semi-infinite crack in an infinite body[END_REF], and superimpose to the primary perturbation δa a proportional yet infinitesimal secondary perturbation δAϕ, where δA is an infinitesimally small parameter. By definition, one has:

K I (a 0 , (A + δA)ϕ ; θ 1 ) -K I (a 0 , Aϕ ; θ 1 ) = ∂K I ∂A (a 0 , Aϕ ; θ 1 )δA + O(δA 2 ) = K 1 I (a 0 , [K 0 I ], [ϕ]; θ 1 ) + 2AK 2 I (a 0 , [K 0 I ], [ϕ]; θ 1 ) + O(A 2 ) δA + O(δA 2 ) (17)
Alternatively, one may consider the front F ′ of reference radius a ′ 0 = a 0 + δAϕ(θ 1 ) perturbed by the quantity δa(θ) = Aϕ(θ), and perturbed it further by a infinitesimal quantity δAϕ(θ) -δAϕ(θ 1 ). Yet, this secondary perturbation is not performed in the direction e r , but rather in the direction normal to the perturbed front F ′ . Given that the error in the resulting front position is of second order in both in δA (see Leblond et al. ( 2012)) 2 , one has [START_REF] Rice | Weight function theory for three-dimensional elastic crack analysis[END_REF]:

K I (a 0 , (A + δA)ϕ ; θ 1 ) -K I (a 0 , Aϕ ; θ 1 ) = K I (a 0 + δAϕ(θ 1 ), Aϕ ; a 0 ; θ 1 ) + PV 2π 0 Z(a 0 + δAϕ(θ 1 ), Aϕ ; θ 1 , θ)K I (a 0 + δAϕ(θ 1 ), Aϕ ; θ) δAϕ(θ) -δAϕ(θ 1 ) ds(θ) -K I (a 0 , Aϕ ; θ 1 ) + O(δA 2 ), ( 18
)
where s is the curvilinear abscissa along the perturbed front F ′ . Equating Eqs. ( 17) and ( 18), one may identify the terms of zero order in A in the expression of ∂K I /∂A:

K 1 I (a 0 , [K 0 I ], [ϕ]; θ 1 ) = ∂K 0 I ∂a (a 0 , θ 1 )ϕ(θ 1 ) + PV 2π 0 Z 0 (a 0 ; θ 1 , θ)K 0 I (a 0 , θ) ϕ(θ) -ϕ(θ 1 ) a 0 dθ, (19)
which is Eq. (2). Identifying the terms of first order in A in the expression of ∂K I /∂A, one has after some simplifications:

K 2 I (a 0 , [K 0 I ], [ϕ]; θ 1 ) = 1 2 ∂ 2 K 0 I ∂a (a 0 , θ 1 )ϕ(θ 1 ) 2 + 1 2 PV 2π 0 Z 0 (a 0 ; θ 1 , θ) ∂K 0 I ∂a (a 0 , θ) ϕ(θ) 2 -ϕ(θ 1 ) 2 a 0 dθ + 1 2 PV 2π 0 Z 0 (a 0 ; θ 1 , θ)K 0 I (a 0 , θ) ϕ(θ) -ϕ(θ 1 ) 2 a 0 dθ + 1 2 PV 2π 0 Z 1 (a 0 , ϕ; θ 1 , θ)K 0 I (a 0 , θ) ϕ(θ) -ϕ(θ 1 ) a 0 dθ + 1 2 PV 2π 0 2π 0 Z 0 (a 0 ; θ 1 , θ)Z 0 (a 0 ; θ, θ ′ )K 0 I (a 0 , θ ′ ) ϕ(θ) -ϕ(θ 1 ) ϕ(θ ′ ) -ϕ(θ) a 2 0 dθdθ ′ . ( 20
)
Details on the derivation of Eq. ( 20) can be found in Appendix B. Equations ( 19) and ( 20) are general equations that prove unfit for analytical calculations, and one has to resort to Fourier series to evaluate them. We shall do it here in the simplified case where the crack is loaded through an axisymmetric system of forces, i.e. when K 0 I is independent of θ. Note however that this is enough to derive second-order analytical solutions for crack propagation in heterogeneous media, even when K 0 I actually depends on θ, as these fluctuations are already of second order in the first-order theory of Gao and Rice (1987b). However, for the sake of completeness, the derivation for a non-axisymmetric

K 0 I is given in Appendix B.
Next, we assume that K 0 I only depends on the radius a 0 of the reference penny-shaped crack. One can decompose ϕ in Fourier series. It reads:

ϕ(θ) = +∞ k=-∞ φk e ikθ , where φk = 1 2π 2π 0 ϕ(θ)e -ikθ dθ. (21)
Then, the first-order contribution of Eq. ( 19) writes as:

K 1 I (a 0 , [K 0 I ], [ϕ]; θ 1 ) = +∞ k=-∞       dK 0 I da (a 0 ) - |k| 2a 0 K 0 I (a 0 )       φk e ikθ 1 , (22)
and the second-order contribution of Eq. ( 20) reads: [START_REF] Leblond | Second-order coplanar perturbation of a semi-infinite crack in an infinite body[END_REF], written in the form simplified by Vasoya et al. (2016b). Note that Eq. ( 20) differs from Eq. ( 8) of [START_REF] Vasoya | A geometrically nonlinear analysis of coplanar crack propagation in some heterogeneous medium[END_REF] only by the term (|k| + |k ′ |)K 0 I (a 0 )/4a 2 0 that typically emerges from the finiteness of the penny-shaped crack geometry. As expected, this term goes to 0 when a 0 → +∞, and one retrieves, in this limit, the results of [START_REF] Leblond | Second-order coplanar perturbation of a semi-infinite crack in an infinite body[END_REF] and [START_REF] Vasoya | A geometrically nonlinear analysis of coplanar crack propagation in some heterogeneous medium[END_REF] for the semi-infinite crack.

K 2 I (a 0 , [K 0 I ], [ϕ]; θ 1 ) = +∞ k=-∞ +∞ k ′ =-∞       1 2 d 2 K 0 I da 2 (a 0 ) - 1 4a 0 dK 0 I da (a 0 )|k + k ′ | + R(k, k ′ ) 16a 2 0 K 0 I (a 0 ) + |k| + |k ′ | 4a 2 0 K 0 I (a 0 )       φk φk ′ e i(k+k ′ )θ1 , ( 23
) where R(k, k ′ ) = 2|k + k ′ |(|k| + |k ′ |) -(k + k ′ ) 2 -k 2 -k ′2 is twice the function F of
One can alternatively express Eq. (16a) in terms of the Fourier coefficient δa k of the perturbation δa:

K I (a 0 , [δa] ; θ 1 ) = K 0 I (a 0 )         1 + k L(k) δa k a 0 e ikθ 1 + k,k ′ H(k, k ′ ) δa k a 0 δa k ′ a 0 e i(k+k ′ )θ 1         + O(A 3 ), ( 24
)
where the first-and second-order convolution kernels L and H read:

                 L(k) = a 0 K 0 I (a 0 ) dK 0 I da (a 0 ) - |k| 2 H(k, k ′ ) = a 2 0 2K 0 I (a 0 ) d 2 K 0 I da 2 (a 0 ) - a 0 K 0 I (a 0 ) dK 0 I da (a 0 ) |k + k ′ | 4 + R(k, k ′ ) 16 + |k| + |k ′ | 4 .
(25)

Validation of the proposed theory

We aim now at validating the second-order expansion of the stress intensity factor given in Eq. ( 24). One may consider that retrieving [START_REF] Leblond | Second-order coplanar perturbation of a semi-infinite crack in an infinite body[END_REF][START_REF] Vasoya | A geometrically nonlinear analysis of coplanar crack propagation in some heterogeneous medium[END_REF]'s results in the limit of very large crack radius a 0 proves enough to validate our calculations. Yet, [START_REF] Leblond | Second-order coplanar perturbation of a semi-infinite crack in an infinite body[END_REF]'s second-order formula has been found conflicting with those previously derived by Adda-Bedia et al. (2006) and [START_REF] Katzav | Roughness of moving elastic lines: Crack and wetting fronts[END_REF], and its validity has been checked through numerical simulations only. As a consequence, one should seek further validation cases to assess the correctness of Eq. ( 24). Here, we follow two different routes: first, (i) we consider elliptical cracks as perturbed circular cracks and test our formula against analytical results obtained for this geometry [START_REF] Irwin | Fracture[END_REF]. It provides a critical test for our second-order theory, and yields a first analytical validation of [START_REF] Leblond | Second-order coplanar perturbation of a semi-infinite crack in an infinite body[END_REF]'s formula for the semi-infinite crack. Next, (ii) we use [START_REF] Lazarus | Brittle fracture and fatigue propagation paths of 3d plane cracks under uniform remote tensile loading[END_REF]'s model that computes the mode I SIF variations at all orders to compare our second-order expansions to numerical results obtained for circular cracks deformed by cosinusoidal perturbations of varying modes and amplitudes. This ultimately allows us to estimate the validity range of our second-order theory.

SIF distribution along the front of an elliptical crack

In contrast with models based on the semi-infinite crack geometry [START_REF] Leblond | Second-order coplanar perturbation of a semi-infinite crack in an infinite body[END_REF][START_REF] Vasoya | A geometrically nonlinear analysis of coplanar crack propagation in some heterogeneous medium[END_REF], our asymptotic theory derived for penny-shaped cracks can be tested against analytical solutions obtained for planar cracks with elliptical front shapes. Here, we consider a crack whose front describes an ellipsis centered in O, of short semi-axis c aligned with e z , and long semi-axis b aligned with e x . In the polar coordinate system of Fig. 1, the front position is parametrized by the function:

a(θ) = bc b 2 sin 2 (θ) + c 2 cos 2 (θ) . ( 26
)
The crack is loaded by a uniform remote tensile stress σ. The distribution of mode I stress intensity factor along the elliptical crack front is known exactly from the work of [START_REF] Irwin | Fracture[END_REF]. It reads:

K e I (θ) = σ √ πc E(κ) c 4 /b 4 cos 2 (θ) + sin 2 (θ) c 2 /b 2 cos 2 (θ) + sin 2 (θ) 1/4 , (27)
where E(κ) is the elliptical integral of the second kind evaluated at κ = 1c 2 /b 2 1/2 . When b = c = a 0 , one retrieves the case of a circular crack front of radius a 0 :

K 0 I (a 0 ) = 2 π σ √ πa 0 , (28)
so that the convolution kernels of Eq. ( 25) read:

             L(k) = - |k| 2 + 1 2 H(k, k ′ ) = R(k, k ′ ) 16 + |k| + |k ′ | 4 - |k + k ′ | 8 - 1 8 . (29)
To derive asymptotic values of K e I from our second-order theory, one needs to choose (i) a small parameter ϵ that quantifies the deviation from circularity of the elliptical crack front, and (ii) the radius a 0 of the reference circular configuration. We follow here the reasoning of Gao and Rice (1987b), and take ϵ = 1c/b. However, they considered a reference radius a 0 equal to the crack front position at the point of SIF evaluation, i.e. a 0 = a(θ 1 ). Consequently, the perturbation δa(θ) = a(θ)a 0 is zero at θ = θ 1 , and the double sum of Eq. ( 23) containing terms in |k| + |k ′ |:

k,k ′ |k| + |k ′ | 4 δa k δa k ′ e i(k+k ′ )θ 1 =        k |k| 2 δa k e i(k+k ′ )θ 1        × δa(θ 1 ) (30)
is also null. Yet, these terms are precisely those we need to check first, as they vanish in the limit of very large crack radii, and cannot be tested against [START_REF] Vasoya | A geometrically nonlinear analysis of coplanar crack propagation in some heterogeneous medium[END_REF]'s formula for the semi-infinite crack. Here, we rather take the length c of the short semi-axis as reference radius a 0 . Following Eq. ( 26), the perturbation δa = aa 0 can be expanded as:

δa(θ)/a 0 = ϵA 1 (θ) + ϵ 2 A 2 (θ) + O(ϵ 3 ) = ϵ 1 2 + 1 2 cos(2θ) + ϵ 2 5 16 + 1 2 cos(2θ) + 3 16 cos(4θ) + O(ϵ 3 ), ( 31
)
where A 1 and A 2 are the first-and second-order expansions of δa/a 0 . Similarly, the exact SIF distribution K e I given in Eq. ( 27) can be expanded to second-order in ϵ:

K e I (θ)/K 0 I (a 0 ) = 1 + ϵ 1 4 - 1 4 cos(2θ) + ϵ 2 11 64 - 3 16 cos(2θ) - 11 64 cos(4θ) + O(ϵ 3 ). ( 32
)
One seeks now to find back Eq. ( 32) from our non-linear theory. Using Eq. ( 32) in Eq. ( 24), the second-order estimate K p I of the mode I SIF reads:

K p I (θ)/K 0 I (a 0 ) = 1 + ϵ        k L(k) Â1 k e ikθ        + ϵ 2         k L(k) Â2 k e ikθ + k,k ′ H(k, k ′ ) Â1 k Â1 k ′ e i(k+k ′ )θ         + O(ϵ 3 ), (33)
where Â1 k and Â2 k are the Fourier coefficients of A 1 and A 2 . Using Eqs. ( 29) and ( 31), one finds:

k L(k) Â1 k e ikθ = 1 4 - 1 4 cos(2θ), k L(k) Â2 k e ikθ = 5 32 - 1 4 cos(2θ) - 9 32 cos(4θ),
and

k,k ′ H(k, k ′ ) Â1 k Â1 k ′ e i(k+k ′ )θ = 1 64 + 1 16 cos(2θ) + 7 64 cos(4θ), (34)
which yields the result expected from Eq. ( 32). Note that we also checked that K e I and our expansion K p I are equal up to second-order when one takes a 0 = a(θ 1 ) as Gao and Rice (1987b) did. Further validations of our second-order theory are provided in Appendix C, even for the case of a non-axisymmetric K 0 I .

The correctness of our Eq. ( 24) in turn provides a first analytical validation of [START_REF] Leblond | Second-order coplanar perturbation of a semi-infinite crack in an infinite body[END_REF]'s formula, which is obtained in the limit a 0 → +∞. Our non-linear perturbative approach yields estimates of the stress intensity factor along elliptical crack fronts for b/c ≃ 1, even when it arises from discontinuous stress distributions (e.g. a point-force loading). It is complementary to other methods available in the literature that provide K I estimates along elliptical crack fronts for arbitrary b/c ratio, but considering polynomial stress distributions [START_REF] Shah | Stress intensity factor for an elliptical crack under arbitrary normal loading[END_REF][START_REF] Kassir | Three-dimensional crack problems : a new selection of crack solutions in three-dimensional elasticity[END_REF][START_REF] Atroshchenko | Stress intensity factor for an embedded elliptical crack under arbitrary normal loading[END_REF]. However, as we will see next, its scope is much wider as it can be used on perturbed fronts of arbitrary shapes.

SIF distribution along a sinusoidal crack front

Now that our second-order theory has been validated considering elliptical crack fronts, one needs to assess its performances in estimating the stress intensity factor distribution along randomly perturbed crack fronts. As elliptical crack fronts are associated with second-order perturbations consisting of the superposition of several deformation modes (see Eq. ( 31)), they do not allow us to quantify the errors made by our second-order model for a given perturbation mode of fixed amplitude. We rather investigate a wavy crack front described by:

a(θ) = a 0 [1 + ϵ cos(kθ)] , (35)
where the integer k ≥ 2 is the mode of the perturbation, and ϵ is a small parameter that quantifies its amplitude. We assume next that the crack parametrized by Eq. ( 35) is loaded by a uniform remote tensile stress σ. As there is no analytical formula available for the mode I stress intensity factor along the front of such wavy cracks, one needs to compute K I numerically. Here, we use the method of [START_REF] Lazarus | Brittle fracture and fatigue propagation paths of 3d plane cracks under uniform remote tensile loading[END_REF] and [START_REF] David | On the key role of crack surface area on the lifetime of arbitrarily shaped flat cracks[END_REF] that computes K I along an arbitrary crack front using the perturbative approach of Rice (1989) iteratively. Namely, (i) we start from a penny-shaped configuration for which the mode I SIF K I and the fundamental kernel W are known. Then, (ii) we discretize the crack front in N equidistant points, and sequentially deform it by a small amount δa(θ) = ∆ϵ cos(kθ) with ∆ϵ = 0.0025, and (iii) update both K I and W using formulae similar to Eqs. ( 2) and ( 5). As ∆ϵ ≪ 1, the first-order assumption of [START_REF] Rice | Weight function theory for three-dimensional elastic crack analysis[END_REF] applies, and K I is computed at all orders in ϵ. The update of the stress intensity factor and the fundamental kernel are costly procedures, but this method still demands less computational resources than e.g. simulations based on the finite element method (FEM) or the boundary element method (BEM), as only the crack front needs to be meshed. This method is also powerful in modeling far-field loading applied at infinity or hypotheses of infinite body. Table 1: Values of perturbation amplitude ϵ (±0.0025) above which the errors in the SIF variation δK I made by the second-order theory exceeds a certain threshold value ∈ {1%, 5%, 25%}. Values of ϵ above 0.5 have not been investigated in this study.

|| • || 0 || • || 2 1% 5% 20% 1% 5% 20% k = 2
Using our second-order model of Eq. ( 24), together with the convolution kernels of Eq. ( 29), the stress intensity factor variations δK p I from the uniform K 0 I reference reads:

δK p I (θ)/K 0 I (a 0 ) = - k -1 2 ϵ cos(kθ) - k 2 -4k + 1 16 ϵ 2 + k 2 + 2k -1 16 ϵ 2 cos(2kθ) + O(ϵ 3 ), ( 36
)
where K 0 I (a 0 ) is defined in Eq. ( 28).

The numerical results using the method of [START_REF] Lazarus | Brittle fracture and fatigue propagation paths of 3d plane cracks under uniform remote tensile loading[END_REF] (in black solid line) are compared to the second-order prediction of Eq. ( 36) (in dashed red line) and to its first-order truncation (in dotted orange line) in Fig. 3. For k = 2 and ϵ = 0.01 (see Fig. 3a), the three quantities are close. When ϵ = 0.1, the first-order assumption is no longer valid, and the linear estimate departs from the numerical results, while the second-order prediction remains accurate. As ϵ grows even larger (ϵ = 0.25), the second-order prediction starts overestimating the stress intensity factor on the regions where the crack is less elongated. Note however that the second-order model do predict an upward shift in the average of the SIF distribution, while the first-order prediction remains symmetric with respect to δK I = 0. For k = 5 (see Fig. 3b), the validity range of our theory shrinks, and departure of the second-order prediction from the numerical results is already noticeable when ϵ = 0.1. When ϵ = 0.25, the discrepancies are fairly large, and our model fails in predicting the strong decrease of SIF amplitude at θ = π/5. We observe that the errors made by our second-order model around this point are even larger than those made by the first-order theory. These preliminary observations are further investigated in Fig. 3c-d where we show the evolution with ϵ of the (maximal) error using the L 0 norm:

|| error || 0 = max θ∈[0,2π] |δK p I (θ) -δK num I (θ)|/K 0 I (a 0 ),
and the (average) error using the L 2 norm:

|| error || 2 = 1 2π/K 0 I (a 0 ) 2π 0 δK p I (θ) -δK num I (θ) 2 dθ.
Here, δK num I is the mode I SIF variation predicted using Lazarus (2003)'s method. Values of ϵ for which the errors made by the first-and second-order models exceeds a threshold value are listed in Table 3.2. Note that we report errors made by our second-order prediction against a numerical estimate of K I , since its exact value is not known. Nonetheless, for the discretization adopted here (N = {252, 378, 680, 850} for k = {2, 3, 4, 5}, and δa = 0.0025), the numerical error made in evaluating δK num I can be estimated to be less than 0.2%. This estimation builds on the error analysis made by [START_REF] Lazarus | Brittle fracture and fatigue propagation paths of 3d plane cracks under uniform remote tensile loading[END_REF] on K I evaluation along elliptical crack fronts with b/c = 3 (i.e. ϵ = 0.66), that necessitates more small perturbations iterations than here where ϵ ≤ 0.25.

We observe nonetheless that our second-order formula generally yields more accurate predictions than the linear theory of Gao and Rice (1987b), and that quantitative agreement between the numerical results and the second-order expansion of Eq. ( 36) is found for ϵ ≲ 0.1 and k ∈ [0, 5]. This is further confirmed by the evolution of the average SIF ⟨K I ⟩ with ϵ (see Fig. 3e): Gao and Rice (1987b) (in dotted orange line) and our second-order prediction (in dashed red line). The errors on the SIF variations between the numerical simulations (in solid line) and the first-order (in dotted line) and second-order (in dashed line) models are compared using (c) L 0 (maximum) or (d) L 2 (Euclidean) norms. They increase with both the perturbation mode k and amplitude ϵ. Deviations of the numerical results from the first-and second-order predictions can be directly observed on (e) the average ⟨K I ⟩ of the stress-intensity factor along the crack front, or (f) its amplitude K max I -K min I , for which first and second-order estimates coincide.

⟨K I ⟩ = 1 2π 2π 0 K I (θ)dθ = K 0 I (a 0 ) × 1 - k 2 -4k + 1 16 ϵ 2 + O(ϵ 3 ), (37)
1 st order 2 nd order Numerics k = 2 k = 3 k = 4 k = 5 1 st order 2 nd order Numerics k = 2 k = 2 k = 2 k = 5 k = 5 k = 5
and the SIF amplitude K max I -K min I (see Fig. 3f):

K max I -K min I = max θ∈[0,2π] K I (θ) -min θ∈[0,2π] K I (θ) = K 0 I (a 0 ) × (k -1)ϵ + O(ϵ 3 ). ( 38
)
The change in average SIF ⟨K I ⟩ appears indeed of second order in ϵ, and is thus overlooked by the first-order theory. Moreover, we observe a change in the sign of the second-order terms in ⟨K I ⟩ between k = 3 and 4. On the contrary, there is no second-order contribution to the SIF amplitude K max I -K min I , and higher-order theories would be required to predict the strong decrease of K I near θ = π/k observed in Fig. 3b. We observed that the abrupt change of slope in the SIF amplitude computed using Lazarus (2003)'s method is associated with a change in the point where K I is minimal, from θ = π/4k (see the third panel of Fig. 3b) to θ = π/k.

We conclude that (i) our second-order theory generally yields better estimates of the SIF variations along perturbed crack fronts than the first-order model of Gao and Rice (1987b), but (ii) its accuracy decreases with increasing perturbation amplitude ϵ and mode k. As we will see next, the latter may be rather inconsequential, as the amplitude of the front perturbations is usually a fast decreasing function (∝ 1/k) of their mode k.

Coplanar crack propagation in weakly disordered materials

The influence of microscopic fluctuations of toughness on the critical load at which a material may fail has been studied quite extensively for the semi-infinite crack geometry since the seminal study of [START_REF] Roux | Effective toughness of heterogeneous brittle materials[END_REF]. Various methods have been used in the literature to address this problem, like e.g. phase-field models [START_REF] Hossain | Effective toughness of heterogeneous media[END_REF][START_REF] Brach | Anisotropy of the effective toughness of layered media[END_REF] or minimal surfaces [START_REF] Ernesti | Computing the effective crack energy of heterogeneous and anisotropic microstructures via anisotropic minimal surfaces[END_REF][START_REF] Michel | Merits and limits of a variational definition of the effective toughness of heterogeneous materials[END_REF]. But perturbative approaches have played an important role in homogenizing fracture properties. Efforts have been made in both proving the existence of a finite effective toughness for random media [START_REF] Dirr | Pinning of interfaces in random media[END_REF][START_REF] Dondl | Threshold phenomenon for homogenized fronts in random elastic media[END_REF], and predicting its value from the statistical [START_REF] Roux | Self-consistent scheme for toughness homogenization[END_REF]Patinet et al., 2013b;Démery et al., 2014b) and spatial (Démery et al., 2014a;[START_REF] Lebihain | Towards brittle materials with tailored fracture properties: the decisive influence of the material disorder and its microstructure[END_REF] distribution of toughness K Ic (z, x). In particular, it has been shown that the material is toughened by the presence of disorder, i.e. the effective toughness K eff Ic is always larger than the average value K 0 Ic of the local toughness field. A major pitfall of these studies is that they predict a toughening that scales as the second-order of the disorder intensity (measured by e.g. the standard deviation σ of the toughness field) from a first-order expansion of the stress-intensity factor. Hence, they neglect potential contributions arising from the non-linearities in K I and may be inaccurate in predicting K eff Ic . Moreover, the semi-infinite crack geometry considered in these works also prevents investigating size-effects emerging from the finite length of the crack.

Here, we investigate crack propagation when K Ic is fluctuating randomly within the failure plane y = 0, and we focus on the effective toughness K eff Ic of such composite materials. We propose a rigorous homogenization framework for weakly disordered materials that correctly predicts the non-linear toughening arising from local fluctuations of toughness from a second-order expansion of the stress intensity factor. We also examine the influence of the finiteness of the crack geometry on material toughening, as previous studies mainly focused on the semi-infinite crack geometry. Our analytical predictions are finally compared to the results of ∼ 1 millions of efficient numerical simulations of crack propagation in heterogeneous brittle media.

Problem statement and field generation

We consider a penny-shape crack of radius a ini centered in O. The crack propagates along the failure plane located at y = 0, and interacts with the non-uniform toughness field:

K Ic (r, θ) = K 0 Ic + σ f (r, θ), (39)
where K 0 Ic corresponds to the spatial average of the toughness field, σ is its standard deviation, and f is a fluctuation term of zero average and unit variance. Fluctuations are either described by the function f in polar coordinates or by f * in Cartesian coordinates. The two functions are linked through the equation:

f (r, θ) = f * (z, x), where z = r cos θ and x = r sin θ. ( 40
)
An example of toughness field is given in Fig. 4a. It is characterized its probability density function P that measures the statistical distribution of the local toughness. Here, we consider that the values of K Ic are uniformly distributed between two extremal values K min Ic and K max Ic (see Fig. 4b), so that:

P(K Ic ) = 1 K max Ic -K min Ic . (41)
Consequently, its mean value K 0 Ic reads:

K 0 Ic = K max Ic + K min Ic 2 , ( 42
)
and the intensity σ = σ/K 0 Ic of the material disorder is equal to:

σ = 1 √ 3 K max Ic -K min Ic K max Ic + K min Ic . ( 43
)
The toughness field is further characterized by its spatial structure, described e.g. by the two-points correlation function of f . In the following, E [•] denotes the expectation of a random variable over multiple disorder realizations.

We consider (i) statistically homogeneous and (ii) ergodic materials, so that ensemble averages (i) do not depend on the observation point M, and (ii) are equivalent to spatial average over large enough surfaces. Using these notations, the two-points correlation function of the fluctuations reads:

E f * (z, x) f * (z ′ , x ′ ) = F z (|z -z ′ |)F x (|x -x ′ |), (44)
where F z and F x are the correlation functions in the z and x direction, respectively. Here we assume that the material is isotropic and that the correlations follows a bell-curve, so that:

F z (u) = F x (u) = F (u) = e -u 2 /d 2 , ( 45
)
where d is the characteristic size of the disorder (see Fig. 4c).

Realizations of K Ic fields are produced following the procedure of [START_REF] Albertini | Stochastic properties of static friction[END_REF], suitably extended to twodimensional surfaces. First, a Gaussian random field of zero average, unit variance, and controlled power-spectrum is constructed using the Python package FyeldGenerator [START_REF] Cadiou | FyeldGenerator[END_REF]. We then apply a non-linear mapping from the Gaussian cumulative distribution function to that of the uniform distribution. This transformation is expected not to affect much the correlation shape F and the correlation length d, as F is positive [START_REF] Grigoriu | Stochastic Calculus: Applications in Science and Engineering[END_REF] (see Fig. 4c).

Computation of the crack front position

The crack is loaded by a pair of symmetric tensile forces of magnitude P applied at the crack center O. It gives rise to a stress intensity K 0 I along the front of a penny-shaped crack of radius a 0 , which reads [START_REF] Tada | The Stress Analysis of Cracks Handbook[END_REF]:

K 0 I (a 0 ) = P (πa 0 ) 3/2 . ( 46
)
The magnitude P of the tensile forces slowly increases in time, so that the crack extends quasi-statically due to the stabilizing nature of the loading. The successive crack front positions are computed from Irwin's criterion:

K I (θ) = K Ic (a(θ), θ). ( 47
)
K I (θ) is usually computed from the front deformation δa from its reference circular configuration of radius a 0 using the first-order expansion of [START_REF] Rice | First-order variation in elastic fields due to variation in location of a planar crack front[END_REF], while K Ic (a(θ), θ) is computed exactly i.e. at all orders in the front deformation.

The non-linear Eq. ( 47) is then solved by building on a viscous regularization of the evolution equation (Patinet et al., 2013b;[START_REF] Lebihain | Towards brittle materials with tailored fracture properties: the decisive influence of the material disorder and its microstructure[END_REF], or by resorting to a comparison theorem to compute efficiently the successive state equilibria [START_REF] Rosso | Roughness at the depinning threshold for a long-range elastic string[END_REF]Démery et al., 2014a). Here, we follow a different approach, and expand the front position to second-order in the disorder intensity σ = σ/K 0 Ic : where a 0 is the radius of a crack propagating on a homogeneous medium of toughness K 0 Ic under the action of the loading P. A 1 and A 2 are fluctuation terms that emerge from the local variations f of the toughness field. One can then expand the local SIF K I (θ) at second-order using our Eq. ( 16):

a(θ) = a 0 1 + σA 1 (a 0 , θ) + σ2 A 2 (a 0 , θ) + O( σ3 ), (48)
K I (θ) = K 0 I (a 0 ) + σK 1 I (a 0 , [K 0 I ], [a 0 A 1 ]; θ) + σ2 K 1 I (a 0 , [K 0 I ], [a 0 A 2 ]; θ) + K 2 I (a 0 , [K 0 I ], [a 0 A 1 ]; θ) + O( σ3 ), ( 49
)
where K 1 I and K 2 I are given by Eq. ( 24), with the following convolution kernels:

             L(k) = - |k| 2 - 3 2
H(k, k ′ ) = R(k, k ′ ) 16 + |k| + |k ′ | 4 + 3|k + k ′ | 8 + 15 8 . (50)
Similarly, the local toughness K Ic (a(θ), θ) can be expanded at second-order following:

K Ic (a(θ), θ) = K 0 Ic 1 + σ f (a 0 , θ) + σ2 ∂ f ∂r (a 0 , θ)a 0 A 1 (a 0 , θ) . ( 51
)
The crack front shape is then computed by solving Irwin's criterion of Eq. ( 47) in cascade, i.e. at successive but increasing orders in σ. Starting with zero order, one gets from Eq. ( 46):

K 0 I (a 0 ) = K 0 Ic ⇒ a 0 = 1 π       P K 0 Ic       2/3 , (52)
which sets the instantaneous value of the loading P. One could have alternatively set the magnitude P of the opening forces, and solve Eq. ( 52) to find the reference radius a 0 . At first order, one has:

K 1 I (a 0 , [K 0 I ], [a 0 A 1 ]; θ) = K 0 Ic f (a 0 , θ) ⇒ k L(k) Â1 k (a 0 )e ikθ = k fk (a 0 )e ikθ , (53)
so that the k-th Fourier coefficient of A 1 reads:

Â1 k (a 0 ) = fk (a 0 ) L(k) , (54)
where L(k) is given in Eq. ( 50). We observe here that the perturbation amplitude is a fast decreasing function of its mode k (∝ 1/k), and that crack fronts are generally stiffer to small wavelength perturbations. Collecting the secondorder terms in Eq. ( 47), one finds:

K 1 I (a 0 , [K 0 I ], [a 0 A 2 ]; θ) + K 2 I (a 0 , [K 0 I ], [a 0 A 1 ]; θ) = K 0 Ic a 0 ∂ f ∂r (a 0 , θ)A 1 (θ) ⇒ k L(k) Â2 k (a 0 )e ikθ = k,k ′ a 0 ∂ f ∂r k ′ (a 0 ) Â1 k-k ′ (a 0 )e ikθ - k,k ′ H(k ′ , k -k ′ ) Â1 k ′ (a 0 ) Â1 k-k ′ (a 0 )e ikθ , (55)
where H(k, k ′ ) is given in Eq. ( 50). We observe that each mode Â2 k of the second-order fluctuations A 2 are linked to all perturbation modes of A 1 . Injecting Eq. ( 54), the k-th Fourier coefficient of A 2 reads:

Â2 k (a 0 ) = 1 L(k) k ′         a 0 L(k -k ′ ) ∂ f ∂r k ′ (a 0 ) fk-k ′ (a 0 ) - H(k ′ , k -k ′ ) L(k ′ )L(k -k ′ ) fk ′ (a 0 ) fk-k ′ (a 0 )         . (56)
In practice, the crack front shape a(θ) and the critical loading P are computed using Eqs. ( 48), ( 52), (54), and (56).

For a given reference crack radius a 0 , the magnitude of the applied force P is computed from Eq. ( 52). The crack front is then discretized in n = 1024 intervals of equal length, and Eqs. (54), and (56) are solved using fast Fourier Transforms. The computation time for one reference front position a 0 scales as O(n 2 ln n). Examples of crack front shapes and loading evolution are given in Fig. 4a and4d.

Our numerical approach differs from those available in the literature, as only second-order contributions of the toughness fluctuations are taken into account in coherence with the second-order expansion of K I . We expect it to break as soon as the front deformation δa(θ) = a(θ)a 0 gets larger than the characteristic length scale d of the toughness variations. Indeed, for larger front distortions, the second-order expansion of the local toughness in Eq. ( 51) may not be valid anymore. This strongly limits (i) the intensity σ of the material disorder, and (ii) the ratio a 0 /d between the crack radius to heterogeneity size considered in the simulations. Namely, we measured from our numerical computations that we could not model crack radii larger than a 0 = 20d for σ = 0.1.

Effective toughness of weakly disordered materials

Due to the strong decrease of K 0 I with crack advance (see Eq. ( 46)), we do not observe any snapback instability in the evolution of P with the average crack radius ⟨a⟩ (see Fig. 4d), which reads:

⟨a⟩ /a 0 = 1 + σ Â1 0 (a 0 ) + σ2 Â2 0 (a 0 ) + O( σ3 ). ( 57
)
As a result, all states computed from Eqs. ( 52), (54), and ( 56) are state equilibria. This suggests that the instability nuclei, often referred to as Larkin domains, are always larger than the crack perimeter P ≃ 2π ⟨a⟩. Accordingly, our simulations should fall down in the weak pinning regime identified by [START_REF] Roux | Effective toughness of heterogeneous brittle materials[END_REF] for small disorder intensity σ. In this regime, the material toughening:

T = K eff Ic /K 0 Ic -1 (58)
is expected to be zero [START_REF] Roux | Effective toughness of heterogeneous brittle materials[END_REF]. Our goals here are (i) to show that this is wrong if one considers a secondorder expansion of both the stress-intensity factor and the toughness field, and (ii) to assess the dependence of the observed toughening on the crack size a 0 .

At the macro-scale, one may overlook the undulations of the crack front and consider the perturbed crack as a circular crack of radius ⟨a⟩ (see Fig. 4a), which departs from its reference value a 0 under the action of the material disorder (see Eq. ( 57)). As the crack propagates under the action of the tensile forces, one may then define an apparent, or effective, toughness of the heterogeneous media from the magnitude P of the tensile forces required to make a circular crack of radius ⟨a⟩ propagate:

K eff Ic = P (π ⟨a⟩) 3/2 , ( 59
)
so that:

K eff Ic = K 0 I (a 0 ) a 0 ⟨a⟩ 3/2 = K 0 Ic a 0 ⟨a⟩ 3/2 . ( 60
)
We observe that the toughening at the macroscale emerges from the difference between the average front position ⟨a⟩ and its reference value a 0 . If the front is pinned by the disorder (⟨a⟩ < a 0 ), the composite appears tougher at the macroscale than the reference homogeneous material of toughness K 0 Ic . On the contrary, if the disorder facilitates crack propagation (⟨a⟩ > a 0 ), the material is weakened by the heterogeneities. Using Eq. ( 57), one finds:

K eff Ic /K 0 Ic = K 0 I (⟨a⟩)/K 0 I (a 0 ) = 1 + σL(0) Â1 0 (a 0 ) + σ2 L(0) Â2 0 (a 0 ) + H(0, 0) Â1 0 (a 0 ) 2 , ( 61
)
where L(0) = a 0 (∂K 0 I /∂a)(a 0 )/K 0 I (a 0 ) and H(0, 0) = a 2 0 (∂ 2 K 0 I /∂a 2 )(a 0 )/2K 0 I (a 0 ). An example of the evolution of K eff Ic with the average crack radius ⟨a⟩ is given in Fig. 4e. We observe in Fig. 4d that the effective toughness is a random variable that seems to fluctuate around K 0 Ic and between the two extremal values K min Ic and K max Ic of the toughness distribution. To have a better perspective on the statistical distribution of the random variable K eff Ic , we perform more than 1 million Monte-Carlo simulations of crack propagation in heterogeneous media, following the procedure described in Section 4.2. Namely, the disorder intensity is varied between 11 different values σ ∈ [0.001, 0.01]. N = 100 ′ 000 independent and identically distributed realizations of the toughness field K Ic (z, x) are generated for each value of σ. For each disorder realization, the crack front position is computed from 16 reference radius a 0 ∈ [0. Eq . ( 62)

Eq . ( 63) from which one can estimate the average front position ⟨a⟩ using Eq. ( 57), and the effective toughness K eff Ic following Eq. ( 61).

The results of our statistical analysis are summarized in Fig. 5. For σ = 0.1, we observe that the effective toughness is indeed bounded by the two extremal values K min Ic and K max Ic of the toughness distribution (see Fig. 5a), as expected from Irwin's criterion. For a small crack radius a 0 = 0.1 d, the effective toughness follows the probability distribution of the local toughness K Ic (see Fig. 4b), as the crack front is often embedded into one single heterogeneity. As the crack grows larger, the distribution of the effective toughness departs from that of the local toughness and strongly resemble a Gaussian distribution. Interestingly, the variance of the distribution looks to decrease with increasing crack to heterogeneity size ratio a 0 /d. One can also notice that the average value of the toughness distribution looks slightly below K 0 Ic for a 0 /d = 1, and closer to it for a 0 /d = 10. Departure of this average value from K 0 Ic seems in general very small with respect to the fluctuations of K eff Ic , so that both quantities may not be of the same order in σ. We also observe in Fig. 5b that the fluctuations of the effective toughness increase with the disorder intensity σ.

In the following, we do not try to predict the exact probability distribution of K eff Ic in Fig. 5a-b, but rather focus on its statistical descriptors, namely its average value and its standard deviation. As the values of K eff Ic seems to gather around K 0 Ic , we concentrate instead on the statistics of the toughening T = K eff Ic /K 0 Ic 1. By performing ensemble averages on Eq. ( 61), one can show that its average reads:

E [T ] =       σ K 0 Ic       2        k 1 L(k) Ĝk (a 0 /d) - H(k, -k) L(k) 2 Fk (a 0 /d) + H(0, 0) L(0) 2 F0 (a 0 /d)        + O              σ K 0 Ic       3        , (62)
while its standard deviation writes as:

std [T ] =       σ K 0 Ic       F0 (a 0 /d) 1/2 + O              σ K 0 Ic       2        , (63)
where F and G are the functions:

F(a 0 /d, θ) = e -4(a 0 /d) 2 sin 2 (θ/2) , and G(a 0 /d, θ) = -4(a 0 /d) 2 sin 2 (θ/2)e -4(a 0 /d) 2 sin 2 (θ/2) , (64) which both relate to the correlation shape F of the fluctuations f . Details on the derivation of Eqs. ( 62) and ( 63) are given in Appendix D. Several comments are in order:

• We observe that the average toughening E [T ] scales as σ2 , while std [T ] increases linearly with σ. As a consequence, is extremely difficult to measure E [T ] from our "naive" Monte-Carlo simulations, as the central limit theorem states that the fluctuations of E [T ] are proportional to σ/ √ N, where N is the number of realizations. To circumvent this issue, we use a variance reduction technique, and resort to antithetic variates. Let us call Ω 1 the ensemble containing all N = 100 ′ 000 realizations generated here for fixed values of σ and a 0 /d. We attribute to each realization an integer i, and note T i 1 the value of T measured at a 0 for a crack propagating in the toughness field K i Ic = K 0 Ic (1 + σ f i ). We now consider the ensemble Ω 2 containing the 100'000 realizations corresponding to the toughness field KIc i = K 0 Ic (1σ f i ) (opposite fluctuations), and note T i 2 the value of T measured at a 0 in this case. T i 2 can be computed at almost no cost, as it corresponds to a change of sign in front of A 1 only. As both f i andf i are identically distributed (zero average, unit variance, uniform probability distribution, same two-points correlation function), one has:

E [T ] = E [T 1 ] + E [T 2 ] 2 ≃ i∈Ω 2 T i 1 + i∈Ω 2 T i 2 2N , (65)
but one expects that the covariance of the random variables T 1 and T 2 is negative, so that the variance of T = (T 1 + T 2 )/2 is greatly reduced. Using antithetic variates, one observe from Fig. 5c that the average toughening E [T ] do scales as σ2 , while its standard deviation std [T ] is only linear in σ (see Fig. 5d). Note that numerical validations of these scalings can only be obtained through the use of variance reduction methods, as no convergence was found even for N = 100 ′ 000. Using our antithetic variates, one observes satisfactory scaling for N as low as N ≃ 250. Our trick may find some use if one wanted to compare the results of perturbative calculations to that of e.g. BEM simulations.

• In contrast with the predictions of [START_REF] Roux | Effective toughness of heterogeneous brittle materials[END_REF], Patinet et al. (2013b), andDémery et al. (2014b), the effective toughness K eff Ic is not equal to the spatial average K 0 Ic of the toughness field while being in the weak pinning regime. In the semi-infinite crack limit (a 0 /d → +∞), the average toughening E [T ] is equal to E [T ] = 0.5(σ/K 0 Ic ) 2 (see Fig. 6a). This toughening contribution is dictated by non-linearities in K I , as only the terms arising from the second-order expansion of K I actually play a role in determining E [T ] (through the terms involving the convolution kernel H in Eq. ( 62)). Moreover, the toughening predicted here in the weak pinning regime is comparable to that measured by Démery et al. (2014b) and Démery et al. (2014a) in the strong pinning regime, where crack propagation is no more smooth and occurs as the succession of depinning instabilities [START_REF] Roux | Effective toughness of heterogeneous brittle materials[END_REF]. Indeed, the predictions E [T ] of Eq. ( 62) for a semi-infinite crack amounts to one fourth of the toughening measured by Démery et al. (2014b) (see their Eq. ( 20)), and half that predicted by Démery et al. (2014a) (see their Eqs. ( 34) and ( 35)). We conclude that one cannot neglect non-linearities in the stress intensity factor variations when estimating the effective toughness of brittle composites.

This calls for a new derivation of the theoretical predictions of Patinet et al. (2013b) and Démery et al. (2014b) for a semi-infinite crack in the strong pinning regime, involving now the second-order terms in the expansion of K I or G derived by [START_REF] Vasoya | A geometrically nonlinear analysis of coplanar crack propagation in some heterogeneous medium[END_REF]. To do so, one can follow the route proposed by Démery et al. (2014a), which is based on a second-order expansion of the front position similar to ours, but taking into account a finite driving velocity. These predictions should then be compared to the results of simulations performed in the strong pinning regime, in which the crack jumps from one equilibrium position to the other. One could build on the approach recently proposed by [START_REF] Sanner | Crack-front model for adhesion of soft elastic spheres with chemical heterogeneity[END_REF], and derive a variational formulation of the problem using a third-order expansion of the elastic energy. The successive state equilibria can then be computed from a trust-region Newton conjugate-gradient algorithm. Another promising route would be to extend the numerical method of [START_REF] Rosso | Roughness at the depinning threshold for a long-range elastic string[END_REF] to second-order, but one would need first to assess the validity of Middleton's comparison theorem for the non-linear convolution integrals of [START_REF] Vasoya | A geometrically nonlinear analysis of coplanar crack propagation in some heterogeneous medium[END_REF].

• We observe in Fig. 5c and Fig. 6a that the presence of material disorder does not always result in an overall toughening of the material, in contrast with the results obtained for the semi-infinite crack (Patinet et al., 2013b;Démery et al., 2014b;[START_REF] Lebihain | Towards brittle materials with tailored fracture properties: the decisive influence of the material disorder and its microstructure[END_REF]. Indeed, for very small crack radii a 0 ≪ d, the average toughening E [T ] is zero, as the crack propagates in a single heterogeneity of uniform toughness, whose probability distribution is given by P of Eq. ( 41) and average value is K 0 Ic . For larger but still small crack radii a 0 ≲ 6d, the average toughening E [T ] is negative. The presence of material disorder leads to an overall weakening of the composite from its average toughness K 0 Ic . On the contrary, the material is reinforced by it for larger crack radii a 0 , in agreement with the results obtained for the semi-infinite crack. These variations are accurately captured by Eq. ( 62). We show in Fig. 6a the normalized average E [T ] / σ2 of the toughening as measured in our numerical simulations. All data points collapse to the master curve described by Eq. ( 62). We further observe from Fig. 6a that the average toughening E [T ] saturates in the limit of very large crack radii a 0 ≫ d. This increase of the macroscopic toughness with crack length up to a saturation value bears striking similarities with the Rcurve behavior observed in a wide variety of materials like e.g. concrete [START_REF] Bažant | R-curve modeling of rate and size effects in quasibrittle fracture[END_REF], ceramics [START_REF] Ebrahimi | Slow crack-growth behavior of alumina ceramics[END_REF], rocks [START_REF] Funatsu | Combined effects of increasing temperature and confining pressure on the fracture toughness of clay bearing rocks[END_REF], bone [START_REF] Nalla | Mechanistic aspects of fracture and r-curve behavior in human cortical bone[END_REF], compacted granular materials [START_REF] Girardot | Disk-shaped compact tension test for fracture analysis on pharmaceutical tablets[END_REF], etc.

A natural consequence of this R-curve behavior is that material disorder may promote micro-cracking below ≃ d. Indeed, the extension of pre-existing microscopic defects is facilitated by the overall decrease of the effective toughness up to a 0 = d, while their growth is further stabilized by the apparent increase of the fracture resistance. Heterogeneous materials may thus display an increased ductility at the macroscale.

• We see in Fig. 6b that the fluctuations of the toughening T (or equivalently that of the effective toughness K eff Ic ) decreases with increasing crack radius a 0 . The dependence of std [T ] in a 0 /d is well captured by our theory, as one can see from Fig. 6b, where numerical estimates of the normalized fluctuations std [T ] / σ all collapse to the master curve described by Eq. ( 63). Again, one notices that std K eff Ic = σ = std [K Ic ] for small crack radii a 0 /d ≪ 1, as the crack propagates in a single heterogeneity only. Moreover, for large crack radii a 0 /d ≫ 1, std [T ] ∝ σ(a 0 /d) -1/2 . This finite size effect may be reminiscent of having Larkin domains smaller than the front perimeter (Démery et al., 2014b). As std K eff Ic → a 0 /d→+∞ 0, it is tempting to think that our theory proves the existence of an intrinsic (deterministic) effective toughness for weakly disordered brittle composites in the limit a 0 /d → +∞. However, one must keep in mind that Eq. ( 63) only provides a first-order expansion of the standard deviation std [T ] of the toughening, so that the second-order fluctuations of the effective toughness are not necessarily zero in that limit. To access these second-order fluctuations, one should derive a third-order expansion of the variance Var [T ] of the toughening. Note that this is feasible here, even within a second-order expansion of K I , but it would involve the three-points correlation function of f , a quantity that is often poorly constrained and difficult to compute [START_REF] Jiao | A superior descriptor of random textures and its predictive capacity[END_REF]. Those calculations are out-of-scope of the present study, and are left for further works.

• The expressions of the average of the toughening E [T ] and its standard deviation std [T ] involve the functions

F and G that relate to the two-points correlation function F of the toughness fluctuations f . This was already noticed by Démery et al. (2014a) for the semi-infinite crack, and it suggests that playing on the spatial structure of the toughness field could potentially lead to a change in material toughening. An interesting route would be to investigate cases where the correlation length d z of the toughness field in the direction z is different from that d x in the direction x. One could then expand the front shape a(θ) with the disorder intensity σ, but also with a small parameter ϵ = d x /d z -1 that quantifies the toughness anisotropy. One could then observe how it modifies both the average front shape, and the overall toughening of the material. Note that such an expansion is only possible within our second-order theory, as terms of first-order in ϵ are necessarily zero (absence of disorder).

• E [T ] involves sums over the Fourier convolution kernels L(k) and H(k, k ′ ). As the latter both depend the reference SIF K 0 I generated by the loading, and its derivatives ∂K 0 I /∂a and ∂ 2 K 0 I /∂a 2 with crack advance (see Eq. ( 25)), one expects that the effective toughness K eff Ic generally depends on the loading conditions. This was already noted by Démery et al. (2014a) and [START_REF] Lebihain | Effective toughness of disordered brittle solids: A homogenization framework[END_REF] in the case of the semi-infinite crack. This dependency breaks down as soon as the structural length scales set by both the crack size a 0 and the loading conditions:

L 1 = K 0 I /(dK 0 I /da) and L 2 = K 0 I /(d 2 K 0 I /da 2 ). ( 66
)
get much larger than the heterogeneity size d. Note however that the standard deviation of the toughening std [T ] given in Eq. ( 63) is independent of the loading conditions.

Overall, our findings shed new lights on the toughening of brittle composites by heterogeneities, as well as on the size effects emerging from the finiteness of the crack geometry. However, one must keep in mind that they are only valid as long as a second-order expansion of K Ic remains valid, i.e. when the front perturbation δa is smaller than the characteristic length scale d of the toughness field (δa ≲ d). As a result, one should expect departure of the numerical results and our theory as soon as the crack enters the so-called strong pinning regime, in which crack propagation articulates as the succession of depinning instabilities (Démery et al., 2014b). But one should also be aware that taking into account size effects and non-linearities in the strong pinning regime may change the conclusions obtained by Patinet et al. (2013b) and Démery et al. (2014b).

Conclusion

This paper address the second-order variation of the local mode I stress intensity factor (SIF) arising from a coplanar perturbation of the front of a penny-shaped crack embedded in some isotropically linearly elastic infinite body. We build next on our model to determine material reinforcement arising from local fluctuations of toughness, and account for both non-linear and finite-size effects.

Section 2 was devoted to the derivation of our theoretical model, following the approach proposed by Rice (1989) and used by [START_REF] Leblond | Second-order coplanar perturbation of a semi-infinite crack in an infinite body[END_REF] for the semi-infinite crack. First, we derived the first-order variations of the fundamental kernel Z arising from the perturbation of the crack front from its reference circular configuration. Next, we determined the first-order variations of the derivative ∂K I /∂A of the mode I SIF K I with respect to the amplitude A of the perturbation, which ultimately yields our second-order formula of Eqs. (24-25) by integration. This was performed in the case where the unperturbed crack is loaded by some axisymmetric system of forces, but the general derivation was proposed in Appendix B.

In Section 3, we validated our second-order theory by comparing its predictions to the second-order expansion of an analytical formula derived by [START_REF] Irwin | Fracture[END_REF] for the SIF distribution along the front of an elliptical crack loaded by a uniform tensile stress. The validity range of our model was then assessed by comparing its output to the local SIF distributions along wavy cracks of varying mode and amplitude, which were computed using the numerical method of [START_REF] Lazarus | Brittle fracture and fatigue propagation paths of 3d plane cracks under uniform remote tensile loading[END_REF]. We showed that the predicting abilities of our second-order theory are improved from the linear model of Gao and Rice (1987b), but that errors are increasing with both the perturbation mode and amplitude.

In Section 4, we use our non-linear model to investigate crack propagation in heterogeneous media, in which the local toughness was fluctuating at random. We only considered moderate fluctuations of material toughness, so that crack propagation was stable and occurred as the continuous succession of equilibrium positions (weak pinning regime). The toughening arising from toughness fluctuations has been widely investigated in the literature for the semi-infinite coplanar crack [START_REF] Roux | Effective toughness of heterogeneous brittle materials[END_REF], even beyond the case of weak pinning, (Patinet et al., 2013b;Démery et al., 2014b;[START_REF] Lebihain | Towards brittle materials with tailored fracture properties: the decisive influence of the material disorder and its microstructure[END_REF]. However, these authors predicted a zero toughening for weak pinning, building a first-order expansion of the local stress intensity factor K I . Here, we showed that this was no more true as soon as one considers non-linearities in the expression of K I . More precisely, we showed that the effective toughness of disordered brittle media is a random variable, whose statistical descriptors, like e.g. its average and its variance, can be predicted rigorously at second order in σ using our non-linear theory. We validated our theoretical predictions on the results of ∼ 1 million numerical simulations of crack propagation in weakly disordered materials. One of our key findings is that material toughening by heterogeneities of typical size d is size-dependent. Namely, cracks of average radius ⟨a⟩ smaller than ≃ 6 d are experiencing an effective toughness K eff Ic that is lower than the average K 0 Ic of the toughness field. On the contrary, larger cracks propagates in a seemingly tougher material (K eff Ic > K 0 Ic ), in accordance with previous results obtained for the semi-infinite crack. It leads to an apparent R-curve behavior of the brittle composite at the macroscale. A natural consequence of this observation is that material disorder may promote the development of micro-cracks at the heterogeneity scale, while preventing them to extend further, leading to an overall increase in material ductility. We also proved that material toughening scales, on average, as the square σ2 of the fluctuations amplitude, and that non-linear terms in the expansion of K I do play a role in the expression of the associated prefactor. Moreover, this contribution in toughening is comparable to that emerging from the intermittent dynamics of crack front at moderate disorder intensity (Démery et al., 2014a). This drive the needs to revisit the results of Patinet et al. (2013b) and Démery et al. (2014b) in a non-linear setting, but appropriate numerical methods must be developed to validate theoretical predictions.

Our work also raises interesting perspectives to model crack propagation in a multi-physics environment. First, we provided here a second-order expansion of the stress intensity factor along quasi-circular crack fronts. However, the perturbative approach of Rice (1989) also yields the first-order variations of local aperture, a quantity that has rarely been used in the past. One could use a reasoning similar to that of [START_REF] Leblond | Second-order coplanar perturbation of a semi-infinite crack in an infinite body[END_REF], and derive its second-order expansion. Combined with that of the crack face weight function, it should provide means to investigate crack propagation in heterogeneous yet cohesive materials, whose strength σ c depends on the local crack opening [START_REF] Lebihain | Quasi-static crack front deformations in cohesive materials[END_REF]. The knowledge of the aperture of perturbed cracks also gives access to the crack volume, which may be relevant to describe 3D hydraulic fracturing [START_REF] Savitski | Propagation of a penny-shaped fluid-driven fracture in an impermeable rock: asymptotic solutions[END_REF] and its interplay with material heterogeneities. Second, one could extend our mode I calculations to mixed mode I+II+III using the perturbative framework of Favier et al. (2006a). Preliminary calculations show that it should allow rationalizing the quasi-elliptical front shapes of fluid-driven shear ruptures propagating along frictional interfaces governed by Coulomb's friction [START_REF] Sáez | Three-dimensional fluid-driven stable frictional ruptures[END_REF]. Third, the knowledge of the perturbation δa * , which has been derived in Section 2.1 to determine the first-order variation of the fundamental kernel, permits us to extend our second-order model to nearly circular connections (Gao and Rice, 1987a) at relatively small cost. One may then apply this theory to investigate adhesion along chemically heterogeneous interfaces [START_REF] Sanner | Crack-front model for adhesion of soft elastic spheres with chemical heterogeneity[END_REF] or contact between a soft material and a rough indenter [START_REF] Argatov | Controlling the adhesive pull-off force via the change of contact geometry[END_REF]. For the purpose of Open Access, a CC-BY public copyright licence has been applied by the authors to the present document and will be applied to all subsequent versions up to the Author Accepted Manuscript arising from this submission.

CRediT authorship contribution statement

Distributed under a Creative Commons Attribution -4.0 International licence Appendix A. Calculation of the first-order variations of the kernel δZ In this Appendix, our aim is to provide useful indications on the derivation of Eq. ( 15) that gives the first-order variations δZ of the kernel, starting from the general formula (8). Following Eqs. ( 4) and ( 9), the formula (8) writes as:

δZ(θ 1 , θ 2 ) = Z 0 (a 0 ; θ 1 , θ 2 )              D(M 1 , M 2 ) D(M 0 1 , M 0 2 )       -2 -1        + PV 2π 0 Z 0 (a 0 ; θ 1 , θ)Z 0 (a 0 ; θ 2 , θ) [δa(θ) -δa * (θ)] a 0 dθ + O(A 2 ) = - δa(θ 1 ) + δa(θ 2 ) 8πa 3 0 sin 2 [(θ 2 -θ 1 )/2] + 1 64π 2 a 3 0 PV 2π 0 1 sin 2 [(θ -θ 1 )/2] sin 2 [(θ -θ 2 )/2] [δa(θ) -δa * (θ)] a 0 dθ + O(A 2 ). (A.1)
Next, we use the first-order expansion of δa * of Eq. ( 13), and the decomposition of the sinus fraction given in Eq. ( 14):

δZ(θ 1 , θ 2 ) = - δa(θ 1 ) + δa(θ 2 ) 8πa 3 0 sin 2 [(θ 2 -θ 1 )/2] + 1 (8π) 2 a 3 0 sin 2 [(θ 2 -θ 1 )/2] ×                PV 2π 0 cos 2 [(θ -θ 1 )/2] sin 2 [(θ -θ 1 )/2] [δa(θ) -δa(θ 1 )] dθ + PV 2π 0 cos 2 [(θ -θ 2 )/2] sin 2 [(θ -θ 2 )/2] [δa(θ) -δa(θ 2 )] dθ - δa(θ 2 ) -δa(θ 1 ) 2 • cos[(θ 2 -θ 1 )/2] sin[(θ 2 -θ 1 )/2] •                  PV 2π 0 cos 2 [(θ -θ 1 )/2] sin 2 [(θ -θ 1 )/2] sin(θ -θ 1 ) odd integrand dθ + PV 2π 0 cos 2 [(θ -θ 2 )/2] sin 2 [(θ -θ 2 )/2] sin(θ -θ 2 ) odd integrand dθ                  - δa(θ 2 ) -δa(θ 1 ) 2 •                2π 0 cos 2 [(θ -θ 1 )/2]dθ - 2π 0 cos 2 [(θ -θ 2 )/2]dθ cancel out                + 2 cos[(θ 2 -θ 1 )/2] sin[(θ 2 -θ 1 )/2] • PV 2π 0 cos 3 [(θ -θ 1 )/2] sin[(θ -θ 1 )/2] - cos 3 [(θ -θ 2 )/2] sin[(θ -θ 2 )/2] δa(θ)dθ -                 PV 2π 0 cos 3 [(θ -θ 1 )/2] sin[(θ -θ 1 )/2] odd integrand dθ                 δa(θ 1 ) -                 PV 2π 0 cos 3 [(θ -θ 2 )/2] sin[(θ -θ 2 )/2] odd integrand dθ                 δa(θ 2 ) - δa(θ 2 ) -δa(θ 1 ) 2 • 2 cos 2 [(θ 2 -θ 1 )/2] sin 2 [(θ 2 -θ 1 )/2] •                2π 0 cos 3 [(θ -θ 1 )/2] sin[(θ -θ 1 )/2] sin(θ -θ 1 )dθ - 2π 0 cos 3 [(θ -θ 2 )/2] sin[(θ -θ 2 )/2] sin(θ -θ 2 )dθ cancel out                - δa(θ 2 ) -δa(θ 1 ) 2 • 2 cos[(θ 2 -θ 1 )/2] sin[(θ 2 -θ 1 )/2] •                 2π 0 cos 3 [(θ -θ 1 )/2] sin[(θ -θ 1 )/2] 2 sin 2 [(θ -θ 1 )/2] odd integrand dθ + 2π 0 cos 3 [(θ -θ 2 )/2] sin[(θ -θ 2 )/2] 2 sin 2 [(θ -θ 2 )/2] odd integrand dθ                 + 2π 0 2 cos 2 [(θ -θ 1 )/2] [δa(θ) -δa(θ 1 )] dθ + 2π 0 2 cos 2 [(θ -θ 2 )/2] [δa(θ) -δa(θ 2 )] dθ - δa(θ 2 ) -δa(θ 1 ) 2 • cos[(θ 2 -θ 1 )/2] sin[(θ 2 -θ 1 )/2]            2π 0 2 cos 2 [(θ -θ 1 )/2] sin(θ -θ 1 ) odd integrand dθ + 2π 0 2 cos 2 [(θ -θ 2 )/2] sin(θ -θ 2 ) odd integrand dθ            - δa(θ 2 ) -δa(θ 1 ) 2 •                2π 0 4 cos 2 [(θ -θ 1 )/2] sin 2 [(θ -θ 1 )/2]dθ - 2π 0 4 cos 2 [(θ -θ 2 )/2] sin 2 [(θ -θ 2 )/2]dθ cancel out                               + O(A 2 ).
(A.2) After simplifications of the integrals that are zero due to the oddness of the integrand, and those that cancel out, one has:

δZ(θ 1 , θ 2 ) = - δa(θ 1 ) + δa(θ 2 ) 8πa 3 0 sin 2 [(θ 2 -θ 1 )/2] + 1 (8π) 2 a 3 0 sin 2 [(θ 2 -θ 1 )/2] × PV 2π 0 cos 2 [(θ -θ 1 )/2] sin 2 [(θ -θ 1 )/2] [δa(θ) -δa(θ 1 )] dθ + PV 2π 0 cos 2 [(θ -θ 2 )/2] sin 2 [(θ -θ 2 )/2] [δa(θ) -δa(θ 2 )] dθ + 2 cos[(θ 2 -θ 1 )/2] sin[(θ 2 -θ 1 )/2] • PV 2π 0 cos 3 [(θ -θ 1 )/2] sin[(θ -θ 1 )/2] - cos 3 [(θ -θ 2 )/2] sin[(θ -θ 2 )/2] δa(θ)dθ + 2π 0 2 cos 2 θ -θ 1 2 + cos 2 θ -θ 2 2 δa(θ)dθ -2π [δa(θ 1 ) + δa(θ 2 )] + O(A 2 ), (A.3)
from which one retrieves Eq. ( 15). Note that the simplifications in Eq. (A.2) arise from the careful choice of δa * in Eq. ( 10) and that of θ * in Eq. ( 11), which allows for the alternate expressions of δa * in Eq. ( 13) that pair well with the decomposition of Eq. ( 14).

Appendix B. Calculation of the second-order variations of the SIF δK

In this Appendix, we provide useful indications on how to calculate the second-order K 2 I perturbations of the mode I stress intensity factor arising from perturbations δa = Aϕ of the circular front of radius a 0 . For the sake of completeness, we do it here in the general case where the reference SIF K 0 I , applied to the unperturbed circular crack, depends on θ. The calculations presented in the main text only address cases where K 0 I is axisymmetric.

As established in Eq. ( 18), the SIF variations arising from an additional infinitesimal perturbation δAϕ reads:

K I (a 0 , (A + δA)ϕ ; θ 1 ) -K I (a 0 , Aϕ ; θ 1 ) = K I (a 0 + δAϕ(θ 1 ), Aϕ ; a 0 ; θ 1 ) -K I (a 0 , Aϕ ; θ 1 )

+ PV 2π 0
Z(a 0 + δAϕ(θ 1 ), Aϕ ; θ 1 , θ)K I (a 0 + δAϕ(θ 1 ), Aϕ ; θ) δAϕ(θ) -δAϕ(θ 1 ) ds(θ) + O(δA 2 ), (B.1) where s is the curvilinear abscissa along the perturbed front F ′ , so that ds reduces here to:

ds(θ) = a 0 + δAϕ(θ 1 ) + Aϕ(θ) + O(A 2 ), (B.2)
Using Eq. ( 2), one has: Z(a 0 + δAϕ(θ 1 ), Aϕ ; θ 1 , θ)K I (a 0 , Aϕ ; θ) δAϕ(θ) -δAϕ(θ 1 ) ds(θ)

K I (a 0 + δAϕ(θ 1 ), Aϕ ; a 0 ; θ 1 ) -K I (a 0 , Aϕ ; θ 1 ) = δA ×        ∂K 0 I ∂a (a 0 , θ 1 )ϕ(θ 1 ) + A ∂ 2 K 0 I ∂a 2 (a 0 , θ 1 )ϕ(θ 1 ) 2 + A PV 2π 0 ∂Z 0 ∂a (a 0 ; θ 1 , θ)K 0 I (a 0 , θ) ϕ(θ) -ϕ(θ 1 ) ϕ(θ 1 )a 0 dθ + A PV
= PV 2π 0 Z(a 0 , Aϕ ; θ 1 , θ)K I (A; a 0 ; θ) δAϕ(θ) -δAϕ(θ 1 ) ds(θ) + O(δA 2 ) = δA × PV 2π 0 Z 0 (a 0 ; θ 1 , θ)K 0 I (a 0 , θ) ϕ(θ) -ϕ(θ 1 ) a 0 dθ + A PV 2π 0 Z 1 (a 0 , ϕ; θ 1 , θ)K 0 I (a 0 , θ) ϕ(θ) -ϕ(θ 1 ) a 0 dθ + A PV 2π 0 Z 0 (a 0 ; θ 1 , θ) ∂K 0 I ∂a (a 0 , θ) ϕ(θ) -ϕ(θ 1 ) ϕ(θ)a 0 dθ + A PV 2π 0 2π 0 Z 0 (a 0 ; θ 1 , θ)Z 0 (a 0 ; θ, θ ′ )K 0 I (a 0 , θ ′ ) ϕ(θ) -ϕ(θ 1 ) ϕ(θ ′ ) -ϕ(θ) a 2 0 dθdθ ′ +A PV 2π 0 Z 0 (a 0 ; θ 1 , θ)K 0 I (a 0 , θ) ϕ(θ) -ϕ(θ 1 ) ϕ(θ)dθ + O(A 2 ) + O(δA 2 ). (B.4)
Identifying the terms in δA in the expansions of Eqs. ( 17) and (B.1), one finds:

K 1 I (a 0 , [K 0 I ], [ϕ]; θ 1 ) = ∂K 0 I ∂a (a 0 , θ 1 )ϕ(θ 1 ) + PV 2π 0
Z 0 (a 0 ; θ 1 , θ)K 0 I (a 0 , θ) ϕ(θ) -ϕ(θ 1 ) a 0 dθ, (B.5) which is Eq. ( 2) of Gao and Rice (1987b). The terms in δA × A yield:

2K 2 I (a 0 , [K 0 I ], [ϕ]; θ 1 ) = ∂ 2 K 0 I ∂a 2 (a 0 , θ 1 )ϕ(θ 1 ) 2 + PV 2π 0 Z 0 (a 0 ; θ 1 , θ) ∂K 0 I ∂a (a 0 , θ) ϕ(θ) 2 -ϕ(θ 1 ) 2 a 0 dθ + PV 2π 0 Z 0 (a 0 ; θ 1 , θ)K 0 I (a 0 , θ) ϕ(θ) -ϕ(θ 1 ) 2 dθ + PV 2π 0 2π 0 Z 0 (a 0 ; θ 1 , θ)Z 0 (a 0 ; θ, θ ′ )K 0 I (a 0 , θ ′ ) ϕ(θ) -ϕ(θ 1 ) ϕ(θ) -ϕ(θ ′ ) a 2 0 dθdθ ′ + PV 2π 0 Z 1 (a 0 , ϕ ; θ 1 , θ)K 0 I (a 0 , θ) ϕ(θ) -ϕ(θ 1 ) a 0 dθ, (B.6)
which is Eq. ( 20). Here we used ∂Z 0 /∂a = -2Z 0 /a 0 from Eq. ( 4).

We need now to calculate Eqs. (B.5) and (B.6) by resorting to Fourier series. The calculations are heavy, and we will show them here only for K 1 I . Eq. (B.5) yields:

K 1 I (a 0 , [K 0 I ], [ϕ]; θ 1 ) = ∂K 0 I ∂a (a 0 , θ 1 )ϕ(θ 1 ) + PV 2π 0 Z 0 (a 0 ; θ 1 , θ)K 0 I (a 0 , θ) ϕ(θ) -ϕ(θ 1 ) a 0 dθ = k,k ′       ∂K 0 I ∂a       k (a 0 ) φk ′ e i(k+k ′ )θ1 + k,k ′ K 0 Ik (a 0 ) φk ′ 1 8a 0 PV 2π 0 e ik ′ (θ-θ1) -1 sin 2 [(θ -θ 1 )/2] e ik(θ-θ1) dθ e i(k+k ′ )θ1 = k,k ′                ∂K 0 I ∂a       k (a 0 ) - |k + k ′ | -|k| 2a 0 K 0 Ik (a 0 )          φk ′ e i(k+k ′ )θ1 , (B.7)
where φk , K 0 Ik (a 0 ), and which is calculated using the residue theorem. The calculations for K 2 I are similar, although more complicated as they involve triple sums. They can be simplified by building on the following expression of the first-order kernel variations Z 1 :

      ∂K 0 I ∂a       k ( a
Z 1 (a 0 , ϕ ; θ 1 , θ 2 ) = - 1 16πa 2 0 sin 2 [(θ 2 -θ 1 )/2] ×        k φk e ikθ2 + e ikθ1 |k| + 2 - ∆ |k|1 4 + i cos [(θ 2 -θ 1 )/2] sin [(θ 2 -θ 1 )/2] k φk e ikθ2 -e ikθ1 1 -∆ k0 - ∆ |k|1 4 sign(k)        , (B.9)
as well as the following integrals:

PV 2π 0 2 cos 3 (u/2) sin(u/2) e iku du = 4iπ sign(k) 1 - 1 4 ∆ |k|1 -∆ k0 PV 2π 0 cos(u/2) sin 3 (u/2) (e iku -1)(e ik ′ u -1)du = -4iπ k(|k + k ′ | -|k|) + k ′ (|k + k ′ | -|k ′ |) , (B.10)
where ∆ i j is the Kronecker symbol that is equal to 1 when i = j and 0 otherwise, and sign(k) is the sign of the integer k. Using Eqs. (B.8), (B.9), and (B.10), Eq. (B.6) yields: (B.11) where the convolution kernel K reads:

K 2 I (a 0 , [K 0 I ], [ϕ]; θ 1 ) = k,k ′ ,k ′′ H K 0 I ; k, k ′ , k ′′ φk ′ φk ′′ e i(k+k ′ +k ′′ )θ 1 ,
H K 0 I ; k, k ′ , k ′′ = 1 2       ∂ 2 K 0 I ∂a 2       k (a 0 ) - 1 4a 0       ∂K 0 I ∂a       k (a 0 ) |k + k ′ + k ′′ | -|k| - 1 4a 2 0 K 0 Ik (a 0 ) |k + k ′ + k ′′ | -|k + k ′ | -|k + k ′′ | + |k| - 1 16a 2 0 K 0 Ik (a 0 ) (k + k ′ ) 2 + (k + k ′′ ) 2 + 1 16a 2 0 K 0 Ik (a 0 ) |k + k ′ | + |k + k ′′ | |k + k ′ + k ′′ | + |k| - 1 8a 2 0 K 0 Ik (a 0 )|k||k + k ′ + k ′′ | + 1 4a 2 0 K 0 Ik (a 0 ) |k + k ′ + k ′′ | -|k| + 1 16a 2 0 K 0 Ik (a 0 ) |k ′ ||k + k ′′ | + |k ′′ ||k + k ′ | - 1 16a 2 0 K 0 Ik (a 0 )|k| |k ′ | + |k ′′ | - 1 16a 2 0 K 0 Ik (a 0 )k |k + k ′ + k ′′ | -|k + k ′ | -|k + k ′′ | + |k| sign(k ′ ) + sign(k ′′ ) - 1 16a 2 0 K 0 Ik (a 0 ) |k + k ′ + k ′′ |(|k ′ | + |k ′′ |) -|k ′ ||k + k ′ | -|k ′′ ||k + k ′′ | sign(k ′ k ′′ ). (B.
12) The set of equations (B.11-12) is validated next in Appendix C on elliptical cracks loaded by a non-axisymmetric loading. Considering the limiting case of a reference SIF K 0 I independent of θ, one has:

H K 0 I ; k, k ′ , k ′′ = 1 2 d 2 K I da 2 (a 0 ) - 1 a 0 dK I da (a 0 ) |k + k ′ | 4 + 1 a 2 0 K 0 I (a 0 ) |k| + |k ′ | 4 + 1 8a 2 0 K 0 I (a 0 )|k||k ′ | + 1 16a 2 0 K 0 I (a 0 ) |k + k ′ |(|k| + |k ′ |) -k 2 -k ′2 1 -sign(kk ′ ) = 1 2 d 2 K I da 2 (a 0 ) - 1 a 0 dK I da (a 0 ) |k + k ′ | 4 + 1 a 2 0 K 0 I (a 0 ) |k| + |k ′ | 4 + 1 16a 2 0 K 0 I (a 0 )R(k, k ′ ), (B.13)
from which one finds back Eq. ( 20).

The general expression (B.11-12) of the second-order variations of stress intensity factor for a perturbed crack front given prove rather complex, and its scope of application seems limited. Indeed, from a theoretical point of view, the influence of a non-uniformity of macroscopic loading on mode I crack propagation is already of secondorder within the linear expansion of Gao and Rice (1987b) (see Eq. ( 2)). Moreover, the triple sum over the Fourier coefficients related to K 0 I , its derivatives with crack advance (∂ n K 0 I /∂a n ), and the perturbation shape ϕ makes it unfit for numerical calculations, as estimating second-order contributions would take a computation time in O(n 3 ln n), where n is the number of discretization points along the crack front. This is already larger than more conventional numerical approaches, like e.g. the boundary element method, that can simulate crack propagation at all orders in "only" O(n 2 ln n).

However, we did not perform here the rigorous calculation for a general K 0 I for the sake of completeness only. Indeed, for shear ruptures, the mode II and mode III stress intensity factors applied to a circular crack depends on θ, even when the applied stress is axisymmetric [START_REF] Gao | Nearly circular shear mode cracks[END_REF]. Propagating internal shear cracks often adopt quasielliptical front shapes [START_REF] Sáez | Three-dimensional fluid-driven stable frictional ruptures[END_REF] that may be described with a perturbative approach to LEFM. We provide thus in this Appendix the basic steps needed to derive an equivalent second-order theory for shear cracks. Future work will be devoted to the construction of this model and its application to the propagation of e.g. fluid-driven shear ruptures.

Appendix C. Validation of the general formula (B.11) on elliptical cracks loaded by non-axisymmetric loadings

This Appendix aims at validating the general expression of the second-order variations of stress intensity factor along a perturbed crack front obtained in Eqs. (B.11-12). We provided in Section 3.1 a validation of Eq. ( 20) or equivalently (B.13) by comparing our second-order predictions to the second-order expansions of the exact stress intensity factor distribution along elliptical crack fronts loaded by a uniform remote stress σ. Fortunately, analytical predictions for the SIF distribution along elliptical crack fronts are available in the literature for more complex loading situations [START_REF] Shah | Stress intensity factor for an elliptical crack under arbitrary normal loading[END_REF][START_REF] Kassir | Three-dimensional crack problems : a new selection of crack solutions in three-dimensional elasticity[END_REF][START_REF] Atroshchenko | Stress intensity factor for an embedded elliptical crack under arbitrary normal loading[END_REF]. Here, a comparison between our model's predictions and second-order expansions of exact analytical solutions is provided for three non-axisymmetric loading cases.

As in Section 3.1, we consider a crack whose front F describes an ellipsis centered in O, of short semi-axis c aligned with e z , and long semi-axis b aligned with e x . The crack front F follows the contour parametrized by the variable φ: (z = b cos φ, x = c sin φ). Its position a(θ) in polar coordinates is given by Eq. ( 26), so that: a(θ) cos θ = b cos φ, and a(θ) sin θ = c sin φ.

(C.1)

We follow here the reasoning of Gao and Rice (1987b), and take as small expansion parameter ϵ = 1c/b, and a reference radius a 0 equal to the front extension at the point θ 1 of SIF evaluation (i.e. a 0 = a(θ 1 )). In that case, the perturbation δa reads:

δa(θ) = a 0 1 + ϵA 1 (θ) + ϵ 2 A 2 (θ) , (C.2) where A 1 (θ) = - 1 2 cos(2θ 1 ) + 1 2 cos(2θ),
and

A 2 (θ) = 1 8 - 1 4 cos(2θ 1 ) - 1 16 cos(4θ 1 ) + 1 4 - 1 4 cos(2θ 1 ) cos(2θ) + 3 16 cos(4θ). (C.3)
Next, we introduce k I , the crack face weight function of a circular crack. k I (a, θ; r, θ ′ ) corresponds to the mode I SIF generated at the point θ along the circular front of radius a by a pair of unitary symmetric tensile forces acting on the point interior to the crack indexed by (r, θ ′ ). It reads [START_REF] Tada | The Stress Analysis of Cracks Handbook[END_REF]:

k I (a, θ; r, θ ′ ) = √ a 2 -r 2 π 3/2 a 1/2 a 2 + r 2 -2ar cos(θ -θ ′ ) . (C.4)
Considering that the crack is loaded by a distribution of tensile stress

σ(r, θ) = k σk (r)e ikθ ,
acting on the failure plane y = 0 of the intact body, the mode I SIF generated along the front of a circular crack of radius a reads:

K 0 I (a, θ) = 2π 0 a 0 σ(r, θ ′ )k I (a, θ; r, θ ′ )rdrdθ ′ = k 2 π √ πa 1 0 u |k|+1 √ 1 -u 2 σk (au)du e ikθ .
(C.5)

We have at our disposal all the basic ingredients to use our second-order model of Eqs. (B.7) and (B.11).

Linear stress field dependent on z.

We first consider the case where the elliptical crack is loaded by a linear stress field dependent on z:

σ(z, x) = σ 0 z b ⇒ σ(r, θ) = σ 0 r b cos(θ). (C.6)
Using Eq. (C.5) the associated SIF distribution along the front of a circular crack of radius a 0 reads:

K 0 I (a, θ) = 4 3 √ π a 3/2 b σ 0 cos(θ). (C.7)
The exact SIF distribution along the elliptical crack front has been calculated by [START_REF] Shah | Stress intensity factor for an elliptical crack under arbitrary normal loading[END_REF]) and [START_REF] Atroshchenko | Stress intensity factor for an embedded elliptical crack under arbitrary normal loading[END_REF]. It reads:

K e I (φ) = - √ πcσ 0 cos(φ) κ 2 (1 -2κ 2 )E(κ) -κ ′2 K(κ) (c/b) 4 cos 2 (φ) + sin 2 (φ) 1/4 ⇒ K e I (θ) = - √ πc 3/2 b σ 0 cos(θ) κ 2 (1 -2κ 2 )E(κ) -κ ′2 K(κ) (c/b) 4 cos 2 (θ) + sin 2 (θ) 1/4 (c/b) 2 cos 2 (θ) + sin 2 (θ) 3/4 . (C.8)
where κ = 1c 2 /b 2 1/2 and κ ′ = c/b. One may then use Eqs. (C.7) and (C.8) to expand K e I /K 0 I to second order in ϵ at ϵ = 0 and θ = θ 1 . One finds:

K e I (A 1 )/K 0 I (a 0 , θ 1 ) = - 3π 4 κ 2 (1 -2κ 2 )E(κ) -κ ′2 K(κ) (c/b) 4 cos 2 (θ 1 ) + sin 2 (θ 1 ) 1/4 = 1 + ϵ 1 4 - 1 2 cos(2θ 1 ) + ϵ 2 5 16 - 3 8 cos(2θ 1 ) - 3 16 cos(4θ 1 ) + O(ϵ 3 ).
(C.9)

Our aim here is to retrieve Eq. (C.9) from our general second-order theory of Eqs. (B.7) and (B.11). The second-order estimate K p I writes as:

K p I (θ 1 ) = K 0 I (a 0 , θ 1 ) + ϵK 1 I (a 0 , [K 0 I ], [A 1 ]; θ 1 ) + ϵ 2 K 1 I (a 0 , [K 0 I ], [A 2 ]; θ 1 ) + K 2 I (a 0 , [K 0 I ], [A 1 ]; θ 1 ) . (C.10)
Combining Eqs. (B.7) and (C.3), one has

K 1 I (a 0 , [K 0 I ], [A 1 ]; θ 1 ) = K 0 I (a 0 , θ 1 ) × 1 4 - 1 2 cos(2θ 1 ) ,
and

K 1 I (a 0 , [K 0 I ], [A 2 ]; θ 1 ) = K 0 I (a 0 , θ 1 ) × 5 32 - 3 16 cos(2θ 1 ) - 1 4 cos(4θ 1 ) . (C.11)
Similarly, the combination of Eqs. (B.11) and (C.3) yields:

K 2 I (a 0 , [K 0 I ], [A 1 ]; θ 1 ) = K 0 I (a 0 , θ 1 ) × 5 32 - 3 16 cos(2θ 1 ) + 1 16 cos(4θ 1 ) . (C.12)
We observe that the two expansions coincide at θ = θ 1 .

Linear stress field dependent on x.

We consider next the case where the elliptical crack is loaded by a linear stress field dependent on x:

σ(z, x) = σ 0 x c ⇒ σ(r, θ) = σ 0 r c sin(θ). (C.13)
As the procedure is similar to the previous case, we should be brief here. One has:

K 0 I (a, θ) = 4 3 √ π a 3/2 c σ 0 sin(θ). (C.14)
The exact SIF distribution along the elliptical crack front has been obtained by [START_REF] Shah | Stress intensity factor for an elliptical crack under arbitrary normal loading[END_REF]) and [START_REF] Atroshchenko | Stress intensity factor for an embedded elliptical crack under arbitrary normal loading[END_REF]. It reads:

K e I (φ) = √ πcσ 0 sin(φ) κ 2 (1 + κ 2 )E(κ) -κ ′2 K(κ) (c/b) 4 cos 2 (φ) + sin 2 (φ) 1/4 ⇒ K e I (θ) = √ πcσ 0 sin(θ) κ 2 (1 + κ 2 )E(κ) -κ ′2 K(κ) (c/b) 4 cos 2 (θ) + sin 2 (θ) 1/4 (c/b) 2 cos 2 (θ) + sin 2 (θ) 3/4 . (C.15)
The expansion of K e I /K 0 I at θ = θ 1 to second order in ϵ at ϵ = 0 and θ = θ 1 yields:

K e I (A 1 )/K 0 I (a 0 , θ 1 ) = 3π 4 κ 2 (1 + κ 2 )E(κ) -κ ′2 K(κ) (c/b) 4 cos 2 (θ 1 ) + sin 2 (θ 1 ) 1/4 = 1 + ϵ - 1 4 - 1 2 cos(2θ 1 ) + ϵ 2 1 16 - 1 8 cos(2θ 1 ) - 3 16 cos(4θ 1 ) + O(ϵ 3 ). (C.16)
Our goal is to retrieve again Eq. (C.16) from Eq. (C.10). One has:

K 1 I (a 0 , [K 0 I ], [A 1 ]; θ 1 ) = K 0 I (a 0 , θ 1 ) × - 1 4 - 1 2 cos(2θ 1 ) , and K 1 I (a 0 , [K 0 I ], [A 2 ]; θ 1 ) = K 0 I (a 0 , θ 1 ) × - 3 32 - 5 16 cos(2θ 1 ) - 1 4 cos(4θ 1 ) ,
(C.17) and:

K 2 I (a 0 , [K 0 I ], [A 1 ]; θ 1 ) = K 0 I (a 0 , θ 1 ) × 5 32 + 3 16 cos(2θ 1 ) + 1 16 cos(4θ 1 ) . (C.18)
The two expansions are found equal once again.

Quadratic stress field odd in z and x.

We conclude this Appendix with a last validation case, in which the stress field is quadratic and odd in both z and x: For which one has:

K 0 I (a, θ) = 8 15
√ π a 5/2 bc σ 0 sin(2θ).

(C.20)

The exact SIF distribution along the elliptical crack front has been obtained by [START_REF] Kassir | Three-dimensional crack problems : a new selection of crack solutions in three-dimensional elasticity[END_REF] and [START_REF] Atroshchenko | Stress intensity factor for an embedded elliptical crack under arbitrary normal loading[END_REF]. It reads: The two expansions are found equal once more, although the shape of the stress field is more complex.

Appendix D. Ensemble averages of crack front propagating in weakly disordered materials

This Appendix is dedicated to the statistical analysis of equation ( 61) defining the effective toughness of a random medium. In the following, we consider a statistical ensemble Ω containing all the possible realizations of the heterogeneous medium considered in Section 4. To each realization, we associate a specific real number ω. We note p : ω → p (ω) the probability function associated to the ensemble Ω. Thus, the probability that the variable ω ′ lies in some neighborhood of ω of measure dω is p(ω)dω. The mathematical expectation E [u(x 1 , x 2 )] of any spatial observable u : (x 1 , x 2 ) → u(x 1 , x 2 ) is defined as: E [u(x 1 , x 2 )] = Ω u(x 1 , x 2 ; ω)p(ω)dω (D.1) Using these notations, the ensemble averaging of Eq. (61) yields: E K eff Ic /K 0 Ic = 1 + σL(0)E Â1 0 (a 0 ) + σ2 L(0)E Â2 0 (a 0 ) + H(0, 0)E Â1 0 (a 0 ) 2 , (D.2) so that its determination requires the prior calculation of E Â1 0 (a 0 ) , E Â2 0 (a 0 ) , and E Â1 0 (a 0 ) 2 . We start here with Â1 0 (a 0 ), for which the calculation is the simplest and illustrates well the reasoning behind the ulterior ones. where we used the expression of Â1 k (a 0 ) given in Eq. ( 54). Even if only the value of E Â1 0 (a 0 ) to estimate the average value of the effective toughness, our calculations show that:

E [A 1 (a 0 , θ)] = E        k Â1 k (a 0 )e ikθ        = k E Â1
k (a 0 ) e ikθ = 0, (D.4) meaning that the first-order contribution to the crack front deformations are on average zero.

We proceed similarly for A 2 , and start by averaging Eq. ( 56):

E Â2 k (a 0 ) = 1 L(k) k ′         a 0 L(k -k ′ ) E         ∂ f ∂r k ′ (a 0 ) fk-k ′ (a 0 )         - H(k ′ , k -k ′ ) L(k ′ )L(k -k ′ ) E fk ′ (a 0 ) fk-k ′ (a 0 )         (D.5)
Let us first concentrate on:

E fk (a 0 ) fk ′ (a 0 ) = 1 4π 2 2π 0 2π 0 E f (a 0 , θ) f (a 0 , θ ′ ) e -i(kθ+k ′ θ ′ ) dθdθ ′ (D.6)

To calculate Eq. (D.6), it is convenient to switch to Cartesian coordinates. The pairs (r, θ) and (r ′ , θ ′ ) are then equivalent to(z = a 0 cos θ, x = a 0 sin θ) and (z ′ = a 0 cos θ ′ , x ′ = a 0 sin θ ′ ). With these notations, one has:

E f (a 0 , θ) f (a 0 , θ ′ ) = E f * (a 0 cos θ, a 0 sin θ) f * (a 0 cos θ ′ , a 0 sin θ ′ e -i(kθ+k ′ θ ′ ) = F z zz ′ F x xx ′ = e -4(a 0 /d) 2 sin 2 [(θ-θ ′ )/2] .

(D.7)

where we used the expressions of F z and F x given in Eq. ( 45), as well as the following identities:

       z -z ′ = -2a 0 sin[(θ -θ ′ )/2] sin[(θ + θ ′ )/2] x -x ′ = 2a 0 sin[(θ -θ ′ )/2] cos[(θ + θ ′ )/2].
Injecting Eq. (D.7) into (D.6), one finds:

E fk (a 0 ) fk ′ (a 0 ) = 1 4π 2 2π 0 2π 0 e -4(a 0 /d) 2 sin 2 [(θ-θ ′ )/2] e -ik(θ-θ ′ ) e -i(k+k ′ )θ ′ dθdθ ′ = ∆ -kk ′ Fk (a 0 /d), (D.8)

where Fk is the k-th Fourier coefficient of the 2π-periodic function:

F(a 0 /d, θ) = e -4(a 0 /d) 2 sin 2 (θ/2) , (D.9) which relates to the correlation shape F of Eq. ( 45). We must now calculate the second average term in the sum of Eq.(D.5). i(m cos θ + n sin θ) F z (m) F x (n)e i(m(z-z ′ )+n(x-x ′ )) dmdndm ′ dn ′ = F ′ z (zz ′ )F x (xx ′ ) cos θ + F z (zz ′ )F ′ x (xx ′ ) sin θ = -4a 0 /d 2 sin 2 [(θ -θ ′ )/2]e -4(a 0 /d) 2 sin 2 [(θ-θ ′ )/2] (D.15) where we used the following identities: -4(a 0 /d) 2 sin 2 [(θ -θ ′ )/2]e -4(a 0 /d) 2 sin 2 [(θ-θ ′ )/2] e -ik(θ-θ ′ ) e -i(k+k ′ )θ ′ dθdθ ′ = ∆ -kk ′ Ĝk (a 0 /d)/a 0 , (D.16)

       F ′ z (z -z ′ )F x (x -x ′ ) = 4a 0 /d 2 sin[(θ -θ ′ )/2] sin[(θ + θ ′ )/2]e -4(a 0 /d) 2 sin 2 [(θ-θ ′ )/2] F z (z -z ′ )F ′ x (x -x ′ ) = -
where Ĝk is the k-th Fourier coefficient of the 2π-periodic function:

G(a 0 /d, θ) = -4(a 0 /d) 2 sin 2 (θ/2)e -4(a 0 /d) 2 sin 2 (θ/2) (D.17)

As a result, all E Â2 k (a 0 ) are zero, except for k = 0, for which one gets:

E Â2 0 (a 0 ) = 1 L(0) k ′ 1 L(k ′ )
Ĝk ′ (a 0 /d) -H(k ′ , -k ′ ) L(k ′ ) 2 Fk ′ (a 0 /d) (D.18) This means that the average second-order contribution to the front deformations is a simple dilation of the reference circular front of radius a 0 by a quantity:

E [A 2 (a 0 , θ)] = 1 L(0) k ′ 1 L(k ′ ) Ĝk ′ (a 0 /d) - H(k ′ , -k ′ ) L(k ′ ) 2 Fk ′ (a 0 /d) (D.19)
We only need to calculate the third and last missing term of Eq. (D.2). Using Eq. (D.8), one has: (D.20) so that the average of the effective toughness finally reads:

E Â1 0 (a 0 ) 2 = E f0 (a 0 ) 2 L(0) 2 = F0 (a 0 /d) L(0) 2 ,
E K eff Ic /K 0 Ic = 1 + σ2        k ′ 1 L(k ′ ) Ĝk ′ (a 0 /d) - H(k ′ , -k ′ ) L(k ′

Figure 3 :

 3 Figure 3: Normalized stress intensity factor variations δK I /ϵK 0 I along a crack front perturbed by δa(θ)/a 0 = ϵ cos(kθ) for (a) k = 2 and (b) k = 5, and ϵ = 0.01 (left panel), ϵ = 0.1 (middle panel), and ϵ = 0.25 (right panel). The SIF variations K num I computed from Lazarus (2003)'s numerical model (in black solid line) are compared to the first-order estimate ofGao and Rice (1987b) (in dotted orange line) and our second-order prediction (in dashed red line). The errors on the SIF variations between the numerical simulations (in solid line) and the first-order (in dotted line) and second-order (in dashed line) models are compared using (c) L 0 (maximum) or (d) L 2 (Euclidean) norms. They increase with both the perturbation mode k and amplitude ϵ. Deviations of the numerical results from the first-and second-order predictions can be directly observed on (e) the average ⟨K I ⟩ of the stress-intensity factor along the crack front, or (f) its amplitude K max

Figure 4 :

 4 Figure 4: (a) Front positions (in solid lines) of a crack propagating in a disordered toughness field K Ic under the action of a point force P applied at the crack center. The fronts are plotted for a constant load increment. (b) The values of K Ic are uniformly distributed around a mean value K 0 Ic

Figure 5 :

 5 Figure 5: (a) At constant disorder intensity σ/K 0 Ic = 10 -1 , the simulated values of effective toughness K Ic are distributed between the two extremal values K min Ic and K max Ic of the local toughness distribution, and its distribution depends on the ratio between the reference crack radius a 0 and the heterogeneity size d. (b) Similarly, at constant size ratio a 0 /d = 0.1, the of K eff Ic depends on the disorder intensity σ/K 0 Ic (the distribution plotted in bold is the same to that shown in 5a). As a result, (c) the average E [T ] and (d) the standard deviation std [T ] of the toughening T = K effIc /K 0 Ic -1 depends on both σ/K 0 Ic and a 0 /d. In particular, E [T ] scales as (σ/K 0 Ic ) 2 and its sign depends on a 0 /d, while std (T ) is linear in σ/K 0 Ic .

Figure 6 :

 6 Figure 6: (a) When normalized by the square of the disorder intensity (σ/K 0 Ic ) 2 , the values of the average E [T ] of the toughening showed in Figure 5c all fall down onto a master curve predicted by Eq. (62). (b) Similarly, when normalized by (σ/K 0 Ic ), the values of the standard deviation std [T ] presented in Figure 5d follow the master curve of Eq. (63) controlled by a 0 /d.

Mathias Lebihain :

 Lebihain Conceptualization, Investigation, Formal analysis, Numerical simulations, Software, Visualization, Writing -Original Draft. Manish Vasoya: Investigation, Numerical simulations, Writing -Review & Editing. Véronique Lazarus: Investigation, Software, Writing -Review & Editing.

  a 0 ; θ 1 , θ)K 0 I (a 0 , θ) ϕ(θ) -ϕ(θ 1 ) ϕ(θ 1 )dθ + O(A 2 ) + O(δA 2 ).(B.3)The calculation of the second term in the right-hand side of Eq. (B.1) is also straightforward. From Eqs. (16) and (B.2), one has: PV 2π 0

  0 ) are the k-th Fourier coefficient of ϕ, K 0 I (a 0 , •), and ∂K I ∂a (a 0 , •) respectively. The terms into brack in Eq. (B.7) reduces to L(k) as given in Eq. (25a) when K 0 I is independent of θ. In deriving (B.7), we used the following integral:

  κ ′4 -κ ′2 + 1)E(κ) -(κ ′4 + κ ′2 )K(κ) (c/b) 4 cos 2 (φ) + sin 2 ′4 -κ ′2 + 1)E(κ) -(κ ′4 + κ ′2 )Kθ = θ 1 to second order in ϵ at ϵ = 0 and θ = θ 1 then yields: ′4 -κ ′2 + 1)E(κ) -(κ ′4 + κ ′2 )K(κ) (c/b) 4 cos 2 (θ 1 ) + sin 2 1 ) + O(ϵ 3 ). (C.22)We need now to compute the three different terms of our second-order equation (C.10). The terms arising from the linearization of K I yields:

  , θ) f (a 0 , θ ′ ) e -i(kθ+k ′ θ ′ ) dθdθ ′ (D.10)The calculation of equation (D.10) proves more difficult than that of Eq. (D.7). In the following, use will be made of Fourier transforms. The definition adopted here for the single Fourier transform g(p) of an arbitrary function g(x) is:g(p) = +∞ -∞ g(x)e -ipx dx ⇐⇒ g(x) = 1 2π +∞ -∞ g(p)e ipx d p. (D.11)and the double (z, x) Fourier transform h(m, n) of an arbitrary function h(z, x) reads:x; ω) cos θ + ∂ f * ∂x (z, x; ω) sin θ f * (z ′ , x ′ ; ω)p(ω,m ′ ,n ′ i(m cos θ + n sin θ)E f * (m, n) f * (m ′ , n ′ ) e i(mz+nx+m ′ z ′ +n ′ x ′ ) dmdndm ′ dn ′ .(D.13) As shown byFavier et al. (2006b) in their Appendix B, one has:E f * (m, n) f * (m ′ , n ′ ) = (2π) 2 δ(m + m ′ )δ(n + n ′ ) F z (m) F x (n), (D.14)where δ is the Dirac function, and F z and F x are the Fourier transforms of the correlation functions F z and F x . Injecting Eq. (D.14) into (D.13), one finds:

  4a 0 /d 2 sin[(θ -θ ′ )/2] sin[(θ + θ ′ )/2]e -4(a 0 /d) 2 sin 2 [(θ-θ ′ )/2] , θ = cos[(θ -θ ′ )/2] cos[(θ + θ ′ )/2]sin[(θ -θ ′ )/2] sin[(θ + θ ′ )/2] sin θ = sin[(θ -θ ′ )/2] cos[(θ + θ ′ )/2] + cos[(θ -θ ′ )/2] sin[(θ + θ ′ )/2].The combination of Eqs. (D.10) and (D.15) finally yields:

  Ic /K 0 Ic = E (K eff Ic /K 0 Ic ) 2 -E K eff Ic /K 0 Ic 2 = σ2 L(0) 2 E Â1 0 (a 0 ) 2 = σ2 F0 (a 0 /d) + O( σ3 ).(D.22)

Note the subtle difference of vocabulary: fundamental kernel standing for the dimensionless quantity W and kernel for Z appearing as the kernel of the integral in Eq.

Note that one must ensure that the secondary perturbation is zero at θ = θ 1 , otherwise the transformation induces a phase-shift that is first-order in δA. This was overlooked by[START_REF] Leblond | Second-order coplanar perturbation of a semi-infinite crack in an infinite body[END_REF] and[START_REF] Vasoya | A geometrically nonlinear analysis of coplanar crack propagation in some heterogeneous medium[END_REF], but their results remain valid as these first-order corrections are zero when K 0 I is independent of the position along the front.

Acknowledgments

The authors gratefully thank Jean-Baptiste Leblond for stimulating discussions on the second-order derivation, Antoine Sanner and Lars Pastewka for fruitful discussions on the extension of our approach to other physical systems and on variational formulations, Sébastien Brisard for suggesting variance reduction methods. We also thank the anonymous reviewer for their constructive and insightful comments. ML thanks Antoine Sanner for a critical reading of the manuscript.

Online access to data

The Python script used to generate the figures and the associated data are available on the French open data plateform recherche.data.gouv at http://dx.doi.org/10.57745/FJFLYB.