## 1. Introduction

In this work we consider finite element discretizations of the problem

(1.1) |

where is a bounded domain and is the integral fractional Laplacian of order ,

(1.2) |

The normalization constant makes the integral in (1.2), calculated in the principal value sense, coincide with the Fourier definition of . It is well understood that, even if the data is smooth (for example, if and ), then the unique solution to (1.1) develops an algebraic singularity near (cf. Example 2.2). This is in stark contrast with the classical Laplacian equation.

Nevertheless, in such a case one expects the solution to be locally smooth in , and thus the discretization error to be smaller in the interior of the domain. Our main result (Theorem 5.3) is a quantitative estimate of the fact that the finite element error is concentrated around .

The fractional Laplacian (1.2) is a nonlocal operator: computing requires the values of at points arbitrarily far away from . Nonlocality is also reflected in the variational formulation of (1.1): the natural space in which the problem is set is the zero-extension fractional Sobolev space , and the norm therein is not subadditive with respect to domain partitions. Furthermore, it is not possible to localize the inner product in , because functions with supports arbitrarily far away from each other may have nonzero -inner product. This is also in stark contrast with the local case (i.e., with the inner product in ), and makes the development of local estimates for such a nonlocal problem a more delicate matter, especially in the case of general shape–regular meshes. This is the main purpose of this paper.

In recent years, there has been significant progress in the numerical analysis and implementation of (1.1) and related fractional-order problems. Finite element discretizations provide naturally the best approximation in the energy norm. A priori convergence rates in the energy norm for approximations using piecewise linear basis functions on either quasi-uniform or graded meshes were derived in [2]; similar results, but regarding convergence in in case , were obtained in [7]. The use of adaptive schemes and a posteriori error estimators has been studied in [3, 20, 23, 32, 36]. A non-conforming discretization, based on a Dunford-Taylor representation was proposed and analyzed in [6]. We refer to [5, 8] for further discussion on these methods. In contrast, the analysis of finite difference schemes typically leads to error estimates in the -norm under regularity assumptions that cannot be guaranteed in general [16, 17, 26].

The rest of the paper is organized as follows. In Section 2, we review the fractional-order spaces and the regularity of solutions to (1.1) in either standard and weighted Sobolev spaces. In Section 3, we describe our finite element discretization, review basic energy based error estimates, and combine such estimates with Aubin-Nitsche techniques to derive novel convergence rates in -norm. In Section 4, we provide a proof of Caccioppoli estimate for the continuous problem. In Section 5, which is the central part of the paper, we combine Caccioppoli estimates and superapproximation techniques, to obtain interior error estimates with respect to -seminorms. At the end of this section we show some applications of our interior error estimates. In particular, we discuss the convergence rates of the finite element error in the interior of the domain with respect to smoothness of the domain and the right hand side in the case of quasi-uniform and graded meshes. The results are summarized in Tables 1 and 2. Finally, several numerical examples at the end of the paper illustrate the theoretical results from Section 5.

## 2. Variational formulation and regularity

In this section, we briefly discuss important features of fractional-order Sobolev spaces that are instrumental for our analysis. Furthermore, we consider regularity properties of the solution to (1.1) and review some negative results that lead to the use of certain weighted spaces, in which the weight compensates the singular behavior of the gradient of the solution near the boundary of the domain. Having regularity estimates in such weighted spaces at hand, we shall be able to increase the convergence rates by constructing a priori graded meshes.

### 2.1. Sobolev spaces

Sobolev spaces of order provide the natural setting for the variational formulation of (1.1). More precisely, we consider to be the set of -functions such that

(2.1) |

where is taken as in (1.2). Clearly, these are Hilbert spaces; we shall denote by the bilinear form that gives rise to the fractional-order seminorms, namely,

(2.2) |

For the variational formulation of (1.1), we need the zero-extension spaces

for which the form becomes an inner product. Moreover, if , then integration in (2.2) takes place in . We shall denote the -norm by , and remark that the -norm of is not needed because a Poincaré inequality holds in the zero-extension Sobolev spaces.

Fractional-order Sobolev spaces can be equivalently defined through interpolation of integer-order spaces; remarkably, if one suitably normalizes the standard

-functional, then the norm equivalence constants can be taken to be independent of [29, Lemma 3.15 and Theorem B.9]. Although the constant in (2.1) is fundamental in terms of continuity of Sobolev seminorms as , we shall omit it whenever is fixed. For simplicity of notation, throughout this paper we shall adopt the convention .Let denote the dual space to , and be their duality pairing. Because of (2.2) it follows that if then and

(2.3) |

This integration by parts formula motivates the following weak formulation of (1.1): given , find such that

(2.4) |

Because this formulation can be cast in the setting of the Lax-Milgram Theorem, existence and uniqueness of weak solutions, and stability of the solution map , are straightforward.

### 2.2. Sobolev regularity

The Lax-Milgram Theorem guarantees the well-posedness of (2.4) in if . A subsequent question is what additional regularity does inherit for smoother . For the sake of finite element analysis, here we shall focus on Sobolev regularity estimates.

By now it is well understood that for smooth domains and data , solutions to (1.1) develop an algebraic singular layer of the form (cf. for example [25])

(2.5) |

that limits the smoothness of solutions. Indeed, if is locally smooth in but behaves as (2.5), then one cannot guarantee that belongs to ; actually, in general (see Example 2.2).

We now quote a recent result [9], that characterizes regularity of solutions in terms of Besov norms. Its proof follows a technique introduced by Savaré [34], that consists in combining the classical Nirenberg difference quotient method with suitably localized translations and exploiting certain convexity properties.

###### Theorem 2.1 (regularity on Lipschitz domains).

Let be a bounded Lipschitz domain, and for some . Then, the solution to (1.1) belongs to the Besov space , where , with

(2.6) |

Consequently, by an elementary embedding, we deduce

(2.7) |

There are two conclusions to be drawn from the previous result. In first place, assuming the domain to be Lipschitz is optimal, in the sense that if was a domain then no further regularity could be inferred. Thus, reentrant corners play no role on the global regularity of solutions: the boundary behavior (2.5) dominates any point singularities that could originate from them; we refer to [22] for further discussion on this point. In second place, in general the smoothness of the right hand side cannot make solutions any smoother than . The expansion (2.5) holds in spite of the smoothness of near . We illustrate these two points with a well-known example [21].

###### Example 2.2 (limited regularity).

We also point out a limitation in the technique of proof in Theorem 2.1 from [9] that is related to the example above. Namely, in case and for some , solutions are expected to be smoother than just ; however, one cannot derive such higher regularity estimates from Theorem 2.1. For smooth domains (i.e., ), the following estimate holds [35]:

(2.9) |

### 2.3. Regularity in weighted Sobolev spaces

By developing a fractional analog of the Krylov boundary Harnack method, Ros-Oton and Serra [33] obtained a fine characterization of boundary behavior of solutions to (1.1) and derived Hölder regularity estimates. In order to exploit these estimates and apply them in a finite element analysis, reference [2] introduced certain weighted Sobolev spaces, where the weight is a power of the distance to . Let

Then, for and , we consider the norm

(2.10) |

and define and as the closures of and , respectively, with respect to the norm (2.10).

Next, for , with and , and , we consider

and the associated space

In analogy with the notation for their unweighted counterparts, we define zero-extension weighted Sobolev spaces by

(2.11) |

We have the following regularity estimate in the scale (2.11) [2, Proposition 3.12], [5, Formula (3.6)].

###### Theorem 2.3 (weighted Sobolev estimate).

Let be a bounded, Lipschitz domain satisfying the exterior ball condition, , for some , , and be the solution of (2.4). Then, it holds that and

###### Remark 2.4 (optimal parameters).

In finite element applications of Theorem 2.3, discussed in Section 3, we will design graded meshes with a grading dictated by . The optimal choice of parameters and depends on both the smoothness of the right hand side and the dimension of the space. We illustrate this now: let , , , and be sufficiently small, and choose and , to obtain the optimal regularity estimate

In contrast, if , we set to be any positive number and take as above to arrive at

###### Remark 2.5 (exterior ball condition).

Taking into account the results from [22], the exterior ball condition could be relaxed. Nevertheless, because the analysis of effects of reentrant corners is beyond our purposes, we leave such an assumption on .

## 3. Finite Element Discretization

We next consider finite element discretizations of (2.4) by using piecewise linear continuous functions. Let ; for , we let denote a triangulation of , i.e., is a partition of into simplices of diameter . We assume the family to be shape-regular, namely,

where and is the diameter of the largest ball contained in . As usual, the subindex denotes the element size, ; moreover, we take elements to be closed sets.

We shall also need a smooth mesh function , which is locally comparable with the element size. Note that shape-regularity yields (cf. [31, Lemma 5.1]), and thus

(3.1) |

Let be the set of interior vertices of , be its cardinality , and the standard piecewise linear Lagrangian basis, with associated to the node . With this notation, the set of discrete functions is

(3.2) |

It is clear that for all and therefore we have a conforming discretization.

### 3.1. Interpolation and inverse estimates

Fractional-order seminorms are not subadditive with respect to domain decompositions; therefore, some caution must be exercised when localizing them. With the goal of deriving interpolation estimates, we define the star (or patch) of a set by

Given , the star of is the first ring of and the star of is the second ring of . The star of the node is .

We have the following localization estimate [18, 19].

(3.3) |

This inequality shows that to estimate fractional seminorms over , it suffices to compute integrals over the set of patches plus local zero-order contributions. In addition, if these contributions have vanishing means over elements –as is often the case whenever is an interpolation error– a Poincaré inequality allows one to estimate them in terms of local -seminorms. Thus, one can prove the following local quasi-interpolation estimates (see, for example, [2, 10, 12]).

###### Proposition 3.1 (local interpolation estimates).

Let , , , and be a suitable quasi-interpolation operator. If , then

(3.4) |

where . Moreover, considering the weighted Sobolev scale (2.11), it holds that for all ,

(3.5) |

For the purpose of this paper, we shall make use of a variant of (3.3). Even though the fractional-order norms can be localized, it is clear that the -inner product of two arbitrary functions cannot: it suffices to consider two positive functions with supports sufficiently far from each other. The following observation is due to Faermann [19, Lemma 3.1]. Since we use it extensively, we reproduce it here for completeness.

###### Lemma 3.2 (symmetry).

For any and bounded, there holds

###### Proof.

Let

be the characteristic function of the set

. Using Fubini’s Theorem, we equivalently writewhere

Let’s fix . Since a.e. provided and otherwise, we deduce

This yields the desired identity. ∎

###### Proposition 3.3 (equivalent fractional inner product).

Let . Then, it holds that

###### Proof.

It suffices to write

and notice that

and

in view of Lemma 3.2 (symmetry) with , where and we recall that is the diameter of the largest ball contained in . This completes the proof. ∎

###### Remark 3.4 (fractional inner product on subdomains).

Proposition 3.3 is also valid for any subdomain , i.e.

Next, we write some inverse estimates that we shall use in what follows. By using standard scaling arguments, one can immediately derive the estimate

(3.6) |

Let be a fixed smooth function. We shall also need the following variant of (3.6) with , whose proof follows immediately because the space is finite dimensional:

(3.7) |

### 3.2. Energy-norm error estimates

The discrete counterpart of (2.4) reads: find such that

(3.8) |

Subtracting (3.8) from (2.4) we get Galerkin orthogonality

(3.9) |

The best approximation property

(3.10) |

follows immediately. Consequently, in view of the regularity estimates of discussed in Section 2, the only ingredient missing to derive convergence rates in the energy norm is some global interpolation estimate. Even though the bilinear form involves integration over , it is possible to prove that the corresponding energy norm is bounded in terms of fractional-order norms on by resorting to fractional Hardy inequalities (see [2]).

Therefore, for quasi-uniform meshes, if one can simply combine (3.3) and (3.4) with a fractional Hardy inequality [24, Theorem 1.4.4.4] to replace by [2, 10] and obtain

(3.11) |

In case , one cannot apply a fractional Hardy inequality. Instead, one may exploit the precise blow-up of the Hardy constant of as to deduce [2, §3.4], [10, Theorem 4.1]

(3.12) |

Alternatively, one could derive either (3.11) or (3.12) by simply interpolating standard global and estimates. However, if we aim to exploit Theorem 2.3 (weighted Sobolev estimate), then we require a suitable mesh refinement near the boundary of . For that purpose, following [24, Section 8.4] we now let the parameter represent the local mesh size in the interior of , and assume that, besides being shape-regular, the family is such that there is a number such that for every ,

(3.13) |

This construction yields a total number of degrees of freedom (see

[4, 10])(3.14) |

Thus, if the interior mesh size and the dimension of satisfy the optimal relation (up to logarithmic factors if ). As anticipated in Remark 2.4 (optimal parameters), the weight in Theorem 2.3 (weighted Sobolev estimate) needs to be related to the parameter , which satisfies (3.13). To do so, we combine (3.3) with either (3.5) or (3.4), depending on whether intersects or not, to find the relation for . If , it suffices to use a fractional Hardy inequality to replace by [2, 10] and obtain

(3.15) |

for all with a constant that depends on and . On the other hand, if , we choose , where is sufficiently small, and exploit the explicit blow-up of the Hardy constant of as , as we did earlier with (3.12), to derive the second estimate in (3.15). We point out that (3.15) does not follow by interpolation of global estimates.

We gather the energy error estimates for quasi-uniform and graded meshes in a single theorem.

###### Theorem 3.5 (global energy-norm convergence rates).

Let be a bounded Lipschitz domain, and denote the solution to (2.4) and denote by the solution of the discrete problem (3.8), computed over a mesh consisting of elements with maximum diameter . If , then we have

(3.16) |

where and if and zero otherwise. Additionally, if satisfies an exterior ball condition, let be such that

(3.17) |

Then, if , and the family satisfies (3.13) with as above, we have

(3.18) |

In terms of the number of degrees of freedom , the estimate above reads

(3.19) |

###### Proof.

If , we combine (3.10), (3.11) with Theorem 2.1 (regularity in Lipschitz domains) with to obtain

(3.20) |

In case , instead of (3.11) we use (3.12) with the same as in (2.7) to get

(3.21) |

Moreover, coupling (3.10), (3.15) and Theorem 2.3 (weighted Sobolev estimate) with and if and and if yields for

(3.22) |

and analogous estimates hold if . Upon taking , we end up with (3.16) and (3.18), as asserted. Inequality (3.19) follows by the choice of and (3.14). ∎

###### Remark 3.6 (exponents of logarithms).

###### Remark 3.7 (optimality).

The convergence rates derived in Theorem 3.5 are theoretically optimal for shape-regular elements. Nevertheless, because we deal with continuous piecewise linear basis functions, one would expect convergence rate with respect to . It is remarkable that such a rate can only be achieved if upon grading meshes according to (3.13). For dimensions , anisotropic meshes are required in order to obtain optimal convergence rates. This limitation stems from the algebraic singular layer (2.5) and becomes more apparent as increases, but comparison of (3.16) and (3.18) shows that in all cases graded meshes improve the convergence rates with respect to .

### 3.3. -norm error estimates

Upon invoking new regularity estimates for data , we now perform a standard Aubin-Nitsche duality argument to derive novel convergence rates in . We distinguish between quasi-uniform and graded meshes.

###### Proposition 3.8 (convergence rates in for quasi-uniform meshes).

Let be a bounded Lipschitz domain. If , then for all we have

(3.23) |

where and if and zero otherwise.

###### Proof.

Let be the error, and let be the solution to (2.4) with instead of the right hand side . Then, the Galerkin orthogonality (3.9) and the Cauchy-Schwarz inequality yield

where is a quasi-interpolation operator satisfying (3.11) if or (3.12) if . Combining these estimates with (2.7), we deduce

(3.24) |

This, in conjunction with the energy error estimates (3.20) and (3.21), implies

Finally, taking gives rise to (3.23). ∎

In Proposition 3.8, the assumption