A comparative study on nonlocal diffusion operators related to the fractional Laplacian

In this paper, we study four nonlocal diffusion operators, including the fractional Laplacian, spectral fractional Laplacian, regional fractional Laplacian, and peridynamic operator. These operators represent the infinitesimal generators of different stochastic processes, and especially their differences on a bounded domain are significant. We provide extensive numerical experiments to understand and compare their differences. We find that these four operators collapse to the classical Laplace operator as \alpha \to 2. The eigenvalues and eigenfunctions of these four operators are different, and the k-th (for k \in N) eigenvalue of the spectral fractional Laplacian is always larger than those of the fractional Laplacian and regional fractional Laplacian. For any \alpha \in (0, 2), the peridynamic operator can provide a good approximation to the fractional Laplacian, if the horizon size \delta is sufficiently large. We find that the solution of the peridynamic model converges to that of the fractional Laplacian model at a rate of O(\delta^{-\alpha}). In contrast, although the regional fractional Laplacian can be used to approximate the fractional Laplacian as \alpha \to 2, it generally provides inconsistent result from that of the fractional Laplacian if \alpha \ll 2. Moreover, some conjectures are made from our numerical results, which could contribute to the mathematics analysis on these operators.


Introduction
In the last couple of decades, nonlocal or fractional differential models have become a powerful tool for modeling challenging phenomena including anomalous transport, long-range interactions, or from local to nonlocal dynamics, which cannot be described properly by integer-order partial differential equations. So far, numerous fractional differential models have been proposed, among which models with the fractional Laplacian have been well applied. The fractional Laplacian operator (−∆) α/2 , representing the infinitesimal generator of a symmetric α-stable Lévy process, has been used to model anomalous diffusion or dispersion [16,13], turbulent flows [6,39], systems of stochastic dynamics [4,11], finance [12], and so on. Nevertheless, the fractional Laplacian is a nonlocal operator defined on the entire space. Moreover, an equation involving the fractional Laplacian has to be enclosed by a nonconventional, nonlocal boundary condition imposed on the complement of the physical domain where the governing equation is defined. The combination of these two facts introduces many significant challenges in the mathematical modeling, numerical simulations, and corresponding mathematical analysis. These issues have not been encountered in the context of integer-order partial differential equations.
To avoid evaluation and analysis over the entire space, one common approach in the literature is to "truncate" and approximate the integral of the fractional Laplacian. Hence, some other nonlocal operators that are closely related to the fractional Laplacian have been proposed in recent years, including the regional fractional Laplacian, the spectral fractional Laplacian, and the peridynamic operator. Similar to the fractional Laplacian, these operators are nonlocal, and on the entire space they are equivalent to the fractional Laplacian. On a bounded domain, the spectral fractional Laplacian is defined via the spectra of the classical Laplace operator on the same domain [11,1]. Thus, local boundary conditions are imposed to the equations involving the spectral fractional Laplacian. In contrast to the fractional Laplacian, the integration domain in the regional fractional Laplacian is reduced from the entire space to a finite domain where the governing equation is defined [23,24,9]. Consequently, the regional fractional Laplacian model considerably simplifies the numerical computations of the original nonlocal problem with the fractional Laplacian.
Peridynamic model was originally proposed as a reformulation of the classical solid mechanics [40]. The classical theory of solid mechanics assumes that all internal forces act through zero distance, and the corresponding mathematical models are expressed in terms of partial differential equations, which cannot describe problems with spontaneous formation of discontinuities and other singularities. The peridynamic model leads to a nonlocal framework that does not explicitly involve the notion of deformation gradients and thus provides a more accurate description of problems with discontinuities and singularities. The peridynamic operator has been recently used to approximate the fractional Laplacian so as to reduce the computational costs [13,21].
In this paper, we study and compare the properties of the fractional Laplacian, spectral fractional Laplacian, regional fractional Laplacian, and the peridynamic operator. We show that the spectral fractional Laplacian and regional fractional Laplacian are significantly different from the Dirichlet fractional Laplacian, although they can be used to approximate each other as the power α → 2. In contrast, the peridynamic operator can provide a consistent approximation to the fractional Laplacian for any α ∈ (0, 2), but a large horizon size is required to obtain a good approximation, especially when α is small. The rest of this paper is organized as follows. In Section 2, we introduce the fractional Laplacian and its related nonlocal diffusion operators, including the spectral fractional Laplacian, regional fractional Laplacian, and peridynamic operator. In Section 3, we numerically compare these four operators by studying their nonlocal effects, eigenvalues and eigenfunctions, solutions to their corresponding Poisson equations, and dynamics of the nonlocal diffusion equations. Finally, we draw some conclusions in Section 4.

Nonlocal diffusion operators
In this section, we introduce the fractional Laplacian and its related nonlocal diffusion operators, including the spectral fractional Laplacian, regional fractional Laplacian, and peridynamic operator. The peridynamic operator with a specially chosen kernel function can be used to approximate the fractional Laplacian [13,21]. If a bounded domain is considered, the spectral fractional Laplacian and regional fractional Laplacian are significantly different from the fractional Laplacian, although they are freely interchanged with the fractional Laplacian in some literature. In the following, we will introduce and compare the properties of these four operators from various aspects. Let Ω ⊂ R n (n = 1, 2, or 3) denote an open bounded domain, and Ω c = R n \Ω represents the complement of Ω.

Fractional Laplacian
On the entire space R n , the fractional Laplacian (−∆) α/2 is defined via a pseudo-differential operator with symbol |ξ| α [29,36]: for α > 0, (2.1) where F represents the Fourier transform, and F −1 is the inverse Fourier transform. From the probabilistic point of view, the fractional Laplacian (−∆) α/2 represents the infinitesimal generator of a symmetric α-stable Lévy process [3,7]. In a special case of α = 2, the definition in (2.1) reduces to the standard Laplace operator −∆. The definition in (2.1) enables one to utilize the fast Fourier transform to efficiently solve problems involving the fractional Laplacian, however, it is suitable only for problems defined either on the whole space R n or on a bounded domain with periodic boundary conditions. In the literature, an equivalent hypersingular integral definition of the fractional Laplacian (−∆) α/2 is given by [36,43,29,15]: where P.V. stands for the principal value, and c n,α is the normalization constant given by with Γ(·) denoting the Gamma function. In contrast to (2.1), the definition in (2.2) can easily incorporate with non-periodic bounded domains. Note that the integral representation in (2.2) is defined for 0 < α < 2, while the pseudo-differential definition in (2.1) is valid for all α > 0. The equivalence of definitions (2.1) and (2.2) for α ∈ (0, 2) are studied in [14,Proposition 3.3]. More discussion can be found in [14,36,28,45] and references therein. The fractional Laplacian on a bounded domain is of great interest, not only from the mathematical point of view, but also in practical applications. Recently, many studies have been carried out on the Dirichlet fractional Laplacian (also known as the restricted fractional Laplacian), i.e., the fractional Laplacian on a bounded domain Ω with extended homogeneous Dirichlet boundary condition (u(x) ≡ 0 for x ∈ Ω c ). However, the current understanding of this topic still remains limited, and the main challenge is from the non-locality of the operator. In the following, we will discuss some fundamental properties of the Dirichlet fractional Laplacian.
Probabilistically, the Dirichlet fractional Laplacian (−∆) α/2 represents the infinitesimal generator of a symmetric α-stable Lévy process that particles are killed upon leaving the domain Ω [7,8,45,42]. One fundamental issue in the study of the Dirichlet fractional Laplacian is its eigenvalues and eigenfunctions. So far, their exact results still remain unknown, and only some estimates and approximations can be found in the literature. In [11], it shows that the k-th eigenvalue λ k (for k ∈ N) of the Dirichlet fractional Laplacian on a convex domain Ω ⊂ R n is bounded by: where µ k represents the k-th eigenvalue of the standard Dirichlet Laplace operator −∆ on the same domain Ω. That is, the eigenvalue of the fractional Laplacian is always smaller than that of the standard Laplacian −∆. If a one-dimensional (i.e., n = 1) domain is considered, the estimates in (2.4) can be improved, and sharper bounds can be found for two special cases, such as k = 1 and α ∈ (0, 2) in [3,19], and α = 1 and k = 1, 2, 3 in [3]. More discussion on the eigenvalue bounds can be found in [26,44,45,20] and references therein. Furthermore, in a one-dimensional interval (−1, 1), the asymptotic approximation of the eigenvalue λ k is given by [27]: It further shows that if α ≥ 1, the eigenvalue λ k (for k ∈ N) is simple, and the corresponding eigenfunction satisfies φ k (−x) = (−1) k−1 φ k (x). Compared to the understanding of eigenvalues, the knowledge of eigenfunctions is even less. As shown in [38], on a smooth bounded domain Ω ⊂ R n , the eigenfunctions φ k (for k ∈ N) are Hölder continuous up to the boundary. Recent numerical results on eigenvalues and eigenfunctions of the Dirichlet fractional Laplacian can be found in [18]. The fractional Poisson equation is one of the building blocks in the study of fractional PDEs. It takes the following form [19,13,2,17]: In (2.7), the extended homogeneous boundary conditions are imposed on the complement Ω c , distinguishing from the classical Poisson problem where boundary conditions are given on ∂Ω. This difference can be explained from probabilistic interpretation of the standard and fractional Laplacian. The standard Laplace operator represents the infinitesimal generator of a Brownian motion with continuous sample paths; thus for a particle in domain Ω, it must leave the domain via the boundary points on ∂Ω. By contrast, the fractional Laplacian is the infinitesimal generator of a symmetric α-stable Lévy process with discontinuous sample paths; particles may "jump" out of the domain without touching any boundary points on ∂Ω. Hence, the solution on Ω can be determined by the values at ∂Ω in the context of classical Poisson equations but not in the context of fractional Poisson equations. The nonlocal problem (2.6)-(2.7) plays an important role in studying stationary behaviors of various problems. Its solutions can be expressed in terms of the Green's function, however, it is challenging to find the explicit form of the Green's function. Some estimates of the Green's function can be found in [5,19,27,8,10]. The regularity of the solution to (2.6)-(2.7) is discussed in [2,34,37,35,33]. As shown in [34], assume that Ω ⊂ R n is a bounded Lipschitz domain which satisfies the exterior ball condition; if the function f ∈ L ∞ (Ω), the solution of the fractional Poisson equation (2.6)-(2.7) satisfies u ∈ C α/2 (R n ), and furthermore u C α/2 (R n ) ≤ c f L ∞ (Ω) with c a constant depending on Ω and α.

Spectral fractional Laplacian
On a bounded domain Ω ⊂ R n , the spectral fractional Laplacian (also known as the fractional power of the Dirichlet Laplacian, or the "Navier" fractional Laplacian) is defined via the spectral decomposition of the standard Laplace operator [38,1,32], i.e., for α > 0, (2.8) where µ k and φ k are the k-th eigenvalue and normalized eigenfunction of the standard Dirichlet Laplace operator −∆ on the domain Ω. From a probabilistic point of view, it represents the infinitesimal generator of a subordinate killed Brownian motion, i.e., the process that first kills Brownian motion in a bounded domain Ω and then subordinates it via a α/2-stable subordinator [41,42]. Here, we include the domain Ω in the notation (−∆ Ω ) α/2 to reflect this process and to distinguish it from the fractional Laplacian (−∆) α/2 that was discussed in Section 2.1. Specially, if α = 2 the definition in (2.8) reduces to the standard Dirichlet Laplace operator −∆ on the domain Ω. The spectral fractional Laplacian is a nonlocal operator, and it is often used in the analysis of (partial) differential equations. The eigenvalues and eigenfunctions of the spectral fractional Laplacian are clearly suggested from its definition in (2.8), that is, the k-th eigenvalue of (−∆ Ω ) α/2 is µ α/2 k , and the corresponding eigenfunction is φ k (x). We remark that the spectral fractional Laplacian and the Dirichlet fractional Laplacian represent generators of different processes, which is also reflected by their eigenvalues and eigenfunctions. The eigenfunctions of the spectral fractional Laplacian are smooth up to the boundary as the boundary allows, while those of the Dirichlet fractional Laplacian are only Hölder continuous up to the boundary [38]. Additionally, it is easy to conclude from (2.4) that the k-th (for k ∈ N) eigenvalue of the Dirichlet fractional Laplacian is always smaller than that of the spectral fractional Laplacian.
The nonlocal Poisson problem with the spectral fractional Laplacian reads [16]: It can be viewed as an analogy to the fractional Poisson equation in (2.6)-(2.7). However, inheriting from the definition of the spectral fractional Laplacian, the boundary condition in (2.10) is defined locally on ∂Ω. From the definition in (2.8), one can formally obtain the solution of (2.9)-(2.10) as: where the coefficient f k is computed by

Regional fractional Laplacian
On a bounded domain Ω ⊂ R n , the regional fractional Laplacian (also known the censored fractional Laplacian) is defined as [4,23,24,22]: with the constant c n,α defined in (2.3). In contrast to the fractional Laplacian, the regional fractional Laplacian (−∆) α/2 Ω represents the infinitesimal generator of a censored α-stable process that is obtained from a symmetric α-stable Lévy process by restricting its measure to Ω. If the domain Ω = R n , the regional fractional Laplacian collapses to the fractional Laplacian (−∆) α/2 . To distinguish it from the fractional Laplacian (−∆) α/2 , we include the subscript 'Ω' in the operator (−∆) α/2 Ω to indicate the restriction of the α-stable Lévy process to the domain Ω. The regional fractional Laplacian is different from the Dirichlet fractional Laplacian, although they are freely interchanged in some literature. In fact, a symmetric α-stable Lévy process killed upon leaving the domain Ω (represented by the Dirichlet fractional Laplacian) is a subprocess of the censored α-stable process (represented by the regional fractional Laplacian) killing inside the domain Ω, i.e., the trajectories may be killed inside Ω through Feynman-Kac transform [42]. Moreover, we will illustrate their difference using a simple example. Consider a one-dimensional interval Ω = (−l, l). Let u be a smooth function satisfying u(x) = 0 for x ∈ Ω c . Then the difference between the regional fractional Laplacian and the Dirichlet fractional Laplacian can be computed as: We find that in the limiting case of α → 2, the difference between the regional fractional Laplacian and the Dirichlet fractional Laplacian vanishes, i.e., Q 1 → 0, due to the constant c 1,α → 0. In other words, the regional fractional Laplacian can be used to approximate the Dirichlet fractional Laplacian as α → 2. While in the limit of α → 0, the difference in (2.13) reduces to Q 1 u → u, i.e., the Dirichlet fractional Laplacian can be written as the summation of the regional fractional Laplacian and an identity operator. Otherwise, if α 0 and α 2, the difference Q 1 u ∼ O 1/(l − |x|) α , which does not tend to zero for any fixed l, and as x → ±l, there is |Q 1 u| → ∞.
In contrast to the fractional Laplacian, the current understanding of the regional fractional Laplacian still remains very limited. Recently, the interior regularity of the regional fractional Laplacian is discussed in [24,31]. It shows that (−∆) 1] or u ∈ C p+1, s (Ω) for s ∈ (α − 1, min(α, 1)]. So far, no results on the eigenvalues or eigenfunctions of the regional fractional Laplacian can be found in the literature. Here, we expect that our numerical results in Section 3.2 could provide insights into the understanding of the spectrum of the regional fractional Laplacian in the future.
The counterpart of the fractional Poisson equation (2.6)-(2.7), but with the regional fractional Laplacian reads: The above nonlocal problem has been used to approximate the fractional Poisson problem, in order to reduce computational complexity in solving (2.6)-(2.7). However, as shown in (2.13), the regional fractional model does not necessarily provide a consistent approximation to the fractional Poisson equation if α 2; see more numerical comparison and discussion in Section 3.3.

Peridynamic operator
The peridynamic models were originally proposed as a reformulation of the classical solid mechanics in [40]. In contrast to the classical models, it properly accounts for the near-field nonlocal interactions so as to effectively model elasticity problems with discontinuity and other singularities. The general form of this nonlocal operator has the following form: where B(x, δ), denoting a ball with its center at point x and radius δ, represents the interaction region of point x. The kernel function K(x, y) = K(|x − y|) describes the interaction strength between points x and y. The constant δ > 0 denotes the size of material horizon, and it is often chosen to be a small number in practical applications. Recently, the operator (2.16) with specially chosen kernel function is used to approximate the fractional Laplacian [13,21]. We will refer it as the peridynamic operator and denote it as i.e., the kernel function in this case is taken as: otherwise.
In other words, K δ (x, y) in the peridynamic operator represents a hard-threshold of the kernel function K(x, y) = c n,α /|x − y| n+α of the fractional Laplacian, which can be viewed as a truncation of K(x, y) in the fractional Laplacian. In the limiting case of δ → ∞, the peridynamic operator (2.17) coincides with the fractional Laplacian (2.2), and thus it is often used to approximate the fractional Laplacian by choosing a sufficiently large δ [13,21]. On the other hand, note that the kernel function K(x, y) has an algebraic decay of order n + α, which presents a heavy tail that accounts for considerable far field interactions. Hence, the cutoff of the kernel function K(x, y) outside of the horizon B(x, δ) may have a significant impact on its approximation to the fractional Laplacian as we shall show next. Similarly, we choose a smooth function u satisfying u(x) = 0 for x ∈ Ω c with Ω = (−l, l) to illustrate the difference between the peridynamic operator and the fractional Laplacian. Here, we assume that the horizon size δ in (2.17) is large enough, such that δ > max{l − x, l + x} for any point x ∈ (−l, l). Then, we can compute their difference as: It shows that the difference of these two operators is of order O(1/δ α ) when u(x) is uniformly bounded on Ω, hence their difference vanishes as δ → ∞. On the other hand, the convergence of the peridynamic operator to the fractional Laplacian as δ → ∞ depends on the power α, and it may degenerate rapidly for small α. Additionally, in the limiting case of α → 2, the difference Q 2 u → 0, because the coefficient c n,α → 0. A peridynamic model for describing the steady-state displacement of a finite microelastic bar can be formulated as follows [25,13,30]: where Ω δ defines the boundary region with width of δ. Note that the peridynamic model (2.19)-(2.20) is enclosed with a nonlocal Dirichlet boundary condition on a nonconventional finite "volume" boundary region of size δ, i.e. Ω δ . Hence, the boundary condition in (2.20) is also referred to as a finite volume constraint in the literature [15,13]. The peridynamic operator (−∆) α/2 δ is often used to approximate the fractional Laplacian (−∆) α/2 , while the nonlocal problem (2.19)- (2.20) is often used to approximate the fractional Poisson equation in (2.6)-(2.7). The well-posedness of (2.19)-(2.20) can be found in the literature [15].
The peridynamic operator in (2.17) can be viewed as an infinitesimal generator of a symmetric α-stable Lévy process by restricting its measure to B(x, δ). In contrast to the regional fractional Laplacian operator, the interaction region of point x in the peridynamic operator is symmetric with respect to itself. Hence, the peridynamic operator is expected to provide a symmetric approximation for a homogeneous elastic material.
In summary, the fractional Laplacian (2.2), spectral fractional Laplacian (2.8), regional fractional Laplacian (2.12), and the peridynamic operator (2.17) are all nonlocal operators in which every point x interacts with other points y over certain long distance. For a point x ∈ Ω, the fractional Laplacian (−∆) α/2 accounts for the interactions between x and y for all y ∈ R n \{x}. By contrast, the interaction region of x in the regional fractional Laplacian (−∆) α/2 Ω is truncated to Ω\{x}, i.e., the same domain of x, while the interaction region of the peridynamic operator (−∆) α/2 δ reduces to B(x, δ)\{x}. We will further compare them in Section 3.

Numerical comparisons
In this section, we further compare these four nonlocal operators by studying their nonlocal effects, eigenvalues and eigenfunctions, and the solution behavior of the corresponding nonlocal problems. In our simulations, the spectral fractional Laplacian is discretized by using the finite difference method combined with matrix transfer techniques introduced in [16], while the other three operators are discretized by the finite difference method based on weighted trapezoidal rules proposed in [17]. Our numerical results provide insights not only to further understand these operators but also to improve the analytical results in the literature.
In the following, we will consider the one-dimensional cases. For notational simplicity, we will also use L h to represent the fractional Laplacian, L s for the spectral fractional Laplacian, L r for the regional fractional Laplacian, and L p for the peridynamic operator.

Nonlocal effects of operators
We compare the nonlocal effects of these four operators by acting them on functions with compact support on the domain Ω = (−1, 1).

Example 1. Consider the function
which is continuous on the whole space R. It is easy to obtain that that is, the function from the spectral fractional Laplacian can be found exactly. While we will numerically compute the functions from the other three operators using the finite difference method proposed in [17].
In Figure 1, we compare the functions L i u for i = s, h, r, or p. The results clearly suggest the difference between these four operators, especially the function L s u from the spectral fractional Laplacian is significantly different from those of the other three operators. It shows that for any α ∈ (0, 2), the function L s u is proportional to the function u on (−1, 1). In contrast, the properties of L i u (for i = h, r, or p) significantly depend on the parameter α. For α ∈ (0, 1), the functions L i u (for i = h, r, or p) exist on the closed domain Ω, but they are very different between operators. The smaller the parameter α, the larger the differences. For α ∈ [1, 2), the functions L i u do not exist at the boundary points, i.e., x = ±1. As α → 2, the functions L i u (for i = h, r, or p) converge to −u xx for x ∈ (−1, 1). Additionally, Figure 1 shows that both the regional fractional Laplacian and peridynamic operator can be used to approximate the fractional Laplacian, if α is close to 2 (see Fig. 1 for α = 1.95). For small α, the results from the regional fractional Laplacian is inconsistent with that from the fractional Laplacian. However, the peridynamic operator can still provide a good approximation to the fractional Laplacian by enlarging the horizon size δ. Figure 2 presents the differences between the functions L p u and L h u for various α and δ. It shows that for a fixed horizon size δ, the difference between these two operators dramatically decreases as α increases. On the other hand, Figure 2 implies that for small α the convergence of the function L p u to L h u could be very slow. For instance, for α = 0.6, the difference in Fig. 2 is around 0.005 for a horizon size δ = 4000. In fact, the nonlocal interactions decay slowly for small α, and thus a large horizon size δ is needed for the peridynamic operator to better approximate the fractional Laplacian.

Example 2. Consider the function
for q ∈ N. For the fractional Laplacian, the analytical value can be found as: where 2 F 1 denotes the Gauss hypergeometric function. Moreover, we can obtain the exact values of (−∆) α/2 Ω u and (−∆) α/2 δ u by using their relation to the fractional Laplacian in (2.13) and (2.18), respectively. For the spectral fractional Laplacian, we numerically computed (−∆ Ω ) α/2 u by the finite difference method proposed in [16]. Figure 3 displays the functions L i u (for i = s, h, r, or p) for various α, where u is defined in (3.2) with q = 2. It shows that the functions L i u exist on the closed domain [−1, 1] for any α ∈ (0, 2), but their values are very different, especially for small α. For the spectral fractional Laplacian, the values of L s u are always zero at boundary points. For the regional fractional Laplacian, the function u in (3.2) with q = 2 satisfies the conditions that u ∈ C 2 ([−1, 1]) and u (±1) = 0, which guarantee the existence of the function L r u for any α ∈ (0, 2) [24]. Since the function u(±1) = 0 and the relations in (2.13) and (2.18), the values of L i u (for i = h, r and p) are the same at boundary points, but they are nonzero. Figure 3 also shows that both the regional fractional Laplacian and the peridynamic operator with relatively small δ could provide a good approximation to the fractional Laplacian, if α is large (see Fig. 3 for α = 1.95). While α is small, although the peridynamic operator can be still used to approximate the fractional Laplacian with a large δ, the regional fractional Laplacian is inconsistent with the fractional Laplacian. Figure 3 additionally shows that as α → 2, the differences between the four operators become insignificant (see Fig. 3 for α = 1.95), and the functions L i u (for i = h, s, r, or p) converge to −∂ xx u, that is, the four operators converge to the standard Dirichlet Laplace operator −∆. To understand their properties as α → 0, Figure 4 shows the functions L i u (for i = s, h, r, or p) for a small value of α = 0.001. The definition of the spectral fractional Laplacian in (2.8) implies that the function L s u converges to u, as α → 0. Fig. 4 shows that the function L h u from the fractional Laplacian converges to u as α → 0, confirming the analytical results in [14]. By contrast, the functions L r u from the regional fractional Laplacian and L p u from the peridynamic operator converge to a zero function. This can be easily obtained from their relation to the fractional Laplacian in (2.13) and (2.18), respectively. Moreover, our numerical results seem to suggest that for α ∈ [1, 2), if the function u ∈ C 1,α/2 (Ω) and u (±1) = 0, then the value L r u from the regional fractional Laplacian exists; see Figure 5. Hence, we conjecture that the regularity results in [24] might be able to improve to u ∈ C 1,α/2 (Ω) at least for one-dimensional case. More analysis needs to be carried out for further understanding of this issue.

Eigenvalues and eigenfunctions
In this section, we compare the four nonlocal operators by studying their eigenvalues and eigenfunctions on a one-dimensional bounded domain Ω = (−l, l). Denote λ i k and φ i k as the k-th (for k ∈ N) eigenvalue and eigenfunction of the nonlocal operator L i on Ω with the corresponding homogeneous Dirichlet boundary conditions, where i = h, s, r, or p. It is well known that the eigenvalues and eigenfunctions of the spectral fractional Laplacian L s  can be found analytically, i.e., for k ∈ N. For the other operators, so far no analytical results can be found in the literature, and thus we will compute their eigenvalues and eigenfunctions numerically.
In Table 1, we present the eigenvalues of the Dirichlet fractional Laplacian L h , spectral fractional Laplacian L s , and regional fractional Laplacian L r , on the domain Ω = (−1, 1). We leave the peridynamic operator L p out of our comparison here, since its spectrum depends on the horizon size δ. From Table 1 and our extensive numerical studies, we find that λ r k < λ h k < λ s k , for α ∈ (0, 2) and k ∈ N, that is, the eigenvalues of the regional fractional Laplacian are much smaller than those of the Dirichlet fractional Laplacian and the spectral fractional Laplacian. However, as α → 2 the eigenvalue λ i k of these three operators converges to µ k = k 2 π 2 /4 -the kth eigenvalue of the standard Dirichlet Laplace operator −∆ on (−1, 1).
In [38], it is proved that the first eigenvalue of the Dirichlet fractional Laplacian is strictly smaller than that of the spectral fractional Laplacian, i.e., λ h 1 < λ s 1 , for α ∈ (0, 2). Our numerical results in Table 1 confirm this conclusion and additionally suggest that the eigenvalue λ h k is strictly smaller than λ s k , for any k ∈ N. Furthermore, we present the difference between the eigenvalues λ s k and λ h k for various α and k in Figure 6. It shows that the difference between the eigenvalues λ s and λ h k depends on both parameters α and k. For a given k ∈ N, there exists a critical value α k,cr where the gap between λ s k and λ h k is maximized. The value of α k,cr increases as k ∈ N increases (see Fig. 6 left). On the other hand, as k → ∞ the relative difference between the eigenvalues λ s k and λ h k decreases quickly (see Fig. 6 right). In Figure 7, we compare the first and second eigenfunctions of the fractional Laplacian, the spectral fractional Laplacian, and the regional fractional Laplacian. For any α ∈ (0, 2), the eigenfunctions for these three operators are all symmetric (for odd k) or antisymmetric (for even k) with respect to the center of the domain Ω. Especially, the eigenfunctions of the spectral fractional Laplacian are independent of the parameter α, which are also the eigenfunctions of the standard Dirichlet Laplace operator −∆. In contrast, the eigenfunctions of the other two operators significantly depend on α, and as α → 2, they converge to sin(kπ(1 + x)/2) -the eigenfunctions of the standard Dirichlet Laplace operator −∆. Our numerical observations in Figure 7 justify the regularity results in [38,Theorem 1], that is, the eigenfunctions of the Dirichlet fractional Laplacian is no better than Hölder continuous up to the boundary, while the eigenfunctions of the spectral fractional Laplacian are smooth up to the boundary as the boundary allows.
From our extensive studies, we find that the eigenvalues of the Dirichlet fractional Laplacian L h , the spectral fractional Laplacian L s , and the regional fractional Laplacian L r reduce as the domain size increases. In particular, if the domain size increases by a ratio of κ, the k-th (for k ∈ N) eigenvalues decreases by a ratio of κ α . Additionally, we explore the eigenvalues of the peridynamic operators for different δ. It shows that the k-th (for k ∈ N) eigenvalue increases as δ decreases.

Solutions to nonlocal problems
In this section, we will further compare the four operators by studying the solutions of their corresponding nonlocal Poisson and diffusion problems. In the following, we choose the one-dimensional domain Ω = (−1, 1) and consider homogeneous Dirichlet boundary conditions, i.e., where Γ = Ω c for the fractional Laplacian, Γ = Ω δ for the peridynamic operator, and Γ = ∂Ω for the spectral fractional Laplacian and the regional fractional Laplacian with α ∈ (1, 2). Example 3. We consider a time-independent problem of the following form: subject to the homogeneous Dirichlet boundary conditions as discussed in (3.3), where i = h, s, r, or p. This problem is often used as a benchmark to test numerical methods for the fractional Laplacian and spectral fractional Laplacian [16,17,13]. In particular, the nonlocal problem (3.4) with the fractional Laplacian L h has diverse applications in various areas [19]. Its solution can be found exactly as: which represents the probability density function of the first exit time of the symmetric α-stable Lévy process from the domain Ω. While if the spectral fractional Laplacian L s is considered, the exact solution of (3.4) is expressed as We will compute the solutions of (3.4) with the regional fractional Laplacian L r or peridynamic operator L p numerically. Figure 8 illustrates the solutions of the nonlocal problem (3.4) with different operators L i . In the cases of the regional fractional Laplacian L r , since the solution does not exist for α ≤ 1, we only present those for α > 1. Generally, the solutions from the four operators are significantly The solution from the regional fractional Laplacian L r is very sensitive to the parameter α, and moreover it is inconsistent with that from the fractional Laplacian L h . The solution of the peridynamic model could serve a good approximation to that of the fractional Poisson problem (2.6)-(2.7), if a proper horizon size δ is chosen. The smaller the parameter α is, the larger the horizon size δ is needed, consistent with our observations in Section 3.1.
In Figure 9, we additionally study the solution of (3.4) with the peridynamic operator L p for various horizon size δ. It shows that the horizon size δ plays an important role in the solution of peridynamic models, especially when α is small (see Fig. 9 left). For the same α, the larger the horizon size δ, the smaller the solution u. Our extensive studies show that as δ → 0, the solution Example 4. Here, we study the following nonlocal diffusion problem: (3.6) subject to the homogeneous Dirichlet boundary conditions as discussed in (3.3), where i = h, s, r, or p. At time t = 0, the initial condition is taken as a step function, i.e., If the spectra of the nonlocal operator L i are known, one can formally express the solution of (3.6) as: , for x ∈ Ω, t ≥ 0, (3.8) where λ i k and φ i k represent the k-th eigenvalue and eigenfunction of the operator L i on the domain Ω, and the coefficient c k is calculated by Hence, the solution of the diffusion equation (3.6) with the spectral fractional Laplacian L s can be obtained analytically as for x ∈ Ω and t ≥ 0. For other cases with L i (i = h, r, p) in (3.6), we will compute their numerical solutions by using the finite difference method proposed in [17].  Figure 10 shows the time evolution of the solution u(x, t) for the cases of the spectral fractional Laplacian L s , fractional Laplacian L h , and regional fractional Laplacian L r , while the results for the peridynamic operator L p are displayed in Figure 11 for various horizon size δ. The similar phenomena to the classical diffusion equation are observed -the solution of nonlocal diffusion equation (3.6) decays over time, and the discontinuous initial condition smooths out quickly during the dynamics. For the same operator L i , the larger the parameter α is, the faster the solution diffuses. Moreover, Fig. 11 shows that for the peridynamic case, the diffusion speed of the solution depends not only on the parameter α, but also on the horizon size δ. For fixed α, the larger the horizon size, the faster the diffusion.
In Figure 12, we further compare the solutions at various time t, where the horizon size δ = 0.5 is used in the peridynamic operator. It shows that the solution quickly diffuses to the boundary if α is large (see Fig. 12 with α = 1.8). For the same parameter α, the solution resulting from the spectral fractional Laplacian diffuses much faster than those from other operators, which can be explained by the solution in (3.8). The solution (3.8) implies that the larger the eigenvalues λ i k , the faster the solution decays over time. On the other hand, Table 1 suggests that the eigenvalues of different operators satisfy λ s k > λ h k > λ r k , for α ∈ (0, 2) and k ∈ N. Hence, it is easy to conclude that the solutions from the spectral fractional Laplacian diffuse the fastest. As mentioned previously, the solution of the peridynamic model depends on the horizon size δ. Example 5. We consider a nonlocal diffusion-reaction equation of the following form: (3.10) subject to the homogeneous Dirichlet boundary conditions as discussed in (3.3), where i = h, s, r, or p. The initial condition is taken as which decays quickly to zero as |x| → ∞. Similarly, the solution of (3.10) can be formulated as: , for x ∈ Ω, t ≥ 0, (3.12) with λ i k and φ i k denoting the k-th eigenvalue and eigenfunction of the operator L i on the domain Ω, and the coefficient c k defined in (3.9). Figure 13 shows the time evolution of the solution u(x, t) to the reaction-diffusion equation (3.10) with the spectral fractional Laplacian L s , fractional Laplacian L h , or regional fractional Laplacian L r , while the results for the peridynamic operator L p are displayed in Figure 14. The results demonstrate the difference between these nonlocal operators, especially when α is small. The solution from the fractional Laplacian L h and the spectral fractional Laplacian L s continuously decay over time for all the α considered here. In contrast, the solution from the regional fractional Laplacian L r or the peridynamic operator L p may increase over time, depending on the competition between the diffusion and reaction terms. For example, for α = 0.3, the solutions increases constantly, but as α increases (e.g., α = 1.8), the diffusion becomes dominant, and the solution decays over time; see Figure 15.

Concluding remarks
We studied four nonlocal diffusion operators, including the fractional Laplacian, spectral fractional Laplacian, regional fractional Laplacian, and peridynamic operator. These four operators are equivalent on the entire space R n , but their differences on a bounded domain are significant. On a bounded domain Ω ∈ R n , they represent the generators of different stochastic processes. The Dirichlet fractional Laplacian represents the infinitesimal generator of a symmetric α-stable Lévy process that particles are killed upon leaving the domain; the regional fractional Laplacian is the generator of a censored α-stable process that is obtained from a symmetric α-stable Lévy process by restricting its measure to Ω, while the spectral fractional Laplacian is the generator of a subordinate killed Brownian motion (i.e. the process that first kills Brownian motion in a bounded domain Ω and then subordinates it via a α-stable subordinator). Our studies clarify the confusion existing in some literature, where the Dirichlet fractional Laplacian, spectral fractional Laplacian, and regional fractional Laplacian are freely interchanged.
We carried out extensive numerical investigations to understand and compare the nonlocal effects of these operators on a bounded domain. Our numerical results suggest that: ii) The eigenvalues and eigenfunctions of these four operators are different, although they all converge to those of the classical Laplace operator as α → 2. For each k ∈ N, the eigenvalues of the spectral fractional Laplacian are always larger than those of the fractional Laplacian and regional fractional Laplacian, which numerically extends the conclusion in the literature [38]. iii) For any α ∈ (0, 2), the peridynamic operator can provide a good approximate to the fractional Laplacian, if the horizon size δ is sufficiently large. We found that the solution of the peridynamic model converges to that of the fractional Laplacian model at a rate of O(δ −α ). Although the regional fractional Laplacian can be used to approximate the fractional Laplacian as α → 2, it generally provides inconsistent results from the fractional Laplacian for α 2.
Moreover, we provided some conjectures from our numerical results, which might contribute to the mathematics analysis on these operators. Due to their nonlocality, numerical simulations of problems involving these four operators are considerably challenging. We refer the readers to [16,17,13,15,21] for more discussions on numerical methods for these nonlocal models.