Exponential convergence of hp FEM for spectral fractional diffusion in polygons

For the spectral fractional diffusion operator of order 2s, s∈(0,1)\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$s \in (0,1)$$\end{document}, in bounded, curvilinear polygonal domains Ω⊂R2\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\varOmega \subset {{\mathbb {R}}}^2$$\end{document} we prove exponential convergence of two classes of hp discretizations under the assumption of analytic data (coefficients and source terms, without any boundary compatibility), in the natural fractional Sobolev norm Hs(Ω)\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$${\mathbb {H}}^s(\varOmega )$$\end{document}. The first hp discretization is based on writing the solution as a co-normal derivative of a 2+1\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$2+1$$\end{document}-dimensional local, linear elliptic boundary value problem, to which an hp-FE discretization is applied. A diagonalization in the extended variable reduces the numerical approximation of the inverse of the spectral fractional diffusion operator to the numerical approximation of a system of local, decoupled, second order reaction-diffusion equations inΩ\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\varOmega $$\end{document}. Leveraging results on robust exponential convergence of hp-FEM for second order, linear reaction diffusion boundary value problems in Ω\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\varOmega $$\end{document}, exponential convergence rates for solutions u∈Hs(Ω)\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$u\in {\mathbb {H}}^s(\varOmega )$$\end{document} of Lsu=f\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\mathcal {L}^s u = f$$\end{document} follow. Key ingredient in this hp-FEM are boundary fitted meshes with geometric mesh refinement towards∂Ω\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\partial \varOmega $$\end{document}. The second discretization is based on exponentially convergent numerical sinc quadrature approximations of the Balakrishnan integral representation of L-s\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\mathcal {L}^{-s}$$\end{document} combined with hp-FE discretizations of a decoupled system of local, linear, singularly perturbed reaction-diffusion equations inΩ\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\varOmega $$\end{document}. The present analysis for either approach extends to (polygonal subsets M~\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$${\widetilde{\mathcal {M}}}$$\end{document} of) analytic, compact 2-manifolds M\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$${\mathcal {M}}$$\end{document}, parametrized by a global, analytic chart χ\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\chi $$\end{document} with polygonal Euclidean parameter domain Ω⊂R2\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\varOmega \subset {{\mathbb {R}}}^2$$\end{document}. Numerical experiments for model problems in nonconvex polygonal domains and with incompatible data confirm the theoretical results. Exponentially small bounds on Kolmogorov n-widths of solution sets for spectral fractional diffusion in curvilinear polygons and for analytic source terms are deduced.


Introduction
In recent years, the mathematical and numerical analysis of initial-boundary value problems for fractional differential operators has received substantial attention.Their numerical treatment has to overcome several challenges.The first challenge arises from their nonlocal nature as integral operators.A direct Galerkin discretization leads to fully populated system matrices, and compression techniques (see, e.g., [18] and the references there) have to be brought to bear to make the discretization computationally tractable.An alternative to a direct Galerkin discretization of an integral operator, which is possible for the presently considered spectral fractional Laplacian, is to realize the nonlocal operator numerically as a Dirichlet-to-Neumann operator for a local (but degenerate) elliptic problem.While this approach, sometime referred to as "Caffarelli-Silvestre" extension ("CS-extension" for short) [14,40] increases the spatial dimension by 1, it permits to use the mathematical and numerical tools that were developed for local, integer order differential operators.In the present paper, we study several hp-FE discretizations of the resulting local (but degenerate) elliptic problem, to which we will refer as Extended hp-FEM.
One alternative to the extension approach is the representation of fractional powers of elliptic operators as Dunford-Taylor integrals proposed in [9,10].Discretizing such an integral leads to a sum of solution operators for local, second order elliptic problems, which turn out to be singularly perturbed, but are amenable to established numerical techniques.In the present paper, we also study this approach under the name sinc Balakrishnan FEM (sinc BK-FEM for short).
A second challenge arises from the fact that the solutions of problems involving fractional operators are typically not smooth, even for smooth input data (cf.the examples and discussion in [8,Sec. 8.4]).Indeed, for the spectral fractional Laplacian, the behavior near a smooth boundary ∂Ω is u ∼ u 0 + O(dist(•, ∂Ω) β ) for some more regular u 0 and a β ∈ N [15, Thm.1.3].Points of non-smoothness of ∂Ω introduce further singularities into the solution.The numerical resolution of both types of singularities requires suitably designed approximation spaces.For the spectral fractional Laplacian in two-dimensional polygonal domains, we present a class of meshes in Ω with anisotropic, geometric refinement towards ∂Ω and with isotropic geometric refinement towards the corners of Ω.We show that spaces of piecewise polynomials on such meshes can lead to exponential convergence.

Spectral Fractional Diffusion
When dealing with fractional operators, care must be exercised in stating the definition of the fractional powers.Here, we consider the so-called spectral fractional diffusion operators as investigated in [14].We refer to the surveys [11,20,34] and the references there for a comparison of the different definitions of fractional powers of the Dirichlet Laplacian.
We consider the linear, elliptic, self-adjoint, second order differential operator w → Lw = −div(A∇w), in a bounded, curvilinear polygon Ω ⊂ R 2 as described in Section 1.1.The diffusion coefficient A ∈ L ∞ (Ω, GL(R 2 )) is assumed symmetric, uniformly positive definite.The data A and f are assumed analytic in Ω.We quantify analyticity of A and f by assuming that there are Here, the notation |D p A| signifies |α|=p |D α A|, with the usual multi-index convention D α denoting mixed weak derivatives of order α ∈ N 2 0 whose total order |α| = α 1 + α 2 .Further, we employ standard notation for (fractional) Sobolev spaces H t (Ω), consistent with the notation and definitions in [22].
We introduce the "energy" inner product a Ω (•, •) on H 1 0 (Ω) associated with the differential operator L by a Ω (w, v) = ˆΩ (A∇w • ∇v) dx . (1. 2) The operator L : H 1 0 (Ω) → H −1 (Ω) induced by this bilinear form is an isomorphism, due to the (assumed) positive definiteness of A. Let {λ k , ϕ k } k∈N ⊂ R + × H 1 0 (Ω) be a sequence of eigenpairs of L, normalized such that {ϕ k } k∈N is an orthonormal basis of L 2 (Ω) and an orthogonal basis of (H 1 0 (Ω), a Ω (•, •)).We introduce, for σ ≥ 0, the domains of fractional powers of L as We denote by H −σ (Ω) the dual space of H σ (Ω).Denoting by •, • the H −σ (Ω)× H σ (Ω) duality pairing that extends the standard L 2 (Ω) inner product, we can identify elements f ∈ H −σ (Ω) with sequences {f k } k (written formally as With this identification, we can extend the definition of the norm in (1.3) to σ < 0. Furthermore, the linear operator k ϕ k is bounded and the Dirichlet problem for the fractional diffusion in Ω may be stated as: given a fractional order s ∈ (0, 1] and f ∈ H −s (Ω), find u ∈ H s (Ω) such that (1.4) The ellipticity estimate w, L s w ≥ λ s 1 w 2 H s (Ω) valid for every w ∈ H s (Ω) implies the unique solvability of (1.4) for every f ∈ H −s (Ω).The hp-FEM approximations of (1.4) developed and analyzed in the present work are not based on explicit or approximated eigenfunctions but instead on the localization of the operator L s in terms of extension discussed in Sec.1.4 and on the Dunford-Taylor integral discussed in Sec.1.5, the so-called Balakrishnan formula.

Contributions
We briefly highlight the principal contributions of this work.For the nonlocal, spectral fractional diffusion problem (1.4) in bounded, curvilinear polygonal domains Ω as described in Section 1.1 and with analytic data A and f as in (1.1), and without any boundary compatibility, we develop two hp-FEMs for (1.4) that converge exponentially in terms of the number of degrees of freedom N DOF in Ω.The setting covers in particular also boundary value problems for fractional surface diffusion on analytic surface pieces as in the setting of Section 8.1.Key insight in our error analysis is that either method, based on the extension of (1.4) combined with a diagonalization procedure as in [8,Sec. 6] or on a contour-integral representation of L s combined with an exponentially converging sinc quadrature, reduce the numerical solution of (1.4) to the numerical solution of local, singularly perturbed second order reaction-diffusion problems in Ω. Drawing on analytic regularity and corresponding hp-FEM in Ω for these reaction-diffusion problems with robust, exponential convergence as developed in [7,[27][28][29], we establish here exponential convergence rate bounds for solutions of (1.4).As we showed in [7,27], the singular perturbation character of the reaction-diffusion problems in Ω mandates both, geometric corner mesh refinement and anisotropic geometric boundary mesh refinement to resolve the algebraic corner and boundary singularities that occur in solutions to (1.4).
Before proceeding to the main part of this paper, we briefly recall the localization due to Caffarelli-Silvestre and the contour integral representation of Balakrishnan [6].

Caffarelli-Silvestre extension
In [14] the (full space) fractional Laplacian L s was localized via a singular elliptic PDE depending on one extra variable and thus represented as Dirichletto-Neumann problem for an elliptic problem in a half-space.Cabré and Tan [13] and Stinga and Torrea [40] extended this to bounded domains Ω and more general operators, thereby obtaining an extension posed on the semi-infinite cylinder C := Ω × (0, ∞).Their extension is given by the local boundary value problem where [14,40].The socalled conormal exterior derivative of U at Ω × {0} is The limit in (1.6) is in the distributional sense [13,14,40].Fractional powers of L s in (1.4) and the Dirichlet-to-Neumann operator of problem (1.5) are related by [14,15] We write x = (x , y) ∈ C with x ∈ Ω and y > 0. For D ⊂ R d × R + , we define L 2 (y α , D) as the Lebesgue space with the measure y α dx and H 1 (y α , D) as the weighted Sobolev space equipped with the norm . (1.9) To investigate (1.5) we include the homogeneous boundary condition on the lateral boundary ∂ L C by setting The bilinear form a C : is continuous and coercive on (1.12) For w ∈ H 1 (y α , C) we denote by tr Ω w its trace on Ω × {0}, which connects the spaces . (1.13) With these definitions at hand, the weak formulation of (1.5) is to find Remark 2 (regularity of U for s = 1/2) In the special case s = 1/2 and A = Id in (1.4) the operator L in the CS extension (1.5) in C coincides with the Laplacian in C. Therefore, the solution U will, in general, exhibit algebraic singularities on ∂Ω, even if ∂Ω is smooth.

Balakrishnan Formula
The second approach we take is via the spectral integral representation of fractional powers of elliptic operators going back to [6].For 0 < s < 1 and L = −div(A∇) with homogeneous Dirichlet boundary conditions the bounded linear operator L −s : H −s (Ω) → H s (Ω) admits the following representation with c B = π −1 sin(πs): (1.15) The representations (1.15) were used in [9,10] in conjunction with an exponentially convergent, so-called sinc quadrature approximation of (1.15) (see, e.g., [39] for details) and an h-version Finite Element projection in Ω to obtain numerical approximations of the fractional diffusion equation (1.4) in Ω.Here, we generalize the results in [9,10] to the hp-FEM, establishing exponential convergence rates in polygonal domains Ω for the resulting sinc BK-FEM for data A and f that are analytic in Ω (cf.(1.1)) without boundary compatibility of f .

Outline
The outline of the remainder of the paper is as follows.The following Section 2 describes the hp-FE spaces and Galerkin methods for (1.5) based on tensor products of discretizations in the x and the y variable.Section 3 develops the diagonalization of the hp-FE semi-discretization in the extended variable.In particular, in Section 3.1 we prove exponential convergence of an hp-FE semidiscretization in (0, ∞).The diagonalization reduces the semidiscrete approximation of the CS-extended, localized problem to a collection of decoupled, linear second order reaction-diffusion problems in Ω.
Section 4 presents the exponential convergence results from [7] of hp-FE discretizations of linear, second order singularly perturbed reaction-diffusion equations in Ω and establishes robust (with respect to the perturbation parameter ε) exponential convergence results for these.
Section 5 completes the proof of exponential convergence for the Extended hp-FEM by applying the hp-FEM from Section 4 in Ω for the reaction-diffusion problems obtained from the diagonalization process in Section 3. Section 5 presents in fact two distinct hp-FE discretizations: a pure Galerkin method (Case B) and a method based on discretizing after diagonalization each decoupled problem separately (Case A).The latter approach features slightly better complexity estimates.
Section 6 is devoted to the analysis of the sinc BK-FEM.There once more the numerical approximation of the fractional Laplacian is reduced to the numerical solution of a sequence of local linear, second order reaction-diffusion problems in Ω. Applying exponential convergence bounds for sinc approximation and for hp-FEM for reaction-diffusion problems in Ω in Section 4 from [7], once again exponential convergence for the resulting sinc BK-FEM for the spectral version of the fractional diffusion operator is established.As in Section 5, we separately discuss the possibilities of approximating the solutions of the decoupled problems from the same space (Case B) or from different spaces (Case A).Section 7 has numerical experiments verifying the theoretical convergence results.Section 8 has a summary and outlines several generalizations and directions for further research.In particular, we address in Section 8.1 the extension to fractional diffusion on manifolds.In Section 8.2, we discuss several exponential bounds on Kolmogoroff n-widths of solution sets for spectral diffusion in polygons that follow from our results.

Notation
Constants C, γ, b may be different in each occurence, but are independent of critical parameters.We denote by S := (0, 1) 2 the reference square and by T := {(ξ 1 , ξ 2 ) |, 0 < ξ 1 < 1, 0 < ξ 2 < ξ 1 } the reference triangle.Sets of the form {x = y}, {x = 0}, {x = y} etc. refer to edges and diagonals of S and analogously {y ≤ x} = {(x, y) ∈ S | y ≤ x}.P q denotes the space of polynomials of total degree q and Q q the tensor product space of polynomial of degree q in each variable separately.

hp-FEM Discretization
In this section, we introduce some hp-FEM space in both the x and the yvariable on which the Extended hp-FEM will be based.In particular, we introduce the geometric meshes G M geo,σ that are used for the discretization in the y-variable.
where P r denotes the space of polynomials of degree r.We will primarily work with the following piecewise polynomial space For constant polynomial degree r i = r ≥ 1, i = 1, . . ., M , we set S r {Y } ((0, Y ), G M ).Henceforth, we abbreviate Of particular interest will be geometric meshes G M geo,σ on [0, Y ], with M elements and grading factor σ ∈ (0, 1): we consider a linear polynomial degree vector r = {r i } M i=1 with slope s > 0 which is defined by For geometric meshes and linear degree vectors we set with constants implied in ∼ depending on s > 0.

hp-FEM in Ω
In the polygon Ω, we consider Lagrangian FEM of uniform1 polynomial degree q ≥ 1 based on regular triangulations of Ω denoted by T .We admit both triangular and quadrilateral elements K ∈ T , but do not assume shape regularity.
In fact, as we shall explain in Section 4 ahead, anisotropic mesh refinement towards ∂Ω will be required to resolve singularities at the singular support ∂Ω that are generically present in solutions of fractional PDEs (cf.Remark 2).We introduce, for a regular (in the sense of [27, Def.2.4.1])triangulation T of Ω comprising curvilinear triangular or quadrilateral elements K ∈ T with associated analytic element maps F K : K → K (where K ∈ { T , S} is either the reference triangle or square depending on whether K is a curvilinear triangle or quadrilateral) the FE space Here, for q ≥ 1, the local polynomial space V q ( K) = P q if K = T and V q ( K) =

Tensor product hp-FE approximation
One hp-FE approximation of the extended problem (1.5) will be based on the finite-dimensional tensor product spaces of the form where T is a regular triangulation of Ω.To analyze this method, we consider semidiscretizations based on the following (infinite-dimensional, closed) Hilbertian tensor product space: Here, the argument C Y indicates that spaces of functions supported in Ω × [0, Y ] are considered.Galerkin projections onto the spaces V q,r h,M (T q , G M ) and V r M (C Y ) with respect to the inner product a C (•, •) are denoted by G q,r h,M and G r M , respectively.For the CS-extension U , i.e., the solution of (1.14), the Galerkin projections G q,r h,M U and G r M U are characterized by (2.9)

Approximation based on semidiscretization in y
A key step in the hp-FE discretization in (0, Y ) is, as in [8], the diagonalization of the semidiscretized, truncated extension problem with solution G r M U given by (2.9).

Exponential
As in [8,23], we exploit the analytic regularity of the extended solution U with respect to the extended variable y.It results in exponential convergence of the hp-semidiscretization error , and consider the geometric mesh G M geo,σ on (0, Y ) and the linear degree vector r with slope s > 0. Let U be given by (1.14) and G r M U be the Galerkin projection onto V r M (T , G M ) given by (2.9).Then there exist C, b > 0 (depending solely on s, L, c 1 , c 2 , σ, ν, s) such that Furthermore, (3.1) also holds for constant polynomial degree r = (r, . . ., r) if Proof The statement is a slight generalization of [8, Lemma 13].In [8, Lemma 13], it is stated that the slope s has to satisfy s ≥ s min for some suitable s min > 0. Inspection of the proof shows that this condition can be removed.Specifically, using [2, Thm.
) provided r is a linear degree vector with suitably chosen slope.
The error bound (3.1) shows that up to an exponentially small (with respect to Y ) error introduced by truncation of (0, ∞) at Y , the solution U can be approximated by the solution of a local problem on the finite cylinder C Y .

Diagonalization
Diagonalization, as introduced in [8], refers to the observation that the solution G r M U of the semidiscrete problem (2.9) can be expressed in terms of M solutions U i ∈ H 1 0 (Ω), of M decoupled, linear local 2nd order reaction-diffusion problems in Ω.As the eigenvalues µ i in the corresponding eigenvalue problem (3.2) ahead govern the length scales in the local reaction-diffusion problems in Ω (3.5) (which, in turn, will be crucial in the mesh-design for the hp-FEM in Ω), it is of interest to know their asymptotic behavior.We investigate this in Lemma 2 below.Diagonalization is based on the explicit representation for the semidiscrete solution U M obtained from the following generalized eigenvalue problem, introduced in [8, Sec.6], and proposed earlier in [21], which reads: find We may expand the semidiscrete approximation G r M U as The coefficient functions U i ∈ H 1 0 (Ω) satisfy a system of M decoupled linear reaction-diffusion equations in Ω: Here v i denotes the i-th eigenfunction of the eigenvalue problem (3.2), (3.3) and with a Ω as introduced in (1.2).Due to the biorthogonality (3.3) of the discrete The following bounds on the µ i were shown in [8, Lemma 14] for the special case of geometric meshes G M geo,σ and linear degree vectors: Lemma 2 (properties of the eigenpairs, [8, Lemma 14]) Let {G M geo,σ } M ≥1 be a sequence of geometric meshes on (0, Y ) and r a linear polynomial degree vector with slope s > 0.

Assume that the truncation parameter
Then, there exists C > 1 (depending on c 1 , c 2 and on σ ∈ (0, 1)) such that there holds for for every M ∈ N for the eigenpairs given by (3.2), (3.3) (3.8)

Fully discrete approximation
The full discretization is obtained by approximating the functions U i of (3.5) from finite-dimensional spaces.Let T i , i = 1, . . ., M, be regular triangulations in Ω and q ∈ N. Let Π q i : H 1 0 (Ω) → S q 0 (Ω, T i ) denote the Ritz projectors for the bilinear forms a µi,Ω , which are characterized by In terms of the projections Π q i we can define the fully discrete approximation By combining (3.5) and (3.9), the functions Π q i U i ∈ S q 0 (Ω, T i ) are explicitly and computably given as the solutions of In view of (3.7), we have the following representation of the difference between the semidiscrete approximation G r M U and the fully discrete approximation U h,M : Lemma 3 Let U h,M be given by (3.10).Then: Concerning the meshes T i , we distinguish two cases in this work: Case A: The meshes T i , i = 1, . . ., M, possibly differ from each other.Case B: The meshes T i , i = 1, . . ., M, coincide.That is, all coefficient functions U i in the semidiscrete solution (3.4) are approximated from one common hp-FE space S q 0 (Ω, T ).In Case B the approximation U h,M actually coincides with the Galerkin projection Lemma 4 (error representation, [8, Lemma 12]) Let (µ i , v i ) M i=1 be the eigenpairs given by (3.2), (3.3).For i = 1, . . ., M, let U i ∈ H 1 0 (Ω) be the solutions to (3.5).Consider Case B and let Π q i : H 1 0 (Ω) → S q 0 (Ω, T ) be the Galerkin projections given as in (3.9), with one common, regular triangulation T of Ω for i = 1, . . ., M. Let G r M U denote the solution to the semidiscrete problem (2.9).Then the tensor product Galerkin approximation Lemma 4 shows that in Case B, the Galerkin projection of U into the tensor product space V q,r h,M (T , G M ) coincides with the approximation U h,M defined in (3.10) in terms of the decoupling procedure.Hence, the decoupling procedure is not essential for numerical purposes in Case B, although it has algorithmic advantages.In contrast, Case A relies on the decoupling in an essential way.In both cases, the exponenial convergence result below will make use of the error estimates of Lemmas 3, 4 obtained by the diagonalization process.
It is advisable to choose the spaces T i in case Case A such that the functions U i can be approximated well from S q 0 (Ω, T i ) in the norm • µi,Ω .Correspondingly in Case B, the commmon space T should be chosen such that each U i can be approximated well from S q 0 (Ω, T ).The bounds (3.8) indicate that, for large M , most of the reaction-diffusion problems (3.5) are singularly perturbed.Hence we design in the following Section 4 hp-FE approximation spaces in Ω which afford exponential convergence rates that are robust with respect to the singular perturbation parameter.

hp-FE Approximation of singular perturbation problems
In the exponential convergence rate analysis of tensorized hp-FEM for the CS extension (Extended hp-FEM) as well as for the ensuing (see Section 6 ahead) sinc BK-FEM approximation, a crucial role is played by robust exponential convergence rate bounds for hp-FEM for singularly perturbed, reaction-diffusion problems in curvilinear polygonal domains Ω. Specifically, we consider the hp-FE approximation of the local reaction-diffusion problem in Ω, where we assume ess inf x ∈Ω c(x ) ≥ c 0 > 0 and A, c, and f are analytic on Ω and A is symmetric, uniformly positive definite.(4.2) We note again that (4.1) does not imply any kind of boundary compatibility of f at ∂Ω (cf.Remark 1).We assume Ω to be scaled so that diam(Ω) = O(1).
Then, for small ε > 0, the boundary value problem (4.1) is a so-called "ellipticelliptic" singular perturbation problem.Under the assumptions (4.2), for every ε > 0 problem (4.1) admits a unique solution u ε ∈ H 1 0 (Ω).In general, u ε exhibits, for small ε > 0, boundary layers near ∂Ω whose robust numerical resolution (i.e., with error bounds whose constants are independent of ε) requires anisotropically refined meshes aligned with ∂Ω (see [17,28,33] and the references there).In addition, the corners of Ω induce point singularities in the (analytic in Ω) solution u ε .In the context of hp-FEM under consideration here, their efficient numerical approximation mandates geometric mesh refinement near the corners.
In the present section, we consider the hp-FEM approximation of u ε that features exponential convergence for two different types of meshes: a) geometric boundary layer meshes in Section 4.2 and b) admissible boundary layer meshes in Section 4.3.In both cases the error estimates are of the form O(e −bq +e −b L ), with the constant hidden in O(•) independent of ε, L, and q and where q is the polynomial degree employed and L measures the number of layers of geometric refinement towards the vertices or edges of Ω.The difference in these two types of meshes is that "admissible boundary layer meshes" are strongly ε-dependent with geometric refinement towards the vertices and only a single layer of thin elements of width O(qε) near ∂Ω to resolve the boundary layer.The number of elements is then O(L) leading to a number of degrees of freedom N = O(Lq 2 ).In contrast, geometric boundary layer meshes are based on geometric, anisotropic refinement towards the edges and corners of Ω.As we show in [7], hp-FEM on such meshes afford exponential convergence for boundary layers with multiple scales.The total number of elements in geometric boundary layer meshes with L layers is O(L 2 ).Combined with local FE spaces of polynomial degree q, this results in a number of degrees of freedom N = O(L 2 q 2 ).Whereas admissible boundary layer meshes are designed to approximate boundary layers of a single, given length scale ε, geometric boundary layer meshes afford concurrent, robust and exponentially convergent approximations of boundary layers with multiple length scales in Ω.These arise, e.g., upon semidiscretization in the extended variable as is evident from (3.5).

Macro triangulation. Geometric boundary layer mesh
We do not consider the most general meshes with anisotropic refinement, but confine the hp-FE approximation theory to meshes generated as push-forwards of a small number of so-called mesh patches.This concept was used in the error analysis of hp-FEM for singular perturbations in [27,Sec. 3.3.3]and in [17].Specifically, we assume given a fixed macro-triangulation analytic patch maps (to be distinguished from the actual element maps) F K M : S = (0, 1) 2 → K M that satisfy the usual compatibility conditions.I.e., T M does not have hanging nodes and, for any two distinct elements K M 1 , K M 2 ∈ T M that share an edge e, their respective element maps induce compatible parametrizations of e (cf., e.g., [27,Def. 2.4.1] for the precise conditions).Each element of the fixed macro-triangulation T M is further subdivided according to one of the refinement patterns in Definition 1 (see also [27,Sec. 3.3.3]or [17]).The actual triangulation is then obtained by transplanting refinement patterns on the square reference patch into the physical domain Ω by means of the element maps F K M of the macro-triangulation.That is, for any element K ∈ T , the element map F K is the concatenation of an affine map-which realizes the mapping from the reference square or triangle to the elements in the patch refinement pattern and will be denoted by A K -and the patch map (which will be denoted by The following refinement patterns were introduced in [7, Def.2.1, 2.3].They are based on geometric refinement towards a vertex and/or an edge; the integer L controls the number of layers of refinement towards an edge whereas n ∈ N measures the refinement towards a vertex. Definition 1 (Catalog P of refinement patterns, [7, Def.2.1]) Given σ ∈ (0, 1), L, n ∈ N 0 with n ≥ L the catalog P consists of the following patterns: 1.The trivial patch: The reference square S = (0, 1) 2 is not further refined.The corresponding triangulation of S consists of the single element: 2. The geometric boundary layer patch Ť BL,L geo,σ : S is refined anisotropically towards { y = 0} into L elements as depicted in Fig. 2 (top left).The mesh Ť BL,L geo,σ is characterized by the nodes (0, 0), (0, σ i ), (1, σ i ), i = 0, . . ., L and the corresponding rectangular elements generated by these nodes.3. The geometric corner patch Ť C,n geo,σ : S is refined isotropically towards (0, 0) as depicted in Fig. 2 (top middle).Specifically, the reference geometric corner patch mesh Ť C,n geo,σ in S with geometric refinement towards (0, 0) and n layers is given by triangles determined by the nodes (0, 0), and (0, σ i ), (σ i , 0), (σ i , σ i ), i = 0, 1, . . ., n. 4. The tensor product patch Ť T,L,n geo,σ : S is triangulated in S 1 := (0, σ L ) 2 and S 2 := S\ S 1 separately as depicted in Fig. 2 (bottom left).The triangulation of S 1 is a scaled version of Ť C,n−L geo,σ characterized by the nodes (0, 0), (0, σ i ), (σ i , 0), (σ i , σ i ), i = L, . . ., n.The triangulation of S 2 is characterized by the nodes (σ i , σ j ), i, j = 0, . . ., L.
Remark 3 We kept the list of possible patch refinement patterns in Definition 1 small in order to reduce the number of cases to be discussed for the hp-FE error bounds.A larger number of refinement patterns could facilitate greater flexibility in mesh generation.In particular, the reference patch meshes do not contain general quadrilaterals but only (axiparallel) rectangles; this restriction is not essential but leads to some simplifications in the hp-FE error analysis in [7].
The addition of the diagonal line in the reference corner, tensor, and mixed patches is done to be able to apply the regularity theory of [27] and probably not necessary in actual computations.We also mention that with additional constraints on the macro triangulation T M the diagonal line could be dispensed with, [7].
The following definition of the geometric boundary layer mesh T L,n geo,σ formalizes the requirement on the meshes that anisotropic refinement towards ∂Ω is needed as well as geometric refinement towards the corners.Given σ ∈ (0, 1), L, n ∈ N 0 with n ≥ L, a mesh T L,n geo,σ is called a geometric boundary layer mesh if the following conditions hold: geo,σ is obtained by refining each element K M ∈ T M according to the finite catalog P of structured patch-refinement patterns specified in Definition 1, governed by the parameters σ, L, and n.

T L,n
geo,σ is a regular triangulation of Ω, i.e., it does not have hanging nodes.Since the element maps for the refinement patterns are assumed to be affine, this requirement ensures that the resulting triangulation satisfies [27, Def.2.4.1].
For each macro-patch K M ∈ T M , exactly one of the following cases is possible: Then the trivial patch is selected as the reference patch.4. K M ∩ ∂Ω is a single point.Then two cases can occur: (a) K M ∩ ∂Ω = {A j } for a vertex A j of Ω.Then the corresponding reference patch is the corner patch Ť C,n geo,σ with n layers of refinement towards the origin O. Additionally, F K M (O) = A j .(b) K M ∩ ∂Ω = {P}, where the boundary point P is not a vertex of Ω.
Then the refinement pattern is the corner patch Ť C,L geo,σ with L layers of geometric mesh refinement towards O. Additionally, it is assumed that F K M (O) = P ∈ ∂Ω. 5. K M ∩ ∂Ω = e for an edge e of K M and neither endpoint of e is a vertex of Ω.Then the refinement pattern is the boundary layer patch Ť BL,L geo,σ and additionally F K M ({ y = 0}) ⊂ ∂Ω. 6.K M ∩∂Ω = e for an edge e of K M and exactly one endpoint of e is a vertex A j of Ω.Then the refinement pattern is the mixed layer patch Ť M,L,n geo,σ and additionally F K M ({ y = 0}) ⊂ ∂Ω as well as F K M (O) = A j .7. Exactly two edges of a macro-element K M are situated on ∂Ω.Then the refinement pattern is the tensor patch Ť T,L,n geo,σ .Additionally, it is assumed that Finally, the following technical condition ensures the existence of certain meshlines: 8.For each vertex A j of Ω, introduce a set of lines Let Γ j , Γ j+1 be the two boundary arcs of Ω that meet at A j .Then there exists a line e ∈ such that the interior angles ∠(e, Γ j ) and ∠(e, Γ j+1 ) are both less than π.Then there are constants C, b > 0, β ∈ [0, 1) (depending solely on the data A, c, f , Ω, on the parameter c 1 , and on the analyticity properties of the patch-maps of the macro-triangulation T M ) such that the following holds: If ε ∈ (0, 1] and L satisfy the (boundary layer) scale resolution condition then, for any q, n ∈ N, the solution u ε ∈ H 1 0 (Ω) of (4.1) can be approximated from S q 0 (Ω, T L,n geo,σ ) such that inf v∈S q 0 (Ω,T L,n geo,σ ) Proposition 1 is restricted to ε ∈ (0, 1].For ε ≥ 1, the solution u ε of (4.1) does not have boundary layer but merely corner singularities.Hence, by Remark 4 meshes with fixed L are appropriate.In particular, the boundary layer scale resolution condition (4.3) is not required: Proposition 2 Assume the hypotheses on Ω and the data A, c, f as in Proposition 1.Let {T L,n geo,σ } L≥0,n≥L be a sequence of geometric boundary layer meshes 2 .
4.3 hp-FE approximation of singularly perturbed problems on admissible meshes T L,q min,λ (ε) in Ω In Proposition 1, the solution u ε is approximated on patchwise geometric meshes.These meshes are able to capture boundary layers (and corner layers) on a whole range of singular perturbation parameters ε: as long as a lower bound for ε is known and provided that geometric mesh refinement resolves all scales, robust exponential convergence is assured.
In particular, for L ∼ q, there are constants b , C > 0 such that ∀ε ≥ 1, q, L ∈ N : inf v∈S q 0 (Ω,T L,q min,λ (ε)) It is worth pointing out the following differences between the approximation on geometric boundary layer meshes T L,n geo,σ and on the minimal admissible boundary layer meshes T L,q min,λ : a) the use of the mesh T L,n geo,σ requires the scale resolution condition (4.3).It requires L | ln ε| so that the approximation result Proposition 1 depends (weakly) on ε. b) Selecting n L q in Proposition 1 yields convergence O(exp(−b 4  √ N )) whereas the choice L q in Proposition 3 yields the better convergence behavior O(exp(−b 3 √ N )).c) The meshes T L,q min,λ are designed to approximate a single scale well whereas the meshes T L,n geo,σ are capable to resolve a range of scales.d) The meshes T L,q min,λ rely on a suitable choice of the parameter λ whereas the geometric boundary layer meshes T L,n geo,σ do not have parameters that need to be suitably chosen.

Exponential Convergence of Extended hp-FEM
Based on the hp semidiscretization in the extended variable combined with the diagonalization in Section 3, we use the hp-approximation results from Section 4 to prove exponential convergence of hp-FEM for the CS-extended problem (1.14).
As is revealed by the diagonalization (3.5), the y-semidiscrete solution G r M U contains M separate length scales associated with the eigenvalues µ i , i = 1, ..., M. The solutions U i ∈ H 1 0 (Ω) of the resulting M many independent, linear second-order reaction-diffusion problems in Ω exhibit both, boundary layers and corner singularities.
In Case A, which we discuss in Section 5.1, we employ for each i a "minimal" hp-FE space in Ω that resolves boundary-and corner layers appearing in the U i due to possibly large/small values of µ i .Mesh design principles for such "minimial" FE spaces that may resolve a single scale of a singularly perturbed problem have already been presented in, e.g., [25,27,28,37,38]; the specific choice T L,q min,λ (ε) has been discussed in Propositions 3, 4 and will be used in our analysis.
In Case B, which we discuss in Section 5.2, we discretize these decoupled, reaction-diffusion problems by one common hp-FEM in the bounded polygon Ω, which employs both, geometric corner refinement as well as geometric boundary refinement, as in [27,28].Due to the need to obtain FE solutions for all µ i in one common FE space in Ω, however (in order that the sum (3.13) belong to a tensor product hp-FE space), our analysis will provide one hp-FE space in Ω which will resolve all boundary and corner layers due to small parameters µ i near ∂Ω.As we shall show, in Case B the total number of DOFs is larger than in Case A.

Exponential Convergence I: Diagonalization and Minimal Meshes
The robust exponential convergence result Proposition 3 allows us to establish, in conjunction with the diagonalization (3.2)-(3.4), a first exponential convergence result in Case A of Section 3. We consider the following numerical scheme, which relies on the "minimal boundary layer meshes" T L,q min,λ (ε) from [27, Sec.2.4.2]already discussed in Proposition 3: (1) Select Y with c 1 M ≤ Y ≤ c 2 M and consider the space S r {Y } ((0, Y ), G M geo,σ ) for the geometric mesh G M geo,σ on (0, Y ) with M elements and a linear degree vector r with slope s > 0.
(3) Select λ > 0. Define U q,L i ∈ S q 0 (T L,q min,λ ( √ µ i )) as the solution of ∀v ∈ S q 0 (T L,q min,λ ( (5.1) (4) Define the approximation U q,L (x, y) := Mgeo i=1 v i (y)U q,L i (x).For the approximation error U − U q,L we have: Theorem 1 Let Ω ⊂ R 2 be a curvilinear polygon with J vertices as described in Section 1.1.Let A, f satisfy (1.1) and let A be uniformly symmetric positive definite on Ω. Fix positive constants c 1 , c 2 , and s.
Then there are constants C, b, b , b , λ 0 > 0 (depending on Ω, A, c, and the parameters characterizing the mesh family T L,q min,λ ) such that for any λ ∈ (0, λ 0 ] there holds for all q, M ≥ 1, L ≥ 0 In particular, for q L M p, denoting U p := U q,L with this choice of q and L, and the total number of degrees of freedom N = i dim S q 0 (Ω, T L,q min,λ ( 3) where the constant b depends additionally on the implied constants in q L M .Remark 5 The approximation result (5.3) still holds if the linear degree vector r in the definition of U q,L is replaced with a constant polynomial degree r ∼ M .

Proof Step 1 (semidiscretization error):
The analyticity of f on Ω implies f ∈ H −s+ν (Ω) for any ν ∈ (0, 1/2 + s).Hence, by (3.1), the semidiscretization error (5.4) Step 2 (representation of G r M U ): The semidiscrete approximation G r M U may be expressed in terms of the eigenbasis where the function U i solve by (3.5) Step 3 : For every i = 1, . . ., M geo , and for every q ∈ N, approximate the U i ∈ H 1 0 (Ω) by its Galerkin approximation U q,L i ∈ S q 0 (Ω, T L,q min,λ ( √ µ i )) ⊂ H 1 0 (Ω).That is, U q,L i = Π q i U i is the a µi,Ω (•, •)-projection of U i given by (3.9).It is the best approximation to U i in the corresponding energy norm and satisfies By linearity of Π q i and the analyticity of f Propositions 3, 4 (depending on whether µ i ≤ 1 or µ i > 1) and Lemma 2 Step 4 (Proof of (5.2)):With the approximations U q,L i the approximation U q,L of U is given by We note that M geo ∼ M 2 .Combining this last estimate with (5.4) yields the second estimate in (5.2).The first estimate in (5.2) expresses the continuity of the trace operator at y = 0.
Step 5 (complexity estimate): ) and the fact that dim S q 0 (Ω, T L,q min,λ ( √ µ i )) ≤ CLq 2 = O(q 3 ) as well as the assumption q L M p, we arrive at a total problem size N = O(q 5 ).Absorbing algebraic factors in the exponentially decaying one in (5.2) then yields (5.3).

Exponential Convergence II: Geometric Boundary Layer Meshes
In this section, we show that exponential convergence of a Galerkin method for (1.14) can be achieved by a suitable choice of meshes T and G M in the tensor product space V q,r h,M (T , G M ) of (2.6).That is, we place ourselves in Case B in Section 3.2.For the discretization in y, we select again the spaces S r {Y } ((0, Y ), G M geo,σ ) with Y ∼ M and the linear degree vector r with slope s.The hp-FE discretization in Ω is based on the space S q 0 (Ω, T L,n geo,σ ) with the geometric boundary layer mesh T L,n geo,σ in Definition 2. Recall that G q,r h,M U denotes the Galerkin projection of the solution ).In Theorem 2 below, we will focus on the case q L M and the corresponding Galerkin projection is denoted U p T P .
Remark 6 (i) In contrast to the procedure of Case A in the preceeding Section 5.1, precise knowledge of the length scales √ µ i is not necessary.
(ii) The diagonalization procedure may be carried out numerically and results in decoupled reaction-diffusion problems, affording parallel numerical solution.(iii) The linear degree vector r could be replaced with a constant degree r ∼ M , and Theorem 2 will still hold.
For the tensor-product hp-FEM in C Y we also have an exponential convergence result: Theorem 2 Let Ω ⊂ R 2 be a curvilinear polygon with J vertices as described in Section 1.1.Let A, f satisfy (1.1) and let A be uniformly symmetric positive definite on Ω. Fix a slope s > 0. Set With these choices, denote by U p T P the Galerkin projection of U onto the tensor product hp-FE space V q,r h,M (T L,n geo,σ , G M geo,σ ) = S q 0 (Ω, T L,n geo,σ )⊗S r {Y } ((0, Y ), G M geo,σ ).
Then N := dim V q,r h,M (T L,n geo,σ , G M geo,σ ) = O(p 6 ) and there are constants C, b, b > 0 depending only on Ω, A, f , the macro triangulation T M underlying the geometric boundary layer meshes T L,n geo,σ , the slope s, and the implied constants in (5.5) such that Proof The proof of this result is structurally along the lines of the proof of Theorem 1.We omit details and merely indicate how the scale resolution condition (4.3) is now accounted for.We note that for fixed Y and M < M , we have that the spaces S r {Y } ((0, Y ), G M ) and ) (if the same slope s for the linear degree vector is chosen).Hence, the Galerkin error for the approximation from the space S q 0 (Ω, ), and we therefore focus on bounding the approximation error for S q 0 (Ω, T L,n geo,σ ) ⊗ S r {Y } ((0, Y ), G M geo,σ ).We select M of the form M = ηM for some η to be chosen below.For ease of notation, we simply set M = ηM .
By Lemma 2 we have that the smallest length scale of the singlarly perturbed problems for the space S q 0 (T L,n geo,σ and that scale resolution condition (4.3) therefore reads Since L M , we see that (5.6) can be satisfied for some fixed c 1 provided η is suitably chosen.The approximation of U from S q 0 (T L,n geo,σ ) ⊗ S r {Y } (G M geo,σ ) then follows by arguments very similar to those of the proof of Theorem 1.

Exponential Convergence of sinc BK-FEM
The hp-since BK FEM is based on exponentially convergence, so-called "sinc" quadratures to the Balakrishnan formula as described in Section 1.5 and (1.15) (see [9][10][11] and the references there).We briefly review the corresponding exponential convergence results in Section 6.1.The numerical realization of the sinc quadrature approximation of (6.1) leads again to the numerical solution of decoupled, local linear reactiondiffusion problems in Ω.These boundary value problems are again singularly perturbed.Accordingly, we discuss two classes of hp-FE approximations for their numerical solution: In Section 6.2.1, we discuss Case A, which is based on the minimal boundary layer meshes T L,q min,λ in Ω.In Section 6.2.2, we detail Case B, where geometric boundary layer meshes in Ω are employed.The latter allow one common hp-FEM for all values of parameters arising from the sinc quadrature approximation of (6.1).

Sinc quadrature approximation
The above integral (6.1) can be discretized by so-called "sinc" quadratures (see, e.g., [9,39]).To that end, we define for K ∈ N For f ∈ L 2 (Ω) ⊂ H −s (Ω) for every 0 < s < 1, the (semidiscrete) sinc quadrature approximation u M of u = L −s f ∈ H s (Ω) as represented in (6.1) reads with ε j := e −yj /2 = e j/(2 We note that for any k, we have that is a bounded linear map for any 0 < s < 1.The semidiscretization error L −s − Q −s k (L) is bound in [9]: Proposition 5 ( [9, Thm.3.2]) For f ∈ L 2 (Ω) and for every 0 ≤ β < s with 0 < s < 1 denoting the exponent of the fractional diffusion operator in (1.4), there exist constants b, C > 0 (depending on β, s, Ω, and L) such that for every k > 0 as in (6.2) holds, with Remark 7 Sinc approximation formulas such as (6.3) have a number of parameters which can be optimized in various ways.The error bound in Proposition 5 is merely one particular choice (the so-called "balanced" choice of parameters), which is sufficient for the exponential sinc error bound (6.4).Other choices yield analogous (exponential) sinc error bounds, with possibly better numerical values for the constants b, C > 0 in (6.4).We point out that we make such a choice in our numerical examples in (7.3) and refer to [9,Rem. 3.1] for details.

hp-FE approximation in Ω
The sinc approximation error bound (6.4) implies exponential convergence of the sinc quadrature sum (6.3), which we write as Here, the w j ∈ H 1 0 (Ω) are solutions of the 2K + 1 reaction-diffusion problems With the bilinear form a ε 2 j ,Ω (•, •) from (3.6), their variational formulations reads: find w j ∈ H 1 0 (Ω) such that ∀v ∈ H 1 0 (Ω) : a ε 2 j ,Ω (w j , v) = (f, v) .(6.7) The reaction diffusion problems (6.6) are again of the type (3.5) for which exponentially convergent hp-FE approximations were presented in Section 5, from [27] and [7].A fully discrete sinc BK-FEM approximation is constructed by replacing w j in (6.5), (6.6) by one of the hp-FE approximations discussed in Section 5.As in the case of the Extended hp-FEM, also for the sinc BK-FEM one can distinguish Case A, in which each problem (6.7) is discretized using a different hp-FE space, and Case B, where all problems (6.7) are discretized by the same hp-FE space in Ω.

Case A
We discretize the singularly perturbed problems (6.7) with length scales ε j using the spaces T L,q min,λ (ε j ).That is, denoting the resulting approximations generically by w hp j , defined by: for |j| ≤ K, find w hp j ∈ S q 0 (Ω, T L,q min,λ (ε j )) such that ∀v ∈ S q 0 (Ω, T L,q min,λ (ε j )) : The hp-FE approximations w hp j are well-defined.Replacing in (6.5) the w j by their hp-FE approximations, we obtain the sinc BK-FEM approximation of the (inverse of) the fractional diffusion operator L s : To bound the error u − Q −s k (L hp )f H s (Ω) , we write . For the first term, the sinc approximation error, we use the error bound (6.4) with β = s/2.Using D(L s/2 ) = H s (Ω) for 0 < s < 1, we obtain from (6.4) and from k K −1/2 (see (6.2)) the bound To bound the second term, definition (6.5) and the triangle inequality imply To invoke the hp-error bound (4.4) with the norm (3.7), we use the interpolation inequality We apply this to each term in (6.11) and, using the definition (3.7) of the norm w ε 2 ,Ω , and for all ε ≥ 0, arrive at We split the summation indices as I + ∪ I − , i.e., I + := {|j| ≤ K} ∩ {j > 0} and With Proposition 4, we estimate the sum over j ∈ I + according to We select L q (i.e., the number L of mesh-layers proportional to the polynomial degree q ≥ 1) with proportionality constant independent of ε j .Furthermore, we note k = 1/ √ K and select K q so that Combining the error bounds (6.10) and (6.15), and suitably adjusting the constant b > 0 in the exponential bounds, we arrive at Given that the approximation u K,hp involves the solution of O(K) = O(q 2 ) reaction-diffusion problems, each of which requires O(q 3 ) DOF, the error bound (6.16) in terms of the total number of degrees of freedom N DOF reads with constants b, C > 0 that are independent of N DOF .We have thus shown: Theorem 3 Let Ω be a curvilinear polygon as defined in Section 1.1, let A, f satisfy (1.1), and let A be uniformly symmetric positive definite on Ω.Let u be the solution to (1.4), and let its fully discrete approximation be given by the sinc BK-FEM approximation (6.5) in conjunction with the hp-FE approximation of w hp j in (6.8) with the hp-FE spaces S q 0 (Ω, T L,q min,λ (ε j )) on the minimal boundary layer meshes T L,q min,λ (ε j ).Choose further the parameters q L K and let N = |j|≤K dim S q 0 (Ω, T L,q min,λ (ε j )) denote the total number of degrees of freedom.
Then there exists a λ 0 (depending on Ω, A, c, and the parameters characterizing the mesh family T L,q min,λ ) such that for any λ ∈ (0, λ 0 ] there are constants C, b > 0 (depending additionally on the implied constants in q L K) such that

Case B
Instead of approximating the problems (6.7) from individual spaces, one may approximate them from the same hp-FE space in Ω.Specifically, we define the approximations w hp j by: Find w hp j ∈ S q 0 (Ω, T L,n geo,σ ) such that ∀v ∈ S q 0 (Ω, T L,n geo,σ ) : Theorem 4 Let Ω be a curvilinear polygon as defined in Section 1.1 and assume that A, f satisfy (1.1) and that A is uniformly symmetric positive definite.Let u be solution to (1.4), and let its discrete approximation Q −s k (L hp )f be given by (6.5) in conjunction with (6.18).Fix c 1 > 0. Let K q L = n and let N = (2K +1) dim S q 0 (Ω, T L,n geo,σ ) ∼ q 6 denote the total number of degrees of freedom.
Then, under the scale resolution condition there are constants C, b > 0 (depending on Ω, A, f , c 1 , σ, and the analyticity properties of the macro triangulation) such that Proof The proof follows the arguments of Theorem 3. Since e −K/2 = min j ε j is the smallest scale, the scale resolution condition (6.19) ensures that Proposition 1 is applicable.

Numerical experiments
We consider the problem (1.4) with diffusion coefficient A = I, i.e., L s = (−∆) s .The domain Ω is chosen as either the unit square Ω 1 = (0, 1) 2 , the so-called L-shaped polygonal domain Ω 2 ⊂ R 2 determined by the vertices {(0, 0), (1, 0), (1, 1), (−1, 1), (−1, −1), (0, −1)}, or the square domain with a slit Ω 3 = (−1, 1) 2 \ (−1, 0] × {0}.As we are in particular interested in smooth, but possibly non-compatible data f in all the numerical examples we take Notice that, in this case, f is analytic on Ω but f ∈ H 1−s (Ω) only for s > 1/2 due to boundary incompatibility (cf.Remark 1).The exact solution is not known, so that the error is estimated numerically with reference to an accurate numerical solution.The error measure is always the functional where u fine is the numerical solution obtained on a fine mesh.Note that for the Galerkin method on the cylinder C (i.e., Extended hp-FEM in Case B) this error measure is equivalent to the energy norm if u fine is replaced by the exact solution u: where U p denotes the discrete solution in C Y .
In Figure 4 we show examples of the meshes used for the three domains.These are constructed using the Netgen/NGSolve package [35].For the square domain Ω = Ω 1 the resulting mesh is the geometric boundary layer mesh T L,L geo,σx with L = 4 and σ x = 1/4.The same parameters are used in Netgen/NGSolve to construct the meshes for the other two domains, with the resulting meshes diverging from the strict definition of T L,L geo,σ2 near the reentrant corners since these meshes are not constructed using mesh patches but instead by applying directly geometric refinement towards edges and vertices.Nevertheless we denote these meshes also by T L,L geo,σx and make use of the finite element spaces S q 0 (Ω, T L,L geo,σx ).Given a polynomial order p ≥ 1, in both approaches the finite element space in Ω is S q 0 (Ω, T L,L geo,σx ) with uniform polynomial degree q = p, number of levels L = p and σ x = 1/4.Next, we describe the parameters used in the hp-FEM on (0, Y ) and the quadrature in the Balakrishnan formula.
For the extended problem, on the geometric mesh G M geo,σy in (0, Y ) as defined in Section 2.1.1 we use FE-spaces S r {Y } ((0, Y ), G M geo,σy ).Given a polynomial degree p ≥ 1, in the definition of these spaces we use Y = 1 2 p, σ y = 1/4 and M = round(0.79p/s)3 , and a uniform degree vector r = (p, . . ., p).
For simplicity in the analysis of the sinc quadrature, we used a symmetric approximation (6.3).For the numerical experiments in order to obtain a more efficient scheme we have followed [9] to define the quadrature as with y = k and the number of quadrature points chosen as For the given polynomial order p ≥ 1, we set k = 4 3 p −1 .We now compare the convergence of the two schemes.We plot the error against the polynomial degree p and against N 1/2 ls , where N ls is the number of linear systems that need to be solved.The convergence curves for the square domain Ω 1 are shown in Figure 5, for the L-shaped domain Ω 2 in Figure 6, and for the slit domain Ω 3 in Figure 7.For all three domains we clearly see exponential convergence as the polynomial order is increased.Also, the Extended hp-FEM requires significantly fewer linear systems to be solved to achieve the same accuracy as the sinc BK-FEM .We should, however, also note that the eigenvalue problem (3.2) becomes ill-conditioned for increasing p and much higher accuracy than the one shown in the above figures cannot be obtained using our approach for the extension problem.No such accuracy limitations could be seen for the sinc approach.
Pulling back the problem (8.1) via χ into the (Euclidean) chart domain Ω ⊂ R 2 , the Dirichlet problem for the fractional power s ∈ (0, 1) of the surface diffusion (8.1) in Γ = χ(Ω) reduces to (1.5) where the bilinear form (1.2) and diffusion coefficient A are given by with G(x ) = Dχ(x ) : Ω → R 3×2 denoting the (assumed analytic in Ω) metric of M in chart χ.The real-analyticity of compositions, sums and product of real-analytic functions implies that A(x ) satisfies (1.1) in Ω, so that the ensuing mathematical results also apply to (1.4) with (8.1).exponential convergence of hp-FEM for reaction-diffusion problems in polygons from [7,27,28], an exponential convergence rate bound C exp(−b 6  √ N DOF ) with respect to the total number of degrees of freedom, N DOF , which are used in the tensor-product hp-FE discretization, is established in the fractional Sobolev norm H s (Ω).We add that the variational semi-discretization in Section 3 with respect to the extruded variable y offers the possibility for residual a posteriori error estimation.
The second discretization is based on the spectral integral representation of L −s due to Balakrishnan [6].A sinc quadrature discretization [39] approximates the spectral integral by an (exponentially convergent [10,39]) finite linear combination of solutions of decoupled elliptic reaction diffusion problems in Ω with analytic input data.Drawing once more on analytic regularity and robust exponential convergence of hp-FEM [7,[27][28][29], we prove exponential convergence also for this approach.A computable a posteriori bound for the semidiscretization error incurred for the sinc BK-FEM approach does not seem to be available currently.
The theoretical convergence rate bounds are verified in a series of numerical experiments.These show, in particular, that exponential convergence is realized in the practical range of discretization parameters.They also indicate a number of practical issues, such as conditioning or algorithmic steering parameter selection, which are beyond the scope of the mathematical convergence analysis.We point out that the proposed algorithms and the exponential convergence results extend in several directions: besides homogeneous Dirichlet boundary conditions, also mixed, Dirichlet-Neumann boundary conditions, and operators with a nonzero first order term could be considered.In either case, the proposed algorithms extend readily.The main result is the construction of hp-FE discretizations with robust exponential convergence rates for spectral fractional diffusion in polygonal domains Ω ⊂ R 2 .Similar results hold in bounded intervals Ω ⊂ R 1 (we refer to [8] for details).In polyhedral Ω ⊂ R 3 , the present line of analysis is also applicable; however, exponential convergence and analytic regularity of hp-FEM for reaction-diffusion problems in space dimension d = 3 does not appear to be available to date.We considered fractional powers only for self-adjoint, second-order elliptic divergenceform differential operators Lw = −div(A∇w) in Ω.The arguments for the sinc BK-FEM extend to non-selfadjoint operators which include first-order terms via [9], provided suitable hp-FEM for advection-reaction-diffusion problems in Ω are available (e.g.[24]).
The present analysis is indicative for achieving high, algebraic rate p+1−s of convergence in H s (Ω) by h-version FEM of fixed order p ≥ 1 in Ω.As in hp-FEM, this will require anisotropic mesh refinement aligned with ∂Ω, ie., so-called "boundary-layer" meshes.Several constructions are available (see, e.g., [38] for so-called "exponential boundary layer meshes" and [33] for socalled "Shiskin meshes").We refrain from developing details for this approach which can be analyzed along the lines of the present paper.

2. 1
Notation and FE spaces 2.1.1Meshes and FE spaces on (0, Y ) Given a truncation parameter Y and a mesh

8 ,
Eqn. (78), Rem.16] (or, referring alternatively to the extended preprint [2, Thm.3.13, Eqn.(3.21), Rem.3.14]) the result holds for any s > 0, with constant b = O(s) as s ↓ 0. The statement about the constant polynomial degree follows from the case of the linear degree vector since a) G r M U is the Galerkin projection of U , b) the minimization property of Galerkin projections, and c) the fact that the space S r {Y

Fig. 2 Fig. 3
Fig.2Catalog P of reference refinement patterns from[7].Top row: reference boundary layer patch Ť BL,L geo,σ with L layers of geometric refinement towards { y = 0}; reference corner patch Ť C,n geo,σ with n layers of geometric refinement towards (0, 0); trivial patch.Bottom row: reference tensor patch Ť T,L,n geo,σ with n layers of refinement towards (0, 0) and L layers of refinement towards { x = 0} and { y = 0}; reference mixed patch Ť M,L,n geo,σ with L layers of refinement towards { y = 0} and n layers of refinement towards (0, 0).Geometric entities shown in boldface indicate parts of ∂ S that are mapped to ∂Ω.These patch meshes are transported into the curvilinear polygon Ω shown in Fig.1via analytic patch maps F K M .

Example 1 Fig. 3 (Remark 4
left and middle) shows an example of an L-shaped domain with macro triangulation and suitable refinement patterns.For fixed L and increasing n, the meshes T L,n geo,σ are geometrically refined towards the vertices of Ω.These meshes are classical geometric meshes for elliptic problems in corner domains as introduced in[4,5] and discussed in [36, Sec.4.4.1].4.2 hp-FE approximation of singularly perturbed problems on geometric boundary layer meshes The principal result [7, Thm.4.1] on robust exponential convergence of hp-FEM for (4.1) reads as follows: Proposition 1 ( [7, Thm.4.1]) Let Ω ⊂ R 2 be a curvilinear polygon with J vertices as described in Section 1.1.Let A, c ≥ c 0 > 0, f satisfy (4.2).Denote by {T L,n geo,σ } L≥0,n≥L a sequence of geometric boundary layer meshes in the sense of Definition 2. Fix c 1 > 0.

Fig. 5
Fig. 5 Error convergence of the Extended hp-FEM and the sinc BK-FEM for the domain Ω 1 , depicted versus the polynomial degree p and N 1/2 ls , where N ls denote the number of linear systems that need to be solved.Solid lines correspond to Extended hp-FEM and dashed lines to the sinc BK-FEM.Results for s = 0.2, s = 0.4 and s = 0.8 are shown.

Fig. 6
Fig.6Error convergence of the Extended hp-FEM and the sinc BK-FEM for L-shaped domain Ω 2 , depicted versus the polynomial degree p and N 1/2 ls , the square root of the number of linear systems to be solved.

Fig. 7
Fig.7Error convergence of the Extended hp-FEM and the sinc BK-FEM for the slit domain Ω 3 , depicted versus the polynomial degree p and N 1/2 ls , the square root of the number of linear systems to be solved.