Regularity and convergence analysis in Sobolev and Hölder spaces for generalized Whittle–Matérn fields

We analyze several types of Galerkin approximations of a Gaussian random field Z:D×Ω→R\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\mathscr {Z}:\mathscr {D}\times \varOmega \rightarrow \mathbb {R}$$\end{document} indexed by a Euclidean domain D⊂Rd\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\mathscr {D}\subset \mathbb {R}^d$$\end{document} whose covariance structure is determined by a negative fractional power L-2β\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$L^{-2\beta }$$\end{document} of a second-order elliptic differential operator L:=-∇·(A∇)+κ2\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$L:= -\nabla \cdot (A\nabla ) + \kappa ^2$$\end{document}. Under minimal assumptions on the domain D\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\mathscr {D}$$\end{document}, the coefficients A:D→Rd×d\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$A:\mathscr {D}\rightarrow \mathbb {R}^{d\times d}$$\end{document}, κ:D→R\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\kappa :\mathscr {D}\rightarrow \mathbb {R}$$\end{document}, and the fractional exponent β>0\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\beta >0$$\end{document}, we prove convergence in Lq(Ω;Hσ(D))\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$L_q(\varOmega ; H^\sigma (\mathscr {D}))$$\end{document} and in Lq(Ω;Cδ(D¯))\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$L_q(\varOmega ; C^\delta (\overline{\mathscr {D}}))$$\end{document} at (essentially) optimal rates for (1) spectral Galerkin methods and (2) finite element approximations. Specifically, our analysis is solely based on H1+α(D)\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$H^{1+\alpha }(\mathscr {D})$$\end{document}-regularity of the differential operator L, where 0<α≤1\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$0<\alpha \le 1$$\end{document}. For this setting, we furthermore provide rigorous estimates for the error in the covariance function of these approximations in L∞(D×D)\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$L_{\infty }(\mathscr {D}\times \mathscr {D})$$\end{document} and in the mixed Sobolev space Hσ,σ(D×D)\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$H^{\sigma ,\sigma }(\mathscr {D}\times \mathscr {D})$$\end{document}, showing convergence which is more than twice as fast compared to the corresponding Lq(Ω;Hσ(D))\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$L_q(\varOmega ; H^\sigma (\mathscr {D}))$$\end{document}-rate. We perform several numerical experiments which validate our theoretical results for (a) the original Whittle–Matérn class, where A≡IdRd\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$A\equiv \mathrm {Id}_{\mathbb {R}^d}$$\end{document} and κ≡const.\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\kappa \equiv {\text {const.}}$$\end{document}, and (b) an example of anisotropic, non-stationary Gaussian random fields in d=2\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$d=2$$\end{document} dimensions, where A:D→R2×2\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$A:\mathscr {D}\rightarrow \mathbb {R}^{2\times 2}$$\end{document} and κ:D→R\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\kappa :\mathscr {D}\rightarrow \mathbb {R}$$\end{document} are spatially varying.

1. Introduction 1.1.Motivation and background.By virtue of their practicality owing to the full characterization by their mean and covariance structure, Gaussian random fields (GRFs for short) are popular models for many applications in spatial statistics and uncertainty quantification, see, e.g., [4,7,26,33,35].As a result, several methodologies in these disciplines require the efficient simulation of GRFs at unstructured locations in various possibly non-convex Euclidean domains, and this topic has been intensively discussed in both areas, spatial statistics and computational mathematics, see, e.g., [2,8,9,13,17,19,24,30].In particular, sampling from non-stationary GRFs, for which methods based on circulant embedding are inapplicable, has become a central topic of current research, see, e.g., [2,9,17].
In order to capture both stationary and non-stationary GRFs, a new class of random fields has been introduced in [26], which is based on the following observation made by P. Whittle [40]: A GRF Z on D := R d with covariance function of Matérn type solves the fractional-order stochastic partial differential equation (SPDE for short) where ∆ denotes the Laplacian, dW is white noise on R d , and κ > 0, β > d /4 are constants which determine the practical correlation length and the smoothness of the field.In [26] this relation has been exploited to formulate generalizations of Matérn fields, the generalized Whittle-Matérn fields, by considering the SPDE (1.1) for non-stationary differential operators L (e.g., by allowing for a spatially varying coefficient κ : D → R) on bounded domains D ⊂ R d , d ∈ {1, 2, 3}.Note that the covariance structure of a GRF is uniquely determined by its covariance operator, in this case given by the negative fractional-order differential operator L −2β .Furthermore, for the case 2β ∈ N, approximations based on a finite element discretization have been proposed in [26].Subsequently, a computational approach which allows for arbitrary fractional exponents β > d /4 has been suggested in [2,3].To this end, a sinc quadrature combined with a Galerkin discretization of the differential operator L is applied to the Balakrishnan integral representation of the fractional-order inverse L −β .
In this work, we investigate Sobolev and Hölder regularity of generalized Whittle-Matérn fields and we perform a rigorous error analysis in these norms for several Galerkin approximations, including the sinc-Galerkin approximations of [2,3].Specifically, we consider a GRF Z β : D × Ω → R, indexed by a Euclidean domain D ⊂ R d , whose covariance operator is given by the negative fractional power L −2β of a second-order elliptic differential operator L : D(L) ⊆ L 2 (D) → L 2 (D) in divergence form with Dirichlet boundary conditions, formally given by Here, we solely assume that D ⊂ R d has a Lipschitz boundary, κ ∈ L ∞ (D), and that A ∈ L ∞ D; R d×d is symmetric and uniformly positive definite.For a sequence Z β N N ∈N of Galerkin approximations for Z β (namely, spectral Galerkin approximations in Section 5 and sinc-Galerkin approximations in Section 6) defined with respect to family (V N ) N ∈N of subspaces V N ⊂ H 1 0 (D) of finite dimension dim(V N ) = N < ∞, we prove convergence at (essentially) optimal rates.More precisely, under minimal regularity conditions on the operator L in (1.2) and for 0 ≤ σ < 2β − d /2, δ ∈ (0, σ), within a suitable parameter range we show that for all ε, q > 0 there exists a constant C > 0 such that, for all N ∈ N, Here, β , β N : D × D → R denote the covariance functions of the Whittle-Matérn field Z β and of the Galerkin approximation Z β N , respectively.For details, see Corollaries 5.1-5.3 for spectral Galerkin approximations, and Theorems 6.18, 6.23 for the sinc-Galerkin approach."Suitable parameter range" refers to the observations that (i) if a finite element method of polynomial degree p ∈ N is used to define the sinc-Galerkin approximation or (ii) if L in (1.2) is H 1+α (D)-regular for 0 < α ≤ 1 maximal (see Definition 6.20), then the convergence rates of the sinc-Galerkin approximation cannot exceed p + 1 − σ or min{1 + α − σ, 2α}, where 0 ≤ σ ≤ 1.
We point out that due to the low regularity of white noise, dW ∈ H − d /2−ε (D), which holds P-almost surely and in L q (Ω) (cf.[2,Prop. 2.3]) the convergence results (1.3)- (1.6) are (essentially, up to ε > 0) optimal and they are also reflected in our numerical experiments, see Section 7 and the discussion in Section 8.Note furthermore that the convergence rates in (1.4), (1.6) of the field with respect to L q (Ω; C δ (D)) and of the covariance function in the C(D × D)-norm, which we obtain via a Kolmogorov-Chentsov argument, are by d /2 better than applying the results (1.3), (1.5) combined with the Sobolev embeddings H δ+ d /2 (D) → C δ (D) and H ε+ d /2, ε+ d /2 (D × D) → C(D × D), respectively.We remark that strong convergence of the sinc-Galerkin approximation with respect to the L 2 (Ω; L 2 (D))-norm, i.e., (1.3) for σ = 0, at the rate 2β − d /2 has already been proven in [2,Thm. 2.10].However, the assumptions made in [2, Ass.2.6 and Eq.(2. 19)] require the differential operator L to be at least H 2 (D)-regular.Thus, our results do not only generalize the analysis of [2] for the strong error to different norms, but also to less regular differential operators.This is of relevance for several practical applications, since the spatial domain, where the GRF is simulated, may be non-convex or the coefficient A may have jumps.For this reason, in Subsection 6.3.2 we work under the assumption that L is H 1+α (D)-regular for some 0 < α ≤ 1 (for instance, α < π /ω if D is a non-convex domain with largest interior angle ω > π).
As an interim result while deriving the error bounds (1.3)-(1.6)for the sinc-Galerkin approximation, we prove a non-trivial extension of one of the main results in [5].Namely, we show that for all β > 0, 0 ≤ σ ≤ min{1, 2β}, −1 ≤ δ ≤ 1 + α, δ = 1 /2 with 2β + δ − σ > 0, and for all ε > 0, there exists a constant C > 0 such that, for N ∈ N and g ∈ H δ (D), Here, L −1 N : H −1 (D) → V N denotes the approximation of the data-to-solution map L −1 : H −1 (D) → H 1 0 (D) with respect to the Galerkin space V N ⊂ H 1 0 (D).For details see Theorem 6.6, Remark 6.7 and Lemmata 6. 21-6.22.This error estimate was proven in [5,Thm. 4.3 & Rem. 4.1] only for β ∈ (0, 1), σ = 0, and δ ≥ 0, see also the comparison in Remark 6.8.1.2.Outline.After specifying the mathematical setting as well as our notation in Subsections 1.3-1.4,we rigorously define the second-order elliptic differential operator L from (1.2) under minimal assumptions on the coefficients A, κ and the domain D ⊂ R d in Section 2; thereby collecting several auxiliary results for this type of operators.Section 3 is devoted to the regularity analysis of a GRF colored by a linear operator T which is bounded on L 2 (D).These results are subsequently applied in Section 4 to the class of generalized Whittle-Matérn fields, where T := L −β with L defined as in Section 2 and β > d /4.In Section 5 we derive the convergence results (1.3)- (1.6) for spectral Galerkin approximations where the finite-dimensional subspace V N is generated by the eigenvectors of the operator L corresponding to the N smallest eigenvalues.We then investigate sinc-Galerkin approximations in Section 6, where we first let V N be an abstract Galerkin space satisfying certain approximation properties, see Subsections 6.1-6.2.Subsequently, in Subsection 6. 3 we show that these properties are indeed satisfied if the Galerkin spaces originate from a quasi-uniform family of finite element discretizations of polynomial degree p ∈ N, and we discuss the convergence behavior for two cases in detail: (i) the coefficients A, κ and the domain D in (1.2) are smooth, and (ii) A, κ, D are such that the differential operator L in (1.2) is only H 1+α (D)-regular for some 0 < α ≤ 1.In Section 7 we perform several numerical experiments for the model example (1.1), d = 1, and sinc-Galerkin discretizations generated with a finite element method of polynomial degree p ∈ {1, 2}.In Section 8 we reflect on our outcomes.If (E, • E ) is a Banach space, then (E * , • E * ) denotes its dual, • , • E * ×E the duality pairing on E * × E, Id E the identity on E, and L(E; F ) the space of bounded linear operators from (E, • E ) to another Banach space (F, If not specified otherwise, ( • , • ) H is the inner product on a Hilbert space H and L 2 (H; U ) ⊆ L(H; U ) denotes the Hilbert space of Hilbert-Schmidt operators between two Hilbert spaces H and U .The adjoint of T ∈ L(H; U ) is identified with T * ∈ L(U ; H) (via the Riesz maps on H and on U ).We write L(E) and L 2 (H) whenever E = F and H = U .The domain of a possibly unbounded operator L is denoted by D(L).
For 1 ≤ q < ∞, L q (D; E) is the space of (equivalence classes of) E-valued, Bochner measurable, q-integrable functions on D and L q (Ω; E) denotes the space of (equivalence classes of) E-valued random variables with finite q-th moment, i.e., The space L ∞ (D; E) consists of all equivalence classes of E-valued, Bochner measurable functions which are essentially bounded on D, i.e., For γ ∈ (0, 1), we furthermore define the mappings on the Banach space Note that the norm • C γ (D;E) renders the subspace of γ-Hölder continuous functions a Banach space.Whenever the functions or random variables are real-valued, we omit the image space and write C(D), C γ (D), L q (D), and L q (Ω), respectively.For σ > 0, the (integer-or fractional-order) Sobolev space is denoted by H σ (D) (see [12,Sec. 2], see also [41,Sec. 1.11.4/5]), and We mark equations which hold almost everywhere or P-almost surely with a.e. and P-a.s., respectively.For two random variables X, Y , we write X d = Y whenever X and Y have the same probability distribution.The Dirac measure at x ∈ D is denoted by ð x .Given a parameter set P and mappings A, B : P → R, we let A(p) B(p) denote the relation that there exists a constant C > 0, independent of p ∈ P, such that A(p) ≤ CB(p) for all p ∈ P. For a further parameter set Q and mappings A, B : P × Q → R, we write A(p, q) q B(p, q) if, for all q ∈ Q, there exists a constant C q > 0, independent of p ∈ P, such that A(p, q) ≤ C q B(p, q) for all p ∈ P and q ∈ Q.Finally, A(p) B(p) indicates that both relations, A(p) B(p) and B(p) A(p), hold simultaneously; and similarly for A(p, q) q B(p, q).

Auxiliary results on second-order elliptic differential operators
As outlined in Subsection 1.1, the overall objective of this article is to study (generalized) Whittle-Matérn fields and Galerkin approximations for them.Here, we call a Gaussian random field a generalized Whittle-Matérn field if its covariance operator is given by a negative fractional power of a second-order elliptic differential operator.The purpose of this section is to present preliminary results on secondorder differential operators which will be of importance for the regularity and error analysis of these fields.
Firstly, we specify the class of differential operators that we consider.We start by formulating assumptions on the coefficients of the operator.).By this we mean that D(L) consists of precisely those u ∈ H 1 0 (D) for which there exists a constant C ≥ 0 such that and, for u ∈ D(L), Lu is the unique element of L 2 (D) which, for all v ∈ H 1 0 (D), satisfies (2. 2) It is well-known that the operator L : 2) is densely defined and self-adjoint (e.g., [31,Prop. 1.22 and Prop. 1.24]).Furthermore, by the Lax-Milgram lemma, its inverse exists and extends to a bounded linear operator [31,Lem. 1.3]).By the Kondrachov compactness theorem L −1 : L 2 (D) → L 2 (D) is compact (e.g., [18,Thm. 7.22]).
For this reason, the spectrum of L consists of a system of only positive eigenvalues (λ j ) j∈N with no accumulation point, whence we can assume them to be in nondecreasing order.The following asymptotic spectral behavior, known as Weyl's law (see, e.g., [11,Thm. 6.3.1]), will be exploited several times in our analysis.Lemma 2.2.Let L be the second-order differential operator in (2.2), defined with respect to the bounded open domain D ⊂ R d , and with coefficients A and κ fulfilling Assumptions 2.1.I-II.Then, the eigenvalues of L (in nondecreasing order) satisfy We let E := {e j } j∈N denote a system of eigenvectors of the operator L in (2.2) which corresponds to the eigenvalues (λ j ) j∈N and which is orthonormal in L 2 (D).Note that, for σ > 0, the fractional power operator is well-defined.Indeed, on the domain the action of L σ is given via the spectral representation is itself a Hilbert space with respect to the inner product and the corresponding induced norm • σ .In what follows, we let Ḣ0 L := L 2 (D) and, for σ > 0, Ḣ−σ L denotes the dual space ( Ḣσ L ) * after identification via the inner product ( • , • ) L2(D) on L 2 (D) which is continuously extended to a duality pairing.
In order to derive regularity and convergence results with respect to the Sobolev space H σ (D) and the space C γ (D) of γ-Hölder continuous functions in (1.9), we wish to relate the various norms involved by well-known results from interpolation theory and by Sobolev embeddings.To this end, we need to consider various assumptions on the spatial domain D, specified below. ) and the norms . We now establish the reverse embedding.By Assumption 2.1.III and, e.g., [16,Thm. 4

Thus (by first approximating
3. General results on Gaussian Random Fields (GRFs) In this section we address different notions of regularity (Hölder and Sobolev) for Gaussian random fields (GRFs) and their covariance functions.We first recall the definition of a GRF and specify then what we mean by a colored GRF.As usually, we work in the setting formulated in Subsection 1.3.= tr(T T * ), where tr( • ) is the trace on L 2 (D).3.1.Hölder regularity of GRFs.We now provide an abstract result on the construction and Hölder regularity of a GRF assuming that the color and, thus, the covariance structure of the field is given.Proposition 3.4.Assume that T ∈ L(L 2 (D); C γ (D)) for some γ ∈ (0, 1).Then there exists a continuous GRF Z colored by T such that Z(x) = W(T * ð x ) P-a.s.∀x ∈ D. (3.2) Furthermore, for q ∈ (0, ∞) and δ ∈ (0, γ), we have Proof.We first define the random field Z 0 : D × Ω → R by Z 0 (x) := W(T * ð x ) for all x ∈ D. By the properties of an isonormal Gaussian process we find, for x, y ∈ D, ) is a real-valued Gaussian random variable, we can apply the Khintchine inequalities (see, e.g., [25,Thm. 4.7 and p. 103]) and conclude with (3.4) that, for all q ∈ (0, ∞), the estimate holds, with a constant C q > 0 depending only on q.
We close the subsection with a brief discussion on (i) the continuity of covariance functions of colored GRFs, and (ii) the L ∞ (D ×D)-distance between two covariance functions of GRFs colored by different operators.
By definition, the covariance function (3.10) We obtain the one-to-one correspondence From this definition it is evident that a GRF Z colored by T (note that E[Z] = 0 by construction, see Definition 3.2) has the covariance operator C = T T * .In the next lemma, this relation is exploited to characterize continuity of the covariance function in terms of the color T of the GRF Z.
Proof.By (3.10), the covariance function of a GRF Z colored by T is given by . Finally, the estimate (3.12) can be shown similarly since, for all x, y ∈ D,

Sobolev regularity of GRFs and their covariances. After having characterized
(i) the Hölder regularity (in L q (Ω)-sense) of a GRF Z, and (ii) continuity of the covariance function in (3.10), in terms of the color of Z, we now proceed with this discussion for Sobolev spaces.Specifically, we investigate the regularity of Z in L q (Ω; H σ (D)) and of the covariance function with respect to the norm on the mixed Sobolev space Here, ⊗ denotes the tensor product of Hilbert spaces.Thus, the inner product on To this end, in the following proposition we first quantify the Ḣσ L -regularity (in L q (Ω)-sense) of a colored GRF in terms of its color, cf.(2.4) and Definition 3.2.In addition, we specify the regularity of the covariance function (3.10) in the Hilbert tensor product space Ḣσ,σ cf. (3.14).Finally, we characterize the distance between two GRFs which are colored by different operators with respect to these norms.

1).
Then Z is square-integrable, i.e., Z ∈ L 2 (D × Ω), if and only if its covariance operator C = T T * has a finite trace on L 2 (D).More generally, for σ ≥ 0 and q ∈ (0, ∞), we have ) Here, tr( • ) is the trace on L 2 (D), L is the differential operator in (2.2) with coefficients A, κ satisfying Assumptions 2.1.I-II, is the covariance function of Z, see (3.10), and , with covariance function and covariance operator C = T T * , we have, for σ ≥ 0 and q ∈ (0, ∞), Proof.Assume first that Z ∈ L 2 (D × Ω).Since Z has mean zero and since it is colored by T ∈ L(L 2 (D)), we obtain C = T T * , i.e., By choosing φ = ψ := λ σ /2 j e j , summing these equalities over j ∈ N, and exchanging the order of summation and expectation, we obtain the identity and the first part of the proposition as well as (3.16) Then we obtain (3.19) from (3.17), since Z − Z is again a GRF, colored by T − T , see (3.1) and Definition 3.2.Furthermore, we find This proves (3.20) and (3.18) follows from this result for Z ≡ 0.

Regularity of Whittle-Matérn fields
In this section we focus on the regularity of (generalized) Whittle-Matérn fields, i.e., of GRFs colored (cf.Definition 3.2) by a negative fractional power of the differential operator L as provided in (2.2).Specifically, we consider for We emphasize the dependence of the covariance structure of Z β on the fractional exponent β > 0 by the index and write β for the covariance function (3.10) of Z β .The first aim of this section is to apply Proposition 3.7 for specifying the regularity of Z β in (4.1) and of its covariance function β with respect to the spaces Ḣσ L and Ḣσ,σ L in (2.4), (3.15).As already pointed out in Remark 3.8, provided that the assumptions of Lemma 2.4 are satisfied, this implies regularity in the Sobolev space H σ (D) and in the mixed Sobolev space H σ,σ (D × D) in (3.14), respectively.
Besides this regularity result with respect to the spaces Ḣσ L and H σ (D), we obtain a stability estimate with respect to the Hölder norm from Corollary 3.5 and continuity of the covariance function from Proposition 3.6.Although we believe that, at least in some specific cases, these results are well-known, for the sake of completeness, we derive them here in our general framework.Lemma 4.1.Let Assumptions 2.1.I-II be fulfilled, β, q ∈ (0, ∞), σ ≥ 0, and Z β be the Whittle-Matérn field in (4.1), with covariance function β .Then, If, in addition, Assumption 2.3.I and 0 ≤ σ ≤ 1 (or Assumptions 2.1.I-III, 2.3.II, and 0 ≤ σ ≤ 2) hold, then the assertions (i)-(ii) remain true if we formulate them with respect to the Sobolev norms Proof.By Proposition 3.7 we have, for any β, q ∈ (0, ∞) and σ ≥ 0, for all x ∈ D, and, for every δ ∈ (0, γ) and q ∈ (0, ∞), the bound for the q-th moment of Z β with respect to the δ-Hölder norm, cf.(1.8), holds.
Proof.Note that by definition of Ḣσ L , see (2.4), for any β > 0, the operator The proof is then completed by applying Corollary 3.5 in both cases (i)/(ii).Proof.By Proposition 3.6(i) we have to show boundedness of is finite, provided that β > d /4.Here, we have used the Cauchy-Schwarz inequality and the spectral behavior (2.3) from Lemma 2.2 in the last estimate.Similarly,  L is the corresponding reproducing kernel Hilbert space, cf.[37].

Spectral Galerkin approximations
In this section we investigate convergence of spectral Galerkin approximations for the Whittle-Matérn field Z β in (4.1).Recall that the covariance structure of the GRF Z β is uniquely determined via its color (3.1) given by the negative fractional power L −β of the second-order differential operator L in (2.2) which is defined with respect to the bounded spatial domain D ⊂ R d .
For N ∈ N, the spectral Galerkin approximation Z β N of Z β is (P-a.s.) defined by i.e., it is a GRF colored by the finite-rank operator mapping to the finite-dimensional subspace V N := span{e 1 , . . ., e N } generated by the first N eigenvectors of L corresponding to the eigenvalues 0 < λ 1 ≤ . . .≤ λ N .The following three corollaries, which provide explicit convergence rates of these approximations and their covariance functions with respect to the truncation parameter N , are consequences of the Propositions 3.4, 3.6 and 3.7.We first formulate the results in the Sobolev norms.
Proof.By Lemma 4.2 there exist continuous random fields and we obtain the convergence result in (5.6) from the stability estimate (4.5) of Lemma 4.2 applied to . Here, we have used the spectral behavior (2.3) from Lemma 2.2 for λ N .

General Galerkin approximations
After having derived error estimates for spectral Galerkin approximations in the previous subsection, we now consider a family of general Galerkin approximations for the Whittle-Matérn field Z β in (4.1) which, for the case β ∈ (0, 1), has been proposed in [2,3].Recall that the random field Z β is indexed by the bounded spatial domain D ⊂ R d .6.1.Sinc-Galerkin approximations.The approximations proposed in [2,3] are based on a Galerkin method for the spatial discretization L h of L and a sinc quadrature for an integral representation of the resulting discrete fractional inverse L −β h .We recall that approach in this subsection, and formulate all assumptions and auxiliary results which are needed for the subsequent error analysis in Subsection 6.2.
6.1.1.Galerkin discretization.We assume that we are given a family (V h ) h>0 of subspaces of H 1 0 (D), with dimension We arrange the eigenvalues of L h in nondecreasing order, and let {e j,h } N h j=1 be a set of corresponding eigenvectors which are orthonormal in L 2 (D).The operator R h : All further assumptions on the finite-dimensional subspaces (V h ) h>0 are summarized below and explicitly referred to, when needed in our error analysis.

Assumption 6.1 (on the Galerkin discretization).
I.There exist θ 1 > θ 0 > 0 and a linear operator I h : H θ1 (D) → V h such that, for all θ 0 < θ ≤ θ 1 , I h : H θ (D) → V h is a continuous extension, and holds for 0 ≤ σ ≤ min{1, θ} and sufficiently small h > 0. II.For all h > 0 sufficiently small and all 0 ≤ σ ≤ 1 the following inverse inequality holds: IV.There exist r, s 0 , t, C 0 , C λ > 0 such that for all h > 0 sufficiently small and for all j ∈ {1, . . ., N h } the following error estimates hold: where {(λ j , e j )} j∈N are the eigenpairs of the operator L in (2.2).
We refer to Subsection 6.3 for explicit examples of finite element spaces (V h ) h>0 , which satisfy these assumptions.Remark 6.2.It is a consequence of the min-max principle that the first inequality in (6.5), λ j ≤ λ j,h , is satisfied for all conforming Galerkin spaces V h ⊂ Ḣ1 L .
In Theorem 6.6 below, we bound the deterministic Galerkin error in the fractional case, i.e., we consider the error between L −β g and L −β h Π h g.This theorem is one of our main results and it will be a crucial ingredient when analyzing general Galerkin approximations of the Whittle-Matérn field Z β from (4.1) in Subsection 6.2.For its derivation, we need the following two lemmata.Lemma 6.3.Let L be as in (2.2) and, for h > 0, let L h , R h be as in (6.1) and (6.2).Suppose Assumptions 2.1.I-II and 2.3.I.Let 0 < α ≤ 1 be such that where Ḣ1+δ L is defined as in in (2.4).Furthermore, let Assumption 6.1.I be satisfied with parameters θ 0 ∈ (0, 1) and θ 1 ≥ 1 + α.Then, for every 0 for sufficiently small h > 0.
Remark 6.9 (p-FEM).Due to the term 2α and 0 < α ≤ 1, (6.13) will be sharp for finite elements of first order, but not for finite elements of polynomial degree p ≥ 2 when β > 1 and the problem is "smooth" such that (6.7) holds with α > 1.
For deriving a bound for (B), we first note that by (6.12) of Lemma 6.5 Next, we define the contour where ω ∈ (0, π) and r := λ1 /2.By, e.g., [32, Ch. 2.6, Eq. ( 6. 3)] we have From the limit ω → π, we then obtain the representation We exploit this integral representation as well as the identity dθ , ( where µ := (1+η+σ) /2, ν := (1+ϑ−δ) /2 and 0 ≤ η ≤ ϑ ≤ α are chosen as follows By (6.12) and (6.9), we find for the term outside of the integral, for h > 0 sufficiently small, where these three cases can be summarized as in (6.13), since 2β It remains to show that the two integrals in (6.14) converge, uniformly in h.To this end, we first note that 0 ≤ µ ≤ 1 and, thus, for any t > 0, By the same argument we find that L ν h (tI + L h ) −1 Π h L(L2(D)) ≤ t ν−1 , for t > 0, since also 0 ≤ ν ≤ 1.Thus, we can bound the first integral arising in (6.14) by Here, we have used that r To estimate the second integral in (6.14), we note that, for any z ∈ C with |z| = λ1 /2, With these observations, we finally can bound the second integral in (6.14), which completes the proof of (6.13) for the case that 0 ≤ δ ≤ 1 + α.Assume now that δ = − σ for some 0 < σ ≤ 1.Then, for g ∈ Ḣδ L , we may exploit (6.13), which has already been proven for 0 ≤ δ ≤ 1 + α, as follows, 0 by assumption.Furthermore, by (6.10) of Lemma 6.5 we have L , for the whole range of parameters σ, δ as stated in the theorem.6.1.2.Sinc quadrature and fully discrete scheme.After the Galerkin discretization (in space), we need a second component to approximate the generalized Whittle-Matérn field Z β in (4.1).Namely, we have to numerically realize a fractional inverse of the Galerkin operator L h in (6.1).To this end, as proposed in [2], we introduce, for β ∈ (0, 1) and k > 0, the sinc quadrature approximation of L −β h from [5], where We also formally define this operator for the case β = 0 by setting Q 0 h,k := Id V h .For a general β = n β + β > 0 as in (4.2), we then consider the approximations ) which are (P-a.s.) defined by Here, the finite-rank operator Π h is given by For β ∈ (0, 1), the construction (6.17) of Z β h,k gives the same approximation as considered in [2,3].Note furthermore that, in contrast to Π h , the operator Π h in (6.18) is neither a projection nor self-adjoint, and its definition depends on the particular choice of the eigenbases {e j } j∈N ⊂ L 2 (D) and {e j,h } N h j=1 ⊂ V h .The reason for why we consider both approximations Z β h,k , Z β h,k will become apparent in the error analysis of Subsection 6.2.Although, in general, they do not coincide in L q (Ω; L 2 (D))-sense, i.e., D) = 0, they have the same Gaussian distribution as shown in the following lemma.Lemma 6.10.Suppose Assumptions 2.1.I-II.Let Π h denote the L 2 (D)-orthogonal projection onto V h , and Π h be the operator in (6.18).Then, if holds for all φ, ψ ∈ L 2 (D).In particular, Z β h,k d = Z β h,k as L 2 (D)-valued random variables, where Z β h,k and Z β h,k are as defined in (6.16)-(6.17).
Proof.Note that (T h Π h ) * (resp.(T h Π h ) * ) denotes the adjoint of T h Π h (resp. of T h Π h ) when interpreted as an operator in L(L 2 (D)).This means, we are identifying (6.19).
By definition of Z β h,k , Z β h,k in (6.16)-(6.17),for M ∈ N and ψ 1 , . . ., ψ M ∈ L 2 (D), the random vectors z, z with entries z j = Z β h,k , ψ j L2(D) and z j = Z β h,k , ψ j L2(D) , 1 ≤ j ≤ M , are multivariate Gaussian distributed.Furthermore, both vanish in expectation and their covariance matrices, C := Cov(z) and C := Cov z , coincide due to (6.19) . This shows that Z β h,k d = Z β h,k as L 2 (D)valued random variables.Remark 6.11 (Simulation in practice).To simulate samples of the in (6.16)-(6.17)abstractly defined (P-a.s.) V h -valued Gaussian random variables Z β h,k and Z β h,k in practice, in both cases, one first has to generate a sample of a multivariate Gaussian random vector b with mean 0 and covariance matrix M, where M is the Gramian with respect to any fixed basis Φ h = {φ j,h } N h j=1 of V h , i.e., M ij := (φ i,h , φ j,h ) L2(D) .This follows from the identical distribution of the GRFs Z 0 h and Z 0 h colored by Π h and Π h , respectively, as well as from the chain of equalities which we obtain from Lemma 6.10 with the random vector Z β k , given by is then the vector of coefficients when expressing the V h -valued sample of Z β h,k (or of Z β h,k ) with respect to the basis Φ h .Here, L ∈ R N h ×N h represents the action of the Galerkin operator L h in (6.1), i.e., L ij := (L h φ j,h , φ i,h ) L2(D) , and, for β ∈ (0, 1), Q β k ∈ R N h ×N h is the matrix analog of the operator Q β h,k from (6.15), i.e., respectively.In order to perform the error analysis for Z β h,k and Z β h,k , we split these operators as follows 2)) which can be estimated with the results from Section 5 on spectral Galerkin approximations.Furthermore, we shall refer to as the Galerkin errors and as the quadrature errors, respectively.
In the following we provide error estimates for both approximations, Z β h,k and Z β h,k in (6.16)-(6.17),with respect to the norm on L q (Ω; H σ (D)) as well as for its covariance functions β h,k , β h,k in the mixed Sobolev norm, cf.(3.14).By exploiting Theorem 6.6 the bounds for Z β h,k and β h,k in Proposition 6.12 below will be sharp if a conforming finite element method with piecewise linear basis functions is used.However, to derive optimal rates for the case of finite elements of higher polynomial degree, a different approach will be necessary, cf.Remark 6.9.To this end, we perform an error analysis for Z β h,k and β h,k based on spectral expansions, see Proposition 6.13.Since these arguments work only if the differential operator L in (2.2) is at least H 2 (D)-regular, both approaches and results are needed for a complete discussion of smooth vs. H 1+α (D)-regular problems in Subsection 6.3.Finally, in Proposition 6.14, we use the approximation Z β h,k from (6.16) to formulate a convergence result with respect to the Hölder norm (1.8) in L q (Ω)-sense and with respect to the L ∞ (D × D)-norm for its covariance function β h,k .We note that, at the cost of other assumptions (e.g., α > 1 /2) on the parameters involved, it is possible to circumvent the additional condition β > 1 (instead of β > 3 /4) needed in the following proposition for the L q (Ω; H σ (D))-estimate if d = 3. Proposition 6.12.Suppose Assumptions 2.1.I-II, 2.3.I, 6.1.II-III, and let Assumption 6.1.I be satisfied with parameters θ 0 ∈ (0, 1) and θ 1 ≥ 1 + α, where 0 < α ≤ 1 is as in (6.7).Assume furthermore that Π h is H 1 (D)-stable, see (6.11), and that d ∈ {1, 2, 3}, β > 0 and 0 ≤ σ ≤ 1 are such that 2β − σ > d /2.Let Z β be the Whittle-Matérn field in (4.1) and, for h, k > 0, let Z β h,k be the sinc-Galerkin approximation in (6.16), with covariance functions β and β h,k , respectively.Then, for every q, ε > 0 and sufficiently small h > 0, where, if d = 3, for (6.24) to hold, we also suppose that β > 1 and α ≥ 1 /2 − σ.
Here Z β h denotes a GRF colored by L −β h Π h , with covariance function β h .Furthermore, we note the following: For m ≥ 0, we have where the observation of Remark 6.2 was used in the last step.Thus, by the spectral asymptotics from Lemma 2.2 and by Assumption 6.1.III we have for m ≥ 0, m For terms (A Z ) and (B Z ), we obtain with the definitions of the Galerkin and quadrature errors For bounding term (A Z ), we let γ ∈ (0, β) and rewrite E β V h from (6.22) as follows, We first bound (A Z ) for d ∈ {1, 2}.To this end, let ε 0 > 0 be chosen sufficiently small such that 2β − σ − d /2 > 4ε 0 and choose γ := d /4 + ε 0 in (6.28).We obtain thus (A Z ) q (A Z ) + (A Z ), where For (A Z ), we find by (6.13) of Theorem 6.6 and by (6.27), applied for the parameters β := β − d /4 − ε 0 , σ := σ, δ := 0, and m = d /4 + ε 0 , respectively, , for any ε > 0 and sufficiently small h > 0.
Finally, we use the estimate as well as the inverse inequality (6.4) to conclude for term (B ) for σ ∈ {0, 1} that Combining the above estimate with (6.30) and stability of the operators which is uniform in h and k for sufficiently small h, k > 0, shows that Interpolation for σ ∈ (0, 1) completes the proof of (6.25).
For the L ∞ (D×D)-estimate (6.38) of the covariance function, fix ε ∈ (0, 2).First, we recall the Sobolev embedding H ε /4+ 1 /2 (D) → C ε /4 (D) as well as the equivalence of the spaces , see Lemma 2.4.We then conclude with (3.12) of Proposition 3.6(ii) that, for σ : . By (6.13) of Theorem 6.6 we have Furthermore, we find, similarly as in (6.31), that , where we have used the inverse inequality (6.4) in the last step.Finally, since , the proof is completed by (6.39) combined with the uniform stability (6.32) of L −β h and Q β h,k .6.3.Application to finite element approximations.We now discuss different scenarios of (i) regularity of the second-order differential L in (2.2), (ii) finite element (FE) discretizations satisfying Assumptions 6.1.I-IV for specific values of 0 < θ 0 < θ 1 and of r, s 0 , t > 0. We then obtain explicit rates of convergence for the FE Galerkin approximations Z β h,k , Z β h,k in (6.16)-(6.17)from Propositions 6.12, 6.13 and 6.14.Assumption 6.15 (FE discretization).Throughout this subsection, we suppose the following setting: • the (minimal) Assumptions 2.1.I-II on the coefficients A, κ of the operator L; • Assumptions 2.3.I, i.e., D ⊂ R d is a bounded Lipschitz domain; • (T h ) h>0 is a quasi-uniform family of triangulations on D, indexed by the mesh width h > 0; • the basis functions of the finite-dimensional space V h ⊂ H 1 0 (D) are continuous on D and piecewise polynomial with respect to T h of degree at most p ∈ N.
All further assumptions on the operator L, on the domain D, and on the FE spaces are explicitly specified for each case.Note that quasi-uniformity of (T h ) h>0 already guarantees that Assumptions 6.1.II and 6.1.III are satisfied (6.1.III is obvious, for the inverse inequality 6.1.II see, e.g., [15,Cor. 1.141]).
In Subsection 6.3.1 we briefly comment on the situation of smooth coefficients and apply Proposition 6.13 to derive optimal convergence rates when p ≥ 1. Afterwards, in Subsection 6.3.2 we focus on less regular problems and p = 1 by using the results from Propositions 6.12 and 6.14.6.3.1.The smooth case.The remaining crucial ingredient in order to derive explicit rates of convergence from Proposition 6.13 is to prove validity of Assumption 6.1.IV for the finite element spaces (V h ) h>0 .For the case of a second-order elliptic differential operator L with smooth coefficients, these results are well-known and we summarize them below.Assumption 6.16 (smooth case).The domain D has a smooth C ∞ -boundary ∂D, and the coefficients of L in (2.2) are smooth, i.e., A ∈ C ∞ (D) d×d and κ ∈ C ∞ (D).Furthermore, the Rayleigh-Ritz projection R h : Lemma 6.17.Suppose Assumptions 6.15 and 6.16.Then, Assumption 6.1.IV is satisfied for r = 2p and s 0 = t = p + 1.
Proof.See, e.g., [38,Thm. 6.1 & Thm.6.2].Theorem 6.18.Suppose Assumptions 6.15 and 6.16.Let d ∈ N, β > 0, and 0 ≤ σ ≤ 1 be such that 2β − σ > d /2, let Z β be the Whittle-Matérn field in (4.1) and, for h, k > 0, let Z β h,k be the sinc-Galerkin approximation in (6.17), and let β , β h,k denote their covariance functions.Then we have, for sufficiently small h > 0, sufficiently small k = k(h) > 0, and all q > 0, where C β,h , C Z β,h and P are as in Proposition 6.13.Proof.By Lemma 6.17 we have r = 2p, s 0 = t = p + 1 and, thus, for γ ∈ {0, 1}, In particular, when the integer part does not vanish, n β ∈ N, a polynomial degree p > 1 is meaningful, since thus higher order convergence rates can be achieved, cf. the numerical experiments in Section 7. 6.3.2.Less regularity.We now discuss convergence of FE discretizations when the operator L in (2.2) has a coefficient A which is not necessarily Lipschitz continuous or the domain D is not convex, i.e., the general case that L is only H 1+α (D)-regular.In the following definition we specify what we mean by this.Definition 6.20.Suppose Assumptions 2.1.I-II, 2.3.I, let 0 < α ≤ 1 and L be the second-order differential operator in (2.2).We say that the elliptic problem associated with L is H 1+α (D)-regular if the restriction of L :
The reference solutions for the field and the covariance function are generated based on an overkill Karhunen-Loève expansion of Z β with N KL = 1000 terms, For d = 1, the operator L does not have multiple eigenvalues and we can assemble the matrix R, for each h ∈ {h 0 , . . ., h 4 }, by computing the discrete eigenfunctions {e j,h } N h j=1 and by adjusting their sign so that e j,h indeed approximates e j for each j ∈ {1, . . ., N h }.Note that we only have to assemble this matrix R to have comparable samples of the sinc-Galerkin approximation and the reference solution needed for the strong error studies.For the simulation practice, one could compute the Cholesky factor of the Gramian M or approximate the matrix square root √ M, e.g., as proposed in [23], in order to sample from b. Since furthermore the dimension of the finite element spaces, even at the highest level = 4, is relatively small, we can assemble the covariance matrices of the sinc-Galerkin approximation directly, without Monte Carlo sampling, as cf. (6.20)-(6.21).
For every of the 100 Monte Carlo samples, we approximate the integrals needed for computing the L 2 (D) and H 1 0 (D)-errors by using MATLAB's built-in function integral with tolerance 1e-6.For the L ∞ -studies we consider the largest error with respect to an equidistant mesh on on D = [0, 1] with N ok = 1001 nodes, i.e., where x j := (j − 1)10 −3 .Furthermore, to compute the L 2 (D × D)-error, we approximate the distance of the covariances by a function which is piecewise constant on a regular lattice with N 2 ok nodes.Finally, the empirical convergence rates, also  2.
shown in Table 2, are obtained via a least-squares affine fit with respect to the data set {(ln h , ln err ) : 2 ≤ ≤ 4}.Here, err denotes the error on level with respect to the norm used in the study and for the respective value of β and p.
The resulting observed errors are displayed in Figure 1 for the fields and in Figure 2 for the covariances.Overall, the empirical results validate our theoretical outcomes fairly well, with a slight deviation for the L ∞ -studies which may be caused by a larger pre-asymptotic range.2.

Conclusion and discussion
We have identified necessary and sufficient conditions for square-integrability, Sobolev regularity, and Hölder continuity (in L q (Ω)-sense) for GRFs in terms of their color, as well as square-integrability, mixed Sobolev regularity, and continuity of their covariance functions, see Propositions 3.4, 3.6 and 3.7.Subsequently, we have applied these findings to generalized Whittle-Matérn fields, see Z β in (4.1),where these conditions become assumptions on the smoothness parameter β > 0, corresponding to the fractional exponent of the color L −β , see Lemmata 4.1-4.3.
While these regularity results readily implied convergence of spectral Galerkin approximations, see Corollaries 5.1-5.3,significantly more work was needed to derive convergence for general Galerkin (such as finite element) approximations, for the following reason: It was unknown, how the deterministic fractional Galerkin error L −β g − L −β h g behaves in the Sobolev space H σ (D), for 0 ≤ σ ≤ 1, all possible exponents β > 0, and sources g ∈ H δ (D) of possibly negative regularity δ < 0. We have identified this behavior in Theorem 6.6 for the general situation that the second-order elliptic differential operator L is H 1+α (D)-regular for some 0 < α ≤ 1.
In conclusion, ) * e ,h = λ Here, we have used the uniform stability of L −β , Q β h,k with respect to h and k, see (6.32), as well as (2.3) from Lemma 2.2 and Assumption 6.1.I.Combining the bounds for (A ), (B ), (C ) completes the proof.

1. 3 .
Setting.Throughout this article, we let (Ω, F, P) be a complete probability space with expectation operator E, and D be a bounded, connected and open subset of R d , d ∈ N, with closure D. In addition, we let W : L 2 (D) → L 2 (Ω) be an L 2 (D)-isonormal Gaussian process in the sense of [29, Def.1.1.1].1.4.Notation.For B ⊆ R d , B(B) denotes the Borel σ-algebra on B (i.e., the σ-algebra generated by the sets that are relatively open in B).For two σ-algebras F and G, F ⊗ G is the σ-algebra generated by F × G.

Proposition 3 . 6 .
Let Z, Z be GRFs colored by T and T , respectively, see (3.1), with covariance functions denoted by and , cf. (3.10).Then, (i) has a continuous representative on D × D (again denoted by ) if and only if T ∈ L(L 2 (D); C(D)).In this case, sup x,y∈D

Figure 2 .
Figure 2. Observed L q (D × D)-error of the covariance function for q ∈ {2, ∞} (top, bottom), polynomial degree p ∈ {1, 2} (left, right), and different values of β, shown in a log-log scale as a function of the mesh width h.The corresponding observed convergence rates are shown in Table2.