Camille Gervais

Denis Roze
email: roze@sb-roscoff.fr

Mutation rate evolution in partially selfing and partially asexual organisms

Keywords: clonality, deleterious mutation, modifier model, multilocus population genetics, self-fertilization

Different factors may influence the evolution of the mutation rate of a species: costs associated with DNA replication fidelity, indirect selection caused by the mutations produced (that should generally favor lower mutation rates, given that most mutations affecting fitness are deleterious) and genetic drift, that may render selection acting on weak mutators inefficient. In this paper, we use a two-locus model to compute the strength of indirect selection acting on a modifier locus that affects the mutation rate towards a deleterious allele at a second, linked locus, in a population undergoing partial selfing or partial clonality. The results show that uniparental reproduction increases the effect of indirect selection for lower mutation rates. Extrapolating to the case of a whole genome with many deleterious alleles and introducing a direct cost a DNA replication fidelity, the results can be used to compute the evolutionarily stable mutation rate U . In the absence of mutational bias towards higher U , the analytical prediction fits well with individual-based, multilocus simulation results. When such a bias is added into the simulations, however, genetic drift may lead to the maintenance of higher mutation rates, and this effect may be amplified in highly selfing or highly clonal populations due to their reduced effective population size.

INTRODUCTION

Rates of spontaneous mutation per nucleotide and per cell division span several orders of magnitudes within eukaryotes (e.g., [START_REF] Sung | Driftbarrier hypothesis and mutation-rate evolution[END_REF][START_REF] Lynch | Genetic drift, selection and the evolution of the mutation rate[END_REF], providing evidence that mutation rates evolve over long timescales. Furthermore, mutation rate variation within a species has been documented in various groups of organisms including bacteria, green algae and fruit flies (e.g., [START_REF] Demerec | Frequency of spontaneous mutations in certain stocks of Drosophila melanogaster[END_REF][START_REF] Woodruff | Variation in spontaneous mutation rates in natural population lines of Drosophila melanogaster[END_REF][START_REF] Miller | Spontaneous mutators in bacteria: insights into pathways of mutagenesis and repair[END_REF][START_REF] Haag-Liautard | Direct estimation of per nucleotide and genomic deleterious mutation rates in Drosophila[END_REF][START_REF] Ness | Extensive de novo mutation rate variation between individuals and across the genome of Chlamydomononas reinhardtii[END_REF], suggesting that mutation rates may possibly change rapidly if natural selection can act upon this genetic variation. Changes in mutation rates have indeed been observed during evolution experiments [START_REF] Sniegowski | The evolution of mutation rates: separating causes from consequences[END_REF]: for example, increased mutation rate in evolving populations of bacteria due to the fixation of mutator genotypes (e.g., [START_REF] Sniegowski | Evolution of high mutation rates in experimental populations of E. coli[END_REF], or decreased mutation rate in populations of Drosophila exposed to Xirradiation during several generations [START_REF] Nöthel | Adaptation of Drosophila melanogaster populations to high mutation pressure: evolutionary adjustment of mutation rates[END_REF].

Starting with [START_REF] Sturtevant | Essays on evolution. I. On the effects of selection on the mutation rate[END_REF], a number of evolutionary forces that may influence the evolution of mutation rates have been identified [START_REF] Drake | Rates of spontaneous mutation[END_REF][START_REF] Sniegowski | The evolution of mutation rates: separating causes from consequences[END_REF][START_REF] Baer | Mutation rate variation in multicellular eukaryotes: causes and consequences[END_REF]Lynch, 2010). Because most mutations affecting fitness are deleterious [START_REF] Eyre-Walker | The ditribution of fitness effects of new mutations[END_REF], alleles coding for higher mutation rates should be associated with less fit genetic backgrounds, thus favoring reduced mutation rates. Using a two-locus modifier model in which one locus affects the mutation rate between alleles at a linked locus directly affecting fitness, [START_REF] Kimura | On the evolutionary adjustment of spontaneous mutation rates[END_REF] showed that the strength of selection to reduce mutation in a panmictic, diploid population is approximately sh δu/ (r + sh), where δu is the change in mutation rate caused by the modifier locus, r the recombination rate between the two loci and sh the heterozygous effect of the deleterious allele, assumed different from zero -see [START_REF] Karlin | Towards a theory of the evolution of modifier genes[END_REF] for the case of a fully recessive deleterious allele. This result was later generalized by [START_REF] Kondrashov | Modifiers of mutation-selection balance: general approach and the evolution of mutation rates[END_REF], [START_REF] Dawson | Evolutionarily stable mutation rates[END_REF]Dawson ( , 1999)), Johnson (1999a) and [START_REF] Lynch | The cellular, developmental and population-genetic determinants of mutation-rate evolution[END_REF] to the case of modifiers changing the deleterious mutation rate over a whole genome.

In asexual populations, this effect may be compensated by the higher rate of production of beneficial alleles by mutator genotypes, which may lead to transient increases in mutation rate when mutators hitchhike with the beneficial alleles they created (e.g., [START_REF] Leigh | Natural selection and mutability[END_REF][START_REF] Eshel | Clone selection and optimal rates of mutation[END_REF][START_REF] Taddei | Role of mutator alleles in adaptive evolution[END_REF][START_REF] Tenaillon | Mutators, population size, adaptive landscape and the adaptation of asexual populations[END_REF][START_REF] André | The evolution of mutation rate in finite asexual populations[END_REF]. In sexual populations, however, recombination destroys the association between mutators and beneficial alleles, and selection for reduced mutation due to the effect of deleterious alleles should generally prevail [START_REF] Leigh | Natural selection and mutability[END_REF]Johnson, 1999b).

The maintenance of nonzero mutation rates is often considered as the result of two opposing forces: selection for reduced mutation rates due to the deleterious effect of most mutations, and the intrinsic cost of DNA replication fidelity (e.g., [START_REF] Kimura | On the evolutionary adjustment of spontaneous mutation rates[END_REF][START_REF] Drake | Rates of spontaneous mutation[END_REF][START_REF] Baer | Mutation rate variation in multicellular eukaryotes: causes and consequences[END_REF][START_REF] Sniegowski | The evolution of mutation rates: separating causes from consequences[END_REF]. More recently, [START_REF] Lynch | The cellular, developmental and population-genetic determinants of mutation-rate evolution[END_REF]Lynch ( , 2011) ) proposed that the equilibrium value of the mutation rate may instead result from a balance between indirect selection and genetic drift: indeed, once the mutation rate has decreased to a very low level, the strength of selection for further increases in replication fidelity may become weaker than genetic drift. The mutation rate would thus reach higher values in populations with lower effective population size N e , due to less efficient selection acting on modifier alleles reducing mutation: this agrees with the observation that the mutation rate is lower in species with larger estimated N e (Lynch, 2010;[START_REF] Sung | Driftbarrier hypothesis and mutation-rate evolution[END_REF][START_REF] Lynch | Genetic drift, selection and the evolution of the mutation rate[END_REF].

Based on [START_REF] Kimura | On the evolutionary adjustment of spontaneous mutation rates[END_REF]'s result mentioned above, reproductive systems that reduce effective recombination rates (such as selfing or clonality) should increase the strength of selection for lower mutation rates (as mutators tend to stay longer associated with the deleterious alleles they produce). In the extreme case of full selfing or full clonality, the strength of selection against a mutator allele becomes equivalent to the increase in mutation load that it causes: neglecting drift, this corresponds to the increase in mutation rate caused by the mutator (e.g., [START_REF] Sturtevant | Essays on evolution. I. On the effects of selection on the mutation rate[END_REF][START_REF] Drake | Rates of spontaneous mutation[END_REF]. Using multilocus simulations incorporating a cost of replication fidelity, [START_REF] Sloan | Evolutionary feedbacks between reproductive mode and mutation rate exacerbate the paradox of sex[END_REF] showed that the equilibrium deleterious mutation rate is indeed lower in asexual than in sexual populations, generating an indirect benefit for asexuality. With selfing, selection for lower mutation rates should be further enhanced by the increased fitness effect of deleterious alleles due to increased homozygosity. However, background selection may strongly reduce the effective size of highly selfing or clonally reproducing populations [START_REF] Nordborg | Structured coalescent processes on different time scales[END_REF][START_REF] Glémin | Adaptation and maladaptation in selfing in outcrossing species: new mutations versus standing variation[END_REF][START_REF] Agrawal | Coalescence with background and balancing selection in systems with bi-and uniparental reproduction: contrasting partial asexuality and selfing[END_REF]Roze, 2016) which, according to Lynch's (2010) hypothesis mentioned above, may possibly increase the equilibrium mutation rate. The overall effect of selfing or clonality on the evolution of mutation rates thus remains unclear, and has been little explored.

In this paper, we extend [START_REF] Kimura | On the evolutionary adjustment of spontaneous mutation rates[END_REF] two-locus model to compute the strength of indirect selection acting on a mutation modifier locus in a partially selfing or partially clonal diploid population. The results confirm that uniparental reproduction increases selection against mutator alleles due to stronger associations with deleterious alleles. Under partial selfing, the strength of indirect selection generated by closely linked loci can be approximated by replacing r and h in [START_REF] Kimura | On the evolutionary adjustment of spontaneous mutation rates[END_REF]'s result by effective recombination and dominance coefficients r (1 -F ) and h (1 -F ) + F (e.g., [START_REF] Glémin | Adaptation and maladaptation in selfing in outcrossing species: new mutations versus standing variation[END_REF]Roze, 2016;[START_REF] Hartfield | Limits to adaptation in partially selfing species[END_REF]. However, this approximation underestimates the effect of more distant loci, which may become important when the selfing rate is high. We then extrapolate from this two-locus model to derive expressions for the genomic deleterious mutation rate at equilibrium between indirect selection generated by deleterious alleles and the cost of replication fidelity, and show that these expressions correctly predict the outcome of individual-based multilocus simulations. Finally, using two different simulation models with different assumptions on the genetic architecture of the mutation rate, we show that consistent with Lynch's (2010) hypothesis, populations with lower effective size may maintain higher mutation rates, provided that mutations increasing replication fidelity (antimutator alleles) occur less frequently than those decreasing it (mutator alleles). In some cases, intermediate rates of outcrossing lead to lower mutation rates than obligate outcrossing or obligate selfing/clonality, due to strong reductions in the effective size of highly selfing or clonal populations caused by background selection.

METHODS

Two-locus model. Our analytical model represents a very large (effectively infinite) population of diploid individuals with discrete generations. As in [START_REF] Kimura | On the evolutionary adjustment of spontaneous mutation rates[END_REF], we consider the evolution of a locus (denoted M ) affecting the mutation rate at a second locus (denoted A), which directly affects fitness. Two alleles (denoted 0 and 1) segregate at each locus; we assume that allele 1 at locus A is deleterious, reducing fitness by a factor 1 -sh in heterozygotes and 1 -s in homozygotes. For simplicity, we assume additivity at the mutation modifier locus (locus M ), the mutation rate at locus A being u 0 , u 0 + δu and u 0 + 2δu in individuals with genotype 00, 01 and 11 at locus M (respectively). We assume that mutations from 0 to 1 and from 1 to 0 occur at the same rate; however this hypothesis should not significantly affect the results, as the effect of back mutations will be negligible as long as the deleterious allele stays at low frequency in the population. We also introduce an intrinsic cost of DNA replication fidelity, so that individuals with lower mutation rates pay a fitness cost. For this, we will assume that the fitness of an individual is multiplied by a function f c that increases with the mutation rate, and will consider different forms of cost function. Individuals contribute to the next generation in proportion to their fitness; under partial selfing a proportion α of juveniles is produced by selfing, while under partial asexuality a proportion γ is produced clonally (the remaining proportion 1 -α or 1 -γ being produced by outcrossing with random union of gametes). Finally, r measures the recombination rate between the two loci. We assume that mutation occurs after selection, before recombination; however, assuming that mutation occurs just after recombination yields the same results (as long as the mutation rate depends on the genotype of the diploid parent).

Following previous works (e.g., [START_REF] Barton | Natural and sexual selection on many loci[END_REF][START_REF] Kirkpatrick | General models of multilocus evolution[END_REF][START_REF] Roze | Effects of interference between selected loci on the mutation load, inbreeding depression and heterosis[END_REF]Roze, , 2016)), genetic associations within and between loci are defined as follows. We define X i,1 and X i,2 as indicator variables that equal 1 if a given individual carries allele 1 at locus i on its first or second haplotype (respectively), and 0 otherwise. The frequency of allele 1 at locus i in the whole population is thus given by p

i = E [(X i,1 + X i,2 ) /2]
, where E stands for the average over all individuals. Defining the centered variables ζ i,1 and ζ i,2 as

ζ i,1 = X i,1 -p i , ζ i,2 = X i,2 -p i , (1)
the genetic association between the sets S and T of loci present on the two haplotypes of the same individual is given by:

D S,T = E [ζ S,T ] (2)
where

ζ S,T = ζ S,1 ζ T,2 + ζ S,2 ζ T,1 2 , ζ S,1 = i∈S ζ i,1 , ζ T,2 = i∈T ζ i,2 (3)
(note that D S,T = D T,S ), and where sets S and T may be the empty set ∅, M , A or MA. Associations between genes present on the same haplotype of an individual (D S,∅ ) will be simply denoted D S . For example,

D M,M = E [(X M,1 -p M ) (X M,2 -p M )]
is a measure of the departure from Hardy-Weinberg equilibrium at locus M , while D MA represents the linkage disequilibrium between loci M and A (genetic association between alleles present on the same haplotype). Similarly, D M,A measures the association between alleles at loci M and A present on different haplotypes of the same individual.

In the following, we assume that both loci have weak effects (s, δu small) and derive an expression for the change in p M (the frequency of allele 1 at locus M ) to the first order in s and δu. We will see that this expression includes different forms of genetic associations. Assuming that the effective recombination rate is large relative to δu, we will then use a quasi-linkage equilibrium approximation (QLE) to express these associations in terms of allele frequencies and of the different parameters of the model. Finally, the result will be extrapolated to compute the overall strength of selection on a modifier allele affecting the mutation rate at a large number of selected loci, assuming that genetic associations between those loci can be neglected.

Multilocus simulations. Our simulation program (written in C++ and available from Dryad) is modified from [START_REF] Roze | Effects of interference between selected loci on the mutation load, inbreeding depression and heterosis[END_REF]Roze ( , 2016)), and represents a finite population of N diploid individuals whose genome consists in a linear chromosome along which deleterious mutations occur every generation. For simplicity, all mutations have the same selection and dominance coefficients (s, h). A mutation modifier locus is located at the mid-point of the chromosome, and controls the deleterious mutation rate (the mutation rate of an individual being the average of the values coded by its two modifier alleles). At the start of each generation, the fitness of every individual is computed as

W = f c (U ) (1 -sh) i (1 -s) j (4)
where U is the deleterious mutation rate of the individual (per haploid genome), f c the function representing the cost of replication fidelity, i, j are the number of heterozygous and homozygous deleterious alleles present in the genome of the individual.

In general, we will use the cost function f c (U ) = e -c U , but different functions will also be considered (as explained in the Results section). To form each of the N juveniles of the next generation, an individual is sampled randomly to serve as a maternal parent.

If the fitness of the individual (divided by the maximal fitness in the population) is higher than a random number sampled from a uniform distribution between 0 and 1, the individual is retained, otherwise another individual is sampled until the test is satisfied. Under partial selfing, the mother self-fertilizes with probability α, in which case the new individual is formed by two recombinant chromosomes from the same parent.

Under partial asexuality, the mother reproduces clonally with probability γ, in which case the genome of the new individual is a copy of the maternal genome. If the mother reproduces by outcrossing (with probability 1 -α or 1 -γ), a second individual is sampled using the same procedure as above to serve as a father, and the genome of the new individual is generated from recombinant chromosomes from both parents.

During meiosis, the number of cross-overs is sampled from a Poisson distribution with parameter R (genome map length, in Morgans), and the position of each cross-over is sampled from a uniform distribution. The parameter R will typically take large values (R = 20) in order to mimic a whole genome with multiple chromosomes. Deleterious mutations occur once the parents have been selected, before recombination (note that different offspring from the same parent will carry different new mutations). The number of new deleterious mutations on each chromosome is sampled from a Poisson distribution whose parameter corresponds to the mutation rate of the parent, and the position of each new mutation is sampled from a uniform distribution along the chromosome. Back mutations do not occur, and any deleterious allele that has reached fixation is removed from the population in order to increase execution speed.

During a number of preliminary generations (usually 2000), the deleterious mutation rate of each individual (per haploid genome) is set to U init = 0.1. Then, new alleles coding for different mutation rates can appear at the modifier locus (at rate µ M = 10 -3 U per generation, where U is the deleterious mutation rate of the individual). When a mutation occurs at the modifier locus, the mutation rate coded by the new allele is sampled from a Gaussian distribution centered on the value of the allele before mutation, with variance σ 2 M = 10 -4 (if the new value is negative, it is set to zero). As explained in the Results section, different mutational models were also considered, including a bias towards higher values of U and scaling of σ 2 M with U .

The program generally runs for 10 6 generations, the equilibrium mutation rate being computed by averaging over the last 8 × 10 5 generations.

RESULTS

Change in frequency at the mutation modifier locus. In the Appendix, we

show that an expression for the change in frequency of allele 1 at locus M , to the first order in δu and s is given by:

∆p M ≈ δu d ln f c du (1 + F ) p M q M -sh (D MA + D M,A ) -s (1 -2h) D MA,A .
(5)

The term on the first line of equation 5 represents the effect of the cost of replication fidelity, favoring alleles that increase the mutation rate. This direct selective pressure increases with the rate of inbreeding (1 + F factor), due to increased homozygosity at locus M . The term on the second line represents the effect of indirect selection disfavoring mutator alleles, as these alleles tend to be more often associated with the deleterious allele at the second locus. Indeed, D MA and D M,A represent the association between allele 1 at locus M and the deleterious allele at locus A on the same or on the other haplotype, while D MA,A represents the association between allele 1 at locus M and homozygotes for the deleterious allele at locus A. We now derive expressions for these associations at QLE.

Expressions for genetic associations. In the following, D sel S,T , D mut S,T and D S,T denote genetic associations measured after selection, mutation and recombination/segregation (respectively). Recursions for genetic associations over each step of the life cycle are computed to the first order in s and δu. We will also assume that the deleterious allele stays at low frequency and neglect terms in p A . Under high effective recombination, it is sufficient to express associations to the first order in δu, neglecting terms in s.

However, these expressions diverge when the effective recombination rate tends to zero (due to terms r (1 -α) or r (1 -γ) in their denominators). In order to obtain more accurate expressions for the case where the effective recombination rate is of order s (assuming that the deleterious allele stays at mutation-selection balance and that δu is sufficiently small, so that the QLE approximation holds), we include terms in s in the recursions for genetic associations, by computing the effect of selection on these associations.

D sel MA ≈ E W W ζ MA,∅ + ζ ∅,MA 2 . (6)
Note that equation 6 is an approximation, as it neglects the change due to selection of allele frequencies p M and p A that appear in ζ MA,∅ and ζ ∅,MA (see equation 1). However, taking this change in allele frequencies into account would introduce a term ∆ sel p M × ∆ sel p A of order δus (e.g., equation 23 in [START_REF] Barton | Natural and sexual selection on many loci[END_REF] that can be neglected here. From equations 6 and A6, neglecting terms in (δu) 2 and terms proportional to p A , one obtains:

D sel MA ≈ (1 -sh) D MA -s (1 -h) D MA,A . (7)
Similarly, one arrives at:

D sel M,A ≈ (1 -sh) D M,A -s (1 -h) D MA,A (8)
D sel MA,A ≈ (1 -s) D MA,A . (9)
Equations 7 -9 show that selection tends to decrease genetic associations between the two loci, as it reduces the frequency of the deleterious allele in the population.

Mutation. The effect of genotype-dependent mutation on genetic associations can be computed as follows. Mutation changes the frequency of allele 1 at locus A to

p mut A = u (1 -p A ) + (1 -u) p A (
where u is the average mutation rate at locus A), while in a given individual, X A,i changes to 1 -X A,i with probability u (the mutation rate of the individual at locus A), and remains unchanged with probability 1 -u. Therefore, D MA after mutation is given by:

D mut MA = 1 2 E (X M,1 -p M ) u (1 -X A,1 ) + (1 -u) X A,1 -p mut A + (X M,2 -p M ) u (1 -X A,2 ) + (1 -u) X A,2 -p mut A ( 10
)
where E is the average over all individuals before mutation. Replacing X i,j by ζ i,j + p i and u by u + δu (ζ M,1 + ζ M,2 ), equation 10 yields (after neglecting terms of order (δu) 2 and terms in p A ):

D mut MA ≈ (1 -2u) D sel MA + δu (p M q M + D M,M ) ≈ (1 -2u) D sel MA + δu (1 + F ) p M q M . (11)
Similarly, one obtains:

D mut M,A ≈ (1 -2u) D sel M,A + δu (1 + F ) p M q M , (12)
D mut MA,A ≈ (1 -4u) D sel MA,A (13)
(indeed, one can show that the term in δu of D mut MA,A is of order p A ). Equations 11 and 12 show that the modifier effect generates an association between the allele increasing mutation and the deleterious allele at the other locus, either on the same or on the other haplotype. The buildup of D MA and D M,A is stronger under inbreeding (term 1+F in equations 11 and 12), due to increased homozygosity at the modifier locus: the probability that a deleterious allele is present in the same genome as an allele coding for higher mutation is increased when the high-mutation allele is also present on the other haplotype.

Recombination/segregation. The effects of recombination and segregation depend on the reproductive system. Under partial selfing, we have:

D MA = (1 -r) D mut MA + r D mut M,A (14)
D M,A = α 2 D mut MA + D mut M,A (15)
D MA,A ≈ α 2 D mut MA,A + (1 -r) D mut MA + r D mut M,A , (16)
equation 16 assuming that p A is small. Under partial asexuality:

D MA = [1 -r (1 -γ)] D mut MA + r (1 -γ) D mut M,A (17)
D M,A = γ D mut M,A , D MA,A = γ D mut MA,A . (18)
QLE expressions. The expressions given above can be used to compute solutions for D MA , D M,A and D MA,A at QLE, corresponding to the equilibrium values of genetic associations for given values of allele frequencies (under the assumption that associations equilibrate fast relative to changes in allele frequencies). For the case of a partially selfing population, one obtains (assuming u s, and thus neglecting terms in u):

D MA ≈ δu (1 + F ) (1 + 2F r) p M q M sh e + r [1 -F -s [h e (1 -2F ) -F (2 -F )]] (19) D M,A ≈ δu (1 + F ) F (1 + 2r) p M q M sh e + r [1 -F -s [h e (1 -2F ) -F (2 -F )]] (20)
D MA,A ≈ F D MA (21)
with h e = h (1 -F ) + F , and F = α/ (2 -α). Equations 19 -21 can be used to obtain the strength of indirect selection acting on the modifier locus, given by

s ind = -s h (D MA + D M,A ) + (1 -2h) D MA,A / (p M q M ) (from equation 5): s ind ≈ - δu (1 + F ) s [h e -2F (F -2h e ) r] sh e + r [1 -F -s [h e (1 -2F ) -F (2 -F )]] (22)
In the absence of selfing (α = 0), s ind simplifies to:

s ind ≈ - δu sh 1 -(1 -sh) (1 -r) (23)
which agrees with the result obtained by Kimura when the modifier effect is weak (equation 2 in [START_REF] Kimura | On the evolutionary adjustment of spontaneous mutation rates[END_REF]. Under complete selfing (α = 1), s ind simplifies to -2δu.

When linkage is tight (small r), a separation of timescales argument yields [START_REF] Nordborg | Structured coalescent processes on different time scales[END_REF]Roze, 2016): this may also be seen from equations 19 and 20, neglecting terms in r in the numerators. From this, one obtains that s ind is equivalent to the result obtained under random mating (equation 23), replacing δu by δu (1 + F ), h by h e and r by r e = r (1 -F ). This can also be obtained directly by neglecting the term in r in the numerator of equation 22, and the term in r s in the denominator. Figure 1A compares the prediction from equation 22 and the prediction obtained by replacing δu, h and r in equation 23 by effective parameters: both yield undistinguishable results when linkage is sufficiently tight, but discrepancies appear when loci are loosely linked.

D M,A ≈ F D MA
Under partial asexuality, genetic associations at QLE are given by:

D MA ≈ δu [1 + γ r] p M q M sh + r [1 -γ + sh (1 + γ)] (24)
D M,A ≈ δu γ p M q M 1 -γ + shγ , D MA,A = 0 . ( 25
)
As shown by Figure 1B, the strength of indirect selection s ind is roughly equivalent to the expression obtained under random mating, replacing r by r (1 -γ), as long as the rate of clonal reproduction γ stays moderate. This is not the case for higher values of γ, however, due to the extra contribution of the association D M,A (mutations generated on the other haplotype, that remain associated to the modifier allele due to clonal reproduction). As under complete selfing, one obtains that s ind = -2δu under full asexuality (γ = 1). In both cases, the strength of indirect selection in the absence of recombination corresponds to the increase in mutation load caused by the modifier (as the modifier allele stays associated with the deleterious alleles it generates).

Under selfing, the increase in load equals the increase in haploid mutation rate (since L = u under full selfing), which is 2δu as individuals quickly become homozygous at the modifier locus. Under clonality, the increase in load is twice the increase in haploid mutation rate (since L = 2u, assuming h is significantly greater than zero), which is δu as individuals stay heterozygous at the modifier locus.

Multilocus extrapolation and simulation results. These two-locus results can be extrapolated to the case of a modifier affecting the mutation rate in the whole genome by integrating s ind over the genetic map. If the map length is sufficiently large, the overall effect can be approximated by assuming free recombination among all loci, replacing r by 1/2 and δu by δU (the increase in haploid mutation rate caused by allele 1 at the modifier locus) in the expressions above. The evolutionarily stable mutation rate (at which indirect selection to reduce the mutation rate exactly balances the cost of replication fidelity) can then be obtained by solving s direct + s ind = 0 for U , where the strength of direct selection is given by s direct ≈ δU (1 + F ) d ln f c /dU , assuming δU is small (see equation 5). Under partial selfing and assuming free recombination, this yields

U ESS ≈ c [2 -α (3 -α) (1 -s)] 2s [α + h (2 -α (1 + α))] (26)
when the cost function is given by f c (U ) = e -c U , so that d ln f c /dU = c/U 2 . Equation 26simplifies to c/ (2sh) when α = 0, and to √ c when α = 1. The equivalent expression for partial asexuality is given by:

U ESS ≈ c [1 -γ + sh γ] [1 -γ + sh (1 + 2γ)] sh (1 + γ) [2 (1 -γ) + 3sh γ] , (27)
simplifying to c/2 under full asexuality (γ = 1).

Figure 2 shows that the predicted value for the evolutionarily stable mutation rate obtained by integrating equations 22 and 24-25 over the genetic map (see Mathematica notebook in Supplementary Material for the integration) generally fits well with the multilocus simulation results (using our first simulation program, with a single modifier locus). With a genome map length of 20 Morgans, the simpler expressions obtained assuming free recombination (equations 26 and 27) also provide accurate predictions: discrepancies appear for lower values of s, due to the fact that deleterious alleles segregate at higher numbers of loci, increasing the number of deleterious alleles that are closely linked to the mutation modifier locus. 3B (see Figure 4). Because the mutation rate evolves to very low levels for these parameter values at sufficiently high selfing or clonality rates, we maintained a minimum mutation rate of µ M = 10 -4 at the modifier locus in the simulation program in order to prevent that the population remains stuck in the absorbing state of perfect replication fidelity (U = 0). Supplementary Figure S2 shows the same simulation results as Figures 3E and3F, displayed on a log scale.

The reason for the limited effect of drift on the average value of U observed in Figure 3 is the absence of mutational bias at the modifier locus: indeed, mutations increase or decrease U with the same probability. Drift may have stronger effects when mutations affecting U tend to occur more often in a particular direction [START_REF] Zhang | The anomalous effects of biased mutation revisited: mean-optimum deviation and apparent directional selection under stabilizing selection[END_REF]Lynch, 2011;[START_REF] Charlesworth | Stabilizing selection, purifying selection, and mutational bias in finite populations[END_REF]: most likely in the direction of increased values of U , as it should be easier to impair DNA replication fidelity than to improve it. Indeed, when such a mutational bias is added in the simulation program (by introducing a parameter β such that a fraction β > 0.5 of mutations at the modifier locus tend to increase U ), the mutation rate evolves towards higher values when population size is sufficiently small. When this is the case, U keeps increasing unless one assumes that the average size of mutational steps δU is proportional to U , so that fewer mutations may fix by drift as U reaches higher values (e.g., Lynch, 2011). Figure 5 shows the results of simulations in which δU is sampled in a half-Gaussian distribution with standard deviation σ M U (for both δU > 0 and δU < 0), where U is the mutation rate coded by the modifier allele before mutation. For the parameter values used in Figure 5, the equilibrium mutation rate is slightly higher for N = 10 4 than for N = 10 5 , but generally remains close to the ESS value (see Supplementary Figure S3 for the same results shown on a log scale). Genetic drift has a much stronger effect on the equilibrium value of U when N = 10 3 , however. As expected, increasing the degree of mutational bias (by increasing β) or decreasing the average size of mutational steps δU (by reducing σ 2 M ) amplifies the effects of drift, causing higher values of U to evolve.

Intermediate selfing or clonality rates bring U closer to its ESS value (by increasing the strength of indirect selection), but U may increase again as α or γ approach 1, due to background selection amplifying the effect of drift. In some simulations, background selection caused a runaway process in which the reduction in N e leads to elevated mutation rate, further reducing N e . When this is the case, U reaches very high values and the program has to be stopped manually: this happened for clonality rates higher than the right-most points in Figure 6B, D, F, and for α = 0.2, α = 0.3 and 0.93 ≤ α ≤ 0.99 in Figure 6E (while an equilibrium was reached for α = 1). As shown by Supplementary Figures S4 -S7, qualitatively similar results were obtained under the different cost functions shown in Figure 4.

Similar results were also obtained using our second simulation program, representing a more realistic genetic architecture in which U is controlled by L biallelic loci.

Since the minimal mutation rate (corresponding to the mutation rate of a genotype carrying allele 0 at each mutation modifier locus) is U min , and since the heterozygous effect of each modifier locus is sampled from an exponential distribution with parameter λ (whose average is 1/λ), the average value of U should thus be close to

U drift = U min + L/λ in
regimes where the evolution of U is mainly controlled by drift (assuming large L, and additivity within and between modifier loci). Below U drift , mutations at modifier loci thus tend to increase U (mutational bias). Figure 6 shows simulation results for L = 1000, λ = 500 and U min = 0.01 (so that U drift ≈ 2). As can be seen on the figure, U becomes closer to its mutation-drift equilibrium value U drift as N decreases. Again, increasing the selfing rate or the clonality rate tends to reduce U by increasing the strength of indirect selection; however, above a given threshold for α or γ (that depends on population size), U increases as the selfing or clonality rate increases (due to background selection effects). Note that we could not obtain estimates for the equilibrium mutation rate under high clonality rates (Figure 6B, on the right of the right-most points) because deleterious alleles tend to accumulate in the heterozygous state in the population, causing the program to become increasingly slow.

DISCUSSION

It is widely accepted that mutation rates are maintained at low levels in order to avoid producing an overly strong burden of deleterious alleles. In this paper, we confirm that this deterministic force favoring lower mutation rates is increased by uniparental reproduction, and compute the strength of this effect in populations with intermediate selfing rates or clonality rates. In agreement with previous separation of timescales arguments [START_REF] Nordborg | Structured coalescent processes on different time scales[END_REF]Roze, 2016), when linkage between loci is sufficiently tight the result obtained under partial selfing becomes equivalent to the expression obtained under random mating [START_REF] Kimura | On the evolutionary adjustment of spontaneous mutation rates[END_REF], replacing the dominance coefficient of deleterious alleles and recombination rates by the effective parameters h (1 -F ) + F and r (1 -F ); however, this expression underestimates the strength of indirect selection generated by loosely linked loci. Introducing a direct fitness cost associated with DNA replication fidelity, we could obtain simple approximations for the evolutionarily stable mutation rate, which were confirmed by multilocus, individualbased simulations, in the absence of mutational bias at the mutation modifier locus.

When a mutational bias towards lower fidelity of DNA replication (i.e., higher mutation rate) is added into the model, the average value of the mutation rate at equilibrium becomes more sensitive to genetic drift, in agreement with general results on the evolution of quantitative traits under mutation, selection and drift [START_REF] Zhang | The anomalous effects of biased mutation revisited: mean-optimum deviation and apparent directional selection under stabilizing selection[END_REF][START_REF] Charlesworth | Stabilizing selection, purifying selection, and mutational bias in finite populations[END_REF]. In that case, the mutation rate U may reach high values when mutations affecting U have a weak effect, so that the effect of indirect selection acting on these mutations becomes weaker than the strength of genetic drift:

this is the essence of the argument proposed by Lynch to explain the observed negative correlation between estimated effective population size N e and the mutation rate (e.g., Lynch, 2010;[START_REF] Sung | Driftbarrier hypothesis and mutation-rate evolution[END_REF][START_REF] Lynch | Genetic drift, selection and the evolution of the mutation rate[END_REF]. However, the mutational bias by itself is not a sufficient condition for U to stabilize around a value that depends on N e . This will occur under the extra assumption that the average size of mutational steps at mutation modifier loci increases with the mutation rate, so that the relative effect of drift at these loci becomes weaker as U increases (Lynch, 2011, Figure 5 of the present article). Alternatively, when U is affected by a sufficiently large number of loci with a distribution of effects, the fraction of loci at which indirect selection is weaker than drift will decrease as N e increases, which may also generate a negative relation between the average value of U and N e (Figure 6). After adding these ingredients into our simulation programs, we observed three possible types of outcomes: either drift has only a limited effect, and U stays close to its deterministic equilibrium (at which indirect selection generated by deleterious alleles exactly balances the cost of replication fidelity), or U stabilizes around a higher value that depends on N e , or a runaway process occurs, under which drift causes the evolution of higher U , in turn

reducing N e through background selection effects, causing further increase in U (which should eventually lead to population extinction through mutational meltdown).

Given that per-nucleotide mutation rates are very low in most species, the hypothesis of a mutational bias towards higher mutation rates seems reasonable. Very little is known on the distribution of the effects of mutations affecting U , however, or on how this distribution may change as U evolves. Interactions between mutations affecting DNA repair pathways have been demonstrated in bacteria and yeast. In some cases, positive epistasis (on the mutation rate) has been shown between mutations in genes with partially redundant effects, such as MutM and MutY in E. coli [START_REF] Michaels | Evidence that MutY and MutM combine to prevent mutations by an oxidatively damaged form of guanine in DNA[END_REF][START_REF] Fowler | Interactions among the Escherichia coli mutT, mutM, and mutY damage prevention pathways[END_REF], or MSH3 and MSH6 in S. cerevisiae [START_REF] Marsischky | Redundancy of Saccharomyces cerevisiae MSH3 and MSH6 in MSH2 -dependent mismatch repair[END_REF]. These examples provide possible scenarios under which the effect of a modifier allele would increase with the baseline mutation rate. However, examples of negative epistasis also exist, for example when a repair pathway involves the combined activity of two proteins (such as MSH2 and MSH6 in yeast, [START_REF] Marsischky | Redundancy of Saccharomyces cerevisiae MSH3 and MSH6 in MSH2 -dependent mismatch repair[END_REF].

Overall, we still lack a clear picture of how the average effect of mutator/antimutator alleles should change with the baseline mutation rate. Furthermore, these studies generally focus on mutators with large effects, which may not be representative of the majority of mutations affecting U . Obtaining more detailed information on the genetic architecture of mutation rate variation within natural populations would represent an important progress, but remains a formidable task.

While our results show that the effects of partial selfing and partial clonality are very similar in regimes where drift has only a limited effect, differences appear when drift is stronger and may lead to mutation accumulation. Muller's ratchet occurred in some of our simulations with partial or complete selfing: for s = 0.01 and α = 1 in Figure 2, and in different cases with N = 1000 in Figure 5 (for α = 1 in 5A, α = 0.98 and 1 in 5B, α = 0.2, 0.3 and 0.93 ≤ α ≤ 1 in 5E). It also occurred for N = 1000, 2000, 3000 and high α (≥ 0.9) in Figure 6. Because our simulation program removes fixed mutations, it could still continue to run and U generally stabilized, except for α = 0.2, 0.3 and 0.93 ≤ α ≤ 0.99 in Figure 5E, where U increased to very high values and the program had to be stopped when U > 30. However, in all these cases any real population would eventually reach extinction due to mutation accumulation. By contrast, at high clonality rates mutations tend to accumulate in the heterozygous state (even when γ < 1), an effect already observed in previous studies by Charlesworth et al. (1993a,b) and [START_REF] Michod | Deleterious mutations and selection for sex in finite, diploid populations[END_REF]. Because these mutations could in principle still be removed from the population by rare segregation events, the simulation program does not eliminate them and becomes increasingly slow. This occurred in nearly all cases with γ = 1 in Figures 2 and3, although the mutation rate reached an equilibrium before the program had to be stopped. It also occurred for values of γ higher than the right-most points in Figures 5 and6 (at a faster rate as γ increased), in which case the program had to be stopped before U had reached equilibrium. Again, in all these situations the population would eventually go extinct by mutational meltdown. These results outline two important differences between partial selfing and partial clonality:

(1) the mutation accumulation regime is reached sooner under partial clonality than under partial selfing as the rate of uniparental reproduction increases (due to the absence of segregation), and (2) a runaway process leading to very high mutation rates may occur at moderate selfing rates (α = 0.2, 0.3 in Figure 5E, see also Figures S2 andS3), while it does not occur at higher selfing rates. This last effect is not observed under partial clonality, and could possibly be due to identity disequilibria between selected loci reducing the efficiency of selection (e.g., [START_REF] Roze | Effects of interference between selected loci on the mutation load, inbreeding depression and heterosis[END_REF].

Provided that mutation rate polymorphism exists within populations, the evo-lution of U could in principle interact with the evolution of different aspects of reproductive systems, such as reproductive modes or mating systems. In regimes where the effect of drift on the evolution of U stays negligible, this should favor the evolution of uniparental reproduction, as it should be associated with lower mutation rates [START_REF] Sloan | Evolutionary feedbacks between reproductive mode and mutation rate exacerbate the paradox of sex[END_REF]. However, the evolution of U may have opposite effects in regimes where it is more strongly affected by drift, and where high rates of uniparental reproduction may trigger the evolution of higher mutation rates, through stronger background selection effects. Whether the evolution of the mutation rate would occur on a sufficiently fast timescale to significantly affect the evolution of reproductive systems should in principle depend on the genetic architecture of U and of the reproductive system; however, this should be explored more rigorously using theoretical approaches. Even if mutation rate evolution is not fast enough to have a significant impact on evolutionary transitions between mating systems, U may change in response to a switch in reproductive system. This may affect the species-level se- This may introduce new effects of reproductive systems on mutation rate evolution, as the reproductive system may affect the proportion of mutations on which natural selection is effective (by affecting N e ). Our model also assumes that all mutations affecting fitness are deleterious: while previous theoretical work has shown that beneficial mutations should only have a minor role on the evolution of mutation rates in sexual, outcrossing populations [START_REF] Leigh | Natural selection and mutability[END_REF]Johnson, 1999b), their effect should become more important in populations undergoing high rates of selfing or clonality (since mutator alleles can stay associated with the beneficial alleles they produced), and may increase the equilibrium mutation rate in such populations. Last, as in most theoretical studies of mutation rate evolution, we have neglected epistatic interactions between selected mutations: in particular, our model does not take into account possible compensatory effects between deleterious alleles (e.g., reciprocal sign epistasis). Allowing the sign of the fitness effect of mutations to depend on the genetic background (which typically occurs in models of directional selection acting on quantitative traits) may affect the selective forces acting on mutation modifier loci. Exploring the evolution of the mutation rate under more realistic assumptions on the genetic architecture of fitness will be the subject of future work.

population, this yields (to the first order in δu and s):

W W ≈ 1 + δu d ln f c du (ζ M,1 + ζ M,2 ) -sh (ζ A,1 + ζ A,2 ) -s (1 -2h) (ζ A,A -D A,A ) -δu d ln f c du sh (ζ MA,∅ + ζ ∅,MA + ζ M,A + ζ A,M -2D MA -2D M,A ) -δu d ln f c du s (1 -2h) [ζ MA,A + ζ A,MA -(ζ M,1 + ζ M,2 ) D A,A -2D MA,A ] . (A6)
The change in frequency of allele 1 at locus M (over one generation) is given by:

∆p M = E W W X M,1 + X M,2 2 -p M = E W W ζ M,1 + ζ M,2 2 (A7)
where E is the average over all individuals just before selection. From equations 2 and A6, and using the fact that repeated indices appearing in genetic associations can be eliminated using the relation D Sii,T = p i q i D S,T + (1 -2p i ) D Si,T (with q i = 1 -p i , e.g., equation 5 in [START_REF] Kirkpatrick | General models of multilocus evolution[END_REF], one arrives at: can be written as F p M q M , where F is the inbreeding coefficient (e.g., [START_REF] Roze | Effects of interference between selected loci on the mutation load, inbreeding depression and heterosis[END_REF], and equation A8 simplifies to: In B, for clonality rates higher than the right-most points (for each population size), deleterious alleles accumulate in the heterozygous state, and the program has to be stopped before U has reached its equilibrium value.

∆p M ≈ δu d ln f c du (p M q M + D M,M ) -sh (D MA + D M,A ) -s (1 -2h) D MA,A -δu d ln f c du sh [(1 -2p M ) (D MA + D M,A ) + 2D MA,M ] -δu d ln f c du s (1 -2h) (D MA,MA -D M,M D A,A ) . (
∆p M ≈ δu d ln f c du (1 + F ) p M q M -sh (D MA + D M,A ) -s (1 -2h) D MA,A . (
A � = ���� � = ���� � = ��� � = � ��� ��� ��� ��� ��� ��� γ ��� ��� ��� ��� � ��� B ��� ��� ��� ��� ��� ��� α ���� ���� ���� ���� � ��� C � = �� � = �� � = � � = � � = � ��� ��� ��� ��� ��� ���

A

  second simulation program considers a different genetic architecture for the mutation rate: instead of being coded by a single locus, U depends on 1000 biallelic loci evenly spaced along the chromosome. Alleles at each of these loci are denoted 0 and 1; genotypes carrying alleles 0 at all loci have mutation rate U min = 0.01. At each locus, allele 1 increases the mutation rate, by an amount that is sampled (independently for each locus) from an exponential distribution with parameter λ. The effects of alleles 1 at the same or at different loci are additive. During the first 2000 generations, all loci affecting the mutation rate are fixed for allele 0; then, during 10 6 generations mutations occur at rate 10 -5 U at each of these loci (mutations and back mutations occur at the same rate). Selection and recombination are implemented as in the previous program. Data availability. The authors state that all data necessary for confirming the conclusions presented in the article are represented fully within the article. Data available from the Dryad Digital Repository: http://dx.doi.org/10.5061/dryad.1j6b0.

  Figures 3E and 3F, yielding a linear selection gradient (d ln f c /dU = a + b U ). Param-

Finally, our model

  makes a number of assumptions on selection against deleterious alleles: in particular, all deleterious alleles have the same selection and dominance coefficient, while drift has no significant effect on their equilibrium frequency. Given the concave shape of the relation between the selection coefficient of deleterious alleles s and the strength of indirect selection acting on the mutation rate, introducing variability in s across loci (while keeping the average constant) should in principle reduce the overall strength of selection for lower mutation rates. However, our deterministic model cannot be used to predict the effect of deleterious alleles for which N e s ∼ 1 (whose frequency is significantly affected by drift), while our infinite sites simulation program cannot deal with very low s values as individuals then carry very large numbers of mutations, causing the program to become extremely slow. Given that an important proportion of mutations may possibly fall in the N e s ∼ 1 parameter region (e.g.,[START_REF] Eyre-Walker | The ditribution of fitness effects of new mutations[END_REF], it would be important to explore the effect of such weakly selected deleterious alleles on the evolution of mutation rate modifiers.

  A8)The first term of equation A8 represents the effect of direct selection acting at locus M (due to the cost of replication fidelity), while the other terms (involving genetic associations) correspond to indirect selection. Associations D MA , D M,A and D MA,A are of order δu, and the term on the second line of equation A8 is thus of order δus.Furthermore, one can show that the association D MA,M is of order s (e.g.,[START_REF] Roze | Effects of interference between selected loci on the mutation load, inbreeding depression and heterosis[END_REF], and the third line of equation A8 can thus be neglected (terms of order (δu) 2 s and δus 2 ). Finally, the term D MA,MA -D M,M D A,A that appears on the last line of equation between loci M and A, generated by partial selfing[START_REF] Roze | Effects of interference between selected loci on the mutation load, inbreeding depression and heterosis[END_REF]. Similarly, one can show that the association D M,M (measuring the excess homozygosity at locus M ) that appears on the first line of equation A8 is affected by the identity disequilibrium through a term proportional to s (1 -2h) (D MA,MA -D M,M D A,A ) -see equation 5 in[START_REF] Roze | Effects of interference between selected loci on the mutation load, inbreeding depression and heterosis[END_REF]. However, we show in the main text that contrarily to D MA,MA -D M,M D A,A , the expressions for D MA , D M,A and D MA,A at QLE do not tend to zero when the frequency of the deleterious allele p A tends to zero. Therefore, assuming that the deleterious allele stays at low frequency (p A small), we may neglect terms proportional to p A , and thus neglect terms involving the identity disequilibrium. In this case, D M,M

Figure 1 .

 1 Figure 1. Strength of indirect selection acting on a mutation modifier locus (scaled

FFigure 2 .FFigure 3 .Figure 4 .FFigure 5 .Figure 6 .

 23456 Figure 2. Evolutionarily stable mutation rate (per haploid genome) as a function of

  lection component acting on the evolution of reproductive systems, for example by accelerating the extinction of selfing or asexual lineages by mutational meltdown if U reaches higher levels due to stronger drift. More generally, it may affect the long-term evolutionary potential of selfing or asexual species, or the relation between the selfing rate and level of inbreeding depression across species. It would thus be of particular interest to obtain mutation rate estimates from different pairs of closely related species with contrasted reproductive systems, in order to see if a general pattern emerges.

Acknowledgements. We thank Jean-Nicolas Jasmin and two anonymous reviewers for helpful discussions and comments, and the bioinformatics and computing service of Roscoff's Biological Station (Abims platform) for computing time. This work was supported by the French Agence Nationale de la Recherche (project TRANS, ANR-11-BSV7-013 and project Clonix, ANR-11-BSV7-007) and by a post-doctoral fellowship from the Conseil Général du Finistère to CG.

APPENDIX: CHANGE IN FREQUENCY AT THE MODIFIER LOCUS

Using the definitions given in the main text, the fitness of an individual can be written as:

where f c represents the cost of replication fidelity and u is the mutation rate of the individual, given by:

Replacing X M,j by ζ M,j + p M and noting that u = u 0 + 2δu p M is the average mutation rate at locus A, we have

), and a Taylor series of f c (u) to the first order in δu yields:

Since f c (u) = f c (u) when δu = 0, equation A4 may also be written as:

Replacing X A,j by ζ A,j + p A in equation A1, one then obtains (to the first order in δu):


Jean-Baptiste Leblond

Pierre Mouro

Crack propagation from a pre-existing flaw at a notch root -II: Detailed form of the stress intensity factors at the initial crack tip and conclusion

Keywords: Notch, pre-existing crack, mixed mode, stress intensity factors, complex potentials, conformal mapping, crack kinking, principle of local symmetry

This paper pursues the study of crack kinking from a pre-existing crack emanating from some notch root. It was shown in Part I that the stress intensity factors at the tip of the small initial crack are given by universal (that is, applicable in all situations, whatever the geometry of the body and the loading) formulae; they depend only on the 'stress intensity factor of the notch' (the multiplicative coefficient of the singular stress field near the apex of the notch in the absence of the crack), the length of the crack, the aperture angle of the notch and the angle between its bisecting line and the direction of the crack. Here we identify the universal functions of the two angles just mentioned which appear in these formulae, by considering the model problem of an infinite body endowed with a notch with straight boundaries and a straight crack of unit length. The treatment uses Muskhelishvili's complex potentials formalism combined with some conformal mapping. The solution is expressed in the form of an infinite series involving an integral operator, which is evaluated numerically. Application of Goldstein and Salganik's principle of local symmetry then leads to prediction of the kink angle of the crack extension. It is found that although the direction of the crack is closer to that of the bisecting line of the notch after kinking than before it, the kink angle is not large enough for the crack tip to get closer to this line after kinking, except perhaps in some special situations.

Introduction

The reader is referred to the Introduction and Fig. 1 of Part I for a general statement of the problem envisaged and notations. It was shown in Section 4 of this first part that the stress intensity factors (SIFs) K p (a) (p = I,II) (where a is recalled to denote the length of the initial crack) at the tip of this crack admit the following asymptotic expression for a → 0 : K p (a) ≡ K (β) p a β + o(a β ).

(1)

In this expression β is an exponent given by

β = α -1 2 , (2)
where α denotes the 'exponent of the displacements' in their asymptotic expression near the notch root in the absence of the crack; α is the smallest positive solution of the equation sin(2αψ) + α sin(2ψ) = 0,

where it is recalled that 2(πψ) represents the aperture angle of the notch. Since α varies from 1 to 1 2 when ψ varies from π/2toπ, β is positive except for ψ = π (case where the notch is itself a crack), so that except in this case, the SIFs K p (a) go to zero with a. Furthermore the quantities K (β) p (p = I,II) are given by the following universal (that is, applicable in all circumstances, whatever the geometry of the body, the notch and the crack and whatever the loading) formula:

K (β) p ≡ L p (ψ, ϕ)κ, (4)
where ϕ is recalled to denote the angle between the bisecting line Ox of the notch and the initial tangent to the crack, the L p (ψ, ϕ) some yet unknown functions of the angles ψ and ϕ and κ the 'mode I stress intensity factor of the notch', that is, the multiplicative coefficient of the singular stress field near the notch root, in the absence of the crack. It is recalled that the definition of the 'stress intensity factor κ of the notch' chosen, implicitly implied by Equations (1) of Part I, possesses the nice properties that κ reduces to the ordinary SIF K I for ψ = π (case of a crack) and to the uniform tensile stress σ yy for ψ = π/2 (case of a regular, corner-free boundary).

It was also shown, in Section V of Part I, that the function L p (ψ, ϕ) (p = I,II) can be identified with the p-th SIF at the crack tip in a 'model problem' involving an infinite body endowed with a notch with straight boundaries and a straight crack originating from the notch root, of unit length (Fig. 1c). The conditions on the notch boundaries, the crack lips and at infinity in this problem read as follows:

                         σ (r, θ =±ψ) • e θ (θ =±ψ) = 0, ∀r>0 (conditions on the notch boundaries) σ (r, θ = ϕ + ).(-e θ (θ = ϕ)) = (r, θ = ϕ).e θ (θ = ϕ)
∀r, 0 <r<1 (condition on the upper lip of the crack)

σ (r, θ = ϕ -).e θ (θ = ϕ) =-(r, θ = ϕ).e θ (θ = ϕ)
∀r, 0 <r<1 (condition on the lower lip of the crack) σ (r, θ) = O(r -α-1 ) for r →+∞, ∀θ,-ψ<θ<+ (condition at infinity).

(5)

In these equations, r and θ are recalled to denote polar coordinates with origin at the notch root, e θ (θ) the unit orthoradial vector and { (r, θ)} the stress field defined by Equations (1) of Part I, except for the omission of the factor κ:

                     rr (r, θ) =- α 2 (2πr) α-1 (3 -α) sin((1 + α)ψ) (1 -α) sin((1 -α)ψ) cos((1 -α)θ) + cos((1 + α)θ) , θθ (r, θ) = α 2 (2πr) α-1 - (1 + α) sin((1 + α)ψ) (1 -α) sin((1 -α)ψ) cos((1 -α)θ) + cos((1 + α)θ) , rθ (r, θ) = α 2 (2πr) α-1 - sin((1 + α)ψ) sin((1 -α)ψ) sin((1 -α)θ) + sin((1 + α)θ) . (6)
This plane elasticity problem is solved here with the aid of [START_REF] Muskhelishvili | Some Basic Problems of the Mathematical Theory of Elasticity[END_REF] complex potentials formalism and conformal mapping. The conformal mapping used maps the lower half-plane of some z-plane onto the physical domain in some Z-plane. To the best of the authors' knowledge, its expression has not previously been available in the literature. It is deduced here from the thesis of [START_REF] Mouchrif | Trajets de propagation de fissures en mécanique linéaire tridimensionnelle de la rupture fragile[END_REF], who has considered the special case where = π . The problem is reduced to an integral equation on Muskhelishvili's potential ϕ ′ (z). The solution is given in the form of an infinite series involving the integral operator appearing in the integral equation. This series is evaluated numerically. Finally, the SIFs and hence the functions L p (ψ, ϕ), are deduced from [START_REF] Andersson | Stress intensity factors at the tips of a star-shaped contour in an infinite tensile sheet[END_REF] formula which relates them to the value of the potential ϕ ′ (z) at that point corresponding to the crack tip in the z-plane.

The results are then compared with those available in the literature. The latter results all pertain to special cases involving a notch with straight boundaries and a straight crack; that is, the universal applicability of the functions L p (ψ, ϕ) has not been perceived. As already mentioned in the Introduction of Part I, four special cases have been previously envisaged. The first one is quite classical and concerns the case where ψ = π/2andϕ = 0 (edge crack originating perpendicularly from a regular free surface). The second one, envisaged by many authors (notably [START_REF] Bilby | The crack with a kinked tip[END_REF][START_REF] Wu | Elasticity problems of a slender Z-crack[END_REF][START_REF] Amestoy | Déviation infinitésimale d'une fissure dans une direction arbitraire[END_REF][START_REF] Amestoy | Crack paths in plane situations -II. Detailed form of the expansion of the stress intensity factors[END_REF], is that where ψ = π but ϕ is arbitrary (semi-infinite crack with an arbitrarily kinked extension). The third one, also considered by several authors [START_REF] Wilson | On Combined Mode Fracture Mechanics[END_REF][START_REF] Sha | Weight function calculations for mixed-mode fracture problems with the virtual crack extension technique[END_REF][START_REF] Aliabadi | Mixed-mode Bueckner weight functions using boundary element analysis[END_REF][START_REF] Fett | Stress Intensity Factors and Weight Functions[END_REF], is that where ψ = π/2butϕ is arbitrary (oblique edge crack emanating from a regular free surface). Finally, the fourth one, which has been treated by [START_REF] Fett | Stress Intensity Factors and Weight Functions[END_REF], using some previous works of [START_REF] Gross | Some Plane Problem Elastostatic Solutions for Plates Having a V-Notch,P h . D .T h e s i s[END_REF], Hasebe and Isida (1978) and Iida (1979), concerns the situation where ϕ = 0butψ is arbitrary (crack lying on the bisecting line of an arbitrary notch). In all of these special cases, we find results in good agreement with those of the literature.

We finally envisage the problem of determining the kink angle χ of the crack extension (see Fig. 1 of Part I). This kink angle is deduced from the values of the SIFs of the initial crack using some suitable criterion, namely [START_REF] Goldstein | Brittle fracture of solids with arbitrary cracks[END_REF] well-known and widely accepted principle of local symmetry. As a consequence of the fact that the notch is always asymptotically predominantly loaded in 'mode I', the angles ϕ and χ are found to be of opposite signs; that is, the direction of the crack is closer to that of the bisecting line of the notch after kinking than before it. It remains to see whether or not the kink angle is large enough for the distance from the crack tip to the bisecting line to decrease after the kink. It is found that it is generally not sufficiently large for this possibility to occur, that is, the crack tip generally still tends to move away from the bisecting line of the notch after the kink, although less markedly than prior to it.

The problem in the physical Z-plane

Let (Z) and (Z) denote Muskhelishvili's potentials in the physical Z-plane. Then, the boundary Z 1 Z 2 Z 3 Z 4 Z 5 (see Fig. 1c) of the body being oriented in that way, the boundary conditions read

i d dS (Z) + Z ′ (Z) + (Z) = T X (Z) + iT Y (Z) on Z 1 Z 2 Z 3 Z 4 Z 5 , (7)
where S denotes the curvilinear length along the arc Z 1 Z 2 Z 3 Z 4 Z 5 and T X (Z), T Y (Z) the components of the traction T(Z) exerted on the boundary1 . Also, since the stresses are O |Z| -α-1 at infinity (Equation (5) 4 ), the conditions at infinity for the potentials (Z), (Z) read

(Z) = O |Z| -α , ( Z ) = O |Z| -α for |Z|→+∞, (8)
which implies that they vanish at infinity (since α is positive). The traction T(Z) being zero on the notch boundaries Z 1 Z 2 and Z 4 Z 5 (Equation (5) 1 ), Equation (7) implies that (Z) + Z ′ (Z) + (Z) is constant along these two half-straight lines. Since both (Z) and (Z) vanish at infinity, one must in fact have

(Z) + Z ′ (Z) + (Z) = 0f o rZ ∈ Z 1 Z 2 ∪ Z 4 Z 5 . (9)
Let us now calculate the traction T(Z) on the lower lip Z 2 Z 3 of the crack. By Equation (5) 3 , one has on this lip

T(Z) =-(r, ϕ).e θ (ϕ) =-rθ (r, ϕ)e r (ϕ) -θθ (r, ϕ)e θ (ϕ) so that T X (Z) + iT Y (Z) = [T r (r, ϕ) + iT θ (r, ϕ)] e iϕ =-[ rθ (r, ϕ) + i θθ (r, ϕ)] e iϕ . A straightforward calculation based on Equations (6) 2 and (6) 3 then leads to T X (Z) + iT Y (Z) = -iαS(2π |Z|) α-1 where S ≡ 1 2 - sin((1 + α)ψ) (1 -α) sin((1 -α)ψ) e iαϕ + αe i(2-α)ϕ + e -iαϕ . (10) Since dS ≡ d|Z| on Z 2 Z 3 , Equation (7) reads d d|Z| (Z) + Z ′ (Z) + (Z) =-αS(2π |Z|) α-1 ,
which yields upon integration

(Z) + Z ′ (Z) + (Z) =- S 2π (2π |Z|) α for Z ∈ Z 2 Z 3 , (11)
the constant of integration being zero since by Equation ( 9), (Z)

+ Z ′ (Z) + (Z) is zero at the point Z 2 ≡ 0 -.
On the upper lip Z 3 Z 4 of the crack, the traction T(Z) is opposite to that on the lower lip Z 2 Z 3 by Equations (5) 2 and (5

) 3 . It follows that T X (Z) + iT Y (Z) = iαS(2π |Z|) α-1 on Z 3 Z 4 .
On the other hand, dS ≡-d|Z| instead of d|Z|, so that Equation ( 7) reads for

Z ∈ Z 3 Z 4 : - d d|Z| (Z) + Z ′ (Z) + (Z) = αS(2π |Z|) α-1 ,
which is the same equation as on Z 2 Z 3 and leads upon integration to the same result:

(Z) + Z ′ (Z) + (Z) =- S 2π (2π |Z|) α for Z ∈ Z 3 Z 4 , (12)
the constant of integration being again zero because (Z) + Z ′ (Z) + (Z) is zero at the point Z 4 ≡ 0 + by Equation ( 9). Equations ( 9), ( 11) and ( 12) can be summarized as follows:

(Z) + Z ′ (Z) + (Z) =- S 2π (2π |Z|) α I(Z) on Z 1 Z 2 Z 3 Z 4 Z 5 , (13)
where I(Z)denotes the characteristic function of the arc Z 2 Z 3 Z 4 :

I(Z) = 1i f Z ∈ Z 2 Z 3 Z 4 , 0i f Z ∈ Z 2 Z 3 Z 4 . ( 14
)
The problem is to find complex potentials (Z), (Z) analytical in the physical domain and satisfying conditions (13) on Z 1 Z 2 Z 3 Z 4 Z 5 and (8) at infinity.

Conformal mapping -the problem in the z-plane

Mouchrif (1994) has provided the expression of the conformal mapping which maps the lower half -of some z-plane onto the exterior of a kinked crack composed of a principal semiinfinite branch and an extension of unit length making an arbitrary angle mπ (-1 <m<+1) with the principal branch in some z-plane (Figs. 1a and1b). The expression of this conformal mapping is

z ≡ ω(z) ≡ A(1 -z 2 ) z + 1 z -1 m ,A ≡ 1 1 -m 2 1 -m 1 + m m . ( 15
) Defining Z ≡ (z) ≡ω(z) ≡ ( ω • ω)(z) ≡z n , ( 16
)
where n is a parameter such that 1 2 ≤ n ≤ 1, we 'open the principal branch of the crack' and thus obtain, in the Z-plane, the geometry we are interested in, with a notch of aperture angle 2(πψ) and a crack of unit length making an angle ϕ with the bisecting line of the notch (Figs. 1b and1c). The a priori free parameters m and n must be chosen in order for the angles and ϕ to take the desired values. It is evident that ψ = nπ; furthermore the point z3 ≡ e imπ in the z-plane becomes z 3 ≡ e imnπ ≡ e iϕ in the z-plane, which implies that ϕ = mnπ. Thus the parameters m and n are given by

n ≡ π ; m ≡ ϕ nπ = ϕ . ( 17
)
We shall need the expression of the quantity (t)/ ′ (t) for t ∈ R in the sequel. After a somewhat lengthy but straightforward calculation, one gets (t)

′ (t) = Q(t)[1 + j(t)], ∀t ∈ R, (18)
where

Q(t) ≡ t 2 -1 2(t -m) (19)
and

j(t) ≡          e -2iψ /n -1i f t<-1 e 2iϕ /n -1i f -1 <t<+1 e +2iψ /n -1i f t>+1. (20)
(Note that the quantity

Q(t)[1 + j(t)
] is continuous at the points t =± 1 in spite of the discontinuities of j(t), since Q(t) vanishes there.)

In the z-plane, Equation (13) takes the form

ϕ(t) + (t) ′ (t) ϕ ′ (t) + ψ(t) =- S 2π (2π | (t)|) α i(t) , ∀t ∈ R
where ϕ(z) ≡ (Z), ψ(z) ≡ (Z) and i(t) denotes the characteristic function of the interval (-1, +1). By Equation ( 18), this equation reads

ϕ(t) + Q(t)ϕ ′ (t) + Q(t)j(t)ϕ ′ (t) + ψ(t) =- S 2π (2π | (t)|) α i(t), ∀t ∈ R. (21)
Furthermore, by Equations ( 8), ( 15) and ( 16),

ϕ(z) = O |z| -2nα ,ψ ( z ) = O |z| -2nα for |z|→+∞. ( 22
)
The problem is reduced to finding some complex potentials ϕ(z), ψ(z) analytical in the lower half-plane -and satisfying conditions (21) on R and ( 22) at infinity.

Solution of the problem

REDUCTION OF THE PROBLEM TO AN INTEGRAL EQUATION

Let us define the function

χ(z) ≡ 1 2iπ +∞ -∞ Q(t)j(t)ϕ ′ (t) dt t -z ≡ 1 4iπ +∞ -∞ (t 2 -1)j (t)ϕ ′ (t)dt (t -m -)(t -z) , ∀z ∈ -(23)
where m -≡ miǫ, ǫ > 0,ǫ → 0. (Note that the integrand behaves like |t| -2nα-1 for t →± ∞ ;hence, since n and α are positive, the integral is convergent at infinity. It appears essential here to have removed the tractions imposed from infinity to the crack lips; otherwise the integrand would behave like |t| 2nα-1 for t →± ∞ , and the integral would diverge at infinity). By Plemelj's formula,

χ(t + ) -χ(t -) = Q(t)j(t)ϕ ′ (t) , ∀t ∈ R
where t ± ≡ t ± iǫ, ǫ > 0,ǫ→ 0. Also, let us define the function

γ(z)≡ 1 2iπ +1 -1 S 2π (2π | (t)|) α dt t -z , ∀z ∈ -. (24)
Then Plemelj's formula again yields

γ(t + ) -γ(t -) = S 2π (2π | (t)|) α i(t), ∀t ∈ R .
These elements and the obvious property Q(t) = Q(t), ∀t ∈ R, imply that Equation ( 21) may be written as

ϕ(t) + Q(t)ϕ ′ (t) + χ(t + ) -χ(t -) + ψ(t) + γ(t + ) -γ(t -) = 0 , ∀t ∈ R. Now, f(z) ≡ f(z) denoting the conjugate of f(z) in the sense of Muskhelishvili, one has χ(t + ) = χ(t -) = χ(t -) and similarly γ(t + ) = γ(t -) = γ(t -)
, so that the preceding equation may be re-written as

ϕ(t) -χ(t -) -γ(t -) =-Q(t)ϕ ′ (t) -χ(t -) -ψ(t) -γ(t -), ∀t ∈ R. (25)
We shall now make use of the following classical result, the (elementary) proof of which is given for instance in [START_REF] Mouchrif | Trajets de propagation de fissures en mécanique linéaire tridimensionnelle de la rupture fragile[END_REF] thesis:

Lemma. Let f(z)and g(z) denote two functions defined and continuous on -∪ R, analytic on -, bounded at infinity and such that f(t) = g(t), ∀t ∈ R. Then f(z) and g(z) are constant and conjugate to each other.

It is straightforward, though somewhat tedious, to check that all the hypotheses of this lemma are verified by the left-and right-hand sides of Equation ( 25). (Here again, it appears as essential to have removed the remote tractions from infinity to the lips of the crack; otherwise ϕ(z) would be O(|z| 2nα ) for |z|→+∞so that the left-hand side of Equation ( 25) would not be bounded at infinity). Hence the lemma yields

ϕ(z) -χ(z) -γ(z)= 0 ⇒ ϕ(z) = χ(z) + γ(z), ∀z ∈ -, (26)
the additive constant being zero in the present case since by Equations (22) 1 , (23) and ( 24), all functions ϕ(z), χ(z), γ (z) vanish at infinity. Differentiating Equation (26) 2 and accounting for the expressions (23) 2 , (24) of the functions χ(z) and γ(z),one gets the following integral equation on the unknown potential ϕ ′ (z) : ϕ

′ (z) = 1 2iπ +1 -1 S 2π (2π | (t)|) α dt (t -z) 2 + 1 4iπ +∞ -∞ (t 2 -1)j (t)ϕ ′ (t)dt (t -m -)(t -z) 2 , ∀z ∈ -, (27)
where, by Equations ( 15) and ( 16),

| (t)|=A n (1 -t 2 ) n 1 + t 1 -t mn , ∀t ∈ (-1, +1). (28)
The function j(t) is constant over each interval (-∞, -1), (-1, +1), (+1, +∞) (see Equation ( 20)). Furthermore, for any real t, ϕ ′ (t) = ϕ ′ (t); hence ϕ ′ (t) may be replaced by ϕ ′ (t) in the second integral of the right-hand side of Equation ( 27). Provided that this replacement is made, the integrand corresponding to each interval (-∞, -1), (-1, +1), (+1, +∞) is then analytic over the upper half-plane + . Hence these intervals may be displaced within + : see Fig. 2. Thus the interval (-∞, -1) may be replaced by the half-straight line + ℓ of equation Re(λ) =-1, Im(λ) > 0, oriented from -1 + i∞ to -1; indeed the integral over the great quarter of a circle which closes the integration path is zero since the integrand is O(|λ| -2-2nα ). Similarly, the interval (+1, +∞) may be replaced by the half-straight line + r : Re(λ) =+1, Im(λ) > 0, oriented from +1to+1 + i∞. Finally, the interval (-1, +1) may be replaced by the semi-circle Ŵ + :| λ|=1, Im(λ) > 0, oriented from -1to+1(since the pole at m -is not crossed in the process). Hence the second integral in the right-hand side of Equation ( 27) may be written as

Lϕ ′ (z) ≡ 1 4iπ C + (λ 2 -1)j (λ)ϕ ′ (λ)dλ (λ -m)(λ -z) 2 ,j ( λ ) ≡        e -2iψ /n -1o n + ℓ e 2iϕ /n -1o nŴ + e +2iψ /n -1o n + r , , (29)
where C + denotes the path + ℓ ∪ Ŵ + ∪ + r , oriented like + ℓ ,Ŵ + and + r . Similarly, taking advantage of the fact that the function (1λ 2 ) nα ((1 + λ) / (1λ)) mnα is analytic on + , one may write the first integral of the right-hand side of Equation ( 27), using Equation ( 28), in the form

ϕ ′ 0 (z) ≡ 1 2iπ Ŵ + S 2π 2πA n (1 -λ 2 ) n 1 + λ 1 -λ mn α dλ (λ -z) 2 . ( 30
)
The integral equation ( 27) on the potential ϕ ′ (z) can therefore be written as

ϕ ′ (z) = ϕ ′ 0 (z) + Lϕ ′ (z) , ∀z ∈ - (31)
where the integral operator L and the function ϕ ′ 0 (z) are defined by Equations ( 29) and (30).

SOLUTION IN THE FORM OF AN INFINITE SERIES

The solution of the integral equation ( 31) is unique and can be written as

ϕ ′ (z) = +∞ p=0 L p ϕ ′ 0 (z), (32)
where L p denotes the p-th iterate of the integral operator L defined by Equation ( 29), provided that this operator is contractant, that is, that it reduces some norm (for instance the ∞ norm) of the function it is applied on. Since it involves the function j(λ) (defined on C + by Equation (29) 2 ), the contractant character will be ensured by the smallness of this function. Now j(λ) is small when ψ is close to π (sharp notch) and ϕ close to 0 (crack lying almost along the bisecting line of the notch); indeed, by Equation (17) 1 ,nisthencloseto1sothat e -2iψ /n -1,e 2iϕ /n -1,e +2iψ /n -1 are all small. Hence the convergence of the series ( 32) is guaranteed under these conditions. Numerical experience (see Section 5 below) confirms that convergence can indeed be achieved except when ψ is close to π/2 ('open notch') and ϕ large (crack lying away from the bisecting line).

In the integral defining the operator L (see Equation (29) 1 ), the only values of the function ϕ ′ (z) it acts on which are involved are those on the conjugate C -of the arc C + (see Fig. 2). Therefore the values of the successive iterates ϕ ′ 0 (z), Lϕ ′ 0 (z), L 2 ϕ ′ 0 (z), etc. can be calculated (numerically) inductively on the sole arc C -, and the values of the unknown potential ϕ ′ (z) on the same arc follow from there and application of Equation ( 32). The value of ϕ ′ (z) at any point z ∈ -can then be obtained from Equation (31).

STRESS INTENSITY FACTORS

The values of the SIFs at the crack tip K p ≡ L p (ψ, ϕ) (p = I,II) can finally be deduced from [START_REF] Andersson | Stress intensity factors at the tips of a star-shaped contour in an infinite tensile sheet[END_REF] formula, which involves the value of the potential ϕ ′ at that point z 3 ≡ m in the z-plane corresponding to the crack tip Z 3 ≡ e iϕ in the physical Z-plane (see Fig. 1):

K I -iK II ≡ L I (ψ, ϕ) -iL II (ψ, ϕ) = 2 √ πe -iϕ/2 [ ′′ (m)] -1/2 ϕ ′ (m).
The value of ϕ ′ (m) here can be obtained as explained above. Calculation of the quantity ′′ (m) leads to the following final expression:

L I (ψ, ϕ) -iL II (ψ, ϕ) =-i 2π n (1 -m 2 )e -iϕ ϕ ′ (m) . ( 33
)
Table 1. Va l u e o f L I as a function of ψ (1st row) and ϕ (1st column)

= 90 • 100 • 110 • 120 • 130 • 140 • 150 • 160 • 170 • ϕ = 0 • 1
= 90 • 100 • 110 • 120 • 130 • 140 • 150 • 160 • 170 • ϕ = 0 • 0 00000000 10 • 0.

Numerical calculation of the functions L p (ψ, ϕ)

In order to perform the numerical integrations, the arc C + is divided into 400 intervals: 100 on each half-straight line + ℓ , + r (numerical integration being 'cut off' on these half-lines at a distance of 100 from the x-axis) and 200 on the semi-circle Ŵ + . The 'mesh' is refined near the singular points ±1 and enlarged near 'infinity' on the half-straight lines + ℓ , + r . Gaussian integration is used, with 2 Gauss points in each interval (which makes a total of 800 Gauss points). The values of the successive iterates ϕ ′ 0 (z), Lϕ ′ 0 (z), L 2 ϕ ′ 0 (z), etc., are calculated at the conjugates (lying on C -) of the Gauss points of the arc C + . The number of terms in the series (32) which is necessary to reach convergence is highly dependent upon the values of the angles ψ and ϕ : it varies from just a few ones when ψ is close to π and ϕ to0toseveral hundreds when ψ is closer to π/2andϕ larger.

Tables 1 and2 show the results derived in that way for the functions L p (ψ, ϕ) (p = I,II), for 90 • ≤ ψ ≤ 170 • and 0 • ≤ ϕ ≤ 70 • . Blank spaces correspond to those cases (involving values of ψ close to 90 • and large values of ϕ) where satisfactory convergence of the series (32) could not be achieved. Results corresponding to negative values of the angle ϕ can also be obtained using the fact that for obvious symmetry reasons, L I (ψ, ϕ) and L II (ψ, ϕ) are even and odd functions, respectively, of this argument. The (absolute) accuracy is better than 10 -2 , and in many cases than 10 -3 .

Of course, these results could also be obtained through more conventional techniques, such as the finite element method. The advantages of the approach used here are three-fold. First, the numerical accuracy which can be achieved is better than that currently encountered when using the finite element method, except in those cases where convergence of the series (32) is difficult to obtain, that is, when ψ is close to π/2andϕ large. Second, the computation time is shorter, except in the same cases as before. Third, changing the values of ψ and ϕ only requires to modify the values of the input parameters of the program instead of demanding a new mesh to be designed for each case. On the other hand, the versatility of the finite element method is obviously superior since it can deal with finite bodies; but consideration of an infinite one is sufficient to the values of the functions L I (ψ, ϕ) and L II (ψ, ϕ) looked for here.

Partial checks on the results

As mentioned in the Introduction, there are four special cases which have been treated by former authors, thus allowing for a comparison with our results.

The first one corresponds to the situation where ψ = π/2andϕ = 0 (edge crack emanating perpendicularly from a regular, corner-free surface). It that case it is well known that K I (a) ≃ 1.12 σ ∞ yy √ πa where σ ∞ yy denotes the remote tensile stress (K II (a) being zero). Now for this value of ψ, the 'exponent of the displacements' α in the absence of the crack is unity; hence, by Equation ( 2), β = 1 2 and Equations ( 1) and ( 4) read K I (a)

≡ L I (ψ = π/2,ϕ = 0) κ √ a + o √ a = L I (ψ = π/2,ϕ = 0) σ ∞ yy √ a + o √
a since the definition of the 'stress intensity factor κ of the notch' in the absence of the crack adopted here is such that it reduces to the uniform tensile stress σ yy for ψ = π/2 : see the Introduction. Comparison between these two expressions of K I (a) yields L I (ψ = π/2,ϕ = 0) ≃ 1.12 √ π ≃ 1.985, which is in excellent agreement with the value of 1.99 found here (see Table 1).

The second special case is that where ψ = π but ϕ is arbitrary (semi-infinite crack with a arbitrarily deviated extension). The comparison can be made only for an initial crack loaded in pure mode I; indeed this initial crack is considered here as the limit of a sharp notch for → π, and any notch (with ψ<π)is always asymptotically loaded in pure mode I (see Section 2 of Part I). [START_REF] Leblond | Crack paths in two-dimensional situations -I. General form of the expansion of the stress intensity factors[END_REF] has shown that in this case

K p (a) ≡ F p,I (ϕ)K I + O( √ a) (p = I,II), (34)
where K I denotes the mode I SIF of the initial crack, in the absence of the extension. This formula is of universal value, that is, it applies in all circumstances, whatever the far geometry of the body and whatever the loading. Furthermore, [START_REF] Amestoy | Crack paths in plane situations -II. Detailed form of the expansion of the stress intensity factors[END_REF] have provided formulae for the functions F p,I (ϕ) which allow for their evaluation with an accuracy better than 10 -6 . Now for ψ = π, α = 1 2 so that, by Equation (2), β = 0. Hence Equations ( 1) and (4) read

K p (a) ≡ L p (ψ = π, ϕ)κ + o(a 0 ) = L p (ψ = π, ϕ)K I + o(a 0 ), (35)
where use has been made of the fact that the 'stress intensity factor κ of the notch' reduces to the ordinary SIF K I in the particular case where the notch becomes a crack (ψ = π) (see the Introduction). Comparison between formulae (34) and ( 35) shows that one must have intensity factor κ of the notch' differing from that adopted here). Table 5 compares our results with those of Fett and Munz; a good agreement can again be observed.

L p (ψ = π, ϕ) ≡ F p,I (ϕ) , ∀ϕ( p = I,II) . ( 36
)

Conclusion: prediction of the kink angle of the crack extension

Let us finally consider the problem of determining the beginning of the propagation path of the crack, and more specifically the kink angle χ of the crack extension (see Fig. 1 of Part I). When this extension is still infinitesimal, the SIFs K * p (p = I,II) at its tip are given by the universal formula (where Einstein's implicit summation convention over repeated indices is used)

K * p ≡ F p,q (χ)K p (a) (p, q = I,II) (34')
which was established by [START_REF] Leblond | Crack paths in two-dimensional situations -I. General form of the expansion of the stress intensity factors[END_REF] and is in fact a mere extension of formula (34) to the case where some mode II is initially present. The kink angle χ can then be deduced from [START_REF] Goldstein | Brittle fracture of solids with arbitrary cracks[END_REF] principle of local symmetry. This criterion is well-known and widely accepted and convincing arguments in its favour have been provided by [START_REF] Amestoy | Propagation de fissures en élasticité plane[END_REF]. It stipulates that

K * II = F II,I (χ)K I (a) + F II,II (χ)K II (a) = 0 ; (37)
the numerical values of the functions F p,q (χ) to be used when solving this equation with respect to χ can be deduced for instance from [START_REF] Amestoy | Crack paths in plane situations -II. Detailed form of the expansion of the stress intensity factors[END_REF] approximate but highly accurate formulae for these functions. If a is small, formulae (1) and (4) for the SIFs K p (a) can be used and the preceding equation reduces to This equation determines the kink angle χ as a function of the sole angles ψ and ϕ, independently of the loading. This feature, which may appear surprising at first sight, arises from the fact that as mentioned above, the notch is always asymptotically loaded in pure mode I, whatever the loading. Precisely because of that property, it is intuitively obvious (and is in fact a rigorous consequence of Equation (37 ′ )) that the angles ϕ and χ are of opposite signs, which means that the direction of the crack is closer to that of the bisecting line of the notch after the kink than before it. However, two possibilities remain (Fig. 3): ϕ being for instance taken as positive, either (a) |χ| <ϕ,which means that the distance between the crack tip and the bisecting line of the notch still increases after the kink, although less markedly than before it; or (b) |χ| >ϕ, so that this distance decreases after the kink.

Figure 4 shows the ratio |χ|/ϕ (ϕ being taken as positive) as a function of the angle ψ for various values of the angle ϕ. It can be observed that in most cases, possibility (a) prevails: the distance from the crack tip to the bisecting line still increases after the kink.H o w e v e r , possibility (b) seems to occur for values of ψ and ϕ close to π/2 and 0, respectively; then this distance decreases after the kink. It may nevertheless be remarked that the ratio |χ|/ϕ is close to unity then, and that conditions for good convergence of the series (32) are not well fulfilled, since ψ is close to π/2. Hence the numerical accuracy may not be sufficient to warrant that the ratio |χ|/ϕ is really greater than 1. Anyway, since it is close to unity then, and since in addition ϕ is small, the decrease of the distance from the crack tip to the bisecting line is marginal: the crack extension almost follows the direction of the bisecting line.

Figure 1 .

 1 Figure 1. Conformal mapping: the z-plane (a); the z-plane (b); the (physical) Z-plane (c).

Figure 2 .

 2 Figure 2. Figure 2. Integration paths in the z-plane.

F

  II,I (χ)L I (ψ, ϕ) + F II,II (χ)L II (ψ, ϕ) = 0 . (37')

Figure 3 .

 3 Figure 3. Figure 3. The two possible situations after the kink: |χ| <ϕ(a); |χ| >ϕ(b).

Figure 4 .

 4 Figure 4.Figure 4. The ratio |χ|/ϕ as a function of the angle ψ, for various values of the angle ϕ.

Figure 4 .

 4 Figure 4.Figure 4. The ratio |χ|/ϕ as a function of the angle ψ, for various values of the angle ϕ.

Table 2 .

 2 Va l u e o f L II as a function of ψ (1st row) and ϕ (1st column)

		.99	1.29	0.99	0.86	0.82	0.83	0.87	0.93	0.98
	10 •	1.95	1.27	0.97	0.85	0.81	0.82	0.86	0.92	0.97
	20 •	1.83	1.20	0.93	0.81	0.78	0.79	0.84	0.89	0.94
	30 •	1.65	1.10	0.86	0.76	0.73	0.75	0.79	0.84	0.89
	40 •	-	0.98	0.78	0.69	0.67	0.69	0.73	0.78	0.82
	50 •	-	-	0.68	0.61	0.60	0.62	0.66	0.70	0.74
	60 •	-	-	-	0.53	0.52	0.54	0.58	0.62	0.65
	70 •	-	-	-	-	0.44	0.46	0.50	0.53	0.56

Table 5 .

 5 Comparison between the functions L p (ψ, ϕ = 0) found here and provided by[START_REF] Fett | Stress Intensity Factors and Weight Functions[END_REF]

		L I (ψ, ϕ = 0)L I (ψ, ϕ = 0) (this work) (Fett and
			Munz, 1997)
	= 120 • 0.856 140 • 0.829	0.868 0.811
	160 •	0.931	0.928
	180 •	11

There is a difference of sign here with respect to Muskhelishvili's classical formula, which arises from the fact that the unit normal exterior vector to the boundary is oriented to the left, instead of to the right as usual, of the unit tangent vector to the boundary.

Acknowledgement

The authors wish to express their sincere thanks to one referee who very rightly pointed out the inaccuracy of the extrapolation made in the first version of the paper of the data contained in Tables 52 and53 of Fett and Munz' (1997) book.

Table 3. Comparison between the functions L p (ψ = π, ϕ) obtained here and the functions F p,I (ϕ) provided by [START_REF] Amestoy | Crack paths in plane situations -II. Detailed form of the expansion of the stress intensity factors[END_REF] L I (ψ = π, ϕ) F I,I (ϕ)

L II (ψ = π, ϕ) F II,I (ϕ) (this work) (Amestoy and (this work) [START_REF] Amestoy | Crack paths in plane situations -II. Detailed form of the expansion of the stress intensity factors[END_REF] Leblond, 1992)

(this work) [START_REF] Fett | Stress Intensity Factors and Weight Functions[END_REF] Munz, 1997) Table 3 compares the functions L p (ψ = π, ϕ) obtained here with the functions F p,I (ϕ) resulting from the work of [START_REF] Amestoy | Crack paths in plane situations -II. Detailed form of the expansion of the stress intensity factors[END_REF]. It can be seen that the discrepancy between the results is smaller than the numerical accuracy in the calculation of the functions L p (ψ = π, ϕ).

The third particular case corresponds to the situation where ψ = π/2b u tϕ is arbitrary (oblique edge crack originating from a regular, corner-free surface). The results of several authors [START_REF] Wilson | On Combined Mode Fracture Mechanics[END_REF][START_REF] Sha | Weight function calculations for mixed-mode fracture problems with the virtual crack extension technique[END_REF][START_REF] Aliabadi | Mixed-mode Bueckner weight functions using boundary element analysis[END_REF] have recently been compiled and critically compared by [START_REF] Fett | Stress Intensity Factors and Weight Functions[END_REF]. These authors provide the values of the functions L p (ψ = π/2,ϕ) in the form of tables (Nos. 52 and 53 in their book). (In fact, the SIFs are given for various crack lengths a, and one must extrapolate the results for a → 0). Table 4 shows the comparison between the functions L p (ψ = π/2,ϕ) obtained here and those provided by these authors; the agreement is quite good.

The last special case is that where ϕ = 0b utψ is arbitrary (crack lying on the bisecting line of an arbitrary notch). This case has again been considered by [START_REF] Fett | Stress Intensity Factors and Weight Functions[END_REF], on the basis of previous works of [START_REF] Gross | Some Plane Problem Elastostatic Solutions for Plates Having a V-Notch,P h . D .T h e s i s[END_REF], Hasebe and Isida (1978) and Iida (1979). Fett and Munz' results are given in the form of an approximate but accurate formula for the function L I (ψ, ϕ = 0)(L II (ψ, ϕ = 0) being zero) (account being taken of a definition of the 'stress