Charles-Edouard Br Éhier
email: brehier@math.univ-lyon1.fr

Shmuel Rakotonirina-Ricquebourg

ON ASYMPTOTIC PRESERVING SCHEMES FOR A CLASS OF STOCHASTIC DIFFERENTIAL EQUATIONS IN AVERAGING AND DIFFUSION APPROXIMATION REGIMES

Keywords: Asymptotic preserving schemes, multiscale methods, slow-fast Stochastic Differential Equations, averaging principle, diffusion approximation, weak approximation AMS subject classifications. 65C30, 60H35

We introduce and study a notion of Asymptotic Preserving schemes, related to convergence in distribution, for a class of slow-fast Stochastic Differential Equations. In some examples, crude schemes fail to capture the correct limiting equation resulting from averaging and diffusion approximation procedures. We propose examples of Asymptotic Preserving schemes: when the time-scale separation vanishes, one obtains a limiting scheme, which is shown to be consistent in distribution with the limiting Stochastic Differential Equation. Numerical experiments illustrate the importance of the proposed Asymptotic Preserving schemes for several examples. In addition, in the averaging regime, error estimates are obtained and the proposed scheme is proved to be uniformly accurate.

1. Introduction. Deterministic and stochastic systems are ubiquitous in science and engineering. Traditional modelling and numerical methods become ineffective when systems evolve at different time scales: see for instance the monographs [START_REF]Principles of multiscale modeling[END_REF][START_REF] Kuehn | Multiple time scale dynamics[END_REF] for comprehensive treatment of multiscale dynamics. Averaging and homogenization [START_REF] Pavliotis | Multiscale methods[END_REF] are two popular techniques which are employed to rigorously derive macroscopic limiting equations, starting from (stochastic) slow-fast systems with separated time-scales.

In the last two decades, constructing efficient numerical methods for multiscale stochastic systems has been a very active research area: let us mention the Heterogeneous Multiscale Method (see [START_REF] Abdulle | The heterogeneous multiscale method[END_REF][START_REF] Bréhier | Analysis of an HMM time-discretization scheme for a system of stochastic PDEs[END_REF][START_REF] Liu | Analysis of multiscale methods for stochastic differential equations[END_REF]), projective integration (see [START_REF] Givon | Strong convergence of projective integration schemes for singularly perturbed stochastic differential systems[END_REF]), equation-free coarse-graining (see [START_REF] Kevrekidis | Equation-free, coarse-grained multiscale computation: enabling microscopic simulators to perform system-level analysis[END_REF]), spectral methods (see [START_REF] Abdulle | Spectral methods for multiscale stochastic differential equations[END_REF]), micro-macro acceleration methods (see [START_REF] Vandecasteele | Efficiency of a micro-macro acceleration method for scale-separated stochastic differential equations[END_REF]), parareal algorithms (see [START_REF] Legoll | Parareal computation of stochastic differential equations with time-scale separation: a numerical convergence study[END_REF]). In the methods mentioned above, the objective is to approximate the limiting model for the slow variables of interest, and only partial but relevant information coming from the fast dynamics is taken into account. As a consequence, these methods may not be appropriate if one wants to approximate simultaneously the original multiscale model and its limit. In this article, we focus on the notion of asymptotic preserving schemes, in order to overcome this issue.

To motivate and illustrate our work, let us introduce simplified versions of the systems of Stochastic Differential Equations (SDE) considered in this article. The time-scale separation parameter is denoted by ∈ (0, 1]. On the one hand, in the averaging regime (see Equation (2.1) in Section 2.1 for the more general version), we consider systems of the type (1.1)

     dX t = b(X t , m t )dt, dm t = - m t dt + √ 2 √ dβ t .
When → 0, the averaging principle (see [START_REF] Pavliotis | Multiscale methods[END_REF]Chapter 10]) states that X converges (at least in distribution) to the solution X of the Ordinary Differential Equation Ẋ = b(X) where b(x) = b(x, m)dν(m) and ν = N (0, 1) is the standard Gaussian random variable. On the other hand, in the diffusion approximation regime (see Equation (2.12) in Section 2.2 for the more general version), we consider systems of the type (1.2)

     dX t = σ(X t )m t dt, dm t = - m t 2 dt + 1 dβ t .
When → 0, the diffusion approximation result (see [START_REF] Pavliotis | Multiscale methods[END_REF]Chapter 11]) states that X converges (in distribution) to the solution X of the SDE

dX t = σ(X t ) • dβ t ,
where the noise is interpreted in the Stratonovich sense. This type of results is related to results known as Wong-Zakai approximation and Smoluchowski-Kramers limits in the literature. In the two SDE systems (1.1) and (1.2), the fast component is an Ornstein-Uhlenbeck process.

In this article, we are interested in the behavior when → 0 of numerical schemes for the SDEs (1.1) and (1.2). To explain the challenge faced and the solutions proposed in this article, we consider the following schemes, which are both consistent for any fixed value of > 0. On the one hand, in the averaging regime one defines

(1.3)      X n+1 = X n + ∆tb(X n , m n+1 ), m n+1 = m n - ∆t m n+1 + 2∆t γ n .
On the other hand, in the diffusion approximation regime one defines (1.4)

       X n+1 = X n + σ(X n ) ∆tm n+1 , m n+1 = m n - ∆t 2 m n+1 + √ ∆t γ n ,
In the schemes (1.3) and (1.4), γ n n≥0 is a sequence of independent standard Gaussian random variables. One may check that X n → X n for all n ≥ 0, in probability, when → 0, where the limiting schemes are given by

X n+1 = X n + ∆tb(X n , 0)
in the averaging regime, and

X n+1 = X n + √ ∆tσ(X n )γ n ,
in the diffusion approximation regime. Note that, in the second case, the limiting scheme is consistent with the Itô interpretation of the noise, instead of the correct Stratonovich one. In the two cases, the limiting scheme is in general not consistent with the limiting equation, and using such a scheme in practice may lead to drawing false conclusions about the limiting system from numerical experiments. We refer to [START_REF] Frank | A note on statistical consistency of numerical integrators for multiscale dynamics[END_REF][START_REF] Li | Effectiveness of implicit methods for stiff stochastic differential equations[END_REF] for other examples of situations where numerical schemes perform badly when applied to multiscale SDE systems. The objective of this article is to design and study Asymptotic Preserving (AP) schemes, such that the following diagram commutes (where convergence is understood in distribution): if T = N ∆t, one has

X N ∆t→0 ----→ X (T )   →0   →0 X N ∆t→0 ----→ X(T )
The two schemes (1.3) and (1.4) described above are not AP. The notion of AP schemes has been introduced in [START_REF] Jin | Efficient asymptotic-preserving (AP) schemes for some multiscale kinetic equations[END_REF], for applications to multiscale kinetic Partial Differential Equations (PDEs), which converge to parabolic diffusion PDEs. We refer to [START_REF] Dimarco | Numerical methods for kinetic equations[END_REF]Section 7], [START_REF] Hu | Asymptotic-preserving schemes for multiscale hyperbolic and kinetic equations[END_REF], [START_REF] Jin | Asymptotic preserving (AP) schemes for multiscale kinetic and hyperbolic equations: a review[END_REF] and [START_REF] Puppo | Kinetic models of bgk type and their numerical integration[END_REF]Section 4] for recent reviews on AP schemes for this type of models. To the best of our knowledge, the design and analysis of Asymptotic Preserving schemes for slow-fast SDEs of the type (1.1) and (1.2) has not been considered so far in the literature. Note that a specific feature (compared with the deterministic case) is the need to consider convergence in distribution. Let us mention related works for Stochastic Partial Differential Equations (SPDEs), in the diffusion approximation regime. First, in [START_REF] Duboscq | Analysis of a splitting scheme for a class of random nonlinear partial differential equations[END_REF][START_REF] Marty | On a splitting scheme for the nonlinear Schrödinger equation in a random medium[END_REF], the authors consider Schrödinger equations and study an abstract asymptotic preserving property. However, they do not propose implementable schemes. In [START_REF] Ayi | Analysis of an asymptotic preserving scheme for stochastic linear kinetic equations in the diffusion limit[END_REF], the authors deal with some multiscale stochastic kinetic PDEs, driven by a Wiener process. However, the structure of the model is different from the one of (1.2). In a future work [START_REF] Bréhier | Asymptotic-preserving schemes for stochastic kinetic equations in the diffusion approximation regime[END_REF], we plan to apply the findings of this article to the SPDE models considered in [START_REF] Ayi | Analysis of an asymptotic preserving scheme for stochastic linear kinetic equations in the diffusion limit[END_REF]. The works mentioned above concerning SPDE models are limited to diffusion coefficients of the type σ(x) = x, for which specific arguments may give a straightforward construction of AP schemes, for appropriate discretization of the fast component. An AP scheme in the case σ(x) = 1 for (1.2) is proposed in [START_REF] Pavliotis | Calculating effective diffusiveness in the limit of vanishing molecular diffusion[END_REF], however the subtlety of the interpretation of the noise at the limit is not relevant in that case. Finally, let us also mention that AP schemes have also been studied for PDEs with random coefficients, see [START_REF] Hu | Uncertainty quantification for kinetic equations[END_REF][START_REF] Jin | Mathematical analysis and numerical methods for multiscale kinetic equations with uncertainties[END_REF][START_REF] Jin | Efficient stochastic asymptotic-preserving implicit-explicit methods for transport equations with diffusive scalings and random inputs[END_REF] or in the context of Monte-Carlo methods for deterministic problems, see [START_REF] Dimarco | Asymptotic-preserving Monte Carlo methods for transport equations in the diffusive limit[END_REF][START_REF] Ren | An asymptotic-preserving Monte Carlo method for the Boltzmann equation[END_REF].

We are now in position to describe the contributions of this article. In Section 3.1, we define the appropriate notion of AP schemes for SDE systems, related to convergence in distribution, and study several general properties.

Our first main result is Theorem 3.7, which exhibits an example of AP scheme in the averaging regime: for the simplified version (1.1), the scheme is given by

(1.5)    X n+1 = X n + ∆tb(X n , m n+1 ), m n+1 = e -∆t m n + 1 -e -2∆t γ n .
The fast component in the scheme above is discretized using a scheme which is exact in distribution.

Our second main result is Theorem 3.8, which states error estimates of the type sup

∈(0,1] E[ϕ(X N ) -E[ϕ(X (T ))] = O √ ∆t ,
for sufficiently smooth real-valued mappings ϕ. This error estimate means that the scheme is Uniformly Accurate.

Finally, our third main result is Theorem 3.9, which exhibits an example of AP scheme in the diffusion approximation regime: for the simplified version (1.2) (see Corollary 3.11), the scheme is given by

(1.6)                m n+1 = m n - ∆t 2 m n+1 + √ ∆t γ n , Y n+1 = X n + σ(X n ) ∆tm n+1 , X n+1 = X n + σ(X n ) + σ(Y n+1 ) 2 ∆tm n+1 .
A prediction-correction method is employed to retrieve the correct interpretation of the noise for the limiting equation: the scheme (1.2) is indeed consistent with the Stratonovich interpretation of the noise.

Let us also mention that another situation is considered in Corollary 3.12: for the model (2.22) taken from [START_REF] Laibe | On the settling of small grains in dusty discs: analysis and formulae[END_REF] (with an application in astrophysics), the limiting equation (2.23) contains a so-called noise-induced drift-term, which is captured only for well-designed AP schemes.

Some numerical experiments (see Section 4) show that the AP schemes (1.5) and (1.6) are effective in all regimes > 0 and → 0, contrary to the schemes (1.3) and (1.4) which fail to capture the correct limiting behavior when → 0.

The article is organized as follows. The general SDE models in the averaging and diffusion approximation regimes are presented in Sections 2.1 and 2.2. The main results of this article are stated in Section 3: the general theory of AP schemes is presented in Section 3.1, and it is applied in the averaging and diffusion approximation regimes in Section 3.2 and 3.3 respectively. Numerical experiments are reported in Section 4. Section 5 is devoted to the proof of the error estimates stated in Theorem 3.8. Finally, Section 6 gives some conclusions and perspectives.

2. Slow-fast SDE models and their limits. Without loss of generality, the time-scale separation parameter satisfies ∈ (0, 1]. The time-step size of the integrators studied in this work is denoted by ∆t. It is assumed that ∆t = T N where T ∈ (0, ∞) is a fixed time and N ∈ N. Without loss of generality, it is assumed that ∆t ∈ (0, 1].

In the slow-fast systems considered in this work, the slow component X takes values in the d-dimensional flat torus T d , where d ∈ N is an arbitrary integer, whereas the fast component m takes values in R. The framework and the models considered in this work may be generalized in many ways to more complex situations, however the arguments and results below are sufficient to illustrate the difficulties of designing asymptotic preserving schemes for stochastic equations.

Let β t t≥0 and B t t≥0 be two independent standard Wiener processes, with values in R and R D respectively, where D ∈ N, defined on a probability space (Ω, F, P) which satisfies the usual conditions.

The following notation for derivatives is used below: ∇ x = ∂ xi 1≤i≤d ∈ R d and ∂ m are the partial gradient and derivative operators with respect to x and m respectively. If σ is a mapping with values in M d,D (R) (the space of d×D matrices with real entries), let σ denote the transpose of σ, and set σσ :

∇ 2 x = d i,j=1 σσ ) i,j ∂ xi ∂ xj . If b is a R d -valued mapping, let b • ∇ x = d i=1 b i ∂ xi . Assumption 2.1.
The initial conditions X 0 ∈ T d and m 0 ∈ R of the processes are deterministic quantities and they satisfy

X 0 = x 0 ---→ →0 x 0 , sup ∈(0,1] |m 0 | < ∞.
2.1. The averaging regime. In the so-called averaging regime, we consider slow-fast SDE systems of the type (2.1)

     dX t = b(X t , m t )dt + σ(X t , m t )dB t , dm t = - m t dt + √ 2h(X t ) √ dβ t .
The coefficients appearing in (2.1) are assumed to satisfy the following conditions.

Assumption 2.2. The functions b :

T d × R → R d and σ : T d × R → M d,D (
R) are assumed to be of class C 4 , and h : T d → R is assumed to be of class C 1 . Moreover, they are all assumed to be bounded and to have bounded derivatives.

Owing to Assumption 2.2, for all initial conditions X 0 ∈ T d and m 0 ∈ R, and for every ∈ (0, 1], there exists a unique global solution X (t), m (t) t≥0 of the SDE system (2.1). Since h is bounded, it is straightforward to check that (2.2) sup

∈(0,1] sup t≥0 E[|m (t)| 2 ] 1 + |m 0 | 2 < ∞.
This estimate will prove useful to prove Proposition 2.4.

The infinitesimal generator L associated with the SDE (2.1) has the following expression:

(2.3) L = 1 L OU + L 0 , where (2.4) L 0 = b(x) • ∇ x + 1 2 σσ : ∇ 2 x , L OU = -m∂ m + h(x) 2 ∂ 2 m ,
Observe that for fixed x ∈ T d , L OU is the generator of an ergodic Ornstein-Uhlenbeck process. The associated invariant distribution is

ν x = N (0, h(x) 2 ).
Define averaged coefficients as follows: for all x ∈ T d

(2.5) b(x) = b(x, m)dν x (m) , a(x) = σ(x, m)σ(x, m) dν x (m).
Note that b :

T d → R d is of class C 4 .
The averaging principle result stated below requires the following condition to be satisfied.

Assumption 2.3.

There exists an integer D ∈ N and a function σ :

T d → M d,D (R) of class C 4 such that for all x ∈ T d (2.6) a(x) = σ(x)σ(x) .
Assumption 2.3 holds if there exists c ∈ (0, ∞) such that a(x) ≥ cI for all x ∈ T d (as symmetric matrices). This condition is satisfied when σ only depends on the slow variable x (∂ m σ(x, m) = 0 for all (x, m) ∈ T d × R), or when σ(x, m)σ(x, m) ≥ cI for all (x, m) ∈ T d × R. In that case, one can choose

D = d. If the diffusion coefficient is of the type σ(x, m) = σ (m)σ † (x), with σ † (x) ∈ R d and σ (m) ∈ R, then one can choose D = D and σ(x) = σ † (x)
σ (m) 2 dν x (m) for all x ∈ T d . We are now in position to state the averaging principle result and to define the limiting process X obtained when → 0.

Proposition 2.4. Let Assumptions 2.1, 2.2 and 2.3 be satisfied. Let T ∈ (0, ∞). When → 0, the C([0, T ], T d )-valued process X (t) 0≤t≤T converges in distribution to the solution X(t) 0≤t≤T of the limiting SDE

(2.7) dX t = b(X t )dt + σ(X t )dB t ,
with initial condition X(0) = x 0 , where the coefficients b and σ are defined by (2.5)-(2.6), and where B t t≥0 is a standard R D -valued Wiener process.

The infinitesimal generator L associated with the limiting SDE (2.7) is given by

(2.8) L = b(x) • ∇ x + 1 2 σ σ : ∇ 2 x ,
and is such that the following property holds: let ϕ ∈ C 4 (T d ), then there exists a function ϕ 1 :

T d × R → R such that ϕ = ϕ + ϕ 1 , (2.9) L ϕ ---→ →0 Lϕ. (2.10)
Finally, let ϕ ∈ C 4 (T d ), then there exists C(T, ϕ) ∈ (0, ∞) such that

(2.11) E[ϕ(X (T ))] -E[ϕ(X(T ))] ≤ C(T, ϕ) .
The averaging principle stated in Proposition 2.4 is a standard result, see for instance [START_REF] Pavliotis | Multiscale methods[END_REF]Chapter 16]. In general the convergence stated in Proposition 2.4 only holds in distribution, however it holds in stronger sense (for instance in mean-square sense) if σ only depends on x. We refer to Appendix A.1 for a sketch of the construction of the perturbed test function ϕ which satisfies (2.9)-(2.10) (see [START_REF] Fouque | Wave propagation and time reversal in randomly layered media[END_REF]Chapter 6] for a detailed description of the perturbed test function method). Note that the perturbed test function appears in Proposition 3.5 below. For the error estimate (2.11), see Lemma 5.2 and its proof below.

The diffusion approximation regime.

2.2.1. General model. In the so-called diffusion approximation regime, we consider slow-fast SDE systems of the type (2.12)

     dX t = b(X t )dt + σ(X t )m t dt, dm t = f (X t ) - m t 2 dt + g(X t ) dt + h(X t ) dβ t .
The coefficients appearing in (2.12) are assumed to satisfy the following conditions.

Assumption 2.5. The functions b : T d → R d and g, h : T d → R are assumed to be of class C 1 . The functions σ : T d → R d and f : T d → R are assumed to be of class C 2 . Moreover, f takes values in (0, ∞): we assume that min

x∈T d f (x) > 0.
Owing to Assumption 2.5, for all initial conditions X 0 ∈ T d and m 0 ∈ R, and for every ∈ (0, 1], there exists a unique global solution X (t), m (t) t≥0 of the SDE system (2.12). The infinitesimal generator L associated with the SDE (2.12) has the following expression:

(2.13) L = 1 2 L OU + 1 L 1 + L 0 , where (2.14)
L 0 = b(x) • ∇ x , L 1 = mσ(x) • ∇ x + f (x)g(x)∂ m , L OU = -f (x)m∂ m + 1 2 f (x) 2 h(x) 2 ∂ 2 m .
Observe that for fixed x ∈ T d , L OU is the generator of an ergodic Ornstein-Uhlenbeck process. The associated invariant distribution is

ν x = N (0, f (x)h(x) 2 2
). We are now in position to state the diffusion approximation result and to define the limiting process X obtained when → 0. Proposition 2.6. Let Assumptions 2.1 and 2.5 be satisfied. Let T ∈ (0, ∞). When → 0, the C([0, T ], T d )-valued process X (t) 0≤t≤T converges in distribution to the solution X(t) 0≤t≤T of the limiting SDE (2.15)

dX t = b + gσ + h 2 2 (σ • ∇ x )σ - h 2 2f σ • ∇ x f σ (X t )dt + h(X t )σ(X t )dW t ,
driven by a standard one-dimensional Wiener process W (t) t≥0 , with initial condition X(0) = x 0 . The infinitesimal generator L associated with the limiting SDE (2.15) is given by

Lϕ = (b + gσ) • ∇ x ϕ + h 2 f σ 2 • ∇ x σ f • ∇ x ϕ (2.16) = (b + gσ) • ∇ x ϕ + h 2 2 σσ : ∇ 2 x ϕ + h 2 2 (σ • ∇ x )σ • ∇ x ϕ - h 2 2f σ • ∇ x f σ • ∇ x ϕ,
and is such that the following property holds: let ϕ ∈ C 3 (T d ), then one constructs two functions ϕ 1 , ϕ 2 :

T d × R → R, such that ϕ = ϕ + ϕ 1 + 2 ϕ 2 , (2.17) L ϕ ---→ →0 Lϕ. (2.18)
Finally, let ϕ ∈ C 3 (T d ), then there exists C(T, ϕ) ∈ (0, ∞) such that

(2.19) E[ϕ(X (T ))] -E[ϕ(X(T ))] ≤ C(T, ϕ) .
The diffusion approximation stated in Proposition 2.6 is a standard result, see for instance [START_REF] Pavliotis | Multiscale methods[END_REF]Chapter 18]. We refer to Appendix A.2 for a sketch of the construction of the perturbed test function ϕ which satisfies (2.17)-(2.18) (see [START_REF] Fouque | Wave propagation and time reversal in randomly layered media[END_REF]Chapter 6] for a detailed description of the perturbed test function method). Since the error estimate (2.19) plays no role in the sequel, the proof is omitted. We refer to [START_REF] Khasminskii | Limit behavior of two-time-scale diffusions revisited[END_REF] for arguments using asymptotic expansions of solutions of Kolmogorov equations leading to (2. [START_REF] Jin | Asymptotic preserving (AP) schemes for multiscale kinetic and hyperbolic equations: a review[END_REF]), (see also [START_REF] Laibe | On the settling of small grains in dusty discs: analysis and formulae[END_REF] for related computations).

2.2.2.

Two examples in the approximation-diffusion regime. The setting described above encompasses several interesting examples of SDE systems. In order to focus on the different possible issues which need to be overcome when constructing asymptotic preserving numerical schemes in the regime → 0, we deal with two examples described below. In addition, the asymptotic preserving numerical schemes will have simpler formulations for these examples than in the general case. In both examples, dimension is set equal to d = 1 to simplify the presentation, and b = 0.

Let us present the first example: consider the system (2.20)

     dX t = σ(X t )m t dt, dm t = - m t 2 dt + 1 dβ t ,
where the coefficients in the fast equation are constant: f (x) = h(x) = 1 and g(x) = 0 for all x ∈ T. Applying Proposition 2.6 in this example yields the following limiting equation

(2.21) dX t = σ(X t ) • dW t ,
where the noise is interpreted using the Stratonovich convention. With the Itô convention, the equation is written as

dX t = 1 2 σ(X t )σ (X t )dt + σ(X t )dW t .
Note that the diffusion approximation result (Proposition 2.6) may be obtained by straightforward arguments in two cases, which will be repeated at the discrete-time levels. Let ζ (t) = 1 t 0 m (s)ds for all t ≥ 0. First, if σ(x) = 1 for all x ∈ T, then one has dX t = dζ t . Therefore passing to the limit yields

X (t) = X 0 + ζ (t) ---→ →0 x 0 + W (t),
and the limiting equation is dX t = dW t . Second, assume that x, X (t) and X(t) take values in the real line R (instead of the torus T) and that σ(x) = x for all x ∈ R. Then (2.21) is written as dX t = X t dζ t . Computing the solution and passing to the limit then yields

X (t) = X 0 exp ζ (t) ---→ →0 x 0 exp(W (t)) = X(t),
and the limiting equation is dX t = X t • dW t .

Note that when the function σ is not constant, the Itô and Stratonovich interpretations differ. Constructing an asymptotic preserving requires to capture the correction term in a limiting scheme (which will naturally be associated with an Itô interpretation of the noise).

Let us now present the second example, taken from [START_REF] Laibe | On the settling of small grains in dusty discs: analysis and formulae[END_REF]. The coefficients f, g, h are allowed to depend on the slow component x, whereas it is assumed that σ(x) = 1 for all x ∈ T. Therefore, the system in the second example has the following expression

(2.22)      dX t = m t dt, dm t = f (X t ) - m t 2 dt + g(X t ) dt + h(X t ) dβ t ,
Applying Proposition 2.6 in this example yields the following limiting equation

(2.23) dX t = g(X t )dt - h(X t ) 2 f (X t ) 2f (X t ) dt + h(X t )dW t .
The noise is interpreted in the Itô sense. Observe that when f is not constant, the noise-induced drift term h 2 f 2f appears. The construction of asymptotic preserving schemes for this problem requires to be careful in order to capture this additional drift term in the limiting scheme.

3. Numerical discretization and asymptotic preserving schemes. The objective of this section is to study the notion of Asymptotic Preserving (AP) schemes for the slow-fast SDE system (2.1) (averaging regime) or (2.12) (diffusion approximation regime) when → 0. The fundamental requirements to have an AP scheme are the following ones: given a consistent discretization scheme for the SDE system,

• for any fixed time-step size ∆t > 0, there exists a limiting scheme when → 0,

• this limiting scheme is consistent with the limiting equation (2.7) given by Proposition 2.4 (averaging regime), or the limiting equation (2.15) given by Proposition 2.6 (diffusion approximation regime). For the SDE considered in this article, consistency is understood in the sense of convergence in distribution. As will be clear below, caution is needed in order to satisfy the second requirement, indeed some standard but naive schemes converge to a limiting scheme which is not consistent with the correct limiting equation. Using such schemes would be dangerous since it could lead to wrong conclusions about the behavior of the SDE system when → 0, hence the need to develop simultaneously the theoretical and numerical analysis.

After discussing general properties of AP schemes, we will provide example of such schemes both for the system (2.1) (averaging regime) and for the system (2.12) (diffusion approximation regime) . We will also study how this scheme applies to the two examples (2.20) and (2.22) described above, and provide a few examples of non AP schemes.

Asymptotic Preserving schemes: definition and properties. Let

T ∈ (0, ∞), and let N ∈ N and ∆t = T N denote the time-step size. Let Γ n 0≤n≤N -1 and γ n 0≤n≤N -1 be two independent families of independent standard R D and Rvalued Gaussian random variables. The initial conditions X 0 and m 0 are assumed to satisfy Assumption 2.1.

On the one hand, a discretization scheme for the SDE (2.1) is defined as

(3.1) (X n+1 , m n+1 ) = Φ ∆t (X n , m n , Γ n , γ n ), n = 0, . . . , N -1.
On the other hand, a discretization scheme for the SDE (2.12) is defined as

(3.2) (X n+1 , m n+1 ) = Φ ∆t (X n , m n , γ n ), n = 0, . . . , N -1.
The presentation is slightly different in the averaging and diffusion approximation regimes. In the remaining of Section 3.1, only the case of schemes of the type (3.1) is considered. This means that if one considers the SDE (2.12) and the scheme (3.2) (approximation diffusion regime) the variable Γ n needs to be omitted -this is also the case if σ = 0 in the SDE (2.1) (averaging regime).

The mapping Φ ∆t appearing in the schemes (3.1) and (3.2) is referred to as the integrator in the sequel.

Let us first discuss stability issues. Due to the presence of factors 1 and 1 2 in the SDE (2.1) and (2.12), using the standard Euler-Maruyama scheme would impose strong stability conditions, of the type ∆t ≤ ∆t 0 ( ) with ∆t 0 ( ) → 0 when → 0. In order to study the behavior of the scheme when → 0 for any fixed time-step size ∆t, it is necessary to avoid such conditions, and we impose the following assumption (which is generally satisfied for some implicit or implicit-explicit methods).

Assumption 3.1. The integrator Φ ∆t is defined for all ∈ (0, 1] and ∆t ∈ (0, ∆t 0 ], where ∆t 0 > 0 is independent of .

We are now in position to study the consistency of the scheme. First, it is assumed that for all ∈ (0, 1], the scheme (3.1) (resp. (3.2)) is consistent with the SDE system (2.1) (resp. (2.12)). When dealing with numerical methods for SDEs, there exist several notions of convergence: in almost sure sense, in probability, in mean-square sense, or in distribution. Since Propositions 2.4 and 2.6 state that X converges in distribution to X when , the relevant notion is consistency in the weak sense, related to convergence in distribution. Assumption 3.2. For all ∈ (0, 1], the numerical scheme (3.1) (resp. (3.2)) is consistent in the weak sense with the SDE system (2.1) (resp. (2.12)): for all bounded continuous functions ϕ :

T d × R, E[ϕ(X N , m N )] ----→ N →∞ E[ϕ(X (T ), m (T ))],
where the time-step size is given by ∆t = T N , for an arbitrary T ∈ (0, ∞). Recall that the consistency in the weak sense of the scheme can be verified using the following equivalent criterion, expressed in terms of the integrator and of the infinitesimal generator: for all

ϕ ∈ C 2 b (T d × R), lim ∆t→0 E[ϕ(Φ ∆t (x, m, Γ, γ))] -ϕ(x, m) ∆t = L ϕ(x, m),
for all (x, m) ∈ T d ×R, where Γ and γ are two independent standard R D and R-valued Gaussian random variables.

The requirements above (Assumptions 3.1 and 3.2) only depend on the behavior of the scheme for fixed ∈ (0, 1]. We are now in position to study the asymptotic behavior as → 0, with fixed time-step size ∆t ∈ (0, ∆t 0 ]. To introduce the notion of Asymptotic Preserving scheme, one first needs to assume the existence of a limiting scheme, as follows.

Assumption 3.3. For every ∆t ∈ (0, ∆t 0 ], there exists a mapping Φ ∆t : T d × R 2 → T d , such that for every (x, m) ∈ T d × R, and every bounded continuous function ϕ :

T d → R, E[ϕ(Φ ∆t (x, m, Γ, γ))] ---→ →0 E[ϕ(Φ ∆t (x, Γ, γ))]
where Γ and γ are two independent standard R D and R valued Gaussian random variables.

Let X n 0≤n≤N be defined by

(3.3) X n+1 = Φ ∆t (X n , Γ n , γ n ), X 0 = x 0 = lim →0 x 0 .
where Γ n 0≤n≤N -1 and γ n 0≤n≤N -1 are two independent families of independent standard R D and R valued Gaussian random variables. By a recursion argument, it is straightforward to check that if Assumptions 2.1 and 3.3 are satisfied, then X n converges in distribution to X n , when → 0, for any fixed ∆t ∈ (0, ∆t 0 ], and 0 ≤ n ≤ N .

We are now in position to introduce the notion of Asymptotic Preserving schemes. As for Assumptions 3.2 and 3.3 above, the consistency is understood in the sense of convergence in distribution. Definition 3.4. Let Assumptions 3.1, 3.2 and 3.3 be satisfied. The scheme (3.1) (resp. (3.2)) is said to be Asymptotic Preserving (AP) if the limiting scheme given by Assumption 3.3 and (3.3) is consistent, in the weak sense, with the limiting equation given by Proposition 2.4 (resp. Proposition 2.6): for every continuous function ϕ :

T d → R, one has E[ϕ(X N )] ----→ N →∞ E[ϕ(X(T ))],
where ∆t = T N , with an arbitrary T ∈ (0, ∞). One of the main contributions of this article is the design of AP schemes in the averaging and in the diffusion approximation regimes, see Sections 3.2 and 3.3 respectively.

To conclude this section, Proposition 3.5 and Corollary 3.6 below are general formulations of the AP property in terms of interverting the limits → 0 and ∆t → 0. As explained above, the result is stated only in the averaging regime to simplify the presentation, however the same result holds also in the diffusion approximation regime with straightforward modifications. Proposition 3.5. Let the setting of Definition 3.4 be satistied. The following statements are equivalent.

(i) The scheme (3.1) is Asymptotic Preserving.

(ii) For any continuous function ϕ :

T d → R, one has lim ∆t→0 lim →0 E[ϕ(X N )] = lim →0 lim ∆t→0 E[ϕ(X N )],
where

T = N ∆t. (iii) For any ϕ ∈ C 3 (T d ), for all (x, m) ∈ T d × R, one has lim ∆t→0 lim →0 E[ϕ (Φ ∆t (x, m, Γ, γ))] -ϕ(x) ∆t = lim →0 lim ∆t→0 E[ϕ (Φ ∆t (x, m, Γ, γ))] -ϕ(x) ∆t ,
where ϕ = ϕ + ϕ 1 is the function introduced by the perturbed test function approach (see (2.9), Proposition 2.4 or (2.17), Proposition 2.6), and Γ and γ are independent R D and R valued standard Gaussian random variables.

Note that using the perturbed test function approach (see Propositions 2.4 and 2.6) is the relevant point of view for the statement (iii) above.

Proof of Proposition 3.5. The equivalence of (i) and (ii) is straightforward. Indeed

lim ∆t→0 lim →0 E[ϕ(X N )] = lim ∆t→0 E[ϕ(X N )], lim →0 lim ∆t→0 E[ϕ(X N )] = lim →0 E[ϕ(X (T ))] = E[ϕ(X(T ))],
using Assumptions 3.2 and 3.3 and Proposition 2.4. The two quantities coincide if and only if the limiting scheme is consistent with the limiting equation.

It remains to prove that (i) and (iii) are equivalent. On the one hand, note that lim

∆t→0 lim →0 E[ϕ (Φ ∆t (x, m, Γ, γ))] -ϕ(x) ∆t = lim ∆t→0 E[ϕ(Φ ∆t (x, Γ, γ))] -ϕ(x) ∆t ,
using the fact that ϕ -ϕ = O( ) and the definition of the limiting scheme from Assumption 3.3. On the other hand, one has lim

→0 lim ∆t→0 E[ϕ (Φ ∆t (x, m, Γ, γ))] -ϕ(x) ∆t = lim →0 L ϕ (x, m) = Lϕ(x),
using the consistency of the scheme for fixed (Assumption 3.2), and the property (2.10), by construction of the perturbed test function ϕ . Then (iii) is equivalent to having

lim ∆t→0 E[ϕ(Φ ∆t (x, γ))] -ϕ(x) ∆t = Lϕ(x),
which means consistency in the weak sense of the limiting scheme with the limiting equation (2.7). This concludes the proof of Proposition 3.5.

The following result is a simple criterion to check whether a scheme satisfies the asymptotic preserving property.

Corollary 3.6. Assume that for all ϕ ∈ C 2 (T d ), one has

Lϕ(x) = lim ∆t→0 E[ϕ(Φ ∆t (x, γ))] -ϕ(x) ∆t
where L is a second-order differential operator. Then the scheme is AP if and only if the property stated in (iii) in Proposition 3.5 holds with ϕ(x) = x i and ϕ(x) = x i x j , with 1 ≤ i, j ≤ d.

The proof of Corollary 3.6 is straightforward and is thus omitted.

3.2. An example of AP scheme in the averaging regime. The objective of this section is to propose an example of AP for the SDE model (2.1), see Theorem 3.7, in the averaging regime. The challenge is to capture the averaged coefficients b and σ, given by (2.5) and (2.6).

Theorem 3.7. Introduce the following numerical scheme:

(3.4)    X n+1 = X n + ∆tb(X n , m n+1 ) + √ ∆tσ(X n , m n+1 )Γ n m n+1 = e -∆t m n + 1 -e -2∆t h(X n )γ n .
This scheme satisfies Assumptions 3.1, 3.2 and 3.3 and is Asymptotic Preserving in the sense of Definition 3.4. Moreover the limiting scheme is given by

(3.5) X n+1 = X n + ∆tb(X n , h(X n )γ n ) + √ ∆tσ(X n , h(X n )γ n )Γ n ,
Let us discuss some properties of the AP scheme (3.4) and of the limiting scheme (3.5).

To simplify the discussion, assume that h(x) = 1. First, assume that σ = 0. Note that even if the limiting equation (2.7) is a deterministic ordinary differential equation, the scheme (3.5) is random. However, in that case, the convergence of X N to X(T ) when ∆t → 0 holds in probability, instead of only in distribution; in that case, the averaging principle result stated in Proposition 2.4 also holds in probability (and even in mean-square sense). The fundamental property to obtain the AP property is that the random quantity appearing in the limiting scheme (3.5) satisfies the property

(3.6) E[b(X n , h(X n )γ n )|X n ] = b(X n ).
In the AP scheme (3.4), the fast component is discretized exactly in distribution (when h(x) = 1): for all n ≥ 0, the Gaussian random variables m n and m (n∆t) are equal in distribution. The fundamental property written above cannot be satisfied if one uses for instance the implicit Euler scheme to discretize the fast component: the scheme defined by

(3.7)      X n+1 = X n + ∆tb(X n , m n+1 ) m n+1 = m n - ∆t m n+1 + 2 ∆t γ n ,
is not asymptotic preserving, since the associated limiting scheme is

X n+1 = X n + ∆tb(X n , 0).
using the identity

m n+1 = 1 1 + ∆t m n + 2 ∆t 1 + ∆t γ n → →0 0,
to pass to the limit.

Second, assume that σ is not equal to 0. Then the convergence of X n to X n only holds in distribution in general. It does not hold in mean-square sense in the following case: assume that d = 1, and that b(x) = 0 and σ(x, m) = m (in that example, the convergence in Proposition 2.4 also does not hold in the mean-square sense). Assume also for simplicity that x 0 = 0, and that m 0 = m 0 ∼ N (0, 1) (and is independent of the Wiener processes β and B). Then one has σ(x) = 1, thus

X n = n-1 k=0 √ ∆tΓ k , X n = n-1 k=0 √ ∆tm k+1 Γ k ,
and one obtains

E|X n -X n | 2 = ∆t n-1 k=0 E|m k+1 -1| 2 = n∆tE|m 0 -1| 2 ,
and the right-hand side does not depend on . It is thus natural to consider convergence in distribution in the notion of asymptotic preserving schemes for SDEs.

Finally, note also that, as above, the scheme

     X n+1 = X n + √ ∆tσ(X n , m n+1 )Γ n m n+1 = m n - ∆t m n+1 + 2 ∆t γ n ,
is not asymptotic preserving, since the associated limiting scheme is

X n+1 = X n + √ ∆tσ(X n , 0)Γ n .
We are now in position to prove Theorem 3.7.

Proof of Theorem 3.7. It is straightforward to check that Assumption 3.1 is satisfied. Let us prove that Assumption 3.3 holds. We have

Φ ∆t (x, m, Γ, γ) = x + ∆tb(x, m ) + √ ∆tσ(x, m )Γ Φ ∆t (x, Γ, γ) = x + ∆tb(x, h(x)γ) + √ ∆tσ(x, h(x)γ)Γ, with m = e -∆t m + 1 -e -2∆t h(x)γ. When → 0,
sup n≥0 E[ m n 2 ] 1 + m 0 2 < ∞.
This estimate will prove useful to prove Lemma 5.4 of Theorem 3.8. Since b and σ are also bounded, we get, in L 1 (Ω), when ∆t → 0

m n+1 -m n = - ∆t m n + 2∆t h(X n )γ n + o(∆t), X n+1 -X n = ∆tb(X n ) + √ ∆tσ(X n )Γ n + o(∆t).
Thus, using that Γ n , γ n and X n are independent, we get the second order Taylor

expansion of ϕ ∈ C 2 b (T d × R), E[ϕ(X n+1 , m n+1 )] -E[ϕ(X n , m n )] = ∆tE[b(X n , m n+1 ) • ∇ x ϕ(X n , m n )] + 1 2 ∆tE[σσ * (X n , m n+1 ) : ∇ 2 x ϕ(X n , m n )] - ∆t E[m n ∂ m ϕ(X n , m n )] + ∆t E[h(X n ) 2 ∂ 2 m ϕ(X n , m n )] + o(∆t) = ∆tE[L ϕ(X n , m n )] + o(∆t).
From there, it is straightforward to check that Assumption 3.2 is satisfied.

Similarly, to prove the consistency of the limiting scheme (3.5) with (2.7), for ϕ ∈ C 2 (T d ), when ∆t → 0, observe that one has

E[ϕ(X n+1 )] -E[ϕ(X n )] = ∆tE[b(X n , h(X n )γ n ) • ∇ x ϕ(X n )] + 1 2 ∆tE[σσ * (X n , h(X n )γ n ) : ∇ 2 x ϕ(X n )] + o(∆t).
The key argument of this proof is the following: by conditioning with respect to X n and the definitions (2.5)-(2.6) of the averaged coefficients, using the fundamental property (3.6) for b and σ σ * = σσ * , yields

E[ϕ(X n+1 )] = E[ϕ(X n )] + ∆tE[Lϕ(X n )] + o(∆t).
The limiting scheme is thus consistent with the limiting equation. This concludes the proof of Theorem 3.7.

Beyond the asymptotic preserving property, it is possible to obtain error estimate, and to prove that the scheme (3.4) given in Theorem 3.7 is uniformly accurate (in distribution).

Theorem 3.8. Let Assumptions 2.1, 2.2 and 2.3 be satisfied. For any T ∈ (0, ∞) and any function ϕ : T d → R of class C 4 , there exists C(T, ϕ) ∈ (0, ∞) such that for all ∆t ∈ (0, ∆t 0 ] and ∈ (0, 1] one has

(3.9) E[ϕ(X N ) -E[ϕ(X (T ))] ≤ C(T, ϕ) min ∆t , ∆t + ,
and the scheme (3.4) is uniformly accurate with the following error estimate: for all ∆t ∈ (0, ∆t 0 ], one has

(3.10) sup ∈(0,1] E[ϕ(X N ) -E[ϕ(X (T ))] ≤ C(T, ϕ) √ ∆t.
The error estimate (3.10) implies that the error E[ϕ(X N ) -E[ϕ(X (T ))] goes to 0 when ∆t → 0 uniformly with respect to → 0. Note that (3.10) is a straightforward consequence of (3.9), considering the cases √ ∆t ≤ and ≤ √ ∆t separately. This argument implies a reduction in the order of convergence appearing in (3.10): it is equal to 1 2 whereas for fixed > 0 (in (3.9)) or when = 0 the order of convergence is equal to 1.

The proof of Theorem 3.8 is long, technical and requires several auxiliary results, it is thus postponed to Section 5.

An example of AP scheme in the diffusion approximation regime.

The objective of this section is to propose an example of AP scheme for the SDE model (2.12), see Theorem 3.9, in the diffusion approximation regime. The challenge is to let the limiting scheme capture the additional drift term appearing in the limiting equation (2.15) when σ or f is not constant.

Theorem 3.9. Let θ ∈ [ 1 2 , 1].
Introduce the following numerical scheme:

(3.11)                                  m n+1 = m n - ∆tf (X n ) m n+θ 2 + ∆tf (X n )g(X n ) + f (X n )h(X n ) √ ∆tγ n , X n+1 = X n + ∆tb(X n ) + σ(X n ) ∆t m n+θ , m n+1 = m n - ∆tf ( X n+1 )m n+θ 2 + ∆tf ( X n+1 )g(X n ) + f (X n )h(X n ) √ ∆tγ n , Y n+1 = X n + ∆tb(X n ) + σ(X n ) ∆tm n+θ , X n+1 = X n + ∆tb(X n ) + σ(X n ) + σ(Y n+1 ) 2 ∆t m n+θ + m n+θ 2 ,
where

m n+θ = (1 -θ)m n + θ m n+1 m n+θ = (1 -θ)m n + θm n+1 .
This scheme satisfies Assumptions 3.1, 3.2 and 3.3 and is Asymptotic Preserving in the sense of Definition 3.4. Moreover the limiting scheme is given by

(3.12)                          Xn+1 = X n + ∆t (b(X n ) + g(X n )σ(X n )) + σ(X n )h(X n ) √ ∆tγ n , Y n+1 = X n + ∆t b(X n ) + g(X n )σ(X n ) + σ(X n )h(X n ) f (X n ) f ( Xn+1 ) √ ∆tγ n , X n+1 = X n + ∆t b(X n ) + g(X n ) σ(X n ) + σ(Y n+1 ) 2 + σ(X n ) + σ(Y n+1 ) 2 1 + f (Xn) f ( Xn+1) 2 h(X n ) √ ∆tγ n .
The design of the scheme 3.11 is based on a carefully chosen prediction-correction procedure. The limiting scheme (3.12) then also contains prediction steps which are the key elements to satisfy the consistency with the limiting SDE (2.15). The choice of the prediction-correction procedure is made clearer looking at the two examples (2.20) and (2.22), see below Corollaries 3.11 and 3.12 respectively. The prediction-correction procedure is crucial to obtain the AP property for the scheme: the following simpler scheme (with θ = 1 to simplify the presentation)

(3.13)        X n+1 = X n + ∆tb(X n ) + σ(X n ) ∆t m n+1 , m n+1 = m n - f (X n )∆t 2 m n+1 + f (X n )g(X n )∆t + f (X n )h(X n ) √ ∆t γ n ,
is not asymptotic preserving, since the associated limiting scheme (see the proof of Theorem 3.9 for the derivation of the limiting scheme) is

X n+1 = X n + ∆t b(X n ) + g(X n )σ(X n ) + h(X n )σ(X n ) √ ∆tγ n .
This limiting scheme is consistent with the SDE dX t = b(X t ) + g(X t )σ(X t ) dt + h(X t )σ(X t )dW t , which differs in general -when σ or f is non constant -from the correct limiting equation (2.15). Observe that in the AP scheme (3.11) the fast component m is discretized using the θ-method. Choosing θ ∈ [ 1 2 , 1] ensures the mean-square stability of the scheme (Assumption 3.1), uniformly with respect to . Note that the same quantity (1 -θ)m n + θ m n+1 appears in the expressions of m n+1 and X n+1 in (3.11). Similarly, the same quantity (1 -θ)m n + θm n+1 appears in the expressions of m n+1 and Y n+1 in (3.11): this highlights the fact that in order to get a limiting scheme, it is fundamental to choose the quadrature rules in this consistent way.

Remark 3.10. There would be no loss of generality to assume that b = 0. Another example of AP scheme would be obtained in the case b = 0, using a splitting technique: combining the scheme (3.11) with b = 0, with a standard explicit Euler scheme to treat the contribution of b. Writing the expression of the resulting scheme is left to the reader.

Proof of Theorem 3.9. It is straightforward to check that Assumption 3.1 is satisfied.

Let us prove that Assumption 3.3 holds, namely that (3.11) converges to (3.12) when → 0. Note that for fixed ∆t > 0 and 0 ≤ n ≤ N , one has (3.14) sup

>0 E[|m n | + | m n |] < +∞.
This is proved by a straightforward recursion argument. As a consequence, one obtains convergence of the quantity,

∆t m n+θ = ∆tg(X n ) + h(X n ) √ ∆tγ n - f (X n ) m n+1 -m n .
Thus one has

Y n+1 ---→ →0 Y n+1 and X n+1 ---→ →0
Xn+1 . Similarly, one obtains the convergence of ∆tm n+θ , which yields

Y n+1 ---→ →0 Y n+1 and X n+1 ---→ →0 X n+1 .
It remains to prove that the scheme satisfies Assumption 3.2 and is asymptotic preserving in the sense of Definition 3.4, namely that the schemes (3.11) and (3.12) are consistent (in the weak sense), with (2.12) and (2.15) respectively.

On the one hand, let > 0 be fixed. To prove that (3.11) is consistent with (2.12), it is sufficient to prove that, for

ϕ ∈ C 2 b (T d × R), when ∆t → 0, (3.15) E[ϕ(X n+1 , m n+1 )] = E[ϕ(X n , m n )] + ∆tE[L ϕ(X n , m n )] + o(∆t). It is straightforward to check that, in L 1 (Ω), m n+1 , X n+1 , m n+1 , Y n+1 = ( m n , X n , m n , X n ) + o(1), hence m n+1 = m n + √ ∆t f (X n )h(X n )γ n + ∆t - f (X n )m n 2 + f (X n )g(X n ) + o(∆t), X n+1 = X n + ∆t b(X n ) + σ(X n )m n + o(∆t).
Since γ n ∼ N (0, 1) and the random variables γ n and X n are independent, one obtains (3.15). On the other hand, it remains to prove that the limiting scheme (3.12) is consistent with (2.15), i.e. that, for ϕ ∈ C 2 (T d ), when ∆t → 0,

(3.16) E[ϕ(X n+1 )] = E[ϕ(X n )] + ∆tE[Lϕ(X n )] + o(∆t).
To simplify the presentation, for any function ψ, the following notation is used below:

ψ n . = ψ(X n ), ∇ψ n . = ∇ψ(X n )
The key argument of this proof is the analysis of the asymptotic behavior of the quantity σ(Xn)+σ(Yn+1)

2 1+ f (Xn ) f ( Xn+1 )

2

, which appears in the scheme in order to capture the drift terms in the limiting equation (2.15).

First, performing expansions at order √ ∆t for Xn+1 and Y n+1 yields

Xn+1 = X n + √ ∆th n γ n σ n + o( √ ∆t), Y n+1 = X n + √ ∆th n γ n σ n + o( √ ∆t).
Second, writing σn+σ(Yn+1)

2 = σ n + σ(Yn+1)-σn 2 and fn f ( Xn+1) = 1 1+ f ( Xn+1 )-fn fn
, one obtains the following expansion at order √ ∆t for the quantity:

σ n + σ(Y n+1 ) 2 1 + fn f ( Xn+1) 2 = σ n + 1 2 √ ∆th n γ n (σ n • ∇)σ n + o( √ ∆t) 1 - 1 2f n √ ∆th n γ n σ n • ∇f n + o( √ ∆t) = σ n + √ ∆t - 1 2f n h n γ n σ n • ∇f n σ n + 1 2 h n γ n (σ n • ∇)σ n + o( √ ∆t).
Finally, one obtains the following asymptotic expansion of X n+1

X n+1 = X n + √ ∆th n γ n σ n + ∆t b n + g n σ n - 1 2f n h 2 n γ 2 n σ n • ∇f n σ n + 1 2 h 2 n γ 2 n (σ n • ∇)σ n + o(∆t).
Since γ n is centered and γ n and X n are independent random variables, one obtains the first order expansion (3.16). This concludes the proof of Theorem 3.9.

The proposed AP scheme given by Theorem 3.9 can be simplified when it is applied to one of the two examples of SDE models introduced in Section 2.2.2. These two examples are employed in the numerical experiments below. To simplify the presentation, we only consider the case θ = 1.

Corollary 3.11. Consider the SDE (2.20). The AP scheme (3.11) given by Theorem 3.9 is written as follows:

(3.17)

               m n+1 = m n - ∆t 2 m n+1 + √ ∆t γ n , Y n+1 = X n + σ(X n ) ∆tm n+1 , X n+1 = X n + σ(X n ) + σ(Y n+1 ) 2 ∆tm n+1 .
The scheme (3.17) is Asymptotic Preserving, and the limiting scheme (3.12) is written as

(3.18)    Y n+1 = X n + √ ∆tσ(X n )γ n , X n+1 = X n + √ ∆t σ(X n ) + σ(Y n+1 ) 2 γ n .
The limiting scheme (3.18) is consistent with the limiting SDE (2.21).

The prediction-correction procedure appearing in the limiting scheme (3.18) allows to recover the Stratonovich interpretation of the noise in the limiting SDE (2.21).

Corollary 3.12. Consider the SDE (2.22). The AP scheme (3.11) given by Theorem 3.9 is written as follows:

(3.19)                        m n+1 = m n - ∆tf (X n ) 2 m n+1 + ∆tf (X n )g(X n ) + f (X n )h(X n ) √ ∆tγ n X n+1 = X n + ∆t m n+1 m n+1 = m n - ∆tf ( X n+1 ) 2 m n+1 + ∆tf ( X n+1 )g(X n ) + f (X n )h(X n ) √ ∆tγ n X n+1 = X n + ∆t m n+1 + m n+1 2
The scheme (3.19) is Asymptotic Preserving, and the limiting scheme (3.12) is written as

(3.20)      Xn+1 = X n + ∆tg(X n ) + h(X n ) √ ∆tγ n , X n+1 = X n + ∆tg(X n ) + 1 + f (Xn) f ( Xn+1) 2 h(X n ) √ ∆tγ n .
The limiting scheme (3.20) is consistent with the limiting SDE (2.23).

The prediction-correction procedure appearing in the limiting scheme (3.18) allows to recover the noise-induced drift term appearing in the limiting SDE (2.23).

Remark 3.13. Consider the first example (2.20), with σ(x) = x and assume that x takes values in the real line R (instead of the torus T). The following scheme (3.21)

       X n+1 = X n exp ∆t [(1 -θ)m n + θm n+1 ] m n+1 = m n - ∆t 2 (1 -θ)m n + θm n+1 + √ ∆t γ n ,
where θ ∈ [ 1 2 , 1], is another example of AP scheme for this problem. The limiting scheme is given by

(3.22) X n+1 = X n exp( √ ∆tγ n )
which is consistent with the limiting equation dX t = X t • dW t . However, the construction is more subtle if the fast component is discretized using an exponential method: let

m n+1 = e -∆t 2 m n + 1 -e -2 ∆t 2 2 γ n , then m n+1 → →0 1 √
2 γ n , and defining

X n+1 = X n exp ∆t m n+1
does not provide an AP scheme, since there exists no limiting scheme when → 0. Inspired by the identity 1 t 0 m (s)ds = (m 0 -m (t)) + β t for all t ≥ 0, one may set

X n+1 = X n exp (m n -m n+1 ) + √ ∆tγ n ,
so that one has

X N = x 0 exp (m 0 -m N ) + √ ∆t N -1 n=0 γ n .
When → 0, one obtains the limiting scheme (3.22) again. However, the generalization of this construction to the case σ(x) = x is not straightforward, whereas the scheme proposed in Theorem 3.9 applies directly to the general case.

Numerical experiments.

To simplify the discussion, the dimension d is set equal to 1.

4.1. Illustration in the averaging regime. The objective of this section is to illustrate qualitatively the superiority of the AP scheme (3.4) proposed in Section 3.2, when the parameter is small, compared with the use of crude integrators which are not AP. In particular, the numerical experiments below confirm that the limiting scheme (3.5) is consistent with the limiting equation (2.7).

We consider the equation (2.1) with a drift given by b(x, m) = cos(2πx)e -m 2 2 , and diffusion coefficient σ(x, m) = 0. Let T = 1, x 0 = 1 and m 0 = 0.

Recall that the AP scheme is given by (3.4), the limiting scheme is given by (3.5) and the limiting equation is given by (2.7), with b(x) = 1 √ 2 cos(2πx) and σ = 0. Let us define X ref n using the standard Euler scheme applied to this limiting equation:

(4.1) X ref n+1 = X ref n + b(X ref n )∆t.
The scheme (4.1) plays the role of a reference scheme to illustrate the consistency of the limiting scheme (3.5) with the limiting equation, and to illustrate the fact that the crude scheme defined by (3.7) fails to capture the correct limit and is not AP.

In Figure 1, we represent the evolution of X n , X n and X ref n as time t n = n∆t evolves, with ∆t = 0.004 and for different values of . In Figure 1(a), X n and X n are computed using the AP scheme (3.4) and the limit scheme (3.5), while in Figure 1(b), X n is computed using the crude scheme (3.7). Observe that, in both case, the scheme converges when → 0 and that the AP scheme (3.4) does capture the correct limiting equation only with AP scheme (3.4), as opposed to the crude scheme (3.7). In Figure 2, we represent the evolution of X n and X ref n as time t n = n∆t evolves, with ∆t = 0.004 and = 0.001, when X n is computed using the AP scheme (3.4) or the crude scheme (3.7). It illustrates the superiority of the AP scheme over the crude scheme for a small .

Illustration in the diffusion approximation regime.

As in the previous section, the objective of this section is to illustrate qualitatively the superiority of the AP scheme (3.11) proposed in Section 3.3, when the parameter is small, compared a not AP scheme.

The two examples described in Sections 2.2.2 are considered below.

First example.

Let us consider the first example, see Equation (2.20). The diffusion coefficient is given by σ(x) = cos(2πx). Let T = 1, x 0 = 1 and m 0 = 0.

Recall that the AP scheme derived from the general case (3.11) in this case is given by (3.17), the limiting scheme is given by (3.18) and the limiting equation is given by (2.21). Let us define X ref n using the standard Euler-Maruyama scheme applied to this limiting equation (rewritten in Itô form):

(4.2) X ref n+1 = X ref n + 1 2 σ(X ref n )σ (X ref n )∆t + σ(X ref n ) √ ∆tγ n .
The scheme (4.2) plays the role of a reference scheme to illustrate the consistency of the limiting scheme (3.18) with the limiting equation, and to illustrate the fact that the crude scheme defined by (3.13) fails to capture the correct limit and is not AP.

In Figure 3, we represent the evolution of X n and X ref n as time t n = n∆t evolves, with ∆t = 0.004 and for different values of . The discretization X n is computed using the AP scheme (3.17) in Figure 3(a) and the crude scheme (3.13) in Figure 3(b). Observe that, in both case, the scheme seem to converge when → 0 but only the AP scheme (3.17) captures the correct limiting. In Figure 4, we represent the evolution of X n and X ref n as time t n = n∆t evolves, with ∆t = 0.004 and = 0.01, when X n is computed using the AP scheme (3.17) or the crude scheme (3.13). Note how the behavior of the crude scheme differs from the reference. It reveals the superiority of the AP scheme for a small . Fig. 4. Evolution of the AP scheme (3.17), the crude scheme (3.13) and the reference scheme (4.2) (diffusion approximation regime, first example), with ∆t = 0.004 and = 0.01. 4.2.2. First example with σ(x) = x. In this section, we illustrate the performance of the AP scheme presented in Remark 3.13, and an important feature of all the AP schemes presented in this article, concerning the consistency of quadrature rules for discretizations of the fast component.

As explained in Remark 3.13, when σ(x) = x, where x ∈ R belongs to the real line instead of imposing periodic conditions, another type of AP scheme (3.21) can be designed. The limiting equation is dX t = X t • dW t or, with an Itô convention, dX t = 1 2 X t dt + X t dW t , and the Euler-Maruyama scheme (used as a reference scheme) for this limiting equation is written as

(4.3) X ref n+1 = X ref n + 1 2 X ref n ∆t + X ref n √ ∆tγ n .
Recall that in (3.21), the quadrature rule used to discretize the integral in the exponential is closely related to the choice of the scheme for the discretization of the fast component. Let us introduce the following scheme where the consistency is not satisfied (scheme (3.21) corresponds to θ = θ below): 5, we represent the evolution of X n and X ref n as time evolves, with ∆t = 0.004 and = 0.01. In Figure 5(a) X n is computed either the specific AP scheme (3.21) or the general AP scheme (3.11), while in Figure 5(b), it is computed using the scheme (4.4) above with θ = 1 = θ = 0.5. It illustrates the AP property of both schemes (3.21) and (3.17) and the non convergence when the quadrature rules are not chosen consistently.

(4.4)        X n+1 = X n exp ∆t [(1 -θ)m n + θm n+1 ] m n+1 = m n - ∆t 2 (1 -θ )m n + θ m n+1 + √ ∆t γ n , with θ, θ ∈ [ 1 2 , 1]. In Figure

Second example.

Let us now consider the second example described in Section 2.2.2, see Equation (2.22). The coefficients are given by f (x) = cos(2πx)+1.5, g(x) = 0 and h(x) = 1. Let T = 1, x 0 = 1 and m 0 = 0. The general case (3.11) gives in this case the AP scheme (3.19) and the limiting scheme (3.20), whereas the limiting equation is given by (2.23). The reference scheme is obtained by using the standard Euler-Maruyama scheme applied to the limiting equation:

(4.5) X ref n+1 = X ref n -∆t h(X ref n ) 2 f (X ref n ) 2f (X ref n ) + √ ∆tγ n .
We represent in Figure 6 the evolution of X n and X ref n as time evolves, with ∆t = 0.004 and = 0.01, where X n is computed using the AP scheme (3.19) (left) and the crude scheme (3.13) (right). Observe that the AP scheme captures the correct limiting equation when → 0, whereas the crude scheme does not. 5. Proof of Theorem 3.8. The objective of this section is to prove the error estimate (3.9). The proof follows from proving the following four auxiliary lemmas. In their statements, let Assumptions 2.1, 2.2 and 2.3 be satisfied. Let T ∈ (0, ∞) be fixed and assume that ϕ : T d → R is of class C 4 . Recall that the identity T = N ∆t is assumed to hold. In addition, recall that X(t) t≥0 and X n n≥0 are defined by the limiting equation (2.7) and the limiting scheme (3.5) respectively. Lemma 5.6. Define u(t, x) = E x [ϕ(X(t))], for all t ≥ 0 and x ∈ T d , where X(t) t≥0 is the solution of the SDE (2.7) and E x means that X(0) = x. One has u ∈ C([0, T ], C 4 b (T d , R)). In addition, there exists C(T, ϕ) ∈ (0, ∞) such that for all j ∈ {1, 2, 3, 4}, one has (5.6) sup

(t,x)∈[0,T ]×T d D j x u(t, x) ≤ C(T, ϕ). Lemma 5.7. Let ∆t ∈ (0, ∆t 0 ]. Define u n (x) = E x [ϕ(X n )]
, for all n ∈ N and x ∈ T d , where X n n≥0 is defined by the limiting scheme (3.5) (see Theorem 3.7), and E x means that X 0 = x. For all n ≥ 0, one has u n ∈ C 2 (T d ). In addition, there exists C(T, ϕ) ∈ (0, ∞) such that for all j ∈ {1, 2}, one has

(5.7) sup 0≤n≤N sup x∈T d D j x u n (x) ≤ C(T, ϕ)
and, for all ∆t ∈ (0, 0 ], one has

(5.8) sup 0≤n≤N -1 sup x∈T d D j x u n+1 (x) -D j x u n (x) ≤ C(T, ϕ)∆t.
Based on the auxiliary results stated above, the proof of Theorem 3.8 is straightforward.

Proof of Theorem 3.8. Note that

E[ϕ(X N ) -E[ϕ(X (T ))] ≤ E[ϕ(X(T )) -E[ϕ(X (T ))] + E[ϕ(X N ) -E[ϕ(X(T ))] + E[ϕ(X N ) -E[ϕ(X N )] ,
thus combining (5.2), (5.3) and (5.4) (with max(∆t, ) ≤ ∆t + ) yields

E[ϕ(X N ) -E[ϕ(X (T ))] ≤ C(T, ϕ) ∆t + .
Combining that error estimate with (5.1) then concludes the proof of the error estimate (3.9). As already explained above, the error estimate (3.10) is a straightforward consequence of (3.9) (considering the cases √ ∆t ≤ and √ ∆t ≥ ). This concludes the proof of Theorem 3.8.

Let us now give proofs of the auxiliary lemmas 5.1, 5.2, 5.3 and 5.4, employing the results of Lemmas 5.5, 5.6 and 5.7 (proofs are given below).

The following notation is used below in the proofs of the auxiliary results: for all λ, µ ≥ 0, λ µ means that there exists C(T, ϕ) ∈ (0, +∞), independent of ∆t, and n, such that λ ≤ C(T, ϕ)µ. In addition, the following notation is used for the infinitesimal generator of the Ornstein-Uhlenbeck process:

L x OU = -m∂ m + h(x) 2 ∂ 2 m ,
in order to let the dependence with respect to x be clear.

Proof of Lemma 5.1. Let us introduce auxiliary continuous-time processes X and m , such that, for all n ∈ {0, . . . , N }, one has X n = X (t n ) and m n = m (t n ) (recall the definition (3.4) of the scheme): for

t n ≤ t ≤ t n+1 X (t) = X n + (t -t n )b(X n , m (t)) + σ(X n , m (t)) (B(t) -B(t n )) d m t = - m t dt + √ 2h(X n ) √ dβ t .
Note that m does satisfy m n = m (t n ), since m n is exact in distribution and m is an Ornstein-Uhlenbeck process with variance

2h(X n ) 2
. The process X satisfies on each subinterval [t n , t n+1 ] the following stochastic differential equation: for all t ∈ [t n , t n+1 ], one has

d X t = b(X n , m t )dt + σ(X n , m t )dB t + (t -t n )∂ m b(X n , m t )d m t + ∂ m σ(X n , m t )(B t -B tn )d m t + h(X n ) 2 (t -t n )∂ 2 m b(X n , m t )dt + h(X n ) 2 ∂ 2 m σ(X n , m t )(B t -B tn )dt.
The expressions for X are complicated due to the fact that in the scheme (3.4), b and σ are evaluated with m = m n+1 , which is required to satisfy the AP property.

Owing to Lemma 5.5, the auxiliary function u is of class C 2 and is solution of the Kolmogorov equation ∂ t u = L u . Using a telescoping sum argument and the definition of the auxiliary processes X and m , the application of Itô's formula yields the following standard expression for the weak error:

E[ϕ(X N )]-E[ϕ(X (T ))] = E[u (0, X N , m (t N ))] -E[u (T, X 0 , m (0))] = N -1 n=0 E[u (T -t n+1 , X (t n+1 ), m (t n+1 )) -u (T -t n , X (t n ), m (t n ))] = N -1 n=0 tn+1 tn E[(-∂ t + L n )u (T -t, X (t), m (t))]dt = N -1 n=0 tn+1 tn E[( L n -L )u (T -t, X (t), m (t))]dt,
where the auxiliary differential operator L n is such that

L n u (T -t, X (t), m (t)) = b(X n , m (t)) • ∇ x u (T -t, X (t), m (t)) + 1 2 σσ * (X n , m (t)) : ∇ 2 x u (T -t, X (t), m (t)) + 1 L X n OU u (T -t, X (t), m (t)) + r n (t),
and where the remainder term r n (t) is given by

r n (t) = 1 (t -t n )L X n OU b(X n , m (t)) • ∇ x u + 1 L X n OU σ(X n , m (t))(B(t) -B(t n )) • ∇ x u + 1 2 D 2 x u • ((t -t n )∂ m b(X n , m (t)) + ∂ m σ(X n , m (t))(B(t) -B(t n ))) 2 ,
where u and its derivatives are evaluated at (T -t, X (t), m (t)), and where we used the notation D 2 φ(x) • (y) 2 = yy * : ∇ 2 x φ(x) to simplify the presentation. Let us first deal with the remainder term r n (t). Observe that the processes B(t) -B(t n ) t∈[tn,tn+1] and X n , m (t) t∈[tn,tn+1] are independent, thus using a conditioning argument and the regularity estimates from Lemma 5.5, one has for all t ∈ [t n , t n+1 ]

E[r n (t)] t -t n .
It remains to deal with (5.9)

( L n -L )u (T -t, X (t), m (t)) -r n (t) = b(X n , m (t)) -b( X (t), m (t)) • ∇ x u (T -t, X (t), m (t)) + 1 2 σσ * (X n , m (t)) -σσ * ( X (t), m (t)) : ∇ 2 x u (T -t, X (t), m (t)) + h(X n ) 2 -h( X (t)) 2 ∂ 2 m u (T -t, X (t), m (t)),
where the expressions (2.3) and (2.4) for the infinitesimal generator L have been used. The three quantities appearing in the right-hand side of (5.9) above are of the type

V ( X (t), m (t)) -V (X n , m (t)) U ( X (t), m (t)), for V = b, σσ * or h 2 and U = ∇ x u (T -t), ∇ 2 x u (T -t) or ∂ 2 m u (T -t).
Using again the independence of (B(t) -B(t n )) tn≤t≤tn+1 and X n , m (t) t∈[tn,tn+1] and regularity properties of u given in Lemma 5.5, applying Itô's formula and conditioning with respect to (X n , ( m (t)) tn≤t≤tn+1 ), one obtains

E V ( X (t), m (t)) -V (X n , m (t)) U (X n , m (t)) t -t n . Moreover, V C 1 (T d ,C(R))
1 and, owing to Lemma 5.5, U C 1 (T d ,C(R)) 1. Therefore, we have

E V ( X (t), m (t)) -V (X n , m (t)) U ( X (t), m (t)) -U (X n , m (t)) ≤ V (•, m (t)) C 1 (T d ) E U (•, m (t)) C 1 (T d ) X (t) -X n 2 E X (t) -X n 2 t -t n ,
where the last inequality comes from the definition of X (t). Gathering the two estimates above gives, for all t ∈ [t n , t n+1 ],

E V ( X (t), m (t)) -V (X n , m (t)) U ( X (t), m (t)) t -t n ,
and using (5.9) finally yields, for all t ∈ [t n , t n+1 ],

E ( L n -L )u (T -t, X (t), m (t)) t -t n + t -t n ∆t .
One then obtains

N -1 n=0 tn+1 tn E[( L n -L )u (T -t, X (t), m (t))]dt ∆t ,
which concludes the proof of (5.1) and of Lemma 5.1.

Proof of Lemma 5.2. For all x ∈ T d , introduce the auxiliary Ornstein-Uhlenbeck process m x solving the SDE

dm x t = -m x t dt + √ 2h(x)dβ t
Let m x (t, m) denote the solution at time t, if the initial condition is given by m x (0, m).

The invariant distribution of the process m x is equal to ν x . Note that for m and m ∈ R, one has

m x (t, m) -m x (t, m ) = (m -m )e -t . Consider V = b or V = σσ , and let
δ(t, x, m) = E[V (x, m (t, m)) -V (x)].
Note that for m and m ∈ R, one has m x (t, m) -m x (t, m ) = (m -m )e -t . As a consequence, we have

δ(t, x, m) -δ(t, x, m ) m -m e -t .
By integrating with respect to m and using the equality

V (x) = V (x, m )dν x (m ), one obtains (5.10) δ(t, x, m) (1 + m )e -t .
Using the fact that δ and its derivatives satisfy (5.10) with V = b, one is able to check that the function ψ b given by, for all x ∈ T d and m ∈ R,

ψ b (x, m) = - ∞ 0 E[b(x, m x (t, m)) -b(x)]dt
is well-defined (by definition of b, see (2.5)) and is of class C 2 (T d × R). Moreover, ψ b and its derivative have at most linear growth in m. In addition,

ψ b (x, •) solves the Poisson equation L x OU ψ b (x, m) = b(x, m) -b(x)
, for all x ∈ T d (indeed ν x is the invariant distribution associated with the generator L x OU of the Ornstein-Uhlenbeck process m x ).

Similarly, define, for all x ∈ T d and m ∈ R,

ψ σ (x, m) = - ∞ 0 E[σσ * (x, m x (t, m)) -σ σ * (x)]dt.
The function ψ σ is well-defined: owing to (2.6) (Assumption 2.2) one has the equality σσ * (x, m)dν x (m) = σ σ * (x) for all x ∈ T d , and using the same arguments as above,

ψ σ (x, •) solves the Poisson equation L x OU ψ σ (x, m) = σσ * (x, m) -σ σ * (x), for all x ∈ T d . Now, for all t ∈ [0, T ], x ∈ T d and m ∈ R, let Φ(t, x, m) = ψ b (x, m) • ∇ x u(T -t, x) + ψ σ (x, m) : ∇ 2 x u(T -t, x),
where u is given by Lemma 5.6 and solves the Kolmogorov equation ∂ t u = Lu. On the one hand, applying Itô's formula yields the following expression for the error term:

E[ϕ(X (T ))]-E[ϕ(X(T ))] = E[u(0, X (T )) -u(T, X (0))] = T 0 E[ b(X (t), m (t)) -b(X (t)) • ∇ x u(T -t, X (t))]dt + 1 2 T 0 E[(σσ * (X (t), m (t)) -σ σ * (X (t))) : ∇ 2 x u(T -t, X (t))]dt = T 0 E[L x OU Φ(t, X (t), m (t))]dt,
by definition of the auxiliary function Φ, since ψ b and ψ σ are solutions of Poisson equations.

On the other hand, applying Itô's formula also gives the identity

E[Φ(T, X (T ), m (T ))] -E[Φ(0, X (0), m (0))] = T 0 E[(∂ t + b • ∇ x + 1 2 σσ * : ∇ 2 x + 1 L x OU )Φ(t, X (t), m (t))]dt.
Combining the two expressions then gives

E[ϕ(X (T ))]-E[ϕ(X(T ))] = T 0 E[L x OU Φ(t, X (t), m (t))]dt = (E[Φ(T, X (T ), m (T ))] -E[Φ(0, X (0), m (0))]) - T 0 E[(∂ t + b • ∇ x + 1 2 σσ * : ∇ 2 x )Φ(t, X (t), m (t))]dt.
Using the regularity estimates from Lemma 5.6 and the identity

∂ t u = b • ∇ x u + σ σ * : ∇ 2 x u,
it is then straightforward to obtain (5.2). This concludes the proof of Lemma 5.2.

Proof of Lemma 5.3. Let us introduce the continuous-time auxiliary process X, such that, for t ∈ [t n , t n+1 ], one has

X(t) = X n + (t -t n )b(X n , h(X n )γ n ) + σ(X n , h(X n )γ n )(B(t) -B(t n )).
Introduce also the second-order differential operator Ln

= b(X n , h(X n )γ n ) • ∇ + 1 2 σσ * (X n , h(X n )γ n ) : ∇ 2 .
With this notation, for any function φ ∈ C 2 (T d ), Itô's formula gives, for all t ∈ [t n , t n+1 ], (5.11)

φ( X(t)) -φ(X n ) = t tn Ln φ( X(s))ds + t tn ∇φ( X(s)) • σ(X n , h(X n )γ n )dB s .
Using the same (standard) arguments as in the proof of Lemma 5.1, one obtains the following decomposition of the error:

E[ϕ(X N )] -E[ϕ(X(T ))] = E[u(0, X N )] -E[ϕ(T, X 0 )] = N -1 n=0 tn+1 tn E[u(T -t n+1 , X n+1 ) -u(T -t n , X n )] = N -1 n=0 tn+1 tn E[(-∂ t + Ln )u(T -t, X(t))]dt = N -1 n=0 tn+1 tn E b(X n , h(X n )γ n ) -b( X(t)) • ∇ x u(T -t, X(t)) dt + 1 2 N -1 n=0 tn+1 tn E σσ * (X n , h(X n )γ n ) -σ σ * ( X(t)) : ∇ 2 x u(T -t, X(t)) dt. The error term E b(X n , h(X n )γ n ) -b( X(t)) • ∇ x u(T -t, X(t)) (with V = b and U = ∇ x u) and E σσ * (X n , h(X n )γ n ) -σ σ * ( X(t)) : ∇ 2 x u(T -t, X(t)) (with V = σσ * and U = ∇ 2
x u) are written as

V (X n , h(X n )γ n ) -V ( X(t)) U (T -t, X(t)).
Note that X(t n ) = X n . As a consequence,

E V (X n , h(X n )γ n ) -V ( X(t)) U (T -t, X(t)) = δ 1 n (t) + δ 2 n (t) + δ 3 n (t) + δ 4 n (t),
where

δ 1 n (t) = E V (X n , h(X n )γ n ) -V (X n ) U (T -t, X n ) δ 2 n (t) = E V (X n ) -V ( X(t)) U (T -t, X n ) δ 3 n (t) = E V (X n , h(X n )γ n ) -V (X n ) U (T -t, X(t)) -U (T -t, X n ) δ 4 n (t) = E V (X n ) -V ( X(t)) U (T -t, X(t)) -U (T -t, X n ) .
It remains to treat the four error terms δ j n (t), j = 1, 2, 3, 4. Let us start with the most important observation: by definition of V (x) = V (x, m)dν x (m), the independence of the random variables X n and γ n yields the identity

E V (X n , h(X n )γ n ) -V (X n ) U (T -t, X n ) = E E V (X n , h(X n )γ n ) -V (X n ) U (T -t, X n ) | X n = 0.
The fact that this term vanishes is fundamental since it justifies the consistency of the scheme (3.5) with the limiting equation (2.7) (see also Theorem 3.7 and its proof), and the AP property.

To treat the second term, observe that X n and ( X(s), B(s)) s>tn are independent random variables, thus conditioning with respect to X n and applying Itô's formula (5.11) gives

E V (X n ) -V ( X(t)) U (T -t, X n ) = E t tn Ln V ( X(s))ds U (T -t, X n ) .
Since U is bounded (owing to the regularity estimates from Lemma 5.6) and since V ∈ C 2 (T d ) (by assumptions on the coefficients b and σ, see Assumption 2.2), one obtains

E V (X n ) -V ( X(t)) U (T -t, X n ) t -t n .
The treatment of the third term uses a conditioning argument, and Itô's formula (5.11): one has

E V (X n , h(X n )γ n ) -V (X n ) U (T -t, X(t)) -U (T -t, X n ) | X n , γ n = E V (X n , h(X n )γ n ) -V (X n ) t tn Ln U (T -t, X(s))ds | X n , γ n .
Using the regularity properties from Lemma 5.6, one obtains

E V (X n , h(X n )γ n ) -V (X n ) U (T -t, X(t)) -U (T -t, X n ) t -t n .
The treatment of the fourth error term is straightforward: since U and V are Lipschitz continuous (owing to Lemma 5.6), one has

E U (T -t, X(t)) -U (T -t, X n ) V (X n ) -V ( X(t)) E X(t) -X n 2 t -t n .
The estimates above are of the type

δ j n (t) t -t n ,
for all t ∈ [t n , t n+1 ] and j = 1, 2, 3, 4. Finally, one obtains

E[ϕ(X N )] -ϕ(X(T )) N -1 n=0 tn+1 tn (t -t n )dt ∆t,
which concludes the proof of Lemma 5.3.

Proof of Lemma 5.4. The idea is to adapt the proof of Lemma 5.2 (see Section A.1) to the discrete-time situation. Let us start with preparatory computations. A telescoping sum argument yields the equality

E[ϕ(X N )] -E[ϕ(X N )] = E[u 0 (X N )] -u N (X 0 ) = N -1 n=0 E[u N -n-1 (X n+1 )] -E[u N -n (X n )] , (5.12)
where the auxiliary function u n is defined in Lemma 5.7. Using the definition of the scheme (3.4), and Markov property combined with the expression of the limiting scheme (3.5), one obtains

u N -n-1 (X n+1 ) = u N -n-1 X n + ∆tb(X n , m n+1 ) + √ ∆tσ(X n , m n+1 )Γ n , E[u N -n (X n )] = E u N -n-1 X n + ∆tb(X n , h(X n )γ n ) + √ ∆tσ(X n , h(X n )γ n )Γ n .
Then, for all m and m ∈ R, and all n ∈ N, one has

m x n (m) -m x n (m ) = e -nτ (m -m ).
Observe that since b is a Lipschitz continuous function, standard arguments give the following upper bound: for all n ∈ N, x ∈ T d and all m, m ∈ R, if δ n (x, m) = (P x τ ) n b(x, m)b(x), then one has

δ n (x, m) = (P x τ ) n b(x, m) -(P x τ ) n b(x, m )dν x (m ) ≤ Eb(x, m x n (m)) -Eb(x, m x n (m ) dν x (m ) e -nτ (1 + |m|).
Similarly, since the derivatives of m x n with respect to x do not depend on m, one can check that the inequality above holds for D x δ and D 2

x δ then concludes the proof of the claim. Since τ ∞ n=0 e -nτ = τ

1-e -τ ≤ max(τ, 1), for all τ ∈ (0, ∞), one obtains inequalities of the type Then ψ τ σ is well-defined and solves the Poisson equation L x τ ψ σ (x, •) = σσ (x, •)σ(x)σ (x), by definition of σ(x), see Assumption 2.3. In addition, ψ τ σ and its derivatives satisfy upper bound of the type (5.14).

Like in the proof of Lemma 5.2 (see Section A.1), introduce the auxiliary function defined by

Φ n (x, m) = ψ τ b (x, m) • ∇ x u N -n-1 (x) + ψ τ σ (x, m) : ∇ 2 x u N -n-1 (x),
for all n ∈ {0, . . . , N -1}, x ∈ T d and m ∈ R, where u n is defined in Lemma 5.7.

Combining the decomposition of the error (5.12) and the identity (5.13), one then obtains the following new expression for the error:

E[ϕ(X N )] -E[ϕ(X N )] = ∆t 2 N -1 n=0 E[R n (∆t)] + ∆t N -1 n=0 L X n τ Φ n (X n , m n+1 ). (5.15)
On the one hand, a telescoping sum argument yields the following expression:

E P X N τ Φ N (X N , m N ) -P X 0 τ Φ 0 (X 0 , m 0 ) = N -1 n=0 E P X n+1 τ Φ n+1 (X n+1 , m n+1 ) -P X n τ Φ n (X n , m n ) = N -1 n=0 E P X n+1 τ Φ n+1 (X n+1 , m n+1 ) -P X n+1 τ Φ n (X n+1 , m n+1 ) + N -1 n=0 E P X n+1 τ Φ n (X n+1 , m n+1 ) -P X n τ Φ n (X n , m n+1 ) + N -1 n=0 E P X n τ Φ n (X n , m n+1 ) -P X n τ Φ n (X n , m n ) .
Note that using Markov property, the first sum on the right-hand side above can be written as

N -1 n=0 E Φ n+1 (X n+1 , m n+2 ) -Φ n (X n+1 , m n+2 ) .
On the other hand, by definition of the operator L x τ with the parameter τ = ∆t , one obtains

∆tE[L X n τ Φ n (X n , m n+1 )] = E[P X n τ Φ n (X n , m n+1 ) -Φ n (X n , m n+1 )] = E[P X n τ Φ n (X n , m n+1 ) -P X n τ Φ n (X n , m n )].
Finally, combining the two identities above, one obtains the following expression for the error:

E[ϕ(X N )] -E[ϕ(X N )] = ∆t 2 N -1 n=0 E[R n ] + E Φ N (X N , m N +1 ) -Φ 0 (X 0 , m 1 ) - N -1 n=0 E Φ n+1 (X n+1 , m n+2 ) -Φ n (X n+1 , m n+2 ) - N -1 n=0 E P X n+1 τ Φ n (X n+1 , m n+1 ) -P X n τ Φ n (X n , m n+1 ) .
It then remains to use auxiliary upper bounds to deduce the result, in particular using (5.14). Note that max(τ, 1) = max(∆t, ).

• as explained above, E[ R n (∆t) ] 1, thus the first term satisfies

E[ ∆t 2 N -1 n=0 R n ] ∆t.
• Using the upper bound (5.7) from Lemma 5.7 and (5.14) with ψ = ψ τ b and ψ = ψ τ σ , the second term satisfies

E Φ N (X N , m N +1 ) -Φ 0 (X 0 , m 1 ) max(∆t, ).
• Observe that one has for all x ∈ T d and x ∈ R,

Φ n+1 (x, m) -Φ n (x, m) = ψ τ b (x, m) • (∇ x u N -n (x) -∇ x u N -n-1 (x)) + ψ τ σ (x, m) : ∇ 2 x u N -n (x) -∇ 2 x u N -n-1 (x) .
Using the upper bound (5.8) from Lemma 5.7 and (5.14), one has

Φ n+1 (x, m) -Φ n (x, m) ∆t max(τ, 1)(1 + |m|).
As a consequence, owing to (3.8) and Assumption 2.1, the third term satisfies

N -1 n=0 E Φ n+1 (X n+1 , m n+2 ) -Φ n (X n+1 , m n+2 ) N -1 n=0 ∆t max(τ, 1)
max(∆t, ).

• Note that Φ n and its derivatives satisfy the upper bound (5.14), owing to (5.7) from Lemma 5.7. Let

f τ n (x, m) = P x τ Φ n (x, m) = E γ∼N (0,1) [Φ n (x, e -τ m + 1 -e -2τ h(x)γ)].
It is straightforward to check that f τ n is twice differentiable. In addition, f τ n and its derivatives satisfy (5.14). Using a second order Taylor expansion, one obtains

P X n+1 τ Φ n (X n+1 , m n+1 ) -P X n τ Φ n (X n , m n+1 ) = f τ n (X n+1 , m n+1 ) -f τ n (X n , m n ) ∆t max(τ, 1)(1 + m ).
Finally, using (3.8) and Assumption 2.1, one obtains

N -1 n=0 E P X n+1 τ Φ n (X n+1 , m n+1 ) -P X n τ Φ n (X n , m n+1 ) max(∆t, ).
Gathering the estimates then concludes the proof of Lemma 5.4.

We refer to [START_REF] Cerrai | Second order PDE's in finite and infinite dimension[END_REF]Theorem 1.3.6] for the proof of Lemma 5.6. It thus remains to provide the arguments for the proofs of Lemmas 5.5 and 5.7. The strategy is standard in the literature, see for instance [START_REF] Cerrai | Second order PDE's in finite and infinite dimension[END_REF]. As a consequence, in order to reduce the length of the manuscript, below the details are only given for the proofs of the estimates for first-order derivatives.

= (k x , k m ) ∈ R d × R, one has D x,m u (t, x, m) • k = E X (0)=x,m (0)=m [Dϕ(X (t)) • η ,k x (t)],
where the process η ,k = (η ,k

x , η ,k m ) is the solution of the first variation equation associated with (2.1): for all t ≥ 0,

     dη ,k x,t = Db(X t , m t ) • η ,k t dt + Dσ(X t , m t ) • η ,k t dB t dη ,k m,t = - η ,k m,t dt + √ 2Dh(X t ) • η ,k x,t √ dβ t , with initial conditions η ,k x (0) = k x and η ,k m (0) = k m .
On the one hand, the component η ,k m satisfies the following equality,

η ,k m (t) = e -t k m + 1 √ t 0 e -t-s √ 2Dh(X (s)) • η ,k x (s)dB s ,
and by means of Itô's isometry formula, one has

E[ η ,k m (t) 2 ] ≤ k 2 m + 2 h C 1 (T d ) t 0 e -2(t-s) E[ η ,k x (s) 2 ]ds.
On the other hand, the functions b and σ are globally Lipschitz continuous (see Assumption 2.2), and one has η

,k 2 = η ,k x 2 + η ,k m
2 ; using Itô's isometry formula, and Minkowski's and Young's inequalities, one obtains

E[ η ,k x (t) 2 ] x 2 + t 0 E[ η ,k x (s) 2 ] + E[ η ,k m (s) 2 ] ds.
Combining the two estimates above then yields

sup s∈[0,t] E[ η ,k x (s) 2 ] k 2 + t 0 E[ η ,k x (s) 2 ] + 1 s 0 e -2 s-r E[ η ,k x (r) 2 ]dr ds k 2 + t 0 sup r∈[0,s] E[ η ,k x (r) 2 ]ds,
where the inequality 1 s 0 e -2 s-r dr ≤ 1 2 has been used. Applying Gronwall's lemma, then inserting the result in the estimate above, one obtains the upper bounds

sup t∈[0,T ] E[ η ,k x (t) 2 ] k 2 , sup t∈[0,T ] E[ η ,k m (t) 2 ] k 2 .
Since ϕ is Lipschitz continuous, using the expression for D x,m u (t, x, m) stated above, one finally obtains (5.5) for the first-order derivative: for all k ∈ R d+1 , one has sup

(t,x,m)∈[0,T ]×T d ×R D x,m u (t, x, m) • k k .
The treatment of higher-order derivatives follows from similar arguments which are omitted. This concludes the proof of Lemma 5.5.

Proof of Lemma 5.7. For all k ∈ R d , one has

Du n (x) • k = E x [Dϕ(X n ) • η k n ]
, where η h 0 = h, and for all n ∈ {0, . . . , N -1}, one has

η k n+1 = η k n + ∆tD x b(X n , h(X n )γ n ) • η k n + ∆t∂ m b(X n , h(X n )γ n )Dh(X n ) • η k n γ n + √ ∆tD x σ(X n , h(X n )γ n ) • η k n Γ n + √ ∆t∂ m σ(X n , h(X n )γ n )Dh(X n ) • η k n Γ n γ n .
The functions b, h and σ are Lipschitz continuous (see Assumption 2.2). Since γ n and Γ are independent centered Gaussian random variables, it is straightforward to obtain the upper bound

E[ η k n+1 2 ] (1 + ∆t)E[ η k n 2 ].
A straightforward recursion argument then gives, for all n ∈ {0, . . . , N },

E[ η k n 2 ] k 2
and one obtains (5.7) for the first-order derivative: for all n ∈ {0, . . . , N }, |Du n (x) • k| k .

The treatment of higher-order derivatives would be similar and is omitted. It thus remains to prove (5.8). On the one hand, by definition (3.5), a second order Taylor expansion yields, for all x ∈ T d and n ∈ {0, . . . , N -1},

u n+1 (x) -u n (x) = E x [ϕ(X n+1 ) -ϕ(X n )] ∆t,
On the other hand, one has

Du n+1 (x) • k -Du n (x) • k = E x [Dϕ(X n+1 ) • η k n+1 -Dϕ(X n ) • η k n ] = E[(Dϕ(X n+1 ) -Dϕ(X n )) • (η k n+1 -η k n )] + E[Dϕ(X n ) • (η k n+1 -η k n )] + E[(Dϕ(X n+1 ) -Dϕ(X n )) • η k n ].
It is straightforward to check that one has the inequalities E[ X n+1 -X n 2 ] ∆t and Using a conditional expectation argument, since γ n , Γ n , X n are independent random variables, one has

E[ η k n+1 -η k n 2 ] ∆t k 2 .
E[Dϕ(X n ) • (η k n+1 -η k n )] = ∆t E[Dϕ(X n ) • D x b(X n , h(X n )γ n ) • η k n ] ∆t k .
Finally, using a second-order Taylor expansion and conditioning arguments, one obtains

E[(Dϕ(X n+1 ) -Dϕ(X n )) • η k n ] ∆t k .
As a consequence, one obtains (5.8) when j = 0 and j = 1. This concludes the proof of Lemma 5.7.

Conclusion.

In this article, we have studied a general notion of Asymptotic Preserving schemes, related to convergence in distribution, for a class of SDE systems in averaging and diffusion approximation regimes. Let us mention that some assumptions made to simplify the setting (the slow component takes values in a compact set T d and the fast component is one-dimensional) may easily be relaxed. Note that when the slow component takes values in R d , it is necessary to also study the stability of the numerical schemes, for instance in mean-square sense.

A limitation of our study is the fact that the fast component is an Ornstein-Uhlenbeck process (when the slow component is frozen): even if the general theory of AP schemes described in Section 3.1 holds in more general settings, the construction of implementable AP schemes (such as the ones described in Sections 3.2 and 3.3) is not straightforward if for instance the fast component is solution of a general ergodic SDE with nonlinear coefficients.

We have also left open the question of obtaining a version of the error estimates stated in Theorem 3.8 in the diffusion approximation case. This question will be studied in future works.

Finally, it would be natural to apply the recipes for the design of AP schemes described in this article to SPDE models. For instance, in a future work [START_REF] Bréhier | Asymptotic-preserving schemes for stochastic kinetic equations in the diffusion approximation regime[END_REF], we plan to design, analyze and test AP schemes for the stochastic kinetic PDE model considered in [START_REF] Rakotonirina-Ricquebourg | Diffusion limit for a stochastic kinetic problem with unbounded driving process[END_REF].

Appendix A. Derivation of the limiting models.

A.1. Sketch of proof of Proposition 2.4 (averaging). Let us first give details concerning the construction of the perturbed test function ϕ given by (2.9), such that (2.10) holds. Recall that this construction is used in the statement of Proposition 3.5.

Owing to the multiscale expansions (2.4) and (2.9) of the generator L and of the perturbed test function ϕ = ϕ + ϕ 1 , one has Note that for m and m ∈ R, one has m x (t, m) -m x (t, m ) = (m -m )e -t . As a consequence, we have δ(t, x, m) -δ(t, x, m ) m -m e -t .

By integrating with respect to m and using the equality V (x) = V (x, m )dν x (m ), one obtains (A.2) δ(t, x, m)

(1 + m )e -t .
Since V is of class C 3 with bounded derivatives, and since the derivatives of m x (t, m) with respect to x do not depend on m, it is straightforward to generalize A.2 to the derivatives of δ. It gives that ϕ 1 ∈ C 2 (T d × R) and that ϕ 1 and its derivatives have at most linear growth in m, hence L 0 ϕ 1 also does. This leads to (2.10) using (2.2). This concludes the identification of the limiting generator L using the perturbed test function method. The remaining ingredients of this strategy to prove the convergence in distribution of the process X to the solution X of the limiting equation associated with the limiting generator L are standard and are thus omitted.

A.2. Sketch of proof of Proposition 2.6 (diffusion approximation). Let us first give details concerning the construction of the perturbed test function ϕ given by (2.17), such that (2.18) holds. Recall that this construction is used in the statement of Proposition 3.5.

Owing to the multiscale expansions (2.14) and (2.17) of the generator L and of the perturbed test function ϕ = ϕ + ϕ 1 + 2 ϕ 2 , one has

(A.3) L ϕ = -2 L OU ϕ + -1 (L 1 ϕ + L OU ϕ 1 ) + (L 0 ϕ + L 1 ϕ 1 + L OU ϕ 2 ) + (L 0 ϕ 1 + L 1 ϕ 2 ) + 2 L 0 ϕ 2 .
Since the test function ϕ does not depend on m, one has L OU ϕ = 0, thus the term of order -2 in (A.3) vanishes. Define (A.4) ϕ 1 (x, m) . = m σ(x) f (x)

• ∇ x ϕ(x).

Then it is straightforward to check that L 1 ϕ + L OU ϕ 1 = 0, thus the term of order -1 in (A.3) vanishes. It remains to construct the function ϕ 2 such that the term of order 1 in (A. for some constant C ϕ depending only on ϕ and on the coefficients of the SDE. It remains to check that Lϕ(x) = R (L 0 ϕ + L 1 ϕ 1 )(x, m)dν x (m) gives the expression (2.16). This concludes the identification of the limiting generator L using the perturbed test function method. The remaining ingredients of this strategy to prove the convergence in distribution of the process X to the solution X of the limiting equation associated with the limiting generator L follows from standard arguments which are omitted.

Fig. 1 .

 1 Fig. 1. Evolution of the AP scheme (3.4) (left), the crude scheme (3.7) (right), and the reference scheme (4.1) (averaging regime), with ∆t = 0.004.

Fig. 2 .

 2 Fig. 2. Evolution of the AP scheme (3.4), the crude scheme (3.7) and the reference scheme (4.1) (averaging regime), with ∆t = 0.004 and = 0.001.

Fig. 3 .

 3 Fig.3. Evolution of the AP scheme (3.17) (left), the crude scheme (3.13) (right), and the reference scheme (4.2) (diffusion approximation regime regime, first example), with ∆t = 0.004.

5 Fig. 5 .

 55 Fig. 5. Evolution of the AP schemes (3.21) and (3.17) (left), the crude scheme (4.4) with θ = θ (right), and the reference scheme (4.3) (diffusion approximation regime, first example with σ(x) = x), ∆t = 0.004 and = 0.01.

Fig. 6 .

 6 Fig. 6. Evolution of the AP scheme (3.19), the crude scheme (3.13), and the reference scheme (4.5) for the second example (2.22) (diffusion approximation regime, second example), with ∆t = 0.004 and = 0.01.

  its derivatives D x ψ τ b and D 2 x ψ τ b . Similarly, define for all x ∈ T d and m ∈ R, ψ τ σ (x, m) = -τ +∞ n=0 ((P x τ ) n σσ * (x, m) -σ σ * (x)).

Proof of Lemma 5 . 5 .

 55 Owing to [7, Proposition 1.3.5], for all k

  Since ϕ is Lipschitz continuous, using Cauchy-Schwarz inequality givesE[(Dϕ(X n+1 (x)) -Dϕ(X n (x))) • (DX n+1 (x) • k -DX n (x) • k)]∆t k .

(A. 1 ) 0 E 1 .

 101 L ϕ = -1 L OU ϕ + (L 0 ϕ + L OU ϕ 1 ) + L 0 ϕ 1 .Since the test function ϕ does not depend on m, one has L OU ϕ = 0, thus the term of order -1 in (A.1) vanishes.Define, for all x ∈ T d and m ∈ R,Lϕ(x) = R L 0 ϕ(x, m)dν x (m) = b(x) • ∇ x ϕ(x) + σ σ * (x) : ∇ 2 x ϕ(x) ϑ(x, m) = L 0 ϕ(x, m) -Lϕ(x) = b(x, m) -b(x) • ∇ x ϕ(x) + (σσ * (x, m) -σ σ * (x)) : ∇ 2 x ϕ(x),where we recall that ν x = N (0, h(x) 2 ) is the invariant distribution of the ergodic Ornstein-Uhlenbeck process m x associated to L OU on R, for any fixed x ∈ T ddm x t = -m x t dt + √ 2h(x)dβ t .Let m x (t, m) denote the solution at time t, if the initial condition is given by m x (0, m). Therefore, the centering condition ϑ(x, m)dν x (m) = 0 is satisfied and the Poisson equation-L OU ϕ 1 (x, •) = ϑ(x, •) admits a solution ϕ 1 (x, m) = ∞ [ϑ(x, m x (t, m))]dt.The multiscale expansion (A.1) becomesL ϕ = Lϕ + L 0 ϕTo prove (2.10), it only remains to get estimates on L 0 ϕ 1 uniformly in . Consider V = b or V = σσ * and let δ(t, x, m) . = E[V (x, m x (t, m)) -V (x)].

2 )

 2 3) is equal to Lϕ. Define, for all x ∈ T d and m ∈ R,Lϕ(x) = R (L 0 ϕ + L 1 ϕ 1 )(x, m)dν x (m), ϑ(x, m) = (L 0 ϕ + L 1 ϕ 1 )(x, m) -Lϕ(x),where we recall that ν x = N (0, f (x)h(x) 2 is the invariant distribution of the ergodic Ornstein-Uhlenbeck process associated toL OU on R, for any fixed x ∈ T d . Let x ∈ T d , then the Poisson equation -L OU ϕ 2 (x, •) = ϑ(x,•) admits a solution ϕ 2 , since the centering condition ϑ(x, m)dν x (m) = 0 is satisfied. Precisely, one has the expressionsϑ(x, m) = -|m| 2x ϕ . (A.5)With the functions ϕ 1 and ϕ 2 constructed above, the multiscale expansion (A.3) is rewritten asL ϕ = Lϕ + (L 0 ϕ 1 + L 1 ϕ 2 ) + 2 L 0 ϕ 2 ,which gives (2.18), more precisely sup x∈T d |L ϕ (x, m) -Lϕ(x)| ≤ C ϕ |m| + 2 |m| 2 ,

  It remains to prove that the scheme satisfies Assumption 3.2 and is asymptotic preserving in the sense of Definition 3.4, namely that the schemes (3.4) and (3.5) are consistent (in the weak sense), with (2.1) and (2.7) respectively.

m converges almost surely to h(x)γ, thus Φ ∆t (x, m, Γ, γ) converges in distribution to Φ ∆t (x, Γ, γ), and Assumption 3.3 is satisfied. Let > 0 be fixed. Since h is bounded, it is straightforward to check that

(3.8)

sup

∈(0,1]

Acknowledgments. The work of C.-E. B. is partially supported by the following projects operated by the French National Research Agency: ADA (ANR-19-CE40-0019-02 ), BORDS (ANR-16-CE40-0027-01) and SIMALIN (ANR-19-CE40-0016).

Lemma 5.1. There exists C(T, ϕ) ∈ (0, ∞) such that for all ∆t ∈ (0, ∆t 0 ] and ∈ (0, 1] one has

Lemma 5.2. There exists C(T, ϕ) ∈ (0, ∞) such that for all ∈ (0, 1] one has

(5.2) E[ϕ(X(T ))] -E[ϕ(X (T ))] ≤ C(T, ϕ) .

Lemma 5.3. There exists C(T, ϕ) ∈ (0, ∞) such that for all ∆t ∈ (0, ∆t 0 ] one has

Lemma 5.4. There exists C(T, ϕ) ∈ (0, ∞) such that for all ∆t ∈ (0, ∆t 0 ] and ∈ (0, 1] one has

The first auxiliary result (Lemma 5.1) states a weak error estimate for the numerical scheme (3.4) for fixed ∈ (0, 1]. Due to the stiffness of the fast component m n , the right-hand side is not uniform with respect to , and it is natural to expect that the upper bound depends on τ = ∆t .

The second auxiliary result (Lemma 5.2) gives an error estimate in the averaging principle (see (2.11) in Proposition 2.4), in the weak sense. This is a standard result in the literature, see for instance [START_REF] Khasminskii | Limit behavior of two-time-scale diffusions revisited[END_REF] for an approach using asymptotic expansions for solutions of Kolmogorov equations. The strategy of the proof provided in Section A.1 is based on the introduction of the solutions of relevant Poisson equations, in the spirit of [START_REF] Pavliotis | Multiscale methods[END_REF]Chapter 17] where strong convergence is studied, see [START_REF] Bréhier | Orders of convergence in the averaging principle for SPDEs: the case of a stochastically forced slow component[END_REF] and [START_REF] Röckner | Strong and weak convergence in the averaging principle for sdes with hö older coefficients[END_REF] for the weak convergence case.

The two remaining auxiliary lemmas and their proofs are more original than the first two. Lemmas 5.3 and 5.4 are quantitative statements concerning two fundamental requirements in the notion of AP scheme (see Definition 3.4). On the one hand, Lemma 5.3 is a quantitative statement of the consistency of the limiting scheme (3.5) with the limiting equation (2.7), since it provides a weak error when ∆t → 0. Since the scheme is not classical (it is not a standard Euler-Maruyama type method, in particular recall that the scheme is random even if X(T ) is deterministic, when σ = 0), a proof is required. On the other hand, Lemma 5.4 is a quantitative statement about the convergence to the limiting scheme, for fixed ∆t ∈ (0, ∆t 0 ] (see Assumption 3.3). In fact, the left-hand side of (5.4) goes to 0 when → 0, however in the right-hand side of (5.4) an additional error term ∆t appears. Proving Lemma 5.4 is the most challenging step towards the proof of Theorem 3.8, whereas a key argument will be identified in the proof of Lemma 5.3 related to the consistency of the limiting scheme with the limiting equation.

The following auxiliary results concerning solutions of Kolmogorov equations are required in order do prove the four auxiliary results stated above. Lemma 5.5. Define u (t, x) = E x,m [ϕ(X (t))], for all t ≥ 0, x ∈ T d and m ∈ R, where X (t), m (t) t≥0 is the solution of the SDE system (2.1), and E x,m means that X (0) = x and m (0) = m. For all ∈ (0, 1], one has u ∈ C([0, T ], C 3 b (T d × R, R)). In addition, there exists C(T, ϕ) ∈ (0, ∞) such that for all j ∈ {1, 2, 3}, one has (5.5) sup

A second order Taylor expansion then gives

Using the regularity estimates from Lemma 5.7, one has R n (∆t)

1. Note that the terms of the orders ∆t 1 2 and ∆t 3 2 in the right-hand side above vanish, since the random variables Γ n and (X n , m n+1 , γ n ) are independent, and

Like in the proofs of Theorem 3.7 and of Lemma 5.3, a conditioning argument allows us to rewrite the expressions above in terms of the functions b and σ: one has

and

We are now in position to employ similar arguments as in the proof of Lemma 5.2 (see Section A.1), with important modifications due to the discrete-time setting. Introduce the auxiliary parameter τ = ∆t . Instead of studying Poisson equations associated with the infinitesimal generator L x OU , one needs to consider the generator L x τ and the transition semigroup of a Markov chain: let