Fractional Modeling in Action:
A Survey of Nonlocal Models for Subsurface Transport, Turbulent Flows, and Anomalous Materials

The connection between Brownian motion $B_{t}$ and the classical diffusion equation was studied in seminal works by Bachelier (1900), Einstein (1905), and Von Smoluchowski (1906). The diffusion equation is posed in an initial value problem,

in which $k^{2}>0$ is the diffusion coefficient and $\Delta u=\partial^{2}u/\partial x^{2}$ denotes the Laplacian. Brownian motion $B_{t}$ is a continuous-time stochastic process defined for $t\geq 0$, which when discretized in time steps of size $\Delta t$ has the property that $B_{t=0}=0$ and

The above notation indicates that the increment $\Delta B$ at each time step of $\Delta t$ is drawn from a normal distribution $\mathcal{N}(0,k\sqrt{\Delta t})$ with mean $\mu=0$ and standard deviation $\sigma=k\sqrt{\Delta t}$. This has probability density function

The rule (2) for sampling a path of $B_{t}$ at times $m\Delta t$, $m=0,1,2,...$ is an example of a discrete stochastic differential equation (SDE), and is referred to as the Euler-Maruyama discretization²²2Another frequently used discrete random walk that leads to Brownian motion simply involves steps of fixed length to the left or right with probability $1/2$ each; see Lawler (2010). In the long-time limit, all such discrete walks that draw increments from a finite-variance distribution lead to Brownian motion, due to the central limit theorem Zaburdaev et al. (2015) of Brownian motion (Kloeden and Platen, 1992).

The discrete process $B_{m\Delta t}$ should be thought of as tracing a path in $\mathbb{R}$ of a particle undergoing “jumps” in a random direction at time intervals of size $\Delta t$. At each time $t$, the position $B_{t}$ of the particle is a random variable. It can be shown that the paths of the continuous-time process $B_{t}$ are almost surely continuous in time Revuz and Yor (2013). From (2) and the central limit theorem, it follows that Brownian motion satisfies the scaling property

where the left-hand side denotes the variance or second moment of the random variable $B_{t}$; see Figure 1. Given an initial distribution of particles $u_{0}(x)$ in $\mathbb{R}$ which then undergo Brownian motion, the distribution $u(x,t)$ of particles in $\mathbb{R}$ is governed by (1). In other words, diffusing particles described at a microscopic scale by Brownian motion, i.e. by (2) in discrete time, have their distribution in space – a macroscopic property – governed by the heat equation (Meerschaert and Sikorskii, 2011, Section 1.1). This is illustrated in Figure 2.

The consistency between this macroscopic description and the microscopic model is illustrated by scaling properties. A necessary property of the Brownian motion model is the second-moment condition (4), which states that on average, particles travel a distance $k\sqrt{t}$ from their initial position after time $t$. This is reflected in the fact that the solution of (1) with initial condition $u_{0}(x)=\delta_{0}(x)$ is

which is the normal density (3) with standard deviation $k\sqrt{t}$. Note that this solution has the property that

Thus, the distribution of plume of particles in this diffusion model spreads out as $(t_{2}/t_{1})^{1/2}$ as time elapses from $t_{1}$ to $t_{2}$, consistent with (4). This scaling of the fundamental solution is illustrated in Figure 1.

The model for normal diffusion reviewed here is also referred to as Fickian diffusion. The heat equation (1) can be derived from the mass conservation with flux term $J$,

under Fick’s law $J=\nabla u$. As discussed by Schumer et al. (2001), the fractional diffusion equations we introduce below follow from mass conservation with non-Fickian fluxes.

Many important systems exhibit diffusive behavior, but do not satisfy the scaling property (4) (Klafter and Sokolov, 2005). This type of diffusion is referred to as anomalous diffusion, as it cannot be described by (2) with normally distributed increments. We desire a microscopic model that generalizes Brownian motion $B_{t}$, and a corresponding macroscopic model that generalizes the diffusion equation (1). The first model we propose remains in the framework of a discrete SDE with independent identically distributed (i.i.d.) increments,

but the increments $\Delta X$ are no longer drawn from a normal distribution. It follows from the central limit theorem that the only way to obtain a microscopic model in this framework that is statistically distinct from $B_{t}$, i.e., not equivalent in distribution, is to draw step sizes from a probability density function with infinite variance (Zaburdaev et al., 2015; Meerschaert and Sikorskii, 2011).

We introduce the isotropic $\alpha$-stable random variable $S_{\alpha}(\gamma,\sigma,\mu)$. This family of random variables is defined³³3Several parametrizations of the $\alpha$-stable characteristic function exist. The parametrization (10) is due to Samorodnitsky and Taqqu (1994). See Nolan (2020, 1998) for discussions of alternate forms. most simply by their characteristic function. For a general random variable $X$, the characteristic function $\varphi_{X}$ is the Fourier transform of the probability density function $p_{X}$ of $X$, i.e.,

Thus, the characteristic function of the normal random variable $\mathcal{N}(\mu,\sigma)$ is $e^{i\xi\mu-\sigma^{2}\xi^{2}/2}$. Generalizing this, the $\alpha$-stable random variable $S_{\alpha}(\gamma,\sigma,\mu)$ has characteristic function

where

The parameter $\alpha\in(0,2]$ is referred to as the stability parameter of the distribution, $\mu\in\mathbb{R}$ as the center, $\gamma\in[-1,1]$ as the skewness, and $\sigma\in(0,\infty)$ as the scale. The isotropic or symmetric $\alpha$-stable distribution $S_{\alpha}(\gamma=0,\sigma,\mu=0)$ therefore has characteristic function

generalizing the characteristic function of the normal distribution with mean $\mu=0$ and standard deviation $\sigma/\sqrt{2}$ and reducing to it when $\alpha=2$. By definition, the probability density function of $S_{\alpha}(\gamma,\sigma,\mu)$ can be written

In general, the $\alpha$-stable density does not admit a closed-form expression⁴⁴4Special cases are $\alpha=2$ corresponding to the normal distribution, $\alpha=1$ and $\gamma=0$ corresponding to the Cauchy distribution, and $\alpha=1/2$ and $\gamma=1$ corresponding to the Lévy distribution., but in the symmetric case where $\gamma=\mu=0$, it has the property that

as discussed in, e.g., Nolan (2020) or Cont and Tankov (2003). In other words, the density exhibits Paretian or power-law tails. This is in contrast to the rapidly decaying square-exponential tails of the normal distribution. In many settings, such tails are informally referred to as being examples of heavy or fat tails (Adler et al., 1998; Haas and Pigorsch, 2009).

Using the isotropic distribution introduced above, we introduce the isotropic $\alpha$-stable Lévy flight $X^{\alpha}_{t}$ by providing the corresponding discrete stochastic process. This is given for $t=k\Delta t$ with integer $k$ by $X^{\alpha}_{0}=0$ and the rule (Meerschaert and Sikorskii, 2011)

The continuous-time stochastic process $X^{\alpha}_{t}$ for $t\geq 0$ can be thought of as a scaling limit as $\Delta t\rightarrow 0$ of the above random walk, and enjoys several theoretical properties such as stability and an extended central limit theorem (Meerschaert and Sikorskii, 2011; Meerschaert and Scheffler, 2001). However, has the property that for $\alpha<2$, the paths of $X^{\alpha}_{t}$ are almost surely discontinuous, in contrast to Brownian motion – hence the name Lévy “flight”. Given an initial distribution $u_{0}(x)$ of particles in $\mathbb{R}$ which undergo $\alpha$-stable Lévy flight, the evolution of the distribution $u(x,t)$ for $t>0$ is governed by the space-fractional diffusion equation (Meerschaert and Sikorskii, 2011, Section 1.2)

as illustrated in Figure 3. The fractional negative Laplacian $(-\Delta)^{\alpha/2}$ is defined for $0<\alpha<2$ and for any dimension $d$ as

with

see Lischke et al. (2020). We have defined this operator in any dimension for future reference, although our present discussion only requires the case $d=1$. Perhaps the simplest characterization of the fractional Laplacian is the Fourier representation,

The simplest case of (16) is the initial condition $u_{0}(x)=\delta_{0}(x)$, in which case the solution is

This is known as the fundamental solution. Although this solution cannot be written in closed form, it satisfies

as shown in Meerschaert and Sikorskii (2011, Section 1.2). This illustrates that a plume of particles undergoing isotropic $\alpha$-stable Lévy flight spreads by a factor of $(t_{2}/t_{1})^{1/\alpha}$ as time elapses from $t_{1}$ to $t_{2}$, a faster rate when $\alpha<2$ than the normal rate $t^{1/2}$ . Thus, $\alpha$-stable Lévy flight is an example of superdiffusion. The dependence of the above solution as well as sample paths on $\alpha$ is shown in Figure 4.

However, since $\alpha>0$, the tail behavior of the isotropic $\alpha$-stable density implies that the second moment of $X^{\alpha}_{t}$ diverges for $\alpha<2$,

with the first moment (the mean) diverging also when $\alpha\leq 1$ Nolan (2020); Zaburdaev et al. (2015). This implies that the variance of $\alpha$-stable motion is not a useful statistic for parameterizing $\alpha$-stable Lévy flight; it bears no useful relationship to $\alpha$. This aspect can be tackled in several ways, motivating the introduction of further fractional-order operators, such as tempered operators and fractional material derivatives discussed below.

We point out several important properties of the fractional Laplacian. From the definition (17), it is clear that $(-\Delta)^{\alpha/2}c=0$, $c$ being a constant. The fractional Laplacian also satisfies the semigroup property $(-\Delta)^{\alpha/2}(-\Delta)^{\beta/2}=(-\Delta)^{(\alpha+\beta)/2}$ (Samko et al., 1993). However, one property that is apparent from the definition is that, unlike integer-order derivatives, the fractional Laplacian is a nonlocal operator, i.e. the value of $(-\Delta)^{\alpha/2}u({{\color[rgb]{0,0,0}\mathbf{x}}})$ depends on the values of $u$ in all of $\mathbb{R}$ (or $\mathbb{R}^{d}$, for $d>1$). In contrast, the value of any integer-order derivative of $u$ at ${{\color[rgb]{0,0,0}\mathbf{x}}}$ depends only on the values of $u$ in an infinitesimal neighborhood of ${{\color[rgb]{0,0,0}\mathbf{x}}}$.

The fractional Laplacian (17) was introduced in the previous section as a symmetric or rotation invariant operator for describing the symmetric or isotropic $\alpha$-stable Lévy flight. That model introduced a stability parameter $0<\alpha\leq 2$ allowing it to generalize normal diffusion, with the scale $\sigma$ and center $\mu$ playing similar roles as the standard deviation and mean of the normal distribution. However, the stable distribution also allows for a skewness parameter $\gamma\in[-1,1]$, with $\beta=0$ in the symmetric case, which has no analogue in the normal distribution or for Brownian motion. This is due to the central limit theorem, which states that the use of any finite-variance distribution for the i.i.d. increments $\Delta X$ in (8), no matter how asymmetric, leads to $X_{t}$ being normally distributed, and therefore necessarily symmetric about the mean. In this section, we introduce the one-sided Riemann-Liouville fractional derivatives as appropriate operators for modeling asymmetric $\alpha$-stably Lévy flights, which are defined by (15) with $\Delta X^{\alpha}\sim S_{\alpha}(\gamma,\sigma=k(\Delta t)^{1/\alpha},\mu=0)$ for nonzero $\beta$.

The left-sided and right-sided Riemann-Liouville derivatives in $\mathbb{R}$ are defined, for $n=\lceil\alpha\rceil$, as

The texts of Podlubny (1998), Oldham and Spanier (1974), and Meerschaert and Sikorskii (2011) discuss these operators in detail. These derivatives are frequently used in models with $a=-\infty$ and $b=\infty$. In connection with initial value problems, the left-sided Riemann-Liouville derivative in time, $\phantom{}{}^{\phantom{}\text{RL}}_{\phantom{R}0}\mathbb{D}^{\alpha}_{t}u(t)$, is sometimes used with $a=0$. We have written the definitions (23) and (24) to avoid ambiguities in notation, and clearly show that substitution of the variable $x$ occurs after integration and differentiation. An alternative approach is to define Riemann-Liouville fractional integrals separately, as in the right-hand sides of (23) and (24); see Samko et al. (1993).

One quirk of the notation for Riemann-Liouville derivatives in (23) and (24) is the writing of the upper and lower limits of integration $[a,x]$ and $[b,x]$, respectively, as subscripts. While this is suggestive, the result is that the variable of evaluation $x$ occurs twice in the notation for each operator. If these derivatives are evaluated at any numerical value of $x$, this value should be substituted in both locations; thus, $\phantom{}{}^{\phantom{}\text{RL}}_{{a}}\mathbb{D}^{\alpha}_{5}u(5)$ represents a valid evaluation of the derivatives, but $\phantom{}{}^{\phantom{}\text{RL}}_{{a}}\mathbb{D}^{\alpha}_{x}u(5)$ and $\phantom{}{}^{\phantom{}\text{RL}}_{{a}}\mathbb{D}^{\alpha}_{5}u(x)$ do not.

With $a=-\infty$ and $b=\infty$, the Riemann-Liouville derivatives can be represented in frequency space by

In one dimension, these can be used in the asymmetric diffusion model

which describes anomalous diffusion of independent particles. Here, the positions of each particle at time steps of $k\Delta t$ for integer $k$ are governed by (8) with increments $\Delta X$ being drawn from the asymmetric $\alpha$-stable distribution

Thus, the skewness ranges from $\gamma=-1$ when $p=0$ to $\gamma=1$ when $p=1$. The fundamental solution of (27) is

cf. equation (20).

Sample paths of the process just described are illustrated in Figure 6. Note that when $p=1/2$, the distribution reverts to the symmetric $\alpha$-stable distribution, and it can be shown in this case that equation (27) reduces to (16); more specifically,

The Fourier representation (25) suggests that the left-sided Riemann-Liouville derivative ${\phantom{}{}^{\phantom{,}\text{RL}}_{{-\infty}}\mathbb{D}^{\alpha}_{x}u}$ should be thought of as a fractional power of the operator $\partial/\partial x$. However, the correspondence between (27) and (28) makes it clear that to obtain a complete description of $\alpha$-stable Lévy flights in one dimension necessitates two operators, a left-sided and a right-sided operator, which agree with one another when $\alpha=2$. Our interest is these models lies in the fact that an extended centralized limit theorems hold for processes with i.i.d. increments drawn from distributions with infinite variance, but for which the tails of the density function satisfy Pareto-type conditions as in (14). For such processes, $\alpha$-stable distributions play an analogous role to the normal distribution in the classical central limit theorem; unlike the classical theorem, for full generality, skewed $\alpha$-stable distributions must be included in such a result. See Meerschaert and Scheffler (2001) or Meerschaert and Sikorskii (2011) for a treatment of these results.

We mention how the Riemann-Liouville derivative can be utilized in dimensions $d>1$. An anisotropic diffusion operator was introduced by Meerschaert et al. (1999) and Benson et al. (2000a) as

Here, $M(d{{\color[rgb]{0,0,0}\bm{\theta}}})$ denotes a nonnegative measure on the angle ${{\color[rgb]{0,0,0}\bm{\theta}}}$ in the unit sphere $\{|{{\color[rgb]{0,0,0}\bm{\theta}}}|=1\}$ in $\mathbb{R}^{d}$, and the Riemann-Liouville directional derivative is given by

Benson et al. (2000a) showed that when the measure $M$ is uniform, the operator (31) reduces to the fractional Laplacian (17). In higher dimensions and for general measures $M$, the operator (31) plays an analogous role to the operator in the right-hand side of (27), which is in fact a special case of it for $d=1$. As such, it is used in models of anistropic multivariate $\alpha$-stable Lévy diffusion.

The superdiffusive model introduced above, in which a plume of particles spreads out in space with higher-order rate $t^{1/\alpha}$ for $0<\alpha<2$, raises the question of whether a process can be constructed which results in diffusion characterized by a lower-order rate than the Brownian rate $t^{1/2}$. In this section, we introduce such a model, constructed as Brownian motion with random waiting times drawn from a skewed stable distribution, supported over positive real numbers with a power-law tail. Here, we step away from the framework of the SDE given by (8). Rather than being defined by a simple time-stepping scheme with i.i.d. increments, the paths of the process are defined by a transformation, or “postprocessing”, of Brownian paths $B_{t}$.

We introduce Brownian motion with waiting times, denoted by $B_{\tau(t)}$. The intuition is that the particle paths traced out in space by a discretization of $B_{\tau(t)}$ are paths of discretized Brownian motion $B_{t}$, but the particles wait at each point of the path for a random time drawn from the totally skewed stable distribution. The operational time $\tau(t)$, which introduces waiting and replaces linear time $t$, is an inverse stable subordinator. This is a stochastic process in the variable $t$, although we write $\tau(t)$ rather than using a subscript for typographical reasons. This process is constructed by first defining the stable subordinator $D(t)$, and defining $\tau(t)$ to be the inverse process⁵⁵5The definition of $\tau$ in terms of $D$ is an example of a right-continuous inverse of an increasing functions. Paths of $D$, thought of as functions of $t$, are nondecreasing, so that each path of $\tau$ constructed in this way is a continuous-from-the-right inverse of the parent path of $D$ used to construct it. of $D(\tau)$. Both $D(t)$ and $\tau(t)$ are nondecreasing processes with units of time. In terms of paths, $\tau(t)$ arises from $D(t)$ as

Intuitively, $D(t)$ represents a cumulative waiting time process, keeping track of the total time waited by a particle throughout a path, while the inverse $\tau(t)$ represents an operational time, i.e., the time spent traveling. The increments of $D(t)$ represent the time waited at each location of a particle before the jump to the next location. More specifically, $D(t)$ is a totally skewed $\beta$-stable Lévy process (28) with stability index $\beta\in(0,1)$, $\gamma=1$, scale $\sigma=\cos(\pi\beta/2)$, and center $\mu=0$; see Meerschaert and Sikorskii (2011), Example 5.14. The construction of sample paths of $B_{\tau(t)}$ is demonstrated in Figure 7.

The resulting probability density function of $D(t)$,

for waiting times is supported in nonneagative real numbers. Due to the nonnegative support of the waiting time density, the characteristic function (10) yields the Laplace transform of the waiting time density as

where the Laplace transform is defined as

See Meerschaert and Sikorskii (2011), p. 108 and p. 156 for a discussion. The variance of the process $B_{\tau(t)}$ is given by

which is the desired subdiffusive property. Note that the finiteness of the variance does not imply that the normal central limit theorem applies to $B_{\tau(t)}$, which is not equal in distribution to Brownian motion nor to any Lévy process. In fact, $B_{\tau(t)}$ is not a Markov processes.

The probability density of Brownian motion with waiting times $B_{\tau(t)}$ is governed by the time-fractional diffusion equation,

Here, the Caputo derivative is defined for $0<\beta<1$ by

For $a=0$, this operator is characterized by the simple Laplace transform representation (see Meerschaert and Sikorskii (2011), page 111)

Higher order Caputo derivative can be defined, although the Laplace transforms of the resulting operators involve initial conditions for derivatives of $u$; see Section 2.3 of Meerschaert and Sikorskii (2011). The Caputo derivative is most frequently utilized as a derivative in time for initial-value problems, with the fractional order $0<\alpha<1$.

Before introducing the fundamental solution to the time-fractional diffusion, we introduce the Mittag-Leffler function (Mainardi et al., 2007; Mainardi, 2020)

This Mittag-Leffer $E_{\theta}(z)$ reduces to the exponential function $e^{z}$ when $\theta=1$, and has Laplace transform property

which immediately implies that $E_{\beta}(-k^{2}t^{\beta})$ solves the fractional ordinary differential equation

Returning to the diffusion equation (38) with initial condition $u(x,t=0)=\delta(x)$, applying the Fourier transform in space implies that

which, as shown by Mainardi et al. (2007), yields a solution that can be written

with

being a special case of the Fox-Wright function. Note that $U(x)=u(x,t=1)$. While the fundamental solution above is transcendental, it has the following properties: for $\alpha=1$, it reduces to the solution (5) of the classical diffusion equation; for $0<\alpha<1$, the solution decays faster than exponential and slower than Gaussian; and the second moment of the solution is

Note that the $t^{\beta}$ scaling of this second moment is consistent with the scaling of the fundamental solution above.

Both the $\alpha$-stable Lévy flight $X^{\alpha,p}_{t}$, which led to the space-fractional diffusion equation discused in Section 2.1.3, and Brownian motion with $\beta$-stable subordinator operational time $B_{T^{\beta}(t)}$, which led to the time-fractional diffusion equation discussed in Section 2.1.4, are examples of continuous time random walks (Meerschaert and Sikorskii, 2011). A continuous time random walk (CTRW) allows for a general family of processes in space to be time-changed by a general family of waiting-time processes. To illustrate this concept, we consider the process $X^{\alpha,p}_{T^{\beta}(t)}$, which is $\alpha$-stable Lévy flight $X^{\alpha,p}_{t}$ defined at the discrete level by (28) time-changed by the $\beta$-stable subordinator process $t\mapsto T^{\beta}(t)$ introduced in Section 2.1.4. This models a particle that performs independent jumps drawn from the $\alpha$-stable process, waiting at each point for a random time drawn independently from the $\beta$-stable subordinator process. As shown by, e.g., Meerschaert and Sikorskii (2011) (Section 4.5), the probability density of this particle position is then governed by a differential equation that is fractional in both time and space,

While intuitive, this result deserves a more detailed outline within the general theory of CTRWs. In the standard CTRW model, particles wait at a location for time drawn from a density function $\psi$, and jump to a new location by an increment drawn from a density function $\phi$. The waiting time and jump samples are assumed to be i.i.d., and uncoupled from each other (Zaburdaev et al., 2015; Metzler and Klafter, 2000a; Scalas et al., 2004; Torrejon and Emelianenko, 2018). Thus, the densities $\psi$ and $\phi$ completely determine the CTRW. From the waiting time density $\psi$, the probability that a particle will remains at any given position for time $t$ is

this is referred to as the survival probability of a CTRW particle. Then, given an initial probability density of a particle $u_{0}(x)=u(x,t=0)$, which can also be thought of as an initial distribution of an ensemble of independent particles, the following equation was derived by Montroll and Weiss (1965) for the density at later times:

This equation is central to the CTRW theory⁶⁶6This equation was also derived by Scher and Lax (1973a, b) and is referred to as the CTRW equation of Scher and Lax by Klafter and Silbey (1980). Other authors, such as Torrejon and Emelianenko (2018) refer to this as the master equation of a CTRW.. Taking the Laplace transform in time, the Fourier transform in space, and solving for $\mathcal{F}\left[\mathcal{L}[u]\right](\xi,s)$ yields the Montroll-Weiss equation (Montroll and Weiss, 1965),

In the case that $\phi$ is the $\alpha$-stable density (28) and $\psi$ is the $\beta$-stable subordinator density (34), then $\mathcal{F}[\phi]$ is given by the analytical formula (10) and $\mathcal{L}[\psi]$ by (35), so that the Montroll-Weiss equation represents a closed-form solution of $u$ in $(\xi,s)$-space. Unsurprisingly, it is impossible to perform inverse transforms and obtain $u$ itself analytically, but $u$ can be shown to satisfy (49) using the representations (41) and (25) (Meerschaert and Sikorskii, 2011).

Superdiffusive $\alpha$-stably Lévy flight exhibits infinite MSD, which is a drawback for certain applications. Related to this is the infinite speed of propagation intrinsic to Lévy flights, i.e., the fact that particles have a nonzero probability of traveling an arbitrary large distance in a unit of time. Brownian motion also suffers from this feature, although this probability of large excursions is so low that MSD remains finite. A prototypical model of superdiffusion that cannot be described by a Lévy flight is ballistic motion, in which particles simply move from an initial configuration in fixed random directions with speed $v$, for all time $t$. A ballistic particle travels a distance $vt$ in time $t$ from an initial position $x_{0}$. If reorientations are allowed, then the positions of these so-called sub-ballistic particles in space-time are confined to a ballistic cone

Because the density function of the particle positions is compactly supported, all moments of the position are finite. Such a process cannot be described by Lévy flights.

To capture such behavior, we introduce the Lévy walk model, following Zaburdaev et al. (2015). Such models are based on continuous-in-time motion of particles, rather than instantaneous jumps. A speed $v$ of particles in a medium is specified; each particle moves with speed $v$ in a chosen direction, before a reorientation event occurs in which the direction changes instantaneously and the particle continues to move with speed $v$ before the next direction. Assuming the direction at reorientation is sampled uniformly on the unit sphere, such a walk is determined by a probability density function for the duration of movement $\psi(\tau)$. This leads to a survival probability $\Psi(t)$ given by (50), with $\psi$ now representing the duration density. Thus, $\Psi(t)$ returns the probability that a particle has persisted in a given direction for time $\tau$, i.e., has not experienced reorientation for time $\tau$. Similar to the CTRW case, a master equation can be derived for the probability density $u(x,t)$ of the location of the particle in Laplace-Fourier space:

Unlike the master equation for CTRWs, this equation exhibits coupling in Fourier and Laplace variables, representing coupling in space-time⁷⁷7A Lévy walk may be compared to a non-standard CTRW in which waiting times prior to jumps are correlated to the jump length, e.g., propertional to the jump length, so that long excursions are penalized by long waiting times. See . This results in governing equations that are considerably more complex than those of a standard CTRW. For a Lévy walk, $\psi$ is taken to be a Pareto-type distribution,

An asymptotic expansion of $\mathcal{L}[\psi]$ and $\mathcal{L}[\Psi]$ substituted in (56) yields the following approximation for the evolution of the density function of a Lévy walk in Fourier-Laplace space:

Given $v$ and $u_{0}$, this equation can be inverted to compute $u(x,t)$, but obtaining a governing equation in $(x,t)$ is less straightforward from this point on, due to space-time coupling. Sokolov and Metzler (2003) suggest defining a fractional material or substantial derivative

in order to obtain a governing equation for $u(x,t)$. Recent works, such as those of Chen and Deng (2015), have explored numerical discretizations for these operators.

Despite the greater mathematical difficulties related to governing equations, as compared to other fractional models, Lévy walks have been widely used due to the physical nature of finite speed of propagation and finite MSD; see Zaburdaev et al. (2015) for a survey. When $1<\gamma<2$, by numerical approximations, it can be seen that $u(x,t)$ evolves from a $\delta$-distribution with “a central part of the profile approximated by the Lévy distribution sandwiched between two ballistic peaks” that propogate at speed $v$, with an MSD and self-similarity property for large $t$ that features a superdiffusive scale factor of $t^{1/\gamma}$ Zaburdaev et al. (2015).

Given the physical meaning within stochastic models of the fractional order $\alpha$ in derivatives such as (17), (23), and (40), it is reasonable to expect that these parameters may vary in space and time. Variable-order fractional models are convenient to describe anomalous diffusion in the case of heterogeneous materials or media, or, more generally, when the nature of the diffusion process (subdiffusive, superdiffusive, and classical) changes with space and time. While models with constant fractional order are the simplest and most widely used, some of the model descriptions we discuss in the following sections are improved by the use of a variable fractional order. In recent years, with the purpose of increasing the descriptive power of fractional operators, new models characterized by a variable fractional order have been introduced for both space- and time-fractional differential operators Antil and Rautenberg (2019a); Darve et al. (2021); D’Elia and Glusa (2021); Razminia et al. (2012); Zheng and Wang (2020b) and several discretization methods have been designed Chen et al. (2015); Schneider et al. (2010); Zeng et al. (2015); Zheng and Wang (2020a); Zhuang et al. (2009). The improved descriptive power of variable-order fractional operators has been demonstrated in some recent works on parameter estimation Pang et al. (2020a, 2019a); Zheng et al. (10 Nov. 2020).

Given a function

i.e., a function $\alpha(\mathbf{x},t)$ of space and time, we define variable-order operators as follows. For a function $u(\mathbf{x},t)$ with $\mathbf{x}\in\mathbb{R}^{d}$ and $t\in\mathbb{R}$, we define the variable-order fractional Laplacian⁸⁸8For more recent works and novel definitions of variable-order fractional Laplacians we refer the reader to Antil and Rautenberg (2019b); D’Elia and Glusa (2021). as

Here, $\alpha(\mathbf{x},t)$ is restricted to take values in $(0,2)$. Note that for constant $\alpha$, $\mathfrak{L}^{\alpha}=(-\Delta)^{\alpha/2}$. For $d=1$ and $\alpha(x,t)$ restricted to $(0,1)$, we define the variable-order left-sided Riemann-Liouville fractional derivative as

The right-sided Riemann-Liouville may be defined for variable order in an analogous way. We define the variable-order Caputo fractional derivative, again for $\alpha(x,t)$ taking values in $(0,1)$, as

To summarize and offer a quick look-up of anomalous diffusion processes, their corresponding fractional models, and applications of each process/model, we have included these relationships in Table 1. This table includes references to the previous sections where each process and model is described, as well as pointers to the applications in the following sections where the models are utilized. We have limited references to applications to only those three areas that we focus on in this article.