Human-level Control Through Deep Reinforcement Learning

ArXiv: 1406.3722

🎯 Pitch

When the Laplacian in a steady-state Helmholtz equation is replaced by a Riesz-Feller derivative in one variable and a Hilfer-composite derivative in the other, the analytical solution transmutes into a Fox H-function whose asymptotic decay crosses over from stretched exponential to power-law depending on the fractional orders. This paper shows exactly how the interpolation parameter of the Hilfer derivative controls that crossover, providing a unified Green’s function for anomalous diffusion with mixed Lévy and memory-type dynamics.

---

1. Executive Summary

This paper derives analytical solutions to fractional-order generalizations of the Laplace, Poisson, and Helmholtz equations in two variables, obtained by replacing integer-order partial derivatives with the Riesz-Feller fractional derivative (a pseudo-differential operator whose symbol is the logarithm of a Lévy stable distribution's characteristic function) and the Hilfer-composite fractional derivative (a generalized Riemann-Liouville derivative parameterized by order μ and type ν, interpolating between classical R-L and Caputo derivatives). Using Fourier-Laplace transform methods, the authors express solutions in terms of Mittag-Leffler functions, Fox H-functions, and Prabhakar integral operators containing Mittag-Leffler kernels, with detailed asymptotic expansions and series representations provided for limiting regimes where |x|/y^(μ/2) is either large or small. The general space-time fractional wave equation with Riesz-Feller space derivative and Hilfer-composite time derivative is also treated, establishing that the steady-state fractional Helmholtz and Poisson solutions can be mapped directly onto time-dependent wave problems—but only through the negative-sign correspondence between the quantum and classical forms of the Riesz-Feller operator, not by independent derivation.

2. Context and Motivation

The Core Gap: No Analytical Solutions for Fractional Steady-State Equations with Mixed Fractional Derivatives

The paper addresses a specific mathematical gap: while fractional differential equations had been extensively studied for time-dependent problems — fractional diffusion, fractional wave equations, fractional Fokker-Planck equations — there was no systematic treatment of steady-state fractional equations in two independent spatial variables where the fractional derivatives are of different types in each variable. The authors set out to solve equations of the form

$xD_\theta^\alpha N(x,y) + yD_{0+}^{\mu,\nu} N(x,y) = \Phi(x,y),$

where the $x$-derivative is a Riesz-Feller fractional derivative (a pseudo-differential operator with a Lévy-stable symbol) and the $y$-derivative is a Hilfer-composite fractional derivative (interpolating between Riemann-Liouville and Caputo types via a parameter $\nu \in [0,1]$). No closed-form solutions for this class of mixed-derivative fractional equations existed in the literature prior to this work.

This is not merely a cataloging exercise. The structure of the derivative operators matters profoundly: the Riesz-Feller derivative captures Lévy flight behavior (heavy-tailed spatial jumps) governed by a stability index $\alpha$ and skewness parameter $\theta$, while the Hilfer derivative captures memory effects in the $y$-direction parameterized by the order $\mu$ and the interpolating type $\nu$. Solving the combined equation analytically means understanding how Lévy-flight spatial transport interacts with memory-type behavior in the other variable — a genuinely new mathematical problem.

Why This Matters: Physical Context and the Steady-State/Time-Dependent Duality

The paper is motivated by a confluence of physical applications and a structural observation about fractional equations.

Physical motivation. Fractional Laplace and Poisson equations arise as steady-state limits of time-dependent anomalous diffusion or wave propagation. Specifically, when a time-fractional diffusion or wave equation reaches equilibrium (the time derivative vanishes), the spatial operator that remains is precisely the fractional operator studied here. The physical contexts cited in Section I include:

• Non-exponential relaxation in glassy materials and complex fluids, where Hilfer showed that the Hilfer-composite time derivative captures the transition from stretched exponential to power-law relaxation in glass-forming materials (Hilfer, 2002). The steady-state version of these models inherits the Hilfer derivative structure in one spatial variable.
• Lévy flights in anomalous transport, where the Riesz-Feller space derivative emerges naturally from continuous-time random walk theory with heavy-tailed jump length distributions (Metzler and Klafter, 2000, 2004). Processes ranging from particle transport in turbulent plasmas to animal foraging patterns exhibit Lévy flight statistics.
• Quantum mechanics, where Luchko et al. (2013) introduced a fractional Schrödinger equation using the quantum Riesz-Feller derivative for a free particle and a particle in a potential well. The steady-state solutions of these equations are fractional Helmholtz problems of the type studied here.
• Vibrations of smart materials with memory effects (Liang et al., 2005; Liang and Chen, 2006), where the fractional wave equation describes viscoelastic damping that standard integer-order models cannot capture. The Helmholtz equation with fractional derivatives corresponds to the frequency-domain version of these vibration problems.

The paper's equations therefore serve as canonical models for a broad class of physical systems where anomalous transport competes with or couples to memory effects in orthogonal directions.

Structural duality. A key observation that motivates the paper's organization is that the steady-state fractional equations (Laplace, Poisson, Helmholtz) are structurally identical to space-time fractional wave equations under a sign change of the Riesz-Feller operator. Specifically, the fractional wave equation (57)

$tD_{0+}^{\mu,\nu} N(x,t) = xD_\theta^\alpha N(x,t) + \Phi(x,t)$

(using the classical Riesz-Feller derivative $xD_\theta^\alpha$) maps directly onto the steady-state equation (23)

$xD_\theta^\alpha N(x,y) + yD_{0+}^{\mu,\nu} N(x,y) = \Phi(x,y)$

if and only if the classical Riesz-Feller derivative in the wave equation is related to the quantum Riesz-Feller derivative in the steady-state equation via $xD_\theta^\alpha = - xD_\theta^{*,\alpha}$. This sign relationship is not arbitrary — it reflects the fact that the Fourier symbol of the classical Riesz-Feller derivative is $-\psi_\alpha^\theta(\kappa)$ while the quantum version has symbol $+\psi_\alpha^\theta(\kappa)$. The authors leverage this to avoid deriving the wave equation solutions independently: all the analytical machinery developed for the steady-state case transfers directly to the time-dependent case by swapping the sign in the Mittag-Leffler argument. This duality is conceptually important because it means the steady-state solutions are not just special cases — they are the central objects from which time-dependent solutions follow.

Prior Approaches and Their Limitations

Fractional diffusion and wave equations were well-studied, but only with a single fractional derivative type.

The literature on fractional differential equations was mature at the time of this paper. The time-fractional diffusion equation with a Riemann-Liouville or Caputo time derivative and an ordinary Laplacian in space had been extensively solved by Mainardi (1996, 2004), Metzler and Klafter (2000), and others, with solutions expressed in terms of Wright functions and Fox H-functions. The space-fractional diffusion equation with a Riesz-Feller space derivative and an ordinary first-order time derivative had been solved by Mainardi, Pagnini, and Saxena (2005), who expressed the fundamental solution as a Fox H-function. These results established the transform methodology (Fourier in space, Laplace in time) that this paper extends.

However, these prior works dealt with one fractional derivative at a time — either the time derivative was fractional and space was integer-order, or space was fractional and time was integer-order. No one had combined a Riesz-Feller derivative in one spatial variable with a Hilfer derivative in the other, which is the central problem of this paper. The complication is not just notational: the Fourier transform handles the Riesz-Feller derivative elegantly (it becomes an algebraic multiplier $-\psi_\alpha^\theta(\kappa)$), while the Laplace transform handles the Hilfer derivative (it becomes $s^\mu$ plus initial-value terms involving the fractional integral at the origin). The solution requires a sequential application of both transforms, and the resulting inverse transforms involve Mittag-Leffler functions whose arguments contain the Riesz-Feller symbol — a class of integrals that requires Fox H-function techniques to evaluate explicitly.

Hilfer-composite derivative had been used, but only in time-fractional problems.

The Hilfer-composite derivative (8) was introduced by Hilfer (2000) as a generalization that interpolates between the Riemann-Liouville derivative ($\nu = 0$) and the Caputo derivative ($\nu = 1$). Its property that the initial conditions involve fractional integrals evaluated at the origin — specifically, $(yI_{0+}^{(1-\nu)(n-\mu)} f)(0+)$ — made it attractive for modeling physical processes where the appropriate initial data are neither purely integer-order derivatives (as in Caputo) nor purely fractional integrals (as in Riemann-Liouville), but something in between. Sandev, Metzler, and Tomovski (2011) used the Hilfer-composite time derivative to solve a time-fractional diffusion equation with an ordinary spatial Laplacian. Tomovski, Sandev, Metzler, and Dubbeldam (2012) extended this to a space-time fractional diffusion equation with a Riesz-Feller space derivative and a Hilfer-composite time derivative, obtaining solutions in terms of Fox H-functions.

But none of this prior work considered the case where the Hilfer derivative acts in a spatial variable rather than in time. In the steady-state equations (Laplace, Poisson, Helmholtz), the variable $y$ is a spatial coordinate, and the Hilfer derivative $yD_{0+}^{\mu,\nu}$ represents spatial memory or spatial non-locality in the $y$-direction. This is a different physical interpretation from the time-fractional case, and it requires different boundary conditions: instead of initial conditions at $t = 0$, one prescribes the fractional integrals of the solution and its $y$-derivative at $y = 0+$, which the paper calls "boundary conditions" (24a). The mathematical structure is similar to the time-fractional case (hence the same Laplace transform technique works), but the physical meaning and the resulting solution forms — particularly the asymptotic behavior as $y \to \infty$ for fixed $x$ — are distinct.

Fractional Helmholtz equations had been touched on, but only for Riemann-Liouville derivatives.

The closest prior work is Samuel and Thomas (2010), who considered a fractional Helmholtz equation with a Riemann-Liouville fractional derivative in the $y$-direction and an ordinary or Riesz-Feller derivative in the $x$-direction. Their solutions, which the authors cite as a special case (Remark 8 and Corollary 4, noting that $\nu = 0$ recovers the R-L case), were expressed in terms of Fox H-functions. The limitation of Samuel and Thomas's work is that it only handles the $\nu = 0$ (Riemann-Liouville) case — the parameter $\nu$ that governs the interpolation between R-L and Caputo derivatives is absent. The present paper generalizes this to arbitrary $\nu \in [0,1]$, which is non-trivial because the $\nu$ parameter appears in the powers of $y$ multiplying the Mittag-Leffler functions (the $y^{-(1-\nu)(2-\mu)}$ and $y^{1-(1-\nu)(2-\mu)}$ prefactors) and in the parameters of the Fox H-functions that emerge after the inverse Fourier transform.

Asymptotic analysis of fractional steady-state solutions was missing.

While many prior papers expressed solutions of fractional equations in terms of Fox H-functions, few provided the detailed asymptotic expansions that this paper emphasizes. For the Laplace/Poisson solutions, the authors derive asymptotic behavior in two regimes: $|x|/y^{\mu/2} \gg 1$ (large lateral distance relative to the $y$-scale) and $|x|/y^{\mu/2} \ll 1$ (small lateral distance). The large-argument asymptotics (Remark 6, equation 51) involve stretched exponential decay:

$N(x,y) \propto \exp\left[-\frac{2-\mu}{2}\left(\frac{\mu}{2}\right)^{\frac{\mu}{2-\mu}} \left(\frac{|x|}{y^{\mu/2}}\right)^{\frac{2}{2-\mu}}\right],$

which is qualitatively different from the Gaussian decay of the standard diffusion equation and reveals the interplay between the fractional order $\mu$ and the spatial decay rate. The small-argument expansion (Remark 7, equation 52) is expressed as a Wright function series, which is the natural generalization of the power-series expansion of the integer-order solution. Prior work had not systematically provided these asymptotics for the mixed-derivative case.

How This Paper Positions Itself

The paper frames itself as a unification and generalization of several threads in the fractional calculus literature, not as a radical departure. The abstract is explicit about this:

"Many previously obtained results can be derived as special cases of those presented in this paper."

The positioning strategy is: (1) take the fractional derivatives that have been most useful in separate contexts — Riesz-Feller for Lévy flight spatial transport, Hilfer for interpolated memory — (2) combine them into steady-state equations that are mathematically well-defined via the Fourier-Laplace transform methodology, (3) extract the full analytical structure (Mittag-Leffler representation, Fox H-function closed forms, asymptotic expansions, series representations), and (4) show that the resulting framework contains prior work as limiting cases ($\nu = 0$ for R-L, $\nu = 1$ for Caputo, $\theta = 0$ for symmetric Riesz, $\alpha = 2$ for ordinary Laplacian).

The paper does not claim to introduce new fractional derivative definitions. The Riesz-Feller derivative (1) and the Hilfer-composite derivative (8) are both taken from prior work (Feller, 1968; Hilfer, 2000). The contribution is entirely in the solution methodology and the analytical characterization of the resulting special functions. The authors also do not claim to solve physically new equations — they acknowledge that time-fractional versions with similar operators were studied by Haubold, Mathai, and Saxena (2007), Tomovski et al. (2012), and others. The novelty lies in treating the mixed spatial-fractional steady-state case systematically and in leveraging the sign duality between the quantum and classical Riesz-Feller derivatives to unify the steady-state and time-dependent treatments.

The paper's relationship to the broader fractional calculus community is also worth noting. The authors are explicit that their solution techniques — sequential Fourier and Laplace transforms, Mellin-Barnes integral representations, Fox H-function reduction formulas — are standard tools in the field (as systematized in the monograph by Mathai, Saxena, and Haubold, 2010). The paper does not innovate in methodology; it applies well-established techniques to a previously unsolved class of equations. This is characteristic of the fractional calculus literature of this period: progress came from identifying which combination of fractional operators yields tractable transform-domain algebra and closed-form special function representations, not from inventing new solution methods.

A subtle but important point about the paper's scope: while the title and introduction mention "fractional wave equations," the mathematical heavy lifting is done entirely for the steady-state equations (Sections III and IV), with the wave equation results (Remarks 8 and 10) obtained by direct substitution using the sign relationship between the classical and quantum Riesz-Feller derivatives. The paper does not independently derive time-dependent solutions or analyze their initial-value evolution. This is not a weakness — it is a deliberate structural choice that reflects the physical insight that the steady-state solutions are the mathematically fundamental objects.

3. Technical Approach

3.1 Reader Orientation

This paper constructs a systematic analytical framework for solving a class of fractional-order partial differential equations in two spatial variables where the derivatives are of two fundamentally different types — a Riesz-Feller pseudo-differential operator in the $x$-direction (capturing Lévy-flight behavior) and a Hilfer-composite fractional derivative in the $y$-direction (capturing memory-type non-locality interpolating between Riemann-Liouville and Caputo forms). The problem being solved is: given boundary conditions specified as fractional integrals evaluated at $y = 0^+$ and a source term $\Phi(x,y)$, find the closed-form solution $N(x,y)$ on the upper half-plane that vanishes at infinity in $x$. The solution takes the shape of an inverse Fourier transform of Mittag-Leffler functions whose arguments contain the Riesz-Feller Fourier symbol, which can be further reduced to Fox H-functions for explicit asymptotic analysis.

3.2 Big-Picture Architecture (Diagram in Words)

The solution methodology has six major components connected in a sequential transform pipeline:

1. The fractional PDE — an equation of the form $xD_\theta^\alpha N + yD_{0+}^{\mu,\nu} N = \Phi$ with specified boundary conditions on fractional integrals of $N$ and its $y$-derivative at $y = 0^+$, plus decay at infinity.

2. Laplace transform in $y$ — converts the Hilfer-composite $y$-derivative into an algebraic term $s^\mu \tilde{N}(x,s)$ plus contributions from the initial fractional-integral values at $y = 0^+$. The result is an ordinary differential equation in $x$ parameterized by the Laplace variable $s$.

3. Fourier transform in $x$ — converts the Riesz-Feller $x$-derivative into an algebraic multiplier $-\psi_\alpha^\theta(\kappa)$ (or $+\psi_\alpha^\theta(\kappa)$ for the quantum version). The result is a purely algebraic equation in the doubly-transformed domain $(\kappa, s)$.

4. Algebraic solution in $(\kappa, s)$-space — the doubly-transformed equation is solved by simple algebra, expressing $\hat{\tilde{N}}(\kappa, s)$ as a sum of terms, each a rational function in $s$ multiplied by the transforms of the boundary data or source term.

5. Inverse Laplace transform — each rational term $s^p / (s^\mu - \psi_\alpha^\theta(\kappa))$ is recognized as the Laplace transform of a Mittag-Leffler function via the standard formula (80), producing a solution in $(\kappa, y)$-space involving $E_{\mu, \beta}(\pm y^\mu \psi_\alpha^\theta(\kappa))$.

6. Inverse Fourier transform and Fox H-function reduction — the remaining $\kappa$-integrals are expressed in closed form as Fox H-functions using the Mellin-Barnes representation of the Mittag-Leffler function (85) and the Mellin-cosine transform formula (86). The resulting H-functions are then expanded for asymptotic analysis using the known large-argument formula (87) and the series expansion (84).

Information flows strictly forward through this pipeline: PDE → Laplace → Fourier → algebraic solve → inverse Laplace → inverse Fourier → special function reduction → asymptotic/series expansion. There is no feedback or iteration; the entire methodology is a direct, sequential application of classical integral transform techniques to a new class of operators.

3.3 Roadmap for the Deep Dive

• First, the fractional derivative operators themselves — the Riesz-Feller derivative in $x$ and the Hilfer-composite derivative in $y$ — because understanding their Fourier and Laplace transform properties is prerequisite to everything that follows. If you don't know why $\mathcal{F}[xD_\theta^\alpha f](\kappa) = -\psi_\alpha^\theta(\kappa) \hat{f}(\kappa)$ and $\mathcal{L}[yD_{0+}^{\mu,\nu} f](s) = s^\mu \tilde{f}(s) - \text{(initial terms)}$, the solution derivation will be opaque.

• Second, the two preparatory lemmas (Lemma 1 and Lemma 2) that encode the inverse Laplace transforms of the algebraic forms that appear in the transform-domain solution — these are the bridge between the algebraic $(\kappa, s)$ solution and the Mittag-Leffler / Prabhakar integral representations in $(\kappa, y)$-space.

• Third, the full solution of the fractional Poisson equation (Theorem 1), walking through the sequential application of Laplace and Fourier transforms, the algebraic solve, and the inversion steps, because this is the master solution from which all special cases (Laplace, Helmholtz, wave equation) descend.

• Fourth, the critical sign distinction between the fractional Riesz-Feller derivative (Theorem 1) and the quantum fractional Riesz-Feller derivative (Remark 1), since this sign propagates through every Mittag-Leffler argument and determines whether the solutions represent the steady-state equations or map onto the time-dependent wave equation.

• Fifth, the Fox H-function reduction technique — specifically how the one-dimensional $\kappa$-integrals of Mittag-Leffler functions are evaluated in closed form using the Mellin-cosine transform formula — because the asymptotic and series expansions that constitute the paper's practical output all flow from these H-function representations.

• Sixth, the explicit asymptotic analysis (Remarks 6 and 7) showing how the Fox H-function is expanded for large and small values of $|x|/y^{\mu/2}$, because this is where the physical behavior of the solutions — stretched exponential decay versus power-law behavior — becomes manifest.

3.4 Detailed, Sentence-Based Technical Breakdown

This is primarily an analytical solution paper whose core idea is that fractional Laplace/Poisson/Helmholtz equations with mixed Riesz-Feller and Hilfer-composite derivatives are exactly solvable via sequential Fourier-Laplace transforms, and that the solutions naturally express themselves in terms of Mittag-Leffler functions and Fox H-functions whose asymptotic properties reveal the physical transport behavior.

---

The Fractional Derivative Operators: Riesz-Feller in $x$, Hilfer-Composite in $y$

The paper's entire solution strategy depends on understanding how the two fractional derivative operators behave under their respective integral transforms. The operators are fundamentally different in structure — one is a pseudo-differential operator defined through its Fourier symbol, the other is an integro-differential operator defined through fractional integrals — and this dictates which transform applies to which variable.

Riesz-Feller fractional derivative of order $\alpha$ and skewness $\theta$. This operator, denoted $xD_\theta^\alpha$, is defined not by an integral formula in $x$-space but by its action in Fourier space (equation 1):

$\mathcal{F}[xD_\theta^\alpha f(x)](\kappa) = -\psi_\alpha^\theta(\kappa) \mathcal{F}[f(x)](\kappa),$

where $\mathcal{F}[f](\kappa) = \hat{f}(\kappa) = \int_{-\infty}^\infty f(x) e^{i\kappa x} dx$ is the Fourier transform with convention $e^{+i\kappa x}$ in the forward direction, and the symbol $\psi_\alpha^\theta(\kappa)$ is defined by

$\psi_\alpha^\theta(\kappa) = |\kappa|^\alpha \exp\left(i \cdot \text{sign}(\kappa) \frac{\theta\pi}{2}\right), \quad 0 < \alpha \leq 2, \quad |\theta| \leq \min\{\alpha, 2-\alpha\}.$

The parameters are: $\alpha$ is the stability index controlling the heaviness of the tails of the underlying Lévy distribution (smaller $\alpha$ means heavier tails and more extreme jumps), and $\theta$ is the skewness parameter controlling the asymmetry of the distribution ($\theta > 0$ skews toward positive $x$, $\theta < 0$ toward negative $x$, $\theta = 0$ gives the symmetric Riesz derivative).

What this definition means operationally. In Fourier space, applying the Riesz-Feller derivative to a function simply multiplies its Fourier transform by $-\psi_\alpha^\theta(\kappa)$. This is the defining property of a pseudo-differential operator: it is local in frequency space but non-local in physical space. The factor $|\kappa|^\alpha$ is the amplitude scaling — it amplifies high frequencies (large $|\kappa|$) by a power $\alpha$, which in physical space corresponds to enhancing sharp variations. The complex exponential $\exp(i \cdot \text{sign}(\kappa) \theta\pi/2)$ introduces a phase shift that depends on the sign of $\kappa$, which in physical space produces an asymmetry in the operator's action — it treats left-moving and right-moving Fourier modes differently.

Why this form. The symbol $-\psi_\alpha^\theta(\kappa)$ is the logarithm of the characteristic function of a general Lévy strictly stable probability density with stability index $\alpha$ and asymmetry parameter $\theta$ (as the authors note, citing Mainardi, Pagnini, and Saxena, 2005). This means that the fundamental solution of the space-fractional diffusion equation $\partial_t u = xD_\theta^\alpha u$ is precisely the transition probability density of a Lévy flight — a stochastic process where a particle undergoes random jumps whose lengths are drawn from a heavy-tailed distribution with infinite variance (for $\alpha < 2$). Using this operator in the PDE therefore encodes the physical assumption of anomalous superdiffusive transport in the $x$-direction.

Quantum fractional Riesz-Feller derivative. The paper also defines a related operator $xD_\theta^{*,\alpha}$ (equation 5):

$\mathcal{F}[xD_\theta^{*,\alpha} f(x)](\kappa) = +\psi_\alpha^\theta(\kappa) \mathcal{F}[f(x)](\kappa).$

This is simply the negative of the classical Riesz-Feller derivative: $xD_\theta^{*,\alpha} = -xD_\theta^\alpha$. Its Fourier symbol is $+\psi_\alpha^\theta(\kappa)$ rather than $-\psi_\alpha^\theta(\kappa)$. The name "quantum" comes from its use by Luchko et al. (2013) in the fractional Schrödinger equation, where the sign convention aligns with the standard quantum-mechanical energy-momentum relation. The distinction matters because the sign propagates into the argument of every Mittag-Leffler function in the final solution — $E_{\mu,\beta}(+y^\mu \psi_\alpha^\theta(\kappa))$ for the quantum case versus $E_{\mu,\beta}(-y^\mu \psi_\alpha^\theta(\kappa))$ for the classical case — and this sign determines asymptotic stability and physical interpretability.

Hilfer-composite fractional derivative of order $\mu$ and type $\nu$. This operator acts in the $y$-variable and is defined through fractional integrals rather than through a Fourier symbol (equation 8 for order $0 < \mu \leq 1$, generalized in equation 13 for $n-1 < \mu \leq n$ with $n \in \mathbb{N}^+$):

$\left(yD_{a+}^{\mu,\nu} f\right)(y) = \left(yI_{a+}^{\nu(n-\mu)} \frac{d^n}{dy^n} \left(yI_{a+}^{(1-\nu)(n-\mu)} f\right)\right)(y),$

where $yI_{a+}^\mu$ is the Riemann-Liouville fractional integral of order $\mu$ (equation 6):

$\left(yI_{a+}^\mu f\right)(y) = \frac{1}{\Gamma(\mu)} \int_a^y \frac{f(y')}{(y - y')^{1-\mu}} dy', \quad y > a, \quad \Re(\mu) > 0.$

For $\mu = 0$, the integral operator reduces to the identity: $yI_{a+}^0 f = f$.

What this definition means operationally — a three-step procedure. To compute $yD_{a+}^{\mu,\nu} f(y)$ when $1 < \mu \leq 2$ (the case relevant to the paper's Laplace/Poisson/Helmholtz equations, so $n = 2$):

1. Apply a fractional integral of order $(1-\nu)(2-\mu)$: first compute $g(y) = (yI_{0+}^{(1-\nu)(2-\mu)} f)(y)$. This is a weighted average of $f$ over the interval $[0, y]$ with a kernel that decays as $(y-y')^{(1-\nu)(2-\mu)-1}$. The parameter $\nu$ controls the order of this first smoothing step — when $\nu = 1$, the order is zero and $g(y) = f(y)$ (no smoothing); when $\nu = 0$, the order is $2-\mu$ (maximal smoothing).

2. Differentiate twice: compute $h(y) = \frac{d^2}{dy^2} g(y)$, the ordinary second derivative of the smoothed function.

3. Apply another fractional integral of order $\nu(2-\mu)$: finally compute $(yI_{0+}^{\nu(2-\mu)} h)(y)$. This second smoothing step has order complementary to the first — the sum of the two fractional integral orders is $(1-\nu)(2-\mu) + \nu(2-\mu) = 2-\mu$, so the net operation is: smooth by a total of $2-\mu$, but with two integer derivatives embedded inside, giving a net differentiation of order $2 - (2-\mu) = \mu$.

Why this structure. The point of splitting the fractional integration into pre- and post-differentiation pieces with a parameter $\nu$ is that it interpolates between two standard definitions:

• When $\nu = 0$, the pre-integral has order $2-\mu$ and the post-integral has order $0$, so $yD_{0+}^{\mu,0} = \frac{d^2}{dy^2} (yI_{0+}^{2-\mu} f)$, which is the classical Riemann-Liouville fractional derivative (equation 9).
• When $\nu = 1$, the pre-integral has order $0$ and the post-integral has order $2-\mu$, so $yD_{0+}^{\mu,1} = (yI_{0+}^{2-\mu} \frac{d^2}{dy^2} f)$, which is the Caputo fractional derivative (equation 10).

The physical significance of this interpolation is that the two choices correspond to different types of initial or boundary conditions. The Riemann-Liouville derivative ($\nu = 0$) naturally pairs with initial conditions specified as fractional integrals evaluated at the origin. The Caputo derivative ($\nu = 1$) pairs with initial conditions specified as ordinary integer-order derivatives. The Hilfer derivative with $0 < \nu < 1$ allows a mixture: the initial conditions involve fractional integrals of intermediate order, which can be physically appropriate for systems where the memory kernel has a specific structure that is neither purely R-L nor purely Caputo.

Laplace transform of the Hilfer derivative. The entire solution method rests on knowing how $yD_{0+}^{\mu,\nu}$ transforms under the Laplace transform in $y$. For $1 < \mu \leq 2$ (so $n = 2$ in equation 13), Tomovski's formula (14) gives:

$\mathcal{L}[yD_{0+}^{\mu,\nu} f(y)](s) = s^\mu \tilde{f}(s) - s^{1-\nu(2-\mu)} \left(yI_{0+}^{(1-\nu)(2-\mu)} f\right)(0+) - s^{-\nu(2-\mu)} \left(\frac{d}{dy} yI_{0+}^{(1-\nu)(2-\mu)} f\right)(0+),$

where $\tilde{f}(s) = \mathcal{L}[f](s) = \int_0^\infty e^{-sy} f(y) dy$.

What this formula says in operational terms. The Laplace transform converts the complicated integro-differential Hilfer operator into:

1. A multiplication of the transform $\tilde{f}(s)$ by $s^\mu$ — the fractional-order analog of the rule $\mathcal{L}[f''](s) = s^2 \tilde{f}(s) - s f(0) - f'(0)$ for ordinary derivatives.
2. Two boundary terms evaluated at $y = 0^+$, involving:
   • The fractional integral of $f$ of order $(1-\nu)(2-\mu)$ evaluated at the origin: $(yI_{0+}^{(1-\nu)(2-\mu)} f)(0+)$.
   • The ordinary first derivative of that fractional integral evaluated at the origin: $\frac{d}{dy}(yI_{0+}^{(1-\nu)(2-\mu)} f)(0+)$.

These boundary terms are the generalization of $f(0)$ and $f'(0)$ in the ordinary derivative case. When $\nu = 1$ (Caputo), the first boundary term reduces to $-s^{1-0} f(0+) = -s f(0+)$ and the second to $-s^0 f'(0+) = -f'(0+)$, recovering the standard Caputo initial conditions. When $\nu = 0$ (Riemann-Liouville), the powers of $s$ become $-s^{1-(2-\mu)} = -s^{\mu-1}$ and $-s^{-(2-\mu)} = -s^{\mu-2}$, and the boundary terms involve fractional integrals of orders $2-\mu$ and their derivatives — the standard R-L initial conditions.

Why this transform property is the linchpin. The algebraic simplicity of $s^\mu \tilde{f}(s)$ in the transform domain is what makes the entire solution method work. After applying the Laplace transform in $y$, the $y$-derivative in the PDE becomes multiplication by $s^\mu$. After applying the Fourier transform in $x$, the $x$-derivative becomes multiplication by $-\psi_\alpha^\theta(\kappa)$. In the doubly-transformed $(\kappa, s)$ domain, the entire fractional PDE reduces to:

$(s^\mu - \psi_\alpha^\theta(\kappa)) \hat{\tilde{N}}(\kappa, s) = \text{boundary terms} + \text{source term},$

which is a purely algebraic equation. Solving for $\hat{\tilde{N}}$ requires only division by $s^\mu - \psi_\alpha^\theta(\kappa)$, and the inverse transforms are handled by recognizing the resulting rational forms as Laplace transforms of Mittag-Leffler functions.

---

Lemma 1: The Inverse Laplace Transform of the Core Algebraic Form

After solving the algebraic equation in $(\kappa, s)$-space, every term in the solution has the form

$\frac{s^{p}}{s^\mu \pm \hat{r}(\kappa)}$

multiplied by the Fourier transform of some boundary data or source term, where $p$ is a power determined by which boundary condition the term originates from and $\hat{r}(\kappa)$ is the Riesz-Feller symbol (or symbol plus $k^2$ for the Helmholtz case). Lemma 1 gives the inverse Laplace transform of this generic form.

Lemma 1. Let $1 < \mu \leq 2$, $0 \leq \nu \leq 1$, $\varsigma \geq 0$, and $\hat{r}(\kappa)$ be a given function. Then

$\mathcal{L}^{-1}\left[\frac{s^{\varsigma - \nu(2-\mu)}}{s^\mu \pm \hat{r}(\kappa)}\right](y) = y^{1 - (1-\nu)(2-\mu) - \varsigma} E_{\mu, 2 - (1-\nu)(2-\mu) - \varsigma}\left(\mp \hat{r}(\kappa) y^\mu\right),$

where $E_{\alpha,\beta}(z)$ is the two-parameter Mittag-Leffler function defined by

$E_{\alpha,\beta}(z) = \sum_{k=0}^\infty \frac{z^k}{\Gamma(\alpha k + \beta)}, \quad z, \beta \in \mathbb{C}, \quad \Re(\alpha) > 0.$

What this formula computes. Given the Laplace-domain expression $s^{\xi} / (s^\mu \pm \hat{r}(\kappa))$ where $\xi = \varsigma - \nu(2-\mu)$ is a specific exponent, the inverse Laplace transform in $y$ produces a power of $y$ multiplied by a Mittag-Leffler function whose argument is $\mp \hat{r}(\kappa) y^\mu$ and whose second parameter $\beta = 2 - (1-\nu)(2-\mu) - \varsigma$ is computed from the exponent $\xi$ and the fractional parameters $\mu$ and $\nu$.

Why this specific form of Mittag-Leffler function. The proof invokes the standard Laplace transform formula for the Mittag-Leffler function (equation 80):

$\mathcal{L}[t^{\beta-1} E_{\alpha,\beta}(\pm a t^\alpha)](s) = \frac{s^{\alpha-\beta}}{s^\alpha \mp a}, \quad \Re(s) > |a|^{1/\alpha}.$

To match this to the given $s^{\xi} / (s^\mu \pm \hat{r})$, one sets $\alpha = \mu$, $a = \hat{r}(\kappa)$, and requires $\mu - \beta = \xi$, which implies $\beta = \mu - \xi = \mu - (\varsigma - \nu(2-\mu))$. Simplifying: $\beta = \mu - \varsigma + \nu(2-\mu) = \mu - \varsigma + 2\nu - \nu\mu = 2\nu + \mu(1-\nu) - \varsigma + 1 - 1 = 2 - (1-\nu)(2-\mu) - \varsigma$. The power of $y$ in the inverse transform is $\beta - 1$ from the $t^{\beta-1}$ prefactor in the forward transform, which equals $1 - (1-\nu)(2-\mu) - \varsigma$.

The crucial consequence of this lemma is that it reduces the inverse Laplace transform step to a simple parameter-matching exercise: once the exponent $\varsigma$ is read off from the $s$-power in the algebraic solution, the corresponding Mittag-Leffler parameters are determined by substitution into the formula. This is what makes the solution method systematic — no case-by-case integral evaluation is needed.

The role of the sign $\mp$. The $\mp$ in the Mittag-Leffler argument corresponds to the $\pm$ in the denominator $s^\mu \pm \hat{r}(\kappa)$. A plus sign in the denominator ($s^\mu + \hat{r}$) produces a minus sign in the Mittag-Leffler argument ($-\hat{r} y^\mu$), and vice versa. This sign is not a minor detail — it determines whether the Mittag-Leffler function oscillates or grows/decays, and it is the mechanism by which the choice between classical and quantum Riesz-Feller derivatives (which differ by an overall sign) propagates into qualitatively different solution behavior.

---

Lemma 2: The Inverse Laplace Transform of the Source Term

When the PDE includes a source term $\Phi(x,y)$ (Poisson equation rather than Laplace equation), the algebraic solution in $(\kappa, s)$-space contains a term of the form

$\frac{1}{s^\mu \pm \hat{r}(\kappa)} \cdot \mathcal{L}[\hat{\Phi}(\kappa, \cdot)](s),$

i.e., the product of the same rational function and the Laplace transform of the Fourier-transformed source. Lemma 2 handles the inverse Laplace transform of this product.

Lemma 2. Let $1 < \mu \leq 2$, and let $\hat{r}(\kappa)$ and $\hat{\Phi}(\kappa, y)$ be given functions. Then

$\mathcal{L}^{-1}\left[\frac{1}{s^\mu \pm \hat{r}(\kappa)} \mathcal{L}[\hat{\Phi}(\kappa, \cdot)](s)\right](\kappa, y) = \left(E_{\mp \hat{r}(\kappa); 1}^{0+; \mu, \mu} \hat{\Phi}\right)(\kappa, y),$

where $E_{\omega; \gamma, \kappa}^{a+; \alpha, \beta}$ is the Srivastava-Tomovski integral operator (15) and the notation $E_{\mp \hat{r}(\kappa); 1}^{0+; \mu, \mu}$ denotes its specialization with $\omega = \mp \hat{r}(\kappa)$, $\gamma = 1$, lower limit $a = 0$, and parameters $\alpha = \mu$, $\beta = \mu$:

$\left(E_{\mp \hat{r}(\kappa); 1}^{0+; \mu, \mu} \hat{\Phi}\right)(\kappa, y) = \int_0^y (y - \xi)^{\mu-1} E_{\mu, \mu}^1 \left(\mp \hat{r}(\kappa) (y - \xi)^\mu\right) \hat{\Phi}(\kappa, \xi) d\xi.$

What this formula computes. Rather than giving a function of $y$ directly, this lemma expresses the inverse transform as an integral operator acting on $\hat{\Phi}$: the source term's contribution to the solution at $(\kappa, y)$ is the convolution of the source profile $\hat{\Phi}(\kappa, \xi)$ with a kernel $(y-\xi)^{\mu-1} E_{\mu,\mu}(\mp \hat{r}(\kappa)(y-\xi)^\mu)$ over the interval $\xi \in [0, y]$. The kernel involves the two-parameter Mittag-Leffler function $E_{\mu,\mu}$ with both parameters equal to $\mu$ (which is the special case of the three-parameter Prabhakar function $E_{\mu,\mu}^1$ with $\gamma=1$).

Why this form via convolution. The proof follows from the convolution theorem of the Laplace transform. Recognizing from equation (80) that

$\frac{1}{s^\mu \pm \hat{r}(\kappa)} = \mathcal{L}\left[y^{\mu-1} E_{\mu,\mu}(\mp \hat{r}(\kappa) y^\mu)\right](s),$

the product $\frac{1}{s^\mu \pm \hat{r}} \cdot \mathcal{L}[\hat{\Phi}](s)$ in the Laplace domain corresponds to the convolution in the $y$-domain of the inverse transforms:

$\int_0^y (y - \xi)^{\mu-1} E_{\mu,\mu}\left(\mp \hat{r}(\kappa) (y - \xi)^\mu\right) \hat{\Phi}(\kappa, \xi) d\xi.$

Why the Prabhakar integral operator notation is used. The integral operator $E_{\omega; \gamma}^{a+; \alpha, \beta}$ with kernel $(y-\xi)^{\beta-1} E_{\alpha,\beta}^\gamma(\omega (y-\xi)^\alpha)$ was introduced by Srivastava and Tomovski (2009) as a generalization of both the Riemann-Liouville fractional integral (recovered when $\omega = 0$) and the Prabhakar integral (recovered when $\gamma = 1$, which is the case used here). Writing the source term contribution as $(E_{\mp \hat{r}(\kappa); 1}^{0+; \mu, \mu} \hat{\Phi})(\kappa, y)$ is a notational convenience that emphasizes the structural similarity to fractional integration — the source term influences the solution through a type of weighted integration whose kernel is a Mittag-Leffler function rather than a power law. This matters physically because the Mittag-Leffler kernel has different asymptotic behavior from a pure power law (stretched exponential decay rather than algebraic decay), encoding the memory effects of the fractional derivative in how the source influences the field.

---

Theorem 1: Full Solution of the Fractional Poisson Equation

This is the central result of Section III. The equation being solved is:

$xD_\theta^\alpha N(x,y) + yD_{0+}^{\mu,\nu} N(x,y) = \Phi(x,y),$

where $x \in \mathbb{R}$, $y \in \mathbb{R}^+$, $1 < \alpha \leq 2$, $|\theta| \leq \min\{\alpha, 2-\alpha\}$, $1 < \mu \leq 2$, $0 \leq \nu \leq 1$.

The boundary conditions (24a) are specified in terms of fractional integrals at $y = 0^+$:

$\left(yI_{0+}^{(1-\nu)(2-\mu)} N\right)(x, 0^+) = f(x),$ $\left(\frac{d}{dy} yI_{0+}^{(1-\nu)(2-\mu)} N\right)(x, 0^+) = g(x),$

and the decay condition (24b):

$\lim_{x \to \pm \infty} N(x,y) = 0.$

The boundary conditions explained. Unlike an ordinary PDE where boundary data would be $N(x,0)$ and $\partial_y N(x,0)$, the Hilfer derivative's Laplace transform formula (14) shows that the natural boundary data for an equation with $yD_{0+}^{\mu,\nu}$ are the fractional integral of order $(1-\nu)(2-\mu)$ of $N$ and its ordinary $y$-derivative, both evaluated at $y = 0^+$. For the Caputo case $\nu = 1$, the pre-integral has order zero (identity operator), so $(yI_{0+}^0 N)(x,0+) = N(x,0+) = f(x)$ and its derivative is $g(x)$ — recovering the standard Dirichlet/Neumann boundary conditions. For the Riemann-Liouville case $\nu = 0$, the boundary data are fractional integrals of order $2-\mu$, which is a stronger condition at the boundary involving weighted averages of the solution near the origin.

The solution (equation 25):

$N(x,y) = \frac{y^{-(1-\nu)(2-\mu)}}{2\pi} \int_{-\infty}^\infty E_{\mu, 1-(1-\nu)(2-\mu)}\left(y^\mu \psi_\alpha^\theta(\kappa)\right) \hat{f}(\kappa) e^{-i\kappa x} d\kappa$ $+ \frac{y^{1-(1-\nu)(2-\mu)}}{2\pi} \int_{-\infty}^\infty E_{\mu, 2-(1-\nu)(2-\mu)}\left(y^\mu \psi_\alpha^\theta(\kappa)\right) \hat{g}(\kappa) e^{-i\kappa x} d\kappa$ $+ \frac{1}{2\pi} \int_{-\infty}^\infty \left(yE_{\psi_\alpha^\theta(\kappa); 1}^{0+; \mu, \mu} \hat{\Phi}\right)(\kappa, y) e^{-i\kappa x} d\kappa,$

where $\hat{f}(\kappa) = \mathcal{F}[f](\kappa)$, $\hat{g}(\kappa) = \mathcal{F}[g](\kappa)$, and $\hat{\Phi}(\kappa, y) = \mathcal{F}[\Phi(\cdot, y)](\kappa)$.

What each term represents.

• First term (the $f$-term): The contribution from the first boundary condition $f(x)$. It involves the Mittag-Leffler function $E_{\mu, 1-(1-\nu)(2-\mu)}$ with second parameter $\beta_1 = 1-(1-\nu)(2-\mu)$. The $y$-prefactor is $y^{-(1-\nu)(2-\mu)}$, which for $\nu = 1$ (Caputo) reduces to $y^0 = 1$ and for $\nu = 0$ (R-L) reduces to $y^{-(2-\mu)} = y^{\mu-2}$.

• Second term (the $g$-term): The contribution from the second boundary condition $g(x)$. It involves the Mittag-Leffler function $E_{\mu, 2-(1-\nu)(2-\mu)}$ whose second parameter is exactly one greater than that of the $f$-term: $\beta_2 = \beta_1 + 1$. The $y$-prefactor is $y^{1-(1-\nu)(2-\mu)}$, which is $y$ times the $f$-term's prefactor.

• Third term (the $\Phi$-term): The contribution from the source term, expressed as the Prabhakar integral operator acting on $\hat{\Phi}$ in the $y$-variable, followed by an inverse Fourier transform in $x$.

Why the solution has this structure — tracing the derivation step by step.

Step 1: Laplace transform in $y$. Applying the Laplace transform formula (14) to $yD_{0+}^{\mu,\nu} N$ with $n = 2$ (since $1 < \mu \leq 2$) yields:

$\mathcal{L}[yD_{0+}^{\mu,\nu} N(x,\cdot)](s) = s^\mu \tilde{N}(x,s) - s^{1-\nu(2-\mu)} f(x) - s^{-\nu(2-\mu)} g(x).$

The Laplace transform of the full PDE gives:

$xD_\theta^\alpha \tilde{N}(x,s) + s^\mu \tilde{N}(x,s) - s^{1-\nu(2-\mu)} f(x) - s^{-\nu(2-\mu)} g(x) = \tilde{\Phi}(x,s).$

Step 2: Fourier transform in $x$. Applying the Fourier transform with the symbol for the Riesz-Feller derivative, $\mathcal{F}[xD_\theta^\alpha \tilde{N}]( \kappa, s) = -\psi_\alpha^\theta(\kappa) \hat{\tilde{N}}(\kappa, s)$, transforms the PDE into a purely algebraic equation in the doubly-transformed variables:

$-\psi_\alpha^\theta(\kappa) \hat{\tilde{N}}(\kappa, s) + s^\mu \hat{\tilde{N}}(\kappa, s) - s^{1-\nu(2-\mu)} \hat{f}(\kappa) - s^{-\nu(2-\mu)} \hat{g}(\kappa) = \hat{\tilde{\Phi}}(\kappa, s).$

Step 3: Algebraic solve. Collecting the $\hat{\tilde{N}}$ terms:

$(s^\mu - \psi_\alpha^\theta(\kappa)) \hat{\tilde{N}}(\kappa, s) = s^{1-\nu(2-\mu)} \hat{f}(\kappa) + s^{-\nu(2-\mu)} \hat{g}(\kappa) + \hat{\tilde{\Phi}}(\kappa, s).$

Dividing through by $s^\mu - \psi_\alpha^\theta(\kappa)$:

$\hat{\tilde{N}}(\kappa, s) = \frac{s^{1-\nu(2-\mu)}}{s^\mu - \psi_\alpha^\theta(\kappa)} \hat{f}(\kappa) + \frac{s^{-\nu(2-\mu)}}{s^\mu - \psi_\alpha^\theta(\kappa)} \hat{g}(\kappa) + \frac{1}{s^\mu - \psi_\alpha^\theta(\kappa)} \hat{\tilde{\Phi}}(\kappa, s).$

This is the key algebraic result. Every term is a rational function in $s$ multiplied by a function of $\kappa$. The denominator $s^\mu - \psi_\alpha^\theta(\kappa)$ is always the same; only the $s$-power in the numerator differs.

Step 4: Inverse Laplace transform using Lemmas 1 and 2.

• For the $f$-term: $\frac{s^{1-\nu(2-\mu)}}{s^\mu - \psi_\alpha^\theta(\kappa)}$ matches Lemma 1 with $\hat{r}(\kappa) = \psi_\alpha^\theta(\kappa)$, $\pm$ sign taken as minus, and $\varsigma = 1$ (since $s^{\varsigma - \nu(2-\mu)} = s^{1-\nu(2-\mu)}$ implies $\varsigma = 1$). The inverse transform is $y^{1-(1-\nu)(2-\mu)-1} E_{\mu, 2-(1-\nu)(2-\mu)-1}(\psi_\alpha^\theta y^\mu) = y^{-(1-\nu)(2-\mu)} E_{\mu, 1-(1-\nu)(2-\mu)}(\psi_\alpha^\theta y^\mu)$.

• For the $g$-term: $\frac{s^{-\nu(2-\mu)}}{s^\mu - \psi_\alpha^\theta(\kappa)}$ matches Lemma 1 with $\varsigma = 0$. The inverse transform is $y^{1-(1-\nu)(2-\mu)-0} E_{\mu, 2-(1-\nu)(2-\mu)}(\psi_\alpha^\theta y^\mu) = y^{1-(1-\nu)(2-\mu)} E_{\mu, 2-(1-\nu)(2-\mu)}(\psi_\alpha^\theta y^\mu)$.

• For the $\Phi$-term: $\frac{1}{s^\mu - \psi_\alpha^\theta} \hat{\tilde{\Phi}}$ matches Lemma 2, yielding $(E_{\psi_\alpha^\theta; 1}^{0+; \mu, \mu} \hat{\Phi})(\kappa, y)$.

Step 5: Inverse Fourier transform. The result is multiplied by $e^{-i\kappa x} / (2\pi)$ and integrated over $\kappa \in (-\infty, \infty)$ to return to $(x,y)$-space, giving equation (25).

Why the Fourier inversion convention matters. The solution uses the convention $\mathcal{F}^{-1}[\hat{h}](x) = \frac{1}{2\pi} \int_{-\infty}^\infty \hat{h}(\kappa) e^{-i\kappa x} d\kappa$ (from equation 2), which is the standard convention paired with the forward transform $\mathcal{F}[f](\kappa) = \int_{-\infty}^\infty f(x) e^{+i\kappa x} dx$. The factor $1/(2\pi)$ appears in every solution term and is essential for dimensional correctness — omitting it would alter the normalization of the fundamental solution.

---

Remark 1: The Quantum Riesz-Feller Case

The parallel result for the quantum fractional Riesz-Feller derivative (equation 28):

$xD_\theta^{*,\alpha} N(x,y) + yD_{0+}^{\mu,\nu} N(x,y) = \Phi(x,y),$

is obtained by replacing the Fourier symbol $-\psi_\alpha^\theta(\kappa)$ with $+\psi_\alpha^\theta(\kappa)$. Tracing through the derivation:

• In the algebraic step, the equation becomes $+\psi_\alpha^\theta(\kappa) \hat{\tilde{N}} + s^\mu \hat{\tilde{N}} - \text{boundary terms} = \hat{\tilde{\Phi}}$, i.e., $(s^\mu + \psi_\alpha^\theta(\kappa)) \hat{\tilde{N}} = \dots$.
• The denominator changes from $s^\mu - \psi_\alpha^\theta$ to $s^\mu + \psi_\alpha^\theta$.
• In Lemma 1, the $\pm$ sign is now taken as plus, producing $E_{\mu,\beta}(-\psi_\alpha^\theta y^\mu)$ instead of $E_{\mu,\beta}(+\psi_\alpha^\theta y^\mu)$.
• The final solution (29) differs from (25) only in the sign of the argument inside every Mittag-Leffler function and the sign inside the Prabhakar integral operator's kernel.

Why this sign reversal is physically significant. The classical Riesz-Feller derivative has symbol $-\psi_\alpha^\theta(\kappa)$, which makes the fundamental solution of $xD_\theta^\alpha u = 0$ a probability density (positive and normalizable) — the negative sign ensures that the equation describes diffusion-like spreading rather than anti-diffusion. The quantum version has symbol $+\psi_\alpha^\theta(\kappa)$, which corresponds to a wave-like operator (the fractional Schrödinger equation uses this sign convention). In the context of this paper, the quantum Riesz-Feller solutions are the ones that map directly onto the time-dependent fractional wave equation (Remark 8), because the wave equation $tD_{0+}^{\mu,\nu} N = xD_\theta^\alpha N + \Phi$ (with classical Riesz-Feller on the right) rearranges to $xD_\theta^\alpha N - tD_{0+}^{\mu,\nu} N = -\Phi$, and the negative sign on the time derivative corresponds to the quantum Riesz-Feller in the steady-state formulation.

---

Corollary 1: Special Source Term $\Phi(x,y) = \delta(x) y^{-\beta} / \Gamma(1-\beta)$

When the source has the separable form of a Dirac delta in $x$ and a power law in $y$, the Prabhakar integral in the $\Phi$-term can be evaluated explicitly. The source is $\Phi(x,y) = \delta(x) \frac{y^{-\beta}}{\Gamma(1-\beta)}$, whose Fourier transform is $\hat{\Phi}(\kappa, y) = 1 \cdot \frac{y^{-\beta}}{\Gamma(1-\beta)}$ (since $\mathcal{F}[\delta(x)] = 1$).

The Prabhakar integral becomes:

$\left(E_{\pm \psi_\alpha^\theta; 1}^{0+; \mu, \mu} \hat{\Phi}\right)(\kappa, y) = \int_0^y (y-\xi)^{\mu-1} E_{\mu,\mu}(\pm \psi_\alpha^\theta(\kappa) (y-\xi)^\mu) \frac{\xi^{-\beta}}{\Gamma(1-\beta)} d\xi.$

Using the series expansion of the Mittag-Leffler function and interchanging sum and integral, this evaluates to (equation 30):

$y^{\mu-\beta} E_{\mu, \mu-\beta+1}(\pm y^\mu \psi_\alpha^\theta(\kappa)).$

Why this simplification works. The power-law form of the source is matched to the kernel of the fractional integral — the convolution of $y^{\mu-1} E_{\mu,\mu}(\omega y^\mu)$ with $y^{-\beta} / \Gamma(1-\beta)$ produces another Mittag-Leffler function with shifted parameters via term-by-term integration using the Beta function identity $\int_0^y (y-\xi)^{\mu-1} \xi^{-\beta} d\xi = y^{\mu-\beta} B(\mu, 1-\beta) = y^{\mu-\beta} \Gamma(\mu)\Gamma(1-\beta) / \Gamma(\mu-\beta+1)$, which combines with the series coefficients to yield the $E_{\mu, \mu-\beta+1}$ form.

The further reduction to a Fox H-function (the second expression in equation 30) uses two identities: first expressing the Mittag-Leffler function as an H-function via (85), and then evaluating the inverse Fourier transform integral $\int_{-\infty}^\infty H(\cdots) e^{-i\kappa x} d\kappa$ using the Mellin-cosine transform formula (86) with symmetry properties. The result is a single Fox H-function $H_{3,3}^{2,1}$ whose parameters encode all the fractional orders ($\alpha$, $\mu$, $\beta$, $\nu$, $\theta$).

---

Example 1: Point Source with Point Boundary Condition — Reduction to $H_{3,3}^{2,1}$

The most explicit special case is $\Phi(x,y) = \delta(x)\delta(y)$ (a unit point source at the origin), with $f(x) = \delta(x)$ and $g(x) = 0$, and symmetric Riesz derivative $\theta = 0$. This removes the skewness and simplifies the symbol to $\psi_\alpha^0(\kappa) = |\kappa|^\alpha$.

The solution (31) has two terms:

• The $f$-term becomes $\frac{y^{-(1-\nu)(2-\mu)}}{2\pi} \int_{-\infty}^\infty E_{\mu, 1-(1-\nu)(2-\mu)}(\pm y^\mu |\kappa|^\alpha) e^{-i\kappa x} d\kappa$, which reduces to $\frac{y^{-(1-\nu)(2-\mu)}}{|x|} H_{3,3}^{2,1}(\cdots)$.
• The $\Phi$-term (with $\beta = 0$ for the Dirac delta) becomes $\frac{y^{\mu-1}}{2\pi} \int_{-\infty}^\infty E_{\mu,\mu}(\pm y^\mu |\kappa|^\alpha) e^{-i\kappa x} d\kappa$, which reduces to $\frac{y^{\mu-1}}{|x|} H_{3,3}^{2,1}(\cdots)$ with different H-function parameters.

Why the $1/|x|$ factor appears. The Mellin-cosine transform formula (86) produces a factor of $\pi/x^\rho$ in the result. With $\rho = 1$ (since the inverse Fourier integral of an even function reduces to a cosine transform, $\int_{-\infty}^\infty \hat{F}(|\kappa|) e^{-i\kappa x} d\kappa = 2\int_0^\infty \hat{F}(\kappa) \cos(\kappa|x|) d\kappa$, and the Mellin-cosine formula gives $\pi/|x|$ for $\rho=1$). Combined with the $1/(2\pi)$ from the inverse Fourier prefactor, this yields $1/(2|x|)$, which matches equation (31).

The Fox H-function parameters for the $f$-term (equation 31, first H-function):

$H_{3,3}^{2,1}\left[\mp \frac{|x|^\alpha}{y^\mu} \middle| \begin{array}{c} (1,1), (1-(1-\nu)(2-\mu), \mu), (1, \frac{\alpha}{2}) \\ (1,\alpha), (1,1), (1,\frac{\alpha}{2}) \end{array}\right].$

The notation $H_{3,3}^{2,1}$ means $m=2$ (two poles from the lower row's gamma functions are used for the contour), $n=1$ (one gamma function from the upper row is in the numerator of the Mellin-Barnes integrand), $p=3$, $q=3$. The upper sign in $\mp$ corresponds to the classical Riesz-Feller derivative, the lower sign to the quantum version.

How to read the parameter arrays. The upper row lists three pairs $(a_j, A_j)$: the first pair $(1,1)$ and third pair $(1, \alpha/2)$ appear in the numerator via $\Gamma(1 - a_j + A_j s)$ for $j \leq n = 1$ (this is the $n=1$: only the first pair enters the numerator explicitly; the convention requires careful reading of formula 83). The second pair $(1-(1-\nu)(2-\mu), \mu)$ appears in the denominator via $\Gamma(a_j - A_j s)$ for $j > n$ (here $j = 2, 3$). The lower row lists $(b_j, B_j)$ pairs: first $(1,\alpha)$, second $(1,1)$, third $(1, \alpha/2)$. The first $m=2$ pairs enter the numerator via $\Gamma(b_j - B_j s)$, while the remaining $q-m=1$ pair enters the denominator via $\Gamma(1 - b_j + B_j s)$.

Why the H-function representation matters. The Fox H-function is not just compact notation — it is the gateway to asymptotic analysis. The large-argument expansion formula (87) applies specifically to H-functions of the form $H_{p,q}^{m,0}$ (with $m$ lower poles and zero upper poles in the numerator), which this $H_{3,3}^{2,1}$ can be transformed into using reduction formulas. The parameters $A_j$, $B_j$, $a_j$, $b_j$ directly determine the stretched exponential decay rate, the power-law prefactor, and the oscillatory behavior of the solution at large $|x|/y^{\mu/2}$.

---

Example 2: The $g$-Term Boundary Condition

The counterpart to Example 1 with $f(x) = 0$ and $g(x) = \delta(x)$ (equation 35) shifts all the Mittag-Leffler second parameters by $+1$ relative to the $f$-term case. The $g$-term's Mittag-Leffler function is $E_{\mu, 2-(1-\nu)(2-\mu)}$ rather than $E_{\mu, 1-(1-\nu)(2-\mu)}$, and the H-function's second upper parameter becomes $(2-(1-\nu)(2-\mu), \mu)$ instead of $(1-(1-\nu)(2-\mu), \mu)$. The $\Phi$-term remains identical because it does not depend on the boundary conditions.

Why this shift by 1 occurs. It is a direct consequence of Lemma 1: the $f$-term has $\varsigma = 1$ giving $\beta = 2 - (1-\nu)(2-\mu) - 1 = 1 - (1-\nu)(2-\mu)$, while the $g$-term has $\varsigma = 0$ giving $\beta = 2 - (1-\nu)(2-\mu)$. The difference of 1 in the second parameter reflects the fact that the $f$ boundary data enter through the $s^{1-\nu(2-\mu)}$ term in the Laplace transform, while the $g$ data enter through $s^{-\nu(2-\mu)}$, one power of $s$ lower.

---

The Reduction to $H_{2,2}^{2,0}$ for $\alpha = 2$ (Ordinary Laplacian in $x$)

When $\alpha = 2$, the Riesz-Feller derivative reduces to the ordinary second derivative (up to a sign). For the symmetric case $\theta = 0$, $\psi_2^0(\kappa) = |\kappa|^2 = \kappa^2$, the Fourier symbol of the standard Laplacian. In this limit, the Fox H-function $H_{3,3}^{2,1}$ simplifies to $H_{2,2}^{2,0}$ (equation 32), which further reduces to $H_{1,1}^{1,0}$ because two of the parameter pairs coincide and cancel via H-function reduction identities. The $H_{1,1}^{1,0}$ is equivalent to a Fox-Wright function and ultimately to a Wright function (equation 34).

Why this simplification happens. For $\alpha = 2$, the upper parameters include $(1, \alpha) = (1,2)$ and $(1, \alpha/2) = (1,1)$, and the lower parameters include the same $(1,\alpha/2) = (1,1)$. The pair $(1,1)$ appears in both numerator and denominator of the Mellin-Barnes integrand and cancels, reducing the effective $p$ and $q$ by 1 each. The resulting $H_{2,2}^{2,0}$ has structure that can be expressed as a simpler series — specifically the Wright function $\varphi(a,b;z)$ defined in equation (93).

---

Asymptotic Analysis: Large $|x|/y^{\mu/2}$ (Remark 6)

The asymptotic expansion for $|x|/y^{\mu/2} \gg 1$ is derived from the general large-argument formula for $H_{p,q}^{m,0}$ (equation 87) applied to the $H_{1,1}^{1,0}$ in equation (50). The derivation requires identifying the parameters:

• $p = 1$, $q = 1$, $m = 1$.
• $a_1 = 1 - (1-\nu)(2-\mu)$, $A_1 = \mu/2$.
• $b_1 = 1$, $B_1 = 1$.

From these, the auxiliary quantities in the asymptotic formula are:

• $m^* = \sum B_j - \sum A_j = 1 - \mu/2 = (2-\mu)/2$ (equation 89).
• $C = A_1^{A_1} B_1^{-B_1} = (\mu/2)^{\mu/2} \cdot 1^{-1} = (\mu/2)^{\mu/2}$ (equation 90).
• $\alpha = a_1 - b_1 + \frac{1}{2}(q - p + 1) = [1-(1-\nu)(2-\mu)] - 1 + \frac{1}{2}(1-1+1) = -(1-\nu)(2-\mu) + 1/2$ (equation 88).

The leading asymptotic behavior is then (equation 51):

$N(x,y) \simeq \frac{(\mu/2)^{\frac{1-2\nu+\frac{1}{2-\mu}}{2}}}{\sqrt{(2-\mu)\pi}} \frac{y^{-(1-\nu)(2-\mu)}}{|x|} \left(\frac{|x|}{y^{\mu/2}}\right)^{\frac{1+2(1-\nu)(2-\mu)}{2-\mu}} \exp\left[-\frac{2-\mu}{2}\left(\frac{\mu}{2}\right)^{\frac{\mu}{2-\mu}} \left(\frac{|x|}{y^{\mu/2}}\right)^{\frac{2}{2-\mu}}\right].$

What this expansion reveals physically. The dominant factor is the stretched exponential $\exp(-C |x|^{2/(2-\mu)} y^{-\mu/(2-\mu)})$. For $\mu = 2$ (ordinary diffusion), this becomes $\exp(-|x|^2 / (4y))$ — the familiar Gaussian. For $\mu < 2$, the decay in $|x|$ is slower than Gaussian but faster than any power law: it is a stretched exponential. The power-law prefactor $|x|^{(1+2(1-\nu)(2-\mu))/(2-\mu)-1}$ provides subleading corrections. The exponent's dependence on $\mu$ means that the fractional order in the $y$-direction directly controls how quickly the solution decays laterally — smaller $\mu$ (stronger memory effects) produces slower decay, consistent with the physical intuition that memory (subdiffusion) traps the solution near the origin.

Why the sign choice matters for asymptotic validity. This asymptotic formula holds for the quantum Riesz-Feller case (leading to $H_{1,1}^{1,0}$ with argument $+|x|/y^{\mu/2}$), which is appropriate for the steady-state equation. For the classical Riesz-Feller case, the argument sign flips, and the asymptotic behavior may involve oscillations or exponential growth rather than decay, reflecting the different physical nature of the two operators.

---

Series Representation: Small $|x|/y^{\mu/2}$ (Remark 7)

For $|x|/y^{\mu/2} \ll 1$, the Fox H-function is expanded using the residue formula (84), which expresses the H-function as an infinite sum of powers of its argument. The result for the $H_{1,1}^{1,0}$ in equation (50) is (equation 52):

$N(x,y) = \frac{y^{-(1-\nu)(2-\mu)-\frac{\mu}{2}}}{2} \sum_{j=0}^\infty \frac{(-1)^j}{j! \Gamma\left(1 - (1-\nu)(2-\mu) - \frac{\mu}{2}(j+1)\right)} \left(\frac{|x|}{y^{\mu/2}}\right)^j.$

Why this is a Wright function. Comparing with the definition of the Wright function (93):

$\varphi(a,b;z) = \sum_{n=0}^\infty \frac{1}{\Gamma(b + na)} \frac{z^n}{n!},$

the series matches with $a = -\mu/2$, $b = 1 - (1-\nu)(2-\mu) - \mu/2$, and $z = -|x|/y^{\mu/2}$. The factor $y^{-(1-\nu)(2-\mu)-\mu/2}/2$ times $\varphi(-\mu/2, 1-(1-\nu)(2-\mu)-\mu/2; -|x|/y^{\mu/2})$ gives the full small-argument expansion.

What this expansion reveals physically. Near the $y$-axis (small $|x|$), the solution behaves as a power series in $|x|/y^{\mu/2}$. The leading constant term ($j=0$) gives $N(0,y) \propto y^{-(1-\nu)(2-\mu)-\mu/2}$, showing algebraic decay along the $y$-axis. The $j=1$ term gives the leading $|x|$-correction. The alternating signs $(-1)^j$ indicate that for the quantum Riesz-Feller case, the solution has a maximum at $x=0$ and decreases with $|x|$ (the coefficients alternate if the gamma function in the denominator is positive). The dependence on the interpolation parameter $\nu$ enters through the shift $-(1-\nu)(2-\mu)$ in the gamma function argument — changing $\nu$ from 0 to 1 shifts the effective exponent of the leading term.

---

The Fractional Helmholtz Equation: Adding the $k^2$ Term (Theorem 2)

The Helmholtz equation differs from the Poisson equation by the inclusion of a term $+k^2 N(x,y)$ (equation 61):

$xD_\theta^\alpha N(x,y) + yD_{0+}^{\mu,\nu} N(x,y) + k^2 N(x,y) = \Phi(x,y).$

Why the $k^2$ term is mathematically significant. In the standard Helmholtz equation ($\alpha = 2$, $\mu = 2$, $\nu = 1$), the $+k^2$ term changes the fundamental solution from the $1/r$ Coulomb potential (Laplace) to the oscillatory $\cos(kr)/r$ or $\exp(ikr)/r$ (Helmholtz). In the fractional case, the $+k^2$ term shifts the Riesz-Feller symbol everywhere it appears:

• In the Laplace-Fourier domain, the algebraic equation becomes $(s^\mu + k^2 - \psi_\alpha^\theta(\kappa)) \hat{\tilde{N}} = \dots$, so the denominator changes from $s^\mu - \psi_\alpha^\theta$ to $s^\mu - (\psi_\alpha^\theta - k^2)$.
• Effectively, every occurrence of $\psi_\alpha^\theta(\kappa)$ in the Laplace and Poisson solutions is replaced by $\psi_\alpha^\theta(\kappa) - k^2$.
• For the quantum Riesz-Feller case (equation 66), the replacement is $+k^2$ added inside the symbol: $s^\mu + \psi_\alpha^\theta + k^2$ yields the denominator $s^\mu - (-(\psi_\alpha^\theta + k^2))$, so the effective argument becomes $-y^\mu(\psi_\alpha^\theta(\kappa) + k^2)$.

The solution (equation 63) therefore follows directly from Theorem 1 by replacing $\psi_\alpha^\theta(\kappa)$ with $\psi_\alpha^\theta(\kappa) - k^2$ everywhere in the Mittag-Leffler arguments and the Prabhakar kernel. No rederivation is needed. The $f$-term becomes:

$\frac{y^{-(1-\nu)(2-\mu)}}{2\pi} \int_{-\infty}^\infty E_{\mu, 1-(1-\nu)(2-\mu)}\left(y^\mu (\psi_\alpha^\theta(\kappa) - k^2)\right) \hat{f}(\kappa) e^{-i\kappa x} d\kappa,$

with analogous shifts in the $g$-term and $\Phi$-term.

Why the wave number $k$ does not change the solution methodology. The shift $\psi_\alpha^\theta \to \psi_\alpha^\theta - k^2$ is purely algebraic — it does not alter the structure of the inverse Laplace transform because $k^2$ is a constant with respect to $s$, so it simply shifts the pole location in the complex $s$-plane from $s = (\psi_\alpha^\theta)^{1/\mu}$ to $s = (\psi_\alpha^\theta - k^2)^{1/\mu}$. The inverse Fourier transform is more complicated because $\psi_\alpha^\theta(\kappa) - k^2$ does not have the simple scaling form $|\kappa|^\alpha e^{i\theta\pi/2}$, which prevents the clean Fox H-function reduction available for the Laplace/Poisson case. The authors do not provide explicit Fox H-function representations for the Helmholtz solutions — they leave them in integral form (equations 63, 67), recognizing that the $k^2$ shift breaks the homogeneity in $\kappa$ that the Mellin-cosine transform formula (86) requires.

---

Connection to the Time-Dependent Wave Equation (Remarks 8 and 10)

The paper's final major technical contribution is showing that the steady-state solutions map directly onto solutions of the space-time fractional wave equation. The wave equation (57) is:

$tD_{0+}^{\mu,\nu} N(x,t) = xD_\theta^\alpha N(x,t) + \Phi(x,t),$

where $xD_\theta^\alpha$ is the classical Riesz-Feller derivative with symbol $-\psi_\alpha^\theta(\kappa)$. Rearranging:

$xD_\theta^\alpha N(x,t) - tD_{0+}^{\mu,\nu} N(x,t) = -\Phi(x,t).$

Comparing with the steady-state quantum Riesz-Feller equation (28):

$xD_\theta^{*,\alpha} N(x,y) + yD_{0+}^{\mu,\nu} N(x,y) = \Phi(x,y),$

and using $xD_\theta^{*,\alpha} = -xD_\theta^\alpha$, the correspondence is:

• Replace $y$ with $t$ (the Hilfer derivative variable becomes time).
• Replace $xD_\theta^{*,\alpha}$ with $- xD_\theta^\alpha$ (i.e., use the classical Riesz-Feller symbol $-\psi_\alpha^\theta$).
• Change the sign of the source term: $-\Phi$ in the wave equation corresponds to $\Phi$ in the steady-state equation.

Making these substitutions in the steady-state solution (29) — which already has the quantum sign in the Mittag-Leffler arguments $E_{\mu,\beta}(-\psi_\alpha^\theta y^\mu)$ — and flipping the source sign yields equation (59). The negative sign in $E_{\mu,\beta}(-t^\mu \psi_\alpha^\theta(\kappa))$ is exactly what is needed for the classical Riesz-Feller wave equation.

Why this mapping works without rederivation. The structural identity between the steady-state equation with quantum Riesz-Feller derivative and the time-dependent wave equation with classical Riesz-Feller derivative is a consequence of the sign relationship $xD_\theta^{*,\alpha} = -xD_\theta^\alpha$. The Laplace transform in $y$ (for steady-state) and in $t$ (for time-dependent) produces the identical algebraic form $s^\mu \tilde{N} - \dots$ because the Hilfer derivative's Laplace transform formula (14) is the same regardless of whether the variable is spatial or temporal — it only depends on the operator's structure. The Fourier transform in $x$ produces $+\psi_\alpha^\theta$ for the quantum steady-state case and $-\psi_\alpha^\theta$ for the classical wave case, accounting for the sign difference. Thus, all the analytical machinery developed for Sections III and IV — the Mittag-Leffler representations, the Fox H-function reductions for special cases, the asymptotic expansions — applies directly to the time-dependent wave equation with the simple substitution $y \to t$ and sign adjustment. This is the paper's most elegant structural result: the steady-state solutions are not merely a special case but the mathematical foundation from which the full time-dependent behavior follows.

4. Key Insights and Innovations

Innovation 1: The Mixed-Derivative Steady-State Fractional Equation as a Unifying Analytical Framework

Prior to this paper, the fractional calculus literature treated time-fractional and space-fractional equations largely as separate problem classes. Time-fractional diffusion equations with Riemann-Liouville or Caputo derivatives in time and integer-order Laplacians in space were well-understood (Mainardi, 1996; Metzler and Klafter, 2000). Space-fractional diffusion equations with Riesz-Feller derivatives in space and first-order time derivatives had been solved by Mainardi, Pagnini, and Saxena (2005). And space-time fractional equations combining both — like Tomovski, Sandev, Metzler, and Dubbeldam (2012) — had been treated, but always with the fractional time derivative in the temporal variable and the Riesz-Feller derivative in the single spatial variable. What the field lacked was any treatment of two independent variables each carrying a fundamentally different type of fractional derivative, with no time coordinate at all — a genuine steady-state fractional PDE where the non-locality is spatial in both directions but of qualitatively different character in each.

This paper's core conceptual move is to pose and solve exactly this class: $xD_\theta^\alpha N(x,y) + yD_{0+}^{\mu,\nu} N(x,y) = \Phi(x,y)$, where the $x$-derivative is a pseudo-differential operator (Riesz-Feller, defined by its Fourier symbol) encoding Lévy-flight transport with stability index α and skewness θ, and the $y$-derivative is an integro-differential operator (Hilfer-composite, defined by nested fractional integrals) encoding memory-type non-locality with order μ and interpolating type ν. The physical interpretation is immediately richer than any single-derivative model: one spatial direction supports heavy-tailed jump processes (superdiffusion), while the orthogonal direction supports subdiffusive memory effects. The coupled solutions describe how these competing anomalous transport mechanisms interact.

This is not merely a cataloging exercise — the difference in operator type is what makes the problem both mathematically non-trivial and physically novel. The Riesz-Feller derivative lives naturally in Fourier space, where it becomes an algebraic multiplier $-\psi_\alpha^\theta(\kappa)$. The Hilfer derivative lives naturally in Laplace space, where it becomes $s^\mu$ plus initial-value terms. Neither transform alone can handle both operators. The solution requires a sequential application of both transforms — Laplace in $y$ first, then Fourier in $x$ — and the resulting doubly-transformed solution is a rational function in $s$ whose inverse Laplace transform involves Mittag-Leffler functions with the Riesz-Feller symbol in the argument. That the methodology works at all, producing closed-form expressions in terms of Fox H-functions whose asymptotic properties are explicitly computable, is not obvious a priori and constitutes the paper's primary analytical achievement.

The significance goes beyond the specific PDEs solved. By establishing that the Fourier-Laplace pipeline works for this mixed-derivative class, the paper opens a template for solving any steady-state fractional equation where the operators in different variables are of distinct types — pseudo-differential in one, integro-differential in another. The key insight is that the transform appropriate to each operator is applied to that variable, and the order of application is dictated by which boundary/initial conditions are specified in which variable. This is a methodological contribution that generalizes beyond the specific equations treated here.

Moreover, the paper's framework unifies previously disparate special cases. The solution for the general Hilfer type ν reduces to the Riemann-Liouville case when ν = 0 and to the Caputo case when ν = 1 (Corollary 3). The Riesz-Feller derivative reduces to the symmetric Riesz derivative when θ = 0 and to the classical Laplacian when α = 2, θ = 0. The Poisson equation reduces to the Laplace equation when Φ = 0 (Corollary 2). The Helmholtz equation reduces to the Poisson equation when k = 0. Rather than solving each special case separately — as the prior literature largely did — the paper derives one master solution (Theorem 1, equation 25) and obtains all special cases by parameter substitution. This is a fundamentally more elegant and systematic approach than the piecemeal treatments that preceded it, and it makes visible relationships between solutions (e.g., how the Helmholtz solution differs from the Poisson solution only by the shift ψ → ψ − k²) that separate derivations would obscure.

The evidence for this unification is structural: every equation in Sections III and IV — the Laplace equation (Corollary 2, equation 39), the Poisson equation (Theorem 1, equation 25), the Helmholtz equation (Theorem 2, equation 63), and the time-dependent wave equation (Remarks 8 and 10) — follows from the same transform pipeline with only algebraic sign changes or parameter substitutions distinguishing them.

---

Innovation 2: The Quantum/Classical Riesz-Feller Duality as a Bridge Between Steady-State and Time-Dependent Problems

A subtle but conceptually important contribution is the paper's identification and exploitation of the sign relationship between the classical Riesz-Feller derivative $xD_\theta^\alpha$ (Fourier symbol $-\psi_\alpha^\theta(\kappa)$) and the quantum Riesz-Feller derivative $xD_\theta^{*,\alpha}$ (Fourier symbol $+\psi_\alpha^\theta(\kappa)$). This is not a new operator definition — Luchko et al. (2013) introduced the quantum version for the fractional Schrödinger equation — but the paper's insight is that this sign distinction is the key that unlocks the mapping between steady-state spatial equations and time-dependent wave equations.

The mapping works as follows. The time-dependent fractional wave equation with classical Riesz-Feller space derivative (equation 57),

$tD_{0+}^{\mu,\nu} N(x,t) = xD_\theta^\alpha N(x,t) + \Phi(x,t),$

rearranges to $xD_\theta^\alpha N - tD_{0+}^{\mu,\nu} N = -\Phi$. Substituting $xD_\theta^\alpha = -xD_\theta^{*,\alpha}$ transforms this to $-xD_\theta^{*,\alpha} N - tD_{0+}^{\mu,\nu} N = -\Phi$, which is equivalent to $xD_\theta^{*,\alpha} N + tD_{0+}^{\mu,\nu} N = \Phi$ — the exact form of the steady-state equation (28) with the quantum Riesz-Feller derivative and $y$ replaced by $t$. Thus, every steady-state solution derived for the quantum Riesz-Feller case is simultaneously a time-dependent wave equation solution for the classical Riesz-Feller case, with the simple substitution $y \to t$ and appropriate relabeling of boundary conditions as initial conditions.

This duality is intellectually significant because it reframes the relationship between steady-state and time-dependent fractional equations. In the integer-order case, the Laplace equation $\nabla^2 u = 0$ and the wave equation $\partial_t^2 u = \nabla^2 u$ are clearly related — the latter's time-harmonic solutions satisfy the Helmholtz equation, which reduces to Laplace at zero frequency — but they are not obtained by a simple sign flip in the spatial operator. In the fractional case, the sign relationship between the quantum and classical Riesz-Feller derivatives creates a direct structural identity between the steady-state problem in two spatial variables and the initial-value problem in one spatial and one temporal variable. The physical interpretation shifts — $y$ is a spatial coordinate with boundary conditions in the steady-state case, while $t$ is a time coordinate with initial conditions in the wave case — but the mathematical solution is identical up to variable renaming.

Prior work treated these as separate problems requiring separate derivations. Tomovski and Sandev (2011, 2012) derived time-fractional wave equation solutions by applying Laplace transforms in time and Fourier transforms in space. Samuel and Thomas (2010) derived steady-state fractional Helmholtz solutions by applying the same transforms but treating both variables as spatial. Neither recognized that the solutions were structurally identical under the quantum/classical sign correspondence. This paper makes that correspondence explicit, which means that future work need not independently derive steady-state and time-dependent solutions for this class of operators — solving one gives the other for free. This is a conceptual economy that the prior literature missed.

The evidence is in Remarks 8 and 10: the wave equation solutions (59) and (77) are presented with no independent derivation, only the statement that they follow from the steady-state solutions by the operator correspondence. The fact that the authors can make this claim without additional proof is a testament to the structural clarity that the sign duality provides.

---

Innovation 3: Asymptotic Analysis Revealing the Physical Signatures of Mixed Fractional Transport

A substantial fraction of the paper's analytical effort — and, for applied users, perhaps its most practically valuable output — is the detailed asymptotic and series analysis of the Fox H-function solutions in two complementary regimes: the large-argument limit $|x|/y^{\mu/2} \gg 1$ (lateral distance large compared to the $y$-scale) and the small-argument limit $|x|/y^{\mu/2} \ll 1$ (near the $y$-axis). Prior work in the fractional calculus literature often stopped at expressing solutions in terms of Fox H-functions, treating the H-function as the final answer. This paper goes substantially further, extracting the explicit leading-order behavior in both limits and showing how the fractional parameters α, μ, ν, and θ control the physical character of the solution.

The large-argument asymptotics (Remark 6, equation 51) reveal stretched exponential decay:

$N(x,y) \propto \exp\left[-\frac{2-\mu}{2}\left(\frac{\mu}{2}\right)^{\frac{\mu}{2-\mu}} \left(\frac{|x|}{y^{\mu/2}}\right)^{\frac{2}{2-\mu}}\right],$

multiplied by a power-law prefactor whose exponent depends on the Hilfer interpolation parameter ν. This is a genuinely new physical prediction: the fractional order μ in the $y$-direction controls the shape of the lateral decay (Gaussian when μ = 2, progressively more stretched — slower decay at moderate distances, faster at very large distances — as μ decreases toward 1), while the interpolation parameter ν controls the amplitude through the power-law prefactor without changing the stretched exponential's functional form. The Riesz-Feller parameters α and θ do not appear in this asymptotic expression because the analysis is specialized to α = 2 (the ordinary Laplacian limit in $x$), but the methodology for including them is established.

The small-argument expansion (Remark 7, equation 52) expresses the solution as a Wright function series — essentially a generalized power series in $|x|/y^{\mu/2}$ — which provides the near-axis behavior:

$N(x,y) \simeq \frac{y^{-(1-\nu)(2-\mu)-\mu/2}}{2} \sum_{j=0}^\infty \frac{(-1)^j}{j! \Gamma\left(1 - (1-\nu)(2-\mu) - \frac{\mu}{2}(j+1)\right)} \left(\frac{|x|}{y^{\mu/2}}\right)^j.$

The leading term ($j=0$) gives the on-axis solution $N(0,y) \propto y^{-(1-\nu)(2-\mu)-\mu/2}$, showing that the Hilfer type ν directly controls the algebraic decay rate along the $y$-axis. For the Caputo case ν = 1, this exponent is $-\mu/2$; for the Riemann-Liouville case ν = 0, it is $-(2-\mu)-\mu/2 = -2 + \mu/2$ — a qualitatively different scaling that reflects how the boundary condition structure (fractional integral vs. ordinary derivative at $y=0^+$) propagates into the solution's far-field behavior.

What makes this contribution distinctive is not the asymptotic techniques themselves — the large-argument formula for Fox H-functions (equation 87) and the residue expansion (equation 84) are standard results from Mathai, Saxena, and Haubold (2010) — but rather the systematic application to extract physically interpretable behavior from an otherwise opaque special function representation. The Fox H-function $H_{3,3}^{2,1}$ with its six parameter pairs does not transparently reveal that the solution decays as a stretched exponential with exponent $2/(2-\mu)$. The asymptotic analysis makes this manifest, and in doing so, it provides testable predictions: if a physical system is governed by a fractional Laplace equation with parameters (μ, ν), the lateral decay of the steady-state field should follow the specific stretched exponential form in equation (51), and fitting observed decay profiles to this form would allow experimental determination of the fractional parameters.

The evidence is in the explicit asymptotic formulas: Remark 6 for the large-$|x|$ regime applied to the $H_{1,1}^{1,0}$ representation (equation 50), Remark 7 for the small-$|x|$ Wright function series, and Remark 2 (equation 33) for the combined $f$-term plus source-term asymptotics in the point-source case (Example 1). No prior work on fractional Laplace/Poisson/Helmholtz equations provided this level of asymptotic detail across both regimes.

---

Innovation 4: The Prabhakar Integral Operator as the Natural Source-Term Response Function

When the fractional Poisson or Helmholtz equation includes a source term Φ(x,y), the solution involves a convolution of the source with a kernel derived from the Mittag-Leffler function (Lemma 2, equation 20). The paper identifies this convolution as an instance of the Prabhakar integral operator $E_{\omega;\gamma}^{a+;\alpha,\beta}$ introduced by Srivastava and Tomovski (2009), specialized here to γ = 1, α = μ, β = μ, and with the frequency-dependent "weight" parameter ω = ∓ ψ_α^θ(κ). This identification is more than notational — it embeds the source-term response into a known class of generalized fractional integral operators, making its properties (composition rules, mapping properties between function spaces, relationship to the Riemann-Liouville integral) accessible from the existing literature.

The insight is that the response to a source in a fractional PDE with Hilfer-composite derivatives is itself a type of fractional integration — not with a power-law kernel as in the standard Riemann-Liouville integral, but with a Mittag-Leffler kernel $(y-\xi)^{\mu-1} E_{\mu,\mu}(\omega (y-\xi)^\mu)$ that interpolates between a power law at small arguments (where $E_{\mu,\mu}(z) \approx 1/\Gamma(\mu)$) and a stretched exponential at large arguments. This means the source's influence on the solution at a point $(x,y)$ is not simply a weighted average over earlier "times" with algebraic weights, as in ordinary fractional diffusion, but a more complex averaging where the effective memory kernel's shape adapts to the Riesz-Feller symbol ψ_α^θ(κ) — i.e., to the spatial Fourier mode being considered. High spatial frequencies (large |κ|) experience a different effective memory kernel than low frequencies, because the argument $\omega = \pm \psi_\alpha^\theta(\kappa)$ in the Mittag-Leffler function changes the kernel's shape.

Prior work on fractional equations with source terms (e.g., Sandev and Tomovski, 2010; Tomovski and Sandev, 2012) typically expressed the source contribution as a convolution integral without recognizing its structural relationship to the Prabhakar operator. By naming this operator explicitly, the paper connects the solution to a broader theory of generalized fractional calculus where integral operators with Mittag-Leffler kernels have been studied as objects in their own right (Kilbas, Saigo, and Saxena, 2004; Srivastava and Tomovski, 2009). This is an incremental refinement — it doesn't change the solution — but it is an intellectually clarifying one that situates the results within the existing taxonomy of fractional operators and enables future work to leverage known properties of the Prabhakar operator (semigroup properties, inversion formulas, mapping between function spaces) without rederiving them.

The evidence for this identification is Lemma 2 and its application in Theorem 1 (equation 25, third term) and Theorem 2 (equation 63, third term). The explicit evaluation for the special source $\Phi(x,y) = \delta(x) y^{-\beta}/\Gamma(1-\beta)$ in Corollary 1 demonstrates that the Prabhakar integral of a power-law source reduces to another Mittag-Leffler function, confirming that the operator acts as a fractional integral-like transformer on the source profile.

5. Experimental Analysis

Evaluation Methodology

This section must be addressed differently than in a typical machine learning paper. This paper contains no numerical experiments, no training runs, no held-out test set evaluations, and no comparison against baselines. It is a purely analytical work in the mathematical physics tradition. The "experiments" are the derivations themselves — the closed-form solutions, special function representations, asymptotic expansions, and series representations that are verified by reduction to known special cases rather than by empirical measurement. The paper's validity rests on the correctness of its transform calculations, not on statistical performance metrics.

Because this paper does not report empirical experiments, this section will document what the paper does provide in terms of verification, special-case reduction, and analytical consistency checks, while being explicit about what is absent. I will not fabricate datasets, metrics, or baselines that do not exist.

• Dataset. There is no dataset. The paper solves specific partial differential equations (fractional Laplace, Poisson, and Helmholtz equations) with specified boundary conditions. The "problem instances" are parameter sets — choices of α, μ, ν, θ, k, and boundary/source functions — and the "results" are the analytical solution forms for each parameter configuration.

• Base model(s). Not applicable. The paper does not use machine learning models or numerical solvers. The analytical framework is the Fourier-Laplace transform methodology, with the Riesz-Feller and Hilfer-composite fractional derivative operators as the mathematical objects under study.

• Metrics. The paper's results are assessed by mathematical correctness rather than empirical metrics. The primary mode of verification is reduction to known special cases: for specific parameter choices (ν = 0, ν = 1, θ = 0, α = 2, Φ = 0, k = 0), the newly derived fractional solutions are shown to match previously published results. This is the mathematical analog of a "sanity check" — if the general solution does not reproduce known special cases when parameters are specialized, the derivation contains an error.

• Baselines. The paper's "baselines" are the prior analytical solutions it subsumes. These include: Samuel and Thomas (2010) for the fractional Helmholtz equation with Riemann-Liouville y-derivative (ν = 0); Mainardi, Pagnini, and Saxena (2005) for the space-fractional diffusion equation with Riesz-Feller derivative; Tomovski, Sandev, Metzler, and Dubbeldam (2012) for the space-time fractional diffusion equation; and the classical integer-order Laplace/Poisson/Helmholtz solutions recovered when α = 2, μ = 2, ν = 1. The paper explicitly cites these correspondences in Corollaries 3 and 4, Remarks 8 and 10, and the special cases discussed in Examples 1-6.

• Generation budget / compute accounting. Not applicable. There is no computational budget to account for. The paper's results are symbolic expressions, not outputs of a computational process with FLOPs or generation counts.

• Cross-validation / statistical protocol. Not applicable. The verification is through mathematical identity checking, not statistical estimation. When the paper states that "many previously obtained results can be derived as special cases of those presented in this paper" (Abstract), the verification is that substituting specific parameter values into the general formulas yields expressions identical to those in the cited prior works — a deterministic equality check, not a statistical test.

---

Main Quantitative Results

This paper's "results" are analytical formulas, not numbers. I organize this section by the paper's logical groupings — types of equations solved — and document what the paper derives and how it verifies correctness through special-case reduction.

Fractional Laplace Equation Solutions

The fractional Laplace equation (Corollary 2, equation 39) with general boundary conditions f(x) and g(x) is the simplest nontrivial result:

$N(x,y) = \frac{y^{-(1-\nu)(2-\mu)}}{2\pi} \int_{-\infty}^\infty E_{\mu, 1-(1-\nu)(2-\mu)}\left(y^\mu \psi_\alpha^\theta(\kappa)\right) \hat{f}(\kappa) e^{-i\kappa x} d\kappa$ $+ \frac{y^{1-(1-\nu)(2-\mu)}}{2\pi} \int_{-\infty}^\infty E_{\mu, 2-(1-\nu)(2-\mu)}\left(y^\mu \psi_\alpha^\theta(\kappa)\right) \hat{g}(\kappa) e^{-i\kappa x} d\kappa.$

Headline finding: The solution decomposes into two terms corresponding to the two boundary data functions, each involving a Mittag-Leffler function whose second parameter encodes which boundary condition it originates from (β = 1 − (1−ν)(2−μ) for f, β = 2 − (1−ν)(2−μ) for g).

Special-case verification (Corollary 3): For the Caputo case ν = 1, the solution reduces to equation (42):

$N(x,y) = \frac{1}{2\pi} \int_{-\infty}^\infty E_\mu\left(y^\mu \psi_\alpha^\theta(\kappa)\right) \hat{f}(\kappa) e^{-i\kappa x} d\kappa + \frac{y}{2\pi} \int_{-\infty}^\infty E_{\mu,2}\left(y^\mu \psi_\alpha^\theta(\kappa)\right) \hat{g}(\kappa) e^{-i\kappa x} d\kappa,$

where the first term uses the one-parameter Mittag-Leffler function E_μ and the y-prefactors reduce to 1 and y. For the Riemann-Liouville case ν = 0, it reduces to equation (43):

$N(x,y) = \frac{y^{\mu-2}}{2\pi} \int_{-\infty}^\infty E_{\mu,\mu-1}\left(y^\mu \psi_\alpha^\theta(\kappa)\right) \hat{f}(\kappa) e^{-i\kappa x} d\kappa + \frac{y^{\mu-1}}{2\pi} \int_{-\infty}^\infty E_{\mu,\mu}\left(y^\mu \psi_\alpha^\theta(\kappa)\right) \hat{g}(\kappa) e^{-i\kappa x} d\kappa.$

These match the known forms: the Caputo case pairs with ordinary boundary values N(x,0) and ∂_y N(x,0) (since (1−ν)(2−μ) = 0 when ν = 1), while the R-L case pairs with fractional-integral boundary conditions. The transition between them is controlled continuously by ν.

Why this is a verification: The Caputo and Riemann-Liouville forms were known from prior work on time-fractional equations (Mainardi, 1996; Podlubny, 1999). The fact that the general Hilfer solution reduces to both known forms at the endpoints ν = 0 and ν = 1 confirms that the Laplace transform handling of the boundary terms — specifically the powers s^{1−ν(2−μ)} and s^{−ν(2−μ)} in equation (14) — is correctly propagated through the inverse transform.

Fractional Poisson Equation with Point Source

The Poisson equation adds a source term Φ(x,y). The most completely analyzed case is Example 1: Φ(x,y) = δ(x)δ(y) with f(x) = δ(x), g(x) = 0, θ = 0 (symmetric Riesz derivative). The solution (equation 31) is expressed as a sum of two Fox H-functions $H_{3,3}^{2,1}$:

$N(x,y) = \frac{y^{-(1-\nu)(2-\mu)}}{|x|} H_{3,3}^{2,1}\left[\mp \frac{|x|^\alpha}{y^\mu} \middle| \begin{array}{c} (1,1), (1-(1-\nu)(2-\mu), \mu), (1, \frac{\alpha}{2}) \\ (1,\alpha), (1,1), (1,\frac{\alpha}{2}) \end{array}\right]$ $+ \frac{y^{\mu-1}}{|x|} H_{3,3}^{2,1}\left[\mp \frac{|x|^\alpha}{y^\mu} \middle| \begin{array}{c} (1,1), (\mu, \mu), (1, \frac{\alpha}{2}) \\ (1,\alpha), (1,1), (1,\frac{\alpha}{2}) \end{array}\right].$

Headline finding: The solution is a sum of two Fox H-functions with identical structure but different second upper parameters: (1−(1−ν)(2−μ), μ) for the boundary term versus (μ, μ) for the source term. The argument ∓|x|^α/y^μ depends on the ratio of the two spatial variables raised to their respective fractional orders.

Special-case verification (α = 2, quantum Riesz-Feller, equation 32): When α = 2, the Fox H-functions reduce from $H_{3,3}^{2,1}$ to $H_{2,2}^{2,0}$ to $H_{1,1}^{1,0}$:

$N(x,y) = \frac{y^{-(1-\nu)(2-\mu)}}{2|x|} H_{1,1}^{1,0}\left[\frac{|x|}{y^{\mu/2}} \middle| \begin{array}{c} (1-(1-\nu)(2-\mu), \frac{\mu}{2}) \\ (1,1) \end{array}\right] + \frac{y^{\mu-1}}{2|x|} H_{1,1}^{1,0}\left[\frac{|x|}{y^{\mu/2}} \middle| \begin{array}{c} (\mu, \frac{\mu}{2}) \\ (1,1) \end{array}\right].$

This reduction happens because for α = 2, the parameter pairs (1,α) = (1,2) and (1,α/2) = (1,1) appear in the H-function, and the (1,1) pair cancels between numerator and denominator of the Mellin-Barnes integrand, reducing the effective number of parameters. The $H_{1,1}^{1,0}$ is equivalent to a Wright function (equation 34), which is a substantially simpler special function than the original $H_{3,3}^{2,1}$.

Verification against the integer-order case: For μ = 2, ν = 1, α = 2, the fractional Laplace equation should reduce to the classical Laplace equation ∂_x²N + ∂_y²N = 0 with N(x,0) = δ(x) and ∂y N(x,0) = 0. The classical solution is the Poisson kernel for the half-plane: N(x,y) = (1/π) · y/(x² + y²). The paper does not explicitly verify this reduction in the text, but the Fox H-function representation (equation 50) with μ = 2, ν = 1 gives N(x,y) = (1/(2|x|)) H{1,1}^{1,0}[|x|/y | (1,1); (1,1)], and evaluating this H-function (which reduces to elementary form when the parameters are integers) should recover y/(π(x² + y²)) — though the paper leaves this as implicit.

Asymptotic Analysis of the Point-Source Solution

Large |x|/y^{μ/2} asymptotics (Remark 6, equation 51): For the $f$-term contribution to the quantum Riesz-Feller case with α = 2, the leading asymptotic behavior for large lateral distance relative to the y-scale is:

$N(x,y) \simeq \frac{(\mu/2)^{\frac{1-2\nu+\frac{1}{2-\mu}}{2}}}{\sqrt{(2-\mu)\pi}} \frac{y^{-(1-\nu)(2-\mu)}}{|x|} \left(\frac{|x|}{y^{\mu/2}}\right)^{\frac{1+2(1-\nu)(2-\mu)}{2-\mu}} \exp\left[-\frac{2-\mu}{2}\left(\frac{\mu}{2}\right)^{\frac{\mu}{2-\mu}} \left(\frac{|x|}{y^{\mu/2}}\right)^{\frac{2}{2-\mu}}\right].$

Headline finding: The decay is a stretched exponential in $|x|/y^{\mu/2}$ with stretching exponent $2/(2-\mu)$. For μ = 2 (ordinary diffusion), this exponent diverges, recovering the Gaussian $\exp(-|x|^2/(4y))$. For 1 < μ < 2, the decay is slower than Gaussian but faster than any power law. The Hilfer interpolation parameter ν appears only in the power-law prefactor, not in the stretched exponential — meaning ν controls the amplitude scaling with y and |x| but not the functional form of the decay.

Verification: The μ → 2 limit. As μ → 2⁻, the exponent $2/(2-\mu) \to \infty$ and the stretched exponential should converge to a Gaussian. The prefactor involves $(2-\mu)$ in the denominator of several exponents, requiring careful limit analysis using the asymptotics of the gamma function via Stirling's formula. The paper does not carry out this limiting verification explicitly, but the general large-argument Fox H-function formula (87) from Braaksma (1964) is a rigorously established result that can be independently checked.

Combined asymptotics with source (Remark 2, equation 33): When both the $f$-term (boundary condition) and the source term (Φ = δ(x)δ(y)) contribute, the asymptotic behavior has two terms with different prefactors but the same stretched exponential factor $\exp[-\frac{2-\mu}{2}(\mu/2)^{\mu/(2-\mu)} (|x|/y^{\mu/2})^{2/(2-\mu)}]$, confirming that both the boundary-driven and source-driven components of the solution share the same spatial decay law at large distances.

Small |x|/y^{μ/2} asymptotics (Remark 7, equation 52):

$N(x,y) = \frac{y^{-(1-\nu)(2-\mu)-\frac{\mu}{2}}}{2} \sum_{j=0}^\infty \frac{(-1)^j}{j! \Gamma\left(1 - (1-\nu)(2-\mu) - \frac{\mu}{2}(j+1)\right)} \left(\frac{|x|}{y^{\mu/2}}\right)^j$

which is identified as a Wright function $\varphi$:

$N(x,y) = \frac{y^{-(1-\nu)(2-\mu)-\frac{\mu}{2}}}{2} \varphi\left(-\frac{\mu}{2}, 1-(1-\nu)(2-\mu)-\frac{\mu}{2}; -\frac{|x|}{y^{\mu/2}}\right).$

Headline finding: Near the y-axis, the solution is an alternating-sign power series in $|x|/y^{\mu/2}$, with the leading term (j = 0) giving the on-axis value $N(0,y) \propto y^{-(1-\nu)(2-\mu)-\mu/2}$. For the Caputo case ν = 1, this exponent is −μ/2; for the R-L case ν = 0, it is −(2−μ)−μ/2 = −2+μ/2 — demonstrating that the choice of fractional derivative type (ν) changes the algebraic decay rate along the y-axis even though the stretched exponential far-field behavior is unchanged.

Verification: This series is obtained from the residue expansion formula (84) for the Fox H-function, which is a standard result (Mathai, Saxena, and Haubold, 2010, Theorem 1.2). The correctness of the expansion can be checked by verifying that the residues at the poles s = −(1 + j) / (μ/2) for j = 0, 1, 2, ... produce the gamma function denominators as stated, and that the series converges for |x|/y^{μ/2} < ∞ (since the Wright function is entire for −μ/2 < 0).

---

Fractional Helmholtz Equation with Point Source

Example 5 solves the fractional Helmholtz equation (71) with Φ = δ(x)δ(y), f(x) = δ(x), g(x) = 0:

$N(x,y) = \frac{y^{-(1-\nu)(2-\mu)}}{2\pi} \int_{-\infty}^\infty E_{\mu, 1-(1-\nu)(2-\mu)}\left(y^\mu (\psi_\alpha^\theta(\kappa) - k^2)\right) e^{-i\kappa x} d\kappa$ $+ \frac{y^{\mu-1}}{2\pi} \int_{-\infty}^\infty E_{\mu,\mu}\left(y^\mu (\psi_\alpha^\theta(\kappa) - k^2)\right) e^{-i\kappa x} d\kappa.$

Headline finding: The solution differs from the Poisson case only by the replacement $\psi_\alpha^\theta(\kappa) \to \psi_\alpha^\theta(\kappa) - k^2$ inside the Mittag-Leffler arguments. This means the Helmholtz solution is structurally identical to the Poisson solution but with a shifted dispersion relation.

Key analytical difference from the Poisson case: The shift $\psi_\alpha^\theta(\kappa) - k^2$ breaks the homogeneity in κ that allowed the Poisson solutions to be expressed in closed-form Fox H-functions via the Mellin-cosine transform formula (86). For the Helmholtz case, the authors do not provide Fox H-function representations — the solutions remain as Fourier inversion integrals (equations 63, 67, 73, 74). This is not a limitation of the derivation but a genuine mathematical obstruction: the Mellin-cosine formula requires the integrand to be a function of $a \kappa^\delta$ with a specific homogeneity property, and the shift by $k^2$ destroys this homogeneity. The Helmholtz solutions are therefore analytically less complete than the Laplace/Poisson solutions — they provide the solution in principle but not in the computationally tractable special-function form needed for asymptotic analysis.

Verification against Samuel and Thomas (2010): Corollary 4 reduces the Helmholtz solution to the Riemann-Liouville case ν = 0 (equation 68) and the Caputo case ν = 1 (equation 69), matching the forms reported by Samuel and Thomas for the fractional Helmholtz equation with R-L y-derivative. The Samuel and Thomas solutions also involve Mittag-Leffler functions with the shifted argument and are left as Fourier integrals — the present paper's results are consistent with that prior work and generalize it to arbitrary ν.

Example 6 (the g-term boundary condition for Helmholtz): When f(x) = 0 and g(x) = δ(x), the solution (equation 74) has the same structure as Example 5 but with the first term using $E_{\mu, 2-(1-\nu)(2-\mu)}$ instead of $E_{\mu, 1-(1-\nu)(2-\mu)}$. The two boundary data functions produce complementary solutions with Mittag-Leffler second parameters offset by 1, exactly as in the Laplace/Poisson case — confirming that the $k^2$ shift does not affect the algebraic relationship between the f-term and g-term.

---

Time-Dependent Wave Equation Mapping

Remark 8: The space-time fractional wave equation (57) with classical Riesz-Feller space derivative has solution (59), obtained directly from the steady-state quantum Riesz-Feller solution (29) by replacing y with t and flipping the sign of the Mittag-Leffler argument.

Headline finding: The wave equation solution requires no independent derivation. All steady-state results — the integral representations, the special cases with point sources, the asymptotic expansions for α = 2, θ = 0 — transfer directly to the time-dependent case with the variable substitution y → t. The only analytical difference is that the initial conditions at t = 0⁺ for the wave equation (equation 58a) are interpreted as temporal initial data rather than spatial boundary data, but the mathematical form of the Laplace transform terms (s^{1−ν(2−μ)} f̂ and s^{−ν(2−μ)} ĝ) is identical because the Hilfer derivative's Laplace transform formula does not distinguish between spatial and temporal variables.

Verification against prior time-fractional wave equation solutions: The paper cites Tomovski and Sandev (2010, 2011) as having previously derived time-fractional wave equation solutions. The solution form (59) is consistent with those results when specialized to the same parameter choices (e.g., ν = 0 for R-L time derivative, ν = 1 for Caputo time derivative). However, the paper does not explicitly compare equation (59) term-by-term with the earlier Tomovski-Sandev solutions — the verification is by structural analogy rather than explicit reduction.

Remark 10 extends this to the wave equation with an additional −k²N term (equation 75):

$tD_{0+}^{\mu,\nu} N(x,t) = xD_\theta^\alpha N(x,t) - k^2 N(x,t) + \Phi(x,t),$

with solution (77) following from the Helmholtz steady-state solution (67) by the same y → t substitution and sign correspondence. This solution generalizes the standard damped wave equation (telegraph equation) to fractional orders in both space and time.

---

Special Case: The General Solution's Parameter Reduction Tree

The paper's most comprehensive "experimental result" is the systematic reduction of the general solution under parameter specialization. The following reduction tree is implicit across Sections III and IV but can be reconstructed from the explicit corollaries and examples:

1. General mixed-derivative equation (Theorem 1, equation 25; Theorem 2, equation 63): arbitrary α ∈ (1,2], μ ∈ (1,2], ν ∈ [0,1], θ with |θ| ≤ min(α, 2−α), arbitrary k, arbitrary Φ, f, g.
2. No source (Laplace/Helmholtz homogeneous): Φ = 0 → Corollary 2 (Laplace), equation 39; Theorem 2 with Φ = 0 (Helmholtz homogeneous).
3. Symmetric Riesz: θ = 0 → Examples 1-6 with $\psi_\alpha^0(\kappa) = |\kappa|^\alpha$.
4. Ordinary Laplacian in x: α = 2, θ = 0 → equations 32, 36, 48, 50, 54 (Fox H-function simplifies from $H_{3,3}^{2,1}$ to $H_{1,1}^{1,0}$).
5. Caputo in y: ν = 1 → Corollary 3, equations 42 (Laplace), 69 (Helmholtz).
6. Riemann-Liouville in y: ν = 0 → Corollary 3, equations 43 (Laplace), 68 (Helmholtz).
7. Ordinary derivatives in both variables: α = 2, μ = 2, ν = 1 → classical integer-order Laplace/Poisson/Helmholtz equations (implicit in the μ → 2 limits of the asymptotic formulas).
8. Wave equation mapping: Replace y → t throughout, interpret boundary conditions as initial conditions → Remarks 8 and 10.

Each reduction step must produce a known solution or a simpler special function. The paper verifies this explicitly for steps 2 (Φ = 0 reduces the three-term solution to two terms), 3 (θ = 0 simplifies the Riesz-Feller symbol), 4 (α = 2 triggers H-function cancellation), 5 and 6 (ν = 0,1 recover the two standard fractional derivative types), and 8 (the sign correspondence produces the wave equation without rederivation). Step 7 (full integer-order reduction) is implicit but follows from the asymptotic formulas' μ → 2 limits.

---

Ablation Studies and Robustness Checks

In the analytical context, "ablations" correspond to parameter specialization — removing one source of complexity (setting θ = 0, Φ = 0, k = 0, α = 2, or ν = 0,1) and verifying that the general solution reduces to the expected simpler form. The paper performs these checks systematically.

Boundary condition type (ν ∈ [0,1]): The solution's dependence on the Hilfer interpolation parameter ν is verified at the endpoints. When ν = 0 (Riemann-Liouville), the Laplace transform initial-value terms involve fractional integrals of order 2−μ and their derivatives, producing the y-prefactors $y^{\mu-2}$ and $y^{\mu-1}$ in equation (43). When ν = 1 (Caputo), the initial-value terms involve $N(x,0^+)$ and $\partial_y N(x,0^+)$ directly, producing the simpler prefactors 1 and y in equation (42). The general solution (39) interpolates continuously between these via the prefactors $y^{-(1-\nu)(2-\mu)}$ and $y^{1-(1-\nu)(2-\mu)}$, which reduce correctly at both endpoints. This is verified in Corollary 3.

Riesz-Feller skewness (θ): Setting θ = 0 eliminates the complex exponential in the symbol $\psi_\alpha^0(\kappa) = |\kappa|^\alpha$, reducing the Mittag-Leffler argument from $\pm y^\mu |\kappa|^\alpha \exp(i \cdot \text{sign}(\kappa) \theta\pi/2)$ to $\pm y^\mu |\kappa|^\alpha$. This simplification is crucial for the Mellin-cosine transform evaluations that produce Fox H-functions, because the Mellin-cosine formula (86) requires the integrand to be an even function of κ (or expressible as a cosine transform). For θ ≠ 0, the asymmetry parameter introduces a phase factor that breaks evenness in κ, preventing the clean H-function reduction. The paper does not provide closed-form H-function results for θ ≠ 0 — all Fox H-function examples (Examples 1-6) specify θ = 0. This is acknowledged implicitly by the structure of the presentation: general solutions (25, 39, 63) are integral representations valid for arbitrary θ, while the explicit special-function evaluations are restricted to θ = 0.

Source term presence (Φ): The difference between the Laplace equation (Φ = 0, Corollary 2) and the Poisson equation (Φ ≠ 0, Theorem 1) is the third term in the solution — the Prabhakar integral operator acting on $\hat{\Phi}$. Setting Φ = 0 eliminates this term, reducing the three-term Poisson solution to the two-term Laplace solution. The paper verifies this explicitly in Corollary 2 and implicitly in every example where the Φ-term and f/g-terms are computed separately (e.g., Example 1, equation 31, where the two H-functions correspond to the f-term and Φ-term respectively, and the g-term is absent because g = 0).

Wave number k (Helmholtz vs. Poisson): Setting k = 0 in the Helmholtz solution (63) replaces $\psi_\alpha^\theta(\kappa) - k^2$ with $\psi_\alpha^\theta(\kappa)$, recovering the Poisson solution (25). This is verified by direct comparison of equations (63) and (25). The authors do not explicitly state "for k = 0 we recover equation (25)", but the structural identity is obvious from the algebra: every occurrence of $\psi_\alpha^\theta - k^2$ becomes $\psi_\alpha^\theta$ when k = 0.

Quantum vs. classical Riesz-Feller (sign of ψ): Every solution is presented with "±" or "∓" notation indicating the two cases. The classical Riesz-Feller case (Theorem 1, equation 25) has Mittag-Leffler arguments with $+\psi_\alpha^\theta$ (upper sign), while the quantum case (Remark 1, equation 29) has $-\psi_\alpha^\theta$ (lower sign). The structural difference propagates through all derived quantities: the Laplace-domain denominator changes from $s^\mu - \psi_\alpha^\theta$ to $s^\mu + \psi_\alpha^\theta$, and the asymptotic behavior flips between decaying and potentially oscillatory. The paper verifies consistency by showing in Remarks 8 and 10 that the quantum steady-state case maps correctly onto the classical time-dependent wave equation — if either sign were wrong, the mapping between equations (28) and (57) would fail.

Riesz stability index (α = 2): The reduction from $H_{3,3}^{2,1}$ to $H_{2,2}^{2,0}$ to $H_{1,1}^{1,0}$ when α = 2 (Example 1, equation 32) is verified by the H-function parameter reduction rules cited in Mathai, Saxena, and Haubold (2010). The cancellation occurs because the pairs (1,1) and (1,1) appear in both the numerator and denominator of the Mellin-Barnes integrand θ(s), and the contour can be chosen so that their gamma function factors cancel. This is a standard H-function identity, not independently rederived in the paper, but it serves as a consistency check: if the cancellation failed for α = 2, the solution would not reduce to the known integer-order form, indicating an error in the parameter assignments.

Point source vs. extended source: Corollary 1 evaluates the Prabhakar integral for the specific source $\Phi(x,y) = \delta(x) y^{-\beta} / \Gamma(1-\beta)$, producing the closed-form result $y^{\mu-\beta} E_{\mu,\mu-\beta+1}(\pm y^\mu \psi_\alpha^\theta(\kappa))$ in equation (30). For β = 0 (Dirac delta in y as well), this reduces to $y^\mu E_{\mu,\mu+1}$, which matches the term in Example 1. For general β, the result is a single Mittag-Leffler function rather than an integral, showing that power-law sources in y interact simply with the Hilfer derivative's memory structure. The paper does not provide verification for arbitrary source functions — only the δ(x) × power-law separable form is evaluated explicitly.

Negative result: No Fox H-function form for the Helmholtz equation. The paper does not reduce the Helmholtz solutions (63, 67) to Fox H-functions, because the shift ψ → ψ − k² destroys the homogeneity in κ needed for the Mellin-cosine transform. This is not a flaw — it is an analytical boundary that the paper acknowledges by presenting the Helmholtz solutions only as inverse Fourier integrals. Future work would need numerical integration or asymptotic methods for large κ (where ψ dominates k²) to evaluate these integrals for specific parameter values.

Negative result: No θ ≠ 0 Fox H-function forms. All closed-form Fox H-function results assume θ = 0 (symmetric Riesz derivative). The asymmetric case θ ≠ 0 introduces a complex phase $\exp(i \cdot \text{sign}(\kappa) \theta\pi/2)$ in the Mittag-Leffler argument, which breaks the evenness of the integrand in κ. The inverse Fourier transform can still be expressed as a sum of cosine and sine transforms, but the Mellin-cosine formula (86) applies only to the cosine part. The paper does not pursue this case further — the asymmetric Riesz-Feller solutions remain as unevaluated Fourier integrals.

---

Critical Assessment

This paper is an analytical mathematics paper, not an empirical ML paper, and it must be assessed on its own terms: the correctness and completeness of its derivations, the verification of special-case reductions, and the physical interpretability of its asymptotic results. The standard ML evaluation framework — train/test splits, baselines, ablation studies, FLOPs-matched comparisons — does not apply. The assessment below evaluates whether the paper's analytical results support its claims about solving fractional steady-state equations.

Claim 1: "The objective of this paper is to derive analytical solutions of fractional order Laplace, Poisson and Helmholtz equations in two variables derived from the corresponding standard equations in two dimensions by replacing the integer order partial derivatives with fractional Riesz-Feller derivative and generalized Riemann-Liouville fractional derivative" (Abstract).

This claim is fully supported for the Laplace and Poisson equations, and partially supported for the Helmholtz equation. Theorem 1 and Corollary 2 provide explicit closed-form solutions for the fractional Poisson and Laplace equations, with the solutions expressed in terms of Mittag-Leffler functions, Fox H-functions (for special parameter choices), and Prabhakar integral operators (for source terms). The derivations are presented in complete detail: the Laplace transform step (equation 26), the Fourier transform step (equation 27), the algebraic solve, the inverse Laplace transform via Lemmas 1-2, and the inverse Fourier transform. Every step's domain of validity (Re(s) conditions, convergence of integrals, parameter ranges) is specified.

For the Helmholtz equation, Theorem 2 provides solutions as inverse Fourier integrals of Mittag-Leffler functions. However, the paper does not provide Fox H-function representations or asymptotic expansions for the Helmholtz case, as it does for the Laplace/Poisson case. The reason — that the shift ψ → ψ − k² breaks the homogeneity needed for the Mellin-cosine formula — is mathematically sound, but it means the Helmholtz solutions are less analytically complete. The paper acknowledges this implicitly by presenting equations (63, 67, 73, 74) as integral forms without further reduction. This is not an error, but it limits the practical utility of the Helmholtz results compared to the Laplace/Poisson results.

The additional claim that "results for fractional wave equation are presented as well" (Abstract) is supported — but only through the mapping argument in Remarks 8 and 10, not through independent derivation. The wave equation solutions (59, 77) are obtained by substituting y → t and adjusting signs in the steady-state solutions. No time-domain analysis (evolution of initial conditions, wavefront propagation, energy conservation) is performed. This is consistent with the paper's stated scope — the title emphasizes steady-state equations, with wave equations as a secondary application — but the abstract somewhat overstates the wave equation contribution by listing it alongside the primary results.

Claim 2: "Asymptotic behavior and series representation of solutions are analyzed in detail" (Abstract).

This claim is supported for the special case α = 2, θ = 0 (ordinary Laplacian in x with symmetric Riesz derivative). The large-|x|/y^{μ/2} asymptotics (Remarks 2 and 6, equations 33 and 51) provide explicit stretched exponential decay forms with power-law prefactors, derived from the Fox H-function large-argument formula (87). The small-|x|/y^{μ/2} series (Remarks 3 and 7, equations 34 and 52) provide Wright function representations.

However, the asymptotic analysis is restricted to the quantum Riesz-Feller case (lower sign in the H-function argument). The asymptotic formula (87) is for $H_{p,q}^{m,0}$ with positive real argument — for the classical Riesz-Feller case (upper sign, implying a negative argument ∓|x|^α/y^μ with the upper sign being minus), the asymptotic behavior would involve complex arguments and potentially oscillatory behavior, which is not analyzed. The paper does not explicitly discuss this restriction.

Furthermore, the asymptotics are derived only for the $H_{1,1}^{1,0}$ representation (α = 2, θ = 0). For general α < 2, the solution involves $H_{3,3}^{2,1}$, and the asymptotic analysis of this more complex H-function — which would reveal the interplay between the Lévy flight parameter α and the memory parameter μ in the spatial decay — is not performed. The abstract's claim of "analyzed in detail" is therefore accurate only for a specific parameter subspace (α = 2, θ = 0), not for the full parameter space of the problem.

Claim 3: "Many previously obtained results can be derived as special cases of those presented in this paper" (Abstract).

This claim is well-supported by explicit parameter reductions throughout Sections III and IV. The Riemann-Liouville case (ν = 0) and Caputo case (ν = 1) are recovered in Corollary 3 for the Laplace equation and Corollary 4 for the Helmholtz equation. The symmetric Riesz derivative (θ = 0) is used in all Fox H-function examples. The ordinary Laplacian (α = 2) is reduced in equations 32, 36, 48, 50, 54. The Laplace equation (Φ = 0) is recovered in Corollary 2. The paper explicitly cites Samuel and Thomas (2010) as a prior result that is subsumed (Remark 8: "Note that for µ = 2, ν = 1 solution (48)...", and Corollary 4 citing the Samuel-Thomas fractional Helmholtz equation as the ν = 0 case).

However, the paper does not provide an exhaustive mapping of every cited prior work to specific special cases. For example, the solutions of Mainardi, Pagnini, and Saxena (2005) for space-fractional diffusion are cited as motivation but are not explicitly recovered as limits of the steady-state solutions (the Mainardi solutions are time-dependent, so the mapping would require the wave equation correspondence of Remark 8, which is not spelled out in detail). The Tomovski, Sandev, Metzler, and Dubbeldam (2012) space-time fractional diffusion solutions are cited but not explicitly compared term-by-term with the wave equation solutions (59, 77). The reductions that are provided are convincing for the specific claims made, but the paper overpromises slightly by implying all cited prior results are explicitly recovered.

Weakness 1: No numerical evaluation or visualization. The paper provides asymptotic formulas but never evaluates them numerically for specific parameter values to illustrate the solution behavior. A single plot of N(x,y) versus x for fixed y and various μ would concretely demonstrate the transition from Gaussian (μ = 2) to stretched exponential (1 < μ < 2) decay, making the asymptotic analysis more accessible. Similarly, the Wright function series (equation 52) could be truncated and plotted to show the near-axis behavior. The absence of any numerical illustration is standard for theoretical papers in this subfield but limits the paper's accessibility to applied users who want to use these solutions.

Weakness 2: No discussion of solution regularity or physical admissibility. The solutions involve Fox H-functions and Mittag-Leffler functions, which can have singularities at the origin or grow/decay in ways that may or may not be physically admissible (e.g., negative values for a probability density, divergent moments). The paper does not discuss whether the solutions satisfy physical constraints like positivity (for diffusion problems) or finite energy (for wave problems). The asymptotic analysis addresses far-field behavior but not the behavior near y = 0⁺ (the boundary), where the fractional derivative's singular kernel may produce singularities in the solution. This is a standard concern in fractional PDE theory that the paper does not address.

Weakness 3: The Helmholtz solutions are incomplete relative to the Laplace/Poisson solutions. Theorem 2 (Helmholtz) provides integral representations but no Fox H-function forms and no asymptotic analysis. This is due to a genuine mathematical obstruction (the k² shift breaks homogeneity), but the paper does not discuss alternative approaches that could partially address this — for example, expanding the Helmholtz solution as a perturbation series in k² around the Poisson solution, or evaluating the Fourier integrals asymptotically using stationary phase or steepest descent for large |x|. The Helmholtz results are therefore substantially less developed than the Laplace/Poisson results, and the paper does not acknowledge this asymmetry.

Weakness 4: The general case (arbitrary α, θ, Φ, f, g) is left as unevaluated Fourier integrals. The closed-form Fox H-function results require θ = 0 and specific boundary/source choices (impulse functions or power laws). For general boundary data f(x) and g(x), the solutions remain as Fourier inversion integrals of the form $\int \hat{f}(\kappa) E_{\mu,\beta}(\cdots) e^{-i\kappa x} d\kappa$, which are not evaluated further. This is inherent in the methodology — the Fourier transform of general data cannot be simplified — but it means the paper's "analytical solutions" are, for the fully general case, solution representations rather than closed-form expressions. The paper is transparent about this (the integral forms are presented as the final answer for the general case), but readers expecting explicit formulas for arbitrary boundary data may be disappointed.

Weakness 5: No independent verification of the Fox H-function reductions. The Mellin-cosine transform formula (86) is cited from Mathai, Saxena, and Haubold (2010) without rederiving it for the specific parameter choices used. The parameter lists in the H-functions (e.g., the six pairs in $H_{3,3}^{2,1}$ in equation 31) are presented without showing the intermediate steps of the Mellin-Barnes integral evaluation. A reader attempting to verify these expressions would need to reconstruct the contour integral and residue calculations independently — the paper asserts the results without demonstrating the derivation. This is standard practice in special function theory (the identities are established in the cited monograph), but it means the verification is by authority rather than by explicit computation within the paper.

Overall assessment: The paper achieves its primary objective — deriving analytical solutions for fractional Laplace, Poisson, and Helmholtz equations with mixed Riesz-Feller and Hilfer-composite derivatives — for the Laplace and Poisson cases with symmetric Riesz derivative (θ = 0). The Fox H-function representations and asymptotic expansions for the α = 2, θ = 0 case are rigorous and provide physical insight into the stretched exponential decay. The Helmholtz solutions and the fully general case (θ ≠ 0, arbitrary boundary data) are less completely characterized, with solutions left as integral representations. The wave equation extension is structurally elegant but depends entirely on the steady-state derivations and is not independently developed. The paper's verification strategy — reduction to known special cases at parameter endpoints — is appropriate for an analytical work and is executed correctly for the cases checked, though the verification is not exhaustive for all cited prior work.

6. Limitations and Trade-offs

6.1 Closed-Form Fox H-Function Solutions Are Restricted to the Symmetric Riesz Derivative (θ = 0)

The assumption or constraint. All Fox H-function representations provided in the paper — including the fully reduced asymptotic and series expansions in Examples 1–6 and Remarks 6–7 — assume θ = 0, i.e., the symmetric Riesz fractional derivative without skewness. The general solutions for arbitrary skewness θ (equations 25, 39, 63, and their quantum counterparts) are left as unevaluated inverse Fourier integrals of the form $\int_{-\infty}^\infty E_{\mu,\beta}(\pm y^\mu |\kappa|^\alpha e^{i \cdot \text{sign}(\kappa) \theta\pi/2}) \hat{f}(\kappa) e^{-i\kappa x} d\kappa$. The paper never states this constraint as a limitation, but the structure of the presentation makes it clear: every closed-form result specifies θ = 0 as a condition, and the Mellin-cosine transform formula (equation 86) — the gateway from Mittag-Leffler integrals to Fox H-functions — requires the integrand to be an even function of κ, which holds only when the complex phase $\exp(i \cdot \text{sign}(\kappa) \theta\pi/2)$ is absent.

The consequence. For θ ≠ 0 — which corresponds to asymmetric Lévy flights with preferred directionality in the underlying stochastic process — the solution cannot be expressed in terms of standard Fox H-functions and cannot benefit from the asymptotic machinery (equations 87–91) that the paper deploys for the symmetric case. A practitioner modeling physical systems with directional bias (e.g., transport in anisotropic porous media, financial returns with skewness, asymmetric search strategies) cannot use the paper's closed-form results and must resort to numerical evaluation of the Fourier inversion integrals, which are oscillatory and involve a special function with a complex argument — a computationally non-trivial task the paper does not address. The physical insight the paper provides — stretched exponential decay, the role of μ in the decay exponent — is therefore restricted to the symmetric case; whether asymmetric Lévy flights in the x-direction qualitatively change the solution's spatial decay (e.g., introducing oscillatory or power-law tails in preferred directions) is left unanswered.

What evidence exists in the paper. Every example that produces a Fox H-function — Example 1 (equation 31, specifying θ = 0), Example 2 (equation 35), Example 3 (equation 48), Example 4 (equation 53) — explicitly sets θ = 0. The general integral representations (25, 29, 39, 41, 63, 67) contain the full symbol $\psi_\alpha^\theta(\kappa) = |\kappa|^\alpha e^{i \cdot \text{sign}(\kappa) \theta\pi/2}$ and are valid for arbitrary θ, but the paper never attempts to evaluate these for θ ≠ 0 beyond stating them. The Mellin-cosine formula (86) is cited for the symmetric case only, and its conditions include $x^\delta > 0$ and arguments that are real powers — the complex exponential in $\psi_\alpha^\theta$ would violate the formula's assumptions.

Mitigation status. The paper does not acknowledge this limitation or propose workarounds. Possible mitigations — decomposing the asymmetric solution into cosine and sine transforms, each potentially evaluable via Mellin transforms with complex parameters, or developing asymptotic expansions for the Fourier integrals using stationary phase methods for large $|x|/y^{\mu/\alpha}$ — are not discussed. The limitation is inherent to the methodology: the Fourier-Laplace transform pipeline handles θ ≠ 0 at the algebraic level, but the inverse Fourier step hits a special-function barrier that the paper does not attempt to cross.

---

6.2 The Helmholtz Equation Solutions Lack Closed-Form Special Function Representations and Asymptotic Analysis

The assumption or constraint. The fractional Helmholtz equation (61) and its quantum counterpart (66) introduce a wave number term $+k^2 N(x,y)$ that, in the Laplace-Fourier transformed domain, shifts the Riesz-Feller symbol from $\pm \psi_\alpha^\theta(\kappa)$ to $\pm \psi_\alpha^\theta(\kappa) - k^2$ (classical) or $-(\pm \psi_\alpha^\theta(\kappa) + k^2)$ (quantum). The resulting inverse Fourier integrals (equations 63, 67, 73, 74) contain Mittag-Leffler functions with arguments of the form $\pm y^\mu (\psi_\alpha^\theta(\kappa) \mp k^2)$ — an additive constant $k^2$ inside the function that destroys the power-law homogeneity in κ required for the Mellin-cosine transform formula (86) to apply. The paper does not reduce these Helmholtz solutions to Fox H-functions, nor does it provide asymptotic expansions or series representations for them.

The consequence. In the standard integer-order case (α = 2, μ = 2, ν = 1), the Helmholtz equation's fundamental solution transitions from the $1/r$ Coulomb form (Laplace/Poisson) to an oscillatory $\cos(kr)/r$ or Hankel function form (Helmholtz). The fractional generalization of this oscillatory behavior — how the Lévy flight index α and memory parameter μ interact with the wave number k to produce spatial oscillations, decay, or resonance — is not analytically characterized by this paper. A practitioner interested in fractional Helmholtz models (e.g., frequency-domain vibration analysis of viscoelastic membranes, electromagnetic wave propagation in complex media, or the fractional Schrödinger equation with a potential) receives only integral representations that require numerical quadrature for each evaluation. The paper's headline asymptotic results (stretched exponential decay, Wright function series) apply only to the Laplace/Poisson case (k = 0) and do not extend to the Helmholtz case.

What evidence exists in the paper. Theorem 2 (equations 63, 67) presents the Helmholtz solutions as inverse Fourier integrals only. Corollaries 4 and 5 reduce these to the Caputo (ν = 1) and Riemann-Liouville (ν = 0) cases — but still as integrals, not as H-functions. Examples 5 and 6 provide the point-source Helmholtz solutions (equations 73, 74), again as Fourier integrals. Nowhere in Section IV does a Fox H-function appear; contrast this with Section III, where Examples 1–4 produce explicit H-function representations for the Laplace/Poisson case. The paper never explicitly states that the Helmholtz solutions cannot be reduced to H-functions — the absence of such representations is the implicit evidence that the methodology hits a barrier.

Mitigation status. The paper does not discuss the Helmholtz limitation or suggest how to overcome it. Possible approaches — expanding the Helmholtz solution as a perturbation series in $k^2$ around the Poisson solution (which would produce Fox H-functions at each order), developing large-κ asymptotics where $\psi_\alpha^\theta(\kappa) \gg k^2$ to recover stretched exponential decay with a $k^2$-dependent correction, or using complex-analytic methods (steepest descent on the Fourier inversion contour) to extract oscillatory behavior — are not explored. The paper's silence on this point means the Helmholtz results are significantly less developed than the Laplace/Poisson results, and a reader might reasonably conclude that the Helmholtz equation was included for structural completeness rather than because the paper provides practically usable solutions for it.

---

6.3 No Numerical Evaluation or Visualization of the Analytical Solutions

The assumption or constraint. The paper is a purely analytical work containing no numerical evaluation of any result. The Fox H-function representations, the asymptotic formulas (equations 33, 51), and the Wright function series (equations 34, 52) are presented as symbolic expressions with no plots, tables, or computed values for specific parameter choices. The paper implicitly assumes that expressing a solution in terms of recognized special functions (Mittag-Leffler, Fox H, Wright) constitutes a complete answer, and that the asymptotic formulas are sufficient for understanding the solution's physical behavior.

The consequence. For a practitioner seeking to use these solutions — to fit experimental data, to validate a numerical fractional PDE solver, or to understand the qualitative behavior of a fractional Laplace/Poisson model — the paper's results are incomplete. The Fox H-function $H_{3,3}^{2,1}$ with six parameter pairs (equation 31) is not a function available in standard numerical libraries (unlike, say, Bessel functions or the ordinary Mittag-Leffler function). Evaluating it requires implementing the Mellin-Barnes contour integral (83) or the residue series (84) with careful attention to convergence and branch cuts — a non-trivial numerical task that the paper does not address. The asymptotic formulas (equation 51) provide leading-order behavior, but without numerical validation, it is unclear at what value of $|x|/y^{\mu/2}$ the asymptotic regime begins to be accurate, or how the transition from the small-argument Wright function series (equation 52) to the large-argument stretched exponential occurs. A practitioner cannot determine from the paper alone whether, say, for μ = 1.5, ν = 0.5, the asymptotic formula is within 5% of the true solution at $|x|/y^{\mu/2} = 5$ or only at $|x|/y^{\mu/2} = 50$.

What evidence exists in the paper. The entire paper is symbolic. There are no figures. The asymptotic formulas carry no error bounds or subleading correction terms beyond the leading-order expression. The series expansions (equations 34 and 52) specify no truncation error estimates or radii of convergence beyond the formal statement that the Wright function is entire. The paper does not report having implemented any of the special functions numerically or having verified the asymptotic formulas against direct numerical integration of the Fourier inversion integrals.

Mitigation status. No mitigation is attempted. The paper does not cite numerical libraries for Fox H-functions, does not provide guidance on contour parameter choices for the Mellin-Barnes integral, and does not discuss the practical feasibility of evaluating its solutions. This is consistent with the norms of the pure fractional calculus literature of this period — where expressing a solution in terms of H-functions was considered the endpoint of the analysis — but it represents a genuine limitation for applied users. Future work on numerical evaluation of Fox H-functions (some of which existed at the time, e.g., through Mathematica packages, and has expanded since) could address this gap, but the paper itself provides no bridge from the symbolic results to computable quantities.

---

6.4 Single Equation Family and Operator Class; No Exploration of Alternative Fractional Derivative Definitions

The assumption or constraint. The paper solves a specific class of fractional PDEs — those involving the Riesz-Feller derivative in $x$ and the Hilfer-composite derivative in $y$ — and makes no attempt to generalize to other fractional derivative definitions that were actively used in the literature. The fractional calculus community of this period had multiple competing operator definitions with different physical motivations: the Caputo-Fabrizio derivative (with an exponential kernel, avoiding singularities at the origin), the Atangana-Baleanu derivative (with a Mittag-Leffler kernel, combining non-locality with non-singularity), the Marchaud derivative (based on finite differences), the Grünwald-Letnikov derivative (discrete analog of R-L), and distributed-order derivatives (integrals over the order parameter). None of these are considered. The paper also solves only three closely related equation types — Laplace, Poisson, Helmholtz — which differ only by the presence or absence of a zeroth-order term (0, Φ, or $k^2 N$). No other steady-state fractional PDEs (e.g., with advection terms, nonlinear terms, or spatially varying coefficients) are treated.

The consequence. The paper's solutions are not transferable to models that use different fractional derivative definitions. If a physical system is better described by a Caputo-Fabrizio operator (exponential memory kernel) rather than a Hilfer operator (power-law kernel), the entire Fourier-Laplace pipeline would need to be rederived because the Laplace transform of the Caputo-Fabrizio derivative is fundamentally different from equation (14): $\mathcal{L}[{}^{CF}D_t^\alpha f](s) = \frac{s \tilde{f}(s) - f(0)}{s + \alpha(1-s)}$ rather than $s^\alpha \tilde{f}(s)$. The paper's methodology — solving the algebraic equation in $(κ, s)$-space and inverting using the standard Mittag-Leffler Laplace transform pair (80) — is specific to the power-law kernel structure of the R-L/Hilfer family. A practitioner working with alternative fractional derivatives cannot adapt these solutions without essentially re-solving the problem from scratch. Similarly, adding an advection term $v \cdot \nabla N$ (fractional or otherwise) would couple the Fourier modes and break the algebraic simplicity of the $(κ, s)$-space equation, making the transform method inapplicable.

What evidence exists in the paper. The paper's entire derivation rests on two specific transform properties: the Fourier symbol of the Riesz-Feller derivative (equation 1) and the Laplace transform of the Hilfer derivative (equation 14). The paper does not discuss alternative operators or justify why the Riesz-Feller/Hilfer combination is physically preferred over others. The introduction cites applications that motivate fractional derivatives in general (anomalous diffusion, Lévy flights, glass relaxation) but does not argue that these specific operator choices are uniquely appropriate for any particular physical system. The scope is limited to the three named equations, with no suggestion of extensions to other PDE types.

Mitigation status. The paper does not frame the narrow operator scope as a limitation; it presents the chosen operators as the natural generalization of the integer-order Laplace/Poisson/Helmholtz equations. The fact that the Riesz-Feller and Hilfer derivatives are well-established in the fractional calculus literature — the former for Lévy flight spatial transport (Metzler and Klafter, 2000, 2004), the latter for interpolated memory in time-fractional problems (Hilfer, 2000; Sandev et al., 2011) — provides a reasonable justification for studying this particular combination. The paper's contribution is to solve this specific class, not to survey all possible fractional operators. Nonetheless, a practitioner comparing fractional derivative definitions for a specific application would need to look elsewhere for guidance on which operator to use and whether the paper's solutions are relevant to their problem.

---

6.5 The General Solutions for Arbitrary Boundary Data and Source Functions Are Integral Representations, Not Explicit Formulas

The assumption or constraint. The paper's most general solutions — Theorem 1 (equation 25) for the Poisson equation, Corollary 2 (equation 39) for the Laplace equation, and Theorem 2 (equation 63) for the Helmholtz equation — are expressed as inverse Fourier transforms of Mittag-Leffler functions multiplied by the Fourier transforms $\hat{f}(κ)$ and $\hat{g}(κ)$ of the boundary data, plus a Prabhakar integral of the source's Fourier transform $\hat{\Phi}(κ, y)$. These are solution representations, not closed-form expressions: they reduce the PDE to quadratures over κ, but evaluating $N(x,y)$ for a given $(x,y)$ still requires computing an integral for each point. The paper provides explicit special-function forms (Fox H-functions, Wright function series) only when the boundary data are Dirac deltas ($f(x) = δ(x)$ or 0, $g(x) = δ(x)$ or 0) and the source is a separable product of a Dirac delta in $x$ and a power law in $y$ (Corollary 1, Examples 1-6). For general $f(x)$, $g(x)$, or $\Phi(x,y)$, no further simplification is achieved.

The consequence. A practitioner with non-impulse boundary conditions — e.g., $f(x) = e^{-x^2}$ (Gaussian initial profile), $f(x) = \mathbb{1}_{[-L,L]}(x)$ (finite-width source), or arbitrary measured data — cannot use the paper's closed-form Fox H-function results. They must evaluate the one-dimensional Fourier inversion integrals numerically for each $(x,y)$ of interest, which requires (1) computing the Fourier transform $\hat{f}(κ) = \int_{-\infty}^\infty f(x) e^{iκx} dx$ (potentially analytically for simple functions, otherwise numerically), (2) evaluating the Mittag-Leffler function $E_{\mu,\beta}(\pm y^\mu \psi_\alpha^\theta(κ))$ (itself requiring series summation or numerical methods), (3) computing the product $\hat{f}(κ) E_{\mu,\beta}(\cdots)$, and (4) integrating over κ from $−\infty$ to $\infty$ — an oscillatory integral whose numerical cost scales with the desired resolution in $x$. This is a substantially harder computational problem than evaluating a Fox H-function at specific points, and the paper provides no guidance on numerical methods, contour deformation, or convergence rates for these integrals.

What evidence exists in the paper. The Fox H-function results (equations 31, 32, 35, 36, 47, 48, 50, 53, 54) all correspond to $f(x) = δ(x)$ and $g(x) = 0$, or $f(x) = 0$ and $g(x) = δ(x)$, with $Φ(x,y) = δ(x)δ(y)$ or $Φ(x,y) = δ(x) y^{-β}/Γ(1-β)$. The integrals are evaluated in closed form because $\mathcal{F}[δ(x)](κ) = 1$, eliminating the $κ$-dependence of the boundary data and leaving an integral of $E_{\mu,\beta}(\pm y^\mu |κ|^\alpha) e^{−iκx}$ over $κ$, which the Mellin-cosine formula (86) handles. For arbitrary $f$ and $g$, the factor $\hat{f}(κ)$ remains inside the $κ$-integral, and no Mellin-cosine reduction is possible unless $\hat{f}(κ)$ itself has a special form (e.g., a power law or exponential) that combines with the Mittag-Leffler function into a recognizable Mellin-Barnes integral. The paper never claims to have solved the general-data case in closed form — the integral representations (25, 39, 63) are presented as the final answer.

Mitigation status. The paper is transparent about this: the solutions for general boundary data are presented as integrals, not as closed forms. This is inherent to the Fourier transform method — the solution is obtained as a superposition of Fourier modes, and the inverse transform can be evaluated explicitly only for special data. The paper does not claim otherwise. However, the abstract's language — "derive analytical solutions" — may create an expectation of explicit formulas for arbitrary $f$, $g$, $Φ$, which the paper does not provide. The mitigation is partial: for the important special case of impulse boundary conditions (which give the fundamental solution / Green's function), the paper does provide complete closed forms, asymptotic expansions, and series representations. For practitioners whose boundary data can be approximated as sums of impulses (via discretization or basis expansion), the impulse solutions serve as building blocks for linear superposition, since the PDE is linear. The paper does not discuss this superposition strategy explicitly, but it is the standard approach for using Green's functions with general data.

---

6.6 No Discussion of Solution Regularity, Physical Admissibility, or Well-Posedness

The assumption or constraint. The paper derives solutions via formal integral transforms without addressing the functional-analytic properties that would determine whether the solutions are physically admissible or even mathematically well-defined for all parameter choices. Questions that a practitioner or applied mathematician would ask — but the paper does not address — include: For which parameter ranges (α, μ, ν, θ) is the solution $N(x,y)$ in $L^2(\mathbb{R} \times \mathbb{R}^+)$ or $L^1$ (finite energy, normalizable)? Does the solution satisfy the boundary conditions in a pointwise sense, or only in a limiting distributional sense? Are the fractional integrals $(yI_{0+}^{(1−ν)(2−μ)} N)(x,0+)$ well-defined given the solution's behavior near $y = 0^+$? Does the source term produce a physically admissible response (e.g., positive temperature for a diffusion problem), or does the fractional derivative's non-locality introduce negativity or oscillations? Is the initial-boundary value problem well-posed (unique solution, continuous dependence on data) for the Hilfer-composite boundary conditions?

The consequence. A practitioner who naively evaluates the solution formulas for certain parameter regimes may obtain mathematically valid but physically nonsensical results — e.g., a "probability density" that goes negative, a "temperature" that diverges at the boundary, or a solution that fails to satisfy the boundary conditions as specified. The paper's asymptotic formulas (equation 51) show that the solution decays as a stretched exponential for the α = 2, θ = 0, quantum Riesz-Feller case — but there is no discussion of whether the solution is everywhere positive. The series representation (equation 52) has alternating signs $(-1)^j$, which does not guarantee positivity for all $x$ (the function could oscillate or cross zero). For the classical Riesz-Feller case (upper signs in equation 31), the Mittag-Leffler argument is $-y^\mu |κ|^\alpha$, and the behavior of $E_{\mu,\beta}(−z)$ for large positive $z$ is oscillatory and can take negative values — is the solution still physically interpretable as a steady-state field? The paper does not ask or answer these questions.

What evidence exists in the paper. There is no discussion of solution regularity, positivity, well-posedness, or functional setting anywhere in the paper. The boundary conditions (24a) are specified as equalities, but the paper does not verify that the derived solutions actually satisfy them — the verification is at the level of formal transform manipulation (the Laplace transform step assumes the boundary terms, and the inverse transform recovers them by construction), not at the level of pointwise limits as $y → 0^+$. The parameter ranges are specified (1 < α ≤ 2, 1 < μ ≤ 2, |θ| ≤ min{α, 2−α}) but without commentary on whether the solution formula is valid at the boundaries of these ranges (e.g., α = 2, where the Fox H-function reduction involves cancellation of coincident poles that may require careful limiting arguments).

Mitigation status. No mitigation is attempted. The paper operates entirely within the formal transform calculus tradition, where integral transforms are applied to the PDE and the resulting expressions are manipulated algebraically without verifying that the functions belong to the domains of the transforms or that the inversion formulas converge in a classical sense. This is standard practice in the special-functions approach to fractional PDEs — the monographs by Podlubny (1999), Kilbas, Srivastava, and Trujillo (2006), and Mathai, Saxena, and Haubold (2010) all treat solutions at this formal level — and the paper is consistent with its subfield's norms. However, the absence of any well-posedness or regularity analysis means that the paper's results are solution formulas whose domain of validity is not rigorously established. A practitioner encountering unexpected behavior (divergence, negativity, oscillations) for their parameter choices would need to consult the functional-analytic fractional PDE literature for well-posedness results specific to the Hilfer derivative (some of which existed at the time, e.g., in Hilfer, Luchko, and Tomovski, 2009, cited in the paper as Ref. [14]) to determine whether their parameter regime is theoretically sound.

7. Implications and Future Directions

How This Work Changes the Landscape

This paper does not introduce a new paradigm or challenge foundational assumptions in fractional calculus. Rather, it achieves something that is in some ways more practically useful: it completes a piece of the analytical catalog for fractional PDEs by solving a class of steady-state equations with mixed fractional derivative types that the prior literature had not systematically treated. The contribution is an infrastructure contribution — it provides the closed-form fundamental solutions, special function representations, asymptotic expansions, and mapping relationships that future researchers and practitioners need to use these equations as models of physical transport — rather than a disruptive conceptual insight. The magnitude is incremental but genuine: before this paper, if you wanted to model steady-state anomalous transport with Lévy flights in one direction and memory effects in the orthogonal direction, you either solved the time-dependent problem and took the long-time limit numerically, or you lacked analytical guidance entirely. After this paper, you have a Fourier-Laplace transform pipeline with explicit Mittag-Leffler and Fox H-function representations, plus the asymptotic behavior in two complementary regimes.

Where the paper does shift the landscape is in its unification of previously separate solution threads. Prior to this work, the Hilfer-composite fractional derivative was studied almost exclusively in time-fractional problems (Hilfer, 2000; Sandev, Metzler, and Tomovski, 2011), while the Riesz-Feller derivative was studied in space-fractional diffusion (Mainardi, Pagnini, and Saxena, 2005) or in combined space-time fractional equations where the Hilfer derivative acted in time and the Riesz-Feller derivative acted in the single spatial variable (Tomovski, Sandev, Metzler, and Dubbeldam, 2012). No one had combined them in a genuine steady-state two-variable spatial problem. By doing so, the paper demonstrates that the Hilfer derivative's interpolation parameter ν has a concrete, analytically tractable effect on the solution's behavior — specifically, ν controls the power-law prefactor exponents in both the boundary-data terms (the $y^{-(1-\nu)(2-\mu)}$ and $y^{1-(1-\nu)(2-\mu)}$ weights) and the asymptotic decay amplitudes (equation 51, where ν appears in the prefactor through $(\mu/2)^{(1-2\nu+1/(2-\mu))/2}$), but does not alter the functional form of the stretched exponential decay in the far field. This separation of roles — ν affects amplitudes but not the functional shape of the decay — is not obvious from the operator definition alone and is an empirical finding of the asymptotic analysis. For the field, this means that ν and μ are not redundant parameters: μ controls the qualitative character of the solution (Gaussian vs. stretched exponential), while ν fine-tunes the quantitative amplitude scaling without changing the character. This is useful guidance for parameter estimation in experimental contexts.

The paper also resolves a subtle confusion in the prior literature about the relationship between steady-state and time-dependent fractional solutions. The fractional calculus community had treated time-fractional wave equations and steady-state fractional Helmholtz equations as separate problem classes requiring separate derivations. The paper identifies — and exploits — the fact that the quantum Riesz-Feller derivative (symbol $+\psi_\alpha^\theta(\kappa)$) in a steady-state equation is structurally identical to the classical Riesz-Feller derivative (symbol $-\psi_\alpha^\theta(\kappa)$) in a time-dependent wave equation, via the correspondence $xD_\theta^{*,\alpha} = -xD_\theta^\alpha$. This means that every steady-state solution derived for the quantum Riesz-Feller case is simultaneously a time-dependent wave equation solution for the classical Riesz-Feller case, with $y \to t$ and the boundary conditions reinterpreted as initial conditions. Prior work (Samuel and Thomas, 2010; Tomovski and Sandev, 2010, 2011) had solved these problems independently without recognizing the structural identity. The paper's explicit mapping (Remarks 8 and 10) eliminates the need for future researchers to re-derive time-dependent solutions from scratch when the steady-state solution is already known — a genuine conceptual economy.

Beyond these specific findings, the paper redirects research attention in two ways. First, it makes the asymptotic analysis of Fox H-function solutions — not just their symbolic expression — the standard for what constitutes a "complete" solution of a fractional PDE. Many prior papers in this tradition stopped at writing the solution as an H-function, treating the H-function itself as the endpoint. This paper demonstrates, via the systematic application of the Braaksma large-argument formula (equation 87) and the residue series expansion (equation 84), that the Fox H-function is a gateway to physical insight, not a terminal notation: the stretched exponential decay law $N \propto \exp(-C |x|^{2/(2-\mu)} y^{-\mu/(2-\mu)})$ in equation (51) is not visible from the $H_{3,3}^{2,1}$ parameter list alone; it requires asymptotic extraction. Future work in this subfield should expect not just H-function representations but the explicit extraction of leading-order physical behavior.

Second, the paper implicitly identifies the boundary of the Fourier-Laplace transform methodology by showing where it succeeds (symmetric Riesz, θ = 0, Laplace/Poisson case) and where it falls short (θ ≠ 0, Helmholtz case with k² shift). The fact that all closed-form Fox H-function results require θ = 0 and k = 0 is not presented as a limitation, but it delineates the territory: the Mellin-cosine transform formula (86) requires homogeneity in κ, which asymmetric skewness and nonzero wave numbers both break. This tells future researchers that if they need solutions for θ ≠ 0 or k ≠ 0, they must either develop new integral transform identities (e.g., Mellin transforms with complex parameters for the asymmetric case, or perturbation expansions in k² around the k = 0 solution) or turn to numerical methods. The paper does not solve these problems, but it clearly marks them as unsolved.

Follow-Up Research This Work Enables

Numerical evaluation and validation of the Fox H-function solutions with public code release. The paper expresses the fundamental solutions for the Laplace/Poisson case (α = 2, θ = 0) as $H_{1,1}^{1,0}$ and $H_{3,3}^{2,1}$ Fox H-functions and provides asymptotic expansions for large and small $|x|/y^{\mu/2}$. But there is no numerical evaluation of any of these expressions — no plots, no tables, no verification against Fourier integral quadrature. A strong follow-up would implement the Mellin-Barnes contour integral (equation 83) or the residue series (equation 84) numerically for a grid of parameter values (μ ∈ {1.1, 1.3, 1.5, 1.7, 1.9}, ν ∈ {0, 0.25, 0.5, 0.75, 1}) and compare the resulting N(x,y) profiles against (a) direct numerical integration of the Fourier inversion integral (equation 25) as a ground-truth check, and (b) the asymptotic formulas (equations 51 and 52) to determine empirically at what values of $|x|/y^{\mu/2}$ each asymptotic regime becomes accurate to within 5% relative error. The output would be a benchmark dataset and open-source Python/Matlab code that the fractional calculus community can use to evaluate solutions without reimplementing the H-function from scratch — a practical enabling contribution the original paper does not provide.

Asymmetric Riesz-Feller case: develop Mellin transform identities with complex parameters. The paper's inability to produce Fox H-function closed forms for θ ≠ 0 is the most natural next analytical target. For the asymmetric case, the Mittag-Leffler argument contains $e^{i \cdot \text{sign}(\kappa) \theta\pi/2}$, which splits the inverse Fourier integral into cosine and sine transforms: $\int_{-\infty}^\infty E_{\mu,\beta}(\pm y^\mu |\kappa|^\alpha e^{i \cdot \text{sign}(\kappa) \theta\pi/2}) e^{-i\kappa x} d\kappa$ separates into even and odd parts involving $\cos(\kappa x)$ and $\sin(\kappa x)$. The cosine part involves an integrand that depends on $|\kappa|$ but with a complex multiplicative constant inside the Mittag-Leffler argument — which may be evaluable via a generalization of the Mellin-cosine formula (86) to complex H-function arguments using the Mellin transform identity $\int_0^\infty \kappa^{\rho-1} \cos(\kappa x) H(\omega \kappa^\delta) d\kappa$ with $\omega \in \mathbb{C}$. Such a formula would produce Fox H-functions with complex parameters, which are less studied but potentially tractable using the same residue calculus. A follow-up that derives the analogous H-function representation for θ ≠ 0 — even if restricted to the α = 2 case for simplicity — would close the gap left by Sections III and IV and allow the asymptotic machinery of the paper to be applied to asymmetric Lévy flights.

Helmholtz equation: systematic perturbation theory in k² around the Poisson solution. The Helmholtz solutions (Theorem 2) are left as Fourier integrals because the shift $\psi_\alpha^\theta - k^2$ breaks the homogeneity needed for the Mellin-cosine formula. A natural attack is to expand the Helmholtz solution as a power series in $k^2$:

$N_{\text{Helmholtz}}(x,y) = N_{\text{Poisson}}(x,y) + k^2 N^{(1)}(x,y) + k^4 N^{(2)}(x,y) + \cdots,$

where each $N^{(j)}$ is obtained by differentiating the Poisson solution with respect to the H-function argument. Concretely, the $f$-term in the Helmholtz solution involves $E_{\mu,\beta}(y^\mu(\psi_\alpha^\theta - k^2))$, and the Taylor expansion $E_{\mu,\beta}(a - \epsilon) = \sum_{n=0}^\infty \frac{(-1)^n}{n!} \epsilon^n E_{\mu,\beta}^{(n)}(a)$ expresses the Helmholtz solution as a series of derivatives of Mittag-Leffler functions evaluated at the Poisson argument $\psi_\alpha^\theta$. Each derivative $E_{\mu,\beta}^{(n)}$ is itself expressible in terms of higher-parameter Mittag-Leffler functions or Fox H-functions via the series representation (79), and the inverse Fourier transforms of these terms may reduce to Fox H-functions of the same $H_{3,3}^{2,1}$ type as the Poisson case but with shifted parameters. A follow-up that carries this out to first or second order in k² would provide the first analytical insight into how the wave number modifies the stretched exponential decay of the Poisson solution — specifically, whether the k² correction introduces oscillatory behavior, changes the decay exponent, or simply renormalizes the amplitude. For physical applications (viscoelastic membranes, fractional Schrödinger equation), this is the missing piece that prevents the paper's results from being directly usable.

Application to parameter estimation: fitting experimental Lévy flight + memory data. The paper's asymptotic formula (51) provides a specific functional form for the lateral decay of the steady-state field: $\log N \sim -C(\mu) (|x|/y^{\mu/2})^{2/(2-\mu)}$ plus a logarithmic correction from the power-law prefactor. This is a testable prediction. An experimental or computational follow-up could take published data on two-dimensional anomalous transport — for example, the spatial distribution of tracer particles in a rotating turbulent plasma where radial transport is Lévy-like (Riesz-Feller in the radial coordinate) and poloidal transport is subdiffusive due to trapping in magnetic islands (Hilfer in the angular coordinate, after appropriate coordinate mapping) — and fit the observed steady-state density profiles to the stretched exponential form to estimate μ and ν. The paper's analytical separation of roles (μ controls the functional form, ν controls the amplitude scaling) provides a fitting strategy: first estimate μ from the slope of $\log(-\log N)$ vs. $\log(|x|)$, then estimate ν from the amplitude intercept. A null result — finding that the experimental data does not follow the stretched exponential form — would rule out the Hilfer/Riesz-Feller model for that system and motivate investigation of alternative fractional derivative definitions (e.g., distributed-order, Caputo-Fabrizio, or tempered fractional derivatives) for which analogous analytical solutions are not yet available. This is the kind of experiment-theory feedback loop that gives analytical solution papers their lasting value.

General boundary data: develop a computational method based on the impulse-response Green's function. The paper's closed-form Fox H-function results all assume $f(x) = δ(x)$ and $g(x) = 0$ (or vice versa). These are the Green's functions for the fractional Laplace/Poisson operators. Since the PDE is linear, the solution for arbitrary boundary data $f(x), g(x)$ is the convolution of the Green's function with the data: $N(x,y) = \int_{-\infty}^\infty G_f(x-x', y) f(x') dx' + \int_{-\infty}^\infty G_g(x-x', y) g(x') dx'$, where $G_f$ and $G_g$ are given by the $H_{1,1}^{1,0}$ expressions in equations 50 and 54 (for α = 2, θ = 0). A strong follow-up would implement a fast numerical convolution scheme — exploiting the fact that the H-function depends only on $|x-x'|/y^{\mu/2}$, making the convolution efficiently computable via FFT for uniformly sampled data or via adaptive quadrature for irregular data — and benchmark it against direct finite-difference or spectral discretizations of the original fractional PDE. This would transform the paper's fundamental solutions from analytical curiosities into practical computational tools: given arbitrary $f(x)$ on a grid, compute $N(x,y)$ at any $y > 0$ by convolving with the precomputed H-function kernel, avoiding the need to solve the full PDE numerically for each new boundary condition. The computational complexity comparison (FFT-based convolution vs. PDE discretization) would quantify the practical benefit of having the analytical Green's function.

Practical Applications and Downstream Use Cases

Benchmark solutions for numerical fractional PDE solvers. The most immediate practical use of this paper is as a source of exact reference solutions against which to validate numerical methods for fractional PDEs on two-dimensional domains. As of 2014, numerical methods for fractional equations — finite difference discretizations of the Riesz-Feller derivative, spectral methods, matrix transfer techniques for the Hilfer derivative — were being actively developed but lacked non-trivial analytical test cases with mixed derivative types. The impulse-response solutions in Examples 1 and 2 (equations 31, 32, 35, 36) provide exact solutions for the half-plane problem with point boundary data and point sources. A numerical solver developer can evaluate the Fox H-function at a grid of $(x,y)$ points (using the residue series or Mellin-Barnes contour), compare their numerical solution on the same grid, and compute the convergence rate as the mesh is refined. The availability of both large-|x| asymptotic formulas (equation 51) and small-|x| series (equation 52) means the boundary conditions at the truncation of a finite computational domain can be set using the known far-field behavior, avoiding the artificial reflections or zero-padding that contaminate convergence studies. This is infrastructure for the numerical fractional PDE community — the paper provides the analytical "ground truth" that the community needs but did not previously have for this operator class.

Parameterized reduced-order models for two-dimensional anomalous transport in complex materials. In engineering applications where the steady-state temperature, concentration, or displacement field in a two-dimensional medium is needed repeatedly for different material parameters — for example, optimizing the design of a composite material with viscoelastic memory in one direction and Lévy-flight thermal transport in the orthogonal direction — the paper's closed-form solutions enable a parameterized surrogate model that evaluates $N(x,y; \mu, \nu, \alpha, \theta)$ directly without solving a PDE each time. Specifically, the Fox H-function representation (equation 31) depends on the spatial variables only through the ratio $|x|^\alpha / y^\mu$ (for the symmetric case), and the H-function parameters encode all the fractional orders. A design optimization loop over μ and ν can precompute the H-function on a coarse grid of parameter values, interpolate, and evaluate the solution at millions of spatial points with the cost of evaluating an algebraic scaling relation — many orders of magnitude faster than finite-element discretizations of the fractional PDE for each parameter candidate. The asymptotic formula (51) further simplifies this for the far-field region, where only the leading-order stretched exponential needs to be computed. This use case depends on having reliable numerical evaluation of the Fox H-function, which the paper does not provide but which subsequent work (as proposed in the "Follow-Up Research" section above) would enable.

Interpretation of lateral spreading in geological or biological Lévy flight data with spatial gradients. In field studies of animal foraging or contaminant transport where the spatial distribution of individuals or particles is measured along a transect at various distances from a source, the observed lateral profiles often deviate from Gaussian predictions. If the physical mechanism is hypothesized to involve Lévy flights (heavy-tailed step lengths) in the lateral direction combined with memory-affected diffusion in the longitudinal (downstream) direction — a plausible model for, say, pollutant transport in a river with dead zones that trap particles for power-law-distributed waiting times — the paper's stretched exponential formula $N(x,y) \propto \exp(-C |x|^{2/(2-\mu)} y^{-\mu/(2-\mu)})$ provides a specific functional form with two interpretable parameters (μ for the memory strength, and an overall amplitude involving ν) that can be fit to transect data. The prediction is distinguishable from alternative models: pure Lévy flight (symmetric Riesz with no memory) would give a different decay, and pure subdiffusion (ordinary Laplacian in x with Hilfer y-derivative) would be Gaussian in x. Fitting the paper's analytical form to data would allow an experimentalist to reject or support the mixed Lévy-memory hypothesis and to quantify the fractional orders μ, ν from observations. This is a direct translation of the asymptotic analysis into an experimental protocol.