# Analysis of an exponential integrator for stochastic PDEs driven by Riesz noise

Charles-Édouard Bréhier, David Cohen, Lluís Quer-Sardanyons, Johan Ulander

► To cite this version:

Charles-Édouard Bréhier, David Cohen, Lluís Quer-Sardanyons, Johan Ulander. Analysis of an exponential integrator for stochastic PDEs driven by Riesz noise. 2026. hal-05520387

HAL Id: hal-05520387

https://hal.science/hal-05520387v1

Preprint submitted on 20 Feb 2026

HAL is a multi-disciplinary open access
archive for the deposit and dissemination of sci-
entific research documents, whether they are pub-
lished or not. The documents may come from
teaching and research institutions in France or
abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire HAL, est
destinée au dépôt et à la diffusion de documents
scientifiques de niveau recherche, publiés ou non,
émanant des établissements d'enseignement et de
recherche français ou étrangers, des laboratoires
publics ou privés.

<figure><img src="image_1.png" /></figure>

Distributed under a Creative Commons CC BY 4.0 - Attribution - International License
# ANALYSIS OF AN EXPONENTIAL INTEGRATOR FOR STOCHASTIC PDEs DRIVEN BY RIESZ NOISE

CHARLES-EDOUARD BRÉHIER, DAVID COHEN, LLUÍS QUER-SARDANYONS, AND JOHAN ULANDER

**ABSTRACT.** We present and study an explicit exponential integrator for parabolic SPDEs in any dimension driven by a Gaussian noise which is white in time and with spatial correlation given by a Riesz kernel. Under assumptions on the coefficients of the SPDE, we prove strong error bounds and exhibit how the rate of convergence depends on the exponent in the Riesz kernel. Finally, numerical experiments in spatial dimensions 1 and 2 are provided in order to confirm our convergence results.

**AMS Classification.** 60H15. 60H35. 60M15. 65C20. 65C30. 65J08.

**Keywords.** Stochastic partial differential equations. Stochastic parabolic equations. Stochastic heat equations. Riesz kernel. Stochastic exponential integrator. $L^p(\Omega)$-convergence.

## 1. INTRODUCTION

The analysis of convergence of numerical schemes for stochastic partial differential equations (SPDEs) started with the seminal paper [32]. In this reference, the authors showed the convergence in probability of an implicit numerical scheme to the solution of a stochastic heat equation with a nonlinear term and driven by additive space-time white noise defined on the interval $[0, 1]$. The literature on (strong) convergence of numerical schemes to solutions of SPDEs, in particular of the parabolic type, in dimension 1 is now well established, see e.g [49, 31, 54, 28, 39, 35, 44, 3, 62, 5, 22, 50, 51, 61, 59, 48, 24, 9, 10, 30, 53, 8, 29, 16, 57] and references therein.

Moving beyond the one dimensional case and focusing on the random field setting, convergence results of numerical methods for SPDEs in dimension larger than 1 are currently lacking. Without being exhaustive, we refer the interested reader to [12, 51, 63]. Indeed, the results presented in our paper are closely related to those obtained in [51]. More precisely, in the latter reference the authors study the speed of convergence of explicit and implicit (semi-) discretization schemes for the following semilinear stochastic heat equation driven by Riesz noise with parameter $\alpha \in (0, 2 \wedge d)$, where $d \ge 1$ is the space dimension:

$$ (1) \qquad \begin{cases} \frac{\partial u}{\partial t}(t, x) = \Delta u(t, x) + b(t, x, u(t, x)) + \sigma(t, x, u(t, x))\dot{F}_{\alpha}(t, x), \\ u(0, x) = u_0(x) \quad \text{for } x \in (0, 1)^d, \end{cases} $$

with homogeneous Dirichlet boundary conditions. Here, $\Delta$ is the Laplacian operator on $[0, 1]^d$ with Dirichlet boundary conditions, $\dot{F}_{\alpha}$ denotes a Gaussian noise which is white in time and has a spatial correlation given by a Riesz kernel (see Section 2.2 for the precise definition), $u_0$ is a continuous function on $[0, 1]^d$ and the functions $b, \sigma$ satisfy some conditions (see Section 2.4).

Our main objective is to analyze the rate of strong convergence of an explicit stochastic exponential integrator for the SPDE (1). The main contributions of the present work are twofold:

First, we prove that the strong rate of convergence of the stochastic exponential integrator is $(\frac{1}{2} - \frac{\alpha}{4})^{-}$ (see Theorem 8 for the precise statement). This rate of convergence is optimal, in the sense that it coincides with the time regularity of the solution $u$ of equation (1) (see Proposition 7). Moreover, the obtained rate of convergence coincides with those of the explicit and implicit Euler-Maruyama schemes
considered in [51, Thm. 3.4. (ii) and (iii)]. On top of that, our error estimates are uniform with respect to time and space (see (35)) and require only the minimal assumption of the initial condition $u_0$ being a continuous function. We also point out that the considered stochastic exponential integrator provides two distinct benefits over the numerical schemes studied in [51]: it avoids the CFL-related step-size restriction of the explicit Euler-Maruyama scheme and it has advantages in implementation over the implicit Euler-Maruyama scheme (see Section 5 for numerical illustrations).

The second important contribution of the present paper is that we numerically confirm the theoretically derived strong convergence rates in dimensions $d = 1$ and $d = 2$ (see Section 5). To the best of our knowledge, no previous references in the literature have shown such numerical results in dimension 2 in the random field setting.

Before elaborating a bit more on the above described contributions, let us mention that exponential integrators have a long history in the numerical analysis of deterministic ordinary and partial differential equations, see the review [36]. It is also not the first time that this kind of time integrators have been applied in the context of SPDEs, see for instance [38, 43, 37, 14, 6, 7, 42, 33] and the references above. In addition, we refer the reader to [3], where the authors consider a one-dimensional stochastic heat equation of the form (1) but driven by space-time white noise. The case of stochastic wave equations has been addressed in [19, 60, 58, 20, 2, 55, 18, 21, 46, 11]. Furthermore, in the context of stochastic Schrödinger equations, there have been advances in, e.g., [52, 17, 1, 34, 41, 15, 23].

The proof of our main result (Theorem 8) on the strong rate of convergence for the stochastic exponential integrator is carried out using two important ingredients. On the one hand, we apply well-known estimates for stochastic integrals with respect to our Riesz noise. Here, we adapt the theory developed by Dalang in the seminal paper [25] to our setting, in which we consider our SPDE in the bounded domain $[0, 1]^d$. On the other hand, most of the technical effort has been focused on proving precise estimates involving the Green function $G_d$ of the stochastic heat equation in $[0, 1]^d$ with Dirichlet boundary conditions (see Lemmas 3 and 4 for the precise statements), and using them sharply throughout the proof of our main theorem. The proofs of some of those results, which involve norms of $G_d$ in the Hilbert space determined by the space covariance of the noise $F_\alpha$, have been sketched in [51]. For the sake of completeness, we have added a detailed proof of all estimates in Appendix A. Finally, we remark that, conveniently used, this kind of bounds allows us to get our rate of convergence uniformly with respect to time and space, with only assuming that $u_0$ is continuous.

The paper is organized as follows. In Section 2, we present the main notation that will be used throughout the paper, we describe the Riesz noise $F_\alpha$ together with the corresponding stochastic integral, we provide some properties and estimates involving the Green function $G_d$ and we establish the main assumptions on the coefficients $b$ and $\sigma$. Section 3 is devoted to defining the mild solution to equation (1) and prove some of its properties. In Section 4, the continuous-time mild formulation of the stochastic exponential integrator is presented and the main result of the paper is stated (Theorem 8). Section 5 is devoted to present the numerical experiments in dimensions $d = 1$ and $d = 2$. Finally, in Section 6 we give the proof of Theorem 8. Appendix A contains detailed proofs of the main properties of the Green function $G_d$, while Appendices B and C contain details on the implementation on the covariance matrix of the noise in dimension $d = 1$ and $d = 2$, respectively.

## 2. SETTING

In this section, we first introduce notation in Subsection 2.1, then we provide the definition of the Riesz noise and properties of the associated stochastic integral in Subsection 2.2. Next, we present properties of the heat kernel in Subsection 2.3. Finally, in Subsection 2.4, we state regularity assumptions on the coefficients $b$ and $\sigma$ appearing in the SPDE (1).

2.1. **Notation.** The set of positive integers is denoted by $\mathbb{N} = \{1, 2, \dots\}$. Let $Q = (0, 1)^d$ and $\bar{Q} = [0, 1]^d$, where $d \in \mathbb{N}$ is the spatial dimension. The Euclidean norm in $\mathbb{R}^d$ is denoted by $|\cdot|$. Let $\mathcal{B}(Q)$ denote the set of Borel subsets of $Q$.
Let $C^0(\overline{Q})$ denote the space of real-valued continuous functions on $\overline{Q}$. For any $v \in C^0(\overline{Q})$, set

$$ \|v\|_{\infty} = \max_{x \in \overline{Q}} |v(x)|. $$

Moreover, let

$$ C_0^0(\overline{Q}) = \{v \in C^0(\overline{Q}); v(x) = 0, \quad \forall x \in \partial Q\}. $$

Let $L^2(Q)$ denote the space of square-integrable real-valued functions on $Q$. In addition, let $D(Q)$ denote the space of real-valued smooth functions defined on $Q$.

In this article, we consider the Laplace operator $\Delta$ on the $d$-dimensional domain $Q$ with homogeneous Dirichlet boundary conditions on $\partial Q = \overline{Q} \setminus Q$. It is well-known that this operator can be diagonalized.

For all $\underline{k} = (k_1, \dots, k_d) \in \mathbb{N}^d$, define $|\underline{k}|^2 = \sum_{j=1}^d |k_j|^2$. For all $j \in \mathbb{N}$, set $\varphi_j(\cdot) = \sqrt{2} \sin(j\pi \cdot)$. For all $\underline{k} \in \mathbb{N}^d$, define

$$ \varphi_{\underline{k}}(x) = \prod_{j=1}^{d} \varphi_{k_j}(x_j), \quad \forall x = (x_1, \dots, x_d) \in \overline{Q}. $$

Then $(\varphi_{\underline{k}})_{\underline{k} \in \mathbb{N}^d}$ is a complete orthonormal system of $L^2(Q)$. Moreover, for any $\underline{k} \in \mathbb{N}^d$, the function $\varphi_{\underline{k}}$ is an eigenfunction of the operator $\Delta$, with associated eigenvalue $-|\underline{k}|^2\pi^2$.

Let $(\Omega, \mathcal{F}, \mathbb{P})$ be a probability space. The expectation operator is denoted by $\mathbb{E}[\cdot]$. For any $p \in [1, \infty)$ and any real-valued random variable $u$, we define $\|u\|_p = (\mathbb{E}[|u|^p])^{\frac{1}{p}} \in [0, \infty]$. Let $L^p(\Omega)$ denote the space of real-valued random variables such that $\|u\|_p < \infty$.

In the sequel, $C$ denotes a constant that may vary from line to line. Subscripts on $C$ may be employed to indicate dependence with respect to parameters.

## 2.2. Riesz noise.
In this subsection, we introduce the Riesz noise $F_\alpha$, depending on the parameter $\alpha \in (0, 2 \wedge d)$, and we provide the main tools to define stochastic integrals with respect to the noise $F_\alpha$.

For all $\alpha \in (0, 2 \wedge d)$, define the inner product $\langle \cdot, \cdot \rangle_\alpha$ and the norms $\| \cdot \|_\alpha$ and $\| \cdot \|_{(\alpha)}$ as follows: for any smooth functions $\phi, \phi_1, \phi_2 \in D(Q)$, set

$$ \langle \phi_1, \phi_2 \rangle_\alpha = \iint_{Q \times Q} \phi_1(x) \phi_2(y) |x - y|^{-\alpha} dx dy, $$

$$ \|\phi\|_{\alpha}^{2} = \iint_{Q \times Q} \phi(x) \phi(y) |x - y|^{-\alpha} dx dy, $$

$$ \|\phi\|_{(\alpha)}^{2} = \iint_{Q \times Q} |\phi(x)||\phi(y)||x - y|^{-\alpha} dx dy. $$

Let us introduce the Hilbert space $\mathcal{H}_\alpha$ as the completion of the space $D(Q)$ with respect to the inner product $\langle \cdot, \cdot \rangle_\alpha$ and the associated norm $\| \cdot \|_\alpha$. Any function $g: \overline{Q} \to \mathbb{R}$ satisfying $\|g\|_{(\alpha)}^2 < \infty$ belongs to $\mathcal{H}_\alpha$ (see, e.g., [26, Prop. 2.6(a)]).

Next, define the Hilbert space $\mathbb{H}_\alpha = L^2([0, \infty); \mathcal{H}_\alpha)$ with inner product $\langle \cdot, \cdot \rangle_{\mathbb{H}_\alpha}$ and norm $\| \cdot \|_{\mathbb{H}_\alpha}$ defined as follows: for all $\varphi, \varphi_1, \varphi_2 \in \mathbb{H}_\alpha$, set

$$ \langle \varphi_1, \varphi_2 \rangle_{\mathbb{H}_\alpha} = \int_0^\infty \langle \varphi_1(t, \cdot), \varphi_2(t, \cdot) \rangle_\alpha dt, $$

$$ \|\varphi\|_{\mathbb{H}_\alpha}^2 = \int_0^\infty \|\varphi(t, \cdot)\|_{\alpha}^2 dt. $$
Given the parameter $\alpha \in (0, 2 \wedge d)$, the Riesz noise $F_\alpha$ is defined as a family $F_\alpha = \{F_\alpha(\varphi); \varphi \in \mathbb{H}_\alpha\}$ of centered Gaussian random variables, characterized by the following covariance structure: for all $\varphi_1, \varphi_2 \in \mathbb{H}_\alpha$, one has

$$ \mathbb{E}[F_{\alpha}(\varphi_1)F_{\alpha}(\varphi_2)] = \langle \varphi_1, \varphi_2 \rangle_{\mathbb{H}_{\alpha}}. $$

For any time $s \ge 0$ and any Borel set $A \in \mathcal{B}(Q)$, the function $\varphi_{s,A} : (t, x) \in \mathbb{R}^+ \times Q \mapsto \mathbb{1}_{[0,s]}(t)\mathbb{1}_A(x)$ belongs to $\mathbb{H}_\alpha$. Set

$$ F_{\alpha}([0, s] \times A) = F_{\alpha}(\varphi_{s,A}). $$

For any $t \ge 0$, we define $\mathcal{F}_t$ as the $\sigma$-algebra

$$ \mathcal{F}_t = \sigma(\{F_\alpha([0, s] \times A); 0 \le s \le t, A \in \mathcal{B}(Q)\}). $$

In the sequel, we consider predictable processes with respect to the filtration $\{\mathcal{F}_t\}_{t \ge 0}$.

The stochastic integral with respect to the noise $F_\alpha$ is understood in the Walsh-Dalang sense (see e.g. [56, Chapter 2], [27, Chapter 2] or [25, 26]). Namely, if $\{X(s, y), s \ge 0, y \in \bar{Q}\}$ is any predictable random field satisfying that $\mathbb{E}[\|X\|_{\mathbb{H}_\alpha}^2] < \infty$, then the stochastic integral $\int_0^\infty \int_Q X(s, y)F_\alpha(ds, dy)$ is a well-defined centered random variable, which satisfies the following isometry property:

$$ (2) \qquad \mathbb{E}\left[\left|\int_0^\infty \int_Q X(s, y)F_\alpha(ds, dy)\right|^2\right] = \mathbb{E}[\|X\|_{\mathbb{H}_\alpha}^2]. $$

If $X$, in addition, is such that $\mathbb{E}[\int_0^\infty \|X(s, \cdot)\|_{(\alpha)}^2 ds] < \infty$, then one has

$$ (3) \qquad \mathbb{E}\left[\left|\int_0^\infty \int_Q X(s, y)F_\alpha(ds, dy)\right|^2\right] \le \mathbb{E}\left[\int_0^\infty \|X(s, \cdot)\|_{(\alpha)}^2 ds\right]. $$

For all $\varphi \in \mathbb{H}_\alpha$, one has

$$ \int_0^\infty \int_Q \varphi(s, y) F_\alpha(ds, dy) = F_\alpha(\varphi). $$

Finally, let us recall the Burkholder-Davis-Gundy inequality, see e.g. [40, Appendix B]: for all $p \ge 1$ and $T \in (0, \infty)$, there exists $C_p(T) \in (0, \infty)$ such that for any predictable process $X$ such that $\int_0^T \|X(s, \cdot)\|_{(\alpha)}^2 ds \in L^p(\Omega)$, the stochastic integral satisfies $\int_0^T \int_Q X(s, y)F_\alpha(ds, dy) \in L^{2p}(\Omega)$, and one has

$$ (4) \qquad \mathbb{E}\left[\left|\int_0^T \int_Q X(s, y)F_\alpha(ds, dy)\right|^{2p}\right] \le C_p(T)\mathbb{E}\left[\left(\int_0^T \|X(s, \cdot)\|_{(\alpha)}^2 ds\right)^p\right], $$

which can be rewritten as

$$ \|\|\int_0^T \int_Q X(s, y) F_\alpha(ds, dy)\|\|_{2p}^{2p} \le C_p(T) \|\|\int_0^T \|X(s, \cdot)\|_{(\alpha)}^2 ds\|\|_p^p. $$

**Remark 1.** The Riesz noise $F_\alpha$ is indeed well-defined for all $\alpha \in (0, d)$. The restriction $\alpha < 2$ is due to the fact that we want the following stochastic convolution to be well-defined, for all $T > 0$:

$$ \int_0^T \int_Q G_d(t-s, x, y) F_\alpha(ds, dy), \quad (t, x) \in [0, T] \times \bar{Q}, $$

where $G_d$ is the heat kernel introduced in the next Section 2.3. More precisely, condition $\alpha < 2$ is used to prove that the function $(s, y) \in [0, t] \times \bar{O} \mapsto G_d(t-s, x, y)$ belongs to $\mathbb{H}_\alpha$ (see (12) in Lemma 4). This is important in order to define the mild solution to equation (1) (see (20)).
2.3. **Heat kernel.** In this subsection, we introduce the heat kernel $G_d$ on $Q$ with Dirichlet homogeneous boundary conditions, and provide several of its fundamental properties. For all $t \in (0, \infty)$ and all $x, y \in \bar{Q}$, set

$$G_d(t, x, y) = \sum_{\underline{k} \in \mathbb{N}^d} \exp(-|\underline{k}|^2 \pi^2 t) \varphi_{\underline{k}}(x) \varphi_{\underline{k}}(y),$$

where for all $\underline{k} \in \mathbb{N}^d$, $|\underline{k}|^2$ and $\varphi_{\underline{k}}$ are defined in Subsection 2.1.

Let us first recall standard properties of the heat kernel $G_d$. First, one has

$$(5) \qquad G_d(t, x, y) \ge 0, \quad \forall t \in (0, \infty), x, y \in \bar{Q}.$$

Moreover, for all $t \in (0, \infty)$ and $x \in \bar{Q}$, one has

$$(6) \qquad \int_Q G_d(t, x, y) \, dy \in [0, 1].$$

Next, let us provide estimates for the heat kernel and its first order spatial and temporal derivatives.

**Lemma 2.** The heat kernel $G_d$, its spatial gradient $\nabla_x G_d$ and its temporal derivative $\partial_t G_d$ satisfy the following upper bounds: there exist $c_d, C_d \in (0, \infty)$ such that for all $(t, x, y) \in (0, \infty) \times \bar{Q} \times \bar{Q}$, one has

$$(7) \qquad G_d(t, x, y) \le C_d t^{-\frac{d}{2}} e^{-c_d \frac{|y-x|^2}{2t}},$$

$$(8) \qquad |\nabla_x G_d(t, x, y)| \le C_d t^{-\frac{d+1}{2}} e^{-c_d \frac{|y-x|^2}{2t}},$$

$$(9) \qquad |\partial_t G_d(t, x, y)| \le C_d t^{-\frac{d+2}{2}} e^{-c_d \frac{|y-x|^2}{2t}}.$$

We refer for instance to [45, page 274 in Chapter 4] for a proof of Lemma 2. Combining (7) and (9) yields the following inequality.

**Lemma 3.** There exists $C_d \in (0, \infty)$ such that for all $t_1, t_2 \in (0, \infty)$ such that $t_1 < t_2$, and all $x, y \in \bar{Q}$, one has

$$(10) \qquad \begin{aligned} |G_d(t_2, x, y) - G_d(t_1, x, y)| &\le C_d \mathbf{1}_{\frac{t_2-t_1}{t_1} < 1} \frac{t_2-t_1}{t_1} t_1^{-\frac{d}{2}} e^{-c_d \frac{|y-x|^2}{2t_1}} \\ &\quad + C_d \mathbf{1}_{\frac{t_2-t_1}{t_1} \ge 1} \left( t_1^{-\frac{d}{2}} e^{-c_d \frac{|y-x|^2}{2t_1}} + t_2^{-\frac{d}{2}} e^{-c_d \frac{|y-x|^2}{2t_2}} \right). \end{aligned}$$

In addition, for all $\gamma \in [0, 1]$, for all $t_1, t_2 \in (0, \infty)$ such that $t_1 < t_2$, and all $x, y \in \bar{Q}$, one has

$$(11) \qquad |G_d(t_2, x, y) - G_d(t_1, x, y)| \le C_d \frac{|t_2 - t_1|^{\gamma}}{t_1^{\gamma}} \left( t_1^{-\frac{d}{2}} e^{-c_d \frac{|y-x|^2}{2t_1}} + t_2^{-\frac{d}{2}} e^{-c_d \frac{|y-x|^2}{2t_2}} \right).$$

Let us state additional properties of the heat kernel $G_d$, which are not standard since they are concerned with the norms $\| \cdot \|_{(\alpha)}$. Some of the estimates below can be found in [51, Appendix A]. For the sake of completeness, we add the whole proof of the lemma below in Appendix A.

**Lemma 4.** Let $d \in \mathbb{N}$ and $\alpha \in (0, 2 \wedge d)$. There exists $C_\alpha \in (0, \infty)$ such that the following inequalities hold:

$$(12) \qquad \sup_{x \in \bar{Q}} \|G_d(t, x, \cdot)\|_{(\alpha)}^2 \le C_\alpha t^{-\frac{\alpha}{2}}, \quad \forall t \in (0, \infty),$$

$$(13) \qquad \sup_{x \in \bar{Q}} \int_{t_1}^{t_2} \|G_d(t_2 - t, x, \cdot)\|_{(\alpha)}^2 dt \le C_\alpha(T) |t_2 - t_1|^{1-\frac{\alpha}{2}}, \quad \forall t_2 \ge t_1 \ge 0,$$

$$(14) \qquad \int_{0}^{+\infty} \|G_d(t, x_2, \cdot) - G_d(t, x_1, \cdot)\|_{(\alpha)}^2 dt \le C_\alpha |x_2 - x_1|^{2-\alpha}, \quad \forall x_1, x_2 \in \bar{Q}.$$

Moreover, for all $\gamma \in (0, \frac{1}{2} - \frac{\alpha}{4})$, there exists $C_{\alpha, \gamma} \in (0, \infty)$ such that

$$(15) \qquad \sup_{x \in \bar{Q}} \int_{0}^{t_1} \|G_d(t_2 - t, x, \cdot) - G_d(t_1 - t, x, \cdot)\|_{(\alpha)}^2 dt \le C_{\alpha, \gamma} t_1^{1-\frac{\alpha}{2}-2\gamma} |t_2 - t_1|^{2\gamma}, \quad \forall t_2 \ge t_1 \ge 0.$$
Let us finally provide a technical lemma on estimates for stochastic integrals of the heat kernel $G_d$ with respect to the noise $F_\alpha$.

**Lemma 5.** Let $d \in \mathbb{N}$, $\alpha \in (0, 2 \wedge d)$, $T \in (0, \infty)$, and $p \in [1, \infty)$. There exists $C_{p,\alpha}(T) \in (0, \infty)$ such that the following holds: if $\kappa: \mathbb{R}_+ \to \mathbb{R}_+$ is a mapping such that $\kappa(s) \le s$ for all $s \in \mathbb{R}_+$, and if $\{X(s, x), (s, x) \in [0, T] \times \overline{Q}\}$ is any square-integrable predictable process taking values in $\mathbb{H}_\alpha$, then one has the upper bounds

$$ \sup_{0 \le t \le T} \sup_{x \in \overline{Q}} \left\| \int_0^t \int_Q G_d(t - \kappa(s), x, y) X(s, y) F_\alpha(ds, dy) \right\|_{2p} $$

$$ (16) \qquad \le C_{p,\alpha}(T) \left( \int_0^t (t - \kappa(s))^{-\frac{\alpha}{2}} \sup_{x \in \overline{Q}} \|X(s, x)\|_{2p}^2 ds \right)^{\frac{1}{2}} $$

$$ (17) \qquad \le C_{p,\alpha}(T) \left( \int_0^t (t - \kappa(s))^{-\frac{\alpha}{2}} \sup_{x \in \overline{Q}} \mathbb{E}[|X(s, x)|^{2p}] ds \right)^{\frac{1}{2p}} $$

$$ (18) \qquad \le C_{p,\alpha}(T) \sup_{0 \le t \le T} \sup_{x \in \overline{Q}} \|X(t, x)\|_{2p}. $$

The proofs of Lemma 3, of Lemma 4 and of Lemma 5 are postponed to Appendix A.

**2.4. Regularity assumptions.** We conclude this section by giving the assumptions on the initial value $u_0$ and on the coefficients $b$ and $\sigma$ which appear in the stochastic evolution equation.

**Assumption 1.** The initial value $u_0$ is a function such that $u_0 \in C_0^0(\overline{Q})$. In other words, it is assumed that $u_0: \overline{Q} \to \mathbb{R}$ is continuous, and that $u_0(x) = 0$ for all $x \in \partial Q = \overline{Q} \setminus Q$.

Note that in this article the initial value $u_0$ is assumed to be deterministic. Considering $\mathcal{F}_0$-measurable random initial values, independent of the Riesz noise $F_\alpha$, would be possible up to minor standard modifications.

Let $T \in (0, \infty)$ denote a fixed time. The mappings $b$ and $\sigma$ are assumed to satisfy the following conditions.

**Assumption 2.** The mappings $b: [0, T] \times \overline{Q} \times \mathbb{R} \to \mathbb{R}$ and $\sigma: [0, T] \times \overline{Q} \times \mathbb{R} \to \mathbb{R}$ satisfy the following assumptions: there exist $C \in (0, \infty)$ such that, for all $s, t \in [0, T]$, $x, y \in \overline{Q}$, $u, v \in \mathbb{R}$, one has

$$ (L) \qquad |b(s, x, u) - b(t, y, v)| + |\sigma(s, x, u) - \sigma(t, y, v)| \le C(|s - t|^{1/2 - \alpha/4} + |x - y|^{1 - \alpha/2} + |u - v|), $$
and such that, for all $t \in [0, T]$, $x \in \overline{Q}$, $u \in \mathbb{R}$, one has

$$ (LG) \qquad |b(t, x, u)| + |\sigma(t, x, u)| \le C(1 + |u|). $$

Observe that the condition (L) means that the two mappings $u \in \mathbb{R} \mapsto b(t, x, u)$ and $u \in \mathbb{R} \mapsto \sigma(t, x, u)$ are globally Lipschitz continuous, uniformly with respect to $t \in [0, T]$ and $x \in \overline{Q}$. Moreover, the temporal and spatial Hölder regularity exponents $1/2 - \alpha/4$ and $1 - \alpha/2$ in the condition (L) are related to the optimal temporal and spatial regularity properties of the solutions, see Section 3 and e.g. [51, Prop. 3.2].

### 3. WELL-POSEDNESS AND PROPERTIES OF THE SPDE DRIVEN BY RIESZ NOISE

The objective of this section is to study the stochastic partial differential equation (1) driven by Riesz noise in dimension $d$, employing the tools presented in Section 2. Proposition 6 provides existence and uniqueness of a mild solution as well as moment bounds. Proposition 7 provides a temporal regularity result.
We consider the stochastic partial differential equation

$$
(19) \quad \begin{cases} \frac{\partial u}{\partial t}(t, x) = \Delta u(t, x) + b(t, x, u(t, x)) + \sigma(t, x, u(t, x)) \frac{\partial^2 F_\alpha}{\partial t \partial x}(t, x), & \forall t > 0, x \in Q, \\ u(t, x) = 0, & \forall x \in \partial Q, \\ u(0, x) = u_0(x), & \forall x \in \bar{Q}, \end{cases}
$$

where the initial value $u_0$ satisfies Assumption 1 and the mappings $b$ and $\sigma$ satisfy Assumption 2. The formulation (19) is formal. We recall that a mild solution to (19) is an adapted process $\{u(t,x), (t,x) \in [0,T] \times \bar{Q}\}$ which satisfies the following identity: almost surely, for any $(t,x) \in (0,T] \times \bar{Q}$, one has

$$
(20) \quad u(t,x) = \int_Q G_d(t,x,y)u_0(y) \, dy + \int_0^t \int_Q G_d(t-s,x,y)b(s,y,u(s,y)) \, dy \, ds \\
\qquad + \int_0^t \int_Q G_d(t-s,x,y)\sigma(s,y,u(s,y))F_\alpha(ds,dy).
$$

The stochastic integral with respect to the Riesz noise $F_\alpha(\mathrm{d}s, \mathrm{d}y)$ appearing in (20) is understood in
the Walsh–Dalang sense, as presented in Subsection 2.2.

Let us state properties of mild solutions of the parabolic SPDE (19).

**Proposition 6.** Let Assumptions 1 and 2 be satisfied. For all $T \in (0, \infty)$, there exists a unique global mild solution $\{u(t,x), (t,x) \in [0, T] \times \bar{Q}\}$ to (19).

In addition, the mild solution satisfies the following moment bounds: for all $p \in [1, \infty)$ and all $T \in (0, \infty)$, there exists $C_p(T) \in (0, \infty)$ such that one has

$$
(21) \quad \sup_{t \in [0,T]} \sup_{x \in \bar{Q}} \mathbb{E}[|u(t,x)|^{2p}] \le C_p(T) (1 + \|u_0\|_{\infty}^{2p}).
$$

Proof. Existence of the unique mild solution $u(t,x)$ to the parabolic SPDE (19) can be found in, for example, [25, 13, 47, 27].

The moment bounds of $u(t,x)$ are established using standard arguments. A detailed proof is given for completeness.

For all $t \in [0, T]$ and all $x \in Q$, one has

$$
(22) \quad \mathbb{E}[|u(t,x)|^{2p}] \le 3^{2p-1} (M_p^1(t,x) + M_p^2(t,x) + M_p^3(t,x)),
$$

where $M_p^1(t, x)$, $M_p^2(t, x)$ and $M_p^3(t, x)$ are defined as

$$
\begin{align*}
M_p^1(t,x) &= \left| \int_Q G_d(t,x,y) u_0(y) \, dy \right|^{2p} \\
M_p^2(t,x) &= \mathbb{E} \left[ \left| \int_0^t \int_Q G_d(t-s,x,y) b(s,y,u(s,y)) \, dy \, ds \right|^{2p} \right] \\
M_p^3(t,x) &= \mathbb{E} \left[ \left| \int_0^t \int_Q G_d(t-s,x,y) \sigma(s,y,u(s,y)) F_\alpha(\mathrm{d}s,\mathrm{d}y) \right|^{2p} \right].
\end{align*}
$$

First, owing to the properties (5) and (6) of the heat kernel $G_d$ (see Subsection 2.3), one has

$$
\left| \int_Q G_d(t, x, y) u_0(y) \, dy \right| \le \int_Q G_d(t, x, y) |u_0(y)| \, dy \le \sup_{y \in Q} |u_0(y)| \int_Q G_d(t, x, y) \, dy \le \|u_0\|_{\infty}.
$$

As a result, one obtains

$$
(23) \quad \sup_{t \in [0,T]} \sup_{x \in \bar{Q}} M_p^1(t,x) \le \|u_0\|_{\infty}^{2p}.
$$
Let us then consider $M_p^2(t, x)$. Applying the Hölder inequality and the properties (5) and (6) of the heat kernel $G_d$ (see Subsection 2.3), one obtains

$$ 
\begin{aligned}
M_p^2(t, x) &\le \int_0^t \int_Q G_d(t-s, x, y) \mathbb{E}[|b(s, y, u(s, y))|^{2p}] dy ds \left( \int_0^t \int_Q G_d(t-s, x, y) dy ds \right)^{2p-1} \\
&\le \int_0^t \sup_{y \in \bar{Q}} \mathbb{E}[|b(s, y, u(s, y))|^{2p}] ds t^{2p-1}.
\end{aligned}
$$

Owing to (LG), the mapping $b$ has at most linear growth, therefore one has

$$ \sup_{y \in \bar{Q}} \mathbb{E}[|b(s, y, u(s, y))|^{2p}] \le C_p (1 + \sup_{y \in \bar{Q}} \mathbb{E}[|u(s, y)|^{2p}]). $$

As a result, one obtains

$$ (24) \qquad M_p^2(t, x) \le C_p(T) \int_0^t \left(1 + \sup_{y \in \bar{Q}} \mathbb{E}[|u(s, y)|^{2p}]\right) ds. $$

Finally, let us consider $M_p^3(t, x)$. Applying the inequality (17) from Lemma 5 (with $\kappa(s) = s$ for all $s \ge 0$), one obtains

$$ M_p^3(t, x) \le C_{p,\alpha}(T) \int_0^t (t-s)^{-\frac{\alpha}{2}} \sup_{y \in \bar{Q}} \mathbb{E}[|\sigma(s, y, u(s, y))|^{2p}] ds. $$

Owing to (LG), the mapping $\sigma$ has at most linear growth, therefore one has

$$ \sup_{y \in \bar{Q}} \mathbb{E}[|\sigma(s, y, u(s, y))|^{2p}] \le C_p (1 + \sup_{y \in \bar{Q}} \mathbb{E}[|u(s, y)|^{2p}]). $$

Hence, one gets that

$$ (25) \qquad M_p^3(t, x) \le C_{p,\alpha}(T) \int_0^t (t-s)^{-\frac{\alpha}{2}} \left(1 + \sup_{y \in \bar{Q}} \mathbb{E}[|u(s, y)|^{2p}]\right) ds. $$

Combining the upper bounds (23), (24) and (25), one obtains the following inequality: there exists $C_{p,\alpha}(T) \in (0, \infty)$ such that for all $t \in [0, T]$ one has

$$ \sup_{x \in \bar{Q}} \mathbb{E}[|u(t, x)|^{2p}] \le C_{p,\alpha}(T) \left( \|u_0\|_{\infty}^{2p} + \int_0^t (1 + (t-s)^{-\frac{\alpha}{2}}) \left(1 + \sup_{y \in \bar{Q}} \mathbb{E}[|u(s, y)|^{2p}]\right) ds \right). $$

Applying a Grönwall inequality (e.g., [25, Lem. 15]) then provides the moment bounds (21). This concludes the proof of Proposition 6. $\square$

Let us now state and prove a result on temporal Hölder regularity of the mild solution (20) to (19). The version presented in Proposition 7 has not previously been obtained in the literature, to the best of our knowledge.

**Proposition 7.** Let Assumptions 1 and 2 be satisfied. Let $T \in (0, \infty)$ and $\{u(t, x), (t, x) \in [0, T] \times \bar{Q}\}$ be the unique mild solution to (19). For all $\alpha \in (0, 2 \wedge d)$, $p \in [1, \infty)$ and $\gamma \in (0, \frac{1}{2} - \frac{\alpha}{4})$, there exists $C_{p,\alpha,\gamma}(T) \in (0, \infty)$ such that for all $0 < t_1 \le t_2 \le T$ one has

$$ (26) \qquad \sup_{x \in \bar{Q}} \left( \mathbb{E}[|u(t_2, x) - u(t_1, x)|^p] \right)^{\frac{1}{p}} \le C_{p,\alpha,\gamma}(T) \frac{|t_2 - t_1|^{\gamma}}{t_1^{\gamma}} (1 + \|u_0\|_{\infty}). $$

*Proof.* Assume that $0 < t_1 \le t_2 \le T$ and let $x \in Q$, then one has the decomposition

$$ u(t_2, x) - u(t_1, x) = \delta^1(t_1, t_2, x) + \delta^2(t_1, t_2, x) + \delta^3(t_1, t_2, x) + \delta^4(t_1, t_2, x) + \delta^5(t_1, t_2, x), $$
where the error terms $\delta^j(t_1, t_2, x)$ for all $j \in \{1, \dots, 5\}$ are given by

$$ \delta^1(t_1, t_2, x) = \int_Q [G_d(t_2, x, y) - G_d(t_1, x, y)] u_0(y) \, dy $$

$$ \delta^2(t_1, t_2, x) = \int_0^{t_1} \int_Q [G_d(t_2 - s, x, y) - G_d(t_1 - s, x, y)] b(s, y, u(s, y)) \, dy \, ds $$

$$ \delta^3(t_1, t_2, x) = \int_{t_1}^{t_2} \int_Q G_d(t_2 - s, x, y) b(s, y, u(s, y)) \, dy \, ds $$

$$ \delta^4(t_1, t_2, x) = \int_0^{t_1} \int_Q [G_d(t_2 - s, x, y) - G_d(t_1 - s, x, y)] \sigma(s, y, u(s, y)) F_\alpha(ds, dy) $$

$$ \delta^5(t_1, t_2, x) = \int_{t_1}^{t_2} \int_Q G_d(t_2 - s, x, y) \sigma(s, y, u(s, y)) F_\alpha(ds, dy). $$

Let us consider the term $\delta^1(t_1, t_2, x)$. Assume that $\gamma \in (0, 1)$. Applying the inequality (11) on the heat kernel (see Subsection 2.3), one obtains the following upper bounds:

$$ \begin{aligned} \left| \int_Q [G_d(t_2, x, y) - G_d(t_1, x, y)] u_0(y) \, dy \right| & \le \int_Q |G_d(t_2, x, y) - G_d(t_1, x, y)| |u_0(y)| \, dy \\ & \le \sup_{y \in \bar{Q}} |u_0(y)| \int_Q |G_d(t_2, x, y) - G_d(t_1, x, y)| \, dy \\ & \le C \frac{|t_2 - t_1|^\gamma}{t_1^\gamma} \|u_0\|_\infty. \end{aligned} $$

Therefore one obtains the following upper bound for $\delta^1(t_1, t_2, x)$: for all $\gamma \in (0, 1)$, there exists $C_\gamma \in (0, \infty)$ such that

$$ (27) \qquad \sup_{x \in \bar{Q}} |\delta^1(t_1, t_2, x)| \le C_\gamma \frac{|t_2 - t_1|^\gamma}{t_1^\gamma} \|u_0\|_\infty, \quad \forall 0 < t_1 < t_2. $$

Let us then consider the term $\delta^2(t_1, t_2, x)$. Let $p \in [1, \infty)$. Applying the Minkowski inequality, one has

$$ \begin{aligned} & \left\| \int_0^{t_1} \int_Q [G_d(t_2 - s, x, y) - G_d(t_1 - s, x, y)] b(s, y, u(s, y)) \, dy \, ds \right\|_p \\ & \le \int_0^{t_1} \int_Q |G_d(t_2 - s, x, y) - G_d(t_1 - s, x, y)| \|b(s, y, u(s, y))\|_p \, dy \, ds \\ & \le \sup_{s \in [0, T]} \sup_{y \in \bar{Q}} \|b(s, y, u(s, y))\|_p \int_0^{t_1} \int_Q |G_d(t_2 - s, x, y) - G_d(t_1 - s, x, y)| \, dy \, ds. \end{aligned} $$

Owing to the condition (LG) from Assumption 2, the mapping $b$ has at most linear growth. Therefore, applying the moment bounds (21) from Proposition (6), one obtains the upper bound

$$ \sup_{s \in [0, T]} \sup_{y \in \bar{Q}} \|b(s, y, u(s, y))\|_p \le C(1 + \sup_{s \in [0, T]} \sup_{y \in \bar{Q}} \|u(s, y)\|_p) \le C_p(T)(1 + \|u_0\|_\infty). $$
Applying the inequality (11) on the heat kernel (see Subsection 2.3), one obtains

$$
\begin{align*}
\int_0^{t_1} \int_Q |G_d(t_2 - s, x, y) - G_d(t_1 - s, x, y)| \, dy \, ds
&\le C \int_0^{t_1} \frac{|t_2 - t_1|^\gamma}{(t_1 - s)^\gamma} \int_Q (t_1 - s)^{-\frac{d}{2}} e^{-c_d \frac{|y-x|^2}{2(t_1-s)}} \, dy \, ds \\
&\quad + C \int_0^{t_1} \frac{|t_2 - t_1|^\gamma}{(t_1 - s)^\gamma} \int_Q (t_2 - s)^{-\frac{d}{2}} e^{-c_d \frac{|y-x|^2}{2(t_2-s)}} \, dy \, ds \\
&\le C \int_0^{t_1} \int_Q \frac{|t_2 - t_1|^\gamma}{(t_1 - s)^\gamma} \, ds \\
&\le C_\gamma(T) |t_2 - t_1|^\gamma.
\end{align*}
$$

Therefore one obtains the following upper bound for $\delta^2(t_1, t_2, x)$: for all $p \in [1, \infty)$ and $\gamma \in (0, 1)$, there exists $C_{p,\gamma}(T) \in (0, \infty)$ such that

$$
(28) \quad \sup_{x \in \bar{Q}} \|\delta^2(t_1, t_2, x)\|_p \le C_{p,\gamma}(T) |t_2 - t_1|^\gamma (1 + \|u_0\|_\infty), \quad \forall 0 < t_1 < t_2 \le T.
$$

The treatment of the term $\delta^3(t_1, t_2, x)$ is similar: applying the Minkowski inequality and the property (6) of the heat kernel (see Subsection 2.3), one has

$$
\begin{align*}
\| \int_{t_1}^{t_2} \int_Q G_d(t_2 - s, x, y) b(s, y, u(s, y)) \, dy \, ds \|_p
&\le \int_{t_1}^{t_2} \int_Q G_d(t_2 - s, x, y) \| b(s, y, u(s, y)) \|_p \, dy \, ds \\
&\le \sup_{s \in [0, T]} \sup_{y \in \bar{Q}} \| b(s, y, u(s, y)) \|_p \int_{t_1}^{t_2} \int_Q G_d(t_2 - s, x, y) \, dy \, ds \\
&\le C_p(T) (1 + \|u_0\|_\infty) |t_2 - t_1|.
\end{align*}
$$

Therefore one obtains the following upper bound for $\delta^3(t_1, t_2, x)$: for all $p \in [1, \infty)$ and $\gamma \in (0, 1)$, there exists $C_{p,\gamma}(T) \in (0, \infty)$ such that

$$
(29) \quad \sup_{x \in \bar{Q}} \|\delta^3(t_1, t_2, x)\|_p \le C_{p,\gamma}(T) |t_2 - t_1|^\gamma (1 + \|u_0\|_\infty), \quad \forall 0 < t_1 < t_2 \le T.
$$

It remains to deal with the terms $\delta^4(t_1, t_2, x)$ and $\delta^5(t_1, t_2, x)$. Applying the Burkholder–Davis–Gundy inequality (4), one obtains the upper bound

$$
\begin{align*}
\mathbb{E}[|\delta^4(t_1, t_2, x)|^{2p}] &= \mathbb{E} \left[ \left| \int_0^{t_1} \int_Q [G_d(t_2 - s, x, y) - G_d(t_1 - s, x, y)] \sigma(s, y, u(s, y)) F_\alpha(\mathrm{d}s, \mathrm{d}y) \right|^{2p} \right] \\
&\le C_p(T) \mathbb{E} \left[ \left| \int_0^{t_1} \left\| [G_d(t_2 - s, x, \cdot) - G_d(t_1 - s, x, \cdot)] \sigma(s, \cdot, u(s, \cdot)) \right\|_{(\alpha)}^2 \right|^p \right].
\end{align*}
$$

For all $(t_1, t_2, x) \in [0, T]^2 \times \bar{Q}$ (with $t_1 < t_2$), introduce the auxiliary positive measure $\nu_{t_1, t_2, x}^{(1)}$ on $[0, t_1] \times Q^2$, which is absolutely continuous with respect to the Lebesgue measure, and with Radon–Nikodym derivative given by

$$
\frac{d\nu_{t_1, t_2, x}^{(1)}(s, y, z)}{ds \, dy \, dz} = |G_d(t_2 - s, x, y) - G_d(t_1 - s, x, y)| |G_d(t_2 - s, x, z) - G_d(t_1 - s, x, z)| |y - z|^{-\alpha}.
$$

Making use of the auxiliary measure $\nu_{t_1, t_2, x}^{(1)}$ introduced above yields the upper bound

$$
\mathbb{E}[|\delta^4(t_1, t_2, x)|^{2p}] \le C_p(T) \mathbb{E} \left[ \left| \int_0^{t_1} \iint_{Q \times Q} |\sigma(s, y, u(s, y))| |\sigma(s, z, u(s, z))| \, d\nu_{t_1, t_2, x}^{(1)}(s, y, z) \right|^p \right].
$$
Applying the Hölder inequality, one then obtains

$$
\begin{aligned}
\mathbb{E}[|\delta^4(t_1, t_2, x)|^{2p}] &\le \int_0^{t_1} \iint_{Q \times Q} \mathbb{E}[|\sigma(s, y, u(s, y))|^p |\sigma(s, z, u(s, z))|^p] d\nu_{t_1, t_2, x}^{(1)}(s, y, z) \nu_{t_1, t_2, x}^{(1)}([0, t_1] \times Q^2)^{p-1} \\
&\le \sup_{s \in [0, T]} \sup_{y, z \in \bar{Q}} \mathbb{E}[|\sigma(s, y, u(s, y))|^p |\sigma(s, z, u(s, z))|^p] \nu_{t_1, t_2, x}^{(1)}([0, t_1] \times Q^2)^p.
\end{aligned}
$$

Owing to the condition (LG) from Assumption 2, the mapping $\sigma$ has at most linear growth. Therefore, applying the moment bounds (21) from Proposition (6), one obtains the upper bounds

$$
\begin{aligned}
\sup_{s \in [0, T]} \sup_{y, z \in \bar{Q}} \mathbb{E}[|\sigma(s, y, u(s, y))|^p |\sigma(s, z, u(s, z))|^p] &\le \sup_{s \in [0, T]} \sup_{y \in \bar{Q}} \mathbb{E}[|\sigma(s, y, u(s, y))|^{2p}] \\
&\le C_p(T) (1 + \|u_0\|_{\infty}^{2p}).
\end{aligned}
$$

Moreover, applying the inequality (15) from Lemma 4, with the condition $\gamma \in (0, \frac{1}{2} - \frac{\alpha}{4})$, if $0 < t_1 < t_2 \le T$, one obtains the upper bound

$$
\nu_{t_1, t_2, x}^{(1)}([0, t_1] \times Q^2) = \int_0^{t_1} \|G_d(t_2 - s, x, \cdot) - G_d(t_1 - s, x, \cdot)\|_{(\alpha)}^2 ds \le C_{\alpha, \mu}(T) |t_2 - t_1|^{2\gamma}.
$$

Gathering the upper bounds, one obtains the following upper bound for $\delta^4(t_1, t_2, x)$: for all $p \in [1, \infty)$ and $\gamma \in (0, \frac{1}{2} - \frac{\alpha}{4})$, there exists $C_{p,\alpha,\gamma}(T) \in (0, \infty)$ such that

$$
(30) \quad \sup_{x \in \bar{Q}} \|\delta^4(t_1, t_2, x)\|_{2p} \le C_{p,\alpha,\gamma}(T) |t_2 - t_1|^\gamma (1 + \|u_0\|_\infty), \quad \forall 0 < t_1 < t_2 \le T.
$$

Finally, the treatment of the term $\delta^5(t_1, t_2, x)$ is similar to the treatment of $\delta^4(t_1, t_2, x)$ above.
Applying the Burkholder–Davis–Gundy inequality (4), one obtains the upper bound

$$
\begin{aligned}
\mathbb{E}[|\delta^5(t_1, t_2, x)|^{2p}] &= \mathbb{E} \left[ \left| \int_{t_1}^{t_2} \int_Q G_d(t_2 - s, x, y) \sigma(s, y, u(s, y)) F_\alpha(ds, dy) \right|^{2p} \right] \\
&\le C_p(T) \mathbb{E} \left[ \left| \int_{t_1}^{t_2} \| G_d(t_2 - s, x, \cdot) \sigma(s, \cdot, u(s, \cdot)) \|_{(\alpha)}^2 \right|^p \right].
\end{aligned}
$$

For all $(t_1, t_2, x) \in [0, T]^2 \times \bar{Q}$ (with $t_1 < t_2$), we introduce the auxiliary positive measure $\nu_{t_1, t_2, x}^{(2)}$ on $[t_1, t_2] \times Q^2$, which is absolutely continuous with respect to the Lebesgue measure, and with Radon-Nikodym derivative given by

$$
\frac{d\nu_{t_1, t_2, x}^{(2)}(s, y, z)}{ds dy dz} = G_d(t_2 - s, x, y)G_d(t_2 - s, x, z) |y - z|^{-\alpha}.
$$

Making use of the auxiliary measure $\nu_{t_1, t_2, x}^{(2)}$ introduced above yields the upper bound

$$
\mathbb{E}[|\delta^5(t_1, t_2, x)|^{2p}] \le C_p(T) \mathbb{E} \left[ \left| \int_{t_1}^{t_2} \iint_{Q \times Q} |\sigma(s, y, u(s, y))| |\sigma(s, z, u(s, z))| d\nu_{t_1, t_2, x}^{(2)}(s, y, z) \right|^p \right].
$$
Applying the Hölder inequality, one then obtains

$$
\begin{align*}
& \mathbb{E}\left[\left|\int_{t_1}^{t_2} \iint_{Q \times Q} |\sigma(s, y, u(s, y))| |\sigma(s, z, u(s, z))| \, d\nu_{t_1, t_2, x}^{(2)}(s, y, z) \right|^p\right] \\
& \le \int_{t_1}^{t_2} \iint_{Q \times Q} \mathbb{E}[|\sigma(s, y, u(s, y))|^p |\sigma(s, z, u(s, z))|^p] \, d\nu_{t_1, t_2, x}^{(2)}(s, y, z) \nu_{t_1, t_2, x}^{(2)}([t_1, t_2] \times Q^2)^{p-1} \\
& \le \sup_{s \in [0, T]} \sup_{y, z \in \bar{Q}} \mathbb{E}[|\sigma(s, y, u(s, y))|^p |\sigma(s, z, u(s, z))|^p] \nu_{t_1, t_2, x}^{(2)}([t_1, t_2] \times Q^2)^p.
\end{align*}
$$

Owing to the condition (LG) from Assumption 2, the mapping σ has at most linear growth. Therefore, applying the moment bounds (21) from Proposition (6), one obtains the upper bounds

$$
\begin{align*}
\sup_{s \in [0,T]} \sup_{y,z \in \bar{Q}} \mathbb{E}[|\sigma(s,y,u(s,y))|^p |\sigma(s,z,u(s,z))|^p] &\le \sup_{s \in [0,T]} \sup_{y \in \bar{Q}} \mathbb{E}[|\sigma(s,y,u(s,y))|^{2p}] \\
&\le C_p(T) (1 + \|u_0\|_{\infty}^{2p}).
\end{align*}
$$

Moreover, applying the inequality (15) from Lemma 4, if $0 < t_1 < t_2 \le T$, one obtains the upper bound

$$
\nu_{t_1, t_2, x}^{(2)}([t_1, t_2] \times Q^2) = \int_{t_1}^{t_2} \|G_d(t_2 - s, x, \cdot)\|_{(\alpha)}^2 ds \le C_\alpha(T) |t_2 - t_1|^{1-\frac{\alpha}{2}}.
$$

Gathering the upper bounds, and taking into account the condition $\gamma \in (0, \frac{1}{2} - \frac{\alpha}{4})$, one obtains the following upper bound for $\delta^5(t_1, t_2, x)$: for all $\gamma \in (0, \frac{1}{2} - \frac{\alpha}{4})$ and $p \in [1, \infty)$, there exists $C_{p,\alpha,\gamma}(T) \in (0, \infty)$ such that

$$
(31) \quad \sup_{x \in \bar{Q}} \|\delta^5(t_1, t_2, x)\|_{2p} \le C_{p,\alpha,\gamma}(T) |t_2 - t_1|^\gamma (1 + \|u_0\|_\infty), \quad \forall 0 < t_1 < t_2 \le T.
$$

Combining the upper bounds (27), (28), (29), (30) and (31) then provides the inequality (26). The proof of Proposition 7 is thus completed. □

4. CONVERGENCE ANALYSIS OF THE STOCHASTIC EXPONENTIAL INTEGRATOR

In this section, we introduce the stochastic exponential integrator given by (32) and (34) applied to the stochastic partial differential equation (19). We then state Theorem 8, which is the main result of this article, and provides strong convergence of the integrator when the time-step size $\tau$ goes to 0, with rate $\gamma$ arbitrarily close to $\frac{1}{2} - \frac{\alpha}{4}$.

Let $T \in (0, \infty)$ and given an integer $m \in \mathbb{N}$ define the time-step size $\tau = \frac{T}{m}$. For all $\ell \in \{0, \dots, m\}$, introduce the grid times $t_\ell = \ell \tau$. In addition, let $[\cdot]$ denote the integer part, and for all $s \in [0, T]$ set

$$
\kappa_m^T(s) = \tau \left\lfloor \frac{s}{\tau} \right\rfloor = \frac{T}{m} \left\lfloor \frac{ms}{T} \right\rfloor = \sup \{t_\ell; t_\ell \le s\}.
$$

Let us first start by giving a continuous-time mild formulation of the stochastic exponential integrator which is inspired by [3] which treats the case of one dimensional SPDE driven by space-time white noise. The random field $\{u^m(t, x), (t, x) \in [0, T] \times \bar{Q}\}$ is defined as the solution to the following system: for all $t \in [0, T]$ and $x \in \bar{Q}$,

$$
\begin{equation} \tag{32}
\begin{split}
u^m(t, x) = \int_Q & G_d(t, x, y) u_0(y) \, dy \\
& + \int_0^t \int_Q G_d(t - \kappa_m^T(s), x, y) b(\kappa_m^T(s), y, u^m(\kappa_m^T(s), y)) \, dy \, ds \\
& + \int_0^t \int_Q G_d(t - \kappa_m^T(s), x, y) \sigma(\kappa_m^T(s), y, u^m(\kappa_m^T(s), y)) F_\alpha(ds, dy).
\end{split}
\end{equation}
$$
The formulation (32) of the stochastic exponential integrator is mainly useful for the theoretical analysis of the scheme. For its practical implementation, it is sufficient to consider approximations $\mathcal{U}_{\ell}^{m} = u^{m}(t_{\ell}, \cdot)$ at the grid times $t_{\ell}$. From the definition (32) of the integrator and the semigroup property of the heat kernel $G_d$, for all $\ell \in \{0, \dots, m-1\}$ and all $x \in \overline{Q}$, one has

$$ (33) \qquad \begin{aligned} \mathcal{U}_{\ell+1}^{m}(x) = & \int_{Q} G_d(\tau, x, y) \mathcal{U}_{\ell}^{m}(y) \, dy \\ & + \tau \int_{Q} G_d(\tau, x, y) b(t_{\ell}, y, \mathcal{U}_{\ell}^{m}(y)) \, dy \\ & + \int_{t_{\ell}}^{t_{\ell+1}} \int_{Q} G_d(\tau, x, y) \sigma(t_{\ell}, y, \mathcal{U}_{\ell}^{m}(y)) F_{\alpha}(ds, dy). \end{aligned} $$

Let the operator $e^{\tau \Delta}$ be defined such that for all $v \in \mathcal{C}_0^0(\overline{Q})$ one has

$$ e^{\tau \Delta} v(x) = \int_Q G_d(\tau, x, y) v(y) \, dy, \quad \forall x \in \overline{Q}. $$

In addition, define Gaussian random variables $(\delta_{\ell}F_{\alpha})_{0 \le \ell \le m-1}$ by

$$ \delta_{\ell}F_{\alpha}(\varphi) = F_{\alpha}(1_{[t_{\ell}, t_{\ell+1}]} \otimes \varphi) = \int_{t_{\ell}}^{t_{\ell+1}} \int_{Q} \varphi(y) F_{\alpha}(ds, dy), \quad \forall \varphi \in \mathcal{H}_{\alpha}. $$

Note that if $\ell_1 \ne \ell_2$ then the Gaussian random variables $\delta_{\ell_1}F_{\alpha}$ and $\delta_{\ell_2}F_{\alpha}$ are independent. Moreover, for all $\ell \in \{0, \dots, m-1\}$, $\delta_{\ell}F_{\alpha}$ is a centered Gaussian random variable, with covariance structure given by

$$ \mathbb{E}[\delta_{\ell}F_{\alpha}(\varphi)\delta_{\ell}F_{\alpha}(\psi)] = \tau\langle\varphi, \psi\rangle_{\alpha} = \tau \iint_{Q \times Q} \varphi(x)\psi(y)|x-y|^{-\alpha} \, dx \, dy, \quad \forall \varphi, \psi \in \mathcal{H}_{\alpha}. $$

Employing the notation introduced above, the stochastic exponential integrator is written as

$$ (34) \qquad \mathcal{U}_{\ell+1}^{m} = e^{\tau \Delta} (\mathcal{U}_{\ell}^{m} + \tau b(t_{\ell}, \cdot, \mathcal{U}_{\ell}^{m}) + \sigma(t_{\ell}, \cdot, \mathcal{U}_{\ell}^{m}) \delta_{\ell} F_{\alpha}), \quad \ell = 0, \dots, m-1. $$

The initial value of the scheme is $\mathcal{U}_{0}^{m} = u_0$ for any value $\tau = T/m$ of the time-step size.

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

**Theorem 8.** Let Assumptions 1 and 2 be satisfied. Let $T \in (0, \infty)$. Let $u$ denote the unique mild solution given by (20) to the stochastic partial differential equation (19). For all $m \in \mathbb{N}$, let $u^m$ denote the solution to the stochastic exponential integrator (32) with time-step size $\tau = T/m$.

For all $p \in [1, \infty)$ and $\gamma \in (0, \frac{1}{2} - \frac{\alpha}{4})$, there exists $C_{p,\alpha,\gamma}(T) \in (0, \infty)$ such that for all $m \in \mathbb{N}$ one has

$$ (35) \qquad \sup_{(t,x) \in [0,T] \times \bar{Q}} \mathbb{E} [|u^m(t,x) - u(t,x)|^{2p}]^{\frac{1}{2p}} \le C_{p,\alpha,\gamma}(T) \tau^{\gamma}. $$

The proof of Theorem 8 is postponed to Section 6.

Note that one obtains the same rate of convergence, arbitrarily close to $\frac{1}{2} - \frac{\alpha}{4}$, as for the explicit Euler-Maruyama scheme and the semi-implicit Euler-Maruyama scheme obtained in [51, Theorem 3.4, part (iii)]. However, in Theorem 8 there is no need to impose regularity conditions on the initial value $u_0$ other than continuity.

## 5. NUMERICAL EXPERIMENTS

In this section, we provide several numerical experiments in dimensions 1 and 2 in order to support and illustrate the theoretical results of this paper. In addition, we shall compare the behavior of the analyzed stochastic exponential integrator (34) with the following classical integrators for SPDEs. Let us recall the notation for all of these integrators:
• the stochastic exponential scheme (34) (denoted by SEXP below):

$$ \mathcal{U}_{\ell+1}^{m} = e^{\tau \Delta} (\mathcal{U}_{\ell}^{m} + \tau b(t_{\ell}, \cdot, \mathcal{U}_{\ell}^{m}) + \sigma(t_{\ell}, \cdot, \mathcal{U}_{\ell}^{m}) \delta_{\ell} F_{\alpha}). $$

• the Euler-Maruyama scheme (denoted EM) from [51]:

$$ \mathcal{U}_{\ell+1}^{m} = \mathcal{U}_{\ell}^{m} + \tau \Delta \mathcal{U}_{\ell+1}^{m} + \tau b(t_{\ell}, \cdot, \mathcal{U}_{\ell}^{m}) + \sigma(t_{\ell}, \cdot, \mathcal{U}_{\ell}^{m}) \delta_{\ell} F_{\alpha}. $$

• the semi-implicit Euler-Maruyama scheme (denoted sSEM) from [51]:

$$ \mathcal{U}_{\ell+1}^{m} = \mathcal{U}_{\ell}^{m} + \tau \Delta \mathcal{U}_{\ell+1}^{m} + \tau b(t_{\ell}, \cdot, \mathcal{U}_{\ell}^{m}) + \sigma(t_{\ell}, \cdot, \mathcal{U}_{\ell}^{m}) \delta_{\ell} F_{\alpha}. $$

Let us recall that the noise in the parabolic SPDE (19) is given by a Riesz potential $f(r) = r^{-\alpha}$ with parameter $\alpha \in (0, 2 \wedge d)$.

We start this section by describing the finite difference discretization of the parabolic SPDE (19) from [51] in Subsection 5.1. We then numerically illustrate the profile in Subsection 5.2, the strong rate of convergence in Subsection 5.3, and the computational costs of the integrators (see Subsection 5.4) applied to the SPDE (19) in spatial dimension $d = 1$. Finally, strong rates of convergence and computational costs of the integrators are presented for our parabolic SPDE in spatial dimension $d = 2$ in Subsection 5.5, resp. Subsection 5.6.

**5.1. Spatial discretization of the parabolic SPDE.** In this subsection, we briefly recall the finite difference scheme for the parabolic SPDE (19) studied in [51].

Fix a positive integer $n$ and consider a uniform grid $\frac{k_j}{n} = (\frac{k_1}{n}, \dots, \frac{k_d}{n})$ of $\overline{Q}$, where $k_j \in \{0, \dots, n\}$ for $1 \le j \le d$. For any non-negative integer $i$, set $x_i := \frac{i}{n}$. In dimension $d=1$, a centered finite difference approximation of the Laplacian is provided by the $(n-1) \times (n-1)$ matrix $n^2 D_n$. The matrix $D_n$ is given by

$$ D_n = \begin{pmatrix} -2 & 1 & 0 & \dots & 0 \\ 1 & -2 & 1 & \dots & \vdots \\ 0 & \dots & \dots & \dots & 0 \\ \vdots & \dots & 1 & -2 & 1 \\ 0 & \dots & 0 & 1 & -2 \end{pmatrix}. $$

For $d \ge 2$, the spatial approximation of the Laplacian in the parabolic SPDE (19) is recursively given by

$$ D_n^{(d)} = \text{diag}(D_n^{(d-1)}) + \begin{pmatrix} -2 \text{Id}_{(n-1)^{d-1}} & \text{Id}_{(n-1)^{d-1}} & 0 & \dots & 0 \\ \text{Id}_{(n-1)^{d-1}} & -2 \text{Id}_{(n-1)^{d-1}} & \text{Id}_{(n-1)^{d-1}} & \dots & \vdots \\ 0 & \dots & \dots & \dots & 0 \\ \vdots & \dots & \text{Id}_{(n-1)^{d-1}} & -2 \text{Id}_{(n-1)^{d-1}} & \text{Id}_{(n-1)^{d-1}} \\ 0 & \dots & 0 & \text{Id}_{(n-1)^{d-1}} & -2 \text{Id}_{(n-1)^{d-1}} \end{pmatrix}, $$

starting with $D_n^{(1)} = D_n$, where $\text{Id}_{(n-1)^{d-1}}$ is the $(n-1)^{d-1} \times (n-1)^{d-1}$ identity matrix and $\text{diag}(D_n^{(d-1)})$ denotes the $(n-1)^d \times (n-1)^d$ matrix with $d$ diagonal blocks equal to the $(n-1)^{d-1} \times (n-1)^{d-1}$ matrix $D_n^{(d-1)}$. We denote by $u^n(t)$ the $(n-1)^d$-dimensional vector approximating the solution at time $t$ to the SPDE (19) on the grid points. For any $k_1, \dots, k_d \in \{1, \dots, n-1\}$, we define the non-negative integer $\mathbf{k} = (k_d - 1)(n-1)^{d-1} + \dots + (k_2 - 1)(n-1) + k_1$ and set $u^n(t, \mathbf{x}_k) := (u^n(t)_{\mathbf{k}})_{\mathbf{k}=1}^{(n-1)^d}$, with $\mathbf{x}_k := (x_{k_1}, \dots, x_{k_d})$, which corresponds to a spatial approximation of $u(t, \mathbf{x}_k)$. Then, the vector $u^n(t)$ satisfies the stochastic differential equation

$$ (36) \qquad \mathrm{d}u^n(t) = n^2 D_n^{(d)} u^n(t) \mathrm{d}t + b(u^n(t)) \mathrm{d}t + \sigma(u^n(t)) \mathrm{d}F_\alpha^n(t), $$
with initial values $u^n(0)_k = u_0(\mathbf{x}_k)$. Here, the following classical convention is used: for any function $h: \mathbb{R} \to \mathbb{R}$ and vector $v \in \mathbb{R}^r$, we write $h(v) := (h(v_1), \dots, h(v_r))$ and interpret the last product in (36) in the elementwise sense. The vector $F_\alpha^n(t)$ is given by

$$F_{\alpha}^{n}(t)_{\mathbf{k}} := n^{d} \int_{0}^{t} \int_{\Box_{\mathbf{x}_{\mathbf{k}}}} F_{\alpha}(ds, dy),$$

where $\Box_{\mathbf{x}_k} = \{x = (x_1, \dots, x_d); x_{k_j} \le x \le x_{k_{j+1}}, \forall 1 \le j \le d\}$.

Note that, in order to apply a time integrator to the SDE (36), one needs to know the covariance matrix of the discretized noise. An explicit formula can be obtained in the 1-dimensional case, see Appendix B for details. To the best of our knowledge, this is not the case in the 2-dimensional case, where we employ quadrature formulas, see Appendix C for details.

## 5.2. Time evolution and profile of the solution in dimension 1.
We consider the semilinear stochastic heat equation (19) on the interval $[0, 1]$ with homogeneous Dirichlet boundary conditions and in the time interval $[0, 0.5]$. For this numerical experiment, we choose $b(u) = \sigma(u) = 1 + 0.5 \cos(u)$ and the initial value $u_0(x) = \sin(\pi x)$. We illustrate the behavior of the solution for three values of the parameter $\alpha$ of the Riesz noise: $\alpha = 0.2, 0.7$ and $\alpha = 1$ (in order to denote a space-time white noise (which is not a Riesz noise) as in [51]). A proof of convergence of the stochastic exponential scheme for space-time white noise is given in [3] for instance. The finite difference method is employed with $n = 2^{10}$ number of grid points and the stochastic exponential integrator with $\tau = 2^{-16}$. We display the time evolution as well as a profile of the numerical solution at $T_{end} = 0.5$ in Figure 1. In this figure, one can observe that the spatial and temporal regularity depends on $\alpha$, see Proposition 7.

## 5.3. Strong convergence in dimension 1.
In this subsection we illustrate the strong convergence of the stochastic exponential integrator (34) as stated in Theorem 8.

We consider the semilinear SPDE (19) for the time interval $[0, 0.5]$. We consider $b(u) = \sigma(u) = 1 + 0.5 \cos(u)$ and the initial value $u_0(x) = \sin(\pi x)$. The parameter for the Riesz noise is set to be $\alpha = 0.7$. In addition, we use 100 samples to approximate the expectations. We have verified that this is enough for the Monte Carlo errors to be negligible. We apply the finite difference method with $n = 2^9$ and the stochastic exponential integrator with $\tau_{ref} = 2^{-20}$ to produce the reference solution $u_{ref}$. The strong errors

$$\sup_{(t_j, x_k) \in [0, 0.5] \times [0, 1]} \mathbb{E} \left[ |u^{j,k} - u_{ref}(t_j, x_k)|^2 \right]$$

of the considered integrators are presented in Figure 2 (left). Here, $u^{j,k}$ denotes a numerical approximation of $u(t_j, x_k)$. An order of convergence $1 - \alpha/2 = 0.65$ for the proposed time integrator (34) is observed. This is in agreement with the results of Theorem 8.

We now illustrate the dependence of the order of convergence of the stochastic exponential integrator with respect to the parameter $\alpha$ and consider a smoother noise with the parameter $\alpha = 0.2$ (the other parameters for this simulation are as above). The results are presented in Figure 2 (right). Strong order of convergence $1 - \alpha/2 = 0.9$ is observed for the stochastic exponential integrator (34). This is in agreement with Theorem 8. Figure 2 also illustrates the rates of convergence of the semi-implicit Euler-Maruyama scheme, as proved in the work [51], as well as the (well-known) fact that the classical Euler-Maruyama scheme has a severe step size restriction when applied to (S)PDEs.

We conclude this section by further illustrating the dependence of the strong order of convergence of the proposed time integrator with respect to the parameter of the noise $\alpha$.
<figure><img src="image_2.png" /><figcaption>FIGURE 1. Time evolution (left) and profile at $T_{end} = 0.5$ (right) for different values of the parameter $\alpha$.</figcaption></figure>

We perform another numerical experiment with the following parameters: $T = 0.5$, $u_0(x) = \sin(\pi x)$, $b(u) = 2 + \sin(u)$, $\sigma(u) = u + 2$, $n = 2^8$, $\tau_{ref} = 2^{-20}$, $\tau = 2^{-9}, \dots, 2^{-16}$, 150 Monte Carlo samples (this is enough for the Monte Carlo errors to be negligible), and $\alpha = 0.2, 0.4, 0.6, 0.8$. The results are presented in Figure 3 and further illustrate Theorem 8.

5.4. **Computational times in dimension 1.** In this subsection we compare the computational cost of the above time integrators. To this end, we consider the SPDE (19). We run 50 samples for each numerical scheme. For each numerical integrator and each sample, we use several time-step sizes ($2^{-\ell}\tau_{ref}$ for $\ell = 1, 2, \dots, 10$) and compare the strong error at the final time $T = 0.75$ with a reference solution computed with the same sample of the noise and by the same numerical scheme for the time-step size $\tau_{ref} = 2^{-20}$. Figure 4 displays the total computational time for all the samples, for each numerical scheme and each time-step size, as a function of the averaged final error. In this figure, one observes better performance for the stochastic exponential integrator (34) than for the two classical time integrators from [51]. In addition, one can observe a time-step restriction for the explicit Euler-Maruyama scheme (top right of the figure).

5.5. **Strong convergence in dimension 2.** In this subsection, we confirm the theoretical result in Theorem 8 in dimension $d = 2$. We consider the semilinear SPDE (19) on the spatial
<figure><img src="image_3.png" /><figcaption>FIGURE 2. Strong errors for the SPDE (19) with Riesz potential $f(r) = r^{-\alpha}$.</figcaption></figure>

<figure><img src="image_4.png" /><figcaption>FIGURE 3. Strong errors for the SPDE (19) with Riesz potential $f(r) = r^{-\alpha}$ for $\alpha = 0.2, 0.4, 0.6, 0.8$.</figcaption></figure>

domain $[0, 1]^2$ with homogeneous Dirichlet boundary conditions and for the time interval $[0, 1]$. We consider $b(u) = \sigma(u) = 1 + \cos(u)$ and the initial value $u_0(x, y) = \sin(2\pi x) \sin(2\pi y)$ for $x, y \in [0, 1]$.

We compute the strong errors

$$ \left( \sup_{(t_k, x_i, x_j) \in [0,1] \times [0,1] \times [0,1]} \mathbb{E} \left[ |u^{k,i,j} - u_{ref}(t_k, x_i, x_j)|^2 \right] \right)^{1/2}, $$

where the reference solution $u_{ref}$ is computed using the stochastic exponential integrator with temporal discretization size $\tau_{ref} = 2^{-13}$. We use 100 samples to approximate the expectation in the strong errors and we have verified that this is enough for the Monte Carlo errors to be negligible.
<figure><img src="image_5.png" /><figcaption>FIGURE 4. Computational time as a function of the averaged final error for the three numerical methods.</figcaption></figure>

In the first numerical experiment, we fix the space discretization parameter to $n = 2^6 = 64$ and take the following values of the noise correlation parameter $\alpha = 0.2, 0.4, 0.6, 0.8, 1, 1.2, 1.4, 1.6, 1.8$. As is seen in Figures 5 and 6, the decay of the strong errors follow the reference lines with slopes $1/2 - \alpha/4$. This confirms the results derived in Theorem 8.

<figure><img src="image_6.png" /><figcaption>FIGURE 5. Root mean square convergence of the stochastic exponential Euler scheme (34) for $\alpha = 0.2, 0.4, 0.6, 0.8$ and for $n = 64$ space grid point in each direction.</figcaption></figure>

In the second numerical experiment, presented in Figure 7, we fix the noise parameter to $\alpha = 0.8$ and use different values of the space discretization parameter $n = 16, 32, 64$ to confirm that the error in Theorem 8 is independent of $n$. In addition, we show that the convergence rate agrees with the derived theoretical rate $1/2 - \alpha/4$.

5.6. **Computational times in dimension 2.** In this subsection we compare the computational costs of the stochastic exponential Euler scheme and of the semi-implicit Euler-Maruyama scheme in dimension $d=2$. We consider the SPDE (19) with the same parameters as in Subsection 5.5. We apply the finite difference discretization with mesh size $n = 2^5 = 32$ and these two integrators with time-step sizes $2^{-l}$ for $l = 2, \dots, 2^{18}$. The strong error is computed using $\tau_{ref} = 2^{-19}$ for the reference solution and using 100 Monte Carlo samples to approximate the expectations.
<figure><img src="image_7.png" /><figcaption>FIGURE 6. Root mean square convergence of the stochastic exponential Euler scheme (34) for $\alpha = 1, 1.2, 1.4, 1.6, 1.8$ and for $n = 64$ space grid points.</figcaption></figure>

<figure><img src="image_8.png" /><figcaption>FIGURE 7. Root mean square convergence of the stochastic exponential Euler scheme (34) for $\alpha = 0.8$ and for $n = 16, 32, 64$ space grid points.</figcaption></figure>

Figure 8 displays the total computational time for the SEXP and sEM schemes for all samples; meaning that, only the parts of the implementations of the SEXP and sEM schemes that are different are included in the elapsed times. In this figure, one can observe better performance for the SEXP scheme compared to the sEM scheme in the regime of small errors.

## 6. PROOF OF THE MAIN RESULT

This section gives the proof of Theorem 8.

*Proof.* Recall that the random field $u$ is the mild solution given by (20) to the stochastic partial differential equation (19), and that, for all $m \in \mathbb{N}$, the random field $u^m$ is the numerical solution given by the stochastic exponential integrator (32) with time-step size $\tau = T/m$.

For all $p \in [1, \infty)$, define

$$ \begin{aligned} \mathcal{E}_{2p}^{m}(t, x) &= \|u(t, x) - u^{m}(t, x)\|_{2p} = (\mathbb{E}[|u(t, x) - u^{m}(t, x)|^{2p}])^{\frac{1}{2p}}, & \forall (t, x) \in [0, T] \times \overline{Q}, \\ \overline{\mathcal{E}}_{2p}^{m}(t) &= \sup_{x \in \overline{Q}} \mathcal{E}_{2p}^{m}(t, x), & \forall t \in [0, T]. \end{aligned} $$
<figure><img src="image_9.png" /><figcaption>FIGURE 8. Computational time as a function of the averaged final error for the two time integrators sEM and SEXP in dimension 2.</figcaption></figure>

We consider the following decomposition of the error $u(t, x) - u^m(t, x)$ for all $t \in [0, T]$ and $x \in \bar{Q}$: one has

$$u(t,x) - u^m(t,x) = A_1^m(t,x) + A_2^m(t,x) + A_3^m(t,x) + A_4^m(t,x) \\ + B_1^m(t,x) + B_2^m(t,x) + B_3^m(t,x) + B_4^m(t,x),$$

where the error terms $A_j^m(t,x)$, for $j \in \{1, \dots, 4\}$, are defined as

$$A_1^m(t,x) = \int_0^t \int_Q [G_d(t-s, x, y) - G_d(t - \kappa_m^T(s), x, y)] b(s, y, u(s, y)) \, dy \, ds,$$

$$A_2^m(t,x) = \int_0^t \int_Q G_d(t - \kappa_m^T(s), x, y) [b(s, y, u(s, y)) - b(\kappa_m^T(s), y, u(s, y))] \, dy \, ds,$$

$$A_3^m(t,x) = \int_0^t \int_Q G_d(t - \kappa_m^T(s), x, y) [b(\kappa_m^T(s), y, u(s, y)) - b(\kappa_m^T(s), y, u(\kappa_m^T(s), y))] \, dy \, ds,$$

$$A_4^m(t,x) = \int_0^t \int_Q G_d(t - \kappa_m^T(s), x, y) [b(\kappa_m^T(s), y, u(\kappa_m^T(s), y)) - b(\kappa_m^T(s), y, u^m(\kappa_m^T(s), y))] \, dy \, ds,$$

and the error terms $B_j^m(t,x)$, for $j \in \{1, \dots, 4\}$, are defined as

$$B_1^m(t,x) = \int_0^t \int_Q [G_d(t-s, x, y) - G_d(t - \kappa_m^T(s), x, y)] \sigma(s, y, u(s, y)) F_\alpha(ds, dy),$$

$$B_2^m(t,x) = \int_0^t \int_Q G_d(t - \kappa_m^T(s), x, y) [\sigma(s, y, u(s, y)) - \sigma(\kappa_m^T(s), y, u(s, y))] F_\alpha(ds, dy),$$

$$B_3^m(t,x) = \int_0^t \int_Q G_d(t - \kappa_m^T(s), x, y) [\sigma(\kappa_m^T(s), y, u(s, y)) - \sigma(\kappa_m^T(s), y, u(\kappa_m^T(s), y))] F_\alpha(ds, dy),$$

$$B_4^m(t,x) = \int_0^t \int_Q G_d(t - \kappa_m^T(s), x, y) [\sigma(\kappa_m^T(s), y, u(\kappa_m^T(s), y)) - \sigma(\kappa_m^T(s), y, u^m(\kappa_m^T(s), y))] F_\alpha(ds, dy).$$
For all $j \in \{1, \dots, 4\}$ and all $t \in [0, T]$, define

$$ \begin{aligned} \overline{\mathcal{A}}_{j,2p}^{m}(t,x) &= \|\overline{\mathcal{A}}_{j}^{m}(t,x)\|_{2p}, \quad \overline{\mathcal{B}}_{j,2p}^{m}(t,x) = \|\overline{\mathcal{B}}_{j}^{m}(t,x)\|_{2p}, & \forall x \in Q, \\ \overline{\mathcal{A}}_{j,2p}^{m}(t) &= \sup_{x \in \overline{Q}} \overline{\mathcal{A}}_{j,2p}^{m}(t,x), \quad \overline{\mathcal{B}}_{j,2p}^{m}(t) = \sup_{x \in \overline{Q}} \overline{\mathcal{B}}_{j,2p}^{m}(t,x). \end{aligned} $$

Applying the Minkowski inequality, one has

$$ \begin{aligned} \mathcal{E}_{2p}^{m}(t,x) &\le \sum_{j=1}^{4} \mathcal{A}_{j,2p}^{m}(t,x) + \sum_{j=1}^{4} \mathcal{B}_{j,2p}^{m}(t,x), & \forall (t,x) \in [0, T] \times Q, \\ \overline{\mathcal{E}}_{2p}^{m}(t) &\le \sum_{j=1}^{4} \overline{\mathcal{A}}_{j,2p}^{m}(t) + \sum_{j=1}^{4} \overline{\mathcal{B}}_{j,2p}^{m}(t), & \forall t \in [0, T]. \end{aligned} $$

The proof proceeds by proving upper bounds for each error term $\overline{\mathcal{A}}_{j,2p}^{m}(t)$ and $\overline{\mathcal{B}}_{j,2p}^{m}(t)$ with $j \in \{1, \dots, 4\}$.

In the sequel, the value of the parameter $\gamma \in (\frac{1}{2} - \frac{\alpha}{4})$ is fixed. Note that for all $s \in [0, T]$ one has $\kappa_m^T(s) \le s$ and $|s - \kappa_m^T(s)| \le \tau$.

• Treatment of the error term $\overline{\mathcal{A}}_{1,2p}^{m}(t)$.

Applying the Minkowski inequality, for all $(t, x) \in [0, T] \times \overline{Q}$, one has

$$ \begin{aligned} \|\mathcal{A}_1^m(t,x)\|_{2p} &\le \int_0^t \int_Q |G_d(t-s, x, y) - G_d(t - \kappa_m^T(s), x, y)| \|b(s, y, u(s, y))\|_{2p} dy ds \\ &\le \sup_{s \in [0, T]} \sup_{y \in \overline{Q}} \|b(s, y, u(s, y))\|_{2p} \int_0^t \int_Q |G_d(t-s, x, y) - G_d(t - \kappa_m^T(s), x, y)| dy ds. \end{aligned} $$

On the one hand, owing to the condition (LG) from Assumption 2, the mapping $b$ has at most linear growth. Applying the moment bounds (21) from Proposition 6 gives the upper bounds

$$ \sup_{s \in [0, T]} \sup_{y \in \overline{Q}} \|b(s, y, u(s, y))\|_{2p} \le C(1 + \sup_{s \in [0, T]} \sup_{y \in \overline{Q}} \|u(s, y)\|_{2p}) \le C_p(T)(1 + \|u_0\|_{\infty}). $$

On the other hand, applying the inequality (11) on the heat kernel (see Subsection 2.3), one obtains the upper bounds

$$ \begin{aligned} \int_0^t \int_Q |G_d(t-s, x, y) - G_d(t - \kappa_m^T(s), x, y)| \, dy \, ds &\le C \int_0^t \frac{|s - \kappa_m^T(s)|^\gamma}{(t-s)^\gamma} \int_Q (t-s)^{-\frac{d}{2}} e^{-c \frac{|y-x|^2}{2(t-s)}} \, dy \, ds \\ &\quad + C \int_0^t \frac{|s - \kappa_m^T(s)|^\gamma}{(t-s)^\gamma} \int_Q (t - \kappa_m^T(s))^{-\frac{d}{2}} e^{-c \frac{|y-x|^2}{2(t-\kappa_m^T(s))}} \, dy \, ds \\ &\le C \int_0^t \frac{\tau^\gamma}{(t-s)^\gamma} \, ds \\ &\le C_\gamma(T) \tau^\gamma. \end{aligned} $$

As a result, for the error term $\overline{\mathcal{A}}_{1,2p}^{m}(t)$, one obtains the following upper bound: there exists $C_{p,\gamma}(T) \in (0, \infty)$ such that one has

$$ (37) \qquad \sup_{t \in [0, T]} \overline{\mathcal{A}}_{1,2p}^{m}(t) \le C_{p,\gamma}(T)(1 + \|u_0\|_{\infty})\tau^{\gamma}. $$
• Treatment of the error term $\overline{\mathcal{A}}_{2,2p}^m(t)$.

Applying the Minkowski inequality, for all $(t, x) \in [0, T] \times \overline{Q}$, one has

$$ \|A_2^m(t, x)\|_{2p} \le \int_0^t \int_Q G_d(t - \kappa_m^T(s), x, y) \|b(s, y, u(s, y)) - b(\kappa_m^T(s), y, u(s, y))\|_{2p} dy ds. $$

Owing to the condition (L) from Assumption 2 and to the assumption $\gamma \in (0, \frac{1}{2} - \frac{\alpha}{4})$, the mapping $s \in [0, T] \mapsto b(s, y, u)$ satisfies a Hölder continuity property with exponent $\frac{1}{2} - \frac{\alpha}{4}$, uniformly with respect to $(y, u) \in \mathbb{R}^2$. Therefore, one has

$$ \begin{aligned} \|A_2^m(t, x)\|_{2p} &\le C \int_0^t \int_Q G_d(t - \kappa_m^T(s), x, y) |s - \kappa_m^T(s)|^{\frac{1}{2}-\frac{\alpha}{4}} dy ds \\ &\le C \tau^{\frac{1}{2}-\frac{\alpha}{4}} \int_0^t \int_Q G_d(t - \kappa_m^T(s), x, y) dy ds. \end{aligned} $$

Applying the property (6) of the heat kernel $G_d$ (see Subsection 2.3), one obtains

$$ \|A_2^m(t, x)\|_{2p} \le C T \tau^{\frac{1}{2}-\frac{\alpha}{4}}. $$

Recall that $\gamma \in (0, \frac{1}{2} - \frac{\alpha}{4})$. As a result, for the error term $\overline{\mathcal{A}}_{2,2p}^m(t)$, one obtains the following upper bound: there exists $C_{p,\gamma}(T) \in (0, \infty)$ such that one has

$$ (38) \qquad \sup_{t \in [0,T]} \overline{\mathcal{A}}_{2,2p}^m(t) \le C_{p,\gamma}(T) \tau^{\gamma}. $$

• Treatment of the error term $\overline{\mathcal{A}}_{3,2p}^m(t)$.

Applying the Minkowski inequality, for all $(t, x) \in [0, T] \times \overline{Q}$, one has

$$ \|A_3^m(t, x)\|_{2p} \le \int_0^t \int_Q G_d(t - \kappa_m^T(s), x, y) \|b(\kappa_m^T(s), y, u(s, y)) - b(\kappa_m^T(s), y, u(\kappa_m^T(s), y))\|_{2p} dy ds. $$

Owing to the condition (L) from Assumption 2, the mapping $u \mapsto b(s, y, u)$ satisfies a global Lipschitz continuity property. Therefore, applying the property (6) of the heat kernel $G_d$ (see Subsection 2.3), one has

$$ \begin{aligned} \|A_3^m(t, x)\|_{2p} &\le C \int_0^t \int_Q G_d(t - \kappa_m^T(s), x, y) \|u(s, y) - u(\kappa_m^T(s), y)\|_{2p} dy ds \\ &\le C \int_0^t \sup_{y \in \overline{Q}} \|u(s, y) - u(\kappa_m^T(s), y)\|_{2p} \int_Q G_d(t - \kappa_m^T(s), x, y) dy ds \\ &\le C \int_0^t \sup_{y \in \overline{Q}} \|u(s, y) - u(\kappa_m^T(s), y)\|_{2p} ds. \end{aligned} $$

To proceed, one needs to apply the regularity property (26) from Proposition 7 of the solution $u$. However, since $\kappa_m^T(s) = 0$ for all $s \in [0, \tau)$, it is not applicable on the time interval $[0, \min(t, \tau))$. The integral is then decomposed as

$$ \int_0^t \sup_{y \in \overline{Q}} \|u(s, y) - u(\kappa_m^T(s), y)\|_{2p} \mathrm{d}s = \int_0^{\min(t, \tau)} \sup_{y \in \overline{Q}} \|u(s, y) - u(\kappa_m^T(s), y)\|_{2p} \mathrm{d}s \\ + \mathbb{1}_{t>\tau} \int_{\tau}^t \sup_{y \in \overline{Q}} \|u(s, y) - u(\kappa_m^T(s), y)\|_{2p} \mathrm{d}s. $$
On the one hand, note that for all $s \in [0, \min(t, \tau))$ one has $u(\kappa_m^T(s), \cdot) = u_0$, thus applying the Minkowski inequality and the moment bounds (21) from Proposition 6, for the first integral one has

$$ \int_0^{\min(t, \tau)} \sup_{y \in \bar{Q}} \|u(s, y) - u(\kappa_m^T(s), y)\|_{2p} \mathrm{d}s \le \int_0^{\min(t, \tau)} \left( \sup_{s \in [0, T]} \sup_{y \in \bar{Q}} \|u(s, y)\|_{2p} + \|u_0\|_{\infty} \right) \mathrm{d}s \\ \le C\tau(1 + \|u_0\|_{\infty}). $$

On the other hand, applying the regularity property (26) from Proposition 7, for all $t \in (\tau, \infty)$ one has

$$ \int_{\tau}^{t} \sup_{y \in Q} \|u(s, y) - u(\kappa_m^T(s), y)\|_{2p} \mathrm{d}s \le C_{p,\gamma}(T)(1 + \|u_0\|_{\infty}) \int_{\tau}^{t} \frac{|s - \kappa_m^T(s)|^{\gamma}}{(\kappa_m^T(s))^{\gamma}} \mathrm{d}s. $$

Using the lower bound $\kappa_m^T(s) \ge s - \tau$, the upper bound $|s - \kappa_m^T(s)| \le \tau$ and a change of variables, one obtains the upper bounds

$$ \int_{\tau}^{t} \frac{|s - \kappa_{m}^{T}(s)|^{\gamma}}{(\kappa_{m}^{T}(s))^{\gamma}} \mathrm{d}s \le \int_{\tau}^{t} \frac{\tau^{\gamma}}{(s - \tau)^{\gamma}} \mathrm{d}s \le \tau^{\gamma} \int_{0}^{t-\tau} \frac{1}{s^{\gamma}} \mathrm{d}s \le C_{\gamma}(T)\tau^{\gamma}. $$

As a result, for the error term $\overline{A}_{3,2p}^m(t)$, one obtains the following upper bound: there exists $C_{p,\gamma}(T) \in (0, \infty)$ such that one has

$$ (39) \qquad \sup_{t \in [0, T]} \overline{A}_{3,2p}^{m}(t) \le C_{p,\gamma}(T)(1 + \|u_0\|_{\infty})\tau^{\gamma}. $$

• Treatment of the error term $\overline{A}_{4,2p}^m(t)$.

Applying the Minkowski inequality, for all $(t, x) \in [0, T] \times \bar{Q}$, one has

$$ \|A_4^m(t, x)\|_{2p} \le \int_0^t \int_Q G_d(t - \kappa_m^T(s), x, y) \|b(\kappa_m^T(s), y, u(\kappa_m^T(s), y)) - b(\kappa_m^T(s), y, u^m(\kappa_m^T(s), y))\|_{2p} \mathrm{d}y \mathrm{d}s. $$

Owing to the condition (L) from Assumption 2, the mapping $u \mapsto b(s, y, u)$ satisfies a global Lipschitz continuity property. Therefore, recalling the definition of $\overline{\mathcal{E}}_{2p}^m(s)$ and applying the property (6) of the heat kernel $G_d$ (see Subsection 2.3), one has

$$ \begin{aligned} \|A_4^m(t, x)\|_{2p} &\le C \int_0^t \int_Q G_d(t - \kappa_m^T(s), x, y) \|u(\kappa_m^T(s), y) - u^m(\kappa_m^T(s), y)\|_{2p} \mathrm{d}y \mathrm{d}s \\ &\le C \int_0^t \int_Q G_d(t - \kappa_m^T(s), x, y) \overline{\mathcal{E}}_{2p}^m(\kappa_m^T(s)) \mathrm{d}y \mathrm{d}s \\ &\le C(T) \int_0^t \overline{\mathcal{E}}_{2p}^m(\kappa_m^T(s)) \mathrm{d}s. \end{aligned} $$

As a result, for the error term $\overline{A}_{4,2p}^m(t)$, one obtains the following upper bound: there exists $C(T) \in (0, \infty)$ such that one has

$$ (40) \qquad \sup_{t \in [0, T]} \overline{A}_{4,2p}^{m}(t) \le C(T) \int_{0}^{t} \overline{\mathcal{E}}_{2p}^{m}(\kappa_{m}^{T}(s)) \mathrm{d}s. $$

This concludes the treatment of the error terms $\overline{A}_{j,2p}^m(t)$ with $j \in \{1, \dots, 4\}$. It remains to deal with the error terms $\overline{B}_{j,2p}^m(t)$ with $j \in \{1, \dots, 4\}$.

• Treatment of the error term $\overline{B}_{1,2p}^m(t)$.
Applying the Burkholder–Davis–Gundy inequality (4), for all $(t, x) \in [0, T] \times \overline{Q}$, one has

$$ \mathbb{E}[|B_1^m(t, x)|^{2p}] \le C_p(T) \mathbb{E} \left[ \left| \int_0^t \left\| [G_d(t-s, x, \cdot) - G_d(t - \kappa_m^T(s), x, \cdot)] \sigma(s, \cdot, u(s, \cdot)) \right\|_{(\alpha)}^2 ds \right|^p \right]. $$

For all $(t, x) \in [0, T] \times \overline{Q}$, introduce the auxiliary positive measure $\nu_{t,x}^m$ on $[0, t] \times Q^2$, which is absolutely continuous with respect to the Lebesgue measure, and with Radon–Nikodym derivative given by

$$ \frac{d\nu_{t,x}^m(s, y, z)}{ds dy dz} = |G_d(t-s, x, y) - G_d(t - \kappa_m^T(s), x, y)| |G_d(t-s, x, z) - G_d(t - \kappa_m^T(s), x, z)| |y-z|^{-\alpha}. $$

Making use of the auxiliary measure $\nu_{t,x}^m$ introduced above yields the upper bound

$$ \mathbb{E}[|B_1^m(t, x)|^{2p}] \le C_p(T) \mathbb{E} \left[ \left| \int_0^t \iint_{Q \times Q} |\sigma(s, y, u(s, y))| |\sigma(s, z, u(s, z))| d\nu_{t,x}^m(s, y, z) \right|^p \right]. $$

Applying the Hölder inequality, one then obtains

$$ \begin{aligned} \mathbb{E}[|B_1^m(t, x)|^{2p}] & \le C_p(T) \int_0^t \iint_{Q \times Q} \mathbb{E}[|\sigma(s, y, u(s, y))|^p |\sigma(s, z, u(s, z))|^p] d\nu_{t,x}^m(s, y, z) \nu_{t,x}^m([0, t] \times Q^2)^{p-1} \\ & \le C_p(T) \sup_{s \in [0, T]} \sup_{y, z \in \overline{Q}} \mathbb{E}[|\sigma(s, y, u(s, y))|^p |\sigma(s, z, u(s, z))|^p] \nu_{t,x}^m([0, t] \times Q^2)^p. \end{aligned} $$

Owing to the condition (LG) from Assumption 2, the mapping $\sigma$ has at most linear growth. Therefore, applying the moment bounds (21) from Proposition 6, one obtains the upper bounds

$$ \begin{aligned} \sup_{s \in [0, T]} \sup_{y, z \in \overline{Q}} \mathbb{E}[|\sigma(s, y, u(s, y))|^p |\sigma(s, z, u(s, z))|^p] & \le \sup_{s \in [0, T]} \sup_{y \in \overline{Q}} \mathbb{E}[|\sigma(s, y, u(s, y))|^{2p}] \\ & \le C_p(T) (1 + \|u_0\|_{\infty}^{2p}). \end{aligned} $$

To proceed, one needs to prove an upper bound for

$$ \nu_{t,x}^m([0, t] \times Q^2) = \int_0^t \|G_d(t-s, x, \cdot) - G_d(t - \kappa_m^T(s), x, \cdot)\|_{(\alpha)}^2 ds. $$

For all $s \in (0, t)$, let $\delta(s) = (t - \kappa_m^T(s)) - (t - s) = s - \kappa_m^T(s)$. In order to apply the inequality (10) on the heat kernel from Lemma 3 (see Subsection 2.3), the cases $\delta(s) < t - s$ and $\delta(s) \ge t - s$ are treated separately.

On the one hand, if $\delta(s) < t - s$, one obtains

$$ \begin{aligned} & \|G_d(t-s, x, \cdot) - G_d(t - \kappa_m^T(s), x, \cdot)\|_{(\alpha)}^2 \\ & \le C \frac{\delta(s)^2}{(t-s)^2} \iint_{Q \times Q} (t-s)^{-d} e^{-c \frac{|y-x|^2}{2(t-s)}} e^{-c \frac{|z-x|^2}{2(t-s)}} |y-z|^{-\alpha} dy dz. \end{aligned} $$

Under the condition $\delta(s) < t - s$, for all $\gamma \in (0, 1)$ one has

$$ \frac{\delta(s)^2}{(t-s)^2} \le \frac{\delta(s)^{2\gamma}}{(t-s)^{2\gamma}}. $$
Applying changes of variable $y' = y-x$ and $z' = z-x$, and $y'' = y'/\sqrt{t-s}$ and $z'' = z'/\sqrt{t-s}$, one has

$$ 
\begin{aligned}
& \iint_{Q \times Q} (t-s)^{-d} e^{-c \frac{|y-x|^2}{2(t-s)}} e^{-c \frac{|z-x|^2}{2(t-s)}} |y-z|^{-\alpha} dy dz \\
& \le \iint_{\mathbb{R}^d \times \mathbb{R}^d} (t-s)^{-d} e^{-c \frac{|y'|^2}{2(t-s)}} e^{-c \frac{|z'|^2}{2(t-s)}} |y'-z'|^{-\alpha} dy' dz' \\
& \le C(t-s)^{-\frac{\alpha}{2}} \iint_{\mathbb{R}^d \times \mathbb{R}^d} e^{-c \frac{|y''|^2}{2}} e^{-c \frac{|z''|^2}{2}} |y''-z''|^{-\alpha} dy'' dz'' \\
& \le C_{\alpha}(t-s)^{-\frac{\alpha}{2}}.
\end{aligned}
$$

On the other hand, if $\delta(s) = s - \kappa_m^T(s) \ge t-s$, applying the inequality (12) from Lemma 4, one obtains

$$ 
\begin{aligned}
\|G_d(t-s, x, \cdot) - G_d(t - \kappa_m^T(s), x, \cdot)\|_{(\alpha)}^2 &\le \|G_d(t-s, x, \cdot)\|_{(\alpha)}^2 + \|G_d(t - \kappa_m^T(s), x, \cdot)\|_{(\alpha)}^2 \\
&\le C_{\alpha}(t-s)^{-\frac{\alpha}{2}} + C_{\alpha}(t - \kappa_m^T(s))^{-\frac{\alpha}{2}} \\
&\le C_{\alpha}(t-s)^{-\frac{\alpha}{2}}.
\end{aligned}
$$

Moreover, under the condition $\delta(s) \ge t-s$, for all $\gamma \in (0, 1)$ one has

$$ 1 \le \frac{\delta(s)^{2\gamma}}{(t-s)^{2\gamma}}. $$

As a result, in both cases one obtains the upper bound

$$ \|G_d(t-s, x, \cdot) - G_d(t - \kappa_m^T(s), x, \cdot)\|_{(\alpha)}^2 \le C_{\alpha} \frac{(s - \kappa_m^T(s))^{2\gamma}}{(t-s)^{2\gamma + \frac{\alpha}{2}}}. $$

Therefore, taking into account the condition $\gamma \in (0, \frac{1}{2} - \frac{\alpha}{4})$, for all $(t, x) \in [0, T] \times Q$, one has

$$ 
\begin{aligned}
\nu_{t,x}^m([0,t] \times Q^2) &= \int_0^t \|G_d(t-s, x, \cdot) - G_d(t - \kappa_m^T(s), x, \cdot)\|_{(\alpha)}^2 ds \\
&\le C_{\alpha} \tau^{2\gamma} \int_0^t (t-s)^{-2\gamma - \frac{\alpha}{2}} ds \le C_{\alpha, \gamma}(T) \tau^{2\gamma}.
\end{aligned}
$$

Finally, for all $(t, x) \in [0, T] \times Q$, one has the upper bound

$$ \|B_1^m(t, x)\|_{2p}^{2p} = \mathbb{E}[|B_1^m(t, x)|^{2p}] \le C_{\alpha, \gamma}(T)(1 + \|u_0\|_{\infty}^{2p})\tau^{2p\gamma}. $$

As a result, for the error term $\overline{B}_{1,2p}^m(t)$, one obtains the following upper bound: there exists $C_{p,\alpha,\gamma}(T) \in (0, \infty)$ such that one has

$$ (41) \qquad \sup_{t \in [0, T]} \overline{B}_{1,2p}^{m}(t) \le C_{p,\alpha,\gamma}(T) (1 + \|u_0\|_{\infty}) \tau^{\gamma}. $$

• Treatment of the error term $\overline{B}_{2,2p}^m(t)$.
Applying the inequality (18) from Lemma 5 (with $\kappa(s) = \kappa_m^T(s)$), for all $(t, x) \in [0, T] \times \overline{Q}$, one has

$$ \mathbb{E}[|B_2^m(t, x)|^{2p}] \le C_{p,\alpha}(T) \sup_{0 \le s \le T} \sup_{y \in \overline{Q}} \mathbb{E}[|\sigma(s, y, u(s, y)) - \sigma(\kappa_m^T(s), y, u(s, y))|^{2p}]. $$

Owing to the condition (L) from Assumption 2, the mapping $s \in [0, T] \mapsto b(s, y, u)$ satisfies a Hölder continuity property with exponent $\frac{1}{2} - \frac{\alpha}{4}$, uniformly with respect to $(y, u) \in \mathbb{R}^2$. Therefore, recalling that $\gamma \in (0, \frac{1}{2} - \frac{\alpha}{4})$, one has

$$ \mathbb{E}[|B_2^m(t, x)|^{2p}] \le C_{p,\alpha}(T) \tau^{2p\gamma}. $$

As a result, for the error term $\overline{B}_{2,2p}^m(t)$, one obtains the following upper bound: there exists $C_{p,\alpha,\gamma}(T) \in (0, \infty)$ such that one has

$$ (42) \qquad \sup_{t \in [0, T]} \overline{B}_{2,2p}^m(t) \le C_{p,\alpha,\gamma}(T) \tau^{\gamma}. $$

• Treatment of the error term $\overline{B}_{3,2p}^m(t)$.

Applying the inequality (16) from Lemma 5 (with $\kappa(s) = \kappa_m^T(s)$), for all $(t, x) \in [0, T] \times \overline{Q}$, one has

$$ \|B_3^m(t, x)\|_{2p}^2 \le C_{p,\alpha}(T) \int_0^t (t - \kappa_m^T(s))^{-\frac{\alpha}{2}} \sup_{y \in \overline{Q}} \|\sigma(\kappa_m^T(s), y, u(s, y)) - \sigma(\kappa_m^T(s), y, u(\kappa_m^T(s), y))\|_{2p}^2 ds. $$

Owing to the condition (L) from Assumption 2, the mapping $\sigma$ satisfies a global Lipschitz continuity property. Therefore, one has

$$ \|B_3^m(t, x)\|_{2p}^2 \le C_{p,\alpha}(T) \int_0^t (t - \kappa_m^T(s))^{-\frac{\alpha}{2}} \sup_{y \in \overline{Q}} \|u(s, y) - u(\kappa_m^T(s), y)\|_{2p}^2 ds. $$

Like in the proof of the upper bound (39) for the error term $\overline{B}_{3,2p}^m(t)$, the regularity property (26) from Proposition 7 is not directly applicable, and it is necessary to decompose the integral as

$$ \begin{aligned} & \int_0^t (t - \kappa_m^T(s))^{-\frac{\alpha}{2}} \sup_{y \in \overline{Q}} \|u(s, y) - u(\kappa_m^T(s), y)\|_{2p}^2 ds \\ & = \int_0^{\min(t, \tau)} (t - \kappa_m^T(s))^{-\frac{\alpha}{2}} \sup_{y \in \overline{Q}} \|u(s, y) - u(\kappa_m^T(s), y)\|_{2p}^2 ds \\ & \quad + 1_{t>\tau} \int_{\tau}^t (t - \kappa_m^T(s))^{-\frac{\alpha}{2}} \sup_{y \in \overline{Q}} \|u(s, y) - u(\kappa_m^T(s), y)\|_{2p}^2 ds. \end{aligned} $$

On the one hand, note that for all $s \in [0, \min(t, \tau))$ one has $u(\kappa_m^T(s), \cdot) = u_0$, thus applying the Minkowski inequality and the moment bounds (21) from Proposition 6, for the first integral
one has

$$
\begin{align*}
& \int_0^{\min(t, \tau)} (t - \kappa_m^T(s))^{-\frac{\alpha}{2}} \sup_{y \in \bar{Q}} \|u(s, y) - u(\kappa_m^T(s), y)\|_{2p}^2 \, ds \\
& \le \int_0^{\min(t, \tau)} (t - \kappa_m^T(s))^{-\frac{\alpha}{2}} \left( \sup_{y \in \bar{Q}} \|u(s, y)\|_{2p}^2 + \|u_0\|_{\infty}^2 \right) \, ds \\
& \le \int_0^{\min(t, \tau)} (t - \kappa_m^T(s))^{-\frac{\alpha}{2}} \, ds (1 + \|u_0\|_{\infty})^2.
\end{align*}
$$

To deal with the integral in the right-hand side above, observe that since $\kappa_m^T(s) = 0$ for all $s \in [0, \min(t, \tau)]$, one has

$$
\int_0^{\min(t, \tau)} (t - \kappa_m^T(s))^{-\frac{\alpha}{2}} \, ds = \min(t, \tau) t^{-\frac{\alpha}{2}} = \min(t, \tau)^{\frac{\alpha}{2}} \min(t, \tau)^{1-\frac{\alpha}{2}} t^{-\frac{\alpha}{2}} \le \tau^{1-\frac{\alpha}{2}} \le C \tau^{2\gamma}.
$$

On the other hand, applying the regularity property (26) from Proposition 7, for all $t \in (\tau, \infty)$, one has

$$
\begin{align*}
& \int_{\tau}^{t} (t - \kappa_{m}^{T}(s))^{-\frac{\alpha}{2}} \sup_{y \in \bar{Q}} \|u(s, y) - u(\kappa_{m}^{T}(s), y)\|_{2p}^{2} \, ds \\
& \le C_{p,\gamma}(T) \int_{\tau}^{t} (t - \kappa_{m}^{T}(s))^{-\frac{\alpha}{2}} \frac{|s - \kappa_{m}^{T}(s)|^{2\gamma}}{\kappa_{m}^{T}(s)^{2\gamma}} \, ds (1 + \|u_0\|_{\infty})^{2}.
\end{align*}
$$

The integral in the right-hand side above can be treated as follows: using the inequalities
$s \ge \kappa_m^T(s)$ and $\kappa_m^T(s) \ge s - \tau$, for all $t > \tau$ one has

$$
\begin{align*}
\int_{\tau}^{t} (t - \kappa_{m}^{T}(s))^{-\frac{\alpha}{2}} \frac{|s - \kappa_{m}^{T}(s)|^{2\gamma}}{\kappa_{m}^{T}(s)^{2\gamma}} \, ds &\le \tau^{2\gamma} \int_{\tau}^{t} (t - \kappa_{m}^{T}(s))^{-\frac{\alpha}{2}} \kappa_{m}^{T}(s)^{-2\gamma} \, ds \\
&\le \tau^{2\gamma} \int_{\tau}^{t} (t-s)^{-\frac{\alpha}{2}} (s-\tau)^{-2\gamma} \, ds \\
&\le \tau^{2\gamma} \int_{0}^{t-\tau} (t-\tau-s)^{-\frac{\alpha}{2}} s^{-2\gamma} \, ds \\
&\le \tau^{2\gamma} (t-\tau)^{1-\frac{\alpha}{2}-2\gamma} \int_{0}^{1} (1-s')^{-\frac{\alpha}{2}} (s')^{-2\gamma} \, ds' \\
&\le C_{\alpha,\gamma}(T) \tau^{2\gamma}.
\end{align*}
$$

As a result, for the error term $\overline{B}_{3,2p}^m(t)$, one obtains the following upper bound: there exists $C_{p,\alpha,\gamma}(T) \in (0, \infty)$ such that one has

$$
(43) \qquad \sup_{t \in [0, T]} \overline{B}_{3,2p}^{m}(t) \le C_{p,\alpha,\gamma}(T) (1 + \|u_0\|_{\infty}) \tau^{\gamma}.
$$

• Treatment of the error term $\overline{B}_{4,2p}^m(t)$.

Applying the inequality (17) from Lemma 5 (with $\kappa(s) = \kappa_m^T(s) \le s$), for all $(t, x) \in [0, T] \times \bar{Q}$, one has

$$
\mathbb{E}[|B_4^m(t, x)|^{2p}] \le C_{p,\alpha}(T) \int_0^t (t-s)^{-\frac{\alpha}{2}} \sup_{y \in \bar{Q}} \mathbb{E}[|\sigma(\kappa_m^T(s), y, u(\kappa_m^T(s), y)) - \sigma(\kappa_m^T(s), y, u^m(\kappa_m^T(s), y))|^{2p}] ds.
$$
Owing to the condition (L) from Assumption 2, the mapping $\sigma$ satisfies a global Lipschitz continuity property. Therefore, recalling the definition of $\overline{\mathcal{E}}_{2p}^m(s)$, one has

$$ \begin{aligned} \mathbb{E}[|B_4^m(t, x)|^{2p}] &\le C_{p,\alpha}(T) \int_0^t (t-s)^{-\frac{\alpha}{2}} \sup_{y \in \bar{Q}} \mathbb{E}[|u(\kappa_m^T(s), y) - u^m(\kappa_m^T(s), y)|^{2p}] ds \\ &\le C_{p,\alpha}(T) \int_0^t (t-s)^{-\frac{\alpha}{2}} \overline{\mathcal{E}}_{2p}^m(\kappa_m^T(s))^{2p} ds. \end{aligned} $$

As a result, for the error term $\overline{B}_{4,2p}^m(t)$, one obtains the following upper bound: there exists $C(T) \in (0, \infty)$ such that one has

$$ (44) \qquad \sup_{t \in [0,T]} \overline{B}_{4,2p}^m(t)^{2p} \le C(T) \int_0^t (t-s)^{-\frac{\alpha}{2}} \overline{\mathcal{E}}_{2p}^m(\kappa_m^T(s))^{2p} ds. $$

This concludes the treatment of the error terms $\overline{B}_{j,2p}^m(t)$ for $j \in \{1, \dots, 4\}$.

### Conclusion.

Gathering the upper bounds (37), (38), (39) and (40) for the error terms $\overline{A}_{j,2p}^m(t)$ for $j \in \{1, \dots, 4\}$, and the upper bounds (41), (42), (43) and (44) for the error terms $\overline{B}_{j,2p}^m(t)$ for $j \in \{1, \dots, 4\}$, one obtains the following inequality: for all $T \in (0, \infty)$, $\alpha \in (0, 2 \wedge d)$, $\gamma \in (\frac{1}{2} - \frac{\alpha}{4})$ and $p \in [1, \infty)$, there exists $C_{p,\alpha,\gamma}(T) \in (0, \infty)$ such that for all $t \in [0, T]$ one has

$$ (45) \qquad \overline{\mathcal{E}}_{2p}^m(t)^{2p} \le C_{p,\alpha,\gamma}(T) \tau^{2p\gamma} (1 + \|u_0\|_{\infty})^{2p} + C_{p,\alpha,\gamma}(T) \int_0^t (1 + (t-s)^{-\frac{\alpha}{2}}) \overline{\mathcal{E}}_{2p}^m(\kappa_m^T(s))^{2p} ds. $$

Applying a version of the Grönwall inequality (c.f [25, Lem. 15]) yields the upper bound

$$ (46) \qquad \sup_{t \in [0,T]} \overline{\mathcal{E}}_{2p}^m(t)^{2p} \le C_{p,\alpha,\gamma}(T) \tau^{2p\gamma} (1 + \|u_0\|_{\infty})^{2p}. $$

This provides the strong error estimates (35) and the proof of Theorem 8 is completed. $\square$

<figure><img src="image_10.png" /></figure>

## APPENDIX A. PROOF OF THE PROPERTIES OF THE HEAT KERNEL

The objective of this section is to provide the proofs of the properties of the heat kernel stated in Subsection 2.3.

Let us first provide the proof of Lemma 3.

*Proof of Lemma 3.* Let $t_1 < t_2$. The proof of the inequality (10) proceeds by considering two cases.

First, assume that $\frac{t_2-t_1}{t_1} < 1$, i. e. that $t_1 < t_2 < 2t_1$. Applying the inequality (9) on the temporal derivative of the heat kernel $G_d$, one obtains

$$ |G_d(t_2, x, y) - G_d(t_1, x, y)| \le \int_{t_1}^{t_2} |\partial_t G_d(t, x, y)| dt \le C_d(t_2 - t_1) t_1^{-\frac{d}{2}-1} e^{-c_d \frac{|y-x|^2}{2t_2}}. $$

Under the condition $t_2 < 2t_1$, one obtains the upper bound

$$ |G_d(t_2, x, y) - G_d(t_1, x, y)| \le C_d \frac{(t_2 - t_1)}{t_1} t_1^{-\frac{d}{2}} e^{-c_d \frac{|y-x|^2}{2t_1}}. $$
Secondly, assume that $\frac{t_2-t_1}{t_1} \ge 1$, i. e. that $t_2 \ge 2t_1$. Applying the triangle inequality and the inequality (7) from Lemma 2, one obtains for all $x, y \in \overline{Q}$

$$|G_d(t_2, x, y) - G_d(t_1, x, y)| \le G_d(t_2, x, y) + G_d(t_1, x, y) \\ \le C_d(t_1^{-\frac{d}{2}}e^{-c_d \frac{|y-x|^2}{2t_1}} + t_2^{-\frac{d}{2}}e^{-c_d \frac{|y-x|^2}{2t_2}}).$$

Combining the results above then provides the inequality (10).

The inequality (11) is a straightforward consequence of the inequality (10). Indeed, one first obtains the inequality

$$|G_d(t_2, x, y) - G_d(t_1, x, y)| \le C_d \min \left( \frac{t_2 - t_1}{t_1}, 1 \right) \left( t_1^{-\frac{d}{2}} e^{-c_d \frac{|y-x|^2}{2t_1}} + t_2^{-\frac{d}{2}} e^{-c_d \frac{|y-x|^2}{2t_2}} \right)$$

and if $\gamma \in (0, 1)$ for all $x \in \mathbb{R}$ one has $\min(x, 1) \le x^\gamma$.

This concludes the proof of Lemma 3.

<figure><img src="image_11.png" /></figure>

Let us now provide the proof of Lemma 4. The inequalities (12), (13), (14) and (15) are established successively. Recall that $d \in \mathbb{N}$ and $\alpha \in (0, 2 \wedge d)$.

*Proof of the inequality* (12). Let $t \in (0, \infty)$ and $x \in \overline{Q}$. Applying the inequality (7), one has

$$\|G_d(t, x, \cdot)\|_{(\alpha)}^2 = \iint_{Q \times Q} G_d(t, x, y) G_d(t, x, z) |y-z|^{-\alpha} dy dz \\ \le Ct^{-d} \iint_{Q \times Q} e^{-c \frac{|y-x|^2 + |z-x|^2}{t}} |y-z|^{-\alpha} dy dz \\ \le G_{t,x}^1 + G_{t,x}^2,$$

where $G_{t,x}^1$ and $G_{t,x}^2$ are defined as

$$G_{t,x}^1 = Ct^{-d} \iint_{Q \times Q} 1_{|x-y| \ge |y-z|} e^{-c \frac{|y-x|^2 + |z-x|^2}{t}} |y-z|^{-\alpha} dy dz \\ G_{t,x}^2 = Ct^{-d} \iint_{Q \times Q} 1_{|x-y| \le |y-z|} e^{-c \frac{|y-x|^2 + |z-x|^2}{t}} |y-z|^{-\alpha} dy dz.$$

Concerning $G_{t,x}^1$, taking into account the condition $|x-y| \ge |y-z|$ and applying the change of variables $(y', z') = (y-z, z-x)$, one obtains

$$G_{t,x}^1 = Ct^{-d} \iint_{Q \times Q} 1_{|x-y| \ge |y-z|} e^{-c \frac{|y-x|^2}{t}} e^{-c \frac{|z-x|^2}{t}} |y-z|^{-\alpha} dy dz \\ \le Ct^{-d} \iint_{Q \times Q} e^{-c \frac{|y-z|^2}{t}} e^{-c \frac{|z-x|^2}{t}} |y-z|^{-\alpha} dy dz \\ \le Ct^{-d} \int_{\mathbb{R}^d} e^{-c \frac{|y'|^2}{t}} |y'|^{-\alpha} dy' \int_{\mathbb{R}^d} e^{-c \frac{|z'|^2}{t}} dz'$$
Concerning $\mathcal{G}_{t,x}^2$, taking into account the condition $|x-y| \le |y-z|$ and applying the change of variables $(y', z') = (y - x, z - x)$, one obtains

$$
\begin{align*}
\mathcal{G}_{t,x}^2 &= Ct^{-d} \iint_{Q \times Q} \mathbf{1}_{|x-y| \le |y-z|} e^{-c \frac{|y-x|}{t}^2} e^{-c \frac{|z-x|}{t}^2} |y-z|^{-\alpha} dy dz \\
&\le Ct^{-d} \iint_{Q \times Q} e^{-c \frac{|y-x|}{t}^2} e^{-c \frac{|z-x|}{t}^2} |y-x|^{-\alpha} dy dz \\
&\le Ct^{-d} \int_{\mathbb{R}^d} e^{-c \frac{|y'|}{t}^2} |y'|^{-\alpha} dy' \int_{\mathbb{R}^d} e^{-c \frac{|z'|}{t}^2} dz'.
\end{align*}
$$

Note that one obtains the same upper bound for the two terms $\mathcal{G}_{t,x}^1$ and $\mathcal{G}_{t,x}^2$. Applying the supplementary change of variables $(y'', z'') = t^{-\frac{1}{2}}(y', z')$, one then has

$$
\begin{align*}
\mathcal{G}_{t,x}^1 + \mathcal{G}_{t,x}^2 &\le Ct^{-d} \int_{\mathbb{R}^d} e^{-c \frac{|y'|^2}{t}} |y'|^{-\alpha} dy' \int_{\mathbb{R}^d} e^{-c \frac{|z'|^2}{t}} dz' \\
&= Ct^{-\frac{\alpha}{2}} \int_{\mathbb{R}^d} |y''|^{-\alpha} e^{-c|y''|^2} dy'' \int_{\mathbb{R}^d} e^{-c|z''|^2} dz'' \\
&\le C_\alpha t^{-\frac{\alpha}{2}}.
\end{align*}
$$

One thus obtains the following upper bound: there exists $C_\alpha \in (0, \infty)$ such that for all $t \in (0, \infty)$ one has

$$
(47) \qquad t^{-d} \iint_{Q \times Q} e^{-c \frac{|y-x|^2 + |z-x|^2}{t}} |y-z|^{-\alpha} dy dz \le C_{\alpha} t^{-\frac{\alpha}{2}}.
$$

As a result, one has

$$
\|G_d(t, x, \cdot)\|_{(\alpha)}^2 \le C_\alpha t^{-\frac{\alpha}{2}}.
$$

The proof of the inequality (12) is thus completed.

<figure><img src="image_12.png" /></figure>

Proof of the inequality (13). Let $t_2 \ge t_1 \ge 0$ and $x \in \bar{Q}$. Applying the change of variable $s = t_2 - t$ and the inequality (12), one has

$$
\int_{t_1}^{t_2} \|G_d(t_2 - t, x, \cdot)\|_{(\alpha)}^2 dt = \int_{0}^{t_2-t_1} \|G_d(s, x, \cdot)\|_{(\alpha)}^2 ds \le C_{\alpha} \int_{0}^{t_2-t_1} s^{-\frac{\alpha}{2}} ds \le C_{\alpha} |t_2 - t_1|^{1-\frac{\alpha}{2}}.
$$

The proof of the inequality (13) is thus completed.

<figure><img src="image_13.png" /></figure>

Proof of the inequality (14). Let $x_1, x_2 \in \bar{Q}$, and set $\delta x = x_2 - x_1$. One has the decomposition

$$
\int_{0}^{+\infty} \|G_d(t, x_2, \cdot) - G_d(t, x_1, \cdot)\|_{(\alpha)}^2 dt = \delta \mathcal{G}^1 + \delta \mathcal{G}^2,
$$

where $\delta \mathcal{G}^1$ and $\delta \mathcal{G}^2$ are defined as

$$
\delta \mathcal{G}^1 = \int_0^{|\delta x|^2} \|G_d(t, x_2, \cdot) - G_d(t, x_1, \cdot)\|_{(\alpha)}^2 dt
$$

$$
\delta \mathcal{G}^2 = \int_{|\delta x|^2}^{\infty} \|G_d(t, x_2, \cdot) - G_d(t, x_1, \cdot)\|_{(\alpha)}^2 dt.
$$
Let us first deal with $\delta\mathcal{G}^1$. One has

$$ \begin{aligned} \int_0^{|\delta x|^2} \|G_d(t, x_2, \cdot) - G_d(t, x_1, \cdot)\|_{(\alpha)}^2 dt &\le 2 \int_0^{|\delta x|^2} \|G_d(t, x_1, \cdot)\|_{(\alpha)}^2 dt + 2 \int_0^{|\delta x|^2} \|G_d(t, x_2, \cdot)\|_{(\alpha)}^2 dt \\ &\le 4 \int_0^{|\delta x|^2} \sup_{x \in Q} \|G_d(t, x, \cdot)\|_{(\alpha)}^2 dt. \end{aligned} $$

Applying the inequality (12), one then gets

$$ \int_0^{|\delta x|^2} \sup_{x \in Q} \|G_d(t, x, \cdot)\|_{(\alpha)}^2 dt \le C_\alpha \int_0^{|\delta x|^2} t^{-\frac{\alpha}{2}} dt = C_\alpha (|\delta x|^2)^{1-\frac{\alpha}{2}} = C_\alpha |\delta x|^{2-\alpha}. $$

As a result, one obtains the upper bound

$$ (48) \qquad \delta\mathcal{G}^1 \le C_\alpha |\delta x|^{2-\alpha}. $$

Let us now deal with $\delta\mathcal{G}^2$. By the definition of the norm $\|\cdot\|_{(\alpha)}$, one has

$$ \begin{aligned} \delta\mathcal{G}^2 &= \int_{|\delta x|^2}^{\infty} \|G_d(t, x_2, \cdot) - G_d(t, x_1, \cdot)\|_{(\alpha)}^2 dt \\ &= \int_{|\delta x|^2}^{\infty} \iint_{Q \times Q} |G_d(t, x_2, y) - G_d(t, x_1, y)| |G_d(t, x_2, z) - G_d(t, x_1, z)| |y-z|^{-\alpha} dy dz dt. \end{aligned} $$

Let $x(\tau) = (1 - \tau)x_1 + \tau x_2$ for all $\tau \in [0, 1]$. For all $t \ge 0$ and all $y, z \in \overline{Q}$, one obtains

$$ \begin{aligned} |G_d(t, x_2, y) - G_d(t, x_1, y)| &= |G_d(t, x(1), y) - G_d(t, x(0), y)| = \left| \int_0^1 \nabla_x G_d(t, x(\theta), y) \cdot \delta x \, d\theta \right| \\ &\le |\delta x| \int_0^1 \left| \nabla_x G_d(t, x(\theta), y) \right| \, d\theta, \\ |G_d(t, x_2, z) - G_d(t, x_1, z)| &= |G_d(t, x(1), z) - G_d(t, x(0), z)| = \left| \int_0^1 \nabla_x G_d(t, x(\eta), z) \cdot \delta x \, d\eta \right| \\ &\le |\delta x| \int_0^1 \left| \nabla_x G_d(t, x(\eta), z) \right| \, d\eta. \end{aligned} $$

As a result, applying the inequality (8), one has the estimates

$$ \begin{aligned} \delta\mathcal{G}^2 &\le |\delta x|^2 \int_{|\delta x|^2}^{\infty} \iint_{Q \times Q} \int_{[0,1]^2} |\nabla_x G_d(t, x(\theta), y)| |\nabla_x G_d(t, x(\eta), y)| |y-z|^{-\alpha} dy dz d\theta d\eta dt \\ &\le C |\delta x|^2 \int_{|\delta x|^2}^{\infty} t^{-(d+1)} \iint_{[0,1]^2} \iint_{Q \times Q} e^{-c \frac{|y-x(\theta)|^2 + |z-x(\eta)|^2}{t}} |y-z|^{-\alpha} dy dz d\theta d\eta dt \\ &\le \delta\mathcal{G}^{2,1} + \delta\mathcal{G}^{2,2}, \end{aligned} $$
where $\delta\mathcal{G}^{2,1}$ and $\delta\mathcal{G}^{2,2}$ are defined as

$$ \delta\mathcal{G}^{2,1} = C|\delta x|^2 \int_{|\delta x|^2}^{\infty} t^{-(d+1)} \iint_{[0,1]^2 Q \times Q} \mathbb{1}_{|y-x(\theta)| \ge |y-z|} e^{-c \frac{|y-x(\theta)|^2 + |z-x(\eta)|^2}{t}} |y-z|^{-\alpha} dy dz d\theta d\eta dt, $$

$$ \delta\mathcal{G}^{2,2} = C|\delta x|^2 \int_{|\delta x|^2}^{\infty} t^{-(d+1)} \iint_{[0,1]^2 Q \times Q} \mathbb{1}_{|y-x(\theta)| \le |y-z|} e^{-c \frac{|y-x(\theta)|^2 + |z-x(\eta)|^2}{t}} |y-z|^{-\alpha} dy dz d\theta d\eta dt. $$

Concerning $\delta\mathcal{G}^{2,1}$, taking into account the condition $|y - x(\theta)| \ge |y - z|$ and applying the change of variable $(y', z') = (y - z, z - x(\eta))$, one obtains

$$ \begin{aligned} \delta\mathcal{G}^{2,1} &\le C|\delta x|^2 \int_{|\delta x|^2}^{\infty} t^{-(d+1)} \iint_{[0,1]^2 Q \times Q} e^{-c \frac{|y-z|^2}{t}} e^{-c \frac{|z-x(\eta)|^2}{t}} |y-z|^{-\alpha} dy dz d\theta d\eta dt \\ &\le C|\delta x|^2 \int_{|\delta x|^2}^{\infty} t^{-(d+1)} \iint_{[0,1]^2 \mathbb{R}^d \times \mathbb{R}^d} e^{-c \frac{|y'|^2}{t}} e^{-c \frac{|z'|^2}{t}} |y'|^{-\alpha} dy' dz' d\theta d\eta dt \\ &\le C|\delta x|^2 \int_{|\delta x|^2}^{\infty} t^{-(d+1)} \iint_{\mathbb{R}^d \times \mathbb{R}^d} e^{-c \frac{|y'|^2}{t}} e^{-c \frac{|z'|^2}{t}} |y'|^{-\alpha} dy' dz' dt. \end{aligned} $$

Concerning $\delta\mathcal{G}^{2,2}$, taking into account the condition $|y - x(\theta)| \le |y - z|$ and applying the change of variable $(y', z') = (y - x(\theta), z - x(\eta))$, one obtains

$$ \begin{aligned} \delta\mathcal{G}^{2,2} &\le C|\delta x|^2 \int_{|\delta x|^2}^{\infty} t^{-(d+1)} \iint_{[0,1]^2 Q \times Q} \mathbb{1}_{|y-x(\theta)| \le |y-z|} e^{-c \frac{|y-x(\theta)|^2}{t}} e^{-c \frac{|z-x(\eta)|^2}{t}} |y-x(\theta)|^{-\alpha} dy dz d\theta d\eta dt \\ &\le C|\delta x|^2 \int_{|\delta x|^2}^{\infty} t^{-(d+1)} \iint_{[0,1]^2 \mathbb{R}^d \times \mathbb{R}^d} e^{-c \frac{|y'|^2}{t}} e^{-c \frac{|z'|^2}{t}} |y'|^{-\alpha} dy' dz' d\theta d\eta dt \\ &\le C|\delta x|^2 \int_{|\delta x|^2}^{\infty} t^{-(d+1)} \iint_{\mathbb{R}^d \times \mathbb{R}^d} e^{-c \frac{|y'|^2}{t}} e^{-c \frac{|z'|^2}{t}} |y'|^{-\alpha} dy' dz' dt. \end{aligned} $$

Observe that one gets the same upper bound for the two terms $\delta\mathcal{G}^{2,1}$ and $\delta\mathcal{G}^{2,2}$. Applying the supplementary change of variables $(y'', z'') = t^{-\frac{1}{2}}(y', z')$, one then obtains the upper bounds

$$ \begin{aligned} \delta\mathcal{G}^2 &\le \delta\mathcal{G}^{2,1} + \delta\mathcal{G}^{2,2} \\ &\le C|\delta x|^2 \int_{|\delta x|^2}^{\infty} t^{-1-\frac{\alpha}{2}} dt \int_{\mathbb{R}^d} |y''|^{-\alpha} e^{-c|y''|^2} dy'' \int_{\mathbb{R}^d} e^{-c|z''|^2} dz'' \\ &\le C_{\alpha}|\delta x|^2 (|\delta x|^2)^{-\frac{\alpha}{2}} = C_{\alpha}|\delta x|^{2-\alpha}. \end{aligned} $$

As a result, one gets the upper bound

$$ (49) \qquad \delta\mathcal{G}^2 \le C_{\alpha}|\delta x|^{2-\alpha}. $$

Combining the upper bounds (48) and (49) yields the inequality

$$ \int_{0}^{+\infty} \|G_d(t, x_2, \cdot) - G_d(t, x_1, \cdot)\|_{(\alpha)}^2 dt = \delta\mathcal{G}^1 + \delta\mathcal{G}^2 \le C_{\alpha}|\delta x|^{2-\alpha}. $$

The proof of the inequality (14) is thus completed. $\square$

<figure><img src="image_14.png" /></figure>
Proof of the inequality (15). Let $t_2 \ge t_1 \ge 0$ and $x \in \overline{Q}$.

Owing to the definition of the norm $\| \cdot \|_{(\alpha)}$, one has

$$ \int_0^{t_1} \|G_d(t_2 - t, x, \cdot) - G_d(t_1 - t, x, \cdot)\|_{(\alpha)}^2 dt \\ = \int_0^{t_1} \iint_{Q \times Q} |G_d(t_2 - t, x, y) - G_d(t_1 - t, x, y)| |G_d(t_2 - t, x, z) - G_d(t_1 - t, x, z)| |y - z|^{-\alpha} dy dz dt. $$

Let $\gamma \in (0, 1)$. Applying the inequality (11) on the heat kernel $G_d$ (see Subsection 2.3), there exists $C_\gamma \in (0, \infty)$ such that for all $y \in \overline{Q}$ and all $t \in (0, t_1)$ one has

$$ \begin{aligned} & |G_d(t_2 - t, x, y) - G_d(t_1 - t, x, y)| \\ & \le (|G_d(t_2 - t, x, y) - G_d(t_1 - t, x, y)|)^\gamma (G_d(t_2 - t, x, y) + G_d(t_1 - t, x, y))^{1-\gamma} \\ & \le C_\gamma (t_1 - t)^{-\frac{d}{2}-\gamma} |t_2 - t_1|^\gamma e^{-c \frac{|y-x|^2}{t_1 - t}}. \end{aligned} $$

Applying the change of variables $(y', z') = (y - x, z - x)$, one has the estimates

$$ \int_0^{t_1} \|G(t_2 - t, x, \cdot) - G(t_1 - t, x, \cdot)\|_{(\alpha)}^2 dt \\ \le C_{\gamma} |t_2 - t_1|^{2\gamma} \int_0^{t_1} (t_1 - t)^{-d-2\gamma} \int_{\mathbb{R}^d \times \mathbb{R}^d} e^{-c\gamma \frac{|y-x|^2}{t_1-t}} e^{-c\gamma \frac{|z-x|^2}{t_1-t}} |y-z|^{-\alpha} dy dz dt. $$

Applying the auxiliary inequality (47) from the proof of the inequality (12), one then obtains

$$ \int_0^{t_1} \|G(t_2 - t, x, \cdot) - G(t_1 - t, x, \cdot)\|_{(\alpha)}^2 dt \le C_{\alpha} |t_2 - t_1|^{2\gamma} \int_0^{t_1} (t_1 - t)^{-2\gamma - \frac{\alpha}{2}} dt. $$

Assuming that $\gamma \in (0, 1 - \frac{\alpha}{2})$, one has $2\gamma + \frac{\alpha}{2} < 1$, and thus one gets the relation

$$ \int_0^{t_1} (t_1 - t)^{-2\gamma - \frac{\alpha}{2}} dt = \int_0^{t_1} t^{-2\gamma - \frac{\alpha}{2}} dt = C_{\alpha, \gamma} t_1^{1-2\gamma - \frac{\alpha}{2}}. $$

As a result, there exists $C_{\alpha, \gamma} \in (0, \infty)$ such that one has the bound

$$ \int_0^{t_1} \|G(t_2 - t, x, \cdot) - G(t_1 - t, x, \cdot)\|_{(\alpha)}^2 dt \le C_{\alpha, \gamma} t_1^{1-2\gamma - \frac{\alpha}{2}} |t_2 - t_1|^{2\gamma}. $$

The proof of the inequality (15) is thus completed. $\square$

It remains to provide the proof of Lemma 5.

Proof of Lemma 5. Applying the Burkholder-Davis-Gundy inequality (4), one has the upper bound

$$ \mathbb{E} \left[ \left| \int_0^t \int_Q G_d(t - \kappa(s), x, y) X(s, y) F_\alpha(\mathrm{d}s, \mathrm{d}y) \right|^{2p} \right] \le C_p(T) \mathbb{E} \left[ \left| \int_0^t \|G_d(t - \kappa(s), x, \cdot) X(s, \cdot)\|_{(\alpha)}^2 \mathrm{d}s \right|^p \right]. $$
Owing to the definition of the norm $||| \cdot |||_{2p}$, one obtains

$$ \begin{aligned} ||| \int_0^t \int_Q G_d(t - \kappa(s), x, y) X(s, y) F_\alpha(\mathrm{d}s, \mathrm{d}y) |||_{2p} &= \left( \mathbb{E} \left[ \left| \int_0^t \int_Q G_d(t - \kappa(s), x, y) X(s, y) F_\alpha(\mathrm{d}s, \mathrm{d}y) \right|^{2p} \right] \right)^{\frac{1}{2p}} \\ &\le C_p(T) \left( ||| \int_0^t ||| G_d(t - \kappa(s), x, \cdot) X(s, \cdot) |||_{(\alpha)}^2 \mathrm{d}s |||_p \right)^{\frac{1}{2}}. \end{aligned} $$

Applying the Minkowski inequality for the norm $||| \cdot |||_p$ and recalling the definition of the norm $|| \cdot ||_{(\alpha)}$, one obtains

$$ \begin{aligned} ||| \int_0^t ||| G_d(t - \kappa(s), x, \cdot) X(s, \cdot) |||_{(\alpha)}^2 \mathrm{d}s |||_p &\le \int_0^t ||| G_d(t - \kappa(s), x, \cdot) X(s, \cdot) |||_{(\alpha)}^2 |||_p \mathrm{d}s \\ &\le \int_0^t \iint_{Q \times Q} G_d(t - \kappa(s), x, y) G_d(t - \kappa(s), x, z) ||| X(s, y) X(s, z) |||_p |y - z|^{-\alpha} \mathrm{d}y \mathrm{d}z \mathrm{d}s. \end{aligned} $$

Note that by the Cauchy–Schwarz inequality, for all $s \ge 0$, one has

$$ \sup_{y,z \in \bar{Q}} ||| X(s,y)X(s,z) |||_p \le \sup_{y \in \bar{Q}} ||| X(s,y) |||_{2p}^2. $$

As a result, applying the inequality (12) from Lemma 4, one obtains

$$ \begin{aligned} ||| \int_0^t ||| G_d(t - \kappa(s), x, \cdot) X(s, \cdot) |||_{(\alpha)}^2 \mathrm{d}s |||_p &\le \int_0^t ||| G_d(t - \kappa(s), x, \cdot) |||_{(\alpha)}^2 \sup_{y \in \bar{Q}} ||| X(s, y) |||_{2p}^2 \mathrm{d}s \\ &\le C_\alpha(T) \int_0^t (t - \kappa(s))^{-\frac{\alpha}{2}} \sup_{y \in \bar{Q}} ||| X(s, y) |||_{2p}^2 \mathrm{d}s. \end{aligned} $$

This shows inequality (16).

To prove the inequality (17), it suffices to apply the Jensen inequality: one has

$$ \begin{aligned} & \left( \int_0^t (t - \kappa(s))^{-\frac{\alpha}{2}} \sup_{x \in \bar{Q}} ||| X(s, x) |||_{2p}^2 \mathrm{d}s \right)^p \\ &\le \left( \int_0^t (t - \kappa(s))^{-\frac{\alpha}{2}} \mathrm{d}s \right)^{p-1} \int_0^t (t - \kappa(s))^{-\frac{\alpha}{2}} \sup_{x \in \bar{Q}} ||| X(s, x) |||_{2p}^2 \mathrm{d}s \\ &\le \left( \int_0^t (t - \kappa(s))^{-\frac{\alpha}{2}} \mathrm{d}s \right)^{p-1} \int_0^t (t - \kappa(s))^{-\frac{\alpha}{2}} \sup_{x \in \bar{Q}} \mathbb{E}[|X(s, x)|^{2p}] \mathrm{d}s. \end{aligned} $$

Recalling that $\alpha < 2$ and that $\kappa(s) \le s$ for all $s \ge 0$, one has the upper bound

$$ \int_0^t (t - \kappa(s))^{-\alpha/2} \mathrm{d}s \le \int_0^t (t - s)^{-\alpha/2} \mathrm{d}s \le C_\alpha(T). $$

It is then straightforward to complete the proof of the inequality (17).
To prove the inequality (18), it suffices to consider the supremum over $s \in [0, T]$; one obtains

$$ \begin{aligned} \int_0^t (t - \kappa(s))^{-\frac{\alpha}{2}} \sup_{x \in Q} \mathbb{E}[|X(s, x)|^{2p}] ds &\le \int_0^t (t - \kappa(s))^{-\frac{\alpha}{2}} ds \sup_{0 \le s \le t} \sup_{x \in Q} \|X(s, x)\|_{2p}^{2p} \\ &\le C_{p,\alpha}(T) \sup_{0 \le t \le T} \sup_{x \in Q} \|X(t, x)\|_{2p}^{2p}. \end{aligned} $$

It is then straightforward to complete the proof of the inequality (18).

This concludes the proof of Lemma 5.

□

## APPENDIX B. IMPLEMENTATION DETAILS ON THE COVARIANCE MATRIX OF THE NOISE IN DIMENSION $d = 1$

In this appendix, we provide some details on the implementation of the discretization of the noise in Section 5.1 in $d = 1$ dimension. We recall that $\alpha \in (0, 1)$ and that the covariance matrix $C \in \mathbb{R}^{(n-1) \times (n-1)}$ has elements given by

$$ C_{ij} = \int_{x_i}^{x_i+\Delta x} \int_{x_j}^{x_j+\Delta x} |x-y|^{-\alpha} dx dy, \text{ for } i,j = 1, \dots, n-1, $$

where we set $\Delta x = \frac{1}{n}$ for the mesh size. These integrals can explicitly be computed as seen in the next lemma.

**Lemma 9.** For $i < j$, one has

$$ C_{ij} = C_{ji} = \frac{\Delta x^{2-\alpha}}{(1-\alpha)(2-\alpha)} \left( (j-i+1)^{2-\alpha} - 2(j-i)^{2-\alpha} + (j-i-1)^{2-\alpha} \right). $$

Furthermore, for all $i \in \{1, \dots, n-1\}$, one has

$$ C_{ii} = \frac{2\Delta x^{2-\alpha}}{(1-\alpha)(2-\alpha)}, \text{ for } i = 1, \dots, n-1. $$

*Proof.* We first consider the case $i < j$. This means that $x_i + \Delta x \le x_j$ and that $y \le x$ for all $y \in [x_i, x_i + \Delta x]$ and for all $x \in [x_j, x_j + \Delta x]$. One then has

$$ \int_{x_j}^{x_j+\Delta x} |x-y|^{-\alpha} dx = \int_{x_j}^{x_j+\Delta x} (x-y)^{-\alpha} dx = \frac{(x_j+\Delta x-y)^{1-\alpha}}{1-\alpha} - \frac{(x_j-y)^{1-\alpha}}{1-\alpha}. $$

Integrating the above with respect to the variable $y$, one thus obtains

$$ C_{ij} = \frac{1}{1-\alpha} \int_{x_i}^{x_i+\Delta x} (x_j + \Delta x - y)^{1-\alpha} dy - \frac{1}{1-\alpha} \int_{x_i}^{x_i+\Delta x} (x_j - y)^{1-\alpha} dy. $$

The first integral gives

$$ \int_{x_i}^{x_i+\Delta x} (x_j + \Delta x - y)^{1-\alpha} dy = \frac{1}{2-\alpha} \left( (x_j + \Delta x - x_i)^{2-\alpha} - (x_j - x_i)^{2-\alpha} \right). $$

The second integral gives

$$ \int_{x_i}^{x_i+\Delta x} (x_j - y)^{1-\alpha} dy = \frac{1}{2-\alpha} \left( (x_j - x_i)^{2-\alpha} - (x_j - x_i - \Delta x)^{2-\alpha} \right). $$
Thus, one obtains

$$C_{ij} = \frac{1}{1-\alpha} \frac{1}{2-\alpha} \left( (x_j + \Delta x - x_i)^{2-\alpha} - 2(x_j - x_i)^{2-\alpha} + (x_j - x_i - \Delta x)^{2-\alpha} \right).$$

Finally, recalling that $x_j = j\Delta x$ and $x_i = i\Delta x$, one obtains the expression

$$C_{ij} = \frac{\Delta x^{2-\alpha}}{(1-\alpha)(2-\alpha)} \left( (j-i+1)^{2-\alpha} - 2(j-i)^{2-\alpha} + (j-i-1)^{2-\alpha} \right).$$

To compute $C_{ii}$, we first compute the integral

$$\int_{x_i}^{x_i+\Delta x} |x-y|^{-\alpha} dx = \int_{x_i}^{y} (y-x)^{-\alpha} dx + \int_{y}^{x_i+\Delta x} (x-y)^{-\alpha} dx = \frac{1}{1-\alpha} \left( (y-x_i)^{1-\alpha} + (x_i + \Delta x - y)^{1-\alpha} \right).$$

Thus, integrating with respect to the variable $y$ gives us that

$$C_{ii} = \frac{1}{1-\alpha} \int_{x_i}^{x_i+\Delta x} (y-x_i)^{1-\alpha} + (x_i + \Delta x - y)^{1-\alpha} dy = \frac{2\Delta x^{2-\alpha}}{(1-\alpha)(2-\alpha)}.$$

This concludes the proof.

<figure><img src="image_15.png" /></figure>

## APPENDIX C. IMPLEMENTATION DETAILS ON THE COVARIANCE MATRIX OF THE NOISE IN DIMENSION $d = 2$

In this appendix, we provide some details on the implementation of the discretization of the noise in Section 5.1 in $d = 2$ dimensions. We recall that $\alpha \in (0, 2)$ and that the covariance matrix $C \in \mathbb{R}^{(n-1)^2 \times (n-1)^2}$ has elements given by

$$C_{ij} = \mathrm{Cov}(F_{\alpha}^{n}(t, x_i), F_{\alpha}^{n}(t, x_j)) = \int_{\Box_{x_i}} \int_{\Box_{x_j}} ||z_1 - z_2||^{-\alpha} dz_1 dz_2, \text{ for } i, j = 1, \dots, (n-1)^2.$$

As we are considering $d = 2$ in this section, $z_1 \in \Box_{x_i} \subset \mathbb{R}^2$ and $z_2 \in \Box_{x_j} \subset \mathbb{R}^2$. For the $i$th space grid point $x_i \in [0, 1]^2$, let us write $x_i = (x_i^1, x_i^2) \in [0, 1]^2$. $x_j^1$ and $x_j^2$ are defined in the analogous way for $x_j \in [0, 1]^2$. The following lemma reduces the computation of the 4-dimensional integrals in $C_{ij}$ to the computation of 2-dimensional integrals.

**Lemma 10.** For $i, j = 1, \dots, (n-1)^2$, one has that

$$C_{ij} = \begin{cases} \Delta x^{4-\alpha} \int_{-1}^{1} \int_{-1}^{1} ((A+q_1)^2 + (B+q_2)^2)^{-\alpha/2} (1-|q_1|)(1-|q_2|) dq_1 dq_2, & i \neq j, \\ \Delta x^{4-\alpha} I(\alpha), & i = j. \end{cases}$$

Here, $A = \frac{I_x(i) - I_x(j)}{\Delta x}$, $B = \frac{I_y(i) - I_y(j)}{\Delta x}$, and $I(\alpha) = \Delta_2(-\alpha)$ is defined by

$$\Delta_2(s) = 8 \frac{(3+s)2^{s/2} + 1}{(s+2)(s+3)(s+4)} + 4B_2(s) - \frac{4(s+4)}{s+2}B_2(s+2), \ s \in \mathbb{C},$$

where

$$B_2(s) = \frac{2}{2+s} {}_2F_1(1/2, -s/2; 3/2; -1), \ s \in \mathbb{C},$$

and ${}_2F_1$ is the hypergeometric function, see [4].
Proof. We write $z_1 = (x_1, y_1) \in [0, 1]^2$ and $z_2 = (x_2, y_2) \in [0, 1]^2$, so that $C_{ij}$ can be written as

$$C_{ij} = \int_{I_x(i)}^{I_x(i)+\Delta x} \int_{I_y(i)}^{I_y(i)+\Delta x} \int_{I_x(j)}^{I_x(j)+\Delta x} \int_{I_y(j)}^{I_y(j)+\Delta x} ((x_1 - x_2)^2 + (y_1 - y_2)^2)^{-\alpha/2} dx_1 dx_2 dy_1 dy_2.$$

Here, we integrate over the following: $x_1 \in [I_x(i), I_x(i) + \Delta x]$, $y_1 \in [I_y(i), I_y(i) + \Delta x]$, $x_2 \in [I_x(j), I_x(j) + \Delta x]$ and $y_2 \in [I_y(j), I_y(j) + \Delta x]$. In other words, $C_{ij}$ is of the more general form

$$(*)= \int_{a}^{a+\Delta x} \int_{b}^{b+\Delta x} \int_{c}^{c+\Delta x} \int_{d}^{d+\Delta x} ((x_1 - x_2)^2 + (y_1 - y_2)^2)^{-\alpha/2} dx_1 dy_1 dx_2 dy_2,$$

for different values of $a, b, c, d \in [0, 1]$. Here we integrate over the following: $x_1 \in [a, a + \Delta x]$, $y_1 \in [b, b + \Delta x]$, $x_2 \in [c, c + \Delta x]$ and $y_2 \in [d, d + \Delta x]$.

Let us now introduce new variables $\tilde{x}_1, \tilde{x}_2, \tilde{y}_1$ and $\tilde{y}_2$ as follows:

$$x_1 = a + \tilde{x}_1 \Delta x, \quad dx_1 = \Delta x \, d\tilde{x}_1$$

$$x_2 = a + \tilde{x}_2 \Delta x, \quad dx_2 = \Delta x \, d\tilde{x}_2$$

$$y_1 = a + \tilde{y}_1 \Delta x, \quad dy_1 = \Delta x \, d\tilde{y}_1$$

$$y_2 = a + \tilde{y}_2 \Delta x, \quad dy_2 = \Delta x \, d\tilde{y}_2.$$

Then, for the above 4 dimensional integral we have that

$$(*)= \Delta x^4 \int_{0}^{1} \int_{0}^{1} \int_{0}^{1} \int_{0}^{1} ((a - c + \Delta x(\tilde{x}_1 - \tilde{x}_2))^2 + (b - d + \Delta x(\tilde{y}_1 - \tilde{y}_2))^2)^{-\alpha/2} d\tilde{x}_1 d\tilde{y}_1 d\tilde{x}_2 d\tilde{y}_2.$$

Relabeling back to $x_1, x_2, y_1, y_2$ and factorizing out $\Delta x^\alpha$, one then obtains

$$(50) \quad (*) = \Delta x^{4-\alpha} \int_{0}^{1} \int_{0}^{1} \int_{0}^{1} \int_{0}^{1} ((A + (x_1 - x_2))^2 + (B + (y_1 - y_2))^2)^{-\alpha/2} dx_1 dy_1 dx_2 dy_2,$$

where $A = \frac{a-c}{\Delta x}$ and $B = \frac{b-d}{\Delta x}$.

If $i = j$ (meaning that $A = 0$ and $B = 0$), then (*) reduces to

$$(*)= \Delta x^{4-\alpha} I(\alpha),$$

where we recall that $I(\alpha)$ is defined in the statement of the lemma. We refer to [4] for more details.

If $i \neq j$ (meaning that $A \neq 0$ or $B \neq 0$), then we first consider

$$(\square) = \int_{0}^{1} \int_{0}^{1} ((A + (x_1 - x_2))^2 + (B + (y_1 - y_2))^2)^{-\alpha/2} dx_1 dx_2.$$

Instead of integration over $[0, 1] \times [0, 1]$ with respect to the measure $dx_1 dx_2$, we integrate along $r \in [-\sqrt{2}/2, \sqrt{2}/2]$ (to match the domain $[0, 1] \times [0, 1]$) with respect to the measure $2(\sqrt{2}/2 - |r|) dr$

$$(\square) = \int_{-\sqrt{2}/2}^{\sqrt{2}/2} ((A + 2r/\sqrt{2})^2 + (B + (y_1 - y_2))^2)^{-\alpha/2} \times 2(\sqrt{2}/2 - |r|) dr.$$

We now change variables according to $q = \frac{2}{\sqrt{2}}r$ and obtain that

$$(\square) = \int_{-1}^{1} ((A + q)^2 + (B + (y_1 - y_2))^2)^{-\alpha/2} (1 - |q|) dq.$$
Repeating the above also for the double integrals involving $y_1, y_2$ in (50), gives us

$$ (*) = \Delta x^{4-\alpha} \int_{-1}^{1} \int_{-1}^{1} ((A+q_1)^2 + (B+q_2)^2)^{-\alpha/2} (1-|q_1|)(1-|q_2|) dq_1 dq_2, $$

where we also added back the factor of $\Delta x^{4-\alpha}$.

<figure><img src="image_16.png" /></figure>

In contrast to the dimension $d=1$ in Appendix B, some integrals in $C_{ij}$ cannot be computed exactly. On the one hand, the diagonal elements $C_{ii}$ involve the integral $I(\alpha)$ which can be computed using matlab's command *hypergeom* to evaluate the hypergeometric function. On the other hand, the other elements of the covariance matrix need numerical approximations, where we used the matlab command *integral2*.

## REFERENCES

[1] R. Anton and D. Cohen. Exponential integrators for stochastic Schrödinger equations driven by Itô noise. *J. Comput. Math.*, 36(2):276-309, 2018.

[2] R. Anton, D. Cohen, S. Larsson, and X. Wang. Full discretization of semilinear stochastic wave equations driven by multiplicative noise. *SIAM J. Numer. Anal.*, 54(2):1093-1119, 2016.

[3] R. Anton, D. Cohen, and L. Quer-Sardanyons. A fully discrete approximation of the one-dimensional stochastic heat equation. *IMA J. Numer. Anal.*, 40(1):247-284, 2020.

[4] D. H. Bailey, J. M. Borwein, and R. E. Crandall. Advances in the theory of box integrals. *Math. Comp.*, 79(271):1839-1866, 2010.

[5] S. Becker, B. Gess, A. Jentzen, and P. E. Kloeden. Lower and upper bounds for strong approximation errors for numerical approximations of stochastic heat equations. *BIT*, 60(4):1057-1073, 2020.

[6] A. Berg, D. Cohen, and G. Dujardin. Approximated exponential integrators for the stochastic Manakov equation. *J. Comput. Dyn.*, 10(2):323-344, 2023.

[7] C.-E. Bréhier. Total variation error bounds for the accelerated exponential Euler scheme approximation of parabolic semilinear SPDEs. *SIAM J. Numer. Anal.*, 62(3):1171-1190, 2024.

[8] C.-E. Bréhier, D. Cohen, and J. Ulander. Analysis of a positivity-preserving splitting scheme for some semilinear stochastic heat equations. *ESAIM Math. Model. Numer. Anal.*, 58(4):1317-1346, 2024.

[9] C.-E. Bréhier and X. Wang. On parareal algorithms for semilinear parabolic stochastic PDEs. *SIAM J. Numer. Anal.*, 58(1):254-278, 2020.

[10] O. Butkovsky, K. Dareiotis, and M. Gerencsér. Optimal rate of convergence for approximations of SPDEs with nonregular drift. *SIAM J. Numer. Anal.*, 61(2):1103-1137, 2023.

[11] M. Cai, D. Cohen, and X. Wang. Strong convergence rates for a full discretization of stochastic wave equation with nonlinear damping. *J. Sci. Comput.*, 102(2):Paper No. 45, 26, 2025.

[12] C. Cardon-Weber. Implicit approximation scheme for the Cahn-Hilliard stochastic equation. *Prépublication 613 du Laboratoire de Probabiliteés et Modèles Aléatoires*, 2000.

[13] C. Cardon-Weber and A. Millet. On strongly Petrovskii's parabolic SPDEs in arbitrary dimension and application to the stochastic Cahn-Hilliard equation. *J. Theoret. Probab.*, 17(1):1-49, 2004.

[14] C. Chen, J. Cui, J. Hong, and D. Sheng. Accelerated exponential Euler scheme for stochastic heat equation: convergence rate of the density. *IMA J. Numer. Anal.*, 43(2):1181-1220, 2023.

[15] C. Chen, T. Dang, and J. Hong. An adaptive time-stepping fully discrete scheme for stochastic NLS equation: strong convergence and numerical asymptotics. *Stochastic Process. Appl.*, 173:Paper No. 104373, 24, 2024.

[16] J. Clausnitzer and A. Kleefeld. A spectral Galerkin exponential Euler time-stepping scheme for parabolic SPDEs on two-dimensional domains with a $C^2$ boundary. *Discrete Contin. Dyn. Syst. Ser. B*, 29(4):1624-1651, 2024.

[17] D. Cohen and G. Dujardin. Exponential integrators for nonlinear Schrödinger equations with white noise dispersion. *Stoch. Partial Differ. Equ. Anal. Comput.*, 5(4):592-613, 2017.

[18] D. Cohen and A. Lang. Numerical approximation and simulation of the stochastic wave equation on the sphere. *Calcolo*, 59(3):Paper No. 32, 32, 2022.

[19] D. Cohen, S. Larsson, and M. Sigg. A trigonometric method for the linear stochastic wave equation. *SIAM J. Numer. Anal.*, 51(1):204-222, 2013.
[20] D. Cohen and L. Quer-Sardanyons. A fully discrete approximation of the one-dimensional stochastic wave equation. *IMA J. Numer. Anal.*, 36(1):400–420, 2016.

[21] S. Cox, A. Jentzen, and F. Lindner. Weak convergence rates for temporal numerical approximations of the semilinear stochastic wave equation with multiplicative noise. *Numer. Math.*, 156(6):2131–2177, 2024.

[22] S. Cox and J. van Neerven. Convergence rates of the splitting scheme for parabolic linear stochastic Cauchy problems. *SIAM J. Numer. Anal.*, 48(2):428–451, 2010.

[23] J. Cui. Explicit approximation for stochastic nonlinear Schrödinger equation. *J. Differential Equations*, 419:1–39, 2025.

[24] J. Cui and J. Hong. Strong and weak convergence rates of a spatial approximation for stochastic partial differential equation with one-sided Lipschitz coefficient. *SIAM J. Numer. Anal.*, 57(4):1815–1841, 2019.

[25] R. C. Dalang. Extending the martingale measure stochastic integral with applications to spatially homogeneous s.p.d.e.'s. *Electron. J. Probab.*, 4:no. 6, 29 pp. 1999.

[26] R. C. Dalang and L. Quer-Sardanyons. Stochastic integrals for spde's: a comparison. *Expo. Math.*, 29(1):67–109, 2011.

[27] R. C. Dalang and M. Sanz-Solé. Stochastic Partial Differential Equations, Space-time White Noise and Random Fields, 2024.

[28] A. M. Davie and J. G. Gaines. Convergence of numerical schemes for the solution of parabolic stochastic partial differential equations. *Math. Comp.*, 70(233):121–134, 2001.

[29] A. Djurdjevac, M. Gerencsér, and H. Kremp. Higher order approximation of nonlinear SPDEs with additive space-time white noise, 2024.

[30] X. Gui, B. Li, and J. Wang. Improved error estimates for a modified exponential Euler method for the semilinear stochastic heat equation with rough initial data. *Sci. China Math.*, 67(12):2873–2898, 2024.

[31] I. Gyöngy. Lattice approximations for stochastic quasi-linear parabolic partial differential equations driven by space-time white noise. *I. Potential Anal.*, 9(1):1–25, 1998.

[32] I. Gyöngy and D. Nualart. Implicit scheme for quasi-linear parabolic partial differential equations perturbed by space-time white noise. *Stochastic Process. Appl.*, 58(1):57–72, 1995.

[33] A.-L. Haji-Ali, H. Hoel, and A. Petersson. The multi-index monte carlo method for semilinear stochastic partial differential equations, 2025.

[34] P. Harms and M. S. Müller. Weak convergence rates for stochastic evolution equations and applications to nonlinear stochastic wave, HJMM, stochastic Schrödinger and linearized stochastic Korteweg–de Vries equations. *Z. Angew. Math. Phys.*, 70(1):Paper No. 16, 28, 2019.

[35] E. Hausenblas. Approximation for semilinear stochastic evolution equations. *Potential Anal.*, 18(2):141–186, 2003.

[36] M. Hochbruck and A. Ostermann. Exponential integrators. *Acta Numer.*, 19:209–286, 2010.

[37] M. Hutzenthaler, A. Jentzen, and D. Salimova. Strong convergence of full-discrete nonlinearity-truncated accelerated exponential Euler-type approximations for stochastic Kuramoto-Sivashinsky equations. *Commun. Math. Sci.*, 16(6):1489–1529, 2018.

[38] A. Jentzen and P. E. Kloeden. The numerical approximation of stochastic partial differential equations. *Milan J. Math.*, 77:205–244, 2009.

[39] A. Jentzen and P. E. Kloeden. Overcoming the order barrier in the numerical approximation of stochastic partial differential equations with additive space-time noise. *Proc. R. Soc. Lond. Ser. A Math. Phys. Eng. Sci.*, 465(2102):649–667, 2009.

[40] D. Khoshnevisan. *Analysis of stochastic partial differential equations*, volume 119. Published for the Conference Board of the Mathematical Sciences, Washington, DC; by the American Mathematical Society, Providence, RI, 2014.

[41] K. Klioba and M. Veraar. Temporal approximation of stochastic evolution equations with irregular nonlinearities. *J. Evol. Equ.*, 24(2):Paper No. 43, 34, 2024.

[42] K. Klioba and M. Veraar. Pathwise uniform convergence of time discretization schemes for SPDEs. *IMA J. Numer. Anal.*, 45(4):2060–2131, 2025.

[43] P. E. Kloeden, G. J. Lord, A. Neuenkirch, and T. Shardlow. The exponential integrator scheme for stochastic partial differential equations: pathwise error bounds. *J. Comput. Appl. Math.*, 235(5):1245–1260, 2011.

[44] M. Kovács, S. Larsson, and F. Lindgren. Strong convergence of the finite element method with truncated noise for semilinear parabolic stochastic equations with additive noise. *Numer. Algorithms*, 53(2-3):309–320, 2010.
[45] O. A. Ladyženskaja. Linear and Quasilinear Equations of Parabolic Type. 1968.

[46] Z. Lei, C.-E. Bréhier, and S. Gan. Numerical approximation of the invariant distribution for a class of stochastic damped wave equations. *J. Comput. Math.*, 43(4):976–1015, 2025.

[47] R. Liu and K. Du. Hölder continuity of solutions to the Dirichlet problem for SPDEs with spatially correlated noise. *Electron. Commun. Probab.*, 26:Paper No. 61, 13, 2021.

[48] Z. Liu and Z. Qiao. Strong approximation of monotone stochastic partial differential equations driven by white noise. *IMA J. Numer. Anal.*, 40(2):1074–1093, 2020.

[49] G. J. Lord, C. E. Powell, and T. Shardlow. *An introduction to computational stochastic PDEs*. Cambridge Texts in Applied Mathematics. Cambridge University Press, New York, 2014.

[50] G. J. Lord and A. Tambue. Stochastic exponential integrators for the finite element discretization of SPDEs for multiplicative and additive noise. *IMA J. Numer. Anal.*, 33(2):515–543, 2013.

[51] A. Millet and P.-L. Morien. On implicit and explicit discretization schemes for parabolic SPDEs in any dimension. *Stochastic Process. Appl.*, 115(7):1073–1106, 2005.

[52] C. M. Mora. Numerical solution of conservative finite-dimensional stochastic Schrödinger equations. *Ann. Appl. Probab.*, 15(3):2144–2171, 2005.

[53] M. Ondreját, A. Prohl, and N. J. Walkington. Numerical approximation of nonlinear SPDE's. *Stoch. Partial Differ. Equ. Anal. Comput.*, 11(4):1553–1634, 2023.

[54] J. Printems. On the discretization in time of parabolic stochastic partial differential equations. *M2AN Math. Model. Numer. Anal.*, 35(6):1055–1078, 2001.

[55] R. Qi and X. Wang. An accelerated exponential time integrator for semi-linear stochastic strongly damped wave equation with additive noise. *J. Math. Anal. Appl.*, 447(2):988–1008, 2017.

[56] J. B. Walsh. An introduction to stochastic partial differential equations. In *École d'été de probabilités de Saint-Flour, XIV–1984*, volume 1180 of *Lecture Notes in Math.*, pages 265–439. Springer, Berlin, 1986.

[57] M. Wang and X. Wang. A linearly implicit finite element full-discretization scheme for SPDEs with nonglobally Lipschitz coefficients. *IMA J. Numer. Anal.*, 45(1):516–579, 2025.

[58] X. Wang. An exponential integrator scheme for time discretization of nonlinear stochastic wave equation. *J. Sci. Comput.*, 64(1):234–263, 2015.

[59] X. Wang. Strong convergence rates of the linear implicit Euler method for the finite element discretization of SPDEs with additive noise. *IMA J. Numer. Anal.*, 37(2):965–984, 2017.

[60] X. Wang, S. Gan, and J. Tang. Higher order strong approximations of semilinear stochastic wave equation with additive space-time white noise. *SIAM J. Sci. Comput.*, 36(6):A2611–A2632, 2014.

[61] X. Wang and R. Qi. A note on an accelerated exponential Euler method for parabolic SPDEs with additive noise. *Appl. Math. Lett.*, 46:31–37, 2015.

[62] Y. Yan. Semidiscrete Galerkin approximation for a linear stochastic parabolic partial differential equation driven by an additive noise. *BIT*, 44(4):829–847, 2004.

[63] Y. Zhang, X. Yang, and R. Qi. Difference approximation of stochastic elastic equation driven by infinite dimensional noise. *Numer. Math. Theory Methods Appl.*, 9(1):123–146, 2016.

## 4. ACKNOWLEDGMENTS

We would like to thank David Krantz (KTH Royal Institute of Technology) for interesting discussions on the numerical approximations of singular quadruple integrals. The work was initiated during a research stay of the authors (DC and LQS) at the Bernoulli Center (EPFL, Switzerland) under a Bernoulli Brainstorm program. We would like to thank the staff of the Bernoulli Center for its support. Part of this work was carried while DC was working for Umeå University. Part of this work was carried thanks to the SFVE-A mobility program between France and Sweden. The work of DC and JU was partially supported by the Swedish Research Council (VR) (projects nr. 2018 – 04443 and 2024 – 04536). LQS is supported by the grant PID2021-123733NB-I00 (Ministerio de Economía y Competitividad, Spain). The computations were performed on resources provided by the Swedish National Infrastructure for Computing (SNIC) at HPC2N, Umeå University and by the National Academic Infrastructure for Supercomputing in Sweden (NAISS) at UPPMAX, Uppsala University, at Dardel,
KTH, and at Vera, Chalmers e-Commons at Chalmers University of Technology and partially funded by the Swedish Research Council through grant agreement no. 2022-06725.

UNIVERSITE DE PAU ET DES PAYS DE L'ADOUR, E2S UPPA, CNRS, LMAP, PAU, FRANCE

Email address: charles-edouard.brehier@univ-pau.fr

DEPARTMENT OF MATHEMATICAL SCIENCES, CHALMERS UNIVERSITY OF TECHNOLOGY AND UNIVERSITY OF GOTHENBURG, 41296 GOTHENBURG, SWEDEN

Email address: david.cohen@chalmers.se

DEPARTMENT OF MATHEMATICS, UNIVERSITAT AUTÒNOMA DE BARCELONA, 08193 BELLATERRA, SPAIN

Email address: lluis.quer@uab.cat

DEPARTMENT OF MATHEMATICAL SCIENCES, CHALMERS UNIVERSITY OF TECHNOLOGY AND UNIVERSITY OF GOTHENBURG, 41296 GOTHENBURG, SWEDEN

Current address: Swedish Defence Research Agency (FOI), 58330 Linköping, Sweden

Email address: johanul@chalmers.se

Email address: johan.ulander@foi.se