Regularity and convergence analysis in Sobolev and H\"older spaces for generalized Whittle-Mat\'ern fields

We analyze several Galerkin approximations of a Gaussian random field $\mathcal{Z}\colon\mathcal{D}\times\Omega\to\mathbb{R}$ indexed by a Euclidean domain $\mathcal{D}\subset\mathbb{R}^d$ whose covariance structure is determined by a negative fractional power $L^{-2\beta}$ of a second-order elliptic differential operator $L:= -\nabla\cdot(A\nabla) + \kappa^2$. Under minimal assumptions on the domain $\mathcal{D}$, the coefficients $A\colon\mathcal{D}\to\mathbb{R}^{d\times d}$, $\kappa\colon\mathcal{D}\to\mathbb{R}$, and the fractional exponent $\beta>0$, we prove convergence in $L_q(\Omega; H^\sigma(\mathcal{D}))$ and in $L_q(\Omega; C^\delta(\overline{\mathcal{D}}))$ at (essentially) optimal rates for (i) spectral Galerkin methods and (ii) finite element approximations. Specifically, our analysis is solely based on $H^{1+\alpha}(\mathcal{D})$-regularity of the differential operator $L$, where $0<\alpha\leq 1$. For this setting, we furthermore provide rigorous estimates for the error in the covariance function of these approximations in $L_{\infty}(\mathcal{D}\times\mathcal{D})$ and in the mixed Sobolev space $H^{\sigma,\sigma}(\mathcal{D}\times\mathcal{D})$, showing convergence which is more than twice as fast compared to the corresponding $L_q(\Omega; H^\sigma(\mathcal{D}))$-rate. For the well-known example of such Gaussian random fields, the original Whittle-Mat\'ern class, where $L=-\Delta + \kappa^2$ and $\kappa \equiv \operatorname{const.}$, we perform several numerical experiments which validate our theoretical results.