FEM-BEM coupling in Fractional Diffusion

We derive and analyze a fully computable discrete scheme for fractional partial differential equations posed on the full space $\mathbb{R}^d$ . Based on a reformulation using the well-known Caffarelli-Silvestre extension, we study a modified variational formulation to obtain well-posedness of the discrete problem. Our scheme is obtained by combining a diagonalization procedure with a reformulation using boundary integral equations and a coupling of finite elements and boundary elements. For our discrete method we present a-priori estimates as well as numerical examples.


Introduction
In this work, we study stationary fractional partial differential equations posed on the full space R d with d = 2, 3 of the form with s ≥ 0, and β ∈ (0, 1).Fractional PDEs of this type are oftentimes used to model non-local effects in physics, finance or image processing, [BV16, SZB + 18].
Regarding the formal definition of non-integer powers L β of differential operators, there are various different descriptions in literature such as Fourier transformation, semigroup approaches, singular integrals or spectral calculus, see [LPG + 20].A distinct advantage of full-space formulations as in (1.1) is that all of these definitions are equivalent, [Kwa17], while there are significant differences in the definitions, if one restricts the problem to a bounded domain.
Nonetheless, there are usually no closed form solutions to these problems available and therefore numerical approximations are used.In order to derive a computable approximation, most numerical methods employ formulations on bounded domains, for which there is a fairly well developed literature.We mention the surveys [BBN + 18, LPG + 20] as well as finite element methods for the integral definition of the fractional Laplacian [AB17, ABH19, FKM22], for the spectral definition [NOS15,NOS16], and semigroup approaches [BP15,BLP19].We especially mention the very influential reformulation using the extension approach by Caffarelli and Silvestre [CS07] (see also [ST10] for a more general setting), which allows to use PDE techniques in the analysis.This approach paired with an hp-FEM approach in the extended direction has proven to be an effective strategy both for elliptic [MPSV18, BMN + 19, BMS23,FMMS22b,FMMS22a] as well as parabolic problems [NOS16,MR21].

Layout
The present paper is structured as follows: In the remainder of Section 1, we introduce our model problem as well as necessary notation and most notably, the Caffarelli-Silvestre extension problem.In Section 2, we formulate our main results: well-posedness of our formulation, the fully-discrete scheme using the diagonalization procedure together with the symmetric FEM-BEM coupling, and, finally, a best-approximation result.Section 3 provides the proofs for the well-posedness and the diagonalization procedure and, most notably, a Poincaré type estimate.Section 4 contains the proofs for the a-priori analysis of the fully discrete formulation using hpfinite elements in the extended variable, which builds upon the regularity and decay properties of [FR23].Finally, Section 5 presents some numerical examples that validate the proposed method.

Notations
Throughout the text we use the symbol a b meaning that a ≤ Cb with a generic constant C > 0 that is independent of any crucial quantities in the analysis.Moreover, we write ≃ to indicate that both estimates and hold.We employ classical integer order Sobolev spaces H k (Ω) on (bounded) Lipschitz domains Ω and the fractional Sobolev spaces H t (R d ) for t ∈ R defined, e.g., via Fourier transformation.We also need Sobolev spaces on the boundary Γ := ∂Ω of a bounded Lipschitz domain Ω ⊂ R d , denoted by H t (Γ) with t ∈ [−1, 1].One way to properly define them is by using local charts, see [SS11] for details.

Assumptions on the model problem
Let d = 2, 3. We consider (1.1) and, for functions u ∈ L 2 (R d ), define the self adjoint operator L β using spectral calculus where E is the spectral measure of L and σ(L) is the spectrum of L. Using standard techniques this definition can be extended to tempered distributions.
For the data, we assume f ∈ L 2 (Ω) and A ∈ L ∞ (R d , R d×d ) is pointwise symmetric and positive definite in the sense that there exists A 0 > 0 such that (A(x)y, y) 2 ≥ A 0 y 2 2 ∀y ∈ R d .
In order to avoid several additional difficulties due to decay conditions at infinity, we assume s ≥ σ 0 > 0 for the case d = 2.
Additionally, in order to be able to apply FEM-BEM coupling techniques, we make the following assumptions on the coefficients in the model problem: There exists a bounded Lipschitz domain This is not the most general setting where our techniques can be applied.For example, also a lowest order term could be included in the definition of L in (1.1); see also Remark 2.10.

Degenerate elliptic extension
In the same way as in our previous work [FR23], we use a reformulation of the fractional PDE as the Dirichlet-to-Neumann mapping for a degenerate elliptic PDE in a half space in R d+1 , the so called Caffarelli-Silvestre extension, [CS07,ST10].We recall the definition of weighted Sobolev spaces used for the extension.For any bounded open subset D ⊂ R d × R, we define L 2 (y α , D) as the space of square integrable functions with respect to the weight y α and the Sobolev space We also employ the spaces L 2 (y α , (0, Y)) and H 1 (y α , (0, Y)) for Y ∈ (0, ∞] defined in an analogous way by omitting the x-integration.For unbounded sets D, we additionally use the weight to take care of the behaviour at infinity.In this case, we define the space H 1 ρ (y α , D) as the space of all square integrable functions U (with respect to the weight function y α ρ −2 ) such that the norm Moreover, we also employ spaces acting only in x.Using the weight 2) by omitting the y-integration.
For functions U in H 1 ρ (y α , R d × R + ), one can give meaning to their trace at y = 0, which we denote by tr 0 U .In fact, by [KM19, Lemma 3.8] and [FR23,Lem. 3.1] we have the trace estimates Then, the extension problem reads as: find where d+1) .By [ST10], the solution to (1.1) is then given by u = tr 0 U .
For the domain Ω with boundary Γ := ∂Ω, we also introduce the usual trace operators γ − Γ (denoting the trace coming from the interior of Ω) and γ + Γ (denoting the trace coming from R d \Ω) and correspondingly the normal derivative operators ∂ ± ν,Γ (see [SS11] for details).The normal vector ν is always assumed to face out of Ω.With theses operators, the jumps across Γ are defined as (1.5) We will apply these operators for functions in , where they are to be understood pointwise with respect to y.This is equivalent to taking the trace and normal derivative along the lateral boundary Γ × R + .

Variational formulation
The weak formulation of (1.4) in In order to also include our discretization scheme, we work in a slightly expanded variational form, inspired by [LS09,Say09].In short, one can formulate an equivalent problem for the solution inside Ω and a function U ⋆ on R d defined in a modified Hilbert space.Definition 2.1.Fix Y ∈ (0, ∞].We consider the space equipped with the norm We note that, by definition, the additional condition of tr 0 U ⋆ being in L 2 (R d ) is only needed for s = 0 as in this case the norm in H Y contains said L 2 -term, which has to be finite.
For f ∈ L 2 (Ω), the weak formulation is given as the problem of finding U ∈ H ∞ such that Problems (1.4) and (2.3) are connected as follows: In order to obtain a computable formulation, we start by cutting the problem from the infinite cylinder R d × R + to a finite cylinder in the y-direction.To do so, we fix a parameter Y > 0 to be chosen later and introduce the truncated bilinear forms The "big" bilinear form is then given by and the cutoff problem reads as: Find (2.4) By the following theorem, we obtain well-posedness of the weak formulation of both variational formulations.
The proof of the theorem is given in Section 3 and, in fact, reduces to the application of suitable Poincaré inequalities.

The discrete scheme
In this section, we describe our discrete scheme to approximate solutions to the truncated variational formulation (2.4).The main idea is to employ a tensor product structure for the approximation by using the diagonalization procedure described in [BMN + 19].
Let V y h be an arbitrary finite dimensional subspace of L 2 (y α , (0, Y)) of dimension N y + 1.Following the ideas of [BMN + 19], we chose an orthonormal basis (ϕ j ) It is easy to see that for s = 0, the constant function is an eigenfunction corresponding to the eigenvalue µ = 0.Moreover, the assumption s > 0 for d = 2 guarantees that there is not a zero eigenvalue.If zero is an eigenvalue (for d = 3), we assume that the eigenvalues are ordered such that µ 0 = 0.
We now give a formal definition of our (semi-)discrete subspace of H Y , which has tensor product structure with respect to the variables x, y.
Definition 2.3.Let V x h ⊂ H 1 (Ω) and V λ h ⊂ H −1/2 (Γ) be finite dimensional spaces and Y ∈ (0, ∞).Additionally, assume that 1 ∈ V λ h .We introduce the closed subspace Then, the semi-discrete problem reads as: (2.7) Using the orthogonal basis for the y-direction, we can actually diagonalize some of the bilinear forms and obtain an equivalent sequence of scalar problems.In fact, functions (U Ω , U ⋆ ) ∈ H h,Y solve (2.7), if and only if they can be written as , where the functions u j,Ω , u j,⋆ solve We refer to Lemma 3.3 for a proof of this statement.
The equation for u j,⋆ is still posed on an unbounded domain.We will replace this with boundary integral equations.Therefore, given µ ∈ C with Re(µ) ≥ 0, we introduce where H (1) 0 denotes the first kind Hankel function of order 0. The single-layer and double-layer potential are then defined as and the corresponding boundary integral operators are given by (2.10) (2.11) We then have the following result, giving a computable approximation of (1.1) that only relies on well-known operators.
The problems (2.12) are standard FEM-BEM coupling problems for what is often called the modified Helmholtz or Yukawa equation.As such, existence and uniqueness of solutions (u j , λ j ) ∈ . Consequently, we also obtain well-posedness of the semi-discrete formulation (2.7) as we have constructed a solution in H h,Y .Uniqueness follows from coercivity of the bilinear form.
, the space V λ h contains the constant functions, and either d = 3 or s > 0.Then, the truncated problem (2.7) has a unique solution Remark 2.6.Due to the construction in Theorem 2.4, we mention that our discrete approximation can very easily be computed with the use of existing FEM/BEM libraries.We refer to Section 5 for a description of the implementation used in the numerical examples therein.

A-priori convergence estimates
In the extended variable y, we employ a hp-FEM discretization.
Let Y > 0 and T y be a geometric grid on (0, Y) with mesh grading factor σ, L-refinement layers towards 0, and M = ⌊ln(Y)/ ln(σ)⌋ levels of growth towards Y.More precisely, we define the grid points as we denote the space of continuous, piecewise polynomials of degree up to p.
The following proposition provides a best-approximation estimate for the hp-semi-discretization in y.We note that in contrast to [BMN + 19], which exploits a known closed form representation of the solution, we only have algebraic convergence of the truncated solution rather than exponential convergence.However, choosing the truncation parameter large enough, the hpsemi-discretization in y allows to recuperate any algebraic convergence rates of the discretization in x.
Remark 2.8.A possible choice for the spatial discretization would be V x h := S 1,1 (T x ), i.e., continuous, piecewise linear polynomials on some (quasi-uniform) mesh T x of Ω.For the operator π Ω one could take the Scott-Zhang projection mapping onto S 1,1 (T x ), see [SZ90].In addition to the required L 2 (Ω)-and H 1 (Ω)-stabilities, the operator has first order approximation properties in H 1 (Ω), provided the input function is sufficiently regular.
Using first-order approximation properties of the Scott-Zhang projection together with bestapproximation of the BEM part (which converges of order h 3/2 assuming sufficient regularity, see [SS11]) and correct choice of the cut-off parameter Y and polynomial degree p, the bestapproximation estimate for the semi-discretization in Theorem 2.7 directly gives first order convergence in h.
Corollary 2.9.Let the assumptions of Theorem 2.7 hold.Assume A ∈ C 1 (R d , R d×d ) and f ∈ H 1 (Ω) and assume Ω has piecewise smooth boundary.Choose V x h := S 1,1 (T x ) with a quasiuniform mesh T x of Ω of maximal mesh-width h and take π Ω to be the Scott-Zhang projection.Let V λ h := S 0,0 (T Γ ) be the space of piecewise constants on the trace mesh T Γ of T x .Moreover, choose p = −c κ,µ,ε ln h with a sufficiently large constant c κ,µ,ε depending only on κ, µ and ε, and Remark 2.10.We note that our main results are valid for more general fractional PDEs as well.Using the same techniques, one obtains the statements also for: 2. operators containing lower order terms, i.e.,

Well-posedness and FEM-BEM formulation
In this section, we provide the proofs of Theorem 2.2 and Theorem 2.4.

Poincaré inequalities
We now show the well-posedness of our variational formulations.The main ingredient is a Poincaré type estimate, which uses the following compactness result.
Proof.We can cover the Lipschitz domain D by a finite number of Lipschitz domains D 1 , . . ., D m which are starshaped with respect to a ball, see for example [Maz11, Sect.1.1.9,Lemma 1].Thus, without loss of generality we may assume that D is starshaped with respect to a ball.With c n := D u n we compute where we used the Poincaré estimate of [NOS15, Corollary 4.4] and the assumed weak convergence.
provided the right-hand side is finite.
Proof.The estimates follow from techniques employed in [AGG94, Theorem 3.3] using a proof by contradiction.In the first step, we show (3.1) (which essentially is covered by [AGG94, Theorem 3.3], we only account for the additional weight y α ) and (3.2) for functions vanishing inside a ball containing the origin.Finally, using a compactness argument this assumption is removed in the second step.
Step 1: First, assume that U ≡ 0 on a sufficiently large (half) ball B R (0) ⊂ R d+1 and has compact support.We employ spherical coordinates in R d × R + , chosen such that y = r cos(ϕ) and collect the remaining d − 1 angles into φ.Using where we denoted by J(ϕ, φ) the angular components of the Jacobian in the transformation theorem.Integration by parts in r and using the assumed support properties of U gives Transforming back to (x, y)-variables and using r µ ≤ ρ µ for µ ≥ 0 this gives the desired bound.By density, we can remove the requirement of compact support of U .
Step 2: In order to get rid of the requirement that U vanishes on the ball B R (0), we use a compactness argument.To keep the notation succinct we set µ = 0 in the following, the general case µ > 0 can be done with the exact same arguments.Assume that (3.1) does not hold, i.e., there exists a sequence Since U n is a bounded sequence in the Hilbert space H 1 ρ (y α , R d \ Γ × (0, ∞)), there exists a weakly convergent subsequence (also denoted by U n ) and we denote the weak limit by U .Since the seminorm is lower semicontinuous, we get |U | H 1 ρ (y α ,R d \Γ×(0,∞)) = 0.A simple calculation (using polar coordinates, similar to the estimate above) shows that -as we are in half-space in R d+1 with d + 1 > 2 and Γ γU ds x = 0 -the space H 1 ρ (y α , R d \ Γ × (0, Y)) does not contain piecewise constant functions except for 0, which means that U = 0. We now show strong convergence of the sequence to U = 0. To that end, fix a ball B R := B R (0) ⊂ R d+1 with sufficiently large R such that Ω × {0} ⊂ B R and consider a smooth cutoff function ψ : R d+1 → R such that ψ ≡ 1 on B R and ψ ≡ 0 on B 2R .We thus decompose U n as From the compactness result of Lemma 3.
) on all bounded half balls B R with sufficiently large R. Since U 2 n vanishes on B R , we can apply step 1 of the proof to determine: Overall, we get that U n → 0 in H 1 ρ (y α , R d \Γ×(0, ∞)), which is a contradiction to the assumption U n H 1 ρ (y α ,R d \Γ×(0,∞)) = 1 for all n ∈ N.
Step 3: Estimate (3.2) for the case d = 3 follows directly from multiplying a full-space Poincaréinequality (see for example [AGG94, Theorem 3.3] for µ = 0 and a similar calculation to step 1 for 0 < µ ≤ 2 with polar coordinates only in x) applied only in x with y α and integrating over (0, Y).
Step 4: It remains to show (3.2) for d = 2, which was also shown in [FR23, Lem.3.2].We write U (x, y) = U (x, 0) + Since Y 0 y α ρ µ−2 dy 1 for sufficiently small µ < µ 0 with µ 0 depending only on α, the first term on the left-hand side can be bounded by For the second term, we employ a weighted Hardy-inequality, see e.g.[Muc72], to obtain which shows the claimed inequality.
We can now look at the well-posedness of our discrete problem.
Proof of Theorem 2.2.Let Y ∈ (0, ∞] and (U Y Ω , U Y ⋆ ) ∈ H Y .On the interior domain Ω, we integrate a standard Poincaré-like estimate to obtain By the conditions on γ − U Y ⋆ and γU Y ⋆ , we observe that the function ⋆ , in R d \Ω has a jump across ∂Ω with vanishing integral mean.Applying Lemma 3.2 to U Y , we get with (3.3) that which shows coercivity.
In order to bound the right-hand side in (2.4), we can directly use the definition of the H Y -norm together with supp f ⊂ Ω for s > 0 to obtain For Y = ∞ and s = 0, which implies d = 3 by assumption, the trace estimate (1.3) gives For the case Y < ∞ and s = 0, we use a cut-off function χ satisfying χ ≡ 1 on (0, Y/2), supp χ ⊂ (0, Y) and ∇χ L ∞ (R + ) Y −1 .As Ω is bounded, this gives with the trace estimate [KM19, Lem.3.7] which finishes the proof.

Diagonalization
We now apply the diagonalization procedure of [BMN + 19] to show that solutions of (2.7) can be written as in (2.8).We recall that (ϕ j ) Ny j=0 is the orthonormal basis of eigenfunctions from (2.5) with corresponding eigenvalues µ j .
), if and only if they can be written as where • ∈ {Ω, ⋆} and Proof.At first, we show unique solvability of (3.4).For that, we consider the weak formulation of (3.4) given by The equivalence between the weak form and the strong form follows from standard arguments, and we refer to [LS09, Sect.7].Coercivity of the weak formulation in H 1 (Ω) × H 1 (R d \Γ) is clear for µ j > 0 as A is positive definite.For µ j = 0, one can employ Poincaré estimates on Ω and R d (with weights) to obtain coercivity in H 1 (Ω) × H 1 ρx (R d \Γ).Therefore, for each j, a unique solution (u j,Ω , u j,⋆ ) ∈

7). By construction we have U Y
h ∈ H h,Y .We next look at the weak formulation of (2.7).First we focus on the ⋆-contribution.Taking V Y ⋆ = v j,⋆ (x)ϕ j (y) with arbitrary v j,⋆ ∈ V x h as test function, we compute For the interior contribution, the same diagonalization procedure gives for Summing up, and using the weak form (3.5) we get that By density we can extend this equality to all test functions V Y h in the space H h,Y and get (2.7).Since the bilinear form B Y (•, •) is coercive we get that the function U Y thus constructed is the only solution to (2.7), which establishes the stated equivalence.
Proof of Theorem 2.4.The statement follows from Lemma 3.3 and [LS09, Section 7], as defining u j,⋆ (x) := V (µ j )λ j (x) − K(µ j )γ − Γ u j,Ω (x) and plugging that into (3.4)gives the stated equations using classical properties of the layer potentials.By definition and decay of the layer potentials, we have that u j,⋆ ∈ H 1 (R d \Γ) for all j ∈ N 0 such that µ j = 0, which gives u j,⋆ ∈ H 1 ρx (R d \Γ) as well.If s > 0, no zero eigenvalue is possible.This matches with the the requirement s tr 0 U ⋆ ∈ L 2 (R d ) in the definition of the space H Y .The case s = 0 is only allowed for d = 3.Here, tr 0 U ⋆ is not required to be in L 2 (R d ), which would not hold.However, in this case the decay property of the layer potentials for the Poisson equation, see e.g.[SS11], give u j,⋆ ∈ H 1 ρx (R d \Γ).In short, we have that our constructed solution is in the semi-discrete space H h,Y .Notably, (3.4) is just the "non-standard transmission problem" corresponding to the standard symmetric FEM-BEM coupling given by Theorem 2.4.

Error analysis
The key to the error analysis are the decay and regularity properties shown in [FR23].In order to make the present paper more accessible, we summarize the key results of [FR23] in the following.

Decay and regularity
The solution to the truncated problem is in fact a weak solution to a Neumann problem.Thus, in this section, we consider solutions U Y to the following truncated problem: Then, the truncation error can be controlled via the following proposition.
The goal in the following is to employ hp-FEM in the extended variable y.Therefore, weighted analytic regularity estimates are the key to the a-priori analysis.
Denoting by L 2 (y α , (0, Y); X) the Bochner spaces of square integrable functions (with respect to the weight y α ) and values in the Banach space X, the regularity results of the previous Proposition can be captured by the solution being in some countably normed space.For constants C, K > 0, we introduce L 2 (y α ,(0,Y);X) Corollary 4.3.Fix Y ∈ (0, ∞] and let U Y solve (4.1).Then, there are constants C, K > 0 such that there holds

Fully discrete analysis
In order to derive error bounds, we employ the reformulation in (2.2) together with the already established decay bounds for the truncation in Y.We will need two quasi-interpolation operators -one for the x-variables and one for the ydirection.Their construction and properties are the subject of the next two lemmas.
(iv) There hold the approximation properties: Proof.We note that a very similar operator is introduced in [MR17, Lemma 4.3].We define: By construction, we have (i), since u h ⋆ = 0 in the interior.Since γ + Ev = γ − v due to the extension property, we get (ii) by The stability estimates follow from the stability of the extension operator and the assumed stabilities of π Ω as The H 1 -estimates follows analogously.
Let ε > 0 be given by Proposition 4.2.Then, choosing L = p, there exists an operator Π y : •) ∈ S p,1 (T y ) for almost all x ∈ R d , and such that the following estimate holds: The constants C, κ > 0 are independent of p, Y.
Proof.We use the hp-interpolation operator from [BMN + 19, Sec.5.5.1] for Π y .This operator is constructed on a geometric mesh in an element-by-element way.On the first element a linear interpolant in σ L /2 and σ L is used, while the remaining elements are mapped to the reference element, on which a polynomial approximation operator that has exponential convergence properties (in the polynomial degree) for analytic functions is used.
For the y-derivatives the situation is a bit more involved, as the same argument can not be made as Π y and ∂ y do not commute.[BMN + 19, Lem.11(ii)] gives an exponentially convergent error bound for the y-derivative provided ).However, in our setting, the requirement , which is enough to regain the exponential estimate as seen in the following.On the first element (0, σ L ) ∈ T y , the definition of the piecewise linear interpolation gives which is nothing else than the L 2 -orthogonal projection of ∂ y v on (σ L /2, σ L ).By choice of the Babȗska-Szabó operator and denoting by Π L 2 p−1 the mapped L 2 -projection onto an element in T y , we have due to the commutator property and the preceding discussion The error estimate for the y-derivative follows from scaling arguments.More precisely, we decompose where K i = (x i , x i+1 ).Using a Hardy inequality, one obtains a bound for the approximation error on the first element using second derivatives only; see [BMN + 19, Lem.15].Together with a scaling argument this leads to . By Corollary 4.3 we can bound the right-hand side.For the remaining elements, we employ a scaling argument from [AM15, Thm.3.13].Denoting by h K i the diameter of K i , we infer y ∼ h K i on K i for i > 0. For any univariate function v satisfying y ℓ−ε v (ℓ+1) L 2 (y α ,(0,Y)) < CK ℓ ℓ! for all ℓ ∈ N 0 there holds where v is the pull-back of v to the reference element.The exponential approximation properties of the Babȗska-Szabó polynomial approximation operator then provides , we can employ (4.4) for v(y) = U (y, •) and square integrate over R d , noting that (4.3) holds due to Corollary 4.3.Summing over i and using i h 2ε shows the claimed estimate.Finally, to show that the operator does indeed map to H 1 ρ (y α , R d × (0, Y)), we note that by the previous considerations we have ∂ y Π y U ∈ L 2 (y α , R d × (0, Y)) as well as Π y U (•, y) = U (•, y) ∈ L 2 ρx (R d ) for certain values y ∈ (0, Y) where it is interpolatory.By the fundamental theorem of calculus, this is sufficient to show that Π y U ∈ L 2 ρ (y α , R d × (0, Y)).
We can now define an interpolation operator acting on both x and y in a tensor product fashion.
In order to keep notation compact, we write with the operators Π x from Lemma 4.4 and Π y from Lemma 4.5.Assume that the operator π Ω in the definition of Π x is both L 2 -and H 1 -stable.Then, the following approximation estimate holds Proof.By the Poincaré inequality (3.2) and the trace inequality (1.3), we only have to estimate the gradient norms.We start with the x-derivatives.Employing the H 1 -stability and approximation properties of Π x from Lemma 4.4 (iii) and (iv) gives Employing again Poincaré inequalities, we can reduce the right-hand side to norms of derivatives only.For the y-derivative, we proceed similarly using the L 2 -stability and approximation which finishes the proof.
In order to obtain a best-approximation estimate for the semi-discretization, we observe that the difference U − U h satisfies some form of Galerkin orthogonality.
Proof.Compared to "standard" Galerkin orthogonality, we observe that V Y h is not an admissible test function in (2.4) due to the weak condition of γ -setting (i.e.working with global functions instead of pairs), the test function However, if we use the pointwise equation (4.1) and integrate back by parts, we get that we can subtract such a term from the right-hand side without changing the equality, which shows the stated Galerkin orthogonality.
Finally, we are in position to show our main result, Theorem 2.7, by combining the decay estimate with the previous two lemmas.
Proof of Theorem 2.7.We start with the triangle inequality For the first term, we use the decay properties of Proposition 4.1 to obtain For the second term, we employ the coercivity of Theorem 2.2, the Galerkin orthogonality of Lemma 4.7, U Y ⋆ | Ω = 0, and a trace inequality for Ω, which gives for arbitrary λ h ∈ V λ h and Taking ε sufficiently small and absorbing the first term in the left-hand side gives with the operator Π of Lemma 4.6.Then, Lemma 4.6 together with the approximation properties of the hpinterpolation in Y gives Combining all estimates gives the stated result.
Finally, we present the proof of Corollary 2.9 that gives first order convergence for a specific choice of discrete spaces.
Proof of Corollary 2.9.Employing [FR23, Pro.2.8] -which with the same techniques also holds for Y < ∞ and a constant independent of Y -together with the assumptions on Ω, A, and f , we obtain control of second order x-derivatives of U Y .As λ := ∂ + ν U Y ⋆ , it is piecewise smooth, depending on the regularity of U Y ⋆ .For m = 0, 1, denoting by π L 2 the L 2 -projection onto S 0,0 (T Γ ) and using a trace estimate, it holds that Interpolating between m = 0 and m = 1 gives Multiplying with y α and integrating with respect to y then controls the second term on the righthand side of Theorem 2.7.For the first term, the approximation properties of the Scott-Zhang projection together with control of the second order x-derivatives gives first order convergence in h.Finally, the last two terms in Theorem 2.7 can also be bounded by Ch by choice of Y and p.

Numerics
In this section, we present two numerical examples to underline the a-priori estimates of Theorem 2.7 and Corollary 2.9.As previously already mentioned, a nice feature of our numerical scheme is that software packages developed for integer order differential operators can be employed directly.As such, we implement our method based on a coupling of the libraries NGSolve ([Sch21], for the FEM-part) and Bempp-cl ([BS21], for the BEM-part) libraries.
In order to validate our numerical method, we consider the case s = 0 and the standard fractional Laplacian, i.e., A = I.In this case a representation formula is available from [CS07].In fact, the fundamental solution for the fractional Laplacian is given by We then calculate u(x) at random sampling points x j using spherical coordinates and Gauss-Jacobi numerical integration to deal with the singularity at r = |x − x j | = 0 as well as standard Gauss-Quadrature for the other coordinate directions.
In order to compute the energy error, we compute the energy differences.For standard FEM with bilinear form a(•, •) and right-hand side f (•), it is well known that one can compute the energy error by the identity u − u h 2 E = a(u, u) − a(u h , u h ) = f (u) − f (u h ).Due to the more complicated form of our method, most notably the presence of the cutoff error, such an identity does not hold exactly.Nevertheless, we expect the following identity to hold approximately We now further replace the unknown value (f, tr 0 U ) L 2 (Ω) by the extrapolation from (f, tr 0 U Y h ) L 2 (Ω) for different refinements using Aitken's ∆ 2 -method.This will be our approximation of the true energy error.For the L 2 -error, we use the approximation U Y h on the finest grid as our standin for the exact solution compare it to the other approximations by computing the L 2 -difference of the traces at y = 0 using Gauss quadrature.For the geometry, we used the unit cube Ω := [−1, 1] 3 .In the bounded domain Ω, we use piecewise linear Lagrangian finite elements on a quasi-uniform mesh of maximal mesh width h.In Figure 5.1, we study the convergence of the proposed fully discrete method as we reduce the mesh size.In order to reduce all the error contributions, we choose the cutoff point Y = h − 2 1+α which gives O(h) for the cutoff error in Proposition 4.1.Since the convergence with respect to the polynomial degree is exponential (but with unknown explicit rate), we use p := round(2m log(m + 1)) where m is the number of uniform h-refinements.This gives a decrease of the y-discretization error which is faster than O(h).Overall we expect the energy error to behave like O(h) by Corollary 2.9.For the pointwise and L 2 -errors we did not establish a rigorous theory in this work.Nonetheless, Figure 5.1 shows convergence rates for these error measures of roughly order O(h 2 ).As a second numerical example, we consider as the domain Ω the unit sphere in R 3 .Instead of using the standard Laplacian with constant coefficients, we consider the following diffusion parameter and right-hand side: (with the slight abuse of notation of making A(x) scalar valued).Since the coefficients are globally continuous, and we are working with lowest order elements, by Corollary 2.9 we expect to obtain first order convergence.Figure 5.2 supports the theoretical results.Since in this case the fundamental solution is not available, we can not compute the pointwise error, but looking at the extrapolated energy and L 2 -errors we get the optimal rates.