In this paper, we propose accurate and efficient finite difference methods to discretize the twoand three-dimensional fractional Laplacian (−∆) α 2 (0 < α < 2) in hypersingular integral form. The proposed finite difference methods provide a fractional analogue of the central difference schemes to the fractional Laplacian, and as α → 2 − , they collapse to the central difference schemes of the classical Laplace operator −∆. We prove that our methods are consistent if u ∈ C α , α− α +ε (R d ), and the local truncation error is O(h ε ), with ε > 0 a small constant and · denoting the floor function. If u ∈ C 2+ α , α− α +ε (R d ), they can achieve the second order of accuracy for any α ∈ (0, 2). These results hold for any dimension d ≥ 1 and thus improve the existing error estimates for the finite difference method of the one-dimensional fractional Laplacian. Extensive numerical experiments are provided and confirm our analytical results. We then apply our method to solve the fractional Poisson problems and the fractional Allen-Cahn equations. Numerical simulations suggest that to achieve the second order of accuracy, the solution of the fractional Poisson problem should at most satisfy u ∈ C 1,1 (R d ). One merit of our methods is that they yield a multilevel Toeplitz stiffness matrix, an appealing property for the development of fast algorithms via the fast Fourier transform (FFT). Our studies of the two-and three-dimensional fractional Allen-Cahn equations demonstrate the efficiency of our methods in solving the high-dimensional fractional problems. In classical partial differential equation (PDE) models, diffusion is described by the classical Laplace operator ∆ = (∂ xx + ∂ yy + ∂ zz ), characterizing the transport mechanics due to Brownian motion. Recently, it has been suggested that many complex (e.g., biological and chemical) systems are instead characterized by Lévy motion, rather than Brownian motion; see [9,17,24,8] and references therein. Hence, the classical PDE models fail to properly describe the phenomena in these systems. To circumvent such issues, the fractional models were proposed, where the classical Laplace operator is replaced by the fractional Laplacian (−∆) α 2 [24,9]. In contrast to the classical diffusion models, the fractional models possess significant advantages for describing problems with long-range interactions, enabling one to describe the power law invasion profiles that have been observed in many applications [7,25,29]. However, the nonlocality of the fractional Laplacian introduces considerable challenges in both mathematical analysis and numerical simulations.
