{'fraction_non_alphanumeric': 0.047597939885173926, 'fraction_numerical': 0.01471209050996285, 'mean_word_length': 4.195416164053076, 'pattern_counts': {'":': 0, '<': 8, '<?xml version=': 0, '>': 18, 'https://': 0, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 7, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
Discretizations of the spectral fractional Laplacian on general domains with Dirichlet, Neumann, and Robin boundary conditions

In this work, we propose novel discretizations of the spectral fractional Laplacian on bounded domains based on the integral formulation of the operator via the heat-semigroup formalism. Specifically, we combine suitable quadrature formulas of the integral with a finite element method for the approximation of the solution of the corresponding heat equation. We derive two families of discretizations with order of convergence depending on the regularity of the domain and the function on which the spectral fractional Laplacian is acting. Our method does not require the computation of the eigenpairs of the Laplacian on the considered domain, can be implemented on possibly irregular bounded domains, and can naturally handle different types of boundary constraints. Various numerical simulations are provided to illustrate performance of the proposed method and support our theoretical results.

Preliminary results of fractional models are promising and encourage the exploration of their use in more realistic settings, so that interesting mathematical models could phase into physically grounded equations, able to predict phenomena observed in experiments. The present study is motivated by providing accessible and effective tools for real world models involving fractional calculus operators.
In order to fix the physical foundation of fractional models, bounded domains and boundary conditions come into play [30].
Typically, experiments are set up so that no effects are expected from outside the domain where the sample is contained and the local boundary conditions should fully determine the evolution of the process under study. For example, the observable under investigation can be zero at the boundaries, mathematically described by homogeneous Dirichlet boundary conditions. However, depending on the phenomena and the domain observed, other boundary conditions involving the specification of a nonzero value of the observable or a particular expression of its flux at the boundaries are required. This last example is the case of an insulated domain where the imposition of reflecting boundaries (i.e., homogeneous Neumann boundary conditions) ensures mass conservation. When the observed process takes place in a domain characterized by high spatial complexity and/or structural heterogeneity we cannot ignore what happens at the interface between different structural compartments nor the conditions imposed at the boundary of the considered domain; hence the need for a consistent definition of nonlocal operators with boundary conditions.
With this in mind, in the present paper we study a family of nonlocal operators defined in a consistent manner with the bounded domain on which they act and the corresponding local boundary conditions.
Since this issue concerns the nonlocality in space, we consider the fractional Laplacian as a prototype operator. For this nonlocal operator many definitions exist (e.g., Fourier multiplier [54], Riesz potential operator [36], hypersingular integral [46], heatsemigroup [56], extension problem [12]) and these definitions are proved to be equivalent for sufficiently smooth functions in R n , n ≥ 1, [35]. However, when the nonlocal operator is considered on bounded domains, this equivalence no longer holds. In the literature, one can find many different approaches for the definition or restriction of the fractional Laplacian in the bounded domain case and often, for each of these definitions, various schemes are proposed for the numerical computation of the considered operator (see e.g., [31,14,18,20,7,5,23,15,60,22,24,1]).
The work presented in this manuscript is based on the spectral definition of the fractional Laplacian on an open set Ω ⊂ R n [11,57,59]. Let L be a linear non-negative second order differential operator on L 2 (Ω), self-adjoint and with a compact inverse. Then there is an orthonormal basis of L 2 (Ω) made of eigenfunctions {ϕ m } m∈N of L, i.e., Lϕ m = λ m ϕ m , such that λ m ≥ 0, and λ m → ∞. Then, Lu = m λ m u, ϕ m ϕ m , and fractional powers of the operator, namely L s , with s ∈ (0, 1), are defined as for all functions u ∈ L 2 (Ω) such that m λ 2s m | u, ϕ m | 2 < ∞ ( ·, · denotes the inner product in L 2 (Ω)). The choice of functional spaces on which L satisfies the aforementioned properties and the particular expression of eigenvalues λ j and eigenfunctions ϕ j depend not only on L but also on Ω and on the boundary conditions considered for this domain. Therefore, here and throughout the rest of the manuscript, we use the subscript B to highlight that a given operator, in the bounded setting, is considered together with some boundary conditions. For example, we denote the Laplacian on Ω by (−∆ B ) rather than simply (−∆).
When L = (−∆ B ), from (1.1) we obtain the definition of the spectral fractional Laplacian (−∆ B ) s on Ω and the notation is consistent with the one used for the Dirichlet and Neumann spectral fractional Laplacians in references [59], [58], and [29]. The fact that different boundary conditions modify the definition of the nonlocal operator on the bounded domain Ω is in complete agreement with our initial physical considerations and motivates our choice of working with the spectral definition rather than any other definition of fractional operator available in the literature.
Regarding the numerical approximation of fractional equations on bounded domains Ω involving the spectral fractional Laplacian, we mention here two strategies. On the one hand, the matrix transfer technique, initially proposed by Ilić et al. [33,34], computes a finite difference approximation of the spectral fractional Laplacian via the direct approximation of the eigenvalues and eigenvectors of the standard discrete Laplacian coupled to some given boundary conditions. This method can be easily adopted on simple geometries but, in the case of general domains, only estimates of the eigenpairs are available, and computations require the solution of a very large number of eigenvalue problems. An efficient method for approximating the eigendecomposition of the Laplacian on more general domains has been recently proposed in [55]. On the other hand, for general Lipschitz domains Ω ⊂ R n , n > 1, Nochetto et al. [47] developed a finite element (FE) technique to compute the solution to the fractional Poisson equation involving the spectral fractional Laplacian with homogeneous Dirichlet boundary conditions and interpreted in light of the Caffarelli-Silvestre extension method [12]. See [18] for a similar approach in the context of finite differences. Bonito and Pasciak in [7] proposed an integral formulation for fractional powers of the inverse of the Laplacian to solve the fractional Poisson equation.
Differently from the approaches mentioned above, we exploit the fact that definition (1.1) is equivalent to the one based on the heat-semigroup [56]. In particular, for the spectral fractional Laplacian we have where e t∆ B u represents the solution of the classical heat equation on the bounded domain Ω with initial condition u and boundary conditions B.
The connection between the nonlocal operator and the local heat equation provided by equation (1.2) allows us to develop an accessible and effective method for the numerical computation of the considered nonlocal operator on general bounded domains and with various boundary constraints. Specifically, in the present paper we compute the value of the operator (−∆ B ) s on a given function by combining a FE approach for the solution of the heat equation and a quadrature rule for the integral in (1.2) with respect to the measure dµ(t) = dt t 1+s . We provide a suitable discretization for the considered fractional operator and prove convergence of the numerical approximation. The cases of Dirichlet, Neumann and Robin boundary conditions are studied in details and several numerical results are provided. The numerical method is also used to compute the numerical approximation to a fractional porous medium type equation on a bounded domain: the sharp boundary behavior predicted in [5] is recovered.
Fractional powers of the Laplacian with homogeneous boundary conditions are our main focus. However, under suitable assumptions, we propose a possible strategy to deal with the case of non-homogeneous boundary conditions, and we outline how our approach can be applied to compute fractional powers of more general second order elliptic operator L and how it can be adapted to deal with the case of more general boundary conditions.
Organization of the paper. In section 2, we introduce our framework, some definitions, the concepts of the FE theory, and provide the general form of the quadrature approximation that will be used in the rest of this work. After stating our main assumptions, in section 3 we provide our first discretization and the corresponding proof of convergence. Under suitable assumptions, a similar approach is used to derive a higher order discretization of the nonlocal operator, as shown in detail in section 4. In section 5 we study the case of non-homogeneous boundary conditions. We define a fractional-type operator corresponding to the Laplacian coupled with non-homogeneous boundary conditions and convergent discretizations of the considered operators are provided. Possible extensions of our main results to more general operators is outlined in section 6. In section 7, we present numerical experiments to support our results and to illustrate the performance of the presented numerical methods. Our conclusions and future research directions are given in section 8.

Definitions and concepts.
2.1. Notations and functional setting. Boundary conditions. Let Ω be a bounded domain in R n , u : Ω → R, κ : ∂Ω → R, and ∂u ∂ν be the derivative of u along the unit outward normal direction on ∂Ω. We denote by B one of the following boundary operators: which are known as Dirichlet, Neumann, and Robin boundary operators, respectively. For now we simply assume that the function u is regular enough so that all of the above makes sense.
Eigenvalue problem. Let Ω be a Lipschitz bounded domain and B one of the the above boundary operators. Then, there exists an orthonormal basis of L 2 (Ω) consisting of eigenfunctions {ϕ m } m∈N of the Laplacian operator −∆ coupled to homogeneous boundary conditions B(·) = 0, that is, non-trivial solutions of Moreover, the corresponding eigenvalues {λ m } m∈N are all non-negative and form a non-decreasing sequence such that λ m → ∞ (see e.g. [53]).

Functional spaces.
Given an open set Ω and r ≥ 0, H r (Ω) denotes the classical Sobolev space W r,2 (Ω) (cf. [21]). We denote by · r the classical norm of H r (Ω). Throughout the manuscript we will always use · 0 to indicate the norm in H 0 (Ω) = L 2 (Ω). H 1 0 (Ω) is the set of functions u ∈ H 1 (Ω) with null trace on ∂Ω. For any given r ≥ 0, let H r B (Ω) be the space of u ∈ L 2 (Ω) such that whereû m := u, ϕ m denotes the L 2 (Ω) inner product of u with the basis function ϕ m previously introduced. Most of the results in this manuscript will be stated in terms of the norm · H r B given by (2.5), for some suitable non-negative integer r. When r is a non-negative integer, H r B (Ω) has the following characterization (see Lemma 2.2 in [9] for the Dirichlet case, with similar proofs for the Neumann and Robin cases): (2.6) H r B (Ω) = {u ∈ H r (Ω) | B(∆ j u) = 0 on ∂Ω, ∀ non-negative integer j < r/2}, and the norm · H r B is equivalent to the usual norm · r in H r (Ω). The conditions imposed on ∆ j u at the boundary are known in the literature as compatibility conditions.

Fractional powers of the Laplacian in bounded domains.
Fractional powers with exponent s ∈ (0, 1) of the Laplacian coupled to the boundary condition B are defined in the spectral sense as for all functions u ∈ L 2 (Ω) such that m λ 2s m |û m | 2 < ∞ (see e.g., [11]). These conditions on u are a quite common requirement in the literature dealing with the spectral fractional Laplacian and we will proceed under these assumptions. However, wanting to work in full generality, the assumptions on the function u could be relaxed.
Let B be one of the boundary operators introduced in (2.1), (2.2), or (2.3). We let w(x, t) := e t∆ B u(x) denote the solution of the heat equation By separation of variables, the function w can be written as the following series: where λ m and ϕ m are the eigenpairs defined in (2.4). Let w ∞ denote the steadystate of (2.11). Using the series expansion of w we can easily see that w ∞ = 0 when λ m > 0 for all m (as in the Dirichlet case and the Robin case with κ > 0), while w ∞ =û 1 ϕ 1 = 1 |Ω| Ω u dx when λ 1 = 0 and λ m > 0 for m ≥ 2 (as in the Neumann case).
On a bounded domain Ω, the definition of (−∆ B ) s given by (2.4)-(2.7) is equivalent to the following one: Let Ω be a bounded domain, s ∈ (0, 1) and u ∈ H 2s B (Ω). The spectral B-fractional Laplacian of u in Ω is given by where e t∆ B u(x) is the solution of (2.8).
In section 5 we will extend definition (2.10) to a fractional-type operator associated with non-homogeneous boundary conditions.
where a(·, ·) is the bilinear form corresponding to the particular set of boundary conditions considered in (2.8) and V is a suitable functional space (i.e., H 1 0 (Ω) for homogeneous Dirichlet conditions or H 1 (Ω) in the homogeneous Neumann and Robin cases). In particular, a(w, v) = Ω ∇w ·∇v for the Dirichlet and Neumann cases, while a(w, v) = Ω ∇w · ∇v + κ w, v L 2 (∂Ω) in the Robin case, with ·, · L 2 (∂Ω) denoting the scalar product in L 2 (∂Ω). Note that the series expansion given by (2.9) still holds for the weak solution w (see e.g. Lemma 7.2-1 in [51]). We introduce here the notation and concepts of the FE theory that we will use throughout this manuscript. Triangulation. Let h > 0 and T h denote a triangulation ofΩ such that the diameter of each element K ∈ T h is smaller than h. Spatial nodes. Let Ω h denote the set {x i } N h i=1 of N h nodes in the considered triangulation. Time points. Let ∆t > 0 and t j = j∆t, for j = 0, 1, . . . , N t , be a uniform grid of time points used to discretize the interval [0, T ], for a given T > 0. Finite dimensional space and FE basis. Let P k , k ≥ 1, denote the space of polynomials of degree less than or equal to k and let X k h be the space of triangular FE, that is, Let u ∈ L 2 (Ω) and P k h be the L 2 (Ω)-orthogonal projection operator over the space X k h , that is, the operator P k h : Stiffness and Mass matrices. The stiffness matrix A and the mass matrix M corresponding to the FE approximation of the problem (2.11) on V h have entries is the FE basis of degree k. 2.4. Definition of the discrete operator. Given u ∈ H 2s B (Ω), we denote by Θ s the following discrete operator in time: where w(x, t j ) is the solution of (2.11) at the time points t j = j∆t, N t ∈ N is a suitable integer, the β j are weights defining a suitable quadrature rule for (2.10), w ∞ is the steady-state of (2.11), and β ∞ is the value of the tail of the integral with respect to the measure dt t 1+s (see Theorem 3.1 and Theorem 4.1 for further details). Combining (2.14) with the FE we get the following approximation of the spectral B-fractional Laplacian (−∆ B ) s of a given u: Definition 2.3. Let w h be the FE solution of degree k of (2.11) defined by equation (2.13), corresponding to the initial datum u. Then, for all x ∈ Ω, the approximation of the B-fractional Laplacian of u is defined as Remark 1. In section 3 and section 4 we propose two different choices of the weights β j leading to approximations of (−∆ B ) s u with different accuracy. The weights β j defined in (3.2) correspond to a composite midpoint quadrature rule while the ones in (4.2) correspond to a quadrature defined via a piecewise linear interpolant.
3.1. Assumptions. In order to prove our main convergence results, we need some assumptions on the domain Ω, the triangulation T h , the function u, the spatial and the temporal steps h and ∆t, respectively. Specifically, we assume: T h , and ρ K := sup{diam(S)|S is a ball contained in K}, there exist two constants, σ ≥ 1 and τ > 0, such where c µ is a positive constant independent from both ∆t and h, derived from a bound for the largest eigenvalue of the bilinear form a(·, ·) on the finite-dimensional space V h (see e.g., Chapter 6 in reference [49] for details).
Remark 2. We have assumed the domain to be a convex polytope so that the technicalities related to FE approximations in the case of curvilinear boundaries and re-entrant corners can be avoided. However, assumption (A Ω ) can be relaxed to Ω being a bounded domain with Lipschitz boundary ∂Ω. The case of Lipschitz domains with curvilinear boundaries can be dealt with, for example, as proposed in [4], while suitable mesh refinements near re-entrant corners allow the recovery of optimal error bounds for non-convex domains as shown in the work by Chatzipantelidis et al. [13] and references therein.

Low order approximation.
We now state our first main result, which concerns homogeneous boundary conditions.
hold, and let Θ s h be defined by (2.15) where: (a) w h is the FE solution of degree k of (2.11) with initial datum u; , λ min being the first nonzero eigenvalue of (−∆ B ); (c) the weights β j , for j = 1, . . . , N t , are given by Then, there exists a constant C > 0 independent of h such that 3.3. Proof of Theorem 3.1. In order to prove Theorem 3.1 we need some intermediate results. We split the integral in (2.10) into three parts: for some . For each of these integrals, in subsection 3.3.1, we define an approximation, that we denote by I D S , I D M , I D ∞ , so that In subsection 3.3.1 we provide a bound to (−∆ B ) s u−Θ s u 0 while in subsection 3.3.2 we estimate Θ s u − Θ s h u 0 . The proof of Theorem 3.1 follows from these preliminary results.
3.3.1. Intermediate results. In the following lemmas, we work under the assumptions (A Ω ), (A B ), and that the regularity of u is as given explicitly by the constants in each error estimate.
Proof. Using the series expansion of u and w = e t∆ B u (see (2.9)), By orthonormality of the functions ϕ m , Therefore, exploiting Minkowski's integral inequality We split the interval [∆t/2, (N t + 1/2)∆t] into N t subintervals of size ∆t, denote w = e t∆ B u, and rewrite I M [u] as the following finite sum of integrals: where t j = j ∆t, for j = 1, . . . , N t . By using the expression (3.7) for I D M [u] we obtain so that, by using the triangular inequality and Minkowski's integral inequality, Using the series expansion of w and the fact that for all λ ≥ 0, e −λt is a convex and monotonically decreasing function of t, one finds that on each subinterval I j = [t j − ∆t/2, t j + ∆t/2], Therefore, w(·, t) − w(·, t j ) 0 ≤ ∆t u H 2 B for all j and all t ∈ I j , which leads to Lemma 3.4. Let s ∈ (0, 1), I ∞ given by (3.6) for some T 2 > 0, and Let λ min be the first nonzero eigenvalue of (−∆ B ). Then there exists a constant Proof. Recall that w ∞ = 0 for the Dirichlet and Robin cases, while w ∞ =û 1 ϕ 1 for Neumann boundary conditions, and let m min be the index corresponding to λ min (i.e., for the boundary operators of interest, m min = 1 in the Dirichlet and Robin cases, while m min = 2 in the Neumann case). Then, for all Hence, Proposition 3.5. Let assumptions of Theorem 3.1 hold. Let (−∆ B ) s and Θ s be defined by (2.10) and (2.14), respectively. Then, there exists a constant c a = c a (s, Proof. The result follows by simply combining Lemma 3.2, Lemma 3.3, and Lemma 3.4. Finally, using the definition of N t and choosing T 2 = (N t + 1/2)∆t we obtain that e −λ min T 2 T s 2 ≤ ∆t 1−s , and a suitable choice of the constant c a completes the proof.

Finite element approximation.
In this section we combine some wellknown estimates of the FE method to obtain a bound for the error made by approximating Θ s u by Θ s h u. As usual, we denote by w the solution of the heat equation with initial condition u, and by w h the corresponding FE approximation of order k.
Recalling that w(·, 0) = u and that w h (·, 0) = P k h u is the L 2 (Ω)-orthogonal projection of the function u, from classical estimates for the FE projection error (e.g., Section 3.5 in [49]) we have for some positive constant c independent from h. In order to bound w(·, t j ) − w h (·, t j ) 0 for t j > 0 we use the fact that there exists a constantc independent from h, ∆t, and t j such that The above result, for t j > 1 can be obtained via classical FE estimates which hold even for rough initial data, u ∈ L 2 (Ω) (see e.g., Remark 11.3.2 in [49] or Theorems 1 and 2 in [39]), while for t j ≤ 1, the bound (3.10) is a consequence of the regularity assumptions made on u, ensuring that the constantc can be controlled by the integral (over [0, 1]) of suitable spatial norms of the heat solution w and its temporal derivatives (see e.g., Theorem 11.3.4 in [49] or Theorem 2.1 in [8]).
Proposition 3.6. Let assumptions of Theorem 3.1 hold and Θ s u given by (2.14).
Proof. Via the triangular inequality we find On one hand, by using (3.9), we can bound the sum of terms in (3.11) involving w(·, 0) and w h (·, 0) as follows: On the other hand, for the remaining terms of (3.11) we use (3.10) and obtain Hence, combining the above estimates we get

Proof of the main result.
Proof of Theorem 3.1. The triangular inequality yields Proposition 3.5 provides the following bound of the first term in (3.12): On the other hand, from Proposition 3.6 we have that Hence, the order of convergence is always determined by h p(1−s) , and by suitably choosing the constant C in (3.3) we conclude the proof.

4.
Higher order discretization. Since s ∈ (0, 1), the rate of convergence presented in Theorem 3.1 can be slow for practical computations when s is close to one. In this section, we present a discretization that allows to obtain a better rate of convergence at the expenses of requiring some extra regularity to the function u, and also to the boundary operator in the Robin case. We adopt the Crank-Nicolson temporal scheme, corresponding to setting θ = 1/2 in the FE solution of the heat equation given by (2.12). We suitably modify the relationship between the discretization parameters ∆t and h, and consider higher order quadrature approximations for both I S and I M .   (c) the weights β j are defined as In the following lemmas, we work under the assumptions (A Ω ), (HA B ), and that the regularity of u is as given explicitly by the constants in each error estimate.
Using the properties of the exponential, for all λ m ≥ 0 we have Hence, Letting w = e t∆ B u, using the series expansion of w and u, and from (4.6), we obtain For the middle part, we consider the piecewise linear interpolant of e t∆ B u(x)−u(x) on a grid of equally spaced points t j = j∆t, j = 1, . . . , N t . Given a function φ : [0, +∞) → R, the linear interpolant of φ is defined as Note that the integrals ρ j are positive since the tent functions P ∆t (t−t j ) are positive for all j. Moreover, following the same procedure used in [31] for the exact computation of quadrature weights we find: where F and F are defined as in (4.3).
Remark 3. Note that the weights β j in (4.2) are equal to ρ j for all j except for j = 1 where an extra positive term (coming from the discretization of the singular part of the integral in the definition of (−∆ B ) s ) is added to ρ 1 .
For notational convenience, we denote g m (t) = e −λmt − 1. Then, using (2.9) it is straightforward to show that and thus, On each subinterval I i , P[g m ] is the Lagrange linear interpolant of g m and hence there exists a constant c independent of ∆t and t j such that for all t ∈ I i , Consequently, for all t ∈ I i ,   = c a (s, u H Proof. We proceed as in the proof of Proposition 3.5 using now Lemma 4.2, Lemma 4.3 and Lemma 3.4 to get Using the definition of N t and choosing T 2 = N t ∆t, we get e −λ min T 2 We follow the exact same procedure and notation used in the proof of Proposition 3.6. The bound (3.9) for the projection error still holds. The FE error for t j > 0 can be now improved exploiting the fact that the second-order Crank-Nicolson discretization scheme is used in time, leading to an error O(h k+1 + ∆t 2 ) (see e.g., Corollary 11.3.1 in [49] and [40]). Thus, we find By definition, β j > 0 for all j. Consequently, Thus, following the proof of Proposition 3.6, we conclude that We can now prove the higher order discretization result.
Proof of Theorem 4.1. Following the same steps done in the proof of Theorem 3.1, now exploiting the estimates of Proposition 4.4 and Proposition 4.5, we obtain , By using the relationship between h and ∆t given by (4.1), we have h k+1 ∆t −s = η −s h k+1−ps and ∆t 2−s = η 2−s h p(2−s) . Moreover, k + 1 − ps ≥ p(2 − s) for p ≤ k+1 2 . Therefore, the order of convergence is determined by ∆t 2−s = h p(2−s) and by suitably choosing the constant C in (4.4) we obtain the desired result.

5.
Nonlocal operators with non-homogeneous boundary conditions. In many practical physical applications, the condition required at the boundary is not necessarily homogeneous. In this section we define and deal with a family of nonlocal operators of fractional type that are compatible with non-homogeneous Dirichlet, Neumann or Robin boundary conditions. 5.1. The nonlocal operator. Consider u : Ω → R, satisfying B(u) = g on ∂Ω for some g : ∂Ω → R. Let v(x, t) be the solution of the following heat equation with initial condition u and non-homogeneous boundary conditions: Generalizing Definition 2.1, we define the nonlocal operator L s g on the given u as We intentionally avoid using the notations e t∆ and (−∆) s when dealing with nonhomogeneous boundary conditions since L s g is not a fractional power of the Laplacian coupled to non-homogeneous boundary conditions. However, as we show in Sec- where z is the harmonic extension of g to Ω. Clearly, the same strategy previously proposed in this work can be followed to obtain numerical approximations of (5.2): the solution of the heat equation (5.1) can be approximated via FE accounting for non-homogeneous boundary conditions, while the quadrature formulas for the approximation of the singular integral (5.2) remain unchanged. Moreover, the above-mentioned characterization of the nonlocal operator L s g in terms of (−∆ B ) s can be exploited as a powerful tool in the numerical approach (see Section 5.3).
After the first submission of this manuscript we came across the work by Antil et al. [2] in which they propose an alternative way to deal with non-homogeneous boundary data. However, in the very recent work by Lischke et al. [38], our definition (given by (5.1)-(5.2)) and the one presented in [2] are shown to be equivalent. Moreover, the authors of [38] also present a nice and detailed comparison of the two approaches from the numerical point of view.

Characterisation of L s
g . The solution of (5.1) can be written as the sum v( is the solution of the heat equation with homogeneous boundary conditions and shifted initial condition, i.e., and z is the steady-state of (5.1), i.e., the solution of the harmonic extension problem The following characterization result provides an alternative definition of L s g .
The condition requiring the solution to the harmonic extension problem (5.4) to be z ∈ H r (Ω) for some r ≥ 0 relies heavily on the geometry of the domain Ω, the specific type of boundary conditions B imposed, and the regularity of the boundary data g. Providing a general theory able to encompass all these regularity considerations is far from the main goal of the present article. However, for illustration purposes, we focus on a particular example and show how the numerical tools derived in the previous sections can be applied in these settings.

Numerical approximation in the non-homogeneous case.
For simplicity, in this section we restrict ourselves to the case of Dirichlet boundary conditions in a two-dimensional convex polygonal domain Ω (thus satisfying (A Ω )).
Consider a function u ∈ H 2 (Ω) and let g = u| ∂Ω be the trace of u in Ω. Classical regularity theory ensures that g ∈ H 3 2 (∂Ω) and hence, the solution z of (5.4) is such that z ∈ H 2 (Ω) (see for example Theorem 1.5 and Theorem 1.8 in [26]). Thus, u − z ∈ H 2 (Ω) and (u − z)| ∂Ω = 0, via (2.6) we obtain that u − z ∈ H 2 B (Ω). In practice, given a function u, z will not be a known analytically. To overcome this issue, we considerz h to be the FE approximation of degree k = 1 of problem (5.4) for a given triangulation. Namely,z h = N h i=1 Z i φ i , according to the notation previously used in this manuscript. Similarly to (5.3), letw be the solution of (5.5) andw h be the corresponding FE approximation of degree k = 1.
The following theorem provides a numerical approximation of L s g [u] and, exploiting the result of Theorem 3.1, gives a bound for the L 2 -norm of the error.
Theorem 5.2. Let s ∈ (0, 1), θ ∈ [0, 1], k = 1 and h, ∆t > 0 such that (3.1) holds. Let Ω ⊂ R 2 , assume that (A Ω ) and (A T h ) hold, and B is the Dirichlet boundary operator. Let u ∈ H 2 (Ω) and g := u| ∂Ω . Letz h be the FE solution of degree k of (5.4) and Θ s h be defined according to (2.15), with N t , β j and β ∞ given by (b) and (c) from Theorem 3.1. Then, Hence, the triangular inequality gives Clearly, the result of Theorem 3.1 (for k = 1) can be applied to obtain On the other hand, To bound the right-hand side, we use the stability of the FE method (see e.g. [49] for details) with respect to the initial data, ensuring that for all j ≥ 0, wherec =c(θ). At this point we use the linearity of the projection operator P k h and the fact that, by definition, P k h is exact on each element of the triangulation T h for polynomial functions of degree k or less. Hence, by construction, P k h (z h ) =z h . Moreover, the accuracy of both the interpolant of z and the FE approximation deriving from the solution of (5.4) are of order h k+1 , so that Proceeding as in the previous sections, we combine all estimates and get The relationship between ∆t and h given by (3.1) concludes the proof.
6. Generalisations. The results of Theorem 3.1 and Theorem 4.1 can be generalized not only to non-homogeneous boundary conditions but also to a wider class of operators coupled to homogeneous boundary conditions. Instead of the Laplacian operator −∆ one can consider a more general positive second order elliptic operator of the form Like in Definition 2.1, the fractional powers of the operator L with the boundary condition B can be expresed as where e −tL B denotes the heat-semigroup associated to L and B. Analogous results to Theorem 3.1 and Theorem 4.1 can be proved for the operator (L B ) s provided that the bilinear form corresponding to the operator L B in the associated variational formulation is weakly coercive (see e.g., [51,53] for precise conditions on L and B), and that the coefficients a ij , b j , c, the domain Ω and the function u ensure a sufficiently smooth e −tL B u(x).
Regarding the boundary conditions, in this work, we have considered the case of Dirichlet, Neumann, and Robin boundary operators, assuming that a single type of condition is applied to the whole boundary ∂Ω. However, more general boundary constraints (e.g., mixed conditions in different parts of the boundary) can be considered, provided that the properties of weak coercivity of the bilinear form and the regularity of e −tL B u(x) needed for the proof can be ensured.
Details on these regularity conditions for both the operator L and more general boundary conditions can be found in reference [37].

Numerical experiments.
In order to illustrate the performance of our numerical methods we perform various numerical experiments and report here our results. The implementations of the methods are done in MATLAB (R2016a, The MathWorks Inc., Natick, Massachusetts, US).
Let h, θ, k be fixed parameters of the FE method, and let ∆t, N t be chosen according to Theorem 3.1 and Theorem 4.1. We solve iteratively the linear system P W (j) = QW (j−1) , where P := (M + θ∆tA) and Q := (M + (θ − 1)∆tA) are the N h × N h square matrices of the FE scheme, we compute the difference vector (W (j) − W (0) ), multiply it by the weights β j , and sum all terms to finally obtain Here W ∞ is the constant vector corresponding to the steady-state of the heat equation (which is known explicitly and depends on the boundary conditions), while W (0) is the vector with coordinates given by L 2 -projection of the function u on the nodes of the spatial grid Ω h .
To test the convergence of the discretizations, we use the fact that when ϕ m is a nomalized eigenfunction of −∆ with the boundary condition B, then the B-fractional Laplacian of ϕ m can be computed analytically according to (2.7) as (−∆ B ) s ϕ m = λ s m ϕ m , where λ m is the eigenvalue corresponding to ϕ m . For Figure 1 we select three values of s ∈ (0, 1), namely s = 0.25, 0.5, 0.75, and compute the error (−∆ B ) s u − Θ s h u 0 as a function of the spatial discretization parameter h, for the Dirichlet, Neumann, and Robin boundary cases. The function u is always a normalized eigenvalue of the corresponding operator, i.e., u = ϕ m /||ϕ m || 0 , for some m ≥ 1. All reported results were obtained with linear FE (k = 1), an implicit temporal discretization scheme (θ = 1), and ∆t = η h k+1 .
In the left column of Figure 1 we compare the expected behavior (dashed lines) predicted by our convergence result Theorem 3.1 with the error decay obtained numerically (solid lines) for different boundary conditions and different values of s. On the right, for each boundary condition, we show the exact value of (−∆ B ) s u for the three values of s considered.
In Figure 2 we show the improved convergence obtained for the same functions considered in the tests of Figure 1 when the higher order discretization method described in section 4 is used. We observe that the convergence decay actually outperforms the trend predicted from our theoretical results for all values of s and that the behavior improves as s increases.
As mentioned various times in this work, the B-fractional Laplacian given by Definition 2.1 is a boundary-dependent nonlocal operator, that is, the operator changes when different boundary conditions B are considered. This can be seen by applying the operator for different choices of B to the same function u. Letting 1 D denote the indicator function of the set D, in Figure 3, we consider a compactly supported function of the form u(x) = e −1/(r 2 −x 2 ) 1 |x|<r , with 0 < r < 1, on Ω = (−1, 1), and apply to it the B-fractional Laplacian with homogeneous Dirichlet, Neumann, and Robin boundary conditions for different values of the fractional parameter s. The fractional Laplacian is computed via the higher order method.
For comparison purposes, in the top left of Figure 3 we plot the standard Laplacian of the considered u, which is boundary independent. As we can see, the standard Laplacian of a compactly supported function is still a compactly supported function, whereas this does no longer hold in the purely fractional cases. The nonlocal character of the operator in fact affects the result and produces functions that are nonzero also on Ω \ [−r, r]. Moreover, differences between the three nonlocal operators considered are clearly visualized for all choices of s < 1.
7.2. Two dimensional case. The spectrum of the classical Laplacian is known analytically also in the two-dimensional case for some simple regular geometries, such as rectangle, disk, or ellipse (see [28] for further details). To prove convergence of our method in the two-dimensional case we exploit these analytical results and use them to compute (−∆ B ) s u − Θ s h u 0 as a function of the spatial discretization parameter h, for the usual three different boundary conditions B, on the unit square Ω = (0, 1) 2 . In this case, eigenfunctions are simply products of the one-dimensional eigenfunctions corresponding to the particular boundary conditions considered, namely ϕ ml (x, y) = ϕ m (x)ϕ l (y), and the corresponding eigenvalues are given by the sum of the one-dimensional counterparts, that is, λ ml = λ m + λ l . The convergence results and the index of the particular eigenfunction ϕ ml used for the test are given in Figure 4.
In Figure 5 we provide an example of the use of our method to compute the action of the fractional operator L s g on a function u satisfying non-homogeneous Dirichlet boundary conditions on the unit square Ω = (0, 1) 2 . Specifically, we set g(x, y) = 0.1 sin(2πx) cos(πy) for (x, y) ∈ ∂Ω. Wanting to study convergence of the method we let ϕ := u − z be a normalized eigenfunction of the homogeneous Dirichlet Laplacian on the unit square so that (−∆ B ) s ϕ is known explicitly.
Not knowing the analytic expresssion of z, we compute its FE approximation z h and hence obtain the FE approximation u h = ϕ + z h (satisfying non-homogeneous boundary conditions) to be used as initial datum for the approximate computation of L s g . The solution is computed with the higher order method and the error once again decays faster than expected for all three values of s considered.  1). However, the analytic expression of the required eigenpairs is not known for general domains and using (1.1) would imply that the solution of a very large number of eigenvalue problems must be computed.
One of the main advantages of our approach and of the use of the numerical method proposed in this work is that with the same strategy we can handle boundary conditions that are not necessarily homogeneous Dirichlet conditions and at the same  Figure 6, for which the eigenpairs are analytically unknown. Although in our theoretical results, the choice of N t depends on the first nonzero eigenvalue of the Laplacian, clearly the same results hold if a lower bound is used for the considered eigenvalue. These bounds are readily available and we refer the reader to Section 4.2 of [28] and references therein for a thorough discussion on these bounds for general domains and homogeneous Dirichlet, Neumann, and Robin boundary conditions. An alternative approach not using bounds for the first nontrivial eigenvalue simply consists in continuing to add terms to the quadrature sum in (7.1) until the difference (W (j) − W (0) ) reaches a steady value, up to a given tolerance. The latter approach is the one we adopted for the computations presented in Figure 6.
Once again the value of the B-fractional Laplacian with s = 0.5 is computed for the compactly supported function with homogeneous Dirichlet, Neumann, and Robin boundary constraints. Different views of (−∆ B ) s u are reported to help in visualizing the differences (especially in the . behavior close to the boundary) between these three cases.

7.
3. An explicit scheme for a fractional porous medium type equation. As a possible application of our method, we compute the numerical approximation of an evolution problem involving the homogeneous Dirichlet fractional Laplacian. Specifically, on a given bounded domain Ω and for τ ≥ 0, we consider the fractional porous medium type equation where m is a suitable integer, m > 1. For a precise presentation of the theoretical aspects regarding (7.2) we refer to [5,6] and references therein. Given an initial condition u(x, 0) = u (0) (x), we use an explicit discretization scheme for the evolution variable τ with uniform step ∆τ > 0 and compute iteratively the numerical approximation of the solution to (7.2) at τ n = (n + 1)∆τ , n ∈ N as As (7.3) is an explicit scheme, it may suffer from numerical instabilities when the time step ∆τ is not properly chosen. We circumvent the issue by imposing on ∆τ a CFL-type condition and setting ∆τ = h 2s /m. More details about these facts can be found in [19] dealing with (7.2) in R N .
In Figure 7 we report the results obtained for two different combinations of the parameters m and s, at the indicated time point τ . In both cases, the spatial domain is Ω = (−1, 1), while the initial condition (plot on the left in Figure 7) is the compactly supported function u (0) (x) = e 4−1/[(0.5−x)(0.5+x)] 1 |x|<0.5 .
For the particular combinations of m and s presented in Figure 7 the theoretical results of Bonforte et al. [5] predict a long-time boundary behavior of the type where φ 1 is the first eigenfunction of the Laplacian with homogeneous Dirichlet boundary condition in Ω. Given a, b ∈ R + we follow the notation used in [5] and write a b whenever there exist universal constants c 0 , c 1 > 0 such that The predicted boundary behaviour presented in [5] is recovered by our simulations and shown in the second and third plots of Figure 7, where the numerical solution obtained with our method (blue) is compared to the expected function v(x, τ ) (red).

Conclusions.
In the present paper we investigate a method to define a nonlocal operator acting on bounded domains that, in the unbounded case, recovers the classical definition of the fractional Laplacian (−∆) s . The studied method is characterized by the fact that the imposition of local boundary conditions is sufficient to obtain a well-posed formulation of the nonlocal operator without requiring additional constraints to be imposed on the complement of the considered finite domain, in analogy with local problems. As pointed out in the introduction, this characteristic is relevant for applied problems and for the physical foundation of fractional models.
We consider for our purposes the spectral definition of the fractional Laplacian and exploit the equivalence between this definition and the heat-semigroup formulation of the operator on bounded domains. The latter definition allows for a much more intuitive and transparent formulation in which the role played by the local boundary conditions (passed on to the local heat problem) is clear and all ambiguity of having local conditions for a nonlocal operator is removed.
This results in a power-law weighted integration of the evolving discrepancy between the solution of the heat equation and its initial condition. Within this framework, we propose a novel discretization scheme. An FE method is adopted to obtain the approximated solution of the heat equation and suitable quadrature formulas are used to approximate the integral.
Unlike other discretization methods, the proposed approach does not rely on the computation (or the approximation) of the spectrum of the Laplacian operator on the considered finite domain. In fact, a key feature of our strategy lies in the flexibility offered by FE in handling well not only simple regular geometries but also general bounded domains for which explicit knowledge of the Laplacian spectrum is missing. Moreover, once the accurate computation of the heat equation solution is available, our method simply requires the approximation of a one-dimensional integral, independently from the dimension in which the fractional Laplacian of the given function is sought.
In this work, we start by considering the case of homogeneous boundary conditions and propose two discretizations of the fractional operator (depending on a discretization parameter h) exhibiting two different rates of convergence. The main difference between the two results lies in the quadrature approximation of the powerlaw weighted integral. In the first discretization, we truncate the integral by ignoring its singular part (near the origin), adopt a naïve mid-point quadrature rule for the middle part, and exploit the exponential convergence to the steady-state of the solution of the heat equation to approximate the tail. Doing so and suitably choosing the parameters of the FE computation, we obtain an error that behaves like O(h p(1−s) ) for some positive p and for values of the fractional order s ∈ (0, 1). This approach produces satisfactory results for values of s close to its lower bound. By adapting the quadrature of the integral representing the operator and ensuring that the FE provide the required higher temporal accuracy for the solution to the heat equation, we derive a second discretization resulting in an error which is O(h p(2−s) ). Although the proof of this result comes at the expenses of requiring some higher regularity to the function on which the action of the fractional Laplacian is computed, the accuracy of this second discretization ensures a decay of the error at least of the order O(h p ), and hence this discretization is suitable to deal with the case of s close to one.
The case of non-homogeneous boundary conditions is considered next. We propose a definition of a fractional-type operator consistent with these inhomogeneous constraints and formulated as the power-law weighted integral of the evolving discrepancy of a heat equation solution, in which non-homogeneous boundary conditions are imposed, and its initial condition satisfying these conditions. We show that such an operator can be rewritten as the fractional Laplacian coupled to homogeneous boundary conditions, for a suitable shifted function, and we use the same two discretizations obtained in the homogeneous case to derive corresponding discretizations of the new nonlocal operator.
Our numerical results show performance of our method for both one and two dimensional domains. The theoretical rate of convergence is very well recovered in all the analyzed cases and in the two-dimensional setting we also illustrate the performance of our method on a general irregular bounded domain. We conclude the numerical section with an example of a parabolic problem for which the computation of the fractional Laplacian could be used. In particular, we consider a fractional porous medium equation in which homogeneous Dirichlet conditions are imposed at the boundary. We compare our numerical solution with the analytic results on the boundary behavior of the solution given in [5] and show that the expected behavior is well captured by our method, for both combinations of parameters considered.
The approach followed in this work provides a natural framework that can be used to obtain either discretizations of the fractional Laplacian defined via the heatsemigroup on the whole R n or to compute discretizations of inverse fractional powers of second order operators (i.e. solutions of elliptic problems associated to fractional operators) on bounded domains. In the unbounded settings, the use of FE to accurately model non rectangular domains is not required anymore due to the absence of a boundary, and the simpler finite difference approach can be adopted. On bounded domains, inverse fractional powers of second order operators can still be represented via integral expressions involving the corresponding heat-semigroups. The kernel of these expressions become locally integrable near the origin, avoiding some of the difficulties we had to deal with in the present work, but integrability is lost at infinity. Nevertheless, once again one can exploit the exponential convergence of the solution to the heat equation towards a steady-state to overcome this issue. Both these generalizations are current research topics of the authors of this paper, and the outlined results will be presented in two forthcoming papers.