François Bignonnet
email: francois.bignonnet@enpc.fr

Luc Dormieux
email: luc.dormieux@enpc.fr

Éric Lemarchand

Eric Lemarchand
email: eric.lemarchand@enpc.fr

Strength of a matrix with elliptic criterion reinforced by rigid inclusions with imperfect interfaces

Keywords: Strength, Homogenization, Interfaces, Micromechanics, Mohr-Coulomb, Tresca

HAL is

Introduction

Imperfect interfaces between constituents of heterogeneous media are increasingly understood to play a major role on the effective strength properties. Such imperfect interfaces can be characterized by a criterion the stress vector acting on the interface must not exceed. Thanks to homogenization methods, the strength of granular geomaterials has been investigated, successively considering the cases of rigid grains interfaced by a Tresca criterion [START_REF] Dormieux | Résistance d'un polycristal avec interfaces intergranulaires imparfaites[END_REF], a frictional criterion [START_REF] Maalej | Micromechanical approach to the failure criterion of granular media[END_REF] or a cohesive frictional criterion [START_REF] He | Cohesive mohr-coulomb interface effects on the strength criterion of materials with granular-based microstructure[END_REF], as well as the competition between Tresca interfacial and Von Mises intra-granular strength [START_REF] Dormieux | Micromechanical models of the strength of a sandstone[END_REF].

Additionally, several types of porous media have been recognized to be governed by an elliptic effective strength criterion in the mean-deviatoric stress plane. For example, a composite made of pores in a matrix with a Von Mises [START_REF] Barthélémy | Approche micromécanique de la rupture et de la fissuration dans les géomatériaux[END_REF] or Drucker-Prager [START_REF] Barthélémy | Approche micromécanique de la rupture et de la fissuration dans les géomatériaux[END_REF][START_REF] Maghous | Micromechanical approach to the strength properties of frictional geomaterials[END_REF] strength criterion has an elliptic macroscopic strength criterion. Above a critical porosity threshold, the granular material with cohesive frictional interfaces considered in [START_REF] He | Cohesive mohr-coulomb interface effects on the strength criterion of materials with granular-based microstructure[END_REF] also proves to follow a similar criterion. Elliptic strength criteria are thus of great interest in geomechanics, for instance to describe the clay matrix of a shale. However, it is worth noting that this class of elliptic strength criteria is obtained only by applying the so-called modified secant modulii approach [START_REF] Suquet | Overall properties of nonlinear composites: a modified secant moduli approach and its link with Ponte-Castaneda's nonlinear variational procedure[END_REF] to the homogenization of ductile porous materials, which is precisely the method we intend to adopt in the present work for a second homogenization step. A Gurson-type analysis [START_REF] Gurson | Continuum theory of ductile rupture by void nucleation and growth. Part I: Yield criteria and flow rules for porous ductile media[END_REF] of the strength of porous materials would have led to the use of a different criterion for the clay matrix, but will not be considered here. micro meso macro rigid grains rigid inclusions r micro r meso

Figure 1: Example of thought-of material. Cohesive frictional granular material reinforced by rigid particles with imperfect interfaces. According to the separation of scale principle, r micro ≪ r meso . This paper deals with upscaling from meso to macro scales.

Driven by the thought-of example of shales for which the clay matrix is reinforced by silica or calcite inclusions (Fig. 1), the macroscopic strength of a material made of a matrix with an elliptic strength criterion reinforced by rigid inclusions is yet to be investigated. Furthermore, the degradation of this reinforcement, coming from matrix-inclusion interface imperfections, is of critical interest. The case of perfect interfaces will be compared to Tresca or Mohr-Coulomb interfaces. The present work addresses the imperfection of the interfaces under the hypothesis of ductility, whereas matrix-inclusion debonding had previously been studied in the context of fracture mechanics (see e.g. [START_REF] Mantič | Interface crack onset at a circular cylindrical inclusion under a remote transverse tension. Application of a coupled stress and energy criterion[END_REF] and [START_REF] Greco | A two-scale failure analysis of composite materials in presence of fiber/matrix crack initiation and propagation[END_REF]).

This paper aims at addressing this strength issue using continuum micromechanics. From a technical point of view, [START_REF] Barthélémy | A micromechanical approach to the strength criterion of drucker-prager materials reinforced by rigid inclusions[END_REF] developed a strength homogenization method in which the macroscopic stress states lying on the boundary of the effective strength criterion domain are obtained by solving a fictitious non linear viscous problem. In turn, the non linear problem can be solved using secant or affine methods [START_REF] Suquet | Overall properties of nonlinear composites: a modified secant moduli approach and its link with Ponte-Castaneda's nonlinear variational procedure[END_REF][START_REF] Suquet | Effective behavior of nonlinear composites[END_REF] which rely on the solution to the associated linear problem. These methods proved successful to predict the strength of heterogeneous material, even in the presence of interface effects [START_REF] Barthélémy | Approche micromécanique de la rupture et de la fissuration dans les géomatériaux[END_REF][START_REF] Barthélémy | A micromechanical approach to the strength criterion of drucker-prager materials reinforced by rigid inclusions[END_REF][START_REF] Dormieux | Micromechanical models of the strength of a sandstone[END_REF][START_REF] Dormieux | Microporomechanics[END_REF][START_REF] Dormieux | Résistance d'un polycristal avec interfaces intergranulaires imparfaites[END_REF][START_REF] He | Cohesive mohr-coulomb interface effects on the strength criterion of materials with granular-based microstructure[END_REF][START_REF] Maalej | Micromechanical approach to the failure criterion of granular media[END_REF][START_REF] Maghous | Micromechanical approach to the strength properties of frictional geomaterials[END_REF][START_REF] Sanahuja | Résistance d'un milieu poreux à phase solide hétérogène[END_REF].

To start with, the homogenization method is briefly recalled and the fictitious non linear problem is derived from the strength properties of the components in section 2.

Next, the non linear homogenization method is presented in section 3 after resolution of the linear problem underlying the fictitious problem arising from the previous section.

Finally, the macroscopic strength criteria are derived in section 4 in the case of perfect, Tresca or Mohr-Coulomb interfaces.

Notations The second and fourth order identity tensors are respectively denoted by 1 and I. The volumic and deviatoric projection tensors J and K are defined as J = 1 3 1 ⊗ 1 and K = I -J.

Limit state equations

Methodology

The aim of this article is to determine the effective strength of a composite made of a matrix reinforced by rigid inclusions with imperfect matrix-inclusion interfaces. A representative elementary volume (rev ) Ω of this composite is introduced. It comprises a matrix (phase Ω m ) and rigid inclusions (phase Ω i ) with volume fraction ρ. Locally, the unit normal to the interface directed outwards from the inclusion is noted n. The volume averages of a function a over the rev Ω, the matrix phase Ω m , the inclusionary phase Ω i , and all the interfaces Γ are respectively denoted a Ω = a, a m , a i and a Γ .

A macroscopic stress state Σ applied to Ω is said to be admissible provided that there exists some microscopic stress field σ(z) defined on the rev Ω which meets the following conditions [START_REF] De Buhan | Approche fondamentale du Calcul à la Rupture des ouvrages en sols renforcés[END_REF], [START_REF] Suquet | Analyse limite et homogénéisation[END_REF]):

divσ = 0 (∀z ∈ Ω) σ(z) ∈ G(z) (∀z ∈ Ω) σ = Σ (1)
where G(z) denotes the domain of admissible microscopic stress states at point z in the rev Ω. The set of admissible macroscopic stress states is denoted by G hom .

The strength properties of the constituents at the microscopic scale therefore need to be characterized. The domain G m of admissible stress states in the matrix is defined by a strength criterion f m (σ) such that:

σ ∈ G m ⇔ f m (σ) 0.
In turn, the inclusions are supposed infinitely resistant. The emphasis of this paper is put on the limited strength of the matrix-inclusions interfaces: The inclusions may be not perfectly bounded to the matrix; instead, the strength of the interface is described by a criterion on the stress vector T acting on the matrixinclusion interface:

T ∈ G Γ ⇔ f Γ (T ) 0.
Γ denotes the set of matrix-inclusion interfaces and G Γ defined above is the set of admissible stress vectors.

Equivalently, the domains G m and G Γ may be characterized by their support functions [START_REF] Dormieux | Microporomechanics[END_REF][START_REF] Salencon | An introduction to the yield design theory and its applications to soil mechanics[END_REF]. This is the so-called dual formulation. The support function π m (d) of the matrix strength is defined as

π m (d) = sup{σ : d, f m (σ) 0},
where d and π m (d) physically represent a virtual strain rate and the associated dissipation. Likewise, the support function π Γ ( v ) of the interface criterion is

π Γ ( v ) = sup{T : v , f Γ (T ) 0},
where v and π Γ ( v ) physically represent a virtual velocity jump across the interface and the associated dissipation.

The direct use of the definition (1) for the determination of G hom is uneasy. Alternatively, as shown in [START_REF] Leblond | Exact results and approximate models for porous viscoplastic solids[END_REF], [START_REF] Barthélémy | A micromechanical approach to the strength criterion of drucker-prager materials reinforced by rigid inclusions[END_REF], the boundary ∂G hom of the macroscopic strength criterion can be retrieved by solving the following boundary value problem defined on the rev

Ω divσ = 0 (Ω) σ = ∂π m ∂d (Ω m ) T = ∂π Γ ∂ v (Γ) σ = C i : d (Ω i ) with C i → ∞ v(z) = D • z (∂Ω) d = grad s v (Ω) (2)
This boundary value problem may be interpreted as a fictitious non linear viscous problem where σ is the local stress field and v the local velocity field. The important feature is that the macroscopic stress state Σ, related to σ by the stress averaging rule Σ = σ lies on the boundary ∂G hom of the macroscopic strength criterion. More precisely it is located at the point of ∂G hom where the outward normal is oriented along the macroscopic strain rate D. Thus, the determination of the boundary ∂G hom of the macroscopic strength criterion amounts to solving the above fictitious non linear viscous problem. The non linear viscous state equations in each phase are defined from the derivatives of their support function. The next parts are dedicated to the characterization of the strength properties of the matrix and the interfaces and to the derivation of the associated fictitious state equations.

Matrix

Using appropriate non linear homogenization techniques, [START_REF] He | Cohesive mohr-coulomb interface effects on the strength criterion of materials with granular-based microstructure[END_REF] recently derived the macroscopic strength of a granular medium made of rigid spherical grains interfaced by cohesive frictional interfaces. The latter follow a Mohr-Coulomb failure criterion. It will be assumed that the result of this analysis provides a suitable description for the strength of the matrix at stake. The strength criterion is isotropic and can be expressed as follows:

f m (σ) = (σ m + c) 2 a + σ 2 d b -1 0 (3)
where a, b, c depend on intergranular interfacial strength properties and on the porosity. The stress invariant appearing in (3) are the mean stress σ m and the deviatoric stress σ d , which are defined as

σ m = 1 3 tr(σ) ; σ d = σ -σ m 1 ; σ d = √ σ d : σ d .
In the (σ m , σ d ) plane, the boundary of G m is an ellipse with semi-axes √ a and √ b. The support function of the failure criterion (3) is the function of the tensorial variable d

π m (d) = ad 2 v + bd 2 d -cd v (4)
where the strain rate invariants are the volume strain rate d v and the deviatoric strain rate d d , defined as

d v = tr(d) ; d d = d - 1 3 d v 1 ; d d = d d : d d .
Differentiating the support function (4), the fictitious state equation in the matrix is written in a secant formulation involving a prestress

σ = ∂π m ∂d = C m : d + σ p , with C m = 3k m J + 2µ m K ; σ p = -c1, k m = a ad 2 v + bd 2 d ; 2µ m = b ad 2 v + bd 2 d .
(5)

Tresca interfaces

First, the case where the strength of the interfaces is governed by a Tresca criterion f tr (T ) is considered. The Tresca criterion is characterized by a threshold k on the tangential component T t of the stress vector:

f tr (T ) = T t -k 0, (6)
with

T t = (1 -n ⊗ n) • T ; T t = T t • T t ; T n = T • n.
The support function of the Tresca criterion [START_REF] Salencon | An introduction to the yield design theory and its applications to soil mechanics[END_REF] is

π tr ( v ) = +∞ if v n = 0, k v t if v n = 0. ( 7
)
This support function is highly singular, thus not suited for the direct implementation from the definition of the fictitious state equation. Following [START_REF] Barthélémy | Approche micromécanique de la rupture et de la fissuration dans les géomatériaux[END_REF]; [START_REF] Sanahuja | Résistance d'un milieu poreux à phase solide hétérogène[END_REF], the Tresca criterion can be viewed as the limit of a sequence of ellipses in the (T n , T t ) plan

f tr (T , L) = T n L 2 + T t k 2 -1 0, (8)
when L/k tends to infinity. The support function of the regularized criterion ( 8) is

π tr ( v , L) = (L v n ) 2 + (k v t ) 2 . ( 9
)
The associated state equation reads

T = ∂π tr ( v , L) ∂ v = K( v , L) • v = K n ( v , L) v n n + K t ( v , L) v t , (10)
where the normal and tangential components of the velocity jump v across the interface are defined similarly to those of the stress vector T . The secant modulii appearing in (10) are:

K n ( v , L) = L 2 (L v n ) 2 + (k v t ) 2 ; K t ( v , L) = k 2 (L v n ) 2 + (k v t ) 2 . ( 11
)
In the limit L/k → ∞, the normal velocity jump v n → 0 while tangential velocity jumps are allowed. Then, the tangent stiffness modulus K t → k/ v t .

Mohr-Coulomb interfaces

In order to account for the effect of the normal stress on the strength of a frictional interface, let now assume that the stress vector at the matrix-inclusion interface is bounded by a Mohr-Coulomb criterion:

f mc (T ) = T t + α(T n -h) 0. ( 12
)
The positive scalars h and α respectively denote the interface tensile strength and the friction coefficient.

The support function of the Mohr-Coulomb criterion [START_REF] Salencon | An introduction to the yield design theory and its applications to soil mechanics[END_REF] is

π mc ( v ) = +∞ if v n < α v t h v n if v n α v t . ( 13
)
The support function is highly singular. Again, in view of deriving the fictitious state equation of this interface, a direct differentiation is not possible. Hence, following [START_REF] He | Cohesive mohr-coulomb interface effects on the strength criterion of materials with granular-based microstructure[END_REF]; [START_REF] Maghous | Micromechanical approach to the strength properties of frictional geomaterials[END_REF], the support function has to be regularized by a sequence of potentials

ψ j ( v ) = f j (Y ) + h v n with Y = v n -α v t . (14)
The positive scalar j will tend to 0 and f j is a convex function of class C 2 defined on ] -j, +∞[ by

• f j (Y ) = 0 if Y 0, • f j decreasing on ] -j, 0], • lim Y →-j + f j (Y ) = +∞.
Thus, when j tends to 0, the sequence of potentials ψ j tends to the support function π mc of the interface criterion. The state equation associated with the regularized potential is

T = ∂ψ j ∂ v = ∂ψ j ∂ v n n + ∂ψ j ∂ v t t. (15)
Introducing the function F (Y ) = ∂ψ j ∂Y which is negative or null, the state equation is

T = F (Y )(n -αt) + hn. (16)
Hence, for T = hn, necessarily F (Y ) < 0 and thus Y ∈] -j; 0[. Asymptotically, the flow rule v n = α v t is met as j tends to 0. This feature of the described technique is very important in so far as it means that the normality rule in the interfaces is correctly accounted for.

In order to have positive tangential and normal modulii, the state equation is written in an affine form:

T = K( v n , v t ) • v + T p ( v n , v t ). ( 17
)
where

K n = F Y ; K t = - αF v t ; T p = h - αF Y v t n. (18)
Asymptotically, the affine modulii meet

lim j→0 K t K n = lim j→0 - αY v t = 0 and lim j→0 T p K n = -αv t . (19)
3. Linear and non linear homogenization

Linear homogenization

As stated before, the dual approach for limit analysis developed in the previous section leads to a fictitious non linear viscous problem. As far as homogenization is concerned, the modified secant modulii method (see e.g. [START_REF] Suquet | Overall properties of nonlinear composites: a modified secant moduli approach and its link with Ponte-Castaneda's nonlinear variational procedure[END_REF][START_REF] Suquet | Effective behavior of nonlinear composites[END_REF]) provides an efficient technique to deal with non linearity. It relies upon the solution of a linear homogenization problem together with the derivation of appropriate phase average quantities. For this reason, this section is devoted to the resolution of the linear viscous problem and to the derivation of average estimates for the strain rate in the matrix as well as for the velocity jump in interfaces.

Definition of the linear problem

To account for the prestressed formulation of the matrix and interface state equations, the linear problem must include three loading parameters: the macroscopic strain rate D, the matrix uniform prestress p and the interface uniform prestress ̟. The linear problem is thus

divσ = 0 (Ω) σ = C m : d -p1 (Ω m ) T = K • v + ̟n (Γ) σ = C i : d (Ω i ) with C i → ∞ v(z) = D • z (∂Ω) d = grad s v (Ω) (20)
The linear problem ( 20) is solved by superposition of the following loading modes

• (P1) : no prestress, macroscopic strain D at the boundary ∂Ω.

• (P2) : prestress -(p + ̟)1 in the matrix only, no prestress in the interfaces, macroscopic strain 0 at the boundary ∂Ω.

• (P3) : prestress ̟1 in the matrix and ̟n in the interfaces, macroscopic strain 0 at the boundary ∂Ω

The local solution of problem (Pi) for i = 1, 2, 3 is denoted (v i , σ i ), the associated strain rate is d i and the macroscopic stress is Σ i . The solution of (P3) is the uniform couple σ 3 = ̟1 and v 3 = 0. The resolution of (P1) will yield the homogenized stiffness C hom such that Σ 1 = C hom : D.

The macroscopic potential energy density Ψ of the linear problem (20) comprises the terms of elastic energy and the work of the prestresses

|Ω|Ψ(D, p, ̟) = 1 2 Ωm d : C : d d V + 1 2 Γ v : K : v d S -(p + ̟) Ωm tr(d) d V + ̟ Ω tr(d) d V. (21)
Akin to poroelasticity where the change in pore volume fraction is the dual macroscopic state variable associated to pore pressure, let's introduce here an additional macroscopic state variable defined as the dual quantity of the matrix prestress. This dual quantity is the matrix volume change normalized by the initial volume of the rev:

v im = (1 -ρ) tr(d i m ) (i = 1, 2) (22)
and by superposition:

v m = (1 -ρ) tr(d m ) (23)
Using the superposition rules v = v 1 + v 2 and d = d 1 + d 2 , the density of potential energy Ψ(D, p, ̟) reads

|Ω|Ψ(D, p, ̟) = 1 2 Ωm d 1 : C : d 1 d V + 1 2 Γ v 1 • K • v 1 d S + 1 2 Ωm d 2 : (σ 2 + (p + ̟)1) d V + 1 2 Γ v 2 • (σ 2 • n) d S + Ωm d 2 : σ 1 d V + Γ v 2 • (σ 1 • n) d S + |Ω| (-(p + ̟)v m + ̟ tr D) .
According to the energy definition of C hom , the sum of the two first integrals of the rhs is equal to

|Ω| 2 D : C hom : D.
The second and the third line can be simplified by application of Hill's lemma extended to interfaces [START_REF] Maalej | Micromechanical approach to the failure criterion of granular media[END_REF] to the couples of kinematically (resp. statically) admissible velocity and stress fields (v 2 , σ 2 ) and (v 2 , σ 1 ). Eventually, one obtains:

Ψ(D, p, ̟) = 1 2 D : C hom : D -(p + ̟) v 1m + 1 2 v 2m + ̟ tr D. (24)

Quadratic effective estimates

In view of the application of the modified secant modulii method to solve the non linear problem, volume and deviatoric estimates of the strain rate in the matrix are needed. In the modified secant modulii method, these effective values are defined as the quadratic averages over the matrix phase

d eff v = d 2 v m ; d eff d = d 2 d m .
Following [START_REF] Dormieux | Micromechanical approach to the behavior of poroelastic materials[END_REF]; [START_REF] He | Cohesive mohr-coulomb interface effects on the strength criterion of materials with granular-based microstructure[END_REF], adapted from [START_REF] Kreher | Residual stresses and stored elastic energy of composites and polycrystals[END_REF], the estimates are

1 2 (1 -ρ) d eff v 2 = ∂Ψ(D, p, ̟) ∂k m , 1 2 (1 -ρ) d eff d 2 = ∂Ψ(D, p, ̟) ∂2µ m . (25)
Similarly, effective jump of velocity at the interface is defined as

v n eff = v n 2 Γ ; v t eff = v t 2 Γ .
The results of [START_REF] Kreher | Residual stresses and stored elastic energy of composites and polycrystals[END_REF] have been extended to interfaces in [START_REF] Maalej | Micromechanical approach to the failure criterion of granular media[END_REF], which yields

1 2 |Γ| |Ω| v n eff 2 = ∂Ψ(D, p, ̟) ∂K n . 1 2 |Γ| |Ω| v t eff 2 = ∂Ψ(D, p, ̟) ∂K t .
(26)

Mori Tanaka estimates

The linear problem (20) will be solved using a generalized Mori-Tanaka scheme in order to take into account the matrix-inclusion morphology. In the case of spherical inclusions, it is well-known that the use of this scheme can be theoretically justified when the spatial distribution of the inclusions is isotropic. This assumption is adopted in the sequel. The radius of the spherical inclusions is r 0 and the volume fraction of the inclusionary phase is ρ. Accordingly, the specific surface of the interfaces is

|Γ| |Ω| = 3ρ r 0 . ( 27
)
Solution of problem (P1). The matrix-inclusion morphology suggests to resort to a Mori-Tanaka homogenization scheme in order to derive an estimate of C hom . However, the morphologically representative pattern is not a homogeneous inclusion but a composite inclusion made up of a rigid sphere surrounded by an elastic interface. The modified Eshelby problem of a sphere surrounded by a weak interface as introduced in Hashin (1991) has thus to be used. Transposing the results of [START_REF] Dormieux | Micromechanical models of the strength of a sandstone[END_REF] derived for a self consistent scheme, the Mori-Tanaka estimate of the homogenized stiffness is C hom = 3k hom J + 2µ hom K where the compressibility and shear modulii are

k hom = 1 3 (3k m + 4ρµ m )r 0 K n + (1 -ρ)12k m µ m (1 -ρ)r 0 K n + 4µ m + 3ρk m , 2µ hom =µ m (((3 + 2ρ)4µ m + (2 + 3ρ)3k m )r 0 K n + 4µ m ((2 + ρ)8µ m + (1 + ρ)9k m ))r 0 K t + ... (3(1 -ρ)(k m + 2µ m )r 0 K n + ((4 -3ρ)4µ m + (3 -2ρ)3k m )2µ m )r 0 K t + ... ...(((3 + ρ)8µ m + (5 + 3ρ)3k m )r 0 K n + (1 -ρ)(8µ m + 9k m )8µ m )2µ m ...(3((5 -2ρ)k m + (2 -ρ)4µ m )r 0 K n + ((2 + 3ρ)4µ m + (3 + 2ρ)3k m )8µ m )µ m . ( 28
)
In the following developments, the limit case K n /K t → ∞ will be of particular interest as seen from ( 11) and ( 19). The homogenized modulii then read

k hom,l = 3k m + 4ρµ m 3(1 -ρ) , 2µ hom,l = µ m (3k m (2 + 3ρ) + 4µ m (2ρ + 3))r 0 K t + 2µ m (8(ρ + 3)µ m + 3(3ρ + 5)k m ) 3 ((1 -ρ)(k m + 2µ m )r 0 K t + µ m ((5 -2ρ)k m + 4(2 -ρ)µ m )) . ( 29
)
In view of the derivation of the potential energy (see ( 24)), it is useful to determine the matrix volume change v 1m as a function of the macroscopic strain rate D. Solving of (P1) using a Mori-Tanaka scheme leads to

v 1m = b m 1 : D (30)
where the coefficient b m is equal to

b m = (1 -ρ) r 0 K n + 4µ m (1 -ρ)r 0 K n + 4µ m + 3ρk m . ( 31
)
Solution of problem (P2). At the microscopic scale, the stress field linearly depends on p + ̟. Hence, there is a second order tensor B 2 such that σ 2 = -B 2 (p+̟). The macroscopic stress is then Σ 2 = -B 2 Ω (p+̟).

In the line of reasoning of the Maxwell-Betti theorem and its application to poroelasticity [START_REF] Dormieux | Micromechanical approach to the behavior of poroelastic materials[END_REF], first note that:

|Ω|d 1 : C : d 2 Ω + |Γ| v 1 • K • v 2 Γ = |Ω|d 2 : C : d 1 Ω + |Γ| v 2 • K • v 1 Γ .
where the symmetry of C and K was used. Careful application of Hill's lemma (in its extended form to interfaces [START_REF] Maalej | Micromechanical approach to the failure criterion of granular media[END_REF]) to the crossed couples (d 1 , σ 2 ) and (d 2 , σ 1 ), with account for the prestress p + ̟ then yields:

d 1 Ω : σ 2 Ω + (p + ̟)1 : (1 -ρ)d 1 m = d 2 Ω : σ 1 Ω .
Next, recalling the boundary conditions in (P1) and (P2) as well as (30), one obtains:

-(p + ̟)B 2 Ω : D + (p + ̟)b m 1 : D = 0.
The macroscopic stress of (P2) is finally derived as

Σ 2 = -(p + ̟)b m 1.
In turn, let us determine the normalized matrix volume change v 2m which is involved in (24). To do so, the previous expression of the macroscopic stress in (P2) is confronted to the stress average rule (see Eqs. ( 35) to (38) for details)

Σ 2 = -(p + ̟)b m 1 = (1 -ρ)(C m : d 2 Ωm -(p + ̟)1) + |Γ| |Ω| r 0 n ⊗ K • v 2 Γ ,
Applying the trace operator, the above equation simplifies into

-3(p + ̟)b m = 3k m v 2m -(1 -ρ)3(p + ̟) + |Γ| |Ω| r 0 K n v 2n Γ .
Since the macroscopic strain rate is null in problem (P2) and recalling that the inclusions are rigid (no strain), the strain average rule reduces to:

0 = v 2m + |Γ| |Ω| v 2n Γ .
Finally, the matrix volume change depends on the prestress -(p + ̟) as

v 2m = - 3(b m -1 + ρ) 3k m -r 0 K n (p + ̟) = 1 N m (p + ̟),
where the modulus N m takes the form

N m = (1 -ρ)r 0 K n + 4µ m + 3ρk m 3ρ(1 -ρ) . ( 32
)
Complete linear problem. Adapting from [START_REF] Dormieux | Microporomechanics[END_REF], the two macroscopic state equations of the complete linear problem are obtained by superposition of the three loading modes

Σ =C hom : D -(p + ̟)b m 1 + ̟1, v m = b m 1 : D + p + ̟ N m , (33)
where C hom is defined by ( 28), b m by ( 31) and N m by (32).

The potential energy of the complete linear viscous problem is finally obtained by injection of ( 33) in ( 24)

Ψ(D, p, ̟) = 1 2 D : C hom : D -(p + ̟)b m tr D - (p + ̟) 2 2N m + ̟ tr D. ( 34
)
For forthcoming use, we also provide an alternative micro-to-macro expression of the macroscopic mean stress Σ m . The starting point is the stress average rule Σ m = σ m which in turn reads:

Σ m = (1 -ρ)σ m m + ρσ m i ( 35
)
We now evaluate the average σ m i from the average T n Γ of the normal component. Indeed, considering a given particle P ⊂ Ω i , one obtains:

1 |P| P σ dV = 1 |P| z ⊗ T dS. (36)
We now consider the case of a spherical particle with radius r 0 . In this case, we note that z is r 0 n, where n is the unit normal vector to ∂P. From the above equation, it is readily seen that:

σ m P = T n ∂P , (37)
which eventually yields

σ m i = T n Γ . ( 38
)
Introducing this result into (35) and recalling (20), we eventually obtain the following expression for Σ m :

Σ m = k m v m -(1 -ρ)p + ρ(K n v n Γ + ̟). (39)

Non linear homogenization

We now consider the non linear problem (2). It is formally identical to the linear one that has just been dealt with, up to the fact that the non linear modulii (5) are substituted for C m and (11) or ( 18) for K in (20).

To account for the non linearity of the matrix stiffness tensor in a simplified manner, the idea is to assume that the non linear stiffness tensors can be viewed as functions of appropriate averages of the strain rate, referred to as "effective strain" (rate). In the modified secant method [START_REF] Suquet | Overall properties of nonlinear composites: a modified secant moduli approach and its link with Ponte-Castaneda's nonlinear variational procedure[END_REF][START_REF] Suquet | Effective behavior of nonlinear composites[END_REF], these effective estimates of the volume and deviatoric strain rates are quadratic averages over the matrix phase of these quantities. Namely, the effective values of the secant modulii (5) in the matrix phase are estimated using (25) as

k eff m ≈ k m (d eff v , d eff d ) ; µ eff m ≈ µ m (d eff v , d eff d ). ( 40
)
Similarly, the non linearity of the interface stiffness is taken into account through the estimates of the jump of velocity components (26) across the interface

K eff n ≈ K n ( v n eff , v t eff ) ; K eff t ≈ K t ( v n eff , v t eff ). ( 41
)

Macroscopic criteria

Perfect interfaces

To start with, perfect interfaces are considered. In this case, no jump of velocity is allowed at the interface, while the tangential and normal stiffnesses of the interface are infinite. Thus, problem (2) is solved based on the results of the linear problem (20) in the limit case where K n and K t are infinite with no prestress in the interfaces (̟ = 0).

As far as the underlying linear problem is concerned, the previous decomposition in the three subproblems becomes trivial. For (P1), the linear homogenized modulii (28) reduce to the classical Mori-Tanaka estimates for a matrix composite with perfectly bounded and infinitely stiff inclusions

k hom = 3k m + 4ρµ m 3(1 -ρ) , 2µ hom = µ m 3k m (2 + 3ρ) + 4µ m (2ρ + 3) 3(1 -ρ)(k m + 2µ m ) . (42)
For (P2), it is readily seen that the solution is (v 2 = 0, σ 2 = -p1) while the solution to (P3) is (v 3 = 0, σ 3 = 0) (since ̟ = 0). Computations are carried out by substituting the homogenized modulii (42) in the potential energy with p = c (see ( 5)) and ̟ = 0. Then estimates of the strain rate are derived using (25) and the non linear equations ( 40) are solved for d eff v and d eff d . The solution is reported in µ eff m :

2µ eff m = (1 -ρ)b D d 6(a + b) 2(a + b)(2ρb + 3a)β 2 + bq , (43)
where the constant q is defined at (48) and β is the strain triaxiality ratio defined as

β = D v /D d . (44)
Then, recalling (5), one obtains k eff m = 2µ eff m a/b. Using these expressions of the effective matrix modulii in (42), the macroscopic stress Σ is eventually obtained from the first equation of (33) which provides parametric equations for Σ m and Σ d as functions of β. Eliminating β yields the equation of the boundary ∂G hom , that is, the macroscopic criterion. As expected, the latter is isotropic. ∂G hom is an ellipse in the (Σ m , Σ d ) plane. The domain G hom is defined as:

(Σ m + c) 2 A + Σ 2 d B -1 0, where A = a 1 + ρ 2b 3a ; B = b 1 + ρ 9a + 4b 6(a + b) . (45)
It is readily seen that G hom contains the matrix criterion (see ( 3)). This result had already been obtained by [START_REF] Shen | A closed-form three scale model for ductile rocks with a plastically compressible porous matrix[END_REF] using the same homogenization techniques.

Tresca interfaces

Next, the case where the interfaces are governed by a Tresca criterion (6) is considered. First, the macroscopic criterion is derived under the assumption of interface activation. Then the conditions for interface activation are derived.

Derivation of the macroscopic criterion when the interfaces are activated

As one can see from ( 5) and ( 10), the considered underlying linear problem has a prestress p = c in the matrix, no prestress in the interface (̟ = 0), and an infinite normal modulus K n in the interface.

The effective value of the tangential jump of velocity v t eff across the interface is expressed as a function of the matrix effective modulii by use of ( 26) and ( 34) where C hom is characterized by ( 29) and K t by (11). Doing so, a second order polynomial equation in v t eff appears. By definition, v t eff is positive, and the equation has only one or no positive root depending on the tangential strength k of the interfaces. Two cases have to be studied 1. Inactive interfaces: when no real positive solution of the equation is available, the interfaces are not active, meaning that the yield condition of the interface criterion is not reached. In this case, v t eff = 0 and K t → ∞ so that the macroscopic criterion for perfect interfaces ( 45) is retrieved. 2. Active interfaces: the equation has a positive root v t eff > 0, which means that the interfaces are active.

The mathematical condition for interface activation will be provided later. The tangential velocity jump v t eff , as well as K eff t through ( 11) are now functions of the effective volume and deviatoric strain rates in the matrix. In turn, these matrix effective strain rates d eff v and d eff d are expressed by use of ( 25), ( 34) and ( 29) in which k m is substituted for k m = 2µ m a/b using ( 5). The expressions of d eff v and d eff d are then injected in µ m , resulting in a second order polynomial equation in µ eff m which has a single positive root. Explicit effective values of k eff m and K eff t are deduced, and substituted in µ hom and k hom to find a parametric equation of the macroscopic stress states lying on the boundary of the macroscopic strength criterion using (33). The macroscopic stress state obtained from (33) lies on a portion of ellipse of equation

(Σ m + c) 2 A + (Σ d -Σ d0 ) 2 B 1. ( 46
)
If the interfaces are not active, the ellipses parameters are defined by ( 45) and Σ d0 = 0. Otherwise, when interfaces mechanisms are active, the parameters of the ellipse are

A = a 1 + 2ρb 3a 1 - k k e 2 B = (1 -ρ)rb 2 18ρ(b + a)k 2 e 1 - k k e 2 Σ d0 = √ 5(2b + 3a) 6(b + a) b k 2 e k ( 47
)
where the following useful notations are introduced:

t = (4 -2ρ)b + (5 -2ρ)a, q = (6 + 4ρ)b + (6 + 9ρ)a, r = (4ρ + 12)b + (9ρ + 15)a, ( 48
)
and k e is a specific value of the tangential strength defined in (49).

Conditions for interface activation

At this stage, the condition for the activation of the interfaces has to be explicited. This condition is shown to depend on the macroscopic strain triaxiality ratio β = D v /D d . First, the existence of roots to the equation on v t eff imposes

k 2 k 2 e = tb 6ρ(b + a) . ( 49
)
Further, one of the roots is positive if

k 2 k 2 c (β) = 5 6 (b(2b + 3a)) 2 2(b + a) 2 (3a + 2bρ)β 2 + (b + a)bq , ( 50
)
where the inequality k 2 c (β) < k 2 e is met irrespective of the value of β. Also note that k c (β) is an even decreasing function of β, thus reaching its maximum for β = 0 (limit case of a purely deviatoric macroscopic strain rate).

To summarize, two cases have to be considered

• If k > k c (0)
, the interfaces are never active because the interfaces are too strong (Fig. 2a,2b): the macroscopic criterion is equal to the criterion with perfect interfaces (45).

• If k k c (0), a critical value of the strain triaxiality ratio β = D v /D d appears

β 2 c = 5 12(3a + 2bρ) (2b + 3a)b (b + a)k 2 1 - k k c (0) 2 . ( 51
)
The discussion concerning the status of the interfaces is controlled by the strain triaxiality ratio (44):

-if β 2 β 2 c , the interfaces are not active. -if β 2 < β 2
c , the interfaces are active. The effective matrix shear modulus is

2µ eff m = (1 -ρ)b D d 3t (t(3a + 2bρ)β 2 + (1 -ρ)br) 1 - k k e 2 (52)
and the tangential velocity jump is

v t eff = √ 5 3 (2b + 3a)r 0 D d t 1 - k k c (β) . ( 53
)
Hence, provided that the tangential strength of the interface is sufficiently low (k k c (0)), the boundary ∂G hom is a combination of two ellipses : for high triaxiality ratios (β 2 β 2 c ), the criterion for perfect interfaces is retrieved while for low triaxiality ratios (β 2 β 2 c ), the criterion is weakened by the activation of the interfaces and branches into a smaller ellipse (Fig. 2c,2d), which center lies on the vertical axis Σ m = -c. Note that the two ellipses are tangent at their intersections for β 2 = β 2 c . At fixed matrix strength properties and inclusions concentration, reduction of the interface tangential strength k results in a lower deviatoric strength, whereas the hydrostatic strengths, either in compression or traction, are unchanged.

Mohr-Coulomb interfaces

Finally, the case where the interfaces strength is governed by a Mohr-Coulomb criterion ( 12) is considered. The underlying linear problem makes full use of the results derived in section 3.1, with prestress in the matrix (p = c) and in the interfaces (̟ = T p ).

Derivation of the macroscopic criterion when the interfaces are activated

Proceeding as for the Tresca interfaces case, the effective tangential jump of velocity v t eff across the interface is expressed as a function of the matrix effective modulii and the interface function F by use of ( 26) and (34) where C hom is characterized by (29) and K t by (18). Recalling from (5), k m = 2µ m a/b, and the only conditionally positive root is The normal velocity jump can be determined either from the normality rule v n eff = α v t eff (see section 2.4) or consistently with the quadratic estimate (26). Moreover, the direct estimate v n Γ proves to be equal to the quadratic estimate v n eff . The interface parameter F (Y ) (see ( 16)) is determined from the comparison of the macroscopic expression of the mean macroscopic stress derived from (33):

v t eff = 6(1 -ρ)(a + b)αF + 2 √ 5µ eff (2b + 3a)D d 6tµ eff r 0 . ( 54
) 0 0.5 1 1.5 2 -4 -3 -2 -1 0 1 2 Σ d Σ m (a) k > ke 0 0.5 1 1.5 2 -4 -3 -2 -1 0 1 2 Σ d Σ m (b) kc(0) < k = 1.5 < ke 0 0.5 1 1.5 2 -4 -3 -2 -1 0 1 2 Σ d Σ m (c) k = 0.8 < kc(0) 0 0.5 1 1.5 2 -4 -3 -2 -1 0 1 2 Σ d Σ m (d) k = 0.3 < kc(0)
Σ m = k hom D v -cb m + (1 -b m )T p
and the direct average of the microscopic mean stress derived from (39):

Σ m = k m v m -(1 -ρ)c + ρ(K n v n eff + T p ).
where we have used the fact that v n eff = v n Γ . The matrix volumic strain rate v m is expressed thanks to the second macroscopic state equation ( 33). Then, according to the interface state equation ( 16), K n v n eff + T p = h + F . Furthermore, in the limit case K n , T p → ∞, the interface fictitious behavior is such that

b m → 1 ; (1 -b m )T p → - ρ(3k m + 4µ m ) (1 -ρ)r 0 v n eff ; T p N m → - 3ρ r 0 v n eff .
From the comparison of the two above expressions of the macroscopic mean stress, we get:

F = 3k eff + 4µ eff 3(1 -ρ) D v - 3ρk eff + 4µ eff r 0 (1 -ρ) v n eff -c -h. (55)
The last step proceeds as follows: An equation on v t eff is to be determined from substitution of ( 55) into (54) with the normality rule v n eff = α v t eff and the identity k eff = 2µ eff a/b. Akin to the Tresca interfaces, two cases have to be studied while solving this equation to get v t eff as a function of µ eff , D v and D d :

1. Inactive interfaces: when this equation has no positive solution, this physically means that the interfaces are not active, so that v t eff = 0, thus K t → ∞ and the macroscopic criterion for perfect interfaces (45) is retrieved. 2. Active interfaces: the equation has a positive root v t eff > 0, which means the interfaces are active.

In the case where interfaces are active, the effective matrix modulii have to be determined assuming a strictly positive value of v t eff exists. For the sake of clarity, the following additional notations are introduced:

u = 2(a + b)(3ρa + 2b)α 2 + bt, w = 2(a + b)qα 2 + 5b(3a + 2bρ). ( 56
)
Assuming v t eff > 0, the effective value v n eff is deduced from the normality rule and K eff t from (18). Then, the matrix strain rate effective values are expressed using (25). Since the interfaces normality rule is asymptotically met, the expression is taken in the limiting case (19) as for the calculation of F . The effective strain rate estimates are cast in (5) to get a second order equation in µ eff m , which has a positive root in the case where (h + c) 2 h 2 4 with

h 2 4 = u 6ρ(a + b)α 2 . ( 57
)
The effective matrix shear modulus µ eff m finally is

2µ eff m = (1 -ρ) 3bu A ′ D 2 v + B ′ D 2 d -2C ′ D v D d 1 - h + c h 4 2 , ( 58
)
where the following notations have been introduced

A ′ = 12a(a + b)(1 -ρ) 2 α 2 + (3a + 2bρ)t, B ′ = (3ρa + 2b)qα 2 + b(1 -ρ)r, C ′ = √ 5ρ(2b + 3a) 2 α. (59)
The tangential velocity jump, whose sign still has to be discussed, is

v t eff = r 0 3u (2b + 3a)( √ 5bD d + 2α(a + b)D v ) - 3(1 -ρ)b(a + b)α(h + c) µ eff m . ( 60
)
The effective modulus k eff m is derived from µ eff m (58) by ( 5) while K eff t (see ( 18)) is obtained from v t eff (60) and F (55). These quantities are finally introduced into (29) in order to get µ hom and k hom . The macroscopic stress states (Σ m , Σ d ) at the boundary of the macroscopic strength criterion are retrieved from (33). The associated support function for the active interface part of the strength criterion is recognized to be

π hom act (D) = 3b u 1 - h + c h 4 2 (A ′ D 2 v + B ′ D 2 d -2C ′ D v D d ) + Σ m0 D v + Σ d0 D d , (61)
with

Σ m0 = (h + c)α(2b + 3a)ρ u 2(a + b)α -c ; Σ d0 = (h + c)α(2b + 3a)ρ u √ 5b. (62)
The active interface part of the macroscopic strength criterion f hom act , of which π hom act (61) is the support function, has for boundary a slant ellipse in the (Σ m , Σ p ) plane, centered at (Σ m0 , Σ p0 ) and characterized by the equation

f hom act (Σ) = (Σ m -Σ m0 ) 2 A + (Σ d -Σ d0 ) 2 B + 2(Σ m -Σ m0 )(Σ d -Σ d0 ) C -1, (63)
with the following ellipse parameters defined using ( 48), ( 59) and (57) as

A = s B ′ ; B = s A ′ ; C = s C ′ , with s = (1 -ρ)b 3 6a(1 -ρ)qα 2 + (3a + 2bρ)r 1 - h + c h 4 2 .
(64)

Conditions for interface activation

Let us now fully explicit the condition for interfaces activation: To do so, the sign of v t eff as a function of h + c and β = D v /D d has now to be studied and the cases of activation of the interfaces discussed. The discussion is a bit technical and requires to introduce 4 parameters denoted by h i (i = 1, . . . , 4). h 2 4 was defined in (57) while the others parameters are:

h 2 1 = (2b + 3a) 2 3(3a + 2bρ) , h 2 2 = 5b(2b + 3a) 2 6(a + b)qα 2 , h 2 3 = w(2b + 3a) 2 6(b + a)(3a + 2ρb)qα 2 . ( 65
)
It can be shown that:

h 1 < h 3 < h 4 ; h 2 < h 3 .
(66)

The sign of h 1 -h 2 depends on the value of α as compared to

α 2 c = 5b(3a + 2ρb) 2(a + b)q . ( 67
) If α > α c , then h 2 < h 1 , else h 2 h 1 .
1. Activation of the interfaces is not possible for (h + c) 2 > h 2 4 defined in (57), that is if their tensile strength is too high (Fig. 3a). The criterion for perfect interfaces is retrieved. 2. Activation is possible only under the additional condition (h + c) 2 h 2 3 which is a more stringent condition than the first one since h 3 < h 4 . When h 2 3 < (h + c) 2 < h 2 4 , a positive root µ eff exists but has no physical relevancy since the equation for v t eff has no positive root. More precisely, the macroscopic criterion based on (63) mathematically exists: it is strictly included in the ellipse defined by (45) (non active interfaces) and does not intersect with it. Nevertheless, it is not part of the boundary of the strength criterion (Fig. 3b) since v t eff < 0. 3. In the intermediate case h 1 < h + c < h 3 , interfaces are active if the strain triaxiality ratio (44)

β ∈ [β c-; β c+ ] where

β c± = √ 5b(2b + 3a) ± (h + c) 3bw 1 -((h + c)/h 3 ) 2 2(a + b)(2b + 3a)α ((h + c)/h 1 ) 2 -1 . (68)
The intersection of the two parts of the macroscopic strength criterion occurs in the two points corresponding to β c± (Fig. 3c).

The sign of the critical value β c-is physically meaningful: When β c-> 0, the corresponding point on the macroscopic criterion lies on the right of the top of the ellipse. In this case, the maximum of the deviatoric stress is reached in a state where the interfaces are not activated. Interface activation is restricted to high values of the macroscopic mean stress. In contrast, if β c-< 0, the corresponding point on the macroscopic criterion lies on the left of the top of the ellipse which implies that the maximum of the deviatoric stress is reached in a state where the interfaces are indeed activated. This is consistent with the fact that β c-is an increasing function of h. Accordingly, the sign of β c-is determined as follows: β c-is negative when (h + c) 2 < h 2 2 and positive otherwise. Note that β c-is not singular for h + c = h 1 since it simplifies into

β c-= 5b 2 (2b + 3a) ((h + c)/h 2 ) 2 -1 2(a + b)α √ 5b(2b + 3a) + (h + c) 3bw 1 -((h + c)/h 3 ) 2 . ( 69
) 4. If h + c < h 1 , β c+ → ∞ and interfaces are always active if β > β c-.
The intersection of the two parts of the macroscopic strength criterion occurs at a single point (Fig. 3d). In this case, interfaces are active for a purely isotropic traction loading mode (β → ∞).

Thus, when the interfaces tensile strength is low enough (h + c < h 3 ), the criterion (45) corresponding to perfect interfaces is truncated and symmetry along the plane Σ m = -c is broken. The ellipse for the active interfaces part of the strength criterion boundary (64) is tangent in β c± to the ellipse for non active interfaces (45), so the criterion has no angular point at the transition. At a given matrix strength and inclusion concentration, the reinforcement effect is seen to decrease with reduction of the interface tensile strength h or friction angle α (Fig. 3). Finally, it is worth noting that the Tresca interfaces case is retrieved from the Mohr-Coulomb case in the limit α → 0 and h → ∞ with h = k/α. For the three cases of interface criteria under study -perfect, Tresca and Mohr-Coulomb -the homogenized strength criteria estimated by the modified secant modulii method appear to depend only on the first and second invariants of the macroscopic stress tensor. However, Tresca and Mohr-Coulomb criteria can generate a dependence to the third invariant of the macroscopic stress tensor. It would be interesting in future works to assess this effect by resorting to a numerical strength homogenization method.

0 0.5 1 1.5 2 -4 -3 -2 -1 0 1 2 Σ d Σ m (a) h + c > h 4 0 0.5 1 1.5 2 -4 -3 -2 -1 0 1 2 Σ d Σ m (b) h 3 < h + c = 5.5 < h 4 0 0.5 1 1.5 2 -4 -3 -2 -1 0 1 2 Σ d Σ m (c) h 1 < h + c = 3.5 < h 3 0 0.5 1 1.5 2 -4 -3 -2 -1 0 1 2 Σ d Σ m (d) h + c = 1.5 < h 1

Conclusion

In this paper, the effective strength of a matrix-inclusion composite has been considered. Specifically, the reinforcement of the matrix strength by rigid inclusions has been studied in the case where the matrix strength domain is delimited by an ellipse in the mean-deviatoric stress plane. The main focus of this work was the matrix-inclusion interface effects.

Three types of interfaces have been considered, characterized by a strength criterion on the stress vector acting upon the interface : perfect interfaces, Tresca interfaces for which the tangential stress vector is uniformly bounded and Mohr-Coulomb interfaces which model cohesive and frictional contact.

The effective strength has been derived within the continuum micromechanics framework. The boundary of the macroscopic strength criterion has been determined by identification of macroscopic stress states corresponding to a fictitious non linear viscous problem defined on a representative elementary volume. The non linear state equations for each constituent have been derived from the different kinds of strength criterion considered, which required a regularization of the interface criteria. The latter involve prestress terms which depend on the interface type at stake. This non linear problem has been solved by an implementation of the so-called modified secant method extended to the case of prestressed phases and interfaces. Estimates for the underlying prestressed linear problem have been provided using a Mori-Tanaka scheme based on a modified Eshelby problem in order to account for interface effects.

Closed form expressions of the effective strength have been proposed for the three interfaces types considered. First, in the case of perfect interfaces, a reinforcement is always observed. The effective criterion is an ellipse in the mean-deviatoric stress plane, concentric to the matrix criterion one, whose size varies linearly with inclusions concentration.

Second, in the case of Tresca interfaces, two situations have been observed. If the tangential strength of the interface is above a threshold which depends on microstructural properties, the interface mechanisms are not active, and the effective criterion is identical to the one for perfect interfaces. In turn, provided that the tangential strength of the Tresca interfaces is low enough, the effective criterion results from the truncation of the one for perfect interfaces by a parallel but smaller ellipse at deviatoric loads, for which interfacial mechanisms are active. The transition between the two parts of the effective criterion is controlled by the strain triaxiality ratio. This hybrid criterion is regular, without angles at the transition points between the two parts and symmetric.

Third, in the case of Mohr-Coulomb interfaces as compared to Tresca interfaces, the part of the effective criterion corresponding to the activation of the interfaces is located on an ellipse with slant axes. The local non symmetry of the Mohr-Coulomb criterion (which accounts for the sign of the normal stress) is transmitted to the macroscopic criterion: The latter has lost its symmetry with respect to the center of the matrix criterion. As expected, interface activation first occurs at higher values of the macroscopic mean stress.

The homogenization strategy adopted in this work is based on the modified secant modulii method. Hence, the homogenized criteria which have been derived are only estimates of the actual macroscopic criteria. To assess the herein presented criteria, numerical simulations should be carried out in future research on representative patterns such as a composite sphere made up of a rigid core surrounded by a matrix, with Tresca or Mohr-Coulomb interfaces between the two phases. We also expect that such simulations could allow to quantify the effect of the third invariant on the homogenized strength.

Figure 2 :

 2 Figure 2: Macroscopic strength criterion for Tresca interfaces (thick blue-red line). The blue line corresponds to inactive interfaces and the red line to active interfaces, for different values of the interfaces tangential strength k. The fine solid black line is the matrix criterion for a = 5, b = 2 and c = 1. The inclusion volume fraction is ρ = 0.5. The dotted blue-red line are the remaining part of the ellipsis which are not part of the criterion boundary. The critical tangential strength values are k e = 1.57 and k c (0) = 1.12.

Figure 3 :

 3 Figure 3: Macroscopic strength criterion for Mohr-Coulomb interfaces (thick blue-red line). The blue line corresponds to inactive interfaces and the red line to active interfaces, for different values of the interfaces tensile strength h. The fine solid black line is the matrix criterion for a = 5, b = 2 and c = 1. The inclusion volume fraction is ρ = 0.5 and the friction angle α = 0.3. The dotted blue-red line are the remaining part of the ellipsis which do not belong to the criterion boundary. The critical interface strength values are h 1 = 2.66, h 2 = 3.73, h 3 = 4.58 and h 4 = 5.93.


Manuel Kauers
email: mkauers@risc.jku.at.

Rika Yatchak
email: ryatchak@risc.jku.at.

Walks in the Quarter Plane with Multiple Steps

Keywords: Lattice Walks, D-finiteness, Computer Algebra

We extend the classification of nearest neighbour walks in the quarter plane to models in which multiplicities are attached to each direction in the step set. Our study leads to a small number of infinite families that completely characterize all the models whose associated group is D4, D6, or D8. These families cover all the models with multiplicites 0, 1, 2, or 3, which were experimentally found to be D-finite -with three noteworthy exceptions.

Résumé. Nous étendons la classification des marches aux plus proches voisins dans le quart de plan à des modèles dans lesquels une multiplicité est attachée à chaque direction de l'ensemble des pas. Notre étude identifie un petit nombre de familles infinies qui caractérisent complétement tous les modèles dont le groupe est D4, D6 ou D8. Ces familles contiennent tous les modèles à multiplicités 0, 1, 2 ou 3 dont il a été prouvé expérimentalement qu'ils étaient D-finis -avec trois exceptions notables.

Introduction

We consider quadrant walk models where step sets may contain several distinguishable steps pointing into the same direction. For example, the step sets {←, ↓, } and {←, ← , ↓, } are considered different, as the latter contains two different ways of going to the left. The objects being counted are then walks in the quarter plane starting at the origin, consisting of n consecutive steps taken from the step set in such a way that the walk never leaves the first quadrant, ending at a point (i, j) ∈ N 2 (with the convention 0 ∈ N), and one of k different colors is attached to each step in the walk whose multiplicity in the step set is k. For each model (viz., for each multiset of admissible directions), we want to know whether the corresponding generating function f (x, y, t) = ∞ n=0 i,j f i,j,n x i y j t n which counts the number f i,j,n of walks of length n ending at (i, j) is D-finite. As usual, a power series in t is D-finite if it satisfies an ordinary linear differential equation with polynomial coefficients. If we let a u,v denote the multiplicity of the direction (u, v) ∈ {-1, 0, 1} 2 \ {(0, 0)}, then the number f i,j,n of walks of length n ending at (i, j) is uniquely determined by the recurrence equation

f i,j,n+1 = u,v a u,v f i-u,j-v,n (n ∈ N, i, j ∈ N)
together with the initial values f 0,0,0 = 1, f i,j,0 = 0 for (i, j) = (0, 0), and the boundary conditions f -1,j,n = f i,-1,n = 0 for all i, j, n. Equivalently, we can say that the generating function f (x, y, t) = ∞ n=0 i,j f i,j,n x i y j t n ∈ Q[x, y][ [t]] satisfies the functional equation

1 -t u,v
a u,v x u y v f (x, y, t)

= 1 -t y u a u,-1 x u f (x, 0, t) -t x v a -1,v y v f (0, y, t) + ta -1,-1 xy f (0, 0, t).

(1)

For the models where all multiplicities a u,v are in {0, 1}, a complete classification is available: among the 2 8 = 256 different models, Bousquet-Mélou and [START_REF] Bousquet-Mélou | Walks with small steps in the quarter plane[END_REF] identified 79 nontrivial cases. For 22 of them they prove that the generating function is D-finite using certain groups G associated to each of the models. For a 23rd model, the notorious Gessel model {←, →, , }, their techniques do not apply but a proof by a different method based on computer algebra was found by [START_REF] Bostan | The complete generating function for Gessel walks is algebraic[END_REF]. A computer-free proof was later found by [START_REF] Bostan | The complete generating function for Gessel walks is algebraic[END_REF]. The remaining 56 models are not D-finite: [START_REF] Mishna | Two non-holonomic lattice walks in the quarter plane[END_REF] and [START_REF] Melczer | Singularity analysis via the iterated kernel method[END_REF] showed that the generating functions of five of these models have infinitely many singularities and therefore are not D-finite. For the remaining models, Bostan et al. (2014b) proved that the counting sequences f 0,0,n for walks returning to the origin have asymptotic behaviour for n → ∞ that D-finite functions cannot possibly have.

The need for a classification of quarter plane models with multiplicities arose in the classification project for octant models in 3D (Bostan et al., 2014a), as it turns out that some models in 3D can be reduced by projection to 2D models with multiplicities. For example, it is easy to see that the generating function for the octant model with step set {

-1 0 0 , -1 0 1 , 0 -1 0 , 1 1 0
} is D-finite if and only if the quadrant model with step set {←, ← , ↓, } is. Bostan et al. (2014a) have classified only the 527 models that they needed for their study, and point out that the classification problem for models with multiplicities is of interest in its own right. If the multiplicities are real numbers in the interval [0, 1] whose sum is 1, they can be interpreted as probabilities of random walks. If the multiplicities a i,j are left symbolic, then f i,j,n is a polynomial in which the coefficient of u,v a

eu,v u,v • • • a e1,1
1,1 indicates how many walks there are in which the direction (u, v) is used e u,v many times, for all u, v.

For the present paper we carried out a systematic search over all the 4 8 = 65536 models where each of the eight directions may have any of the four multiplicities 0, 1, 2, 3. Of these, 30307 are nontrivial and essentially different. Of these nontrivial models, 1457 turn out to be D-finite. 79 of the D-finite models are even algebraic. Going one step further, we have identified families of D-finite models in which some or all of the "multiplicities" are arbitrary complex numbers. Rather than asking for a fixed model what the corresponding group is, we ask for a fixed group what all the models leading to this group are. In this way we obtain a small number of families that completely characterize all the models which lead to groups with at most eight elements. This characterization covers 1454 of the 1457 D-finite cases we discovered for multiplicities in {0, 1, 2, 3}. The remaining three models have a group of order 10, which was too hard for us to analyze in full generality. In view of the fact that all models previously considered had either a finite group of order at most eight or an infinite group, the appearance of these models was a surprise to us. We were less surprised to find, after spending some 6.5 years of computation time, that none of the models with a (probably) infinite group appears to be D-finite based on the inspection of the first 5000 terms.

Models of Interest

Our reasoning largely follows that of [START_REF] Bousquet-Mélou | Walks with small steps in the quarter plane[END_REF]. Their first step is to identify the interesting models. By a model, we understand here a particular choice of multiplicities a u,v ∈ C (not necessarily integers). For each such model, there is a corresponding generating function f (x, y, t) ∈ C[x, y][[t]], and we want to identify the models whose generating functions are D-finite.

A model is uninteresting if a 1,-1 = a 1,0 = a 1,1 = 0 or a -1,1 = a 0,1 = a 1,1 = 0 or a -1,-1 = a -1,0 = a -1,1 = 0 or a -1,-1 = a 0,-1 = a 1,-1 = 0, because in either of these cases the corresponding generating function is algebraic and it is well-understood why (Flajolet and Sedgewick, 2009, Section VII.8). Secondly, if two models can be obtained from one another by reflecting the step set about the diagonal x = y, then the corresponding generating functions can be obtained from one another by exchanging the variables x ↔ y, and therefore either both are D-finite or neither is. Similarly, if one model can be obtained from another by multiplying all multiplicities by a nonzero constant λ, then its generating function can be obtained from the generating function of the other by sending t to λt, and therefore again either both are D-finite or neither is.

Applying all these filters to the 4 8 = 65536 models with possible multiplicities a u,v ∈ {0, 1, 2, 3} leaves us with 30307 nontrivial models (including, for the sake of completeness, the 79 interesting models with a u,v ∈ {0, 1} that have already been completely classified).

3 The Group of the Model For a fixed model, i.e., for a fixed choice of multiplicities a u,v ∈ C, consider the functional equation ( 1). The group associated to the model acts on this equation. Its elements map the variables x and y to certain rational functions in x and y, which are chosen in such a way that all the group elements leave the kernel polynomial K(x, y, t)

:= 1 -t u,v a u,v x u y v
fixed. It is easy to check that the two particular transformations Φ, Ψ : C(x, y) → C(x, y) defined by

Φ : (x, y) → 1 x v a -1,v y v v a 1,v y v , y , Ψ : (x, y) → x, 1 y u a u,-1 x u u a u,1 x u
have this property. It is also easy to check that Φ and Ψ are involutions, i.e., Φ 2 = Ψ 2 = id.

The group G is defined as the group generated by Φ and Ψ under composition. Note that we do not need to worry that one of the denominators u a u,1 x u or v a 1,v y v is identically zero, because this only happens for models that are uninteresting in the sense of the previous section. For the same reason, we may also assume that the numerators u a u,-1 x u and v a -1,v y v , respectively, are nonzero polynomials. In order to argue that the composition of rational functions into the power series of equation ( 1) is algebraically meaningful, recall that the series in question belong to Q[x, y][[t]], so the result of the composition can be naturally interpreted as an element of C(x, y)

[[t]].
The group G is finite if and only if (ΦΨ) n = id for some n ∈ N, and this is the case if and only if G is the dihedral group with 2n elements. The sign sgn(g) of an element g ∈ G is defined to be 1 if g = (ΦΨ) k for some k, and -1 otherwise.

Our definition for G corresponds to the notion of the "group of the walk" in Bousquet-Mélou and [START_REF] Bousquet-Mélou | Walks with small steps in the quarter plane[END_REF]. This notion was in turn inspired by the "group of the random walk", a well-known tool in the probabilistic context that was introduced by Malyshev and further detailed in [START_REF] Fayolle | Random walks in the quarter-plane[END_REF]. The two definitions are not identical, but the group of the random walk can be interpreted as a quotient group of G. Recent results regarding the finiteness of the group of the random walk can be found in [START_REF] Fayolle | Random walks in the quarter-plane with zero drift: an explicit criterion for the finiteness of the associated group. Markov Process[END_REF]; [START_REF] Fayolle | Random walks in the quarter-plane: Advances in explicit criterions for the finiteness of the associated group in the genus 1 case[END_REF].

Models with Group D4

As there is obviously no way to choose a u,v such that ΦΨ = ( 1x r(y), 1 y s(x)) = (x, y) = id, the smallest possible n ∈ N with (ΦΨ) n = id is 2. The group with (ΦΨ) 2 = id is the dihedral group D4 with four elements. In order to determine the models which lead to this group, regard the a u,v as variables and compute (p, q) := Φ(Ψ(x, y)) -Ψ(Φ(x, y)). This is a pair of rational functions in x, y whose coefficients are rational functions in the a u,v over the rational numbers. Write p, q as quotients of polynomials in x, y whose coefficients are polynomials in a u,v with integer coefficients. We want to know the possible choices of a u,v for which p and q become zero. (Note that ΦΨ = ΨΦ ⇐⇒ (ΦΨ) 2 = id because Φ and Ψ are involutions.) In order to find these a u,v , consider the ideal in Q[a -1,-1 , . . . , a 1,1 ] generated by the coefficients of all monomials x i y j in the numerator of p and the coefficients of all monomials x i y j in the numerator of q. This ideal basis consists of 36 homogeneous polynomials of degree 4, which we don't reproduce here because of its length. Using Gröbner basis techniques [START_REF] Becker | Gröbner Bases[END_REF], we can determine the irreducible components of the radical of this ideal. We have used the commands facstd and minAssGTZ of the software package Singular [START_REF] Greuel | A Singular Introduction to Commutative Algebra[END_REF] for this step. It turns out that the two irreducible components are generated by

{ a 0,1 a 1,-1 -a 0,-1 a 1,1 , a -1,1 a 1,-1 -a -1,-1 a 1,1 , a -1,1 a 0,-1 -a -1,-1 a 0,1 }, and { a 1,0 a -1,1 -a -1,0 a 1,1 , a 1,-1 a -1,1 -a -1,-1 a 1,1 , a 1,-1 a -1,0 -a -1,-1 a 1,0 }.
As the latter is obtained from the former by replacing all a u,v by a v,u , it suffices to consider one of the two components, say the first. The equations in this component are equivalent to saying that the vectors (a -1,-1 , a 0,-1 , a 1,-1 ) and (a -1,1 , a 0,1 , a 1,1 ) are linearly dependent. Since the models where one or both of these vectors are zero are uninteresting, the interesting models leading to the group D4 are precisely those for which there exists a constant λ = 0 such that a -1,v = λa 1,v for v = -1, 0, 1. We then have

Φ(x, y) = λ x , y and Ψ(x, y) = x, 1 y λa 1,-1 x -1 + a 0,-1 + a 1,-1 x λa 1,1 x -1 + a 0,1 + a 1,1 x .
At this point, we can proceed analogously to Bousquet-Mélou and Mishna (cf. their Proposition 5): multiplying (1) on both sides by xy/K(x, y, t) and forming the orbit sum gives the general relation

g∈G sgn(g) g(xy f (x, y, t)) = 1 K(x, y, t) g∈G sgn(g) g(xy),
which holds whenever the group is finite. For the special case under consideration, the right hand side evaluates to

1 K(x, y, t) xy - yλ x - x y λa 1,-1 x -1 + a 0,-1 + a 1,-1 x λa 1,1 x -1 + a 0,1 + a 1,1 x + λ xy λa 1,-1 x -1 + a 0,-1 + a 1,-1 x λa 1,1 x -1 + a 0,1 + a 1,1 x = (x 2 -λ)(a 0,1 xy 2 -a 0,-1 x -(λ + x 2 )(a 1,-1 -a 1,1 y 2 )) xy(a 1,1 (λ + x 2 ) + a 0,1 x)K(x, y, t) .
For the left hand side, we have

xy f (x, y, t) - λy x f λ x , y, t - x y s(x) f x, 1 y s(x), t + λ xy s(x) f λ x , 1 y s(x), t ,
where we abbreviate s(x) = λa1,-1x -1 +a0,-1+a1,-1x λa1,1x -1 +a0,1+a1,1x . The identity holds in Q(x, y) [[t]], but it can be seen that all quantities actually belong to

Q(x)[y, y -1 ][[t]].
The last two terms of the equation involve only negative exponents with respect to y, so taking the positive part [y > ] will kill them. The remaining terms happen to belong to Q

[x, x -1 ][[t]],
and since the second term only has negative exponents with respect to x, taking the positive part [x > ] will eliminate it and only leave the first. It follows that

f (x, y, t) = 1 xy [x > ][y > ] (x 2 -λ)(a 0,1 xy 2 -a 0,-1 x -(λ + x 2 )(a 1,-1 -a 1,1 y 2 )) xy(a 1,1 (λ + x 2 ) + a 0,1 x)K(x, y, t) .
Alternatively, we could interpret the elements of Q(x, y) [[t]] as elements of multivariate formal Laurent series field Q ≤ ((x, y, t)) for a term order ≤ with x, y ≤ 1 ≤ t and do the positive part extraction with respect to x and y simultaneously. See Aparicio Monforte and Kauers (2013) for a discussion of formal Laurent series in several variables. In any case, we can summarize the result of this section as the following theorem. A similar characterization already appears as Lemma 4.1.1 in [START_REF] Fayolle | Random walks in the quarter-plane[END_REF].

Theorem 1 The interesting quarter plane models whose group is D4 are precisely those where a -1,v = λa 1,v for v = -1, 0, 1 and some λ = 0. All these models are D-finite.

Family 0

Defining equations:

a 1,0 a -1,1 = a -1,0 a 1,1 , a 1,-1 a -1,1 = a -1,-1 a 1,1 , a -1,1 a -1,0 = a -1,-1 a 1,0
Example:

2 -3 5 -7 13 6 -9 15
5 Models with Group D6

We now determine all the choices for a u,v such that (ΦΨ) 3 = id. As before, we compute (p, q) := Ψ(Φ(Ψ(x, y))) -Φ(Ψ(Φ(x, y))) and consider the ideal generated by the coefficients of the numerators with respect to x, y. The basis consists of 210 homogeneous polynomials of degree 9. The ideal has 34 irreducible components, 18 of which turn out to contain only uninteresting models. Of the remaining 16 components, 6 can be discarded because their solution sets are properly contained in the solution set of others. Of the remaining 10 components, 4 can be discarded because they are reflections of others. This leaves us with the following 6 families:

Family 1a

Defining equations:

a 1,1 = a -1,-1 = 0, a -1,1 a 1,-1 = a -1,0 a 1,0 = a 0,1 a 0,-1
Example:

1/2 1/3 1/5 2 3 5

Family 1b

Defining equations:

a 1,-1 = a -1,1 = 0, a -1,0 a 1,0 = a -1,-1 a 1,1 = a 0,-1 a 0,1 Example: 1/2 1/3 1/5 2 3 5

Family 2a

Defining equations:

a 1,0 = a 1,1 = 0, a 0,-1 a -1,1 = 2a 0,1 a -1,-1 , a 2 0,-1 = 4a 1,-1 a -1,-1 , a 0,-1 a 0,1 = 2a -1,1 a 1,-1
Example:

7 7 5 1 2 1

Family 2b

Defining equations:

a 1,0 = a 1,-1 = 0, a 0,1 a -1,-1 = 2a 0,-1 a -1,1 , a 2 0,1 = 4a 1,1 a -1,1 , a 0,1 a 0,-1 = 2a -1,-1 a 1,1
Example:

7 7 5 1 2 1

Family 3a

Defining equations:

a -1,0 = a -1,-1 = 0, a 0,1 a 1,-1 = 2a 0,-1 a 1,1 , a 2 0,1 = 4a -1,1 a 1,1 , a 0,1 a 0,-1 = 2a 1,-1 a -1,1
Example:

7 7 5 1 2 1

Family 3b

Defining equations:

a -1,0 = a -1,1 = 0, a 0,-1 a 1,1 = 2a 0,1 a 1,-1 , a 2 0,-1 = 4a -1,-1 a 1,-1 , a 0,-1 a 0,1 = 2a 1,1 a -1,-1
Example:

7 7 5 1 2 1
Note that the families on the right can be obtained from those on the left by reflection about the horizontal axis and the families in the third row can be obtained from those in the second row by reversing all arrows. The families in the first row are closed under reversing arrows.

Theorem 2 The interesting quarter plane models whose group is D6 are precisely those that belong to one or more of the families described in the table above. All these models are D-finite.

The remainder of this section is devoted to the D-finiteness claim of this theorem.

Families 1a, 2a, 3a

These families can be handled very much like the family in Section 4 above. Without going into further details, we just report the resulting formulas for the generating functions.

For family 1a, let λ = a -1,1 a 1,-1 = a -1,0 a 1,0 = a 0,1 a 0,-1 . If λ = 0, then all the a u,v are nonzero (except a 1,1 and a -1,-1 of course). In this case, the resulting formula for the generating function is

f (x, y, t) = 1 xy [x > y > ] (a -1,1 -a 0,-1 xy -2 )(a 1,-1 -a -1,0 yx -2 )(λxy -a -1,0 a 0,-1 ) λ 2 K(x, y, t) .
Otherwise, if λ = 0 and a -1,1 = 0, then a -1,0 = 0 and a 0,1 = 0 (otherwise the model is not interesting), but then a 1,0 = 0 and a 0,-1 = 0 (by the defining equations), and then a 1,-1 = 0 (otherwise again the model is not interesting). In this case, the resulting formula for the generating function is

f (x, y, t) = 1 xy [x > y > ] (a 0,1 -a 1,-1 xy -2 )(a 1,-1 -a -1,0 x -2 y)(a 0,1 xy -a -1,0 ) a 2 0,1 a 1,-1 K(x, y, t)
.

Finally, if λ = 0 and a -1,1 = 0, then a 1,-1 = 0 and the only interesting cases have a 1,0 = 0, a 0,-1 = 0, and a -1,0 = a 0,1 = 0. This case is symmetric to the previous case and therefore not interesting.

For family 2a, we may assume that a 1,-1 = 0, because otherwise the model is uninteresting. Then we can also assume a 0,1 = 0, because if a 0,1 = 0, then the last defining equation would imply a -1,1 = 0, which together with a 1,1 = 0 would also render the model not interesting. Under the assumption a 1,-1 = 0, a 0,1 = 0, the generating function can be expressed as

f (x, y, t) = 1 xy [x > y > ] P (x, y) a -1,-1 -a 1,-1 x 2 + a -1,1 y 2 + a -1,0 y 2a 0,1 y 2 -2a 1,-1 x -a 0,-1 4x 2 y 3 a 2 0,1 a 1,-1 K(x, y, t)
where P (x, y) = 2a -1,-1 -2a 0,1 xy 2 + a 0,-1 x + 2a -1,1 y 2 + 2a -1,0 y. For family 3a, models are interesting only when a -1,1 = 0 and a 0,-1 = 0 and (a 1,-1 , a 1,0 , a 1,1 ) = (0, 0, 0). Under these assumptions, we obtain the following expression for the generating function:

f (x, y, t) = 1 xy [x > y > ] Q(x, y) a 0,1 y 2 -2a 0,-1 +2a -1,1 y 2 /x a 1,1 x 2 y+a 1,0 x 2 +a 1,-1 x 2 /y-a -1,1 y a 1,-1 + a 1,0 y + a 1,1 y 2 4a -1,1 a 0,-1 + (2a -1,1 + a 0,1 x Q(x, y))K(x, y, t)
where Q(x, y) = 2a 1,1 xy 2 +2a 1,0 xy +2a 1,-1 x+a 0,1 y 2 -2a 0,-1 . Note that in this case the denominator contains nontrivial factors involving both x and y, so the ad-hoc reasoning used in Section 4, which also works for the families 1a and 2a, does not work here. However, there is no problem if we take the viewpoint of multivariate Laurent series (Aparicio Monforte and Kauers, 2013), because all that is needed for the argument to go through is the property that there exists a term order ≤ so that for all g ∈ G \ {id} and all positive integers i, j the expansion of g(x) i g(y) j ∈ C(x, y) in the multivariate Laurent series field C ≤ ((x, y)) contains no terms x k y where both k and are positive. This turns out to be the case.

Family 1b

For the family 1b there are three cases to distinguish. First, when a -1,-1 = a 1,0 = a 0,1 = 0, then a 1,1 , a -1,0 , a 0,-1 all must be nonzero in order for the model to be interesting. In this case, the generating function is

f (x, y, t) = k 3 a 0,-1 a 1,1 /a 2 -1,0 x, 3 a -1,0 a 1,1 /a 2 0,-1 y, 3 a -1,0 a 0,-1 a 1,1 t ,
where k(x, y, t) is the generating function for classical Kreweras walks (i.e., a 1,1 = a -1,0 = a 0,-1 = 1), which is known to be algebraic [START_REF] Kreweras | Sur une classe de problèmes liés au treillis des partitions d'entiers[END_REF][START_REF] Bousquet-Mélou | Walks in the quarter plane: Kreweras' algebraic model[END_REF]. Secondly, when a 1,1 = a -1,0 = a 0,-1 = 0, algebraicity of the generating function can be established by a similar argument. The third case is when a 1,1 , a -1,-1 , a 1,0 , a -1,0 , a 0,1 , a 0,-1 are all nonzero. In this case it is impossible to express the generating function in terms of the generating function for the corresponding model without multiplicities, known as the double Kreweras model. However, if we let f λ (x, y, t) be the generating function for the family where a

-1,-1 = a -1,0 = a 0,-1 = 1 and a 1,1 = a 1,0 = a 0,1 = λ = 0, then f (x, y, t) = f a0,1a 2 1,0 /(a0,-1a 2 1,1 ) a 1,1 a 0,1 x, a 1,1 a 1,0 y, a -1,-1 a 2 1,1 a 0,1 a 1,0 t
is the generating function of an arbitrary model of family 1b with a 1,1 a -1,-1 = 0. It therefore suffices to show that f λ (x, y, t) is D-finite. We will show that it is in fact algebraic, following the treatment in Section 6.3 of Bousquet-Mélou and [START_REF] Bousquet-Mélou | Walks with small steps in the quarter plane[END_REF] with the added parameter λ. Since the orbit sum is zero, we sum (1) over only half the orbit to obtain a nonzero expression on both sides. This new expression will be more complicated than in the orbit sum case: it will involve the unknown series f λ (x, y, t), f λ (x, 0, t) and f λ (0, 0, t). Using careful coefficient extraction, we will obtain the algebraicity result.

Writing A v = u a u,v x u for v = -1, 0, 1, the half-orbit sum equation reads

xyf λ (x, y, t)- 1 λx f λ 1 λxy , y + 1 λy f λ 1 λxy , x = xy -1 λx + 1 λy -2txA -1 f λ (x, 0, t) + tf λ (0, 0, t) K(x, y, t) .
Next we extract the coefficient of y 0 using Lemma 7 from Bousquet-Mélou and [START_REF] Bousquet-Mélou | Walks with small steps in the quarter plane[END_REF]. We solve K(x, y, t) = 0 for y in terms of x and t: writing

∆(x) := t 2 x -2 -2(t+2λt 2 )x -1 +(1-6λt 2 )-2λt(1+ 2t)x + λ 2 t 2 x 2 for the discriminant of K(x, y, t), the two solutions are Y 0 = 1 -tA 0 -∆(x) /(2tA 0 )
and Y 1 = 1/(λxY 0 ). The coefficient of y n in 1/K(x, y, t) can be expressed in terms of Y 0 , Y 1 , and ∆(x):

[y n ] 1 K(x, y, t) = 1 ∆(x) × Y -n 0 if n ≤ 0 Y -n 1 if n ≥ 0 .
Using these facts, extracting the coefficient of y 0 on both sides of the half-orbit sum equation leads to

- 1 λx d λ 1 λx , t = 1 ∆(x) xY 0 - 1 λx + 1 λY 1 -2txA -1 f λ (x, 0, t) + tf λ (0, 0, t) , (2)
where d λ (x, t) := i,n (f λ ) i,i,n x i t n is the generating function for walks ending on the diagonal.

Now we write

∆(x) = t 2 Z 2 ∆ -(x)∆ + (x)
, where

∆ + (x) = 1 - 2λZ(1 + 2Z + 2λZ 2 + 2λ 2 Z 3 + λ 2 Z 4 ) (1 -λZ 2 ) 2 x + λ 2 Z 2 x 2 , ∆ -(x) = ∆ + 1 x ,
and where

Z ∈ Q[λ][[t]] is defined through Z = t(1+3λZ 2 +4λ(1+λ)Z 3 +3λ 2 Z 4 +λ 3 Z 6 ) (1-λZ 2 ) 2
and Z(0) = 0.

Multiplying (2) by A 1 ∆ -(x) and using the explicit expressions for Y 0 and Y 1 given above, we obtain

∆ -(x) x t - 1 λx A 1 d λ 1 λx , t = ZA 1 t ∆ + (x) x (1 -tA 0 ) tA 1 - 1 xλ -2tA -1 f λ (x, 0, t) + tf λ (0, 0, t) .
From this equation, we extract the coefficient of

x 0 . Using [x 0 ]d λ ( 1 λx , t) = [x 0 ]f λ (x, 0, t) = f λ (0, 0, t), we find f λ (0, 0, t) = Z -4λZ 3 -2λZ 4 -2λ 2 Z 4 -λ 2 Z 5 t(1 -λZ 2 ) 2 .
We can now extract the positive part in x on both sides of the same equation to obtain

f λ (x, 0, t) = x 2 (λZ 2 -1) + 2xZ(λZ + 1) -λZ 3 + Z 2λtx(x + 1) 2 Z(1 -λZ 2 ) ∆ + (x) - Z 2t(1 + x) λtx 3 + 2tx + t -x 2 λtx(x + 1)Z + 2(λ 2 Z 3 + λ(Z + 3)Z 2 -1) (1 -λZ 2 ) 2 + 1 .
Noting that f λ (0, y, t) = f λ (y, 0, t), we conclude from equation (1) that f λ (x, y, t) is algebraic.

Families 2b, 3b

The orbit sum argument also fails for these families. For the models in family 2b the orbit sum is zero, while in family 3b the orbit sum is nonzero but the desired term f (x, y, t) cannot be isolated by taking the positive part because there are group elements g = id for which f (g(x), g(y), t) also contributes terms with positive exponents to the orbit sum. Because of the lack of symmetry, the half orbit sum argument used for family 1b does not seem to apply either. One model from each of these two families were already encountered by Bostan et al. (2014a), and computer proofs have been given there that the generating function for the model belonging to family 2b is algebraic and the model belonging to family 3b is (transcendental) D-finite. The models considered by Bostan et al. (2014a) are

a 1,0 = a 1,-1 = a -1,0 = 0, a -1,1 = 1 2 a 0,1 = a 1,1 = a -1,-1 = a 0,-1 = 1 (case 2b), and its reverse a -1,0 = a -1,1 = a 1,0 = 0, a 1,-1 = 1 2 a 0,-1 = a -1,-1 = a 1,1 = a 0,1 = 1 (case 3b).
We were able to extend these computer proofs to the more general cases where a -1,0 = λ (case 2b), and a 1,0 = λ (case 3b), respectively, are formal parameters. From here, every other model of the respective family can be reached by an appropriate algebraic substitution: if f λ (x, y, t) is the generating function for the model

a 1,0 = a 1,-1 = a -1,0 = 0, a -1,1 = 1 2 a 0,1 = a 1,1 = a -1,-1 = a 0,-1 = 1, a -1,0 = λ, then f (x, y, t) = f a-1,0/ √ a-1,-1a-1,1 a 0,-1 a -1,-1 x, a -1,1 a -1,-1 y, 1 2 a -1,-1 a -1,1 t
is the generating function for an arbitrary model of family 2b, and likewise for family 3b. The computational techniques we used were introduced by [START_REF] Kauers | The quasi-holonomic ansatz and restricted lattice walks[END_REF]; [START_REF] Kauers | Proof of Ira Gessel's lattice path conjecture[END_REF]; [START_REF] Bostan | The complete generating function for Gessel walks is algebraic[END_REF], and they have been described for the cases λ = 0 in the paper of Bostan et al. (2014a). We do not repeat these explanations again but only remark that the additional symbolic parameter λ has made the calculations considerably more expensive. The computations were done using software of [START_REF] Kauers | Guessing handbook[END_REF] and [START_REF] Koutschan | HolonomicFunctions (User's Guide)[END_REF]. The bottleneck was the construction of a certified recurrence for (f λ ) 0,0,n . The (nonminimal) recurrence we found has order 14 and degrees 30, 26 in n, λ, respectively; the certificate for this recurrence is 16 gigabytes long! From this recurrence it can be deduced that f λ (0, 0, t) is the unique formal power series T ∈ Q[λ][[t]] with T (0) = 1 and t 4 T 2 + (2tλ + 1)t 2 T + t(4t + 1) -(3t 2 (λ -4) + 3t + 1)Z + t(6t + 1)(λ + 2)Z 2 = 0,

where Z ∈ Q[λ][[t]
] is the unique formal power series with Z(0) = 0 and t = Z(4Z+1) 1+6Z+12Z 2 +4(2+λ)Z 3 . Using this equation and the functional equation (1) (with f λ in place of f ), we could then prove the correctness of guessed polynomial equations P (x, t, λ, f λ (x, 0, t)) = Q(y, t, λ, f λ (0, y, t)) = 0, which in turn can be used to deduce that f λ (x, 0, t) is the unique formal power series

U ∈ Q[x, λ][[t]] with U (0) = 1 and (x + 1) 2 t 4 U 2 + (2tλ -x + 1)t 2 U + t(t(x + 4) + 1) -(3t 2 (λ -4) + 3t + 1)Z + t(6t + 1)(λ + 2)Z 2 = 0 and that f λ (0, y, t) is -1 + 1 + tyV /(ty) where V is the unique formal power series V ∈ Q[y, λ][[t]] with V (0) = 1 and (λy + y 2 + 1) 2 t 4 V 2 + 4t 2 (6t + 1)(λ + 2)yZ 2 -4ty(3t 2 λ -12t 2 + 3t + 1)Z + t(6t 2 λy 2 + 4t 2 λ + 2t 2 y 3 + 18t 2 y -2tλy + 2ty 2 + 4ty + 2t -y) V + t(4tλy + ty 2 + 16t + 2y + 4) -4(3t 2 λ -12t 2 + 3t + 1)Z + 4(6t + 1)t(λ + 2)Z 2 = 0.
Together with the functional equation, it finally follows that f λ (x, y, t) is algebraic.

For the generating function fλ (x, y, t) of the model with a -1,0 = a -1,1 = a 0,1 = 0, a -1,-1 = 1 2 a 0,-1 = a 1,-1 = a 1,1 = a 0,1 = 1, a 1,0 = λ from family 3b, we have that fλ (0, 0, t) = f λ (0, 0, t) (for combinatorial reasons), and we can use this and the functional equation to certify guessed systems of partial linear differential equations for fλ (x, 0, t) and fλ (0, y, t) which then together with the function equation ( 1) (now with fλ in place of f ) imply that fλ (x, y, t) is D-finite. The equations are somewhat too large to be included here: fλ (x, 0, t) satisfies a differential equation of order 11 with respect to t with polynomial coefficients of respective degrees 82, 90, 110 in x, λ and t, while fλ (0, y, t) satisfies a differential equation of order 11 with respect to t with polynomial coefficients of respective degrees 70, 58, 90 in y, λ and t.

6 Models with Group D8

For the possible values of a u,v such that (ΦΨ) 4 = id, we obtain three essentially different prime ideals after discarding uninteresting or redundant components. One of them is the ideal from Section 4, which appears again because (ΦΨ) 2 = id implies (ΦΨ) 4 = id. The other two define the following families:

Family 4a

Defining equations:

a 1,-1 a -1,1 = a 1,0 a -1,0 , a 1,1 = a 0,1 = a 0,-1 = a -1,-1 = 0 Example: 6 3 2 4

Family 4b

Defining equations:

a 1,1 a -1,-1 = a 1,0 a -1,0 , a 1,-1 = a 0,1 = a 0,-1 = a -1,1 = 0 Example: 2 4 6 3
In family 4a, we must have a 1,-1 = 0 and a -1,1 = 0 for a model to be interesting. But then a 1,-1 a -1,1 = 0 implies also a 1,0 = 0 and a -1,0 = 0 through the first defining equation. Similarly, we can assume for the models in family 4b that a 1,0 , a -1,0 , a 1,1 a -1,-1 all are nonzero.

For family 4a, the orbit sum argument applies and yields

f (x, y, t) = 1 xy [x > y > ] (a 1,-1 x/y -a -1,1 y/x)(a 1,-1 /y -a 1,0 )(a 1,0 x -a -1,0 /x)(a 1,0 x -a -1,1 y/x) a -1,1 a 3 1,0 K(x, y, t)
For family 4b the orbit sum is zero, but it was pointed out by Bostan et al. (2014a) in their Section 6.2 that its D-finiteness can be deduced from the D-finiteness of the corresponding model without multiplicities. If g(x, y, t) denotes the generating function for the Gessel model without multiplicities, we have f (x, y, t) = g a 1,0 a -1,0 x, a 1,1 a 1,0 y, a 1,0 a -1,0 t for the general generating function of models of family 4b. Since g(x, y, t) is known to be algebraic [START_REF] Bostan | The complete generating function for Gessel walks is algebraic[END_REF][START_REF] Bostan | The complete generating function for Gessel walks is algebraic[END_REF], it follows that all the models of family 4b are algebraic.

Theorem 3 The interesting quarter plane models whose group is D8 are precisely those that belong to one of the families described in the table above. All these models are D-finite.

Models with Larger Groups

For n ≥ 5 we failed to compute the prime decomposition of the ideal of relations among the a u,v that ensures (ΦΨ) n = id, as the required calculations become too expensive. However, in our search over all the 30307 quarter plane models with multiplicites in {0, 1, 2, 3} we did encounter, much to our surprise, the following three models that do not belong to any of the families discussed so far. Their group is D10.

1 1 2 1 2 1 1 1 2 1 1 1 1 2 1 2 1 2 1 1 1
The orbit sum is zero for all three models. By guessing, we found that they seem to be algebraic. According to Raschel (2015), they indeed are.

The models on the left and in the middle can be obtained from one another by reversing arrows, so these two models have the same number of walks returning to the origin. We previously conjectured that these two models could be proven to be algebraic using half-orbit sum techniques. [START_REF] Bousquet-Mélou | An elementary solution of Gessel's walks in the quadrant[END_REF] confirms this.

For the model on the right, f (x, 0, t) seems to satisfy an algebraic equation P (x, t, f (x, 0, t)) = 0 for some irreducible polynomial P ∈ Z[x, t, T ] of respective degrees 40, 45, 24 in x, t, T , and f (0, y, t) seems to satisfy an algebraic equation Q(y, t, f (0, y, t)) = 0 for some irreducible polynomial Q ∈ Z[y, t, T ] of respective degrees 64, 45, 24 in y, t, T . We expect that these algebraic equations can be proven by computer algebra in a similar way as the models of family 2b in Section 5.3 above, but this would require immense calculations which we have not carried out.

Using substitutions like in earlier sections, the three models can be used to generate three families of models, which we call Family 5a, 5b, and 5c. The corresponding ideals of defining relations for the a u,v have dimension three. We do not know whether these families completely characterize all the interesting models whose group is D10, nor do we know anything about models for even larger groups. Does there exist for every n ≥ 2 a quarter plane model with multiplicities whose group is D2n?

Partially supported by the Austria FWF grants Y464-N18 and F50-04. Partially supported by the Austria FWF grant F50-04 1365-8050 c 2015