LMPA Marc Prévost
email: marc.prevost@univ-littoral.fr

Tanguy Rivoal

Marc Pr

REMAINDER PAD É APPROXIMANTS FOR HYPERGEOMETRIC SERIES

Keywords: 2010 Mathematics Subject Classification. Primary 41A21; Secondary 41A28 Remainder Padé approximants, hypergeometric function, multiple Charlier polynomials

1. Introduction 1.1. Remainder Padé Approximants. In [START_REF] Prévost | A new proof of the irrationality of ζ(2) and ζ(3) using Padé approximants[END_REF], the first author introduced a new kind of rational approximation, the Remainder Padé Approximants (RPA). They are defined as follows. Let

f (z) = ∞ k=0 a k z k ∈ C[[z]]
; for simplicity of the exposition, in the discussion below f is assumed to be convergent in a disk D of positive radius (1 ). For any integer n ≥ 0, we assume the existence of certain functions b n (z) ∈ C(z) such that the normalized remainder series Φ(z, n) := 1 bn(z) ∞ k=n a k z k admits, for any fixed z ∈ D \ {zeros of all the b n }, an asymptotic expansion as n → +∞ of the form

Φ(z, n) ∼ Φ z 1 n := ∞ k=0 φ k (z) n k ∈ C(z) 1 n .
Let [p/q] Φz(t) denote the ordinary Padé approximant of degrees p and q of the formal power series Φ z (t) at t = 0; it belongs to C(z, t). We thus have

f (z) = n-1 k=0 a k z k + b n (z)Φ n (z) ≈ n k=0 a k z k + b n (z) [p/q] Φz(t) t=1/u(n) . (1.1)
At this stage, the ≈ sign does not have any particular meaning, and the sequence of complex numbers u(n) is n + o(n) but not necessarily equal to n. We say that the rational fraction A(z) B(z)

:= n-1 k=0 a k z k + b n (z) [p/q] Φz(t) t=1/u(n)
is a Remainder Padé Approximant of f (z), where A(z), B(z) ∈ C[z] are of lowest degrees. In certain cases, there exists a sequence u(n) such that the corresponding RPA of f (z) coincides with an ordinary Padé approximant of f (z) at z = 0, say [N/D] f (z) . It is not yet understood for which class of functions f this remarkable fact happens. When it does happen, the ≈ sign in (1.1) can then be understood as follows: the order at z = 0 of the remainder B(z)f (z) -A(z) is larger than or equal to N + D + 1, the order of Padé approximation.

More generally, it also happens that, instead of Padé approximants, we recover some previsously known functional approximations of f (z), or even numerical approximations of certain of its values. For instance, the first author constructed certain sequences of

RPA for ζ(2) = ∞ k=1 1/k 2 and ζ(3) = ∞ k=1 1/k 3 (both can be viewed as values of the polylogarithms ∞ k=1 z k /k s at z = 1
) and proved that they are exactly the sequences of rationals numbers used by Apéry [START_REF] Apéry | Irrationality of ζ(2) and ζ(3)[END_REF] to prove the irrationality of ζ(2) and ζ(3). In [START_REF]Nombres d'Euler, approximants de Padé et constante de Catalan[END_REF], the second author used a modification of these RPA to produce a sequence of fast converging rational approximations for Catalan's constant G = ∞ k=0 (-1) k /(2k + 1) 2 . In [START_REF] Prévost | Remainder Padé approximants for the exponential function[END_REF], using certain discrete multiple orthogonal polynomials introduced in [START_REF] Arvesú | Some discrete multiple orthogonal polynomials[END_REF], which generalize the classical Charlier orthogonal polynomials, we proved that certain RPA of e z provide the Padé approximants of the function e z .

All these examples provide various meanings to the ≈ sign in (1.1) and we use the vocable RPA phenomenon to cover such non-trivial properties.

We recall the notion of type II (diagonal) Padé approximants. For a given family (F j (X)) j=1,...,r of formal series in C[[X]] and any integers p and q such that p ≥ (r-1)q ≥ 0, there exist (by linear algebra) some polynomials P 1 (x), . . . , P r (x) and Q(X) in C[X], not all 0, such that deg(P

j ) ≤ p, deg(Q) ≤ r • q and Q(X)F j (X) -P j (X) = O(X p+q+1 ), 1 ≤ j ≤ r.
The fractions P j /Q are unique and are by definition the type II Padé approximants of (F j (X)) j=1,...,r of parameters (p, q). When r = 1, we recover the ordinary Padé approximants [p/q] of F 1 (X).

1.2. The results. The goal of this paper is to show that the RPA phenomenon occurs (in various generality) for the hypergeometric series

E a (z) := 1 F 1 (1; a; z) = ∞ k=0 z k (a) k and 2 F 0 (b, 1; z) = ∞ k=0 (b) k z k .
We assume that a ∈ C \ Z ≤0 . The former series generalizes the exponential function

1 F 1 (1; 1; z) = exp(z), while the latter generalizes Euler's series 2 F 0 (1, 1; z) = ∞ k=0 k!z k .
The variable z is a formal one in Theorems 1 and 2 below. However, it is important to have in mind that, though the series E a (z) is an entire function of z, the 2 F 0 (b, 1; z) diverges for every z. To overcome this divergence, we shall first work the function

E b (z) := 1 Γ(b) ∞ 0 u b-1 1 -uz e -u du
defined for (b) > 0 and analytic for z ∈ C \ [0, +∞). It admits an asymptotic expansion E b (z) as z → 0 in any angular sector centered at 0 that does not contain [0, +∞). It turns out that E b (z) = 2 F 0 (b, 1; z) in such sectors.

To state our results, we need to introduce a certain number of notations. The parameters α j , j = 1, . . . , r, are pairwise distinct non-zero complex numbers, and we set α := (α j ) j=1,...,r . For any z ∈ C and any t ∈ C such that 1/t / ∈ N, we define

Φ z (t) := ∞ k=0 z k (1 -1/t) k .
For any j = 1, . . . , r and any n ≥ 0, we have the trivial identity

∞ k=0 (α j z) k (a) k = n-1 k=0 (α j z) k (a) k + (α j z) n (a) n Φ α j z -1 a -1 + n . (1.2)
In [START_REF] Prévost | Remainder Padé approximants for the exponential function[END_REF], it has been proved that, for any fixed z ∈ C, Φ z (t) admits an asymptotic expansion Φ z (t) := ∞ k=0 ϕ k (-z)t k as t → 0, t < 0, where the coefficients ϕ k (z) are Touchard exponential polynomials of degree k. The latter are defined by the exponential generating function e z(e X -1) = ∞ k=0 ϕ k (z) X k k! (see [START_REF] Touchard | Nombres exponentiels et nombres de Bernoulli[END_REF]). We provided in [START_REF] Prévost | Remainder Padé approximants for the exponential function[END_REF] explicit expressions for the type II Padé approximants at t = 0 of parameters (rp -1, p) for the formal power series Φ α j z (t) in C[z][[t]] for j = 1, . . . , r: we denote them by P j,α,p (t, z)/Q α,p (t, z), and the expressions of these polynomials in t are given in §2. It turns out that they are also polynomials in z.

We now "replace" Φ α j z -1 a-1+n in (1.2) by the corresponding simultaneous Padé approximant of parameters (rp -1, p) of Φ α j z (t) evaluated at t = -1/(a -1 + n) (this substitution is possible because Q α,p (-1/(a -1 + n), z) does not vanish identically as a function of z). This alteration of (1.2) provides the formal Laurent series

R j,n,α,p (z) := ∞ k=0 (α j z) k (a) k - n-1 k=0 (α j z) k (a) k - (α j z) n (a) n P j,α,p ( -1 a-1+n , z) Q α,p ( -1 a-1+n , z) , j = 1, . . . , r,
and we want to know in which sense they are ≈ 0. By construction, for every j, the rational fraction

n-1 k=0 (α j z) k (a) k + (α j z) n (a) n P j,α,p -1 a-1+n , z Q α,p -1 a-1+n , z (1.3) is a RPA of E a (α j z).
The degrees of its numerator and denominator are a priori bounded by n + rp -1 and rp respectively. A better bound holds for the numerator ((i) below) and it is crucial to prove (ii). (Throughout the paper, O(z c ) denotes a Laurent series in z with order equal to c at z = 0.)

Theorem 1. Let us fix the integers r ≥ 1, n ≥ 1, p ≥ 0, such that n ≥ p, a ∈ C \ Z ≤0 , α 1 , . . . , α r ∈ C pairwise distinct.
(i) The degrees of the numerator and denominator of the RPA in (1.3) are bounded by n + (r -1)p -1 and rp respectively, and moreover R j,n,α,p (z) = O(z n+rp ).

(ii) The collection of the RPA in (1.3), j = 1, . . . , r, is the type II Padé approximant of (E a (α j z)) j=1,...,r of parameters (n + (r -1)p -1, p).

(iii) For r = 1 and α 1 = 1, (ii) gives the following equality of rational fractions:

n-1 k=0 z k (a) k + z n (a) n [p -1/p] Φz(t) t=-1/(a-1+n) = [n -1/p] Ea(z) ,
For a generic value of a, our proof needs the assumption that n ≥ p. However, if a = 1, it can be removed in the statement of Theorem 1, and we simply recover the results from [START_REF] Prévost | Remainder Padé approximants for the exponential function[END_REF]. The explanation of this fact is Identity (2.5) in §2.

Let us now describe our results for the function E b (z). For any n ≥ 0, b such that (b) > 0 and z ∈ C \ [0, +∞), we have the identity (proved in §3)

E b (z) = n-1 k=0 (b) k z k + (b) n z n Ψ z 1 b -1 + n (1.4)
where

Ψ z (t) = 1 Γ(1 + 1/t) ∞ 0 u 1/t 1 -uz e -u du
is defined for any z ∈ C \ (0, +∞) and any t such that (1/t) > -1 (in particular, for any t > 0). We shall prove in §3 that when z < 0, Ψ z (t) admits an asymptotic expansion as t → 0, t > 0, given by ( 2)

Ψ z (t) := - 1 z ∞ k=0 ϕ k (1/z)t k+1 = - t z Φ -1/z (t). (1.5)
It then immediately follows that the Padé approximant [p/p] at t = 0 of the formal power series Ψ z (t

) ∈ C[1/z][[t]] is -t z P j,1,p (t,-1/z) Q 1,p (t,-1/z)
where the polynomials are those used to construct the Padé approximants [p -1/p] of Φ z (t) when r = 1.

In (1.4), we now "replace" E b (z) by E b (z), and Ψ z

1 b-1+n by the Padé approximant [p/p] of Ψ z (t) evaluated at t = 1/(b -1 + n) (this subtitution is possible because Q 1,p (1/(b -1 + n), 1/z)
does not vanish identically as a function of z). This alteration of (1.4) provides the formal Laurent series

S n,p (z) := E b (z) - n-1 k=0 (b) k z k + (b) n z n-1 (b -1 + n) P j,1,p ( 1 b-1+n , -1 z ) Q 1,p ( 1 b-1+n , - 1 z )
, 2 This identity is reminiscent of André's duality between solutions of E-operators at z = 0 and z = ∞;

see [START_REF] André | Séries Gevrey de type arithmétique, I. Théorèmes de pureté et de dualité[END_REF]. The former involve function like E a (z), while the latter involve functions like E b (1/z), with a, b ∈ Q.

and we want to know in which sense it is ≈ 0. Because of (1.4), the rational fraction

n-1 k=0 (b) k z k - (b) n z n-1 (b -1 + n) P j,1,p ( 1 b-1+n , -1 z ) Q 1,p ( 1 b-1+n , -1 z ) (1.6)
can be viewed as a RPA of E b (z) in an extended sense. The degrees of its numerator and denominator are a priori bounded by n + p -1 and p respectively. Again, a better bound holds for the numerator ((i) below) and again it is crucial to prove (ii).

) k z k + (b) n z n [p/p] Ψz(t) t=1/(b-1+n) = [n -1/p] E b (z) .
It seems that a RPA phenomenom similar to (ii) in Theorem 1 does not hold for any family ( E b (β j )) j=1,...,r with r ≥ 2.

The proofs of Theorems 1 and 2 are given in §2 and §3 respectively.

Proof of Theorem 1

Proof of (i). We recall the expressions of the type II approximants of the functions Φ α j z (t) (obtained in [START_REF] Prévost | Remainder Padé approximants for the exponential function[END_REF]):

P j,α,p (t, z) = (-t) rp-1 p k 1 ,...,kr=0 r i=1 p k i (-α i z) p-k i (-1/t) k 1 +•••+kr k 1 +•••+kr i=1 (α j z) i-1 (-1/t) i , (2.1)
and

Q α,p (t, z) = (-t) rp p k 1 ,...,kr=0 r i=1 p k i (-α i z) p-k i (-1/t) k 1 +•••+kr . (2.
2)

The P 's have degree rp -1 in both t and z, while the Q's have degree rp in both t and z. Moreover, Q α,p (0, 1) = 1.

For the sake of simplicity, we set

K r = k 1 + k 2 + • • • + k r . Substituting -1/(a -1 + n) for t (which is possible because a /
∈ Z ≤0 and n ≥ 1), we get

Q α,p -1 a -1 + n , z = (a + n -1) -rp p k 1 ,...,kr=0 r i=1 p k i (-α i z) p-k i (a -1 + n) Kr and P j,α,p -1 a -1 + n , z = (a -1 + n) 1-rp (a) n-1 (α j z) n × p k 1 ,...,kr=0 r i=1 p k i (-α i z) p-k i (a -1 + n) Kr S n+Kr-1 (α j z) -S n-1 (α j z)
where

S n (z) = n k=0 z k /(a) k . (Note that Q α,p -1
a-1+n , z is clearly not identically zero, as stated in the Introduction.)

Hence, after some simplifications, the RPAs of the functions E α j (z), j = 1, . . . , r, are given by

n-1 k=0 (α j z) k (a) k + (α j z) n (a) n P j,α,p ( -1 a-1+n , z) Q α,p ( -1 a-1+n , z) = p k 1 ,...,kr=0 r i=1 p k i (-α i z) p-k i (a -1 + n) Kr S n+Kr-1 (α j z) p k 1 ,...,kr=0 r i=1 p k i (-α i z) p-k i (a -1 + n) Kr . (2.3)
At first sight, the numerator and denominator of this rational fraction have degrees bounded by n + rp -1 and rp in z respectively, but a better estimate holds for the numerator. Indeed, we shall now prove that its degree is bounded by n + (r -1)p -1. Indeed, with

K j,r = k 1 + • • • + k j-1 + k j+1 + • • • + k r and k j,r = (k 1 , .
. . , k j-1 , k j+1 , . . . , k r ), the numerator of (2.3) can be expressed as

p k 1 ,...,k j-1 ,k j+1 ,...,kr=0 r i=1,i =j p k i (-α i z) p-k i (a -1 + n) K j,r A j,k j,r (z) (2.4) with A j,k j,r (z) = p k j =0 p k j (-α j z) p-k j (a -1 + n + K j,r ) k j S n+Kr-1 (α j z).
We claim that the degree in z of the polynomial A j,k j,r (z) is at most max(p, n + K j,r ) -1.

To prove this observation, which is non-trivial, we first observe that, for any integers m, q such that m ≥ q ≥ 0,

t a-1 S m (xt) (q) |t=1 = x q S m-q (x) + U q (x) (2.5)
where

U q (x) := q-1 j=0
(-1) q-j (1 -a) q-j x j .

(Above and below in this section, all the differentiations are with respect to t.) Observe that U q (x) ≡ 0 if a = 1, but in general U q (x) has degree q -1 in x. Since (a

-1 + m) k = (-1) k (t -a+1-m ) (k)
|t=1 , it follows that

A j,k j,r (z) = (-1) p p k=0 p k t -a+1-n-K j,r (k) |t=1 t a-1 S n+K j,r +p-1 (α j zt) (p-k) |t=1 + p k=0 (-1) p-k p k (a -1 + n + K j,r ) k U p-k (α j z) = (-1) p t -n-K j,r S n+K j,r +p-1 (α j zt) (p) |t=1
(2.6)

+ p k=0 (-1) p-k p k (a -1 + n + K j,r ) k U p-k (α j z). (2.7)
The polynomial in (2.7) obviously has degree at most p -1 in z. On the other hand

t -n-K j,r S n+K j,r +p-1 (α j zt) (p) |t=1 = n+K j,r +p-1 j=0 (α j z) j (a) j t j-n-K j,r (p) |t=1 = (-1) p n+K j,r +p-1 j=0 (α j z) j (a) j (n + K j,r -j) p .
Since (n + K j,r -j) p = 0 for any j ∈ {n + K j,r , . . . , n + K j,r + p -1}, it follows that the polynomial in (2.6) has degree at most n + K j,r -1, so that the degree of A j,k j,r (z) is at most max(p, n + K j,r ) -1 as claimed. Since n ≥ p by assumption, we deduce that the degree of the numerator of the RPAs given by (2.4) is at most n + (r -1)p -1, as expected.

To prove that R j,n,α,p (z) = O(z n+rp ), we have to multiply the denominator of (2.3) by the series ∞ k=0 (α j z) k (a) k and prove that the product is the sum of a polynomial of degree n + (r -1)p -1 and a function which is O(z n+rp ). This amounts to proving that the coefficient c m of z m in the Taylor expansion of

p k 1 ,...,kr=0 r i=1 p k i (-α i z) p-k i (a -1 + n) Kr • ∞ k=0 (α j z) k (a) k is zero for m ∈ {n + (r -1)p, . . . , n + rp -1}. Let us assume that m ≥ n + (r -1)p. Then c m = p k 1 ,...,kr=0 r i=1 p k i (-α i ) p-k i (a -1 + n) Kr α k j (a) k (2.8)
where in (2.8) the integer k :

= m + K r -rp is ≥ 0 (because m ≥ n + (r -1)p ≥ rp).
We now also assume that m ≤ n + rp -1. In that case, we use the identity (a

-1 + n) Kr (a) k = 1 (a) n-1 t a+n+Kr-2 (n+rp-m-1)
|t=1 in (2.8) and get

c m = 1 (a) n-1 p k 1 ,...,kr=0 r i=1 p k i (-α i ) p-k i α m-rp+Kr j t a+n+Kr-2 (n+rp-m-1) |t=1 = 1 (a) n-1 t a-1+n-rp r i=1 (α j -α i t) p (n+rp-m-1)

|t=1

.

Since t = 1 is a zero of order p of the product r i=1 (α j -α i t) p , it follows that c m = 0 when n + rp -m -1 ≤ p -1, i.e when m ≥ n + (r -1)p. This complete the proof of (i).

Proof of (ii). By unicity of type II Padé approximants, the bounds for the degrees of the simultaneous RPAs and the order of R j,n,α,p at 0 ensure that (2.3) provides the expression of the type II Padé approximants of the family (E a (α j z)) j=1,...,r of parameters (n + (r -1)p -1, p) when n ≥ p.

Proof of Theorem 2

Let b such that (b) > 0. We first prove (1.4) of the Introduction, i.e. that for any n ≥ 0 and any z ∈ C \ [0, +∞),

E b (z) = n-1 k=0 (b) k z k + (b) n z n Ψ z 1 b -1 + n .
This follows from

(b) n z n Ψ z 1 b -1 + n = (b) n z n Γ(b + n) ∞ 0 u b+n-1 1 -uz e -u du = (b) n Γ(b + n) ∞ 0 u b-1 1 -uz e -u du + (b) n Γ(b + n) ∞ 0 u b-1 (uz) n -1 1 -uz e -u du = E b (z) - (b) n Γ(b + n) n-1 k=0 z k ∞ 0 u k+b-1 e -u du = E b (z) - n-1 k=0 (b) n Γ(b + k) Γ(b + n) z k = E b (z) - n-1 k=0 (b) k z k .
We now prove the following lemma.

Lemma 1. Let us fix z ∈ (-∞, 0). The function Ψ z (t) admits an asymptotic expansion as t → 0, t > 0, given by

Ψ z (t) := - 1 z ∞ k=0 ϕ k (1/z)t k+1 = - t z Φ -1/z (t),
where the coefficients ϕ k (X) are Touchard exponential polynomials of degree k in X.

Proof. We first use Proposition 2 in [START_REF]On the arithmetic nature of the values of the Gamma function, Euler's constant and Gompertz's constant[END_REF] to get the following alternative expression for Ψ z (t): for any z ∈ (-∞, 0) and t > 0,

Ψ z (t) = ∞ 0 e -u (1 -zu) 1+1/t du. We now set v = ln(1 -uz), so that Ψ z (t) = - 1 z ∞ 0 e (e v -1)/z e -v/t dv (3.1)
for z ∈ (-∞, 0) and t > 0. Since v → e (e v -1)/z is C ∞ at v = 0, non-zero at v = 0, and

∞ 0 e (e v -1)/z e -v/
Ψ z (t) ∼ - 1 z ∞ k=0 d k dv k e (e v -1)/z |v=0 t k+1 = - t z ∞ k=0 ϕ k (1/z)t k , t → 0, t > 0.
This complete the proof of the lemma.

From now on, we view Ψ z (t) as a formal series in C

[1/z][[t]].
As observed in the Introduction, the identity Ψ z (t) = -t z Φ -1/z (t) enables us to immediately get the Padé approximant [p/p] of Ψ z (t) at t = 0. More generally, let β := (β j ) j=1,...,r denote a family of pairwise distinct complex numbers. Then, without further efforts, we even obtain the type II Padé approximants at t = 0 of parameters (rp, p) of the formal power series ( Ψ β j z (t)) j=1,...,r for j = 1, . . . , r. They are given by t β j z P j,1/β,p (t, -1/z)

Q 1/β,p (t, -1/z) , j = 1, . . . , r where
P j,1/β,p (t, -1/z) = (-t) rp-1 p k 1 ,...,kr=0 r i=1 p k i (β i z) k i -p (-1/t) k 1 +•••+kr k 1 +•••+kr i=1 (-β j z) 1-i (-1/t) i , Q 1/β,p (t, -1/z) = (-t) rp p k 1 ,...,kr=0 r i=1 p k i (β i z) k i -p (-1/t) k 1 +•••+kr .
We now substitute 1/(b -1 + n) for t in these expressions, which is possible because (b) > 0 and n ≥ 1. Since Q 1/β,p (1/(b -1 + n), -1/z) is clearly not identically zero, we get simultaneous RPAs of the family (E b (β j z)) j=1,...,r given by the formulas: where T n (z) := n k=0 (b) k z k . We are now in position to complete the proof of Theorem 2. We take r = 1 and β 1 = 1 in (3.2).

n-1 k=0 (b) k (β j z) k -(b) n (β j z) n-1 (b -1 + n) P j,1/β,p ( 1 b-1+n , -1 z ) Q 1/β,p
Proof of (i). It is immediate that the degrees of the numerator and denominator of the rational function in (3.2) are bounded by n -1 and p respectively. We now want to prove that, formally,

E b (z) - n-1 k=0 (b) k z k + (b) n z n-1 (b -1 + n) P j,1,p ( 1 b-1+n , -1 z ) Q 1,p ( 1 b-1+n , - 1 z )
= O(z n+p ).

(

The denominator of the right-hand side of (3. is the sum of a polynomial P (z) of degree ≤ n -1 and a formal power series of order ≥ n + p at z = 0. This in fact means that we want to prove that P (z)/Q(z) is the [n -1/p] Padé approximant of E b (z). This is indeed the case because we recognize Q(z) in (3.4) as being of the form given in [START_REF] Sidi | A new method for deriving Padé approximants for some hypergeometric functions[END_REF] and unicity of Padé approximant enables us to conclude.

Proof of (ii). This is simply a summary of what was proved in (i). We have

E b (z) =

Theorem 2 .

 2 Let us fix b ∈ C \ Z ≤0 , and the integers n ≥ 1, p ≥ 0 such that n ≥ p. (i) The degrees of the numerator and denominator of the RPA in (1.6) are bounded by n -1 and p respectively, and moreover S n,p (z) = O(z n+p ) exists and is finite. (ii) The RPA in (1.6) coincides with the Padé approximant [n -1/p] of E b (z). In other words,

  z) k i -p (-1) Kr T n-Kr-1 (β j z)/(b) z) k i -p (-1) Kr /(b) n-Kr . (3.2)

  prove(3.3) it is enough to prove that the formal series expansion of n-k z k E b (z)

  k z k + (b) n z n ([p/p] Ψz(t) ) |t=1/(b-1+n) + O(z n+p ) (3.5)and the bounds for the degrees show thatn-1 k=0 (b) k z k + (b) n z n ([p/p] Ψz(t) ) |t=1/(b-1+n) = [n -1/p] E b (z) .

It is possible to extend the notion of RPA when f diverges everywhere but is summable in some sense; we present an exemple of this situation in Theorem

1

The above results are a priori proved under the assumption that (b) > 0. They hold under the weaker assumption that b ∈ C \ Z ≤0 because all these identities between rational fractions extend to this case.


Emmanuel Soubies

Laure Blanc-Féraud

Sébastien Schaub

Ellen Van Obberghen-Schilling

IMPROVING 3D MA-TIRF RECONSTRUCTION WITH DECONVOLUTION AND BACKGROUND ESTIMATION

Keywords: Total-Internal Reflection Fluorescence Microscopy, Inverse Problems, SDMM

Total internal reflection fluorescence microscopy (TIRF) produces 2D images of the fluorescent activity integrated over a very thin layer adjacent to the glass coverslip. By varying the illumination angle (multi-angle TIRF), a stack of 2D images is acquired from which it is possible to estimate the axial position of the observed biological structures. Due to its unique optical sectioning capability, this technique is ideal to observe and study biological processes at the vicinity of the cell membrane. In this paper, we propose an efficient reconstruction algorithm for multi-angle TIRF microscopy which accounts for both the PSF of the acquisition system (diffraction) and the background signal (e.g., autofluorescence). It jointly performs volume reconstruction, deconvolution, and background estimation. This algorithm, based on the simultaneous-direction method of multipliers (SDMM), relies on a suitable splitting of the optimization problem which allows to obtain closed form solutions at each step of the algorithm. Finally, numerical experiments reveal the importance of considering the background signal into the reconstruction process, which reinforces the relevance of the proposed approach.

INTRODUCTION

Due to its capability to limit the observed region to a thin layer of a few hundred nanometers (below 100nm up to 1µm), total internal reflection fluorescence (TIRF) microscopy [START_REF] Axelrod | Cell-substrate contacts illuminated by total internal reflection fluorescence[END_REF][START_REF] Axelrod | Total internal reflection fluorescence microscopy[END_REF][START_REF] Poulter | The physical basis of total internal reflection fluorescence (TIRF) microscopy and its cellular applications[END_REF][START_REF] Ml Martin-Fernandez | A "pocket guide" to total internal reflection fluorescence[END_REF] is ideally suited to the observation of biological activities near the cell membrane. This unique optical sectioning is made possible by using an evanescent excitation produced in the total internal reflection regime. Moreover, given a set of TIRF acquisitions with different incident angles (multi-angle TIRF, MA-TIRF) it is possible to reconstruct a 3D volume through dedicated algorithms [START_REF] Reichert | Total internal reflection fluorescence (tirf) microscopy. i. modelling cell contact region fluorescence[END_REF][START_REF] Jeffrey S Burmeister | Quantitative analysis of variable-angle total internal reflection fluorescence microscopy (va-tirfm) of cell/substrate contacts[END_REF][START_REF] Bence P Olveczky | Mapping fluorophore distributions in three dimensions by quantitative multiple angle-total internal reflection fluorescence microscopy[END_REF][START_REF] Rohrbach | Observing secretory granules with a multiangle evanescent wave microscope[END_REF][START_REF] Loerke | Quantifying axial secretory-granule motion with variable-angle evanescent-field excitation[END_REF][START_REF] Saffarian | Differential evanescence nanometry: live-cell fluorescence measurements with 10-nm axial resolution on the plasma membrane[END_REF][START_REF] Liang | A Bayesian method for 3D estimation of subcellular particle features in multi-angle tirf microscopy[END_REF][START_REF] Yang | 3-d reconstruction of microtubules from multi-angle total internal reflection fluorescence microscopy using bayesian framework[END_REF][START_REF] Cardoso | Axial nanoscale localization by normalized total internal reflection fluorescence microscopy[END_REF][START_REF] Soubies | A 3D model with shape prior information for biological structures reconstruction using Multiple-Angle Total Internal Reflection Fluorescence Microscopy[END_REF][START_REF] Boulanger | Fast high-resolution 3D total internal reflection fluorescence microscopy by incidence angle scanning and azimuthal averaging[END_REF][START_REF] Daniel R Stabley | Real-time fluorescence imaging with 20 nm axial resolution[END_REF][START_REF] Cardoso | Topography of cells revealed by variable-angle total internal reflection fluorescence microscopy[END_REF][START_REF] Jin | High-resolution 3d reconstruction of microtubule structures by quantitative multi-angle total internal reflection fluorescence microscopy[END_REF][START_REF] Zheng | Threedimensional super-resolved live cell imaging through polarized multiangle tirf[END_REF]. This gives access to a precise axial localization of the molecules within the observed layer, with a precision of the order of 30nm. To further improve the quality of the reconstruction, it is important to consider both the PSF of the system as well as a background signal which corrupt the acquired data. In a prior work [START_REF] Soubies | Nanometric axial resolution of fibronectin assembly units achieved with an efficient reconstruction approach for multi-angle-tirf microscopy[END_REF], we proposed a joint reconstruction and deconvolution algorithm for MA-TIRF microscopy. However, the background sig-E.S. work was funded by an ERC grant (ERC-692726-GlobalBioIm). Work on TIRF was supported by the French Government (National Research Agency) through the Investments for the Future LABEX SIGNALIFE (ANR-11-LABX-0028-01) and the ANR-FNS Collaborative Research Programme (ANR-16-CE93-0005-01). nal is estimated and subtracted from the acquired data within a preprocessing step, as proposed by the authors in [START_REF] Boulanger | Fast high-resolution 3D total internal reflection fluorescence microscopy by incidence angle scanning and azimuthal averaging[END_REF]. In the present communication, we extend the algorithm in [START_REF] Soubies | Nanometric axial resolution of fibronectin assembly units achieved with an efficient reconstruction approach for multi-angle-tirf microscopy[END_REF] to jointly perform 3D reconstruction, deconvolution, and background estimation. A reconstruction obtained with the proposed method is presented in Fig- ure 1. Moreover, we analyze the effect of the background signal on the reconstructed volume and show that it is crucial to estimate it in order to produce relevant reconstructions.

With deconvolution

Without deconvolution 5µm

IMAGE FORMATION MODEL

A TIRF acquisition (2D) g ∈ R M follows the generic fluorescent microscopy model

gm = (h * (wf ))(xm, z fp ) + bm ( 1
)
where gm is the mth component of the discrete vector g, w ∈ L2(R 3 ) is the illumination, h ∈ L2(R 3 ) denotes the point-spread function (PSF) (including the integration over the mth pixel), and f ∈ L2(R 3 ) is the unknown 3D fluorophores density map. Model (1) maps the continuous object f to the discrete measurements g. This is achieved by sampling the continuous quantity (h * (wf )) at the lateral positions {xm ∈ R 2 } M m=1 of the camera pixels and the axial position z fp of the objective focal plane. Fi-nally b ∈ R M is a vector modeling the background signal (e.g., auto-fluorescence) and noise (discrete setting).

The specificity of TIRF microscopy lies in the illumination process. It proceeds by illuminating the sample with an incident angle greater than the critical angle of reflection, so that the light is totally reflected by the coverslip (total internal reflection regime). This phenomenon produces the evanescent excitation [START_REF] Axelrod | Cell-substrate contacts illuminated by total internal reflection fluorescence[END_REF][START_REF] Axelrod | Total internal reflection fluorescence microscopy[END_REF][START_REF] Ml Martin-Fernandez | A "pocket guide" to total internal reflection fluorescence[END_REF] w(x, z) = w(z) = w0 exp (-zp) ,

for z > 0, where w0 is the intensity at the interface (i.e., z = 0) and p characterizes the inverse of the penetration depth of the evanescent wave. It has the expression

p = 4πni λexc sin 2 (α) -sin 2 (αc) , (3)
where ni (nt, respectively) is the refractive index of the incident medium (the transmitted medium, respectively), λexc is the wavelength of the incident beam, αc = asin(nt/ni) is the critical angle from which there is total internal reflection of the light, and α > αc is the incident angle of the illumination beam. The intensity w0 depends on the angle α as well as the polarization of the incident beam and is assumed spatially constant. A theoretical expression can be derived from Fresnel's equations [START_REF] Ml Martin-Fernandez | A "pocket guide" to total internal reflection fluorescence[END_REF][START_REF] Boulanger | Fast high-resolution 3D total internal reflection fluorescence microscopy by incidence angle scanning and azimuthal averaging[END_REF]. Because of the fast decay of the excitation field (2) in depth (with respect to the z variable), only a thin layer of the sample at the vicinity of the glass coverslip is imaged (between 100nm to 1µm, depending on the incident angle). Within this thin region, the variations of the PSF in depth are small and we can thus assume that h is constant along z (i.e., h(x, z) = h(x)). This allows to decouple the spatial convolution and the TIRF excitation in Model (1),

gm = ∞ 0 R 2 h(xm -x) w(z)f (x, z) dxdz + bm = ∞ 0 w(z) R 2 h(xm -x)f (x, z) dxdz + bm, (4)
which is the key ingredient to derive the SDMM that we propose in Section 3.

A single TIRF image does not provide a quantitative information on the location of the fluorescent probes within the observed layer. However, given a set of Q TIRF images {gq ∈ R M } Q q=1 acquired with the incident angles {αq > αc} Q q=1 , it is possible to numerically reconstruct the 3D biological sample f . To that end, we define the discrete version of Model (4) for the incident angle αq as

gq = TqHf + b, (5)
where f ∈ R N 0 denotes the 3D discrete fluorophore density map, H ∈ R N ×N is a convolution operator that convolves each z-slice of the volume f with a sampled version of the 2D kernel h, and Tq ∈ R N ×M is the discrete TIRF operator associated to the incident angle αq. It performs a weighted sum of the slices of Hf , where the weights are samples of w associated to αq. Note that, in contrast, the background signal b is considered independent from the incident angle (i.e., independent of the TIRF excitation process).

JOINT DECONVOLUTION, RECONSTRUCTION AND BACKGROUND ESTIMATION

We aim at estimating both the 3D fluorophore density f and the (2D) background signal b. We formulate the reconstruction problem as the joint optimization problem

( f , b) ∈ arg min f ∈R N b∈R M Q q=1 1 2 TqHf + b -gq 2 2 + µ0 ∇b 2 2 + µ1R1(Lf ) + µ2R2(f ) . (6)
Because the background has smooth variations, we consider the first order Tikhonov regularizer ∇b 2 2 for b. For f , the problem is regularized using the combination of two terms: R1(L • ) and R2( • ). The only requirement we make is that R1 and R2 are simple in the sense that their proximal operators can be computed efficiently. In this work, we will consider the Hessian Schatten-norm of order one [START_REF] Lefkimmiatis | Poisson image reconstruction with hessian schatten-norm regularization[END_REF], defined by R1(L • ) = He • ,1, as well as the sum of a nonnegativity constraint and a sparsity promoting functional

R2 = i ≥0 + • 1, where i ≥0 (f ) = {0 if f ∈ R N ≥0 ; +∞ otherwise}.
It is noteworthy that the proposed approach can also deal with the popular total variation (TV) [START_REF] Leonid I Rudin | Nonlinear total variation based noise removal algorithms[END_REF] regularizer (i.e., R1 = • 2,1 and L = ∇).

One can easily check that the objective function in [START_REF] Jeffrey S Burmeister | Quantitative analysis of variable-angle total internal reflection fluorescence microscopy (va-tirfm) of cell/substrate contacts[END_REF], with the regularizers considered in this work, is convex in (f , b). Hence, we benefit from the variety of existing convex optimization algorithms to tackle this minimization problem. We deploy the SDMM1 [START_REF] Combettes | Proximal Splitting Methods in Signal Processing[END_REF] to solve the augmented problem

( f , b, {ûi} 2 i=0 ) ∈ arg min f ,b {û i } 2 i=0 Q q=1 1 2 Tqu0 + b -gq 2 2 + µ0 ∇b 2 2 + µ1R1(u1) + µ2R2(u2) such that u0 = Hf , u1 = Lf , u2 = f . (7)
The iterates are summarized in Algorithm 1. The parameters {ρi} 2 i=0 are the Lagrangian multipliers that act as preconditioners. In practice, we tune them empirically to optimize the convergence speed. Lines 7 to 9 require to compute different proximity operators. For R1 = • * ,1, it admits a closed form expression which can be efficiently computed through singular value decomposition [START_REF] Lefkimmiatis | Poisson image reconstruction with hessian schatten-norm regularization[END_REF]. The proximal operator of R2 = i ≥0 + • 1 has a very simple component-wise expression given by

prox µ 2 ρ 2 R 2 (v) n = max (0, vn -µ2/ρ2) . (8)
Finally, for

D = 1 2ρ 0 Q q=1 Tq • + b k -gq 2 2 , we have prox D (v) = Q q=1 T T q Tq + ρ0I -1 Q q=1 T T q (gq -b k ) + ρ0v
Because the TIRF excitation (2) is independant from the lateral variable x, the operator Q q=1 T T q Tq is block diagonal with identical blocks B ∈ R Nz ×Nz defined by

[B]ij = Q q=1 w 2 0,q exp(-zipq) exp(-zjpq), (9)
where {zi} Nz i=1 are the axial discretization points and pq (w0,q, respectively) is given by (3) ((2), respectively) with α = αq. Hence Algorithm 1 SDMM for solving Problem [START_REF] Bence P Olveczky | Mapping fluorophore distributions in three dimensions by quantitative multiple angle-total internal reflection fluorescence microscopy[END_REF].

Require: f 0 ∈ R N , b 0 ∈ R M , ρ0 > 0, ρ1 > 0, ρ2 > 0 1: u 0 0 = Hf 0 2: u 0 1 = Lf 0 , u 0 2 = f 0 3: w 0 i = u 0 i , ∀i ∈ {0, 1, 2} 4: C = ρ0H T H + ρ1L T L + ρ2I 5: k = 0 6: while (not converged) do 7: u k+1 0 = prox 1 2ρ 0 Q q=1 Tq • +b k -gq 2 2 Hf k + w k 0 /ρ0 8: u k+1 1 = prox µ 1 ρ 1 R 1 Lf k + w k 1 /ρ1 9: u k+1 2 = prox µ 2 ρ 2 R 2 f k + w k 2 /ρ2
10:

f k+1 = C -1 H T ρ0u k+1 0 -w k 0 +L T ρ1u k+1 1 -w k 1 + ρ2u k+1 2 -w k
w k+1 1 = w k 1 + ρ1(Lf k+1 -u k+1 1 )
13:

w k+1 2 = w k 2 + ρ2(f k+1 -u k+1 2 )
14:

b k+1 = (LI + µ0∇ T ∇) -1 L l=1 g l -T l u k+1 0 15: k = k + 1 16: end while ( Q q=1 T T q Tq + ρ0I
) is also block diagonal and its inversion requires to compute once and store the inverse of the small size bloc (B + ρ0I) (typically smaller than 50 × 50).

Regarding the linear steps in Lines 10 and 14, they can be easily solved in the Fourier domain. Indeed, considering periodic boundary conditions, both He T He (i.e., L T L in the present paper) and ∇ T ∇ are convolution operators. Hence, (ρ0H T H + ρ1L T L + ρ2I) and (LI + µ0∇ T ∇) are also convolution operators which can be inverted in the Fourier domain.

IMPORTANCE OF BACKGROUND ESTIMATION

In this section we study the influence of the background signal on the reconstructed volume. Our main result is stated in Proposition 4.1.

Proposition 4.1. Let T : L2(R) → R be the 1D TIRF operator such that T {f }(α) = w0(α) ∞ 0 f (z) exp(-zp(α)) dz. Let A = {αq > αc} Q
q=1 be a set of incident angles and g ∈ R Q be the constant vector such that gq = bg > 0 for all q ∈ {1, . . . , Q}. Then, the solution of the optimization problem

f = arg min f ≥0 1 2 Q q=1 (T {f }(αq) -gq) 2 ( 10
)
is given by

f = Q q=1 w0(αq) Q q=1 w 2 0 (αq) bg δ ( 11
)
where δ denotes the Dirac impulse.

Proof. The proof can be found in [START_REF] Soubies | Sur quelques problèmes de reconstruction en imagerie MA-TIRF et en optimisation parcimonieuse par relaxation continue exacte de critères pénalisés en norme-l0[END_REF]Proposition 3.1].

From Proposition 4.1, we get that the best interpretation (in the sense of positive least-squares) of a constant signal by the TIRF model (i.e., (4) without the convolution operator) is a Dirac mass localized at z = 0. It is noteworthy to mention that the TIRF operator (i.e., the first integral in (4)) only depends on z, and that the result of Proposition 4.1 is thus applicable in 3D. In practice, a TIRF acquisition contains the signal resulting from the excitation of the sample plus a constant signal (background). Hence, trying to reconstruct the 3D sample from such data without considering the background signal will produce a strong signal at z = 0 and will introduce errors in the axial localization of the fluorescent probes. One solution to this problem, adopted by Boulanger et al [START_REF] Boulanger | Fast high-resolution 3D total internal reflection fluorescence microscopy by incidence angle scanning and azimuthal averaging[END_REF], consists in estimating the background signal using an additional "dark image" and subtracting this background value to the acquired data before reconstruction. Note that, to account for the uncertainty of this estimation, the authors in [START_REF] Boulanger | Fast high-resolution 3D total internal reflection fluorescence microscopy by incidence angle scanning and azimuthal averaging[END_REF] relax the nonnegativity constraint i ≥0 by the distance to the set of nonnegative vectors. In contrast, we propose in this work to estimate the background signal jointly with the reconstruction of the fluorophore density map.

NUMERICAL EXPERIMENTS Simulation Settings

We consider the synthetic sample depicted in Figure 2 (top left). It has been obtained from the list of molecules positions of the MT0 dataset created for the SMLM challenge [START_REF] Sage | Super-resolution fight club: A broad assessment of 2d & 3d single-molecule localization microscopy software[END_REF]. More precisely, these molecules have been binned into a 3D discrete Reconstruction Results and Discussion From these simulated data, we performed several reconstructions by considering the background estimation (B) or not (NB), as well as the deconvolution (D) or not (ND). The reconstruction algorithm has been implemented using the GlobalBioIm Library [START_REF] Unser | Global-BioIm: A unifying computational framework for solving inverse problems[END_REF]. Regularization parameters have been empirically tuned so as to maximize the performances.

No background estimation Background estimation 0 -100 nm Color-coded depth maps are depicted in Figure 3. One can clearly observe the improvement obtained when the deconvolution is activated, in particular with the line profiles presented in the bottom graph of Figure 3. Then, in order to analyze the effect of the background estimation on the reconstructed volume, we have computed mean axial error maps between the reconstructed volumes and the ground-truth. They are provided in Figure 4. For structures that are close to the glass coverslip (i.e., z = 0), the mean axial error is reasonably low, even without background estimation (thin white arrow). However, for deeper structures, estimating the background signal results in reconstructions which are more accurate (thick magenta arrows). These observations are corroborated with the (y,z) projections depicted in Figure 5. These projections reveal that the axial localization of the observed structures is significantly improved when the background signal is estimated. Moreover, when the background signal is not estimated, the reconstructed volume presents a strong signal at the interface (white arrow), which illustrates the statement of Proposition 4.1. Finally, an example of reconstruction with real data is presented in Figure 1. Here also, one can appreciate the effect of the deconvolution. We refer the reader to the companion paper [START_REF] Soubies | Nanometric axial resolution of fibronectin assembly units achieved with an efficient reconstruction approach for multi-angle-tirf microscopy[END_REF] for further experiments on real data.

CONCLUSION

We proposed an efficient reconstruction algorithm for MA-TIRF microscopy which performs at the same time reconstruction, deconvolution, and background estimation. Moreover, we demonstrated that it is crucial to estimate the background signal present in the TIRF acquisitions in order to produce relevant reconstructions.

Fig. 1 :

 1 Fig. 1: MA-TIRF reconstruction (proposed method) of actin filaments. Results with and without considering the convolution operator in the model are presented. The colormap encodes the axial position of the reconstructed structures from 0 to 800nm.

Fig. 2 :

 2 Fig. 2: Simulated data. Top row: color-coded depth representation of the ground-truth (left) and background signal (right). Middle row: zoom corresponding to the yellow square (left) and the corresponding (y,z) projection (average intensity). The yellow line in the (x,y) zoom is the one used for the profiles that are displayed in the bottom graph of Figure 3. Bottom row: examples of TIRF acquisitions for αc < α 1 < α 2 < α 3 < α 4 < αmax. When not specified the scale bar represents 2λexc.

Fig. 3 :

 3 Fig. 3: Color-coded depth representations of the reconstructed volumes within the yellow region in the top-left image of Figure 2. The colormap is the one displayed in Figure 2. Profiles of the (x,y) projections (average intensity) along the yellow line drawn in Figure 2.

Fig. 4 :Fig. 5 :

 45 Fig. 4: Mean axial error maps for the reconstructed volumes with deconvolution.

Extension of the popular alternating direction method of multipliers (ADMM) to minimize the sum of more than two functionals.


Pascale Royer

Pierre Recho

Claude Verdier

On the quasi-static effective behaviour of poroelastic media containing elastic inclusions

The aim of the present study is to derive the effective quasi-static behaviour of a composite medium, made of a poroelastic matrix containing elastic impervious inclusions. For this purpose, the asymptotic homogenisation method is used. On the local scale, the governing equations include Biot's model of poroelasticity in the porous matrix and Navier equations in the inclusions, with elastic properties of the same order of magnitude. Biot's diphasic model of poroelasticity is obtained on the macroscopic scale, but with effective parameters that are strongly impacted by the distribution of inclusions, even at low volume fraction. The impact on fluid flow is strictly geometrical, showing that the inclusions do not play the role of a porous network.

I. Introduction

Composites made of a porous matrix reinforced with solid impervious inclusions occur in several engineering disciplines, involving natural media such as geomaterials [START_REF] Rice | Deformation of spherical cavities and inclusions in fluid-infiltrated elastic materials[END_REF], biological tissues [START_REF] Rauch | A coupled approach for fluid saturated poroelastic media and immersed solids for modeling cell-tissue interactions[END_REF], [START_REF] Loboa | Mechanobiology of soft skeletal tissue differentiation -a computational approach of a fiber-reinforced poroelastic model based on homogeneous and isotropic simplifications[END_REF], or tumors [START_REF] Xue | A non-equilibrium thermodynamic model for tumor extracellular matrix with enzymatic degradation[END_REF], as well as man-made structures as, cementbased [START_REF] Lemarchand | Effect of Inclusions on Friction Coefficient of Highly Filled Composite Materials[END_REF] or biomimetic materials Rajkhowa et al. (2010). A comprehensive understanding of the overall behaviour of these composites, on the basis of their microstructure, can enhance the knowledge concerning physical scenarios, with respect to key physical properties. In this way, a generalisation of Eshelby's formula is proposed in [START_REF] Berryman | Generalization of Eshelby's Formula for a Single Ellipsoidal Elastic Inclusion to Poroelasticity and Thermoelasticity[END_REF], to give the response of a single ellipsoidal elastic inclusion, in a poroelastic whole space, to a uniform strain imposed at infinity. Some other works focus on fluid flow or solute transport, such as the mathematical model developed in [START_REF] Federico | On the permeability of fibre-reinforced porous materials[END_REF], of the effect of fibre arrangement on the permeability of a porous fibre-reinforced composite, or the macroscopic models obtained in [START_REF] Baltean | Diffusion-Convection in a Porous Medium with Impervious Inclusions at Low Flow Rates[END_REF] by asymptotic homogenisation, for passive solute transport in a rigid medium made of a porous matrix with impervious inclusions.

The present work is aimed at deriving the effective quasi-static mechanical behaviour of a saturated poroelastic medium containing elastic impervious inclusions and is focused on constituents with elastic properties of the same order of magnitude. For this purpose, the method of asymptotic homogenisation is used and Biot's model of poroelasticity is obtained on the macroscopic scale, but with effective parameters that are strongly affected by the distribution of inclusions. The paper is organised as follows. Section II presents a brief description of the homogenisation methodology. Then, homogenisation of a poroelastic medium with elastic inclusions is detailed in Section III, and the derived macroscopic description is commented in Section IV. Finally, Section V presents a summary of the main theoretical results contained in this work and highlights conclusive remarks.

II. Homogenisation method i. Medium under consideration

We consider a periodic medium, of characteristic size L, and made of a fluid saturated elastic porous matrix which contains isolated elastic impervious inclusions. We further denote the spacing between two inclusions by l, and and we formulate the condition of separation of scales by ε = l/L 1. Within the periodic cell Ω, we denote by Ω p the fluid-saturated porous matrix domain, by Ω c the volume occupied by the inclusion, and by Γ their common interface, as depicted in Fig. 1. Using the two characteristic lengths, l and L, and the physical space variable, X, we define two dimensionless space variables: y = X/l, x = X/L, and y and x describe variations on the microscopic and the macroscopic scales, respectively. Invoking the differentiation rule of multiple variables, the gradient operator with respect to X is written as

∇ X = 1 l ∇ y + 1 L ∇ x . (1)
We further introduce the following cell averages

< . > Ω =< . > Ω p + < . > Ω c , < . > Ω ff = 1 | Ω | Ω α . dΩ (α = p, c).
ii. Homogenisation procedure

The methodology firstly consists in writing, in dimensionless form, the governing equations which describe the problem on the periodic unit cell. This dimensionless writing of the equations requires the choice of a characteristic length for the dimensionless writing of space derivatives. We arbitrarily choose L as the reference characteristic length. The dimensionless gradient operator is thus L ∇ X , which by Eq.( 1) is given by

∇ = L ∇ X = ε -1 ∇ y + ∇ x . (2)
The homogenisation method being used is based upon the fundamental assumption that the unknown fields can be written in the form of asymptotic expansions in powers of ε ψ = ψ 0 ( y, x) + εψ1 ( y, x) + ε 2 ψ 2 ( y, x) + ..., (3) in which functions ψ i are Ω-periodic in variable y. The method consists in incorporating the asymptotic expansions in the dimensionless local description, while taking into account the expression of the dimensionless gradient operator Eq.( 2). This leads to approximate governing equations and boundary conditions at the successive orders, which together with the condition of periodicity define well posed boundary value problems within the periodic unit cell, from which functions ψ i can be determined. Existence of solutions requires that volume averaged equations be satisfied. The latter ones actually describe the macroscopic behaviour at successive orders.

III. Quasi-static homogenisation in a poroelastic medium with elastic inclusions i. Dimensionless governing equations on the local scale

The poroelastic matrix (Ω p ) is made of a linear elastic skeleton saturated with a viscous incompressible Newtonian fluid, and its behaviour is described by quasi-static Biot's model [START_REF] Biot | Theory of elasticity and consolidation for a porous anisotropic solid[END_REF] 1 :

                   ∇ • σp = 0, (4) σp = cp : ẽ( u s ) -αp p f , (5)
∇ • v p = -α p : ẽ( ∂ u s ∂t ) -β p ∂p f ∂t , (6)
v p = φ p ( v f - ∂ u s ∂t ) = - Kp µ ∇p f . ( 7
)
The four above equations express the momentum balance, the poroelastic constitutive law, the conservation of fluid mass and Darcy's law, respectively. The distinct quantities involved in the model are the following: σp and p f denote the total stress and the interstitial fluid pressure, respectively; u s is the solid displacement of the porous matrix, while v f and v p stand for the mean fluid velocity within the volume of the micropores and the mean fluid relative velocity within the porous matrix; ẽ( u s ) = 1/2( ∇ u s + ∇ T u s ) is the strain tensor, while cp , αp , β p > 0, and Kp represent the fourth order elastic tensor of the drained porous matrix, the second order symmetric and positive Biot coupling tensor, Biot's bulk modulus and the second order tensor of permeability of the porous matrix, respectively; φ p and µ denote the porosity of the porous matrix, and the fluid viscosity, respectively. The inclusion (Ω c ) is linear elastic and satisfies the Navier equations

∇ • σc = 0, (8) σc = cc : ẽ( u c ), (9)
where σc and u c stand for the solid stress tensor and displacement, respectively, and where cc represents the elastic tensor.

The appropriate conditions over the interface Γ between the porous matrix and the inclusion include the continuity of normal stresses and displacements and the normal mean fluid relative velocity within the porous matrix must be set to zero [START_REF] Mikelic | On the interface law between a deformable porous medium containing a viscous fluid and an elastic body[END_REF]:

     σp • n = σc • n over Γ, (10)
u s = u c over Γ, (11)
v p • n = 0 over Γ, (12)
where n denotes the unit vector giving the normal to Γ exterior to Ω p .

ii. Homogenisation

We consider equations Eqs.( 4)-( 9) and boundary conditions Eqs.( 10)-( 12), and we look for solutions in the form of Eq.(3) for σp , σc , u s , u c , p f and v p . Note that, due to Eq.( 2), the expansion of v p starts with a term in ε -1 . Furthermore, the strain tensors read

ẽ( u α ) = ε -1 ẽy ( u α ) + ẽx ( u α ), (α = p, c), (13)
and consequently the expansions of both stress tensors also start with a ε -1 term. Incorporating the asymptotic expansions and the expressions of the dimensionless gradient operator Eq.( 2) and of the strain tensors Eq.( 13) into Eqs.( 4)-( 12), then identifying terms of same power of ε, leads to boundary value problems at the successive orders.

ii.1 Boundary value problem for

σ-1 p , σ-1 c , u 0 s and u 0 c
Considering the leading order of Eqs.( 4)-( 5) and Eqs.( 8)-( 11), we deduce the following boundary value problem of unknowns u 0 s and

u 0 c :                                    ∂ ∂y j c p ijlm e y lm ( u 0 s ) = 0 within Ω p , (14) ∂ ∂y j c c ijlm e y lm ( u 0 s ) = 0 within Ω c , ( 15
)
c p ijlm e y lm ( u 0 s ) n j = c c ijlm e y lm ( u 0 c ) n j over Γ, (16)
u 0 s i = u 0 c i over Γ, (17) σ-1 p , σ-1 c , u 0 s , u 0 c : periodic in y, (18)
from which it is clear that the displacements u 0 s and u 0 c are constant over the period

u 0 s = u 0 c = u 0 (x). ( 19
)
Since by Eqs.( 5) and ( 9) at O(ε -1 )

σ -1 α ij = c p ijlm e y lm ( u 0 ) (α = p, s), (20)
we consequently get

σ-1 p = σ-1 c = 0. ( 21
)
ii.2 Boundary value problem for v -1 p and p 0 f At the lowest order, Eqs.( 6)-( 7) lead to

                 ∂ ∂y i K p ij µ ∂p 0 f ∂y j = 0 within Ω p , (22)
K p ij µ ∂p 0 f ∂y j n j = 0 over Γ, (23)
p 0 f : periodic in y. (24)
Consequently, we get:

p 0 f = p 0 f ( x), (25)
v -1 p = 0. ( 26
)
ii.3 Boundary value problem for σ0 p , σ0 c , u 1 s and u 1 c

We now consider the second order of Eqs.( 4)-( 5) and Eqs.( 8)-( 11), from which we deduce the following system of unknowns u 1 s and u 1 c :

                                         ∂ ∂y j c p ijlm [e y lm ( u 1 s ) + e x lm ( u 0 s )] -α p ij p 0 f = 0 within Ω p , ( 27
) ∂ ∂y j c c ijlm [e y lm ( u 1 c ) + e x lm ( u 0 c )] = 0 within Ω c , ( 28
) (c p ijlm [e y lm ( u 1 s ) + e x lm ( u 0 s )] -α p ij p 0 f ) n j = (c c ijlm (e y lm ( u 1 c ) + e x lm ( u 0 c ))) n j over Γ, (29) u 1 s i = u 1 c i over Γ, (30)
u 1 s , u 1 c : periodic in y. ( 31
)
The above set of equations constitutes a wellposed problem for u 1 s and u 1 c , and by virtue of linearity, the solutions read (Cf. A):

u 1 s i = ω kh p i e x kh ( u 0 ) -π p i p 0 f + ū1 s i ( x), (32) u 1 c i = ω kh c i e x kh ( u 0 ) + ū1 c i ( x), (33)
where ū1 s i ( x) and ū1 c i ( x) are arbitrary functions. Note that, to render the solution unique, we impose that ωp , π p and ωc are with zero average [START_REF] Bensoussan | Asymptotic Analysis for Periodic Structures[END_REF], Sanchez-Palencia (1980):

< ω kh p i > Ω p = 0, < ω kh c i > Ω c = 0, < π p i > Ω p = 0.
By Eqs.( 5) and (9) at O(ε 0 ), we obtain

       σ 0 p ij = c p ijlm [e y lm ( u 1 s ) + e x lm ( u 0 s )] -α p ij p 0 f , ( 34
)
σ 0 c ij = c c ijlm [e y lm ( u 1 c ) + e x lm ( u 0 c )], (35)
and then employing Eqs.( 32)-( 33), we deduce

         σ 0 p ij = (c p ijlm e y lm ( ω kh ) + c p ijkh )e x kh ( u 0 ), -(c p ijkh e y kh ( π) + α p ij )p 0 f , ( 36
)
σ 0 c ij = (c c ijlm e y lm ( ω kh ) + c c ijkh )e x kh ( u 0 ), . (37)
ii.4 Boundary value problem for v 0 p and p 1 f We now consider Eq.( 6) at O(ε -1 ), Eq.( 7) at O(ε O ), and boundary condition Eq.( 12) at O(ε 0 ), from which we deduce the differential system

                 ∂ ∂y i K p ij µ ( ∂p 1 f ∂y j + ∂p 0 f ∂x j ) = 0 in Ω p , ( 38
)
K p ij µ ( ∂p 1 f ∂y j + ∂p 0 f ∂x j ) n i = 0 over Γ, (39)
p 1 f : periodic in y.

(40)

The above set of equations is a well-posed boundary value problem of unknown p 1 f , from which it appears that p 1 f is a linear function of ∇ x p 0 f :

p 1 f = χ p i ∂p 0 f ∂x i + p1 f ( x, t), (41)
where p1 f ( x, t) is an arbitrary function and where < χ p >

Ω p Ω = 0. ( 42
)
Vector χ p k is the specific solution for p 1 f , corresponding to ∂p 0 f /∂x j = δ jk . Then, by Eq.( 7) at O(ε -1 )

v 0 p i = - K p ij µ ( ∂p 1 f ∂y j + ∂p 0 f ∂x j ), (43)
and by Eq.( 41), we deduce the following expression for v 0 p :

v 0 p i = - K p ij µ ( ∂χ k ∂y j + δ jk ) ∂p 0 f ∂x k , ( 44
)
where δ denotes Kronecker's symbol.

ii.5 Macroscopic momentum balance

Let now consider Eqs.( 4), ( 8), (10) at the third order:

                 ∂σ 1 p ij ∂y j + ∂σ 0 p ij ∂x j = 0 within Ω p , ( 45
)
∂σ 1 c ij ∂y j + ∂σ 0 c ij ∂x j = 0 within Ω c , (46) σ 1 p ij n j = σ 1 c ij n j over Γ. (47)
The homogenisation procedure consists now in integrating Eqs.( 45) and ( 46) over Ω p and Ω c , respectively. This leads to a compatibility condition, i.e. a necessary and sufficient condition for the existence of solutions for u 2 s and u 2 c , which further represents the first order macroscopic momentum balance. Invoking Gauss' theorem, integration of Eq.( 45) over Ω p yields

1 | Ω | Γ σ 1 p ij n j dS + 1 | Ω | Ω p ∂σ 0 p ij ∂x j dΩ = 0,
(48) where the contribution over the cell boundaries, δΩ ∩ δΩ p , cancel due to y-periodicity. Then, employing Eq.( 47) and Gauss' theorem, we get

1 | Ω | Γ σ 1 p ij n j dΩ = 1 | Ω | Ω c ∂σ 0 c ij ∂x j dΩ. (49)
Finally, Eq.( 48) becomes

∂ < σ 0 p ij > Ω p ∂x i + ∂ < σ 0 c ij > Ω c ∂x i = 0. ( 50
)
Let us define the total stress σT as

σT = σp in Ω p , σc in Ω c .
Thus, Eq.( 50) is rewritten as

       ∂ < σ 0 T ij > Ω ∂x i = 0, (51)
< σ 0 T ij > Ω =< σ 0 p ij > Ω p + < σ 0 c ij > Ω c . ( 52
)
Finally, using Eqs.( 36)-( 37), we get

< σ 0 T ij > Ω = C eff ijkh e x kh ( u 0 ) -A eff ij p 0 f , (53)
with

         C eff ijkh =< c p ijlm e y lm ( ω kh p ) + c p ijkh > Ω p Ω + < c c ijlm e y lm ( ω kh c ) + c c ijkh > Ω c Ω , ( 54
)
A eff ij =< c p ijkh e y kh ( π p ) + α p ij > Ω p Ω . ( 55
)
The first order momentum balance is thus described by Eqs.( 51), ( 53).

ii.6 Macroscopic mass balance

At the second order, Eqs.( 6) and ( 12) yield

               ∂v 1 p i ∂y i + ∂v 0 p i ∂x i = -α p ij e y ij ( ∂ u 1 s ∂t ) + e x ij ( ∂ u 0 ∂t ) -β p ∂p 0 f ∂t within Ω p , ( 56
)
v 1 p i n i = 0 over Γ. ( 57
)
Integrating Eq.( 56) and invoking Gauss' theorem, while taking boundary condition Eq.( 57) into account, together with the condition of periodicity, yields

∂ < v 0 p i > Ω p ∂x i = -< α p ij (e y ij ( ∂ u 1 s ∂t ) + e x ij ( ∂ u 0 ∂t ) > Ω p Ω -< β p > Ω p Ω ∂p 0 f ∂t . (58)
Using Eq.( 32), the above equation can be written as

∂ < v 0 p i > Ω p ∂x i = -G eff lm e x lm ( ∂ u 0 ∂t ) -B eff ∂p 0 f ∂t , (59) where
G eff lm =< α p ij e y ij ( ω lm p ) + α p lm > Ω p , (60) B eff =< β p -α p ij e y ij ( π p ) > Ω p . (61)
Now, by Eq.( 44), we get

         < v 0 p i > Ω p = - K eff ij µ ∂p 0 f ∂x j , (62)
K eff ik =< K p ij ( ∂χ k ∂y j + δ jk ) > Ω p . ( 63
)
The first order macroscopic mass balance is thus given by Eqs.( 59) and (62).

IV. Macroscopic description

The first order macroscopic description thus consists of by Eqs.( 51), ( 53), ( 59) and ( 62), with the effective properties defined by Eqs.( 54), ( 55), ( 60), ( 61) and ( 63). From its definition, it is clear that tensor C eff ijkh is the effective elastic tensor of the whole empty medium, made of the empty porous matrix and the inclusions. It therefore coincides with the effective elasticity tensor that would be obtained for a two-phase elastic composite and thus possesses all the required symmetry properties that characterise an elastic tensor (e.g. Penta and Gerish (2017)). Now, from the variational formulation associated with the definitions of u 1 s and u 1 c , we show that the coupling tensors Ãeff and Geff are equal (Cf. B), and that the bulk modulus B eff is positive (Cf. C). Furthermore, from the symmetries of cp and of αp , it follows that Ãeff is symmetric. With the above mentioned properties, the macroscopic description, Eqs.( 51), ( 53), ( 59) and ( 62), is a Biot diphasic model of poroelasticity, but in which the effective properties are strongly impacted by the local distribution of inclusions. An illustration of this appears when considering homogeneous materials. The definitions of the effective Biot parameters and permeability then reduce to

           A eff ij = c p ijkh < e y kh ( π p ) > Ω p +(1 -n c )α p ij , B eff = (1 -n c )β p -α p ij < e y ij ( π p ) > Ω p , K eff ij = K p ik < ∂χ j ∂y k + δ kj > Ω p ,
where n c = | Ω c |/| Ω | denotes the inclusion volume fraction. We firstly note that at low inclusion concentration, i.e. when n c ≈ 0, the inclusions still have an impact since Ãeff = αp and B eff = β p . Furthermore, when the porous matrix is incompressible, i.e when αeff p = Ĩ and β eff p = 0, the whole poroelastic composite remains compressible since Ãeff = Ĩ and B eff = 0, even at low inclusion volume fraction. Finally, we see that the permeability is such that K eff ij = K p ik T p kj , where

T p kj =< ∂χ j ∂y k + δ kj > Ω p .
Tensor Tp is a purely geometrical parameter, and we note that φ -1 Tp actually represents the tortuosity (e.g. see [START_REF] Royer | Homogenisation of advective-diffusive transport in poroelastic Mech[END_REF]) associated with the distribution of inclusions.

V. Conclusion

We have thus shown that the first order macroscopic behaviour, i.e. with precision in the order of O(ε), of a poroelastic matrix containing elastic inclusions is described by Biot's diphasic model of poroelasticity

                           ∂ < σ T ij > Ω ∂x i = 0, < σ T ij > Ω = C eff ijkh e x kh ( u) -A eff ij p f , ∂ < v p i > Ω p ∂x i = -A eff lm e x lm ( ∂ u ∂t ) -B eff ∂p f ∂t , < v p i > Ω p = - K eff ij µ ∂p f ∂x j .
The effective parameters are strongly affected by the distribution of inclusions, even at low volume fraction. In the above developments, it is implicitely assumed that the inclusion size, l c , is of same order of magnitude as the distance between two inclusions. Considering the case of low inclusion concentration, l c l, would not modify the macroscopic behaviour. But, since this introduces the additional small parameter l c /l 1, simplified formulas can be obtained for the effective parameters. While sometimes qualified as a multiporous medium and although three distinct scales are actually involved, this composite is distinct from a double porosity microstructure (e.g. [START_REF] Auriault | Gas flow through a double-porosity porous media[END_REF], [START_REF] Royer | Time analysis of the three characteristic behaviours of dualporosity media. I: fluid flow and solute transport[END_REF], [START_REF] Boutin | On models of double porosity poroelastic media[END_REF]), as the distribution of the inclusions does not play the role of a porous network. Indeed, the impact of inclusions on fluid flow transfer is characterised by a purely geometrical parameter, with no reference to their fluid conductivity. Finally, we shall underline that all the above results are valid for elastic properties of both constituents in the same order of magnitude and for perfect interface bonding. A particular case of the model is derived and solved in [START_REF] Chen | Multiscale modelling and homogenisation of fibre-reinforced hydrogels for tissue engineering[END_REF], in which an incompressible istropic porous matrix reinforced with isotropic fibers is considered to model a construct for tissue engineering. Further work should include numerical simulations on specific geometries, so as to analyse the sensitivity of inclusion concentration.

A. Boundary value problem for u 1 s and u 1 c Let us multiply the system Eqs.( 27)-( 31) by a vectorial test function γ, and then let us integrate over Ω. We obtain the following variational formulation:

from which it appears that u 1 s is a linear vectorial function of ẽx ( u 0 ) and p 0 f , Eq.( 32), and that u 1 c is as a linear vectorial function of ẽx ( u 0 ), Eq.(33). Third-order tensors ω kh p i and ω kh c i are the specific solutions, u 1 s = ω kh p , u 1 c = ω kh c , to system Eqs.( 27)-(31), corresponding to    e x lm ( u 0 ) = 1 2 (δ lk δ mh + δ mk δ lh ) , p 0 f = 0.

As for vector π p , it is the specific solution for u 1 s when p 0 f = -1 and e x lm ( u 0 s ) = 0.

B. Equality of coupling tensors Ãeff and Geff

By taking u 1 s = ω kh p , u 1 c = ω kh c , γ = π p in Ω p , and γ = 0 in Ω c , in the variational formulation Eq.( 64 The left hand side of the above equation is positive, due to the positiveness of the local strain energy. Then, from the positiveness of the right hand side, it follows by Eq.( 61), that B eff > 0.

Figure 1 :

 1 Figure 1: Porous medium :(a) Macroscopic sample; (b) Periodic unit cell.

  Ω p c p ijlm e y lm ( u 1 s )e y ij ( γ) dΩ + Ω c c c ijlm e y lm ( u 1 c )e y ij ( γ) dΩ = -Ω p c p ijlm e y ij ( γ) dΩ e x lm ( u 0 ) -Ω c c c ijlm e y ij ( γ) dΩ e x lm ( u 0 ) + Ω p α p ij e y ij ( γ) dΩ p 0 f ,

  ), we getΩ p c p ijlm e y lm ( ω kh p )e y ij ( π p ) dΩ = -Ω p c p ijkh e y ij ( π p ) dΩ.Next, we take u 1 s = π p and γ = ω kh p in Eq.(64):Ω p c p ijlm e y lm ( π p )e y ij ( ω kh p ) dΩ = -Ω p α p ij e y ij ( ω kh p ) dΩ.The left hand sides of Eqs.(65) and (66) are equal. Then, from the equality of both right hand sides, it appears by Eqs.(55), (60), that Ãeff = Geff . C. Positiveness of the effective Biot bulk modulus B eff Considering u 1 s = π p , γ = π p , in Eq.(64), we get: Ω p c p ijlm e y lm ( π p )e y ij ( π p )dΩ = -Ω p α p ij e y ij ( π p )dΩ.

	(66)
	(65)

This is possible provided that the matrix pore size, l p be greatly smaller than the spacing between two inclusions, l.

Acknowledgements

This research is supported by CNRS (AAP "Osez l'Interdisciplinarité 2018", MoTiV Project). P. Recho and C. Verdier are members of LabeX Tec 21 (PIA : Grant No ANR-11-LABEX-0030).