high-dimensional problems, sum-of-Gaussians approximation, Fokker–Planck equation, sparse grids, tensor neural networks, fast algorithms

1 Introduction↩︎

The Fokker–Planck equation (FPE) provides a deterministic description of the time evolution of probability density functions for stochastic systems [1], [2], with applications in statistical mechanics, stochastic processes, mathematical finance, information theory, and machine learning [3]–[7]. In the classical regime, the underlying dynamics are typically driven by Gaussian white noise, leading to the well-known linear growth of mean-squared displacement, i.e., \(\langle |x|^2 \rangle \propto t\) [8]. However, many anomalous-transport models require non-Gaussian jump statistics with algebraic tails and a self-similar length scale that differs from the Brownian scale. In the spatially fractional model considered here, the Fourier symbol \(|\boldsymbol{k}|^{2\alpha}\) with \(0<\alpha<1\) corresponds to a symmetric stable process of index \(2\alpha\): the characteristic length grows like \(t^{1/(2\alpha)}\), while moments of order \(q\ge 2\alpha\) are infinite. Such heavy-tailed anomalous diffusion is modeled by the fractional Fokker–Planck equation (FFPE), where the classical Laplacian is replaced by a fractional Laplacian operator \((-\Delta)^\alpha\) [9], defined for \(\boldsymbol{x}\in \mathbb{R}^{d}\) through the Cauchy principal value \[\label{eq::fractional95operator} (-\Delta)^{\alpha}p(\boldsymbol{x})=\frac{2^{2\alpha}\Gamma(\alpha+d/2)}{\pi^{d/2}|\Gamma(-\alpha)|}\text{P.V.}\int_{\mathbb{R}^{d}}\frac{p(\boldsymbol{x})-p(\boldsymbol{y})}{|\boldsymbol{x}-\boldsymbol{y}|^{d+2\alpha}}\mathrm{d}\boldsymbol{y},\tag{1}\] where \(\Gamma(z)\) denotes the Gamma function.

The fractional Laplacian \((-\Delta)^\alpha\) makes the equation nonlocal and captures long-range jumps associated with Lévy stable processes [10]–[13]. We consider the following initial value problem for the high-dimensional FFPE: \[\label{eq::FFPE95formula} \begin{cases} \frac{\partial }{\partial t}p(\boldsymbol{x},t) = -\boldsymbol{b} \cdot \nabla p(\boldsymbol{x},t) + D_{o}\Delta p(\boldsymbol{x},t) - D_{f}(-\Delta)^{\alpha}p(\boldsymbol{x},t), \\ p(\boldsymbol{x},0) = \delta(\boldsymbol{x}-\boldsymbol{x}_{0}), \quad \boldsymbol{x} \in \mathbb{R}^{d}, \end{cases}\tag{2}\] where \(\boldsymbol{b} \in \mathbb{R}^d\) is the drift vector, and \(D_o \ge 0\) and \(D_f > 0\) are the ordinary and fractional diffusion coefficients. The solution of 2 is the fundamental solution (Green’s function) of the FFPE; by linearity it determines, through convolution, the solution for a general initial datum. We therefore take the Dirac-delta case as the basic building block and use the same representation for separated initial data in high dimension.

Numerical treatment of Eq. 2 in high dimensions poses substantial mathematical and computational challenges. First, traditional grid-based methods, such as finite difference or finite element schemes, suffer from the curse of dimensionality: the degrees of freedom grow exponentially with the dimension \(d\) [14], [15]. Second, the nonlocality of \((-\Delta)^\alpha\) leads to dense discretized operators, which are costly for large-scale problems. Monte Carlo sampling and deep-learning-based solvers [16], [17] avoid full grids, but the Dirac-delta initial condition and slow convergence can still be limiting factors [15]. Recent methods based on functional hierarchical tensors [18] and fundamental-solution integrals [19] improve this situation. In particular, Ye et al. [19] reduce the free-space FFPE with Dirac-delta initial data to a one-dimensional radial integral that is evaluated to high precision in low to moderate dimensions. That approach relies on radial quadrature: as \(d\) increases, the Bessel order \((d-2)/2\) and the power \(r^{d/2}\) in the integrand both grow, making the integrand more oscillatory and increasing its dynamic range; the method is demonstrated up to \(d=29\). This motivates the separable representation developed here, which is aimed at higher dimensions and at separated initial data. The two approaches are complementary: the radial integral remains effective in low to moderate dimension, whereas the present method is designed for high-dimensional separated representations.

To address these challenges, we develop a fast algorithm based on a sum-of-Gaussians (SOG) approximation of the fundamental solution. The fractional operator enters the Fourier-space solution only through the stretched-exponential factor \(\exp(-D_ft|\boldsymbol{k}|^{2\alpha})\), the one piece that does not factorize across coordinates. Because this factor is completely monotone [20], it is the Laplace transform of a one-sided \(\alpha\)-stable density, and a trapezoidal discretization of that representation turns it into a sum of Gaussians – each an ordinary heat kernel that factorizes across dimensions. The FFPE thus reduces to a short separated sum of decoupled heat solutions. On a tensor-product grid, its one-dimensional factors are assembled in \(O(MdN)\) work and storage rather than forming \(O(N^d)\) dense grid values, where \(M\) is the number of Gaussians and \(N\) the number of points per dimension.

This construction is accurate, admits a rigorous error analysis, and extends naturally to low-rank summation of separated initial data. We give a rigorous a priori error analysis for the pure-fractional kernel that fixes the quadrature step and truncation for any prescribed tolerance, with convergence governed by a complex-plane bound on the stable density. For prescribed physical windows, the same scaled formulation gives a domain-adapted parameter choice for all \(\alpha\in(0,1)\) and includes the ordinary-diffusion case \(D_o>0\); when \(\alpha=1/2\), a closed-form identity for the trapezoidal error provides an additional analytic reference. In numerical experiments the method attains more than ten digits of relative accuracy with \(M\) growing only logarithmically in the inverse tolerance, and sustains this accuracy up to \(d=10^{5}\), well beyond the dimensions reachable by radial-quadrature methods. Because the relative error obeys a self-similar scaling, a single approximation sized at the smallest time serves an entire space-time window and avoids the small-time, high-dimensional loss of accuracy observed for direct quadrature [19]. Finally, since the method approximates the fundamental solution by a separated sum of Gaussians, any initial datum written as a sum of tensor products evolves in closed form through one-dimensional convolutions alone – a class of high-dimensional functions that is computable and amenable to error analysis [21]–[23]. Sparse-grid approximation methods [24], [25] and tensor neural networks [26]–[28] provide complementary ways to construct reduced representations for more general data.

Gaussian-sum approximations also arise in other high-dimensional PDE contexts. For example, in time-independent many-electron Schrödinger eigenvalue problems, mixed-derivative regularity [29] provides analytic support for sparse-grid approximations, and sparse-grid methods have been developed for the Schrödinger equation [30]. The pairwise Coulomb kernel \(1/|\boldsymbol{r}_i-\boldsymbol{r}_j|\), which is nonseparable in the electronic coordinates, admits accurate SOG approximations; after expansion into Gaussian factors, it is compatible with tensor-product integration and tensor neural network representations [31], [32]. For high-dimensional evolution problems, the same separability mechanism is relevant whenever a Fourier-space propagator, or a linear subproblem arising from time discretization, admits an accurate Gaussian-sum representation. For many nonlinear evolution equations, an unconditionally energy-stable scalar auxiliary variable (SAV) temporal discretization reduces each time step to a linear problem with known source terms [33]. If the resulting linear subproblem has a constant-coefficient solution operator that admits an accurate and efficient Gaussian-sum representation, the framework developed here extends naturally to such problems. These connections motivate SOG approximations as building blocks for separable representations in high-dimensional PDEs, although the analysis below is restricted to the FFPE fundamental solution.

The remainder of this paper is organized as follows. In 2, we review the mathematical preliminaries, including Fourier transforms, the theory of completely monotone functions, and properties of the stretched exponential function. In 3, we detail the SOG algorithm and provide its rigorous error estimate. Numerical experiments assessing the performance of the proposed solver are presented in 4, followed by concluding remarks in 5.

2 Preliminaries↩︎

2.1 Fourier transform↩︎

For a function \(f\in L^{1}(\mathbb{R}^{d})\cap L^{2}(\mathbb{R}^{d})\) we define its Fourier transform by \[\widehat{f}(\boldsymbol{k})=\mathcal{F}[f](\boldsymbol{k})=\int_{\mathbb{R}^{d}} f(\boldsymbol{x})\,e^{-i\,\boldsymbol{k}\cdot\boldsymbol{x}}\,\mathrm{d}\boldsymbol{x}, \qquad \boldsymbol{k}\in\mathbb{R}^{d},\] and the inverse Fourier transform by \[f(\boldsymbol{x})=\mathcal{F}^{-1}[\widehat{f}](\boldsymbol{x})=\frac{1}{(2\pi)^d}\int_{\mathbb{R}^{d}} \widehat{f}(\boldsymbol{k})\,e^{i\,\boldsymbol{k}\cdot\boldsymbol{x}}\,\mathrm{d}\boldsymbol{k}, \qquad \boldsymbol{x}\in\mathbb{R}^{d}.\] If \(f\) is radial, i.e.\(f(\boldsymbol{x})=f(r)\) with \(r=\sqrt{x_{1}^{2}+\dots+x_{d}^{2}}\), its Fourier transform is again radial. In this case, the radial Fourier transform (i.e., the Hankel transform) pairs are \[\begin{align} \widehat{f}(k)&=\frac{(2\pi)^{d/2}}{k^{(d-2)/2}} \int_{0}^{\infty} f(r)\,r^{\frac{d}{2}} J_{\frac{d-2}{2}}(kr)\,\mathrm{d}r,\\ f(r) &=\frac{1}{(2\pi)^{d/2} r^{(d-2)/2}} \int_{0}^{\infty} \widehat{f}(k)\,k^{\frac{d}{2}} J_{\frac{d-2}{2}}(kr)\,\mathrm{d}k, \end{align}\] where \(J_{\nu}\) is the Bessel function of the first kind of order \(\nu\). Finally, for a sufficiently regular and rapidly decaying \(f\) the Poisson summation formula [34] links sums over the integer lattice to sums over its dual: \[h\sum_{n=-\infty}^{\infty} f(hn + a) = \sum_{m=-\infty}^{\infty} e^{i \frac{2\pi m a}{h}} \widehat{f}\left(\frac{2\pi m}{h}\right), \label{eq:poissonsf}\tag{3}\] providing a powerful bridge between spatial and frequency-domain information.

2.2 Completely monotone functions↩︎

The following result, Bernstein’s theorem, provides a crucial integral representation that is often used as an alternative definition (see, for example, [35]).

2.3 Properties of the stretched exponential and its inverse Laplace transform↩︎

The function \[f(x)=e^{-x^\alpha}, \qquad 0<\alpha<1, \label{eq:kwwfun}\tag{4}\] also known as the stretched exponential or Kohlrausch–Williams–Watts (KWW) function, possesses several important properties.

The following properties of \(\rho_{\alpha}\) can be found in [36].

3 A fast sum-of-Gaussians FFPE solver↩︎

In this section, we present a fast algorithm for solving the FFPE 2 . By combining the radial Fourier-integral representation of the FFPE solution with a sum-of-Gaussians (SOG) approximation of the stretched exponential, the proposed method reduces the anomalous-diffusion solution to a sum of heat-equation solutions with closed-form expressions, thereby achieving a computational cost that grows linearly with the dimension.

3.1 SOE approximation of the stretched exponential function↩︎

We approximate the stretched exponential by a sum of exponentials (SOE), \[e^{-x^{\alpha}} \approx \sum_{\ell=-M_1}^{M_2} w_\ell e^{-s_\ell x}, \qquad x\in [\delta,R]. \label{eq:soeappr}\tag{5}\] Our starting point is the integral representation ?? . Applying the change of variables \(t=e^{u}\) to ?? , we obtain \[e^{-x^{\alpha}} = \int_{-\infty}^{\infty} e^{-x e^u+u}\rho_\alpha(e^u)\mathrm{d}u. \label{eq:kwwfunrep2}\tag{6}\] The integrand decays rapidly to zero as \(u\to\pm\infty\), so the trapezoidal rule converges exponentially fast [37], and the nodes \(s_\ell\) and weights \(w_\ell\) in 5 are given by \[s_\ell = e^{h\ell}, \qquad w_\ell = h s_\ell \rho_\alpha(s_\ell),\] where \(h>0\) is the step size in the trapezoidal rule.

We analyze the approximation error of 5 . We record the asymptotic magnitude of the Gamma function for a complex argument (Lemma 4, used in the parameter selection of 3.4); its proof is given in the appendix of [38].

We first analyze the discretization error of the (infinite) trapezoidal rule \[\label{eq::SOE} \begin{align} e^{-x^{\alpha}}&\approx \sum_{\ell=-\infty}^{+\infty}w_\ell e^{-s_\ell x}. \end{align}\tag{7}\]

Theorem 1 bounds the discretization error of the infinite rule uniformly on \([0,+\infty)\). Truncating the series to \(\ell\in[-M_1,M_2]\), as in 5 , restricts the accurate range to the finite window \([\delta,R]\), whose endpoints are governed by the extreme retained nodes: since a single term \(e^{-s_\ell x}\) acts on the scale \(x\sim s_\ell^{-1}\), the largest node \(s_{M_2}=e^{hM_2}\) sets the lower limit \(\delta\sim s_{M_2}^{-1}\) and the smallest node \(s_{-M_1}=e^{-hM_1}\) the upper limit \(R\sim s_{-M_1}^{-1}\). The truncation indices \(M_1,M_2\) are chosen in 3.3 to meet the target tolerance.

Theorem 2 provides a uniform bound of \(I_\alpha\).

The detailed proof is provided in Appendices 6 and 7. The estimate in ?? shows only mild endpoint singularities. The large factor that appears in a direct pointwise bound for \(|\rho_\alpha(z)|\) near \(\alpha=1\) is not intrinsic to \(I_\alpha\); in the proof of 7 the compact part is integrated in the radial variable before the saddle-contour integral is estimated, which preserves the cancellation in the phase. At \(\alpha=1\) the representing measure degenerates to a Dirac mass at \(z=1\), while the limit \(\alpha\to0^+\) is also degenerate and is no longer described by a regular probability density on \((0,\infty)\). Accordingly, endpoint regimes still require care in numerical density evaluation, but the contour constant entering the trapezoidal error does not grow like \(c^{1/(1-\alpha)}\).

3.2 The SOG approximation of the solution↩︎

Taking the Fourier transform of 2 with respect to the space variable \(\boldsymbol{x}\in \mathbb{R}^d\), we obtain \[\label{eq::FFPE95Fourier} \left\{\begin{align} &\frac{\partial}{\partial t}\widehat{p}(\boldsymbol{k}, t) =-i(\boldsymbol{k}\cdot\boldsymbol{b}) \widehat{p}(\boldsymbol{k}, t)-D_{\mathrm{o}}|\boldsymbol{k}|^2 \widehat{p}(\boldsymbol{k}, t)-D_{\mathrm{f}}|\boldsymbol{k}|^{2 \alpha} \widehat{p}(\boldsymbol{k}, t) \\ &\widehat{p}(\boldsymbol{k}, 0) =\exp \left(-i\boldsymbol{k}\cdot\boldsymbol{x}_0\right), \end{align}\right.\tag{13}\] which can be explicitly solved with \[\label{eq::p95hat95explicit} \widehat{p}(\boldsymbol{k}, t)=e^{-i\boldsymbol{k}\cdot\boldsymbol{x}_0^t}e^{-(D_o|\boldsymbol{k}|^2+D_f|\boldsymbol{k}|^{2\alpha})t},\quad \boldsymbol{x}_0^t:=\boldsymbol{x}_0+\boldsymbol{b}t.\tag{14}\] By the inverse Fourier transform, we have \[\label{eq::p95intergral} p(\boldsymbol{x},t)=\frac{1}{(2\pi)^d}\int_{\mathbb{R}^d}e^{i\boldsymbol{k}\cdot(\boldsymbol{x}-\boldsymbol{x}_0^t)}e^{-(D_o|\boldsymbol{k}|^2+D_f|\boldsymbol{k}|^{2\alpha})t}\mathrm{d}\boldsymbol{k}.\tag{15}\] Since the factor \(\exp(-(D_o|\boldsymbol{k}|^2+D_f|\boldsymbol{k}|^{2\alpha})t)\) is radially symmetric and hence invariant under coordinate rotations, we choose a new orthogonal coordinate frame \((\boldsymbol{e}_1,\cdots,\boldsymbol{e}_d)\) such that \(\boldsymbol{x}-\boldsymbol{x}_0^{t}\) lies along the first coordinate axis, that is, \(\boldsymbol{x}-\boldsymbol{x}_0^{t}=y\boldsymbol{e}_1\) with the direction vector \(|\boldsymbol{e}_1|=1\). For \(d\ge2\) and \(y>0\), by performing a \(d\)-dimensional spherical coordinate transformation \((k_1,k_2,\cdots,k_d)\rightarrow (r,\theta_1,\theta_2,\cdots,\theta_{d-1})\), the integral in 15 is equivalent to \[\label{eq::p95int95equiv} p(\boldsymbol{x},t)=\frac{\int_0^\pi \sin ^{d-2}(\theta) \int_0^{\infty} r^{d-1} \cos (\cos (\theta) y r) e^{-(D_or^2+D_f r^{2\alpha})t} \mathrm{d}r \mathrm{d} \theta}{2^{d-1}\pi^{\frac{d+1}{2}}\Gamma(\frac{d-1}{2})}\tag{16}\] where \(\Gamma(\cdot)\) denotes the Gamma function, and \(\theta\) abbreviates the first angular coordinate \(\theta_1\). Performing the angular integration for \(y\neq0\), and evaluating 15 directly for \(y=0\), gives the explicit radial integral expression [19] \[\label{eq::int95expression} p(\boldsymbol{x},t)= \left\{ \begin{align} &\frac{1}{y^{\frac{d-2}{2}}}\int_0^{\infty}\left(\frac{r}{2\pi}\right)^{\frac{d}{2}}J_{\frac{d-2}{2}}(yr)\exp(-(D_or^2+D_fr^{2\alpha})t)\mathrm{d}r,\;y\neq 0\\ &\frac{2^{1-d}}{\pi^{d/2}\Gamma\left(d/2\right)}\int_0^{\infty}r^{d-1}\exp(-(D_or^2+D_fr^{2\alpha})t)\mathrm{d}r,\;y=0, \end{align} \right.\tag{17}\] where \(J_\nu(\cdot)\) represents the \(\nu\)-th Bessel function of the first kind \[\label{eq::Bessel} J_\nu(z)=\frac{(z/2)^{\nu}}{\pi^{1/2}\Gamma(\nu+1/2)}\int_0^{\pi}\sin^{2\nu}(\theta)\cos(\cos(\theta)z)\mathrm{d}\theta.\tag{18}\] For a fixed terminal time \(T > 0\) and \(t\in[0,T]\), we apply the SOE expansion 5 to the fractional factor in \(\widehat{p}(\boldsymbol{k},t)\), namely \(e^{-D_ft|\boldsymbol{k}|^{2\alpha}}=e^{-x^{\alpha}}\) with \(x=(D_ft)^{1/\alpha}|\boldsymbol{k}|^2\), giving the sum-of-Gaussians (SOG) approximation \[\label{eq::SOG95approx} \widehat{p}(\boldsymbol{k}, t)\approx\widehat{p}_{\text{SOG}}(\boldsymbol{k},t)=\sum_{\ell=-M_1}^{M_2}w_\ell e^{-i\boldsymbol{k}\cdot\boldsymbol{x}_0^t}e^{-[D_ot+s_\ell(D_ft)^{1/\alpha}]|\boldsymbol{k}|^2}\tag{19}\] In real space, this procedure is equivalent to approximating the true solution \(p(\boldsymbol{x}, T)\) through a linear superposition of ordinary diffusion kernels. Specifically, this approximation is formulated as \[\label{eq::p95approx} p(\boldsymbol{x},T)\approx p_{\text{SOG}}(\boldsymbol{x},T):=\sum_{\ell=-M_1}^{M_2}w_\ell p_\ell(\boldsymbol{x},T),\tag{20}\] where \(p_\ell(\boldsymbol{x},T)\) is the solution to the following heat equation \[\label{eq::heat} \left\{ \begin{align} \frac{\partial}{\partial t}p_\ell(\boldsymbol{x},t)&=-\boldsymbol{b}\cdot\nabla p_\ell(\boldsymbol{x},t)+\left(D_{o}+\frac{s_\ell D_f}{\alpha} (D_ft)^{\frac{1}{\alpha}-1}\right)\Delta p_\ell(\boldsymbol{x},t)\\ p_\ell(\boldsymbol{x},0)&=\delta_{\boldsymbol{x}_0}(\boldsymbol{x}). \end{align} \right.\tag{21}\] The exact solution of 21 reads \[\label{eq::Heat95explicit} p_\ell(\boldsymbol{x},T)=\frac{1}{(4\pi C_\ell^{T})^{d/2}}\exp\left(-\frac{y^2}{4C_\ell^{T}}\right),\tag{22}\] where \(C_\ell^{T}:=D_oT+s_\ell (D_fT)^{1/\alpha}\) is the \(\ell\)-th ordinary diffusion constant in the SOG approximation.

Therefore, on a tensor-product observation grid \(\boldsymbol{X}=\otimes_{j=1}^d X_j\), each kernel \(p_\ell\) in 22 is represented by its one-dimensional Gaussian factors. Assembling the factors for one term costs \(O(\sum_{j=1}^d n_j)\) work and storage, and assembling all factors costs \(O(M\sum_{j=1}^d n_j)\), or \(O(MdN)\) when \(n_j=N\) for all \(j\). This is the cost of the separated representation; forming and storing all dense values on the full tensor grid would still require \(O(N^d)\) entries.

3.3 Error estimate of the SOG approximation↩︎

We now analyze the SOG approximation error of the FFPE solution given by 20 . Denote \[\mathcal{E}_T=\|p(\boldsymbol{x},T)-p_{\text{SOG}}(\boldsymbol{x},T)\|_{\infty}\] as the \(L_\infty\) error at the fixed terminal time \(T\). By the inverse Fourier transform, \(\mathcal{E}_T\) admits the simple bound \[\label{eq::IFT95error} \mathcal{E}_T\le \frac{1}{(2\pi)^d}\int_{\mathbb{R}^d}|\widehat{e}(\boldsymbol{k},T)|\mathrm{d}\boldsymbol{k},\quad \widehat{e}(\boldsymbol{k},T):=\widehat{p}(\boldsymbol{k},T)-\widehat{p}_{\text{SOG}}(\boldsymbol{k},T).\tag{23}\] Using the explicit expression in 14 and the invariance of Gaussians under the Fourier transform, we obtain \[\label{eq::IFT95error95explicit} \mathcal{E}_T\le \frac{1}{(2\pi)^d}\int_{\mathbb{R}^d}\left|e^{-D_oT|\boldsymbol{k}|^2}\left(e^{-D_fT|\boldsymbol{k}|^{2\alpha}}-\sum_{\ell=-M_1}^{M_2}w_\ell e^{-s_\ell(D_fT)^{1/\alpha}|\boldsymbol{k}|^2}\right)\right|\mathrm{d}\boldsymbol{k}.\tag{24}\] In fact, the error \(\mathcal{E}_T\) stems from the trapezoidal discretization of the kernel integral and the truncation of the infinite series at both ends; hence, we decompose it into three parts \[\label{eq::error95split} \mathcal{E}_T\le \mathcal{E}_\infty+\mathcal{E}_{\text{up}}+\mathcal{E}_{\text{down}},\tag{25}\] where \[\begin{align} \mathcal{E}_\infty&=\frac{1}{(2\pi)^d}\int_{\mathbb{R}^d}\left|e^{-D_oT|\boldsymbol{k}|^2}\left(e^{-D_fT|\boldsymbol{k}|^{2\alpha}}-\sum_{\ell=-\infty}^{\infty}w_\ell e^{-s_\ell(D_fT)^{1/\alpha}|\boldsymbol{k}|^2}\right)\right|\mathrm{d}\boldsymbol{k},\\ \mathcal{E}_{\text{up}}&=\frac{1}{(2\pi)^d}\int_{\mathbb{R}^d}\left|e^{-D_oT|\boldsymbol{k}|^2}\sum_{\ell=M_2+1}^{\infty}w_\ell e^{-s_\ell(D_fT)^{1/\alpha}|\boldsymbol{k}|^2}\right|\mathrm{d}\boldsymbol{k},\\ \mathcal{E}_{\text{down}}&=\frac{1}{(2\pi)^d}\int_{\mathbb{R}^d}\left|e^{-D_oT|\boldsymbol{k}|^2}\sum_{\ell=-\infty}^{-M_1-1}w_\ell e^{-s_\ell(D_fT)^{1/\alpha}|\boldsymbol{k}|^2}\right|\mathrm{d}\boldsymbol{k}.\\ \end{align}\] Clearly, when the ordinary diffusion is absent, i.e., \(D_o=0\), controlling the error is most challenging because there is no additional Gaussian damping. We therefore prove the conservative a priori estimates in this worst case. When \(D_o>0\), the same scaled representation is used, with the extra parameter \(\lambda(t)\) retained over the prescribed time window in 3.4.

For \(\mathcal{E}_\infty\), since the integrand is radially symmetric, we switch to \(d\)-dimensional spherical coordinates and set \(u=(D_fT)^{1/\alpha}|\boldsymbol{k}|^2\). Then the integral in 24 can be reduced to a one-dimensional integral with respect to \(u\), namely, \[\label{IFT951d} \mathcal{E}_{\infty}= \frac{\Omega_d}{2(2\pi)^d}(D_fT)^{-\frac{d}{2\alpha}}\int_{0}^{\infty}|\epsilon_\infty(u)|u^{\frac{d}{2}-1}\mathrm{d}u,\tag{26}\] where \(\Omega_d=2\pi^{d/2}/\Gamma(d/2)\) denotes the surface area of the unit sphere in \(\mathbb{R}^d\), and \[\epsilon_\infty(u):=e^{-u^{\alpha}}-\sum_{\ell=-\infty}^{+\infty}w_\ell e^{-s_\ell u}\] denotes the approximation error of the kernel function for \(u\in [0,+\infty)\), which is bounded pointwise by \(|\epsilon_\infty(u)|\le \epsilon_{\text{SOE}}\) under the parameter choice of Theorem 1. This uniform bound is useful only on a finite interval; it cannot by itself be integrated over \([U,\infty)\) against the growing measure factor \(u^{d/2-1}\). We therefore split the integral in 26 over \([0,U]\) and \([U,+\infty)\), and bound the tail directly. Let \[\label{eq::tail95infty} \begin{align} \mathcal{T}_h(U)={}&\int_U^\infty e^{-u^\alpha}u^{\frac{d}{2}-1}\,\mathrm du\\ &+\sum_{\ell=-\infty}^{\infty}w_\ell \int_U^\infty e^{-s_\ell u}u^{\frac{d}{2}-1}\,\mathrm du \\ ={}&\frac{1}{\alpha}\Gamma\!\left(\frac{d}{2\alpha},U^\alpha\right) +\sum_{\ell=-\infty}^{\infty}w_\ell s_\ell^{-\frac{d}{2}} \Gamma\!\left(\frac{d}{2},s_\ell U\right), \end{align}\tag{27}\] where \(\Gamma(a,z)\) is the upper incomplete Gamma function. Then \[\label{eq::E95infty95estimate} \mathcal{E}_{\infty}\le \frac{1}{2^d\pi^{d/2}\Gamma(d/2)}(D_fT)^{-\frac{d}{2\alpha}}\left[\frac{2\epsilon_{\text{SOE}}}{d}U^{\frac{d}{2}}+\mathcal{T}_h(U)\right].\tag{28}\]

Next, we analyze the error bounds for \(\mathcal{E}_{\text{up}}\) and \(\mathcal{E}_{\text{down}}\). By applying the same change of variables as for \(\mathcal{E}_\infty\), the exponential integrals can be evaluated in closed form. Indeed, let \[\begin{align} I_\ell&:=w_\ell \int_{0}^{\infty}e^{-s_\ell u}u^{\frac{d}{2}-1}\mathrm{d}u=h\Gamma\left(d/2\right)s_\ell^{1-\frac{d}{2}}\rho_\alpha(s_\ell), \end{align}\] and \(\mathcal{E}_{\text{up}}\) is then estimated using property 2 of Lemma 3, which gives \[\label{eq::err95up} \begin{align} \mathcal{E}_{\text{up}}&=\frac{1}{2^d\pi^{d/2}\Gamma(d/2)}(D_fT)^{-\frac{d}{2\alpha}}\sum_{\ell=M_2+1}^{\infty}I_\ell\\ &= \frac{1}{2^d\pi^{d/2}\Gamma(d/2)}(D_fT)^{-\frac{d}{2\alpha}}\sum_{\ell=M_2+1}^{\infty}h\Gamma\left(d/2\right)s_\ell^{1-\frac{d}{2}}\rho_\alpha(s_\ell)\\ &\le \sum_{\ell=M_2+1}^{\infty}\frac{h}{2^d\pi^{d/2}}(D_fT)^{-\frac{d}{2\alpha}}s_\ell^{1-\frac{d}{2}}\cdot C_\alpha s_\ell^{-(1+\alpha)}\\ &=\frac{hC_\alpha (D_fT)^{-\frac{d}{2\alpha}}}{2^d\pi^{d/2}(e^{h(\alpha+d/2)}-1)}e^{-h(\alpha+d/2)M_2}. \end{align}\tag{29}\] For \(\mathcal{E}_{\text{down}}\) we invoke the asymptotic estimate of \(\rho_\alpha(t)\) as \(t\to 0^{+}\) (property 3 of Lemma 3): \[\label{eq::err95down} \begin{align} \mathcal{E}_{\text{down}}&=\frac{1}{2^d\pi^{d/2}\Gamma(d/2)}(D_fT)^{-\frac{d}{2\alpha}}\sum_{\ell=-\infty}^{-M_1-1}I_\ell\\ &=\frac{1}{2^d\pi^{d/2}\Gamma(d/2)}(D_fT)^{-\frac{d}{2\alpha}}\sum_{\ell=-\infty}^{-M_1-1}h\Gamma\left(d/2\right)s_\ell^{1-\frac{d}{2}}\rho_\alpha(s_\ell)\\ &\le \sum_{\ell=M_1+1}^{\infty}\frac{h (D_fT)^{-\frac{d}{2\alpha}}}{2^d\pi^{d/2}}\cdot A_\alpha e^{h\ell(\gamma+\frac{d}{2}-1)}\exp\big(-D e^{h\ell\frac{\alpha}{1-\alpha}}\big)\\ &\le \frac{C_{0}(D_fT)^{-\frac{d}{2\alpha}}}{2^d\pi^{d/2}}\cdot hA_\alpha e^{hM_1(\gamma+\frac{d}{2}-1)}\exp\big(-D e^{hM_1\frac{\alpha}{1-\alpha}}\big),\\ \end{align}\tag{30}\] Here \(C_0\) denotes a tail-majorization constant after \(M_1\) is chosen beyond the maximizer of \(e^{h\ell(\gamma+d/2-1)}\exp\{-D e^{h\ell\alpha/(1-\alpha)}\}\); from that point the summand is decreasing, and the remaining lower tail is bounded by a fixed multiple of the displayed value at \(M_1\).

Next, based on the error bounds in 28 29 30 , we provide a rigorous parameter-selection strategy. Given a prescribed target tolerance \(\epsilon\), split it into budgets \(\epsilon_\infty+\epsilon_{\mathrm{up}}+\epsilon_{\mathrm{down}}\le\epsilon\). For the infinite-rule part, choose \(U\) and \(\epsilon_{\text{SOE}}\) so that \[\label{eq::para95U} \mathcal{T}_h(U)\le 2^{d-1}\pi^{\frac{d}{2}}\Gamma(d/2)(D_fT)^{\frac{d}{2\alpha}}\epsilon_\infty .\tag{31}\] The SOE tolerance \(\epsilon_{\text{SOE}}\) is chosen to satisfy \[\label{eq::para95SOE} \epsilon_{\text{SOE}}\le \frac{d\epsilon_\infty}{4}\cdot(D_fT)^{\frac{d}{2\alpha}}\cdot 2^d\pi^{\frac{d}{2}}\Gamma(d/2)\cdot U^{-\frac{d}{2}},\tag{32}\] which in turn fixes the trapezoidal step size \(h\) via ?? . The tail condition 31 is then checked with this \(h\); if necessary, \(U\) and \(\epsilon_{\text{SOE}}\) are adjusted iteratively. Finally, we determine the SOG truncation parameters \(M_1\) and \(M_2\). For the upper truncation \(M_2\), the selection becomes \[\label{eq::para95M2} M_2\ge \frac{1}{h(\alpha+\frac{d}{2})}\log\left[\frac{hC_\alpha (D_fT)^{-\frac{d}{2\alpha}}}{2^d\pi^{\frac{d}{2}}\epsilon_{\mathrm{up}}(e^{h(\alpha+\frac{d}{2})}-1)}\right].\tag{33}\] For the lower truncation we take \(M_1=\log X/h\), where \(X\) satisfies \[\label{eq::para95M1} X^{\gamma+\frac{d}{2}-1}\exp\left(-DX^{\frac{\alpha}{1-\alpha}}\right)\le \frac{2^d\pi^{\frac{d}{2}}\epsilon_{\mathrm{down}}\cdot (D_fT)^{\frac{d}{2\alpha}}}{C_0 A_\alpha\cdot h},\tag{34}\] which provides the asymptotic bound \[\label{eq::para95M195asymp} M_1\ge \frac{1}{h}\log\left[X_0+\frac{(1-\alpha)(\gamma+\frac{d}{2}-1)}{D\alpha}X_0^{\frac{1-2\alpha}{1-\alpha}}\log X_0\right]\tag{35}\] with \[\label{eq::X950} X_0=\left[\frac{1}{D}\log\left(\frac{C_0A_\alpha h }{2^d\pi^{\frac{d}{2}}\epsilon_{\mathrm{down}}\cdot (D_fT)^{\frac{d}{2\alpha}}}\right)\right]^{\frac{1-\alpha}{\alpha}}.\tag{36}\]

3.4 Parameter selection for a prescribed tolerance and evaluation window↩︎

The estimates in [sec::SOE_approx] [sec::Error] give an all-\(\alpha\) a priori construction for the scalar multiplier \(e^{-x^\alpha}\) and its induced Green’s-function approximation. This subsection formulates the corresponding parameter choice directly at the level of the scaled Green’s function. The construction applies for every \(\alpha\in(0,1)\); the special value \(\alpha=1/2\) enters only through the closed-form identities recorded at the end of the subsection.

Scaled Green’s-function representation. Let \(y=|\boldsymbol{x}-\boldsymbol{x}_0-\boldsymbol{b}t|\) and define \[\label{eq::scaled95profile} \eta=\frac{y}{(D_ft)^{1/(2\alpha)}},\qquad \lambda(t)=\frac{D_ot}{(D_ft)^{1/\alpha}}.\tag{37}\] The following proposition separates the self-similar scaling from the finite Gaussian quadrature and transfers relative-error estimates from the scaled profile to the physical Green’s function.

The relevant range of the kernel argument. The SOE approximates \(e^{-x^\alpha}\), and in the solution its Fourier-side argument is \(x=(D_ft)^{1/\alpha}|\boldsymbol{k}|^2\). Writing the solution radially and changing variables from \(|\boldsymbol{k}|\) to \(x\) (the substitution used for \(\mathcal{E}_\infty\) in 3.3), \[\label{eq::radial95x} p(Y,t)\;\propto\;(D_ft)^{-\frac{d}{2\alpha}}\int_0^\infty x^{\,d/2-1}\, \Lambda_d\!\big(\eta\sqrt{x}\big)\,e^{-x^\alpha}\,\mathrm{d}x,\qquad \eta=\frac{|Y|}{(D_ft)^{1/(2\alpha)}},\tag{41}\] where \[\Lambda_d(z)=2^{(d-2)/2}\Gamma\!\left(\frac{d}{2}\right) z^{-(d-2)/2}J_{(d-2)/2}(z)\] is the normalized radial kernel (\(\Lambda_d(0)=1\)). The solution therefore samples \(e^{-x^\alpha}\) only through the weight \(x^{d/2-1}\Lambda_d(\eta\sqrt{x})\). For the pure-fractional case \(D_o=0\), a physical window \(|Y|\le R\), \(t\in[t_{\min},t_{\max}]\) enters this scaled description through the largest self-similar displacement \[\eta_{\max}=\frac{R}{(D_ft_{\min})^{1/(2\alpha)}}.\] The active scalar interval \(x\in[x_{\min},x_{\max}]\) is chosen so that this weight is retained on the region where its envelope exceeds an \(\epsilon\)-level threshold. The two endpoints are governed by different data:

• the upper end \(x_{\max}\) is independent of \(Y\) and \(t\): for \(d>2\) the radial envelope \(x^{d/2-1}e^{-x^\alpha}\) peaks at \(x_\star=\big(\tfrac{d-2}{2\alpha}\big)^{1/\alpha}\), and \(x_{\max}\) is the largest \(x\) with \((\tfrac d2-1)\log\tfrac{x}{x_\star}-\big(x^\alpha-x_\star^\alpha\big)=\log\epsilon\). For \(d\le2\) the same upper-tail scale is obtained directly from \(e^{-x^\alpha}\lesssim\epsilon\). Thus \(x_{\max}\sim\max\!\big(x_\star,\,(\log\tfrac1\epsilon)^{1/\alpha}\big)\) for \(d>2\), and \(x_{\max}\sim(\log\tfrac1\epsilon)^{1/\alpha}\) for \(d\le2\);

• the lower end \(x_{\min}\) is the far-field, low-frequency cutoff set by \(\eta_{\max}\): as \(\Lambda_d(\eta_{\max}\sqrt{x})\) departs from unity only once \(\eta_{\max}\sqrt{x}\gtrsim\sqrt{d}\), the window extends down to \(x_{\min}\sim x_\star/(1+\eta_{\max}^2)\).

Thus \(x_{\max}\) (from \(\alpha,d,\epsilon\)) fixes the smallest node and hence \(M_1\), \(x_{\min}\) (from \(\eta_{\max}\)) the largest node and hence \(M_2\), and the step \(h\) controls the discretization error. When \(D_o>0\), the Fourier-side factor \(e^{-D_ot|\boldsymbol{k}|^2}\) gives additional Gaussian damping; the comparison domain in 1 nevertheless retains the full dependence on \(\lambda(t)\).

A domain-adapted selection procedure. When \(D_o=0\), \(\lambda=0\) and a fixed physical window \(0\le y\le R\), \(t\in[t_{\min},t_{\max}]\) reduces to \[\label{eq::etamax95selector} 0\le\eta\le\eta_{\max},\qquad \eta_{\max}=\frac{R}{(D_ft_{\min})^{1/(2\alpha)}},\tag{42}\] which formalizes Remark [rmk::selfsim]: the selection depends on \(t_{\min}\), through \(\eta_{\max}\), but not on \(t_{\max}\). When \(D_o>0\), \(\lambda(t)\) must be retained over the whole interval in the comparison domain. Except in the closed-form case \(\alpha=1/2\), \(D_o=0\), the profile \(Q_{\alpha,d}\) is computed from 39 ; the finite set \(\mathcal{A}\subset\mathcal{D}_{\mathcal{W}}\) in 1 approximates the supremum in Proposition [prop::Q95sog].

The intermediate tolerances in 1 enter only in the preliminary determination of the step size and truncation interval. The final SOG parameters are required to satisfy \(E_{\mathcal{A}}\le\epsilon\) on the adaptively refined set \(\mathcal{A}\). For \(\alpha=1/2\) and \(D_o=0\), \(Q_{\rm ref}\) is the closed-form Cauchy profile; in the remaining cases, it is computed from the scaled integral representation by high-accuracy quadrature. This finite-set comparison is used for the domain-adapted parameter choice, while the continuum a priori estimate of 3.3 provides the rigorous error bound for the pure-fractional case \(D_o=0\). In high dimensions, both \(Q_{\rm ref}\) and \(Q_{h,L,U}\) are evaluated in logarithmic form, and the largest term is factored out of the finite SOG sum to preserve numerical stability.

Closed form at \(\alpha=1/2\). The preceding construction does not rely on a closed form for \(\rho_\alpha\). When \(\alpha=1/2\) and \(D_o=0\), however, the one-sided stable density is elementary, \(\rho_{1/2}(s)=(2\sqrt{\pi})^{-1}s^{-3/2}e^{-1/(4s)}\), and the discretization error of the infinite trapezoidal rule can be written explicitly. These formulas give a closed-form reference for the general parameter-selection criterion.

In this closed-form case, a sufficient all-space choice of \(h\) is the largest value for which \[\label{eq::h95sharp} 2\sum_{m\ge1}\frac{\big|\Gamma\!\big(\tfrac{d+1}{2}+i\,2\pi m/h\big)\big|}{\Gamma\!\big(\tfrac{d+1}{2}\big)} \le\epsilon\tag{45}\] holds. By Lemma 4, \(|\Gamma(\nu+i\tau)|\sim\sqrt{2\pi}\,(\nu^2+\tau^2)^{(2\nu-1)/4}e^{-\pi|\tau|/2}\) decays exponentially in \(|\tau|\). The leading-mode approximation, combined with the Gaussian approximation of the Gamma ratio for \(2\pi/h\ll (d+1)/2\), gives the large-\(d\) estimate \[\label{eq::h95sharp95asymp} h\approx\frac{2\pi}{\sqrt{(d+1)\ln(2/\epsilon)}}.\tag{46}\] The integrand peaks at \(s_\star=\beta_Y/(2(d+1))\) and decays as \(\exp(-\tfrac{d+1}2\psi(v))\), where \(\psi(v)=v+e^{-v}-1\) and \(v=\log(s/s_\star)\). Equating this decay factor to \(\epsilon\) gives the retained band \[\label{eq::band} s_\ell=e^{h\ell}\in\Big[\frac{e^{v_{\mathrm{lo}}}}{2(d+1)},\; \frac{(1+\eta_{\max}^2)\,e^{v_{\mathrm{hi}}}}{2(d+1)}\Big],\qquad \eta_{\max}=\frac{R}{(D_fT_{\min})^{1/(2\alpha)}},\tag{47}\] where the endpoints \(v_{\mathrm{lo}}<0\) and \(v_{\mathrm{hi}}>0\) solve \[\label{eq::vlo95vhi} \psi(v_{\mathrm{lo}})=\psi(v_{\mathrm{hi}})=q_d:=\frac{2\log(1/\epsilon)}{d+1}.\tag{48}\] The corresponding term count is \[\label{eq::nexp} M=M_1+M_2+1=\frac{\ln(1+\eta_{\max}^2)+(v_{\mathrm{hi}}-v_{\mathrm{lo}})}{h}.\tag{49}\]

The complete FFPE solver is summarized in 2.

4 Numerical results↩︎

In this section we assess the fast SOG solver of 2 and study its accuracy, robustness, and computational cost. All experiments are carried out in MATLAB R2025b on a laptop with an Intel Core Ultra 7 255H CPU and 64 GB of memory, with a serial implementation. As an exact reference we use the only nontrivial case for which the fundamental solution is known in closed form in every dimension: the pure fractional case \(D_o=0\), \(\alpha=1/2\), for which 17 evaluates to the \(d\)-dimensional Cauchy distribution [19] \[\label{eq::cauchy95ref} p(\boldsymbol{x},T)=\frac{\Gamma\!\left(\frac{d+1}{2}\right)}{\pi^{\frac{d+1}{2}}} \frac{D_fT}{\big[(D_fT)^2+y^2\big]^{\frac{d+1}{2}}}, \qquad y=|\boldsymbol{x}-\boldsymbol{x}_0^T|.\tag{51}\] This reference is exact to machine precision for arbitrary \(d\), and is therefore well suited to assessing the solver in the high-dimensional and small-time regimes that are otherwise the most difficult to compare against closed-form references.

All reported term counts are the consecutive bands returned by 1. A run is fully specified by \((\alpha,d,D_o,D_f,\epsilon)\), the physical window in \(y\) and \(t\), and the resulting parameters \((h,L,U)\). The SOG weights are then \(w_\ell=h e^{h\ell}\rho_\alpha(e^{h\ell})\), and all high-dimensional sums and references are evaluated through logarithms using 52 and log-Gamma functions. For \(\alpha=1/2\) errors are measured against the exact Cauchy density. For \(\alpha\ne1/2\) the reported errors are the solution-level self-convergence estimates described in 4.2.

4.1 High-order accuracy and convergence↩︎

We first verify that the SOG solver attains a prescribed accuracy in moderate and high dimension. Fixing \(D_f=1\), \(D_o=0\), \(\alpha=1/2\), and \(T=1\), we prescribe a relative accuracy \(\epsilon\), select the SOG parameters by 1 (with the step \(h\) initialized by 45 and the band by 47 49 ), and measure the relative \(L_\infty\) error of \(p_{\text{SOG}}\) against the exact Cauchy density 51 over \(y\in[0,2]\). The dimensions in 3 are \(d=7,15,29,100\). Panel (a) shows that the measured error tracks the prescribed tolerance with slope one across eleven orders of magnitude, reaching \(7.8\times10^{-14}\) at \(d=100\) when \(\epsilon=10^{-13}\). The accuracy is therefore limited by the prescribed tolerance and, ultimately, by the arithmetic of the double-precision evaluation. Two further observations confirm the solution-level parameter selection of 3.4:

• The number of Gaussians \(M\) grows only logarithmically in \(1/\epsilon\) (3 (b)). Across all dimensions and tolerances shown, \(M\) lies between \(8\) and \(56\); at \(\epsilon=10^{-13}\) the term counts are \(56,44,39,39\) for \(d=7,15,29,100\), respectively. The dimension enters through both the step \(h\propto(d+1)^{-1/2}\) (asymptotically, 46 ) and the truncation band 47 49 . The decrease over this range is the low-to-moderate-dimensional branch of the non-monotonic dependence described in Remark [rmk::nexp95d]; in very high dimension the displacement-band term \(\ln(1+\eta_{\max}^2)/h\propto\sqrt d\) dominates and \(M\) grows again (4.4).

• Replacing the evaluation of \(\rho_\alpha\) by the exact closed form of the one-sided stable (Lévy) density \(\rho_{1/2}(s)=\tfrac{1}{2\sqrt{\pi}} s^{-3/2}e^{-1/(4s)}\) reproduces the errors of 3 to every displayed digit. The density evaluation therefore does not limit the accuracy; the residual error is due entirely to the SOE truncation.

For \(\alpha=1/2\), Proposition [prop::disc] controls the relative discretization error independently of the evaluation point and of \(T\), and 1 enforces \(E_{\mathcal{A}}\le\epsilon\) for the finite retained band on the prescribed spatial window; by Proposition [prop::Q95sog], this scaled relative error equals that of the physical Green’s function, so the guarantee transfers directly to \(p_{\text{SOG}}\) and the measured error tracks \(\epsilon\) uniformly in \(d\).

4.2 Self-convergence study↩︎

The exact Cauchy reference 51 is available only for \(\alpha=1/2\). For other fractional orders we assess the assembled solution by self-convergence: the chosen finite Gaussian sum is compared, on the scaled radial interval, with a refined reference at half the step and a wider band. Since the inverse Fourier transform of each Gaussian is exact, this probes the only numerical approximation in the method, the SOE approximation of \(e^{-x^\alpha}\) on the active solution window. These entries should therefore be interpreted as finite-domain self-convergence estimates; the independent exact-reference checks are the \(\alpha=1/2\) row and the Cauchy tests in the other tables.

1 reports the result at \(d=1000\) for \(\alpha=0.1,0.3,0.5,0.7,0.9\), with \(D_o=0\), \(D_f=8\), \(t=0.04\), and \(y\in[0,2]\), at prescribed tolerance \(\epsilon=10^{-10}\). It also lists the scaled endpoint \(\eta_{\max}=R/(D_ft)^{1/(2\alpha)}\), since this – not \(t\) alone – sets the active range of Gaussian scales. The largest term count occurs at \(\alpha=0.1\), where the fixed window maps to the largest scaled radius; the \(\alpha=0.3\) row, in which a saddle-point approximation of \(\rho_\alpha\) is blended smoothly with the direct evaluation at small \(s\), indicates that the density evaluation is not the limiting factor in this test. The \(\alpha=1/2\) row, where the true error against 51 is available, also meets the tolerance and anchors the self-convergence estimates at the other orders.

4.3 Robustness across time and dimension↩︎

We now fix a relative tolerance \(\epsilon=10^{-12}\) and, for each \((t,d)\), select the SOG parameters from 1 (with the \(\alpha=1/2\) step from 45 and band from 47 49 ). We evaluate \(p_{\text{SOG}}\) for \(D_o=0\), \(D_f=8\), \(\alpha=1/2\) on the higher-dimensional part of the \(t\times d\) grid used in [19], reporting the maximum relative error over \(y\in[0,2]\).

Two features stand out. First, the time dependence is handled by self-similar scaling (Proposition [prop::disc] and Remark [rmk::selfsim]): a single approximation sized at \(t_{\min}=0.004\), using \(M=71\)–\(82\) terms over \(d=5\)–\(29\), holds relative error about \(10^{-12}\) across the whole time grid, with no degradation as \(t\to0^+\) (the regime in which the solution concentrates toward the Dirac measure). Second, the error is uniform in dimension, staying at the prescribed \(10^{-12}\) throughout. 2 contrasts this with the integral solver of [19]: at small time (\(t=0.004\)) its relative error grows rapidly over the tested dimensions, exceeding unity for \(d\ge17\) and reaching \(\mathcal{O}(10^{9})\) at \(d=29\), whereas at moderate time (\(t=0.2\)) it attains machine precision. The two approaches are thus complementary: the integral solver of [19] is effective for high-precision evaluation in low to moderate dimension, while the SOG holds the prescribed accuracy uniformly across the tested \((t,d)\) range – including the small-\(t\), high-\(d\) corner where the integral solver loses accuracy, as shown in 4 – and scales to much higher dimensions (4.4).

4.4 Scaling to very high dimension↩︎

Unlike the oscillatory radial integrand \(r^{d/2}J_{(d-2)/2}(yr)\) that limits [19] to moderate \(d\), the SOG solution is a positive sum of separable Gaussians, so its logarithm can be evaluated stably. Setting \(a_\ell(\boldsymbol{x})=\log w_\ell-\tfrac{d}{2}\log(4\pi C_\ell^{T}) -|\boldsymbol{x}-\boldsymbol{x}_0^T|^2/(4C_\ell^{T})\) and \(a_\star=\max_\ell a_\ell(\boldsymbol{x})\), \[\label{eq::logsumexp} \log p_{\text{SOG}}(\boldsymbol{x},T)=a_\star+\log\sum_{\ell} \exp\!\big(a_\ell(\boldsymbol{x})-a_\star\big),\tag{52}\] so that the per-term magnitudes \((4\pi C_\ell^{T})^{-d/2}\), which overflow or underflow for large \(d\), never appear explicitly. Each evaluation requires \(O(M)\) operations for the radial profile.

These dimensions exceed those reached in the radial-quadrature experiments of [19]. We again use the exact Cauchy solution 51 : although its magnitude over- or underflows, its logarithm is computable to full relative precision in any dimension via the log-Gamma function. 3 reports the relative error and evaluation time per point for \(d\) up to \(10^{5}\). The SOG attains ten-digit relative accuracy at \(d=10^{5}\) in about \(4\) μs per point, with the term count \(M\) rising only from \(36\) at \(d=100\) to \(451\) at \(d=10^5\) (the slow large-\(d\) growth is analyzed in Remark [rmk::nexp95d]). These tests give deterministic, high-accuracy evaluations at dimensions beyond the reach of radial-quadrature methods.

4.5 General initial data: tensor-product representations↩︎

A key advantage of the SOG representation is that it provides a reusable building block that maps separated (tensor-product) data to separated data. The SOG approximation \(p_{\text{SOG}}\) is itself a rank-\(M\) separated function: by 22 each term \(p_\ell\) is a tensor product \(\prod_{j=1}^{d} g_\ell(x_j)\) of one-dimensional Gaussians. Consequently, for any initial datum written as a linear combination of tensor products, \[\label{eq::tensor95ic} p(\boldsymbol{x},0)=\sum_{r=1}^{R}\prod_{j=1}^{d}\phi_{r,j}(x_j),\tag{53}\] linearity and the separability of the heat kernel give the solution in closed form as \[\label{eq::tensor95sol} p(\boldsymbol{x},T)=\sum_{r=1}^{R}\sum_{\ell=-M_1}^{M_2}w_\ell\prod_{j=1}^{d} \big(g_\ell * \phi_{r,j}\big)(x_j-b_jT),\tag{54}\] which is again a sum of tensor products. The dimension \(d\) enters only through the one-dimensional convolutions \(g_\ell*\phi_{r,j}\), and the rank grows from \(R\) to \(RM\), recompressible by standard tensor-rank truncation. Such separated, low-rank formats – the canonical, Tucker, and tensor-train decompositions [21]–[23], [39], [40] and the related tensor networks [41]–[43] – provide a natural high-dimensional function class for this solver. Tensor neural networks [26]–[28], [44] can be used to construct separated approximations for more general initial data.

Gaussian initial data are the canonical example: the one-dimensional convolutions remain Gaussian and are available analytically. Since each \(p_\ell\) in 22 has covariance \(2C_\ell^{T}I\), the solution for the initial condition \[p(\boldsymbol{x},0)=\sum_{j=1}^{N_g}a_j\, \mathcal{N}(\boldsymbol{x};\boldsymbol{c}_j,\sigma_j^2 I)\] is \[\label{eq::sog95ic} p(\boldsymbol{x},T)=\sum_{j=1}^{N_g}\sum_{\ell=-M_1}^{M_2}a_j\,w_\ell\, \mathcal{N}\!\big(\boldsymbol{x};\,\boldsymbol{c}_j+\boldsymbol{b}T,\,(\sigma_j^2+2C_\ell^{T})I\big),\tag{55}\] a sum of \(N_gM\) Gaussians, each a tensor product across the \(d\) coordinates, stored and manipulated in the same factored \(O(MdN)\) representation as the fundamental solution.

We test 55 directly in high dimension. For a single Gaussian source, the value at the advected center has the independent non-oscillatory reference \[\frac{\Omega_d}{(2\pi)^d}\int_0^\infty r^{d-1}\exp\!\left(-\frac{\sigma^2r^2}{2}-D_fT r^{2\alpha}\right)\mathrm{d}r,\] which is evaluated after a saddle-point change of variables and does not involve the oscillatory Bessel integral. With \(\alpha=1/2\), \(\sigma=0.5\), \(D_f=1\), and a broad \(M=3251\) term SOE that preserves \(\sum_\ell w_\ell=1\) to the displayed digits, the relative errors at \(T=0.1\) are \(1.1\times10^{-13}\), \(2.7\times10^{-12}\), and \(6.2\times10^{-11}\) for \(d=10^3,10^4,10^5\), respectively; at \(T=1\) they are \(3.4\times10^{-13}\), \(9.1\times10^{-13}\), and \(9.8\times10^{-11}\).

5 (a) shows a two-dimensional slice of the \(d=1000\) solution for a non-radial, three-Gaussian initial condition under drift and pure fractional diffusion, computed from the separated representation directly rather than by reducing the data to a radial profile.

Finally, 6 confirms the linear-in-\(d\) cost. Holding the SOG approximation fixed at \(M=1012\) terms (\(P=N_gM=8096\) Gaussians for \(N_g=8\) sources) and varying only the dimension, the factored assembly time is indistinguishable from a straight line through the origin, reaching \(0.38\) s at \(d=60\). In \(d=10\) the same \(P=8096\) factored Gaussians are assembled in \(0.07\) s, and a factored evaluation agrees with the explicit double sum over \(j\) and \(\ell\) to \(8\times10^{-16}\) – whereas the corresponding dense tensor has \(N^{10}\approx10^{18}\) entries and is not a feasible object to store.

Combining this building block with the logarithmic evaluation of 4.4 lets us evolve sum-of-Gaussians initial data in very high dimension; 5 (b) shows a six-Gaussian initial condition evolved in \(d=1000\) with drift along the line of centers. Spanning more than eighty orders of magnitude, the solution is representable only through its logarithm 52 . The mass \(\int_{\mathbb{R}^d}p\,\mathrm{d}\boldsymbol{x} =\big(\sum_j a_j\big)\big(\sum_\ell w_\ell\big)\) is conserved to \(1.000000\) at every time, in any dimension, and evaluating one 600-point time slice costs about \(0.09\) s at \(d=1000\).

5 Conclusion↩︎

We have developed a sum-of-Gaussians solver for the high-dimensional fractional Fokker–Planck equation with high-order accuracy and error control. The construction uses the complete monotonicity of the fractional symbol \(e^{-D_ft|\boldsymbol{k}|^{2\alpha}}\), which represents it as a continuous superposition of Gaussians. Discretizing that superposition gives a finite sum of ordinary heat flows that act independently in each coordinate, with storage and assembly cost linear in the dimension and quadrature parameters fixed a priori for any target accuracy. In our experiments the solver attains more than ten digits of accuracy, with \(M\) growing only logarithmically as the tolerance tightens and the accuracy uniform across the time windows tested; evaluating the positive Gaussian form through its logarithm carries the computation to dimension \(10^{5}\).

Because the method approximates the fundamental solution by a separated sum of Gaussians, it propagates any tensor-product initial datum in closed form, so the only additional ingredient needed for more general data is a separated, low-rank representation of it. Sparse-grid and tensor neural network approximation frameworks are compatible with this requirement and can be combined directly with the present solver. When a high-dimensional Fourier symbol or interaction kernel admits an accurate Gaussian-sum expansion, the expansion expresses nonseparable terms as sums of separable Gaussian contributions. The FFPE results presented here therefore support the use of SOG approximations as building blocks for separable representations of kernels and solution operators in high-dimensional PDEs.

The approach is also relevant beyond the particular FFPE fundamental solution studied here. For a linear constant-coefficient evolution equation in high dimension, if the Fourier-space propagator admits an accurate Gaussian-sum approximation, the same construction expresses the solution operator as a sum of separable Gaussian evolution operators. For many nonlinear evolution equations, applying an unconditionally energy-stable scalar auxiliary variable (SAV) temporal discretization reduces each time step to a linear problem with known source terms [33]. When this linear problem has a constant-coefficient solution operator that admits an accurate and efficient Gaussian-sum representation, the present solver can serve as a building block for high-dimensional nonlinear evolution PDEs.

Acknowledgments↩︎

The Flatiron Institute is a division of the Simons Foundation. The work of Q. Zhou was supported by the National Natural Science Foundation of China (Grant No. 125B2023).

6 A bound for \(|\rho_\alpha(z)|\) on \(|z|\le 1\)↩︎

In this section, we provide some estimates on \(|\rho_\alpha(z)|\), where \(\rho_\alpha(z)\) is the one-sided stable PDF with \(0<\alpha<1\), analytically continued to the complex plane. Denote \(z=re^{i\theta}\), where the modulus \(r\in [0,1]\) and the angle \(|\theta|<\theta_{*}(1-\alpha)\), with \(0<\theta_{*}<\pi/6\) a critical angle to be determined later. The inverse Laplace transform gives \[\label{eq::ILT} \rho_\alpha(z)=\frac{1}{2\pi i}\int_{C}e^{\Phi_{z}(s)}\mathrm{d}s,\tag{56}\] where \[\label{eq::phase95original} \Phi_{z}(s)=sz-s^{\alpha},\tag{57}\] and \(C\) is an appropriate inversion contour. To estimate \(|\rho_\alpha(z)|\), we proceed in the following steps:

1. Find the saddle point of \(\Phi_z(s)\). Transform \(C\) into the (approximated) steepest descent contour. Specifically, a parabolic contour derived from the local and global properties of \(\Phi_z(s)\).

2. Analyze properties of the phase function on the parabolic contour.

3. Based on the results above, provide an upper bound of \(|\rho_\alpha(z)|\).

We proceed according to these steps in sequence.

6.1 The parabolic contour↩︎

From 57 , we take derivative of \(\Phi_z(s)\), then \[\label{eq::Phi95deri} \Phi_z'(s)=z-\alpha s^{\alpha-1},\tag{58}\] which shows the saddle point is \(s_c=(\alpha/z)^{1/(1-\alpha)}\). The polar representation of \(s_c\) is provided by \(s_c=\sigma e^{i\psi_0}\), where \[\label{eq::sigma95psi} \sigma=\alpha^{1/(1-\alpha)}r^{-1/(1-\alpha)},\quad \psi_0=-\frac{\theta}{1-\alpha}.\tag{59}\] Thus one has \(\Phi_z(s_c)=-(1-\alpha)s_c^{\alpha}\) and \(\Phi_z''(s_c)=\lambda s_c^{\alpha-2}\), \(\lambda=\alpha(1-\alpha)\). Based on these results, one selects the parameterized contour as \(s(u)=s_cw(u)\), where \[\label{eq::sd95quad} w(u)=1+iu-\eta u^2,\quad \eta=\frac{2-\alpha}{6}\in(\frac{1}{6},\frac{1}{3}),\quad u\in\mathbb{R}.\tag{60}\] Correspondingly, the phase and the parameter derivative have the expression as \[\label{eq::phase95para} \Phi_z(s(u))=\sigma^\alpha e^{i\phi} (\alpha w(u) - w(u)^\alpha), \quad s'(u)=\sigma e^{i\psi_0}(i-2\eta u).\tag{61}\] where \(\phi=\alpha\psi_0\). Hence the integral in 56 becomes \[\label{eq::integral95transform} \rho_\alpha(z)=\frac{1}{2\pi i}\int_{-\infty}^{\infty} \exp\Big(\sigma^\alpha e^{i\phi}(\alpha w(u) - w(u)^\alpha)\Big) e^{i\psi_0}\sigma\bigl(i-2\eta u\bigr)\mathrm{d}u.\tag{62}\] Define the main phase function as \[\label{eq::main95phase} F(u):=\operatorname{Re}\{e^{i\phi}(\alpha w(u)-w(u)^{\alpha})\},\tag{63}\] and one readily derives the upper bound \[\label{eq::rho95bound} |\rho_\alpha(z)|\le\frac{1}{2\pi}\int_{-\infty}^{\infty} \exp\Big(\sigma^\alpha F(u)\Big)\cdot \sigma(1+2\eta |u|)\mathrm{d}u.\tag{64}\]

6.2 The global maximum property of \(F(u)\)↩︎

For the main phase function \(F(u)\) defined in 63 , Proposition [prop::global95maximum] characterizes its most essential property.

By the symmetry \(F(-u;\phi)=F(u;-\phi)\), it suffices to prove \(F'(u)<0\) for \(u>0\) and all admissible \(\phi\). The case \(u<0\) then follows by replacing \(\phi\) with \(-\phi\). Before proving the global maximum property, we first consider the behavior of \(F(u)\) as \(|u|\rightarrow0\) and \(|u|\rightarrow \infty\), described in Lemmas 5 and 6.

Now we introduce the polar representation \(w(u):=\varrho(u)e^{i\delta(u)}\), with \(\varrho(u)=|w(u)|\) and \(\delta(u)=\arg w(u)\in (0,\pi)\). Lemma 7 shows that the argument of \(w(u)\) is monotonically increasing for \(u>0\).

Furthermore, with the polar representation of \(w(u)\), the derivative of \(F(u)\) is then \[\label{eq::F95derivative} F'(u)=\operatorname{Re}\Big\{\alpha e^{i\phi}(1-w(u)^{\alpha})w'(u)\Big\}.\tag{65}\] Hence we denote \[\label{eq::tangent} k(u)=2\eta u>0,\quad \tau(u)=\arctan(k(u))\in[0,\frac{\pi}{2}),\quad S(u)=\phi+\tau(u)\tag{66}\] as the positive tangent slope, tangent angle, and the transport angle of \(w'(u)\), respectively, since \(w'(u)=i-2\eta u\). We also introduce auxiliary exponents of modulus and argument that \[\label{eq::auxiliary} \mu(u)=\varrho(u)^{\alpha-1},\quad \beta(u)=(1-\alpha)\delta(u)\in (0,\pi),\tag{67}\] which indicates that \[\label{eq::F95expression} \frac{F'(u)}{\alpha\sqrt{1+k(u)^2}} =\mu(u)\sin\big(S(u)-\beta(u)\big)-\sin S(u).\tag{68}\] Subsequent proofs concerning monotonicity will rely on the expression provided in 68 , requiring a detailed analysis of the ranges and quantitative relationships of \(\mu(u)\), \(S(u)\), and \(\beta(u)\). Lemma 8 provides the bound of \(\varrho(u)\) and \(\mu(u)\).

For the transport angle \(S(u)\), Lemmas 9 and 10 show the geometry of the case when \(S(u)\ge \pi/2\).

Based on all the properties discussed above, we refine the angular constraints of Proposition [prop::global95maximum] and provide a proof that \(F(u)\) has a unique global maximum at \(F(0)\). The detailed description is given in Theorem 3.

Proving Theorem 3 requires a detailed discussion of the phase relationship between the two angles \(S(u)\) and \(\beta(u)\). Recall that \(S(u)=\phi+\tau(u)\in[-\theta_{*},\pi/2+\theta_{*})\) and \(\beta(u)=(1-\alpha)\delta(u)\in (0,\pi)\). For convenience, we split the whole phase diagram into four parts, which is clearly shown in 7. Specifically, the four cases are:

• Case A. \(S(u)< 0\).

• Case B. \(0\le S(u)\le \pi/2\).

• Case C. \(S(u)>\pi/2\), and \(\beta(u)\ge\pi/2\).

• Case D. \(S(u)>\pi/2\), and \(0<\beta(u)<\pi/2\).

The following Propositions [prop::A]-[prop::D] provide the whole proof process of Theorem 3.

6.3 An upper bound of \(|\rho_\alpha(z)|\)↩︎

In 6.2, we have shown that \(F(0)\) becomes the strict global maximum of the main phase function \(F(u)\). Here we provide an upper bound of \(|\rho_\alpha(z)|\) based on estimate 64 , described in Theorem 4.

Straightforward manipulations lead to Corollary 1, which gives the estimate on the contour selected in 3.1.

7 Integral bound of \(I_\alpha\)↩︎

To estimate the error of SOE approximation in Theorem 1, an upper bound of \(I_\alpha\) is required with the expression of \[\label{eq::I95alpha95new} I_\alpha=\int_{0}^{\infty}|\rho_\alpha(re^{i\theta})|\mathrm{d}r,\tag{96}\] where \(\theta=-(1-\alpha)\theta_{*}\). The case \(\theta=(1-\alpha)\theta_{*}\) has the same value by the reflection property \(\rho_\alpha(\bar z)=\overline{\rho_\alpha(z)}\). 6 provides estimates for \(|\rho_\alpha(z)|\) within \(|z|\le 1\). We therefore split the integral into \[\label{eq::I95alpha95split} I_\alpha=\int_{0}^{1}|\rho_\alpha(re^{i\theta})|\mathrm{d}r+\int_{1}^{\infty}|\rho_\alpha(re^{i\theta})|\mathrm{d}r:=I_\alpha^1+I_\alpha^2,\tag{97}\] and we provide the bound for each component.

For the compact part \(I_\alpha^{1}\), the pointwise bound of Corollary 1 is too conservative when \(\alpha\) is close to one. It controls \(|\rho_\alpha(z)|\) for each fixed \(z\) by forcing all lower-order terms in the phase to be dominated by the quadratic term, thereby introducing the spurious factor \(c^{1/(1-\alpha)}\). For the integral defining \(I_\alpha^1\), a sharper route is to integrate first in \(r\) and only then estimate the saddle-contour integral. We use the pointwise estimate only for the complementary range \(0<\alpha\le 1/2\), where its constants grow only algebraically.

For the infinite part \(I_\alpha^{2}\), we shall use the asymptotic properties of \(\rho_\alpha(z)\) as \(|z|\ge 1\). Indeed, we have the following estimate on \(|\rho_{\alpha}(z)|\) when \(|z|>1\), as shown in Lemma 12.

The following Corollary 2 follows from integrating both sides of ?? , which provides the error bound of \(I_{\alpha}^2\).

Since \(\Gamma(\alpha+1)\) is bounded on \((0,1)\) and \(2^{1-\alpha}-1\ge c(1-\alpha)\) for \(0<\alpha<1\), ?? implies \[I_\alpha^2\le C\left(\frac{1}{\alpha}+\frac{1}{1-\alpha}\right).\] Combining this estimate with Theorem 5 proves ?? .

References↩︎

---

1. Center for Computational Mathematics, Flatiron Institute, Simons Foundation, New York, NY 10010, USA ().↩︎

2. School of Science and Engineering, The Chinese University of Hong Kong, Shenzhen, Shenzhen, Guangdong 518172, P. R. China; Shenzhen International Center for Industrial and Applied Mathematics, Shenzhen Research Institute of Big Data, Shenzhen, Guangdong 518172, P. R. China ().↩︎

3. School of Mathematical Sciences, Shanghai Jiao Tong University, Shanghai 200240, P. R. China ().↩︎