V Lemaire

S Menozzi
email: menozzi@math.jussieu.fr

ON SOME NON ASYMPTOTIC BOUNDS FOR THE EULER SCHEME

Keywords: 1991 Mathematics Subject Classification. 60H35, 65C30, 65C05, 60E15. January 8, 2010 Non asymptotic Monte Carlo bounds, Discretization schemes, Gaussian concentration

We obtain non asymptotic bounds for the Monte Carlo algorithm associated to the Euler discretization of some diffusion processes. The key tool is the Gaussian concentration satisfied by the density of the discretization scheme. This Gaussian concentration is derived from a Gaussian upper bound of the density of the scheme and a modification of the so-called "Herbst argument" used to prove Logarithmic Sobolev inequalities. We eventually establish a Gaussian lower bound for the density of the scheme that emphasizes the concentration is sharp.

Introduction

Statement of the problem

Let the R d -valued process (X t ) t 0 satisfy the dynamics

X t = x + t 0 b(s, X s )ds + t 0 Bσ(s, X s )dW s , (1.1)
where (W t ) t 0 is a d ′ -dimensional (d ′ d) standard Brownian motion defined on a filtered probability space (Ω, F , (F t ) t 0 , P) satisfying the usual assumptions. The matrix

B = I d ′ ×d ′ 0 (d-d ′ )×d ′ is the embedding matrix from R d ′ into R d . The coefficients b : R + × R d → R d , σ : R + × R d → R d ′ ⊗ R d ′
are assumed to be Lipschitz continuous in space, 1/2-Hölder continuous in time so that there exists a unique strong solution to (1.1).

Let us fix T > 0 and introduce for (t, x)

∈ [0, T ] × R d , Q(t, x) := E[f (T, X t,x T )],
where f is a measurable function, bounded in time and with polynomial growth in space. The numerical approximation of Q(t, x) appears in many applicative fields. In mathematical finance, Q(t, x) can be related to the price of an option when the underlying asset follows the dynamics (1.1). In this framework we consider two important cases: (a) If d = d ′ , Q(t, x) corresponds to the price at time t when X t = x of the vanilla option with maturity T and pay-off f .

(b) If d ′ = d/2, b(x) = b 1 (x) b 2 (x) where b 1 (x) ∈ R d ′ , b 2 (x) = (x 1 , • • • , x d ′ ) * , Q(t, x
) corresponds to the price of an Asian option.

It is also well known, see e.g. Friedman [START_REF] Friedman | Stochastic differential equations[END_REF], that Q(t, x) is the Feynman-Kac representation of the solution of the parabolic PDE ∂ t Q(t, x) + LQ(t, x) = 0, (t, x)

∈ [0, T ) × R d , Q(T, x) = f (T, x), x ∈ R d , (1.2)
where L stands for the infinitesimal generator of (1.1). Hence, the quantity Q(t, x) can also be related to problems of heat diffusion with Cauchy boundary conditions (case (a)) or to kinetic systems (case (b)).

The natural probabilistic approximation of Q(t, x) consists in considering the Monte Carlo algorithm. This approach is particularly relevant compared to deterministic methods if the dimension d is large. To this end we introduce some discretization schemes. For case (a) we consider the Euler scheme with time step ∆ := T /N, N ∈ N * . Set ∀i ∈ N, t i = i∆ and for t 0, define φ(t) = t i for t i t < t i+1 . The Euler scheme writes

X ∆ t = x + t 0 b(φ(s), X ∆ φ(s) )ds + t 0
σ(φ(s), X ∆ φ(s) )dW s .

(1.3)

For case (b) we define

X ∆ t = x + t 0 b 1 (φ(s), X ∆ φ(s) ) (X ∆ s ) 1,d ′ ds + t 0 Bσ(φ(s), X ∆ φ(s) )dW s , (1.4)
where

(X ∆ s ) 1,d ′ := (X ∆ s ) 1 , • • • , (X ∆ s ) d ′ *
. Equation (1.4) defines a completely simulatable scheme with Gaussian increments. On every time step, the last d ′ components are the integral of a Gaussian process.

The weak error for the above problems has been widely investigated in the literature. Under suitable assumptions on the coefficients b, σ and f (namely smoothness) it is shown in Talay and Tubaro [START_REF] Talay | Expansion of the global error for numerical schemes solving stochastic differential equations[END_REF] that

E D (∆) := E x [f (T, X ∆ T )] -E x [f (T, X T )] = C∆ + O(∆ 2
). Bally and Talay [START_REF] Bally | The law of the Euler scheme for stochastic differential equations: I. Convergence rate of the distribution function[END_REF] then extended this result to the case of bounded measurable functions f in a hypoelliptic setting for time homogeneous coefficients b, σ. Also, still for time homogeneous coefficients, similar expansions have been derived for the difference of the densities of the process and the discretization scheme, see Konakov and Mammen [START_REF] Konakov | Edgeworth type expansions for euler schemes for stochastic differential equations[END_REF] in case (a), Konakov et al. [START_REF] Konakov | Explicit parametrix and local limit theorems for some degenerate diffusion processes[END_REF] in case (b) for a uniformly elliptic diffusion coefficient σσ * , and eventually Bally and Talay [START_REF] Bally | The law of the Euler scheme for stochastic differential equations, II. Convergence rate of the density[END_REF] for a hypoelliptic diffusion and a slight modification of the Euler scheme. The constant C in the above development involves the derivatives of Q and therefore depends on f, b, σ, x.

The expansion of E D (∆) gives a good control on the impact of the discretization procedure of the initial diffusion, and also permits to improve the convergence rate using e.g. Richardson-Romberg extrapolation (see [START_REF] Talay | Expansion of the global error for numerical schemes solving stochastic differential equations[END_REF]). Anyhow, to have a global sharp control of the numerical procedure it remains to consider the quantities

E MC (M, ∆) = 1 M M i=1 f (T, (X ∆ T ) i ) -E x f (T, X ∆ T ) .
(1.5)

In the previous quantities M stands for the number of independent samples in the Monte Carlo algorithm and (X ∆ t ) i t 0 i∈[[1,M]] are independent sample paths. Indeed, the global error associated to the Monte Carlo algorithm writes:

E(M, ∆) = E D (∆) + E MC (M, ∆),
where E D (∆) is the discretization error and E MC (M, ∆) is the pure Monte Carlo error. The convergence of E MC (M, ∆), to 0 when M → ∞ is ensured under the above assumptions on f by the strong law of large numbers. A speed of convergence can also be derived from the central limit theorem, but these results are asymptotic, i.e. they hold for a sufficiently large M . On the other hand, a non asymptotic result is provided by the Berry-Esseen Theorem that compares the distribution function of the normalized Monte Carlo error to the distribution function of the normal law at order O(M -1/2 ).

In the current work we are interested in giving, for Lipschitz continuous in space functions f , non asymptotic error bounds for the quantity E MC (M, ∆). Similar issues had previously been studied by Malrieu and Talay [START_REF] Malrieu | Concentration inequalities for euler schemes[END_REF]. In that work, the authors investigated the concentration properties of the Euler scheme and obtained Logarithmic Sobolev inequalities, that imply Gaussian concentration see e.g. Ledoux [Led99], for multi-dimensional Euler schemes with constant diffusion coefficients. Their goal was in some sense different than ours since they were mainly interested in ergodic simulations. In that framework we also mention the recent work of Joulin and Ollivier for Markov chains [START_REF] Joulin | Curvature, concentration, and error estimates for Markov chain Monte Carlo[END_REF].

Our strategy is here different. We are interested in the approximation of Q(t, x), t T where T > 0 is fixed. It turns out that the log-Sobolev machinery is in some sense too rigid and too ergodic oriented. Also, as far as approximation schemes are concerned it seems really difficult to obtain log-Sobolev inequalities in dimension greater or equal than two without the constant diffusion assumption, see [START_REF] Malrieu | Concentration inequalities for euler schemes[END_REF]. Anyhow, under suitable assumptions on b, σ (namely uniform ellipticity of σσ * and mild space regularity), the discretization schemes (1.3), (1.4) can be shown to have a density admitting a Gaussian upper bound. From this a priori control we can modify Herbst's argument to obtain an expected Gaussian concentration as well as the tensorization property (see [START_REF] Ledoux | Concentration of measure and logarithmic sobolev inequalities[END_REF]) that will yield for r > 0 and a Lipschitz continuous in space f ,

P |E MC (M, ∆)| r + δ 2 exp(-M α(T ) r 2 ) for α(T ) > 0 independent of M uniformly in ∆ = T /N .
Here δ 0 is a bias term (independent of M ) depending on the constants appearing in the Gaussian domination (see Theorem 2.1) and on the Wasserstein distance between the law of the discretization scheme and the Gaussian upper bound. We also prove that a Gaussian lower bound holds true for the density of the scheme. Hence, the Gaussian concentration is sharp, i.e. for a function f with suitable non vanishing behavior at infinity, the concentration is at most Gaussian, i.e.

P |E MC (M, ∆)| r -δ 2 exp(-M ᾱ(T ) r 2
), for r large enough, δ depending on f , and the Gaussian upper and lower bounds, ᾱ(T ) > 0 independent of M uniformly in ∆ = T /N .

The paper is organized as follows, we first give our standing assumptions and some notations in Section 1.2. We state our main results in Section 2. Section 3 is dedicated to concentration properties and non asymptotic Monte Carlo bounds for random variables whose law admits a density dominated by a probability density satisfying a log-Sobolev inequality. We prove our main deviations results at the end of that section as well.

In Section 4 we show how to obtain the previously mentioned Gaussian bounds in the two cases introduced above. The main tool for the upper bound is a discrete parametrix representation of Mc Kean-Singer type for the density of the scheme, see [START_REF] Mckean | Curvature and the eigenvalues of the Laplacian[END_REF] and Konakov and Mammen [START_REF] Konakov | Local limit theorems for transition densities of Markov chains converging to diffusions[END_REF] or [START_REF] Konakov | Edgeworth type expansions for euler schemes for stochastic differential equations[END_REF]. The lower bound is then derived through suitable chaining arguments adapted to our non Markovian setting.

Assumptions and Notations

We first specify some assumptions on the coefficients. Namely, we assume: (UE) The diffusion coefficient is uniformly elliptic. There exists

λ 0 > 0 s.t. for (t, x, ξ) ∈ [0, T ] × R d × R d ′ we have λ -1 0 |ξ| 2 a(t, x)ξ, ξ λ 0 |ξ| 2
where a(t, x) := σσ * (t, x), and |.| stands for the Euclidean norm. (SB) The diffusion matrix a is uniformly η-Hölder continuous in space, η > 0, uniformly in time, and the drift b is bounded. That is there exists

L 0 > 0 s.t. sup (t,x)∈[0,T ]×R d |b(t, x)| + sup t∈[0,T ],(x,y)∈R 2d , x =y |a(t, x) -a(t, y)| |x -y| η L 0 .
Throughout the paper we assume that (UE), (SB) are in force.

In the following we will denote by C a generic positive constant that can depend on L 0 , λ 0 , η, d, T . We reserve the notation c for constant depending on L 0 , λ 0 , η, d but not on T . In particular the constants c, C are uniform w.r.t the discretization parameter ∆ = T /N and eventually the value of both c, C may change from line to line.

To establish concentration properties, we will work with the class of Lipschitz continuous functions Denote now by S d-1 the unit sphere of R d . For z ∈ R d \{0}, π S d-1 (z) stands for the uniquely defined projection on S d-1 . For given ρ 0 > 0, β > 0, we introduce the following growth assumption in space for F in the above class of functions:

F : R d → R satisfying |∇F | ∞ = esssup x |∇F (x)
(G ρ0,β ) There exists A ⊂ S d-1 such that ∀y ∈ R d \B(ρ 0 ), π S d-1 (y) ∈ A, y 0 := ρ 0 π S d-1 (y), F (y) -F (y 0 ) β|y -y 0 |,
with A of non empty interior and |A| ε > 0 for d 2 (|.| standing here for the Lebesgue measure of S d-1 ), and A ⊂ {-1, 1} for d = 1. In the above equation B(ρ 0 ) stands for the Euclidean ball of R d of radius ρ 0 , and ., . denotes the scalar product in R d .

Remark 1.1. The above assumption simply means that for |y| ρ 0 the graph of F stays above a given hyperplane. In particular, for all z ∈ A, F (rz) → r→+∞ +∞.

The bounds of the quantities E MC (M, ∆) will be established for real valued functions f that are uniformly Lipschitz continuous in space and measurable bounded in time, such that for a fixed T , F (.) := f (T, .) will be Lispchitz continuous satisfying |∇F | ∞ 1. Moreover, for the lower bounds, we will suppose that the above F satisfies (G ρ0,β ).

Results

Let us first justify that under the assumptions (UE), (SB), the discretization schemes admit a density. For all

x ∈ R d , 0 j < j ′ N , A ∈ B(R d ) (where B(R d ) stands for the Borel σ-field of R d ) we get P X ∆ t j ′ ∈ A X ∆ tj = x = (R d ) j ′ -j-1 ×A p ∆ (t j , t j+1 , x, x j+1 )p ∆ (t j+1 , t j+2 , x j+1 , x j+2 ) • • • p ∆ (t j ′ -1 , t j ′ , x j ′ -1 , x j ′ ) × dx j+1 dx j+2 • • • dx j ′ , (2.1)
where the notation p

∆ (t i , t i+1 , x i , x i+1 ), i ∈ [[0, N -1]] stands in case (a)
for the density at point x i+1 of a Gaussian random variable with mean x i +b(t i , x i )∆ and non degenerated covariance matrix a(t i , x i )∆, whereas in case (b) it stands for the density of a Gaussian random variable with mean

x 1,d ′ i + b 1 (t i , x i )∆, x d ′ +1,d i + x 1,d ′ i ∆ + b 1 (t i , x i )∆ 2 /2
and non degenerated as well covariance matrix

a(t i , x i )∆ a(t i , x i )∆ 2 /2 a(t i , x i )∆ 2 /2 a(t i , x i )∆ 3 /3
, where ∀y ∈ R d , y 1,d ′ = (y 1 , . . . , y d ′ ) * and y d ′ +1,d = (y d ′ +1 , . . . , y d ) * . Equation (2.1) therefore guarantees the existence of the density for the discretization schemes. From now on, we denote by p ∆ (t j , t j ′ , x, •) the transition densities between times t j and t j ′ , 0 j < j ′ N , of the discretization schemes (1.3), (1.4). Let us denote by P x (resp. P tj ,x , 0 j < N ) the conditional probability given X ∆ 0 = x (resp. {X ∆ tj = x}), so that in particular P x X ∆ T ∈ A = A p ∆ (0, T, x, x ′ )dx ′ . We have the following Gaussian estimates for the densities of the schemes. Theorem 2.1 ("Aronson" Gaussian estimates for the discrete Euler scheme). Assume (UE), (SB). There exist constants c > 0, C 1, s.t. for every 0 j < j ′ N :

C -1 p c -1 (t j ′ -t j , x, x ′ ) p ∆ (t j , t j ′ , x, x ′ ) Cp c (t j ′ -t j , x, x ′ ), (2.2)
where for all 0 s < t T , in case (a),

p c (t -s, x, x ′ ) := c 2π(t-s) d/2 exp -c |x ′ -x| 2 2(t-s)
and in case (b)

p c (t-s, x, x ′ ) := √ 3c 2π(t -s) 2 d/2 exp -c |(x ′ ) 1,d ′ -x 1,d ′ | 2 4(t -s) +3 |(x ′ ) d ′ +1,d -x d ′ +1,d -x 1,d ′ +(x ′ ) 1,d ′ 2 (t -s)| 2 (t -s) 3 .
Note that p c enjoys the semigroup property, i.e. ∀0 < s < t, R d p c (ts, x, u)p c (s, u, x ′ )du = p c (t, x, x ′ ) (see Kolmogorov [START_REF] Kolmogorov | Zufällige Bewegungen (zur Theorie der Brownschen Bewegung)[END_REF] or [START_REF] Konakov | Explicit parametrix and local limit theorems for some degenerate diffusion processes[END_REF] for case (b)).

Remark 2.1. The above upper bound can be found in [START_REF] Konakov | Edgeworth type expansions for euler schemes for stochastic differential equations[END_REF] in the case of time homogeneous Lipschitz continuous coefficients. Both bounds can be derived (for time dependent coefficients) from the work of Gobet and Labart [START_REF] Gobet | Sharp estimates for the convergence of the density of the euler scheme in small time[END_REF] under stronger smoothness assumptions. Here, our framework is the one of the "standard" PDE assumptions to derive Aronson's estimates for the fundamental solution of non degenerated non-divergence form second order operators, see e.g. Sheu [START_REF] Sheu | Some estimates of the transition density of a nondegenerate diffusion Markov process[END_REF] or [START_REF] Delarue | Density estimates for a random noise propagating through a chain of differential equations[END_REF]. In particular no regularity in time is needed.

Our second result is the Gaussian concentration of the Monte Carlo error E MC (M, ∆) defined in (1.5) for a fixed M uniformly in ∆ = T /N, N 1.

Theorem 2.2 (Gaussian concentration). Assume (UE), (SB). For the constants c and C of Theorem 2.1, we have for every ∆ = T /N, N 1, and every Lipschitz continuous function in space and measurable bounded in time f :

R d → R satisfying |∇f (T, .)| ∞ 1 in (1.5), ∀r > 0, ∀M 1, P x E MC (M, ∆) r + δ C,α(T ) 2e -M α(T ) r 2 , (2.3) with 1 α(T ) =    c 2T in case (a), c 2T 1 + 3 T 2 1 -1 + T 2 3 + T 4 9 in case (b), (2.4) and δ C,α(T ) = 2 α(T ) log C.
Moreover, if F (.) := f (T, .) 0 satisfies for a given ρ 0 > 0 and β > 0, the growth assumption (G ρ0,β ),

∀r > 0, ∀M 1, P x E MC (M, ∆) r -δc,C,T,f 2 exp - M ᾱ(T ) r β ∨ ρ 0 2 , (2.5)
where δc,C,T,f = (1 + √ 2) α(T ) log C + γ c -1 ,T (F ) + ρ 0 β -F , γ c -1 ,T (dx ′ ) = p c -1 (T, x, x ′ )dx ′ , and F := inf s∈S d-1 F (sρ 0 ). The constant ᾱ(T ) -1 appearing in (2.5) writes in case (a)

1 ᾱ(T ) = Λ + χ :=    c -1 2T + 1 ρ 2 0 log π d/2 C K(d,A) + for d even, c -1 θ 2T + 1 ρ 2 0 log π d/2 C arccos(θ -1/2 )K(d,A) + for d odd, θ ∈ (1, +∞),
where for all d ∈ N * , A ⊂ S d-1 appearing in (G ρ0,β ),

K(d, A) =    |A|(d/2-1)! 2 , d even, |A| Q d-1 2 j=1 (j-1/2) π 1/2
, d odd.

(2.6)

In case (b), d is even and

1 ᾱ(T ) = Λ + χ := c -1 2T 1 + 3 T 2 1 + 1 + T 2 3 + T 4 9 + 1 ρ 2 0 log   π T d/2 [T 2 + 3(1 + 1 + T 2 3 + T 4 9 )] d/2 C K(d, A)   + .
From Theorem 2.1 and our current assumptions on f , we can deduce from the central limit theorem that

M 1/2 E MC (M, ∆) (law) → M N (0, σ 2 (f, ∆)), σ 2 (f, ∆) := E x [f (X ∆ T ) 2 ] -E x [f (X ∆ T )] 2 .
From this asymptotic regime, we thus derive that for large M the typical deviation rate r (i.e. the size of the confidence interval) in (2.3) has order cσ(f, ∆)M -1/2 where for a given threshold α ∈ (0, 1), c := c(α) can be deduced from the inverse of the Gaussian distribution function. In other words, r is typically small for large M . On the other hand, we have a systematic bias δ C,α(T ) , independently of M . In whole generality, this bias is inherent to the concentration arguments used to derive the above bounds, see Section 3, and cannot be avoided. Hence, those bounds turn out to be particularly relevant to derive non asymptotic confidence intervals when r and δ C,α(T ) have the same order. In particular, the parameter M is not meant to go to infinity. This kind of result can be useful if for instance it is particularly heavy to simulate the underlying Euler scheme and that only a relatively small number M of samples is reasonably allowed. On the other hand, the smaller T is the bigger M can be. Precisely, one can prove that the constant C of Theorem 2.1 is bounded by c exp(L 0 T ) (see Section 4). Hence from (2.4), we have δ C,α(T ) = O(T ) for T small. Remark 2.2. For the lower bound, the "expected" value for ᾱ(T ) -1 would be λ corresponding to the largest eigenvalue of one half the inverse of the covariance matrix of the random variable with density p c -1 (T, x, .) appearing in the lower bound of Theorem 2.1. There are two corrections with respect to this intuitive approach. First, there is in case (a) an additional multiplicative term θ > 1 (that can be optimized) when d is odd. This correction is unavoidable for d = 1, anyhow for odd d > 1, it can be avoided up to an additional additive factor like the above χ (see the proof of Proposition 3.3 for details). We kept this presentation to be homogeneous for all odd dimensions. Also, an additive correction (or penalty) factor χ appears. It is mainly due to our growth assumption (G ρ0,β ).Observe anyhow that, for given T > 0, C 1, ε > 0 s.t. |A| ε, if the dimension d is large enough, by definition of K(d, A), we have χ = 0. Still, for d = 1 (which can only occur in case (a)) we cannot avoid the correction factor χ.

Remark 2.3. Let us also specify that in the above definition of χ, ρ 0 is not meant to go to zero, even though some useful functions like |.| satisfy (G ρ0,1 ) with any ρ 0 > 0. Actually, the bound is particularly relevant in 'large regimes", that is when r/β is not assumed to be small. Also, we could replace in the above definition of χ, ρ 0 by R > 0 as soon as r/β R. In particular, if F satisfies (G ρ0,β ), for R ρ 0 it also satisfies (G R,β ). We gave the statement with ρ 0 in order to be uniform w.r.t. the threshold ρ 0 appearing in the growth assumption of F but the correction term can be improved in function of the deviation factor r/β. Remark 2.4. Note that under (UE), (SB), in case (a), the martingale problem in the sense of Stroock and Varadhan is well posed for equation (1.1), see Theorem 7.2.1 in [START_REF] Stroock | Multidimensional diffusion processes[END_REF]. Also, from Theorem 2.1 and the estimates of Section 4, one can deduce that the unique weak solution of the martingale problem has a smooth density that satisfies Aronson like bounds. Furthermore, a careful reading of [START_REF] Konakov | Edgeworth type expansions for euler schemes for stochastic differential equations[END_REF] emphasizes that the discretization error analysis carried therein can be extended to our current framework, i.e. we only need boundedness of the drift and uniform spatial Hölder continuity of the (non-degenerated) diffusion coefficient to control E D (∆). Hence, the above concentration result gives that in case (a), one can control the global error

E(M, ∆) := E D (∆) + E MC (M,

∆). The well-posedness of the martingale problem in case (b) remains to our best knowledge an open question and will concern further research.

Remark 2.5. In case (b), the concentration regime in the above bounds highly depends on T . Since the two components do not have the same scale we have that, in short time, the concentration regime is the one of the non degenerated component in the upper bound (resp. of the degenerated component in the lower bound). For large T , it is the contrary.

We now consider an important case for applications in case (b). Namely, in kinetic models (resp. in financial mathematics) it is often useful to evaluate the expectation of functions that involve the difference of the first component and its normalized average (which corresponds to a time normalization of the second component). This allows to compare the velocity (resp. the price) at a given time T and the averaged velocity (resp. averaged price) on the associated time interval. Obviously, the normalization is made so that the two components have time-homogeneous scales. We have the following result.

Corollary 2.1. In case (b), if f in (1.5) writes f (T, x) = g(T, T -1 T x) where T -1 T = I d ′ ×d ′ 0 d ′ ×d ′ 0 d ′ ×d ′ T -1 I d ′ ×d ′
and g is a Lipschitz continuous function in space and measurable bounded in time satisfying |∇g(T, .)| ∞ 1 then we have for every ∆ = T /N, N 1,

∀M 1, P x E MC (M, ∆) r + δ C,α(T ) 2e -M α(T ) r 2 , with α(T ) = (4 - √ 13) c T and δ C,α(T ) = 2 α(T ) log C.
A lower bound could be derived similarly to Theorem 2.2. The proof of Theorems 2.1 and 2.2 (as well as Corollary 2.1) are respectively postponed to Sections 4.2 and 3.3.

Gaussian concentration and non asymptotic Monte Carlo bounds

Gaussian concentration -Upper bound

We recall that a probability measure γ on R d satisfies a logarithmic Sobolev inequality with constant α > 0 if for all f ∈ H 1 (dγ

) := {g ∈ L 2 (dγ) : |∇g| 2 dγ < +∞} such that f 0, one has Ent γ (f 2 ) α |∇f | 2 dγ, (LSI α )
where Ent γ (φ) = φ log(φ)dγφdγ log φdγ denotes the entropy of the measure γ. In particular, we have the following result (see [START_REF] Ledoux | Concentration of measure and logarithmic sobolev inequalities[END_REF] Section 2.2 eq. (2.17)).

Proposition 3.1. Let V be a C 2 convex function on R d with HessV λI d×d , λ > 0 and such that e -V is integrable with respect to the Lebesgue measure. Let γ(dx) = 1 Z e -V (x) dx be a probability measure (Gibbs measure). Then γ satisfies a logarithmic Sobolev inequality with constant α = 2 λ . Throughout this section we consider a probability measure µ with density m with respect to the Lebesgue measure λ K on R K (here we have in mind K = d or K = M d, M being the number of Monte Carlo paths). We assume that µ is dominated by a probability measure γ in the following sense γ(dx) = q(x)dx satisfies (LSI α ) and ∃κ 1, ∀x ∈ R K , m(x) κq(x).

(H κ,α ) Proposition 3.2. Assume that µ and γ satisfy (H κ,α ). Then for all Lipschitz continuous function

F : R K → R s.t. |∇F | ∞ 1, ∀r > 0, P µ F (Y ) -µ(F ) r + W 1 (µ, γ) κe -r 2 α ,
where

W 1 (µ, γ) = sup |∇F | ∞ 1 |µ(F ) -γ(F )| (Wasserstein distance W 1 between µ and γ).
Proof. By the Markov inequality, one has for every λ > 0,

P µ F (Y ) -µ(F ) r + W 1 (µ, γ) e -λ(µ(F )+r+W1(µ,γ)) E µ e λF (Y ) , (3.1)
and by (H κ,α ), E µ e λF (Y ) κE γ e λF (Y ) . Since γ satisfies a logarithmic Sobolev inequality with constant α > 0, the Herbst argument (see e.g. Ledoux [START_REF] Ledoux | Concentration of measure and logarithmic sobolev inequalities[END_REF] section 2.3) gives

E γ e λF (Y ) e λγ(F )+ α 4 λ 2 , so that E µ e λF (Y ) κe λµ(F )+ α 4 λ 2 +λ(γ(F )-µ(F ))
and

E µ e λF (Y ) κe λµ(F )+ α 4 λ 2 +λW1(µ,γ) , (3.2)
since owing to the definition of W 1 one has W 1 (µ, γ) γ(F )µ(F ). Plugging the above control (3.2) into (3.1) yields

P µ F (Y ) -µ(F ) r + W 1 (µ, γ) κe -λr+ α 4 λ 2 .
An optimization on λ gives the result.

Lemma 3.1. Assume that µ with density m and γ with density q satisfy the domination condition

∃κ 1, ∀x ∈ R d , m(x) κq(x)
and that there exist (α,

β 1 , β 2 ) ∈ (R + ) 3 such that for all Lipschitz continuous function F satisfying |∇F | ∞ 1 and for all λ > 0, E γ e λF (Y ) e λγ(F )+ α 4 λ 2 +β1λ+β2 . Then we have, W 1 (µ, γ) β 1 + α (β 2 + log(κ)).
Proof. Recall first that for a non-negative function f , we have the following variational formulation of the entropy:

Ent γ (f ) = sup E γ [f h] ; E γ e h 1 . (3.3) W.l.o.g. we consider F such that µ(F ) γ(F ). Let λ > 0 and h := λF -λγ(F ) -α 4 λ 2 -β 1 λ -β 2 so that E γ e h
1 and

E µ [h] = E γ m q h = λ (µ(F ) -γ(F )) - α 4 λ 2 -β 1 λ -β 2 .
We then have

µ(F ) -γ(F ) = α 4 λ + β 1 + 1 λ β 2 + E γ m q h , (3.3) α 4 λ + β 1 + 1 λ β 2 + Ent γ m q .
An optimization in λ yields

µ(F ) -γ(F ) β 1 + α β 2 + Ent γ m q (3.4)
Now using the domination condition, one has Ent γ m q = m q log m q dγ log(κ) and the results follows.

Remark 3.1. Note that if γ satisfies an (LSI α ) we have β 1 = β 2 = 0 and the result (3.4) of Lemma 3.1

reads W 1 (µ, γ) α Ent γ m q
α log(κ). For similar controls concerning the W 2 Wasserstein distance see Theorem 1 of Otto and Villani [START_REF] Otto | Generalization of an inequality by Talagrand and links with the logarithmic Sobolev inequality[END_REF] or Bobkov et al. [START_REF] Sergey | Hypercontractivity of Hamilton-Jacobi equations[END_REF].

Using the tensorization property of the logarithmic Sobolev inequality we derive the following Corollary. Note that the term δ κ,α can be seen as a penalty term due on the one hand to the transport between µ and γ, and on the other hand to the explosion of the domination constant κ M between µ ⊗M and γ ⊗M when M tends to infinity. We emphasize that the bias δ κ,α is independent of M . Hence, the result below is especially relevant when r and δ κ,α have the same order. In particular, the non-asymptotic confidence interval given by (3.5) cannot be compared to the asymptotic confidence interval deriving from the central limit theorem whose size has order O(M -1/2 ).

Corollary 3.1. Let Y 1 , . . . , Y M be i.i.d. R d -valued random variables with law µ. Assume there exist α > 0, κ 1 and γ such that (H κ,α ) holds on R d . Then, for all Lipschitz continuous function f :

R d → R satisfying |∇f | ∞ 1, we have ∀r > 0, M 1, P 1 M M f (Y k ) -E f (Y 1 ) r + δ κ,α 2e -M r 2 α , (3.5)
with δ κ,α = 2 α log(κ) 0.

Proof. Let r > 0 and M 1. Clearly, changing f into -f , it suffices to prove that

P 1 M M k=1 f (Y k ) -E µ f (Y 1 ) r + δ κ,α e -M r 2 α .
By tensorization, the measure γ ⊗M satisfies an (LSI α ) with the same constant α as γ, and then the probabilities

µ ⊗M and γ ⊗M satisfy (H κ M ,α ) on R K , K = M d. In this case, Lemma 3.1 gives √ M δ κ,α W 1 (µ ⊗M , γ ⊗M ) + M α log(κ) and then P 1 M M k=1 f (Y k ) -E µ f (Y 1 ) r + δ κ,α P 1 √ M M k=1 f (Y k ) - √ M E µ f (Y 1 ) √ M r + α log(κ) + W 1 (µ ⊗M , γ ⊗M ) .
Applying Proposition 3.2 with the measures µ ⊗M and γ ⊗M , the function

F (x 1 , . . . , x M ) = 1 √ M M k=1 f (x k ) (which satisfies |∇F | ∞ 1) and r = √ M (r + α log(κ)) we obtain P 1 M M k=1 f (Y k ) -E µ f (Y 1 ) r + δ κ,α κ M e -M ( r+ √ α log(κ) ) 2 α ,
and we easily conclude.

Remark 3.2. Note that to obtain the non-asymptotic bounds of the Monte Carlo procedure (3.5), we successively used the concentration properties of the reference measure γ, the control of the distance W 1 (µ, γ) given by the variational formulation of the entropy (see Lemma 3.1) and the tensorization property of the functional inequality satisfied by γ. The same arguments can therefore be applied to a reference measure γ satisfying a Poincaré inequality.

Gaussian concentration -Lower bound

Concerning the previous deviation rate of Proposition 3.2, a natural question consists in understanding whether it is sharp or not. Namely, for a given function f satisfying suitable growth conditions at infinity, otherwise we cannot see the asymptotic growth, do we have a lower bound of the same order, i.e. with Gaussian decay at infinity? The next proposition gives a positive answer to that question.

Proposition 3.3. Let f : R d → R + be a Lipschitz continuous function satisfying |∇f | ∞ 1 and assumption (G ρ0,β ) for given ρ 0 , β > 0.
For a C 2 function V on R d such that e -V is integrable with respect to λ d and s.t. ∃ λ 1, λI d×d Hess(V ) 0, let γ(dx) = e -V (x) Z -1 dx be the associated Gibbs probability measure. We assume that ∃κ 1 s.t. for |x| ρ 0 the measures µ(dx) = m(x)dx and γ(dx) satisfy

m(x) κ -1 e -V (x) Z -1 . Let Λ := λ 2 + sup s∈S d-1 |V (sρ0)| ρ 2 0 + sup s∈S d-1 |∇V (sρ0)| ρ0
.

We have

∀r > 0, P µ [f (Y ) -µ(f ) r -(W 1 (µ, γ) + δ(f, γ))]        K(d,A) Z Λd/2 κ exp -Λ r β ∨ ρ 0 2 , d even, arccos(θ -1/2 )K(d,A) Z Λd/2 κ exp -θ Λ r β ∨ ρ 0 2 , ∀θ > 1, d odd, with δ(f, γ) = γ(f ) + βρ 0 -f , f := inf s∈S d-1 f (sρ 0 ), and K(d, A) defined in (2.6) where A ⊂ S d-1 appears in (G ρ0,β ). Proof. Set E := {Y ∈ A × [ρ 0 , ∞)}.
Here we use the convention that for

d = 1, A × [ρ 0 , +∞) ⊂ (-∞, -ρ 0 ] ∪ [ρ 0 , +∞). Write now P µ [f (Y ) -µ(f ) r -(W 1 (µ, γ) + δ(f, γ))] P µ [f (Y ) -µ(f ) r -(W 1 (µ, γ) + δ(f, γ)), E] κ -1 P γ f (Y ) r -βρ 0 + f , E := κ -1 P. (3.6) Denoting Y 0 = ρ 0 π S d-1 (Y )
, we have

P P γ f (Y 0 ) + f (Y ) -f (Y 0 ) r -βρ 0 + f , E , (G ρ 0 ,β ) P γ β|Y -Y 0 | r -βρ 0 + f -f (Y 0 ), E P γ |Y | r -βρ 0 β + |Y 0 |, E P γ |Y | r β ∨ ρ 0 , π S d-1 (Y ) ∈ A . (3.7) Write P A σ(ds) +∞ ρ0∨ r β dρρ d-1 exp(-V (sρ))Z -1 ,
where σ(ds) stands for the Lebesgue measure of S d-1 . Now, Hess(V ) λI d×d yields ∀ρ

ρ 0 ∨ r β , |V (sρ)|/ρ 2 Λ, Λ := λ 2 + sup s∈S d-1 |V (sρ0)| ρ 2 0 + sup s∈S d-1 |∇V (sρ0)| ρ0
and therefore

P |A| +∞ ρ0∨ r β dρρ d-1 exp(-Λρ 2 )Z -1 |A| Z(2 Λ) d/2 +∞ (ρ0∨ r β )(2 Λ) 1/2 dρρ d-1 exp(- ρ 2 2 ) = |A| Z(2 Λ) d/2 Q d (ρ 0 ∨ r β )(2 Λ) 1/2 . (3.8)
We now have the following explicit expression:

∀x > 0, Q d (x) := exp(- x 2 2 )M (d, x), M (d, x) :=                d 2 -1 i=0 x 2i d 2 -1 j=i+1 2j, d even, d-1 2 -1 i=0 x 2i+1 d-1 2 -1 j=i (2j + 1) + d-1 2 -1 j=0 (2j + 1) exp( x 2 2 ) +∞ x exp(- ρ 2 2 )dρ, d odd,
with the convention that

-1 i=0 = 0, ∀k ∈ N, k-1 j=k j = 1. Observe now that ∞ x exp(-ρ 2 /2)dρ = (2π) 1/2 P[N (0, 1) x] (2π) 1/2 P[Y ∈ K, |Y | x/ cos( θ)] := (2π) 1/2 Q(x)
, where Y ∼ N (0 2×1 , I 2×2 ) is a standard bidimensional Gaussian vector and K := {z ∈ R 2 , z, e 1 cos( θ)|z|}, θ ∈ (0, π 2 ), e 1 = (1, 0). Since Q(x) = θ π exp(-x 2 2 cos 2 ( θ) ), we derive that

Q d (x)    2 d/2-1 (d/2 -1)! exp(-x 2
2 ), d even,

θ2 d/2 π 1/2 d-1 2 j=1 (j -1 2 ) exp(-x 2 2 cos 2 ( θ) ), d odd,
which plugged into (3.8) yields:

P        K(d,A) Z Λd/2 exp -Λ r β ∨ ρ 0 2 , d even, θK(d,A) Z Λd/2 exp -Λ cos 2 ( θ) r β ∨ ρ 0 2 , d odd.
Corollary 3.2. Under the assumptions of Proposition 3.3, let Y 1 , • • • , Y M be i.i.d. R d -valued random variables with law µ. We have ∀r > 0, ∀M 1,

P 1 M M k=1 f (Y k ) -E µ f (Y 1 ) r -(W 1 (µ, γ) + δ(f, γ)) 2 ×        exp -M ( Λ + χ) r β ∨ ρ 0 2 , χ = 1 ρ 2 0 log Z Λd/2 κ K(d,A) + , d even, exp -M (θ Λ + χ) r β ∨ ρ 0 2 , χ = 1 ρ 2 0 log Z Λd/2 κ K(d,A) arccos(θ -1/2 ) + , θ ∈ (1, +∞), d odd, with K(d, A) defined in (2.6).
Proof. We only consider d even. By independence of the ((

Y ) k ) k∈[[1,M]] , exploiting M k=1 {f (Y k )-E µ f (Y 1 ) r -(W 1 (µ, γ) + δ(f, γ))} ⊂ { 1 M M k=1 f (Y k ) -E µ f (Y 1 ) r -(W 1 (µ, γ) + δ(f, γ))}, we have P 1 M M k=1 f (Y k ) -E µ f (Y 1 ) r -(W 1 (µ, γ) + δ(f, γ)) K(d, A) Z Λd/2 κ M exp -M Λ r β ∨ ρ 0 2 . For χ = 1 ρ 2 0 log Z Λd/2 κ K(d,A) +
, we thus obtain

P 1 M M k=1 f (Y k ) -E µ f (Y 1 ) r -(W 1 (µ, γ) + δ(f, γ)) exp -M ( Λ + χ) r β ∨ ρ 0 2
, which completes the proof. .

Proofs of Theorem 2.2 and Corollary 2.1

-Theorem 2.2 -Upper bound (2.3).

In case (a), the Gaussian probability γ c,T with density p c (T, x, .) defined in Theorem 2.1 satisfies a logarithmic Sobolev inequality with constant α(T ) = 2T c . The result then follows from Theorem 2.1 and Corollary 3.1.

In case (b), γ c,T (dx ′ ) = p c (T, x, x ′ )dx ′ = Z -1 e -VT,x(x ′ ) dx ′ where

V T,x (x ′ ) = c |(x ′ ) 1,d ′ -x 1,d ′ | 2 4T + 3 |(x ′ ) d ′ +1,d -x d ′ +1,d -x 1,d ′ +(x ′ ) 1,d ′ 2 T | 2 T 3 .
(3.9)

The Hessian matrix of V T,x satisfies

∀x ′ ∈ R d , HessV T,x (x ′ ) = 2c T I d ′ ×d ′ -3c T 2 I d ′ ×d ′ -3c T 2 I d ′ ×d ′ 6c T 3 I d ′ ×d ′ λI d×d , with λ = c T + 3c T 3 1 -1 + T 2 3 + T 4 9
> 0. By Proposition 3.1, the probability γ c,T satisfies a logarithmic Sobolev inequality with constant α(T

) = 2T c 1 1+ 3 T 2 " 1- q 1+ T 2 3 + T 4 9
« . We still conclude by Theorem 2.1 and Corollary 3.1. -Theorem 2.2 -Lower bound (2.5).

With the notation p c -1 (ts, x, x ′ ) = Z -1 e -Vt-s,x(x ′ ) , the Hessian of the potential

V T,x satisfies ∀x ′ ∈ R d , HessV T,x (x ′ ) λI d×d where λ = c -1 T in case (a) and λ = c -1 T + 3c -1 T 3 1 + 1 + T 2 3 + T 4 9 in case (b). Set γ c -1 ,T (dx ′ ) = p c -1 (T, x, x ′ )dx ′ and µ T (dx ′ ) = p ∆ (0, T, x, x ′ )dx ′ .
Since µ T and γ c,T satisfy (H κ,α ) with κ = C and α = α(T ) defined in (2.4), the probability µ T satisfies (3.2), and Lemma 3.1 yields W 1 (µ T , γ c,T ) α(T ) log(C). Now, γ c -1 ,T and γ c,T satisfy (H κ,α ) with κ = C 2 and α = α(T ). We therefore get from Lemma 3.1, W 1 (γ c -1 ,T , γ c,T ) 2α(T ) log(C). Hence,

W 1 (γ c -1 ,T , µ T ) W 1 (µ T , γ c,T ) + W 1 (γ c -1 ,T , γ c,T ) (1 + √
2) α(T ) log(C). Now, by definition of δc,C,T,f we have δc,C,T,f W 1 (γ c -1 ,T , µ T ) + δ(f, γ c -1 ,T ), (δ(f, γ c -1 ,T ) introduced in Proposition 3.3) and Corollary 3.2 yields

P x 1 M M k=1 f (T, (X ∆ T ) k ) -E x f (T, X ∆ T ) r -δc,C,T,f exp - M ᾱ(T ) r β ∨ ρ 0 2 ,
where observing that for our Gaussian bounds Λ = λ 2 , and andK(d, A) defined in (2.6).

1 ᾱ(T ) =    λ 2 + χ, χ = 1 ρ 2 0 log 2 -d/2 Z λd/2 C K(d,A) + for d even, θ λ 2 + χ, χ = 1 ρ 2 0 log 2 -d/2 Z λd/2 C K(d,A) arccos(θ -1/2 ) + for d odd, θ > 1,
Observe now that in case (a), the normalization factor Z = Z(T, d) associated to p c -1 (T, x, .) writes Z = (2πcT ) d/2 . Hence, recalling that λ = (cT ) -1 , we obtain in this case

χ =    1 ρ 2 0 log π d/2 C K(d,A) + for d even, 1 ρ 2 0 log π d/2 C K(d,A) arccos(θ -1/2 ) for d odd, θ > 1. In case (b), we have Z = (2πc) d/2 T d , λ = 1 cT 1 + 3 T 2 1 + 1 + T 2 3 + T 4 9 so that 2 -d/2 Z λd/2 = π T d/2 [T 2 + 3(1 + 1 + T 2 3 + T 4 9 )] d/2
. Eventually, since in case (b) we always have d even, the correction writes

χ = 1 ρ 2 0 log   π T d/2 [T 2 + 3(1 + 1 + T 2 3 + T 4 9 )] d/2 C K(d, A)   + .
This completes the proof.

-Proof of Corollary 2.1.

Note that the random variable

Y ∆ T = T -1 T X ∆ T admits the density p ∆ Y (T, y, y ′ ) = T d ′ p ∆ (0, T, T T y, T T y ′
) with respect to λ d (dy ′ ). By Theorem 2.1 this density is dominated by (ZT d ′ ) -1 e -V T ,T T y (TT y ′ ) where V T,x is defined in (3.9). The Hessian of y ′ → V T,TT y (T T y ′ ) satisfies

∀y ′ ∈ R d , HessV T,TT y (T T y ′ ) = 2c T I d ′ ×d ′ -3c T I d ′ ×d ′ -3c T I d ′ ×d ′ 6c T I d ′ ×d ′ λI d×d , with λ = c T (4 - √ 13
). We still conclude by Proposition 3.1, and Corollary 3.1.

Derivation of the Gaussian bounds for the discretization schemes

Parametrix representation of the densities

We first derive a parametrix representation of the densities of the schemes. The key idea is to express this density in terms of iterated convolutions of the density of a scheme with frozen coefficients, that therefore admits a Gaussian density, and a suitable kernel, that has an integrable singularity. These representations have previously been obtained in Konakov and Mammen [KM00] and Konakov et al. [START_REF] Konakov | Explicit parametrix and local limit theorems for some degenerate diffusion processes[END_REF].

We first need to introduce some objects and notations. Let us begin with the "frozen" inhomogeneous scheme. For fixed x, x ′ ∈ R d , 0 j < j ′ N , we define

X ∆ ti i∈[[j,j ′ ]] ≡ X ∆,x ′ ti i∈[[j,j ′ ]] by X ∆ tj = x, ∀i ∈ [[j, j ′ ), X ∆ ti+1 = X ∆ ti + b(t i , x ′ )∆ + σ(t i , x ′ )(W ti+1 -W ti ) (4.1)
for case (a). Note that in the above definition the coefficients of the process are frozen at x ′ , but we omit this dependence for notational convenience. In case (b) we define

X ∆ ti i∈[[j,j ′ ]] = X ∆,x ′ ,j ′ ti i∈[[j,j ′ ]] by X ∆ tj = x, and ∀i ∈ [[j, j ′ ), X ∆ ti+1 = X ∆ ti +   b 1 (t i , x ′ )∆ ti+1 ti ( X ∆ s ) 1,d ′ ds   + Bσ t i , x ′ - 0 d ′ ×1 (x ′ ) 1,d ′ (t j ′ -t i ) (W ti+1 -W ti ). (4.2)
That is, in case (b) the frozen process also depends on j ′ through an additional term in the diffusion coefficient. This correction term is needed, in order to have good continuity properties w.r.t. the underlying metric associated to p c when performing differences of the form a(t j , x)a(t j , x ′ -

0 d ′ ×1 (x ′ ) 1,d ′ (t j ′ -t i ))
, see the definition (4.7) and Sections 4.2 and 4.3 for details. From now on, p ∆ (t j , t j ′ , x, •) and p ∆,t j ′ ,x ′ (t j , t j ′ , x, •) denote the transition densities between times t j and t j ′ of the discretization schemes (1.3), (1.4) and the "frozen" schemes (4.1), (4.2) respectively. Let us introduce a discrete "analogue" to the inhomogeneous infinitesimal generators of the continuous objects from which we derive the kernel of the discrete parametrix representation. For a sufficiently smooth function ψ : R d → R and fixed x ′ ∈ R d , j ′ ∈ (0, N ]], define the family of operators (L ∆ tj ) j∈[[0,j ′ ) and (

L ∆ tj ) j∈[[0,j ′ ) = ( L ∆,t j ′ ,x ′ tj ) j∈[[0,j ′ ) by L ∆ tj ψ(x) = E ψ(X ∆ tj +∆ ) X ∆ tj = x -ψ(x) ∆ , and
L ∆ tj ψ(x) = E ψ( X ∆ tj +∆ ) X ∆ tj = x -ψ(x) ∆ .
Using the notation p ∆ (t j , t j ′ , x, x ′ ) = p ∆,t j ′ ,x ′ (t j , t j ′ , x, x ′ ), we now define the discrete kernel H ∆ by

H ∆ (t j , t j ′ , x, x ′ ) = L ∆ tj -L ∆ tj p ∆ (t j + ∆, t j ′ , x, x ′ ), 0 j < j ′ N. (4.3)
Note carefully that the fixed variable x ′ appears here twice: as the final point where we consider the density and as freezing point in the previous schemes (4.1), (4.2). Note also that if j ′ = j + 1 i.e. t j ′ = t j + ∆, the transition probability p ∆,t j ′ ,x ′ (t j+1 , t j+1 , ., x ′ ) is the Dirac measure δ x ′ so that

H ∆ (t j , t j+1 , x, x ′ ) = ∆ -1 E δ x ′ (X ∆ tj+1 ) X ∆ tj = x -E δ x ′ ( X ∆ tj+1 ) X ∆ tj = x , = ∆ -1 p ∆ (t j , t j+1 , x, x ′ ) -p ∆,t j ′ ,x ′ (t j , t j+1 , x, x ′ ) .
From the previous definition (4.3), for all 0 j < j ′ N ,

H ∆ (t j , t j ′ , x, x ′ ) = ∆ -1 R d p ∆ -p ∆,t j ′ ,x ′ (t j , t j+1 , x, u) p ∆,t j ′ ,x ′ (t j+1 , t j ′ , u, x ′ )du.
Analogously to Lemma 3.6 in [START_REF] Konakov | Local limit theorems for transition densities of Markov chains converging to diffusions[END_REF] we obtain the following result.

Proposition 4.1 (Parametrix for the density of the Euler scheme). Assume (UE), (SB) are in force. Then, for 0 t j < t j ′ T ,

p ∆ (t j , t j ′ , x, x ′ ) = j ′ -j r=0 p ∆ ⊗ ∆ H ∆,(r) (t j , t j ′ , x, x ′ ), (4.4)
where the discrete time convolution type operator ⊗ ∆ is defined by

(g ⊗ ∆ f )(t j , t j ′ , x, x ′ ) = j ′ -j-1 k=0 ∆ R d g(t j , t j+k , x, u)f (t j+k , t j ′ , u, x ′ )du,
where g ⊗ ∆ H ∆,(0) = g and for all r 1, H ∆,(r) = H ∆ ⊗ ∆ H ∆,(r-1) denotes the r-fold discrete convolution of the kernel H ∆ . W.r.t. the above definition, we use the convention that p ∆ ⊗ ∆ H ∆,(r) (t j , t j , x, x ′ ) = 0, r 1.

Proof of the Gaussian estimates of Theorem 2.1

The key argument for the proof is given in the following lemma whose proof is postponed to Section 4.3.

Lemma 4.1. There exists c > 0, C 1 s.t. for all 0 j < j ′ N , for all r ∈ [[0, j

′ -j]], ∀(x, x ′ ) ∈ R d , | p ∆ ⊗ ∆ H ∆,(r) (t j , t j ′ , x, x ′ )| C r+1 (t j ′ -t j ) rη/2 r+1 i=1 B 1 + (i -1)η 2 , η 2 p c (t j ′ -t j , x, x ′ ). (4.5)
In the above equation B(m, n) := 1 0 s m-1 (1s) n-1 ds stands for the β function. The upper bound in (2.2) then follows from Proposition 4.1 and the asymptotics of the β function. It is also useful to achieve the first step of the lower bound.

Proof of the lower bound. We provide in this section the global lower bound in short time. W.l.o.g. we assume that T 1. This allows to substitute the constant C appearing in (4.5) by a constant c 0 c exp(|b| ∞ ) uniformly for t j ′t j T . From the upper bound, we derive the lower bound in short time, on the compact sets of the underlying metric, see (4.7) below. This gives the diagonal decay. To get the whole bound in short time it remains to obtain the "off-diagonal" bound. To this end a chaining argument is needed. In case (a) it is quite standard in the Markovian framework, see Chapter VII of Bass [START_REF] Bass | Diffusions and Elliptic Operators[END_REF] or Kusuoka and Stroock [START_REF] Kusuoka | Applications of the Malliavin calculus[END_REF]. In case (b), the chaining in the appendix of [START_REF] Delarue | Density estimates for a random noise propagating through a chain of differential equations[END_REF] can be adapted to our discrete framework. We adapt below these arguments to our non Markovian setting for the sake of completeness.

Eventually, to derive the lower bound for an arbitrary fixed T > 0 it suffices to use the bound in short time and the semigroup property of p c -1 . Naturally, the biggest is T , the worse is the constant in the global lower bound.

From Proposition 4.1 we have

p ∆ (t j , t j ′ , x, x ′ ) p ∆ (t j , t j ′ , x, x ′ ) - j ′ -j r=1 | p ∆ ⊗ ∆ H ∆,(r) (t j , t j ′ , x, x ′ )| c -1 0 p c -1 (t j ′ -t j , x, x ′ ) -c 0 (t j ′ -t j ) η/2 p c (t j ′ -t j , x, x ′ ), (4.6) exploiting p ∆ (t j , t j ′ , x, x ′ ) c -1 0 p c -1 (t j ′ -t j , x, x ′ ) (cf. Lemma 3.1 of [KMM09] in case (b)
) and (4.5) (replacing C by c 0 ) for the last inequality. Equation (4.6) provides a lower bound on compact sets provided that T is small enough. Precisely, denoting

d 2 t j ′ -tj (x, x ′ ) =        |x-x ′ | 2 t j ′ -tj
in case (a),

˛(x ′ ) 1,d ′ -x 1,d ′ ˛2 2(t j ′ -tj ) + 6 ˛(x ′ ) d ′ +1,d -x d ′ +1,d -x 1,d ′ +(x ′ ) 1,d ′ 2 (t j ′ -tj ) ˛2 (t j ′ -tj ) 3
in case (b), (4.7)

we have that, for a given

R 0 1/2, if d 2 t j ′ -tj (x, x ′ ) 2R 0 and (t j ′ -t j ) T c -2 0 exp(-c -1 R 0 )/2 2/η , p ∆ (t j , t j ′ , x, x ′ ) 1 (t j ′ -t j ) S (c -1 0 exp(-c -1 R 0 ) -c 0 T η/2 ) c -1 0 2(t j ′ -t j ) S exp(-c -1 R 0 )
where the parameter S is the intrinsic scale of the scheme. In case (a) S = d/2, in case (b) S = d. Hence, up to a modification of c -1 0 we have that

∃c 0 1, ∀0 j < j ′ N, ∀(x, x ′ ) ∈ (R d ) 2 , d 2 t j ′ -tj (x, x ′ ) 2R 0 , p ∆ (t j , t j ′ , x, x ′ ) c -1 0 (t j ′ -t j ) S . (4.8) In particular ∃c > 0, c 0 1, ∀0 j < j ′ N, ∀(x, x ′ ) ∈ (R d ) 2 , d 2 t j ′ -tj (x, x ′ ) 2R 0 , p ∆ (t j , t j ′ , x, x ′ ) c -1 0 p c -1 (t j ′ -t j , x, x ′ ).
Chaining in case (a). Let us introduce: ∀0 s < t T, (x, x ′ , y)

∈ (R d ) 3 , p ∆,y (s, t, x, x ′ )dx ′ := P[X ∆ t ∈ dx ′ |X ∆ s = x, X ∆ φ(s) = y]
. Equation (4.8) provides a lower bound for the density of the scheme when s, t correspond to discretization times. For the chaining the first step consists in extending this result to arbitrary times 0 s < t T . Precisely, if d 2 t-s (x, x ′ ) R 0 /12 we prove that ∃c 0 1, ∀0 s < t T, ∀y, p ∆,y (s, t, x, x ′ ) c -1 0 (ts) -d/2 . (4.9)

If φ(t) = φ(s), the above density is Gaussian and (4.9) holds. If φ(t) = (φ(s) + ∆), equation (4.9) directly follows from a convolution argument between two Gaussian random variables. Note anyhow carefully that the "crude" convolution argument cannot be iterated L times for an arbitrary large L. Indeed, in that case the constants would have a geometric decay. Thus, for φ(t) -(φ(s) + ∆) ∆ we write

p ∆,y (s, t, x, x ′ ) = (R d ) 2 p ∆,y (s, φ(s) + ∆, x, x 1 )p ∆ (φ(s) + ∆, φ(t), x 1 , x 2 )p ∆ (φ(t), t, x 2 , x ′ )dx 1 dx 2 BR(s,t,x,x ′ ) p ∆,y (s, φ(s) + ∆, x, x 1 )p ∆ (φ(s) + ∆, φ(t), x 1 , x 2 )p ∆ (φ(t), t, x 2 , x ′ )dx 1 dx 2 (4.10) where B R (s, t, x, x ′ ) := {x 1 ∈ R d : d 2 φ(s)+∆-s (x, x 1 ) R} × {x 2 ∈ R d : d 2 t-φ(t) (x 2 , x ′ )
R} for R > 0 to be specified later on. Now, for (x 1 , x 2 ) ∈ B R (s, t, x, x ′ ),

d 2 φ(t)-(φ(s)+∆) (x 1 , x 2 ) = |x 1 -x 2 | 2 φ(t) -(φ(s) + ∆) 2|x 1 -x| 2 + 4|x -x ′ | 2 + 4|x 2 -x ′ | 2 φ(t) -(φ(s) + ∆) 6R + R 0 ,
where we used that for φ(t) -(φ(s) + ∆) ∆,

1 φ(t)-(φ(s)+∆) 3 t-s in the last inequality. Taking R = R 0 /6 we obtain that ∀(x 1 , x 2 ) ∈ B R (s, t, x, x ′ ), d 2 φ(t)-(φ(s)+∆) (x 1 , x 2 ) 2R 0 .
We therefore derive from (4.8) and (4.10) that ∃c 0 > 0,

p ∆,y (s, t, x, x ′ ) c -1 0 (φ(s) + ∆ -s) -d/2 (t -φ(t)) -d/2 (φ(t) -(φ(s) + ∆)) -d/2 (R d ) 2 I (x1,x2)∈BR(s,t,x,x ′ ) dx 1 dx 2 .
Since φ(t) -(φ(s) + ∆) ts and there exists c > 0 s.t.

|{x 1 ∈ R d : d 2 φ(s)+∆-s (x, x 1 ) R}| c(φ(s) + ∆ - s) d/2 , |{x 2 ∈ R d : d 2 t-φ(t) (x 2 , x ′ ) R}| c(t -φ(t)) d/2
where |.| stands for the Lebesgue measure of a given set in R d , we derive (4.9) from the above equation up to a modification of c 0 .

It now remains to do the chaining when for 0 j < j ′ N, (x, x ′ ) ∈ (R d ) 2 we have d 2 t j ′ -tj (x, x ′ ) 2R 0 1. Set L = ⌈Kd 2 t j ′ -tj (x, x ′ )⌉, for K 1 to be specified later on and h := (t j ′t j )/L. Note that L 1. For all i ∈ [[0, L]] we denote s i = t j + ih,

y i = x + i L (x ′ -x) so that s 0 = t j , s L = t j ′ , y 0 = x, y L = x ′ . Introduce now ρ := d t j ′ -tj (x, x ′ )(t j ′ -t j ) 1/2 /L = |x ′ -x|/L and for all i ∈ [[1, L -1]], B i := {x ∈ R d : |x -y i | ρ}. Note that with the previous definitions ∀i ∈ [[0, L -1]], |y i+1 -y i | = |x ′ -x|/L = ρ. Thus, ∀x 1 ∈ B 1 , |x -x 1 | 2ρ, ∀i ∈ [[1, L -2]], (x i , x i+1 ) ∈ B i × B i+1 , |x i -x i+1 | 3ρ, ∀x L-1 ∈ B L-1 , |x L-1 -x ′ | 2ρ.
(4.11)

We can now choose K large enough s.t.

3ρ/ √ h = 3d t j ′ -tj (x, x ′ )/ √ L (R 0 /12) 1/2 (4.12)
so that according to (4.9), denoting

x 0 = x, x L = x ′ , for all i ∈ [[0, L -1]], ∀y ∈ R d , (x i , x i+1 ) ∈ B i × B i+1 , p ∆,y (s i , s i+1 , x i , x i+1 ) c -1 0 h -d/2
(with the slight abuse of notation B 0 = {x}, B L = {x ′ } and p ∆,y (0, h, x 0 , x 1 ) = p ∆ (0, h, x, x 1 )).

We have

p ∆ (t j , t j ′ , x, x ′ ) E tj ,x I ∩ L-1 i=1 X ∆ s i ∈Bi p ∆,X ∆ φ(s L-1 ) (s L-1 , t j ′ , X ∆ sL-1 , x ′ ) . (4.13)
To proceed we have to distinguish two cases: h ∆ and h < ∆.

-If h ∆, write from (4.13),

p ∆ (t j , t j ′ , x, x ′ ) E tj ,x I ∩ L-1 i=1 X ∆ s i ∈Bi E[p ∆,X ∆ φ(s L-1 ) (s L-1 , t j ′ , X ∆ sL-1 , x ′ )|X ∆ sL-1 , X ∆ φ(sL-1) ] .
Since we consider the events X ∆ sL-1 ∈ B L-1 , we derive from (4.11), (4.12) that

|X ∆ sL-1 -x ′ |/ √ h 2ρ/ √ h 3d t j ′ -tj (x, x ′ )/ √ L (R 0 /12) 1/2
. Hence, from (4.9)

p ∆ (t j , t j ′ , x, x ′ ) c -1 0 h -d/2 E tj ,x I ∩ L-1 i=1 X ∆ s i ∈Bi = c -1 0 h -d/2 E tj ,x I ∩ L-2 i=1 X ∆ s i ∈Bi P[X ∆ sL-1 ∈ B L-1 |X ∆ sL-2 , X ∆ φ(sL-2) ] . Now P[X ∆ sL-1 ∈ B L-1 |X ∆ sL-2 , X ∆ φ(sL-2) ] = BL-1 p ∆,X ∆ φ(s L-2 ) (s L-2 , s L-1 , X ∆ sL-2 , y)dy, but since we restrict to X ∆ sL-2 ∈ B L-2 , according to (4.11), we have for all y ∈ B L-1 , |X ∆ sL-2 -y|/ √ h 3ρ/ √ h (R 0 /12) 1/2
for the previous R and therefore (4.9) yields

p ∆ (t j , t j ′ , x, x ′ ) (c -1 0 h -d/2 ) 2 |B L-1 |E tj ,x [I ∩ L-2 i=1 X ∆ s i ∈Bi ].
Iterating the process we finally get

p ∆ (t j , t j ′ , x, x ′ ) (c -1 0 h -d/2 ) L L-1 i=1 |B i |.
Observing that

∃c > 0, ∀i ∈ [[1, L -1]], |B i | cρ d , (4.14)
we obtain from the previous definition of h and L:

p ∆ (t j , t j ′ , x, x ′ ) (c -1 0 h -d/2 ) L (cρ d ) L-1 (4.15) c -1 0 (t j ′ -t j ) -d/2 exp((L -1) log(c -1 0 c(ρ/ √ h) d )) c -1 0 (t j ′ -t j ) -d/2 exp(-cd 2 t j ′ -tj (x, x ′ ))
for a suitable c up to a modification of c 0 . -If h < ∆. We have to introduce for all k ∈ [[j, j ′ ),

I k := {l ∈ [[0, L -1]], s l ∈ [t k , t k+1 [}. Rewrite from (4.13) p ∆ (t j , t j ′ , x, x ′ ) E tj ,x [I ∩ j ′ -1 k=j ∩ i∈I k X ∆ s i ∈Bi p ∆,X ∆ φ(s L-1 ) (s L-1 , t j ′ , X ∆ sL-1 , x ′ )]. Define for all k ∈ [[j, j ′ ), i ∈ [[1, ♯I k ]], I i k ∈ I k and t k s I 1 k < s I 2 k < • • • < s I ♯I k k < t k+1 . In particular, for all i ∈ [[1, ♯I k -1]], s I i+1 k -s I i k = h. Rewrite now, p ∆ (t j , t j ′ , x, x ′ ) E tj ,x [I ∩ j ′ -2 k=j ∩ i∈I k X ∆ s i ∈Bi E[I ∩ i∈I j ′ -1 X ∆ s i ∈Bi p ∆,X ∆ t j ′ -1 (s L-1 , t j ′ , X ∆ sL-1 , x ′ )|F s I ♯I j ′ -2 j ′ -2 ]]. (4.16) Introducing P j ′ -1,j := E[I ∩ i∈I j ′ -1 X ∆ s i ∈Bi p ∆,X ∆ t j ′ -1 (s L-1 , t j ′ , X ∆ sL-1 , x ′ )|F s I ♯I j ′ -2 j ′ -2 ] = E[I X ∆ s I 1 j ′ -1 ∈B I 1 j ′ -1 Q♯I j ′ -1 i=2 B I i j ′ -1 p ∆,X ∆ t j ′ -1 (s I 1 j ′ -1 , s I 2 j ′ -1 , X ∆ s I 1 j ′ -1 , x 2 ) × ♯I j ′ -1 -1 i=2 p ∆,X ∆ t j ′ -1 (s I i j ′ -1 , s I i+1 j ′ -1 , x i , x i+1 )p ∆,X ∆ t j ′ -1 (s I ♯I j ′ -1 j ′ -1 , t j ′ , x ♯I j ′ -1 , x ′ ) ♯I j ′ -1 i=2 dx i |F s I ♯I j ′ -2 j ′ -2 ],
we derive from (4.11), (4.12) and (4.9)

P j ′ -1,j (c -1 0 h -d/2 ) ♯I j ′ -1 ♯I j ′ -1 i=2 |B I i j ′ -1 | B I 1 j ′ -1 p ∆,X ∆ t j ′ -2 (s I ♯I j ′ -2 j ′ -2 , s I 1 j ′ -1 , X ∆ s I ♯I j ′ -2 j ′ -2 , x 1 )dx 1 (4.14) (c -1 0 h -d/2 ) ♯I j ′ -1 (cρ d ) ♯I j ′ -1 -1 B I 1 j ′ -1 p ∆,X ∆ t j ′ -2 (s I ♯I j ′ -2 j ′ -2 , s I 1 j ′ -1 , X ∆ s I ♯I j ′ -2 j ′ -2 , x 1 )dx 1 .
Plugging this estimate in (4.16) we obtain

p ∆ (t j , t j ′ , x, x ′ ) (c -1 0 h -d/2 ) ♯I j ′ -1 (cρ d ) ♯I j ′ -1 -1 ×E tj ,x [I ∩ j ′ -2 k=j ∩ i∈I k X ∆ s i ∈Bi B I 1 j ′ -1 p ∆,X ∆ t j ′ -2 (s I ♯I j ′ -2 j ′ -2 , s I 1 j ′ -1 , X ∆ s I ♯I j ′ -2 j ′ -2 , x 1 )dx 1 ] (c -1 0 h -d/2 ) ♯I j ′ -1 +1 (cρ d ) ♯I j ′ -1 E tj ,x [I ∩ j ′ -2 k=j ∩ i∈I k X ∆ s i ∈Bi ]
using once again (4.12), (4.9) for the last inequality. Iterating this procedure we still obtain (4.15) and can conclude as in the previous case.

Chaining in case (b). If d 2 t-s (x, x ′ ) c-1 R 0 , for c large enough, we derive similarly to case (a) that

∃c 0 > 0, ∀0 s < t T, ∀y, p ∆,y (s, t, x, x ′ ) c -1 0 (t -s) -d . (4.17)
Similarly to the previous paragraph we reduce to the case φ(t) -(φ(s) + ∆) ∆. Then, equation (4.10) still holds and for the previous set B R with the current definition of d 2 . (., .). From standard computations, we derive taking a suitable R that ∀

(x 1 , x 2 ) ∈ B R (s, t, x, x ′ ), d 2 φ(t)-(φ(s)+∆) (x 1 , x 2 ) 2R 0 . Therefore, p ∆,y (s, t, x, x ′ ) c -1 0 (φ(s) + ∆ -s) -d (t -φ(t)) -d (φ(t) -(φ(s) + ∆)) -d (R d ) 2 I (x1,x2)∈BR(s,t,x,x ′ ) dx 1 dx 2 . (4.18) Define now ∀(u, y) ∈ (0, T ] × R d , R > 0, BR (u, y) := {z ∈ R d : |z 1,d ′ -y 1,d ′ | 2 /u R/7, |z d ′ +1,d -y d ′ +1,d -y 1,d ′ u| 2 /u 3 R/24}.
We have that ∀z ∈ BR (u, y):

d 2 u (y, z) := |z 1,d ′ -y 1,d ′ | 2 2u + 6 |z d ′ +1,d -y d ′ +1,d -y 1,d ′ +z 1,d ′ 2 u| 2 u 3 12 |z d ′ +1,d -y d ′ +1,d -y 1,d ′ u| 2 u 3 + 7 |z 1,d ′ -y 1,d ′ | 2 2u R.
Hence BR (φ(s) + ∆s, x) × BR (tφ(t), x ′ ) ⊂ B R (s, t, x, x ′ ) and therefore ∃c > 0, |B R (s, t, x, x ′ )| c(tφ(t)) d (φ(s) + ∆s) d which plugged into (4.18) yields (4.17).

It now remains to do the chaining when d 2 t j ′ -tj (x, x ′ ) 2R 0 . The crucial point is to choose a "good" path between x and x ′ . In the non degenerated case it was naturally the straight line between the two points (Euclidean geodesic). In our current framework we can relate d 2 t j ′ -tj (x, x ′ ) to a deterministic control problem. Introduce:

I(t j ′ -t j , x, x ′ ) = inf{ t j ′ -tj 0 |ϕ(s)| 2 ds, φ(0) = x, φ(t j ′ -t j ) = x ′ }, • φ t = Aφ t + Bϕ t , (CD) with A = 0 d ′ ×d ′ 0 d ′ ×d ′ I d ′ ×d ′ 0 d ′ ×d ′ , B = 1 d ′ ×d ′ 0 d ′ ×d ′ , ϕ ∈ L 2 ([0, t j ′ -t j ], R d ′ ). Problem (CD) is a linear determin-
istic controllability problem that has a unique solution reached for

ϕ s = B * [R(t j ′ -t j , s)] * [Q -1 t j ′ -tj ](x ′ -R(t j ′ -t j , 0)x), (4.19)
where R stands for the resolvent, i.e. ∀0 t, t 0

t j ′ -t j , ∂ t R(t, t 0 ) = AR(t, t 0 ), R(t 0 , t 0 ) = I d×d and Q t j ′ -tj = t j ′ -tj 0 R(t j ′ -t j , s)BB * R(t j ′ -t j ,
s) * ds is the Gram matrix, see e.g. Theorem 1.11 Chapter 1 in Coron [START_REF] Coron | Control and nonlinearity[END_REF]. For (CD) the resolvent writes R(t, t 0 ) =

I d ′ ×d ′ 0 d ′ ×d ′ (t -t 0 )I d ′ ×d ′ I d ′ ×d ′
and therefore the Gram matrix of the control problem corresponds to the covariance matrix of the process

X t = x + t 0 AX s ds + BW t at time t j ′ -t j , that is Q t j ′ -tj = (t j ′ -t j )I d ′ ×d ′ (t j ′ -t j ) 2 /2I d ′ ×d ′ (t j ′ -t j ) 2 /2I d ′ ×d ′ (t j ′ -t j ) 3 /3I d ′ ×d ′
. Hence, explicit computations give:

∀s ∈ [0, t j ′ -t j ], ϕ s = (x ′ ) 1,d ′ -x 1,d ′ (t j ′ -t j ) 2 [6s -2(t j ′ -t j )] + 6 (x ′ ) d ′ +1,d -x d ′ +1,d -(x ′ ) 1,d ′ (t j ′ -t j ) (t j ′ -t j ) 3 [t j ′ -t j -2s],
(4.20) and thus, 1 2 I(t j ′t j , x, x ′ ) = d 2 t j ′ -tj (x, x ′ ) defined in (4.7). Now we have a candidate for a deterministic curve around which we can do the chaining. It is simply the deterministic curve (φ s ) s∈[0,t j ′ -tj ] solution of (CD) for the above control (ϕ s ) s∈[0,t j ′ -tj ] .

To complete the proof of the chaining it remains to specify how to define the (s i ) i 1 , (y i ) i 1 and the associated sets. Recall that 2R 0 1. We set here L := ⌈Kd 2 t j ′ -tj (x, x ′ )⌉ for an integer K 3 to be specified later on. In term of the new distance, L is similar in its definition to the one of the previous paragraph. Define s 0 = 0,

s i := inf{t ∈ [s i-1 , t j ′ -t j ] : t si-1 |ϕ s | 2 ds = I(t j ′ -t j , x, x ′ )/L} ∧ (s i-1 + (t j ′ -t j )/L)I si-1<(t j ′ -tj )(1-2 L ) + (t j ′ - t j )I si-1 (t j ′ -tj )(1-2 L )
, i 1. The previous conditions on R 0 , K give the well posedness of this definition. Lemma 4.2 (Controls on the time step). Set for all i 0, ε i := s i+1s i . There exist a constant c 1 1 and an integer L ∈

[L -1, L/c 1 ], s.t. sL = t j ′ -t j and ∀i ∈ [[0, L -2]], c 1 t j ′ -t j L ε i t j ′ -t j L , t j ′ -t j L ε L-1 2 t j ′ -t j L . (4.21)
Proof. We first set L = inf{k ≥ 1 : s k = t j ′t j }. The set {k ≥ 1 : s k = t j ′t j } is clearly nonempty. The upper bound in (4.21) then follows from the definition of the family (s i ) i 1 . Suppose now that s i < (t j ′t j )(1 -2/L) for a given 0 ≤ i ≤ L -2. Assume also that s i+1s i < (t j ′t j )/L (otherwise ε i = (t j ′t j )/L). Then, si+1 si |ϕ s | 2 ds = I(t j ′t j , x, x ′ )/L. From (4.19), (4.20), we deduce that ∃c 2 > 0, sup

0≤s≤t j ′ -tj |ϕ s | ≤ c 2 (t j ′ -t j ) -1/2 d t j ′ -tj (x, x ′ ).
Hence, we obtain

si+1 si |ϕ s | 2 ds := I(t j ′ -t j , x, x ′ ) L c 2 2 ε i d 2 t j ′ -tj (x, x ′ ) (t j ′ -t j ) .
Recalling that I(t j ′t j , x, x ′ ) = 2d 2 t j ′ -tj (x, x ′ ), the lower bound in (4.21) follows for all i s.t. s i < T (1 -2/L). The bound for L and the last time step are then easily derived.

Define now for all i ∈ [[0, L]], y i = φ si (in particular y 0 = x and y L = x ′ ), and for all i ∈ [[1, L -1]],

B i := {z ∈ R d : |Q -1/2 Kρ 2 (R(s i , s i-1 )y i-1 -z)| + |Q -1/2 Kρ 2 (z -R(s i , s i+1 )y i+1 )| 2R 0 K -1/2 },
where ρ := d t j ′ -tj (x, x ′ )(t j ′t j ) 1/2 /L. Because of the transport term, we are led to consider sets that involve the forward transport from the previous point on the optimal curve and the backward transport of the next point in the above definition. Equation (4.13) still holds with L replaced by L. Following the strategy of the previous paragraph concerning the conditioning, the end of the proof relies on the following Lemma 4.3 (Controls for the chaining). With the previous assumptions and definitions we have that for K large enough: The key estimate is the following control of the convolution kernel H ∆ . There exist c > 0, C 1, s.t. for all 0 j < j ′ N, x, x ′ ∈ R d , |H ∆ (t j , t j ′ , x, x ′ )| C(t j ′t j ) -1+η/2 p c (t j ′t j , x, x ′ ).

∀i ∈ [[1, L -2]], ∀(x i , x i+1 ) ∈ B i × B i+1 , d 2 εi (x i , x i+1 ) 2R 0 , ∀x 1 ∈ B 1 , d 2 s1 (x, x 1 ) 2R 0 , ∀xL -1 ∈ BL -1 , d 2 ε L-1 (xL -1 , x ′ ) 2R 0 .
(4.24) Indeed this bound yields that for all 0 j < j

′ N, x, x ′ ∈ R d | p ∆ ⊗ ∆ H ∆ (t j , t j ′ , x, x ′ )| ∆ j ′ -j-1 k=0 R d p ∆ (t j , t j+k , x, u)|H ∆ (t j+k , t j ′ , u, x ′ )|du C 2 ∆ j ′ -j-1 k=0 (t j ′ -t j+k ) -1+η/2 p c (t j ′ -t j , x, x ′ ) C 2 (t j ′ -t j ) η/2 B 1, η 2 p c (t j ′ -t j , x, x ′ )
using the inequality p ∆ (t j+kt j , x, u) Cp c (t j , t j+k , x, u) (cf. Lemma 3.1 of [START_REF] Konakov | Explicit parametrix and local limit theorems for some degenerate diffusion processes[END_REF] in case (b)) and the semigroup property of p c for the last but one inequality. The bound (4.5) then follows from the above control and (4.24) by induction.

Proof of (4.24) We consider two cases.

j ′ = j + 1. From (4.3) we have in this case, ∀x,

x ′ ∈ R d , H ∆ (t j , t j ′ , x, x ′ ) = ∆ -1 (p ∆ -p ∆ )(t j , t j ′ , x, x ′ )
which are Gaussian densities. In case (a) we have

H ∆ (t j , t j ′ , x, x ′ ) = ∆ -1 G ( √ ∆σ(t j , x)) -1 (x ′ -x -b(t j , x)∆) ∆ d det(a(t j , x)) 1/2 - G ( √ ∆σ(t j , x ′ )) -1 (x ′ -x -b(t j , x ′ )∆) ∆ d det(a(t j , x ′ )) 1/2
, where ∀z ∈ R d , G(z) = exp(-|z| 2 /2)(2π) -d/2 stands for the density of the standard Gaussian vector of R d . In case (b) we get

H ∆ (t j , t j ′ , x, x ′ ) = ∆ -1 (2 √ 3) d ′ ×       G (∆ 1/2 σ(t j , x)) -1 ((x ′ ) 1,d ′ -x 1,d ′ -b 1 (t j , x)∆) 2 √ 3(∆ 3/2 σ(t j , x)) -1 ((x ′ ) d ′ +1,d -x d ′ +1,d -x 1,d ′ +(x ′ ) 1,d ′ 2 ∆) ∆ d det(a(t j , x)) - G (∆ 1/2 σ(t j , (x ′ ) ∆ )) -1 ((x ′ ) 1,d ′ -x 1,d ′ -b 1 (t j , x ′ )∆) 2 √ 3(∆ 3/2 σ(t j , (x ′ ) ∆ )) -1 ((x ′ ) d ′ +1,d -x d ′ +1,d -x 1,d ′ +(x ′ ) 1,d ′ 2 ∆) ∆ d det(a(t j , (x ′ ) ∆ ))       , where (x ′ ) ∆ := x ′ - 0 d ′ ×1 (x ′ ) 1,d ′ ∆
allows to have good continuity properties to equilibrate the singularities coming from the difference |x -

(x ′ ) ∆ | |(x ′ ) 1,d ′ -x 1,d ′ |(1 + ∆ 2 ) + |(x ′ ) d ′ +1,d -x d ′ +1,d -x 1,d ′ +(x ′ ) 1,d ′ 2
∆| with the terms appearing in the exponential. In all cases, tedious but elementary computations involving the mean value theorem yield that ∃c > 0, C 1 s.t.

|H ∆ (t j , t j ′ , x, x ′ )| C∆ -1+η/2 p c (∆, x, x ′ ).

j ′ > j + 1. We write in case (a)

H ∆ (t j , t j ′ , x, x ′ ) = ∆ -1 R d G(z) p ∆ (t j+1 , t j ′ , x + b(t j , x)∆ + ∆ 1/2 σ(t j , x)z, x ′ ) -p ∆ (t j+1 , t j ′ , x, x ′ ) -p ∆ (t j+1 , t j ′ , x + b(t j , x ′ )∆ + ∆ 1/2 σ(t j , x ′ )z, x ′ ) -p ∆ (t j+1 , t j ′ , x, x ′ ) dz := T (a) 1 -T (a) 2 . Now exploiting that R d G(z)zdz = 0, a Taylor expansion at order 3 of T (a) 1 , T (a) 2 yields H ∆ (t j , t j ′ , x, x ′ ) = b(t j , x) -b(t j , x ′ ), D x p ∆ (t j+1 , t j ′ , x, x ′ ) + 1 2 Tr a(t j , x) -a(t j , x ′ ) D 2 x p ∆ (t j+1 , t j ′ , x, x ′ ) +R ∆ (t j , t j ′ , x, x ′ ) := (H + R ∆ )(t j , t j ′ , x, x ′ ). ( 4
R ∆ (t j , t j ′ , x, x ′ ) = ∆ 2 Tr bb * (t j , x) -bb * (t j , x ′ ) D 2 x p ∆ (t j+1 , t j ′ , x, x ′ ) + 3∆ -1 |ν|=3 R d dzG(z) 1 0 dδ(1 -δ) 2 D ν x p ∆ (t j+1 , t j ′ , x + δ(b(t j , x)∆ + σ(t j , x)∆ 1/2 z), x ′ ) (b(t j , x)∆ + σ(t j , x)∆ 1/2 z) ν ν! -D ν x p ∆ (t j+1 , t j ′ , x + δ(b(t j , x ′ )∆ + σ(t j , x ′ )∆ 1/2 z), x ′ ) (b(t j , x ′ )∆ + σ(t j , x ′ )∆ 1/2 z) ν ν!
using the following notations for multi-indices and powers. For ν = (ν 1 , ..., ν

d ) ∈ N d , x = (x 1 , ..., x d ) * set |ν| = ν 1 + ... + ν d , ν! = ν 1 !...ν d !, (x) ν = x ν1 1 ... x ν d d , D ν x = D ν1 x1 ...D ν d x d . Recalling the standard control ∃c > 0, C 1, ∀ν, |ν| 4, ∀0 j < j ′ N, (x, x ′ ) ∈ (R d ) 2 , |D ν x p ∆ (t j , t j ′ , x, x ′ )| C(t j ′ -t j ) -|ν|/2 p c (t j ′ -t j , x, x ′ ) ( 4
.26) for the derivatives of Gaussian densities, we obtain:

|R ∆ (t j , t j ′ , x, x ′ )| C∆|b| 2 ∞ (t j ′ -t j+1 ) -1 p c (t j ′ -t j+1 , x, x ′ ) + 3∆ -1 |ν|=3 R d dzG(z) 1 0 dδ(1 -δ) 2 × D ν x p ∆ (t j+1 , t j ′ , x + δ(b(t j , x)∆ + σ(t j , x)∆ 1/2 z), x ′ ) (b(t j , x)∆ + σ(t j , x)∆ 1/2 z) ν ν! - (b(t j , x ′ )∆ + σ(t j , x ′ )∆ 1/2 z) ν ν! -D ν x p ∆ (t j+1 , t j ′ , x + δ(b(t j , x ′ )∆ + σ(t j , x ′ )∆ 1/2 z), x ′ ) -D ν x p ∆ (t j , t j ′ , x + δ(b(t j , x)∆ + σ(t j , x)∆ 1/2 z), x ′ ) × (b(t j , x ′ )∆ + σ(t j , x ′ )∆ 1/2 z) ν ν! C p c (t j ′ -t j+1 , x, x ′ ) + ∆ 1/2 |x -x ′ | η (t j ′ -t j+1 ) 3/2 R d dzG(z) 1 0 dδ(1 -δ) 2 |z| 3 × p c (t j ′ -t j+1 , x + δ(b(t j , x)∆ + σ(t j , x)∆ 1/2 z), x ′ ) + ∆|x -x ′ | η (t j ′ -t j+1 ) 2 R d dzG(z) [0,1] 2 dδdγ(1 -δ) 2 |z| 4 × p c (t j ′ -t j+1 , x + γδ[(b(t j , x) -b(t j , x ′ ))∆ + (σ(t j , x) -σ(t j , x ′ ))∆ 1/2 z], x ′ ). Now, using the inequality ∀ε ∈ (0, 1), |x -x ′ + ρ| 2 |x -x ′ | 2 (1 -ε) + |ρ| 2 (1 -ε -1 ), ∀ρ ∈ R d , taking ρ = δ(b(t j , x)∆ + σ(t j , x)∆ 1/2 z) and ρ = γδ[(b(t j , x) -b(t j , x ′ ))∆ + (σ(t j , x) -σ(t j , x ′ ))∆ 1/2 z
] respectively in the first and second integral we get

|R ∆ (t j , t j ′ , x, x ′ )| C (1 -ε) d/2 p (1-ε)c (t j ′ -t j+1 , x, x ′ ) 1 + ∆ 1/2 |x -x ′ | η (t j ′ -t j+1 ) 3/2 R d dzG(z) 1 0 dδ(1 -δ) 2 |z| 3 × exp c |σ| 2 ∞ |z| 2 ∆ t j ′ -t j+1 (ε -1 -1) + ∆|x -x ′ | η (t j ′ -t j+1 ) 2 R d dzG(z) [0,1] 2 dδdγ(1 -δ) 2 |z| 4 exp c |σ| 2 ∞ |z| 2 ∆ t j ′ -t j+1 (ε -1 -1) .
Choosing ε sufficiently close to one the above integrals are finite and therefore for different c, C depending on ε as well, we have

|R ∆ (t j , t j ′ , x, x ′ )| Cp c (t j ′ -t j+1 , x, x ′ ) 1 + |x -x ′ | η (t j ′ -t j+1 ) C(t j ′ -t j+1 ) -1+η/2 p c (t j ′ -t j+1 , x, x ′ ). (4.27)
Now with the definitions of (4.25) we also have from (4.26)

|H(t j , t j ′ , x, x ′ )| C(t j ′ -t j+1 ) -1+η/2 p c (t j ′ -t j+1 , x, x ′ ).
Plugging this last estimate and (4.27) in (4.25) we derive

|H ∆ (t j , t j ′ , x, x ′ )| C(t j ′ -t j+1 ) -1+η/2 p c (t j ′ -t j+1 , x, x ′ ) C(t j ′ -t j ) -1+η/2 p c (t j ′ -t j , x, x ′ ).
For case (b) we have

H ∆ (t j , t j ′ , x, x ′ ) = ∆ -1 R d G(z) p ∆ (t j+1 , t j ′ , x ∆ + B ∆ (t j , x) + Σ ∆ (t j , x)z, x ′ ) -p ∆ (t j+1 , t j ′ , x ∆ , x ′ ) -p ∆ (t j+1 , t j ′ , x ∆ + B ∆ (t j , x ′ ) + Σ ∆ (t j , (x ′ ) ∆,j,j ′ ) z, x ′ ) -p ∆ (t j+1 , t j ′ , x ∆ , x ′ ) dz := T (b) 1 -T (b) 2 ,
where

x ∆ := x + 0 d ′ ×1 x 1,d ′ ∆ , (x ′ ) ∆,j,j ′ := x ′ - 0 d ′ ×1 (x ′ ) 1,d ′ (t j ′ -t j ), and we set ∀y ∈ R d , Σ ∆ (t j , y) := ∆ 1/2 σ(t j , y) 0 ∆ 3/2 σ(t j , y)/2 ∆ 3/2 σ(t j , y)/(2 √ 3) , B ∆ (t j , y) := b 1 (t j , y)∆ b 1 (t j , y)∆ 2 /2
. The strategy now relies as in case (a) on Taylor expansions. Actually the d ′ first component in the above expression can be handled exactly in the same way. Anyhow, in order to deal with the d ′ last components we perform the expansion of T

(b) 1 , T (b) 2
around the point x ∆ . We obtain:

H ∆ (t j , t j ′ , x, x ′ ) = b 1 (t j , x) -b 1 (t j , x ′ ), D x 1,d ′ p ∆ (t j+1 , t j ′ , x ∆ , x ′ ) (4.28) + 1 2 Tr (a(t j , x) -a (t j , (x ′ ) ∆,j,j ′ )) D 2 x 1,d ′ p ∆ (t j+1 , t j ′ , x ∆ , x ′ ) +R ∆ (t j , t j ′ , x ∆ , x ′ ) := (H + R ∆ )(t j , t j ′ , x ∆ , x ′ ),
where D x 1,d ′ denotes the differentiation w.r.t. the first d ′ components and similarly to (4.25), H is the difference of the generators at time t j of the processes (X t ) t 0 satisfying (1.1) and the Gaussian process Xt -D x d ′ +1,d p ∆ (t j+1 , t ′ j , x ∆ + γ(B ∆ (t j , x ′ ) + Σ ∆ (t j , (x ′ ) ∆,j,j ′ ) z), x ′ ), B ∆ (t j , x ′ ) + Σ ∆ (t j , (x ′ ) ∆,j,j ′ ) z

= x + t tj b 1 (s, x ′ ) ( Xs ) 1,d ′ ds + t tj Bσ s, x ′ - 0 d ′ ×1 (x ′ ) 1,d ′ (t j ′ -s) dW s , t ∈ [t j , t j ′ ] (continuous version of ( X∆ ti ) i∈[[j,j ′ ]]
d ′ +1,d +    2∆ -1 |θ|=2 R d dzG(z) [0,1] 2
dγdδ(1δ)γ 2 × D θ D x 1,d p ∆ (t j+1 , t j ′ , x ∆ + δγ(B ∆ (t j , x) + Σ ∆ (t j , x)z), x ′ ) (B ∆ (t j , x) + Σ ∆ (t j , x)z) θ θ! (B ∆ (t j , x) + Σ ∆ (t j , x)z) 1,d ′ -D θ D ν x 1,d p ∆ (t j+1 , t j ′ , x ∆ + δγ(B ∆ (t j , x ′ ) + Σ ∆ (t j , (x ′ ) ∆,j,j ′ )z), x ′ ) (B ∆ (t j , x ′ ) + Σ ∆ (t j , (x ′ ) ∆,j,j ′ )z) θ θ! × (B ∆ (t j , x ′ ) + Σ ∆ (t j , (x ′ ) ∆,j,j ′ )z) 1,d ′ := (R 1,∆ + R 2,∆ + R 3,∆ )(t j , t j ′ , x, x ′ ). The above control yields |R 1,∆ (t j , t j ′ , x, x ′ )| Cp c (t j ′t j+1 , x ∆ , x ′ ) Cp c (t j ′t j , x, x ′ ), (4.30) up to a modification of c, C in the last inequality, observing that

Let µ = (µ 1 , • • • , µ d ′ ) ∈ N d ′ , ν = (ν 1 , • • • , ν d ′ ) ∈ N d ′ be
|(x ′ ) d ′ +1,d -x d ′ +1,d ∆ - (x ′ ) 1,d ′ +x 1,d ′ ∆ 2 (t j ′ -t j+1 )| 2 (t j ′ -t j+1 ) 3 = |(x ′ ) d ′ +1,d -x d ′ +1,d -(x ′ ) 1,d ′ +x 1,d ′ 2 (t j ′ -t j+1 ) -x 1,d ′ ∆| 2 (t j ′ -t j+1 ) 3 = |(x ′ ) d ′ +1,d -x d ′ +1,d -(x ′ ) 1,d ′ +x 1,d ′ 2 (t j ′ -t j ) -(x 1,d ′ -(x ′ ) 1,d ′ ) 2 ∆| 2 (t j ′ -t j+1 ) 3 (1 -ε) |(x ′ ) d ′ +1,d -x d ′ +1,d -(x ′ ) 1,d ′ +x 1,d ′ 2 (t j ′ -t j )| 2 (t j ′ -t j ) 3 + (1 -ε -1 ) |(x ′ ) 1,d ′ -x 1,d ′ | 2 ∆ 2 4(t j ′ -t j ) 3
, ∀ε ∈ (0, 1).

Let us now turn to R 2,∆ (t j , t j ′ , x, x ′ ). From (4.29) we obtain:

|R 2,∆ (t j , t j ′ , x, x ′ )| C (1 -ε) d/2 p (1-ε)c (t j ′ -t j+1 , x ∆ , x ′ ) 1 (t j ′ -t j+1 ) 1/2 + ∆ 1/2 |x -(x ′ ) ∆,j,j ′ | η (t j ′ -t j+1 ) 3/2 R d dzG(z)|z|× exp c|σ| 2 ∞ |z| 2 (ε -1 -1) + ∆ 2 |x -(x ′ ) ∆,j,j ′ | η (t j ′ -t j+1 ) 3 R d dzG(z)|z| 2 exp c|σ| 2 ∞ |z| 2 (ε -1 -1) .
Taking ε sufficiently close to 1 the above integrals are finite. Also

|x -(x ′ ) ∆,j,j ′ | η C(|x 1,d ′ -(x ′ ) 1,d ′ | η + |x d ′ +1,d -(x ′ ) d ′ +1,d -(x ′ ) 1,d (t j ′ -t j )| η ) C(|x 1,d ′ -(x ′ ) 1,d ′ | η + |x d ′ +1,d -(x ′ ) d ′ +1,d - x 1,d + (x ′ ) 1,d 2 (t j ′ -t j )| η + |x 1,d ′ -(x ′ ) 1,d ′ | η t j ′ -t j 2 η
),(4.31) and similarly to (4.30),

C

(1-ε) d/2 p (1-ε)c (t j ′t j+1 , x ∆ , x ′ ) Cp c (t j ′t j , x, x ′ ) up to a modification of c, C in the previous r.h.s. Hence, from the definition of p c in case (b) we deduce:

|R 2,∆ (t j , t j ′ , x, x ′ )| C(t j ′t j ) -1+η/2 p c (t j ′t j , x, x ′ ). Thus, proceeding as in case (a) using (4.31), (4.29) to handle R 3,∆ (t j , t j ′ , x, x ′ ), we eventually derive that |R ∆ (t j , t j ′ , x, x ′ )| C(t j ′t j ) -1+η/2 p c (t j ′t j , x, x ′ ). From the definition in equation (4.28) we also obtain |H(t j , t j ′ , x, x ′ )| C(t j ′t j ) -1+η/2 p c (t j ′t j , x, x ′ ) using (4.31), (4.29) and the statement (4.24) follows.

Remark 4.1. Note that the time dependence in the frozen dynamics (4.2) somehow corresponds to the backward transport of the terminal condition. It is crucial in order to allow from (4.31) the compensation of the exploding terms associated to derivatives in x 1,d ′ of order greater than 2 and derivatives in x d ′ +1,d of order greater than 1 appearing in the kernel H ∆ . A similar construction was used in [START_REF] Konakov | Explicit parametrix and local limit theorems for some degenerate diffusion processes[END_REF].

= Q 1 i + Q 2 i + Q 3 i .
One has

Q 1 i c 2 j=1 ε 1/2-j i |(x i -R(s i , s i-1 )y i-1 ) j | c 2 j=1 ε i Kρ 2 1/2-j
(K 1/2 ρ) 1-2j |(x i -R(s i , s i-1 )y i-1 ) j |, denoting for all z ∈ R d , z 1 := z 1,d ′ , z 2 := z d ′ +1,d with a slight abuse of notation. Now, from (4.21), ε i /(Kρ 2 ) c 1 ((t j ′t j )/L)/(Kd 2 t j ′ -tj (x, x ′ )(t j ′t j )/L 2 ) = c 1 L Kd 2 t j ′ -t j (x,x ′ ) . Thus, recalling L = ⌈Kd 2 t j ′ -tj (x, x ′ )⌉, ∃c >

0, ∀j ∈ [[1, 2]], εi Kρ 2
1/2-j c and

Q 1 i c 2 j=1 (K 1/2 ρ) 1-2j |(x i -R(s i , s i-1 )y i-1 ) j | c|Q -1/2
Kρ 2 (x i -R(s i , s i-1 )y i-1 )| cR 0 K -1/2 , exploiting x i ∈ B i for the last identity. The term Q 3 i could be handled in a similar way so that 2R 0 for K large enough independent of t j ′t j . Eventually, for x 1 ∈ B 1 , xL -1 ∈ BL -1 the terms Q 0 := |Q -1/2 ε0 (R(s 1 , 0)xx 1 )| and Q L-1 := |Q -1/2 ε L-1 (R(t j ′t j , sL -1 )xL -1x ′ )|

Q 1 i + Q 3 i cR 0 K -1/2 .

  | 1 where |∇F | denotes the Euclidean norm of the gradient ∇F , defined almost everywhere, of F .

(4. 22 )

 22 For the same c 1 as in Lemma 4.2,∀i ∈ [[1, L -1]], |B i | c 1 ρ 2d , (4.23) where |B i | stands for the Lebesgue measure of the set B i . Indeed, exploiting, (4.21), (4.22) (resp. (4.23)) instead of (4.11), (4.12) (resp. (4.14)), the proof remains unchanged. The proof of Lemma 4.3 is postponed to Section 4.3. 4.3. Proof of the technical Lemmas 4.3.1. Proof of Lemma (4.1).

  multi-indices. Similarly to (4.26) we have,∃c > 0, C 1, ∀(µ, ν), |µ| 3, |ν| 4,∀0 j < j ′ N, (x, x ′ ) ∈ R d × R d , |D ν x 1,d ′ D µ x d ′ +1,d p ∆ (t j , t j ′ , x, x ′ )| C(t j ′t j ) -(|ν|/2+3/2|µ|) p c (t j ′t j , x, x ′ ). (4.29)

4. 4 .

 4 Proof of Lemma 4.3Let us first prove (4.22). We begin with (x i , x i+1 ) ∈ B i × B i+1 , i ∈ [[1, L -2]]. From Section 4.2, one can check that |Q -1/2 εi (R(s i+1 , s i )x ix i+1 )| 2 = 2d 2 εi (x i , x i+1 ).Hence,Q i := d εi (x i , x i+1 ) s i+1 , s i )x ix i+1 )| c|Q -1/2 εi (x i -R(s i , s i+1 )x i+1 )| c |Q -1/2 εi (x i -R(s i , s i-1 )y i-1 )| + |Q -1/2 εi (R(s i , s i-1 )y i-1y i )| + |Q -1/2εi (y i -R(s i , s i+1 )x i+1 )| :

  .25)In the above equation H is the difference of the infinitesimal generators at time t j of the processes (X t ) t 0 satisfying (1.1) and the Gaussian process Xt = x + ′ )dW s , t t j , which can be seen as the continuous version of the frozen Markov chain ( X∆ ti ) i∈[[j,N ]] introduced in (4.1), applied to the Gaussian density p ∆ (t j+1 , t j ′ , •, x ′ ) at point x. The remainder term writes

t tj b(s, x ′ )ds + t tj σ(s, x

  introduced in (4.2)), applied to the Gaussian density p ∆ (t j+1 , t j ′ , •, x ′ ) at point x ∆ . The remainder term writesR ∆ (t j , t j ′ , x, x ′ ) = ∆ 2 Tr b 1 b * 1 (t j , x)b 1 b * 1 (t j , x ′ ) D 2 x 1,d ′ p ∆ (t j+1 , t j ′ , x ∆ , x ′ ) + ∆ 4 Tr b 1 b * 1 (t j , x)b 1 b * 1 (t j , x ′ ) ∆ + a(t j , x)a(t j , (x ′ ) ∆,j,j ′ ) D 2 x 1,d ′ ,x d ′ +1,d p ∆ (t j+1 , t j ′ , x ∆ , x ′ ) dγ D x d ′ +1,d p ∆ (t j+1 , t j ′ , x ∆ + γ(B ∆ (t j , x) + Σ ∆ (t j , x)z), x ′), B ∆ (t j , x) + Σ ∆ (t j , x)z

	∆ -1	dzG(z)	1	d ′ +1,d
	R d	0

  Now Q 2 i := √ 2d εi (y i , y i+1 ) I(s i , s i+1 , y i , y i+1 ) 1/2 c

	Hence, for all i ∈ [[1, L -2]], Q i	si+1 si	|ϕ s | 2 ds	1/2	c	dt j ′ -t j (x,x ′ ) L 1/2	c K 1/2 .

c|Q

-1/2 ε L-1 (xL -1 -R(sL -1 , t j ′t j )x ′ )| can be controlled as the previous

] from the definitions of B 1 , BL -1 , so that Q i 2R 0 , i ∈ {0, L -1} as well. This proves (4.22).

It now remains to control the Lebesgue measure of the sets

The previous definition of

3K 1/2 . Now, arguments similar to those used to control the above

Since the term R 3

i could be handled in the same way we deduce that for K large enough and R 0 large enough w.r.t. the above c, R i 2R 0 K -1/2 . Hence E i ⊂ B i which completes the proof.


Rosana El Jurdi
email: eljurdi.rosana@gmail.com

Lucia Nichelli

Agusti Alentorn

Ghislain Vaillant

Guanghui Fu

Khê Hoang-Xuan

Caroline Houillier

Stéphane Lehéricy

Olivier Colliot

Rosana El

Border irregularity loss for automated segmentation of primary brain lymphomas on post-contrast MRI

Keywords: Lymphoma, Brain, CNS, MRI, segmentation, deep learning, loss function, constraint

Unlike for other brain tumors, there has been little work on the automatic segmentation of primary central nervous system (CNS) lymphomas. This is a challenging task due the highly variable pattern of the tumor and its boundaries. In this work, we propose a new loss function that controls border irregularity for deep learning-based automatic segmentation of primary CNS lymphomas. We introduce a border irregularity loss which is based on the comparison of the segmentation and it smoothed version. The border irregularity loss is combined with a previously proposed topological loss to better control the different connected components. The approach is general and can be used with any segmentation network. We studied a population of 99 patients with primary CNS lymphoma. 40 patients were isolated from the very beginning and formed the independent test set. The segmentations were performed on post-contrast T1-weighted MRI. The MRI were acquired in clinical routine and were highly heterogeneous. The proposed approach substantially outperformed the baseline across the various evaluation metrics (by 6 percent points of Dice, 40mm of Hausdorff distance and 6mm of mean average surface distance). However, the overall performance was moderate, highlighting that automatic segmentation of primary CNS lymphomas is a difficult task, especially when dealing with clinical routine MRI. The code is publicly available here: https://github.com/rosanajurdi/LymphSeg.

INTRODUCTION

Primary central nervous system (CNS) lymphoma (PCNSL) is a rare and aggressive type of cancer that primarily affects the brain and spinal cord. It accounts for approximately 4% of newly diagnosed primary CNS tumors and 1% of all non-Hodgkin lymphomas (NHL). Magnetic resonance imaging (MRI) suggests its diagnosis, is needed for tumor biopsy planification, and has a pivotal role in PCNSL post-treatment assessment. Precise tumor quantification is highly desirable, and one first needed step to promote better patient care. It could aid in pre-surgical planning and in objective tumor response evaluation. [START_REF] Razek | Clinical applications of artificial intelligence and radiomics in neuro-oncology imaging[END_REF] Manual segmentation is time-consuming and is partly subjective.

A large number of methods have been proposed for automated segmentation of brain tumors. [START_REF] Menze | The multimodal brain tumor image segmentation benchmark (BRATS)[END_REF][START_REF] Cabral | Interactive segmentation of brain tumors in MR images using 3D region growing[END_REF][START_REF] Liu | Brain tumor segmentation in MRI by using the fuzzy connectedness method[END_REF][START_REF] Prastawa | A brain tumor segmentation framework based on outlier detection[END_REF][START_REF] Khotanlou | 3D brain tumor segmentation in MRI using fuzzy classification, symmetry analysis and spatially constrained deformable models[END_REF][START_REF] Bauer | A survey of MRI-based medical image analysis for brain tumor studies[END_REF][START_REF] Havaei | Brain tumor segmentation with deep neural networks[END_REF][START_REF] Pereira | Brain tumor segmentation using convolutional neural networks in MRI images[END_REF][START_REF] Rahimpour | Improving T1w MRI-based brain tumor segmentation using cross-modal distillation[END_REF][START_REF] Yerachmiel | Weakly supervised brain tumor segmentation via semantic affinity deep neural network[END_REF] However, while many papers have dealt with other types of tumors such as gliomas, [START_REF] Chow | Updates on deep learning and glioma: use of convolutional neural networks to image glioma heterogeneity[END_REF] only a limited number of studies have focused on brain lymphomas. Some works focused on classification [START_REF] Miao | The value of convolutional neural networks-based deep learning model in differential diagnosis of space-occupying brain diseases[END_REF][START_REF] Hsu | Differentiation of primary central nervous system lymphoma from glioblastoma using optical coherence tomography based on attention ResNet[END_REF] or survival analysis. [START_REF] She | Deep learning-based overall survival prediction model in patients with rare cancer: a case study for primary central nervous system lymphoma[END_REF] Automatic segmentation of PCNSL from MRI data has been performed in a few publications. [START_REF] Pennig | Primary central nervous system lymphoma: Clinical evaluation of automated segmentation on multiparametric MRI using deep learning[END_REF][START_REF] Conte | Generative adversarial networks to synthesize missing t1 and flair mri sequences for use in a multisequence brain tumor segmentation model[END_REF] However, these works were not specific to lymphoma and included also gliomas for training or validation. Moreover, they did not report distance-based metrics for evaluation which is critical when dealing with lymphomas which have complex boundaries and multiple components.

PCNSL has an extremely heterogenous MRI appearance. Lesions can be single or multiple. Several morphologies, topographies, and mass effect patterns can (co)exist. As a result, tumor size, shape and boundary can considerably vary between patients.

These variations pose significant challenges in accurately segmenting lymphomas. Convolutional neural network (CNN)-based segmentation approaches are widely used but often produce errors near boundaries. To address this issue, prior knowledge integration methods at the level of the loss function have been explored to enhance the plausibility of automatic segmentations. [START_REF] El Jurdi | High-level prior-based loss functions for medical image segmentation: A survey[END_REF] In this paper, we propose a novel prior-based loss that integrates border irregularity attributes of the tumor in order to improve segmentation performance. The proposed method was trained and validated for segmentation of PCNSL on post-contrast T1-weighted MR images using a dataset of 99 patients.

The rest of the paper is organized as follows. Section 2 describes the dataset and pre-processing. In Section 3, we introduce the proposed loss. Section 4 is devoted to experiments and results. The discussion is provided in Section 5.

MATERIALS

Participants and MRI data

We studied 118 patients with primary CNS lymphoma. The study was approved by the Institutional Ethical Committee (Pitié Salpêtrière Hospital, Ile-de-France VI, n°DC-2009-957) and by the French Data Protection Authority (CNIL, Commission Nationale de l'Informatique et des Libertés, DR 2013-279). According to French regulation, consent was waived as these images were acquired as part of the routine clinical care of the patients. Each patient had a T1-weighted MRI after injection of gadolinium. The images were acquired as part of clinical routine and were thus not harmonized. They were acquired on different scanners and at different field strengths (57 at 3T, 54 at 1.5T, and 7 at 1T). The MRI scan was either 3D or 2D. Manual segmentations were performed by a trained radiologist (L.N.) who also rated the difficulty of the segmentation process. The segmentation was considered difficult when the lymphoma tissue was less extended, when the lesions boundary were difficult to visualize, and/or, in rarer cases, when there were hemorrhagic remnants that are spontaneously T1 hyperintense and that can therefore mimic tumors. Furthermore, the radiologist noted whether the images presented substantial artifacts. This led to partition the dataset into four subsets: D1, D2, D3, D4 according to these two criteria (D1: easy, no artifact; D2: easy, artifacts; D3: difficult, no artifact; D4: difficult, artifacts).

We applied the following preprocessing. All images were converted from DICOM to NIfTI using dicom2niix 19 * and organized according to the Brain Imaging Data Structure (BIDS) standard. [START_REF] Gorgolewski | The brain imaging data structure, a format for organizing and describing outputs of neuroimaging experiments[END_REF] Using FSL FLIRT, [START_REF] Jenkinson | Improved optimization for the robust and accurate linear registration and motion correction of brain images[END_REF] we linearly registered each image to the MNI-152 template which is 1mm isotropic and of dimensions 181 × 217 × 181. [START_REF] Ciric | Templateflow: Fair-sharing of multi-scale, multi-species brain models[END_REF][START_REF] Mazziotta | A probabilistic atlas of the human brain: Theory and rationale for its development: The international consortium for brain mapping (icbm)[END_REF] We applied a brain mask in order to remove unnecessary information like the skull, nose, and ears. Pydra was used to implement the preprocessing steps. [START_REF] Jarecka | Pydra-a flexible and lightweight dataflow engine for scientific analyses[END_REF] We visually checked the preprocessing results and found that preprocessing failed for 19 patients (9 in D1, 3 in D2, 3 in D3 and 4 in D4). These were excluded from the study. The characteristics of the studied patients and the corresponding subsets are presented in Table 1.

Splits

From D1 and D3, 26 participants (30%) were isolated for testing (21 participants from D1, 5 from D3). The remainder of D1 and D3 were separated into training (80% of the remaining participants) and validation sets (20%). This resulted in 47 participants for training (39 from D1, 8 from D3) and 12 for validation (10 from D1, 2 from D3). D2 and D4 were used only for testing (n = 14

PROPOSED METHOD

Proposed border irregularity loss

In order to include it in a loss function, one needs to be able to quantify the border irregularity of an object. One approach is to smooth the segmentation map, via a Gaussian filter of a kernel size s and standard deviation σ, until it becomes more uniform in shape, and then compare the smoothed segmentation and the original segmentation. [START_REF] Lee | New approach to measure border irregularity for melanocytic lesions[END_REF] A segmentation with a greater degree of irregularity would require stronger smoothing, resulting in a higher index of border irregularity. Conversely, a smoother segmentation would yield a lower index of irregularity. Applying a Gaussian filter with to a segmentation map is a well-established technique in computer vision for achieving smoothness and extracting the global structure of a segmentation map without generating new irregularities, indentations, or protrusions. The irregularity index is defined as

BI = |A * ∪A|-|A * ∩A| |A * | = |A * |+|A|-2|A * ∩A| |A * |
where A * and A are the smoothed and non-smoothed segmentation maps. The border irregularity index represents the level of dissimilarity between the smoothed and non-smoothed segmentations. Another way to express dissimilarity between two segmentations is via the complement of the Dice coefficient. Specifically, the Dice coefficient can be used to derive a border irregularity index as

I b = 1 -2 A * ∩A |A * |+|A| . If in BI, we normalize by |A + A * | instead of |A * |, it is equivalent to I b .
In this paper, we will use I b as a measure of border irregularity at the loss function level, since it leads to a direct implementation of a differentiable and smooth gradient loss function. We optimize the difference between the ground-truth border irregularity index I b and the predicted border irregularity index Îb . The impact of the smoothing procedure on the segmentation is determined by two factors: the smoothing level (σ) and the border irregularity. When dealing with a segmentation map that exhibits severe irregularity, a more intense smoothing is needed, leading to a higher σ value, and conversely. Conversely, a smoother segmentation would result in a lower σ measure of irregularity. In this work, the smoothing maps for the ground-truth segmentations are obtained statically prior to the training process and fed to the framework in order to compute the border irregularity of the prediction maps.

The process can be summarized as follows: 1) irregularities are gradually smoothed out in a systematic manner by applying a variable σ. Smaller irregularities are first eliminated, followed by the larger ones. As some indentations or protrusions are smoothed, they may reveal the presence of a larger irregularity in their respective locations. This larger irregularity is considered global irregularity, while the smaller ones are regarded as "local" irregularities. Consequently, a hierarchical structure of irregularities is established given varying σ values.

Segmentation smoothing module

Implementation of the smoothing module was carried out using a Gaussian kernel of size k = {5, 10} and a standard deviation of σ. A variable standard deviation σ = 2 x for x ∈ {0, 1, 2, 3, ...10} was iteratively applied. The smoothed segmentation map is obtained when the value obtained for I b saturates reaching a fixed value over two iterations. This means that the smoothing process has reached a stable state, and further iterations may not improve the result. The stopping conditions ensure that the smoothing process is performed until a desirable outcome is achieved or until further iterations do not significantly affect the result. They help control the iterations and prevent unnecessary computations, improving the efficiency of the smoothing process.

Model and implementation details

We zero-padded the data to a size of (184, 220, 184). [START_REF] Pérez-García | Torchio: A python library for efficient loading, preprocessing, augmentation and patch-based sampling of medical images in deep learning[END_REF] The MRI scan intensities were normalized between 0 and 1. The border irregularity loss was either used in conjunction with the Dice loss only (corresponding results are denoted as BIL) or in combination with our previously proposed topological loss [START_REF] Fu | Introducing Soft Topology Constraints in Deep Learning-based Segmentation using Projected Pooling Loss[END_REF] and the Dice loss (results denoted as BIL-Topo). The proposed approaches were compared to a baseline which used the Dice loss (corresponding results are denoted as Baseline). Each 3D volume processed as a stack of independent 2D images. The network architecture was a 2D U-net [START_REF] Ronneberger | U-Net: Convolutional networks for biomedical image segmentation[END_REF] which architecture has been used in previous publications. [START_REF] Jurdi | How precise are performance estimates for typical medical image segmentation tasks?[END_REF][START_REF] El Jurdi | A Surprisingly Effective Perimeter-based Loss for Medical Image Segmentation[END_REF] The architecture is a 3-stage structure composed of convolutional, de-convolutional blocks, bottleneck and skip connections. The encoder part is composed of an ensemble of convolutional and batch normalization layers, whereas the decoder part is composed of 2 consecutive convolutional blocks and an upsampling layer in each of the 3 stages. The bottleneck is composed of 2 convolutional blocks separated by a residual block. [START_REF] Zhang | Road extraction by deep residual U-Net[END_REF] The optimizer was Adam and the learning rate was 0.001. The learning rate was halved if the validation performances did not improve over 20 epochs as proposed by. [START_REF] Kervadec | Boundary loss for highly unbalanced segmentation[END_REF] We used batches of 8.

At inference, we predict for each slice independently and then stack the slices belonging to the same patient to form a 3D prediction. Table 2. Results on the whole test set and separately for D1, D2, D3 and D4 test sets. Results are presented as mean and confidence interval (computing using bootstrapping with 15,000 resamplings). D1+D2+D3+D4 refers to the union of the test sets of D1 and D3 together with D2 and D4. HD: 95% 3D Hausdorff distance. MASD: mean average surface distance. n is the number of samples. Best result in each case is in bold. * indicates that the improvement over the baseline is statistically significant.

EXPERIMENTS AND RESULTS

Evaluation framework

We chose the following performance metrics, based on the recommendations of Reinke et al, [START_REF] Reinke | Metrics reloaded -a new recommendation framework for biomedical image analysis validation[END_REF] the 3D Dice coefficient, the 95% 3D Hausdorf distance, and the mean average surface distance (MASD) † . The results are reported as the mean values along with their corresponding confidence intervals, which were computed using bootstrapping with 15,000 resamplings.

Since the dataset is constituted of 3D and 2D MRI scans, we report the mean and the 95 % confidence interval obtained via bootstrapping via 15000 resamples for each of the 2D and 3D patients both combined and separately. † Performance metrics were computed using this code: https://github.com/deepmind/surface-distance

Results

Table 2 presents the results on the whole test set as well as on the different subsets. The proposed method (BIL-Topo) outperformed the baseline across all performance metrics (improvement of about 6 percent points of Dice, 40mm of HD and 6mm of MASD). Confidence intervals are quite wide due to the relatively small size of the test set (n=40) [START_REF] Jurdi | How precise are performance estimates for typical medical image segmentation tasks?[END_REF] but the difference was statistically significant for the boundary-based metrics (HD, MASD).

For the BIL alone, the mean Dice over the entire test set was substantially higher than for the baseline (5 points) and similar to that of the BIL-Topo. On the other hand, the boundary-based metrics were substantially better with the BIL-Topo, demonstrating the added value of the topological loss. Results in the different subsets are consistent with those on the entire set. One can observe that performances tend to be lower in subsets for which the segmentation was considered difficult by the rater (D3 and D4). Some examples of segmentation are shown on Figure 1.

Table 3 presents separately the results obtained on 2D vs 3D acquisitions. For both 2D and 3D acquisitions, the BIL-Topo method achieved higher Dice, lower HD and MASD a compared to the baseline. The improvement was statistically significant for HD and MASD.

Subset

DISCUSSION

The paper introduces a novel border irregularity loss for automatic segmentation of brain lymphomas. It considers border information and combines it with a topological loss to better handle multiple connected components.

Results demonstrate the usefulness of the proposed approaches. The BIL and topological losses capture border characteristics, improve boundary delineation, and enhance segmentation performance. They both resulted in improvement in terms of Dice score over the baseline. The topological loss provided additional improvements in boundary metrics. Nevertheless, the overall performances remain moderate. The average Dice is around 65% which corresponds to a moderate spatial agreement. The relatively high Hausdorff distance mainly reflects the fact that in several cases, some tumor components are missed while erroneous connected components are detected. This highlights that automatic segmentation of brain lymphomas is a very difficult task, in particular when dealing with clinical routine MRI data of heterogeneous quality. Therefore, further work in needed on this application. The present work remains preliminary and has several limitations. Firstly, the BIL exhibits sensitivity towards small connected components. Also, the smoothing parameters are chosen in an ad-hoc manner. Future work should propose more general ways to set these parameters. Finally, it will be necessary to assess the impact of the BIL when associated with other segmentation architectures.

In future work, the aim is to address these limitations. Specifically, efforts will be made to decouple the BIL loss from its reliance on the topological loss, enabling it to effectively handle very small objects and multiconnected components.

Figure 1 .

 1 Figure 1. Examples of segmentations (ground truth, baseline and proposed methods). Red boxes are tumor parts that have been missed by the models.

Table 1 .

 1 ). Characteristics of the study population. Age (in years) is reported as mean±standard-deviation. The age was missing for 8 patients (4 of them were part of D1, 2 of D2, 1 of D3, and 1 of D4). The sex was missing for one patient from D4. The table reports the number of patients for which a 3D acquisition was available (the others had a 2D acquisition). 1T, 1.5T or 3T indicates the MRI magnetic field strength.

	Subset Total	Sex	Age	3D acquisition 1T 1.5T 3T
	D1	70	40 F (57.14%) 61.64 ± 20.28	50	3	36	31
	D2	10	4 F (40.00%) 52.60 ± 29.21	9	-	1	9
	D3	15	7 F (46.66%) 63.20 ± 20.81	11	1	8	6
	D4	4	2 F (66.66%) 47.75 ± 32.37	3	-	2	2
	Total	99	53 F (54.08%) 60.40 ± 21.84	73	4	47	48

Table 3 .

 3 Results on test set depending on whether the acquisition was 2D or 3D. Results are presented as mean and confidence interval (computing using bootstrapping with 15,000 resamplings). HD: 95% 3D Hausdorff distance. MASD: mean average surface distance. n is the number of samples. Best result in each case is in bold. * indicates that the improvement over the baseline is statistically significant.

ACKNOWLEDGMENTS

The research leading to these results has received funding from the French government under the management of Agence Nationale de la Recherche as part of the "Investissements d'avenir" program, reference ANR-19-P3IA-0001 (PRAIRIE 3IA Institute) and reference ANR-10-IAIHU-06 (Agence Nationale de la Recherche-10-IA Institut Hospitalo-Universitaire-6).