text:
\section{Introduction}

 Let $D$ be a bounded open subset of $\mathbb{R}^d$ (satisfying certain regularity conditions) and consider the following parabolic stochastic PDE on $[0,\infty)\times D$ with zero Dirichlet boundary conditions:
\begin{equation}\label{SHE}
\left\{\begin{array}{ll} \frac{\partial u}{\partial t}(t,x) = -(-\Delta)^{\gamma} u(t,x)+\sigma(u(t,x)) \dot W(t, x), & (t,x) \in [0,\infty)\times D,\\ u(t,x)=0, &  (t,x)\in[0,\infty)\times\partial D, \\
	u(0,x)=0, & x\in D. \end{array}\right.\end{equation}
Here, $\dot W$ is a Gaussian space-time white noise on $[0,\infty)\times D$, ${\sigma}\colon \mathbb{R}\to\mathbb{R}$ is a Lipschitz function, and $(-\Delta)^{\gamma}$ is a spectral power of $-\Delta$ (see Section~\ref{sec2} for details). The purpose of this article is study the regularity of $t\mapsto u(t,\cdot)$ as a stochastic process taking values in a scale of Sobolev spaces $H_r=H_r(D)$ indexed by $r\in\mathbb{R}$ (to be defined in Section~\ref{sec2}). More precisely,  we are interested in the asymptotic behavior as $n\to\infty$ of the \emph{(normalized) power variation of order $p$} given by
\begin{equation}\label{Vnp} V^{n,r}_p(u,t):= \Delta_n\sum_{i=1}^{[t/\Delta_n]} \biggl( \frac{\|u(i\Delta_n,\cdot)-u((i-1)\Delta_n,\cdot)\|_{H_r}}{\tau_n(r)} \biggr)^p,\quad t\in[0,\infty),\quad n\geq1. \end{equation}
In this paper, $\Delta_n$ is a strictly positive sequence decreasing to $0$ (e.g., $\Delta_n=\frac1n$), $p>0$ is a fixed but arbitrary power, and $\tau_n(r)$ is a normalizing factor depending on $r$  and $n$ (and ${\gamma}$), chosen if possible in such a way that $V^{n,r}_p(u,t)$ converges to a limit $V^r_p(u,t)$, say, uniformly on compact sets in probability.

In order to describe the flavor of our results, let us specialize to the case where $d=1$, $D=(0,\pi)$, ${\gamma}=1$, and ${\sigma}\equiv1$ (i.e., to the stochastic heat equation on an interval with additive noise) in this introductory part. As Proposition~\ref{uinHr} below shows, in order that $u(t,\cdot)\in H_r$ for $t>0$, the smoothness parameter $r$ must be taken in the range $(-\infty,\frac12)$.
\begin{theoremA}\label{thmA} Let $u$ be the solution of \eqref{SHE} with $d=1$, $D=(0,\pi)$, ${\gamma}=1$, and ${\sigma}\equiv1$. Assume $r<\frac12$ and define
\begin{equation}\label{taunr0} \tau_n(r):=\begin{cases} \Delta_n^{\frac12} &\text{if } r<-\frac12,\\ (\Delta_n \lvert\log\Delta_n\rvert)^{\frac12} & \text{if } r=-\frac12,\\ \Delta_n^{\frac14-\frac r2} &\text{if } -\frac12<r<\frac12. \end{cases} \end{equation}
If $p\geq1$ is an integer, then
\begin{equation}\label{conv} \lim_{n\to\infty} \mathbb{E}\Biggl[ \sup_{t\in[0,T]} \Bigl|V^{n,r}_{2p}(u,t)-K(r,p)t\Bigr|^q\Biggr] =0 \end{equation}
for every $q,T\in(0,\infty)$, where the constant $K(r,p)$ is given by
\begin{equation}\label{const} K(r,p):=\begin{cases} 2^p B_p\Bigl( \frac{(1-1)!}{2}\zeta(-2\cdot  1 \cdot r),\dots,\frac{(p-1)!}{2}\zeta(-2\cdot  p \cdot r) \Bigr)&\text{if } r<-\frac12,\\ 2^{-p}& \text{if } r=-\frac12,\\ \biggl( \frac{{\Gamma}(r+\frac12)}{2(\frac12-r)} \biggr)^p &\text{if } -\frac12<r<\frac12. \end{cases} \end{equation}
In the last formula, $B_p$ is the complete Bell polynomial in $p$ variables (see \cite[Definition 2.4.1]{Peccati11}), $\zeta(z):= \sum_{k=1}^\infty k^{-z}$, for $z>1$, is the Riemann zeta function, and ${\Gamma}(z):=\int_0^\infty y^{z-1}\mathrm{e}^{-y}\,\mathrm{d} y$ is Euler's gamma function.
\end{theoremA}

The appearance of the Riemann zeta function in formula \eqref{const} is somewhat unexpected. It is related to the fact that in spatial dimension $d=1$, the $k$th eigenvalue of the Laplacian is proportional to $k^2$ and the norm on $H_r$ is defined using these eigenvalues (see \eqref{Hr}). Somewhat surprisingly, for $d\geq 1$ and a wide class of bounded open sets $D$,  conclusions similar to those of Theorem \ref{thmA} remain valid, after replacing the Riemann zeta function by the so-called spectral zeta function (see Corollaries \ref{cor1} and \ref{cor2} and Remarks \ref{rem2.7} and \ref{rem2.10}).

As we can see from both the formula for the normalizing sequence in \eqref{taunr0} and the formula for the limiting constant in \eqref{const}, the behavior of $V^{n,r}_p(u,t)$ changes at the critical value $r=-\frac12$. In fact, our analysis shows that for $r<-\frac12$, the terms in the sum \eqref{Vnp} become nearly independent and identically distributed as $n$ gets large, so that the convergence in \eqref{conv} can be interpreted as a \emph{(weak) law of large numbers}. By contrast, for $r\in[-\frac12, \frac12)$, they become \emph{deterministic} as $n\to\infty$, so the law of large numbers  is degenerate in this case.
The existence of a critical value of $r$ for the limiting behavior of $V^{n,r}_p(u,t)$ persists if we consider the more general equation \eqref{SHE}. In fact, Theorem~\ref{thmA} is a special case of our main Theorems~\ref{LLN} and \ref{LLN2} and their Corollaries~\ref{cor1} and \ref{cor2}.

Various authors have previously studied the  regularity of solutions to equations like \eqref{SHE}. For instance, the joint space-time Hölder regularity of the solution to the stochastic heat equation on $\mathbb{R}^d$ was established in \cite{SanzSole02}. Hölder regularity of the solution   as a process in $H_r$ (or related spaces) was investigated, for example,   in \cite[Theorem~11.8]{Peszat07} and \cite{vanNeerven08}. When $d=1$ and $D$ is an interval or $\mathbb{R}$, the power variations of $t\mapsto u(t,x)$, for fixed $x\in\mathbb{R}$, were analyzed  in \cite{Bibinger19b, Bibinger19, Cialenco20, Pospisil07,Swanson07} for $p\in\{2,4\}$ and in \cite{Chong20} for general powers. In the same setting, for fixed $t>0$, the power variations of $x\mapsto u(t,x)$ were discussed in \cite{Cialenco20,Hildebrandt19,Pospisil07} for $p=2$ and in \cite{Foondun15} for ${\gamma}\in(\frac12,1]$ and $p=2/(2{\gamma}-1)$. In the context of one-parameter stochastic processes, power variations have been investigated for semimartingales \cite{Jacod12}, fractional Brownian motion and related processes \cite{Corcuera06,Corcuera09, Nourdin08,Nourdin10, Nualart18}, and moving average processes \cite{BN11,Basse17, Corcuera13}, just to name a few.

To our best knowledge, power variations in Sobolev-type spaces have not been considered in the literature before. Even the existence of an $L^2$-continuous random field solution to \eqref{SHE} seems not yet to have been considered (our Proposition \ref{exist}), since it relies on fairly strong estimates concerning series of eigenfunctions of the Laplacian (see e.g.~Lemma \ref{A1}). Also, as Theorem~\ref{thmA} or the more general Theorems~\ref{LLN} and \ref{LLN2} below reveal, power variations in Sobolev-type spaces already show a rich behavior in first order. In particular, there is a phase transition in the law of large numbers that is not present in, for example, the power variations of fractional Brownian motion (cf.\ \cite{BN11,Corcuera06,Corcuera13}) or the power variations of the stochastic heat equation in time at a fixed spatial point (cf.\ \cite{Bibinger19,Chong20,Cialenco20,Swanson07}).

In the remainder of this paper, we write $A\lesssim B$ if there exists a finite constant $C>0$ (that is independent of all quantities of interest) such that $A\leq CB$. Moreover, we use the notation $\mathbb{N}=\{1,2,\dots\}$ and $\mathbb{N}_0=\{0,1,\dots\}$.

\section{Definitions and results}\label{sec2}

We first recall some basic definitions.
Let $0<{\lambda}_1\leq{\lambda}_2\leq \cdots$ and $(\phi_k)_{k\in\mathbb{N}}$ be the eigenvalues and corresponding normalized eigenfunctions of $-\Delta$, where $\Delta$ is the Dirichlet Laplacian on $D$ with zero boundary conditions, such that $(\phi_k)_{k\in\mathbb{N}}$ forms a complete orthonormal basis of $L^2(D)$, and each $\phi_k$ is smooth \cite[Corollary 8.11]{Gilbarg01} and bounded: see Lemma \ref{A1} (ii). As in
\cite[Chapter~IV, Example~3]{Walsh86}), let $E_0$ be the set of $f$ of the form $f= \sum_{j=1}^N a_k \phi_k$ and define, for $r\in\mathbb{R}$,
$$
    \|f\|_{H_r} := \Biggl( \sum_{k=1}^{\infty} {\lambda}_k^r a_k^2 \Biggr)^{1/2}.
$$
Let $H_r:=H_r(D)$ be the completion of $E_0$ with respect to $\|\cdot \|_{H_r}$. Then $H_r$ is a Hilbert space, the \emph{Sobolev space of order $r$}, each element $\Phi$ of which can be identified with a series of the form
\begin{equation}\label{Hr}  \Phi=\sum_{k=1}^\infty a_k(\Phi)\phi_k, \qquad\text{where } a_k(\Phi) \in \mathbb{R} \text{ and } \|\Phi\|_{H_r} := \Biggl( \sum_{k=1}^{\infty} {\lambda}_k^r a_k(\Phi)^2 \Biggr)^{1/2}<\infty . \end{equation}
As noted in \cite[Remark~2.8]{Sheffield07}, the series defining $\Phi$ converges in the topology of the space of distributions on $D$ and in the $H_r$-norm. The inner product on $H_r$ is given by
\begin{equation}\label{Hr-ip} \langle \Phi_1,\Phi_2\rangle_{H_r}:=\sum_{k=1}^\infty {\lambda}_k^r a_k(\Phi)a_k(\Phi'),\qquad \Phi_1,\Phi_2\in H_r. \end{equation}
Moreover, $H_r\subseteq H_s$ for $s\leq r$, $H_0=L^2(D)$ with $\|\cdot\|_0=\|\cdot\|_{L^2(D)}$ by Parseval's identity, and the evaluation
$$
\langle \Phi_1, \Phi_2 \rangle := \sum_{k=1}^\infty a_k(\Phi_1) a_k(\Phi_2),\quad \Phi_1\in H_{-r},\quad \Phi_2\in H_r,
$$
puts $H_{-r}$ and $H_r$ in duality.

The \emph{spectral power of $-\Delta$ of order ${\gamma}\in(0,\infty)$} is now defined   via
\begin{equation}\label{power} (-\Delta)^{\gamma}\colon \bigcup_{r\in\mathbb{R}} H_r \to \bigcup_{r\in\mathbb{R}} H_r,\qquad (-\Delta)^{\gamma}\Phi := \sum_{k=1}^\infty {\lambda}_k^{\gamma} a_k(\Phi)\phi_k.  \end{equation}

It is standard to interpret the SPDE \eqref{SHE} using the notion of random field solution (see e.g.~\cite[Chapter 4]{Dalang20}), in which one gives a suitable meaning to the integral equation
\begin{equation}\label{eq:mild} u(t,x)  = \int_0^t \int_D g(t-s; x,y){\sigma}(u(s,y))\, W( \mathrm{d} s, \mathrm{d} y)\qquad\text{a.s.,}  \end{equation}
which should be satisfied for all $(t,x)\in[0,\infty)\times D$, where  the Dirichlet Green's kernel $g$ takes the form
\begin{equation} \label{eq:def-heat}
g(t;x,y) :=   \sum_{k=1}^\infty \phi_k(x)\phi_k(y) \mathrm{e}^{-{\lambda}^{\gamma}_k t}\mathds 1_{t> 0},\qquad (t,x,y)\in [0,\infty)\times D^2.
\end{equation}

The existence and uniqueness of a mild solution to \eqref{SHE} can be established under appropriate conditions on ${\gamma}$ and $D$. The proof of the next Proposition is given in the \hyperref[appn]{Appendix} after Lemma~\ref{A1}.

\begin{proposition}\label{exist} Let $D$ be a bounded open subset of $\mathbb{R}^d$ that satisfies the {\em cone property} (see \cite[Section~2]{Agmon65}). Let ${\sigma}$ be a globally Lipschitz function and let ${\gamma}\in(\frac d2,\infty)$. Then there is a predictable random field $(t,x)\mapsto u(t,x)$, called the \emph{random field solution to \eqref{SHE}} that satisfies \eqref{eq:mild} for all $(t,x)\in(0,\infty)\times D$ and is such that
\begin{equation}\label{eq:Lp} \sup_{(t,x)\in[0,T]\times D} \mathbb{E}[|u(t,x)|^p] < \infty \end{equation}
for all $p,T>0$. In addition, $u$ is $L^p({\Omega})$-continuous for all $p>0$.  Up to versions, $u$ is unique among all predictable random fields satisfying \eqref{eq:Lp} for $p=2$.
\end{proposition}

\begin{remark}
A hyperbolic analogue of \eqref{SHE} (with $D=\mathbb{R}^d$ and a spatially homogeneous noise) was considered by the second author together with M.\ Sanz-Solé in \cite{Dalang05}.
\end{remark}

In fact, our methods apply to a more general situation. Instead of studying \eqref{SHE} with a nonlinearity of a specific functional form, we shall consider
\begin{equation}\label{SHE-2}
\left\{\begin{array}{ll} \frac{\partial u}{\partial t}(t,x) = -(-\Delta)^{\gamma} u(t,x)+\sigma(t,x) \dot W(t, x), & (t,x) \in [0,\infty)\times D,\\ u(t,x)=0, &  (t,x)\in[0,\infty)\times\partial D, \\
	u(0,x)=0, & x\in D, \end{array}\right.\end{equation}
with a general predictable random field ${\sigma}(t,x)$ satisfying mild regularity assumptions (so that \eqref{SHE} is included as a special case).
When ${\gamma} >\frac d2$, the random field solution to \eqref{SHE-2} is given  by
\[
 u(t,x)= \int_0^t \int_D g(t-s;x,y) {\sigma}(s,y)\,W(\mathrm{d} s,\mathrm{d} y)
\]
for $t\in(0,\infty)$ and $x\in D$.
By the stochastic Fubini theorem \cite[Theorem 2.6]{Walsh86}, whose assumption is satisfied by \eqref{est1} because $\gamma > \tfrac{d}{2}$,
\begin{equation}
u(t,\cdot) =\sum_{k=1}^\infty  a_k(t)\phi_k \qquad\text{with}\qquad a_k(t) :=\int_0^t\int_D \mathrm{e}^{-{\lambda}_k^{\gamma}(t-s)}\phi_k(y){\sigma}(s,y)\,W(\mathrm{d} s,\mathrm{d} y)
\label{akt}
\end{equation}
for $k\in\mathbb{N}$, where the series converges in $L^2(\Omega)$.

Another advantage of studying \eqref{SHE-2} is that we can drop the condition ${\gamma}>\frac d2$ (which guarantees a random field solution and is needed for \eqref{SHE} in the case of multiplicative noise) if we allow $t\mapsto u(t,\cdot)$ to be a distribution-valued process. Indeed, by a direct calculation and Weyl's law (see \eqref{eq:Weyl}), we obtain for $t>0$ that
\begin{equation}\label{eq:L2finite}
\mathbb{E}\Biggl[\sum_{k=1}^\infty {\lambda}_k^r (a_k(t))^2 \Biggr] = \frac12\sum_{k=1}^\infty {\lambda}_k^{r-{\gamma}} (1-\mathrm{e}^{-2{\lambda}_k^{\gamma} t})< \infty \iff r<{\gamma}-\tfrac d2.
\end{equation}
Therefore, if $r<{\gamma}-\tfrac d2$ and $\sigma$ satisfies \eqref{mom} below with $q=2$, then, for every $t\geq0$, the series in \eqref{akt} converges in $H_r$ almost surely and hence \emph{defines} a random variable $u(t,\cdot)$ in $H_r$ with $\mathbb{E}[\|u(t,\cdot)\|^2_{H_r}]< \infty$. As shown in \cite[Proposition~5.3]{Walsh86}, the $H_r$-valued process $t\mapsto u(t,\cdot)$ obtained in this way is  the \emph{weak solution} to \eqref{SHE-2} in the sense of \cite[Equation~(5.4)]{Walsh86} (see also \cite[Chapter~3]{Dalang20}). Let us remark that if ${\gamma}=1$ and ${\sigma}\equiv 1$, then the weak solution to \eqref{SHE-2} is a dynamical analogue of the \emph{Gaussian free field} considered, for example, in \cite{Giacomin01, Sheffield07}.   The restriction  $r<{\gamma}-\tfrac d2$ cannot be removed as the following result shows (see the \hyperref[appn]{Appendix} for a proof):
\begin{proposition}\label{uinHr} If ${\sigma}\equiv1$, then for all $t>0$, we have $u(t,\cdot)\in H_r$ a.s.~if and only if $r<{\gamma}-\frac d2$.
\end{proposition}

It is therefore natural to consider power variations in the spaces $H_r$ with $r< {\gamma}-\tfrac d 2$. In fact,  more general functionals can be considered. Indeed, if $r<-\frac d2$, we shall investigate the \emph{(normalized) $F$-variations} of the solution $u$ to \eqref{SHE-2}, that is,
\begin{equation}\label{VnF} V^{n,r}_F(u,t):= \Delta_n\sum_{i=1}^{[t/\Delta_n]} F\biggl( \frac{u(i\Delta_n,\cdot)-u((i-1)\Delta_n,\cdot)}{\tau_n(r)} \biggr),\quad t\in[0,\infty),\quad n\geq1, \end{equation}
where $F\colon H_r\to\mathbb{R}$ is a continuous functional on $H_r$ satisfying certain regularity properties and the $\tau_n(r)$ are normalizing constants that we will define later in \eqref{taunr}.

If $-\frac d2\leq r<{\gamma}-\frac d2$, we consider
\begin{equation}\label{Vnf} V^{n,r}_f(u,t):= \Delta_n\sum_{i=1}^{[t/\Delta_n]} f\biggl( \frac{\|u(i\Delta_n,\cdot)-u((i-1)\Delta_n,\cdot)\|_{H_r}}{\tau_n(r)} \biggr),\quad t\in[0,\infty),\quad n\geq1, \end{equation}
where $f\colon[0,\infty)\to\mathbb{R}$ is a continuous function satisfying certain regularity properties. For the reasons given in Remarks~\ref{rem-tight} and \ref{difference}, we do not consider functionals as general as \eqref{VnF} when $-\frac d2\leq r<{\gamma}-\frac d2$.

For the case where $r<-\frac d2$, we   recall some facts about Gaussian measures on Hilbert spaces. Let ${\cal L}^{+,\mathrm{sym}}_1(H_r)$ be the space of all bounded linear operators $H_r\to H_r$ that are symmetric, nonnegative (meaning ``nonnegative definite'') and of finite trace (see \cite[Appendix~B]{Prevot07}). By \cite[Theorem~2.1.2]{Prevot07},  for every $Q\in {\cal L}^{+,\mathrm{sym}}_1(H_r)$, there is a Gaussian measure $N_r(0,Q)$ on $H_r$ with mean zero and covariance operator $Q$ (i.e., if $H\sim N_r(0,Q)$, then for every $h\in H_r$, the real-valued random variable $\langle H,h\rangle_{H_r}$ is normally distributed with mean $0$ and variance $\langle Qh,h\rangle_{H_r}$).

In the following, given a functional $F\colon H_r\to \mathbb{R}$ and some $Q\in{\cal L}^{+,\mathrm{sym}}_1(H_r)$, we write $\mu_F(Q) = \mathbb{E}[F(H)]$, where $H\sim N_r(0,Q)$, whenever this expectation exists. Still for $r<-\tfrac d2$, with a slight abuse of notation, given a nonnegative function $w\in L^\infty(D)$ (or a nonnegative random field $w$ with uniformly bounded expectation on $D$), we shall use the abbreviation $\mu_{r,F}(w)=\mu_F(Q_r(w))$ where the operator $Q_r(w)\colon H_r\to H_r$ is defined, for $h\in H_r$, by
\begin{equation}\label{Qphi1}
   Q_r(w) h := \sum_{k=1}^\infty \phi_k \sum_{\ell=1}^\infty {\lambda}_\ell^r a_\ell(h)\int_D \phi_k(y)\phi_\ell(y)w(y)\,\mathrm{d} y,
\end{equation}
or, equivalently, for $h_1\in H_r$ and $h_2\in H_{-r}$, by
\begin{equation}\label{Qphi
(Q_r(w) h_1)(h_2):= \sum_{k,\ell=1}^\infty  {\lambda}_\ell^r a_k(h_2)a_\ell(h_1)\int_D \phi_k(y)\phi_\ell(y)w(y)\,\mathrm{d} y .\end{equation}
One easily checks that $Q_r(w) \in {\cal L}^{+,\mathrm{sym}}_1(H_r)$ since $r <-\tfrac d2$. A more tractable formula for the law $N_r(0,Q_r(w))$ is given in Lemma \ref{eNoQ}.

We first formulate our main result for $r < -\tfrac d2$. Given processes $(X^n(t))_{t\geq0}$ and $(X(t))_{t\geq0}$, we write $X^n\stackrel{L^q}{\Longrightarrow} X$ (or $X^n(t)\stackrel{L^q}{\Longrightarrow} X(t)$) if $\sup_{t\in[0,T]} |X^n(t)-X(t)|\stackrel{L^q}{\longrightarrow}0$  for all $T>0$ as $n\to\infty$.

\begin{theorem}\label{LLN} Let $D$ be a bounded open subset of $\mathbb{R}^d$ with the cone property (see \cite[Section~2]{Agmon65}) and ${\gamma}\in(0,\infty)$. Let $u$ be the solution to \eqref{SHE-2} given by \eqref{akt}. Suppose that $r<-\frac d2$ and let $\tau_n(r)=\sqrt{\Delta_n}$.
Assume that $F\colon H_r \to\mathbb{R}$ is continuous with at most polynomial growth (i.e., there is $p>0$ such that $|F(h)|\lesssim 1+\|h\|_{H_r}^p$ for all $h\in H_r$). If ${\sigma}$ is a predictable random field that is continuous in probability and satisfies
\begin{equation}\label{mom} \sup_{(t,x)\in[0,T]\times D} \mathbb{E}[|{\sigma}(t,x)|^q] < \infty
 \qquad\text{for all } q>0\text{ and } T>0,
\end{equation}
then for every $q>0$, as $n\to\infty$, $V^r_F(u,t) := L^q\text{-}\lim_{n\to\infty} V^{n,r}_F(u,t)$ exists,
\begin{equation}\label{eq:LLN} V^{n,r}_F(u,t) \stackrel{L^q}{\Longrightarrow} V^r_F(u,t)\qquad\text{and}\qquad V^r_F(u,t)= \int_0^t \mu_{r,F} \Bigl( {\sigma}^2(s,\cdot)\Bigr)\,\mathrm{d} s. \end{equation}
\end{theorem}
Theorem \ref{LLN} is proved in Section \ref{sec3.1}.

\begin{remark}\label{L2-rem} By well-known results on uniform integrability (combine Proposition~4.12 with the discussion before Lemma~4.10 in \cite{Kallenberg02}), if ${\sigma}$ is continuous in probability and satisfies \eqref{mom}, then it is automatically $L^p({\Omega})$-continuous for every $p\geq1$. In particular, for every $p\geq1$ and $T>0$,
\begin{equation}\label{L2-cont} w_p({\varepsilon};T):= \sup\Bigl\{(\mathbb{E}[|{\sigma}(t,x)-{\sigma}(s,y)|^p])^{\frac1p}: |t-s|+|x-y|<{\varepsilon},~s,t\in[0,T]\Bigr\} \to 0
\end{equation}
as ${\varepsilon}\to0$.
\end{remark}

Let us apply Theorem~\ref{LLN} to the power variations in \eqref{Vnp} of \emph{even} order, which correspond to the functionals $F_{2p}(h) = \|h\|_{H_r}^{2p}$ for some $p\in\mathbb{N}$. In this case, $\mu_{F_{2p}}$ can be determined explicitly. To this end, we introduce
for a symmetric matrix $A\in\mathbb{R}^{p\times p}$, its \emph{${\alpha}$-permanent} (see \cite{VereJones97}),  defined as
\begin{equation}\label{perm} \mathrm{per}_{\alpha}(A):=\sum_{\Sigma\in\mathrm{Sym}_p} {\alpha}^{\#\Sigma} \prod_{i=1}^p A_{i\Sigma(i)},\qquad {\alpha}\in\mathbb{R}, \end{equation}
where $\mathrm{Sym}_p$ is the symmetric group of all permutations $\Sigma$ on $\{1,\ldots,p\}$ and $\#\Sigma$ is the number of cycles in $\Sigma$. If ${\alpha}=1$, this is the ordinary permanent of $A$, while for ${\alpha}=-1$, this is $(-1)^n$ times the determinant of $A$. For ${\alpha}=\frac12$, this matrix function has an important connection to moments of Gaussian random variables, see  \cite[Equation~(3.23)]{McCullagh18}: If $X_1,\ldots,X_p$ are jointly Gaussian with mean $0$ and covariance matrix $C$, then
\begin{equation}\label{prodmom} \mathbb{E}[X_1^2\cdots X_p^2] = 2^p\mathrm{per}_{\frac12}(C). \end{equation}
In the following, we will apply the permanent function to matrices of the form
\[ D_r(y_1,\ldots,y_p) := \begin{pmatrix} \langle \delta_{y_1},\delta_{y_1}\rangle_{H_r} &\dots& \langle \delta_{y_1},\delta_{y_p}\rangle_{H_r} \\ \vdots & \ddots & \vdots \\ \langle \delta_{y_p},\delta_{y_1}\rangle_{H_r} &\dots& \langle \delta_{y_p},\delta_{y_p}\rangle_{H_r}  \end{pmatrix}\]
and define
\begin{equation}\label{perm2} \mathrm{per}^{(r)}_{\alpha}(y_1,\ldots,y_p) := \mathrm{per}_{{\alpha}}(D_r(y_1,\ldots,y_p)), \end{equation}
where $r<-\frac d2$, $y_1,\dots,y_p\in D$ are given spatial points and $\delta_y$ is the Dirac delta function at $y$. Notice that by \eqref{est1}, $\delta_{y} \in H_r$ for $r < - \frac d2$, and
\[
   \langle \delta_{y_i},\delta_{y_j}\rangle_{H_r} =  \sum_{k=1}^\infty {\lambda}_k^r \phi_k(y_i)\phi_k(y_j).
\]

\begin{corollary}\label{cor1} Assume the hypotheses of Theorem~\ref{LLN} (in particular, $r<-\frac d2$).
\begin{longlist}
\item The process $t\mapsto u(t,\cdot)$ has a locally finite quadratic variation in $H_r$: for any $q>0$,
\begin{equation}\label{2var}\begin{aligned} V^{n,r}_2(u,t)&=\sum_{i=1}^{[t/\Delta_n]} \|u(i\Delta_n,\cdot)-u((i-1)\Delta_n,\cdot)\|^2_{H_r}\\ & \stackrel{L^q}{\Longrightarrow} V^r_2(u,t)=\int_0^t \int_D {\sigma}^2(s,y)\|\delta_y\|_{H_r}^2\,\mathrm{d} y\,\mathrm{d} s. \end{aligned}  \end{equation}
\item More generally, if $p\in\mathbb{N}$, then for any $q>0$, $ V^{n,r}_{2p}(u,t) \stackrel{L^q}{\Longrightarrow}  V^r_{2p}(u,t)$, where
$$
V^{n,r}_{2p}(u,t)=\Delta_n^{1-p} \sum_{i=1}^{[t/\Delta_n]} \|u(i\Delta_n,\cdot)-u((i-1)\Delta_n,\cdot)\|^{2p}_{H_r}
$$
and
\begin{equation}\label{2pvar}
\begin{split}
 V^r_{2p}(u,t)&=\int_0^t \sum_{\Sigma\in\mathrm{Sym}_p} 2^{p-\#\Sigma} \int_D \cdots\int_D \langle\delta_{y_1},\delta_{y_{\Sigma^{-1}(1)}}\rangle_{H_r}\cdots\langle\delta_{y_p},\delta_{y_{\Sigma^{-1}(p)}}\rangle_{H_r}\\
	&\qquad\times {\sigma}^2(s,y_1)\cdots{\sigma}^2(s,y_p)\,\mathrm{d} y_1\cdots\,\mathrm{d} y_p\,\mathrm{d} s\\
	&=2^p\int_0^t \int_D \cdots\int_D \mathrm{per}^{(r)}_{\frac12}(y_1,\ldots,y_p)\,  {\sigma}^2(s,y_1)\cdots{\sigma}^2(s,y_p)\,\mathrm{d} y_1\cdots\,\mathrm{d} y_p \,\mathrm{d} s.
	\end{split}
\end{equation}
\item In particular, in the case of additive noise where ${\sigma}(s,y)\equiv {\sigma}$ for some ${\sigma}\in\mathbb{R}$, we have
\begin{equation}\label{2pvar-add}
 V^r_{2p}(u,t)={\sigma}^{2p} 2^p B_p(x_1,\dots,x_p)t,

\end{equation}
where $B_p$ is the complete Bell polynomial in $p$ variables (see \cite[Definition 2.4.1]{Peccati11}),
\begin{equation}\label{xl}
   x_\ell := \frac{(\ell-1)!}{2} \zeta_D(-\ell r),\qquad \ell=1,\dots,p,
\end{equation}
and $\zeta_D(z):= \sum_{k=1}^\infty {\lambda}_k^{-z}$. In particular, for $p=1$, the quadratic variation in $H_r$ of $t\mapsto u(t,\cdot)$ is $V^r_{2}(u,t) = {\sigma}^2 \zeta_D(-r) t$.
\end{longlist}
\end{corollary}

Corollary \ref{cor1} is proved at the end of Section \ref{sec3.1}.

\begin{remark}\label{rem2.7} The function $\zeta_D(z):= \sum_{k=1}^\infty {\lambda}_k^{-z}$ is known as the  {\em spectral zeta function} or as the \emph{Minakshisundaram--Pleijel zeta function} in the mathematical physics literature (see \cite{Elizalde12,Hawking77,Voros87,Voros92}). By Weyl's law (see \eqref{eq:Weyl}), $\zeta_D(z)$ well defined when $\mathrm{Re}(z)>\frac d2$. In the special case where $d=1$ and $D=[0,\pi]$, then ${\lambda}_k = k^2$ and $\zeta_D(z) = \zeta(2z)$, where $\zeta(\cdot)$ is the classical Riemann zeta function. When ${\sigma}(s,y) \equiv \sigma$, the formula $V^r_{2}(u,t) = {\sigma}^2 \zeta_D(-r) t$ shows that knowing the  quadratic variation of $u$ in $H_r$ for $r<-\tfrac{d}{2}$ yields, via the function $\zeta_D$, information about the eigenvalues $({\lambda}_k)_{k\in\mathbb{N}}$ and hence about the domain $D$.
\end{remark}

Next, in the case where $-\frac d2\leq r<{\gamma}-\frac d2$, we shall only consider functionals of the form $F(h)=f(\|h\|_{H_r})$ with a real-valued function $f\colon [0,\infty)\to\mathbb{R}$; see \eqref{Vnf}.

\begin{theorem}\label{LLN2} Let ${\gamma}\in(0,\infty)$ and $D$ be a bounded connected open subset of $\mathbb{R}^d$ with a piecewise smooth boundary in the sense of \cite[Definition~1.17]{Hezari18} and the cone property (see \cite[Section~2]{Agmon65}). Suppose that $-\frac d2\leq r<{\gamma}-\frac d2$ and let
\begin{equation}\label{taunr} \tau_n(r):= \begin{cases} \Delta_n^{\frac1{2{\gamma}}({\gamma}-\frac{d}{2}-r)}&\text{if } -\frac d2<r<{\gamma}-\frac d2,\\ (\Delta_n\lvert\log\Delta_n\rvert)^{\frac12}&\text{if } r=-\frac d2.\end{cases}\end{equation}
Assume that $f\colon[0,\infty)\to\mathbb{R}$ is continuous with at most polynomial growth (i.e., there is $p>0$ such that $|f(x)|\lesssim 1+x^p$ for all $x\geq0$).

Let $u$ be the solution to \eqref{SHE-2} given by \eqref{akt}, and assume that ${\sigma}$ is a predictable random field that is continuous in probability and satisfies \eqref{mom}.
Then, as $n\to\infty$, $V^r_f(u,t) := L^q\text{-}\lim_{n\to\infty} V^{n,r}_f(u,t)$ exists,
\begin{equation}\label{eq:LLN-2} V^{n,r}_f(u,t) \stackrel{L^q}{\Longrightarrow} V^r_f(u,t)\qquad\text{and}\qquad V^r_f(u,t)= \int_0^t f\Biggl( \sqrt{\frac{K_r}{|D|}\int_D {\sigma}^2(s,y)\,\mathrm{d} y}\Biggr)\,\mathrm{d} s \end{equation}
 for all $q>0$,
where
\begin{equation}\label{Kr} K_r :=
\begin{cases} \frac{|D|}{(4\pi)^{d/2} {\Gamma}(\frac d2)({\gamma}-\frac{d}{2}-r)} {\Gamma}\Bigl(\frac1{\gamma}( r  +\frac{d}{2})\Bigr
	&\text{if } -\frac d2<r<{\gamma}-\frac d2,\\ \frac{|D|}{{\gamma}(4\pi)^{d/2} {\Gamma}(\frac d2)} &\text{if } r=-\frac d2,\end{cases}
\end{equation}
and where $|D|$ denotes the Lebesgue measure of the set $D$.
\end{theorem}

Theorem \ref{LLN2} is proved in Section \ref{nsec3.3}. Using the functions $f(x) = x^p$, $p>0$, we immediately get the following corollary.

\begin{corollary}\label{cor2} Assume the hypotheses of Theorem~\ref{LLN2} and let $-\frac d2< r<{\gamma}-\frac d2$.
\begin{longlist}
	\item The process $t\mapsto u(t,\cdot)$ has a locally finite $\frac{2{\gamma}}{{\gamma}- d/2 -r}$-variation in $H_r$. More precisely, for any $q>0$,
	\begin{equation}\label{exact-var}\begin{aligned} V^{n,r}_{\frac{2{\gamma}}{{\gamma}- d/2 -r}}(u,t)&=\sum_{i=1}^{[t/\Delta_n]} \|u(i\Delta_n,\cdot)-u((i-1)\Delta_n,\cdot)\|^{\frac{2{\gamma}}{{\gamma}- d/2 -r}}_{H_r}\\ & \stackrel{L^q}{\Longrightarrow} V^r_{\frac{2{\gamma}}{{\gamma}- d/2 -r}}(u,t)=\biggl( \frac{K_r}{|D|}\biggr)^{\frac{{\gamma}}{{\gamma}- d/2 -r}}\int_0^t \biggl( \int_D {\sigma}^2(s,y)\,\mathrm{d} y\biggr)^{\frac{{\gamma}}{{\gamma}- d/2 -r}}\,\mathrm{d} s. \end{aligned}  \end{equation}
	\item More generally, for any $p>0$,
	\begin{equation}\label{p-var}\begin{aligned} V^{n,r}_p(u,t)&=\Delta_n^{1-\frac p{2{\gamma}}({\gamma}-\frac{d}{2}-r)}\sum_{i=1}^{[t/\Delta_n]} \|u(i\Delta_n,\cdot)-u((i-1)\Delta_n,\cdot)\|^p_{H_r}\\ & \stackrel{L^q}{\Longrightarrow} V^r_p(u,t)=\biggl( \frac{K_r}{|D|}\biggr)^{\frac p2}\int_0^t \biggl( \int_D{\sigma}^2(s,y)\,\mathrm{d} y\biggr)^{\frac p2}\,\mathrm{d} s. \end{aligned}  \end{equation}
	\item In particular, in the case of additive noise where ${\sigma}(s,y)\equiv {\sigma}$ for some ${\sigma}\in\mathbb{R}$, we have
		\begin{equation}\label{p-var-add}V^r_p(u,t)=|{\sigma}|^p K_r^{\frac p2}t. \end{equation}
	\item Part (ii) and part (iii) remain valid for $r=-\frac d2$ when we replace the first line in \eqref{p-var} by
	\begin{equation}\label{p-var-12} V^{n,r}_p(u,t)=\frac{\Delta_n^{1-\frac{p}{2}}}{\lvert\log \Delta_n\rvert^{\frac p2}}\sum_{i=1}^{[t/\Delta_n]} \|u(i\Delta_n,\cdot)-u((i-1)\Delta_n,\cdot)\|^p_{H_r}.  \end{equation}
\end{longlist}
\end{corollary}

\begin{remark}\label{rem2.10} Contrary to the case $r < -\tfrac{d}{2}$, in the case where $-\tfrac{d}{2} \leq r < \gamma - \tfrac{d}{2}$, the quadratic variation of $u$ is infinite and the $\frac{2{\gamma}}{{\gamma}- d/2 -r}$-variation, which is finite, contains very little information about the domain $D$, other than its volume $|D|$, via the constant $K_r$.
\end{remark}
Our proofs also yield the optimal Hölder regularity of the sample paths of $u$ (as a process with values in $H_r$); cf.\ also \cite[Theorem~11.8]{Peszat07} and \cite{vanNeerven08}. Similar results in the hyperbolic case (with $D=\mathbb{R}^d$) were obtained in \cite{Dalang05}. The proof of the following Corollary \ref{Holder} is given at the end of Section~\ref{nsec3.4}.
\begin{corollary}\label{Holder} Let $r<{\gamma}-\frac d2$ and
assume the hypotheses of Theorem~\ref{LLN} (if $r<-\frac d2$) or Theorem~\ref{LLN2} (if $r\geq -\frac d2$). Then $t\mapsto u(t,\cdot)$, viewed as a stochastic process with values in $H_r$ has a version that is almost surely locally ${\alpha}$-Hölder continuous for all ${\alpha}<{\alpha}(r)$,
where ${\alpha}(r):=\frac12$ if $ r\leq -\frac d2$ and ${\alpha}(r):=\frac{1}{2{\gamma}}({\gamma}-\frac d2 -r)$ if $-\frac d2<r<{\gamma}-\frac d2$.
Moreover, if ${\sigma}\equiv1$ and ${\alpha}>{\alpha}(r)$, then $t\mapsto u(t,\cdot)$ does not have a version that is almost surely locally ${\alpha}$-Hölder continuous.
\end{corollary}

The next section contains the proofs of all the results that we have stated above. In Section~\ref{compare}, we also explain in Remarks \ref{rem-tight} and \ref{difference} why in Theorem~\ref{LLN2} (where $-\frac d2\leq r<{\gamma}-\frac d2$), we did not study functionals as general as in Theorem~\ref{LLN} (where $r<-\frac d2$). In the \hyperref[appn]{Appendix}, we gather several important properties and estimates concerning eigenvalues and eigenfunctions of the Laplacian on $D$, which are used in these proofs.

\section{Proofs}

We first show that if ${\sigma}\equiv 1$, then the $H_r$-norm of an increment $u(i\Delta_n,\cdot)-u((i-1)\Delta_n,\cdot)$ is typically of order $\tau_n(r)$ as defined in Theorem \ref{LLN} or \eqref{taunr} (depending on the value of $r$).
\begin{lemma}\label{size-incr} Let
\begin{equation}\label{Kr-2} K_r := \sum_{k=1}^\infty {\lambda}_k^r = \zeta_D(-r), \qquad\text{if } r<-\frac d2, \end{equation}
and let $K_r$ be defined as in \eqref{Kr} if $-\frac d2\leq r<{\gamma}-\frac d2$. If $u$ is the solution to \eqref{SHE-2} with ${\sigma}\equiv1$, then for every ${\varepsilon}>0$,
\begin{equation}\label{limit} \lim_{n\to\infty} \sup_{i\colon i\Delta_n\geq{\varepsilon}} \Biggl\lvert\frac{\mathbb{E}[\|u(i\Delta_n,\cdot)-u((i-1)\Delta_n,\cdot)\|_{H_r}^2]}{\tau_n(r)^2} -K_r\Biggr\rvert = 0.\end{equation}
\end{lemma}
\begin{proof}  By standard calculations, and using the notation $\iint_a^b = \int_a^b \int_D$,
\begin{equation}\label{calc}
\begin{aligned}
&\mathbb{E}[\|u(i\Delta_n,\cdot)-u((i-1)\Delta_n,\cdot)\|_{H_r}^2]\\
&\quad= \sum_{k=1}^\infty {\lambda}_k^{r}\mathbb{E}[(a_k(i\Delta_n)-a_k((i-1)\Delta_n))^2] \\
&\quad=\sum_{k=1}^\infty {\lambda}_k^{r} \Biggl( \iint_0^{i\Delta_n} \mathrm{e}^{-2{\lambda}_k^{\gamma} (i\Delta_n-s)}\phi_k(y)^2\,\mathrm{d} s\,\mathrm{d} y+\iint_0^{(i-1)\Delta_n} \mathrm{e}^{-2{\lambda}_k^{\gamma}((i-1)\Delta_n-s)}\phi_k(y)^2\,\mathrm{d} s\,\mathrm{d} y\\
&\quad\qquad-2\iint_0^{(i-1)\Delta_n} \mathrm{e}^{-{\lambda}_k^{\gamma}(i\Delta_n-s)-{\lambda}_k^{\gamma}((i-1)\Delta_n-s)}\phi_k(y)^2\,\mathrm{d} s\,\mathrm{d} y \Biggr)\\
&\quad=\sum_{k=1}^\infty {\lambda}_k^{r}\Biggl( \int_0^{i\Delta_n} \mathrm{e}^{-2{\lambda}_k^{\gamma} s}\,\mathrm{d} s -2\int_0^{(i-1)\Delta_n} \mathrm{e}^{-{\lambda}_k^{\gamma} (s+\Delta_n)-{\lambda}_k^{\gamma} s}\,\mathrm{d} s+\int_0^{(i-1)\Delta_n} \mathrm{e}^{-2{\lambda}_k^{\gamma} s}\,\mathrm{d} s   \Biggr)\\

\quad= \displaystyle \sum_{k=1}^\infty {\lambda}_k^{r} \frac{1-\mathrm{e}^{-2{\lambda}_k^{\gamma} i\Delta_n}-2\mathrm{e}^{-{\lambda}_k^{\gamma} \Delta_n}+2\mathrm{e}^{-{\lambda}_k^{\gamma} \Delta_n}\mathrm{e}^{-2{\lambda}_k^{\gamma} (i-1)\Delta_n}+1-\mathrm{e}^{-2{\lambda}_k^{\gamma} (i-1)\Delta_n}}{2{\lambda}_k^{\gamma}}\\
&\quad=\sum_{k=1}^\infty {\lambda}_k^{r-{\gamma}}(1-\mathrm{e}^{-{\lambda}_k^{\gamma}\Delta_n}) -\frac12 \sum_{k=1}^\infty {\lambda}_k^{r-{\gamma}} \mathrm{e}^{-2{\lambda}_k^{\gamma} i\Delta_n}  (\mathrm{e}^{{\lambda}_k^{\gamma}\Delta_n}-1)^2.
\end{aligned}\end{equation}
Let us denote the first and the second series by $A^{n,1}_r$ and $A^{n,i,2}_r$, respectively. By dominated convergence, we have
\begin{equation}\label{second-term} \limsup_{n\to\infty} \sup_{i\colon i\Delta_n\geq{\varepsilon}} \Delta_n^{-2}A^{n,i,2}_r \leq \lim_{n\to\infty}  \sum_{k=1}^\infty \frac{{\lambda}_k^{r-{\gamma}}}{2} \mathrm{e}^{-2{\lambda}_k^{\gamma} {\varepsilon}}  \Biggl(\frac{\mathrm{e}^{{\lambda}_k^{\gamma}\Delta_n}-1}{\Delta_n} \Biggr)^2 = \sum_{k=1}^\infty \frac{{\lambda}_k^{r+{\gamma}}}{2}\mathrm{e}^{-2{\lambda}_k^{\gamma} {\varepsilon}}<\infty, \end{equation}
which shows that the contribution of $A^{n,i,2}_r$ is negligible for the limit taken in \eqref{limit}.

Concerning $A^{n,1}_r$, if $r<-\frac d2$, then we immediately deduce from dominated convergence that
\begin{equation}\label{convKr} \lim_{n\to\infty} \Delta_n^{-1}A^{n,1}_r = \sum_{k=1}^\infty {\lambda}_k^r = K_r,  \end{equation}
which implies \eqref{limit}.
\begin{comment}
If $-\frac d2<r<{\gamma}-\frac d2$, Lemma~\ref{A1} (i) shows that for every $\delta\in(0,1)$, there is $N\in \mathbb{N}$ such that $(1-\delta)C_D k^{2/d}\leq {\lambda}_k\leq (1+\delta)C_D k^{2/d}$ for all $k\geq N$, where $C_D$ is the constant on the right-hand side of \eqref{eq:Weyl}. Because $\Delta_n^{-1}\sum_{k=1}^N {\lambda}_k^{r-{\gamma}}(1-\mathrm{e}^{-{\lambda}_k^{\gamma} \Delta_n}) \to \sum_{k=1}^N {\lambda}_k^r$ by the same argument as in \eqref{convKr} and $\Delta_n/\tau_n(r)^2 \to 0$, the term $A^{n,1}_r/\tau_n(r)^2$ has the same limiting behavior as $A^{n,1}_r(N)/\tau_n(r)^2$ when $n\to\infty$, where
\[ A^{n,1}_r(N)=\sum_{k=N+1}^\infty {\lambda}_k^{r-{\gamma}}(1-\mathrm{e}^{-{\lambda}_k^{\gamma} \Delta_n}). \]
By our choice of $N$ and \eqref{taunr}, we have
\begin{equation}\label{help}\begin{split} \frac{A^{n,1}_r(N)}{\tau_n(r)^2} &\leq ((1-\delta)C_D)^{r-{\gamma}}\tau_n(r)^{-2}\sum_{k=N+1}^\infty k^{\frac 2d(r-{\gamma})}(1-\mathrm{e}^{-((1+\delta)C_Dk^{2/d})^{\gamma}\Delta_n})\\
&=((1-\delta)C_D)^{r-{\gamma}}\sum_{k=N+1}^\infty(k\Delta_n^{\frac{d}{2{\gamma}}})^{\frac 2d (r-{\gamma})} (1-\mathrm{e}^{-((1+\delta)C_D)^{\gamma} (k\Delta_n^{d/(2{\gamma})})^{2{\gamma}/d}})\Delta_n^{\frac{d}{2{\gamma}}}
\end{split}\end{equation}
Since $\Delta_n^{d/(2{\gamma})}\to 0$,
the last sum is in fact a Riemann sum. Thus, as we shall explain below,
\begin{equation}\label{riemann} \limsup_{n\to\infty} \frac{A^{n,1}_r}{\tau_n(r)^2}=\limsup_{n\to\infty} \frac{A^{n,1}_r(N)}{\tau_n(r)^2} \leq ((1-\delta)C_D)^{r-{\gamma}} \int_0^\infty z^{\frac2d(r-{\gamma})}(1-\mathrm{e}^{-((1+\delta)C_D)^{\gamma} z^{2{\gamma}/d}})\,\mathrm{d} z.
\end{equation}
Substituting $((1+\delta)C_D)^{\gamma} z^{2{\gamma}/d}\mapsto y$, integrating by parts, and employing the relation $z{\Gamma}(z) = {\Gamma}(z+1)$ for $z\notin\{0,-1,-2,\dots\}$, the integral on the right-hand side of \eqref{riemann} is equal to
\begin{align*}
&\frac{(1-\delta)^{r-{\gamma}}}{(1+\delta)^{r-{\gamma}+d/2}}\,\frac{d}{2{\gamma} C_D^{d/2}} \int_0^\infty y^{\frac r{\gamma} -2+\frac{d}{2{\gamma}}}(1-\mathrm{e}^{-y})\,\mathrm{d} y\\
&\quad =-\frac{(1-\delta)^{r-{\gamma}}}{(1+\delta)^{r-{\gamma}+d/2}}\frac{d}{2{\gamma} C_D^{d/2}} \,\frac{1}{\frac r {\gamma}-1+\frac{d}{2{\gamma}}} \int_0^\infty y^{\frac r {\gamma} -1+\frac{d}{2{\gamma}}} \mathrm{e}^{-y}\, \mathrm{d} y\\
&\quad= -\frac{(1-\delta)^{r-{\gamma}}}{(1+\delta)^{r-{\gamma}+d/2}}\,\frac{d}{2{\gamma} C_D^{d/2}} \,\frac{1}{\frac r {\gamma}-1+\frac{d}{2{\gamma}}} {\Gamma}\Bigl(\tfrac r {\gamma} +\tfrac{d}{2{\gamma}}\Bigr) = \frac{(1-\delta)^{r-{\gamma}}}{(1+\delta)^{r-{\gamma}+d/2}} K_r. \end{align*}
Replacing $\delta$ by $-\delta$ in the above argument, one proves in an analogous fashion that
\[ \liminf_{n\to\infty} \frac{A^{n,1}_r}{\tau_n^2(r)} = \liminf_{n\to\infty} \frac{A^{n,1}_r(N)}{\tau_n^2(r)}\geq \frac{(1+\delta)^{r-{\gamma}}}{(1-\delta)^{r-{\gamma}+d/2}}K_r. \]
Sending $\delta\to0$ completes the proof of Lemma~\ref{size-incr} for $-\frac d2< r<{\gamma}-\frac d2$.
We note that \eqref{riemann} requires an explanation. To this end, we first write the last sum in \eqref{help} as
\begin{equation}\label{help2} \int_0^\infty \sum_{k=N+1}^\infty (k\Delta_n^{\frac{d}{2{\gamma}}})^{\frac 2d (r-{\gamma})} (1-\mathrm{e}^{-((1+\delta)C_D)^{\gamma} (k\Delta_n^{d/(2{\gamma})})^{2{\gamma}/d}})\mathds 1_{(k\Delta_n^{d/(2{\gamma})},(k+1)\Delta_n^{d/(2{\gamma})}]}(z)\,\mathrm{d} z \end{equation}
and notice that the integrand converges pointwise in $z$ to the integrand in \eqref{riemann} as $\Delta_n\to0$ (because $N$ is fixed). Furthermore, for any $a,b,c>0$ with $a<b$, there is a unique $z_0=z_0(a,b,c)\in(0,\infty)$ such that the function $z\mapsto h(z)=z^{a-b}(1-\mathrm{e}^{-cz^b})$ is increasing on $(0,z_0)$ and decreasing on $(z_0,\infty)$. Therefore, if $a=\frac2d r$, $b=\frac2d {\gamma}$, and $c=((1+\delta)C_D)^{\gamma}$, then the integrand in \eqref{help2} is uniformly bounded in $n$ by
 \begin{equation}\label{help3} h(z)\mathds 1_{(0,z_0]}(z)+h(z_0)\mathds 1_{(z_0,z_0+{\varepsilon}]}(z)+h(z-{\varepsilon})\mathds 1_{(z_0+{\varepsilon},\infty)}(z)\end{equation}
for any ${\varepsilon}>\Delta_1$. For our choices of $a$, $b$, and $c$, the function in \eqref{help3} is integrable on $(0,\infty)$, so \eqref{riemann} follows by dominated convergence.
 %
 %
It remains to consider the term $A^{n,1}_r/\tau_n(r)^2$ for $r=-\frac d2$. It can be decomposed as
\[ \frac{A^{n,1}_{-d/2}}{\tau_n(-\frac d2)^2}= \frac{1}{\Delta_n\lvert\log\Delta_n\rvert} \sum_{k=1}^{[\Delta_n^{-d/(2{\gamma})}]} \frac{1-\mathrm{e}^{-{\lambda}_k^{\gamma}\Delta_n}}{{\lambda}_k^{d/2+{\gamma}}} + \frac{1}{\Delta_n\lvert\log\Delta_n\rvert} \sum_{k=[\Delta_n^{-d/(2{\gamma})}]+1}^\infty \frac{1-\mathrm{e}^{-{\lambda}_k^{\gamma}\Delta_n}}{{\lambda}_k^{d/2+{\gamma}}}.   \]
Let us denote the two terms on the right-hand side by $B^{n,1}$ and $B^{n,2}$, respectively. Clearly, by \eqref{eq:Weyl} and integral approximation,
\[ |B^{n,2}| \lesssim \frac{1}{\Delta_n\lvert\log\Delta_n\rvert} \sum_{k=[\Delta_n^{-d/(2{\gamma})}]+1}^\infty \frac{1}{k^{1+2{\gamma}/d}} \leq \frac{1}{\Delta_n \lvert\log\Delta_n\rvert} \int_{[\Delta_n^{-d/(2{\gamma})}]}^\infty x^{-1-2{\gamma}/d}\,\mathrm{d} x \lesssim \frac{1}{\lvert\log \Delta_n\rvert} \to 0 \]
as $n\to\infty$.

Next, the term $B^{n,1}$ can further be written as $B^{n,11}+B^{n,12}$ where
\begin{align*} B^{n,11}&= \frac{1}{\Delta_n\lvert\log\Delta_n\rvert} \sum_{k=1}^{[\Delta_n^{-d/(2{\gamma})}]} \frac{{\lambda}_k^{\gamma}\Delta_n}{{\lambda}_k^{d/2+{\gamma}}} = \frac{1}{\lvert\log\Delta_n\rvert} \sum_{k=1}^{[\Delta_n^{-d/(2{\gamma})}]} \frac{1}{{\lambda}_k^{d/2}},\\
 B^{n,12}&=\frac{1}{\Delta_n\lvert\log\Delta_n\rvert} \sum_{k=1}^{[\Delta_n^{-d/(2{\gamma})}]} \frac{1-\mathrm{e}^{-{\lambda}_k^{\gamma}\Delta_n}-{\lambda}_k^{\gamma}\Delta_n}{{\lambda}_k^{d/2+{\gamma}}}. \end{align*}
 Using \eqref{eq:Weyl} and the fact that $(\log n)^{-1}\sum_{k=1}^n k^{-1} \to 1$ as $n\to\infty$, it is easy to verify that $B^{n,11} \to\frac{|D|}{{\gamma}(4\pi)^{d/2} {\Gamma}(\frac d2)}$. Moreover, as $|1-\mathrm{e}^{-x}-x|\leq \frac{x^2}{2}$ for all $x\geq0$, \eqref{eq:Weyl} leads to
 \[ |B^{n,12}|\leq \frac{\Delta_n}{2\lvert\log\Delta_n\rvert}\sum_{k=1}^{[\Delta_n^{-d/(2{\gamma})}]} {\lambda}_k^{{\gamma}-\frac d2} \lesssim \frac{1}{\lvert\log\Delta_n\rvert}\to0. \]
 Altogether we obtain $A^{n,1}_r/\tau_n(r)^2\to\frac{|D|}{{\gamma}(4\pi)^{d/2} {\Gamma}(\frac d2)} = K_r$ for $r=-\frac d2$, which completes the proof of \eqref{limit}.
 \end{comment}
 If $-\frac d2\leq r<{\gamma}-\frac d2$, then Lemma~\ref{convdens} (i) shows $\tau(n)^{-2}A^{n,1}_r=H_n\to K_r$.
\end{proof}

\subsection{Proof of Theorem~\ref{LLN} and Corollary~\ref{cor1} ($r < -\tfrac d2$)}\label{sec3.1}

\begin{proof}[Proof of Theorem~\ref{LLN}]
We only show \eqref{eq:LLN} for $q=1$. By \cite[Proposition~4.12]{Kallenberg02} and the discussion before Lemma~4.10 in \cite{Kallenberg02}, the statement for all larger values of $q$ follows easily from the hypothesis that $F$ has polynomial growth and the fact that $u$ has locally uniformly bounded moments of all orders, by using the fact that convergence in probability together with bounded $L^{q+{\varepsilon}}$-moments implies convergence in $L^q$.

Since $V^{n,r}_F$ and $V^r_F$ in \eqref{eq:LLN} are linear in $F$, by decomposing $F$ into its positive and negative parts, we may further assume without loss of generality that $F$ is nonnegative. In this case, both $V^{n,r}_F(u,t)$ and $V^r_F(u,t)$ are nonnegative and increasing in $t$, so uniform convergence on compacts is equivalent to pointwise convergence in $t$ (see \cite[Chapter VI, Theorem 2.15 c)]{Jacod03}). Hence, in the following, we consider a fixed time point $t>0$. Furthermore, for brevity, we use the notation
\begin{equation}\label{abbr}\begin{aligned}
	\Delta^n_i u &:= u(i\Delta_n,\cdot)-u((i-1)\Delta_n,\cdot),\quad \Delta^n_i a_k := a_k(i\Delta_n)-a_k((i-1)\Delta_n),\\
	\Delta^n_i g(s;x,y)&:=g(i\Delta_n-s;x,y)-g((i-1)\Delta_n-s;x,y)\mathds 1_{[0,(i-1)\Delta_n]}(s),\\
	\Delta^n_i \mathrm{e}_k(s)&:=\mathrm{e}^{-{\lambda}_k^{\gamma}(i\Delta_n-s)}\mathds 1_{[0,i\Delta_n]}(s)-\mathrm{e}^{-{\lambda}_k^{\gamma}((i-1)\Delta_n-s)}\mathds 1_{[0,(i-1)\Delta_n]}(s).\end{aligned}
\end{equation}

Recall from \eqref{akt} that the coefficients $a_k(t)$ of $u(t,\cdot)$ are in fact Ornstein--Uhlenbeck-type processes,
that is, for every $k\in\mathbb{N}$, $a_k$ satisfies the stochastic differential equation
\begin{equation}\label{SDE} \mathrm{d} a_k(t)= -{\lambda}_k^{\gamma} a_k(t)\,\mathrm{d} t + \mathrm{d} X_k(t),\qquad a_k(0)=0,\end{equation}
where
\begin{equation}\label{Xkt} X_k(t):=\int_0^t\int_D \phi_k(y){\sigma}(s,y)\,W(\mathrm{d} s,\mathrm{d} y).\end{equation}
As a first step, we will show that the drift part is asymptotically negligible, leading us to approximate $u(t,\cdot)$ by $u'(t,\cdot)$, where
\begin{equation}\label{uprime}
   u'(t,\cdot):= \sum_{k=1}^\infty \phi_k X_k(t),\qquad t\in[0,\infty),
   \end{equation}
and the limit in \eqref{uprime} is taken in $H_r$.

As a second step, we locally freeze the coefficient $\sigma$, leading us to approximate $V^{n,r}_F(u',t)$ by
$$
\Delta_n\sum_{i=1}^{[t/\Delta_n]} F\biggl( \frac{h^n_i}{\sqrt{\Delta_n}}\biggr),
$$
where, for $n\in\mathbb{N}$ and $i=1,\ldots,[t/\Delta_n]$, $h^n_i$ is the $H_r$-valued random variable
\begin{equation}\label{hni} h^n_i := \sum_{k=1}^\infty \phi_k\iint_{(i-1)\Delta_n}^{i\Delta_n}\phi_k(y){\sigma}((i-1)\Delta_n,y)\,W(\mathrm{d} s,\mathrm{d} y).  \end{equation}
With this notation, we will prove that \eqref{eq:LLN} holds by writing
\begin{equation}\label{eq3.33}
   V^{n,r}_F(u,t) - \int_0^t \mu_F( {\sigma}^2(s,\cdot))\,\mathrm{d} s = I_{n,1} + \cdots + I_{n,4},
\end{equation}
where
\begin{align*}
   I_{n,1} &:= V^{n,r}_F(u,t)-V^{n,r}_F(u',t),
  \qquad I_{n,2} := V^{n,r}_F(u',t) - \Delta_n\sum_{i=1}^{[t/\Delta_n]} F\biggl( \frac{h^n_i}{\sqrt{\Delta_n}}\biggr),\\
   I_{n,3} &:= \Delta_n\sum_{i=1}^{[t/\Delta_n]} F\biggl( \frac{h^n_i}{\sqrt{\Delta_n}}\biggr) - \Delta_n\sum_{i=1}^{[t/\Delta_n]} \mu_F(Q^{n,i}_r),\\
  I_{n,4} &:= \Delta_n\sum_{i=1}^{[t/\Delta_n]} \mu_F(Q^{n,i}_r) - \int_0^t \mu_F( {\sigma}^2(s,\cdot))\,\mathrm{d} s,
\end{align*}
and the operator $Q^{n,i}_r$ is defined in \eqref{Qni} below. The terms $I_{n,1}$ and $I_{n,2}$ are exactly the errors incurred by the two approximations described above. Their convergence to $0$ in $L^1({\Omega})$ will be proved in Lemmas \ref{remove-drift} and \ref{discr-sigma}, respectively, based on a tightness argument (see Lemma~\ref{tight}). Next, we establish the law of large numbers for $\Delta_n\sum_{i=1}^{[t/\Delta_n]} F( {h^n_i}/{\sqrt{\Delta_n}})$ in Lemma~\ref{LLN-core-1}, which, together with a calculation of the conditional expectation that appears there (see Lemma~\ref{comp-condexp}), shows that $\mathbb{E}[|I_{n,3}|]\to0$ as $n\to\infty$. Finally, Lemma~\ref{cex-conv} shows that the discrete sum that we obtain from the previous step approximates $V^r_F(u,t)=\int_0^t\mu_F({\sigma}^2(s,\cdot))\,\mathrm{d} s$ in $L^1(\Omega)$, that is, $\mathbb{E}[|I_{n,4}|]\to 0$. To summarize, once Lemmas~\ref{remove-drift}--\ref{cex-conv} are proved, the proof of Theorem \ref{LLN} will be complete.
\end{proof}

\begin{lemma}\label{eNoQ}
Let $r < -\tfrac d2$. If $H$ is a random vector with values in $H_r$ and law $N_r(0,Q_r(w))$, where $Q_r(w)$ is defined in \eqref{Qphi1} and \eqref{Qphi}, then $H$ has the same law as $\sum_{k=1}^\infty X_k\, {\lambda}_k^{- r/2} \phi_k$, where the $X_k$ are jointly Gaussian centered random variables with covariances
\begin{equation}\label{covkl}
{\mathrm{Cov}}(X_k,X_\ell)  =  {\lambda}_k^{\frac r2}{\lambda}_\ell^{\frac r2} \int_D \phi_k(y)\phi_\ell(y)w(y)\,\mathrm{d} y.
\end{equation}
\end{lemma}
\begin{proof}
Consider the orthonormal basis of $H_r$ given by
\begin{equation}\label{basis_b}
\Bigl(b_k:={\lambda}_k^{-\frac r2}\phi_k: k\in\mathbb{N}\Bigr)
\end{equation}
 and note that $H =\sum_{k=1}^\infty a_k(H) {\lambda}_k^{r/2}b_k$. Define $X_k = {\lambda}_k^{r/2}a_k(H)$. Since $a_k(b_\ell) = {\lambda}_k^{-r/2}\delta_{k,\ell}$, we have $X_k = \langle H,b_k\rangle_{H_r}$, and by \cite[Theorem 2.1.2]{Prevot07},
\[
{\mathrm{Cov}}(X_k,X_\ell) = \mathbb{E}[\langle H,b_k\rangle_{H_r}\, \langle H,b_\ell\rangle_{H_r}]= \langle Q_r(w)b_k, b_\ell \rangle_{H_r}.
\]
By \eqref{Qphi1},
\begin{align*}
Q_r(w)b_k &=  \sum_{k_1=1}^\infty \phi_{k_1}\, {\lambda}_k^r a_k(b_k)\int_D \phi_{k_1}(y)\phi_k(y)w(y)\,\mathrm{d} y\\
&=  \sum_{k_1=1}^\infty b_{k_1} \, {\lambda}_{k_1}^{\frac r2} \, {\lambda}_k^{\frac r2} \int_D \phi_{k_1}(y)\phi_k(y)w(y)\,\mathrm{d} y.
\end{align*}
Therefore,
\begin{align*}
\langle Q_r(w)b_k, b_\ell \rangle_{H_r}
& =   {\lambda}_\ell^{\frac r2}  {\lambda}_k^{\frac r2}   \int_D \phi_{\ell}(y)\phi_k(y)w(y)\,\mathrm{d} y.
\end{align*}
This completes the proof.
\end{proof}

\begin{proof}[Proof of Corollary~\ref{cor1}] We first prove part (ii). The $2p$-variation obviously corresponds to the functional $F\colon H_r\to \mathbb{R}$ given by $F(h)=\|h\|_{H_r}^{2p}$, which satisfies the hypotheses of Theorem~\ref{LLN}. For a deterministic function $w \in L^2(D)$, let $H$ be a random vector in $H_r$ with law $N_r(0,Q_r(w))$. By definition and Lemma \ref{eNoQ},
\begin{equation}\label{mu2p}
\begin{aligned}\mu_{r,F}(w) &= \mathbb{E}[F(H)] = \mathbb{E}[\|H\|_{H_r}^{2p}] = \mathbb{E}\Biggl[ \Biggl(\sum_{k=1}^\infty X_k^2 \Biggl)^p \Biggr]
	=\sum_{k_1,\ldots,k_p=1}^\infty \mathbb{E}\Big[ X_{k_1}^2\cdots X_{k_p}^2\Big],
\end{aligned}
\end{equation}
where
$(X_k)_{k\in\mathbb{N}}$ is a sequence of jointly Gaussian centered random variables with covariances  given by \eqref{covkl}.
Using \eqref{prodmom} to evaluate the expectations in the last line of \eqref{mu2p}, we obtain
\begin{align*}
\mu_{r,F}(w) &= 2^p\sum_{k_1,\ldots,k_p=1}^\infty \mathrm{per}_{\frac12}\Bigl({\mathrm{Cov}}(X_{k_i},X_{k_j})_{i,j=1}^p\Bigr)\\
&=\sum_{\Sigma\in\mathrm{Sym}_p} 2^{p-\#\Sigma} \int_D \cdots\int_D \sum_{k_1,\ldots,k_p=1}^\infty {\lambda}_{k_1}^{\frac r2}\, {\lambda}^{\frac r 2}_{k_{\Sigma(1)}}\cdots {\lambda}_{k_p}^{\frac r2} \, {\lambda}_{k_{\Sigma(p)}}^{\frac r2}\\
&\qquad\times \phi_{k_1}(y_1)\phi_{k_{\Sigma(1)}}(y_1)\cdots \phi_{k_p}(y_p) \phi_{k_{\Sigma(p)}}(y_p)  w(y_1)\cdots w(y_p)\,\mathrm{d} y_1\cdots\,\mathrm{d} y_p.
\end{align*}
Since $\Sigma$ is one-to-one, each factor ${\lambda}_{k_i}^{r/2}$, for $i=1,\ldots,p$, appears exactly twice, and the same holds for the function $\phi_{k_i}$, but once with argument $y_i$ and once with argument $y_{\Sigma^{-1}(i)}$. Thus,
\begin{equation}\label{muFsigma}\begin{aligned}
	\mu_{r,F}({\sigma}^2(s,\cdot)) &= \sum_{\Sigma\in\mathrm{Sym}_p} 2^{p-\#\Sigma} \int_D \cdots\int_D \sum_{k_1=1}^\infty {\lambda}_{k_1}^{r} \phi_{k_1}(y_1)\phi_{k_1}(y_{\Sigma^{-1}(1)}) \\
	&\qquad\times\cdots\times \sum_{k_p=1}^\infty {\lambda}_{k_p}^{r} \phi_{k_p}(y_p)\phi_{k_p}(y_{\Sigma^{-1}(p)})  {\sigma}^2(s,y_1)\cdots{\sigma}^2(s,y_p)\,\mathrm{d} y_1\cdots\,\mathrm{d} y_p\\
	&=\sum_{\Sigma\in\mathrm{Sym}_p} 2^{p-\#\Sigma} \int_D \cdots\int_D \langle\delta_{y_1},\delta_{y_{\Sigma^{-1}(1)}}\rangle_{H_r}\cdots\langle\delta_{y_p},\delta_{y_{\Sigma^{-1}(p)}}\rangle_{H_r}\\
	&\qquad\times {\sigma}^2(s,y_1)\cdots{\sigma}^2(s,y_p)\,\mathrm{d} y_1\cdots\,\mathrm{d} y_p,
\end{aligned}\end{equation}
and part (ii) of Corollary~\ref{cor1} follows from Theorem~\ref{LLN}.

The other two parts of the corollary are special cases of the result that we have just proved. Indeed, for $p=1$, the permanent of a number $A\in\mathbb{R}$ is simply $\mathrm{per}_{\alpha}(A)={\alpha} A$. And because $D_r(y)=\|\delta_y\|_{H_r}^2$, part (i) follows.

Regarding part (iii), if ${\sigma}(s,y)\equiv {\sigma}$, then \eqref{muFsigma} becomes
\begin{equation}\label{3.37}\begin{aligned}
	\mu_{r,F}({\sigma}^2) &= {\sigma}^{2p}\sum_{\Sigma\in\mathrm{Sym}_p} 2^{p-\#\Sigma} \int_D \mathrm{d} y_1 \cdots\int_D \mathrm{d} y_p \Biggl[\sum_{k_1=1}^\infty {\lambda}_{k_1}^{r} \phi_{k_1}(y_1)\phi_{k_1}(y_{\Sigma^{-1}(1)}) \Biggr] \\
	&\qquad \times\cdots\times \Biggl[\sum_{k_p=1}^\infty {\lambda}_{k_p}^{r} \phi_{k_p}(y_p)\phi_{k_p}(y_{\Sigma^{-1}(p)}) \Biggr] \\
	&= {\sigma}^{2p}\sum_{\Sigma\in\mathrm{Sym}_p} 2^{p-\#\Sigma} \sum_{k_1=1}^\infty \cdots \sum_{k_p=1}^\infty  {\lambda}_{k_1}^{r} \cdots {\lambda}_{k_p}^{r}  \int_D \mathrm{d} y_1 \cdots\int_D \mathrm{d} y_p \, \phi_{k_1}(y_1)\phi_{k_1}(y_{\Sigma^{-1}(1)}) \\
	&\qquad  \times \cdots \times \phi_{k_p}(y_p)\phi_{k_p}(y_{\Sigma^{-1}(p)}).
\end{aligned}\end{equation}

Suppose that $\#\Sigma = \ell$ and these $\ell$ cycles have respective lengths $q_1,\dots, q_\ell$, with $q_1 + \cdots + q_\ell = p$. From the orthogonality property of the functions $\phi_{k_i}$, we see that the integral in \eqref{3.37} will vanish unless
\[
k_1 = k_{\Sigma(1)} = k_{\Sigma(\Sigma(1))} = \cdots,
\]
that is, the indices $k_i$ in each cycle of $\Sigma$ coincide.

The term in \eqref{3.37} corresponding to such a $\Sigma$ is
\begin{equation}\label{3.38}
2^{p-\ell} \Biggl[\sum_{i_1=1}^\infty {\lambda}_{i_1}^{r q_1} \Biggr] \cdots \Biggl[\sum_{i_\ell=1}^\infty {\lambda}_{i_\ell}^{r q_\ell} \Biggr].
\end{equation}
We would like to determine the number of permutations $\Sigma$ in $\mathrm{Sym}_p$ which consist of $\ell$ cycles  with respective lengths $q_1,\dots, q_\ell$. In order to obtain such a permutation of $\{1,\dots,p\}$, we first build, using the notation of \cite[Section~2.1]{Peccati11}, a partition of $\{1,\dots,p\}$ with $r_i$ blocks of size $i$ for every $i=1,\dots,p$, where $r_1,\dots,r_p \in \mathbb{N}_0$ are such that $1 r_1 + 2 r_2 + \cdots + p r_p = p$. According to  \cite[(2.3.8)]{Peccati11}, the number of such partitions is
\[
\Biggl[\begin{array}{c} p \\ r_1,\dots,r_p \end{array} \Biggr] := \frac{p!}{(1!)^{r_1}\, r_1! \, (2!)^{r_2}\, r_2!\,  \cdots\, (p!)^{r_p}\, r_p!}.
\]
For each of the $r_i$ blocks of size $i$, there are $(i-1)!$ possible cycles formed with the elements of this block. For each choice of these cycles, we obtain a permutation $\Sigma \in \mathrm{Sym}_p$ consisting of $r_1$ cycles of length $1,\dots, r_p$ cycles of length $p$. For this permutation, the term \eqref{3.38} is equal to
\begin{align*}
& 2^{p-(r_1+ \cdots + r_p)} \Biggl[\sum_{i_1 = 1}^\infty {\lambda}_{i_1}^{r \cdot 1} \Biggr]^{r_1} \Biggl[\sum_{i_2 = 1}^\infty {\lambda}_{i_2}^{r \cdot 2} \Biggr]^{r_2} \cdots \Biggl[\sum_{i_p = 1}^\infty {\lambda}_{i_p}^{r \cdot p} \Biggr]^{r_p} \\
& \quad = 2^p \Bigl[\tfrac12 \zeta_D (-r\cdot 1) \Bigr]^{r_1 } \Bigl[\tfrac12 \zeta_D (-r\cdot 2) \Bigr]^{r_2} \cdots \Bigl[\tfrac12 \zeta_D (-r\cdot p) \Bigr]^{r_p} .
\end{align*}
Therefore, with the sums below taken over nonnegative integers $r_1,\dots,r_p$ such that $1 r_1 + 2 r_2 + \cdots + p r_p = p$,
\begin{align*}
\mu_{r,F}({\sigma}^2) &= \sigma^{2p} 2^p \sum_{} \Biggl[\begin{array}{c} p \\ r_1,\dots,r_p \end{array} \Biggr] [(1-1)!]^{r_1} [(2-1)!]^{r_2} \cdots [(p-1)!]^{r_p} \prod_{\ell=1}^p [\tfrac12 \zeta_D(-r\ell)]^{r_\ell} \\
&= \sigma^{2p} 2^p \sum_{} \Biggl[\begin{array}{c} p \\ r_1,\dots,r_p \end{array} \Biggr] \prod_{\ell=1}^p \Biggl[ \frac{(\ell-1)!}{2} \zeta_D (-r\ell) \Biggr]^{r_\ell} \\
&= \sigma^{2p} 2^p B_p(x_1,\dots,x_p),
\end{align*}
where $B_p$ is the complete Bell polynomial in $p$ variables  and $x_\ell$ is defined in \eqref{xl}.
\end{proof}

\subsection{Lemmas used in the proof of Theorem \ref{LLN}, and their proofs}\label{nsec3.2}

   We now prove the lemmas that were quoted in the proof of Theorem \ref{LLN}.

\begin{lemma}\label{remove-drift}
Let $r<-\frac d2$ and the assumptions of Theorem~\ref{LLN} hold. Let $u'(t,\cdot)$ be defined in \eqref{uprime}, where $X_k$ is given by \eqref{Xkt}. Then
\[ V^{n,r}_F(u,t)-V^{n,r}_F(u',t) \stackrel{L^1}{\longrightarrow} 0\qquad\text{as}\qquad n\to\infty. \]
\end{lemma}

For the proof of Lemma \ref{remove-drift}, we need an auxiliary result.
\begin{lemma}\label{tight}
The family
\begin{equation}\label{calu} {\cal U} :=\Biggl\{ \frac{\Delta^n_i u}{\sqrt{\Delta_n}}, \frac{\Delta^n_i u'}{\sqrt{\Delta_n}}: n\in\mathbb{N},~ i=1,\ldots,[T/\Delta_n]\Biggr\}  \end{equation}
of random elements in $H_r$ is tight.
\end{lemma}
\begin{proof} We use the tightness criterion given in \cite[Theorem~1]{Suquet99} and consider the orthonormal basis in $H_r$ given by \eqref{basis_b}.
Then we have for any $1\leq m_1\leq m_2\leq \infty$,
\begin{align*}
\sum_{k=m_1}^{m_2} \langle \Delta^n_i u, b_k\rangle_{H_r}^2 &= \sum_{k=m_1}^{m_2} \Biggl( \sum_{\ell=1}^\infty {\lambda}_\ell^r  a_\ell(\Delta^n_i u)a_\ell(b_k)\Biggr)^2 = \sum_{k=m_1}^{m_2}{\lambda}_k^r \left[a_k(\Delta^n_i u)\right]^2 \\
&= \sum_{k=m_1}^{m_2}{\lambda}_k^r  \Biggl( \iint_0^{i\Delta_n} \Delta^n_i \mathrm{e}_k(s)\phi_k(y){\sigma}(s,y)\,W(\mathrm{d} s,\mathrm{d} y)\Biggr)^2,
\end{align*}
and therefore by \eqref{mom} and calculations similar to those in \eqref{calc},
\begin{equation}\label{series-calc}\begin{split}
\mathbb{E}\Biggl[ \sum_{k=m_1}^{m_2} \langle \Delta^n_i u, b_k\rangle_{H_r}^2\Biggr] &= \sum_{k=m_1}^{m_2}{\lambda}_k^r \iint_0^{i\Delta_n} (\Delta^n_i \mathrm{e}_k(s))^2\phi_k(y)^2\mathbb{E}[{\sigma}^2(s,y)]\,\mathrm{d} s\,\mathrm{d} y\\
&\lesssim \sum_{k=m_1}^{m_2}{\lambda}_k^r \int_0^{i\Delta_n} (\Delta^n_i \mathrm{e}_k(s))^2\,\mathrm{d} s\\
&=\sum_{k=m_1}^{m_2} {\lambda}_k^{r-{\gamma}} (1-\mathrm{e}^{-{\lambda}_k^{\gamma}\Delta_n}) - \frac12\sum_{k=m_1}^{m_2} {\lambda}_k^{r -{\gamma}} \mathrm{e}^{-2{\lambda}_k^{\gamma} i\Delta_n}  (\mathrm{e}^{{\lambda}_k^{\gamma}\Delta_n}-1)^2 \\
&\leq \Delta_n\sum_{k=m_1}^{m_2} {\lambda}_k^r.
\end{split}\end{equation}
Hence, we obtain
\begin{equation}\label{m1m2}
\mathbb{E}\Biggl[ \sum_{k=m_1}^{m_2} \biggl\langle \frac{\Delta^n_i u}{\sqrt{\Delta_n}}, b_k\biggr\rangle_{H_r}^2\Biggr] \lesssim \sum_{k=m_1}^{m_2} {\lambda}_k^r.
\end{equation}
The same estimate holds true when $\Delta^n_i u$ is replaced by $\Delta^n_i u'$ as the reader may quickly verify. We obtain on the one hand, for any $m\geq1$,
\[ \lim_{M\to\infty} \sup_{h \in {\cal U}} \mathbb{P}\Biggl( \sum_{k=1}^m \langle h, b_k\rangle_{H_r}^2 > M \Biggr) \lesssim \lim_{M\to\infty}\frac1M \sum_{k=1}^{m} {\lambda}_k^r = 0, \]
and on the other hand, as a consequence of Lemma~\ref{A1} (i), for any $\delta>0$,
 \[ \lim_{m\to\infty} \sup_{h \in {\cal U}} \mathbb{P}\Biggl( \sum_{k=m}^\infty \langle h, b_k\rangle_{H_r}^2 > \delta \Biggr) \lesssim \lim_{m\to\infty}\frac1\delta \sum_{k=m}^{\infty} {\lambda}_k^r = 0. \]
 Hence, the claim follows from \cite[Theorem~1]{Suquet99}.
\end{proof}

\begin{proof}[Proof of Lemma~\ref{remove-drift}]
By Lemma~\ref{tight}, we can find for every $\theta>0$, a compact subset $K_\theta$ of $H_r$ such that
\begin{equation}\label{theta}
\sup_{n\in\mathbb{N}} \sup_{i=1,\ldots,[t/\Delta_n]} \Biggl\{ \mathbb{P}\biggl( \frac{\Delta^n_i u}{\sqrt{\Delta_n}} \notin K_\theta\biggr) +\mathbb{P}\biggl( \frac{\Delta^n_i u'}{\sqrt{\Delta_n}} \notin K_\theta\biggr) \Biggr\} <\theta.
\end{equation}
Moreover, for all $\theta>0$, since $F$ is uniformly continuous on $K_\theta$ by the Heine--Cantor theorem, there exist for given $\delta>0$, numbers ${\varepsilon}_\theta(\delta)>0$ with the property ${\varepsilon}_\theta(\delta)\to0$ as $\delta\to0$, such that
\begin{equation}\label{unifcont} h,h'\in K_\theta,~\|h-h'\|_{H_r} \leq \delta \implies |F(h)-F(h')|\leq {\varepsilon}_\theta(\delta). \end{equation}
If we further let $M(\theta)=\sup_{h\in K_\theta} |F(h)|$, which is finite for all $\theta>0$, and recall that $|F(h)|\lesssim 1+|h|^p$ for some $p\geq2$, we derive
\begin{equation}\label{Fdiff} |F(h)-F(h')| \lesssim {\varepsilon}_\theta(\delta)+M(\theta)\frac{\|h-h'\|_{H_r}}{\delta} + \mathds 1_{\{h \text{ or } h'\notin K_\theta\}}(1+\|h\|_{H_r}^p + \|h'\|_{H_r}^p). \end{equation}

Suppose now that we can show the following:
\begin{align}
\label{toshow-1} \lim_{n\to\infty} \sup_{i=1,\ldots,[t/\Delta_n]} \frac{\mathbb{E}[\|\Delta^n_i u-\Delta^n_i u'\|_{H_r}]}{\sqrt{\Delta_n}} &= 0,\\
\label{toshow-2} \limsup_{n\to\infty} \sup_{i=1,\ldots,[t/\Delta_n]} \frac{(\mathbb{E}[\|\Delta^n_i u\|_{H_r}^{2p}])^{1/2}+(\mathbb{E}[\|\Delta^n_i u'\|_{H_r}^{2p}])^{1/2}}{\Delta_n^{p/2}}&<\infty.
\end{align}
Then using \eqref{Fdiff} and applying the Cauchy--Schwarz inequality to the third term results in
\begin{align*}
&\mathbb{E}[|V^{n,r}_F(u,t)-V^{n,r}_F(u',t)|]\\
&\quad\leq \Delta_n\, [t/\Delta_n]\sup_{i=1,\ldots,[t/\Delta_n]} \mathbb{E}\Biggl[\biggl\lvert F\biggl(\frac{\Delta^n_i u}{\sqrt{\Delta_n}}\biggr) -F\biggl( \frac{\Delta^n_i u'}{\sqrt{\Delta_n}}\biggr) \biggr\rvert\Biggr]\\
&\quad\lesssim {\varepsilon}_{\theta}(\delta) +M(\theta)\sup_{i=1,\ldots,[t/\Delta_n]}\frac{\mathbb{E}[\|\Delta^n_i u-\Delta^n_i u'\|_{H_r}]}{\delta\sqrt{\Delta_n}}+ \mathbb{P}\biggl(\frac{\Delta^n_i u}{\sqrt{\Delta_n}}\text{ or }\frac{\Delta^n_i u'}{\sqrt{\Delta_n}}\notin K_\theta\biggr)^{1/2}\\
&\qquad\quad\times\sup_{i=1,\ldots,[t/\Delta_n]} \Biggl(1+\frac{(\mathbb{E}[\|\Delta^n_i u\|_{H_r}^{2p}])^{1/2}}{\Delta_n^{p/2}} + \frac{(\mathbb{E}[\|\Delta^n_i u'\|_{H_r}^{2p}])^{1/2}}{\Delta_n^{p/2}}\Biggr).
\end{align*}
Hence, \eqref{theta}, \eqref{toshow-1}, and \eqref{toshow-2} imply
\begin{align*}
\limsup_{n\to\infty} \mathbb{E}[|V^{n,r}_F(u,t)-V^{n,r}_F(u',t)|]\lesssim {\varepsilon}_{\theta}(\delta)+\theta^{\frac12}.
\end{align*}
Letting $\delta\to0$ and then $\theta\to0$ yields the assertion of the lemma.

Therefore, it remains to verify \eqref{toshow-1} and \eqref{toshow-2}. For \eqref{toshow-2},
by the Minkowski and the Burkholder--Davis--Gundy inequalities and  calculations that are similar to \eqref{series-calc}, we have
\begin{align}\nonumber
(\mathbb{E}[\|\Delta^n_i u\|_{H_r}^{2p}])^{\frac1p} &= \Biggl(\mathbb{E}\Biggl[ \Biggl( \sum_{k=1}^\infty {\lambda}_k^r(\Delta^n_i a_k)^2\Biggr)^p\Biggr]\Biggr)^{1/p}\leq \sum_{k=1}^\infty {\lambda}_k^r\Bigl(\mathbb{E}[|\Delta^n_i a_k|^{2p}]\Bigr)^{\frac1p}\\
&\leq \sum_{k=1}^\infty {\lambda}_k^r \iint_0^{i\Delta_n} (\Delta^n_i \mathrm{e}_k(s))^2\phi_k(y)^2(\mathbb{E}[|{\sigma}(s,y)|^{2p}])^{\frac1p}\,\mathrm{d} s\,\mathrm{d} y\lesssim \Delta_n,
\label{rd3.22}
\end{align}
from which the desired bound for $\Delta^n_i u$ follows. The arguments for $\Delta^n_i u'$ are completely analogous. This proves \eqref{toshow-2}.

For \eqref{toshow-1}, we first compute the variance of an increment of the drift process in \eqref{SDE}. By standard stochastic calculus and straightforward manipulations, we have
\begin{align*}
&\mathbb{E}[|\Delta^n_i a_k-\Delta^n_i X_k|^2]\\
&\quad= 2{\lambda}_k^{2{\gamma}} \int_{(i-1)\Delta_n}^{i\Delta_n}\int_{(i-1)\Delta_n}^s\mathbb{E}[a_k(s)a_k(r)]\,\mathrm{d} r\,\mathrm{d} s\\
&\quad=2 {\lambda}_k^{2{\gamma}} \int_{(i-1)\Delta_n}^{i\Delta_n}\int_{(i-1)\Delta_n}^s\iint_0^r \mathrm{e}^{-{\lambda}_k^{\gamma}(s-v)-{\lambda}_k^{\gamma}(r-v)}\phi_k(y)^2\mathbb{E}[{\sigma}^2(v,y)]\,\mathrm{d} v\,\mathrm{d} y\,\mathrm{d} r\,\mathrm{d} s\\
&\quad\lesssim  2{\lambda}_k^{2{\gamma}} \int_{(i-1)\Delta_n}^{i\Delta_n} \mathrm{e}^{-{\lambda}_k^{\gamma} s} \int_{(i-1)\Delta_n}^s\mathrm{e}^{-{\lambda}_k^{\gamma} r}\int_0^r \mathrm{e}^{2{\lambda}_k^{\gamma} v}\,\mathrm{d} v\,\mathrm{d} r\,\mathrm{d} s\\
& \quad={\lambda}_k^{\gamma}  \int_{(i-1)\Delta_n}^{i\Delta_n} \mathrm{e}^{-{\lambda}_k^{\gamma} s} \int_{(i-1)\Delta_n}^s(\mathrm{e}^{{\lambda}_k^{\gamma} r}-\mathrm{e}^{-{\lambda}_k^{\gamma} r})\,\mathrm{d} r\,\mathrm{d} s\\
&\quad ={\lambda}_k^{\gamma}  \int_{(i-1)\Delta_n}^{i\Delta_n} \mathrm{e}^{-{\lambda}_k^{\gamma} (s-(i-1)\Delta_n)} \int_{(i-1)\Delta_n}^s\mathrm{e}^{-{\lambda}_k^{\gamma} (i-1)\Delta_n)}(\mathrm{e}^{{\lambda}_k^{\gamma}  r}-\mathrm{e}^{-{\lambda}_k^{\gamma} r})\,\mathrm{d} r\,\mathrm{d} s\\
&\quad ={\lambda}_k^{\gamma}  \int_{(i-1)\Delta_n}^{i\Delta_n} \mathrm{e}^{-{\lambda}_k^{\gamma} (s-(i-1)\Delta_n)} \int_0^{s-(i-1)\Delta_n} (\mathrm{e}^{{\lambda}_k^{\gamma} r}-\mathrm{e}^{-2{\lambda}_k^{\gamma} (i-1)\Delta_n)}\mathrm{e}^{-{\lambda}_k^{\gamma} r})\,\mathrm{d} r\,\mathrm{d} s\\
&\quad={\lambda}_k^{\gamma}  \int_0^{\Delta_n} \mathrm{e}^{-{\lambda}_k^{\gamma} s} \int_0^{s} (\mathrm{e}^{{\lambda}_k^{\gamma} r}-\mathrm{e}^{-2{\lambda}_k^{\gamma} (i-1)\Delta_n)}\mathrm{e}^{-{\lambda}_k^{\gamma} r})\,\mathrm{d} r\,\mathrm{d} s\\
&\quad=\int_0^{\Delta_n} (1-\mathrm{e}^{-{\lambda}_k^{\gamma} s})\,\mathrm{d} s - \mathrm{e}^{-2{\lambda}_k^{\gamma} (i-1)\Delta_n}\int_0^{\Delta_n} \mathrm{e}^{-{\lambda}_k^{\gamma} s}(1-\mathrm{e}^{-{\lambda}_k^{\gamma} s})\,\mathrm{d} s\\
&\quad\leq \int_0^{\Delta_n}(1-\mathrm{e}^{-{\lambda}_k^{\gamma} s})\,\mathrm{d} s.
\end{align*}
Now let us choose a small $\rho\in(0,1)$ such that $r+{\gamma}\rho<-\frac d2$. Then, with the estimate $1-\mathrm{e}^{-{\lambda}_k^{\gamma}  s}\leq ({\lambda}_k^{\gamma}  s)^{\rho}(1-\mathrm{e}^{-{\lambda}_k^{\gamma}  s})^{1-\rho} \leq ({\lambda}_k^{\gamma}  \Delta_n)^{\rho}$ for all $s\leq\Delta_n$, we obtain $\mathbb{E}[|\Delta^n_i a_k-\Delta^n_i X_k|^2] \lesssim \Delta_n^{1+\rho}{\lambda}_k^{{\gamma} \rho}$, and therefore,
\begin{align*}
\mathbb{E}[\|\Delta^n_i u - \Delta^n_i u'\|^2_{H_r}] &= \sum_{k=1}^\infty {\lambda}_k^r  \mathbb{E}[|\Delta^n_i a_k-\Delta^n_i X_k|^2]\lesssim \Delta_n^{1+\rho}\sum_{k=1}^\infty {\lambda}_k^{r+{\gamma}\rho} .
\end{align*}
By Lemma~\ref{A1} (i), the series converges with our choice of $\rho$, which proves \eqref{toshow-1}.
\end{proof}

Due to the martingale property of $X_k$, the increments of $u'$ are uncorrelated to each other. However, we would like them to be more, namely, conditionally independent given each other. For this purpose, we have to discretize the random field ${\sigma}$.
\begin{lemma}\label{discr-sigma}
Let $h^n_i$ be as defined in \eqref{hni}. As $n\to\infty$,
\[ V^{n,r}_F(u',t) - \Delta_n\sum_{i=1}^{[t/\Delta_n]} F\biggl( \frac{h^n_i}{\sqrt{\Delta_n}}\biggr) \stackrel{L^1}{\longrightarrow} 0. \]
\end{lemma}

\begin{proof}
The proof follows the same scheme as the proof of Lemma~\ref{remove-drift}. The main tool is the estimate \eqref{Fdiff}, this time applied to $\Delta^n_i u'$ and $h^{n}_i$ (note that ${\cal U}$ in \eqref{calu} remains tight if we include also the vectors $h^n_i/\sqrt{\Delta_n}$). Therefore, the lemma is proved once we can show an analogue of \eqref{toshow-1} and \eqref{toshow-2}. We shall only give the details for the modified version of \eqref{toshow-1}.
Notice that, with the notation of \eqref{L2-cont},
\begin{align*}
\mathbb{E}[\|\Delta^n_i u'-h^n_i\|_{H_r}^2]&= \sum_{k=1}^\infty {\lambda}_k^r  \iint_{(i-1)\Delta_n}^{i\Delta_n} \phi_k(y)^2\mathbb{E}[|{\sigma}(s,y)-{\sigma}((i-1)\Delta_n,y)|^2]\,\mathrm{d} s\,\mathrm{d} y\\
&\leq \Delta_n\, w_2(\Delta_n;t)^2 \sum_{k=1}^\infty {\lambda}_k^r.
\end{align*}
Since $w_2(\Delta_n;t)\to0$ as $n\to\infty$, the statement corresponding to \eqref{toshow-1} follows.
\end{proof}

The law of large numbers for $h^n_i$ can now be shown via covariance analysis.
\begin{lemma}\label{LLN-core-1} We have
\begin{equation}\label{LLN-core-eq-1}   \lim_{n\to\infty} \mathbb{E}\Biggl[\Biggl\lvert \Delta_n\sum_{i=1}^{[t/\Delta_n]} \Biggl\{  F\biggl( \frac{h^n_i}{\sqrt{\Delta_n}}\biggr)-\mathbb{E}\biggl[ F\biggl( \frac{h^n_i}{\sqrt{\Delta_n}}\biggr)\mathrel{\bigg|} {\cal F}_{(i-1)\Delta_n}\biggr]\Biggr\}\Biggr\rvert^2\Biggr] = 0. \end{equation}
\end{lemma}
\begin{proof} By construction, if $j>i$, then $h^{n}_i$ and $h^{n}_j$ are conditionally independent given ${\cal F}_{(j-1)\Delta_n}$. Moreover, since $F$ has at most polynomial growth, the variance of $F(h^n_i/\sqrt{\Delta_n})$ can be bounded independently of $n$ and $i$, cf.\ \eqref{toshow-2}. Therefore, the second moment on the left-hand side of \eqref{LLN-core-eq-1} is bounded by a constant times $\Delta_n^2\, [t/\Delta_n]\leq t\Delta_n$, which converges to $0$ as $n\to\infty$.
\end{proof}

The next step is to compute the conditional expectation in \eqref{LLN-core-eq-1}. To this end, we need to recall and prove some results concerning Gaussian measures on Hilbert spaces. By \cite[Theorem~2.1.2]{Prevot07}, the second moment of a mean-zero Gaussian random vector on $H_r$ with covariance operator $Q$  is given by the \emph{trace-class norm} of $Q$:
\begin{equation}\label{trace} \|Q\|_1 := \sum_{k=1}^\infty \langle Qb_k,b_k\rangle_{H_r},  \end{equation}
where $(b_k)_{k\geq1}$ is an arbitrary orthonormal basis of $H_r$. Another norm on ${\cal L}^{+,\mathrm{sym}}_1(H_r)$ is given by the \emph{Hilbert--Schmidt norm}
\begin{equation}\label{HS} \|Q\|_2 := \Biggl(\sum_{k=1}^\infty \|Qb_k\|_{H_r}^2\Biggr)^{1/2}.\end{equation}
By \cite[Remark~B.0.4 and B.0.6]{Prevot07}, neither $\|Q\|_1$ nor $\|Q\|_2$ depends on the choice of $(b_k)_{k\geq1}$. We introduce a third metric on ${\cal L}^{+,\mathrm{sym}}_1(H_r)$ by setting
\begin{equation}\label{dr} d_r(Q,Q') := \|Q-Q'\|_{{\cal L}(H_r)} + | \|Q\|_1 - \|Q'\|_1|,\qquad Q,Q'\in {\cal L}^{+,\mathrm{sym}}_1(H_r),\end{equation}
where $\|\cdot\|_{{\cal L}(H_r)}$ is the usual operator norm for bounded linear operators $H_r\to H_r$.

\begin{lemma}\label{comp-condexp}
We have
\[ \mathbb{E}\Biggl[ F\biggl( \frac{h^n_i}{\sqrt{\Delta_n}}\biggr)\mathrel{\bigg|} {\cal F}_{(i-1)\Delta_n}\Biggr] = \mu_F(Q^{n,i}_r), \]
where $Q^{n,i}_r\colon H_r \to H_r$ is given by
\begin{equation}\label{Qni
(Q^{n,i}_r h)(h'):= \sum_{k,\ell=1}^\infty  {\lambda}_\ell^r  a_k(h')a_\ell(h)\int_D \phi_k(y)\phi_\ell(y) {\sigma}^2((i-1)\Delta_n,y)\,\mathrm{d} y\end{equation}
for $h\in H_r$ and $h'\in H_{-r}$.
\end{lemma}
\begin{proof} For any $h\in H_r$, we have
\[
\langle h^n_i, h\rangle_{H_r} = \sum_{k=1}^\infty {\lambda}_k^r  a_k(h^n_i)a_k(h) = \sum_{k=1}^\infty {\lambda}_k^r a_k(h) \iint_{(i-1)\Delta_n}^{i\Delta_n}\phi_k(y){\sigma}((i-1)\Delta_n,y)\,W(\mathrm{d} s,\mathrm{d} y),
\]
which is, conditionally on ${\cal F}_{(i-1)\Delta_n}$, normally distributed with mean $0$ and variance
\begin{equation}\label{vni} v^n_i(h):= \Delta_n\sum_{k,\ell=1}^\infty {\lambda}_k^r{\lambda}_\ell^r a_k(h)a_\ell(h)c^{n,i}_{k,\ell}, \end{equation}
where $c^{n,i}_{k,\ell}$ is the integral in \eqref{Qni}. Hence, by \cite[Definition~2.1.1]{Prevot07}, $h^n_i/\sqrt{\Delta_n}$ is conditionally Gaussian given ${\cal F}_{(i-1)\Delta_n}$ with zero mean. It remains to show that its conditional covariance operator is given by $Q^{n,i}_r$. But this follows
immediately from \eqref{Qni} because
\begin{equation}\label{ak}\begin{aligned}
	a_k(Q^{n,i}_r h) = (Q^{n,i}_r h)(\phi_k)= \sum_{\ell=1}^\infty{\lambda}_\ell^r a_\ell(h) c^{n,i}_{k,\ell}
\end{aligned}\end{equation}
and consequently,
\begin{equation}\label{QHH} \langle Q^{n,i}_r h, h\rangle_{H_r} = \sum_{k=1}^\infty {\lambda}_k^r a_k(Q^{n,i}_r h)a_k(h)=\frac{1}{\Delta_n} v^n_i(h). \end{equation}
\end{proof}

\begin{lemma}\label{muf-prop}
If $F\colon H_r\to\mathbb{R}$ is continuous with $|F(h)|\lesssim 1+\|h\|_{H_r}^p$ for some $p>0$, then $\mu_F$ is continuous on ${\cal L}^{+,\mathrm{sym}}_1(H_r)$ with respect to $d_r$, and $|\mu_F(Q)| \lesssim 1+\|Q\|_{1}^{p/2}$.
\end{lemma}
\begin{proof} Since $\|Q\|_1$ is precisely the second moment of $N_r(0,Q)$, the bound on $\mu_F(Q)$ follows immediately from the growth assumption on $F$ and \cite[Corollary~3.2]{Ledoux91}. For the continuity of $\mu_F$, let
$Q_n$ and $Q$ be elements of ${\cal L}^{+,\mathrm{sym}}_1(H_r)$ such that $d_r(Q_n,Q)\to0$.
Then, the assumptions of \cite[Example~3.8.15]{Bogachev98} are satisfied, and we have $N_r(0,Q_n)\stackrel{\rm w}{\longrightarrow} N_r(0,Q)$ as $n\to\infty$, where $\stackrel{\rm w}{\longrightarrow}$ denotes weak convergence of measures.
Since $F$ is of polynomial growth and the moments of any order of $N_r(0,Q_n)$ are uniformly bounded in $n$ (as the first part of the proof has shown), we obtain $\mu_F(Q_n)\to \mu_F(Q)$ from \cite[Lemma 3.8.7]{Bogachev98}.
\end{proof}

\begin{comment}
\begin{lemma}\label{tight-Q} Let $Q_r(u)$ be defined in the same way as $Q^{n,i}_r$ in \eqref{Qni}, but with $(i-1)\Delta_n$ replaced by $u$. For every $r<-\frac12$ and $t>0$, the collection of random operators given by
\begin{equation}\label{calq} {\cal Q}_r =\{Q^{n}_r(u)\colon n\in\mathbb{N},  \, u\in[0,t]\} \end{equation}
is tight in ${\cal L}^{+,\rm{sym}}_1(H_r)$ with respect to the metric $d_r$.
\end{lemma}
\begin{proof}
By the definition of $d_r$, the statement is equivalent to the tightness of ${\cal Q}_r^{1/2} = \{Q^n_r(u)^{1/2}\colon n\in\mathbb{N},\, u\in[0,t]\}$ in the Hilbert--Schmidt norm. By \cite[Proposition~B.0.7]{Prevot07}, the space $L_2(H_r)$ of all bounded linear operators $H_r\to H_r$ for which \eqref{HS} is finite becomes a separable Hilbert space with the inner product
\[ \langle Q_1,Q_2\rangle_2 = \sum_{k=1}^\infty \langle Q_1b_k,Q_2b_k\rangle_{H_r},\quad Q_1,Q_2\in L_2(H_r). \]
In the following, we take $(b_k=\sqrt{\frac2\pi}(1+k^2)^{-r/2}\sin(k\cdot)\colon k\geq1)$ as orthonormal basis for $H_r)$. Then it is not difficult to see that the operators $B_{jk}\colon H_r\to H_r$ given by $B_{jk}H=\langle H,b_j\rangle_{H_r}b_k$, where $j,k\geq1$, form an orthonormal basis of $L_2(H_r)$.

So for the proof of tightness, we can again apply \cite[Theorem~1]{Suquet99}. For  $1\leq j_1\leq j_2\leq \infty$ and $1\leq k_1\leq k_2\leq \infty$, we have by Parseval's identity,
\begin{equation}\label{exp}\begin{aligned}
 \sum_{j=j_1}^{j_2}\sum_{k=k_1}^{k_2}  \langle Q^{n}_r(u)^{1/2},B_{jk}\rangle_{2}^2   &= \sum_{j=j_1}^{j_2}\sum_{k=k_1}^{k_2} \left( \sum_{m=1}^\infty \langle Q^n_r(u)^{1/2}b_m, B_{jk}b_m\rangle_{H_r} \right)^2 \\
	&=\sum_{j=j_1}^{j_2}\sum_{k=k_1}^{k_2} \langle Q^n_r(u)^{1/2}b_j,b_k\rangle_{H_r}^2\\
	&\leq \sum_{j=j_1}^{j_2} \| Q^n_r(u)^{1/2}b_j\|^2_{H_r} = \sum_{j=j_1}^{j_2} \langle Q^n_r(u)b_j,b_j\rangle_{H_r},
\end{aligned}\end{equation}
where the last step is valid because $Q^n_r(u)^{1/2}$ is symmetric. If we let $c^{n}_{k,l}(u)$ be the double integral in \eqref{Qni} with $i\Delta_n$ (resp., $(i-1)\Delta_n$) replaced by $u$ (resp., $u-\Delta_n$), then by \eqref{vni} and \eqref{QHH},
\begin{align*}
\mathbb{E}\left[ \sum_{j=j_1}^{j_2}\sum_{k=k_1}^{k_2}  \langle Q^{n}_r(u)^{1/2},B_{jk}\rangle_{2}^2\right] &\leq \frac2{\pi\Delta_n} \sum_{j=j_1}^{j_2} \sum_{l,m=1}^\infty (1+l^2)^r(1+m^2)^r a_l(b_j)a_m(b_j)\mathbb{E}[c^{n,i}_{l,m}]\\
&\lesssim  \sum_{j=j_1}^{j_2} (1+j^2)^r.
\end{align*}
Since $\langle Q^n_r(u)^{1/2}b_j,b_k\rangle_{H_r}^2\rangle_{H_r} = \langle Q^n_r(u)^{1/2}b_k,b_j\rangle_{H_r}^2\rangle_{H_r}$ in the second line of \eqref{exp} by the symmetry of $Q^n_r(u)^{1/2}$, the last calculation actually implies
\begin{equation}\label{min} \mathbb{E}\left[ \sum_{j=j_1}^{j_2}\sum_{k=k_1}^{k_2}  \langle Q^{n}_r(u)^{1/2},B_{jk}\rangle_{2}^2\right] \lesssim \min\left( \sum_{j=j_1}^{j_2} (1+j^2)^r, \sum_{k=k_1}^{k_2} (1+k^2)^r\right). \end{equation}
From this, by similar arguments as in the last part of the proof of Lemma~\ref{tight}, the two conditions in \cite[Theorem~1]{Suquet99} hold, and ${\cal Q}_r$ is tight as claimed.
\end{proof}
\end{comment}

The next lemma is the final ingredient in the proof of Theorem~\ref{LLN} for $r<-\frac d2$.
\begin{lemma}\label{cex-conv} We have
\begin{equation}\label{eq1} \Delta_n\sum_{i=1}^{[t/\Delta_n]} \mu_F(Q^{n,i}_r) - \int_0^t \mu_F({\sigma}^2(s,\cdot))\,\mathrm{d} s \stackrel{L^1}{\longrightarrow}0\qquad\text{as} \qquad n\to\infty. \end{equation}
\end{lemma}
\begin{proof} Let $Q_r(v)$ (resp., $c_{k,\ell}(v)$) be defined in the same way as $Q^{n,i}_r$ in \eqref{Qni} (resp., the integral in \eqref{Qni}), but with $(i-1)\Delta_n$ replaced by $v$. By definition (see \eqref{Qphi1}), we have $\mu_F(Q_r(v))=\mu_{r,F}({\sigma}^2(v,\cdot))$ and $Q^{n,i}_r=Q_r((i-1)\Delta_n)$. Therefore, we can write the difference in \eqref{eq1} as $C^n_1+C^n_2$ where
\begin{align*}
C^n_1 &:=\int_0^t \Biggl\{\sum_{i=1}^{[t/\Delta_n]} \Bigl(\mu_F(Q_r((i-1)\Delta_n)) - \mu_F(Q_r(v))\Bigr)\mathds 1_{((i-1)\Delta_n,i\Delta_n]}(v)\Biggr\}\, \mathrm{d} v,\\
C^n_2&:=- \int_{[t/\Delta_n]\Delta_n}^{t} \mu_F(Q_r(v))\,\mathrm{d} v.
\end{align*}
Suppose we show that
\begin{equation}\label{Q-conv} \sup_{v,w\in[0,t]:\,|v-w|\leq\Delta_n}\mathbb{E}[d_r(Q_r(v),Q_r(w))^2] \lesssim w_2(\Delta_n;t)^2. \end{equation}
Then, together with the $L^2({\Omega})$-continuity of ${\sigma}$ and the continuity of $Q\mapsto \mu_F(Q)$ (see Lemma~\ref{muf-prop}), this will imply that for every $v\in[0,t]$, the integrand in the definition of $C^n_1$ converges to $0$ in probability as $n\to\infty$. In fact, this convergence also takes place in $L^1({\Omega})$ by uniform integrability: indeed, by Lemma~\ref{muf-prop}, \eqref{mom} and equations \eqref{trace}, \eqref{vni}, and \eqref{QHH},
\begin{equation}\label{muf-2}\begin{aligned}
\mathbb{E}[\mu_F(Q_r(v))^2] &\lesssim 1+ \mathbb{E}[\|Q_r(v)\|_1^p] \leq 1+\mathbb{E}\Biggl[ \Biggl(  \sum_{k=1}^\infty {\lambda}_k^rc_{k,k}(v)\Biggr)^p\Biggr]\\
&\leq 1+\Biggl(  \sum_{k=1}^\infty {\lambda}_k^r \mathbb{E}[|c_{k,k}(v)|^p]^{\frac1p}\Biggr)^p \lesssim 1+\Biggl(\sum_{k=1}^\infty {\lambda}_k^r\Biggr)^p,
\end{aligned}\end{equation}
which is finite and independent of $v$. Therefore, for $v \in [0,t]$,
\begin{align*}
    &  \mathbb{E}\Biggl[ \Biggl(\sum_{i=1}^{[t/\Delta_n]} \Bigl(\mu_F(Q_r((i-1)\Delta_n) - \mu_F(Q_r(v))\Bigr)\mathds 1_{((i-1)\Delta_n,i\Delta_n]}(v)\Biggr)^2 \Biggr] \\
   &\quad =  \mathbb{E} \Biggl[\sum_{i=1}^{[t/\Delta_n]} \Bigl(\mu_F(Q_r((i-1)\Delta_n) - \mu_F(Q_r(v))\Bigr)^2\mathds 1_{((i-1)\Delta_n,i\Delta_n]}(v)\Biggr]  \\
   &\quad\lesssim 2 \sum_{i=1}^{[t/\Delta_n]} \Biggl( 1+\Biggl(\sum_{k=1}^\infty {\lambda}_k^r\Biggr)^p \Biggr) \mathds 1_{((i-1)\Delta_n,i\Delta_n]}(v)   \leq 2 \Biggl( 1+\Biggl(\sum_{k=1}^\infty {\lambda}_k^r\Biggr)^p\Biggr)  < \infty.
\end{align*}
This establishes the uniform integrability of the integrand in the definition of $C^n_1$. Applying dominated convergence to the $\mathrm{d} v$-integral, we obtain $\mathbb{E}[|C^n_1|]\to0$. The same estimate \eqref{muf-2} also shows $\mathbb{E}[|C^n_2|]\to0$.

Hence, it remains to prove \eqref{Q-conv}, and we shall consider the two parts defining the metric $d_r$ in \eqref{dr} separately. Using \cite[Remark~B.0.6(ii)]{Prevot07}, we obtain
\begin{equation}\label{eq2}
\|Q_r(v)-Q_r(w)\|_{{\cal L}(H_r)}^2 \leq \|Q_r(v)-Q_r(w)\|_2^2 = \sum_{k=1}^\infty \| Q_r(v)b_k-Q_r(w)b_k\|_{H_r}^2.
\end{equation}
Recalling the formula for $b_k$ from \eqref{basis_b},
\[
\langle Q_r(v)b_k,Q_r(w)b_k\rangle_{H_r} = \sum_{j=1}^\infty {\lambda}_j^r a_j(Q_r(v)b_k)a_j(Q_r(w)b_k)= {\lambda}_k^r \sum_{j=1}^\infty {\lambda}_j^r c_{j,k}(v)c_{j,k}(w)
\]
by \eqref{ak}, so expanding the square in the last term in \eqref{eq2} yields for $v,w \in [0,1]$ with $\lvert v - w \rvert \leq \Delta_n$,
\begin{align*}
\mathbb{E}[\|Q_r(v)-Q_r(w)\|_{{\cal L}(H_r)}^2]\leq   \sum_{k=1}^\infty {\lambda}_k^r \sum_{j=1}^\infty {\lambda}_j^r \mathbb{E}[ (c_{j,k}(v)-c_{j,k}(w))^2] \lesssim w_2(\Delta_n;t)^2
\end{align*}
where the last inequality uses the Cauchy-Schwartz inequality and the fact that $\sigma$ has bounded moments by \eqref{mom}, as well as Lemma~\ref{A1} (i)  and the fact that $r < -\tfrac d2$. For the second part of the metric $d_r$, observe from \eqref{vni} and \eqref{QHH} that
\begin{align*}
\|Q_r(v)\|_1-\|Q_r(w)\|_1 &= \sum_{k=1}^\infty \Big( \langle Q_r(v)b_k,b_k\rangle_{H_r} -  \langle Q_r(w)b_k,b_k\rangle_{H_r}\Big)\\
&= \sum_{k=1}^\infty {\lambda}_k^r \Big( c_{k,k}(v) - c_{k,k}(w)\Big).
\end{align*}
This in turn gives
\begin{align*}
\mathbb{E}\biggl[\Big|\|Q_r(v)\|_1-\|Q_r(w)\|_1\Big|^2\biggr] \leq \Biggl(  \sum_{k=1}^\infty {\lambda}_k^r \mathbb{E}[( c_{k,k}(v) - c_{k,k}(w))^2]^{\frac12} \Biggr)^2 \lesssim w_2(\Delta_n;t)^2,
\end{align*}
which completes the proof of \eqref{Q-conv}.
\end{proof}

\subsection{Proof of Theorem~\ref{LLN2} ($-\tfrac{d}{2} \leq r <{\gamma}- \tfrac d2$)}\label{nsec3.3}
\begin{proof}[Proof of Theorem~\ref{LLN2}]
As in the proof of Theorem~\ref{LLN}, we will approximate $\|\Delta^n_i u\|_{H_r}$ by simpler expressions: For ${\varepsilon}>0$, $i= 1,\dots, [t/\Delta_n]$, $k \geq 1$, and $-\frac d2\leq r<{\gamma}-\frac d2$, we consider, in a first truncation step,
\begin{equation}\label{hnieps} a^{n,i,{\varepsilon}}_k:=\iint_{(i\Delta_n-{\varepsilon})\vee 0}^{i\Delta_n} \Delta^n_i \mathrm{e}_k(s)\phi_k(y){\sigma}(s,y
\,W(\mathrm{d} s,\mathrm{d} y),\qquad
h^{n,i,{\varepsilon}}_r :=\Biggl(\sum_{k=1}^\infty {\lambda}_k^r(a^{n,i,{\varepsilon}}_k)^2\Biggr)^{1/2}\end{equation}
as approximations for $\Delta^n_i a_k$ (as defined in \eqref{abbr}) and $\|\Delta^n_i u\|_{H_r}$, respectively. In a second step, we discretize the random field ${\sigma}$ by introducing the variables
\begin{equation}\label{wthnieps}\begin{split} \widetilde a^{n,i,{\varepsilon}}_k&:=\iint_{(i\Delta_n-{\varepsilon})\vee 0}^{i\Delta_n} \Delta^n_i \mathrm{e}_k(s)\phi_k(y){\sigma}(i\Delta_n-{\varepsilon},y
\,W(\mathrm{d} s,\mathrm{d} y),\\
\widetilde h^{n,i,{\varepsilon}}_r  &:=\Biggl(\sum_{k=1}^\infty {\lambda}_k^r(\widetilde a^{n,i,{\varepsilon}}_k)^2\Biggr)^{1/2},  \end{split}\end{equation}
where  ${\sigma}$ is extended to $\mathbb{R}\times D$ by setting ${\sigma}(s,y):={\sigma}(0,y)$ for all $s<0$ and $y\in D$.

In this set-up, we consider the decomposition
\[ V^{n,r}_f(u,t)-V^r_f(u,t)=J^{\varepsilon}_{n,1}+J^{\varepsilon}_{n,2}+J^{\varepsilon}_{n,3},  \]
where
\begin{align*}
J^{\varepsilon}_{n,1}&:=V^{n,r}_f(u,t)-\Delta_n\sum_{i=1}^{[t/\Delta_n]} f\biggl( \frac{h^{n,i,{\varepsilon}}_r}{\tau_n(r)}\biggr), \\
J^{\varepsilon}_{n,2}&:=\Delta_n\sum_{i=1}^{[t/\Delta_n]}   f\biggl( \frac{h^{n,i,{\varepsilon}}_r}{\tau_n(r)}\biggr)-\Delta_n\sum_{i=1}^{[t/\Delta_n]} f\biggl( \frac{\widetilde h^{n,i,{\varepsilon}}_r}{\tau_n(r)}\biggr),\\
J^{\varepsilon}_{n,3}&:=\Delta_n\sum_{i=1}^{[t/\Delta_n]}  f\biggl( \frac{\widetilde h^{n,i,{\varepsilon}}_r}{\tau_n(r)}\biggr) -V^r_f(u,t)
\end{align*}
and prove \eqref{eq:LLN-2} by showing
\[ \lim_{{\varepsilon}\to0}\limsup_{n\to\infty} \mathbb{E}[|J^{\varepsilon}_{n,1}|+|J^{\varepsilon}_{n,2}|+|J^{\varepsilon}_{n,3}|] =0 \]
in Lemmas~\ref{lemma-eps}, \ref{LLN-discr}, and \ref{LLN-condexp} below. Indeed, as mentioned at the beginning of Section \ref{sec3.1}, having $\mathbb{E}[|V^{n,r}_f(u,t)-V^r_f(u,t)|]\to0$ for all $t\in[0,T]$  is sufficient to obtain the uniform convergence to $0$ on $[0,T]$ in $L^q({\Omega})$ that is claimed in \eqref{eq:LLN-2}. This establishes Theorem \ref{LLN2}.
\end{proof}

While Lemmas~\ref{lemma-eps} and \ref{LLN-discr} (needed for the proof of Theorem~\ref{LLN}) have similar proofs to Lemmas~\ref{remove-drift} and \ref{discr-sigma} (needed for the proof of Theorem~\ref{LLN2}), both the assertion and the proof of Lemma~\ref{LLN-condexp} are significantly different from what we saw in Lemma~\ref{LLN-core-1}. In fact, the main step in proving Lemma~\ref{LLN-condexp} is Lemma~\ref{convP}, where we show that, in the case $-\frac d2\leq r<{\gamma}-\frac d2$, not only the  \emph{average} but actually \emph{each} of the variables $\widetilde h^{n,i,{\varepsilon}}_r/\tau_n(r)$ converges in $L^2({\Omega})$; cf.\ Remark~\ref{difference} (ii).

\subsection{Lemmas used in the proof of Theorem \ref{LLN2}, and their proofs}\label{nsec3.4}
\begin{lemma}\label{lemma-eps} Let $-\frac d2\leq r<{\gamma}-\frac d2$ and consider the variables $h^{n,i,{\varepsilon}}_r$ from \eqref{hnieps}.
Then
\begin{equation}\label{diff} \lim_{{\varepsilon}\to0} \limsup_{n\to\infty} \mathbb{E}\Biggl[\biggl\lvert V^{n,r}_f(u,t)-\Delta_n\sum_{i=1}^{[t/\Delta_n]} f\biggl( \frac{h^{n,i,{\varepsilon}}_r}{\tau_n(r)}\biggr)\biggr\rvert\Biggr] = 0. \end{equation}
\end{lemma}
\begin{proof} The difference in \eqref{diff} is
\begin{align*}
D^{n,{\varepsilon}}_r &:=\Delta_n\sum_{i=1}^{[t/\Delta_n]} \Biggl\{f\biggl(\frac{\|\Delta^n_i u\|_{H_r}}{\tau_n(r)}\biggr) -f\biggl( \frac{h^{n,i,{\varepsilon}}_r}{\tau_n(r)}\biggr)\Biggr\}
\end{align*}
(note that for $i=1,\dots, [{\varepsilon}/\Delta_n]$, the terms in this sum vanish).
Since $f$ is continuous with polynomial growth, we can find $p\geq2$ such that $|f(x)|/x^p\to 0$ as $x\to \infty$. Hence, similarly to \cite[Equation~(3.4.16)]{Jacod12}, we have
\begin{equation}\label{fdiff} |f(x+y)-f(x)| \leq C\biggl( \Phi'_A(\delta) + \Phi(A) \frac{|y|}{\delta} + \Phi''(A)(x^p+y^p)\biggr),\qquad x,y\geq0,\end{equation}
where $C>0$ only depends on $p$, and $\Phi(A)$, $\Phi'_A(\delta)$, and $\Phi''(A)$ are finite numbers satisfying $\Phi'_A(\delta)\to0$ as $\delta\to0$ for given $A>0$ and $\Phi''(A)\to0$ as $A\to\infty$. If we can show that
\begin{align}
\label{toshow1} &\lim_{n\to\infty} \sup_{i=1,\ldots,[t/\Delta_n]} \frac{\mathbb{E}[|\|\Delta^n_i u\|_{H_r}-h^{n,i,{\varepsilon}}_r|]}{\tau_n(r)} = 0,\\
\label{toshow2} &\limsup_{{\varepsilon}\to0} \limsup_{n\to\infty} \sup_{i=1,\ldots,[t/\Delta_n]} \frac{\mathbb{E}[\|\Delta^n_i u\|_{H_r}^p]+\mathbb{E}[(h^{n,i,{\varepsilon}}_r)^p]}{\tau_n(r)^p}<\infty,
\end{align}
then \eqref{fdiff} implies
\begin{align*}
\limsup_{n\to\infty} \mathbb{E}[|D^{n,{\varepsilon}}_r|] &\leq \limsup_{n\to\infty} \Delta_n\, [t/\Delta_n]\sup_{i=1,\ldots,[t/\Delta_n]} \mathbb{E}\Biggl[\biggl\lvert f\biggl(\frac{\|\Delta^n_i u\|_{H_r}}{\tau_n(r)}\biggr) -f\biggl( \frac{h^{n,i,{\varepsilon}}_r}{\tau_n(r)}\biggr) \biggr\rvert\Biggr]\\
&\lesssim \Phi'_A(\delta)+\Phi''(A),
\end{align*}
which goes to $0$ if we first let $\delta\to0$ and then $A\to\infty$.
This establishes \eqref{diff}.

So it remains to show \eqref{toshow1} and \eqref{toshow2}.
For \eqref{toshow2}, we note that
\begin{align*}
   \mathbb{E}[\|\Delta^n_i u\|_{H_r}^p] &= \mathbb{E}\Biggl[\Biggl(\sum_{k=1}^\infty {\lambda}_k^r (\Delta^n_i a_k)^2 \Biggr)^{p/2} \Biggr] = \Biggl\lVert \sum_{k=1}^\infty {\lambda}_k^r(\Delta^n_i a_k)^2 \Biggr\rVert_{L^{p/2}(\Omega)}^{p/2}  \\
   &\leq \Biggl(\sum_{k=1}^\infty {\lambda}_k^r \lVert \Delta^n_i a_k \rVert_{L^p(\Omega)}^2 \Biggr)^{p/2}
\end{align*}
by Minkowski's inequality. Now we apply the Burkholder--Davis--Gundy inequality, and again the Minkowski's inequality and hypothesis \eqref{mom} to see, as in the \eqref{calc}, that
\begin{align*}
   \lVert \Delta^n_i a_k \rVert_{L^p(\Omega)}^2 &\leq \Biggl\lVert \iint_0^{i\Delta_n} (\Delta^n_i \mathrm{e}_k(s))^2\phi_k(y)^2{\sigma}^{2}(s,y)\,\mathrm{d} s\,\mathrm{d} y \Biggr\rVert_{L^{p/2}(\Omega)} \\
    &\leq \iint_0^{i\Delta_n} (\Delta^n_i \mathrm{e}_k(s))^2\phi_k(y)^2(\mathbb{E}[|{\sigma}(s,y)|^p])^{\frac 2p}\,\mathrm{d} s\,\mathrm{d} y \lesssim \frac{1-\mathrm{e}^{-{\lambda}_k^{\gamma}\Delta_n}}{{\lambda}_k^{\gamma}},
\end{align*}
so
$
   \mathbb{E}[\|\Delta^n_i u\|_{H_r}^p] \lesssim  (A^{n,1}_r)^{p/2},
$
where $A^{n,1}_r/\tau_n(r)^2 \to K_r$ as seen in the proof of Lemma~\ref{size-incr}.
This show \eqref{toshow2} for the first term.  The proof for the second term is similar and is left to the reader.

We now prove \eqref{toshow1}. To this end, we use the reverse triangle inequality to deduce for $i\Delta_n\geq{\varepsilon}$,
\begin{equation}\label{diff-mom}\begin{aligned}
	\mathbb{E}[(\|\Delta^n_i u\|_{H_r}-h^{n,i,{\varepsilon}}_r)^2]&\leq \sum_{k=1}^\infty {\lambda}_k^r \mathbb{E}[(\Delta^n_i a_k- a^{n,i,{\varepsilon}}_k)^2]\\
	&=\sum_{k=1}^\infty {\lambda}_k^r \iint_0^{i\Delta_n-{\varepsilon}} (\Delta^n_i \mathrm{e}_k(s))^2\phi_k(y)^2\mathbb{E}[{\sigma}^2(s,y)]\,\mathrm{d} s\,\mathrm{d} y\\
	&\lesssim \sum_{k=1}^\infty {\lambda}_k^r \int_0^{i\Delta_n-{\varepsilon}} (\mathrm{e}^{-{\lambda}_k^{\gamma}(i\Delta_n-s)}-\mathrm{e}^{-{\lambda}_k^{\gamma}((i-1)\Delta_n-s)})^2\,\mathrm{d} s\\
	&=\frac12\sum_{k=1}^\infty {\lambda}_k^{r-{\gamma}}(1-\mathrm{e}^{{\lambda}_k^{\gamma}\Delta_n})^2(\mathrm{e}^{-2{\lambda}_k^{\gamma}{\varepsilon}}-\mathrm{e}^{-2{\lambda}_k^{\gamma} i\Delta_n})\\
	&\leq\frac12 \sum_{k=1}^\infty {\lambda}_k^{r-{\gamma}}(1-\mathrm{e}^{{\lambda}_k^{\gamma}\Delta_n})^2\mathrm{e}^{-2{\lambda}_k^{\gamma}{\varepsilon}}.
\end{aligned}\end{equation}
Due to the presence of the exponential term $\mathrm{e}^{-2{\lambda}_k^{\gamma}{\varepsilon}}$,
the right-hand side converges to a finite limit after dividing by $\Delta_n^2$, which shows \eqref{toshow1} by Jensen's inequality because $\Delta_n/\tau_n(r)\to0$.
\end{proof}

\begin{lemma}\label{LLN-discr}
With the variables $h^{n,i,{\varepsilon}}_r$ and $\widetilde h^{n,i,{\varepsilon}}_r$ from \eqref{hnieps} and \eqref{wthnieps}, we have,  as $n\to\infty$, that
\begin{equation}\label{diff2} \lim_{{\varepsilon}\to0} \limsup_{n\to\infty} \Delta_n\sum_{i=1}^{[t/\Delta_n]} \mathbb{E}\Biggl[\biggl\lvert  f\biggl( \frac{h^{n,i,{\varepsilon}}_r}{\tau_n(r)}\biggr)-f\biggl( \frac{\widetilde h^{n,i,{\varepsilon}}_r}{\tau_n(r)}\biggr)\biggr\rvert\Biggr] = 0. \end{equation}
\end{lemma}
\begin{proof}
As in the proof of Lemma~\ref{lemma-eps}, it suffices to show \eqref{toshow1} and \eqref{toshow2} with $\|\Delta^n_i u\|_{H_r}$ replaced by $\widetilde h^{n,i,{\varepsilon}}_r$. We shall only show the modified version of \eqref{toshow1}. With the same calculations as in \eqref{diff-mom}, we obtain for all $i\geq 1$,
\begin{align*}
&\mathbb{E}[(h^{n,i,{\varepsilon}}_r-\widetilde h^{n,i,{\varepsilon}}_r)^2] \\
&\quad\leq \sum_{k=1}^\infty  {\lambda}_k^r \iint_{(i\Delta_n-{\varepsilon})\vee 0}^{i\Delta_n} (\Delta^n_i \mathrm{e}_k(s))^2\phi_k(y)^2\mathbb{E}[({\sigma}(s,y)-{\sigma}(i\Delta_n-{\varepsilon},y))^2]\,\mathrm{d} s\,\mathrm{d} y\\
&\quad\lesssim w_2({\varepsilon};t)^2 \sum_{k=1}^\infty  {\lambda}_k^r \Biggl( \int_0^{i\Delta_n} \mathrm{e}^{-2 {\lambda}_k^{\gamma} s} \,\mathrm{d} s+ \int_0^{(i-1)\Delta_n}\mathrm{e}^{-2 {\lambda}_k^{\gamma} s} \,\mathrm{d} s\Biggr) \lesssim w_2({\varepsilon};t)^2 \sum_{k=1}^\infty {\lambda}_k^{r-{\gamma}},
\end{align*}
so the $L^2({\Omega})$-continuity of ${\sigma}$ completes the proof since $r -{\gamma}< -\frac d2$.
\end{proof}

\begin{lemma}\label{Pinr} For $n\in\mathbb{N}$ and $-\frac d2\leq r<{\gamma}-\frac d2$, consider the measure on $[0,\infty)\times D$ defined by
\begin{equation}\label{Pin} \Pi^{n}_{r,{\gamma}}(\mathrm{d} s,\mathrm{d} y) := \frac1{\tau_n(r)^2} \sum_{k=1}^\infty {\lambda}_k^r (\mathrm{e}^{-{\lambda}_k^{\gamma} s}-\mathrm{e}^{-{\lambda}_k^{\gamma}(s-\Delta_n)}\mathds 1_{s\geq\Delta_n})^2\phi_k(y)^2\,\mathrm{d} s\,\mathrm{d} y. \end{equation}
Then, as $n\to\infty$,
\begin{equation}\label{weak} \Pi^{n}_{r,{\gamma}} \stackrel{\rm w}{\longrightarrow} \frac{K_r}{|D|}\,\delta_0(\mathrm{d} s)\, \mathds 1_D(y)\,\mathrm{d} y, \end{equation}
where $K_r$ is the number from \eqref{Kr}.
\end{lemma}
\begin{proof} Clearly, \eqref{weak} is equivalent to showing
\begin{equation}\label{weak-2} \begin{split} &\lim_{n\to\infty} \Pi^n_{r,{\gamma}}([0,t]\times D_z) \\
	&\quad= \lim_{n\to\infty} \int_0^t \int_{D_z} \frac1{\tau_n(r)^2} \sum_{k=1}^\infty {\lambda}_k^r (\mathrm{e}^{-{\lambda}_k^{\gamma} s}-\mathrm{e}^{-{\lambda}_k^{\gamma} (s-\Delta_n)}\mathds 1_{s\geq\Delta_n})^2\phi_k(y)^2\,\mathrm{d} s\,\mathrm{d} y \\
	&\quad= K_r\frac{|D_z|}{|D|} = \frac{K_r}{|D|}\iint \mathds 1_{[0,t]\times D_z}(s,y)\,\delta_0(\mathrm{d} s)\,\mathrm{d} y\end{split}\end{equation}
for all $t> 0$ and $z\in \mathbb{R}^d$, where $D_z =\{y\in D: y_i\leq z_i \text{ for all } i=1,\dots,d\}$.

Consider the double integral in the first line of \eqref{weak-2}. With similar calculations as in \eqref{calc}, if we first integrate with respect to $s$, it is equal to
\begin{equation}\label{weak-3}  \frac1{\tau_n(r)^2}  \int_{D_z} \sum_{k=1}^\infty {\lambda}_k^{r-{\gamma}}\Biggl( (1-\mathrm{e}^{-{\lambda}_k^{\gamma}\Delta_n})-\frac{\mathrm{e}^{-2{\lambda}_k^{\gamma} t}}{2}(\mathrm{e}^{{\lambda}_k^{\gamma} \Delta_n}-1)^2\Biggr)\phi_k(y)^2\,\mathrm{d} y.  \end{equation}
Thus, \eqref{weak-2} is equivalent to the convergence of the previous line to $\frac{K_r}{|D|}\int_{D_z} 1\, \mathrm{d} y$ for all $t>0$ and $z\in\mathbb{R}^d$. If $z$ is such that $D_z=D$, then \eqref{weak-3} is exactly the last line in \eqref{calc} (with $i\Delta_n$ replaced by $t$) divided by $\tau_n(r)^2$, which, as seen in the proof of Lemma~\ref{size-incr}, converges to $K_r$. As a result, by Scheffé's theorem (see \cite[Theorem~16.12]{Billingsley95}), the lemma is proved if we can show that
\[ \frac1{\tau_n(r)^2}   \sum_{k=1}^\infty {\lambda}_k^{r-{\gamma}}\Biggl( (1-\mathrm{e}^{-{\lambda}_k^{\gamma}\Delta_n})-\frac{\mathrm{e}^{-2{\lambda}_k^{\gamma} t}}{2}(\mathrm{e}^{{\lambda}_k^{\gamma} \Delta_n}-1)^2\Biggr)\phi_k(y)^2 \to \frac{K_r}{|D|} \]
as $n\to\infty$, pointwise for all $y\in D$. But this follows from Lemma~\ref{convdens}
because the second term inside the parentheses is negligible as the calculations for $A^{n,i,2}_r$ in \eqref{second-term} demonstrate.
\end{proof}

The next lemma is the main step in the proof of Theorem~\ref{LLN2}.
\begin{lemma}\label{convP}
For $n\in\mathbb{N}$, $k\in\mathbb{N}$, ${\varepsilon}>0$, $-\frac d2\leq r<{\gamma}-\frac d2$, and $v\geq{\varepsilon}$, define
\begin{equation}\label{tildeahu}\begin{aligned} \widetilde a^{n,{\varepsilon}}_k(v)&:= \iint_{v-{\varepsilon}}^{v} (\mathrm{e}^{-{\lambda}_k^{\gamma}(v-s)}-\mathrm{e}^{-{\lambda}_k^{\gamma}(v-\Delta_n-s)}\mathds 1_{s\leq v-\Delta_n})\phi_k(y){\sigma}(v-{\varepsilon},y)\,W(\mathrm{d} s,\mathrm{d} y),\\
	\widetilde h^{n,{\varepsilon}}_r(v) &:=\Biggl(\sum_{k=1}^\infty {\lambda}_k^r(\widetilde a^{n,{\varepsilon}}_k(v))^2\Biggr)^{1/2}. \end{aligned}\end{equation}
Then, for any $q>0$,
\begin{equation}\label{convP-eq} \sup_{v\in[{\varepsilon},t]}\mathbb{E}\Biggl[\Biggl\lvert \frac{\widetilde h^{n,{\varepsilon}}_r(v)}{\tau_n(r)} - \biggl(\frac{K_r}{|D|}\int_D {\sigma}^2(v-{\varepsilon},y)\,\mathrm{d} y\biggr)^{1/2}\Biggr\rvert^q\Biggr] \to0\quad\text{as}\quad n\to\infty. \end{equation}
\end{lemma}
\begin{proof} For simplicity, we adopt the notation $\mathbb{E}_{v-{\varepsilon}}[\cdot] = \mathbb{E}[\cdot\,|\,{\cal F}_{v-{\varepsilon}}]$, and similarly ${\mathrm{Var}}_{v-{\varepsilon}}$ and ${\mathrm{Cov}}_{v-{\varepsilon}}$. As the moments of any order of $\widetilde h^{n,{\varepsilon}}_r(v) /\tau_n(r)$ are uniformly bounded in $n$ and $v$ (cf.\ \eqref{toshow2}) and ${\sigma}^2(v-{\varepsilon},y)$ has uniformly bounded moments of all orders by \eqref{mom}, it suffices to prove the statement for $q=1$. Therefore, the result will follow from Lemma \ref{lemapp}  (with ${\cal G} = {\cal F}_{v-{\varepsilon}}$) if we can show that, as $n\to\infty$,
\begin{align}
\label{conv-mean} \sup_{v\in[{\varepsilon},t]} \mathbb{E}\Biggl[\biggl\lvert\tau_n(r)^{-2}\mathbb{E}_{v-{\varepsilon}}[(\widetilde h^{n,{\varepsilon}}_r(v))^2] - \frac{K_r}{|D|}\int_D {\sigma}^2(v-{\varepsilon},y)\,\mathrm{d} y\biggr\rvert\Biggr]&\to0,\\
\sup_{v\in[{\varepsilon},t]} \mathbb{E}\Bigl[\tau_n(r)^{-4}{\mathrm{Var}}_{v-{\varepsilon}}[(\widetilde h^{n,{\varepsilon}}_r(v))^2]\Bigr]&\to0.
\label{conv-var}
\end{align}

Given a large number $L>0$, consider a partition $D=\bigcup_{\ell=1}^{L} D^L_\ell$ of $D$ into  pairwise disjoint measurable sets $D^L_l$ such that
\[ \lim_{L\to\infty} \sup_{\ell=1,\dots,L} \sup_{x,y\in D^L_\ell} |x-y| =0, \]
where $|\cdot|$ is an arbitrary norm on $\mathbb{R}^d$.
Picking an arbitrary point $y^L_\ell$ from each $D^L_\ell$, we can then
decompose the difference in \eqref{conv-mean}   into $E^{n,{\varepsilon},L,1}_r(v)+E^{n,{\varepsilon},L,2}_r(v)+E^{{\varepsilon},L,3}_r(v)$, where
\begin{equation}\label{D}\begin{aligned}
E^{n,{\varepsilon},L,1}_r(v)&:= \iint_0^v  \sum_{\ell=1}^L ({\sigma}^2(v-{\varepsilon},y)-{\sigma}^2(v-{\varepsilon},y^L_\ell))\mathds 1_{D^L_\ell}(y)\mathds 1_{s<{\varepsilon}}\,\Pi^n_{r,{\gamma}}(\mathrm{d} s,\mathrm{d} y),\\
E^{n,{\varepsilon},L,2}_r(v)&:=\sum_{\ell=1}^L {\sigma}^2(v-{\varepsilon},y^L_\ell) (\Pi^{n}_{r,{\gamma}}([0,{\varepsilon}]\times D^L_\ell)-\tfrac{K_r}{|D|}|D^L_\ell|),\\
E^{{\varepsilon},L,3}_r(v)&:=\frac{K_r}{|D|}\sum_{\ell=1}^L \int_D ({\sigma}^2(v-{\varepsilon},y^L_\ell) -{\sigma}^2(v-{\varepsilon},y))\mathds 1_{D^L_\ell}(y)\,\mathrm{d} y.
\end{aligned}\end{equation}
By \eqref{mom} and Lemma~\ref{Pinr}, we have $\sup_{v\in[{\varepsilon},t]}\mathbb{E}[|E^{n,{\varepsilon},L,2}_r(v)|]\to0$  as $n\to\infty$. Using the $L^2({\Omega})$-continuity of ${\sigma}$, one can further show that
\[ \sup_{v\in[{\varepsilon},t]}\mathbb{E}[|E^{{\varepsilon},L,3}_r(v)|]\to 0\qquad\text{and}\qquad \limsup_{n\to\infty} \sup_{v\in[{\varepsilon},t]} \mathbb{E}[|E^{n,{\varepsilon},L,1}_r(v)|]\to 0\] as $L\to\infty$, proving \eqref{conv-mean}.

Next, we turn to \eqref{conv-var}. Using the notation $\Delta^n\mathrm{e}_k(s)=\mathrm{e}^{-{\lambda}_k^{\gamma} s}-\mathrm{e}^{-{\lambda}_k^{\gamma}(s-\Delta_n)}\mathds 1_{s\geq\Delta_n}$, and the formula ${\mathrm{Cov}}(X^2,Y^2)=2{\mathrm{Cov}}(X,Y)^2$ for mean-zero bivariate Gaussian variables $X$ and $Y$ (which follows, for example, from \eqref{prodmom}), we derive
\begin{equation}\label{var-line1}\begin{aligned}
\frac{1}{\tau_n(r)^4}{\mathrm{Var}}_{v-{\varepsilon}}[(\widetilde h^{n,{\varepsilon}}_r(v))^2] &=\frac{1}{\tau_n(r)^4}\sum_{k,m=1}^\infty {\lambda}_k^r {\lambda}_m^r{\mathrm{Cov}}_{v-{\varepsilon}}\Bigl( \widetilde a^{n,{\varepsilon}}_k(v)^2,\widetilde a^{n,{\varepsilon}}_m(v)^2\Bigr)\\
&=\frac{2}{\tau_n(r)^4}\sum_{k,m=1}^\infty {\lambda}_k^r {\lambda}_m^r\biggl(\iint_0^{\varepsilon} \Delta^n\mathrm{e}_k(s)\Delta^n\mathrm{e}_m(s) \\
&\qquad\times \phi_k(y)\phi_m(y){\sigma}^2(v-{\varepsilon},y)\,\mathrm{d} s\,\mathrm{d} y\biggr)^2.
\end{aligned}\end{equation}

As $D$ is a bounded open set, we know
from the Vitali covering theorem (see \cite[Theorem~5.5.2]{Bogachev98}) that for small $\delta>0$, there are $L(\delta)\in\mathbb{N}$ and pairwise disjoint open balls $(B^\delta_\ell)_{\ell=1,\ldots, L(\delta)}$ of radius $\leq\delta$ such that $B^\delta_\ell \subseteq D$ and $\dist(B^\delta_\ell, \partial D)>\delta$ for every $i$ and such that $|D\setminus\bigcup_{\ell=1}^{L(\delta)} B^\delta_\ell| \downarrow 0$ as $\delta\downarrow 0$.
Denoting the center of $B^\delta_\ell$ by $z^\delta_\ell$
and using the elementary estimate $(x_1+x_2+x_3)^2\leq 3(x_1^2+x_2^2+x_3^2)$, we bound the terms in  \eqref{var-line1} by $3\sum_{i=1}^3 F^{n,{\varepsilon},\delta,i}_r(v)$, where
\begin{align*}
F^{n,{\varepsilon},\delta,1}_r(v)  &:=\frac{2}{\tau_n(r)^4}\sum_{k,m=1}^\infty {\lambda}_k^{r}{\lambda}_m^{r} \Biggl(\int_0^{{\varepsilon}}  \Delta^n\mathrm{e}_k(s)\Delta^n\mathrm{e}_m(s) \,\mathrm{d} s\\
&\qquad\times \int_D \phi_k(y)\phi_m(y)  \sum_{\ell=1}^{L(\delta)}\mathds 1_{B^\delta_\ell}(y)  \Bigl({\sigma}^2(v-{\varepsilon},y) -{\sigma}^2(v-{\varepsilon},z^\delta_\ell) \Bigr) \,\mathrm{d} y\Biggr)^2,\\
F^{n,{\varepsilon},\delta,2}_r(v) &:=\frac{2}{\tau_n(r)^4}\sum_{k,m=1}^\infty {\lambda}_k^{r}{\lambda}_m^{r}\Biggl( \int_0^{{\varepsilon}}  \Delta^n\mathrm{e}_k(s)\Delta^n\mathrm{e}_m(s)  \,\mathrm{d} s \\
&\qquad\times \sum_{\ell=1}^{L(\delta)}{\sigma}^2(v-{\varepsilon},z^\delta_\ell) \int_{B^\delta_\ell}\phi_k(y)\phi_m(y)\,\mathrm{d} y\Biggr)^2,\\
F^{n,{\varepsilon},\delta,3}_r(v) &:=\frac{2}{\tau_n(r)^4}\sum_{k,m=1}^\infty {\lambda}_k^{r}{\lambda}_m^{r}  \Biggl(\int_0^{{\varepsilon}} \Delta^n\mathrm{e}_k(s)\Delta^n\mathrm{e}_m(s) \,\mathrm{d} s \\
&\qquad\times\int_{D\setminus\bigcup_{\ell=1}^{L(\delta)} B^\delta_\ell} \phi_k(y)\phi_m(y)  {\sigma}^2(v-{\varepsilon},y) \,\mathrm{d} y\Biggr)^2.
\end{align*}
Because ${\sigma}$ has uniformly bounded moments of all orders,
\[ \sup_{t\in[0,T]}\sup_{\ell =1,\dots,L(\delta)}\sup_{y\in B^\delta_\ell} \mathbb{E}[|{\sigma}^2(t,y) -{\sigma}^2(t,z^\delta_i) |^2]^{\frac12} \lesssim w_4(\delta;T) \]
with the notation from \eqref{L2-cont}. So following the calculations in the proof of Lemma~\ref{size-incr}, we deduce that for every $t\in[0,T]$,
\[  \lim_{\delta\to0}\limsup_{n\to\infty} \sup_{v\in[{\varepsilon},t]} \mathbb{E}[|F^{n,{\varepsilon},\delta,1}_r(v)|]= 0.\]
By a similar argument and the Cauchy--Schwarz inequality,
\begin{align*} &\limsup_{n\to\infty} \sup_{v\in[{\varepsilon},t]}  \mathbb{E}[|F^{n,{\varepsilon},\delta,3}_r(v)|]\\
&\quad\lesssim \limsup_{n\to\infty} \Biggl( \frac{1}{\tau_n(r)^2} \sum_{k=1}^\infty {\lambda}_k^r \int_0^{{\varepsilon}} (\Delta^n\mathrm{e}_k(s))^2\,\mathrm{d} s\int_{D\setminus\bigcup_{\ell=1}^{L(\delta)} B^\delta_\ell} \phi_k(y)^2\,\mathrm{d} y \Biggr)^2\\
&\quad\lesssim \limsup_{n\to\infty}  \Biggl( \frac{1}{\tau_n(r)^2} \sum_{k=1}^\infty {\lambda}_k^{r-{\gamma}}(1-\mathrm{e}^{-{\lambda}_k^{\gamma}\Delta_n})  \int_{D\setminus\bigcup_{\ell=1}^{L(\delta)} B^\delta_\ell} \phi_k(y)^2\,\mathrm{d} y \Biggr)^2. \end{align*}
By Lemma~\ref{convdens} (ii), this is $\lesssim |D\setminus\bigcup_{\ell=1}^{L(\delta)} B^\delta_\ell|^2$, which converges to $0$ as $\delta\to0$.

For $F^{n,{\varepsilon},\delta,2}_r(v)$, with small $\delta >0$ fixed, we consider the terms where $k=m$ and the terms where $k\neq m$ separately. By \eqref{mom} and \eqref{calc}, the sum of the terms where $k=m$ has a first moment that is bounded by
\begin{equation}\label{eq3}\begin{aligned}
	&\frac{2}{\tau_n(r)^4}\sum_{k=1}^\infty {\lambda}_k^{2r} \biggl( \int_0^{\varepsilon} \Delta^n\mathrm{e}_k(s)^2\,\mathrm{d} s\biggr)^2\\
	&\quad\leq \frac{2}{\tau_n(r)^4}\sum_{k=1}^\infty {\lambda}_k^{2r-2{\gamma}}  \Bigl( 1-\mathrm{e}^{-{\lambda}_k^{\gamma}\Delta_n}  - \tfrac12\mathrm{e}^{-2{\lambda}_k^{\gamma}{\varepsilon}} (\mathrm{e}^{{\lambda}_k^{\gamma}\Delta_n}-1)^2 \Bigr)^2.
\end{aligned} \end{equation}
As in the proof of Lemma~\ref{size-incr}, the second part is of order $\tau_n(r)^{-4}\Delta_n^4$, which goes to $0$ as $n \to \infty$. Regarding the first part, we have
\begin{equation}\label{sum}
\frac{1}{\tau_n(r)^4}\sum_{k=1}^\infty {\lambda}_k^{2r-2{\gamma}} (1-\mathrm{e}^{-{\lambda}_k^{\gamma}\Delta_n})^2\leq \Biggl(\sup_{k\in\mathbb{N}} \frac{{\lambda}_k^{r-{\gamma}} (1-\mathrm{e}^{-{\lambda}_k^{\gamma}\Delta_n})}{\tau_n(r)^2}\Biggr)\frac{1}{\tau_n(r)^2}\sum_{k=1}^\infty {\lambda}_k^{r-{\gamma}} (1-\mathrm{e}^{-{\lambda}_k^{\gamma}\Delta_n}).
\end{equation}
The term after the supremum converges to $K_r$ by Lemma~\ref{convdens} (i). At the same time, if $r\neq -\frac d2$, then $0< \frac1{\gamma}({\gamma}-\frac d2-r)<1$. So choosing ${\varepsilon}>0$ small enough, we obtain
\begin{align*}
\sup_{k\in\mathbb{N}} \frac{{\lambda}_k^{r-{\gamma}} (1-\mathrm{e}^{-{\lambda}_k^{\gamma}\Delta_n})}{\tau_n(r)^2} \leq \sup_{k\in\mathbb{N}} \frac{{\lambda}_k^{r-{\gamma}}({\lambda}_k^{\gamma}\Delta_n)^{({\gamma}-d/2-r+{\varepsilon})/{\gamma}}}{\Delta_n^{({\gamma}-d/2-r)/{\gamma}}} = \Delta_n^{\frac{{\varepsilon}}{{\gamma}}} \sup_{k\in\mathbb{N}}{\lambda}_k^{-\frac d2+{\varepsilon}} = {\lambda}_1^{-\frac d2+{\varepsilon}}\Delta_n^{\frac {\varepsilon} {\gamma}} \to 0
\end{align*}
as $n\to\infty$. Similarly, if $r=-\frac d2$,
\begin{align*}
\sup_{k\in\mathbb{N}} \frac{{\lambda}_k^{-\frac d2-{\gamma}} (1-\mathrm{e}^{-{\lambda}_k^{\gamma}\Delta_n})}{\tau_n(-\frac d2)^2} \leq \sup_{k\in\mathbb{N}} \frac{{\lambda}_k^{-\frac d2-{\gamma}} {\lambda}_k^{\gamma}\Delta_n}{\Delta_n\lvert\log \Delta_n\rvert} = \frac1{\lvert\log \Delta_n\rvert} \sup_{k\in\mathbb{N}}{\lambda}_k^{-\frac d2} = {\lambda}_1^{-\frac d2}\frac1{\lvert\log \Delta_n\rvert} \to 0.
\end{align*}
We conclude that the right-hand side of \eqref{sum} converges to $0$ as $n\to\infty$.

Finally, let us discuss the terms in $F^{n,{\varepsilon},\delta,2}_r(v)$ when $k\neq m$. Applying the Cauchy--Schwarz inequality to the $\mathrm{d} s$-integral and realizing that, as before, the asymptotically relevant contribution of $\int_0^{\varepsilon} (\Delta^n\mathrm{e}_{j}(s))^2\,\mathrm{d} s$ is ${\lambda}_{j}^{-{\gamma}}(1-\mathrm{e}^{-{\lambda}_{j}^{\gamma}\Delta_n})$, where $j \in \{k,m \}$, we only have to analyze further the following expression:
\begin{align*} &\frac{4}{\tau_n(r)^4}\sum_{k=1}^\infty\sum_{m=k+1}^\infty {\lambda}_k^{r-{\gamma}}{\lambda}_m^{r-{\gamma}}(1-\mathrm{e}^{-{\lambda}_k^{\gamma}\Delta_n})(1-\mathrm{e}^{-{\lambda}_m^{\gamma}\Delta_n})\\
&\qquad\times\Biggl(    \sum_{\ell=1}^{L(\delta)}{\sigma}^2(v-{\varepsilon},z^\delta_\ell) \int_{B^\delta_\ell}\phi_k(y)\phi_m(y)\,\mathrm{d} y\Biggr)^2. \end{align*}
For any given $K\in\mathbb{N}$, we can further neglect all terms where $k\leq K$. Indeed, their sum is bounded in $L^1({\Omega})$ by a constant times
\[ \frac{\Delta_n}{\tau_n(r)^4}\sum_{k=1}^K \sum_{m=k+1}^\infty {\lambda}_k^r {\lambda}_m^{r-{\gamma}}(1-\mathrm{e}^{-{\lambda}_m^{\gamma}\Delta_n})\leq \frac{\Delta_n}{\tau_n(r)^2}\sum_{k=1}^K {\lambda}_k^r \cdot \frac{1}{\tau_n(r)^2}\sum_{m=1}^\infty {\lambda}_m^{r-{\gamma}}(1-\mathrm{e}^{-{\lambda}_m^{\gamma}\Delta_n}).\]
We have $\Delta_n/\tau_n(r)^2\to0$ and the second factor tends to $K_r$ by Lemma~\ref{convdens} (i). Thus, the previous display
converges to $0$ as $n\to\infty$, uniformly in $v\in[{\varepsilon},t]$.  As a consequence, only
\begin{equation}\label{ovF}\begin{aligned} \overline F^{n,{\varepsilon},\delta}_r(v)&:=\frac{4}{\tau_n(r)^4}\sum_{k=K+1}^\infty \sum_{m=k+1} {\lambda}_k^{r-{\gamma}}{\lambda}_m^{r-{\gamma}}(1-\mathrm{e}^{-{\lambda}_k^{\gamma}\Delta_n})(1-\mathrm{e}^{-{\lambda}_m^{\gamma}\Delta_n})\\
&\qquad\times\Biggl(    \sum_{\ell=1}^{L(\delta)}{\sigma}^2(v-{\varepsilon},z^\delta_\ell) \int_{B^\delta_\ell}\phi_k(y)\phi_m(y)\,\mathrm{d} y\Biggr)^2 \end{aligned}\end{equation}
remains to be considered, where $K$ will be chosen shortly.

If $r>-\frac d2$, we decompose $\overline F^{n,{\varepsilon},\delta}_r(v)$ into $\overline F^{n,{\varepsilon},\delta,\theta,1}_r(v)+\overline F^{n,{\varepsilon},\delta,\theta,2}_r(v)$, where the first (resp., second) term comprises the indices $(k,m)$ satisfying $K<k<m\leq (1+\theta)k$ (resp., $k>K$ and $m>(1+\theta)k$) and where $\theta\in(0,1)$ is some fixed number for the moment. We can bound the $\mathrm{d} y$-integral of $\phi_k\phi_m$ over $B^\delta_\ell$ in the previous display either by $1$ (using the Cauchy--Schwarz inequality) or by the estimate in Lemma~\ref{lem:innerprod} (recall that $\delta$ and the $B^\delta_\ell$ are fixed). Therefore, letting $\rho\in[0,1]$ be another parameter that we choose later, we obtain
\begin{align*} &\sup_{v\in[{\varepsilon},t]}\mathbb{E}[|\overline F^{n,{\varepsilon},\delta,\theta,2}_r(v)|] \\
&\quad\leq C_\delta \frac{4}{\tau_n(r)^4}\sum_{k=K+1}^\infty\sum_{m>(1+\theta)k}{\lambda}_k^{r-{\gamma}}{\lambda}_m^{r-{\gamma}}(1-\mathrm{e}^{-{\lambda}_k^{\gamma}\Delta_n})(1-\mathrm{e}^{-{\lambda}_m^{\gamma}\Delta_n})  \frac{{\lambda}_m^\rho {\lambda}_k^{\frac14 \rho}}{({\lambda}_m-{\lambda}_k)^{2\rho}} \end{align*}
for some constant $C_\delta$ that is independent of $n$. Thanks to \eqref{eq:Weyl}, we can write ${\lambda}_k=C_D(1+\eta_k)k^{2/d}$ for all $k\in\mathbb{N}$ with a sequence  $(\eta_k)_{k\in\mathbb{N}}$   that converges to $0$. In particular, if $m>(1+\theta)k$, then
\[
  {\lambda}_m \geq \frac{1+\eta_m}{1+\eta_k} (1+\theta)^{\frac{2}{d}} {\lambda}_k,
  \qquad\text{therefore}\qquad  {\lambda}_m - {\lambda}_k \geq K(m,k,\theta) \, {\lambda}_k,
\]
where
\[
    K(m,k,\theta):= \frac{1+\eta_m}{1+\eta_k} (1+\theta)^{\frac{2}{d}} - 1 = \frac{1+\eta_m}{1+\eta_k} \left(\left[ (1+\theta)^{\frac{2}{d}} - 1\right] + 1 - \frac{1+\eta_k}{1+\eta_m}\right).
\]
Since $(1+\theta)^{\frac{2}{d}} - 1 \geq \frac{2\theta}{d}$, we have $K(m,k,\theta) \geq \frac{\theta}{2 d}$, for $m>k\geq K=K(\theta)\in \mathbb{N}$, where $K(\theta)$ is large enough so that
\[
   \frac{1+\eta_m}{1+\eta_k} \geq \frac12,\qquad
   1 - \frac{1+\eta_k}{1+\eta_m} \geq -\frac{\theta}{d}.
\]
Therefore,
 \begin{align*} &\sup_{v\in[{\varepsilon},t]}\mathbb{E}[|\overline F^{n,{\varepsilon},\delta,\theta,2}_r(v)|] \\
 &\quad\lesssim  \frac{1}{\tau_n(r)^4}\sum_{k=K+1}^\infty\sum_{m>(1+\theta)k}{\lambda}_k^{r-{\gamma}}{\lambda}_m^{r-{\gamma}}(1-\mathrm{e}^{-{\lambda}_k^{\gamma}\Delta_n})(1-\mathrm{e}^{-{\lambda}_m^{\gamma}\Delta_n})   {\lambda}_m^\rho {\lambda}_k^{\frac14 \rho}{\lambda}_k^{-2 \rho}\\
&\quad\lesssim \frac{1}{\tau_n(r)^4}\sum_{k=1}^\infty {\lambda}_k^{r-\frac74 \rho-{\gamma}} (1-\mathrm{e}^{-{\lambda}_k^{\gamma}\Delta_n}) \sum_{m=1}^\infty {\lambda}_m^{ r-{\gamma}+\rho}(1-\mathrm{e}^{-{\lambda}_m^{\gamma}\Delta_n}). \end{align*}
Let $\rho>0$ be small enough such that $r_1=r-\frac 74 \rho>-\frac d2$ and $r_2=r+\rho<{\gamma}-\frac d2$. Then, by Lemma~\ref{convdens} (i), the sum over $k$ is of order $\tau_n(r_1)^2$, while the sum over $m$ is of order $\tau_n(r_2)^2$. As a consequence,
the formula \eqref{taunr} for $\tau_n(r)$ shows
\begin{equation}\label{F2-2}   \lim_{n\to\infty} \sup_{v\in[{\varepsilon},t]}\mathbb{E}[|\overline F^{n,{\varepsilon},\delta,\theta,2}_r(v)|] \lesssim \lim_{n\to\infty} \frac{\tau_n(r_1)^2\tau_n(r_2)^2}{\tau_n(r)^4} = \lim_{n\to\infty} \Delta_n^{\frac1{\gamma}(\frac74\rho-\rho)}= 0. \end{equation}

Next, we consider $\overline F^{n,{\varepsilon},\delta,\theta,1}_r(v)$. Bounding the second moment of the sum in parentheses in \eqref{ovF} by a constant and using $1-\mathrm{e}^{-x}\leq x^{\frac1{\gamma}({\gamma}-\frac d2-r)}$, we obtain
\begin{equation}\label{F1}\begin{aligned}
\sup_{v\in[{\varepsilon},t]} \mathbb{E}[|\overline F^{n,{\varepsilon},\delta,\theta,1}_r(v)|] &\lesssim \frac{1}{\tau_n(r)^4}\sum_{k=1}^\infty \sum_{k<m\leq (1+\theta)k}  {\lambda}_k^{r-{\gamma}}{\lambda}_m^{r-{\gamma}}(1-\mathrm{e}^{-{\lambda}_k^{\gamma}\Delta_n})(1-\mathrm{e}^{-{\lambda}_m^{\gamma}\Delta_n})\\
&\lesssim \frac{1}{\tau_n(r)^2}\sum_{k=1}^\infty {\lambda}_k^{r-{\gamma}} (1-\mathrm{e}^{-{\lambda}_k^{\gamma}\Delta_n})  \sum_{k<m\leq (1+\theta)k} m^{-1}\\
&\leq \frac{1}{\tau_n(r)^2}\sum_{k=1}^\infty {\lambda}_k^{r-{\gamma}} (1-\mathrm{e}^{-{\lambda}_k^{\gamma}\Delta_n}) (\log((1+\theta)k)-\log k)\\
&=\log(1+\theta)\frac{1}{\tau_n(r)^2}\sum_{k=1}^\infty {\lambda}_k^{r-{\gamma}} (1-\mathrm{e}^{-{\lambda}_k^{\gamma}\Delta_n})  \\
&\to K_r \log(1+\theta)
\end{aligned}\end{equation}
by Lemma~\ref{convdens} (i). As a consequence,
\[ \lim_{\theta\to0}\limsup_{n\to\infty}\sup_{v\in[{\varepsilon},t]} \mathbb{E}[|\overline F^{n,{\varepsilon},\delta,\theta,1}_r(v)|]  =0,  \]
which completes the proof of \eqref{conv-var} for $r>-\frac d2$.

If $r=-\frac d2$, we go back to \eqref{ovF} and let $\overline F^{n,{\varepsilon},\delta,\theta,1}_{-d/2}(v)$ (resp., $\overline F^{n,{\varepsilon},\delta,\theta,2}_{-d/2}(v)$) contain the indices $(k,m)$ with $K<k<m\leq k^\theta$ (resp., $k>K$ and $m>k^\theta$), where $\theta>1$ and $K=K(\theta)>0$ is chosen such that $|\eta_k|\leq \frac12$ for all $k>K$ and $m^{\frac2d}\geq 6 m^{\frac 2{d\theta}}$ for all $m>K^\theta$. Then $k>K$ and $m>k^\theta$ imply
\begin{align*} {\lambda}_m-{\lambda}_k &= C_D(1+\eta_m)m^{\frac 2d} - C_D(1+\eta_k)k^{\frac 2d} \geq \tfrac12C_Dm^{\frac 2d} -\tfrac32 C_Dm^{\frac2{d\theta}}\geq\tfrac14 C_Dm^{\frac 2d}\\
&\begin{cases} \gtrsim {\lambda}_m,  \\  \geq\frac14 C_Dk^{\frac 2d \theta}\gtrsim {\lambda}_k^\theta. \end{cases}\end{align*}

So arguing as above (with $\rho=1$) and bounding $({\lambda}_m-{\lambda}_k)^2\gtrsim {\lambda}_m{\lambda}_k^\theta$, we derive
 \begin{align*} &\sup_{v\in[{\varepsilon},t]}\mathbb{E}[|\overline F^{n,{\varepsilon},\delta,\theta,2}_{-d/2}(v)|]\\
  &\quad\lesssim  \frac{1}{\tau_n(-\frac d2)^4}\sum_{k=K+1}^\infty\sum_{m>k^\theta}{\lambda}_k^{-\frac d2-{\gamma}}{\lambda}_m^{-\frac d2-{\gamma}}(1-\mathrm{e}^{-{\lambda}_k^{\gamma}\Delta_n})(1-\mathrm{e}^{-{\lambda}_m^{\gamma}\Delta_n})     {\lambda}_k^{\frac14 -\theta} \\
&\quad\lesssim  \frac{1}{\tau_n(-\frac d2)^4}\sum_{k= 1}^\infty{\lambda}_k^{-\frac d2+\frac14 -\theta-{\gamma}}(1-\mathrm{e}^{-{\lambda}_k^{\gamma}\Delta_n})\sum_{m=1}^\infty{\lambda}_m^{-\frac d2-{\gamma}}(1-\mathrm{e}^{-{\lambda}_m^{\gamma}\Delta_n}).
\end{align*}
By Lemma \ref{convdens}, the sum over $m$ is ${\cal O}(\tau_n(-\frac d2)^2)$, whereas, as seen in the proof of Lemma \ref{size-incr}, the sum over $k$ is ${\cal O}(\tau_n(-\frac d2+\frac14-\theta))={\cal O}(\Delta_n)$. Therefore,
 \begin{align*} \lim_{n\to\infty}\sup_{v\in[{\varepsilon},t]}\mathbb{E}[|\overline F^{n,{\varepsilon},\delta,\theta,2}_{-d/2}(v)|] &\lesssim\lim_{n\to\infty} \frac{\Delta_n}{\tau_n(-\frac d2)^2} = \lim_{n\to\infty} \frac1{\lvert\log \Delta_n\rvert} = 0. \end{align*}

Analogously to \eqref{F1}, we further have
\begin{align*} &\sup_{v\in[{\varepsilon},t]} \mathbb{E}[|\overline F^{n,{\varepsilon},\delta,\theta,1}_{-d/2}(v)|] \\
&\quad\lesssim \frac{\theta-1}{\tau_n(-\frac d2)^2\lvert\log\Delta_n\rvert}\sum_{k=1}^\infty {\lambda}_k^{-\frac d2-{\gamma}} (1-\mathrm{e}^{-{\lambda}_k^{\gamma}\Delta_n})\log k\\
&\quad\lesssim \frac{\theta-1}{\lvert\log\Delta_n\rvert^2}\sum_{k=1}^{[\Delta_n^{-d/(2{\gamma})}]} k^{-1} \log k +\frac{\theta-1}{\Delta_n\lvert\log\Delta_n\rvert^2}\sum_{k=[\Delta_n^{-d/(2{\gamma})}]+1}^\infty k^{-1-\frac2d{\gamma}}  \log k. \end{align*}
Bounding $k^{-\frac 2d {\gamma}}\leq(1+[\Delta_n^{-d/(2{\gamma})}])^{-\frac 2d {\gamma}}\leq \Delta_n$, we conclude that
\begin{align*} \limsup_{n\to\infty} \sup_{v\in[{\varepsilon},t]} \mathbb{E}[|\overline F^{n,{\varepsilon},\delta,\theta,1}_{-d/2}(v)|] &\leq \lim_{n\to\infty} \frac{2(\theta-1)}{\lvert\log\Delta_n\rvert^2}\log \Delta_n^{-\frac{d}{2{\gamma}}}(1+\log \Delta_n^{-\frac{d}{2{\gamma}}})\lesssim \theta-1. \end{align*}
The proof of \eqref{conv-var} when $r=-\frac d2$ is completed by sending $\theta\downarrow 1$.
\end{proof}

\begin{comment}
To complete the proof of \eqref{conv-var},  we must show that the integral in \eqref{var-line2}, with ${\sigma}^2(v-{\varepsilon},y){\sigma}^2(v-{\varepsilon},w)$ replaced by $\mathds 1_{[y_1,y_2]}(y)\mathds 1_{[w_1,w_2]}(w)$, converges to $0$ for any $y_1,y_2,w_1,w_2\in[0,\pi]$. By additivity of $v \mapsto \int_0^v \cdots$, if suffices to consider $y_1=w_1=0$. Then, because of the estimate
$$\left|\int_0^{y_2} \sin(ky)\sin(my) \,\mathrm{d} y\right|=\left|\frac{k\sin(my_2)\cos(ky_2)-m\sin(ky_2)\cos(my_2)}{m^2-k^2}\right| \leq \frac{1}{m-k},$$
the term \eqref{var-line2}, modified as described above, is bounded by
\begin{align*}
\frac{1}{\tau_n(r)^4} \sum_{(k,m)\in\mathbb{N}^2\colon m>k} \frac{k^{2r}m^{2r}}{(m-k)^2} \left(\int_0^{\varepsilon} \Delta^n\mathrm{e}_k(s)\Delta^n \mathrm{e}_m(s)\,\mathrm{d} s\right)^2.
\end{align*}
Using the Cauchy--Schwarz inequality on the $\mathrm{d} s$-integral, we can further bound this term by (cf.\ \eqref{eq3})
\begin{align*}
&\frac{1}{\tau_n(r)^4}\sum_{(k,m)\in\mathbb{N}^2\colon m>k}   \frac{k^{2r}m^{2r}}{(m-k)^2} \left( \frac{1-\mathrm{e}^{-k^2\Delta_n}}{k^2} - \mathrm{e}^{-2k^2{\varepsilon}}\frac{(\mathrm{e}^{k^2\Delta_n}-1)^2}{2k^2}\right)\\
&\quad\times\left( \frac{1-\mathrm{e}^{-m^2\Delta_n}}{m^2} - \mathrm{e}^{-2m^2{\varepsilon}}\frac{(\mathrm{e}^{m^2\Delta_n}-1)^2}{2m^2}\right).
\end{align*}
When we expand the product, as soon as we take at least one of the second terms in each of the parentheses as a factor, the resulting expression converges to $0$ as $n\to\infty$. The last remaining term is given by (define $j=m-k$ in the following)
$$
\frac{1}{\tau_n(r)^2} \sum_{k=1}^\infty k^{2r} \frac{1-\mathrm{e}^{-k^2\Delta_n}}{k^2}\times \frac{1}{\tau_n(r)^2} \sum_{j=1}^\infty \frac{(k+j)^{2r}}{j^2(k+j)^2}(1-\mathrm{e}^{-(k+j)^2\Delta_n}).
$$
Let $\theta\in(\frac12-r,\min(1,1-r))$ if $r>-\frac12$ and let $\theta=1$ if $r=-\frac12$. Then the previous line is
$$
  \leq \frac{1}{\tau_n(r)^2} \sum_{k=1}^\infty k^{2r} \frac{1-\mathrm{e}^{-k^2\Delta_n}}{k^2} \times \frac{1}{\tau_n(r)^2} \sum_{j=1}^\infty \frac{1}{j^2} (k+j)^{2r-2 + 2\theta} \Delta_n^\theta.

$$
Since $2r-2 + 2\theta \leq 0$ by the choice of $\theta$, we replace $k+j$ by $j$ in the second series to see that this is bounded above by
$$
  \frac{1}{\tau_n(r)^2} \sum_{k=1}^\infty k^{2r} \frac{1-\mathrm{e}^{-k^2\Delta_n}}{k^2} \times \frac{\Delta_n^\theta}{\tau_n(r)^2} \sum_{j=1}^\infty j^{2r - 4 +2\theta}.
$$
As shown in Lemma~\ref{size-incr}, the first factor, which is equal to $\tau_n(r)^{-2}A^{n,1}_r$, converges to $K_r$ as $n\to\infty$. Concerning the second factor, since $2r - 4 +2\theta \leq -2$, the series in $j$ converges, and $\Delta_n^\theta \tau_n(r)^{-2}  \to 0$ as $n\to\infty$ by our choice of $\theta$.

\end{comment}

\begin{lemma}\label{LLN-condexp}
For every $t>0$,
\[
\lim_{{\varepsilon}\to0} \limsup_{n\to\infty} \mathbb{E}\Biggl[\Biggl\lvert \Delta_n\sum_{i=1}^{[t/\Delta_n]}  f\biggl( \frac{\widetilde h^{n,i,{\varepsilon}}_r}{\tau_n(r)}\biggr) -V^r_f(u,t)\Biggr\rvert\Biggr] = 0.\]
\end{lemma}
\begin{proof}
We write the difference in the previous line as a sum of four terms:
\begin{align*}
G^{n,{\varepsilon},1}_r &:= \Delta_n \sum_{i=[{\varepsilon}/\Delta_n]+1}^{[t/\Delta_n]} \Biggl\{  f\Biggl(\frac{\widetilde h^{n,i,{\varepsilon}}_r}{\tau_n(r)}\Biggr)  - f\Biggl( \sqrt{\frac{K_r}{|D|} \int_D{\sigma}^2(i\Delta_n-{\varepsilon},y)\,\mathrm{d} y}\Biggr) \Biggr\},\\
G^{n,{\varepsilon},2}_r &:=  \int_{\Delta_n[{\varepsilon}/\Delta_n]}^{\Delta_n[t/\Delta_n]}  \sum_{i=[{\varepsilon}/\Delta_n]+1}^{[t/\Delta_n]} \Biggl\{ f\Biggl( \sqrt{\frac{K_r}{|D|} \int_D {\sigma}^2(i\Delta_n-{\varepsilon},y)\,\mathrm{d} y}\Biggr)\\
&\qquad- f\Biggl( \sqrt{\frac{K_r}{|D|} \int_D {\sigma}^2(v-{\varepsilon},y)\,\mathrm{d} y}\Biggr) \Biggr\}\mathds 1_{[(i-1)\Delta_n,i\Delta_n)}(v) \,\mathrm{d} v,\\
G^{n,{\varepsilon},3}_r &:= \int_{\Delta_n[{\varepsilon}/\Delta_n]}^{\Delta_n[t/\Delta_n]} \Biggl\{  f\Biggl( \sqrt{\frac{K_r}{|D|} \int_D{\sigma}^2(v-{\varepsilon},y)\,\mathrm{d} y}\Biggr)- f\Biggl( \sqrt{\frac{K_r}{|D|} \int_D {\sigma}^2(v,y)\,\mathrm{d} y}\Biggr)\Biggr\}\,\mathrm{d} v,\\
G^{n,{\varepsilon},4}_r &:=\Delta_n\sum_{i=1}^{[{\varepsilon}/\Delta_n]}  f\Biggl( \frac{\widetilde h^{n,i,{\varepsilon}}_r}{\tau_n(r)}\Biggr) -\int_0^{\Delta_n[{\varepsilon}/\Delta_n]} f\Biggl( \sqrt{\frac{K_r}{|D|} \int_D{\sigma}^2(v,y)\,\mathrm{d} y}\Biggr)\,\mathrm{d} v\\
&\qquad-\int_{\Delta_n[t/\Delta_n]}^{t} f\Biggl( \sqrt{\frac{K_r}{|D|} \int_D {\sigma}^2(v,y)\,\mathrm{d} y}\Biggr)\,\mathrm{d} v.
\end{align*}

Employing the estimate \eqref{fdiff} and Lemma~\ref{convP} and arguing as in the beginning of the proof of Lemma~\ref{LLN-discr} to get \eqref{toshow1} for $\widetilde h^{n,i,{\varepsilon}}_r$, we can show that $\mathbb{E}[|G^{n,{\varepsilon},1}_r|]\to0$ as $n\to\infty$. Using the $L^2({\Omega})$-continuity of ${\sigma}$ in addition (see Remark~\ref{L2-rem}), it follows similarly to the proof of Lemma~\ref{LLN-discr}, that $\mathbb{E}[|G^{n,{\varepsilon},2}_r|]\to0$ as $n\to\infty$ and $\limsup_{n\to\infty} \mathbb{E}[|G^{n,{\varepsilon},3}_r|] \to 0$ as ${\varepsilon}\to0$. Furthermore, as the reader may easily verify by using \eqref{toshow2} for $\widetilde h^{n,i,{\varepsilon}}_r$ and \eqref{mom},  also $\lim_{{\varepsilon}\to0}\limsup_{n\to\infty} \mathbb{E}[|G^{n,{\varepsilon},4}_r|]= 0$.
\end{proof}

\begin{proof}[Proof of Corollary~\ref{Holder}] With the same arguments as for \eqref{toshow-2} and \eqref{toshow2}, one can show that for all $p,T>0$, there is a constant $C=C(p,r,{\gamma},{\sigma},T)>0$ such that for all $s,t\in[0,T]$
\[ \mathbb{E}[\|u(t,\cdot)-u(s,\cdot)\|^p_{H_r}]^{\frac1p} \leq C|t-s|^{{\alpha}(r)}  \]
($\leq C|t-s|^{{\alpha}(r)-{\varepsilon}}$ when $r=-\frac d2$, where ${\varepsilon}$ can be taken arbitrarily small). So the first part of the corollary follows from Kolmogorov's continuity theorem (see \cite[Theorem~3.3]{DaPrato14}). For the second part, if ${\sigma}\equiv1$, observe from Lemma~\ref{size-incr} that for all $0<T_1<T_2<\infty$, there is $C'=C'(r,{\gamma},T_1,T_2)$ such that for all $s,t\in[T_1,T_2]$,
\[ \mathbb{E}[\|u(t,\cdot)-u(s,\cdot)\|^2_{H_r}]\geq C'|t-s|^{2{\alpha}(r)}.  \]
Hence, the second assertion  can be obtained with the same arguments as in \cite[Theorem~4]{Dalang05}.
\end{proof}

\subsection{Comparison between Theorem~\ref{LLN}  and Theorem~\ref{LLN2}}\label{compare}

We conclude this paper by comparing Theorem~\ref{LLN} (where $r<-\frac d2$) and Theorem~\ref{LLN2} (where $-\frac d2 \leq r<{\gamma}-\frac d2$) in two remarks.
\begin{remark}\label{rem-tight} Let us explain why in Theorem~\ref{LLN2}, we did not study functionals as general as in Theorem~\ref{LLN}. As Lemma~\ref{tight} shows, the family
\[ {\cal U}' := \Biggl\{ \frac{u(i\Delta_n,\cdot)-u((i-1)\Delta_n,\cdot)}{\tau_n(r)}: \quad n\in\mathbb{N},~i=1,\dots,[T/\Delta_n]\Biggr\} \]
is tight for $r<-\frac d2$ and $\tau_n(r)=\sqrt{\Delta_n}$, which makes it possible to prove Lemma \ref{remove-drift}.

By contrast, if $-\frac d2\leq r<{\gamma}-\frac d2$ and ${\sigma}$ is, say, identically $1$, then there is no normalizing sequence $\tau_n(r)$ such that ${\cal U}'$ becomes tight in $H_r$ without all limit points being zero. Indeed,
it follows from Lemma~\ref{size-incr} that for $-\tfrac d2 \leq r < {\gamma}-\tfrac d2$, the normalization $\tau_n(r)$ must be chosen as in \eqref{taunr} (or at least, of the same order) to avoid degenerate limits. Since ${\cal U}'$ is a collection of centered Gaussian vectors on $H_r$, we know from \cite[Example~3.8.13~(iv)]{Bogachev98} that if ${\cal U}'$ were tight, then the series $\sum_{k=1}^\infty \mathbb{E}[\langle \Delta^n_i u/\tau_n(r), b_k\rangle_{H_r}^2]$ would converge uniformly in $n$ and $i$, where $b_k$ are the orthonormal basis functions from \eqref{basis_b} and $\Delta^n_i u= u(i\Delta_n,\cdot)-u((i-1)\Delta_n,\cdot)$. In other words, we would have
\begin{equation}\label{false} \lim_{m\to\infty} \sup_{n\in\mathbb{N}}\sup_{i=1,\ldots,[T/\Delta_n]}\sum_{k=m}^\infty \mathbb{E}\Biggl[\biggl\langle \frac{\Delta^n_i u}{\tau_n(r)}, b_k\biggr\rangle_{H_r}^2\Biggr] = 0. \end{equation}
But this is not true because from the calculations in \eqref{series-calc}, we deduce that
\begin{align*}&\sum_{k=m}^\infty \mathbb{E}\Biggl[\biggl\langle \frac{\Delta^n_i u}{\tau_n(r)}, b_k\biggr\rangle_{H_r}^2\Biggr]\\
&\quad= \frac{1}{\tau_n(r)^2}\Biggl(\sum_{k=m}^\infty {\lambda}_k^{r-{\gamma}}(1-\mathrm{e}^{-{\lambda}_k^{\gamma}\Delta_n})  - \frac12\sum_{k=m}^\infty {\lambda}_k^{r-{\gamma}} \mathrm{e}^{-2{\lambda}_k^{\gamma} i\Delta_n}  (\mathrm{e}^{{\lambda}_k^{\gamma}\Delta_n}-1)^2\Biggr). \end{align*}
The expression within the parentheses is exactly the last line of \eqref{calc}, except that summation starts at $k=m$. By \eqref{second-term}, if $i\Delta_n\geq{\varepsilon}$ for some ${\varepsilon}>0$, then the second part, divided by $\tau_n(r)^2$, vanishes as $n\to\infty$. In the notation of the proof of Lemma~\ref{size-incr}, the first part equals
\begin{equation}\label{eq4} \frac{1}{\tau_n(r)^2}A^{n,1}_r- \frac{1}{\tau_n(r)^2}\sum_{k=1}^{m-1}{\lambda}_k^{r-{\gamma}}(1-\mathrm{e}^{-{\lambda}_k^{\gamma}\Delta_n}). \end{equation}
For any fixed value of $m\in\mathbb{N}$, the second term is bounded by $(\sum_{k=1}^{m-1} {\lambda}_k^r)\Delta_n/\tau_n(r)^2$, which converges to $0$ as $n\to\infty$. At the same time, the first term converges to the constant $K_r$ of \eqref{Kr}, as seen in the proof of the lemma. This is a nonzero limit that is independent of the value of $m$, so \eqref{false} must be false and therefore ${\cal U}'$ cannot be tight.
\end{remark}

\begin{remark}\label{difference} There are further important differences between Theorem~\ref{LLN} and Theorem~\ref{LLN2}. Let us take ${\sigma}\equiv1$  and consider quadratic variations (i.e., $p=2$ in \eqref{Vnp}, so $f(x) = x^2$) in the following to simplify the discussion.
\begin{longlist}
	\item If $r<-\frac d2$, Theorem~\ref{LLN} gives a genuine law of large numbers in the following sense: As $n$ tends to infinity, most of the squared normalized increments $\|\Delta^n_i u\|^2_{H_r}/\tau_n(r)^2$ (where ``most'' means that  $i\Delta_n\geq{\varepsilon}$ for some fixed but arbitrary ${\varepsilon}>0$) have equal mean (namely $K_r$ by Lemma~\ref{size-incr}) and bounded variance (see \eqref{toshow-2}). Moreover, they are only weakly correlated to each other, so a classical $L^2({\Omega})$-argument gives the law of large numbers; cf.\ Lemma~\ref{LLN-core-1}.
	\item By contrast, if $-\frac d2\leq r<{\gamma}-\frac d2$, Theorem~\ref{LLN2} is  a \emph{degenerate} law of large numbers. Here, not only do we have $\mathbb{E}[\|\Delta^n_i u\|_{H_r}^2]/\tau_n(r)^2 \to K_r$ (i.e., convergence of means) for most values of $i$, but we actually have $\mathbb{E}[|\|\Delta^n_i u\|_{H_r}^2/\tau_n(r)^2 - K_r|]$ (i.e., convergence in mean); cf.\ Lemma~\ref{convP}. In the case $r<-\frac d2$, only the \emph{sum} $\Delta_n\sum_{i=1}^{[t/\Delta_n]} \|\Delta^n_i u\|_{H_r}^2/\tau_n(r)^2$ converges in mean, not the individual summands. In particular, no averaging or cancellation argument as in Lemma~\ref{LLN-core-1}, which would be typical of a genuine law of large numbers, is needed for $-\frac d2\leq r<{\gamma}-\frac d2$. This also explains why if we compare \eqref{2pvar-add} and \eqref{p-var-add}, then only in \eqref{2pvar-add} is there a factor (namely $2^pB_p(x_1,\ldots,x_p)$) that is related to moment formulas.
\end{longlist}
\end{remark}

\begin{appendix}
	\setcounter{equation}{0}
\setcounter{theorem}{0}
\section*{}\label{appn}
We start with some analytic preliminaries needed for the proof of Proposition~\ref{exist}.
\begin{lemma}\label{A1} Let $D$ be an open bounded subset of $\mathbb{R}^d$. Then the eigenvalues and eigenfunctions of $-\Delta$ have the following properties:
\begin{longlist}
\item The eigenvalues $({\lambda}_n)_{n\in\mathbb{N}}$ satisfy \emph{Weyl's law}, that is,
\begin{equation}\label{eq:Weyl} \lim_{n\to\infty} \frac{{\lambda}_n}{n^{2/d}} =C_D:=\frac{4\pi {\Gamma}(1+\frac d2)^{2/d}}{|D|^{2/d}}. \end{equation}
\item For all $n\in\mathbb{N}$,
\begin{equation}\label{eq:Linfty} \|\phi_n\|_{L^\infty(D)} \leq \biggl(\frac{2\mathrm{e} {\lambda}_n}{\pi d}\biggr)^{d/4}. \end{equation}
\item For every $x\in D$,  the following asymptotics are valid as ${\lambda}\to\infty$ (where $f(x)\sim g(x)$ means $f(x)/g(x)\to1$):
\begin{itemize}
\item for $r<-\frac d 2$,
\begin{equation}\label{est1} \sum_{k:\,{\lambda}_k>{\lambda}} {\lambda}_k^r |\phi_k(x)|^{2} \sim -\frac{1}{(4\pi)^{d/2}{\Gamma}(1+\frac d2)(1+\frac 2d r)} {\lambda}^{\frac d2+r}; \end{equation}
\item for $r>-\frac d2$,
\begin{equation}\label{est2} \sum_{k:\,{\lambda}_k\leq{\lambda}} {\lambda}_k^r |\phi_k(x)|^{2} \sim \frac{1}{(4\pi)^{d/2}{\Gamma}(1+\frac d2)(1+\frac 2d r)} {\lambda}^{\frac d2+r}; \end{equation}
\item for $r=-\frac d2$,
\begin{equation}\label{est3} \sum_{k:\,{\lambda}_k\leq{\lambda}} {\lambda}_k^r |\phi_k(x)|^{2} \sim \frac{1}{(4\pi)^{d/2}{\Gamma}(\frac d2)}\log {\lambda}. \end{equation}
\end{itemize}
\item If $D$ additionally has the cone property, then, for any $r<-\frac d2$, there is ${\varepsilon}>0$ such that
\begin{equation}\label{summable} \sup_{x\in D} \sum_{k=1}^\infty {\lambda}_k^{r} |\phi_k(x)|^{2+{\varepsilon}}<\infty. \end{equation}
\end{longlist}
\end{lemma}
\begin{proof} Part (i) is well known; see \cite[Theorem~1.11]{Birman80} or \cite{Geisinger14}. Part (ii) is shown in \cite[Lemma~3.1]{Davies74}. For (iii), we use \cite[Equation~(0.7)]{Garding53}, which shows, for every $x\in D$, that
\[ V(t):=V(t;x):=\sum_{1\leq k\leq t} |\phi_k(x)|^2 \sim \frac{t}{|D|} \]
as $t\to\infty$.
Therefore, by Abel's summation formula (see \cite[Theorem~4.2]{Apostol76}),
\begin{align*} \sum_{1\leq k\leq t} k^{\frac 2d r} |\phi_k(x)|^2 &= t^{\frac 2d r}V(t) -\frac{2}{d}r \int_1^t s^{\frac2d r-1} V(s)\,\mathrm{d} s \sim  \frac{t^{\frac2d r+1}}{|D|} - \frac{2r}{d\lvert D\rvert}\int_1^t s^{\frac2d r}\,\mathrm{d} s\\
&\sim \begin{cases} \frac{1}{(\frac 2d r+1)\lvert D\rvert}t^{\frac 2d r+1}&\text{if } r>-\frac d2, \\ \frac{1}{|D|}\log t &\text{if } r=-\frac d2.\end{cases} \end{align*}
As ${\lambda}_k \sim C_Dk^{2/d}$ by \eqref{eq:Weyl}, we derive
\[ \sum_{k:\, {\lambda}_k\leq{\lambda}} {\lambda}_k^r |\phi_k(x)|^{2} \sim C_D^r \sum_{1\leq k\leq C_D^{-d/2}{\lambda}^{d/2}} k^{\frac 2d r}|\phi_k(x)|^2 \sim \begin{cases} \frac{C_D^{- d/2}}{(\frac 2d r+1)\lvert D\rvert}{\lambda}^{r+\frac d2}&\text{if } r>-\frac d2, \\ \frac{dC_D^{-d/2}}{2|D|}\log {\lambda} &\text{if } r=-\frac d2,\end{cases} \]
as ${\lambda}\to\infty$, which becomes \eqref{est2} and \eqref{est3} after simplifying the constants. If $r<-\frac d2$, we obtain, with similar arguments,
\begin{align*} \sum_{k:\, k> t} k^{\frac 2d r} |\phi_k(x)|^2 &= \lim_{K\to\infty} \Biggl( K^{\frac 2d r}V(K) - t^{\frac2d r} V(t) -\frac{2}{d}r \int_t^K s^{\frac2d r-1} V(s)\,\mathrm{d} s\Biggr)\\
&= - t^{\frac2d r} V(t) -\frac{2}{d}r \int_t^\infty s^{\frac2d r-1} V(s)\,\mathrm{d} s\\
&\sim  -\frac{t^{\frac2d r+1}}{|D|} - \frac{2r}{d\lvert D\rvert}\int_t^\infty s^{\frac2d r}\,\mathrm{d} s= - \frac{1}{(\frac 2d r+1)\lvert D\rvert}t^{\frac 2d r+1} \end{align*}
and consequently,
\[ \sum_{{\lambda}_k>{\lambda}} {\lambda}_k^r |\phi_k(x)|^{2} \sim C_D^r \sum_{k> C_D^{-d/2}{\lambda}^{d/2}} k^{\frac 2d r}|\phi_k(x)|^2 \sim  - \frac{C_D^{- d/2}}{(\frac 2d r+1)\lvert D\rvert}{\lambda}^{r+\frac d2},   \]
which is \eqref{est1}.

For (iv), we use (ii), \cite[Theorem~8.2]{Agmon65}, and the hypothesis $r <-\frac d2$ to deduce
\begin{align*}
\sum_{{\lambda}_k\geq{\lambda}} {\lambda}_k^{r} |\phi_k(x)|^{2+{\varepsilon}} &=\sum_{k=0}^\infty \sum_{{\lambda}_k\in [2^n{\lambda}, 2^{n+1}{\lambda})} {\lambda}_k^{r} |\phi_k(x)|^{2+{\varepsilon}} \lesssim \sum_{n=0}^\infty \sum_{{\lambda}_k\in [2^n{\lambda}, 2^{n+1}{\lambda})} {\lambda}_k^{r+ \frac d 4 {\varepsilon}} \phi_k(x)^2 \\
&\lesssim {\lambda}^{r+ \frac d 4 {\varepsilon}}\sum_{n=0}^\infty 2^{nr} \sum_{{\lambda}_k\leq 2^{n+1}{\lambda}}\phi_k(x)^2 \lesssim {\lambda}^{r+ \frac d 4 {\varepsilon}} \sum_{n=0}^\infty 2^{nr} (2^{n+1}{\lambda})^{\frac d2} \lesssim {\lambda}^{\frac d2 +r+ \frac d 4 {\varepsilon}},
\end{align*}
which converges to $0$ as ${\lambda}\to\infty$, uniformly in $x\in D$, if ${\varepsilon}$ is small enough. This proves \eqref{summable}.
\end{proof}

\begin{proof}[Proof of Proposition~\ref{exist}] The existence of a jointly measurable and adapted solution of \eqref{SHE} will be a consequence of Theorem 4.2.1 in \cite{Dalang20} once we show that
\begin{equation}\label{L2kernel}   \int_0^t \mathrm{d} s \, \sup_{x\in\mathbb{R}^d} \int_D \mathrm{d} y\, g^2(s;x,y) < \infty.
\end{equation}

According to \cite[Theorem 8.2]{Agmon65}, because $D$ satisfies the cone property, there is a finite constant $C>0$ such that for all $x \in D$ and all $t \geq 0$,
\begin{equation}\label{eA.8}
   \widetilde V(t,x) := \sum_{k:\, \lambda_k \leq t} |\phi_k(x)|^2 \leq C t^{d/2}.
\end{equation}

Using the orthonormal property of $(\phi_k)_{k\in\mathbb{N}}$, we see that
\begin{align*} \int_D g(s;x,y)^2\,\mathrm{d} y
&= \sum_{k=1}^\infty \phi_k(x)^2\mathrm{e}^{-2{\lambda}_k^{\gamma} s}
= \int_{\lambda_1-}^{+\infty} \mathrm{e}^{-2u^\gamma s} \, \widetilde V(\mathrm{d} u,x).
\end{align*}
By Abel's summation formula (see \cite[Theorem~4.2]{Apostol76}), this is equal to
\[
   - \int_{\lambda_1}^{+\infty} \widetilde V(u,x) (-2\gamma u^{\gamma-1} s) \mathrm{e}^{-2u^\gamma s}\, \mathrm{d} u \leq 2\gamma C \int_{\lambda_1}^{+\infty} u^{\frac{d}{2} + \gamma -1} s \mathrm{e}^{2 u^\gamma s} \, \mathrm{d} u.
\]
Therefore,
\[
   \int_0^t \mathrm{d} s \, \sup_{x\in\mathbb{R}^d} \int_D \mathrm{d} y\, g^2(s;x,y) \leq 2\gamma C \int_{\lambda_1}^{+\infty} \mathrm{d} u\, u^{\frac{d}{2} + \gamma -1}\int_0^t \mathrm{d} s \, s\, \mathrm{e}^{2 u^\gamma s}.
\]
Since $\int_0^\infty se^{-as}\mathrm{d} s = a^{-2}$ ($a>0$), this is bounded above by
\[
   2\gamma C \int_{\lambda_1-}^{+\infty} u^{\frac{d}{2} + \gamma -1} (2 u^\gamma)^{-2}\, \mathrm{d} u = \tfrac{\gamma C}{2} \int_{\lambda_1-}^{+\infty}  u^{\frac{d}{2} - \gamma -1}\, \mathrm{d} u < \infty
\]
because $\tfrac{d}{2} - \gamma <0$. This implies \eqref{L2kernel}.

In order to establish the $L^p$-continuity of $u$, according to \cite[Section 4.2.2]{Dalang20}, it suffices to show that
\begin{align}\label{L2cont}
\lim_{(\tau,h)\to0}  \int_0^t \int_D (g(t+\tau-s;x+h,y)-g(t-s;x,y))^2\,\mathrm{d} y\,\mathrm{d} s &=0,\\
\lim_{\tau \downarrow 0}  \int_t^{t+\tau}\int_D g^2(t+\tau-s;x,y)\,\mathrm{d} y\,\mathrm{d} s &= 0, \label{L2cont2}\end{align}
for all $t>0$ and $x\in D$. Once $L^2$-continuity is established, predictability follows immediately \cite[Proposition 2]{Dalang98}.

In order to establish \eqref{L2cont}, because $\int_0^t (\mathrm{e}^{-{\lambda}_k^{\gamma} (s+\tau)}-\mathrm{e}^{-{\lambda}_k^{\gamma} s})^2\,\mathrm{d} s = (1-\mathrm{e}^{-{\lambda}_k^{\gamma} \tau})^2(1-\mathrm{e}^{-2{\lambda}_k^{\gamma} t})/(2{\lambda}_k^{\gamma})$, one easily checks that the integral in \eqref{L2cont} is bounded by
\begin{align*} &2\int_0^t\int_D (g(s+\tau;x+h,y)-g(s;x+h,y))^2\,\mathrm{d} y\,\mathrm{d} s \\
&\quad\quad+2\int_0^t\int_D (g(s;x+h,y)-g(s;x,y))^2\,\mathrm{d} y\,\mathrm{d} s\\
&\quad \leq \sum_{k=1}^\infty \frac{(1-\mathrm{e}^{-{\lambda}_k^{\gamma} \tau})^2}{{\lambda}_k^{\gamma}}\phi_k(x+h)^2 + \sum_{k=1}^\infty \frac{(\phi_k(x+h)-\phi_k(x))^2}{{\lambda}_k^{\gamma}}. \end{align*}
As $(\tau,h)\to0$, the terms within both sums converge to $0$ for each $k\in\mathbb{N}$. Moreover, $\sum_{k=1}^\infty {\lambda}_k^{-{\gamma}}\delta_k$ is a finite measure on $\mathbb{N}$ by Lemma~\ref{A1} (i) and the hypothesis ${\gamma}>\frac d2$. So by Proposition~4.12 and the discussion before Lemma~4.10 in \cite{Kallenberg02}, \eqref{L2cont} is implied by $\sup_{x\in D} \sum_{k=1}^\infty (2{\lambda}_k^{\gamma})^{-1}|\phi_k(x)|^{2+{\varepsilon}}<\infty$, which was shown in \eqref{summable} for some small ${\varepsilon}>0$.

Concerning \eqref{L2cont2}, the same arguments show that
\[ \int_t^{t+\tau}\int_D g^2(t+\tau-s;x,y)\,\mathrm{d} y\,\mathrm{d} s = \int_0^\tau \int_D g^2(s;x,y)\,\mathrm{d} y\,\mathrm{d} s = \sum_{k=1}^\infty \frac{1-\mathrm{e}^{-{\lambda}_k^{\gamma} \tau}}{2{\lambda}_k^{\gamma}} \phi_k(x)^2 \to 0 \]
as $\tau \downarrow 0$.
\end{proof}

\begin{proof}[Proof of Proposition~\ref{uinHr}] Observe from \eqref{akt} and the orthogonality of $(\phi_k)_{k\in\mathbb{N}}$ that $(a_k(t))_{k\in\mathbb{N}}$ is a sequence of independent centered Gaussian random variables. Thus, by \cite[Theorem~4.17]{Kallenberg02} and \cite[Theorem~22.7]{Billingsley95}, the almost-sure finiteness of $\|u(t,\cdot)\|^2_{H_r} = \sum_{k=1}^\infty {\lambda}_k^r (a_k(t))^2$ is equivalent to the convergence in probability of $\sum_{k=1}^n {\lambda}_k^{r/2} a_k(t)$ as $n\to\infty$. By Gaussianity, this is in turn equivalent to the summability of the variances, that is, to $r<{\gamma}-\frac d2$ by \eqref{eq:L2finite}.
\end{proof}

\begin{lemma}\label{convdens} Suppose that $D\subseteq\mathbb{R}^d$ is open and bounded and that $-\frac d2\leq r<{\gamma}-\frac d2$. Let
\begin{align*}
H_n&:=\frac1{\tau_n(r)^2}   \sum_{k=1}^\infty {\lambda}_k^{r-{\gamma}} (1-\mathrm{e}^{-{\lambda}_k^{\gamma}\Delta_n}),\\
H_n(y)&:= \frac1{\tau_n(r)^2}   \sum_{k=1}^\infty {\lambda}_k^{r-{\gamma}} (1-\mathrm{e}^{-{\lambda}_k^{\gamma}\Delta_n})\phi_k(y)^2,\qquad y\in D, \end{align*}
with $\tau_n(r)$ from \eqref{taunr}. Further recall the definition of $K_r$ from \eqref{Kr}.
\begin{longlist}
\item As $n\to\infty$, we have $H_n\to K_r$ and $H_n(y)\to K_r/|D|$ for every $y\in D$.
\item If $D$ additionally satisfies the cone property and $D_0$ is a measurable subset of $D$, then
\begin{equation}\label{limsup} \limsup_{n\to\infty} \int_{D_0} H_n(y)\,\mathrm{d} y \leq C|D_0|\end{equation}
with a constant $C\in(0,\infty)$ that does not depend on $D_0$.
\end{longlist}
\end{lemma}
\begin{proof} Consider (i) first. We only show the claim for $H_n(y)$; the assertion for $H_n$ follows analogously once we replace $\phi_k(y)^2$ by $1$ in the subsequent arguments (in this case, there would be an additional factor $|D|$ in \eqref{eA.12}, which explains that the limit of $H_n$ differs from the limit of $H_n(y)$ also by a factor of $|D|$). If $r\neq -\frac d2$, we fix $y\in D$ and
let
\begin{align*} U({\lambda})&:=U({\lambda};y):=\sum_{k:\,{\lambda}_k\leq {\lambda}} \phi_k(y)^2,\qquad {\lambda}>0,\\
k(z)&:=z^{{\gamma}-r}(1-\mathrm{e}^{-z^{-{\gamma}}}),\qquad z>0,\\
M(z):=M(z;y)&:= \int_0^\infty k(z/{\lambda})\, \mathrm{d} U({\lambda}) = z^{{\gamma}-r}\sum_{k=1}^\infty {\lambda}_k^{r-{\gamma}}(1-\mathrm{e}^{-{\lambda}_k^{\gamma}/z^{{\gamma}}})\phi_k(y)^2 ,\qquad z>0.
\end{align*}
so that $M(z)$ is the Mellin--Stieltjes transform of $U$ \cite[(4.0.3)]{Bingham87}. Then
\begin{equation}\label{Hn} H_n(y)=\Delta_n^{-\frac1{\gamma}({\gamma}-\frac d2 -r)}\sum_{k=1}^\infty {\lambda}_k^{r-{\gamma}} (1-\mathrm{e}^{-{\lambda}_k^{\gamma}\Delta_n})\phi_k(y)^2 = \Delta_n^{\frac {d}{2{\gamma}}}M(\Delta_n^{-\frac1{\gamma}}). \end{equation}

We use an Abelian theorem from \cite{Bingham87} to determine the behavior of $M(z)$ as $z\to\infty$. To this end, observe from \cite[Equation~(0.7)]{Garding53} and \eqref{eq:Weyl} that
\begin{equation}\label{eA.12} U({\lambda})\sim \sum_{1\leq k\leq C_D^{-d/2}{\lambda}^{d/2}} \phi_k(y)^2\sim (2\pi)^{-d}\frac{\pi^{\frac d2}}{{\Gamma}(1+\frac d2)}{\lambda}^{\frac d2},\qquad {\lambda}\to\infty. \end{equation}
 Moreover, since ${\lambda}_1>0$, $U({\lambda})$ is identically $0$ for $0<{\lambda}<{\lambda}_1$. Next, pick ${\alpha},\beta>0$ such that $(-r)\vee 0<{\alpha}<\frac d2 <\beta<{\gamma} -r$, which is possible because $-\frac d 2<r<{\gamma}-\frac d2$ by assumption. Because $z\mapsto k(z)$ is continuous on $(0,\infty)$, increasing at $0$ and eventually decreasing for large $z$, there is $n_0\in\mathbb{N}$ such that
\begin{align*} &\sum_{n\in\mathbb{Z}} (\mathrm{e}^{-{\alpha} n}\vee\mathrm{e}^{-\beta n})\sup_{\mathrm{e}^n \leq z\leq \mathrm{e}^{n+1}} k(z)\\
&\quad \leq \sum_{n= -\infty}^{-n_0-1} \mathrm{e}^{-\beta n}\mathrm{e}^{(n+1)({\gamma}-r)} + \sum_{n=-n_0}^{n_0} (\mathrm{e}^{-{\alpha} n}\vee\mathrm{e}^{-\beta n}) \sup_{z\in(0,\infty)}k(z)+\sum_{n=n_0+1}^\infty \mathrm{e}^{-{\alpha} n}\mathrm{e}^{-nr}\\
&\quad<\infty.\end{align*}
We conclude from \cite[Theorem~4.4.2]{Bingham87} that
\[ M(z)\sim \frac d2(2\pi)^{-d}\frac{\pi^{\frac d2}}{{\Gamma}(1+\frac d2)}\check k(\tfrac d2)z^{\frac d2},\qquad z\to\infty, \]
where, by a substitution $z^{\gamma}\mapsto y$ and integration by parts,
\begin{align*} \check k(\tfrac d2)&:=\int_0^\infty z^{\frac d2-1}k(z^{-1})\,\mathrm{d} z = \int_0^\infty z^{\frac d2-1}z^{r-{\gamma}}(1-\mathrm{e}^{-z^{\gamma}})\,\mathrm{d} z=\frac1{\gamma} \int_0^\infty y^{\frac{r}{{\gamma}}-2+\frac {d}{2{\gamma}}}(1-\mathrm{e}^{-y}) \,\mathrm{d} y \\
&=\frac{-1}{r-{\gamma}+\frac {d}{2}}\int_0^\infty  y^{\frac{r}{{\gamma}}-1+\frac {d}{2{\gamma}}}\mathrm{e}^{-y}\,\mathrm{d} y=\frac{-1}{r-{\gamma}+\frac {d}{2}}{\Gamma}(\tfrac{r}{{\gamma}}+\tfrac {d}{2{\gamma}})
\end{align*}
 Recalling the formula for $K_r$ from \eqref{Kr}, we obtain
\[ M(z)\sim \frac d2(2\pi)^{-d}\frac{\pi^{\frac d2}}{{\Gamma}(1+\frac d2)}\frac1{{\gamma}-\frac d2-r}{\Gamma}(\tfrac{r}{{\gamma}}+\tfrac {d}{2{\gamma}})z^{\frac d2} = \frac{K_r}{|D|}z^{\frac d2},\qquad z\to\infty. \]
Inserting this into \eqref{Hn} yields the desired result.

If $r=-\frac d2$, we split $H_n(y)$ into $H_{n,1}(y)+H_{n,2}(y)+H_{n,3}(y)$, where
\begin{align*}
H_{n,1}(y)&:=\frac1{\lvert\log \Delta_n\rvert}   \sum_{{\lambda}_k\leq \Delta_n^{-1/{\gamma}}} {\lambda}_k^{-\frac d2} \phi_k(y)^2,\\
H_{n,2}(y)&:=\frac1{\Delta_n\lvert\log \Delta_n\rvert}  \sum_{{\lambda}_k\leq \Delta_n^{-1/{\gamma}}} {\lambda}_k^{-\frac d2-{\gamma}} (1-\mathrm{e}^{-{\lambda}_k^{\gamma}\Delta_n}-{\lambda}_k^{\gamma}\Delta_n)\phi_k(y)^2,\\
H_{n,3}(y)&:=\frac1{\Delta_n\lvert\log \Delta_n\rvert}  \sum_{{\lambda}_k> \Delta_n^{-1/{\gamma}}} {\lambda}_k^{-\frac d2-{\gamma}} (1-\mathrm{e}^{-{\lambda}_k^{\gamma}\Delta_n})\phi_k(y)^2.
\end{align*}
Because $|1-\mathrm{e}^{-x}-x|\leq \frac{x^2}{2}$ for all $x\geq0$, we have from \eqref{est2},
\begin{align*}
H_{n,2}(y)&\lesssim \frac{\Delta_n}{\lvert\log \Delta_n\rvert} \sum_{{\lambda}_k\leq \Delta_n^{-1/{\gamma}}} {\lambda}_k^{{\gamma}-\frac d2} \phi_k(y)^2 \lesssim \frac{1}{\lvert\log \Delta_n\rvert} \to 0
\end{align*}
as $n\to\infty$. Similarly, using \eqref{est1} and the simple bound $1-\mathrm{e}^{-x}\leq 1$, we obtain
\begin{align*}
H_{n,3}(y)&\lesssim \frac{1}{\Delta_n\lvert\log \Delta_n\rvert} \sum_{{\lambda}_k> \Delta_n^{-1/{\gamma}}} {\lambda}_k^{-\frac d2-{\gamma}} \phi_k(y)^2 \lesssim \frac{1}{\lvert\log \Delta_n\rvert} \to 0.
\end{align*}
Finally, \eqref{est3} and \eqref{Kr} show $H_{n,1}(y)\to ({\gamma}(4\pi)^{\frac d2}{\Gamma}(\frac d2))^{-1}=\frac{K_r}{D}$, which completes the proof of (i).

For (ii), let $\widetilde V(t,x)$ be defined in \eqref{eA.8}, so that
\[
   H_n(y) = \frac{1}{\tau_n(r)^2} \int_{{\lambda}_1-}^\infty t^{r-{\gamma}} (1- \mathrm{e}^{-t^{\gamma} \Delta_n})) \, \widetilde V(\mathrm{d} t,y).
\]
By Abel summation (see \cite[Theorem~4.2]{Apostol76}), we have
\begin{align*}
   H_n(y) &= \frac{1}{\tau_n(r)^2} \lim_{t \to \infty} \Bigl[\widetilde V(t,y) t^{r-{\gamma}} (1- \mathrm{e}^{-t^{\gamma} \Delta_n})\Bigr] - \frac{1}{\tau_n(r)^2} \int_{\lambda_1}^\infty \widetilde V(t,y) \Bigl( t^{r-{\gamma}} (1- \mathrm{e}^{-t^{\gamma} \Delta_n})\Bigr)'\, \mathrm{d} t \\
   &\lesssim \frac{1}{\tau_n(r)^2} \int_{\lambda_1}^\infty (f_{n,1}(t) + f_{n,2}(t))\, \mathrm{d} t,
\end{align*}
where
\[
   f_{n,1}(t):= (\gamma - r)\, t^{\frac{d}{2}+r-{\gamma} -1} (1- \mathrm{e}^{-t^{\gamma} \Delta_n}) ,\qquad f_{n,2}(t):= -{\gamma}\, t^{\frac{d}{2}+r-1}\, \mathrm{e}^{-t^{\gamma} \Delta_n},
\]
and we have used \eqref{eA.8}. The right-hand side no longer depends on $y$, so the integral over $D_0$ contributes a factor $|D_0|$. Using the change of variables $u = t^{\gamma} \Delta_n$, we see that for $i=1,2$,
\[
   \lim_{n \to \infty } \frac{1}{\tau_n(r)^2} \int_{\lambda_1}^\infty f_{n,i}(t)\, \mathrm{d} t
\]
exists and is finite. This completes the proof of (ii).
\end{proof}

\begin{lemma}\label{lem:innerprod}
\begin{longlist}
\item Fix $d \geq 2$. Let $D\subseteq\mathbb{R}^d$ be a bounded connected open set with a piecewise smooth boundary (in the sense of \cite[Definition~1.17]{Hezari18}) and let $D_0$ be an open set with a smooth boundary satisfying $\overline D_0 \subseteq D$. Then there is $C\in(0,\infty)$, only depending on $D$ and $D_0$, such that for all $k> \ell$,
\[ \Biggl\lvert\int_{D_0} \phi_k(x)\phi_\ell(x) \,\mathrm{d} x\Biggr\rvert \leq C\frac{{\lambda}_k^{\frac12}{\lambda}_\ell^{\frac18}}{{\lambda}_k-{\lambda}_\ell}.\]
\item Let $d=1$ and let $D \subseteq\mathbb{R}$ be an interval and let $D_0$ be a subinterval of $D$. Then there is $C\in(0,\infty)$, only depending on $D$ and $D_0$, such that for all $k> \ell$,
\[
  \Biggl\lvert\int_{D_0} \phi_k(x)\phi_\ell(x) \,\mathrm{d} x\Biggr\rvert \leq C\frac{{\lambda}_k^{\frac12}}{{\lambda}_k-{\lambda}_\ell}.
\]
\end{longlist}
\end{lemma}
\begin{proof} For (i),   we can use the eigenfunction properties of $\phi_k$ and $\phi_\ell$ and  Green's second identity to get
\begin{align*} ({\lambda}_k-{\lambda}_\ell)\int_{D_0} \phi_k(x)\phi_\ell(x)\,\mathrm{d} x &= \int_{D_0} (\phi_k(x)\Delta\phi_\ell(x)-\phi_\ell(x) \Delta \phi_k(x))\,\mathrm{d} x\\
&=\int_{\partial D_0} (\phi_k(x)\partial_\nu\phi_\ell(x)-\phi_\ell(x) \partial_\nu \phi_k(x))\,S(\mathrm{d} x),
\end{align*}
where $S$ is the surface measure and $\partial_\nu$ is the derivative in direction of the outward pointing unit normal field of the boundary $\partial D_0$. Applying the Cauchy--Schwarz inequality, we obtain
\begin{align*} ({\lambda}_k-{\lambda}_\ell) \Biggl\lvert\int_{D_0} \phi_k(x)\phi_\ell(x)\,\mathrm{d} x  \Biggl\rvert&\leq \Biggl(\int_{\partial D_0} \phi_k(x)^2\,S(\mathrm{d} x)\Biggr)^{\frac12}\Biggr(\int_{\partial D_0} (\partial_\nu\phi_\ell(x))^2\,S(\mathrm{d} x)\Biggr)^{\frac12}\\
&\qquad+\Biggl(\int_{\partial D_0} \phi_\ell(x)^2\,S(\mathrm{d} x)\Biggr)^{\frac12}\Biggr(\int_{\partial D_0} (\partial_\nu\phi_k(x))^2\,S(\mathrm{d} x)\Biggr)^{\frac12}.
\end{align*}
Since $d \geq 2$, by \cite[Corollary~1.8]{Hezari18} and \cite[Theorem~1.1]{Christianson15}, it follows that
\begin{align*} ({\lambda}_k-{\lambda}_\ell) \Biggl\lvert\int_{D_0} \phi_k(x)\phi_\ell(x)\,\mathrm{d} x  \Biggl\rvert&\lesssim{\lambda}_k^{\frac18}{\lambda}_\ell^{\frac12}+{\lambda}_\ell^{\frac18}{\lambda}_k^{\frac12}\lesssim {\lambda}_\ell^{\frac18}{\lambda}_k^{\frac12},
\end{align*}
which implies the assertion of statement (i).

In the case of (ii), since $d=1$, we may as well assume that $D = (0,\pi)$. In this case, $\lambda_k = k^2$ and $\phi_k(x) = \sin(kx)$. Assuming that $D_0 = (y_1,y_2)$, we consider
\[
   \int_{y_1}^{y_2} \sin(k x) \sin(\ell x)\, \mathrm{d} x.
\]
By additivity of $v \mapsto \int_0^v \cdots$, it suffices to consider $y_1=w_1=0$. The estimate
\[\Biggl\lvert\int_0^{y_2} \sin(ky)\sin(\ell y) \,\mathrm{d} y\Biggr\rvert=\Biggl\lvert\frac{k\sin(\ell y_2)\cos(ky_2)-\ell\sin(ky_2)\cos(my_2)}{k^2-\ell^2}\Biggr\rvert \leq \frac{2\lambda_k^{\frac12}}{\lambda_k-\lambda_\ell}
\]
completes the proof of the lemma.
\end{proof}

\begin{lemma} \label{lemapp}
Let ${\cal G}$ be a sub-$\sigma$-field of ${\cal F}$, $I$ be an index set, and for $v\in I$, let $(Y_n(v))_{n \in \mathbb{N}}$ and $Z(v)$ be nonnegative random variables such that
\begin{equation} \label{e0app}
   \sup_{v\in I} \sup_{n\in \mathbb{N}} \mathbb{E}[Y_n(v)^6] < \infty \qquad\mbox{and}\qquad \sup_{v\in I} \mathbb{E}[Z(v)^6] < \infty.
\end{equation}
Suppose that
\begin{equation} \label{e1app}
   \lim_{n \to \infty}\sup_{v\in I} \mathbb{E}\biggl[ \Bigl\lvert \mathbb{E}[Y_n(v)^2 \mid {\cal G}] - Z(v)\Bigr\rvert\biggr] = 0\qquad{and}\qquad
   \lim_{n \to \infty} \sup_{v\in I} \mathbb{E}\Bigl[{\mathrm{Var}}(Y_n(v)^2 \mid {\cal G})\Bigr] = 0.
\end{equation}
Then
\begin{equation} \label{e21app}
   \lim_{n \to \infty} \sup_{v\in I} \mathbb{E}\biggl[\Bigl\lvert Y_n(v)  - \sqrt{Z(v)} \Bigr\rvert\biggr] = 0.
\end{equation}
\end{lemma}

\begin{proof} Recalling that
$$
   \mathbb{E}\Bigl[{\mathrm{Var}}(Y_n(v)^2 \mid {\cal G})\Bigr] = \mathbb{E}\biggl[\Bigl(Y_n(v)^2 - \mathbb{E}[Y_n(v)^2 \mid {\cal G}]\Bigr)^2\biggr],
$$
we deduce from \eqref{e1app} that
\begin{equation} \label{e3app}
   \lim_{n\to \infty} \sup_{v\in I} \Vert Y_n(v)^2 - \mathbb{E}[Y_n(v)^2 \mid {\cal G}] \Vert_{L^2({\Omega})} = 0.
\end{equation}
By \eqref{e1app},
$$
   \lim_{n\to \infty} \sup_{v\in I} \Vert \mathbb{E}[Y_n(v)^2 \mid {\cal G}] - Z(v) \Vert_{L^1({\Omega})} = 0.
$$
Because of the bounded moments assumption \eqref{e0app}, we get
\begin{equation} \label{e5app}
   \lim_{n\to \infty} \sup_{v\in I} \Vert \mathbb{E}[Y_n(v)^2 \mid {\cal G}] - Z(v) \Vert_{L^2({\Omega})}= 0.
\end{equation}
From \eqref{e3app} and \eqref{e5app}, we see that
\begin{equation}\label{e6app}
   \lim_{n\to \infty} \sup_{v\in I} \Vert Y_n(v)^2 - Z(v) \Vert_{L^2({\Omega})} = 0.
\end{equation}
Since $Y_n(v)$ and $Z(v)$ are nonnegative,
$$
   \mathbb{E}\biggl[\Bigl\lvert Y_n(v) - \sqrt{Z(v)}\Bigr\rvert 1_{\{Y_n(v) + \sqrt{Z(v)} \leq {\varepsilon}\}}\biggr] \leq {\varepsilon}.
$$
Therefore,
\begin{align*}
   \mathbb{E}\biggl[\Bigl\lvert Y_n(v) - \sqrt{Z(v)}\Bigr\rvert\biggr] &\leq \mathbb{E}\biggl[\Bigl\lvert Y_n(v) - \sqrt{Z(v)}\Bigr\rvert 1_{\{Y_n(v) + \sqrt{Z(v)} > {\varepsilon}\}} \biggr] + {\varepsilon}\\
   &\leq \mathbb{E}\Biggl[\Bigl\lvert Y_n(v) - \sqrt{Z(v)}\Bigr\rvert \frac{Y_n(v) + \sqrt{Z(v)}}{{\varepsilon}} 1_{\{Y_n(v) + \sqrt{Z(v)} > {\varepsilon}\}} \Biggr] + {\varepsilon} \\
   &\leq \frac{1}{{\varepsilon}} \mathbb{E}\biggl[\Bigl\lvert Y_n(v)^2 - Z(v)\Bigr\rvert\biggr]  + {\varepsilon}.
\end{align*}
From this bound and \eqref{e6app}, we immediately deduce \eqref{e21app}.
\end{proof}

\end{appendix}

\section*{Acknowledgements}

The second author was supported in part by the Swiss National Foundation for Scientific Research.