Multi-D Riesz Space Fractional Diffusion

• Multi-dimensional Riesz space fractional diffusion equations generalize classical diffusion by replacing the Laplacian with the Riesz fractional operator, enabling the modeling of anomalous, heavy-tailed transport.
• High-order discretization methods, including fourth-order compact differences and finite element schemes, are developed to ensure accurate spatial and temporal convergence with mesh-independent spectral preconditioning.
• These equations are pivotal in capturing jump-diffusion and Lévy-flight processes, with applications in physics, finance, and biology, demonstrating robust stability and convergence properties.

A multi-dimensional Riesz space fractional diffusion equation generalizes the classical diffusion equation by replacing the Laplacian with the Riesz fractional Laplacian, introducing spatial nonlocality and heavy-tailed propagation that models anomalous transport. Such equations are pivotal in capturing superdiffusive processes in physics, finance, and biology where jump-diffusion or Lévy-flight phenomena are significant. The analysis and simulation of these equations require high-precision finite difference or finite element discretizations, fast iterative solvers, and rigorous spectral analysis, all scaled to multi-dimensional domains.

1. Definition and Problem Setting

A typical multi-dimensional Riesz space-fractional diffusion(-reaction) equation on a rectangular domain $\Omega=\prod_{i=1}^d (a_i,b_i)$, $d\geq2$, is formulated as

$$\frac{\partial u(\mathbf{x},t)}{\partial t} = \sum_{i=1}^d K_i\, \frac{\partial^{\alpha_i}u(\mathbf{x},t)}{\partial|x_i|^{\alpha_i}} + f(\mathbf{x},t,u),$$

where $1<\alpha_i<2$, $K_i>0$, and $f$ is a possibly nonlinear source (e.g., reaction term) subject to initial data $u(\mathbf{x},0)=u_0(\mathbf{x})$ and homogeneous Dirichlet conditions $u(\mathbf{x},t)=0$ for $\mathbf{x} \in \partial\Omega$. The Riesz derivative in coordinate $x_i$ is given by

$$\frac{\partial^{\alpha_i}u}{\partial|x_i|^{\alpha_i}} = -\frac{1}{2\cos(\pi\alpha_i/2)}\left({}_{a_i}D_{x_i}^{\alpha_i}u + {}_{x_i}D_{b_i}^{\alpha_i}u\right),$$

where the left/right Riemann–Liouville operators of order $\alpha_i$ are

$${}_{a_i}D_{x_i}^{\alpha_i}u(\mathbf{x}) = \frac{1}{\Gamma(2-\alpha_i)}\frac{\partial^2}{\partial x_i^2} \int_{a_i}^{x_i} (x_i - \xi)^{1-\alpha_i} u(\ldots,\xi,\ldots)\,d\xi.$$

This operator is symmetric and positive definite under Dirichlet boundary conditions. The case $f\equiv 0$ reduces to the pure Riesz diffusion equation; variable coefficients and nonlinearity are common in modeling real systems (Qu et al., 2024, Huang et al., 31 Jul 2025).

2. Analytical Properties and Fundamental Solutions

The Riesz fractional Laplacian $(-\Delta)^{\alpha/2}$ has a Fourier symbol $|\mathbf{k}|^{\alpha}$, leading to the Cauchy problem

$$\partial_t u(\mathbf{x},t) = -D(-\Delta)^{\alpha/2}u,$$

whose fundamental solution is expressed as

$$G(\mathbf{x},t) = (2\pi)^{-d} \int_{\mathbb{R}^d} e^{i\mathbf{k}\cdot\mathbf{x} - D t|\mathbf{k}|^{\alpha}}d^d\mathbf{k},$$

which can be recast using Fox $H$ functions or Wright functions (Parisis et al., 2018). This kernel exhibits heavy tails $G(\mathbf{x},t)\sim |\mathbf{x}|^{-(d+\alpha)}$ for large $|\mathbf{x}|$, in contrast to the Gaussian decay for $\alpha=2$. The solution family interpolates between Gaussian diffusion and non-spreading decay as $\alpha$ transitions from $2$ to $0^{+}$.

For space-time-fractional equations, the fundamental solution $G_{\alpha,\beta,d}(\mathbf{x},t)$ (with time-fractional Caputo order $\beta$) is given by

$$G_{\alpha,\beta,d}(\mathbf{x},t) = (2\pi)^{-d} \int_{\mathbb{R}^d} e^{i\mathbf{k}\cdot\mathbf{x}} E_\beta(-|\mathbf{k}|^\alpha t^\beta)d^d\mathbf{k},$$

where $E_\beta$ is the Mittag-Leffler function. Subordination principles connect solutions for different $\alpha$ and $\beta$ (Luchko, 2018).

3. High-Order Numerical Discretizations

Spatial discretization for multi-dimensional Riesz fractional operators has progressed from first-order Grünwald–Letnikov-type stencils to high-order compact and quasi-compact schemes. For instance, the fourth-order spatial discretization is realized through the formula

$$(-\Delta_x)^{\alpha/2} u(x_i) = h_x^{-\alpha} \sum_{k=-N_x}^{N_x} s_k^{(\alpha)} u(x_{i-k}) + O(h_x^4),$$

with coefficients $s_k^{(\alpha)}$ specified via the Fourier transform, ensuring symmetric and positive-definite difference matrices in each coordinate, applicable dimension-wise (Qu et al., 2024, Huang et al., 31 Jul 2025, Hu et al., 2019, She et al., 2024). Second-order schemes are also well studied using weighted and shifted Grünwald operators (Tian et al., 2012, Chen et al., 2013).

Time discretization is commonly Crank–Nicolson (second order), often combined with explicit or semi-implicit linearizations for nonlinear terms: $$\frac{u^{m+1}-u^m}{\Delta t} = \frac12[\mathcal{L}_h u^{m+1} + \mathcal{L}_h u^m] + f_{lin},$$ where $\mathcal{L}_h$ is the discrete fractional Laplacian (Qu et al., 2024, Huang et al., 31 Jul 2025).

In multi-D, the global system at each time step is a linear system with a symmetric positive definite “diagonal plus multi-level Toeplitz” structure, exploiting Kronecker tensor assembly (Huang et al., 31 Jul 2025, She et al., 2024). ADI and LOD splitting are often employed to reduce solution of multi-dimensional systems to sequences of unidirectional or block-tridiagonal solves (Deng et al., 2013, Valizadeh et al., 2018, Hu et al., 2019).

4. Stability, Convergence, and Error Estimates

Analysis across the high-order schemes establishes unconditional stability and convergence:

• For the Crank–Nicolson–fourth-order-centered difference (CN–4FCD) schemes,

$$\|U_{\rm exact}^m - U^m\|_2 \leq C(\Delta t^2 + h_x^4 + h_y^4)$$

for uniform meshes (discrete $L_2$ norm) (Qu et al., 2024, Huang et al., 31 Jul 2025, Hu et al., 2019, She et al., 2024).

• Second-order difference schemes with WSGD and standard CN achieve

$\|e^n\| = O(\Delta t^2 + h^2)$

global error (Tian et al., 2012, Chen et al., 2013, Jian et al., 2018).

• The convergence rates are independent of mesh ratio or step sizes, i.e., no CFL restriction.

For nonlinearities satisfying a global Lipschitz bound, stability and convergence are proved using discrete energy methods and Grönwall-type inequalities. Compact ADI-type schemes preserve unconditional stability under mild constraints on time step $\tau$ with respect to the nonlinearity (Hu et al., 2019).

5. Fast Iterative Solvers and Spectral Preconditioning

Discretization yields multi-level symmetric Toeplitz systems that are highly ill-conditioned for fine meshes. Optimal preconditioners—specifically, the sine-transform-based “$\tau$-matrix” preconditioners—have been rigorously analyzed and shown to cluster the spectrum of the preconditioned matrix in $[3/8,2]$ or $[1/2,3/2]$ independent of the mesh parameters (Qu et al., 2024, Huang et al., 31 Jul 2025, Huang et al., 2021, She et al., 2024, Huang et al., 2021). For example,

$$\lambda_{\min}(P^{-1}A) \geq \frac{3}{8},\qquad \lambda_{\max}(P^{-1}A) \leq 2.$$

Consequently, the PCG iteration count is invariant to grid refinement, typically requiring $\approx$8–12 iterations for two and three-dimensional problems, even with variable coefficients or nontrivial nonlinearities.

The typical structure:

Scheme	Spatial Order	Temporal Order	Preconditioner Type	PCG Iteration Count
CN–4FCD	$O(h^4)$	$O(\Delta t^2)$	Sine-transform ($\tau$)	8–12
Compact ADI	$O(h^4)$	$O(\Delta t^2)$	N/A (tri/banded ADI)	O(1) per substep
Second-order FD	$O(h^2)$	$O(\Delta t^2)$	Circulant/sine	O(1)

This scalability is confirmed across both theory and numerical tests (Qu et al., 2024, Huang et al., 31 Jul 2025, Huang et al., 2021).

6. Generalizations and High-Dimensional Extension

All fourth-order compact difference and quasi-compact schemes, as well as the block-Toeplitz preconditioning logic, extend naturally to three and higher spatial dimensions. The Kronecker structure is preserved, ensuring that fast sine/cosine transforms in $O(N^d\log N)$ complexity per iteration remain valid, and spectral clustering is maintained (Qu et al., 2024, Huang et al., 31 Jul 2025, She et al., 2024).

Finite element approaches on irregular and convex domains require assembly of fractional stiffness matrices on unstructured meshes, with accurate quadrature over fractional paths; these are essential for preserving accuracy in non-rectangular geometries (Yang et al., 2016, Partohaghighi et al., 2023). ETD-RDP-FEM schemes deliver second-order time accuracy and optimal spatial accuracy dictated by the finite element basis (Partohaghighi et al., 2023).

7. Applications, Limitations, and Representative Results

Multi-dimensional Riesz space fractional diffusion models underpin the simulation of superdiffusive phenomena, anomalous transport, and nonlocal reaction–diffusion systems. State-of-the-art solvers now allow simulation on grids $N^d$ with $N\sim 10^3$ per axis, leveraging FFT-enabled Toeplitz/arithmetic and mesh-independent PCG convergence (Qu et al., 2024, Huang et al., 31 Jul 2025, Huang et al., 2021, She et al., 2024).

Key validations include:

• $L_2$ and $L_\infty$ error reductions by factors of $16$ per halving of $h$ in spatial fourth-order schemes (Qu et al., 2024, Huang et al., 31 Jul 2025).
• PCG iteration counts for T-matrix/sine preconditioners remain constant as system sizes increase to $>10^6$ unknowns, with CPU time scaling as $O(N\log N)$ (Qu et al., 2024, Huang et al., 2021, Huang et al., 31 Jul 2025).
• For convex or irregular domains, penalty or finite element frameworks achieve similar accuracy and scalability (Yang et al., 2016, Huang et al., 2021, Partohaghighi et al., 2023).

Limitations:

• The proofs for spectrum clustering for inexact (e.g., penalized or generalized) preconditioners rely partly on numerical evidence in non-rectangular domains (Huang et al., 2021).
• Assembly of multi-dimensional fractional stiffnesses for unstructured meshes may become computationally demanding without meshless/H-matrix acceleration.

In conclusion, the synthesis of high-order difference or finite element discretization, mesh-independent spectral preconditioning, and advanced ADI/LOD factorization constitutes the mature numerical analysis toolkit for multi-dimensional Riesz space fractional diffusion equations (Qu et al., 2024, Huang et al., 31 Jul 2025, Huang et al., 2021, She et al., 2024, Yang et al., 2016, Partohaghighi et al., 2023).