Integral Fractional Laplacian Overview

Updated 2 February 2026

Integral Fractional Laplacian is a nonlocal operator defined as a hypersingular integral with a rotationally invariant kernel, modeling phenomena like anomalous diffusion.
It operates within specialized Sobolev spaces with nonlocal exterior conditions and serves as the generator of symmetric stable Lévy processes.
Numerical methods such as FEM, FD, and Monte Carlo leverage its matrix structures to efficiently handle singularities and boundary challenges in nonlocal PDEs.

The integral fractional Laplacian—sometimes called the Riesz or regional fractional Laplacian—is a nonlocal operator arising in analysis, probability, and mathematical physics, defined as a hypersingular integral with a rotationally invariant kernel. On bounded domains, this operator is supplementally prescribed with nonlocal exterior conditions and serves as the generator for symmetric stable Lévy processes. Its mathematical properties, regularity, numerical discretization, and application have been the subject of extensive research in PDE theory, computational mathematics, and statistical physics.

1. Mathematical Definition and Normalization

The integral fractional Laplacian of order $s \in (0,1)$ for a sufficiently smooth function $u:\mathbb{R}^d\to\mathbb{R}$ is given by the hypersingular integral

$(-\Delta)^s u(x) = C_{d,s}\;\mathrm{P.V.}\!\int_{\mathbb{R}^d} \frac{u(x)-u(y)}{|x-y|^{d+2s}}\,dy,$

where “P.V.” denotes the Cauchy principal value. The normalization constant

$C_{d,s} = \frac{2^{2s} s\,\Gamma(s + d/2)}{\pi^{d/2} \Gamma(1-s)}$

ensures that the operator matches its Fourier-multiplier characterization: $\widehat{(-\Delta)^s u}(\xi) = |\xi|^{2s}\widehat{u}(\xi).$ The operator extends to bounded Lipschitz domains $\Omega$ by prescribing $u=0$ on $\mathbb{R}^d \setminus \Omega$ and restricting the integral to functions supported in $\Omega$ (Lischke et al., 2018, Bonito et al., 2017).

2. Functional Framework and PDEs

On domains $\Omega$ , the associated Sobolev space is the zero-extension space: $u:\mathbb{R}^d\to\mathbb{R}$ 0 with seminorm induced by the energy bilinear form

$u:\mathbb{R}^d\to\mathbb{R}$ 1

The weak Dirichlet problem is: find $u:\mathbb{R}^d\to\mathbb{R}$ 2 such that

$u:\mathbb{R}^d\to\mathbb{R}$ 3

Well-posedness follows from Lax–Milgram. The operator also admits a probabilistic interpretation as generator of the killed symmetric $u:\mathbb{R}^d\to\mathbb{R}$ 4-stable Lévy process, and solutions admit the Feynman–Kac formula (Lischke et al., 2018, Sheng et al., 2022).

3. Regularity and Weighted Estimates

Even with smooth data, solutions of the integral fractional Laplacian exhibit boundary singularities; for example, $u:\mathbb{R}^d\to\mathbb{R}$ 5 near $u:\mathbb{R}^d\to\mathbb{R}$ 6, so $u:\mathbb{R}^d\to\mathbb{R}$ 7 is not in $u:\mathbb{R}^d\to\mathbb{R}$ 8 generally (Lischke et al., 2018, Borthagaray et al., 2021). Weighted Sobolev and Hölder regularity results for solutions are available:

For $u:\mathbb{R}^d\to\mathbb{R}$ 9, $(-\Delta)^s u(x) = C_{d,s}\;\mathrm{P.V.}\!\int_{\mathbb{R}^d} \frac{u(x)-u(y)}{|x-y|^{d+2s}}\,dy,$ 0 for $(-\Delta)^s u(x) = C_{d,s}\;\mathrm{P.V.}\!\int_{\mathbb{R}^d} \frac{u(x)-u(y)}{|x-y|^{d+2s}}\,dy,$ 1 and any $(-\Delta)^s u(x) = C_{d,s}\;\mathrm{P.V.}\!\int_{\mathbb{R}^d} \frac{u(x)-u(y)}{|x-y|^{d+2s}}\,dy,$ 2 (Borthagaray et al., 2018).
Analytic regularity up to boundary vertices/edges: in polygons and polyhedra, $(-\Delta)^s u(x) = C_{d,s}\;\mathrm{P.V.}\!\int_{\mathbb{R}^d} \frac{u(x)-u(y)}{|x-y|^{d+2s}}\,dy,$ 3 and its Caffarelli–Silvestre extension admit weighted $(-\Delta)^s u(x) = C_{d,s}\;\mathrm{P.V.}\!\int_{\mathbb{R}^d} \frac{u(x)-u(y)}{|x-y|^{d+2s}}\,dy,$ 4-estimates with factorial growth, precisely characterizing singular layers via vertex-edge-face weights and bootstrapping of Caccioppoli inequalities (Faustmann et al., 2021, Faustmann et al., 2023).
For piecewise regularity $(-\Delta)^s u(x) = C_{d,s}\;\mathrm{P.V.}\!\int_{\mathbb{R}^d} \frac{u(x)-u(y)}{|x-y|^{d+2s}}\,dy,$ 5, discretization errors saturate at $(-\Delta)^s u(x) = C_{d,s}\;\mathrm{P.V.}\!\int_{\mathbb{R}^d} \frac{u(x)-u(y)}{|x-y|^{d+2s}}\,dy,$ 6 or higher (Wu et al., 2021).

4. Numerical Discretization Methods

Finite Element Methods (FEM):

Direct assembly of the dense stiffness matrix (double integrals over $(-\Delta)^s u(x) = C_{d,s}\;\mathrm{P.V.}\!\int_{\mathbb{R}^d} \frac{u(x)-u(y)}{|x-y|^{d+2s}}\,dy,$ 7) with explicit quadrature handling near and far-field singularities, panel clustering, and hierarchical matrices for efficient matvec (Ainsworth et al., 2017).
Weighted-residual a posteriori error estimators and adaptive mesh refinement, with skeleton-weighted indicators for $(-\Delta)^s u(x) = C_{d,s}\;\mathrm{P.V.}\!\int_{\mathbb{R}^d} \frac{u(x)-u(y)}{|x-y|^{d+2s}}\,dy,$ 8 (Faustmann et al., 2019, Faustmann et al., 2022).
Greedy mesh grading to match local regularity for quasi-optimal convergence (Borthagaray et al., 2021, Borthagaray et al., 2018).

Finite Difference Methods (FD):

Central difference analogues, operator factorization, and quadrature splitting schemes to obtain block-Toeplitz structure (Wu et al., 2021, Duo et al., 2018, Hao et al., 2024, Huang et al., 2013).
Singularity subtraction to regularize integrands for high-order schemes (Minden et al., 2018).
Fast implicit schemes for time-space fractional diffusion using Toeplitz matrices, M-matrix theory, SOE convolution acceleration, and circulant preconditioning (Gu et al., 2020).

Meshless and RBF Methods:

Meshless collocation based on order-dependent generalized multiquadric RBFs, leveraging analytic pseudo-spectral identities to avoid hypersingular integration, and using low-rank tail correction with rapid convergence (Hao et al., 2024).

Monte Carlo, Walk-on-Spheres:

Monte Carlo algorithms based on spatial Green’s function and Poisson kernel for the fractional Laplacian, extending the classical walk-on-spheres to multiple dimensions and arbitrary domains, with rigorous error and step-count estimates (Sheng et al., 2022).

Discontinuous Galerkin (DG) Methods:

Local DG and LDG schemes using mixed formulations with nonlocal fluxes, stabilization on graded meshes, and analysis of Riesz potentials (Nie et al., 2021, Han et al., 13 Dec 2025).

5. Matrix Structure and Fast Solvers

Discretization schemes exploit the translation-invariance and Toeplitz/block-Toeplitz structure of the discrete operator when possible, both in FD, RBF, and factorization approaches (Duo et al., 2018, Wu et al., 2021, Gu et al., 2020, Minden et al., 2018). This enables $(-\Delta)^s u(x) = C_{d,s}\;\mathrm{P.V.}\!\int_{\mathbb{R}^d} \frac{u(x)-u(y)}{|x-y|^{d+2s}}\,dy,$ 9 applications via FFTs and efficient iterative solutions (CG, PCG). Panel clustering and hierarchical matrices in FEM and BEM settings further reduce complexity for global matvec (Ainsworth et al., 2017).

Method	Matrix Structure	Solver Complexity
FD (Toeplitz, circulant)	Block-Toeplitz/Toeplitz	$C_{d,s} = \frac{2^{2s} s\,\Gamma(s + d/2)}{\pi^{d/2} \Gamma(1-s)}$ 0 via FFT, PCG
FEM (dense, H-matrix)	Dense/H-matrix	$C_{d,s} = \frac{2^{2s} s\,\Gamma(s + d/2)}{\pi^{d/2} \Gamma(1-s)}$ 1, multigrid, CG
Meshless GMQ RBF	Dense/low-rank	$C_{d,s} = \frac{2^{2s} s\,\Gamma(s + d/2)}{\pi^{d/2} \Gamma(1-s)}$ 2, analytic + quadrature
Monte Carlo (WOS)	N/A	$C_{d,s} = \frac{2^{2s} s\,\Gamma(s + d/2)}{\pi^{d/2} \Gamma(1-s)}$ 3, parallel
LDG (mixed flux)	Dense/structured	Sparse local + nonlocal; fast BEM

6. Error Estimates and Convergence

For $C_{d,s} = \frac{2^{2s} s\,\Gamma(s + d/2)}{\pi^{d/2} \Gamma(1-s)}$ 4: FEM error $C_{d,s} = \frac{2^{2s} s\,\Gamma(s + d/2)}{\pi^{d/2} \Gamma(1-s)}$ 5 (Ainsworth et al., 2017), optimal rates when $C_{d,s} = \frac{2^{2s} s\,\Gamma(s + d/2)}{\pi^{d/2} \Gamma(1-s)}$ 6 is sufficiently smooth.
For piecewise-linear FE and adaptive refinement: energy norm error $C_{d,s} = \frac{2^{2s} s\,\Gamma(s + d/2)}{\pi^{d/2} \Gamma(1-s)}$ 7, $C_{d,s} = \frac{2^{2s} s\,\Gamma(s + d/2)}{\pi^{d/2} \Gamma(1-s)}$ 8 error $C_{d,s} = \frac{2^{2s} s\,\Gamma(s + d/2)}{\pi^{d/2} \Gamma(1-s)}$ 9 in 2D (Lischke et al., 2018, Faustmann et al., 2022).
FD schemes: $\widehat{(-\Delta)^s u}(\xi) = |\xi|^{2s}\widehat{u}(\xi).$ 0 for central schemes when splitting parameter is chosen optimally and $\widehat{(-\Delta)^s u}(\xi) = |\xi|^{2s}\widehat{u}(\xi).$ 1 is regular, $\widehat{(-\Delta)^s u}(\xi) = |\xi|^{2s}\widehat{u}(\xi).$ 2 for splitting/interpolation approaches (Wu et al., 2021, Huang et al., 2013).
GMQ RBF meshless approximation: algebraic rate $\widehat{(-\Delta)^s u}(\xi) = |\xi|^{2s}\widehat{u}(\xi).$ 3 if $\widehat{(-\Delta)^s u}(\xi) = |\xi|^{2s}\widehat{u}(\xi).$ 4, spectral if $\widehat{(-\Delta)^s u}(\xi) = |\xi|^{2s}\widehat{u}(\xi).$ 5 is RBF-native (Hao et al., 2024).
LDG methods: $\widehat{(-\Delta)^s u}(\xi) = |\xi|^{2s}\widehat{u}(\xi).$ 6 for $\widehat{(-\Delta)^s u}(\xi) = |\xi|^{2s}\widehat{u}(\xi).$ 7th-order polynomial basis, with numerical stability proven (Nie et al., 2021), and optimal rates on graded meshes (Han et al., 13 Dec 2025).
Walk-on-spheres: $\widehat{(-\Delta)^s u}(\xi) = |\xi|^{2s}\widehat{u}(\xi).$ 8 convergence, expected number of steps per path $\widehat{(-\Delta)^s u}(\xi) = |\xi|^{2s}\widehat{u}(\xi).$ 9 (Sheng et al., 2022).

7. Applications and Interpretive Remarks

The integral fractional Laplacian finds use in modeling nonlocal phenomena:

Anomalous diffusion and superdiffusive processes via Lévy flights (Lischke et al., 2018, Sheng et al., 2022).
Nonlocal elliptic and parabolic PDEs: fractional Poisson, Allen–Cahn, Gray–Scott, and reaction–diffusion systems (Duo et al., 2018, Hao et al., 2024).
Nonlocal obstacle problems with weighted boundary behavior and regularity estimates (Borthagaray et al., 2018).
PDEs with variable-order fractional Laplacians, efficiently handled by quasi-linear solvers and FFT-based methods (Hao et al., 2024).
Fractional Laplacian-driven pattern formation and Turing instability in reaction–diffusion systems (Ainsworth et al., 2017).

Boundary conditions are intrinsically nonlocal: exterior data must be prescribed for the Riesz definition, affecting solution regularity and boundary behaviors distinct from the spectral definition (Lischke et al., 2018). Mesh grading and weighted regularity theory are essential for accurate resolution of singular layers, especially in numerical implementations.

The operator is recommended whenever the modeling requires infinite-range interaction, genuine jumps, or physically motivated Lévy processes; for subordinated Brownian motion or local boundary data, spectral definitions may be more appropriate. Horizon-based truncated approximations are useful when restricting nonlocality to a finite interaction band (Lischke et al., 2018).