Jean-Baptiste Leblond

Alain Karma

Véronique Lazarus

Theoretical analysis of crack front instability in mode I+III

Keywords: Crack front, mode I+III, instability, elliptic helix, Griffith criterion, principle of local symmetry

This paper focusses on the theoretical prediction of the widely observed crack front instability in mode I+III, that causes both the crack surface and crack front to deviate from planar and straight shapes, respectively. This problem is addressed within the classical framework of fracture mechanics where the crack front evolution is governed by conditions of constant energyrelease-rate (Griffith criterion) and vanishing stress intensity factor of mode II (principle of local symmetry) along the front. The formulation of the linear stability problem for the evolution of small perturbations of the crack front exploits previous results of Movchan et al. (1998) (suitably extended) and Gao and Rice (1986), which are used to derive expressions for the variations of the stress intensity factors along the front resulting from both in-plane and out-of-plane perturbations. We find exact eigenmode solutions to this problem, which correspond to perturbations of the crack front that are shaped as elliptic helices with their axis coinciding with the unperturbed straight front and an amplitude exponentially growing or decaying along the propagation direction. Exponential growth corresponding to unstable propagation occurs when the ratio of the unperturbed mode III to mode I stress intensity factors exceeds some "threshold" depending on Poisson's ratio. Moreover, the growth rate of helical perturbations is inversely proportional to their wavelength along the front. This growth rate therefore diverges when this wavelength goes to zero, which emphasizes the need for some "regularization" of crack propagation laws at very short scales. This divergence also reveals an interesting similarity between crack front instability in mode I+III and well-known growth front instabilities of interfaces governed by a Laplacian or diffusion field.

Introduction

The phenomenon of crack front instability in mode I+III, leading to growing deviations of the crack from coplanarity, has been known for a long time and studied experimentally quite extensively [START_REF] Sommer | Formation of fracture "lances" in glass[END_REF][START_REF] Knauss | An observation of crack propagation in antiplane shear[END_REF][START_REF] Palaniswamy | Crack extension in brittle solids[END_REF][START_REF] Hourlier | Fissuration par fatigue sous sollicitations polymodales (mode I ondulé + mode III permanent) d'un acier pour rotors 26NCDV14[END_REF][START_REF] Pollard | Formation and interpretation of dilatant echelon cracks[END_REF][START_REF] Suresh | Combined mode I -mode III fracture of fatigueprecracked alumina[END_REF][START_REF] Pollard | Progress in understanding jointing over the past century[END_REF][START_REF] Yates | Mixed-mode (I+III) fatigue thresholds in a forging steel[END_REF][START_REF] Hull | Tilting cracks: the evolution of fracture surface topology in brittle solids[END_REF][START_REF] Yates | Crack propagation under mixed mode (I+III) loading[END_REF][START_REF] Hubbard | Determining the fatigue threshold and crack propagation properties of a steel under mixed mode I and III loading[END_REF][START_REF] Hull | The effect of mixed mode I/III on crack evolution in brittle solids[END_REF][START_REF] Cooke | Fracture propagation paths under mixed mode loading within rectangular blocks of polymethyl methacrylate[END_REF][START_REF]Quelques problèmes tridimensionnels de mécanique de la rupture fragile[END_REF]Lazarus et al., 2001b[START_REF] Buchholz | Comparison of predictions by mode II or mode III criteria on crack front twisting in three or four point bending experiments[END_REF][START_REF] Lin | Criterion for initiation of cracks under mixed-mode I+III loading[END_REF]. The fact that a qualitatively similar instability is observed in diverse materials such as glass [START_REF] Sommer | Formation of fracture "lances" in glass[END_REF], alumina [START_REF] Suresh | Combined mode I -mode III fracture of fatigueprecracked alumina[END_REF], steels [START_REF] Hourlier | Fissuration par fatigue sous sollicitations polymodales (mode I ondulé + mode III permanent) d'un acier pour rotors 26NCDV14[END_REF][START_REF] Yates | Mixed-mode (I+III) fatigue thresholds in a forging steel[END_REF][START_REF]Quelques problèmes tridimensionnels de mécanique de la rupture fragile[END_REF], rocks [START_REF] Pollard | Formation and interpretation of dilatant echelon cracks[END_REF][START_REF] Pollard | Progress in understanding jointing over the past century[END_REF][START_REF] Cooke | Fracture propagation paths under mixed mode loading within rectangular blocks of polymethyl methacrylate[END_REF] and PMMA [START_REF] Buchholz | Comparison of predictions by mode II or mode III criteria on crack front twisting in three or four point bending experiments[END_REF], which have different atomic scale structures, suggests that the standard tools of macroscopic fracture mechanics should be able to capture some of its basic aspects.

In qualitative terms, the evolution of the crack seems to be governed by a general tendency of its surface to rotate around the direction of propagation, so as to eliminate mode III and reach a situation of pure mode I. When mode III is initially predominant, the crack is observed to discontinuously tilt around the direction of propagation [START_REF] Knauss | An observation of crack propagation in antiplane shear[END_REF][START_REF] Palaniswamy | Crack extension in brittle solids[END_REF]. In this case the original and new crack surfaces do not intersect along the original crack front. Thus the new crack surface cannot develop smoothly from this original front; in practice "lenses" or "plates" inclined with respect to the original front grow from a number of discrete points, implying a very complex geometry. When mode I is initially significant, on the other hand, the crack often develops smoothly from its original front, twisting continuously around its direction of propagation [START_REF] Sommer | Formation of fracture "lances" in glass[END_REF]. The resultant geometric situation is then simpler. Because of the complexity of the crack geometry in the first case, we shall essentially concentrate in this paper on the sole second case of continuous twisting.

The thickness of the body in the direction of the crack front also plays an important role. In massive bodies, the front commonly splits into multiple "facets" (or "lances" in [START_REF] Sommer | Formation of fracture "lances" in glass[END_REF]'s terminology). In the less common case of thin bodies, the crack surface appears to rotate macroscopically as a whole, thus more or less quickly achieving a state of pure mode I along the whole front [START_REF]Quelques problèmes tridimensionnels de mécanique de la rupture fragile[END_REF]. However careful examination of the fracture surface in the vicinity of the initial front reveals that even in this case, crack front segmentation occurs on a smaller scale than the sample thickness.

It has been observed by [START_REF] Hourlier | Fissuration par fatigue sous sollicitations polymodales (mode I ondulé + mode III permanent) d'un acier pour rotors 26NCDV14[END_REF] that front segmentation leads to two types of facets: "type A" facets of large width, rotating in such a way that the local stress intensity factor of mode I tends to increase with the distance of propagation while that of mode III tends to decrease; and "type B" facets rotating oppositely, so that the local stress intensity factor of mode I tends to decrease while that of mode III tends to increase. [START_REF] Hourlier | Fissuration par fatigue sous sollicitations polymodales (mode I ondulé + mode III permanent) d'un acier pour rotors 26NCDV14[END_REF] also noted that the crack propagates more along type A facets than along type B ones. This observation suggests that the crack would ideally like to develop exclusively along type A facets, and that type B facets are present only because they are "geometrically necessary" to join type A ones.

Another common observation is that the facets tend to "coarsen" in time. In Lazarus (1997)'s and [START_REF] Yates | Crack propagation under mixed mode (I+III) loading[END_REF]'s experiments of 3-point bending of beams containing an inclined crack, for instance, the facets merged into the mean surface of the crack, itself only a larger facet, after a certain distance of propagation.

Theoretical studies of the phenomenon are much more scarce than experimental ones. In the case of discontinuous tilting, [START_REF] Pollard | Formation and interpretation of dilatant echelon cracks[END_REF] and [START_REF] Cooke | Fracture propagation paths under mixed mode loading within rectangular blocks of polymethyl methacrylate[END_REF] suggested to estimate the angle of tilt from the hypothesis that the stress intensity factor of mode III after the tilt (evaluated in a very approximate way) should be zero. Lazarus et al. (2001b) showed that this leads to overestimating this angle by a factor of about 2 in practice, and proposed to modify [START_REF] Pollard | Formation and interpretation of dilatant echelon cracks[END_REF]'s and [START_REF] Cooke | Fracture propagation paths under mixed mode loading within rectangular blocks of polymethyl methacrylate[END_REF]'s formula accordingly. But this formula, whether in its original or modified form, is purely heuristic anyway and lacks any theoretical justification. Lazarus et al. (2001b) tried to rationalize experimental observations in terms of the standard energetic theory of fracture. Their main contribution was to show, using theoretical estimates of the stress intensity factors after some continuous twisting (Lazarus et al., 2001a), that the energy-releaserate is larger at the center of type A facets than at that of type B ones1 , implying that propagation along type A facets is "energetically favored". This provided a qualitative explanation to [START_REF] Hourlier | Fissuration par fatigue sous sollicitations polymodales (mode I ondulé + mode III permanent) d'un acier pour rotors 26NCDV14[END_REF]'s observation that the length of type A facets is larger than that of type B ones. Also, Lazarus et al. (2001b) provided semi-heuristic estimates of the "rate of rotation" of the crack around its direction of propagation, agreeing tolerably well with experimental measurements performed on 3-point bent beams containing an inclined crack. These estimates were recently improved by [START_REF] Buchholz | Comparison of predictions by mode II or mode III criteria on crack front twisting in three or four point bending experiments[END_REF] by accounting for the influence of mode II near the free surface of such beams. Even more recently, [START_REF] Lin | Criterion for initiation of cracks under mixed-mode I+III loading[END_REF] introduced the reasonable notion that the energy dissipated by fracture should be larger on type B facets than on type A ones, due to different physical mechanisms, and showed how the difference in width of the two types of facets could be related to this feature.

In spite of their qualities, none of these theoretical studies has succeeded in providing a simple, clear and rigorous explanation of crack front instability in mode I+III, based on the standard tools of fracture mechanics.

However, a new avenue to the problem was very recently opened by [START_REF] Pons | Helical crack-front instability in mixed-mode fracture[END_REF]. These authors performed numerical simulations of crack propagation in mode I+III using a previously developed "phase field" model of brittle fracture [START_REF] Kessler | Phase-field model of mode III dynamic fracture[END_REF], which includes a phenomenological description of short scale failure in the process zone around the crack front. At least for values of the ratio of the initial mode III to mode I stress intensity factors larger than some threshold, simulations reproduced salient features of the crack front instability, including the tendency of the crack to deviate out of its original plane through formation of an array of inclined facets and the tendency of the energetically-favored A facets to propagate ahead of B facets. In addition, simulations led to two key findings. Firstly, the initial development of the crack front instability, in what can be considered the "linear regime" of the instability, is associated with the amplification in time of small helically shaped perturbations, where the helix axis is parallel to the straight unperturbed front and the amplification rate is maximum for a preferred wavelength. Secondly, although facets form initially at this preferred wavelength, in the subsequent "nonlinear regime" of instability, they tend to coarsen in time with the mean facet spacing increasing roughly proportionally to the mean distance of crack advance. This finding strongly suggests that facet coarsening is generally the dominant mechanism that allows a macroscopic rotation of the average crack front position (i.e. the envelope of the tips of A facets) in both thin and thick bodies.

The relationship of the phase-field model to conventional fracture mechanics has been explored by [START_REF] Hakim | Laws of crack motion and phase-field models of fracture[END_REF], who showed that quasi-static crack propagation in this model is approximately governed by a Griffith-type condition of constant energyrelease-rate plus, in isotropic media, [START_REF] Goldstein | Brittle fracture of solids with arbitrary cracks[END_REF]'s condition of zero stress intensity factor of mode II (principle of local symmetry) all along the crack front. This suggests that the problem of crack front instability in mode I+III might be treatable analytically by a linear stability analysis based on a double criterion consisting of a Griffith condition and the principle of local symmetry, combined with a first-order calculation of the variation of the stress intensity factors resulting from the perturbation of the crack.

The purpose of this paper is to describe such an analysis, apparently the first of its kind. It is organized as follows:

• Section 2 explains general hypotheses and notations, and contains in particular a geometric description of the perturbations of the crack considered. • Section 3, based on previous works of [START_REF] Gao | Shear stress intensity factors for planar crack with slightly curved front[END_REF] and [START_REF] Movchan | On perturbations of plane cracks[END_REF] (suitably extended), provides first-order formulae for the variations of the stress intensity factors along the crack front resulting from these perturbations. • Section 4 expounds the linear stability analysis of crack propagation in mode I+III announced, based on these formulae. • Section 5 finally provides a summary of, and some comments on, the results obtained.

Hypotheses and notations

In the reference unperturbed configuration, the geometry consists of an infinite body made of some isotropic elastic material, containing a planar, semi-infinite crack with a straight front (Figure 1). The definition of the axes Ox, Oy, Oz follows the usual convention: Ox is oriented along the direction of propagation of the crack, Oy is perpendicular to the crack plane and Oz is parallel to the crack front. The body is loaded through some system of body forces or surface tractions independent of the coordinate z, so that the problem is in reality 2D in the Oxy plane and the stress intensity factors K 0 I , K 0 II , K 0 III are independent of the position along the crack front. Some general formulae to follow will be provided for arbitrary values of K 0 I , K 0 II , K 0 III , but in the ensuing application K 0 II will be assumed to be zero.

In the true, perturbed configuration, the crack front takes the shape of an elliptic helix at every instant (Figure 2).2 This helix is described by the following parametric equations, where x represents the current position of the unperturbed front, z that of the current point on this front, and x , y , z the coordinates of the current point on the perturbed one:

⎧ ⎪ ⎪ ⎪ ⎪ ⎪ ⎨ ⎪ ⎪ ⎪ ⎪ ⎪ ⎩ x = x + A x e x/a cos(kz) y = A y e x/a sin(kz) z = z.
(1)

In these expressions, is a small parameter, A x and A y are unknown "amplitudes", a is a parameter homogeneous to a distance and k is an unknown "wavenumber" which can safely be assumed to be positive.3 Equations (1) mean that the helix grows (for a > 0) or decays (for a < 0) in proportion in the directions Ox and Oy and exponentially with the distance of propagation; the dimensionless quantity 1/(ka) represents its "normalized" rate of growth or decay.

The perturbed crack surface is generated by this family of helices, with x varying from -∞ to +∞. To find its equation, it suffices to note that since x = x + O( ) and z = z, equation (1) 2 may be rewritten in the form y = A y e x /a sin(kz ) to first order in . Hence the perturbed crack surface is obtained from the unperturbed one through application of a small displacement φ y (x, z) in the direction Oy given by φ y (x, z) = A y e x/a sin(kz).

( In addition, the perturbed crack front is obtained from the unperturbed one through application of a small displacement φ x (x, z) in the direction Ox of the tangent plane given by φ x (x, z) = A x e x/a cos(kz).

(3)

The perturbations of the crack surface and crack front envisaged involve two characteristic distances, |a| in the direction Ox and the wavelength λ ≡ 2π/k in the direction Oz. Both of these distances are assumed to be small in comparison to those defined by:

• the characteristic dimensions of the body and the crack -this is the key hypothesis that allows to schematize the unperturbed crack as semi-infinite and the unperturbed front as straight; • the ratios of the various constants (stress intensity factors, non-singular stresses, ...) characterizing the unperturbed stress field near the crack front -this will allow to express this stress field in terms of the sole unperturbed stress intensity factors in the perturbation analysis to follow.

3 The stress intensity factors along the perturbed crack front

Preliminary comments

It was noted above that the perturbed crack geometry is obtained from the unperturbed one through application of an out-of-plane displacement to the crack surface, φ y (x, z), plus an in-plane displacement to the crack front, φ x (x, z). Formulae for the variations δK α (z) (α = I, II, III) of the stress intensity factors accurate to first order in may be obtained by just adding the contributions arising from these two perturbations, δ 1 K α (z) and δ 2 K α (z).

The situation is somewhat complex with regard to the variations δ 1 K α (z) arising from the out-of-plane perturbation of the crack surface, and requires comments.

Early attempts at the calculation of the δ 1 K α (z) [START_REF] Gao | Three-dimensional slightly nonplanar cracks[END_REF][START_REF] Xu | An analysis of non-planar crack growth under mixed mode loading[END_REF] yielded results in conflict and contradicted by later ones. This was due to errors evidenced and analyzed in later papers (Ball and Larralde, 1995;[START_REF] Movchan | On perturbations of plane cracks[END_REF][START_REF] Leblond | Crack paths in three-dimensional solids. II: Three-term expansion of the stress intensity factors -Applications and perspectives[END_REF]Lazarus et al., 2001a). Ball and Larralde (1995) re-calculated δ 1 K I (z), but only for a mode I loading. [START_REF] Movchan | On perturbations of plane cracks[END_REF] then calculated all variations δ 1 K α (z) for arbitrary mixed-mode loadings. Their result for δ 1 K I (z) fully confirmed that of Ball and Larralde (1995). Finally Leblond and coworkers [START_REF] Leblond | Crack paths in three-dimensional solids. I: Two-term expansion of the stress intensity factors -Application to crack path stability in hydraulic fracturing[END_REF][START_REF] Leblond | Crack paths in three-dimensional solids. II: Three-term expansion of the stress intensity factors -Applications and perspectives[END_REF]Lazarus et al., 2001a) again calculated all the δ 1 K α (z), now for a crack perturbation confined to the immediate vicinity of the front, but without the assumption of almost identical orientations of the old and new tangent planes, that is for arbitrary values of the spatial derivatives of the perturbation. With the additional assumption that these derivatives were small, their results confirmed all terms in Movchan et al. (1998)'s formulae "accessible" by this approach, that is, pertaining to the influence of the near-front perturbation.

The conclusion of this literature survey is that [START_REF] Movchan | On perturbations of plane cracks[END_REF]'s work is the only one which provides the expression of the δ 1 K α (z) with enough generality for our present needs. That it is in addition reliable is strongly suggested by the satisfactory comparisons made with partial results of other authors.

This question being settled, another difficulty arises. Indeed there are three types of terms in [START_REF] Movchan | On perturbations of plane cracks[END_REF]'s formulae for the δ 1 K α (z):

• local terms involving only the values of the perturbation and its two derivatives at that point where the perturbed stress intensity factors are evaluated; • semi-local terms involving values of the perturbation along the whole crack front;

• non-local terms involving values of the perturbation on the whole crack surface. Now, whereas [START_REF] Movchan | On perturbations of plane cracks[END_REF] calculated all local and semi-local terms in detail, they provided explicit expressions of the non-local ones only for perturbations independent of the x coordinate, which is insufficient for our purposes. Additional, non-trivial work is needed to derive formulae for the non-local terms in full generality.

With regard to the variations δ 2 K α (z) arising from the in-plane perturbation of the crack front, the situation is much more favorable since fully explicit formulae for the specific crack geometry considered were derived 25 years ago by [START_REF] Gao | Shear stress intensity factors for planar crack with slightly curved front[END_REF].

The Movchan-Gao-Willis formulae

With the hypotheses made, [START_REF] Movchan | On perturbations of plane cracks[END_REF]'s formulae for the δ 1 K α (z), when the unperturbed front lies at the position x = 04 , reduce to:

⎧ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎨ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎩ δ 1 K I (z) = - 3 2 K 0 II ∂φ y ∂x (0, z) -2K 0 III ∂φ y ∂z (0, z) + δ 1 K skew I (z); δ 1 K II (z) = K 0 I 2 ∂φ y ∂x (0, z) - 2 -3ν 2(2 -ν) K 0 I π P V +∞ -∞ φ y (0, z ) -φ y (0, z) (z -z) 2 dz ; δ 1 K III (z) = 2(1 -ν) 2 2 -ν K 0 I ∂φ y ∂z (0, z) (4)
where ν denotes Poisson's ratio. The only non-local term here is δK skew I (z), the precise expression of which is established in Subsection 3.3 below. The explanation of the notation is that this term is connected to Bueckner's skew-symmetric crack-face weight functions, which provide the stress intensity factors along the crack front arising from application of a pair of identical (not opposite!) point forces on the crack faces; see [START_REF] Movchan | On perturbations of plane cracks[END_REF] for details.

Several terms present in [START_REF] Movchan | On perturbations of plane cracks[END_REF]'s formulae have been omitted here:

• terms arising from the variation of the initial unperturbed stress field in the direction of the crack front, which vanish here since this stress field is assumed to be independent of z; • terms proportional to the initial unperturbed non-singular stresses T 0 xx , T 0 zz , T 0 xz , which are negligible here because the characteristic distances |a| and λ ≡ 2π/k of the perturbation are assumed to be small in comparison to the various characteristic distances (K 0 α /T 0 ij ) 2 defined by the loading; • terms proportional to the coefficients characterizing the third term of the expansion of the unperturbed stresses near the initial crack front, which are negligible for the same reason as those proportional to the unperturbed non-singular stresses.

For the perturbation given by equation ( 2), the integral in equation ( 4) 2 may be calculated by standard methods. One thus gets

⎧ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎨ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎩ δ 1 K I (z) = - 3 2 K 0 II A y a sin(kz) -2K 0 III A y k cos(kz) + δ 1 K skew I (z); δ 1 K II (z) = K 0 I 2 A y a sin(kz) + 2 -3ν 2(2 -ν) K 0 I A y k sin(kz); δ 1 K III (z) = 2(1 -ν) 2 2 -ν K 0 I A y k cos(kz).
(5)

3.3 Calculation of δ 1 K skew I
The detailed expression of δ 1 K skew I was established in Appendix B of [START_REF] Movchan | On perturbations of plane cracks[END_REF]'s paper in the special case of a perturbation independent of x; the first task is to extend it to the general case. Adapting elements of this Appendix, one gets after a few integrations by parts, for an arbitrary perturbation φ y (x, z), when the unperturbed front lies at the position x = 0:

δ 1 K skew I (z) = √ 2 4π 1 -2ν 1 -ν Re 0 -∞ dx +∞ -∞ K 0 III -i(1 -ν)K 0 II (-x + i(z -z )) 3/2 ∂φ y /∂z(x , z ) √ -x dz . (6)
We shall now integrate on z first, unlike [START_REF] Movchan | On perturbations of plane cracks[END_REF] who integrated on x first. 5

For a perturbation of type (2), equation ( 6) yields

δ 1 K skew I (z) = √ 2 8π 1 -2ν 1 -ν A y k × ×Re K 0 III -i(1 -ν)K 0 II 0 -∞ e x /a √ -x [I(k; x , z) + I(-k; x , z)] dx (7)
where

I(μ; x , z) ≡ +∞ -∞ e iμz dz [-x + i(z -z )] 3/2 (x < 0). (8)
The integrals I(k; x , z) and I(-k; x , z) are calculated in Appendix A, and the results are given by equations (A.1) and (A.8) of this Appendix. It follows from these equations that

δ 1 K skew I (z) = 1 -2ν 1 -ν k 3/2 √ 2π A y K 0 III cos(kz) -(1 -ν)K 0 II sin(kz) 0 -∞ e kx e x /a √ -x dx . (9)
Note that the integral here converges only provided that the condition

k + 1 a > 0 ⇔ 1 a > -k (10)
is satisfied. In the unstable case where perturbations grow as the crack advances (a > 0), this condition is always satisfied since k > 0. In the stable case where perturbations decay (a < 0), this condition is satisfied only if 1/(ka) > -1. However, we shall find that the linear stability spectrum, i.e. the relationship between k and a obtained finally, automatically satisfies this condition. Therefore, our analysis also describes stable perturbations without the necessity to introduce, when a < 0, an arbitrary cut-off of the exponentially increasing amplitude in the limit x → -∞.

5 One advantage of the order of integration adopted here is that it in fact allows to calculate δ 1 K skew I for perturbations varying arbitrarily in the direction x, not necessarily as a growing exponential.

The value of the integral in the right-hand side of equation ( 9) is provided by [START_REF] Gradshteyn | Table of Integrals, Series, and Products[END_REF]'s formula (3.361.1):

0 -∞ e kx e x /a √ -x dx ≡ π k F (ka) (11)
where

F (u) ≡ u 1 + u . ( 12
)
The expression (9) of δ 1 K skew I (z) then becomes

δ 1 K skew I (z) = 1 -2ν √ 2(1 -ν) A y k K 0 III cos(kz) -(1 -ν)K 0 II sin(kz) F (ka). ( 13
)

Final formulae for the perturbed stress intensity factors

In the case considered, Gao and Rice (1986)'s formulae for the variations δ 2 K α (z) of the stress intensity factors arising from the in-plane perturbation φ x (x, z) of the crack front, when the unperturbed front lies at the position x = 0, read as follows:

⎧ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎨ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎩ δ 2 K I (z) = K 0 I 2π P V +∞ -∞ φ x (0, z ) -φ x (0, z) (z -z) 2 dz ; δ 2 K II (z) = - 2 2 -ν K 0 III ∂φ x ∂z (0, z) + 2 -3ν 2 -ν K 0 II 2π P V +∞ -∞ φ x (0, z ) -φ x (0, z) (z -z) 2 dz ; δ 2 K III (z) = 2(1 -ν) 2 -ν K 0 II ∂φ x ∂z (0, z) + 2 + ν 2 -ν K 0 III 2π P V +∞ -∞ φ x (0, z ) -φ x (0, z) (z -z) 2 dz ,
(14) that is, for the perturbation defined by equation (3):

⎧ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎨ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎩ δ 2 K I (z) = - K 0 I 2 A x k cos(kz); δ 2 K II (z) = - 2 -3ν 2(2 -ν) K 0 II A x k cos(kz) + 2 2 -ν K 0 III A x k sin(kz); δ 2 K III (z) = - 2(1 -ν) 2 -ν K 0 II A x k sin(kz) - 2 + ν 2(2 -ν) K 0 III A x k cos(kz). ( 15
)
Combining equations ( 5), ( 13) and ( 15), one obtains the total variations δK α (z) ≡ δ 1 K α (z) + δ 2 K α (z) of the stress intensity factors arising from both perturbations, when the unperturbed front lies at the position x = 0:

⎧ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎨ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎩ δK I (z) = - K 0 I 2 A x k cos(kz) -K 0 II A y a sin(kz) 3 2 + 1 -2ν √ 2 ka F (ka) +K 0 III A y k cos(kz) -2 + 1 -2ν √ 2(1 -ν) F (ka) ; δK II (z) = K 0 I 2 A y a sin(kz) 1 + 2 -3ν 2 -ν ka - 2 -3ν 2(2 -ν) K 0 II A x k cos(kz) + 2 2 -ν K 0 III A x k sin(kz); δK III (z) = 2(1 -ν) 2 2 -ν K 0 I A y k cos(kz) - 2(1 -ν) 2 -ν K 0 II A x k sin(kz) - 2 + ν 2(2 -ν) K 0 III A x k cos(kz) (16)
where the function F is given by equation ( 12).

It is easy to establish formulae similar to ( 16) but providing the variations of the stress intensity factors when the unperturbed front lies at an arbitrary position x = 0 instead of x = 0. The conclusion is that equations ( 16) are still applicable provided that all terms are multiplied by the same factor e x/a .

Elements of verification

Since formulae ( 16) play a fundamental role in the sequel, it is important to warrant that they are free of any error. That this is indeed very likely results from the following elements:

• All local and semi-local terms, that is all terms except δ 1 K skew I in equations (4), taken from [START_REF] Movchan | On perturbations of plane cracks[END_REF]'s paper, have been confirmed by Leblond and coworkers [START_REF] Leblond | Crack paths in three-dimensional solids. I: Two-term expansion of the stress intensity factors -Application to crack path stability in hydraulic fracturing[END_REF][START_REF] Leblond | Crack paths in three-dimensional solids. II: Three-term expansion of the stress intensity factors -Applications and perspectives[END_REF]Lazarus et al., 2001a) by a completely different method. The same is true of all terms in equations ( 14), taken from Gao and Rice (1986)'s paper.

• On the other hand [START_REF] Movchan | On perturbations of plane cracks[END_REF]'s explicit formula for δ 1 K skew I was never confirmed in later papers, and furthermore had to be extended here to perturbations of the crack plane depending on x. However, as already mentioned, this term is connected to Bueckner's skew-symmetric crack-face weight functions, the expressions of which are available from other sources, notably [START_REF] Bueckner | Weight functions and fundamental fields for the penny-shaped and the half-plane crack in three-space[END_REF] and [START_REF] Kuo | Stress intensity factors for a semi-infinite plane crack under a pair of point forces on the faces[END_REF]. Using the convergent expressions provided by these authors, we have been able to check the correctness of formula (6) for δ 1 K skew I by a reasoning completely independent of what precedes (not expounded here for space reasons).

• Finally all reasonings and calculations leading from the original formula (6) for δ 1 K skew I to its final expression (13) have been checked very carefully. Also, this final expression coincides (with the required changes of notations) with Movchan et al. (1998)'s result (3.42) in the special case of a perturbation of the crack surface independent of x (a = +∞, F (ka) = 1).

Linear stability analysis

Generalities

We shall now study crack propagation in mode I+III, and thus assume the unperturbed stress intensity factor K 0 II of mode II to be zero henceforward. It will also be harmless to assume a positive K 0 III .

In this section, we investigate the following problem: given some suitable propagation criterion, is there a configuration of the crack and its front, other than the trivial one consisting of a perfectly planar crack with a perfectly straight front, satisfying this criterion? It will be shown that some small perturbations of the type defined in Section 2 are exact solutions of the linearized equations of crack propagation for all successive positions of the front. These solutions depart from, or approach more and more in time, the trivial one. The trivial configuration is necessarily stable versus small arbitrary in-plane or outof-plane perturbations of the crack front or crack surface if all linear modes decay with the distance of propagation, and unstable if one of them at least grows.

The criterion adopted here, inspired from Pons and Karma (2010)'s recent numerical study based on a phase-field model, will consist of a combination of a Griffith-type condition and [START_REF] Goldstein | Brittle fracture of solids with arbitrary cracks[END_REF]'s principle of local symmetry. The first condition implies, with the hypothesis of a uniform fracture toughness, that the energy-release-rate must be uniform along the crack front, whereas the second one implies that the stress intensity factor of mode II must be uniformly zero.

It is worth noting that the idea of combining these two conditions to predict out-of-plane propagation of cracks loaded in mode I+III could already be found in the work of [START_REF] Xu | An analysis of non-planar crack growth under mixed mode loading[END_REF]. However the results of these authors were unfortunately corrupted by errors in the formulae used for the perturbed stress intensity factors (see Subsection 3.1 above). In fact, the idea of using the principle of local symmetry to study crack propagation in mode I+III goes back to the earlier work of [START_REF] Gao | Shear stress intensity factors for planar crack with slightly curved front[END_REF], even though these authors could not pursue it since their analysis was limited to in-plane perturbations of the crack front. The idea can even be traced back to the work of [START_REF] Goldstein | Brittle fracture of solids with arbitrary cracks[END_REF]; but it could even less be put to test then because of total lack of studies of perturbations of 3D cracks at the time.

Clearly, the conditions to be written of cancellation of the variations of the stress intensity factor of mode II and the energy-release-rate along the perturbed front, are insensitive to the presence of the factor e x/a multiplying all terms in the expressions (16) of the variations of the stress intensity factors, when the unperturbed front lies at some arbitrary position x. Hence writing these conditions for the sole position x = 0 of this front will be sufficient to ensure their continuous satisfaction along the whole propagation path.

Analysis

The expression of the energy-release-rate after perturbation of the surface and front of the crack is required in order to apply Griffith's condition. By equations ( 16) (with K 0 II = 0), the variation δG(z) of this energy-release-rate is given, when the unperturbed front lies at the position x = 0, by

δG(z) 2(1 -ν 2 )/E = K 0 I δK I (z) + K 0 III δK III (z) 1 -ν = - (K 0 I ) 2 2 A x k cos(kz) - 2 + ν 2(1 -ν)(2 -ν) (K 0 III ) 2 A x k cos(kz) +K 0 I K 0 III A y k cos(kz) - 2 2 -ν + 1 -2ν √ 2(1 -ν) F (ka) (17)
where E denotes Young's modulus.

Let us now write the two conditions of the criterion adopted. Using first the principle of local symmetry to equate δK II (z) to zero, one gets upon use of equation ( 16) 2 :

A y A x = -4 K 0 III K 0 I ka 2 -ν + (2 -3ν)ka . ( 18
)
This equation defines the elliptic shape of the projection of the helix onto the Oxy plane.

Inserting then this result into equation ( 17), one gets

δG(z) 2(1 -ν 2 )/E = A x k cos(kz) - (K 0 I ) 2 2 + (K 0 III ) 2 - 2 + ν 2(1 -ν)(2 -ν) + 4ka 2 -ν + (2 -3ν)ka 2 2 -ν - 1 -2ν √ 2(1 -ν) F (ka) . ( 19
)
Using now Griffith's condition to equate δG(z) to zero, one gets from there, after some algebraic manipulations:

ka = (2 + ν) K 0 III K 0 I 2 + (1 -ν)(2 -ν) 3(2 -ν) -4 √ 2 (1 -2ν)F (ka) K 0 III K 0 I 2 -(1 -ν)(2 -3ν) . ( 20
)
This equation defines the ratio of the characteristic distances |a| and λ ≡ 2π/k of the crack perturbation in the directions Ox and Oz, as a function of the ratio K 0 III /K 0 I .

Our main interest is of course in positive values of the quantity ka, corresponding to exponentially growing perturbations of the crack. For such values, the predictions of equation ( 20) may be qualitatively understood by a simple reasoning, in spite of the fact that it does not provide a truly explicit expression of ka since this quantity unfortunately enters its right-hand side via the function F (ka) given by equation ( 12). Indeed 0 < F (ka) < 1 for ka > 0, and F (ka) is multiplied by the factor 1-2ν in the term [...] of the denominator of the right-hand side of equation ( 20), so that for values of Poisson's ratio close to 1/2, the effect of this function is modest. For ν = 0.4 for instance, the term [...] varies only from 4.80 to 3.67 when ka varies from 0 to +∞. Thus it may approximately be considered as a constant. The quantity ka then becomes a simple homographic function of the square of the ratio K 0 III /K 0 I , the numerator of which is always positive whereas the denominator may be positive or negative. When K 0 III /K 0 I is small, the denominator is negative and so is consequently the predicted value of ka. When K 0 III /K 0 I becomes larger than some critical value (K 0 III /K 0 I ) cr , however, the predicted value of ka becomes positive, and goes to some positive limit when K 0 III /K 0 I further increases to infinity.

Positive solutions in ka of equation ( 20) may also of course be found numerically. Figure 3 shows the results obtained for the normalized growth rate 1/(ka), for the values ν = 0.2, 0.3 and 0.4 of Poisson's ratio, typical for concrete, steel and plexiglass respectively. (The curves were simply obtained by considering given values of ka and then deducing those of K 0 III /K 0 I from equation ( 20), using the fact that it may be solved explicitly with respect to this ratio). This figure fully confirms the above qualitative considerations, even though the values ν = 0.2 and 0.3 are not so close to 1/2. 1/(ka)

K 0 III /K 0 I ν = 0.2 ν = 0.3 ν = 0.4
Fig. 3. Plot of 1/(ka) versus K 0 III /K 0 I for ν = 0.2, 0.3 and 0.4.

Equations ( 18) and ( 20) also allow to calculate the ratio A y /A x as a function of the sole quantity K 0 III /K 0 I . Figure 4 shows the results obtained for the same values of Poisson's ratio as before.

Quite remarkably, the critical value of the ratio K 0 III /K 0 I may be calculated exactly and analytically. Indeed this ratio takes its critical value when the denominator of the righthand side of equation ( 20) becomes zero, and the predicted value of ka infinite. In this case F (ka) = F (+∞) = 1 so that

K 0 III K 0 I cr = (1 -ν)(2 -3ν) 3(2 -ν) -4 √ 2 (1 -2ν) . ( 21
)
This critical ratio is represented versus Poisson's ratio in Figure 5.

1 10 0.1 1 10 -A y /A x K 0 III /K 0 I ν = 0.2 ν = 0.3 ν = 0.4 Fig. 4. Plot of -A y /A x versus K 0 III /K 0 I for ν = 0.2,
(K 0 III /K 0 I ) cr ν Fig. 5. Plot of (K 0 III /K 0 I ) cr versus ν.
Although only positive solutions in ka of equation ( 20) were considered up to now, one may also discuss negative solutions (satisfying the necessary condition ( 10)). Such solutions correspond to perturbations of the crack decaying exponentially with the distance of propagation. In fact a truly complete stability analysis must even consider possibly complex solutions of equation ( 20), corresponding to perturbations of the surface and front of the crack of the form

⎧ ⎪ ⎨ ⎪ ⎩ φ x (x, z) = Re A x e x/a cos(kz) φ y (x, z) = Re A y e x/a sin(kz), (22)
for complex parameters A x , A y and a. The simple qualitative reasoning presented above is no longer applicable to such general solutions, for which ka may be close to -1 and F (ka) large. But the study of these solutions remains possible and even relatively easy, equation ( 20) being basically equivalent to a polynomial equation of the third degree on the variable ζ ≡ F (ka). Such a study is briefly presented in Appendix B. The conclusion is that for values of Poisson's ratio larger than a very small threshold ν c 0.03, and thus for most physically interesting values of ν like 0.2, 0.3, 0.4, equation ( 20) admits only two real solutions in ka. 6 For values of the ratio K 0 III /K 0 I larger than the critical value (K 0 III /K 0 I ) cr defined by equation ( 21), one of these solutions is positive while the other is negative (and satisfies condition ( 10)). For values of K 0 III /K 0 I smaller than (K 0 III /K 0 I ) cr , both solutions are negative (and satisfy condition ( 10)). Figure 6 shows these solutions, in the form of the normalized growth or decay rate 1/(ka), as a function of K 0 III /K 0 I , for the typical value ν = 0.3. (This figure was drawn by the same method as Figure 3). As explained in Appendix B, for very small values of Poisson's ratio (ν < ν c ), equation (20) on ka admits complex solutions and the instability threshold is no longer strictly defined by equation ( 21). We do not elaborate this case in detail here given that values of ν close to zero are of little practical interest.

Effect of the propagation law

Crack propagation does not always obey a Griffith-type criterion. In subcritical crack growth and fatigue, for instance, it is rather governed by some evolution law of the crack front providing its local velocity in terms of the local energy-release-rate. It is therefore interesting to repeat the stability analysis for such an evolution law. This is done in Appendix C, and the main conclusion is the same as for Griffith's criterion: for usual values of Poisson's ratio, coplanar propagation in mode I+III is unstable for values of the ratio K 0 III /K 0 I larger than some threshold. Furthermore the change of propagation law leaves the value of this threshold strictly unmodified, although it does alter that of the normalized growth rate 1/(ka) of the instability modes.

Synthesis of results and discussion

We have presented, apparently for the first time, a rigorous linear stability analysis of the problem of crack propagation in mode I+III. This analysis is based on formulae for the variations of the stress intensity factors resulting from small but otherwise arbitrary deviations of the crack surface from flatness and the crack front from straightness, which are exact to first order in these perturbations. It is self-consistent in the sense that the criterion used, which combines a Griffith-type condition and [START_REF] Goldstein | Brittle fracture of solids with arbitrary cracks[END_REF]'s principle of local symmetry, is continuously satisfied, for the perturbations considered, all along the front and for all successive positions of this front.

The major output of the analysis is the existence, for values of the ratio K 0 III /K 0 I of the mode III to mode I unperturbed stress intensity factors larger than some threshold depending on Poisson's ratio, of a continuous spectrum of linear instability modes. The deviations of the crack surface from flatness and the crack front from straightness corresponding to these modes oscillate sinusoidally in the direction Oz of the front, and grow exponentially in the direction Ox of propagation. The parameter characterizing the rate of growth of each mode in the direction Ox (1/a in the notation introduced above) is inversely proportional to its wavelength (λ ≡ 2π/k) in the direction Oz; the constant of proportionality depends on the values of the ratio K 0 III /K 0 I and Poisson's ratio.

As a consequence, the rate of growth of a given linear instability mode diverges to infinity when its wavelength goes to zero. This raises a basic problem of physical soundness of the model employed, since this wavelength may take arbitrarily small values. It is clear that a "regularization" of the propagation criterion of some sort is needed at short scales, in order to prevent formation of disturbances of the crack geometry varying over vanishingly small distances in the direction of the front and growing arbitrarily quickly in the direction of propagation. The need for such a regularization has already been noted by other authors in various situations, when interest was focussed on small scales, for which the predictions of standard fracture mechanics break down. But it seems more essential here, since the physically absurd prediction of arbitrarily large growth of some disturbances of the crack geometry precludes any study of the evolution of cracks loaded in mode I+III for values of K 0 III /K 0 I larger than the threshold, whatever the scale over which it is envisaged. The role of the regularization will be to predict the preferred initial scale of instability observed experimentally.

It is worth noting that the need for a short scale regularization is well-known in the context of growth front instabilities ranging from viscous fingering to crystal growth via flame propagation (for a collection of articles on this class of instabilities, see [START_REF] Pelcé | Dynamics of Curved Fronts[END_REF]).

In the example of the two-dimensional Saffman-Taylor viscous fingering problem [START_REF] Saffman | The penetration of a fluid into a porous medium or Hele-Shaw cell containing a more viscous fluid[END_REF], the growth rate of sinusoidal perturbations of the front between two fluids of different viscosities diverges inversely proportionally to the wavelength in the absence of surface tension, like in the present mode I+III crack front instability. The cure to this divergence is well-known: the inclusion of surface tension stabilizes short wavelength perturbations and leads to a preferred wavelength, as observed in the phasefield simulations of mode I+III crack front instability [START_REF] Pons | Helical crack-front instability in mixed-mode fracture[END_REF]. Based on this analogy, Pons and Karma have recently proposed a regularization of the fracture problem based on a modification of the principle of local symmetry where the mode II stress intensity factor along the front is required to be proportional to the local out-ofplane curvature and a phenomenological cut-off length on the process zone scale, instead of simply zero [START_REF] Pons | Helical crack-front instability in mixed-mode fracture[END_REF]. An extension of the present linear stability analysis to include this regularization is presently underway.

In addition, in a future paper, we shall envisage interpretation of some experimental results in terms of the results obtained here (and possible extensions). One major challenge in this study will be to reconcile the threshold value of the ratio K 0 III /K 0 I predicted theoretically with that observed experimentally. For ν = 0.3 for instance, equation ( 21) yields a threshold value of 0.52; in contrast, [START_REF] Sommer | Formation of fracture "lances" in glass[END_REF]'s experiments on glass tubes yielded a value of only 0.06, and some experimenters even claim that there is no threshold at all (Ravi-Chandar, 2010). One possible explanation of the discrepancy, among others, could be a very strong influence of initial imperfections which would lower the theoretical threshold. Additional work is obviously needed in order to settle this difficult question.

J(k; 0 + ) = - 1 + i √ 2 0 -∞ e -ikζ -1 dζ (-ζ) 3/2 + -1 + i √ 2 +∞ 0 e -ikζ -1 dζ ζ 3/2 = - 1 + i √ 2 +∞ 0 e ikζ -1 dζ ζ 3/2 + -1 + i √ 2 +∞ 0 e -ikζ -1 dζ ζ 3/2 = √ 2 +∞ 0 [1 -cos (kζ) + sin (kζ)] dζ ζ 3/2 = 2 √ 2 k +∞ 0 sin (kζ) dζ √ ζ + √ 2 +∞ 0 sin (kζ) dζ ζ 3/2 .
The values of these integrals are finally deduced from [START_REF] Gradshteyn | Table of Integrals, Series, and Products[END_REF]'s formulae (3.761.4), (8.338.2) and (8.338.3), and one gets

J(k; 0 + ) = 4 √ πk (k > 0). (A.6)
Then, by equation (A.4),

J(k; ξ) = 4 √ πk e -kξ (k > 0, ξ > 0) (A.7)
so that, by equation (A.2),

I(-k; x , z) = 4 √ πk e kx e -ikz (k > 0, x < 0). (A.8)
B Appendix : study of complex solutions of equation ( 20)

In this Appendix, we study general complex solutions of equation ( 20). The determination of the square root in the expression (12) of the function F (u) is the standard one in the complex plane: that is, its argument u/(1 + u) cannot be a negative real, and its values lie in the half-plane of complex numbers with a positive real part.

The study of solutions of equation ( 20) is made easier by the change of unknown defined by

ζ ≡ F (ka) ≡ ka 1 + ka , (B.1)
where the variable ζ lies in the half-plane Re(ζ) > 0. Expressing ka as ζ 2 /(1ζ 2 ) in equation ( 20), one obtains after some manipulations the following polynomial equation of the third degree in ζ: 20), which is straightforward using the method employed to draw Figures 3 and6. The results of this study are as follows for standard values of Poisson's ratio like 0.2, 0.3, 0.4:

P (ζ) ≡ 4 √ 2 (1 -2ν) K 0 III K 0 I 2 ζ 3 - ⎡ ⎣ 2(4 -ν) K 0 III K 0 I 2 + 2ν(1 -ν) ⎤ ⎦ ζ 2 +(2 + ν) K 0 III K 0 I 2 + (1 -ν)(2 -ν) = 0. (B.
• For values of the ratio K 0 III /K 0 I larger than the critical value (K 0 III /K 0 I ) cr defined by equation ( 21), equation ( 20) has two real solutions in ka, one of which is positive and the other smaller than -1; equivalently, equation (B.2) admits, in addition to its negative solution ζ 1 , two real and positive solutions ζ 2 and ζ 3 . The set of roots of the polynomial P (ζ) is thus complete, there are no additional complex solutions.

• For values of K 0 III /K 0 I smaller than (K 0 III /K 0 I ) cr , equation ( 20) again has two real solutions in ka, both of which are now smaller than -1; equivalently, equation (B.2) again admits two real and positive solutions ζ 2 and ζ 3 in addition to its negative solution ζ 1 . Again, the set of roots of polynomial P (ζ) is thus complete and no additional complex solutions exist.

For exceptional values of Poisson's ratio very close to 0, the situation is different: a detailed study show that for ν ≤ ν c 0.03 and some values of K 0 III /K 0 I , equation (B.2) admits complex conjugate solutions, implying existence of complex conjugate solutions of equation (20). But this case is of little practical interest.

C Appendix : stability analysis for an arbitrary evolution law of the crack front

In this Appendix, we study the possible instability of a crack loaded in mode I+III, propagation being assumed to be governed by an arbitrary evolution law of the crack front connecting its local velocity to the local energy-release-rate.

C.1 Perturbations considered

The in-plane and out-of-plane perturbations considered are of the form

⎧ ⎪ ⎨ ⎪ ⎩ φ x (x, z) = f (x) cos(kz) φ y (x, z) = g(x) sin(kz) (C.1)
where the functions f (x) and g(x) are arbitrary and k denotes a positive wavenumber. These formulae apply to the initial perturbation from which the crack is assumed to evolve, as well as to that which develops from it.

Clearly, the initial perturbations thus considered are not generic since, for instance, the phase difference in the direction Oz between their in-plane and out-of-plane components amounts to precisely π/2, and the phases of these components are independent of x.

Attention is restricted to such perturbations because considering completely arbitrary ones would make the treatment formally much more complex, without changing its main conclusions.

C.2 Calculation of the stress intensity factors and the energy-release-rate

Gao and Rice (1986)'s and [START_REF] Movchan | On perturbations of plane cracks[END_REF]'s formulae for the variations of the stress intensity factors resulting from the in-plane and out-of-plane perturbations are applied, and the treatment of Subsection 3.3 is adapted to re-calculate the quantity δ 1 K skew I with the new, more general hypotheses made here. The resulting expressions of the variations of the stress intensity factors read, when the unperturbed crack front lies at an arbitrary position x:

⎧ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎨ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎩ δK I (x, z) = - K 0 I 2 f (x)k cos(kz) -2K 0 III g(x)k cos(kz) + 1 -2ν 1 -ν k 3/2 √ 2π K 0 III cos(kz) 0 -∞ g(x + x )e kx √ -x dx ; δK II (x, z) = K 0 I 2 g (x) sin(kz) + 2 -3ν 2(2 -ν) K 0 I g(x)k sin(kz) + 2 2 -ν K 0 III f (x)k sin(kz); δK III (x, z) = 2(1 -ν) 2 2 -ν K 0 I g(x)k cos(kz) - 2 + ν 2(2 -ν) K 0 III f (x)k cos(kz), (C.2)
where K 0 II has been assumed to be zero. The expression of the variation of the energyrelease-rate follows from there:

δG(x, z) 2(1 -ν 2 )/E = - (K 0 I ) 2 2 f (x)k cos(kz) - 2 + ν 2(1 -ν)(2 -ν) (K 0 III ) 2 f (x)k cos(kz) - 2 2 -ν K 0 I K 0 III g(x)k cos(kz) + 1 -2ν 1 -ν k 3/2 √ 2π K 0 I K 0 III cos(kz) 0 -∞ g(x + x )e kx √ -x dx .
(C.3)

C.3 Application of the principle of local symmetry and the evolution law of the crack front

We introduce some generic evolution law of the crack front in the form

V (x, z) ≡ V [G(x, z)] (C.4)
where V (x, z) denotes the local normal velocity of the front in the tangent plane to the crack, and G(x, z) the local energy-release-rate. Linearization of this law yields

V (x, z) = V [G 0 + δG(x, z)] = V (G 0 ) 1 + N 2 δG(x, z) G 0 , (C.5)
where

G 0 ≡ 1 -ν 2 E K 0 I 2 + (K 0 III ) 2 1 -ν (C.6)
is the unperturbed energy-release-rate and N is the dimensionless parameter defined by

N ≡ 2 V (G 0 ) V (G 0 ) G 0 . (C.7)
The advantage of the notation introduced in equation (C.7) is that for a Paris-type law, N simply reduces to Paris's exponent, which is independent of G 0 ; note however that N does depend on G 0 for a more general law.

Application of the evolution law of the crack front in its linearized form (C.5) yields

∂ ∂t [x(t) + φ x (x(t), z)] = dx dt (t) + ∂φ x ∂x (x(t), z) dx dt (t) = V (G 0 ) 1 + N 2 δG (x(t), z) G 0 ,
which implies, upon identification of terms of order 0 and 1 in , that

⎧ ⎪ ⎪ ⎨ ⎪ ⎪ ⎩ dx dt = V (G 0 ) ∂φ x ∂x (x, z) dx dt = V (G 0 ) N 2 δG (x, z) G 0 ⇒ ∂φ x ∂x (x, z) = N 2 δG (x, z) G 0 .
Using now equations (C.1) 1 , (C.3) and (C.6), one gets from there the equation

f (x) = Nk 1 + (K 0 III /K 0 I ) 2 1 -ν ⎡ ⎣ - f (x) 2 - 2 + ν 2(1 -ν)(2 -ν) K 0 III K 0 I 2 f (x) - 2 2 -ν K 0 III K 0 I g(x) + 1 -2ν 1 -ν K 0 III K 0 I k 2π 0 -∞ g(x + x )e kx √ -x dx ⎤ ⎦ . (C.8)
Furthermore, the principle of local symmetry implies that δK II (x, z) = 0, that is, by equation

(C.2) 2 , g (x) = - k 2 -ν (2 -3ν)g(x) + 4 K 0 III K 0 I f (x) . (C.9)
Equations (C.8) and (C.9) allow for the calculation of the evolution of the crack from some arbitrary initial configuration defined by the values of f (x 0 ) and g(x), x ≤ x 0 . Note that (C.8) is not a simple differential equation but an integro-differential equation; this is because at a given stage of propagation, the current stress intensity factors depend on the whole previous geometry of the crack, via the non-local term δ 1 K skew I .

C.4 Asymptotic behavior of the perturbation for large times

Although equations (C.8) and (C.9) may be integrated numerically for all possible initial conditions, it is more interesting to use them to analytically study the evolution of the perturbation of the crack for long times, which governs the stability or instability of its unperturbed configuration. More specifically, we shall examine the possibility of some asymptotic exponential divergence of the in-plane and out-of-plane components of the perturbation, that is, of asymptotic expressions of the functions f (x) and g(x) of the form

⎧ ⎪ ⎨ ⎪ ⎩ f (x) ∼ A x e x/a
g(x) ∼ A y e x/a for x → +∞ (C.10)

where A x and A y are unknown "amplitudes" and a a positive parameter homogeneous to a distance.

Combined with these expressions, equation (C.9) yields in the limit x → +∞: .11) this result is exactly identical to (18).

A y A x = -4 K 0 III K 0 I ka 2 -ν + (2 -3ν)ka ; (C
To similarly write the asymptotic form of equation (C.8), it is necessary to determine the asymptotic behavior of the integral 0 -∞ g(x+x )e kx √ -x dx . Now when x moves away from the origin in the negative direction, both g(x + x ) and e kx / √ -x quickly decrease; hence the behavior of the integral is determined by that of the function g in the vicinity of the point x, which means that this function may safely be replaced by its asymptotic expression (C.10) 2 . One thus gets where F (u) is the function defined by equation ( 12). Inserting this result into equation (C.8), one gets in the limit x → +∞:

- 1 2 + K 0 III K 0 I 2 - 2 + ν 2(1 -ν)(2 -ν) + 4ka 2 -ν + (2 -3ν)ka × × 2 2 -ν - 1 -2ν √ 2(1 -ν) F (ka) = 1 Nka 1 + (K 0 III /K 0 I ) 2 1 -ν .
(C.13)

One may check that when the right-hand side of condition (C.13) is zero, that is when the parameter N is infinite, this condition is equivalent to equation ( 20) on ka, applicable to propagation governed by Griffith's criterion. This is not surprising since an evolution law of the front of the form (C.5) predicts, with such a value of N , a very steep increase of the crack front velocity when G(x, z) reaches G 0 , and thus "simulates" a Griffith-type criterion. When N takes finite values, however, the nonzero right-hand side of equation (C.13) brings a correction depending on the evolution law adopted.

Figures C.1 and C.2 show the values of the normalized growth rate 1/(ka) and the ratio A y /A x characterizing the elliptic shape of the helix versus that of the ratio K 0 III /K 0 I , for ν = 0.4 and various values of N ; the curves corresponding to propagation governed by Griffith's criterion, marked "G c ", are also given for reference. (Again, the figures have been obtained by taking various values of ka and then deducing those of K 0 III /K 0 I and A y /A x from equations (C.13) and (C.11), using the fact that the former equation may be solved explicitly with respect to K 0 III /K 0 I ). One clearly sees in Figures C.1 and C.2 that an asymptotic diverging behavior of type (C.10) of the functions f (x) and g(x) is possible only for values of the ratio K 0 III /K 0 I larger than some threshold. Quite remarkably, the value of this threshold can again be calculated analytically and is exactly identical to that for propagation governed by Griffith's criterion, provided by equation ( 21). Indeed 1/(ka) can be observed in Figure C.1 to go to zero when K 0 III /K 0 I gets close to the threshold, so that in this limit, the right-hand side of equation (C.13) vanishes and this equation becomes equivalent to (20).

  Fig. 1. Unperturbed geometry and loading

  Fig. 2. Perturbed geometry

IFig. 6 .

 6 Fig.6. Plot of 1/(ka) (positive and negative solutions) versus K 0 III /K 0 I for ν = 0.3.

  a F (ka) for x → +∞ (C.12)

  Fig. C.1. Plot of 1/(ka) versus K 0 III /K 0 I for ν = 0.4 and various values of N .

  Fig. C.2. Plot of -A y /A x versus K 0 III /K 0 I for ν = 0.4 and various values of N .

  0.3 and 0.4.

	1.2
	1.1
	1
	0.9
	0.8
	0.7
	0.6
	0.5
	0.4
	0.3
	0.2
	0.1	0.15	0.2	0.25	0.3	0.35	0.4	0.45	0.5

  2) It is obvious that over the real interval (-∞, 0], P (ζ) varies from -∞ to some positive value. It follows that equation (B.2) necessarily has a real solution, say ζ 1 , in the interval (-∞, 0). But this solution is inadmissible since it violates the condition Re(ζ 1 ) > 0. (In other words, the equation ka 1+ka = ζ 1 in ka cannot be solved). Equation (B.2) necessarily admits two other solutions, ζ 2 and ζ 3 , which are the only ones of interest. Since the coefficients of ζ 3 and ζ 0 of the polynomial P (ζ) are both positive, the product ζ 1 ζ 2 ζ 3 is negative. Since ζ 1 < 0, it follows that ζ 2 ζ 3 > 0. Hence the solutions ζ 2 and ζ 3 may be real of identical signs, or complex and conjugate. It is easy to study the real solutions in ζ of equation (B.2), since it is equivalent to studying the real solutions in ka of equation (

This property is not a triviality since the stress intensity factors of mode I and III evolve oppositely upon propagation on both type A and B facets.

Only the initial crack configuration and the helical perturbed front are shown here; the perturbed crack surface is not represented because the figure would become illegible.

3 This only requires to suitably choose the sign of A y .

This momentary restrictive hypothesis is made for simplicity and will be relaxed at the end of Subsection 3.4.

The equivalent polynomial equation of the third degree on ζ of course admits three solutions, but one of these is inadmissible since for this value of ζ, the equation F (ka) = ζ in ka cannot be solved; see Appendix B for details.

Acknowledgements

We thank Vincent Hakim for valuable discussions. A.K. was supported by US DOE grant DE-FG02-07ER46400.

A Appendix : calculation of the integrals I(k; x , z) and I (-k; x , z) To calculate the integral I(k; x , z) defined by equation (8) (with k > 0 and x < 0), consider the integral of e ikz [-x + i(zz )] -3/2 over a closed path composed of a real interval [-A, A] completed by a semi-circle of radius A in the upper half of the complex plane. This integral is zero since the pole of the integrand is at z = z + ix , that is in the lower half-plane. Now in the limit A → +∞, the integral from -A to A goes to I(k; x , z), and the integral over the semi-circle goes to zero. It follows that

To calculate I(-k; x , z), write it in the form

where

the problem is thus reduced to calculating the integral J(k; ξ).

To do so, differentiate it with respect to ξ and integrate the result by parts to get

integration of this differential equation yields

To calculate J(k; 0 + ), note first that since obviously +∞ -∞ (0 +iζ) -3/2 dζ = 0, this integral may be expressed as

where the integral is now absolutely convergent, i.e. exists in the standard Lebesgue sense. Now note that since the cut of the power function is along the line of negative reals:

< 0, then 0 + -iζ = e iπ/2 (-ζ) so that (0 + -iζ) -3/2 = e -3iπ/4 (-ζ) -3/2 = -1+i √ 2 (-ζ) -3/2 .

It follows from there, using equation (A.5), that