Multi-Dimensional Riesz Diffusion Equations

• Multi-dimensional Riesz space fractional diffusion equations are generalizations that capture nonlocal, anisotropic, and asymmetric anomalous transport phenomena across high-dimensional domains.
• They replace the classical Laplacian with fractional Riesz or Riesz–Feller operators to model long-range jumps, heavy-tailed propagation, and complex stochastic processes.
• Advanced numerical schemes, including high-order finite differences and mesh-independent preconditioners, provide unconditional stability and scalable solutions for these complex models.

Multi-dimensional Riesz space fractional diffusion equations constitute a major generalization of classical diffusion models, enabling the representation of spatially nonlocal, possibly anisotropic, and/or asymmetric anomalous transport phenomena in high-dimensional domains. These equations replace the standard Laplacian with a (potentially vector-valued) Riesz or Riesz–Feller fractional differential operator, providing a natural framework for modeling long-range jumps, heavy-tailed spatial propagation, and related behaviors in physics, biology, finance, and materials science. Modern theory addresses both the analytic structure (existence, probabilistic representation, and fundamental solution properties) and a rapidly evolving landscape of numerical methods, favored by advances in discrete high-order schemes and mesh-independent preconditioning.

1. Mathematical Structure and Operator Definition

The multi-dimensional Riesz space fractional diffusion equation is generically written as: $$\partial_t u(\mathbf{x}, t) = \sum_{i=1}^d K_i\, \frac{\partial^{\alpha_i} u(\mathbf{x}, t)}{\partial |x_i|^{\alpha_i}} + f(\mathbf{x}, t)$$ where $1 < \alpha_i < 2$ and $K_i > 0$ are spatially varying or constant diffusion coefficients. The operator $$\frac{\partial^{\alpha_i}}{\partial |x_i|^{\alpha_i}}$$ denotes the Riesz fractional derivative, defined through its Fourier symbol: $$\mathcal{F}\left\{\, \frac{\partial^{\alpha_i}}{\partial |x_i|^{\alpha_i}} u(\mathbf{x}, t)\,; k_i\right\} = -|k_i|^{\alpha_i} \hat{u}(k_1, ..., k_d, t)$$ The Riesz–Feller derivative generalizes this by introducing an asymmetry parameter $\theta_i$: $$\mathcal{F}\left\{\, {}_x D^{\alpha_i}_{\theta_i} u(\mathbf{x}, t);\, k_i \right\} = -|k_i|^{\alpha_i} \exp[i\, \operatorname{sign}(k_i)\, \theta_i \pi/2]\, \hat{u}(k_1,...,k_d, t)$$ where $|\theta_i| \leq \min\{\alpha_i, 2-\alpha_i\}$ and $\hat{u}$ denotes the spatial Fourier transform (1109.4841, Saxena et al., 2014, Saxena et al., 2015). The operator is pseudo-differential in nature; in multi-dimensional settings, one often considers superpositions of such operators along different coordinate axes or in directions associated with anisotropy (Prodanov, 2018).

2. Analytic Solution Structures and Special Functions

Closed-form analytic representations of the fundamental solution (Green's function) are available via Laplace and Fourier transform techniques. Solutions often express the time evolution, for suitable initial data, as: $$u(\mathbf{x}, t) = \int_{\mathbb{R}^d} e^{-i\mathbf{k} \cdot \mathbf{x}}\, f^*(\mathbf{k})\, E_{\alpha, \beta}\left( -t^\alpha \, \eta \, \psi_{\boldsymbol{\alpha}, \boldsymbol{\theta}}(\mathbf{k}) \right) \, d\mathbf{k}$$ where $E_{\alpha,\beta}(\cdot)$ denotes a (generalized) Mittag–Leffler function, $f^*(\mathbf{k})$ encodes initial conditions, and $$\psi_{\boldsymbol{\alpha}, \boldsymbol{\theta}}(\mathbf{k})$$ generalizes the symbol for joint multi-directional Riesz–Feller derivatives (Saxena et al., 2014, Saxena et al., 2015). Compact expressions are also obtained in terms of the Fox $H$-function, with the scalar propagator (fundamental solution) written as: $u(\mathbf{x}, t) = H\mbox{-function expression}$ These solutions generalize previous symmetric Riesz diffusion models by encoding both memory (via fractional time derivatives, e.g., Caputo, Hilfer) and spatial asymmetry or anisotropy (via multiple Riesz–Feller derivatives) (1109.4841, Saxena et al., 2015). For the space-time-fractional diffusion-wave equation, the Mellin–Barnes integral and subordination principles connect the multi-dimensional fundamental solutions to convolutions of “classical” (local) solutions via kernels constructed from Wright and generalized Wright functions, enforcing completely monotone, normalized pdf properties (Luchko, 2018).

An essential concern is the convergence of these analytic representations. Convergence of double and higher-series (e.g., arising from the generalized Mittag–Leffler expansion or Srivastava–Daoust hypergeometric functions) is rigorously addressed using parameter constraints, ensuring the well-posedness of the H-function and series representations (1109.4841).

3. Probabilistic and Physical Interpretation

The Riesz space fractional diffusion operator models Lévy flights—a class of spatially nonlocal Markov processes (stable processes) with heavy-tailed jump statistics. When generalized to the Riesz–Feller form, asymmetric jump distributions and directed anomalous transport naturally emerge. The solution may be interpreted as a spatial probability density function evolving in time according to a generalized, possibly anisotropic, Lévy process (1109.4841, Beghin, 2016, Prodanov, 2018).

Extensions involving exponential or logarithmic pseudo-differential operators (as functional calculus of Riesz–Feller derivatives) correspond to subordinate time-changed stable processes—yielding additional stochastic interpretations involving Poisson processes with drift and gamma subordinators, relevant, for example, to physical and financial applications (Beghin, 2016). The subordination paradigm in space-time–fractional equations further bridges classical and fractional models: solutions of the fractional model are written as integrals (convolutions) over solutions of the classical equation with respect to completely monotone probability kernels, often represented by the Mittag–Leffler or Wright function (Luchko, 2018).

4. Numerical Schemes and High-Dimensional Solver Framework

Modern numerical discretizations for multi-dimensional Riesz space fractional diffusion equations exploit the Toeplitz (or block-Toeplitz) structure induced by the fractional finite difference (or quasi-compact) stencils. High-order schemes, such as fourth-order fractional centered difference (4FCD) or quasi-compact finite difference methods, are now developed for multi-dimensional domains and variable coefficients (Huang et al., 31 Jul 2025, Qu et al., 6 May 2024, She et al., 16 Apr 2024). Typical time-stepping is Crank–Nicolson, backward Euler, BDF2, or exponential time differencing, with methods (e.g., ADI, LOD splitting, trapezoidal, leapfrog) rigorously optimized for unconditional stability and error order (Chen et al., 2013, Deng et al., 2013, Valizadeh et al., 2018, Hu et al., 2019, Partohaghighi et al., 2023).

The discrete space operator is a sum of $d$-level (symmetric) Toeplitz matrices, yielding linear systems of the form: $D + T\,,$ where $D$ is diagonal (e.g., variable coefficient mass-matrix) and $T$ is multilevel Toeplitz, typically dense due to long-range nonlocal couplings (Huang et al., 31 Jul 2025, She et al., 16 Apr 2024).

Efficient solvers are a focus; preconditioned Krylov subspace methods (CG, GMRES) are employed, with mesh-independent convergence facilitated by preconditioners based on the discrete sine transform ($\tau$-matrices, T-matrix) (Huang et al., 2021, Huang et al., 2021, Qu et al., 6 May 2024, Huang et al., 31 Jul 2025). It is rigorously proven that when preconditioning with an appropriate sine transform-based matrix for high-order Toeplitz-like systems, the spectrum of the preconditioned matrix remains uniformly bounded in an interval (e.g., $(3/8, 2)$ or $(1/2, 3/2)$), yielding convergence rates for iterative solution that do not deteriorate with mesh refinement or problem dimension (Huang et al., 2021, Qu et al., 6 May 2024, Huang et al., 31 Jul 2025).

Key features can be tabulated as:

Numerical method	Spatial accuracy	Stability	Preconditioner	Spectral bound
CN-4FCD + sine-transform CG	O($h^4$)	Unconditional	Sine-transform ($\tau$-mat)	$(3/8, 2)$ (Huang et al., 31 Jul 2025)
Quasi-compact FD + GMRES	O($h^4$)	Unconditional	Sine-transform ($\tau$-mat)	mesh-independent (She et al., 16 Apr 2024)
ADI 4th order (compact)	O($h^{2\alpha}, h^{2\beta}$)	Unconditional	N/A	N/A
FE (unstructured mesh)	O($h^{2-\max{\alpha,\beta}}$)	Proven	N/A	N/A

Spectral analyses show that, with properly constructed preconditioners, eigenvalues of the preconditioned matrices do not cluster towards zero as the mesh is refined, granting scalability even for $d\geq3$.

5. Convergence, Stability, and Error Analysis

Theoretical convergence and stability analyses are rigorously developed for both analytic and numerical schemes:

• All high-order Crank–Nicolson and quasi-compact FD schemes are unconditionally stable and convergent. Error estimates $$\|u_{numerical} - u_{exact}\| = O(\Delta t^2 + h^4)$$ in the $L_2$-norm are proven for the CN-4FCD scheme (Huang et al., 31 Jul 2025).
• The fully discrete FE schemes with Caputo or Riesz derivatives achieve saturation error orders $\mathcal{O}(\tau^2 + h^2)$ (Yue et al., 2017).
• With explicit consideration of non-smooth initial data and stiff nonlinear source terms, exponential time differencing FEM schemes demonstrate L-stability and sharper error bounds for both linear and nonlinear, even nonsmooth, problems (Partohaghighi et al., 2023).
• Convergence of the analytic series (e.g., those involving the Srivastava–Daoust function) relies on precise restrictions for the real parameters associated with the time and space fractional orders (1109.4841).

6. Extensions: Anisotropy, Reaction Dynamics, and Irregular Domains

Recent works extend the framework in several significant ways:

• Anisotropy is incorporated by combining Riesz (-Feller) derivatives along arbitrary directions or by explicit use of directional vectors (see (Prodanov, 2018)), leading to Green’s functions parameterized by the direction of propagation.
• Reaction-diffusion systems are unified by combining time-fractional derivatives (Caputo, Riemann–Liouville, Hilfer) of distributed or generalized orders with Riesz–Feller spatial operators. The resulting solutions, in closed forms involving generalized Mittag–Leffler or $H$-functions, subsume classic diffusion, anomalous reaction–diffusion, and even telegraph and memory-telegraph models (by varying operator parameters) (1109.4841, Saxena et al., 2015).
• Robust numerical treatment of nonlinearities, nonsmooth data, and irregular geometries is realized through the use of FEM and penalization strategies (for convex/irregular domains) in both stationary and time-dependent settings. Efficient assembly algorithms are designed for stiffness matrices on unstructured meshes, enabling accurate simulation on arbitrary convex domains (Yang et al., 2016, Partohaghighi et al., 2023).

7. Impact, Practical Applications, and Future Directions

Multi-dimensional Riesz space fractional diffusion models and their high-order numerics underpin the modeling of non-Gaussian transport, anomalously diffusive processes, and nonlocal information propagation in complex systems. Successful applications range from hydrology and turbulent transport to image processing and quantitative finance. The development of mesh-independent preconditioners for high-order discretizations and their rigorous spectral analysis represents a foundational computational advance, affording scalability to large-scale, variable-coefficient, multi-dimensional problems (Huang et al., 31 Jul 2025, Qu et al., 6 May 2024, Huang et al., 2021).

Future research directions include increasing order-of-accuracy for non-symmetric operators, further improving solver robustness for extreme anisotropy/heterogeneity, developing stochastic solution representations for extended reaction–diffusion systems, and extending spectral preconditioning theory to even more general classes of structured operators (She et al., 16 Apr 2024, Partohaghighi et al., 2023). The core methodologies and analytic frameworks for multi-dimensional Riesz space fractional diffusion equations continue to evolve, providing essential mathematical infrastructure to the paper of nonlocal and fractional phenomena in science and engineering.