Summary: In this paper, we propose accurate and efficient finite difference methods to discretize the two- and three-dimensional fractional Laplacian \((- \Delta)^{\frac{\alpha}{2}}\) (\(0 < \alpha < 2\)) in hypersingular integral form. The proposed methods provide a fractional analogue of the central difference schemes to the fractional Laplacian. As \(\alpha \rightarrow 2^-\), they collapse to the central difference schemes of the classical Laplace operator \(- \Delta \). We prove that our methods are consistent if \(u \in C^{\lfloor \alpha \rfloor, \alpha - \lfloor \alpha \rfloor + \varepsilon}(\mathbb{R}^d)\), and the local truncation error is \(\mathcal{O}(h^\varepsilon)\), with \(\varepsilon > 0\) a small constant and \(\lfloor \cdot \rfloor\) denoting the floor function. If \(u \in C^{2 + \lfloor \alpha \rfloor, \alpha - \lfloor \alpha \rfloor + \varepsilon}(\mathbb{R}^d)\), they can achieve the second order of accuracy for any \(\alpha \in(0, 2)\). These results hold for any dimension \(d \geq 1\) and thus improve the existing error estimates of the one-dimensional cases in the literature. 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 \in C^{1, 1}(\mathbb{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.