Slices of the Anomalous Phase Cube Depict Regions of Sub-and Super-Diffusion in the Fractional Diffusion Equation

: Fractional-order time and space derivatives are one way to augment the classical diffusion equation so that it accounts for the non-Gaussian processes often observed in heterogeneous materials. Two-dimensional phase diagrams—plots whose axes represent the fractional derivative order— typically display: (i) points corresponding to distinct diffusion propagators (Gaussian, Cauchy), (ii) lines along which speciﬁc stochastic models apply (L é vy process, subordinated Brownian motion), and (iii) regions of super- and sub-diffusion where the mean squared displacement grows faster or slower than a linear function of diffusion time (i.e., anomalous diffusion). Three-dimensional phase cubes are a convenient way to classify models of anomalous diffusion (continuous time random walk, fractional motion, fractal derivative). Speciﬁcally, each type of fractional derivative when combined with an assumed power law behavior in the diffusion coefﬁcient renders a characteristic picture of the underlying particle motion. The corresponding phase diagrams, like pages in a sketch book, provide a portfolio of representations of anomalous diffusion. The anomalous diffusion phase cube employs lines of super-diffusion (L é vy process), sub-diffusion (subordinated Brownian motion), and quasi-Gaussian behavior to stitch together equivalent regions.


Introduction
Fractional calculus (non-integer order integration and differentiation) extends the local, memory-free concepts of Newton and Leibniz [1,2], but this growth entails a mosaic of noninteger derivatives whose diverse properties can overwhelm new users. Given the long history and distinguished pedigree of fractional calculus [3,4], one would think that by now it would have an established niche in mathematical physics (as have generalized functions and pseudo-differential operators [5]). Nevertheless, while widely used in electrochemistry, viscoelasticity, dielectric relaxation and anomalous diffusion, the tools of fractional calculus are often overlooked. This paradox suggests that we search for new ways to introduce the multiple forms of fractional-order derivatives (Grünwald-Letnikov, Riemann-Liouville, Caputo, Riesz-Weyl, etc.) to new investigators seeking to apply fractional calculus in their research.
In the field of magnetic resonance imaging (MRI), this question can be focused on the Bloch-Torrey equation [6], and on the selection of the proper fractional-order tools for generalizing relaxation and diffusion phenomena that occur in complex biological tissues. Incorporation of a fractional time derivative into a first order relaxation equation yields a single parameter Mittag-Leffler function that enhances the integer-order exponential decay model [7], while employing a fractional space derivative in the classical diffusion equation provides a stretched exponential solution that expands the Gaussian model to one of Lévy motion [8]. Selection of a particular fractional derivative for each case is based on anticipating how material composition, heterogeneity and, in some cases, fractal structure might be incorporated into the analysis [9,10].
What can guide the researcher in choosing one fractional derivative over another? When should the material parameters in the models be constants or functions of position and time? Why use only a linear model? These questions quickly take us beyond the scope of the current paper but provide context for this investigation. Specifically, here we will examine fractional-order models that generalize the classical Gaussian solution to the diffusion equation as it is encountered in MRI. The goal is not to apply a particular model to a given set of data, but to compare and contrast different models (e.g., integerorder, fractional-order; linear, non-linear; constant diffusion coefficients, variable diffusion coefficients). The comparison tool we have chosen will be a scaling analysis of the governing one-dimensional diffusion equation. We will illustrate the results using a phase diagram as a means of classifying the time-dependence of the mean-squared displacement of the underlying random motion. We will illustrate the behavior of each model by presenting the phase diagram as an orthogonal slice of a three-dimensional anomalous diffusion phase cube. Lines of super-diffusion (Lévy process) and sub-diffusion (subordinated Brownian motion, fractional Brownian motion, etc.) for each two-dimensional phase diagram will be seen to link different diffusion models-displaying corresponding regions of sub-and super-diffusion.

Background
There are now available more than one hundred books on fractional calculus, each describing specific aspects of the formulation, properties and applications of fractional derivatives in physics, chemistry, biology, and engineering (https://mechatronics.ucmerced.edu/ fcbooks) (accessed on 20 May 2021). An annotated bibliography of books that we have found useful can be found in Chapter XII of the book, Fractional Calculus in Bioengineering [4]. In addition, for reference, the definitions and basic properties of the fractional operators used in this paper are listed in the Supplementary Materials (Tables S1 and S2) of this paper. These tables and Appendices A and B of this paper should help researchers who are not familiar with specific properties of different forms of fractional integration and differentiation.
In addition, it is an inconvenient truth that the Greek letters often used in fractional calculus to indicate, non-integer, fractional-order operations (α, β, γ, . . . ) are not applied in a consistent manner for time and space derivatives but shift with the author's preference and historical precedent. Since Jeremy Bentham's 'calculus of felicity' (https://plato. stanford.edu/entries/) (accessed on 20 May 2021) is not part of fractional calculus, the operator notation in fractional calculus was not formulated to keep all users happy; hence, we must be ever wary of the local terminology. It should also be kept in mind that interchanging operations between fractional and integer-order calculus should be viewed with suspicion; one should not mix the tools and notation together any more carelessly than trying to assemble a motor using parts with both metric and British Standard Whitworth threads [11]. Consequently, in the results below we will employ a restricted set of fractional derivatives and identify the local notation in a consistent manner for each case.

Results
In this section two groups of fractional order diffusion equations, Sections 3.1 and 3.2, are considered where β, denotes the order of the fractional time derivative, and α, indicates the order of the fractional space derivative. In the first group, Section 3.1, Table 1, we assume a linear diffusion equation with a diffusion coefficient that is a power law function of time, D β,α,ν (t ν ), where |ν| < 1. The general solution for each of these equations can be expressed in terms of the Fox, H-function [10], and its Fourier transform has been used to model diffusion studies in MRI [12]. In the second group, Section 3.2, Table 2, we assume both linear and non-linear fractional order models of diffusion with a diffusion coefficient that is a power law function of both space and time, D β,α,θ, ν (|x| −θ t ν ), where θ > 0 and |ν| < 1. The general solution for this group is not available in a closed form and, only a few special cases have been used in MRI. For both groups we use scaling principles [13] to determine the mean squared displacement (MSD) for a multi-dimensional, Euclidean phase space, {β, α, ν, θ, . . .}. Two-dimensional phase diagrams of each model are displayed, for convenience, as a slice of a three-dimensional anomalous phase cube, e.g., {β, α, ν} or {β, α, θ}. The MSD appears as a line of quasi-Gaussian behavior on the phase diagram separating regions of sub-and super-diffusion. Table 1. Anomalous diffusion with power law diffusion coefficient, D β,α,ν (t ν ).
In this table, α is the order of symmetric Riesz space derivative is (0 < α < 2), β is the order of the Caputo time derivative (0 < β < 1), υ represents the power law behavior of the diffusion coefficient (−1 < υ < 1), and P(x,t) is particle density. FTS is fractional time and space diffusion, FT is fractional time diffusion, and FS is fractional space diffusion. Also, x 2 , is the mean squared displacement, and D β,α,ν (t ν ) = D β,α,ν (1 + ν)t ν with the units, mm α /s υ+β . Table 2. Anomalous diffusion with power law diffusion coefficient, In this table, FTS is fractional time and space diffusion with a power law space diffusion coefficient, NFS is non-linear fractional space diffusion with a power law space diffusion coefficient, FTD is fractional time diffusion with a power law space diffusion coefficient, and P(x,t) is the particle density. Also, x 2 , is the mean squared displacement, and D β,α,ν,θ,ρ has the units of mm (α+θ) /s (β+ν) P ρ−1 .
In the first group of equations, Section 3.1, we describe linear fractional time and space diffusion where the diffusion coefficient is assumed to vary as function of time, D β,α,ν (t ν ). The first diffusion example (Section 3.1.1) describes the floor of the phase cube as a two-dimensional phase diagram {α, β}, for fractional time and space derivatives with a constant diffusion constant, D β,α , with dimensions, mm α /s β . The next two cases describe faces of the phase cube {β, α, ν} that intersect the fractional derivative floor along the line β = 1 (fractional motion, Section 3.1.2), and along the line α = 2 (fractional Kilbas-Saigo decay, Section 3.1.3). In both cases the diffusion coefficient varies with time as a power law.
In the second group of equations, Section 3.2, we describe linear fractional time diffusion with linear and non-linear space diffusion of degree, {ρ}, where the diffusion coefficient is assumed to vary as function of space and time, D β,α,θ, ν (|x| −θ t ν ). The first example (Section 3.2.1) is linear, ρ = 1, with ν = 0, and describes the phase cube {β, α, θ}. The next example (Section 3.2.2) is a non-linear model of diffusion with a power law space diffusion coefficient, D 1,α,θ,ρ (|x| −θ ), and describes the phase cube {α, θ, ρ}. The last example (Section 3.2.3) is a linear diffusion model with a fractional order time and an integer order space derivative, and a diffusion coefficient that varies as a power law in both space and time, D β,α,ν,θ t ν |x| −θ , so that corresponding to the phase polytope is a tesseract {β, α, θ, ν}.
In the following, the diffusion equation will be described in terms of the particle density, P(x,t), with the units of particles/mm 3 , where the particle can be a molecule a microbe, or any object whose physical size is small compared to its surroundings.

Fractional Time and Space Derivatives with Power Law D(t)
Consider the one-dimensional fractional diffusion equation of P(x,t) that includes a time-varying power law diffusion coefficient D(t) = D β,α,ν (1 + ν)t ν with the Caputo time and symmetric Riesz space derivatives: Here, the order of symmetric Riesz space derivative is α, (0 < α < 2), the order of the Caputo time derivative is β, (0 < β < 1), and the Greek letter ν, (−1 < υ < 1) is assigned to the assumed power law dependence of diffusion coefficient. The diffusion coefficient in this case is a constant, D β,α,ν (1 + ν), multiplying t υ . The units of D β,α,ν are mm α /s υ+β . The solution to Equation (1) for an initial value problem with a Dirac delta function concentration of material at the origin of an unbounded domain, P(x, 0) = δ(x), is given in Appendix A.2. The result can be written as an inverse Fourier transformation of the Kilbas-Saigo function, If we assume that Equation (2) can be written in the scaled form, with ξ = (ν + β)/α, [14][15][16], then the corresponding mean squared displacement (MSD) is: where I β,α,ν may not exist for all values of β, α and ν (See Appendix A.1 for details). Analyzing the time behavior of the mean squared displacement in Equation (4) leads us to recognize sub-diffusion, normal diffusion, and super-diffusion when: 2(υ + β) < α, 2(υ + β) = α, and 2(υ + β) > α, respectively. The three fractional parameters of the anomalous phase cube {α, β, υ} describe: (i) Gaussian diffusion with a power law D(t ν ), (ii) fractional time and space diffusion (FTS) with a power law D(t ν ), (iii) fractional space diffusion (FS) with a power law D(t ν ), and (iv) fractional time diffusion (FT) with a power law D(t ν ), models of anomalous diffusion-all subsumed by Equation (1)-and displayed with the corresponding power law time-dependence of the MSD in Table 1.

Diffusion Equation with Fractional Time and Space Derivatives
The continuous time random walk (CTRW) model of anomalous diffusion is widely used because it provides a direct connection between the statistical properties of particle trapping and jump increments (inverse power laws) with the order of the time and space fractional derivatives in the associated diffusion equation [17,18].
Consider the one-dimensional diffusion equation for the CTRW model of P(x,t) expressed in terms of fractional time and space derivatives: ( , ) = , ( , ) | | (5) Here, the Greek letter , (0 < < 1), is used to designate the order of Caputo fractional time derivative, while the Greek letter , (0 < < 2), is used to designate the order of the symmetric Riesz fractional space derivative (both derivatives are defined in Supplementary Tables S1 and S2). The diffusion coefficient in this case, , , is a constant, and has the units of / .
The solution to Equation (5) when all particles are initially concentrated at the origin of an unbounded domain, ( , 0) = ( ), is a Fox, H-function [19,20], which can be written as: where for = / / we have, The corresponding MSD is formally given by: where , = , ( ) , with the understanding that the MSD does not exist for the entire range of the fractional order parameters. Analyzing the behavior of the 〈 〉 in

Diffusion Equation with Fractional Time and Space Derivatives
The continuous time random walk (CTRW) model of anomalous diffusion is widely used because it provides a direct connection between the statistical properties of particle trapping and jump increments (inverse power laws) with the order of the time and space fractional derivatives in the associated diffusion equation [17,18].
Consider the one-dimensional diffusion equation for the CTRW model of P(x,t) expressed in terms of fractional time and space derivatives: Here, the Greek letter β, (0 < β < 1), is used to designate the order of Caputo fractional time derivative, while the Greek letter α, (0 < α < 2), is used to designate the order of the symmetric Riesz fractional space derivative (both derivatives are defined in Supplementary Tables S1 and S2). The diffusion coefficient in this case, D β,α , is a constant, and has the units of mm α /s β .
The solution to Equation (5) when all particles are initially concentrated at the origin of an unbounded domain, P(x, 0) = δ(x), is a Fox, H-function [19,20], which can be written as: where for z = x/t β/α we have, The corresponding MSD is formally given by: where I β,α = ∞ −∞ P α,β (z)z 2 dz, with the understanding that the MSD does not exist for the entire range of the fractional order parameters. Analyzing the behavior of the x 2 in Equation (8) leads us to recognize sub-diffusion, normal diffusion, and super-diffusion when: 2β < α, 2β = α, and 2β > α, respectively.
The Fourier transformation of Equation (6), p(k, t) can be expressed as a one-parameter Mittag-Leffler function [20]: This solution, which is the characteristic function of the P(x,t), also describes the time decay of each spatial frequency of the associated random motion. It yields particular cases for selected values of the fractional derivatives: Gaussian (β = 1, α = 2), Lévy process (β = 1, α), and subordinated Brownian motion for (β, α = 1). These special cases of anomalous diffusion for the CTRW model are displayed in Figure 1 as a two-dimensional {β, α} phase diagram. This phase diagram corresponds with the ν = 0 plane of the anomalous phase cube {β, α, ν} in Figure 1.

Integer Time Derivative and Fractional Space Derivative with Power Law D(t)
This anomalous diffusion model was defined by Eliazar and Shlesinger as 'fractional motions' [21]. It differs from the CRTW model by allowing persistence and anti-persistence in the jump increments of particle motion. This behavior is incorporated into a modified fractional diffusion equation by setting the order of the time derivative, β = 1, and introducing a time-dependent diffusion coefficient, which is assumed to be a power law. Eliazar and Shlesinger provide a complete description of this extended model of Brownian motion in [21] where the rare extremes of large jumps is described as the Noah exponent (α in our notation), while persistence, or anti-persistence of motion is characterized as the Joseph exponent (ν in our notation). They also show its connection to Lévy motion and the Hurst exponent (H). The utility of this diffusion model in MRI is illustrated by Fan and Gao [22] and by Karaman et al. [23,24].
Consider the one-dimensional diffusion equation for the fractional motions model of P(x,t) expressed in terms of the integer time derivative and fractional space derivative: Here, the Greek letter α, (0 < α < 2) is used to designate the order of the symmetric Riesz fractional space derivative, while the Greek letter ν, (−1 < ν < 1) is the assumed power law dependence of the time-dependent diffusion coefficient, D α,ν (t). The diffusion coefficient in this case includes a constant, D α,ν (1 + ν) multiplying t ν , where D α,ν has the units of mm α /s (1+ν) .
The solution to Equation (10) for an initial concentration of material at the origin of an unbounded domain is a fractional Lévy motion [21] can be written as: where for z = x/t 2(1+ν)//α we have, The corresponding mean squared displacement is: with I α,ν = ∞ −∞ P α,ν (z)z 2 dz, with the understanding that the MSD does not exist for the entire range of the fractional order parameters The Fourier transformation of Equation (11), p(k, t) has a simple stretched exponential solution: This solution describes the time decay of each spatial frequency component. It yields particular cases for selected values of the fractional derivatives: Gaussian (α = 2, ν = 0), Cauchy (α = 1, ν = 0), Levy process (α, ν = 0), and fractional Brownian motion for (α = 2, ν). The governing diffusion equations for these cases are shown in the Supplementary Table S3.

Fractional Time Derivative and Integer Space Derivative with Power Law D(t)
The fractional power law model complements the 'fractional motions' model by assuming a fractional time derivative, β, 0 < β < 1, but with the order of the space derivative, α = 2. In the context of the CTRW model, particle trapping is proscribed, while jump increments must have a finite second moment. Introducing a time-dependent diffusion coefficient provides an independent time-scaling parameter, one that has been associated with fractal properties of the medium [fractal/fractional model]. A complete description of this extended model of Brownian motion and its connection to CTRW motion is provided by Gorenflo et al. in [7]. The utility of this diffusion model in MRI is illustrated by Fan and Gao [22] and by Magin et al. [12].
Consider the one-dimensional diffusion equation of P(x,t) with a power law D(t) expressed in terms of a fractional time and an integer space derivative: Here, the order of space derivative, α, is set to two, the order of the Caputo time derivative is β, (0 < β < 1), and the Greek letter ν, (−1 < υ < 1) is assigned to the assumed power law dependence of diffusion coefficient. The diffusion coefficient in this case is a constant, D β,ν (1 + ν), multiplying t υ . The units of D β,ν are mm 2 /s υ+β . The solution, P(x,t), for an initial concentration of material at the origin of an unbounded domain is a Kilbas-Saigo function [10] (p. 21): which may also be written as: with, The corresponding mean squared displacement is: with the understanding that the MSD does not exist for the entire range of the fractional order parameters.
Classification of this behavior in terms of sub-and super-diffusion in the {β, ν} plane is displayed in Figure 1. The mean squared diffusion for this model of anomalous diffusion has a time dependence of the form t υ+β . This expression corresponds to quasi-diffusion when υ + β = 1, super-diffusion when υ + β > 1, and sub-diffusion when υ + β < 1. This phase diagram corresponds with the α = 2 plane on the right side of the anomalous phase cube {α, β, υ} in Figure 1.

Fractional Time and Space Derivatives with Power Law D(x,t) and Non-Linear in P(x,t)
The classical derivation of the diffusion equation is a combination of a flux defined as the negative gradient of particle concentration, [25]. Extending this derivation using fractional calculus suggests-in addition to introducing fractional operators for the time and space derivatives-that one might have reason to modify either the flux or the continuity equation (or both) [26,27]. This approach also introduces the notion of sequential fractional space derivatives [1,3], and the possibility of inserting a space or time dependent diffusion coefficient in the definition of the diffusion coefficient [10], as well as the possibility of a non-linear flux. Although the last step might seem arbitrary, it is consistent with the porous media model of the diffusion equation [28], the extension of diffusion from one to n dimensions [29], the characterization of diffusion across fractal interfaces [30], and the description of porous materials in terms of fractal dimensions [31]. In this section we will investigate a linear/non-linear model that describes fractional diffusion combined with an inverse power law time and space dependence diffusion coefficient.
The general solution to Equation (21) for an initial value problem with a Dirac delta function concentration of material at the origin of an unbounded domain, P(x, 0) = δ(x), is not available. If we assume that the normalized solution to Equation (21) can be written in the scaled form: then by the approach used in [32,33] we find (Appendix B.1), ξ = (β + ν)/(ρ + α + θ).
The corresponding mean squared displacement (MSD) is: where I ξ = ∞ −∞ z 2 p(z)dz for z = x/t ξ and with the understanding that the MSD does not exist for the entire range of the fractional order parameters (see Appendix B.1 and [16] (Chapter 16) for details).

Fractional Time and Space Derivatives with Power Law, D(x), and Linear in P(x,t)
Consider the one-dimensional fractional diffusion equation of P(x,t) that includes a space-varying power law diffusion coefficient D(|x| −θ ) with the Caputo time and symmetric Riesz space derivatives: The corresponding mean squared displacement is: Here, the Greek letter, β, (0 < β < 1), is the order of the Caputo fractional time derivative, and the Greek letter α, (0 < α < 1), is the order of the symmetric Riesz fractional space derivative. The diffusion coefficient in this case, D β,α,θ , is a constant and has the units of mm (α+θ) /s β . The general solution to Equation (24) is not known. Solutions of this equation for an initial concentration of material at the origin of an unbounded domain are presented for specific cases in [10] and Appendix B.2.
For the α = 1 case, when all particles are initially concentrated at the origin, P(x, 0) = δ(x), the solution to Equation (24) in an unbounded domain is a Fox, H-function [36] (Appendix B.2), which can be written as: where, with z = x/t β/(2+θ) , and, In the general case, this model fills the anomalous phase cube {β, α, θ} with regions of sub-and super-diffusion that correspond to the situations where the exponent of the MSD, 2β/(1 + α + θ), is less than or greater than one. Classification of this behavior in terms of sub-and super-diffusion is displayed in Figure 2. When θ = 0, the super-diffusion domain of the β = 1 plane merges with the fractional diffusion model (Lévy motion), while also for θ = 0 the sub-diffusion domain of the α = 1 plane merges with the fractional diffusion model (subordinated Brownian motion). When α = β = 1, the vertical line, z = θ reflects the stretched exponential solution for D(x) = D θ |x| −θ , which is Gaussian for θ = 0, and exhibits sub-diffusion when θ is positive and super-diffusion when θ is negative.

Integer Time Derivative with a Non-Linear Space Derivative and a Power Law
Diffusion in complex materials has been modeled using a non-linear version of the classical diffusion equation (porous media equation) where P(x,t) is raised to a power, ρ [37][38][39]. In this section we will outline the fractional order generalization of this non-linear model, largely following the work described in [37,39], and in Chapter 6 of the book on fractional diffusion equations and anomalous diffusion by Evangelista and Lenzi [10].
Consider the one-dimensional non-linear fractional diffusion equation for P(x,t) expressed in terms of the symmetric Riesz space derivative and a power law diffusion coefficient D(|x| −θ ): The corresponding mean squared displacement is derived from a scaling argument in Appendix B.3: with , and with the understanding that the MSD does not exist for the entire range of the fractional order parameters. Here the order of the fractional space derivative, α is (0 < α < 1), θ is a positive real number, and the range of the fractional power ρ is (−1 < ρ < 2). The diffusion constant, D α,θ,ρ has the units mm (α+θ) /sP ρ−1 .
Using the scaling argument in Appendix B.3, and expressing P(x, t) in terms of z = x/t ξ with ξ = 1/(θ + α + ρ), as: it is possible to simplify Equation (31) as: For an initial concentration of material at the origin of an unbounded domain particular solutions to Equation (34) are [38]: For the case of α = 1, and in the case of α = 1, with where N and N are a normalization factors, and a and b are chosen such that Equation (36) satisfies Equation (34). The parameter b is chosen ±1, depending on the behavior of the solution. In fact, for a compact behavior, we have b = −1 and for a long-tailed behavior b = 1. Note that the solutions shown for the non-linear case are given in terms of power laws and may be related to the Levy distributions as discussed in [10,[37][38][39].
In the general case, this model fills the anomalous phase cube {α, θ, ρ} with regions of sub-and super-diffusion that correspond to the situations where the exponent of the MSD, 2/(θ + α + ρ), is less than or greater than one. Classification of this behavior in terms of sub-and super-diffusion is displayed in Figure 3, which is one of the eight phase cubes that comprise the {β, α, θ, ρ} tesseract.
The corresponding mean squared displacement is derived from a scaling argument in Appendix B.4: where , , = ( ) for = / ( )/( ) with the understanding that the MSD In Figure 3, the linear model of anomalous diffusion is displayed as the top of a phase cube (ρ = 1), as a two-dimensional {α, θ}, phase diagram that is divided into regions of super-and sub-diffusion by the line, θ = 1 − α, with the Gaussian diffusion case at the point {α = 1, β = 1, θ = 0, ρ = 1}. The ρ = 1 plane is intersected by the θ = 0 and the θ = 1 planes illustrating complementary phase behavior for the non-linear model when θ = 0. For the α = 1 plane, the line of quasi-diffusion (ρ = 1 − θ, where the MSD is proportional to diffusion time), also has a negative slope when D(x) = D θ |x| −θ . Finally, the ρ = 1, θ = 0, and α = 1 planes slice off a corner of subdiffusion from the phase cube {α, θ, ρ}.

Fractional Time Derivative and Integer Space Derivatives with a Time and Space-Dependent Diffusion Coefficient, and Linear in P(x,t)
Consider the one-dimensional linear diffusion equation for P(x,t), with integer space derivatives and a Caputo time derivative, β, (0 < β < 1): The corresponding mean squared displacement is derived from a scaling argument in Appendix B.4: x 2 = I β,ν,θ t 2(β+ν)/(2+θ) (39) where I β,ν,θ = ∞ −∞ P(z)z 2 dz for z = x/t (β+ν)/(2+θ) with the understanding that the MSD does not exist for the entire range of the fractional order parameters.
Here, the diffusion coefficient in this case, D β,ν,θ is a constant and has the units of mm (2+θ) /s (β+ν) . Solutions of this equation for an initial concentration of material at the origin of an unbounded domain are presented for specific cases (additional cases can be found in [10]).
The general solution to Equation (38) for an initial concentration of material at the origin of an unbounded domain [40] is derived in Appendix B.4 using: with In this case, the model reflects the anomalous phase cube {β, θ, ν} with regions of suband super-diffusion that correspond to the situations where the exponent of the MSD, (β + ν)/(2 + θ), is less than or greater than one. Classification of this behavior in terms of sub-and super-diffusion is displayed in Figure 4, which is one of the eight phase cubes that comprise the {β, α, θ, ν} tesseract. The general solution to Equation (38) for an initial concentration of material at the origin of an unbounded domain [40] is derived in Appendix B.4 using: In this case, the model reflects the anomalous phase cube { , , } with regions of suband super-diffusion that correspond to the situations where the exponent of the MSD, ( + )/(2 + ), is less than or greater than one. Classification of this behavior in terms of sub-and super-diffusion is displayed in Figure 4, which is one of the eight phase cubes that comprise the { , , , } tesseract.  In Figure 4, the linear model of anomalous diffusion is displayed on the top face of a phase cube (ν = 0), as a two-dimensional {β, θ}, phase diagram that is divided into regions of super-and sub-diffusion by the line, β = 1 + θ/2, with the Gaussian diffusion case appearing at the point {β = 1, θ = 0, ν = 0}. The ν = 0 plane is intersected by the θ = 0 plane and the θ = −2 plane (where the MSD does not exist). For the β = 1 plane, the line of quasi-diffusion (ν = 1 + θ/2, where the MSD is proportional to diffusion time), shows a negative slope when D(x) = D β,ν,θ t ν |x| −θ . Finally, the ν = 0, and θ = 0, planes slice off a wedge of sub-diffusion from the phase cube {β, θ, ν}, while the ν = 0, and β = 1, planes slice off a wedge of super-diffusion.

Discussion
The phase diagrams displayed in the results section (Figures 1-4) depict models of anomalous diffusion (Tables 1 and 2) as a slice of a three or n-dimensional solid (anomalous phase hypercube). The intent of this representation is to provide a perspective that connects the underlying mathematical models (fractional diffusion equations) with the expected regions of sub-and super-diffusion. Our overall goal is to present these diagrams as an aid in the interpretation of diffusion data collected from porous, heterogeneous material, such as, that acquired in diffusion-weighted MRI experiments.
Anomalous diffusion, typically characterized by asymptotic, power law tails for P(x, t) and a power law growth in time of the mean squared displacement (MSD) [41], sprouts from the deep root of Brownian motion in its many variegated forms. Here we classify these forms by comparing the power law exponent (R) of the MSD for each model (t 2R ) with the classical result for Brownian motion ( x 2 = 2Dt, where R = 1 2 ). In Section 3.1, we combine the CTRW model of anomalous diffusion with a power-law form of the time-dependent diffusion coefficient, D(t) = D β,α,ν (1 + ν)t ν . For this model, Equation (1), we show that the MSD can be obtained directly from the analytical solution of specific cases (using Fox, H-functions) as well as from a scaling analysis of the governing fractional order differential equation (details in Appendix A). The results include the special cases of Lévy motion, fractional motion and subordinated fractional diffusion.
In Section 3.2, we consider a modified CTRW model of anomalous diffusion derived by assuming a fractional flux gradient, J α (x, t) = −D α ∂ α |x| P(x, t), acting in concert with an integer order divergence term in the continuity equation. The resulting diffusion equation is also allowed to be non-linear, and to include a power-law time-and space-dependent diffusion coefficient, D(x, t) = D β,α,ν,θ,ρ t υ |x| −θ . For this model, Equation (21), we do not have an analytical solution for all cases but use scaling analysis (Appendix B.1) to obtain the MSD for three situations: first, a linear anomalous diffusion equation with a fractional time and a fractional space derivative, and D(x, t) = D α,θ,ρ |x| −θ (Section 3.2.1); second a non-linear anomalous diffusion equation with an integer time derivative, a fractional space derivative, and D(x, t) = D α,θ,ρ |x| −θ (Section 3.2.2); and finally a linear anomalous diffusion equation that uses a fractional time derivative, an integer order space derivative, and includes a time-and space-dependent power law diffusion coefficient D(x, t) = D β,α,ν,θ,ρ t υ |x| −θ (Section 3.2.3). The so-called porous medium equation appears under the umbrella of the non-linear case when both the time and space derivatives are of integer order [42].
Although each model springs from a different physical or mathematical perspective [43,44] (linear/non-linear, fractional/fractal, deterministic/stochastic, and constant/varying diffusion coefficient), its parameter space (order of the time and space derivatives, power-law exponents of the diffusion coefficients, degree of non-linearity and fractal dimensions) is divided into zones of sub-and super-diffusion by considering the behavior of the mean squared displacement for each model. In the phase cube representation, these zones coalesce to a point for the Gaussian behavior of Brownian motion, a line along which Lévy motion or subordinated Brownian motion appear, and a plane where regions of sub-and super-diffusion diffusion are displayed.
Phase diagrams are used throughout physics and engineering to map the behavior of physical systems [45]. Pressure-temperature, pressure-volume, temperature-volume (gases) and material composition (liquids and solids) phase diagrams are widely used in metallurgy and chemistry to guide our understanding of phase changes, optimization of fractional distillation, and chemical extraction procedures. General principles such as the Gibb's phase rule (F = C − P + 2; the degrees of freedom of a system at equilibrium (F) is equal to the number of components (C) minus the number of phases (P) plus 2) guide the expected changes of state when the temperature, pressure or composition are varied. The idea is easily generalized to binary systems at constant pressure and to eutectic mixtures as well as to alloys and amalgams. In such cases slices of a multidimensional phase volume provide a means of organizing the complex material using multiple perspectives.
Chemical phase diagrams of single and multiple components (e.g., water, salt-water, lead-tin, iron-iron carbide) capture reams of experimental data (melting, vaporization, sublimation) in regions representing solid, liquid, and gaseous phases [46]. The molecular configuration in each phase is verified by microscopy and spectroscopy, and the diagrams modified to account for new phenomena (polymorphic solids, glasses, compressible liquids, and supercritical fluids). The result is not just an equilibrium description of the material at a fixed specific volume, but a road map showing how the terrain will change when experiments alter the temperature, pressure or composition. Implied, but not displayed in the phase diagram are the consequences of changing the temperature or pressure on the mobility and configurations of molecules.
Mathematical phase diagrams of fractional derivative orders (β time; and α, space) in the diffusion equation describe a myriad of analytical, numerical and Monte Carlo simulation results for the normal (Gaussian) and anomalous (Levy, fractional Brownian, subordinated, ballistic) motion of ideal particles [27,28]. The particle motion in each phase reflects an underlying material configuration where the mobility is characterized by sub-diffusion, super-diffusion, or quasi-Gaussian motion. Association of experimental measurements of 'anomalous' diffusion in a complex, heterogeneous material with a specific model using fitted values of β and α is not a one-to-one map of fractional derivative order to molecular configuration; it is doubly degenerate. First, at least in nuclear magnetic resonance (NMR) and MRI [47], because there are multiple sets of parameters that can with high fidelity match the observed data, and second because there are integer and fractional-order models of diffusion, some employing a power-law diffusion coefficient D(|x| −θ , t ν ), which also fit the data.
Here, we would like to suggest several ways to combine these models. One idea is to overlap chemical and mathematical phase diagrams, beginning with simple materials. As an example, in Figure 5 we have superimposed a portion Figure 5 (e.g., the {β, α} plane for fractional time and space diffusion with a constant coefficient) on the phase diagram of pure water [45] (triple point, tp, to critical point, cp,). Taking inspiration from the Clausius-Clapeyron equation [46], we associate α = 2 [log(P/P tp )]/[ log(P cp /P tp )] with the pressure ratio, and β = T − T tp /T tp / T( cp −T tp )/T tp with the temperature ratio so that the critical point of the {P, T} diagram corresponds with the Gaussian point of the anomalous diffusion phase diagram. We then see that overlap of the phase diagrams divides the liquid and vapor phases of water into regions of putative sub-and super-diffusion, respectively. The combined diagram also suggests and that the supercritical fluid phase corresponds with the domain of Brownian motion in the mathematical phase diagram. This inchoate connection can perhaps be extended to other single-and multi-phase situations where anomalous diffusion or behavior appears. It also has the potential of linking chemical composition and material structure-for a given material-with the most appropriate slice of the phase cube and thus could provide a basis for selection of one model versus another. This correspondence might be extended by applying fractional calculus to theoretical models of the vapor-liquid equilibria (Chapter IX in [5,48]). Alternatively, we can probe the underlying assumptions of kinetic, fractional order and stochastic models to establish how β and α reflect changes in kinetic theory concepts such as mean free path and the time between collisions with the CTRW probabilities of particle jumps and waiting times. The simplest compartmental representation of fractional derivatives is expressed in terms of a distribution of first order exponential decay processes. This approach, used by Berberan-Santos in luminescence studies [49], links lower order fractional derivatives with a broader distribution of compartments-a single compartment as the order approaches one. The CTRW model assumes particle jump | | ( ) and waiting time ( ) probabilities that by convolution embed material properties in the fractional order derivative operators. This correspondence was applied by Schumer and colleagues (see Figures 7 and 9 in [50]) to problems in erosion and groundwater contamination.
Classical diffusion theory links the time rate of change of particle concentration with the divergence of the flux, which is assumed to be proportional to the gradient of concentration. Here, we used the Caputo fractional time derivative to generalize the time rate of change and the symmetric Reisz fractional space derivative to generalize the flux (in Section 3.2, we have not modified the gradient operator in the continuity equation). The fractional time derivative affects particle dynamics by introducing memory in the time derivative, ( ) * ( ) / (− ), while the fractional space derivative introduces non-locality, ( ) * | | ( ) / (− ), [5,43]. Thus, memory and non-locality are introduced by convolution into mathematical physics when fractional order operators are employed.
Thinking of the classical diffusion problem as it might be affected by Maxwell's demon-as William Thompson did in 1875 (see pages 57-60 in Reference [51])-the demon, is viewed as a batsman in cricket who sorts out pitches/particles entering a 'control' volume based on speed and direction. In the fractional calculus version of cricket, the bowler is also a 'demon', one who can release balls with multiple delays, ( − ) all weighted by  Alternatively, we can probe the underlying assumptions of kinetic, fractional order and stochastic models to establish how β and α reflect changes in kinetic theory concepts such as mean free path and the time between collisions with the CTRW probabilities of particle jumps and waiting times. The simplest compartmental representation of fractional derivatives is expressed in terms of a distribution of first order exponential decay processes. This approach, used by Berberan-Santos in luminescence studies [49], links lower order fractional derivatives with a broader distribution of compartments-a single compartment as the order approaches one. The CTRW model assumes particle jump |x| −(1+α) and waiting time t −(1+β) probabilities that by convolution embed material properties in the fractional order derivative operators. This correspondence was applied by Schumer and colleagues (see Figures 7 and 9 in [50]) to problems in erosion and groundwater contamination.
Classical diffusion theory links the time rate of change of particle concentration with the divergence of the flux, which is assumed to be proportional to the gradient of concentration. Here, we used the Caputo fractional time derivative to generalize the time rate of change and the symmetric Reisz fractional space derivative to generalize the flux (in Section 3.2, we have not modified the gradient operator in the continuity equation). The fractional time derivative affects particle dynamics by introducing memory in the time derivative, f (t) * t −(1+β) + /Γ(−β), while the fractional space derivative introduces non-locality, g(x) * |x| −(1+α) /Γ(−α), [5,43]. Thus, memory and non-locality are introduced by convolution into mathematical physics when fractional order operators are employed.
Thinking of the classical diffusion problem as it might be affected by Maxwell's demon-as William Thompson did in 1875 (see pages 57-60 in Reference [51])-the demon, is viewed as a batsman in cricket who sorts out pitches/particles entering a 'control' volume based on speed and direction. In the fractional calculus version of cricket, the bowler is also a 'demon', one who can release balls with multiple delays, f (t − τ) all weighted by t −(1+β) + /Γ(−β), and from different distances, g(x − ξ), with each one weighted by |x| −(1+α) /Γ(−α).
Mathematically, the fractional time derivative can be thought of either as a convolution integral [52] or as a power law weighted sum of delay derivatives [53]. Rewriting, Equation (5), we can describe anomalous diffusion by: Physically, this changes the complete memory of ordinary integration into a fractional integral [13] (Chapter 12, page 217). In a similar way, the fractional space derivative allows particles to jump into the control volume from any distance, with a likelihood decreasing as an inverse power law of distance [50].
The CTRW phase diagram is used in this paper as the basis for comparison between fractional-derivative and power law models. It was selected for convenience and its prior appearance in the literature [27,28,30]. There are other options (e.g., fractal derivative, Lévy walks). Even within the CTRW context, different choices of the waiting time and jump increment probabilities (independent, coupled, truncated, symmetric or anti-symmetric) introduce new sub-regions (ballistic, Gaussian, sub-diffusion, super-diffusion). The CTRW study by Weeks et al. [54], which uses separate decoupled and truncated particle flight and sticking time distributions shows 5 and 6 regions, for the symmetric and anti-symmetric models, each with a different MSD diffusion exponent.
Fractional-order time and space derivatives extend the classical diffusion equation in a way that accounts for the non-Gaussian diffusion often observed in heterogeneous materials. One criticism of fractional-order models is their apparent ad hoc nature that interpolates between the well understood integer order operations of calculus. This weakness is actually a strength as illustrated in two-dimensional phase diagrams-plots whose axes represent the fractional derivative order-and typically display points corresponding to distinct diffusion propagators (Gaussian, Cauchy), lines along which specific stochastic models apply (Lévy, subordinated Brownian motion, quasi-diffusion), and regions of superand sub-diffusion where the mean squared displacement grows faster or slower than a linear function of diffusion time.
In addition to depicting the integer-and fractional-order operators of mathematical physics, anomalous diffusion phase diagrams provide a convenient way to classify different models of anomalous diffusion (continuous time random walk, fractal derivative). Specifically, each type of fractional derivative or assumed functional form of the diffusion coefficient corresponds to a difference picture of the underlying particle motion. The corresponding phase diagrams provide a portfolio of representations of anomalous diffusion. The anomalous phase cube employs lines of super-diffusion (Lévy process) and sub-diffusion (subordinated Brownian motion) to stitch together the different diffusion models displaying corresponding regions. In addition to the examples described in Section 3, the anomalous phase cube provides a framework for displaying other models and phenomena. For example, anomalous (non-linear growth of the MSD) behavior often exists only for a specific time window. This is also true for Brownian motion, which assumes that particle displacement is observed at times much greater than the time between collisions. Hence, in anomalous diffusion the MSD can be linear, sub-linear, and then linear again when plotted over 4 to 5 orders of magnitude in time [10]. Experimental examples and mathematical models typically involve one or more of the cases illustrated in Sections 3.1 and 3.2, and therefore are amenable to phase cube representation.
In addition, the complete Lévy description of anomalous diffusion (e.g., the Feller-Takayasu representation of an asymmetric stable law where α is the Lévy index and β is not the fractional or fractal time derivative, but a symmetry index) can be inserted into the base of the anomalous phase cube to replace the CTRW model. In a similar manner the fractional/fractal extension of Hamiltonian dynamics on a fractal as described by Zaslavsky [55] has the same general form as case the CTRW model (Sections 3.1 and 3.2) but for this model the order of the fractional time and space derivatives are directly connected with a fractal structure in time and distance in the underlying phase space (p, q) of the Hamiltonian system [43]. Another example, the generalized diffusion coefficient D(x, t) = D|x| a t b , discussed by Lenzi gives rise to anomalous diffusion depicted by stretched Bessel and Mittag-Leffler functions [10]. Finally, the power law diffusion coefficient model (Section 3.1), which can also be described as a stretched exponential is equivalent to with a time-dependent power law (in time) diffusion coefficient, and hence is similar to the Novikov model: D(t) = D 0 + D 1 /t 1+α , [56,57].
Additional studies brought to our attention during the review of this paper include the work modeling plasma dynamics using fractional calculus by Moradi et al. [58] and by del Castillo-Negrete [59]. Further information concerning the Kilbas-Saigo function can be found in the papers by Boudabsa and Simon [60] and by de Oliveira et al. [61]. Also, a recent mini review of anomalous diffusion in MRI is available in the work of Capuani and Palombo [62]. Furthermore, (i) details needed to evaluate the inverse Fourier transform of the Mittag-Leffler function can be found in [63], (ii) connections between anomalous diffusion and non-linear fractional diffusion equations are described in [64], and (iii) representations of the solution to fractional order anomalous diffusion equations in terms of the M-Wright function are given in [65].
The examples described in this paper do not exhaust the potential association of the phase cube with the wide variety of models for anomalous diffusion. Extensions using the Langevin, Fokker-Planck, and generalized Master equations are certainly possible, perhaps, as suggest by Bruce West under a new 'complexity' hypothesis [66]. As with all visual aids, the significance resides with the fidelity of the perspective with the experimental observations. The contribution of the phase cube is through its panoramic view of the disparate mathematical representations of anomalous diffusion, a view that can perhaps serve as a map for new explorations.

Conclusions
Three factors support the phase cube description of anomalous diffusion. First, is its successful portrayal of the variety of models (e.g., fractional, fractal, statistical) derived over the past 50 years. Second, is the visual comparison of the generalized diffusion equations with Gaussian and non-Gaussian cases as points, lines, and slices of the phase cube. Third, is the graphic presentation of regions of sub-and super-diffusion on each slice, which may suggest theoretical connections and experimental (or numerical) tests of the different models. The movement of material, such as exhaust gases, Li ions and anti-cancer agents in a porous, heterogeneous matrix is very likely to be non-Gaussian. Therefore, the continued refinement of anomalous diffusion models that capture this behavior is needed to advance our understanding of catalysis, battery storage capacity and drug delivery. Such insight can be placed into a physical and mathematical context using the anomalous phase cube.

Conflicts of Interest:
The authors declare no conflict of interest.

Appendix A. Solution to Fractional Diffusion Equation
Appendix A. 1

. Fractional Time and Space Derivatives with Power Law D(t)
In this appendix we provide the analytical solution of Equation (1), in Section 3.1, following the development presented in Garra et al. [14], and also show how a scaled form of the solution can be used to obtain the mean squared displacement. The fractional time and space diffusion equation with a time-dependent diffusion coefficient can be written as: with the diffusion coefficient given by D(t) = D ν t ν , where D ν = D β, α,ν (1 + ν) and |ν| ≤ 1.
We may use the Fourier transform, which when applied to Equation (A1) yields, The solution of Equation (A2) may be found in terms of the Kilbas-Saigo function [1], which is solution of the following Cauchy problem [14]: where f (t) = p(k, t), and λ = D ν |k| α . It is given by: Thus, the solution for the fractional diffusion Equation (A2) with the time dependent coefficient D(t) = D ν t ν , in the Fourier space, can be written as: By applying the inverse of Fourier transform in Equation (A5), we obtain We may perform the following change of variable k = D ν t ν+β 1/α k in Equation (A6), in order to write the solution as follows: It is worth mentioning that the Equation (A7) is written in a scaled form, i.e., with z = x/t ν+β α and P β,α,ν (z) = 1 2π Notice that Equation (A8) has the same form of Equation (3) in Section 3.1 of the text with ξ = (ν + β)/α.
The time dependence of the mean square displacement, which for this case corresponds to the second moment for the previous initial condition, can be obtained by using Equation (A8) with a change of variables from x to z as follows: and Now, we consider the fractional diffusion equation: This is a particular case of the Equation (A1) with ν = 0, i.e., a constant diffusion coefficient. Similar to the previous case, we can use the Fourier tranform to simplify the equation to obtain: whose solution is given in terms of the Mittag-Leffler function, By performing the inverse of Fourier transform, we obtain, which can be written as with P α,β (z) = 1 where z = x/t β/α . Note that to obtain the inverse of Fourier transform, we used the following result [15]: and some properties of the Fox functions. Another case that emerges from the Equation (A1) is obtained for β = 1, which yields: The solution for the previous equation with an initial concentration of material at the origin of an unbounded domain is a fractional Lévy motion, and may be obtained using the previous procedure by applying the Fourier transform. In this case we find the equation: whose solution is a stretched exponential, i.e., By performing the inverse of Fourier transform, we have, which may also be written as: where z = z/t (1+ν)/α . Another particular case of Equation (A1) is obtained for the case of α = 2, yielding, with the diffusion coefficient given by D(t) = D β,ν t ν , where D β,ν = D β,ν (1 + ν) and |ν| ≤ 1.
We use the Fourier transform as in the previous cases. After applying the Fourier transform to Equation (A25), we obtain: The solution for the fractional diffusion Equation (A26), in the Fourier space, can be written as: By applying the inverse of Fourier transform in Equation (A27), we have, We may perform the following change of variable k = D β,ν t ν+β 1/2 k in Equation (A6), in order to write the solution as follows: It is worth mentioning that Equation (A8) is written in a scaled form, i.e., and, consequently, where I ξ = ∞ −∞ P ξ (z)z 2 dz, z = x/t ξ , assuming I ξ exists. It is worth mentioning that the same time dependence for the second moment is obtained from the following equation: with the diffusion coefficient given by D(x) = D β,θ |x| −θ . The solution for this case may be found by using the eigenfunctions of the spatial operator present in Equation (A43). They are obtained from the Sturm-Liouville problem: subject to the conditions ψ(±∞) = 0, which implies in: In terms of these eigenfunctions, the solution can be written as: where, By using the orthogonality condition of the eigenfunctions in the fractional diffusion Equation (A43), it is possible to show that: whose solution is given by in terms of the Mittag-Leffler function as follows: Thus, for an arbitrary initial condition, we have: which can be written as follows: By using some properties and identities [15,36], Equation (A51) may be written in terms of a generalized H Fox function, i.e., with And For the initial condition P(x, 0) = δ(x), the previous solution may be simplified to: Equation (A55) has a scaled form, i.e., it can be written as follows: with P β,θ (z) = and z = x/t β/(2+θ) . Similar to the previous cases, for example, Equation (A40), we may perform the same analysis to obtain the time behavior for the second moment related to the Equation (A56). In particular, after some calculations, it is possible to show that: where I β,θ = ∞ −∞ P β,θ (z)z 2 dz, assuming I β,θ exists. Note that for β = 1, Equation (A55) may be simplified to: with the mean square displacement given by: Now, we consider the non-linear version of Equation (A32) with β = 1 and ν = 0. Following [37][38][39] and the above results, it is possible to simplify the partial differential Equation (A32) to an ordinary differential equation and to find solutions for some cases. In particular, by substituting Equation (A39) in Equation (A32), we obtain: with z = x/t ξ and ξ = 1/(θ + α + ρ). In this manner, we simplify our non-linear partial differential equation to an ordinary differential equation as mentioned before. Note that as discussed in [37,39], we will use the Riemann-Liouville operator to obtain the exact solutions. In this case, we will work with the positive axis and, later on, we will use symmetry to extend the results to the entire real axis. Particular solutions of this equation can be found as shown in [38]. Two of them are: for the case α = 1, and P α,θ,ρ (z) = N z δ (a + bz) ζ (A63) for α = 1. Here a and b are chosen so that Equation (A63) satisfies Equation (A61), and δ = (α − 1)(1 + α + θ) where N and N are normalization factors, which should be found for each case. The parameter b is chosen ±1, depending on the behavior of the solution with a = 1. In fact, for a compact behavior, we have b = −1 and for a long-tailed behavior b = 1. Note that the solutions shown for the non-linear case are given in terms of power laws and may be related to the Levy distributions as discussed in [10,[37][38][39]. Let us consider some linear cases of the Equation (A32) for ρ = 1. The first one involves a fractional time diffusion equation with a power law, time-and space-dependent diffusion coefficient as follows: with D(x, t) = D β,θ,ν t ν |x| −θ . Note that Equation (A65) combines some of the cases worked out above, i.e., Equations (A32) and (A43). The solution for this case may be found by using the Kilbas-Saigo function given in Equation (A35) and the eigenfunctions of the spatial operator present in Equation (A56). In particular, they are defined in Equation (A45).
Thus, the solution may be formally written as: for an arbitrary initial condition. By considering the initial condition P(x, 0) = δ(x), Equation (A70) may be simplified to: P(x, t) = 1 The mean square displacement was derived in Equation (A41), and is given by: where I β,θ,ν = ∞ −∞ P β,θ,ν (z)z 2 dz, assuming I β,θ,ν exists. Now, we consider the Equation (A32) with β = 1 and ρ = 1. For this case, it can be simplified to: ∂ ∂t P(x, t) = ∂ ∂x D(x, t) ∂ α ∂x α P(x, t) , with the diffusion coefficient given by D(x, t) = D α,θ,ν (1 + ν)|x| −θ t ν . By using the scaling arguments, it is possible to simplify Equation (A76) to an ordinary differential equation as follows: (1 + ν)D α,θ,ν |z| −θ d α dz α P ν,θ,α (z) = −ξzP ν,θ,α (z), where z = x/t ξ with ξ = (1 + ν)/(1 + θ + α). The solutions for this equation can be found by using the symmetries arguments used to find the solutions of the non-linear case and the Kilbas-Saigo function defined in Equation (A34). In particular, the solution for Equation (A77) is given by: where N is a normalization constant. A similar analysis has been performed in [40] by considering reaction terms. By using Equation (A78), it is possible to write the solution as follows: By using the previous equation, we obtain the time behavior for the second moment, or the mean square displacement, for this case. In this direction, it is interesting to note that Equation (A79) may be written as follows: Using Equation (A73), it is possible to show that x 2 = t 2(1+ν) 1+α+θ I α,θ,ν , by assuming that I α,θ,ν = ∞ −∞ z 2 P α,θ,ν (z)dz exists.