From Diffusion to Fractional Diffusion: A CTRW Approach to Anomalous Transport

Introduction

Since the earliest tokamak experiments, measured cross-field transport has exceeded both classical and neoclassical predictions, often by one to two orders of magnitude [3, 4]. The standard picture attributes this to turbulence. There is also a mathematical issue. The local diffusive model for a scalar field $\phi$,

$$\partial_t \phi = \partial_x \bigl[\chi\, \partial_x \phi\bigr] + S,$$

follows from a conservation law and Fick’s law. Its derivation from a Brownian random walk requires Gaussian statistics, Markovian dynamics, and a well-defined transport scale [1]. Turbulent trapping and large radial flights violate all three.

We ask what equation replaces ordinary diffusion when the underlying stochastic process has heavy-tailed statistics. One such equation, developed by Montroll and Weiss [5], Metzler and Klafter [1], and others, is a fractional diffusion equation

$${}^{C}_{0}D_t^{\beta} P = \chi\, D_{|x|}^{\alpha} P,$$

where ${}^{C}_{0}D_t^{\beta}$ is the Caputo fractional time derivative of order $0 < \beta \leq 1$, and $D_{|x|}^{\alpha}$ is the Riesz fractional space derivative of order $0 < \alpha \leq 2$. For $\alpha=2$ and $\beta=1$ the equation reduces to the ordinary diffusion equation. The fractional orders are set by the tail exponents of the jump and waiting-time distributions in the underlying CTRW.

An application to plasma turbulence was given by del-Castillo-Negrete, Carreras, and Lynch [2, 7]. They simulated tracer particles in pressure-gradient-driven plasma turbulence and found that the radial displacement pdf had algebraic tails decaying as $|x|^{-(1+\alpha)}$ with $\alpha \approx 3/4$, and that the similarity exponent was $\nu = 2/3$. Equation [eq:fde] with $\alpha = 3/4$ and $\beta = 1/2$ reproduced these results quantitatively.

The thesis is organized as follows. Section 2 derives the fractional diffusion equation from the CTRW, treating the Brownian and heavy-tailed cases side by side. Section 3 defines the fractional operators and shows how they invert the Fourier-Laplace transform of the CTRW solution. Section 4 discusses the physical content of the model and its limitations.

From Random Walks to Fractional Diffusion

The Montroll—Weiss Equation

Consider a particle that waits a random time $\tau_i$ drawn from a pdf $\psi(\tau)$, then jumps by a random displacement $\zeta_i$ drawn from a pdf $\lambda(\zeta)$. We assume waiting times and jump lengths are independent. The probability $P(x,t)$ of finding the particle at position $x$ at time $t$ satisfies the Montroll—Weiss equation [1]:

$$P(x,t) = \delta(x)\int_t^{\infty} \psi(t')\,dt' + \int_0^t \psi(t-t') \left[\int_{-\infty}^{\infty} \lambda(x-x')\,P(x',t')\,dx'\right] dt'.$$

The first term accounts for particles that have not yet jumped. The second accounts for particles that arrived at $x'$ at time $t'$ and then jumped to $x$. Introducing the Fourier transform $\hat{f}(k) = \int e^{ikx} f(x)\,dx$ and the Laplace transform $\tilde{f}(s) = \int_0^{\infty} e^{-st} f(t)\,dt$, equation [eq:MW_real] becomes algebraic [1]:

$$\hat{\tilde{P}}(k,s) = \frac{1 - \tilde{\psi}(s)}{s} \cdot \frac{1}{1 - \tilde{\psi}(s)\,\hat{\lambda}(k)}.$$

Apply the Laplace transform to [eq:MW_real]. The survival probability $\int_t^{\infty}\psi(t')\,dt'$ in the first term has Laplace transform

$$\mathcal{L}\!\left[\int_t^{\infty}\psi(t')\,dt'\right] = \int_0^{\infty}e^{-st}\int_t^{\infty}\psi(t')\,dt'\,dt = \int_0^{\infty}\psi(t')\int_0^{t'}e^{-st}\,dt\,dt' = \frac{1-\tilde{\psi}(s)}{s},$$

where Fubini’s theorem was used to exchange the order of integration. The second term in [eq:MW_real] is a time convolution of $\psi$ with the spatial convolution $\int\lambda(x-x')P(x',t')\,dx'$. By the Laplace convolution theorem it transforms to $\tilde{\psi}(s)$ times the Laplace transform of the spatial term; the Fourier convolution theorem then converts the spatial integral to $\hat{\lambda}(k)\,\hat{\tilde{P}}(k,s)$. Taking the Fourier transform of the first term gives $\hat{\delta}(k)\cdot(1-\tilde{\psi}(s))/s = (1-\tilde{\psi}(s))/s$. Assembling both sides:

$$\hat{\tilde{P}} = \frac{1-\tilde{\psi}(s)}{s} + \tilde{\psi}(s)\,\hat{\lambda}(k)\,\hat{\tilde{P}}.$$

Solving for $\hat{\tilde{P}}$ gives [eq:MW] directly.

Recovering Ordinary Diffusion

Take an exponential waiting-time pdf $\psi_M(\tau) = \mu\, e^{-\mu\tau}$ with finite mean $\langle \tau \rangle = 1/\mu$, and a Gaussian jump pdf $\lambda_F(\zeta) = (2\pi\sigma)^{-1/2}\exp(-\zeta^2/2\sigma)$ with finite variance $\langle \zeta^2 \rangle = \sigma$. To take the continuum limit, we rescale waiting times by a factor $r$ and jumps by a factor $h$, then send $r, h \to 0$. The small-parameter expansions of the transforms are

$$\begin{aligned} \tilde{\psi}_M(rs) &= \frac{1}{1 + rs\langle\tau\rangle} \approx 1 - rs\langle\tau\rangle + \cdots,\\[4pt] \hat{\lambda}_F(hk) &= e^{-\langle\zeta^2\rangle h^2 k^2/2} \approx 1 - \langle\zeta^2\rangle h^2\,\frac{k^2}{2} + \cdots. \end{aligned}$$

Substituting [eq:psi_gauss]—[eq:lam_gauss] into [eq:MW] and expanding to leading order in $r$ and $h$, the numerator becomes $1-\tilde{\psi}(rs)\approx rs\langle\tau\rangle$, while the denominator factor satisfies

$$1 - \tilde{\psi}(rs)\,\hat{\lambda}(hk) \approx 1 - \bigl(1 - rs\langle\tau\rangle\bigr)\!\bigl(1 - \tfrac{1}{2}\langle\zeta^2\rangle h^2 k^2\bigr) \approx rs\langle\tau\rangle + \tfrac{1}{2}\langle\zeta^2\rangle h^2 k^2,$$

where the cross-term $rs\langle\tau\rangle\cdot\frac{1}{2}\langle\zeta^2\rangle h^2k^2$ is of higher order in $r$ and $h$ and is dropped. To see this consistently: the diffusivity $\chi = h^2\langle\zeta^2\rangle/(2r\langle\tau\rangle)$ is held finite as $r,h\to 0$, which requires $h^2 \sim r$. Under this coupling the cross-term scales as $r\cdot h^2 \sim r^2$, while the two retained terms each scale as $r$, so the cross-term is sub-leading and the truncation is self-consistent. Therefore

$$\hat{\tilde{P}} \approx \frac{rs\langle\tau\rangle}{s\bigl(rs\langle\tau\rangle + \frac{1}{2}\langle\zeta^2\rangle h^2 k^2\bigr)} = \frac{1}{s + \dfrac{\langle\zeta^2\rangle h^2}{2r\langle\tau\rangle}k^2} = \frac{1}{s + \chi k^2},$$

where $\chi = h^2\langle\zeta^2\rangle/(2r\langle\tau\rangle)$ is held finite as $r,h\to 0$. Rearranging gives

$$s\,\hat{\tilde{P}} - 1 = -\chi\, k^2\, \hat{\tilde{P}}.$$

Using $P(x,0)=\delta(x)$, the standard transform identities

$$\begin{aligned} \mathcal{L}[\partial_t P] &= s\tilde{P}(s) - \delta(x), \\ \mathcal{F}[\partial_x^2 P] &= -k^2 \hat{P}(k), \end{aligned}$$

invert [eq:diff_FL] to give $\partial_t P = \chi\,\partial_x^2 P$. This holds whenever both $\langle \tau \rangle$ and $\langle \zeta^2 \rangle$ are finite [1].

Heavy Tails and the Breakdown of the CLT

Now suppose the dynamics produces power-law tails [2]:

$$\psi(\tau) \sim \tau^{-(1+\beta)}, \quad 0 < \beta < 1, \qquad \lambda(\zeta) \sim |\zeta|^{-(1+\alpha)}, \quad 0 < \alpha < 2.$$

The mean waiting time $\langle \tau \rangle = \int \tau\,\psi(\tau)\,d\tau$ diverges because $\beta < 1$. The jump variance $\langle \zeta^2 \rangle = \int \zeta^2\,\lambda(\zeta)\,d\zeta$ diverges because $\alpha < 2$. More generally, the moment $\langle |\zeta|^n \rangle$ diverges for any $n \geq \alpha$ [1]. Neither $\langle \tau \rangle$ nor $\langle \zeta^2 \rangle$ exists, so the standard CLT does not apply and the continuum limit cannot produce the ordinary diffusion equation.

What happens instead is governed by a Tauberian theorem: a pdf with algebraic tail $\psi(\tau) \sim \tau^{-(1+\beta)}$ has a Laplace transform that behaves, for small $s$, as $\tilde{\psi}(s) \approx 1 - c_1 s^{\beta}$. The singularity structure of the transform at small argument encodes the slow decay of the pdf at large argument. Similarly, $\hat{\lambda}(k) \approx 1 - c_2 |k|^{\alpha}$ for small $k$ [1]. Rescaling as before and substituting into [eq:MW]:

$$\begin{aligned} \tilde{\psi}(rs) &\approx 1 - c_1(rs)^{\beta} + \cdots,\\[4pt] \hat{\lambda}(hk) &\approx 1 - c_2(h|k|)^{\alpha} + \cdots. \end{aligned}$$

Substituting [eq:psi_heavy]—[eq:lam_heavy] into [eq:MW], the numerator is $1-\tilde{\psi}(rs)\approx c_1 r^{\beta}s^{\beta}$, and the denominator factor becomes

$$1 - \tilde{\psi}(rs)\,\hat{\lambda}(hk) \approx c_1 r^{\beta}s^{\beta} + c_2 h^{\alpha}|k|^{\alpha},$$

again dropping the higher-order cross-term. Assembling:

$$\hat{\tilde{P}} \approx \frac{c_1 r^{\beta} s^{\beta}}{s\bigl(c_1 r^{\beta} s^{\beta} + c_2 h^{\alpha}|k|^{\alpha}\bigr)} = \frac{s^{\beta-1}}{s^{\beta} + \chi|k|^{\alpha}},$$

where $\chi = c_2 h^{\alpha}/(c_1 r^{\beta})$ is held finite. To obtain a nontrivial continuum limit, the two denominator terms must remain of the same order as $r,h\to 0$. This scaling keeps the spatial and temporal terms at the same order in the continuum limit. Rearranging gives

$$\boxed{s^{\beta}\,\hat{\tilde{P}} - s^{\beta - 1} = -\chi\,|k|^{\alpha}\,\hat{\tilde{P}}.}$$

Compare this with [eq:diff_FL]: $s$ has been replaced by $s^{\beta}$, and $k^2$ by $|k|^{\alpha}$. The tail exponents of the underlying random walk now appear as the orders of the operators in transform space. The next question is which real-space operators correspond to $s^{\beta}$ and $|k|^{\alpha}$.

Fractional Operators and the Fractional Diffusion Equation

The Riemann—Liouville Fractional Integral and Derivative

The fractional integral is a generalization of the Cauchy formula for repeated integration. For integer $n$, the $n$-fold integral of $\phi$ from $a$ to $x$ can be written as a single convolution integral [1]:

$${}_a D_x^{-n}\,\phi = \frac{1}{(n-1)!}\int_a^x (x-y)^{n-1}\,\phi(y)\,dy.$$

Replacing $(n-1)!$ with $\Gamma(\nu)$ extends this to non-integer order $\nu > 0$:

$${}_a D_x^{-\nu}\,\phi = \frac{1}{\Gamma(\nu)}\int_a^x (x-y)^{\nu-1}\,\phi(y)\,dy.$$

The Riemann—Liouville fractional derivative of order $\alpha$ is then defined as an integer differentiation of a fractional integral [2]:

$${}_a D_x^{\alpha}\,\phi = \frac{d^m}{dx^m}\bigl[{}_a D_x^{-(m-\alpha)}\phi\bigr], \qquad m = \lceil \alpha \rceil,$$

where $m$ is the smallest integer greater than $\alpha$. For $1 < \alpha \leq 2$, we have $m=2$ and the derivative takes the explicit form

$${}_a D_x^{\alpha}\,\phi = \frac{1}{\Gamma(2-\alpha)}\,\partial_x^2 \int_a^x \frac{\phi(y)}{(x-y)^{\alpha-1}}\,dy.$$

For $0 < \alpha \leq 1$, one has $m=1$ and the second derivative above is replaced by a first derivative. In either case the operator is non-local: it depends on the values of $\phi$ over the entire interval $(a,x)$, not just at $x$.

A key property for Fourier analysis: on the semi-infinite domain $(-\infty, x)$, the Weyl fractional derivative satisfies ${}_{-\infty}D_x^{\mu}\,e^{ikx} = (ik)^{\mu}\,e^{ikx}$ [1].

Left, Right, and Symmetric Derivatives

The RL derivative [eq:RL_deriv] integrates to the left of $x$. A right derivative, integrating over $(x,b)$, is defined analogously [2]:

$${}_x D_b^{\alpha}\,\phi = \frac{(-1)^m}{\Gamma(m-\alpha)}\,\partial_x^m \int_x^b \frac{\phi(y)}{(y-x)^{\alpha+1-m}}\,dy.$$

The general spatial fractional operator is a weighted sum $D_x^{\alpha} = w^{-}\,{}_a D_x^{\alpha} + w^{+}\,{}_x D_b^{\alpha}$, where $w^{\pm}$ control left-right asymmetry. The symmetric Riesz derivative on $(-\infty,\infty)$ is defined as [2]:

$$D_{|x|}^{\alpha} = \frac{-1}{2\cos(\pi\alpha/2)}\bigl[{}_{-\infty}D_x^{\alpha} + {}_x D_{\infty}^{\alpha}\bigr],$$

and its Fourier transform, for $0 < \alpha < 2$, gives

$$\mathcal{F}\bigl[D_{|x|}^{\alpha}\,\phi\bigr] = -|k|^{\alpha}\,\hat{\phi}(k,t).$$

To derive [eq:riesz_FT], apply the Fourier transform directly using the Weyl eigenvalue relations ${}_{-\infty}D_x^{\alpha}e^{ikx}=(ik)^{\alpha}e^{ikx}$ and ${}_xD_{\infty}^{\alpha}e^{ikx}=(-ik)^{\alpha}e^{ikx}$ [1]. For any real $k$, writing $ik=|k|e^{i(\pi/2)\operatorname{sgn}(k)}$ and $-ik=|k|e^{-i(\pi/2)\operatorname{sgn}(k)}$, one has

$$(ik)^{\alpha}+(-ik)^{\alpha} = |k|^{\alpha}\!\left(e^{i\pi\alpha\operatorname{sgn}(k)/2}+e^{-i\pi\alpha\operatorname{sgn}(k)/2}\right) = 2|k|^{\alpha}\cos\!\left(\frac{\pi\alpha}{2}\right).$$

Applying the Fourier transform to the Riesz operator therefore gives

$$\mathcal{F}\bigl[D_{|x|}^{\alpha}\phi\bigr] = \frac{-1}{2\cos(\pi\alpha/2)}\bigl[(ik)^{\alpha}+(-ik)^{\alpha}\bigr]\hat{\phi}(k) = \frac{-1}{2\cos(\pi\alpha/2)}\cdot 2|k|^{\alpha}\cos\!\left(\frac{\pi\alpha}{2}\right)\hat{\phi}(k) = -|k|^{\alpha}\hat{\phi}(k),$$

confirming [eq:riesz_FT]. The prefactor $-1/(2\cos(\pi\alpha/2))$ in the definition of the Riesz operator is the constant needed to cancel the cosine and yield the Fourier symbol $-|k|^{\alpha}$. This holds for the full range $0 < \alpha < 2$, which is important because the plasma case study of del-Castillo-Negrete et al. uses $\alpha = 3/4 < 1$.

The Caputo Time Derivative

For the time direction, the RL time derivative has the Laplace transform $\mathcal{L}$, which involves the initial value of a fractional integral of $\phi$ rather than $\phi$ itself. This is not suitable for initial-value problems. The Caputo derivative resolves this by switching the order of differentiation and integration [2]:

$${}^{C}_{0}D_t^{\beta}\,\phi = \frac{1}{\Gamma(1-\beta)}\int_0^t \frac{\partial_{\tau}\phi(x,\tau)}{(t-\tau)^{\beta}}\,d\tau, \qquad 0 < \beta < 1.$$

Its Laplace transform is

$$\mathcal{L}\bigl[{}^{C}_{0}D_t^{\beta}\,\phi\bigr] = s^{\beta}\,\tilde{\phi}(s) - s^{\beta-1}\,\phi(0),$$

To derive [eq:caputo_LT], write the Caputo derivative [eq:caputo] as the time convolution ${}^{C}_{0}D_t^{\beta}\phi = ({\partial_t\phi})*g$ where $g(t)=t^{-\beta}/\Gamma(1-\beta)$. By the Laplace convolution theorem,

$$\mathcal{L}\bigl[{}^{C}_{0}D_t^{\beta}\phi\bigr] = \mathcal{L}[\partial_t\phi]\cdot\mathcal{L}[g].$$

The first factor is $\mathcal{L}[\partial_t\phi]=s\tilde{\phi}(s)-\phi(0)$. For the second, the substitution $u=st$ gives

$$\mathcal{L}[t^{-\beta}] = \int_0^{\infty}e^{-st}t^{-\beta}\,dt = s^{\beta-1}\int_0^{\infty}e^{-u}u^{-\beta}\,du = \Gamma(1-\beta)\,s^{\beta-1},$$

so $\mathcal{L}[g]=s^{\beta-1}$. Combining and expanding:

$$\mathcal{L}\bigl[{}^{C}_{0}D_t^{\beta}\phi\bigr] = s^{\beta-1}\bigl(s\tilde{\phi}(s)-\phi(0)\bigr) = s^{\beta}\tilde{\phi}(s) - s^{\beta-1}\phi(0),$$

which is [eq:caputo_LT], and depends directly on $\phi(0)$. The Caputo derivative of a constant is zero, unlike the RL derivative. This is necessary for well-defined steady-state solutions.

Assembling the Equation

We had from the CTRW continuum limit, equation [eq:fde_FL]:

$$s^{\beta}\,\hat{\tilde{P}} - s^{\beta-1} = -\chi\,|k|^{\alpha}\,\hat{\tilde{P}}.$$

From [eq:caputo_LT], the left side is the Fourier-Laplace transform of ${}^{C}_{0}D_t^{\beta} P$ with $P(x,0)=\delta(x)$. From [eq:riesz_FT], the right side is $\chi$ times the transform of $D_{|x|}^{\alpha} P$. Inverting:

$$\boxed{{}^{C}_{0}D_t^{\beta}\,P = \chi\,D_{|x|}^{\alpha}\,P.}$$

This is the fractional diffusion equation. At $\alpha=2$, $\beta=1$ the Caputo derivative reduces to $\partial_t$ and the Riesz derivative to $\partial_x^2$, recovering ordinary diffusion.

Properties of the Solution

The Fourier—Laplace transform of [eq:fde_final], with initial condition $G(x,0)=\delta(x)$ so that $\hat{G}(k,0)=1$, is precisely [eq:fde_FL]. Solving algebraically for $\hat{\tilde{G}}$:

$$\hat{\tilde{G}}(k,s) = \frac{s^{\beta-1}}{s^{\beta}+\chi|k|^{\alpha}}.$$

Inverting the Laplace transform uses the Mittag-Leffler identity $\mathcal{L}[E_{\beta}(-ct^{\beta})]=s^{\beta-1}/(s^{\beta}+c)$, which follows from the series $E_{\beta}(-ct^{\beta})=\sum_{n=0}^{\infty}(-c)^n t^{\beta n}/\Gamma(\beta n+1)$ and the standard transform $\mathcal{L}[t^{\beta n}]=\Gamma(\beta n+1)/s^{\beta n+1}$:

$$\mathcal{L}[E_{\beta}(-ct^{\beta})] = \sum_{n=0}^{\infty}\frac{(-c)^n}{\Gamma(\beta n+1)}\cdot\frac{\Gamma(\beta n+1)}{s^{\beta n+1}} = \frac{1}{s}\sum_{n=0}^{\infty}\!\left(\!\frac{-c}{s^{\beta}}\right)^{\!n} = \frac{s^{\beta-1}}{s^{\beta}+c}.$$

Setting $c=\chi|k|^{\alpha}$ in [eq:green_FL] therefore gives

$$\hat{G}(k,t) = E_{\beta}\!\left(-\chi|k|^{\alpha}t^{\beta}\right).$$

For the parameter choice used here, $\alpha>\beta$, and the corresponding Green’s function is a proper probability density [2]; in the plasma case, $\alpha=3/4>\beta=1/2$. The self-similar form of the Green’s function of [eq:fde_final] now follows. To see the scaling explicitly, suppose $G(x,t) = t^{-\nu} K(x/t^{\nu})$ for some exponent $\nu$ and profile $K$. Substituting into the Fourier transform and changing variables $\eta = x/t^{\nu}$ gives $\hat{G}(k,t) = \hat{K}(kt^{\nu})$. For this to equal $E_{\beta}(-\chi|k|^{\alpha}t^{\beta})$ from [eq:green_Fourier], the argument $kt^{\nu}$ must absorb all $t$-dependence; this requires $|kt^{\nu}|^{\alpha} \propto |k|^{\alpha}t^{\beta}$, which forces $\nu\alpha = \beta$, i.e. $\nu = \beta/\alpha$. The full scaling factor follows from normalisation $\int G\,dx = 1$. Since $\hat{G}$ depends on $k$ and $t$ only through the combination $\chi|k|^{\alpha}t^{\beta} = |k(\chi^{1/\beta}t)^{\beta/\alpha}|^{\alpha}$, the inverse Fourier transform must satisfy $G(x,t)=(\chi^{1/\beta}t)^{-\beta/\alpha}K(\eta)$ with similarity variable

$$G(x,t) = (\chi^{1/\beta}\,t)^{-\beta/\alpha}\,K(\eta), \qquad \eta = x\,(\chi^{1/\beta}\,t)^{-\beta/\alpha},$$

where the profile $K$ is the inverse Fourier transform of $\hat{K}(q)=E_{\beta}(-|q|^{\alpha})$, i.e., $K(\eta)=(1/\pi)\int_0^{\infty}\cos(\eta z)\,E_{\beta}(-z^{\alpha})\,dz$ [2]. Two asymptotic regimes characterize $K$: near the origin, $K(\eta) \sim A/\eta^{1-\alpha} + B$ (a weak cusp for $\alpha < 1$), and for large $|\eta|$,

$$K(\eta) \sim \frac{C}{\eta^{1+\alpha}}, \qquad C = \frac{1}{\pi}\frac{\Gamma(1+\alpha)}{\Gamma(1+\beta)}\sin\frac{\pi\alpha}{2}.$$

The large-$|\eta|$ tail of $K$ follows from the small-$|q|$ behaviour of $\hat{K}$. Using the power-series definition

$$E_{\beta}(-z)=\sum_{n=0}^{\infty}\frac{(-z)^n}{\Gamma(\beta n+1)},$$

we obtain, for small $z$,

$$E_{\beta}(-z)=1-\frac{z}{\Gamma(1+\beta)}+O(z^2).$$

Substituting $z=|q|^\alpha$ gives

$$\hat{K}(q)=E_{\beta}(-|q|^\alpha) =1-\frac{|q|^\alpha}{\Gamma(1+\beta)}+O(|q|^{2\alpha}),$$

so the leading non-analytic term is $-|q|^\alpha/\Gamma(1+\beta)$. The Tauberian theorem for Fourier transforms [1] states that a function whose Fourier transform behaves as $-A|q|^{\alpha}$ near $q=0$ has algebraic tails $A\Gamma(1+\alpha)\sin(\pi\alpha/2)/(\pi|\eta|^{1+\alpha})$; setting $A=1/\Gamma(1+\beta)$ gives $C=\Gamma(1+\alpha)\sin(\pi\alpha/2)/(\pi\Gamma(1+\beta))$, confirming [eq:tail_decay].

An important caveat concerns the moments. Formally, the $n$-th moment of the self-similar propagator scales as

$$\langle |x|^n \rangle \sim t^{n\beta/\alpha},$$

so the similarity exponent $\nu = \beta/\alpha$ determines the transport regime: $\nu = 1/2$ for ordinary diffusion, $\nu > 1/2$ for superdiffusion, $\nu < 1/2$ for subdiffusion. However, because $K(\eta) \sim \eta^{-(1+\alpha)}$ for large $\eta$, the integrand $\eta^n K(\eta) \sim \eta^{n-1-\alpha}$ diverges for $n \geq \alpha$ when integrated over $(-\infty, \infty)$ [2]. In particular, for $\alpha < 2$ the mean squared displacement is not defined on the infinite line. Del-Castillo-Negrete et al. handle this by introducing a finite spatial cutoff, which is physically natural since real systems always have a bounded domain [2]. The moments that do remain well-defined are the fractional moments $\langle |x|^{\delta} \rangle$ for $\delta < \alpha$, which are finite and scale as $t^{\delta\beta/\alpha}$.

Discussion and Conclusion

The Plasma Case Study

Del-Castillo-Negrete et al. [2] simulated 25,000 tracer particles in 3D resistive pressure-gradient-driven turbulence. The turbulence produces $E\times B$ eddies that trap particles, and avalanche-like instabilities that launch them on long radial flights. Both effects are visible in the tracer orbits (their Figure 2). The measured pdf of radial displacements has tails decaying as $|x|^{-1.75}$, giving $\alpha \approx 3/4$. The similarity exponent, measured from the collapse of the pdf under rescaling by $t^{\nu}$, is $\nu \approx 2/3$, which with $\nu = \beta/\alpha$ requires $\beta = 1/2$. With $\chi = 0.09$, the fractional model [eq:fde_final] reproduces both the pdf shape and the self-similar collapse. Note that because $\alpha = 3/4 < 1$, the mean squared displacement of the fractional Green’s function is not finite on the full line, and the moment scaling reported in [2] relies on the finite spatial extent of the simulation domain acting as a natural cutoff.

This case is interesting because $\alpha < 1$ places it outside the regime where the RL derivative takes the more familiar $m=2$ form [eq:RL_explicit]. The Riesz derivative is still well-defined for this value through its Fourier transform property [eq:riesz_FT], and the derivation from the CTRW in Section 2 does not require $\alpha > 1$.

The subdiffusion case $\alpha = 2$, $\beta < 1$ visible in Figure 5 (right column) is worth isolating, as it clarifies the role of the temporal fractional order in the absence of any spatial anomaly. Setting $\alpha = 2$ in [eq:fde_final] gives ${}^{C}_{0}D_t^{\beta} P = \chi\,\partial_x^2 P$: the spatial operator is the ordinary Laplacian, but the Caputo time derivative introduces non-Markovian memory through the convolution in [eq:caputo]. The Green’s function [eq:green] reduces to $G(x,t) = t^{-\beta/2} K(\eta)$ with $\eta = x/t^{\beta/2}$, and the MSD grows as $\langle x^2\rangle \sim t^{\beta} < t$ for $\beta < 1$, that is, slower than ordinary diffusion. The connection to Figure 4 is direct: the Fourier transform $\hat{G}(k,t) = E_{\beta}(-\chi k^2 t^{\beta})$ decays as a Mittag-Leffler function in time at fixed $k$. At large $t$, the asymptotic $E_{\beta}(-\chi k^2 t^{\beta}) \sim (\chi k^2 t^{\beta})^{-1}/\Gamma(1-\beta)$ corresponds to the algebraic relaxation visible in Figure 4 for $\beta < 1$. The slowdown arises entirely from the divergent mean waiting time $\langle\tau\rangle$: the spatial jump distribution remains Gaussian, and only the waiting-time distribution is anomalous.

Connection to Experimental Diagnostics

Transfer-entropy studies have also been discussed in CTRW language. Van Milligen et al. [8] describe “slow” and “fast” transport channels, with the latter consistent with jump-like propagation across minor transport barriers. This is qualitatively compatible with a heavy-tailed jump distribution, but it is not a derivation of the fractional model.

The NLCC method of Ding et al. [9] addresses a different question: directional asymmetry in nonlinear coupling between radially separated probe signals. Such asymmetry may be compared with an imbalance between left and right spatial fractional derivatives, but I did not find such a connection in the literature. The NLCC construction itself is based on phase-space reconstruction rather than fractional calculus.

They are discussed further in the author’s physics thesis.

Limitations of the Fractional Approximation

The fractional diffusion equation is an asymptotic result: it describes the large-scale, long-time behaviour of the CTRW, and can converge slowly at finite times or near the origin. Barkai [6] showed that the exact CTRW propagator at $x = 0$ may differ from the fractional Green’s function over a long transient, because the fractional equation captures only the leading term of the asymptotic expansions [eq:psi_heavy]—[eq:lam_heavy]. Higher-order corrections in $s$ and $k$ are discarded in taking the continuum limit, and these corrections matter most near the origin and at early times.

Separately, Rodriguez-Fernandez et al. [10] showed that cold-pulse phenomena in tokamaks, which have been cited as evidence of nonlocal transport, can be reproduced by local quasilinear turbulent transport models (specifically the TGLF-SAT1 saturation rule). This does not invalidate the fractional framework, but it does mean that apparently nonlocal transport signatures in experiments are not by themselves enough to establish genuinely nonlocal transport. The fractional model is a reduced model: it compactly captures certain statistical features of anomalous transport, not the underlying first-principles dynamics.

Conclusion

We derived the fractional diffusion equation from the CTRW by showing that heavy-tailed jump and waiting-time distributions produce, in the continuum limit, transform-space operators $|k|^{\alpha}$ and $s^{\beta}$ that correspond to the Riesz and Caputo fractional derivatives respectively. The fractional orders are not free parameters: they are fixed by the tail exponents of the underlying stochastic process. The plasma-turbulence results of del-Castillo-Negrete et al. provide a case where this framework is tested against simulation data, and the CTRW simulations presented here reproduce the main qualitative signatures: algebraic tails, self-similar collapse, and, for moments of order below $\alpha$, anomalous scaling. The natural extension is to bounded domains, where the truncated RL derivatives become singular at the boundaries and require Caputo-style regularization [11]. That problem is directly relevant to radial transport modelling in confined plasmas.