Optimal explicit stabilized integrator of weak order one for stiff and ergodic stochastic differential equations

A new explicit stabilized scheme of weak order one for stiff and ergodic stochastic differential equations (SDEs) is introduced. In the absence of noise, the new method coincides with the classical deterministic stabilized scheme (or Chebyshev method) for diffusion dominated advection-diffusion problems and it inherits its optimal stability domain size that grows quadratically with the number of internal stages of the method. For mean-square stable stiff stochastic problems, the scheme has an optimal extended mean-square stability domain that grows at the same quadratic rate as the deterministic stability domain size in contrast to known existing methods for stiff SDEs [A. Abdulle and T. Li. Commun. Math. Sci., 6(4), 2008, A. Abdulle, G. Vilmart, and K. C. Zygalakis, SIAM J. Sci. Comput., 35(4), 2013]. Combined with postprocessing techniques, the new methods achieve a convergence rate of order two for sampling the invariant measure of a class of ergodic SDEs, achieving a stabilized version of the non-Markovian scheme introduced in [B. Leimkuhler, C. Matthews, and M. V. Tretyakov, Proc. R. Soc. A, 470, 2014].


Introduction
We consider Itô systems of stochastic differential equations of the form dX(t) = f (X(t))dt + m r=1 g r (X(t))dW r (t), where X(t) is a stochastic process with values in R d , f : R d → R d is the drift term, g r : R d → R d , r = 1, . . . , m are the diffusion terms, and W r (t), r = 1, . . . , m, are independent one-dimensional Weiner processes fulfilling the usual assumptions. We assume that the drift and diffusion functions are smooth enough and Lipschitz continuous to ensure the existence and uniqueness of a solution of (1) on a given time interval (0, T ). We consider autonomous problems to simplify the presentation, but we emphasise that the scheme can also be extended to non-autonomous SDEs. A one step numerical integrator for the approximation of (1) at time t = nh is a discrete dynamical system of the form X n+1 = Ψ(X n , h, ξ n ) (2) where h denotes the stepsize and ξ n are independent random vectors. Analogously to the deterministic case, standard explicit numerical schemes for stiff stochastic problems, such as the simplest Euler-Maruyama method defined as X n+1 = X n + hf (X n ) + m r=1 g r (X n )∆W n,r , where ∆W n,r = W r (t n+1 ) − W r (t n ) are the Brownian increments, face a severe timestep restriction [17,16,18], and one can use an implicit or semi-implicit scheme with favourable stability properties. In particular, it is shown in [17] that the implicit θ-method of weak order one is mean-square A-stable if and only if θ ≥ 1/2, while weak order two mean-square A-stable are constructed in [7]. An alternative approach is to consider explicit stabilized schemes with extended stability domains, as proposed in [4,5]. In [5] the deterministic Chebyshev method is extended to the context of mean-square stiff stochastic differential equations with Itô noise, while the Stratonovitch noise case is treated in [4]. In place of a standard small damping, the main idea in [4,5] is to use a large damping parameter η optimized for each number s of stages to stabilize the noise term. This yields a family of Runge-Kutta type schemes with extended stability domain with size L s 0.33s 2 . This stability domain size was improved to L s 0.42s 2 in [8] where a family of weak second order stabilized schemes (and strong order one under suitable assumptions) is constructed based on the deterministic orthogonal Runge-Kutta-Chebyshev method of order 2 (ROCK2) [6].
For ergodic SDEs, i.e., when (1) has a unique invariant measure µ satisfying for each test function φ and for any deterministic initial condition X 0 = x, one is interested in approximating numerically the long-time dynamics and to design numerical scheme with a unique invariant measure such that where C is independent of h small enough and X 0 . In such a situation, we say that the numerical scheme has order r with respect to the invariant measure. For instance, the Euler-Maruyama method has order 1 with respect to the invariant measure. In [20] the following non-Markovian scheme with the same cost as the Euler-Maruyama method was proposed for Brownian dynamics, i.e where the vector field is a gradient f (x) = −∇V (x) and the noise is additive (g(x) = σ), X n+1 = X n + hf (X n ) + σ ∆W n,j + ∆W n+1,j 2 , and it was shown in [21] that (6) has order 2 with respect to the invariant measure for Brownian dynamics. However, the admissible stepsizes for such an explicit method to be stable may face a severe restriction and alternatively to switching to drift-implicit methods, one may ask if a stabilized version of such an attractive non-Markovian scheme exists. In this paper we introduce a new family of explicit stabilized schemes with optimal meansquare stability domain of size L s = Cs 2 , where C ≥ 2 − 4 3 η and η ≥ 0 is a small parameter. We emphasize that in the deterministic case, L s = 2s 2 is the largest, i.e. optimal, stability domain along the negative real axis for an explicit s-stage Runge-Kutta method [16]. We note that the Chebyshev method (8) (with η = 0) realizes such an optimal stability domain. The new schemes have strong order 1/2 and weak order 1. The main ingredient for the design of the new schemes is to consider second kind Chebyshev polynomials, in addition to the usual first kind Chebyshev polynomials involved in the deterministic Chebychev method and stochastic extensions [5,4]. For stiff stochastic problems, the stability domain sizes are close to the optimal value 2s 2 and in the deterministic setting the method coincide with the optimal first order explicit stabilized method. Thus these methods are more efficient than previously introduced stochastic stabilized methods [5,8]. For ergodic dynamical systems, in the context of the ergodic Brownian dynamics, the new family of explicit stabilized schemes allows for a postprocessing [29] (see also [9,19] in the context of Runge-Kutta methods) to achieve order two of accuracy for sampling the invariant measure. In this context, our new methods can be seen as a stabilized version of the non-Markovian scheme (6) introduced in [20,21]. This paper is organized as follows. In Section 2, we introduce the new family of schemes with optimal stability domain and we recall the main tools for the study of stiff integrators in the mean-square sense. We then analyze its mean-square stability properties (Section 3), and convergence properties (Section 4). In Section 5, using a postprocessor we present a modification with negligible overcost that yields order two of accuracy for the invariant measure of a class of ergodic overdamped Langevin equation. Finally, Section 6 is dedicated to the numerical experiments that confirm our theoretical analysis and illustrate the efficiency of the new schemes.

New second kind Chebyshev methods
In this section we introduce our new stabilized stochastic method. We first briefly recall the concept of stabilized methods. In the context of ordinary differential equations (ODEs), and the Euler method X 1 = X 0 + hf (X(0)), a stabilization procedure based on recurrence formula has been introduced in [27]. Its construction relies on Chebyshev polynomials (hence the alternative name "Chebyshev methods"), T s (cos x) = cos(sx) and it is based on the explicit s-stage Runge-Kutta method where and for all i = 2, . . . , s, One can easily check that the (family of) methods (8) has the same first order accuracy as the Euler method (recovered for s = 1). In addition, the scheme (8) has a low memory requirement (only two stages should be stored when applying the recurrence formula) and it has a good internal stability with respect to round-off errors [27]. The attractive feature of such a scheme comes from its stability behavior. Indeed, the method (8) applied to the linear test problem dX(t)/dt = λX(t) yields, using the recurrence relation  where T 0 (p) = 1, T 1 (p) = p, with p = λh, where the dependence of the stability function R s,η on the parameters s and η is emphasized with a corresponding subscript. The real negative interval (−C s (η) · s 2 , 0) is included in the stability domain of the method S := {p ∈ C; |R s,η (p)| ≤ 1}.
The constant C s (η) 2 − 4/3 η depends on the so-called damping parameter η and for η = 0, it reaches the maximal value C s (0) = 2. Hence, given the stepsize h, for systems with a Jacobian having large real negative eigenvalues (such as diffusion problems) with spectral radius λ max at X n , the parameter s for the next step X n+1 can be chosen adaptively as 1 see [1] in the context of deterministically stabilized schemes of order two with adaptative stepsizes. The method (8) is much more efficient as its stability domain increases quadratically with the number s of function evaluations while a composition of s explicit Euler steps (same cost) has a stability domain that only increases linearly with s. In Figure 1(a) we plot the complex stability domain {p ∈ C ; |R s,η (p)| ≤ 1} for s = 7 stages and different values η = 0, η = 0.05 and η = 3.98, respectively. We also plot in Figure 1(b) the corresponding stability function R s,η (p) as a function of p real, to illustrate that the stability domain along the negative real axis corresponds to the values for which |R s,η (p)| ≤ 1. We observe that in the absence of damping (η = 0), the stability domain includes the large real interval [−2 · s 2 , 0] of width 2 · 7 2 = 98. However for all p that are a local extrema of the stability function, where |R s,η (p)| = 1, the stability domain is very thin and does not include a neighbourhood close to the negative real axis. To make the scheme robust with respect to small perturbations of the eigenvalues, it is therefore needed to add some damping and a typical value is η = 0.05, see for instance the reviews [28,2]. The advantage is that the stability domain now includes a neighbourhood of the negative real axis portion. The price of this improvement is a slight reduction of the stability domain size C η s 2 , where C η 2 − 4 3 η. Chebyshev methods have been first generalized for Itô SDEs in [5] (see [4] for Stratonovitch SDEs) with the following stochastic orthogonal Runge-Kutta-Chebyshev method (S-ROCK): 2 . . , s, where the coefficients µ i , ν i , κ i are defined in (9), (10). In contrast to the deterministic method (8), where η is chosen small and fixed (typically η = 0.05), in stochastic case for the classical S-ROCK method [5], the damping η = η s is not small and chosen as an increasing function of s that plays a crucial role in stabilizing the noise and in obtaining an increasing portion of the true stability domain (19) as s increases.
In the context of stiff SDEs, a relevant stability concept is that of mean-square stability. A test problem widely used in the literature is [25,17,11,26] , in dimensions d = m = 1 with fixed complex parameters λ, µ. Note that other stability test problem in multiple dimensions are also be considered in [10] and references therein. The exact solution of (16) is called mean-square stable if lim t→∞ E |X(t)| 2 = 0 and this holds if and only if (λ, µ) ∈ S M S , where Indeed, the exact solution of (16) is given by X(t) = exp((λ + 1 2 µ 2 )t + µW (t)), and an application of the the Itô formula yields E(|X(t)| 2 ) = exp(( (λ) + 1 2 µ 2 )t) which tends to zero at infinity if and only if (λ) + 1 2 µ 2 < 0. We say that a numerical scheme {X n } for the test problem (16) is mean-square stable if and only if lim n→∞ E(|X n | 2 ) = 0. For a one-step integrator applied to the test SDE (16), we obtain in general a induction of the form where p = λh, q = µ √ h, and ξ n is a random variable (e.g. a Gaussian ξ n ∼ N (0, 1) or a discrete random variable). Using E(|X n+1 | 2 ) = E(|R(p, q, ξ n )| 2 )E(|X n | 2 ), we obtain the mean-square stability condition [25,17] where we define S num = (p, q) ∈ C 2 ; E|R(p, q, ξ)| 2 < 1 . The function R(p, q, ξ n ) is called the stability function of the one-step integrator. For instance, the stability function of the Euler-Maruyama method (3) reads R(p, q, ξ) = 1 + p + qξ and we have E(|R(p, q, ξ)| 2 ) = (1 + p) 2 + q 2 .
We say that a numerical integrator is mean-square A-stable if S M S ⊂ S num . This means that the numerical scheme applied to (16)      for which the exact solution of (16) is mean-square stable. An explicit Runge-Kutta type scheme cannot however be mean-square stable because its stability domain S num is necessary bounded along the p-axis. Following [4,5], we consider the following portion of the true mean-square stability domain and define for a given method We search for explicit schemes for which the length L of the stability domain is large. For example, for the classical S-ROCK method [5], the value η = 3.98 is the optimal damping maximising L for s = 7 stages and we can see in Figure 1 that this damping reduces significantly the stability domain compared to the optimal deterministic domain.
The new S-ROCK method, denoted SK-ROCK (for stochastic second kind orthogonal Runge-Kutta-Chebyshev method) introduced in this paper is defined as . . , s.
where Q = m r=1 g r (X 0 )∆W r , and µ 1 = ω 1 /ω 0 , ν 1 = sω 1 /2, κ 1 = sω 1 /ω 0 and µ i , ν i , κ i , i = 2, . . . , s are given by (10), with a fixed small damping parameter η. In the absence of noise (g r = 0, r = 1, . . . , m, deterministic case), this method coincides with the standard deterministic order 1 Chebychev method, see the review [2]. We observe that the new class of methods (21) is closely related to the standard S-ROCK method (15). Comparing the two schemes (21) and (15), the two differences are on the one hand that the noise term is computed at the first internal stage K 1 for (21), whereas it is computed at the final stage in (15), and on the other hand, for the new method (21) the damping parameter η involved in (9) is fixed and small independently of s (typically η = 0.05), whereas for the standard method (15), the damping η is an increasing function of s, optimized numerically for each number of stages s.
If we apply the above scheme (21) to the linear test equation (16), we obtain and correspond to the drift and diffusion contributions, respectively. The above stability function (see Lemma 3.1 in Section 3) is obtained by using the recurrence relation for the first kind Chebyshev polynomials (11) and the similar recurrence relation for the second kind Chebyshev polynomials where U 0 (p) = 1, U 1 (p) = 2p. Notice that the relation T s (p) = sU s−1 (p) between first and second kind Chebyshev polynomials will be repeatedly used in our analysis. In Figure 2(b)(d), we plot the mean-square stability domain of the SK-ROCK method for s = 7 and s = 20 stages, respectively and the same small damping η = 0.05 as for the deterministic Chebyshev method. We observe that the stability domain has length L s (2 − 4 3 η)s 2 . For comparison, we also include in Figure 2(a)(c) the mean-square stability domain of the standard S-ROCK method with smaller stability domain size L s 0.33 · s 2 .
Noticing that E(|R(p, q, ξ)| 2 ) is an increasing function of q, the case q 2 = −2p represented in Figure 3(c) corresponds to the upper border of the stability domain S L defined in (19) (note that this is the stability function value along the dashed boundary in Figure 2), while the scaling q 2 = −p in Figure 3(c) is an intermediate regime. In Figures 3(b)(c), we also include the drift function A(p) 2 (red dotted lines) and diffusion function B(p) 2 q 2 (blue dashed lines), and it can be observed that their oscillations alternate, which means that any local maxima of one function is close to a zero of the other function. This is not surprising because A(p) and B(p) are related to the first kind and second kind Chebyshev polynomials, respectively, corresponding to the cosine and sine functions. This also explains how a large mean-square stability domain can be achieved by the new SK-ROCK method (21) with a small damping parameter η, in contrast to the standard S-ROCK method (15) from [5] that uses a large and s-dependent damping parameter η with smaller stability domain size L s 0.33 · s 2 (see Figures 3(a)(c)).

Mean-square stability analysis
In this section, we prove asymptotically that the new SK-ROCK methods have an extended mean-square stability domain with size Cs 2 growing quadratically as a function of the number of internal stages s, where the constant C ≥ 2 − 4 3 η is the same as the optimal constant of the standard Chebyshev method in the deterministic case, using a fixed and small damping parameter η.
is a Gaussian variable and the stability function given by Proof. Indeed, we take advantage that T j and U j have the same recurrence relations (11), (23), and only the initialization changes with T 1 (x) = x and U 1 (x) = 2x, we deduce Q = X 0 µ √ hξ, and we obtain by induction on i ≥ 1, For a positive damping η, we prove the following main result of this section, showing that a quadratic growth L ≥ (2−4/3 η)s 2 of the mean-square stability domain defined in (20) is achieved for all η small enough and all stage number s large enough.
Before we prove Theorem 3.2, we have the following lemma, see [27] for analogous results.
Lemma 3.4. We have the following convergences as s → ∞ to analytic functions 3 uniformly for z in any bounded set of the complex plan, 3 Note that for z < 0, we can use √ 2z = i √ −2z and obtain Ts(1 + z/s 2 ) → cos( √ −2z) for s → ∞ and similarly Proof. We prove the uniform convergence of the first limit only, since it will be useful in the proof of the next theorem. The others can be proved in a similar way. First, let us write the two functions of η in Taylor series, where we used the formula sU Noticing that 2s−1 converges to zero as s → ∞ and is bounded by n! for all integers s, n, which is the general term of the convergent series of exp(η 0 ) = ∞ n=0 η n 0 n! , the Lebesgue dominated convergence theorem implies that (25) converges to zero as s → ∞, which concludes the proof.
Proof of Lemma 3.5. Using the Lemma 3.4 we have for s → ∞, uniformly for all η ∈ [0, η 0 ], Now if we expand both functions in Taylor series we get: and this implies that for all s large enough and all η ≤ η 0 , which is less than 1 for η 0 small enough.
Remark 3.6. Numerical evidence suggests that the result of Theorem 3.2 holds for all s ≥ 1 and all η ≥ 0. Indeed, it can be checked numerically that (26) holds for all η ∈ (0, 1) and all s ≥ 1.
Proof of Theorem 3.2. Setting x = w 0 + w 1 p, a calculation yields The proof is conducted in two steps, where we treat separately the cases p ∈ [−2ω −1 where we denote is an increasing function of x, smaller than its value at Using Lemma 3.5 we obtain |Q s (x)| ≤ 1. This yields E(|R(p, q, ξ) Using Lemma 3.4, we get for s → ∞, where the above convergence is uniform for p ∈ [0, 1],η ≤ η 0 . Using the fact that By Taylor series in the neighbourhood of zero we have α(p) 2 = 1 +

Convergence analysis
We show in this section that the proposed scheme (21) has strong order 1/2 and weak order 1 for general systems of SDEs of the form (1) with Lipschitz and smooth vector fields, analogously to the simplest Euler-Maruyama method.
We denote by C 4 P (R d , R d ) the set of functions from R d to R d that are 4 times continuously differentiable with all derivatives with at most polynomial growth. The following theorem shows that the proposed SK-ROCK has strong order 1/2 and weak order 1 for general SDEs.
For the proof the Theorem 4.1, the following lemma will be useful. It relies on the linear stability analysis of Lemma 3.1.
Lemma 4.2. The scheme (21) has the following Taylor expansion after one timestep, where all the moments of R h (X 0 ) are bounded uniformly with respect to h assumed small enough, with a polynomial growth with respect to X 0 .
Proof. Using the definition (21) of the scheme and the recurrence relations (11),(23), we obtain by induction on i = 1, . . . , s, and R i,h (X 0 ) has the properties claimed on R h (X 0 ). Using ω 1 = T s (ω 0 )/T s (ω 0 ), this yields the result for X 1 = K s .
for all initial value X(0) = X 0 and where C has at most polynomial growth with respect to X 0 , to deduce the weak convergence estimate (30). For the strong convergence (30), using the classical result from [22], it is sufficient to show in addition the local error estimate for all initial value X(0) = X 0 and where C has at most polynomial growth with respect to X 0 . These later two local estimates are an immediate consequence of Lemma 4.2.
To conclude the proof of the global error estimates, it remains to check that for all r ∈ N the moments E(|X n | 2r ) are bounded uniformly with respect to all h small enough for all 0 ≤ nh ≤ T . We use here the approach of [24, Lemma 2.2, p. 102] which states that it is sufficient to show where C is independent of h and M n is a random variable with moments of all orders bounded uniformly with respect to all h small enough. These estimates are a straightforward consequence of the definition (21)   , 0] and all q with p + |q| 2 2 ≤ 0, as it can be check numerically. The idea is to modify the first stages of the scheme such that the stability function (24) becomes We refer to [8, Remark 3.2] for details.

Long term accuracy for Brownian dynamics
In this section we discuss the long-time accuracy of the SK-ROCK for Brownian dynamics (also called overdamped Langevin dynamics). We will see that using postprocessing techniques we can derive an SK-ROCK method that captures the invariant measure of Brownian dynamics with second order accuracy. In doing so, we do not need our stabilized method to be of weak order 2 on bounded time intervals and we obtain a method that is cheaper than the stochastic orthogonal Runge-Kutta-Chebyshev method of weak order 2 (S-ROCK2) proposed in [8], as S-ROCK2 uses many more function evaluations per time-step and a smaller stability domain.

An exact SK-ROCK method for the Orstein-Uhlenbeck process
We consider the 1-dimensional Orstein-Uhlenbeck problem with 1-dimensional noise with constants δ, σ > 0, that is ergodic and has a Gaussian invariant measure with mean zero and variance given by lim t→∞ E(X(t) 2 ) = σ 2 /(2δ). Applying the SK-ROCK method to the above system we obtain where p = −δh, ξ n ∼ N (0, 1) is a Gaussian variable and similarly as for (24) we have A simple calculation (using that |A(p)| < 1) gives From the above equation, we see that the SK-ROCK method has order r for the invariant measure of (33) if and only if R(p) = 1 + O(p r ) and a short calculation using (35) reveals that R(p) = 1 + O(p), it has order one for the invariant measure (this is of course not surprising because the SK-ROCK has weak order one). We next apply the techniques of postprocessed integrators popular in the deterministic literature [12] and proposed in the stochastic context in [29]. The idea is to consider a postprocessed dynamics X n = G n (X n ) (of negligible cost) such that the process X n approximates the invariant measure of the dynamical system with higher order. For the process (33), we consider the postprocessor which yields lim n→∞ E(X 2 n ) = σ 2 2δ (R(p) − 2c 2 p). In the case of the SK-ROCK method with η = 0 (zero damping), we have A(p) = T s (1 + p/s 2 ), B(p) = U s−1 (1 + p/s 2 )(1 + p/(2s 2 ))/s. Setting c = 1/(2s) and using the identity ( Hence the postprocessed SK-ROCK method (that will be denoted PSK-ROCK) captures exactly the invariant measure of the 1-dimensional Orstein-Uhlenbeck problem (33). Such a behavior is known for the drift-implicit θ method with θ = 1/2 (see [13] in the context of the stochastic heat equation) and has recently also been shown for the non-Markovian Euler scheme [20]. In [29] an interpretation of the scheme [20] as an Euler-Maruyama method with postprocessing (36) with c = 1/2 has been proposed and we observe that this is exactly the same postprocessor as for the PSK-ROCK method (with s = 1, η = 0). As the SK-ROCK method with zero damping is mean-square stable (see Remark 3.3 for η = 0), it can be seen as a stabilized version of the scheme [20]. However, the PSK-ROCK method with s > 1 and zero damping is not robust to use as its stability domain along the drift axis does not allow for any imaginary perturbation at the points where |T s (1 + p/s 2 )| = 1 and it is not ergodic (see Remark 5.2 below).
Stability analysis for Orstein-Uhlenbeck Let M ∈ R d×d denote a symmetric matrix with eigenvalues −λ d ≤ . . . ≤ −λ 1 < 0, and consider the d-dimensional Orstein-Uhlenbeck problem where W (t) denotes a d-dimensional standard Wiener process. The following theorem shows that the damping parameter η > 0 plays an essential role to warranty the convergence to the numerical invariant measure ρ h ∞ (x)dx at an exponentially fast rate.
Proof. It is sufficient to show the estimate for all h ≤ h 0 , where we denote A(z) = T s (ω 0 + ω 1 z)/T s (ω 0 ). Indeed, considering two initial conditions X 1 0 , X 2 0 for (21) and the corresponding numerical solutions X 1 n , X 2 n (obtained for the same realizations of {ξ n }) with postprocessors X 1 n , X 2 n , we obtain X 1 n − X 2 n = A(hM )(X 1 n−1 − X 2 n−1 ) and using the matrix 2-norm A(hM ) = max j |A(−λ j h)| and (39), we deduce by induction on n, X and taking X Using ω 1 s 2 ≥ 1 and T s (ω 0 ) ≥ 1 + η, we obtain
Remark 5.2. Note that η > 0 is a crucial assumption in Theorem 5.1. Indeed, the estimate of Theorem 5.1 is false for η = 0 already in dimension d = 1 for all s > 1: for a stepsize h such that 1 − hλ 1 /s 2 = cos(π/s) we obtain A(−λ 1 h) = −1 and B(−λ 1 h) = 0 in (35) (corresponding to the local extrema p = −λ 1 h of A(p) closest to zero) and X n = (−1) n X 0 for all n, and the scheme is not ergodic. In addition, notice that Theorem 5.1 allows to use an h-dependent value of η such as η = h λ 1 where λ 1 ≥ λ 1 is an upper bound for λ 1 . In this case, the exponential convergence of Theorem 5.1 holds for all stepsize h ≤ 1.
We end this section by noting that being exact for the invariant measure of Brownian dynamics (40) is only true for the PSK-ROCK method (or the method in [20]) in the linear case, i.e. for a quadratic potential V . Second order accuracy for the invariant measure has been shown in [21] (see also [29]) for the method in (6) (equivalent to PSK-ROCK with s = 1 stage) for general nonlinear Brownian dynamics (40). This will also be shown for the nonlinear PSK-ROCK method in the next section.

PSK-ROCK: a second order postprocessed SK-ROCK method for nonlinear Brownian dynamics
We consider the overdamped Langevin equation, where the stochastic process X(t) takes values in R d and W (t) is a d-dimensional Wiener process. We assume that the potential V : R d → R has class C ∞ and satisfies the at least quadratic growth assumption for two constants C 1 , C 2 > 0 independent of x ∈ R d . The above assumptions warranty that the system (40) is ergodic with exponential convergence to a unique invariant measure with Gibbs for test function φ and all initial condition X 0 , where C, λ are independent of t. We propose to modify the internal stage K 1 = X 0 + µ 1 hf (X 0 + ν 1 Q) + κ 1 Q of the method (21) as follows: where α is a parameter depending on s and η given in Theorem 5.4 below. where and r s is defined by induction as r 0 = 0, r 1 = s 2 ω 3 1 4ω 0 := ∆ 1 and Then, X n yields order two for the invariant measure, i.e. (5) holds with r = 2, and in addition for all t n = nh, φ ∈ C ∞ P (R d , R), where C 1 , C 2 are independent of h assumed small enough, and C 2 is independent of the initial condition X 0 .
The proof of Theorem 5.4 relies on the following postprocessing analysis from [29]. Consider a scheme (2) with bounded moments and assumed ergodic when applied to (40), where the potential V satisfies the above ergodicity assumptions. Assume that the scheme has a weak Taylor expansion after one time step of the form and consider a postprocessor of the form X n = G n (X n ) where where the constants in O in (45),(46) have at most a polynomial growth with respect to x. Here Lφ = φ f + σ 2 /2∆φ denotes generator of the SDE and A 1 , A 1 are linear differential operators with smooth coefficients. Note that A 1 = L 2 /2 in general (otherwise the scheme has weak order 2). If the condition (A 1 + [L, A 1 ]) * ρ ∞ = 0 holds, equivalently, for all test function φ, where we define φ = R d φρ ∞ dx, then it is shown in [29, Theorem 4.1] that X n has order two for the invariant measure, i.e. the convergence estimates (5) with r = 2 and (44) hold. Before we can apply the above result, the following lemma allow to compute the weak Taylor expansion of the modified scheme.
Lemma 5.5. Consider the scheme (21) with modified stage (42) and assume the hypotheses of Theorem 5.4. Then (45) holds where the linear differential operator A 1 is given by where f = −∇V (x) and Proof. Adapting the proof of Lemma 4.2, the internal stage K i defined in (21) (and (42) for i = 1) satisfies (31) where h 2 R i,h (X 0 ) can be replaced by where E( R i ) = 0 and all the moments of R i , R i,h (X 0 ) are bounded with polynomial growth with respect to X 0 . Here,r i is defined by induction asr 0 = 0,r 1 = ∆ 1 + α, and We have E( R i ) = 0 because R i is a linear combination of f (X 0 )f (X 0 )ξ n , f (X 0 )(f (X 0 ), ξ n ), and f (X 0 )(ξ n , ξ n , ξ n ) with zero mean values (recall that odd moments of ξ n vanish). Next, observing that the difference d i =r i − r i satisfies d 0 = 0, d 1 = α, and d i = ν i d i−1 + κ i d i−2 , i = 2, . . . , s, we deducer In particular, taking i = s in (31),(50), and expanding (45), we deduce that (48) holds with c 2 , c 3 , c 4 defined in (49) where we note that c 3 = r s = r s + d s .   where we apply repeatedly integration by parts for the integral in (47), using Lemma 5.5 for the expression of A 1 , we deduce that the quantity in (47) satisfies where we use [L, For the values of α, c defined in (43), we obtain that (51) indeed holds and we deduce that the order two condition (47) for the invariant measure is satisfied. This concludes the proof.

Numerical experiments
In this Section, we illustrate numerically our theoretical analysis and we show the performance of the proposed SK-ROCK method and its postprocessed modification PSK-ROCK.

A nonlinear nonstiff problem
We first consider the following non-stiff nonlinear SDE, whose exact solution is X(t) = sinh( t 2 + W (t) √ 2 ). In Figure 4, we consider the SK-ROCK method (21) and plot the strong error E(|X(T ) − X N | and the weak error |E(arcsinh(X(T )) − E(arcsinh(X N ))| at the final time T = N h = 1 using 10 4 samples and number of stages s = 1, 5, 10, 100. We obtain convergence slopes 1 and 1/2, respectively, which confirms Theorem 4.1 stating the strong   order 1/2 and weak order 1 of the proposed scheme. Note that s = 1 stage is sufficient for the stability of the scheme in the non-stiff case. The results for s = 5, 10, 100 yield nearly identical curves which illustrates that the error constants of the method are nearly independent of the stage number of the scheme.

Nonlinear nonglobally Lipschitz stiff problems
Consider the following nonlinear SDE in dimensions d = 2 with a one-dimensional noise (d = 2, m = 1). This is a modification of a one-dimensional population dynamics model from [15,Chap. 6.2] considered in [7,8,3] for testing stiff integrator performances, Observe that linearizing (53) close to the equilibrium (X, Y ) = (1, 1), we recover for ν = 0 the scalar test problem (16). In Figure 5 we consider the SK-ROCK method applied to (53) with parameters that are identical to those used in [3,Sect. 4.2]. We take the initial condition X(0) = Y (0) = 0.95 close to this steady state and use the parameters ν = 2, µ 2 = 0.5,λ 2 = −1. In a nonstiff regime (−λ 1 = µ 1 = 1 in Figure 5(a)), we observe a convergence slope 1 for the second moment E(X(T ) 2 ) which illustrates the weak order one of the scheme, although our analysis in Theorem 4.1 applies only for globally Lipschitz vector fields. The stage number s = 1 is sufficient for stability, but we also include for comparison the results for s = 10 (note that the results for s = 50, 100 not displayed here are nearly identical to the case s = 10). The convergence curves are obtained as an average over 10 6 samples. In a stiff regime (−λ 1 = µ 2 1 = 100 in Figure 5(b)), we observe for the standard small damping η = 0.05 a stable but not very accurate convergence, due to the severe nonlinear stiffness. However, considering a slightly larger damping η = 4, in the spirit of the S-ROCK method, yields a stable integration for all considered timesteps and all trajectories and we observe a line with slope one for the SK-ROCK method. Here, given the timesteps h, the numbers of stages s are adjusted as proposed in (54) where λ max = |λ 1 | = 100.   Remark 6.1. For severely stiff problems, alternatively to switching to drift-implicit schemes [17,7], one can consider in SK-ROCK a slightly larger damping η and the corresponding stage parameter s below, similar to (14) and chosen such that the mean-square stability domain length (20) satisfies L > hλ max , where Ω(η) is given in Lemma 3.4. In Figure 6, we consider the SK-ROCK and PSK-ROCK methods with s = 1, 5, 10, 100 stages, respectively. For a short time T = 0.5 ( Fig. 6(a)(b)), we observe weak convergence slopes one for both SK-ROCK and PSK-ROCK (second moment E(X(T ) 2 )) as predicted by Theorem 4.1, and the postprocessor has nearly no effect of the errors. For a long time T = 10 where the solution of this ergodic SDE is close to equilibrium, we observe that the weak order one of SK-ROCK (Fig. 6(c)) is improved to order two using the postprocessor in PSK-ROCK (Fig. 6(d)), which confirms the statement of Theorem 5.4 that the postprocessed scheme has order two of accuracy for the invariant measure. For comparison, in Figure 7, we also include the results of PSK-ROCK without damping (η = 0) using M = 10 8 samples. We recall that for the scalar linear Orstein-Uhlenbeck process, the PSK-ROCK method with zero damping is exact for the invariant measure (see Section 5.1). We observe only Monte-Carlo errors with size M −1/2 = 10 −4 , which confirms that the PSK-ROCK method has no bias at equilibrium for the invariant measure in the absence of damping, as shown in (37). We emphasise however that this exactness results holds only for linear problems, and a positive damping parameter η should be used for nonlinear SDEs for stabilization, as shown in Sections 3 and 5.1.

Nonglobally Lipschitz Brownian dynamics
To illustrate the advantage of the PSK-ROCK method applied to nonglobally Lipschitz ergodic Brownian dynamics, we next consider the following double well potential V (x) = (1 − x 2 ) 2 /4 and the corresponding one-dimensional Brownian dynamics problem In Figure 8, we compare the performances of S-ROCK, S-ROCK2 considered in [8] (a method with weak order 2 for general SDEs), and the new SK-ROCK and PSK-ROCK methods at short time T = 0.5 (Figures 8(a)(b)) and long time T = 10 (Figures 8(c)(d)). As we focus on invariant measure convergence and not on strong convergence, we consider here discrete random increments with P(ξ n = ± √ 3) = 1/6, P(ξ n = 0) = 2/3, which has the correct moments so that Theorem 5.4 remains valid. Our numerical tests indicate that it makes PSK-ROCK with modified stage (42) more stable. For a fair comparison, we use the same discrete random increments for all schemes. We plot the second moment error versus the time stepsize h and versus the average cost which is the total number of function evaluations during the time integration divided by the total number number of samples. Indeed, the number of function evaluations depends on the    trajectories because the stage parameter s is adaptive at each time step. For short time, we can see that the S-ROCK and the SK-ROCK method have order 1 (Figure 8(a)) and exhibit similar performance with nearly identical error versus cost curves in Figure 8(b), while PSK-ROCK is less advantageous for short time. This illustrates that the postprocessing has no advantage for short times. The S-ROCK2 method is the most accurate for small time steps, and it has order 2 as shown in Figures 8(a)(c), but at the same time it has a larger average cost as observed in Figures 8(b)(d) due to its smaller stability domain with size 0.42 · s 2 . For long time, the SK-ROCK and S-ROCK both exhibit order 1 of accuracy (Figure 8(c)), with an advantage in terms of error versus cost for the SK-ROCK method that is about 10 times more accurate for large time steps. In contrast, the postprocessed scheme PSK-ROCK exhibits order 2 of convergence ( Figure 8(c)) which corroborates Theorem 5.4. Since the postprocessing overcost is negligible (two additional vector field evaluations per timestep due to the modified stage K 1 in (42)), this makes PSK-ROCK the most efficient in terms of error versus cost, as shown in Figure 8(d). The S-ROCK2 method has order 2 here but with poor accuracy compared to the PSK-ROCK method with approximately the same cost. Note that typically the SK-ROCK method used s = 1, 2, 3 stages in contrast to the S-ROCK method using s = 2, . . . , 6 stages per timesteps.

Stochastic heat equation with multiplicative space-time noise
Although our analysis applies only to finite dimensional systems of SDEs, we consider the following stochastic partial differential equation (SPDE) obtained by adding multiplicative noise to the heat equation, whereẆ (t, x) denotes a space-time white noise that we discretize together with the Laplace operator with a standard finite difference formula [14]. We obtain the following stiff system of SDEs where u(x i , t) ≈ u i (t), with x i = i∆x, ∆x = 1/N , where the Dirichlet and the Neumann conditions impose u 0 = 5 and u N +1 = u N −1 , respectively. Here, w 1 , . . . , w N are independent standard Wiener processes and dw i indicates Itô noise. In Figure 9(a), we plot one realization of the SPDE using space stepsize ∆x = 1/100 and timestep size ∆t = 1/50. Note that the Lipchitz constant associated to the space-discretization of (57) has size ρ = 4∆x −2 , and the stability condition is fulfilled for s = 22 stages. For comparison, the standard S-ROCK method would require s = 46 stages, while applying the standard Euler-Maruyama with a smaller stable timestep ∆t/s would require s ≥ ∆tρ/2 = 400 intermediate steps.
Notice that the initial condition in (57) satisfies the boundary conditions, which permits a smooth solution close to time t = 0. Taking alternatively an initial condition that does not satisfy the boundary conditions (for instance u(x, 0) = 1) yields an inaccurate numerical solution with large oscillations close to the boundary x = 0. A simple remedy in such a case is to consider a larger damping parameter η, as described in Remark 6.1.
In Figure 9(b), we compare the number of vector field evaluations of the standard S-ROCK and new SK-ROCK methods when applied to the SPDE (57) with finite difference discretization   with parameter ∆x = 1/100. The better performance of SK-ROCK with damping η = 0.05 is due to its larger stability domain with size 1.94 · s 2 compared to the size 0.33 · s 2 for S-ROCK. Observing the ratio of the two costs in Figure 9(b), we see that the new SK-ROCK methods has a reduced cost for stabilization by an asymptotic factor of about 1.94/0.33 2.4 for large s and large stepsizes, which confirms the stability analysis of Section 3. The convergence analysis of the SK-ROCK method for the stochastic heat equation is the topic of future work. Remark 6.2. Notice that SK-ROCK with s = 1 stage has the optimal mean-square stability length (L = 2 for η = 0) as defined in (20). In contrast, the S-ROCK method with s = 1 has the smaller stability length L = 3/2, while the standard Euler-Maruyama has L = 0. This explains why for the smallest considered stepsize ∆t = 2 −15 in Figure 9(b), we have s = 1 for SK-ROCK while S-ROCK uses s = 2 stages.
In Figure 10 we consider again one realization with SK-ROCK of the SPDE problem (57) but with a different initial condition u(0, x) = 1 not fulfilling the boundary conditions, i.e. that is outside the domain of the Laplace operator, as considered in [3]. We compare the result for the same sets of random numbers but for for different values of the damping parameter η. We observe numerically high oscillations in time and space for the small damping value η = 0.05 in Figure (10a) while the larger damping η = 10 yields a smoother solution in Figure (10b). This illustrates again Remark 6.1 showing that the damping parameter η can be increased in the case of severely stiff problems, adjusting the stage parameter accordingly with (54).