Assyr Abdulle

Charles-Edouard Bréhier

Gilles Vilmart

Convergence analysis

Keywords: explicit stabilized methods, second kind Chebyshev polynomials, stochastic partial differential equations, finite element methods AMS subject classification (2020): 65C30, 60H35, 65L20

come

Introduction

In this paper, we consider semilinear parabolic stochastic partial differential equations (SPDEs), in the framework of [START_REF] Da Prato | Stochastic equations in infinite dimensions[END_REF], of the form

du(t) = Λu(t)dt + F (u(t))dt + σdW Q (t), u(0) = u 0 , (1)
where u 0 is a given initial condition, Λ is a symmetric diffusion operator, and F is a smooth nonlinearity, σ > 0, W Q (t) t≥0 is a Q-Wiener process. The semi-discrete approximation obtained by spatial discretization using the finite element method with piecewise linear elements is defined by the following equation: for a spatial mesh size h ∈ (0, 1],

du h (t) = Λ h u h (t)dt + P h F (u h (t))dt + σP h dW Q (t), u h (0) = P h u 0 , (2)
where P h denotes the L2 -orthogonal projection operator onto the finite element space. The stochastic evolution problem ( 2) is driven by a Q h -Wiener process, with Q h = P h QP h . The spatial discretization Λ h of the diffusion operator Λ typically yields large eigenvalues, with spectral radius of size O(h -2 ) for a Laplace operator and a finite element mesh with size h, which yields a severe timestep restriction τ = O(h 2 ) 1 for standard explicit integrators. To avoid such prohibitive time-step size restrictions, a natural approach is to consider implicit or semi-implicit schemes, such as the linear implicit Euler method, for n ≥ 0,

u h n+1 = (I -τ Λ h ) -1 u h n + P h τ F (u h n ) + σ∆W Q n , (3)
where h ∈ (0, 1] and τ ∈ (0, 1), and the Wiener increments are defined by ∆W Q n = W Q (t n+1 ) -W Q (t n ) with t n = nτ . Alternatively to implicit methods, in this paper we consider families of explicit stabilized schemes of the form

u h n+1 = A s (τ Λ h )u h n + B s (τ Λ h )P h τ F (u h n ) + σ∆W Q n , (4)
where A s and B s are polynomials of degree s = s(τ, h) which are chosen to satisfy a suitable stability condition depending on τ and h. For well-chosen polynomials, the stability domain can be very large, which allow to choose the time-step size τ independently of the mesh size h. In particular, we shall consider variants of the SK-ROCK method from [START_REF] Abdulle | Optimal explicit stabilized integrator of weak order 1 for stiff and ergodic stochastic differential equations[END_REF]. Explicit stabilized integrators are well known efficient time integrators for stiff problems, in particular arising from diffusion PDE problems, both in the deterministic and stochastic settings, see [START_REF] Hairer | Solving ordinary differential equations II. Stiff and differential-algebraic problems[END_REF]Sect. IV.2] and the review [START_REF] Abdulle | Explicit Stabilized Runge-Kutta Methods[END_REF]. Compared to the explicit stabilized integrators first proposed in the stochastic context in [START_REF] Abdulle | Chebyshev methods for stiff stochastic differential equations[END_REF][START_REF] Abdulle | S-ROCK methods for stiff Ito SDEs[END_REF] using only first kind Chebyshev polynomials and a large damping parameter η, a new optimal family of explicit stabilized method involving second kind Chebyshev polynomials, analogous to [START_REF] Abdulle | S-ROCK methods for stiff Ito SDEs[END_REF], was introduced in [START_REF] Abdulle | Optimal explicit stabilized integrator of weak order 1 for stiff and ergodic stochastic differential equations[END_REF]. It permits the efficient integration of stiff systems of stochastic differential equations, with favourable mean-square stability properties and high-order for sampling the invariant measure of ergodic problems such as the overdamped Langevin equation.

While the convergence analysis of stochastic explicit stabilized methods is known in finite dimension, see e.g. [START_REF] Abdulle | Chebyshev methods for stiff stochastic differential equations[END_REF][START_REF] Abdulle | S-ROCK methods for stiff Ito SDEs[END_REF][START_REF] Abdulle | PIROCK: a swiss-knife partitioned implicit-explicit orthogonal Runge-Kutta Chebyshev integrator for stiff diffusion-advection-reaction problems with or without noise[END_REF][START_REF] Abdulle | Weak second order explicit stabilized methods for stiff stochastic differential equations[END_REF], applied to a semi-discrete in space problem [START_REF] Abdulle | Optimal explicit stabilized integrator of weak order 1 for stiff and ergodic stochastic differential equations[END_REF], in general one obtains convergence estimates with error constants that depend on the spatial mesh size h. The aim of this paper is to remove this dimension dependency and to provide the first strong convergence analysis of families of explicit stabilized methods in the context of semilinear parabolic SPDEs such as the stochastic heat equation, with error constants that are independent of the space and time mesh parameters h and τ . Under natural assumptions, for a fixed time interval size T > 0, we prove that the families of explicit stabilized schemes (4) satisfy the following space-time strong error estimate. For all α ∈ [0, 1] chosen small enough such that (-Λ) (α-1)/2 Q 1/2 is a bounded Hilbert-Schmidt operator (corresponding to α ∈ [0, 1/2) for the one-dimensional stochastic heat equation with space time-white noise), we show that the mean-square strong error satisfies

E|u(t n ) -u h n | 2 1 2 ≤ C(1 + |u 0 |t -α 2 n )(τ α 2 + h α ),
for all τ, h ∈ (0, 1] with t n = nτ ≤ T , where C is independent of n and both the space and time mesh parameters h, τ and the degree s = s(τ, h). The constant C depends however on T and α. The parameter α is related to the spatial and temporal regularity of the process. Note that in the deterministic case (σ = 0 in (1)), we obtain an order one of convergence in time and order two in space, corresponding formally to α → 2. Here, | • | denotes the L 2 (D) norm in space on a bounded, open and convex polyhedral domain D in dimension d.

While the scheme (4) achieves the same strong convergence rate as the classical linear implicit Euler method (3), as first studied in [START_REF] Printems | On the discretization in time of parabolic stochastic partial differential equations[END_REF], we emphasize that the analysis in this paper is not a straightforward generalization of that in [START_REF] Printems | On the discretization in time of parabolic stochastic partial differential equations[END_REF] due to lower regularization properties of the explicit numerical flow for the diffusion A s (τ Λ h ) in (4) compared to the linear implicit Euler method [START_REF] Abdulle | Chebyshev methods for stiff stochastic differential equations[END_REF]. Furthermore, the convergence analysis of an explicit method for parabolic SPDEs with space and time mesh independent constants is a new result, to the best of our knowledge. In addition, note that the convergence analysis is performed for a general class of explicit stabilized methods satisfying suitable regularization conditions, and it provides a unified natural framework for Runge-Kutta type schemes, including the linear implicit Euler method (3), applied to semilinear parabolic SPDEs.

This paper is organized as follows. In Section 2 we describe the classical Hilbert space setting and assumptions used for the analysis of SPDEs and of their numerical approximations. Section 3 is devoted to the definition of the explicit stabilized method in time, coupled with a standard finite element method in space and presents our main convergence results. Section 4 is dedicated to the strong convergence analysis of the method. Finally Section 5 is dedicated to numerical experiments which confirm the theoretical findings.

Setting

Let H be the separable Hilbert space L The assumptions stated in this section are standard in the literature on SPDEs [START_REF] Da Prato | Stochastic equations in infinite dimensions[END_REF] and their numerical approximations [START_REF] Kruse | Optimal error estimates of Galerkin finite element methods for stochastic partial differential equations with multiplicative noise[END_REF], [START_REF] Lord | An introduction to computational stochastic PDEs[END_REF].

Linear operator Λ

The linear operator Λ is defined as the linear second-order elliptic differential operator Λu(x) = div a(x)∇u(x) ,

with domain D(Λ) = H 2 (D) ∩ H 1 0 (D)
, where we consider for simplicity homogeneous Dirichlet boundary conditions on the boundary of the domain (u = 0 on ∂D). We assume that the field a : D → R is continuous on the closure D of D, of class C ∞ on D, and satisfies min x∈D a(x) > 0. The following result is standard.

Proposition 2.1. The linear operator Λ is unbounded and self-adjoint. There exists a complete orthonormal system e m m≥1 of H, and a non-decreasing sequence λ m m≥1 , such that one has

λ 1 > 0, lim m→∞ m -2 d λ m ∈ (0, ∞)
, and for all m ≥ 1,

Λe m = -λ m e m .
The linear operator Λ generates a strongly continuous semi-group on H, denoted by e tΛ t≥0 . For all u ∈ H and all t ≥ 0, one has

e tΛ u = m≥1 e -λmt u, e m e m .
In addition, for any α ∈ [-1, 1], the linear operator (-Λ) α is defined as follows: for all u ∈ H, one has

(-Λ) α u = m≥1 λ α m u, e m e m .
The semi-group e tΛ t≥0 satisfies the following smoothing and regularity properties: for any γ ∈ [0, 1],

sup t≥0 t γ (-Λ) γ e tΛ L(H) < ∞ , sup t>0 t -γ e tΛ -I)(-Λ) -γ L(H) < ∞. (5)

Nonlinear operator F

The nonlinear operator F in the stochastic evolution equation ( 1) is assumed to be globally Lipschitz continuous.

Assumption 2.2. The nonlinear operator F is a globally Lipschitz continuous mapping from H to H. Precisely there exists a constant C such that for all u, v ∈ H,

|F (u) -F (v)| ≤ C|u -v|.
For instance, let f : R → R be a globally Lipschitz continuous real-valued function. If F is defined as the associated so-called Nemytskii operator, precisely F (u) = f (u(•)) for all u ∈ H, then Assumption 2.2 is satisfied. The local Lipschitz continuous case is out of the scope of this article, and would require to modify the scheme (using for instance implicit methods, splitting [START_REF] Bréhier | Analysis of some splitting schemes for the stochastic Allen-Cahn equation[END_REF] or tamed or truncated explicit [START_REF] Hutzenthaler | Numerical approximations of stochastic differential equations with non-globally Lipschitz continuous coefficients[END_REF] schemes, adaptive timestepping [START_REF] Kelly | Adaptive time-stepping strategies for nonlinear stochastic systems[END_REF] techniques).

Q-Wiener process and mild solution of the SPDE (1)

Let us first introduce the Q-Wiener process W Q (t) t≥0 . Let m m≥1 be a complete orthonormal system of the Hilbert space H, and consider a bounded sequence q m m≥1 of nonnegative real numbers. Let Q and Q 1 2 be the bounded, linear, self-adjoint operators on H, defined by

Qu = m≥1 q m u, e m e m , Q 1 2 u = m≥1 √ q m u, e m e m .
Definition 2.3. Let β m m≥1 be a family of independent standard real-valued Wiener processes, β m = β m (t) t∈R + , defined on a probability space which satisfies the usual conditions. The Q-Wiener process W Q (t) t∈R + is defined as follows: for all t ≥ 0,

W Q (t) = m≥1 √ q m β m (t) m .
Note that the Q-Wiener process W Q takes values in H if and only if m≥1 q m < ∞ (which means that Q is a trace-class operator, or equivalently that Q 1 2 is an Hilbert-Schmidt operator). Using the regularization properties [START_REF] Abdulle | Coupling heterogeneous multiscale FEM with Runge-Kutta methods for parabolic homogenization problems: a fully discrete space-time analysis[END_REF] for the semi-group e tΛ t≥0 generated by Λ, it is possible to define solutions of stochastic evolution equations driven by a Q-Wiener process without requiring that Q is trace-class. More precisely, the wellposedness of the stochastic evolution equation ( 1) is ensured by the assumption that the parameter α defined in Assumption 2.4 below is positive. We refer for details to [START_REF] Da Prato | Stochastic equations in infinite dimensions[END_REF]Chap. 5] in the context of linear SPDEs with additive noise and to [START_REF] Da Prato | Stochastic equations in infinite dimensions[END_REF]Sect. 7.1] in the context of semilinear SPDEs. Assumption 2.4. Assume that there exists α > 0 such that (-Λ)

α-1 2 Q 1 2 L 2 (H) < ∞, and define α = sup α ∈ [0, 1]; (-Λ) α-1 2 Q 1 2 L 2 (H) < ∞ .
In the sequel, the parameter α plays a key role, indeed it determines the spatial and temporal regularity properties of the solutions, as well as strong rates of convergence of the numerical methods. We observe that α ∈ (0, 1] and for all α ∈ [0, α), (-Λ)

α-1 2 Q 1 2 L 2 (H) < ∞.
We recall that for an equation driven by space-time white noise in dimension d = 1 (Q = I), one has α = 1 2 . However, if d ≥ 2 and Q = I, then (-Λ)

α-1 2 2 L 2 (H) = m≥1 λ α-1 m = ∞ for all α ≥ 0 (see Proposition 2.1), hence the need to consider colored noise (Q = I) if d ≥ 2.
Under Assumption 2.2 and 2.4, for any initial condition u 0 ∈ H and any time T ∈ (0, ∞), the stochastic evolution equation (1) admits a unique global mild solution, which satisfies

u(t) = e tA u 0 + t 0 e (t-s)Λ F (u(s))ds + σ t 0 e (t-s)Λ dW Q (s) , for all t ∈ [0, T ], (6) see

Numerical methods

The aim of this section is to define the semi-discrete and full-discrete approximations of the process u(t) t≥0 . We first describe the spatial discretization, performed using a finite element method. We then detail the fully-discrete scheme, which is the main topic of this article: the temporal discretization is performed using explicit-stabilized integrators. We provide several concrete examples of such schemes using Chebyshev polynomials, and state a general assumption for the analysis. Finally, we state our main convergence results.

Spatial discretization: finite element method

For every mesh size h > 0 (without loss of generality, assume that h ∈ (0, 1]), let T h be a triangulation of the bounded, convex, polyhedral domain D. Assume that the triangulation is quasi-uniform. The finite element spaces V h are finite dimensional linear subspaces of V = H 1 0 (D). Let N h denote the dimension of V h . For simplicity, we consider V h to be the space of continuous functions on D, which are piecewise linear on T h and zero at the boundary ∂D.

For all h ∈ (0, 1], define the linear operator Λ h on V h , equipped with the inner product •, • inherited from H = L 2 (D), as follows: for all

u h , v h ∈ V h , Λ h u h , v h = - D a(x)∇u h (x) • ∇v h (x)dx.
Observe that Λ h is a negative and self-adjoint linear operator on V h , thus there exists an orthonormal basis e m,h 1≤m≤N h of V h and positive real numbers λ m,h 1≤m≤N h such that

Λ h e m,h = -λ m,h e m,h , 1 ≤ m ≤ N h , and λ 1,h ≤ • • • ≤ λ m,h ≤ λ m+1,h ≤ • • • ≤ λ N h ,h .
In addition, for all h > 0, the spectral gap

λ 1,h = inf u h ∈V h -Λ h u h ,u h |u h | 2
of -Λ h is larger or equal than the spectral gap

λ 1 = inf u∈H -Λu,u |u| 2
of -Λ (this is an immediate consequence of the inclusion V h ⊂ H). Define P h : H → V h as the L 2 -orthogonal projection operator onto V h (for the inner product •, • ), for all h ∈ (0, 1]. Then, for all u ∈ H, one has

P h u = N h m=1
u, e m,h e m,h .

In the sequel, the following notation is extensively used. For any function φ : (-∞, 0] → R, and all h ∈ (0, 1], define the self-adjoint linear operator φ(Λ h ) as follows: for all

u h ∈ V h , φ(Λ h )u h = N h m=1 φ(-λ m,h ) u h , e m,h e m,h . (7)
In particular, the linear operators (-Λ h ) α , for all α ∈ [-1, 1], and e tΛ h , for all t ≥ 0, are defined by the expression [START_REF] Abdulle | Weak second order explicit stabilized methods for stiff stochastic differential equations[END_REF] with φ(z) = (-z) α and φ(z) = e tz respectively. Let us state the conditions required for the analysis below. Those conditions are satisfied for piecewise linear finite element methods on a quasi-uniform triangulation of the bounded, convex, polygonal domain D. We refer for instance to [17, Section 3.2] and references therein for details. Assumption 3.1. For all α ∈ [- 1 2 , 1 2 ], there exists C α ∈ (0, ∞) such that for all h ∈ (0, 1) and all

u h ∈ V h , C -1 α |(-Λ h ) α u h | ≤ |(-Λ) α u h | ≤ C α |(-Λ h ) α u h |, and for all u ∈ D (-Λ) α , |(-Λ h ) α P h u| ≤ C α |(-Λ) α u|. Moreover, for every α 1 ∈ [0, 1], α 2 ∈ [α 1 , 2], there exists C α 1 ,α 2 ∈ (0, ∞) such that for all h ∈ (0, 1], (-Λ) α 1 2 I -P h (-Λ) -α 2 2 L(H) ≤ C α 1 ,α 2 h α 2 -α 1 .
It is then straightforward to prove the following result.

Proposition 3.2. Let Assumptions 2.4 and 3.1 be satisfied. Then for all α ∈ [0, α),

sup h∈(0,1] (-Λ h ) α-1 2 P h Q 1 2 L 2 (H) < ∞.
The semi-discrete approximation of (1) obtained by spatial discretization using the finite element method is defined by the following equation: for h ∈ (0, 1],

du h (t) = Λ h u h (t)dt + P h F (u h (t))dt + σP h dW Q (t) , u h (0) = P h u 0 . (8)
This stochastic evolution equation is driven by a Q h -Wiener process, with Q h = P h QP h . The initial condition is given by u h 0 = P h u 0 . The problem is globally well-posed, and its unique solution u h (t) t≥0 takes values in the finite dimensional space V h . It satisfies the following mild formulation

u h (t) = e tΛ h u h 0 + t 0 e (t-s)Λ h P h F (u h (s))ds + σ t 0 e (t-s)Λ h P h dW Q (s), (9)
where we recall that the semi-group is defined by e tΛ h u h = N h m=1 e -tλ m,h u h , e m,h e m,h for all t ≥ 0 and u h ∈ V h using [START_REF] Abdulle | Weak second order explicit stabilized methods for stiff stochastic differential equations[END_REF]. The following smoothing and regularity properties (which are uniform in h ∈ (0, 1]) are used in the sequel: for any γ ∈ [0, 1],

sup h∈(0,1] sup t≥0 t γ (-Λ h ) γ e tΛ h L(V h ) < ∞. ( 10
)
and sup

h∈(0,1] sup t>0 t -γ e tΛ h -I)(-Λ h ) -γ L(V h ) < ∞, (11)
where • L(V h ) denotes the operator norm on L(V h ).

Space-time discretization by explicit stabilized methods

The objective of this section is to describe the temporal discretization of (8) using explicitstabilized integrators. More precisely, we consider families of functions A s s≥1 and B s s≥1 , indexed by the integer parameter s ≥ 1 which is used to describe the stability region of the associated integrators. Let τ > 0 denote the time-step size (without loss of generality, assume that τ ∈ (0, 1]). The integrators studied in this article are defined by

u h n+1 = A s (τ Λ h )u h n + B s (τ Λ h )P h τ F (u h n ) + σ∆W Q n , (12)
for all n ≥ 0, with the initial condition u h 0 = P h u 0 = u h (0). In [START_REF] Eftekhari | Explicit stabilised gradient descent for faster strongly convex optimisation[END_REF], the Wiener increments are defined by ∆W

Q n = W Q (t n+1 ) -W Q (t n ), with t n = nτ .
The self-adjoint linear operators A s (τ Λ h ) and B s (τ Λ h ) are defined using [START_REF] Abdulle | Weak second order explicit stabilized methods for stiff stochastic differential equations[END_REF]:

A s (τ Λ h ) = N h m=1 A s (-τ λ m,h ) •, e m,h e m,h , B s (τ Λ h ) = N h m=1 B s (-τ λ m,h ) •, e m,h e m,h . (13)
From [START_REF] Eftekhari | Explicit stabilised gradient descent for faster strongly convex optimisation[END_REF], one gets a discrete-time mild formulation, similar to ( 6) and ( 9): for all n ≥ 0,

u h n = A s (τ Λ h ) n u h 0 + n-1 k=0 A s (τ Λ h ) n-1-k B s (τ Λ h )P h τ F (u h k ) + σ∆W Q k . ( 14
)
For given time-step size τ and the mesh size h, the parameter s is chosen large enough such that the following stability condition is satisfied:

τ λ N h ,h ≤ L s , (15)
where the size L s ∈ (0, ∞] of the stability domain depends on the choice of A s . One of the minimal requirements (see Assumption 3.3 below for the full list of conditions) is that sup 19) of the explicit stabilized methods (normal and logarithmic scales), with degree s = 7, damping parameter η = 0.05.

z∈[-Ls,0] |A s (z)| ≤ 1. (16)
Note that if A s (z) is a (non-constant) polynomial, then |A s (z)| → ∞ as |z| → ∞
, which yields that L s < ∞ necessarily has a finite value, whereas L s = +∞ for the rational function A(z) = (1-z) -1 related to the implicit Euler method. Several choices are possible for the definition of A s s≥1 and B s s≥1 and will be discussed in the next section.

Chebyshev explicit stabilized integrators

Let us provide examples of explicit stabilized integrators inspired from the literature and which fit in this framework, and discuss implementation details. The first choice is inspired from [START_REF] Abdulle | Optimal explicit stabilized integrator of weak order 1 for stiff and ergodic stochastic differential equations[END_REF] in the context of stiff stochastic integrators with favourable mean-square stability properties, while the second choice in inspired from [START_REF] Eftekhari | Explicit stabilised gradient descent for faster strongly convex optimisation[END_REF] in the context of deterministic partitioned stabilized Runge-Kutta methods for convex optimisation. The idea of explicit stabilized methods is to consider suitably chosen polynomials A s such that the stability domain size L s in ( 16) grows rapidly with s. Under the consistency condition

A s (z) = 1 + z + O(z 2 ) (z → 0), it is known [13, Sect. IV.2]
that the optimal choice maximizing the value of L s is given by A s = T s (1 + z/s 2 ), where T s is the first kind Chebyshev polynomial of order s satisfying T s (cos θ) = cos(sθ), givingL s = 2s 2 . This quadratic growth of L s with respect to s permits to achieve much larger stability domains, in contrast to standard explicit Runge-Kutta methods. However, it turns out, both in the deterministic and stochastic settings [START_REF] Abdulle | Explicit Stabilized Runge-Kutta Methods[END_REF], that the condition ( 16) is not sufficient in general for the design of a robust integrator, and some damping should be introduced to obtain a stronger bound sup

z∈[-Ls,-δ] |A s (z)| ≤ 1 -,
where > 0 and δ > 0 are fixed constants independent of s and the mesh sizes τ, h. A natural choice, as considered in this article, is to consider a fixed damping parameter η > 0 and define A s (z) as

A s (z) = T s (ω 0 + ω 1 z) T s (ω 0 ) , (17)
where T s is the first-kind Chebyshev polynomial of degree s, and the parameters ω 0 , ω 1 are defined by

ω 0 = 1 + η s 2 , ω 1 = T s (ω 0 ) T s (ω 0 ) . ( 18
)
where η is referred to as the damping parameter. The real numbers ω 0 and ω 1 depend on s and η, but the dependence is omitted to simplify notation. Definition ( 17)-( 18) yields a stability domain with length 16), and it can be shown that L s ≥ (2 -4 3 η)s 2 , with a small and fixed damping parameter η > 0. Observe that L s grows quadratically as a function of s, and is arbitrarily close to the optimal value L s = 2s 2 , see [START_REF] Abdulle | Explicit Stabilized Runge-Kutta Methods[END_REF].

L s ≥ 2 ω 1 in (
The choice for the associated function B s is not unique, and also not trivial to be able to take into account the noise source term, which needs to be regularized. We present two families of explicit stabilized integrators.

For the first family, as proposed in [START_REF] Abdulle | Optimal explicit stabilized integrator of weak order 1 for stiff and ergodic stochastic differential equations[END_REF], we use the second kind Chebyshev polynomials defined by T s (x) = sU s-1 (x) and sin(sθ) = sin(θ)U s-1 (cos θ) for all s ≥ 0. Then the SK-ROCK method [START_REF] Abdulle | Optimal explicit stabilized integrator of weak order 1 for stiff and ergodic stochastic differential equations[END_REF] is obtained with choosing A s , B s as follows: for s ≥ 1, set

A s (z) = T s (ω 0 + ω 1 z) T s (ω 0 ) , B s (z) = U s-1 (ω 0 + ω 1 z) U s-1 (ω 0 ) 1 + ω 1 2 z , (19)
where ω 0 , ω 1 are defined in [START_REF] Lord | An introduction to computational stochastic PDEs[END_REF]. It is shown in [START_REF] Abdulle | Optimal explicit stabilized integrator of weak order 1 for stiff and ergodic stochastic differential equations[END_REF] that the nearly optimal stability domain size L s ≥ (2 -4 3 η)s 2 also persists for the size of the mean-square stability domain of SK-ROCK in the context of stiff SDEs, taking advantage of the classical identity on first and second kind Chebyshev polynomials,

T s (x) 2 + U s-1 (x) 2 (1 -x 2 ) = 1.
To illustrate the behavior of these polynomials, in Figure 1, we plot the polynomials A s (z) and B s (z) in [START_REF] Printems | On the discretization in time of parabolic stochastic partial differential equations[END_REF] as a function of real negative z, and we observe that for s = 7, η = 0.05, the corresponding stability domain length L s in ( 16) is close to the optimal value 2s 2 = 98. In addition, the polynomial A s (z) oscillates between the values ±1/T s (ω 0 ) where 1/T s (ω 0 )

1 -η = 0.95, while B s (z) oscillates with a very small amplitude. We also include for comparison the polynomial (A s (z) -1)/z discussed below as an alternative definition of B s (z), and the stability function 1/(1 -z) of the implicit Euler method.

Since A s and B s are polynomials (with degree depending on s), the integrators can be implemented explicitly, in particular no diagonalization or preconditioning techniques of the operator Λ h are required in practice, in constrat to implicit methods. In practice, to avoid a dramatic accumulation of round-off errors, computing the operators A s (τ Λ h ) and B s (τ Λ h ) should note be performed naively as high degree polynomials (see e.g. [START_REF] Abdulle | Explicit Stabilized Runge-Kutta Methods[END_REF]).

Alternatively, taking advantage that first and second kind Chebyshev polynomials have satisfy the same recursion relation: for all s ≥ 1

T s+1 (x) = 2xT s (x) -T s-1 (x), U s+1 (x) = 2xU s (x) -U s-1 (x),
an efficient implementation inspired by [START_REF] Abdulle | Optimal explicit stabilized integrator of weak order 1 for stiff and ergodic stochastic differential equations[END_REF] of the method [START_REF] Eftekhari | Explicit stabilised gradient descent for faster strongly convex optimisation[END_REF] with stability polynomials A s , B s defined in [START_REF] Printems | On the discretization in time of parabolic stochastic partial differential equations[END_REF] can then be achieved as follows. The first method considered in the paper then writes: given u h n , compute u h n+1 by induction as

           K h n,0 = u h n , G h n = P h τ F (u h n ) + σ∆W Q n , K h n,1 = K h n,0 + τ µ 1 Λ h K h n,0 + ν 1 G h n + κ 1 G h n , K h n,i = µ i τ Λ h K h n,i-1 + ν i K h n,i-1 + κ i K h n,i-2 , i = 2, . . . , s, u h n+1 = K h n,s , (20)
with

µ 1 = ω 1 /ω 0 , ν 1 = sω 1 /2, κ 1 = sµ 1
, and for all i = 2, . . . , s,

µ i = 2ω 1 T i-1 (ω 0 ) T i (ω 0 ) , ν i = 2ω 0 T i-1 (ω 0 ) T i (ω 0 ) , κ i = - T i-2 (ω 0 ) T i (ω 0 ) = 1 -ν i . ( 21
)
We now describe the second family of explicit stabilized integrators considered in this article: the polynomials A s and B s are given by

A s (z) = T s (ω 0 + ω 1 z) T s (ω 0 ) , B s (z) = A s (z) -1 z , (22)
a choice considered in [START_REF] Eftekhari | Explicit stabilised gradient descent for faster strongly convex optimisation[END_REF] in the context of convex optimisation. Note that the two methods use the same definition of A s . The second method considered in this paper can be implemented as follows: given u h n , compute u h n+1 by induction as

           K h n,0 = u h n , G h n = P h τ F (u h n ) + σ∆W Q n , K h n,1 = K h n,0 + µ 1 τ Λ h K h n,0 + G h n , K h n,i = µ i τ Λ h K h n,i-1 + G h n + ν i K h n,i-1 + κ i K h n,i-2 , i = 2, . . . , s, u h n+1 = K h n,s , (23)
where µ 1 = ω 1 /ω 0 and µ i , ν i , κ i , i = 2, . . . , s are defined in (21).

General assumptions and main results

We are now in position to state the general assumptions and convergence results of this article. The analysis of convergence of the integrators ( 12) is performed under the following abstract conditions for the family of functions A s s≥1 and B s s≥1 .

Assumption 3.3. For all s ≥ 1, A s and B s are meromorphic functions, and for some sequence L s s≥1 of positive real numbers, the conditions below are satisfied.

A s (0) = A s (0) = B s (0) = 1, for all s ≥ 1, (24)
In addition sup

s≥1,z∈[-Ls,0] (1 + |z|)|B s (z)| 2 < ∞, (25)
and for all δ ∈ (0, 1],

sup s≥1,z∈[-Ls,-δ] |A s (z)| < 1. (26)
Finally, there exists δ ∈ (0, 1] such that sup s≥1,z∈D(0,δ)

|A s (z)| + |B s (z)| < ∞, (27)
where D(0, δ) = {z ∈ C; |z| < δ} denotes the open disc of radius d and center 0 in C.

The strong order of convergence of the explicit stabilized methods [START_REF] Eftekhari | Explicit stabilised gradient descent for faster strongly convex optimisation[END_REF] for the temporal discretization with respect to τ is equal to α/2, analogously to the implicit Euler method (3). The strong order of convergence α/2 is related the Hölder regularity in time of the solution. For all α ∈ [0, α) and T ∈ (0, ∞), there exists C α,T ∈ (0, ∞) such that for any initial condition u 0 ∈ H and all h ∈ (0, 1], τ ∈ (0, 1] and s ≥ 1 such that the stability condition (15) holds, then for all t n = nτ ≤ T one has

E|u h (t n ) -u h n | 2 1 2 ≤ C α,T (1 + |u 0 |t -α 2 n )τ α 2 . ( 28
)
We deduce the strong convergence rates both in time and space, taking into account the finite element spatial discretization. This is an immediate consequence of Theorem 3.4, classical finite element discretization estimates (see Proposition A.2 in Appendix) and the triangular inequality

E|u(t n ) -u h n | 2 1 2 ≤ E|u h (t n ) -u h n | 2 1 2 + E|u(t n ) -u h (t n )| 2 1 2 .
Corollary 3.5. (Strong convergence in time and space) Let Assumptions 2.4 and 3.3 be satisfied. For all α ∈ [0, α) and T ∈ (0, ∞), there exists C α,T ∈ (0, ∞) such that for any initial condition u 0 ∈ H and all h ∈ (0, 1], τ ∈ (0, 1] and s ≥ 1 such that the stability condition (15) holds, then for all t n = nτ ≤ T one has

E|u(t n ) -u h n | 2 1 2 ≤ C α,T (1 + |u 0 |t -α 2 n ) τ α 2 + h α , (29)
for all t n = nτ ≤ T .

Remark 3.6. The framework of our analysis includes not only explicit stabilized methods but also more general explicit or implicit Runge-Kutta type linearized methods. Indeed, observe that the conditions stated in Assumption 3.3 are satisfied if one considers the explicit Euler method (for all s ≥ 1, set A s (z) = 1 + z, B s (z) = 1), with the constraint L s < 2, or the implicit Euler method (for all s ≥ 1, set A s (z) = B s (z) = (1 -z) -1 ), with L s = ∞. However, note that the Crank-Nicolson method (for all s ≥ 1, set A s (z) =

1+ z 2 1-z 2 , B s (z) = 1 1-z 2
) satisfies (24), ( 25) and (27), with L s = ∞, however the condition (26) is not satisfied since |A s (∞)| = 1 (the method is A-stable but not L-stable, see [START_REF] Hairer | Solving ordinary differential equations II. Stiff and differential-algebraic problems[END_REF]).

The convergence estimate of Theorem 3.4 can be generalized in several directions. The proofs would follow the same strategy as for proving Theorem 3.4, with a few modifications described below and omitting some technical details.

In the case of a spatially regular noise, the following strong order one estimate can be obtained, stated for F = 0 for simplicity (the case of a non-zero F is discussed in Remark 4.7). Note that this corresponds formally to α = 2 in Assumption 2.4 and the statement of Theorem 3.7, where τ α/2 can be replaced by τ 1-for arbitrarily small > 0. Note that proving this result requires different techniques from the proof of Theorem 3.4.

Theorem 3.7. (Higher-order of convergence for spatially regular additive noise) Under the assumptions of Theorem 3.4, assume in addition that

sup s≥1,z∈[-Ls,0] min(1, |z|) 1 -A s (z) 2 < ∞. ( 30
)
Consider (1) with F = 0 and assume

(-Λ) 1 2 Q 1 2 L 2 < ∞. ( 31
)
For all ∈ (0, 1) and T ∈ (0, ∞), there exists C ,T ∈ (0, ∞) such that for any initial condition u 0 ∈ H and all h ∈ (0, 1], τ ∈ (0, 1] and s ≥ 1 such that the stability condition (15) holds, then for all t n = nτ ≤ T one has

E|u h (t n ) -u h n | 2 1 2 ≤ C ,T (1 + |u 0 |t -1+ n )τ 1-. (32)
Finally, in the deterministic case (σ = 0), the explicit-stabilized scheme has order of convergence in time equal to 1 and order of convergence in space equal to 2 in the deterministic case σ = 0. The proof is omitted since it follows from straightforward modifications of the proof of Theorem 3.4, using the Hölder regularity of the solution with order 1 -when σ = 0. Theorem 3.8. (Deterministic case) Let Assumption 3.3 be satisfied.

For all ∈ (0, 1) and T ∈ (0, ∞), there exists C ,T ∈ (0, ∞) such that for any initial condition u 0 ∈ H and all h ∈ (0, 1], τ ∈ (0, 1] and s ≥ 1 such that the stability condition (15) holds, then for all t n = nτ ≤ T one has

|u(t n ) -u h n | ≤ C ,T (1 + |u 0 |t -1+ n ) τ 1-+ h 2-. ( 33
)

Convergence analysis

This section is dedicated to the convergence analysis of the considered class of explicit stabilized methods for semilinear parabolic SPDEs. We first prove that the abstract conditions (Assumption 3.3) are indeed satisfied by the considered SK-ROCK method and its variant both described in Section 3.3. We then derive the convergence analysis, based on general consistency results, moment estimates, spatial discretization estimates.

Verification of the abstract conditions for the Chebyshev methods

The goal of this section is to prove Proposition 4.2 below, which states that the two explicitstabilized methods (20) and ( 23), with stability functions ( 19) and ( 22), respectively, satisfy the conditions given in Assumption 3.3. First, we recall some useful asymptotic results used to analyze properties of the explicitstabilized integrators based on the first and second-kind Chebyshev polynomials, see [START_REF] Abdulle | Optimal explicit stabilized integrator of weak order 1 for stiff and ergodic stochastic differential equations[END_REF]. Recall that η is the damping parameter, see [START_REF] Lord | An introduction to computational stochastic PDEs[END_REF]. Proposition 4.1. Let η > 0. Then, when s → ∞, one has

ω 1 s 2 → s→∞ tanh( √ 2η) √ 2η =: Ω(η),
and for all δ ∈ [0, η/Ω(η)],

T s (ω 0 -ω 1 δ) → s→∞ cosh 2(η -δΩ(η)) > 1.
Let η 0 = inf{η > 0; η/Ω(η) = 1}. Observe that η 0 > 0, and that η/Ω(η) < 1 for all η ∈ (0, η 0 ). The value of η 0 has been estimated numerically: η 0 0.7001.

We are now in position to state the main result of this section.

Proposition 4.2. Let A s s=0,1,... be defined by [START_REF] Kruse | Strong and weak approximation of semilinear stochastic evolution equations[END_REF] and B s s=0,1,... be defined either by [START_REF] Printems | On the discretization in time of parabolic stochastic partial differential equations[END_REF] or by (22), with the parameters given by [START_REF] Lord | An introduction to computational stochastic PDEs[END_REF]. Assume that η ∈ (0, η 0 ). Then Assumption 3.3 is satisfied, with L s = 2ω -1 1 .

Proof of Proposition 4.2. Note that A s and B s are polynomial, thus they are analytic functions on C. Proof of (24). This follows from straightforward computations. Proof of (27). We recall the following formula for all s ≥ 1, k ≥ 0,

T (k) s (1) = k-1 j=0 s 2 -j 2 2j + 1 ≤ s 2k . (34)
Moreover, for all x ∈ [1, +∞), one has T (k) s (x) ≥ 0, because all the roots of the polynomial T s (and of its derivatives T (k) s ) belong to the interval (-1, 1). Choose δ = 1. Let s ≥ 1, then T s is a polynomial of degree s, and the Taylor formula yields, for all z ∈ C,

A s (z) = T s (ω 0 + ω 1 z) T s (ω 0 ) = s k=0 1 k!T s (ω 0 ) T (k) s (1)(η/s 2 + ω 1 z) k .
Since the function T s is increasing on [1, +∞) and ω 0 ≥ 1, one has T s (ω 0 ) ≥ T s (1) = 1.

In addition, T s is increasing on [1, ∞), thus using T s (ω 0 ) ≥ T s (1) = s 2 , we obtain ω 1 s 2 ≤ ω 1 T s (ω) = T s (ω 0 ). Finally, using (34) then yields sup

s≥1 ω 1 s 2 ≤ sup s≥1 s k=0 T (k) s (1) s 2k η k k! ≤ e η . ( 35
)
These properties and the estimate (34) yield sup z∈D(0,1)

|A s (z)| ≤ s k=0 1 k! η + ω 1 s 2 k ≤ e η+ω 1 s 2 ≤ e η+e η .
To obtain an upper bound for sup z∈D(0,1)

|B s (z)|, the two versions need to be treated separately.

First, assume that B s is defined by [START_REF] Printems | On the discretization in time of parabolic stochastic partial differential equations[END_REF]. Using the relation T s = sU s-1 , the inequality T s (ω 0 ) ≥ T s (1) = s 2 , and the upper bound

|1 + ω 1 2 z| ≤ 1 + ω 1 2 ≤ 1 + e η
2 ≤ e η when z ∈ D(0, 1), one obtains for all s ≥ 1 sup z∈D(0,1)

|B s (z)| ≤ e η s-1 k=0 1 k!T s (ω 0 ) T (k+1) s (1)(η/s 2 + ω 1 ) k ≤ e η s-1 k=0 1 k! (η + ω 1 s 2 ) k ≤ e 2η+e η ,
where we used a Taylor expansion for the polynomial U s-1 (ω 0 + ω 1 z).

Second, assume that B s is defined by ( 22). Observe that B s is an analytic function (since A s (0) = 1), thus using the maximum principle yields for all s ≥ 1 sup z∈D(0,1)

|B s (z)| = sup |z|=1 A(z) -1 z ≤ 1 + sup z∈D(0,1) |A s (z)| ≤ 1 + e η+e η .
This concludes the proof of (27) for the choice δ = 1.

Proof of (26). For a fixed η, observe that it is sufficient to prove (26) for all δ ∈ (0, ηe -η ]. Indeed, the result then immediately follows as the left-hand side of ( 26) is by definition a decreasing function of δ. We now assume δ ≤ ηe -η . Using (35) yields δ ≤ η/(ω 1 s 2 ) = (ω 0 -1)/ω 1 . On the one hand, if z ∈ [-L s , (1

-ω 0 )/ω 1 ], then x = ω 0 + ω 1 z ∈ [-1, 1], hence |T s (ω 0 + ω 1 z)| ≤ 1. On the other hand, if z ∈ [(1 -ω 0 )/ω 1 , -δ], then x = ω 0 + ω 1 z ≥ 1, hence 0 ≤ T s (ω 0 + ω 1 z) ≤ T s (ω 0 -ω 1 δ) < T s (ω 0 ), since T s is increasing on [1, +∞). Thus, sup z∈[-Ls,-δ] |A s (z)| ≤ max(1, T s (ω 0 -ω 1 δ)) T s (ω 0 ) = A s (-δ) ∈ (0, 1
), for all s ≥ 1.

Moreover, using Proposition 4.1,

A s (-δ) → s→∞ cosh 2(η -δΩ(η)) cosh √ 2η ∈ (0, 1).

One thus obtains the inequality sup s≥1

A s (-δ) ∈ (0, 1), which concludes the proof of (26).

Proof of (25). Owing to (27), which holds with δ = 1, it is sufficient to prove sup

s≥1,z∈[-Ls,-1] |z||B s (z)| 2 < ∞.
It is again necessary to treat separately the two versions of the integrator. First, assume that B s is defined by [START_REF] Printems | On the discretization in time of parabolic stochastic partial differential equations[END_REF]. Recall that

L s = 2ω -1 1 , thus if z ∈ [-2ω -1 1 , -1], then x = ω 0 + ω 1 z ∈ [-1 + η/s 2 , ω 0 -ω 1 ]. Define Q s (x) = s 2 ω 1 T s (ω 0 ) 2 1 + x -η/s 2 2 1 -(x -η/s 2 ) 2 1 -x 2 ,
then a straightforward computation implies the equality

-2zB s (z) 2 = U s-1 (x) 2 (1 -x 2 )Q s (x).
Let us first prove the following claim: for all x ∈ [-1 + η/s 2 , ω 0 -ω 1 ], one has Q s (x) ≥ 0, and sup

s≥1 sup x∈[-1+η/s 2 ,ω 0 -ω 1 ] Q s (x) < ∞. (36)
Indeed, owing to (35), s 2 ω 1 Ts(ω 0 ) 2 ≤ e η . In addition, for all z ∈ [-2ω -1 1 , -1],

1 + x -η/s 2 2 = 1 + ω 1 z 2 ∈ [0, 1].
To treat the last term, note that for all x ∈ R,

d dx 1 -(x -η/s 2 ) 2 1 -x 2 = 2η 1 + x 2 -η/s 2 x s 2 (1 -x 2 ) 2 ≥ 0,
where we used η/s 2 ≤ η 0 ≤ 2. As a consequence, for -1

+ η s 2 ≤ x ≤ ω 0 -ω 1 = 1 + η s 2 -ω 1 , one obtains for all s ≥ 1 0 ≤ Q s (x) ≤ e η 1 -(x -η/s 2 ) 2 1 -x 2 ≤ sup s≥1 e η 1 -(1 -ω 1 ) 2 1 -(ω 0 -ω 1 ) 2 < ∞.
Here, we used the property

e η 1 -(1 -ω 1 ) 2 1 -(ω 0 -ω 1 ) 2 → s→∞ e η 1 -η/Ω(η) ,
and an asymptotic expansion with ω 1 ∼ s→∞ Ω(η)s -2 (owing to Proposition 4.1), ω 0 = 1+ η s 2 , and the property η/Ω(η) < 1 for η < η 0 . This concludes the proof of (36).

Since Q s (x) ≥ 0 for x ∈ [-1 + η/s 2 , ω 0 -ω 1 ], and since z ≤ 0, one obtains that U s-1 (x) 2 (1 -x 2 ) ≥ 0. In addition, properties of the first and second-kind Chebyshev polynomials yield the equality

U s-1 (x) 2 (1 -x 2 ) = 1 -T s (x) 2 ≤ 1 for all x ∈ R. As a consequence, sup s≥1 sup z∈[-2ω -1 1 ,-1] -zB s (z) 2 ≤ sup s≥1 sup x∈[-1+η/s 2 ,ω 0 -ω 1 ] Q s (x) < ∞.
This concludes the proof of (25) for the first version [START_REF] Printems | On the discretization in time of parabolic stochastic partial differential equations[END_REF].

Let us now focus on the second version: assume that B s is defined by (22). Then, using (A s (z) -1)

2 ≤ 2 + 2A s (z) 2 , sup s≥1 sup z∈[-2ω -1 1 ,-1] -zB s (z) 2 = sup s≥1 sup z∈[-2ω -1 1 ,-1] |A s (z) -1| 2 |z| ≤ 2+2sup s≥1 sup z∈[-2ω -1 1 ,-1] |A s (z)| 2 < ∞,
owing to (26) with δ = 1. This concludes the proof of (25) for the two versions. The proof of Proposition 4.2 is thus completed.

Finally, we verify that the considered Chebyshev methods satisfy the additional assumption (30) needed in Theorem 3.7. Proposition 4.3. Let A s s=0,1,... be defined by [START_REF] Kruse | Strong and weak approximation of semilinear stochastic evolution equations[END_REF] with η ∈ (0, η 0 ). Then (30) holds.

Proof. Let δ ∈ (0, ηe -η ). First, (26) yields sup 35)), then one has x ≥ ω 0 -ω 1 δ ≥ 1. Thus, A s (z) and A s (z) are positive and increasing on the interval [-δ, 0]. In addition, using

s≥1,z∈[-Ls,-δ] min(1, |z|) 1 -A s (z) 2 ≤ 1 1 -sup s≥1 sup z∈[-Ls,-δ] |A s (z)| 2 < ∞. Second, let z ∈ [-δ, 0] and x = ω 0 + ω 1 z. Using δ ≤ ηe -η ≤ inf s≥1 ω 0 -1 ω 1 (owing to (
A s (z) -1 z = 1 0 A s (tz)dt ≥ 1 0 A s (-tδ)dt = 1 -A s (-δ) δ ,
and noting 1

-A s (z) 2 ≥ 1 -A s (z), we obtain sup s≥1,z∈[-δ,0] min(1, |z|) 1 -A s (z) 2 ≤ δ 1 -sup s≥1 A s (-δ) < ∞,
where we used (26). This concludes the proof of Proposition 4.3.

General consistency analysis

The objective of this section is to state and prove some consequences of the abstract conditions given in Assumption 3.3. We start with the following lemma, which states that B s (z) = 1 + O(z) and A s (z) n = e nz + O(n -1 ), uniformly for z ∈ [-L s , 0] and s ≥ 1.

Lemma 4.4. Let Assumption 3.3 be satisfied. Then

sup s≥1,z∈[-Ls,0] |1 -B s (z)| |z| < ∞, (37)
and sup

n≥1,s≥1,z∈[-Ls,0] n A s (z) n -e nz < ∞. ( 38
)
The estimate (38) is inspired from [START_REF] Crouzeix | Approximation of parabolic equations[END_REF]Theorem 9] for the time discretization of parabolic PDEs and from [START_REF] Abdulle | Coupling heterogeneous multiscale FEM with Runge-Kutta methods for parabolic homogenization problems: a fully discrete space-time analysis[END_REF]Lemma 5.2] in the context of explicit stabilized methods for parabolic homogenization problems.

|r s (z)| ≤ δ -1 1 + sup z∈[-Ls,-δ] |B s (z)| < ∞.
This concludes the proof of (37). It remains to prove (38). For n ≥ 1 and s ≥ 1, set ϕ n,s (z) = A s (z) n -e nz , and let r 1,s (z) = ϕ 1,s (z)/z 2 . Since A s , and thus ϕ 1,s , are meromorphic functions (by Assumption 3.3), with ϕ 1,s (0) = ϕ 1,s (0) = 0 (owing to (24)), then r 1,s is an analytic function in the open disc D(0, δ), where δ is given by (27) (again, δ does not depend on s). Owing to the maximum principle for analytic functions, one has

C 0 = sup s≥1,|z|≤δ |r 1,s (z)| = sup s≥1 sup |z|=δ |r 1,s (z)| ≤ 1 δ 2 sup s≥1,|z|≤δ (e δ + |A s (z)|) < ∞.
Let γ ∈ (0, δ) be chosen sufficiently small such that γ 8 + C 0 γe γ/2 < 1 2 . Then, for all z ∈ [-γ, 0],

|A s (z)| = e z/2 |e z/2 -e -z/2 ϕ 1,s (z)| ≤ e z/2 (1 -|z|/2 + |z| 2 /8 + C 0 |z| 2 e -z/2 ) ≤ e z/2 .
Using the identity

|ϕ n,s (z)| = |A s (z) n -e nz | = |z -2 ϕ 1,s (z)||z 2 n-1 k=0 e kz A(z) n-k-1 |, one obtains for all n ≥ 1, sup s≥1,z∈[-γ,0] |ϕ n,s (z)| ≤ C 0 |z| 2 ne (n-1)z/2 ≤ 1 n C 0 e γ/
n|ϕ n,s (z)| ≤ sup n≥1 n e -nγ + ρ n < ∞.
Gathering the upper bounds for z ∈ [-L s , -γ] and z ∈ [-γ, 0] concludes the proof of (38). The proof of Lemma 4.4 is thus completed.

Let us now apply the conditions given in Assumption 3.3 and obtained in Lemma 4.4, to deduce properties for the operators A s (τ Λ h ) and B s (τ Λ h ). Recall that • L(V h ) denotes the operator norm on L(V h ). Proposition 4.5. Let Assumption 3.3 be satisfied. For all γ ∈ [0, 1], there exists C ∈ (0, ∞) such that for all τ, h, and for s ≥ 1 satisfying the stability condition (15), then

A(τ Λ h ) L(V h ) ≤ 1, ( 39
) (-Λ h ) γ 2 B(τ Λ h ) L(V h ) ≤ C τ γ 2 , ( 40
)
I -B(τ Λ h ) (-Λ h ) -γ L(V h ) ≤ Cτ γ , (41)
A(τ Λ h ) n -e nτ Λ h L(V h ) ≤ C n , for all n ≥ 1. ( 42
)
Proof. The linear operators A(τ Λ h ), B(τ Λ h ), (-Λ h ) γ , (-Λ h ) -γ , and e nτ Λ h are self-adjoint operators, recall the notation [START_REF] Abdulle | Weak second order explicit stabilized methods for stiff stochastic differential equations[END_REF]. Owing to the stability condition [START_REF] Kelly | Adaptive time-stepping strategies for nonlinear stochastic systems[END_REF], one has -τ λ m,h ∈ [-L s , 0] for all m ∈ {1, . . . , N h }. Inequality (39) is then a straightforward consequence of (26) (which holds for all arbitrarily small δ > 0). Inequality ( 40) is a consequence of (25) (in the cases γ = 0 and γ = 1), and of the following interpolation argument: using the inequalities B s (z) 2 ≤ C 0 and |z|B s (z 41) is a consequence of (37) (in the case γ = 1). Finally, Inequality (42) follows from (38).

) 2 ≤ C 1 , one obtains |z| γ B s (z) 2 ≤ C 1-γ 0 C γ 1 , for γ ∈ [0, 1] and z = -τ λ m,h . Similarly, Inequality (

Proof of the main convergence estimates

We have obtained in the previous section all the ingredients for the proof of Theorem 3.4. The structure of the proof follows the same approach as for the convergence of the linear implicit Euler scheme applied to parabolic semilinear SPDE discretization analysis with Lipschitz continuous nonlinearities, see for instance [START_REF] Printems | On the discretization in time of parabolic stochastic partial differential equations[END_REF]. Our proof however illustrates the behavior of explicit-stabilized integrators in this context.

In the proofs, the value of the constant C ∈ (0, ∞) may change from line to line.

The following estimate on the moments of the fully-discrete scheme is a key ingredient for the proof of the convergence estimates. Proposition 4.6. Let Assumptions 2.2, 2.4 and 3.1 be satisfied.

Let u h n n≥0 be defined by [START_REF] Eftekhari | Explicit stabilised gradient descent for faster strongly convex optimisation[END_REF], where A s s≥1 and B s s≥1 satisfy Assumption 3.3. For all α ∈ [0, α) and T ∈ (0, ∞), there exists C α,T ∈ (0, ∞) such that for all u 0 ∈ H, and all h ∈ (0, 1], τ ∈ (0, 1] and s ≥ 1 such that the stability condition (15) is satisfied, one has max

nτ ≤T E|(-Λ h ) α 2 u h n | 2 ≤ C α,T (1 + |(-Λ) α 2 u 0 | 2 ).
In the proof of Theorem 3.4, only the result with α = 0 is used. Note that Proposition 4.6 shows that the spatial regularity is preserved by the temporal discretization, uniformly in the admissible parameters h, τ, s.

Proof of Proposition 4.6. By the Minkowski inequality, owing to the formulation ( 14) of the fully-discrete scheme, it is sufficient to deal with the three following contributions.

(i) Owing to (39) and using Assumption 3.1, one has

|(-Λ h ) α 2 A s (τ Λ h ) n u h 0 | ≤ |(-Λ h ) α 2 P h u 0 | ≤ C|(-Λ) α 2 u 0 |.
(ii) Owing to the Minkowski inequality, and using the linear growth

|F (x)| ≤ C(1 + |x|) (a consequence of the Lipschitz continuity of F ), one has E (-Λ h ) α 2 τ n-1 k=0 A s (τ Λ h ) n-1-k B s (τ Λ h )P h F (u h k ) 2 1 2 ≤ τ n-1 k=0 E (-Λ h ) α 2 A s (τ Λ h ) n-1-k B s (τ Λ h )P h F (u h k ) 2 1 2 ≤ Cτ n-1 k=0 (-Λ h ) α 2 A s (τ Λ h ) n-1-k B s (τ Λ h ) L(V h ) 1 + E|u h k | 2 1 2 ≤ Cτ n-1 k=0 (-Λ h ) α 2 e (n-1-k)τ Λ h B s (τ Λ h ) L(V h ) 1 + E|u h k | 2 1 2 + Cτ n-1 k=0 A s (τ Λ h ) n-1-k -e (n-1-k)τ Λ h (-Λ h ) α 2 B s (τ Λ h ) L(V h ) 1 + E|u h k | 2 1 2 .
On the one hand, owing to (40) (applied with γ = α), and to [START_REF] Crouzeix | Approximation of parabolic equations[END_REF], one has

τ n-1 k=0 (-Λ h ) α 2 e (n-1-k)τ Λ h B s (τ Λ h ) L(V h ) 1 + E|u h k | 2 1 2 = τ (-Λ h ) α 2 B s (τ Λ h ) L(V h ) 1 + E|u h n-1 | 2 1 2 + τ n-2 k=0 (-Λ h ) α 2 e (n-1-k)τ Λ h B s (τ Λ h ) L(V h ) 1 + E|u h k | 2 1 2 ≤ Cτ 1-α 2 1 + E|u h n-1 | 2 1 2 + Cτ n-2 k=0 1 (n -1 -k)τ α 2 1 + E|u h k | 2 1 2 .
On the other hand, using first (42), then (

A s (τ Λ h ) n-1-k -e (n-1-k)τ Λ h (-Λ h ) α 2 B s (τ Λ h ) L(V h ) 1 + E|u h k | 2 1 2 ≤ Cτ n-2 k=0 1 (n -1 -k) (-Λ h ) α 2 B s (τ Λ h L(V h ) 1 + E|u h k | 2 1 2 ≤ Cτ n-2 k=0 1 (n -1 -k)τ α 2 1 + E|u h k | 2 1 2 . 40), τ n-1 k=0
(iii) Using the Itô isometry formula, the condition (-Λ h )

α+ -1 2 P h Q 1 2 L 2 (H) ≤ C, with ∈ (0, α-α) (see Proposition 3.
2), and the inequalities (10), ( 40) and (42), combined with the same arguments as above, one obtains

E (-Λ h ) α 2 n-1 k=0 A s (τ Λ h ) n-1-k B s (τ Λ h )P h ∆W Q k | 2 = τ n-1 k=0 (-Λ h ) α 2 A s (τ Λ h ) n-1-k B s (τ Λ h )P h Q 1 2 2 L 2 (H) ≤ τ n-1 k=0 (-Λ h ) 1- 2 A s (τ Λ h ) n-1-k B s (τ Λ h ) 2 L(V h ) (-Λ h ) α+ -1 2 P h Q 1 2 2 L 2 (H) ≤ Cτ n-1 k=0 (-Λ h ) 1- 2 e (n-1-k)τ Λ h B s (τ Λ h ) 2 L(V h ) + Cτ n-1 k=0 (-Λ h ) 1- 2 A s (τ Λ h ) n-1-k -e (n-1-k)τ )Λ h B s (τ Λ h ) 2 L(V h ) ≤ C τ + τ n-2 k=0 1 (n -1 -k)τ 1- + Cτ n-2 k=0 1 (n -1 -k)τ 1-(-τ Λ h ) 1- 2 B s (τ Λ h ) 2 L(V h ) ≤ C.
To conclude, first assume that α = 0. The estimate then follows from the application of the discrete Gronwall lemma. The case α ∈ (0, α) then follows from the calculations above. This concludes the proof of Proposition 4.6.

Proof of Theorem 3.4. Introduce the notation t n = nτ , n ∈ N 0 , and assume that

T = N τ , with N ∈ N. Set F h = P h F . Let also n = E|u h (nτ ) -u h n | 2 1 2 .
Using the mild formulations ( 9) and ( 14), one obtains the decomposition

u h (t n ) -u h n = e nτ Λ h -A s (τ Λ h ) n u h 0 + n-1 k=0 t k+1 t k e (tn-t)Λ h -A s (τ Λ h ) n-1-k B s (τ Λ h ) P h dW Q (t) + n-1 k=0 t k+1 t k e (tn-t)Λ h P h F (u h (t)) -A s (τ Λ h ) n-1-k B s (τ Λ h )F h (u h k ) dt.
Using Minkowski inequality, one obtains n ≤ 1 n + 2 n + 3 n , where

1 n = e nτ Λ h -A s (τ Λ h ) n u h 0 2 n = n-1 k=0 E t k+1 t k e (tn-t)Λ h -A s (τ Λ h ) n-1-k B s (τ Λ h ) P h dW Q (t) 2 1 2 3 n = n-1 k=0 t k+1 t k E e (tn-t)Λ h F h (u h (t)) -A s (τ Λ h ) n-1-k B s (τ Λ h )F h (u h k ) 2 1 2 dt.
It remains to prove the three claims below: for all α ∈ [0, α), there exists a constant C such that

1 n ≤ C|u 0 |τ α 2 t α 2 n , ( 43
) 2 n ≤ Cτ α 2 , ( 44
) 3 n ≤ Cτ n-1 k=0 k + Cτ α 2 (1 + |u 0 |). (45)
Proof of (43). This claim follows from (42), indeed this inequality yields for all n ≥ 1

| 1 n | ≤ C n |u h 0 | ≤ C|u 0 |τ α 2 (nτ ) α 2
.

Proof of (44). Using the Itô isometry formula,

( 2 n ) 2 = n-1 k=0 t k+1 t k e (tn-t)Λ h -A s (τ Λ h ) n-1-k B s (τ Λ h ) P h Q 1 2 2 L 2 (H) dt ≤ 3 ( 2,1 n ) 2 + ( 2,2 n ) 2 + ( 2,3 n ) 2 ,
where

( 2,1 n ) 2 = tn 0 e (tn-t)Λ h (I -B s (τ Λ h ))P h Q 1 2 2 L 2 (H) dt ( 2,2 n ) 2 = n-1 k=0 t k+1 t k e (tn-t)Λ h -e (tn-t k+1 )Λ h B s (τ Λ h )P h Q 1 2 2 L 2 (H) dt ( 2,3 n ) 2 = τ n-1 k=0 e (tn-t k+1 )Λ h -A s (τ Λ h ) n-1-k B s (τ Λ h ))P h Q 1 2 2 L 2 (H)
We next prove upper bounds for the quantities ( 2,j n ) 2 , j = 1, 2, 3 as follows.

• Estimate of 2,1 n . Owing to Proposition 3.2, [START_REF] Crouzeix | Approximation of parabolic equations[END_REF], and to (41),

( 2,1 n ) 2 ≤ tn 0 e (tn-t)Λ h (I -B s (τ Λ h ))(-Λ h ) 1-α- 2 2 L(V h ) dt (-Λ h ) α+ -1 2 P h Q 1 2 2 L 2 (H) ≤ C tn 0 (-Λ h ) 1- 2 e (tn-t)Λ h 2 L(V h ) dt (I -B(τ Λ h ))(-Λ h ) -α 2 2 L(V h ) ≤ C tn 0 t -1 dtτ α ≤ Cτ α .
• Estimate of 2,2 n . Owing to Proposition 3.2, [START_REF] Da Prato | Stochastic equations in infinite dimensions[END_REF], and to (40),

( 2,2 n ) 2 ≤ C n-2 k=0 t k+1 t k e (t k+1 -t)Λ h -I)(-Λ h ) -α 2 2 L(V h ) dt e (tn-t k+1 )Λ h (-Λ h ) 1- 2 2 L(V h ) + C tn t n-1 e (tn-t)Λ h -I)(-Λ h ) -α 2 2 L(V h ) dt B s (τ Λ h )(-Λ h ) 1- 2 2 L(V h ) ≤ Cτ α τ n-2 k=0 1 (n -1 -k)τ 1-+ τ ≤ Cτ α .
• Estimate of 2,3 n . Owing to Assumption 2.4, (42), and to (40),

( 2,3 n ) 2 ≤ Cτ n-2 k=0 e (tn-t k+1 )Λ h -A s (τ Λ h ) n-1-k 2 L(V h ) B s (τ Λ h )(-Λ h ) 1--α 2 2 L(V h ) ≤ Cτ n-2 k=0 1 (n -1 -k) 1-τ α+ -1 ≤ Cτ α τ n-2 k=0 1 (n -1 -k)τ 1- ≤ Cτ α .
Proof of (45). The error term 3 n is decomposed as follows:

3 n ≤ 3,1 n + 3,2 n + 3,3 n + 3,4 n ,
where

3,1 n = n-1 k=0 t k+1 t k E e (tn-t)Λ h F h (u h (t)) -F h (u h k ) 2 1 2 dt 3,2 n = n-1 k=0 t k+1 t k E e (tn-t)Λ h I -B s (τ Λ h ) F h (u h k ) 2 1 2 dt 3,3 n = n-1 k=0 t k+1 t k E e (tn-t)Λ h -e (tn-t k+1 )Λ h B s (τ Λ h )F h (u h k ) 2 1 2 dt 3,4 n = τ n-1 k=0 E e (tn-t k+1 )Λ h -A s (τ Λ h ) n-1-k B s (τ Λ h )F h (u h k ) 2 1 2 .
We next estimate the quantities 3,1 n , j = 1, 2, 3, 4 as follows.

• Estimate of 3,1 n . By the Lipschitz continuity of F h (uniformly with respect to the parameter h ∈ (0, 1]), and to the temporal regularity estimate from Proposition A.2, one obtains

3,1 n ≤ Cτ n-1 k=0 E|u h (t k ) -u h k | 2 1 2 + C n-1 k=0 t k+1 t k E|u h (t) -u h (t k )| 2 1 2 dt ≤ Cτ n-1 k=0 k + Cτ α 2 τ 1-α 2 + τ n-1 k=1 (1 + |u h 0 | (kτ ) α 2
) .

• Estimate of 3,2 n . Using the linear growth property of F , Proposition 4.6, and the inequalities (41) and ( 10), one has

3,2 n ≤ C (-Λ h ) -α 2 (I -B s (τ Λ h )) L(V h ) n-1 k=0 t k+1 t k e (tn-t)Λ h (-Λ h ) α 2 L(V h ) dt 1 + E|u h k | 2 1 2 ≤ Cτ α 2 tn 0 dt (t n -t) α 2 (1 + |u 0 |).
• Estimate of 3,3 n . Using the linear growth property of F , Proposition 4.6, inequality (40), and using a combination of ( 10) and [START_REF] Da Prato | Stochastic equations in infinite dimensions[END_REF], then one obtains

3,3 n ≤ C(1 + |u 0 |) n-1 k=0 t k+1 t k e (tn-t)Λ h -e (tn-t k+1 )Λ h L(V h ) dt ≤ C(1 + |u 0 |) n-1 k=0 t k+1 t k e (t k+1 -t)Λ h -I e (tn-t k+1 )Λ h L(V h ) dt ≤ C(1 + |u 0 |)τ α 2 τ 1-α 2 + τ n-2 k=0 C (n -1 -k)τ α 2
.

• Estimate of 3,4 n . Using the linear growth property of F , and the inequality (40), and then using (42), one obtains

3,3 n ≤ C(1 + |u 0 |)τ n-1 k=0 e (tn-t k+1 )Λ h -A s (τ Λ h ) n-1-k L(V h ) ≤ C(1 + |u 0 |)τ α 2 τ 1-α 2 + τ n-2 k=0 C (n -1 -k)τ α 2
.

Combining the error estimates (43), ( 44) and ( 45), the error satisfies 0 = 0 and

n ≤ Cτ n-1 k=0 k + Cτ α 2 (1 + t -α 2 n |u 0 |).
Applying the discrete Gronwall lemma yields the result and this concludes the proof of Theorem 3.4.

The proof of Theorem 3.7 requires a different approach from the proof above of Theorem 3.4 to obtain order 1 of strong convergence for the time discretization of the SPDE (1) driven by additive noise (instead of order at most 1/2 in Theorem 3.4).

Proof of Theorem 3.7. Let us first establish that there exists C ∈ (0, ∞) such that for all s ≥ 1, z ∈ [-L s , 0] and n ≥ 1, one has

|A s (z) n -e nz | ≤ C min(1, n|z| 2 )(|A(z)| n-1 + e (n-1)z ). ( 46
)
On the one hand, sup

s≥1,z∈[-Ls,0] |A(z)| < ∞ by (26), thus |A s (z) n -e nz | ≤ |A s (z)| n + e nz ≤ C|A s (z)| n-1 + e (n-1)z .
On the other hand, using b m -a m ≤ mb m-1 (b -a) for all 0 ≤ a ≤ b and all m ≥ 1, we deduce

|A s (z) n -e nz | ≤ n|A s (z) -e z |(|A s (z)| n-1 + e (n-1)z ) ≤ n|z| 2 |r 1,s (z)|(|A s (z)| n-1 + e (n-1)z ),
where r 1,s (z) = As(z)-e z z 2

. We use the same techniques as in the proof of Lemma 4.4. Let δ ∈ (0, 1). Using (24) and the maximum principle one has

(-Λ h ) 1 2 P h Q 1 2 L 2 (H) ≤ (-Λ) 1 2 Q 1 2 L 2 (H) < ∞,
where we have used the assumption (31).

We then have the following decomposition of the error, which is analogous to that in the proof of Theorem 3.4, but with 3 n = 0 due to F = 0,

E|u h (nτ ) -u h n | 2 1 2 ≤ 1 n + 2 n , with 2 n ≤ 3 ( 2,1 n ) 2 + ( 2,2 n ) 2 + ( 2,3 n ) 2 .

The term 1

n is treated using (42):

| 1 n | ≤ C n |u h 0 | ≤ Cτ |u 0 | t n .
In addition, using the inequalities from Proposition 4.5, the terms 2,1 n and 2,2 n can be treated similarly:

( 2,1 n ) 2 = tn 0 e (tn-t)Λ h (I -B s (τ Λ h ))P h Q 1 2 2 L 2 (H) dt ≤ tn 0 e (tn-t)Λ h (I -B s (τ Λ h ))(-Λ h ) -1 2 2 L(V h ) dt (-Λ h ) 1 2 P h Q 1 2 2 L 2 (H) ≤ C tn 0 (-Λ h ) 1 2 -e tΛ h 2 L(V h ) dt (I -B s (τ Λ h ))(-Λ h ) -1+ 2 L(V h ) ≤ C(T )τ 2(1-) ,
and analogously,

( 2,2 n ) 2 = n-1 k=0 t k+1 t k e (tn-t)Λ h -e (tn-t k+1 )Λ h B s (τ Λ h )P h Q 1 2 2 L 2 (H) dt ≤ C n-2 k=0 t k+1 t k e (t k+1 -t)Λ h -I)(-Λ h ) -1+ 2 L(V h ) dt e (tn-t k+1 )Λ h (-Λ h ) 1 2 -2 L(V h ) + C tn t n-1 e (tn-t)Λ h -I)(-Λ h ) -1+ 2 L(V h ) dt B s (τ Λ h )(-Λ h ) 1 2 -2 L(V h ) ≤ C τ n-2 k=0 τ 2(1-) (n -1 -k)τ 1-+ τ 2 .
It remains to deal with 2,3 n , using different arguments from the proof of Theorem 3.4.

(

2,3 n ) 2 = τ n-1 k=0 e (tn-t k+1 )Λ h -A s (τ Λ h ) n-1-k B s (τ Λ h ))P h Q 1 2 2 L 2 (H) ≤ τ ∞ k=0 N h m=1 |Q 1 2 e m,h | 2 e -kτ λ m,h -A s (-τ λ m,h ) k 2 |B s (-τ λ m,h )| 2 ≤ Cτ N h m=1 |Q 1 2 e m,h | 2 min(1, τ 4 λ 4 m,h ) ∞ k=0 k 2 e -2kτ λ m,h + |A s (-τ λ m,h )| 2k ,
where we have used the claim (46). Since sup

x∈[0,1) (1 -x) 3 ∞ k=0 k 2
x k < ∞ and min(1, y 4 ) ≤ y min(1, y 3 ) for all y ≥ 0, one obtains

( 2,3 n ) 2 ≤ Cτ N h m=1 |Q 1 2 e m,h | 2 min(1, τ 4 λ 4 m,h ) 1 (1 -e -2τ λ m,h ) 3 + 1 (1 -|A s (-τ λ m,h )| 2 ) 3 ≤ Cτ 2 N h m=1 |Q 1 2 (-Λ h ) 1 2 e m,h | 2 min(1, τ λ m,h ) 3 (1 -e -2τ λ m,h ) 3 + min(1, τ λ m,h ) 3 (1 -|A s (-τ λ m,h | 2 ) 3 .
In addition, from assumption (30), we deduce sup

s≥1,z∈[-Ls,0] min(1, |z|) 1 -e -2z + min(1, |z|) 1 -|A s (z)| 2 < ∞.
Finally, this yields

( 2,3 n ) 2 ≤ Cτ 2 N h m=1 |Q 1 2 (-Λ h ) 1 2 e m,h | 2 ≤ Cτ 2 (-Λ) 1 2 Q 1 2 2
L 2 (H) .

Gathering the above estimates concludes the proof of Theorem 3.7.

Remark 4.7. The statement of Theorem 3.7 could be extended for a non zero F , assuming that F is of class C 2 from H to itself with bounded first and second derivatives. In the proof of the main Theorem 3.4 one would need to change the way 3,1 n is treated (with this approach one can not overcome the order 1/2). Performing a second-order Taylor expansion and using the stochastic Fubini theorem would give the result. The details are standard and are omitted for brevity.

Numerical experiments

In this section, we illustrate our convergence analysis on the following stochastic heat equation with additive space-time white noise in dimension d = 1, on the domain D = (0, 1) where we consider homogeneous Dirichlet boundary conditions u(0, t) = u(1, t) = 0 and the initial condition u(x, 0) = sin(2πx). We also consider the stochastic heat equation with multiplicative space-time white noise in dimension d = 1

∂ t u(x, t) = ∂ xx u(x, t) + f (u(x, t)) + g(u(x, t))∂ t W (x, t), (47)
∂ t u(x, t) = ∂ xx u(x, t) + f (u(x, t)) + u(x, t)∂ t W (x, t), (48)
We use a standard finite difference method1 with constant mesh size h = 1/100 for discretizing the Laplacian and the white noise, and obtain the following system where u(x i , t) u i (t) with x i = ih, i = 1, . . . N , h = 1/N ,

du i = u i+1 -2u i + u i-1 h 2 dt + f (u i )dt + 1 √ h g(u i )dw i , i = 1, . . . , N -1,
where w i , i = 1, . . . N are independent Wiener processes, and u 0 = u N = 0 to take into account the homogeneous Dirichlet boundary conditions, with g(u) = 1 for (47) and g(u) = u for (48).

In Figure 2, we plot the convergence curves for the strong error E|u

h (nτ ) -u h n | 2 1 2
both for the additive noise case (47) (see Fig. 2(a)) and the multiplicative noise case (48) (see Fig. 2(b)) using the nonlinearity f (u) = -u -sin(u). We used averages over 10 4 trajectories and a reference solution was computed with the small timestep τ = 0.1 • 2 -14 . We observe in Figure 2 lines of slope 1/4, this corroborates the strong convergence estimate of Theorem 3.4 in the additive noise case and suggests that the strong convergence estimate persists in the multiplicative noise case (recall that α = 1/2 in this case). We see that the error constant for the SK-ROCK method (20) is better than that of the variant (23). We also observe that the linear implicit Euler method has a slightly better error constant in the additive noise case, but not in the multiplicative noise case. Note however that this experiment is meant to illustrate the theory, rather than the performance of explicit stabilized methods which typically reveal efficient in larger spatial dimensions, see e.g. [START_REF] Abdulle | PIROCK: a swiss-knife partitioned implicit-explicit orthogonal Runge-Kutta Chebyshev integrator for stiff diffusion-advection-reaction problems with or without noise[END_REF]. In Figure 3, we plot one realization of the stochastic nonlinear heat equation (47) with f (u) = -u -sin(u) using the two explicit stabilized methods (20), (23) and the linear implicit Euler method (3), respectively. For comparison, we used the same random realizations for sampling the noise. We also plot in Figure 3(d)the corresponding profiles at final time T = 1. It can be observed that compared to the SK-ROCK method (20) the variant method (23) and the linear implicit Euler method (3) exhibit an increased regularity. Note that increasing the damping parameter η of the explicit stabilized methods would increase the regularity of the numerical solutions (this is not illustrated here for brevity). We mention that the question of the spatial regularity of the numerical solution is addressed in [START_REF] Bréhier | High order integrator for sampling the invariant distribution of a class of parabolic stochastic PDEs with additive space-time noise[END_REF]Prop. 3.9], where it is shown that the same regularity as the exact solution can be recovered by introducing a suitable postprocessor for the linear implicit Euler method applied to the stochastic heat equation.

Research Agency. The work of GV was partially supported by the Swiss National Science Foundation, projects No. 200020 184614 and No. 200020 178752. The computations were performed at the University of Geneva on the Baobab cluster.

A Appendix

In this appendix, we first provide spatial and temporal regularity properties on the process u(t) t≥0 and its semi-discrete approximation u h (t) t≥0 . The proofs are omitted since they are quite standard and can be found for instance in [START_REF] Da Prato | Stochastic equations in infinite dimensions[END_REF] or [START_REF] Kruse | Optimal error estimates of Galerkin finite element methods for stochastic partial differential equations with multiplicative noise[END_REF] For any initial condition u 0 ∈ H, there exists a unique global mild solution u(t) t≥0 of (1): the process u is continuous with values in H and satisfies (6) for all t ≥ 0. Moreover, assume that |(-Λ) α 0 2 u 0 | < ∞, with α 0 ∈ [0, α) (recall that α is defined in Assumption 2.4). For every α ∈ [α 0 , α) and T ∈ (0, ∞), there exists a constant C ∈ (0, ∞) such that for all t ∈ (0, T ],

E|(-Λ) α 2 u(t)| 2 ≤ C 1 + |(-Λ) α 0 2 u 0 | 2 t α-α 0 ,
and for all 0 < t 1 ≤ t 2 ≤ T ,

E|u(t 2 ) -u(t 1 )| 2 ≤ C|t 2 -t 1 | α 1 + |(-Λ) α 0 2 u 0 | 2 t α-α 0 1 .
The well-posedness part of the result is obtained by a standard fixed point argument and the computation below with α = 0. For the spatial regularity estimate, note that the stochastic contribution is treated as follows: by the Itô isometry formula, E|(-Λ) L(H) ds (-Λ)

α+ -1 2 Q 1 2 2 L 2 (H) ≤ C t 0 (t -s) -1 ds (-Λ) α+ -1 2 Q 1 2 2 L 2 (H) ,
where ∈ (0, α -α). In addition, the smoothing properties (5) of the semi-group yields (-Λ) For the temporal regularity estimate, similar arguments are used to prove that, first, We now state without proof the following estimate concerning the regularity properties of u h defined by [START_REF] Bréhier | Analysis of some splitting schemes for the stochastic Allen-Cahn equation[END_REF], and the spatial discretization error. We refer for instance to [17, Chap. 2, Thm. 2.27], [START_REF] Lord | An introduction to computational stochastic PDEs[END_REF]Chap. 10,Thm. 10.28], and [START_REF] Kruse | Optimal error estimates of Galerkin finite element methods for stochastic partial differential equations with multiplicative noise[END_REF].

Proposition A.2. Let Assumptions 2.2, 2.4 and 3.1 be satisfied. Assume |(-Λ) α 0 u 0 | < ∞, with α 0 ∈ [0, α). Then, for every α ∈ [α 0 , α) and T ∈ (0, ∞), there exists C ∈ (0, ∞) such that, for all t ∈ (0, T ], sup h∈(0,1] E|(-Λ)

α 2 h u h (t)| 2 ≤ C 1 + |(-Λ) α 0 2 u 0 | 2 t α-α 0 ,
and for all 0 < t 1 ≤ t 2 ≤ T ,

sup h∈(0,1] E|u h (t 2 ) -u h (t 1 )| 2 ≤ C|t 2 -t 1 | α 1 + |(-Λ) α 0 2 u 0 | 2 t α-α 0 1 .
For all α ∈ [α 0 , α) and T ∈ (0, ∞), there exists C ∈ (0, ∞), such that for all t ∈ (0, T ], E u h (t) -u(t)

2 1 2 ≤ Ch α 1 + |(-Λ) α 0 2 u 0 | t α-α 0 2 .
Proposition A.2 states that the strong order of convergence in space is α, which is twice the order α/2 of convergence in time (Theorem 3.4). This is also consistent with the spatial regularity for the process u stated in Proposition A.1.

  2 (D), where D ⊂ R d is an open and convex polyhedral bounded domain in dimension d ∈ {1, 2, 3}. The inner product in H is denoted by •, • . The norm in H is denoted by | • |. The operator norm on the space L(H) of bounded linear operators from H to H is denoted by • L(H) . The Hilbert-Schmidt norm on the space L 2 (H) of Hilbert-Schmidt operators on H is denoted by • L 2 (H) .

Figure 1 :

 1 Figure 1: Stability functions A s (z) and B s (z) in (19) of the explicit stabilized methods (normal and logarithmic scales), with degree s = 7, damping parameter η = 0.05.

  Theorem 3.4. (Strong convergence in time) Let Assumptions 2.4 and 3.3 be satisfied.

Proof of Lemma 4 . 4 .

 44 Let r s (z) = (1 -B s (z))/z. Since B s is a meromorphic function (by Assumption 3.3) with B s (0) = 1 (owing to (24)) and which is bounded on D(0, δ), then r s is an analytic function in the open disc D(0, δ), where δ is given by (27) (recall that δ does not depend on s). Owing to the maximum principle for analytic functions, sup |z|≤δ |r s (z)| = max θ∈[0,2π] |r s (δe iθ )| ≤ C δ In addition, using (25), we deduce sup z∈[-Ls,-δ]

  sup s≥1,z∈[-δ,0] |r 1,s (z)| < ∞.In addition, using (26), one has sups≥1,z∈[-Ls,-δ] |r 1,s (z)| ≤ 1 δ 2 ( sup s≥1,z∈[-Ls,-δ] |A s (z)| + e -δ ) < ∞.This concludes the proof of (46). It remains to prove the error estimate (32). Note that, due to Assumption 3.1, one has the following variant of the result of Proposition 3.2, sup h∈(0,1]

  Additive noise ((47), g(u) = 1).

  Multiplicative noise ((48), g(u) = u).

Figure 2 :

 2 Figure 2: Strong convergence plots using SK-ROCK (20), the variant (23), and the implicit Euler method (3) for the stochastic heat equation with space-time white noise, with final time T = 0.1. Averages over 10 4 samples.

  (a) SK-ROCK method (20). (b) Variant method (23). (c) Linear implicit Euler method (3).

  Solution profile at final time T = 1.

Figure 3 :

 3 Figure 3: Samples of realisation of the stochastic nonlinear heat equation (47) with nonlinearity f (u) = -u -sin(u) and additive noise g(u) = 1 using explicit stabilized methods (20), (23). With time and spatial stepsizes τ = h = 1/100.

α 2 t 0 e

 0 (t-s)Λ dW Q (s)| 2 =

1 - 2 e

 12 (t-s)Λ 2

α 2 e

 2 tΛ u 0 ≤ (-Λ)

t 2 0e (t 2

 22 -s)Λ dW Q (s) -

t 1 0e (t 1 2 ≤ 2 L 2 tα-α 0 2 1 .

 112221 -s)Λ dW Q (s) C|t 2 -t 1 | α ,and, second, using the smoothing properties[START_REF] Abdulle | Coupling heterogeneous multiscale FEM with Runge-Kutta methods for parabolic homogenization problems: a fully discrete space-time analysis[END_REF],| e t 2 Λ -e t 1 Λ u 0 | ≤ e (t 2 -t 1 )Λ -I e t 1 Λ (-Λ) -α 0 (H) |(-Λ) α 0 2 u 0 | ≤ C |t 2 -t 1 | α

  . Let us first recall the following well-posedness and time regularity results, see for instance [17, Chap. 2, Thm. 2.31] and [18, Chap. 10, Thm. 10.26, Thm. 10.27], and [11]. Proposition A.1. Let Assumptions 2.2 and 2.4 be satisfied.

Mathematics Section, École Polytechnique Fédérale de Lausanne, Station 8, 1015 Lausanne, Switzerland, Assyr.Abdulle@epfl.ch

[START_REF] Abdulle | Optimal explicit stabilized integrator of weak order 1 for stiff and ergodic stochastic differential equations[END_REF] Univ Lyon, CNRS, Université Claude Bernard Lyon 1, UMR5208, Institut Camille Jordan, F-69622 Villeurbanne, France. brehier@math.univ-lyon1.

fr[START_REF] Abdulle | Chebyshev methods for stiff stochastic differential equations[END_REF] Université de Genève, Section de mathématiques, 2-4 rue du Lièvre, CP 64, 1211

Genève 4, Switzerland, Gilles.Vilmart@unige.ch

Note that such a standard finite difference method on a uniform mesh can be seen as the simplest finite element method which is thus covered by our analysis.

Acknowledgements. The work of AA was partially supported by the Swiss National Science Foundation, project No. 200020 172710. The work of C.-E. B. was partially supported by the project SIMALIN (ANR-19-CE40-0016) operated by the French National