Zeina Al Masry
email: zeina.almasry@univ-pau.fr

Sophie Mercier
email: sophie.mercier@univ-pau.fr

Ghislain Verdier
email: ghislain.verdier@univ-pau.fr

Approximate simulation techniques and distribution of an extended Gamma process

Keywords: Reliability, Degradation, Process with independent increments, Series representation, Inverse Laplace transform, Post-Widder formula, Weighted Gamma process

In reliability theory, many papers use a standard Gamma process to model the evolution of the cumulative deterioration of a system over time. When the variance-to-mean ratio of the system deterioration level varies over time, the standard Gamma process is not convenient any more because it provides a constant ratio. A way to overcome this restriction is to consider the extended version of a Gamma process proposed by Cinlar (1980). However, based on its technicality, the use of such a process for applicative purpose requires the preliminary development of technical tools for simulating its paths and for the numerical assessment of its distribution. This paper is devoted to these two points.

Introduction

Safety and dependability is a crucial issue in many industries, which has lead to the development of a huge literature devoted to the so-called reliability theory. In the oldest literature, the lifetimes of industrial systems or components were usually directly modeled through random variables, see e.g. [START_REF] Barlow | On a generalization of gamma processes[END_REF] for a pioneer work on the subject. Based on the development of on-line monitoring which allows the effective measurement of a system deterioration, numerous papers nowadays model the degradation in itself, which often is considered to be accumulating over time. This is done through the use of a non decreasing stochastic process (X t ) t≥0 (say), where the system lifetime τ typically corresponds to the hitting time of a given failure threshold L: τ = inf(t > 0 : X t > L).

A major quantity of interest for such a system is its reliability function R(t), which corresponds to the probability that, without any repair, the system is still in order at time t:

R(t) = P(τ > t) = P(X t ≤ L) = F X t (L) = L 0 f X t (s) ds,
where F X t and f X t stand for the cumulative distribution function (cdf) and probability density function (pdf) of X t , respectively. Another important indicator for prediction purpose is the system Mean Time Before Failure (MTBF), with

M T BF = E(τ ) = ∞ 0 P(τ > t) dt = ∞ 0 F X t (L) dt.
Most of other quantities of interest (Mean Residual Lifetime, cost/profit functions, . . . ) may also be expressed with respect to the pdf or cdf of X t , e.g. see Rausand and Høyland (2004). Before being able to use a stochastic process (X t ) t≥0 for deterioration modeling purpose and make prediction over the system future behavior, there hence is a crucial need to develop tools for computing its pointwise pdf and cdf at time t and/or for simulating trajectories of (X t ) t≥0 , which provides another way to compute the aimed quantities through Monte-Carlo simulations.

Since seminal papers [START_REF] Verdier | A gamma wear process[END_REF] and (C ¸inlar et al 1977), the most widely used process to model the phenomena of cumulative degradation is the (standard) Gamma process, see van [START_REF] Van Noortwijk | A survey of the application of gamma processes in maintenance[END_REF] for a comprehensive presentation and overview of applications of the Gamma process to reliability theory, including more than 100 references and both scientific and empirical justification of its use. However, a notable restriction of a standard Gamma process is that its variance-to-mean ratio is constant over time, which may be restrictive within an applicative context, see [START_REF] Guida | A time-discrete extended gamma process for time-dependent degradation phenomena[END_REF] for a real data set consisting of "sliding wear data of four metal alloy specimens" and where there is "empirical evidence that the variance-to-mean ratio is not a constant but varies with [time]". To overcome this drawback, we propose to use an Extended Gamma Process (EGP), which was introduced mostly simultaneously by C ¸inlar (1980) and [START_REF] Dykstra | A bayesian nonparametric approach to reliability[END_REF]. Note that EGPs are also called weighted Gamma processes in the literature [START_REF] Ishwaran | Computational methods for multiplicative intensity models using weighted gamma processes: proportional hazards, marked point processes, and panel count data[END_REF]). An EGP is a non decreasing process with independent increments, which can be constructed as a stochastic integral with respect to a standard Gamma process. If the EGP has been used for Bayesian modeling of the hazard function [START_REF] Dykstra | A bayesian nonparametric approach to reliability[END_REF][START_REF] Ishwaran | Computational methods for multiplicative intensity models using weighted gamma processes: proportional hazards, marked point processes, and panel count data[END_REF][START_REF] Laud | Monte Carlo methods for approximating a posterior hazard rate process[END_REF], up to our knowledge, it has not been much studied for cumulative degradation modeling, except in [START_REF] Guida | A time-discrete extended gamma process for time-dependent degradation phenomena[END_REF] in a simplified setting.

A standard Gamma process is characterized by a shape function and a constant scale parameter. For an EGP, the scale parameter may vary over time. This allows for more flexibility than its standard version, for modeling purpose. However, there is a cost and the use of an EGP presents some technical difficulties. Firstly, except for specific cases, the exact simulation of such a process is generally impossible. Secondly, there is no explicit formula for the probability distribution of an EGP (which is not Gamma). These technical difficulties have lead [START_REF] Guida | A time-discrete extended gamma process for time-dependent degradation phenomena[END_REF] to use a discrete version of an EGP. We here propose to deal with the original continuous time version.

The aim of the present paper is to develop technical tools necessary to the practical use of an EGP for cumulative degradation modeling (or for any other purpose). More precisely, we focus on the simulation of approximate sample paths and on the numerical assessment of both its pdf and cdf. Our contribution is threefold. Firstly, the series representations provided by [START_REF] Rosinski | Series representations of Lévy processes from the perspective of point processes[END_REF] for standard Gamma processes are extended to EGPs. Such representations are used to generate approximate sample paths. The quality of the approximation is studied and compared to the method proposed by [START_REF] Ishwaran | Computational methods for multiplicative intensity models using weighted gamma processes: proportional hazards, marked point processes, and panel count data[END_REF], which is based on an alternate approximate representation of an EGP. Secondly, following the results of [START_REF] Veillette | A technique for computing the pdfs and cdfs of nonnegative infinitely divisible random variables[END_REF] for Laplace transform inversion through Post-Widder formulas, explicit asymptotic expressions are provided for both pdf and cdf of an EGP. Thirdly, a discretization method is proposed, based on the approximation of a general EGP by a specific EGP with a piecewise constant scale function, which is easy to deal with. The discretization method allows both to simulate approximate sample paths of a general EGP and to compute approximations of its pdf and cdf. Convergence results are obtained for the approximate discretized EGP towards the initial general EGP and the quality of the approximation is studied. The method allows to compute the cdf of a general EGP at a known and adjustable precision. Up to our knowledge, a similar result was not available in the previous literature.

The paper is organized as follows. Section 2 gives an overview of the EGP. In Section 3, the approximate representation of [START_REF] Ishwaran | Computational methods for multiplicative intensity models using weighted gamma processes: proportional hazards, marked point processes, and panel count data[END_REF] is presented as well as four series representations of an EGP. Section 4 is devoted to the approximation of the pdf and cdf through Laplace transform inversion. Section 5 introduces and studies the discretization method. The different methods are illustrated through numerical experiments in Section 6. Conclusions are formulated in Section 7.

Definition of an EGP and first properties

Let a : R + → R + be a measurable, increasing and right-continuous function with a(0) = 0 and let b 0 > 0. Recall that a standard (non homogeneous) Gamma process Γ 0 (a(t), b 0 ) with a(.) as shape function and b 0 as (constant) scale parameter is a stochastic process with independent, non-negative and Gamma distributed increments. Its pdf at a specified time point t is given by

f t (x) = b a(t) 0 Γ (a(t))
x a(t)-1 exp(-b 0 x), ∀x ≥ 0, e.g. see [START_REF] Verdier | A gamma wear process[END_REF]. Now, let b : R * + → R * + be a measurable positive function such that, for all t > 0:

(0,t] da(s) b(s) < ∞. (1)
Following (C ¸inlar 1980;[START_REF] Dykstra | A bayesian nonparametric approach to reliability[END_REF], the process X = (X t ) t≥0 is said to be an EGP with shape function a(.) and scale function b(.) (written X ∼ Γ (a(t), b(t)) in the sequel) if it can be represented as a stochastic integral with respect to a standard Gamma process (Y t ) t≥0 ∼ Γ 0 (a(t), 1):

X t = (0,t] dY s b(s) , ∀t > 0 (2)
and X 0 = 0. If b(.) is constant and equal to b 0 , the EGP simply reduces to a standard Gamma process Γ 0 (a(t), b 0 ). An EGP can be proved to have independent increments and its distribution to be infinitely divisible (C ¸inlar 1980). Also, an explicit formula is available for the Laplace transform of an increment, with

L X t+h -X t (λ) := E e -λ(X t+h -X t ) = exp - (t,t+h] log 1 + λ b(s) da(s) ,
(3) for all t, λ ≥ 0 and h > 0.

Conversely, a stochastic process with independent increments and Laplace transform of an increment provided by (3) can be proved to be an EGP ∼ Γ (a(t), b(t)).

Finally, the mean and variance of an EGP are given by

E(X t ) = (0,t] da(s) b(s) and V(X t ) = (0,t] da(s) b(s) 2 . ( 4
)
Remark 1 Let X ∼ Γ (a(t), b(t)). So far, the shape function a was supposed to be rightcontinuous. Then, following C ¸inlar (1980), the function a can be written as a = a c + a d , where a c is the continuous part of a and a d its jump part. The EGP X is then the sum of two independent EGPs X c ∼ Γ (a c (t), b(t)) and X d ∼ Γ (a d (t), b(t)). This second EGP simply reduces, for each t ≥ 0, to a sum of independent Gamma variables, for which the distribution is known in full form (see Section 5). This second EGP is hence easy to deal with and can be omitted without loss of generality. Thus, the shape function a is assumed to be continuous in the sequel.

For a better understanding of the possible evolution of an EGP over time, we now look at its asymptotic behavior and we set

X ∞ = lim t→∞ X t (≤ ∞),
which exists, due to the non decreasingness property of an EGP.

Proposition 1 Let X = (X t ) t≥0 be an EGP with X ∼ Γ (a(t), b(t)). We then have:

(0,∞) log 1 + 1 b(s) da(s) ≤ E(X ∞ ) = (0,∞) da(s) b(s) ≤ ∞ (5)
and the following results:

1

. If E(X ∞ ) < ∞, then X ∞ is almost surely finite. 2. If (0,∞) log 1 + 1 b(s) da(s) = ∞, (6)
then X ∞ is almost surely infinite.

Proof Inequation (5) and the first point are clear. For the second one, let us note that

E(e -X ∞ ) = L X ∞ (1) = exp - (0,∞) log 1 + 1 b(s) da(s) .
Under assumption (6), we consequently have E(e -X ∞ ) = 0. This implies that e -X ∞ = 0 almost surely and the result.

Remark 2 In the case where

(0,∞) log 1 + 1 b(s) da(s) < ∞ = (0,∞) da(s) b(s) , (7)
the dichotomy P (X ∞ = ∞) ∈ {0, 1} is not valid any more and both P (X ∞ = ∞) > 0 and P (X ∞ < ∞) > 0 are possible, as is illustrated in the following example (third case).

Example 1 The different types of behavior are illustrated in Figure 1, 2 and 3 (left) where a few sample paths are plotted for three different EGPs (using the discretized rate method from Section 5):

X (1) ∼ Γ (t, (t + 1) 2 ), X (2) ∼ Γ t, 1 (t+1) 2
and

X (3) ∼ Γ 1 2 1 -1 (t+1) 2 , 1 exp(1+t)-1 . The process X (1) (resp. X (2)
) corresponds to the first (resp. second) case of Proposition 1 whereas X (3) corresponds to case (7). As expected, all trajectories of Figure 1 are stabilizing whereas all trajectories of Figure 2 are exploding. In Figure 3, some trajectories are exploding and others are stabilizing. For a better insight into their behavior, the variance-to-mean ratios of the X (i) 's are also plotted in the same figures (right). These figures illustrate the flexibility of EGPs for modeling purpose.

In all the following, X = (X t ) t≥0 stands for an EGP ∼ Γ (a(t), b(t)) and Y = (Y t ) t≥0 stands for a standard Gamma process Γ 0 (a(t), 1), without any further notification. We recall that the shape function a(t) is assumed to be continuous.

Simulation techniques of an EGP

The goal of this section is to present methods for generating approximate sample paths of an EGP on a given compact set [0, T ]. The section is divided into two parts. The first part quickly recalls the method proposed by [START_REF] Ishwaran | Computational methods for multiplicative intensity models using weighted gamma processes: proportional hazards, marked point processes, and panel count data[END_REF]. In the second part, the series representations of a standard Gamma process provided in [START_REF] Rosinski | Series representations of Lévy processes from the perspective of point processes[END_REF] are extended to the case of an EGP. Algorithms based on these series representations are further provided in Section 6 for approximate sample paths generation.

Ishwaran and James's method

The method proposed by [START_REF] Ishwaran | Computational methods for multiplicative intensity models using weighted gamma processes: proportional hazards, marked point processes, and panel count data[END_REF] roughly boils down to some discretization of the stochastic integral (2) for t ∈ [0, T ], where T > 0. More specifically, for each t ∈ [0; T ] and each K * ∈ N, they consider:

G (K * ) t = K * k=1 1 [0,t] (V k ) 1 b(V k ) W (K * ) k , (8)
where

(V k ) k=1,...K * are i.i.d random variables (r.v.s) with distribution da(v) a(T ) 1 [0,T ] (v) and (W (K * ) k ) k=1,...K * are i.i.d r.v.s with Gamma distribution Γ 0 ( a(T )
K * , 1), independent from the V k 's. [START_REF] Ishwaran | Computational methods for multiplicative intensity models using weighted gamma processes: proportional hazards, marked point processes, and panel count data[END_REF]

proved that G (K * ) t t≥0 -→ (X t ) t≥0 weakly when K * → ∞
in the Skorokhod space D (0, ∞) of right-continuous functions with left-side limits. (To be V(X

0 0.2 0.4 0.6 0.8 0.5 0.6 0.7 0.8 0.9 1 V(X (1) t )/E(X (1) t ) t 0 2 4 6 8 10 0 0.5 1 1.5 2 X (1) t t Fig. 1 Sample paths (left) and ratio V(X (1) t ) E(X (1) t ) (right) as a function of time 0 2 4 6 8 10 0 200 400 600 800 X (2) t t 0 0.2 0.4 0.6 0.8 1 1.5 2 2.5 3 V(X (2) t )/E(X (2) t ) t Fig. 2 Sample paths (left) and ratio V(X (2) t ) E(X (2) t ) (right)
(3) t )/E(X (3) t ) t Fig. 3 Sample paths (left) and ratio V(X (3) t ) E(X (3) t )
(right) as a function of time more specific, they show the weak convergence of the random measure with corresponding cumulative process defined by (8) towards a weighted gamma measure with shape measure a(dt) and scale function b(t). The weak convergence of the cumulative process (8) towards an EGP is then a consequence of (Daley and Vere-Jones 2007, Lemma 11.1.XI)).

One can easily check that

E G (K * ) t = E[X t ], V G (K * ) t = V[X t ] + 1 K * a(T )V[X t ] -(E[X t ]) 2
so that the approximation of X t by G (K * ) t is unbiased. Also, the rate of convergence of the variance is O 1 K * .

Series representations

The following results extend those from [START_REF] Rosinski | Series representations of Lévy processes from the perspective of point processes[END_REF], Section 6) devoted to standard Gamma processes, see [START_REF] Rosinski | Series representations of Lévy processes from the perspective of point processes[END_REF] for more references and historical remarks.

Proposition 2 Let T > 0 and let (U n ) n≥1 be the points of a homogeneous Poisson process M with parameter a(T ). Let also

(V n ) n≥1 be a sequence of i.i.d. r.v.s with distribution H(dv) = da(v) a(T ) 1 [0,T ] (v)
, independent of M and let (W n ) n≥1 be a sequence of i.i.d r.v.s with distribution P W (dw), independent of M and of the V n 's. Then we have the four following series representations of an EGP, where D = means "is identically distributed as":

(1) Bondesson's series representation (Bond):

X t D = n≥1 1 b(V n ) exp(-U n ) W n 1 [0,t] (V n ), for 0 ≤ t ≤ T, (9)
where {W n } n≥1 is a sequence of i.i.d exponential r.v.s with mean 1, (2) Rejection's series representation (Rej):

X t D = n≥1 1 b(V n ) 1 exp(U n ) -1 1 E (U n , W n ) 1 [0,t] (V n ), for 0 ≤ t ≤ T, (10)
where

E = (u, w) ∈ R + × [0, 1] | exp(u) exp(u)-1 exp(-(exp(u) -1) -1 ) ≥ w and where {W n } n≥1 is a sequence of i.i.d uniform r.v.s on [0, 1],
(3) Thinning's series representation (Thin):

X t D = n≥1 1 b(V n ) W n 1 (W n U n ≤1) 1 [0,t] (V n ), for 0 ≤ t ≤ T, (11)
where {W n } n≥1 is a sequence of i.i.d exponential r.v.s with mean 1, (4) Inversion of Lévy measure's series representation (ILM):

X t D = n≥1 1 b(V n ) E -1 1 (U n )1 [0,t] (V n ), for 0 ≤ t ≤ T, (12
)

with E 1 (x) = ∞ x u -1 exp(-u)du, the exponential integral function.
Proof The four series representations can be written in the same manner:

X t = n≥1 f t (U n , V n , W n ) = n≥1 1 b(V n ) H(U n , W n )1 [0,t] (V n ), (13)
with f t (u, v, w) = 1 b(v) H(u, w)1 [0,t] (v)
and an appropriate choice for H : R 2 + → R + . As already mentioned in Section 2, (X t ) 0≤t≤T is an EGP if its increments are independent and the Laplace transform of an increment can be written as in (3). Let us first note that (U n , V n , W n ) n≥1 can be seen as the points of a Poisson random measure N on R 3 + with intensity

µ(du, dv, dw) = du × da(v)1 [0,T ] (v) × P W (dw)
(see for example (C ¸ınlar 2011, Corollary 3.5, p. 265)).

Considering

t 1 < • • • < t n , we then have X t i+1 -X t i = N (f t i+1 -f t i ) for all 1 ≤ i ≤ n -1, where f t i+1 -f t i takes range in R + × (t i , t i+1 ] × R + .
As the ranges of the f t i+1 -f t i are non overlapping, the increments X t i+1 -X t i are independent. Now, based on the Laplace functional theorem for a Poisson random measure (C ¸ınlar 2011, Theorem 2.9, p. 252) for the third line, we have for all t ∈ [0, T ] and all h, λ > 0:

L X t+h -X t (λ) = E[exp(-λ(X t+h -X t ))] = E[exp(-N (λ(f t+h -f t )))] = exp - R 3 1 -e -(λ(f t+h -f t )) dµ = exp - (t,t+h] Q(λ, v)da(v) (14)
where

Q(λ, v) = R 2 1 -exp -λ 1 b(v) H(u, w) du × P W (dw).
From ( 3) and ( 14), it then suffices to check that Q(λ, v) = log 1 + λ b(v) for each of the four series representations.

1) Bondesson:

In this case, H(u, w) = exp(-u)w and P W (dw) = exp(-w)dw. Then

Q(λ, v) = ∞ 0 ∞ 0 1 -exp - λ b(v) exp(-u)w exp(-w)dw du = ∞ 0 1 - b(v) b(v) + λ exp(-u) du = ∞ 0 λ exp(-u) b(v) + λ exp(-u) du = log 1 + λ b(v)
.

2) Rejection:

Here, H(u, w)

= 1 exp(u)-1 1 ( exp(u) exp(u)-1 exp(- 1 exp(u)-1 )≥w) and P W (dw) = dw. Accord- ingly, Q(λ, v) = ∞ 0 exp(u) exp(u)-1 exp(- 1 exp(u)-1 ) 0 1 -exp - λ b(v) 1 exp(u) -1 dw du = ∞ 0 1 -exp - λ b(v) 1 exp(u) -1 exp(u) exp(u) -1 exp - 1 exp(u) -1 du = ∞ 0 1 -exp - λ b(v) z exp(-z) z dz (15)
setting z = 1 exp(u)-1 in the last line. (15) is a Frullani's integral, which can classically be computed by differentiating with respect of λ:

Setting o(λ, v, z) = 1 -exp -λ b(v) z exp(-z) z
, we get by dominated convergence that

∂Q(λ, v) ∂λ = ∞ 0 ∂o(λ, v, z) ∂λ dz = 1 b(v) + λ ,
which provides the result.

3) Thinning: H(u, w) = w1 (uw≤1) and P W (dw) = exp(-w)dw.

Q(λ, v) = ∞ 0 1 w 0 1 -exp - λ b(v) w exp(-w)du dw = ∞ 0 1 -exp - λ b(v) w exp(-w) w dw.
We recognize (15), which allows to conclude.

4) Inversion of Lévy measure:

H(u, w) = E -1 1 (u) is independent on w. Then Q(λ, v) = ∞ 0 1 -exp - λ b(v) E -1 1 (u) du = ∞ 0 1 -exp - λ b(v) z exp(-z) z dz setting z = E -1 1 (u).
We recognize (15) again, which allows to conclude.

Approximate simulation of (X t ) 0≤t≤T is done by truncating the infinite series. As is classically done in the Lévy case, we retain only the points of the Poisson process M which belong to a compact set [0, B], where B > 0:

X (B) t = n≥1 1 b(V n ) H(U n , W n )1 [0,t] (V n )1 [0,B] (U n ) for all t ∈ [0, T ],
where we use the notations of (13). This allows to better control the truncation error than simply retaining a fixed number of terms in the series, see e.g. [START_REF] Imai | Numerical inverse Lévy measure for infinite shot noise series representation[END_REF] for further details in the Lévy case. We also set

X(B) t = X t -X (B) t = n≥1 1 b(V n ) H(U n , W n )1 [0,t] (V n )1 (B,∞) (U n ) (16)
for all t ∈ [0, T ], to be the remainder of the truncated series. By construction, both X(B)

t t∈[0,T ]
and

X (B) t t∈[0,T ]
have independent increments. The moments of X(B) t are given in the following proposition.

Proposition 3 Recall that X(B)

t
is the remainder of the truncated series (16). We have

(1) Bondesson's series representation:

X(B) t ∼ Γ a(t), b(t) exp(-B)
;

E[ X(B) t ] = exp(-B)E[X t ]; V[ X(B) t ] = exp(-2B)V[X t ],
(2) Rejection's series representation:

E[ X(B) t ] = E[X t ] 1 -exp - 1 exp(B) -1 ; V[ X(B) t ] = V [X t ] 1 -1 + 1 exp(B) -1 exp - 1 exp(B) -1 ,
(3) Thinning's series representation:

E[ X(B) t ] = E[X t ] 1 -B + B exp - 1 B ; V[ X(B) t ] = V[X t ] 1 + exp - 1 B + 2B exp - 1 B -2B ,
(4) Inversion of Lévy measure's series representation:

E[ X(B) t ] = E[X t ](1 -E 1 (B)); V[ X(B) t ] = V[X t ](1 -exp(-E 1 (B)) -E 1 (B) exp(-E 1 (B))).
Proof Starting from (16), the Laplace transform of X(B) t is computed using the Laplace functional theorem for Poisson random measures. In the Bondesson case, we recognize the Laplace transform of an EGP ∼ Γ a(t), b(t) exp(-B) . For the three other representations, the mean and variance are obtained through the first and second order derivatives of the Laplace transform at λ = 0.

Proposition 3 allows to determine the minimal value of B which ensures a given relative precision on both mean and variance of X t for all t ∈ [0, T ] (analytically for Bondesson's series representation and numerically for the three other procedures). Asymptotic equivalents of the series remainders are provided in Table 1 when B tends to infinity, which are obtained through Taylor expansions.

Method E X(B) t /E(Xt) V X(B) t /V(Xt) Bondesson exp(-B) exp(-2B) Rejection exp(-B) exp(-2B) 2 Thinning 1 2B 1 6B 2 Inversion of Lévy measure E 1 (B) E 1 (B) 2 2 Table 1 Asymptotic equivalents of E X(B) t /E(Xt) and V X(B) t /V(Xt) when B → ∞
From Table 1, we can note that the Inversion of Lévy measure's method converges very fast (E 1 (B) = o (exp(-B)) when B → ∞). However, the computing time of the Inversion of Lévy measure's method may be long because of the numerical computation of the inverse of the exponential integral function, for which no explicit formula is available. As far as the Thinning's method is concerned, it clearly converges more slowly than the others. Thus, Rejection's and Bondesson's methods seem to be the most suitable ones among the four series representations to approximate the paths of an EGP.

Post-Widder formulas

As already noted in Section 2, the distribution of an EGP is infinitely divisible and its Laplace transform is available in full form (see (3)). This allows to use the method from [START_REF] Veillette | A technique for computing the pdfs and cdfs of nonnegative infinitely divisible random variables[END_REF] to compute the pdf and cdf of an EGP, by inverting its Laplace transform through the Post-Widder method. We use the following formulas from (Masol and Teugels 2010, Section 4) for the pdf and cdf (see also [START_REF] Veillette | A technique for computing the pdfs and cdfs of nonnegative infinitely divisible random variables[END_REF] for the cdf): for all x ≥ 0,

F X t (x) = lim H→∞ H k=0 - H x k 1 k! L (k) X t H x , (17)
f X t (x) = lim H→∞ - 1 H! - H x H+1 L (H) X t H x , (18)
where L

(k)

X t denotes the k-th derivative of the Laplace transform L X t for k ≥ 1 and L (0)

X t = L X t .
Let φ t be the Laplace exponent of X t , with

φ t (λ) = -ln (L X t (λ)) , ∀λ ≥ 0.
Following Veillette and Taqqu (2011), Leibniz's formula provides:

L (k+1) X t (λ) = - k j=0 k j L (j) X t (λ)φ (k+1-j) t (λ), ∀λ, t, k ≥ 0, (19)
with φ (j) t the j-th derivative of φ t . An expression of φ (j) t is given in the following Lemma in the case of an EGP.

Lemma 1 Let λ ≥ 0, j ≥ 1 and let X ∼ Γ (a(t), b(t)). Then φ (j) t (λ) = (-1) j-1 (j -1)! m (j) t (λ) where m (j) t (λ) = t 0 da(s) (b(s) + λ) j .
Proof Based on (3) and (15) and setting u = z/b(s) in the third line we have:

φ t (λ) = (0,t] log 1 + λ b(s) da(s) = (0,t] ∞ 0 1 -exp - λ b(s) z exp(-z) z dz da(s) = ∞ 0 1 -e -λu µ t (du)
where [START_REF] Veillette | A technique for computing the pdfs and cdfs of nonnegative infinitely divisible random variables[END_REF], φ

µ t (du) = 1 u (0,t] e -b(s)u da(s) du is a measure on R + such that R + (1 ∧ x)µ t (dx) < ∞. From
(j) t (λ) is directly calculated by φ (j) t (λ) = (-1) j+1 ∞ 0 x j e -λx µ t (dx) (20)
for all j ≥ 1. This provides:

φ (j) t (λ) = (-1) j+1 t 0 ∞ 0 x j-1 e -(b(s)+λ)x dx da(s) = (-1) j+1 t 0 Γ (j) (b(s) + λ) j da(s),
renormalizing the pdf of Γ 0 (j, b(s) + λ) to get the last line.

Based on (19) and Lemma 1, L (k) X t (λ) can be recursively computed through

L (k+1) X t (λ) = -k! k j=0 L (j) X t (λ) (-1) k-j j! m (k+1-j) t (λ) .
This allows to compute an approximation of the pdf and cdf of an EGP through (17 -18), taking H large enough.

Discretized rate function method

In case of a piecewise constant scale function b(.), the process (X t ) t≥0 can easily be constructed from standard Gamma processes. The simulation of its paths is hence immediate (see e.g. van Noortwijk ( 2009)). Also, the random variable X t simply is the sum of standard Gamma variables, and different tools are available in the literature to compute both its pdf and cdf (see [START_REF] Nadarajah | A review of results on sums of random variables[END_REF] for a review). Based on this, the aim of this section is to propose an approximation of an EGP with a general scale function by another EGP with a piecewise constant scale function.

Construction of the approximate process X (ε)

Let T > 0. The following assumption is considered: b(.) is continuous on (0, T ] and ∃m > 0 such that ∀t ∈ (0, T ], b(t) ≥ m.

Let ε > 0. The following piecewise constant approximation of 1 b(.) on (0, T ], denoted by

1 b (ε) (.) , is considered: ∀t ∈ (0, T ], 1 b (ε) (t) = n(ε) i=0 1 b i 1 [l i ,l i+1 ) (t) (22)
where n(ε) is such that l n(ε) ≤ T < l n(ε)+1 , l 0 = 0, and the l i 's, for i = 1, . . . , n(ε) + 1, are recursively defined by:

l i+1 = sup l ∈ (l i , T ] : ∀l ∈ [l i , l], 1 b(l i ) - 1 b(l ) < ε . ( 23
)
Note that assumption (21) insures that n(ε) is finite. The constants b i , i = 0, . . . , n(ε) are next defined by:

∀i = 0, . . . , n(ε) -1, 1 b i = 1 l i+1 -l i l i+1 l i 1 b(s) ds, (24)
and

1 b n(ε) = 1 T -l n(ε) T l n(ε) 1 b(s) ds. (25)
Now that b (ε) (.) is constructed, we set

X (ε) t = (0,t] dY s b (ε) (s) , ∀t ∈ [0, T ], (26)
where (Y t ) t≥0 is the same standard Gamma process Γ 0 (a(t), 1) as in Equation (2) defining

X. Then X (ε) t t∈[0,T ]
is the restriction to [0, T ] of an EGP ∼ Γ a(t), b (ε) (t) with the same shape function a(.) as X and with b (ε) (.) as scale function.

Remark 3 Based of the definition (2) of an EGP, we prefer to construct an approximation of 1 b(.) instead of b (.). Also, by construction, for all t in (0, T ],

1 b(t) - 1 b (ε) (t) ≤ ε. (27)
5.2 Quality of the approximation of X by X (ε) We first look at the moments of the residual part in the approximation of X t by X

(ε) t , for t ∈ [0, T ]. Proposition 4 Let X(ε) t = X t -X (ε) t = (0,t] 1 b(s) - 1 b (ε) (s) dY s , ∀t ∈ [0, T ]. (28)
Then:

E[| X(ε) t |] ≤ ε a(t), (29)
V[ X(ε) t ] ≤ ε 2 a(t) (30)
for all t ∈ [0, T ]. Proof Setting b(ε) (s) = 1/b(s) -1/b (ε) (s) -1 if b(s) = b (ε) (s) and b (ε) (s) = ∞ otherwise, we have X(ε) t = t 0 1 b(ε) (s) dY s where 1 | b(ε) (s)| ≤ , based on (27). We easily derive that | X(ε) t | ≤ t 0 1 b(ε) (s) dY s ≤ ε Y t (31)
from where we derive (29). Also, one can check that (4) is still valid for X(ε) t with b(t) substituted by b(ε) , so that

V[ X(ε) t ] = (0,t] 1 b(ε) (s) 2 da(s) ≤ ε 2 a(t).
We now provide some convergence results.

Proposition 5 Let X ∼ Γ (a(t), b(t)), (ε n ) n≥1 be a sequence of positive real numbers and let (X (ε n ) ) n≥1 be a sequence of EGPs with shape function a(.) and scale function b (ε n ) (.). Suppose that assumption (21) is satisfied.

1) If n≥1 √ ε n < ∞, then sup t∈[0,T ] |X (ε n ) t -X t | ----→ n→∞ 0 almost surely. 2) If ε n ----→ n→∞ 0, then sup t∈[0,T ] |X (ε n ) t -X t | ----→ n→∞ 0 in quadratic mean.
Proof 1) According to (31), we have:

P sup t∈[0,T ] |X (ε n ) t -X t | > √ ε n ≤ P ε n sup t∈[0,T ] Y t > √ ε n = P Y T > 1 √ ε n ≤ E (Y T ) √ ε n ,
based on Markov inequality for the last line. We derive that

n≥1 P sup t∈[0,T ] |X (ε n ) t -X t | > √ ε n < ∞,
which concludes the first point, based on e.g. (C ¸ınlar 2011, Proposition 2.7, p. 98). 2) In the same way, we have:

E sup t∈[0,T ] |X (ε n ) t -X t | 2 ≤ E ε n sup t∈[0,T ] Y t 2 = ε 2 n E Y 2 T ,
which allows to conclude.

Remark 4 Note that a consequence of point 2 is that, under the same assumptions, ] weakly in the Skorokhod space D([0, T ]).

X (ε) t t∈[0,T ] -→ (X t ) t∈[0,T

Approximation of the pdf and cdf of X

(ε) t

(and of X t )

As previously mentioned, X (ε) is an EGP with a piecewise constant scale function. For each t, X (ε) t can be written as a sum of independent Gamma distributed random variables with different scale parameters. Several methods are available in the literature to evaluate the pdf of such a sum, see [START_REF] Nadarajah | A review of results on sums of random variables[END_REF]. In [START_REF] Guida | A time-discrete extended gamma process for time-dependent degradation phenomena[END_REF], the authors use the method from [START_REF] Moschopoulos | The distribution of the sum of independent gamma random variables[END_REF] to evaluate the pdf of their discrete time EGP. After performing some numerical comparisons (not provided here), we finally chose to use the method from (Peppas 2011, Sub. 4.2), which provides an expression of both pdf and cdf of X (ε) t in terms of an infinite integral. For all t in [0, T ] and all x ≥ 0, this writes:

f X (ε) t (x) = lim K→∞ f (K) X (ε) t (x) and F X (ε) t (x) = lim K→∞ F (K) X (ε) t (x), where f (K) X (ε) t (x) = 1 π K 0 cos( P p=0 α p arctan(u/b p ) -xu) P p=0 (1 + (u/b p ) 2 ) α p /2 du, (32)
F (K) X (ε) t (x) = 1 2 - 1 π K 0 sin( P p=0 α p arctan(u/b p ) -xu) u P p=0 (1 + (u/b p ) 2 ) α p /2 du ( 33
)
with P = P ( , t) the single integer such that l P < t ≤ l P +1 , α p = a(l p+1 ) -a(l p ) for p = 0, . . . , P -1 and α P = a(t) -a(l P ).

Taking K large enough, (32 -33) provide approximations for the pdf and cdf of X (ε) t , and consequently of X t . Using that | sin(.)| ≤ 1 and that (1 + (u/b p ) 2 ) α p /2 ≥ (u/b p ) α p , we easily get an upper bound for the approximation error, with:

E (K) X (ε) t (x) = |F X (ε) t (x) -F (K) X (ε) t (x)| ≤ 1 π ∞ K 1 u P p=0 (u/b p ) α p du ≤ P p=0 b α p p π ∞ K 1 u a(t)+1 du ≤ P p=0 b α p p π a(t)K a(t) (34)
for all K ∈ N, due to P p=0 α p = a(t) for the second-to-last line.

Bounds for the cdf of X t

Consider two other piecewise constant approximations of 1 b(t) defined by

1 b (ε,-) (t) = n(ε) i=0 1 b - i 1 [l i ,l i+1 ) (t) and 1 b (ε,+) (t) = n(ε) i=0 1 b + i 1 [l i ,l i+1 ) (t)
for all t ∈ (0, T ], where

1 b - i = inf 1 b(l) , l ∈ [l i , min{l i+1 , T }) and 1 b + i = sup 1 b(l) , l ∈ [l i , min{l i+1 , T }) .
The following bounds are obtained:

1 b (ε,-) (t) ≤ 1 b(t) ≤ 1 b (ε,+) (t)
, for all t ≤ T.

(

) 35
The induced EGPs are denoted by X (ε,-) and X (ε,+) , respectively. They satisfy

X (ε,-) t ≤ X t ≤ X (ε,+) t
, for all t ≤ T so that

F X (ε,+) t (x) ≤ F X t (x) ≤ F X (ε,-) t (x), ∀t ≤ T, ∀x ≥ 0. ( 36
)
Theorem 1 Using the notations of Subsection 5.3 substituting X (ε) by X (ε,-) and by X (ε,+) , we set:

m t (x, ε, K) = F (K) X (ε,+) t (x) -E (K) X (ε,+) t (x), M t (x, ε, K) = F (K) X (ε,-) t (x) + E (K) X (ε,-) t (x)
for all K ∈ N * and all ε > 0. Then,

m t (x, ε, K) ≤ F X t (x) ≤ M t (x, ε, K) (37)
and

F X t (x) - m t (x, ε, K) + M t (x, ε, K) 2 ≤ M t (x, ε, K) -m t (x, ε, K) 2 . ( 38
)
Proof The bounds (37) are easily obtained starting from Equation (36) and using that

F X (ε,-) t (x) ≤ F (K) X (ε,-) t (x) + E (K) X (ε,-) t (x) = M t (x, ε, K)
for the upper bound, and similar arguments for the lower one. Inequality (38) is a direct consequence of (37).

The previous theorem provides computable bounds for the cdf F X t of an EGP. These bounds can be made as tight as necessary, taking ε small enough and K large enough. This method hence provides a way to numerically assess the cdf F X t at a known precision. This is used in Section 6 to get numerical reference results. Remark also that the two piecewise constant approximations b (ε,-) (.) and b (ε,+) (.) satisfy ( 27). The theoretical results of Subsection 5.2 are then still valid for X (ε,-) and X (ε,+) . Finally, remark that for any fixed s > 0, the increment process (X t+s -X s ) t≥0 still is an EGP ∼ Γ (a(t + s) -a(s), b(t + s)) (based on its independent increments and on its pointwise Laplace transform provided by (3)). Thus, all the results of the paper can be used to compute both pdf and cdf of an EGP increment, or to simulate it. (For simulation purpose, one can also only retain the terms for which V n ∈ (s, s + t] in the series representations from Proposition 2).

When assumption (21) is not satisfied

Assumption (21) allows to construct the piecewise constant function b (ε) (.) specified in Subsection 5.1. When this assumption is not satisfied, alternatives for the construction of b (ε) (.) have to be considered.

• An interesting case is when function b(.) satisfies lim

t→0 + b(t) = 0. If the following as- sumption is checked, ∃η > 0 such as η 0 ds b(s) < ∞,
then the piecewise constant approximation can be constructed from l 0 = 0, l 1 such that:

l 1 0 ds b(s) ≤ ε and l 1 0 da(s) b(s) ≤ ε,
and the l i 's, for i = 2, . . . , n(ε) + 1, recursively defined according to (23). Note that in this case, inequality ( 27) is not satisfied any more. • If function b(.) is not continuous on (0, T ] but only piecewise continuous, the procedure described in Subsection 5.1 can be applied on each continuous section and inequality ( 27) is still valid.

6 Numerical Experiments

Summary and Algorithms

To sum up the above sections, three main methods are considered to simulate an EGP:

-Method 1: Discretization of stochastic integral (DSI) from Subsection 3.1, -Method 2: Series representations from Subsection 3.2 : ILM, Bond, Thin, Rej, -Method 3 : Discretized rate function (DRF) from Subsection 5.1.

Furthermore, two approaches to approximate the pdf and cdf of an EGP are proposed:

-Approach 1: Post-Widder formula from Section 4, -Approach 2: DRF method + approximation of the cdf and pdf through the infinite integral formulas from Subsection 5.3. This method is referred to as DRFI in the sequel.

For the DRFI method, the cdf of X t is approximated by

F (K) X (ε) t
as provided by (33), where we recall that X (ε) t is defined by (26). There hence are two parameters to define (ε and K).

All these approaches are investigated in this section. We set X ∼ Γ (a(t), b(t)) to be an EGP, where a(t) is assumed to be one-to-one. We first state the algorithms resulting from the three methods previously mentioned to simulate N approximate sample paths of (X t ) t∈[0,T ] . Note that for Method 2, only Bondesson's series representation algorithm is given.

Algorithm 1 Method 1: Ishwaran and James's approach 1. Fix K * . For i in 1 : N , repeat the following steps : 2. Simulate K * realizations Z 1 , . . . , Z K * according to the uniform distribution on [0, 1]. 3. Compute V j = a -1 (Z j a(T )) for j = 1, . . . , K * (inverse transform sampling).

Simulate

K * realizations W (K * ) 1 , . . . , W (K * ) K * according to the Gamma distribution Γ 0 a(T ) K * , 1 . 5. Compute G (K * ) t,i = K * j=1 1 [0,t] (V j ) 1 b(V j ) W (K * ) j .
Algorithm 2 Method 2: Bondesson's series representation

[0, 1]. 6. Compute V j = a -1 (Z j a(T )) for j = 1, ..., N B . 7. Compute X (B) t,i = N B j=1 1 b(V j ) exp(-U j )W j 1 [0,t] (V j ).
Algorithm 3 Method 3: Discretized rate method 1. Fix ε. 2. Construct the sequence (l j ) j=1,...,n(ε)+1 using (23). 3. Approximate the rate function using (22). For i in 1 : N , repeat the following steps : 4. Simulate Y j according to the Gamma distribution Γ 0 a(l j+1 ) -a(l j ), b (ε) (l j ) for j = 0, . . . , n(ε) -1 and Y n(ε) according to the distribution Γ 0 a(T ) -a(n(ε)), b (ε) (n(ε)) .

5. Compute X (ε)

t,i = n(ε) j=1 Y j 1 {t≤l j+1 } .

Comparison of the simulation methods

The different approaches for the sample paths generation of an EGP are compared. For each method, N = 10 5 approximate sample paths are generated and the relative errors between the theoretical and empirical mean, variance and Laplace transform are computed. Let us recall that the Laplace transform of X t fully characterizes its distribution. In order to scan a good part of this distribution, we introduce λ j = L -1 X t (0.01 * j), j = 1, . . . , 100 and we define the relative error on the Laplace transform by:

E L (t) = 1 100 100 j=1 |L X t (λ j ) -LX (ε) t (λ j )| L X t (λ j ) , (39)
where

LX (ε) t (λ j ) is the empirical Laplace transform of X (ε) t
at point λ j . The previous simulation is repeated 500 times and for each method, we compute the mean and the standard deviations of the tree relative errors (mean, variance and Laplace transform) based on these 500 sets of N = 10 5 approximate sample paths.

As a first step, the four series representations of Method 2 are compared and the results provided in Table 2. As expected from Table 1, the thinning method seems less accurate than the others. Moreover, the computing time for the inverse Lévy measure method is thrice the computing time of Bondesson and Rejection methods for a similar precision. Thus, only Bondesson and Rejection approaches are maintained to continue the comparison. Tables 3, 4 and 5 provide the results for the DSI method from Subsection 3.1, the two selected series representations (Bondesson and Rejection) and for the proposed approach of the discretized rate function for different shape and rate functions (see the legends of the tables). The parameters of the four methods were adjusted to have similar computing times. The results for the discretized rate function are quite similar to the series representations methods, while the DSI method is sometimes less accurate for the variance or for the Laplace transform (see Table 4). The discretized rate method hence seems to behave as well as Bondesson's and Rejection methods for simulation purpose, whereas the DSI method seems a little below.

Comparison of the approximation methods for the cdf and pdf

We here compare the quality of the numerical assessment of the cdf/pdf of an EGP through both Post-Widder's formula and the DRFI method. We begin with the cdf, for which we recall that we are able to get reference results, with an exact control of the error (see Subsection 5.4). An EGP X ∼ Γ (a(t), b(t)) is considered with a(t) = t and b(t) = (t + 1) 2 . Reference results are first computed for F X t (x) with x ∈ [0.01 : 0.1 : 4] and t = 2, up to a precision lower than 3.10 -6 , obtained for K = 800 and ε = 10 -6 . The results are provided in Table 6 for the mean absolute error on F X t (x) with respect to the reference results for x ∈ [0.01 : 0.1 : 4] and both Post-Widder's and DRFI methods. For H 100 in Post-Widder's method and ε = 10 -2 , K = 12 for the DRFI approximation, both methods provide similar results with similar computation times. However, the DRFI method easily yields more accurate results by decreasing ε and increasing K whereas the highest precision by Post-Widder's method is obtained for H 100 in the present example, with a clearly lower precision for H = 200. This problem of poor convergence when H → ∞ for Post-Widder's method has already been observed in [START_REF] Masol | Numerical accuracy of real inversion formulas for the Laplace transform[END_REF]. 15 DRFI (ε = 10 -2 , K = 800) 2.10 -5 15 DRFI (ε = 10 -3 , K = 800) < precision 30

Keeping the same a(t) and b(t), three approximations of the cdf of X 10 are next plotted in Figure 4: a non-parametric estimation obtained from a sample of size 10 5 generated by Rejection method, Post-Widder's estimation (H = 100) and the DRFI approximation (ε = 0.001 and K = 100). We observe the good superposition of the DRFI approximation and of the non-parametric estimation. These results are in concordance with the convergence results obtained in Subsection 5.2, which imply the convergence of the cdf F (ε) X 10 at each continuity point of F X 10 . In the zoomed part of Figure 4 (right), one can observe that Post-Widder's method seems a little less accurate than the proposed DRFI approximation.

Figure 5 represents the corresponding approximations for the pdf of X 10 . Even if there is no theoretical convergence result for the DRFI approximation, the closeness of the new approach with the non-parametric one gives us good reason to believe that our approximation is efficient, keeping in mind, however, that the non-parametric estimation is obtained from approximate realizations of X t . Note that, here again, the zoomed part of Figure 5 (right) shows that Post-Widder's method seems a little less accurate than the DRFI approximation.

Conclusion

Different tools based on some discretization of the rate function of an EGP have been presented here, firstly, for the approximate simulation of an EGP, and secondly, for the numerical assessment of the cdf/pdf of an EGP. Three discretization schemes have been proposed: one provides the best approximation among the three (b (ε) (t)) and the other two (b (ε,+) (t), b (ε,-) (t)) provide bounds for X t and for its cdf F X t .

As far as simulation procedures are concerned, it seems that our approximate simulation scheme behaves as well as two of the most usual ones (Bondesson's and rejection methods), previously developed in the context of subordinators. Also, our simulation scheme seems to behave better than the one proposed by [START_REF] Ishwaran | Computational methods for multiplicative intensity models using weighted gamma processes: proportional hazards, marked point processes, and panel count data[END_REF], specifically developed for the simulation of an EGP.

As for the numerical assessment of the cdf of an EGP, we could not find in the literature any available procedure with a similar control on the precision (except from [START_REF] Veillette | A technique for computing the pdfs and cdfs of nonnegative infinitely divisible random variables[END_REF], which provides a refinement of the Post-Widder's method as well as some asymptotic bounds for the error when H → ∞). Beyond this control on the precision, the results provided by the proposed method have been compared to those obtained by Post-Widder's formula, at the advantage of our method. Even if we do not have any similar control on the precision for the pdf, a good behavior of the proposed approximation has been numerically observed, here again at its advantage when compared to Post-Widder.

To sum up, the discretized rate function method seems to behave well, both for simulating approximate paths and for the numerical assessment of the cdf/pdf of an EGP.

Before being able to use the proposed tools in an applied situation (for e.g. deterioration modelling in an industrial reliability context), another important issue that requires further study concerns the development of statistical estimation procedures for an EGP. In his seminal paper, C ¸inlar (1980) proposes an iterative procedure that seems difficult to use in practice (since it needs notably a test for deciding whether a sample path comes from an ordinary Gamma process), whereas [START_REF] Guida | A time-discrete extended gamma process for time-dependent degradation phenomena[END_REF], in a parametric context, apply an approximate maximum likelihood method. The study of estimation procedures for an EGP is in progress and will be the subject of a future paper.

Fig. 4 Fig. 5

 45 Fig.4The cdf of Xt as a function of x for a(t) = t, b(t) = (t + 1) 2 and t = 10

  Fix B. For i in 1 : N , repeat the following steps : 2. Simulate the total number of jumps N B according to the Poisson distribution P (a(T )B). 3. Simulate N B realizations U 1 , . . . , U N B according to the uniform distribution on [0, B]. 4. Simulate N B realizations W 1 , . . . , W N B according to the exponential distribution with parameter 1. 5. Simulate N B realizations Z 1 , ..., Z N B according to the uniform distribution on

1.

Table 2

 2 Mean (standard deviation) of the relative errors for the four series representations of Method 2 for a(t) = t, b(t) = (t + 1) 2 and t = 1

		ILM	Bond	Thin	Rej
	simulation parameter B = 10	B = 10	B = 10	B = 10
	cpu time (s)	3	1	1	1
	mean	0.0027 (0.0022) 0.0028 (0.0020) 0.0246 (0.0033) 0.0026 (0.0021)
	variance	0.0087 (0.0066) 0.0089 (0.0064) 0.0084 (0.0064) 0.0086 (0.0062)
	Laplace transform	0.0017 (0.0011) 0.0021 (0.0008) 0.0511 (0.0022) 0.0021 (0.0010)

Table 3

 3 Mean (standard deviation) of the relative errors for the different methods of simulation of Xt for a(t) = t, b(t) = (t + 1) 2 and t = 1

		DSI	Bond	Rej	DRF
	simulation parameters K * = 20	B = 10	B = 10	ε = 0.035
	cpu time (s)	1	1	1	1
	mean	0.0026 (0.0020) 0.0028 (0.0020) 0.0026 (0.0021) 0.0028 (0.0022)
	variance	0.0108 (0.0077) 0.0089 (0.0064) 0.0086 (0.0062) 0.0088 (0.0070)
	Laplace transform	0.0020 (0.0013) 0.0021 (0.0008) 0.0021 (0.0010) 0.0017 (0.0011)

Table 4

 4 Mean (standard deviation) of the relative errors for the different methods of simulation of Xt for a(t) = 2t, b(t) =

	1 (t+1) 2 and t = 2
		DSI	Bond	Rej	DRF
	simulation parameters K * = 30	B = 10	B = 10	ε = 0.25
	cpu time (s)	1	1	1	1
	mean	0.0014 (0.0011) 0.0014 (0.0011) 0.0014 (0.0011) 0.0015 (0.0011)
	variance	0.0301 (0.0070) 0.0056 (0.0040) 0.0056 (0.0041) 0.0055 (0.0041)
	Laplace transform	0.0038 (0.0012) 0.0010 (0.0007) 0.0011 (0.0007) 0.0011 (0.0007)

Table 5

 5 Mean (standard deviation) of the relative errors for the different methods of simulation of Xt for

		DSI	Bond	Rej	DRF
	simulation parameters K * = 20	B = 30	B = 30	ε = 0.021
	cpu time (s)	1	1	1	1
	mean	0.0030 (0.0024) 0.0032 (0.0024) 0.0031 (0.0025) 0.0030 (0.0023)
	variance	0.0085 (0.0066) 0.0092 (0.0070) 0.0086 (0.0065) 0.0088 (0.0066)
	Laplace transform	0.0019 (0.0012) 0.0019 (0.0011) 0.0019 (0.0013) 0.0019 (0.0011)

a(t) = 1 -exp(-t), b(t) = 1 -exp(-t -1) and t = 1

Table 6

 6 Mean absolute error on F X t (x) for a(t) = t, b(t) = (t + 1) 2 , t = 2 and x ∈ [0.01 : 0.1 : 4]

	Method	Mean absolute error cpu time
	Post-Widder (H = 50)	0.0033	2
	Post-Widder (H = 100)	0.0015	4
	Post-Widder (H = 110)	0.0016	4
	Post-Widder (H = 200)	0.0645	7
	DRFI (ε = 10 -2 , K = 12)	0.0014	3
	DRFI (ε = 10 -2 , K = 200) 2.10 -5


T T Trang Nghiêm

P.-Olivier Chapuis

Phononic thermal resistance due to a finite periodic array of nano-scatterers

The wave property of phonons is employed to explore the thermal transport across a finite periodic array of nano-scatterers such as circular and triangular holes. As thermal phonons are generated in all directions, we study their transmission through a single array for both normal and oblique incidences, using a linear dispersionless time-dependent acoustic frame in a two-dimensional system. Roughness effects can be directly considered within the computations without relying on approximate analytical formulae. Analysis by spatio-temporal Fourier transform allows us to observe the diffraction effects and the conversion of polarization. Frequency-dependent energy transmission coefficients are computed for symmetric and asymmetric objects. We demonstrate that the phononic array acts as an efficient thermal barrier by applying the theory of thermal boundary (Kapitza) resistances to arrays of smooth scattering holes in silicon for an exemplifying periodicity of 10 nm in the [5-100 K] temperature range. It is observed that the associated thermal conductance has the same temperature dependence than that without phononic filtering.

I.

INTRODUCTION

Characterizing heat flow at interfaces [START_REF] Cahill | Nanoscale thermal transport[END_REF][START_REF] Hopkins | Thermal transport across solid interfaces with nanoscale imperfections: effects of roughness, disorder, dislocations, and bonding on thermal boundary conductance[END_REF], where thermal boundary resistances (TBR, also called Kapitza resistances) can impede the heat transfer, is of critical importance at the sub-micron length scales due to the high surface-to-volume ratio. Adding carbon nanotubes to the metal/metal interface [START_REF] Xu | Enhancement of thermal interface materials with carbon nanotube arrays[END_REF][START_REF] Mcnamara | Characterization of nanostructured thermal interface materials-a review[END_REF], or including self-assembled nanoparticles embedded in the material [START_REF] Shukla | Thermal transport in composites of self-assembled nickel nanoparticles embedded in yttria stabilized zirconia[END_REF] [START_REF] Losego | Effects of chemical bonding on heat transport across interfaces[END_REF] are some experimental examples that can be used to tune the TBR, as well as chemical etching [START_REF] Duda | Systematically controlling Kapitza conductance via chemical etching[END_REF][START_REF] Hopkins | Controlling thermal conductance through quantum dot roughening at interfaces[END_REF][START_REF] Hopkins | Effects of surface roughness and oxide layer on the thermal boundary conductance at aluminum/silicon interfaces[END_REF]. Other strategies have been highlighted, such as the addition of a new material to the solid/solid coupling or modulating the roughness [START_REF] Merabia | Thermal boundary conductance across rough interfaces probed by molecular dynamics[END_REF]. Lattice dynamics [START_REF] Young | Lattice-dynamical calculation of the Kapitza resistance between fcc lattices[END_REF][START_REF] Zhao | Lattice-dynamical calculation of phonon scattering at ideal Si-Ge interfaces[END_REF], Green"s functions [START_REF] Zhang | The atomistic Green's function method: An efficient simulation approach for nanoscale phonon transport[END_REF] and molecular dynamics are often employed for the calculations [START_REF] Stevens | Effects of temperature and disorder on thermal boundary conductance at solid-solid interfaces: Nonequilibrium molecular dynamics simulations[END_REF][START_REF] Landry | Thermal boundary resistance predictions from molecular dynamics simulations and theoretical calculations[END_REF][START_REF] Han | Phonon interference and thermal conductance reduction in atomic-scale metamaterials[END_REF]. We observe that the two commonly applied models, the acoustic mismatch model (AMM) and diffuse mismatch model (DMM), do not necessarily lead to values that are comparable to the available experimental data [START_REF] Cahill | Nanoscale thermal transport[END_REF][START_REF] Landry | Thermal boundary resistance predictions from molecular dynamics simulations and theoretical calculations[END_REF][START_REF] Stoner | Kapitza conductance and heat flow between solids at temperatures from 50 to 300[END_REF][START_REF] Lyeo | Thermal conductance of interfaces between highly dissimilar materials[END_REF]. In all cases, it is clear that the shapes of embedded elements and/or interfaces have an important impact on the TBR [START_REF] Shukla | Thermal transport in composites of self-assembled nickel nanoparticles embedded in yttria stabilized zirconia[END_REF], [START_REF] Losego | Effects of chemical bonding on heat transport across interfaces[END_REF], [START_REF] Merabia | Thermal boundary conductance across rough interfaces probed by molecular dynamics[END_REF] [START_REF] Zhou | Relationship of thermal boundary conductance to structure from an analytical model plus molecular dynamics simulations[END_REF], [START_REF] Gotsmann | Quantized thermal transport across contacts of rough surfaces[END_REF].

Phonon coherence effects may provide a new way to control heat transfer properties at boundaries. Such effects linked to the wave nature of the phonons have recently been experimentally evidenced by the Chen group [START_REF] Luckyanova | Coherent phonon heat conduction in superlattices[END_REF] and the California group [START_REF] Ravichandran | Crossover from incoherent to coherent phonon scattering in epitaxial oxide superlattices[END_REF], showing in particular that peculiar effects may take place for nanometer sizes at low temperatures and up to temperatures close to ambient. In addition, phononic configurations, which involve periodic arrays of holes [START_REF] Hopkins | Reduction in the thermal conductivity of single crystalline silicon by phononic crystal patterning[END_REF] [24] [START_REF] Zen | Engineering thermal conductance using a two-dimensional phononic crystal[END_REF] or embedded particles [START_REF] Kim | Thermal conductivity reduction and thermoelectric figure of merit increase by embedding nanoparticles in crystalline semiconductors[END_REF], have shown a strong reduction of the effective thermal conductivity. It is debated if the reason is due to coherence in these particular experimental works or due to the involved sizes which appear larger than the thermal wavelengths, but many interesting proposals have been highlighted, from the GHz to the THz [START_REF] Maldovan | Narrow low-frequency spectrum and heat management by thermocrystals[END_REF][START_REF] Davis | Nanophononic metamaterial: Thermal conductivity reduction by local resonance[END_REF] range where periodicity is expected to play a key role for heat transport [START_REF] Hopkins | Reduction in the thermal conductivity of single crystalline silicon by phononic crystal patterning[END_REF][START_REF] Maldovan | Sound and heat revolutions in phononics[END_REF]. Most of the theoretical suggestions deal with crystals involving many periods and derive an effective thermal conductivity based on the infinitecrystal approach. Since boundaries can already decrease strongly the heat transfer, it may be possible to block the heat transfer with smaller structures.

Here, we study the transmission and the thermal resistance of a finite phononic structure consisting of a single array of periodic holes, by solving the elastic wave equation. The acoustic frame is particularly suitable to reproduce the low-temperature phonon behaviors [START_REF] Cleland | Thermal conductance of nanostructured phononic crystals[END_REF] and may help to disentangle the elastic and inelastic contributions to thermal resistances [START_REF] Hopkins | Relative contributions of inelastic and elastic diffuse phonon scattering to thermal boundary conductance across solid interfaces[END_REF]. Various hole shapes are considered, such as disks, equilateral and isosceles triangles. We note that roughness can be directly included in the model by designing associated shapes. In contrast to previous works that analyzed phonon transmission through such single array with the goal of highlighting non-symmetrical acoustic transmission (sometimes improperly called "acoustic rectification" [START_REF] Maznev | Reciprocity in reflection and transmission: what is a 'phonon diode'?[END_REF]), we analyze the acoustic transmission not only for the normal incidence but also for oblique cases. This is required because thermal phonons are generated randomly in all directions. It allows observing reciprocity also for the case of asymmetric phononic structures. We perform an analysis of the displacement fields in order to determine if phonons of particular wavelengths and direction of propagation are especially filtered by the single array. Finally, we calculate the thermal conductances of the phononic array with the help of a Landauer-based approach which is similar to the AMM theory of thermal boundary resistance, both at equilibrium and out of equilibrium. A 10 nm-periodicity is considered for this example in two dimensions, showing that the structure indeed blocks a large portion of heat depending on the geometrical parameters associated to the considered shapes.

The article is organized as follows: In Sec. II, we describe (i) the structure, (ii) how we simulate phonon propagation with the elastic wave equation, (iii) the implementation and the computational domain, and (iv) present briefly the derivation of the thermal conductance based on the frequency-dependent transmission coefficients. In Sec. III, we (i) analyze the results of the spatial filtering resulting in diffraction effects and (ii) compute the transmission coefficients and the thermal conductance in the presence of the hole array. Finally, we present a summary and the consequences of this work in Sec. IV.

II. Method

A. Simulated structure and computational domain

The illustration of one simulated structure in two dimensions (2D) is shown in Figure 1. The periodicity in the y-direction is a and we simulate often more than one period, as will be explained in Sec. II.C. Empty holes with specific shapes are located in the middle of the simulated domain. Periodic boundary conditions are applied at the bottom and the top of structure. Two absorbing layers are created to avoid the reflection at left and right of the system [START_REF] Berenger | A perfectly matched layer for the absorption of electromagnetic waves[END_REF][START_REF] Hastings | Application of the perfectly matched layer (PML) absorbing boundary condition to elastic wave propagation[END_REF][START_REF] Johnson | Notes on perfectly matched layers (PMLs)[END_REF]. The size of these layers is large enough to ensure that the fluxes at the end walls of the system are null. We recall the approach that we use to derive the elastic equation in the absorbing domain in Appendix A. Acoustic waves are generated by the source at the left end of the lossless computational domain (left red line) and detected at its right end (right red line), and vice versa. The comparison between the propagation from left to right and right to left is particularly useful for the analysis of asymmetric objects. The distances between the single array and the source/detector are equal to L source = 10, where  is the wavelength. The mesh is chosen so that both x and y are always smaller than /5. The wavelengths can be compared to the periodicity by introducing a non-dimensioned positive number N  = a/. Hence, the circular frequency is also related to the medium velocity v and the periodicity a:

2 v N a     , ( 1
)
where v is the velocity for the considered polarization. Three shapes of holes are analyzed in this study: (a) circular, (b) equilateral triangle, (c) right isosceles triangle. To compare the area of these holes, we introduce a "filling factor" f defined by the ratio between the hole area and the a-side square area:

h o le s q u a r e S f S  . ( 2
)
This can be seen as the filling factor of the corresponding phononic crystal (of 2D periodicity). In addition, we consider what we will call the "blocking ratio" defined as the ratio of the projected length l to the periodicity a: l r a  .

(3) Fig. 2 presents three hole shapes with the same filling factor f = 0.2, with blocking ratios which are different. The values of the blocking ratios increase with the following order: disk, equilateral triangle, and then isosceles triangle. FIG 2. Three different hole shapes studied in this work: disk, equilateral triangle and right isosceles triangle, respectively from left to right. The periodicity is a and the filling factor f = 0.2 is identical for these three cases.

B. Phonons as elastic waves

As phonons are propagative acoustic waves, we solve the linear elastic equation in two dimensions (2D) to compute the phonon transmission experiment through an array:

      2 2 , , u r t r T r t t        , ( 4
)
where  is the mass density, u  is displacement field and T is the stress tensor. The stress tensor relates to the elastic constant C ijkl by

        , , 1 , 2 k l ij ijk l l k u r t u r t T r t C r x x                
. In the following, we use the usual abbreviated subscripts for elastic constants [START_REF] Auld | Acoustic fields and waves in solids[END_REF] and require only three constants C 11 , C

      * , , , 2
v r t T r t P r t       , ( 5
)
where v  is the velocity calculated by derivation of the displacement field / d u d t  . In the following, we consider a bi-dimensional system and the Poynting vector is reduced to a vector with two components. The projection of the flux on the x-direction allows monitoring the energy propagation across the array:

        * * , , , , , 2 x
        . ( 6
)
We note that this method can be used to study phonon transmission for many types of structures, such as boundaries between two dissimilar solids in direct contact (see Appendix B). It may provide a mesoscopic alternative to atomistic techniques such as molecular dynamics, especially at low temperatures where the latter does not behave well. Note that it could also be extended to nonlinear media.

C. Condition linking oblique waves and the simulated domain

As mentioned in Sec. I, a thermal source emits phonons to its surrounding environment in all directions. Hence, to calculate the heat flux generated and transmitted through a structure, all directions have to be included. We excite not only the normal acoustic waves that are perpendicular to the periodic direction, but also oblique waves (see Fig. 3). As longitudinal and transverse wave propagations can be separated [START_REF] Auld | Acoustic fields and waves in solids[END_REF] (see more in Sec. III.A), we show here how to simulate the longitudinal waves.

The acoustic wave displacement is expressed as

    0 0 x y i t k r i t k x k y u u e u e            , (7)
where k

  

. The periodic condition in the y-direction requires that at one given position x 0 , the displacement field at the top and the bottom are the same. Fig. 3b shows two typical points y 1 and y 2 at the boundaries and illustrates the condition on the incident angle  

, ux   
and the number of simulated rows N rows . The condition is written as

1 0 2 0 , , 2 2
r o w s r o w s a N a N u y x u y x                 (8)
Taking Eq. 8 into account, we have

0 0 0 0 e x p s in c o s 2 e x p s in c o s 2 r o w s r o w s aN u i t k k x aN u i t k k x                                    (9)
Eq. 9 leads easily to the relation

2 s in r o w s v n aN     , ( 10
)
where n is an integer that allows to satisfy 0 s in 1    . By taking Eq. 1 into account, we obtain the final condition:

s in r o w s n NN    . ( 11
)
As we consider angles  from 0 to π/2, n takes values between 0 and N  N rows . Consequently, for a given frequency N  and a given number of simulated rows N rows , we can simulate certain waves with incident angles satisfying the condition defined by Eq. 11.

Fig. 4 illustrates the propagation of the waves by presenting a snapshot of the u x component of the displacement field for two incident angles in the case N  = 2.5 (as an example, for a = 10 nm,  = a/2.5 = 4 nm): (a) normal incident angle ( = 0), and (b) oblique incident angle with sin = 0.4. In this example the stationary regime has not been reached. The scattering objects separate the displacement field into two regions: a first one at left where incident and reflected waves are observed; a second one at right that contains the fields associated to the transmitted waves, where the propagation has been obviously modified. In addition, the interference between incident and reflected waves is observed in the first region. We note in particular that a u y component is generated immediately after the objects even though only the u x component is excited in the normal incident case (not shown here). This is due to the fact that the propagation direction is changed and is not anymore only in the normal incidence for longitudinal waves; the projection of the displacement field on the yaxis is non-null. Animations related to normal and oblique incidence of phonons on the single array during 313 wave periods can be found in [38]. As waves travel through a grating-like periodic structure, they can be diffracted according to Bragg"s law expressed as

s in n a n    , ( 12
)
where n is an integer which corresponds to the diffraction order characterized by the angle  n (see also Appendix C).

Considering the wave shown in Fig. 4 with N  =2.5, the diffraction angles are  1 = 1/2.5 = 23.6° and  2 = 2/2.5 = 53.1°. In the following, the magnitude of these waves will be determined by the spatio-temporal Fourier transform.

D. Equilibrium thermal conductance in 2D and frequency-dependent transmission coefficients

The TBR (Kapitza) measures the boundary resistance to the propagation of the thermal flux. The theory has been applied to predict the thermal resistance of different types of junction [START_REF] Landry | Thermal boundary resistance predictions from molecular dynamics simulations and theoretical calculations[END_REF][START_REF] Han | Phonon interference and thermal conductance reduction in atomic-scale metamaterials[END_REF][START_REF] Merabia | Thermal conductance at the interface between crystals using equilibrium and nonequilibrium molecular dynamics[END_REF][START_REF] Liang | Finite-size effects on molecular dynamics interfacial thermal-resistance predictions[END_REF]. Here, we consider the hole array as a barrier between two parts of same material. In the following, we recall how to apply the formula for the 2D case. The thermal conductance G between two media is defined as:

2 1 q G TT   ( 13
)
where q is the heat flux across the junction, and T i is the temperatures at the lead i. The heat flux across the structure relates to the transmitted phonons. At steadystate, the net heat flux in 2D is

                2 1 1 2 1 2 2 2 2 1 2 2 d , , , 2 d , , , 2 x k x k k q v k t k f k T k v k t k f k T                           ( 14
)
where ћ is the reduced Plank constant, v ix denotes the phonon velocity in medium i projected along the direction x normal to the array,  

, i f k T  is the phonon distribution function at the medium temperature T i ,   , ij t k  
is the wave-vector dependent transmission coefficient from the medium i to the medium j. The signs + andindicate that the integrals deal with k x > 0 and k x < 0, respectively. When the thermal transport is close to the equilibrium state, the phonon distribution can be assimilated to the equilibrium Bose-Einstein distribution that does not anymore depend on the direction:

  1 , e x p 1 BE B f T kT                   , ( 15
)
where k B is the Boltzmann constant. In a lossless medium, by invoking the balance principle when the heat flux is zero at thermal equilibrium [START_REF] Landry | Thermal boundary resistance predictions from molecular dynamics simulations and theoretical calculations[END_REF], Eq. 14 can be simplified to

          2 1 2 2 1 2 d , , , , 2 x B E B E k k q v k t k f T f T                    . ( 16
)
In 2D coordinates, 2

2 d d . d d . d k k k v       
, where  is the angle between k  and the x axis, , 2 2

        
and therefore v x = v cos. We can perform the integration on 0, 2

        
and the 2D thermal conductance at equilibrium

d d q G
T  associated to one polarization state is finally calculated as

            3 2 12 2 2 0 /2 3 12 2 2 0 d , 1 2 , c o s d d d 2 , 1 2 , c o s d . 2 BE eq BE f T G t v T d f T d t v d T                                     (17)
We introduce the frequency-dependent transmission coefficient  12 (ω) as the transmission coefficient including all wavevector directions:

    /2 1 2 1 2 0 2 , c o s d t             , (18)
so that e 2D thermal conductance writes:

      3 12 2 2 d , 1 d d 2 BE eq f T G v T           . ( 19
)
The maximal monochromatic thermal conductance (

  ) for one polarization state is then given by

    3 , e q , m a x 2 2 d , 1 d 2 BE f T g v T       . ( 20
)
The expression of Eq. 17 is established at equilibrium. Out of equilibrium, a discontinuity of temperature occurs, which can be accounted for by using the local distribution () fr  [START_REF] Merabia | Thermal conductance at the interface between crystals using equilibrium and nonequilibrium molecular dynamics[END_REF], which obeys the Boltzmann transport equation

(BTE) with         BE f r f T r f r       . ( 24
)
Including the deviation of the distribution function f in the relaxation time approximation of the BTE, the nonequilibrium thermal conductance can finally be written as [START_REF] Merabia | Thermal conductance at the interface between crystals using equilibrium and nonequilibrium molecular dynamics[END_REF] 1 2 2 1

1 eq n e q G G      (25)
involving the quantity

2 1 2 1 1 1 2 1 d/ eq s c a tt x k f v t k T             (26)
where  scatt,1 is the relaxation time and   is the thermal conductivity in material 1.  21 has a similar definition. To be consistent with the 2D conduction, these fractions are also calculated in 2D. By plugging the AMM expression into Eq. 25, the conductance G neq of this model has been shown to compare well with molecular dynamics simulations [START_REF] Merabia | Thermal conductance at the interface between crystals using equilibrium and nonequilibrium molecular dynamics[END_REF].

E. Numerical implementation

The simulation allows us to compute the average flux density at a given position x 0 including the contributions of all waves:

    m a x m in 0 0 m a x m in 1 ,d y x x y P x P x y y y y    (23)
with y max -y min being a multiple of a, and P x (x 0 ,y) is defined in Eq. 5. The frequency-dependent transmission coefficient can be calculated as the ratio between the sum of transmitted heat flux intensities and the sum of the incident heat flux intensities 0 P :

      0 , , c o s d d c o s d i t i t i t t t i i i P P                           (
        
, where the index x stands for the projection on the x axis. The frequency-dependent coefficient is computed by discretizing Eq. 21:

      , 0 , i tx i i P P          . ( 22
)

III. Results

A. Diffraction effect -Spatial filtering

To analyze the transmission and reflection processes, we use the spatio-temporal Fourier Transform (FT) [START_REF] Alleyne | A two-dimensional Fourier transform method for the measurement of propagating multimode signals[END_REF]. It has been used to investigate the dispersion relation of surface acoustic modes [START_REF] Sugawara | Direct access to the dispersion relations of multiple anisotropic surface acoustic modes by Fourier image analysis[END_REF], acoustic band structure [START_REF] Profunser | Imaging ripples on phononic crystals reveals acoustic band structure and Bloch harmonics[END_REF] and to observe the diffraction angles [START_REF] Danworaphong | Real-time imaging of acoustic rectification[END_REF]. In infinite space and time domains, the FT is written as

      , , d d i t k r U k t u r t e r t               . ( 27
)
By limiting to a given spatial domain and integrating over one period T = 2π/ω, the spatio-temporal FT becomes

      0 0 1 1 0 0 /2 /2 / 2 / 2 1 1 1 , , , e d d d y x x y x y y l x l t T i t k x k y x y t x l y l U k t u x y t x y t l l T              ; (28)
[x 0 -l x /2, x 0 + l x /2] and [y 0 -l y /2, y 0 + l y /2] are the bounds of the analyzed domain in x and y directions. We consider a time t 1 for which the stationary regime is already established, l y proportional to a, x 0 at left or right of the scattering elements with l

x such that the addressed area is in the far field of the scattering elements. Both

              ( 29
)
where  is the angle between x-axis and the wavevector k  . Fig. 5 shows the maps of spatio-temporal FT of the displacement field in reciprocal space with N  =2.5, i.e. =a/2.5. The independent propagation of longitudinal u l and transverse u t modes are for two cases: (a-b-e-f) circular holes, (c-d-g-h) triangular holes. The two dotted circles represent the iso-wavevector curves for longitudinal waves (k l = ω/v l ) and for transverse waves (k

t = ω/v t ).
As the longitudinal velocity v l is larger than the transverse one v t , the inner circle corresponds to k l and the outer one corresponds to k t . In this figure, the wave vectors are normalized by the maximal longitudinal one.

In Fig. 5a and5c, the incident waves in normal incidence ( = 0) are represented by the points on the longitudinal circle, with k

x > 0 and k y = 0; this shows the propagation in the positive direction. On the same circle, the points of negative k

x represent the diffracted waves: the centered point is associated to the reflected wave (order 0), the next two symmetric points, which have diffraction. These ratios are exactly the sine values of diffracted angles  1 and  2 , as predicted by Bragg"s law.

In addition, diffracted waves are also represented on the transverse circle. This shows the effect of conversion of polarization due to the arrays. In the transmitted region, diffracted waves that travel through the single array are represented with k x >0 in Fig. 5b and5d. According to Snell"s law, these transverse waves have the same k y component than the longitudinal wave of identical diffraction order. For an oblique wave with incidence such that sin = 0.4, the incident waves in Fig. 5e and5g are found on the iso-longitudinal wavevector circle according to k y /k l = 0.4. As a consequence, the zerothorder diffracted waves also have this ratio. In addition, higher orders are shifted of n/a on the k y -axis. The analysis in the two regions remains similar to the normal incidence case. While the amplitude of each diffracted wave depends on the hole shape, Bragg"s law and Snell"s law apply for all structures. Diffraction plays an important role on the propagation direction, leading to spatial filtering effect. Note that for long wavelengths ( > a  N  <1) the propagation direction of the transmitted wave is the same as the one of the incident wave. First, we analyze the frequency-dependent transmission coefficients as a function of the number of incident angles. As the coefficient depends on cosine values, the large angles, especially those near π/2, are less important than the small ones. In our simulations incident angles are characterized by sine values in [0-0.89]. Fig. 6 shows the frequency-dependent coefficients as a function of number of incident angles for the case N  = 3.0. In this example, the sine values included are 0, 0.11, 0.25, 0.4, 0.53, 0.67 and 0.85. The convergence of the coefficient value appears to be obtained for a rather small number of angles. We have verified that increasing this number up to 18 does not significantly improve the calculations in our case. In the following, we calculate the total transmission coefficients and the thermal conductance with 7 angles for each frequency. We observe in Fig. 6 (and Fig. 7b) that the frequencydependent transmission coefficients associated to heat fluxes impinging the bases or the vertices of the asymmetric triangles are equal when all propagation directions are included. This is a manifestation of reciprocity, which is fulfilled in our lossless and linear acoustic system (see more in [START_REF] Maznev | Reciprocity in reflection and transmission: what is a 'phonon diode'?[END_REF]). Indeed no rectification can be observed in absence of non-linear mechanism. This is different to the case of acoustic waves only excited and observed in normal incidence [START_REF] Danworaphong | Real-time imaging of acoustic rectification[END_REF][START_REF] Zanjani | One-way phonon isolation in acoustic waveguides[END_REF]. We also note in Fig. 6 that the frequency-dependent transmission coefficients depend strongly on the filling factor f: as an example, the lower the filling factor ratio, the larger the transmission.

B. Transmission coefficient and 2D thermal conductance

Fig. 7 shows the frequency-dependent transmission coefficients as a function of the frequency 2 v/a) = N  for (a) circular holes and (b) equilateral and right isosceles triangular holes. For all the presented longitudinal cases, we observe that the curves have same trend. The transmission coefficients reach unity at low frequency, and then reduce to a certain value when increasing the frequency. Two peaks at N  = 1.5 and 3.5 are observed, while the non-zero minimum values take place for N  = 2.5. However, the values of the coefficient remain quite close for each filling factor/blocking ratio for frequencies larger than N  = 0.5 and we do not observe sharp features in the spectrum. This may be due to the "thermal" averaging due to the integration over all the excited angles. For the circular holes of Fig. 7a, the transmission coefficients are around 65% for f = 0.1 (r = 0.36), then around 55-60% for f = 0.2 (r = 0.50), and finally around 30% for f = 0.1 (r = 0.48). The same order is obtained for equilateral triangular holes, the transmission coefficients are the largest for f = 0.1 (r = 0.48), being close to only 60%. This means that phonons are already efficiently blocked and diffracted for a modest density of scatterers. The values of the coefficients are around 40% for f = 0.2 (r = 0.68), then less than 10% for f = 0.4 (r = 0.96). The same trend is observed for transverse waves (not shown here). Comparing the total coefficients with f = 0.2, they are different for each shape. The values of the frequency-dependent transmission coefficients decrease with the following order: disk, equilateral and right isosceles holes, in agreement with the values of the associated blocking ratios. Considering two cases with close values of r but different values of f -(i) circular array with f = 0.2, r = 0.50, and (ii) triangular array with f = 0.1, r = 0.48 -we observe close transmission coefficients through these arrays, around 60%. Despite the fact that wave scattering is often related to the area associated to the scatterer, the blocking ratio may be also a convenient way to describe the transmission. Let us note that we have considered ideal shapes until now, neglecting possible roughness on the walls of the holes. In contrast to many other works, such roughness can be accounted for by designing directly complex shapes without relying on approximate analytical expressions. For example, a roughness of =0.3 nm was introduced to the case f = 0.2 studied in Fig. 7a by extruding half disks and grafting them between the extruded area (see inset in Fig. 7a). Note that such roughness is not random but identical on each hole if no supercell is considered. A strong decrease of the transmission can be observed, reaching 17%. Roughness especially impacts the transmission coefficients at high frequencies when short wavelengths become comparable to the roughness characteristic size. In the following, we restrict our study to smooth shapes. The single array of periodic holes acts as a thermal barrier to which a thermal conductance can be associated. In order to calculate the thermal conductance as described in Section II.D, we consider temperatures exciting thermal wavelengths which are commensurate with the geometric parameters of the array such as its periodicity.

Here, a periodicity a = 10 nm is chosen as an example. The frequency dependence of the maximal monochromatic thermal conductance defined in Eq. 20 (section II. D) is considered. The thermal frequencies associated to temperatures ranging from 5 K to 100 K are well located in the simulated frequency band 2v/a) = [0.025 -4.5]. We have checked that the summation of monochromatic thermal conductances g ,max with our discretization fills over 90% of the values of the integral for each temperature, as is customary for bosons which possess broad Note that the time-domain simulations are computer intensive (few hours of CPU time on a 12-core processor with RAM 64GB for one simulation) and that the final computation requires summation over frequency and angle.

The equilibrium thermal conductance G eq is shown in Fig. 8a for circular and triangular hole-based single arrays. The 2D conductance without the presence of hole G 0 (integral of Eq. 20) is also plotted for comparison The thermal conductance decreases when the filling factor increases, as the transmission coefficients are reduced (see Eq. [START_REF] Zhou | Relationship of thermal boundary conductance to structure from an analytical model plus molecular dynamics simulations[END_REF]. Moreover, for the same filling factor the thermal conductance of the triangular array is always smaller than that of the circular one. We remark that the temperature dependence of the thermal conductance in all cases is quadratic, i. e. G eq ∝ 𝑇 2 : the modulation of the Bose-Einstein factor by the transmission coefficients does not lead to a different thermal behavior. The reduction of thermal conduction in the presence of single arrays normalized to G 0 is plotted in Fig. 8b. With f = 0.1, the relative conductance reduces to 65%-62% for the two hole types, while the conduction drops to 30% and 10% with f = 0.4.

To obtain the non-equilibrium conductance (Eq. 25), we take into account phonon-phonon Umklapp scattering.

Our purpose is more to analyze how non-equilibrium affects the values of the phononic thermal conductance than to properly account for phonon volume scattering. Slack and Galginaitis [START_REF] Slack | Thermal conductivity and phonon scattering by magnetic impurities in CdTe[END_REF] suggested the following form for the Umklapp process:

  /3 1 2 D T U U B T e        (30)
where B U is a fitted parameter and  D is the Debye temperature. We consider the values of parameters calibrated to reproduce the experimental thermal conductivity of silicon in Ref. [START_REF] Morelli | Estimation of the isotope effect on the lattice thermal conductivity of group IV and group III-V semiconductors[END_REF]. They are  L = 586 K and B U L = 5.5×10 -20 s -1 K -3 .

The results are shown in Fig. 9 for the case of longitudinal waves propagating through an array of triangular holes. As the same material is present on the two sides of the array, one can define only one fraction  =  12 =  21 . Due to reduction of the transmission coefficient when widening the holes, the same trend is obtained for  for temperatures from 5K to 100K. The non-equilibrium conductances are always higher than the equilibrium ones for all cases (see Eq. 25). The shift of the non-equilibrium conductance with respect to the equilibrium one increases with the transmission coefficient, i.e. when reducing the element size. The main conclusions of the preceding sections stay valid.

IV. Conclusion

In summary, we have investigated the transmission of heat though a finite periodic array of scatterers, here holes, by solving the elastic wave equation in finite geometries. The acoustic waves, which model thermal phonons, have been excited in various directions to model the thermal "emission" by a heat source. By analyzing the spatio-temporal FT of the displacement fields, we observed (i) that the periodicity follows Bragg"s law of diffraction, and (ii) that the conversion of polarization takes place according to Snell"s law. The amplitude of the transmission of each phonon mode varies as a function of the scattering object shape, and we found that there is strong filtering effect even for modest size of the elements, such as for an equivalent filling ratio f = 0.1. In particular, we observed that rough scatterers lead to lower transmission than smooth ones.

The frequency-dependent transmission coefficients of three types of arrays were found to have similar trends when varying the frequency. We noticed that the "blocking ratio" appears useful to determine if phonons can be transmitted: despite of the different shapes, the frequency-dependent transmission coefficients through objects of similar ratios are very close. Furthermore, for triangular (asymmetric) objects, the frequency-dependent coefficients associated to the heat flux impinging (i) the bases and (ii) the vertices have the same values. This is a consequence of reciprocity. We also observed that there is no particular frequency at which sharp spectral features such as a total gap could appear.

We highlighted that the single array acts as a thermal barrier, in a way very-closely related to thermal (Kapitza) boundary conductances. The phononic thermal conduction in the 2D case was characterized for these structures for temperatures ranging from 5K to 100K, in the example of isotropic silicon. Due to the impact on the frequency-dependent transmission coefficient, we observed that the hole size strongly influences the thermal conductance. A reduction of 90% could be reached for a blocking ratio r = 0.71. We also noticed that the presence of the array does not change the temperature dependence of the equilibrium thermal conductance: it merely filters the whole spectrum once a certain frequency is reached. We finally considered the impact of non-equilibrium close to the array.

This work is an important step related to the study of the wake-like phonon scattering properties through artificial interfaces, which are created by adding a periodic structure which is finite in at least one dimension. It may provide a basis for future investigations dealing also with non-linear mechanisms [START_REF] Liang | Acoustic diode: Rectification of acoustic energy flux in one-dimensional systems[END_REF][START_REF] Chen | A photon thermal diode[END_REF] expected to exhibit thermal rectification effects.

Assuming that no inelastic scattering takes place at the interface, the transmission coefficient t 12 through a perfect interface between two media 1 and 2 is given in the AMM framework by     . The vertical dashed line shows the critical limit of 35.7° the waves come from the Ge medium.

Appendix C: Periodic conditions lead to diffraction

The periodic condition for wave propagation in the simulated domain, in particular for the oblique ones, is expressed in Eq. 11:

C.1. Propagation through a perfect interface

According to Snell"s law, the refraction angle at the interface between two different materials is defined from the condition that can be rewritten This shows that the sine value of the refraction angle satisfies the periodic condition in medium 2. This is due to the well-known condition that the tangential component is conserved when crossing an interface.

C.2. Propagation through a periodic array

The diffracted waves satisfy Bragg"s law as expressed in Eq. 12: ). By replacing (1) in [START_REF] Zhao | Lattice-dynamical calculation of phonon scattering at ideal Si-Ge interfaces[END_REF], we obtain: The numerator of Eq. C.6 is obviously an integer. By comparing Eq. 11 and Eq. C.6, we verify that the diffracted waves are satisfying the periodic conditions.

FIG 1 .

 1 FIG 1. Illustration of a simulated structure involving two rows with triangular shapes. Absorbing conditions are applied in the two blue regions (see more in Appendix); phonon source and detector are highlighted with the red lines.

,

  cubic crystal case. To simplify, the material is assumed isotropic, so that C 12 = Eq. 1 with the finite element method[37]. We obtain the time evolution of the displacement field   at each point in the system. This allows computing the acoustic Poynting vector   , P r t   that carries the energy flux:

FIG 3 .

 3 FIG 3. (a) Illustration of waves impinging the scatterer from different incident angles. (b) Periodic condition for the oblique incidence.



  is the wave vector, kx and k y are the projection of k on the x and y axis, respectively.

FIG 4 .

 4 FIG 4. u x displacement field close to a single array of equilateral triangular holes: (a) wave in normal incidence, and (b) wave in normal incidence with sin = 0.4. Two regions are defined: the one with incident and reflected (IR) waves, and the one with transmitted (T) waves.

  u x and u y components are analyzed with FT. In order to examine the conversion of polarization, we combine these two components to determine the relative contributions of the longitudinal and the transverse modes. The longitudinal and transverse amplitudes are obviously calculated with

  FIG 5. Maps of spatio-temporal TF in reciprocal space for N  = 2.5 in the case of circular holes (a-b-e-f) and in the case of equilateral triangular holes (c-d-g-h) : (a-c) incident region and (b-d) transmitted region for a wave in normal incidence ( = 0); (e-f) incident region and (g-h) transmitted region for a wave in oblique incidence with sin = 0.4.The dashed lines represent the iso-wavevector curves: the inner circle for the longitudinal one, the outer circle for the transverse one. The orders of diffraction are also shown in the maps (0 for the zeroth order, 1 for the first order, 2 for the second order).

1 .FIG 6 .

 16 FIG 6. Frequency-dependent transmission coefficient  as a function of the number of incident angles for N  = 3.0 in two cases: (i) incidence toward the triangle bases (solid lines) and (ii) toward the vertices (dashed lines), both for a single array of equilateral triangles.

FIG 7 .

 7 FIG 7. Frequency-dependent transmission coefficients for arrays of (a) circular holes and (b) equilateral and right isosceles triangular holes (plain lines). In (a), the dashed line represents

FIG 8 .

 8 FIG 8. (a) Equilibrium phononic thermal conductance as a function of temperature in absence of hole array G 0 (solid black line), in presence of an array with circular holes (solid lines with symbols) and with equilateral triangular holes (dashed lines with symbols). (b) Ratio between the thermal conductance in

FIG 9 .

 9 FIG 9. Phononic thermal conductance as a function of temperature: in absence of hole array G 0 (solid black lineequilibrium case for reference), in presence of a triangular array at equilibrium (dashed line with symbols) and out of equilibrium (solid lines with symbols). Blue lines for filling factor f = 0.1, green lines for f = 0.2, red lines for f = 0.4.

1 ,

 1 FIG B.1. Transmission coefficients as a function of the angle of incidence for waves crossing the Si/Ge interface (a) from Si to Ge, (b) from Ge to Si. The blue solid lines present the AMM calculation and symbols are for the results obtained from

  m is an integer which corresponds to the diffraction order characterized by the angle  m (   0 mN  

Acknowledgements

This work has been supported by ANR project RPDOC NanoHeat and INSA BQR MaNaTherm. We acknowledge useful discussions with S. Merabia. P.O.C. also thanks F. Alzina, E. Chavez, J. Gomis and C. Sotomayor.

Appendix A: Wave equation in absorbing region

We recall the derivation of the equation used for the absorbing zones [START_REF] Johnson | Notes on perfectly matched layers (PMLs)[END_REF]. In 1D, the elastic wave equation is written as

This equation can be expressed in the following form:

In the absorbing zones, we introduce a new mass density and new elastic constants while keeping the same

Then, we look for a solution with a plane wave form

in the absorbing regions. We can take  

By combining Eqs. A2 and A3, we obtain

Finally, Eq. A. 4 can be rewritten in the real space as a third-order partial differential equation, or also as follows:

This linear equation, which mixes frequency and time, may appear as unusual, but it leads exactly to the same displacement field that the one of the third-order equation when a plane wave is excited at . It can be useful if one prefers only to solve second-order partial differential equations in the computational domain to avoid potential additional discretization requirements.

Appendix B: Transmission through a perfect Si-Ge interface

In this section, we compute the transmission through a perfect Si-Ge interface and compare the results with those obtained semi-analytically within the acoustic-mismatch model (AMM) [START_REF] Little | The transport of heat between dissimilar solids at low temperatures[END_REF]. This model captures the impedance mismatch effect of phonon transmission. At the interface, one part of wave is reflected and the other part is refracted at the other side of the interface following Snell"s law:

 1 and  2 are the incident and refraction angles, respectively. It is required that the incident angle be smaller than the critical angle  c = asin(v 2 /v 1 ). We consider longitudinal wave with velocity ratio v Si /v Ge  1.7. As a consequence, the critical angle of waves propagating from Ge towards Si is  c  35.°, while there is no angle limit in the opposite direction.