Time-Fractional Black–Scholes Equation

• TFBSE is a fractional generalization of the classical Black–Scholes PDE that replaces the time derivative with a Caputo derivative to incorporate history-dependent dynamics.
• It is derived through probabilistic subordination, linking the option price to inverse subordinator processes and accommodating nonlocal temporal effects.
• Numerical methods like finite differences, compact schemes, spline collocation, and neural solvers are used to efficiently tackle the nonlocal memory and weak initial singularities in TFBSE.

Searching arXiv for recent and foundational papers on the Time-Fractional Black–Scholes Equation and related variants. The Time-Fractional Black–Scholes Equation (TFBSE) is a fractional-in-time generalization of the Black–Scholes option-pricing PDE in which the first-order time derivative is replaced by a nonlocal fractional derivative, most commonly of Caputo type or an equivalent modified Riemann–Liouville form under suitable regularity. In the TFBSE, the spatial Black–Scholes operator is typically retained while temporal evolution acquires memory through a convolution kernel, so current option values depend on the full past history of the solution rather than only its instantaneous state. Across the literature, the TFBSE appears in several mathematically distinct but related forms: as an ad hoc time-fractional extension of the classical pricing PDE, as the pricing equation associated with subdiffusive or tempered-subdiffusive asset dynamics, and as the target of a wide range of numerical schemes including finite differences, compact methods, spline collocation, meshless RBF methods, neural solvers, and transform-based series representations (Krzyżanowski et al., 2019, Dimitrov et al., 2016, Singh et al., 2022, Song et al., 2021, Zhou et al., 2023, Garg et al., 30 Jan 2026). A persistent point of clarification is that not every “fractional Black–Scholes” model is a TFBSE: some papers introduce fractionality through stochastic volatility, long memory, or roughness in the volatility factor rather than through a time-fractional pricing operator (Garnier et al., 2015).

1. Definition and model variants

In its standard one-asset form, the TFBSE replaces the classical Black–Scholes time derivative by a fractional derivative of order in $(0,1)$. A representative formulation is

$${}_{0}^{c}D_t^\alpha v(z,t) = -\frac{1}{2}\sigma^2 z^2 \frac{\partial^2 v(z,t)}{\partial z^2} -rz\frac{\partial v(z,t)}{\partial z} +r v(z,t),$$

with European call terminal and boundary conditions (Krzyżanowski et al., 2019). After time reversal and log-price transformation, this becomes a constant-coefficient convection–diffusion–reaction equation,

$${}_{0}^{c}D_t^\alpha u(x,t) = \frac{1}{2}\sigma^2 u_{xx}(x,t)+\left(r-\frac{1}{2}\sigma^2\right)u_x(x,t)-ru(x,t),$$

which is the form used by many numerical papers (Krzyżanowski et al., 2019, Song et al., 2021, Singh et al., 2022, Garg et al., 30 Jan 2026).

A closely related form appears when the original problem is posed backward from maturity using a modified right Riemann–Liouville derivative and then transformed into a forward problem. In such settings, the transformed derivative is treated as equivalent to the Caputo derivative under $C^{(1)}$ temporal regularity (V et al., 9 Aug 2025, Garg et al., 30 Jan 2026, Singh et al., 2022). This equivalence is important because Caputo derivatives accommodate classical payoff-based initial conditions.

Several extensions of the TFBSE are represented in the literature. A two-asset time-fractional Black–Scholes PDE with Caputo derivative and mixed derivative term has been written as

$\prescript{C}{}{D}_t^\alpha c(S_1,S_2,t) +\frac{1}{2}\sigma_1^2 S_1^2 c_{S_1S_1} +\frac{1}{2}\sigma_2^2 S_2^2 c_{S_2S_2} +\rho \sigma_1 \sigma_2 S_1S_2 c_{S_1S_2} +rS_1 c_{S_1} +rS_2 c_{S_2} -r c=0,$

with $0<\alpha\le 1$ (Zakaria et al., 2020). Tempered variants replace the standard Caputo derivative by a tempered operator, typically to weaken the heaviest memory tails while retaining nonlocality (Zhou et al., 2023, Krzyżanowski et al., 2021). Time-space-fractional Black–Scholes-type models generalize further by using Caputo time derivatives together with radial Riemann–Liouville spatial derivatives, so the pure TFBSE is recovered only as a specialization (Torres-Hernandez et al., 2020).

Not all “fractional Black–Scholes” papers belong to this class. In particular, a model in which volatility is driven by a fractional Ornstein–Uhlenbeck process built from fractional Brownian motion yields fractional-power implied-volatility corrections to classical Black–Scholes, but it does not replace the time derivative in the pricing PDE by a fractional operator (Garnier et al., 2015). That distinction is foundational for the taxonomy of the field.

2. Fractional derivatives and probabilistic interpretation

The dominant operator in the TFBSE literature is the Caputo derivative. In the one-dimensional setting, it is written as

$${}_{0}^{c}D_t^\alpha g(t) = \frac{1}{\Gamma(1-\alpha)} \int_0^t g'(s)(t-s)^{-\alpha}\,ds, \qquad 0<\alpha<1,$$

or in equivalent modified Riemann–Liouville form after appropriate time reversal (Krzyżanowski et al., 2019, V et al., 9 Aug 2025, Singh et al., 2022). The Caputo choice is standard because it preserves the interpretation of initial conditions in terms of option payoffs.

A more structural interpretation arises in subdiffusive and tempered-subdiffusive Black–Scholes models. In the subdiffusive case, the asset process is obtained by replacing the physical clock of geometric Brownian motion with an inverse $\alpha$-stable subordinator,

$Z_\alpha(t)=Z(S_\alpha(t)),$

where $S_\alpha$ is independent of Brownian motion (Krzyżanowski et al., 2019). The resulting option price can be expressed as

$${}_{0}^{c}D_t^\alpha v(z,t) = -\frac{1}{2}\sigma^2 z^2 \frac{\partial^2 v(z,t)}{\partial z^2} -rz\frac{\partial v(z,t)}{\partial z} +r v(z,t),$$0

where $${}_{0}^{c}D_t^\alpha v(z,t) = -\frac{1}{2}\sigma^2 z^2 \frac{\partial^2 v(z,t)}{\partial z^2} -rz\frac{\partial v(z,t)}{\partial z} +r v(z,t),$$1 is the classical Black–Scholes price. Laplace-transform arguments then yield the Caputo TFBSE as the governing pricing equation (Krzyżanowski et al., 2019). In this reading, the fractional time derivative is not merely phenomenological; it is the PDE manifestation of random waiting times and market stagnation encoded by inverse-subordinator time change.

A more general construction uses the inverse of a strictly increasing subordinator with Laplace exponent

$${}_{0}^{c}D_t^\alpha v(z,t) = -\frac{1}{2}\sigma^2 z^2 \frac{\partial^2 v(z,t)}{\partial z^2} -rz\frac{\partial v(z,t)}{\partial z} +r v(z,t),$$2

leading to a generalized time-fractional operator

$${}_{0}^{c}D_t^\alpha v(z,t) = -\frac{1}{2}\sigma^2 z^2 \frac{\partial^2 v(z,t)}{\partial z^2} -rz\frac{\partial v(z,t)}{\partial z} +r v(z,t),$$3

(Zhang et al., 13 Nov 2025). In this framework, the option-pricing kernel solves

$${}_{0}^{c}D_t^\alpha v(z,t) = -\frac{1}{2}\sigma^2 z^2 \frac{\partial^2 v(z,t)}{\partial z^2} -rz\frac{\partial v(z,t)}{\partial z} +r v(z,t),$$4

and the TFBSE is derived from a rigorously specified arbitrage-free but generally incomplete sub-diffusive market model rather than postulated directly (Zhang et al., 13 Nov 2025). This suggests that the standard Caputo TFBSE is a special case within a broader inverse-subordinator theory.

Tempered subdiffusive models modify the inverse stable clock by using an inverse tempered stable subordinator, leading to the tempered operator

$${}_{0}^{c}D_t^\alpha v(z,t) = -\frac{1}{2}\sigma^2 z^2 \frac{\partial^2 v(z,t)}{\partial z^2} -rz\frac{\partial v(z,t)}{\partial z} +r v(z,t),$$5

which reduces to the ordinary TFBSE operator when $${}_{0}^{c}D_t^\alpha v(z,t) = -\frac{1}{2}\sigma^2 z^2 \frac{\partial^2 v(z,t)}{\partial z^2} -rz\frac{\partial v(z,t)}{\partial z} +r v(z,t),$$6 (Krzyżanowski et al., 2021). Here tempering weakens the most extreme heavy-tail effects while preserving temporal nonlocality.

3. Transformations and canonical computational forms

The standard computational treatment of the TFBSE begins with time reversal and logarithmic transformation. Starting from a terminal-value problem in $${}_{0}^{c}D_t^\alpha v(z,t) = -\frac{1}{2}\sigma^2 z^2 \frac{\partial^2 v(z,t)}{\partial z^2} -rz\frac{\partial v(z,t)}{\partial z} +r v(z,t),$$7, one sets

$${}_{0}^{c}D_t^\alpha v(z,t) = -\frac{1}{2}\sigma^2 z^2 \frac{\partial^2 v(z,t)}{\partial z^2} -rz\frac{\partial v(z,t)}{\partial z} +r v(z,t),$$8

or equivalent notation, thereby converting the problem into an initial-value equation on a log-price domain (Krzyżanowski et al., 2019, V et al., 9 Aug 2025, Singh et al., 2022, Song et al., 2021, Garg et al., 30 Jan 2026). This removes the variable coefficients $${}_{0}^{c}D_t^\alpha v(z,t) = -\frac{1}{2}\sigma^2 z^2 \frac{\partial^2 v(z,t)}{\partial z^2} -rz\frac{\partial v(z,t)}{\partial z} +r v(z,t),$$9 and $${}_{0}^{c}D_t^\alpha u(x,t) = \frac{1}{2}\sigma^2 u_{xx}(x,t)+\left(r-\frac{1}{2}\sigma^2\right)u_x(x,t)-ru(x,t),$$0 from the Black–Scholes operator and yields a constant-coefficient PDE in $${}_{0}^{c}D_t^\alpha u(x,t) = \frac{1}{2}\sigma^2 u_{xx}(x,t)+\left(r-\frac{1}{2}\sigma^2\right)u_x(x,t)-ru(x,t),$$1.

For bounded-domain computations, the infinite log-price interval is truncated to $${}_{0}^{c}D_t^\alpha u(x,t) = \frac{1}{2}\sigma^2 u_{xx}(x,t)+\left(r-\frac{1}{2}\sigma^2\right)u_x(x,t)-ru(x,t),$$2 or $${}_{0}^{c}D_t^\alpha u(x,t) = \frac{1}{2}\sigma^2 u_{xx}(x,t)+\left(r-\frac{1}{2}\sigma^2\right)u_x(x,t)-ru(x,t),$$3, and Dirichlet conditions are imposed based on asymptotic option behavior (Krzyżanowski et al., 2019, V et al., 9 Aug 2025, Singh et al., 2022). In some approaches, nonhomogeneous boundary values are eliminated by an affine lifting. For example, one may define

$${}_{0}^{c}D_t^\alpha u(x,t) = \frac{1}{2}\sigma^2 u_{xx}(x,t)+\left(r-\frac{1}{2}\sigma^2\right)u_x(x,t)-ru(x,t),$$4

followed by a transformed unknown $${}_{0}^{c}D_t^\alpha u(x,t) = \frac{1}{2}\sigma^2 u_{xx}(x,t)+\left(r-\frac{1}{2}\sigma^2\right)u_x(x,t)-ru(x,t),$$5 satisfying homogeneous boundary conditions (Zhou et al., 2023, Song et al., 2021). In the tempered case, an additional exponential integrating factor is used: $${}_{0}^{c}D_t^\alpha u(x,t) = \frac{1}{2}\sigma^2 u_{xx}(x,t)+\left(r-\frac{1}{2}\sigma^2\right)u_x(x,t)-ru(x,t),$$6 which transforms the Black–Scholes operator into a diffusion–reaction equation with homogeneous boundaries (Zhou et al., 2023).

An older high-order compact treatment uses a different route for a European put. After time reversal and boundary homogenization,

$${}_{0}^{c}D_t^\alpha u(x,t) = \frac{1}{2}\sigma^2 u_{xx}(x,t)+\left(r-\frac{1}{2}\sigma^2\right)u_x(x,t)-ru(x,t),$$7

and then

$${}_{0}^{c}D_t^\alpha u(x,t) = \frac{1}{2}\sigma^2 u_{xx}(x,t)+\left(r-\frac{1}{2}\sigma^2\right)u_x(x,t)-ru(x,t),$$8

to eliminate the first derivative term and obtain a diffusion-form TFBSE involving $${}_{0}^{c}D_t^\alpha u(x,t) = \frac{1}{2}\sigma^2 u_{xx}(x,t)+\left(r-\frac{1}{2}\sigma^2\right)u_x(x,t)-ru(x,t),$$9 and a reaction term (Dimitrov et al., 2016). This diffusion form is particularly convenient for nonuniform spatial meshes tailored to the degenerate left boundary and strike-region nonsmoothness.

These transformations are not mere algebraic conveniences. They determine the structure of the discrete operator, the appropriate boundary treatment, and the attainable order of accuracy, especially when combined with graded temporal meshes designed to handle startup singularities (Song et al., 2021, Zhou et al., 2023).

4. Analytical solutions, transform methods, and formal series

Closed-form solutions to the TFBSE are generally unavailable except through subordinated semigroup representations, Laplace-domain formulas, or formal series. In the subdiffusive setting, the option price can be written as the classical Black–Scholes price averaged over the law of the inverse subordinator: $C^{(1)}$0 which directly links the TFBSE to the classical model (Krzyżanowski et al., 2019). The more general sub-diffusive formulation yields

$C^{(1)}$1

where $C^{(1)}$2 is geometric Brownian motion in operational time and $C^{(1)}$3 is the inverse subordinator, and this function is proved to be the unique solution of the generalized TFBSE (Zhang et al., 13 Nov 2025).

For a European call in the sub-diffusive model, the value at time zero can be written as

$C^{(1)}$4

which is an explicit mixture of Black–Scholes call values over the inverse-subordinator law (Zhang et al., 13 Nov 2025). A plausible implication is that generalized TFBSEs can often be solved more naturally by probabilistic subordination than by direct PDE methods.

Transform-based formal series also appear. A two-asset Caputo TFBSE has been treated using the Sumudu transform, leading to an integral equation and a recursive infinite-series solution

$C^{(1)}$5

with

$C^{(1)}$6

(Zakaria et al., 2020). The same recurrence can be compactly interpreted as an operator series, although the paper itself does not formalize that viewpoint (Zakaria et al., 2020). Its main limitation is that convergence is asserted but not rigorously proved.

ANN-based work also treats the TFBSE from a semidiscrete perspective. There the Caputo derivative is approximated by a history-weighted scheme, and each time step is solved as a spatial problem by a two-layer feed-forward network trained by Adam (Bajalan et al., 2021). However, the paper’s explicit fractional Black–Scholes example is a manufactured problem with internal typographical inconsistencies, so its value is methodological rather than canonical (Bajalan et al., 2021).

5. Numerical methods

A large fraction of TFBSE research is numerical, and the field is organized around the tradeoff between temporal nonlocality, spatial accuracy, boundary treatment, and the weak initial singularity typical of fractional evolution equations.

The most classical family uses $C^{(1)}$7-type approximations for the Caputo derivative on uniform meshes. In the subdiffusive Black–Scholes model, a weighted finite difference method interpolates between implicit, explicit, and fractional Crank–Nicolson analogues through a parameter $C^{(1)}$8, with the fully implicit case $C^{(1)}$9 unconditionally stable for all $\prescript{C}{}{D}_t^\alpha c(S_1,S_2,t) +\frac{1}{2}\sigma_1^2 S_1^2 c_{S_1S_1} +\frac{1}{2}\sigma_2^2 S_2^2 c_{S_2S_2} +\rho \sigma_1 \sigma_2 S_1S_2 c_{S_1S_2} +rS_1 c_{S_1} +rS_2 c_{S_2} -r c=0,$0 and the fractional Crank–Nicolson case conditionally stable (Krzyżanowski et al., 2019). The same paper proves convergence with order $\prescript{C}{}{D}_t^\alpha c(S_1,S_2,t) +\frac{1}{2}\sigma_1^2 S_1^2 c_{S_1S_1} +\frac{1}{2}\sigma_2^2 S_2^2 c_{S_2S_2} +\rho \sigma_1 \sigma_2 S_1S_2 c_{S_1S_2} +rS_1 c_{S_1} +rS_2 c_{S_2} -r c=0,$1 in time and $\prescript{C}{}{D}_t^\alpha c(S_1,S_2,t) +\frac{1}{2}\sigma_1^2 S_1^2 c_{S_1S_1} +\frac{1}{2}\sigma_2^2 S_2^2 c_{S_2S_2} +\rho \sigma_1 \sigma_2 S_1S_2 c_{S_1S_2} +rS_1 c_{S_1} +rS_2 c_{S_2} -r c=0,$2 in space under its stability conditions (Krzyżanowski et al., 2019). A tempered extension modifies the time operator and the convolution weights but retains the same general finite-difference architecture (Krzyżanowski et al., 2021).

High-order compact finite differences on nonuniform meshes constitute another major line. A three-point compact approximation on smooth nonuniform meshes generated by Tavella–Randall or quadratic mappings has fourth-order local spatial accuracy, and when applied to the transformed TFBSE yields observed $\prescript{C}{}{D}_t^\alpha c(S_1,S_2,t) +\frac{1}{2}\sigma_1^2 S_1^2 c_{S_1S_1} +\frac{1}{2}\sigma_2^2 S_2^2 c_{S_2S_2} +\rho \sigma_1 \sigma_2 S_1S_2 c_{S_1S_2} +rS_1 c_{S_1} +rS_2 c_{S_2} -r c=0,$3 behavior, with fourth-order spatial convergence and roughly first-order temporal accuracy in experiments due to the startup singularity (Dimitrov et al., 2016). This work is important because it ties mesh design directly to two structural features of Black–Scholes operators: degeneracy near $\prescript{C}{}{D}_t^\alpha c(S_1,S_2,t) +\frac{1}{2}\sigma_1^2 S_1^2 c_{S_1S_1} +\frac{1}{2}\sigma_2^2 S_2^2 c_{S_2S_2} +\rho \sigma_1 \sigma_2 S_1S_2 c_{S_1S_2} +rS_1 c_{S_1} +rS_2 c_{S_2} -r c=0,$4 and reduced regularity near the strike (Dimitrov et al., 2016).

Later high-order methods explicitly address the weak initial singularity. A nonuniform Alikhanov formula combined with a fourth-order average approximation in space and sum-of-exponentials acceleration produces a fast variable-step solver with unconditional stability and global convergence

$\prescript{C}{}{D}_t^\alpha c(S_1,S_2,t) +\frac{1}{2}\sigma_1^2 S_1^2 c_{S_1S_1} +\frac{1}{2}\sigma_2^2 S_2^2 c_{S_2S_2} +\rho \sigma_1 \sigma_2 S_1S_2 c_{S_1S_2} +rS_1 c_{S_1} +rS_2 c_{S_2} -r c=0,$5

so second-order temporal accuracy is recovered when the grading parameter satisfies $\prescript{C}{}{D}_t^\alpha c(S_1,S_2,t) +\frac{1}{2}\sigma_1^2 S_1^2 c_{S_1S_1} +\frac{1}{2}\sigma_2^2 S_2^2 c_{S_2S_2} +\rho \sigma_1 \sigma_2 S_1S_2 c_{S_1S_2} +rS_1 c_{S_1} +rS_2 c_{S_2} -r c=0,$6 (Song et al., 2021). A tempered analogue combines a compact fourth-order spatial operator, graded time mesh $\prescript{C}{}{D}_t^\alpha c(S_1,S_2,t) +\frac{1}{2}\sigma_1^2 S_1^2 c_{S_1S_1} +\frac{1}{2}\sigma_2^2 S_2^2 c_{S_2S_2} +\rho \sigma_1 \sigma_2 S_1S_2 c_{S_1S_2} +rS_1 c_{S_1} +rS_2 c_{S_2} -r c=0,$7, nonuniform tempered $\prescript{C}{}{D}_t^\alpha c(S_1,S_2,t) +\frac{1}{2}\sigma_1^2 S_1^2 c_{S_1S_1} +\frac{1}{2}\sigma_2^2 S_2^2 c_{S_2S_2} +\rho \sigma_1 \sigma_2 S_1S_2 c_{S_1S_2} +rS_1 c_{S_1} +rS_2 c_{S_2} -r c=0,$8 approximation, and SOE history compression, achieving temporal order $\prescript{C}{}{D}_t^\alpha c(S_1,S_2,t) +\frac{1}{2}\sigma_1^2 S_1^2 c_{S_1S_1} +\frac{1}{2}\sigma_2^2 S_2^2 c_{S_2S_2} +\rho \sigma_1 \sigma_2 S_1S_2 c_{S_1S_2} +rS_1 c_{S_1} +rS_2 c_{S_2} -r c=0,$9, fourth-order space, and unconditional stability for both direct and fast schemes (Zhou et al., 2023).

Spline-based collocation methods form a distinct class. Exponential B-spline collocation combined with an $0<\alpha\le 1$0-type time discretization yields a tridiagonal system at each time step, unconditional stability by von Neumann analysis, and convergence of order $0<\alpha\le 1$1 (Singh et al., 2022). The method has been applied not only to manufactured tests but also to European call, European put, and double barrier knock-out call options (Singh et al., 2022). A more recent exponential B-spline method uses a Crank–Nicolson-type approximation for the fractional derivative, claims unconditional stability, and reports temporal order $0<\alpha\le 1$2 and second-order spatial convergence on smooth test problems (Garg et al., 30 Jan 2026). Modified cubic B-spline differential quadrature has also been used together with $0<\alpha\le 1$3 time discretization, with reported fourth-order spatial convergence and order $0<\alpha\le 1$4 in time under a matrix inverse bound derived via the Neumann series theorem (V et al., 9 Aug 2025).

Meshless methods are represented by a global RBF collocation approach for time-space-fractional Black–Scholes-type PDEs. Specializing to the pure time-fractional case by taking $0<\alpha\le 1$5 yields a TFBSE solver using Caputo $0<\alpha\le 1$6-type time stepping, cubic polyharmonic spline interpolation, and a preconditioned global collocation matrix (Torres-Hernandez et al., 2020). The paper emphasizes flexibility in multidimensional settings and condition-number reduction, but it does not provide a full convergence theorem for the fully discrete RBF method (Torres-Hernandez et al., 2020).

ANN-based solvers, finally, replace the spatial linear solver by optimization over a neural network ansatz at each time step, with domain mapping for semi-infinite $0<\alpha\le 1$7-domains and fine-tuning across adjacent time levels (Bajalan et al., 2021). This suggests a broadening of TFBSE computation beyond classical discretizations, though the validation remains limited (Bajalan et al., 2021).

6. Applications, interpretation, and conceptual boundaries

The TFBSE has been used to price a range of European-style derivatives. Standard applications include European calls and puts under time-fractional dynamics (Krzyżanowski et al., 2019, Singh et al., 2022). More specialized applications include double barrier knock-out call options (Song et al., 2021, Singh et al., 2022) and two-asset basket-type payoffs in a Caputo two-dimensional Black–Scholes PDE (Zakaria et al., 2020). In the generalized sub-diffusive setting, the pricing theory is extended to European-style contingent claims under a sub-diffusion equivalent martingale measure, with the call price represented by a shifted TFBSE and an explicit inverse-subordinator mixture formula (Zhang et al., 13 Nov 2025).

A recurring interpretive theme is that the fractional order governs the strength of temporal memory or subdiffusion. In the subdiffusive Black–Scholes model, smaller $0<\alpha\le 1$8 corresponds to more pronounced waiting times or stagnation, interpreted as illiquidity or inactivity (Krzyżanowski et al., 2019). In the tempered subdiffusive extension, $0<\alpha\le 1$9 weakens the heaviest waiting-time tails, but the effect of a numerically small $${}_{0}^{c}D_t^\alpha g(t) = \frac{1}{\Gamma(1-\alpha)} \int_0^t g'(s)(t-s)^{-\alpha}\,ds, \qquad 0<\alpha<1,$$0 can remain significant when $${}_{0}^{c}D_t^\alpha g(t) = \frac{1}{\Gamma(1-\alpha)} \int_0^t g'(s)(t-s)^{-\alpha}\,ds, \qquad 0<\alpha<1,$$1 is small because $${}_{0}^{c}D_t^\alpha g(t) = \frac{1}{\Gamma(1-\alpha)} \int_0^t g'(s)(t-s)^{-\alpha}\,ds, \qquad 0<\alpha<1,$$2 need not be small (Krzyżanowski et al., 2021). This matters for both modeling and numerical interpretation.

At the level of computed prices, several papers report that the influence of $${}_{0}^{c}D_t^\alpha g(t) = \frac{1}{\Gamma(1-\alpha)} \int_0^t g'(s)(t-s)^{-\alpha}\,ds, \qquad 0<\alpha<1,$$3 is strongest near the money or in barrier-sensitive regions. In the exponential B-spline collocation study, the impact of $${}_{0}^{c}D_t^\alpha g(t) = \frac{1}{\Gamma(1-\alpha)} \int_0^t g'(s)(t-s)^{-\alpha}\,ds, \qquad 0<\alpha<1,$$4 on European call prices is mild far from the strike but significant near $${}_{0}^{c}D_t^\alpha g(t) = \frac{1}{\Gamma(1-\alpha)} \int_0^t g'(s)(t-s)^{-\alpha}\,ds, \qquad 0<\alpha<1,$$5, while for double barrier knock-out calls the price is reported to be inversely proportional to $${}_{0}^{c}D_t^\alpha g(t) = \frac{1}{\Gamma(1-\alpha)} \int_0^t g'(s)(t-s)^{-\alpha}\,ds, \qquad 0<\alpha<1,$$6 when the stock price is near or above strike (Singh et al., 2022). A plausible implication is that fractional memory effects manifest most clearly where payoff geometry and boundary interaction already make the pricing surface sensitive.

The main conceptual boundary in the literature concerns what counts as a TFBSE. A paper on fractional stochastic volatility driven by an fOU process with fBm input derives $${}_{0}^{c}D_t^\alpha g(t) = \frac{1}{\Gamma(1-\alpha)} \int_0^t g'(s)(t-s)^{-\alpha}\,ds, \qquad 0<\alpha<1,$$7 corrections to Black–Scholes prices and implied volatilities, including fractional-power maturity effects, but the pricing operator itself remains the classical Black–Scholes PDE plus perturbative corrections (Garnier et al., 2015). It belongs to the broader fractional-finance literature, but not to the TFBSE in the PDE sense. This distinction is necessary to avoid conflating memory in volatility dynamics with memory in the pricing equation.

Another internal boundary separates direct fractional-PDE modeling from probabilistically derived TFBSEs. The ad hoc Caputo replacement

$${}_{0}^{c}D_t^\alpha g(t) = \frac{1}{\Gamma(1-\alpha)} \int_0^t g'(s)(t-s)^{-\alpha}\,ds, \qquad 0<\alpha<1,$$8

is common and computationally tractable (Singh et al., 2022, Garg et al., 30 Jan 2026), but the subdiffusive and generalized inverse-subordinator works provide a more rigorous foundation by deriving the pricing PDE from time-changed processes and martingale pricing (Krzyżanowski et al., 2019, Zhang et al., 13 Nov 2025). This suggests two parallel strands of the literature: phenomenological fractionalization and stochastic-clock derivation.

7. Open issues and research directions

The literature identifies several persistent technical issues. The first is the weak initial singularity of TFBSE solutions. Uniform-step schemes with formal order $${}_{0}^{c}D_t^\alpha g(t) = \frac{1}{\Gamma(1-\alpha)} \int_0^t g'(s)(t-s)^{-\alpha}\,ds, \qquad 0<\alpha<1,$$9 or higher often lose their advertised convergence unless sufficient regularity is assumed. This has driven the development of graded meshes, variable-step Alikhanov formulas, and nonuniform tempered $\alpha$0 methods (Song et al., 2021, Zhou et al., 2023, Dimitrov et al., 2016). The dominance of this issue suggests that temporal regularity, not merely fractional order, is central to the design of reliable solvers.

The second issue is the history cost of Caputo-type derivatives. Direct convolution implies $\alpha$1 temporal work in straightforward implementations. Some papers leave this unaddressed (V et al., 9 Aug 2025, Singh et al., 2022), while others introduce sum-of-exponentials compression or recursive auxiliary variables to reduce cost to essentially linear growth in the number of time steps up to the SOE rank (Song et al., 2021, Zhou et al., 2023). Fast memory evaluation is therefore one of the main algorithmic frontiers of TFBSE computation.

A third issue is the relationship between high spatial order and payoff nonsmoothness. Fourth-order compact and average schemes perform strongly on manufactured smooth solutions (Dimitrov et al., 2016, Song et al., 2021, Zhou et al., 2023), but practical option payoffs contain kinks at the strike. Some papers explicitly note that this can degrade the full formal order and motivate local mesh refinement or piecewise-uniform meshes (Zhou et al., 2023). This suggests that future high-order work must treat payoff regularization and adaptive resolution as first-class concerns.

A fourth issue is financial derivation and arbitrage interpretation. Many PDE papers simply replace the time derivative by a fractional operator, without deriving the modified pricing equation from self-financing arguments, anomalous diffusion, or time-changed martingale models (Zakaria et al., 2020, Bajalan et al., 2021). By contrast, the subdiffusive and generalized sub-diffusive works explicitly derive the TFBSE from time-changed stochastic models and construct appropriate martingale measures (Krzyżanowski et al., 2019, Zhang et al., 13 Nov 2025). A plausible implication is that future TFBSE research may increasingly distinguish between numerically useful fractional surrogates and fully specified no-arbitrage market models.

Finally, multidimensional and nonstandard-contract TFBSEs remain comparatively underdeveloped. Two-asset Caputo PDEs exist in formal series form (Zakaria et al., 2020), and time-space-fractional meshless frameworks can in principle handle higher-dimensional asset spaces (Torres-Hernandez et al., 2020), but rigorous, financially calibrated multidimensional TFBSE pricing remains sparse. The same is true for American-style products, where free boundaries would interact nontrivially with temporal memory.

Taken together, these works position the TFBSE as a mature but still evolving subject at the intersection of anomalous diffusion, inverse-subordinator stochastic modeling, and high-order nonlocal PDE numerics. The core idea is stable across the literature: the pricing operator retains Black–Scholes spatial structure while time evolution is made nonlocal. What varies is the derivative class, the probabilistic foundation, the contract geometry, and the numerical strategy used to resolve long memory efficiently and accurately (Krzyżanowski et al., 2019, Song et al., 2021, Zhou et al., 2023, Zhang et al., 13 Nov 2025).