Fractional Laplace operator in two dimensions, approximating matrices, and related spectral analysis

In this work we review some proposals to define the fractional Laplace operator in two or more spatial variables and we provide their approximations using finite differences or the so-called Matrix Transfer Technique. We study the structure of the resulting large matrices from the spectral viewpoint. In particular, by considering the matrix-sequences involved, we analyze the extreme eigenvalues, we give estimates on conditioning, and we study the spectral distribution in the Weyl sense using the tools of the theory of Generalized Locally Toeplitz matrix-sequences. Furthermore, we give a concise description of the spectral properties when non-constant coefficients come into play. Several numerical experiments are reported and critically discussed.


Introduction
The study of the fractional Laplace operator was initiated in 1938 by Marcel Riesz in his seminal work [37]. However, since then numerous definitions for such operator have been proposed due to the need of describing several distinct types of applications (for a recent survey, refer to [24] and to the references therein). In the present paper, we review two main proposals to define the fractional Laplace operator which enable us to construct approximations in two or more spatial variables. As it is wellknown, the standard Laplacian in ℝ d can be formulated as For generalizing such operator towards the fractional case two different strategies can be followed: (a) compute a fractional power of the d-dimensional Laplacian in (1); (b) 'fractionalize' each integer order partial derivative in (1).
Obviously, these two strategies coincide in the 1D case, while they lead to different operators in more than one dimension. Formulation (a) is particularly useful when imposing boundary conditions which differ from the standard Dirichlet ones. Indeed, in [21] Ilić et al. showed that the problems involving the fractional Laplacian of the form (a) are well defined on finite domains with homogeneous boundary conditions of Dirichlet, Neumann and Robin type. This is due to the fact that the fractional Laplacian as in (a) has the same interpretation as the standard Laplacian in terms of spectral decomposition for such boundary conditions. Formulation (b) should be preferred instead when modeling phenomena that show a directiondependent character and then require a different fractional order along each spatial dimension. In this regard, we mention the application of the fractional paradigm to cardiac electrophysiology (see, e.g., [26]). Formulation (b) suits better than (a) here because in the heart the myocyte fibers tend to align in a particular direction and therefore the propagation of the electrical wave is different orthogonally, longitudinally, and transversally along the fibers.
Discretizing (1) by means of finite differences with uniform step-size in each direction and considering an hypercube as domain, we obtain a d-level Toeplitz matrix associated to a generating function [10], also known as symbol [17,18], which provides an asymptotic description of the spectrum of the discretization matrices. In the light of formulations (a) and (b), for generalizing the symbol to the fractional case we can: (1) take a fractional power of the symbol associated to the discretization of the d-dimensional Laplacian in (1); (2) 'fractionalize' the symbol associated to the discretization of each integer order partial derivative in (1).
Suitable (finite difference) approximations of the fractional Laplacian are then obtained by considering the corresponding d-level Toeplitz matrices (we refer the reader to [30,34,35] for the unidimensional case). It is worth mentioning that a similar idea was introduced by Lubich in [27,28] for constructing the so-called Fractional Linear Multistep Methods (FLMMs). In that case, the resulting schemes for a one-sided fractional operator (in the Riemann-Liouville or in the Caputo sense) were derived by considering a fractional power of the generating function of the coefficients of a classical linear multistep method for ordinary differential equations.
An alternative matrix representation of the fractional Laplacian can be obtained by means of the so-called Matrix Transfer Technique (MTT) proposed in [21,22]. The approximation to the fractional Laplacian in this case is given by a fractional power of the matrix representation obtained by using any reasonable method for discretizing the standard Laplacian. As an application, we mention that in [45] the MTT has been used for solving fractional Bloch-Torrey models in 2D with fractional Laplacian following both formulations (a) and (b). The interest in these models is due to the fact that they are used to study anomalous diffusion in the human brain. We point out that in order to face the high computational efforts that may result from this approach, rational approximations to the MTT and rational Krylov methods have been recently introduced in [1-4, 8, 20].
The aim of this paper is twofold: first, we review the spectral properties of finite difference discretizations, which are of Toeplitz type associated with symbols, as in (i) and (ii); second, we investigate the spectral properties of the matrices obtained by the MTT as fractional powers of both finite difference and finite element approximations of the standard Laplacian. Aside from the study of the spectrum/conditioning of the above matrices, we show that the corresponding matrix-sequences belong to the Generalized Locally Toeplitz (GLT) class.
One reason for characterizing all the aforementioned discretization matrixsequences as GLT matrix-sequences relates to the solution of the corresponding linear systems. Iterative schemes such as multigrid and preconditioned Krylov methods-able to exploit the GLT structure-can be found in [13,23,25,29,32].
The rest of this paper is organized as follows. In Sect. 2 we briefly recall the main results of the classical theory on the spectral behavior of Toeplitz matrices generated by integrable functions, as well as the notion of symbol for a generic matrixsequence and the main properties of the GLT class. In Sect. 3 we collect various definitions and results concerning the 1D fractional Laplacian. Section 4 is devoted to the construction of matrix approximations to the fractional Laplacian both in one and in two dimensions. Section 5 contains a detailed spectral analysis of the underlying matrices. Few numerical tests that confirm our theoretical findings are discussed in Sect. 6. Finally, in Sect. 7 we draw some conclusions.

Preliminaries on Toeplitz and GLT matrix-sequences
For an independent reading, in this section we formalize the definition of (d-level) Toeplitz matrix associated to a Lebesgue integrable function and we recall few key analytical features of the generating function, which provide information on the spectral properties of the associated Toeplitz matrix (Sect. 2.1). Moreover, we introduce the GLT class, a * -algebra of matrix-sequences containing Toeplitz sequences, diagonal sampling matrix-sequences, and zero-distributed matrix-sequences (Sect. 2.2).

Spectral results for Toeplitz matrices
Let f be an integrable function defined on [− , ] and extended periodically to the whole real axis ℝ. We denote by the Toeplitz matrix related to f (called generating funtion) whose coefficients a j , j = 0, ±1, ±2, … , ±(n − 1), are the Fourier coefficients of f, i.e., In the case where f is an even, real-valued and nonnegative function, the related T n (f ) is a real, symmetric and positive semidefinite matrix. If f is only real-valued, then the matrix T n (f ) is Hermitian.
The following result allows to localize the spectrum of T n (f ) in terms of ess inf f and ess sup f , the essential infimum and the essential supremum of f (that is the infimum and the supremum of f, up to zero-measure sets); see, e.g., [17] for a proof. More generally, when f and g are real-valued and g is nonnegative and not identically zero, the range of f/g gives precise indications regarding the spectrum of T −1 n (g)T n (f ). with r = ess inf(f ∕g) and R = ess sup(f ∕g). If g ≥ 0 (not identically zero) and r < R , then the spectrum of T −1 n (g)T n (f ) is contained in the interval (r, R) independently of n ∈ ℕ . If r = R then we have the trivial case where T −1 n (g)T n (f ) = rI n , with I n denoting the identity of order n.
To provide a result concerning the behavior of the minimum eigenvalue of T n (f ) the following definitions must be introduced first.
Definition 1 Let f, g be two nonnegative integrable functions (not identically zero) on ⊂ ℝ d , d ≥ 1 . We write f ∼ g if and only if there exist two positive constants 0 < r ≤ R < ∞ such that uniformly on .
Definition 2 Let f be a nonnegative integrable function on ⊂ ℝ d , d ≥ 1 . We say that f has a zero of order in =̄ if and only if there exist two positive constants 0 < c ≤ĉ < ∞ such that holds in a suitable neighborhood I ⊂ of ̄ . If I = , then ̄ is also the unique zero of f and in this case the definition is equivalent to the relation f ∼ ‖ −̄‖ 2 , according to Definition 1.
Using the notations introduced in Definitions 1 and 2, the following result holds (see [9,39]). The previous results can be extended to two-level Toeplitz matrices which are defined as follows. Let f ( 1 , 2 ) be an integrable function defined on the square [− , ] 2 and extended periodically on ℝ 2 . Denoting by r min (T n (g)) ≤ min (T n (f )) ≤ R min (T n (g)).
the Fourier coefficients of f, we define the two-level Toeplitz matrix T n,m (f ) generated by f as follows: with By the construction given in (2)-(4), we can easily deduce that when the generating function is a multiplicatively separable function, i.e., f ( 1 , 2 ) = s 1 ( 1 )s 2 ( 2 ) , we have with ⊗ being the standard Kronecker product [7]. The latter and the linearity of the Toeplitz operator imply that for a generic integrable function f where J r,1 ( J r,−1 ) denotes the lower (upper) shift matrix of order r. The last representation suggests the following compact notation in the case of a d-variate generating function and of a d-level Toeplitz matrix. Let f ( ) , , and the ordering considered in 1 − n ≤ j ≤ n − 1 is the lexicographical one employed for j and with s ≤ t if and , see [42].
The generating function f provides a description of the spectrum of T n (f ) , for n large enough in the sense of the following definition.
Definition 3 Let f ∶ ⊂ ℝ d → ℂ be a measurable function defined on a set , measurable in the sense of Lebesgue, with finite measure ̄( ). Moreover, let us assume that {A n } n is a sequence of matrices with eigenvalues j (A n ) and singular values j (A n ) , j = 1, … , d n , and that d n is the size of A n with d n → ∞ , as n → ∞ , i.e., n i → ∞, i = 1, … , d. • We say that {A n } n is distributed as f over in the sense of the eigenvalues, and we write {A n } n ∼ (f , ), if for every continuous function F with compact support. In this case, we say that f is the symbol of {A n } n . • We say that {A n } n is distributed as f over in the sense of the singular values, and we write {A n } n ∼ (f , ), if for every continuous function F with compact support. In particular, when f ≡ 0 we say that {A n } n is a zero-distributed matrix-sequence.
In the following, for the sake of simplicity we will often write {A n } n ∼ f , Remark 1 If f is smooth enough, an informal interpretation of the limit relation 6 (resp. 7) is that when n is sufficiently large, then the eigenvalues (resp. singular values) of A n can be approximated by a sampling of f (resp. |f|) on a uniform equispaced grid of the domain .

GLT class
The formal definition of GLT matrix-sequences is rather technical and involves somewhat cumbersome notation. Therefore, in this subsection we briefly discuss some properties of the GLT class which are sufficient for our aims and we refer the reader to [17] for more details. Throughout, we use the following notation to say that the matrix-sequence {A n } n is a GLT matrix-sequence with GLT symbol (x, ) , where x = (x 1 , … , x d ) are the physical variables and = ( 1 , … , d ) are the Fourier variables. First we introduce a certain type of diagonal matrices which, as we shall see, come into play when dealing with differential problems with variable coefficients.
We list below key features of GLT matrix-sequences which represent a characterization of GLT matrix-sequences. In other words the * -algebra of matrix-sequences satisfying items GLT1-GLT6 coincides with the whole class of GLT matrixsequences. We have chosen this way, not only for the sake of brevity and notational simplicity, but also because items GLT1-GLT6 furnish a constructive way for verifying whether a concrete matrix-sequence belongs to the GLT class.
then it also holds that {A n } n ∼ ( , G). GLT2 The set of GLT matrix-sequences forms a * -algebra that is closed under linear combinations, products, inversion, conjugation. In formulas, let {A n } n ∼ GLT 1 and {B n } n ∼ GLT 2 , then GLT4 Every diagonal sampling matrix-sequence {D n ( )} n is a GLT matrixsequence whose symbol is given by (x, ) = (x). GLT5 Every zero-distributed matrix-sequence is a GLT matrix-sequence with symbol 0 and viceversa, i.e., According to Definition 3, in the presence of a zero-distributed matrix-sequence the singular values of the n th matrix (weakly) cluster around 0. This is formalized in the following result.

The 1D fractional Laplacian
In this section we review two different definitions of the fractional Laplacian in one dimension and we report their basic properties useful in the sequel. We shall always assume that the fractional Laplacian is applied to a suitable smooth function, say v (for instance, v ∈ L 1 (ℝ) ). In addition, for simplicity, we denote the space variable by x.
Consider the Riesz potential operator that in the 1D case can be formulated as where denotes the Gamma function. It verifies the following properties: Such operator allows to introduce the first definition of fractional Laplacian to which we refer, also known as Riesz fractional derivative operator, cf. [11,38]. In [11,Lemma 2.2], under the hypotheses that v ∈ C 5 (ℝ) and all derivatives up to order five belong to L 1 (ℝ), Çelik and Duman proved that where c 1 denotes a generalization of the integer-order centered differences to the fractional case proposed by Ortigueira in [30, Definition 3.5], namely Taking the limit as h approaches 0 of both sides of (9) gives another form for the fractional Laplacian (see (8)) which in [30, Definition 4.1] was called type 1 fractional centred derivative. Therefore, setting in (10) we have It can be easily verified that the coefficients k satisfy the following properties [11, Property 2.1]: where z = + m, m < z < m + 1, ∈ (0, 1), for all nonnegative values of k we obtain Consequently, for k ≥ 0 we can write (see [46]) or (10) By using the ratio expansion of two Gamma functions at infinity (see, e.g., [38, p. 17]) one deduces For this reason, the series in (10) converges absolutely for a bounded function v ∈ L 1 (ℝ). By exploiting the decay of the coefficients k , a "short memory" version of the scheme proposed by Ortigueira has been implemented in [36]. Finally, using the Fourier definition of the fractional Laplacian, in [31] was proved that In this case, the Fourier coefficients are given by We conclude this section by pointing out that several applications require the solution of fractional differential equations expressed in terms of fractional Laplacian multiplied by a certain non-constant coefficient. In such cases, the spatial operator has the following form with D a non-negative function. Later we also give spectral insights on this operator under the additional hypothesis that D is Riemann integrable.

Matrix forms for Laplace operators
Aiming to discretize the fractional Laplace operator in two (or more) dimensions, in the following we always impose absorbing boundary conditions (see [12] for further details). In the one-dimensional case this is equivalent to applying the fractional Laplace operator to functions of the form Different boundary conditions can also be imposed and the spectral analysis in terms of distribution is quite simplified using the GLT machinery. In fact, the perturbation .
induced by a different choice of boundary conditions can be represented as a zerodistributed matrix-sequence and hence GLT5 combined with items GLT1-GLT2 represents the key for the related spectral analysis. We finally emphasize that from a matrix viewpoint the method in Sect. 4.1 amounts in taking the Toeplitz structure generated by a proper 'fractionalization' of the symbol of the discrete Laplacian, while the method in Sect. 4.2 amounts in taking the same fractional power applied to a (possibly Toeplitz) structure. It is evident that the first method is computationally much more convenient: the interesting fact is that the two resulting GLT sequences possess the same symbol (so that they differ for a zero-distributed matrix-sequence), even if they look quite different in terms of entries and of global algebraic structure (for instance in the second form the Toeplitz structure is completely lost).

Toeplitz matrix form
Consider the uniform mesh of the domain [a, b] with step-size h = (b − a)∕(n + 1), i.e., From (11), when = 2, it is clear that the standard 1D Laplacian is expressed in terms of centered differences as where we abbreviate w j = w( j ). Setting all these formulas can be written simultaneously in matrix form as −h −2 T n (f 2 ).
We remark that T n (f 2 ) is the Toeplitz matrix whose entries on the kth diagonal are the Fourier coefficients of the generating function Raising f 2 to the appropriate power yields Therefore, from (12) we deduce that k are the Fourier coefficients of f . Denoting by (14) j = a + jh, j = 0, 1, … , n + 1.
the Toeplitz matrix related to f , and taking into account (11) it is a simple matter to verify that −h − T n (f ) represents the matrix form for the 1D fractional Laplacian described in Sect. 3, cf. [11,46]. Now we turn to the Laplace operator in two dimensions. Discretizing in space with a uniform mesh of step-size h in each domain direction and using the 5-point finite-difference stencil leads to the block tridiagonal matrix having the following form: with I n denoting the identity matrix of size n and T n (f 2 ) defined according to (15). T n,n ( 2 ) is a 2-level Toeplitz matrix which corresponds to the bivariate symbol (see (2)) Bearing in mind the idea presented above for approximating the 1D fractional Laplacian, in the two-dimensional case we may approximate the fractional Laplacian by considering the matrix T n,n ( ) related to the generating function From (16), we can immediately rewrite 2 ( 1 , 2 ) = f 2 ( 1 ) + f 2 ( 2 ). Consequently, Another way for constructing an approximation of the 2D fractional Laplacian is based on the 2-level Toeplitz matrix generated by the separable function hereafter we denote this matrix by T n,n ( ). Using (5), we obtain that where T n (f ) is given by (17). Obviously, T n,n ( 2 ) = T n,n ( 2 ).

MTT matrix form
As already mentioned in the introduction, the Matrix Transfer Technique (MTT) is a strategy for discretizing the fractional Laplace operator. Originally introduced in [21,22] for the 1D case, the MTT has been also generalized to the two-dimensional case in ( 1 , 2 ) = 4 − 2 cos 1 − 2 cos 2 ∕2 .
(20) T n,n ( ) = T n (f ( 1 )) ⊗ I n + I n ⊗ T n (f ( 2 )), [44]. The approximation provided by this method is given by a fractional power of any reasonable discretization of the standard Laplacian. Here we consider two particular discretizations in the one-dimensional case: the one obtained by using centered finite differences and the one corresponding to a generic finite element approximation. We denote by L FD ,n and L FE ,n , the corresponding discrete fractional operators provided by the MTT. In formulas (see Sect. 4.1) where M n and K n are the usual finite element mass and stiffness matrices. In particular, M n = hT n ((2 + cos )∕3) and K n = h −1 T n (f 2 ) when piecewise linear basis functions are adopted.
Concerning the 2D case, in this work we will focus our spectral analysis on the following MTT matrix forms: where M n,n and K n,n denote the mass and stiffness matrices in the 2D case.

Matrix form for the case of non-constant coefficient
Once we have a matrix representation for the 1D fractional Laplacian, say L ,n , obtained as explained in the last two subsections (see (17), (21) and (22)), we can easily recover a discretization for the operator (13) on the uniform mesh (14) just taking where D n = diag j=1,…,n (D( j )) . Similar arguments apply to the d-dimensional operator In particular, a discretization of this operator when d = 2 takes the following form where D n,n is the diagonal matrix whose non-zero entries are the values of D(x) at the grid points and L ,n,n denotes a matrix form of the 2D fractional Laplacian.

Spectral analysis
In this section, divided into three parts, the spectral analysis is performed following the results reported in Sect. 2 concerning the Toeplitz technology and the quite powerful GLT theory. The interesting message we can give is the following: independently of the used method, despite the substantial difference in the structure of the involved matrices, the smallest eigenvalue scales like n − and all the related matrix-sequences present the same global eigenvalue distribution.

Toeplitz case
First we focus our attention on the spectral properties of T n (f ). Since its generating function f is an even, real-valued and nonnegative function, the matrix T n (f ) is real, symmetric and positive semidefinite. Consequently, where n, = diag( 1, , 2, , … , n, ) is the diagonal matrix containing the eigenvalues of T n (f ) sorted in nondecreasing order and Q n, is the orthogonal matrix whose columns are the eigenvectors associated to the eigenvalues j, , j = 1, 2, … , n. In addition, considering that ess inf f = 0, ess sup f > 0 and f ∼ | | , from Theorem 3 we deduce that • T n (f ) is positive definite for each n ∈ ℕ; • the minimum eigenvalue of T n (f ) behaves as follows: 1, ∼ 1∕n . Now, using (20) and (24) we obtain that The matrix ( n, ⊗ I n + I n ⊗ n, ) is block diagonal with diagonal blocks, say D j , having the following form: In particular, this implies that the minimum and the maximum eigenvalues of T n,n ( ) are, respectively, It is worth to observe that since 1, and n, are the minimum and the maximum eigenvalues of T n (f ), the condition number in Euclidean norm of the matrix T n,n ( ) coincides with the one of the matrix T n (f ), i.e., In addition, the behavior of 1, leads to deduce that (24) T n (f ) = Q n, n, Q T n, , T n,n ( ) = (Q n, n, Q T n, ) ⊗ I n + I n ⊗ (Q n, n, Q T n, ) = (Q n, ⊗ Q n, )( n, ⊗ I n + I n ⊗ n, )(Q n, ⊗ Q n, ) T .
By virtue of the fact that f 2 ∼ | | 2 (see (16)), the two generating functions defined in (18) and (19) satisfy, respectively, Using the polar coordinates ( , ) we can write Therefore, from the previous relations we have Since (|cos | + |sin | ) is never identically zero, we conclude that By applying the block version of Theorem 4 (see Theorem 3.9 in [40] for a more general result), we have Consequently, Notice that because of (25) and (26), both and have a zero at the origin of order in the sense of Definition 2.

MTT case
Finite Differences It is well-known that the discretization T n (f 2 ) given in (15) can be diagonalized by means of the discrete sine transform. In formulas, T n (f 2 ) = Q n,2 n,2 Q n,2 , where the entries of the symmetric orthogonal matrix Q n,2 are given by and, for j = 1, … , n, the non-zero entries of the diagonal matrix n,2 are the eignevalues j,2 = 2 − 2 cos(j ∕(n + 1)). This enables to write directly the spectral decomposition of L FD ,n by means of Q n,2 ∕2 n,2 Q n,2 and to conclude that min (T n,n ( )) ∼ 1∕n .
Similarly, in two dimensions the matrix T n,n ( 2 ) = T n (f 2 ) ⊗ I n + I n ⊗ T n (f 2 ) can be diagonalized by means of Q n,2 ⊗ Q n,2 , and its diagonal form is I n ⊗ n,2 + n,2 ⊗ I n . Therefore, its eigenvalues are all the possible sums i,2 + j,2 , with i, j = 1, … , n and we can compute L FD ,n,n by means of the decomposition (Q n,2 ⊗ Q n,2 )(I n ⊗ n,2 + n,2 ⊗ I n ) ∕2 (Q n,2 ⊗ Q n,2 ) which implies that and 2 (L FD ,n ) = 2 (L FD ,n,n ). GLT considerations Since T n (f 2 ) is a symmetric positive definite matrix and ,n is still real symmetric and by GLT1 we also have {−h L FD ,n } n ∼ f . In other words, {−h L FD ,n } n is a GLT matrix-sequence that shares the same symbol of the Toeplitz sequence {T n (f )} n obtained by a direct centered finite difference discretization of the fractional Laplacian. As a consequence, by GLT2 which thanks to GLT5 and Proposition 1 is equivalent to say that the matrixsequence {−h L FD ,n − T n (f )} n is weakly clustered at zero. A similar reasoning proves that {−h L FD ,n,n } n ∼ GLT, and that {−h L FD ,n,n − T n,n ( )} n ∼ GLT, 0. Finite Elements The matrix representation of the 1D standard Laplacian resulting from the finite element method is given by C n = M −1 n K n . Both M n and K n are real and symmetric and M n is positive definite. Considering that C n = M −1∕2 n K n M −1∕2 n is a real symmetric matrix and it follows that C n is diagonalizable. In particular, when piecewise linear basis functions are adopted, that is M n = hT n ((2 + cos )∕3) and K n = h −1 T n (f 2 ), then C n can be still diagonalized by means of Q n,2 and the eigenvalues of both C n and L FE ,n are known exactly. Specifically, we have Similar considerations can be made in the two-dimensional setting. The general case of high-order finite elements is less trivial. For example, in [14] it has been shown that the matrices M n and K n do not belong to the same matrix algebra when considering finite elements with maximum regularity of order greater than or equal to 3. In this context, an approximation of the eigenvalues of L FE ,n and L FE ,n,n can be obtained using the notion of symbol, as outlined in the next paragraph.
GLT considerations Here we discuss the case of high-order finite elements with maximum regularity and we exploit the results in [15] where it has been shown that {h −1 M n } n and {hK n } n are both GLT matrix-sequences. We denote the corresponding symbols as m and , respectively. Since both M n and K n are real symmetric, by GLT1 the matrix-sequences {h −1 M n } n and {hK n } n distribute as m, also in the eigenvalue sense. In formulas, Since h −1 M n is positive definite with minimal eigenvalue well separated from zero (see again [15]), the corresponding function m is invertible everywhere. For instance, when using piecewise linear elements, h −1 M n = T n ((2 + cos )∕3) and hence on the whole definition domain. Consequently, item GLT2 implies that {h 2 C n } n ∼ GLT m −1 =∶ . In addition, from (27) and by recalling that C n is real symmetric, by using items GLT1, GLT2, and GLT6 we get Furthermore, since K n is positive semidefinite, also C n is positive semidefinite and then (C n ) ∕2 is well-defined. Taking into account that (C n ) ∕2 is still real symmetric, by GLT6 we infer that To prove that also {−h L FE ,n } n is a GLT matrix-sequence we first observe that (see (22) and (27)) The last equality in (30) comes directly from the definition of matrix function and from the fact that C n is symmetric positive semidefinite. More precisely, if we write the diagonal form of C n as C n = ZJZ T , then Starting from (30) and putting together (28), (29), and by using items GLT1, GLT2, and GLT6, we deduce that Notice that the above reasoning still works in the two-dimensional setting. When replacing M n and K n with M n,n and K n,n we have with = m −1 where now m, are the symbols of the properly scaled two-dimensional mass and stiffness matrices.
As an example, we consider the finite element method constructed by using C 2 piecewise cubic basis functions. Since and in this case while Remark 2 In the case of high-order finite element with a certain fixed regularity the resulting matrix-sequences are (sequences of submatrices of) block GLTs. Analogous results to (31) and (32) follow using the multilevel block GLT theory [5,6,16].

Non-constant coefficient case
In this section we review how the GLT theory can be used for studying the spectral behavior of the matrix-sequence {−h M ,n } n , with M ,n given by (23). In particular, we focus on the case of the finite differences which has already been analyzed in [13,Proposition 5] and in [29,Proposition 3.11] and on the finite element one. Everything works in the two-dimensional setting as well, so we skip the details for this case.
Let h (x, ) = D(x)g ( ) be the symbol associated to −h M ,n = −h D n L ,n , where L ,n is one of the discretizations of the 1D fractional Laplacian discussed in the previous sections whose symbol is denoted by g . Noticing that D n is a positive definite diagonal matrix, then that is M ,n is similar to a real symmetric matrix. Therefore, all the eigenvalues of   a)x, −h M ,n n ∼ĥ (x, ) is equivalent to −h M ,n n ∼ h (x, ). Therefore, From the limit relation which is true simply because it is clear that the symbol h has a zero of order at = 0 , whenever the diffusion coefficient D is bounded and positive.

Summary of spectral results
We now summarize all the spectral results obtained by collecting them in two groups, those relating to the one-dimensional case and those relating to the two-dimensional case. Furthermore, we assign to each operator an acronymous that will be useful in the next section.

Numerical results
In this section we numerically check some of the relations we have listed at the end of the previous section. In all our experiments we used the Matlab function eig to compute the eigenvalues. In addition, when finite elements are involved, we opt for C 2 piecewise cubic basis functions. We start by verifying the spectral relations 1FD-T and 1FE-MTT. Since the symbols involved are even functions, we can restrict to the interval [0, ] and consider a uniform mesh with step-size h = ∕(n + 1). As shown in Fig. 1, setting = 1.3 and n = 100, there is a good matching between the eigenvalues of the matrices T n (f ) and −h L FE ,n , and a sampling of f on the points grid. Similar results can be obtained when checking 1FD-NC with g = f , = 1.2, n = 100, and D(x) = (3 − )x or D(x) = exp(8x) (refer to Fig. 2).
Consider now the square domain [0, ] 2 and define a uniform mesh with step-size h = ∕(n + 1) in each domain direction. In Fig. 3 we check the validity of 2FD-T-S and 2FE-MTT comparing the spectrum of the matrices T n,n ( ) , −h L FE ,n,n with a sampling of , ∕2 , when = 1.7 and n = 10. Also in this two-dimensional example the symbol is describing quite well the spectrum of the considered discretization matrices.
As a final test, in Fig. 4 we compare the contour plots of all the symbols corresponding to two-dimensional discretization matrices, i.e., , ∕2 , and , when = 1.4 . What can be immediately noticed is that the symbol shows sharper patterns than and especially ∕2 , when approaching zero. Moreover, we note that the 2FE-MTT symbol is close to zero on a slightly larger portion of the domain than the other two symbols. This results in a worse ill-conditioning of the corresponding discretization matrices and refers to the degree of the approximation space.

Conclusions
Starting from two different ways of defining the fractional Laplace operator in more than one dimension we have studied the spectral properties of related finite difference and finite element discretization matrices. Precisely, we have analyzed the spectrum/conditioning of • the Toeplitz matrices whose symbol is a proper fraction of the symbol of the standard Laplacian operator; • the matrices obtained by the MTT as fractional powers of both finite difference and finite element approximations of the standard Laplacian.
Furthermore, we have shown that the corresponding matrix-sequences belong to the GLT class. The characterization of all the aforementioned discretization matrixsequences as GLT matrix-sequences relates to the solution of the corresponding linear systems. Inspired by the GLT-based algorithmic proposals in [29] for finite difference discretizations of the 2D fractional Laplacian and due to the ill-conditioning observed when considering the 2FE-MTT symbol, in a future work it could be interesting to exploit the GLT structure for the definition of preconditioners suitable for finite element discretizations.