 # Crossover dynamics from superdiffusion to subdiffusion: Models and solutions

## Abstract

The Cattaneo or telegrapher’s equation describes the crossover from initial ballistic to normal diffusion. Here we study and survey time-fractional generalisations of this equation that are shown to produce the crossover of the mean squared displacement from superdiffusion to subdiffusion. Conditional solutions are derived in terms of Fox H-functions and the δth-order moments as well as the diffusive flux of the different models are derived. Moreover, the concept of the distribution-like is proposed as an alternative to the probability density function.

## 1 Introduction

Brownian diffusion processes [18, 81, 37] are characterised by the linear growth in time of the mean squared displacement (MSD), xϱ2(t)t. This process is governed by the coupled system of partial differential equations,

Jx,t=xϱx,t,xJx,t=tϱx,t,(1.1)

here used in dimensionless units, respectively representing the constitutive equation (Fick’s first law) and the continuity equation, where J(x, t) is the diffusive flux and ϱ(x, t) is the distribution of the diffusing substance at position x and time t. Combining both equations, we obtain the diffusion equation (Fick’s second law) tϱ(x,t)=x2ϱ(x,t). Here, ζn denotes the nth partial derivative with respect to the variable ζ. Equations (1.1) subject to the initial condition ϱ(x, 0+) = δ(x), x ∈ ℝ, in terms of the Dirac delta function (unit impulse function), and the “natural” boundary conditions lim|x|→∞ϱ(x, t) = 0, produce the famed Gaussian shape of ϱ(x, t) with the linear growth of the MSD of the diffusing particles.

In this formulation Brownian motion is overdamped, i.e., any inertial effects of the diffusing tracer particle are considered to have decayed. In the original Langevin equation formulation, the particle motion at times shorter than the relaxation time (given in terms of the particle mobility) is inertia-dominated and thus ballistic [37, 34]. One way to incorporate the crossover from initial ballistic motion 〈x2(t)〉 ∼ t2 to final diffusive motion 〈x2(t)〉 ∼ t was achieved in the discussion of heat transport  in the form of the Cattaneo or telegrapher’s equation discussed in Section 2. Today, this field has adopted more general stochastic processes to account for the ever better resolved patterns involved in heat transport, . In this paper we study and survey the mathematical properties of generalised Cattaneo-type equations, in which the integer scaling exponents of ballistic and diffusive motion are replaced by given real-valued exponents.

Already in 1926 Richardson reported a process which deviates from equations (1.1) in the form of the cubic time dependence of the relative MSD of two particles in a turbulent flow . In a wide variety of systems, the MSD of the diffusing particles follows the power-law form [6, 48, 49]

xϱ2(t)tβ,β>0,(1.2)

distinguishing subdiffusive processes with 0 < β < 1, the normal Brownian process with β = 1, superdiffusive processes with 1 < β < 2, the ballistic process with β = 2, and superballistic processes, β > 2.

To describe such anomalous diffusion processes, fractional kinetic equations as proposed by Schneider and Wyss [88, 79], Zaslavsky , and Nigmatullin , are found to be outstanding tools. The time-fractional diffusion equation, for instance, given in the natural (Caputo) form [40, 41, 42]

JC(x,t)=xϱ(x,t),xJC(x,t)=0+CDtβϱ(x,t),(1.3)

or in the modified (Riemann-Liouville) form

JRLx,t=0+RLDt1βxϱx,t,xJRLx,t=tϱx,t,(1.4)

capture well the subdiffusive regime for ϱ(x, 0+) = δ(x), x ∈ ℝ, and lim|x|→∞ϱ(x, t) = 0. Here 0+CDtβ and 0+RLDtβ, 0 < β < 1, are respectively the Caputo and Riemann-Liouville fractional derivatives (see Appendix A). Similarly, the time-fractional wave equation in the natural form

0+CDtβJCx,t=xϱx,t,xJCx,t=0+CDtβϱx,t,(1.5)

or in the modified form

0+RLDtβJRLx,t=0+RLDt1βxϱx,t,xJRLx,t=tϱx,t,(1.6)

0 < β < 1, cover the superdiffusive regime when ϱ(x, 0+) = δ(x), x ∈ ℝ, and lim|x|→∞ϱ(x, t) = 0.

Several early studies suggested that the scaling behaviour in terms of a single exponent, equation (1.2), is unable to describe many physical processes, such as the Sinai diffusion model , truncated Lévy-flights , or subdiffusive motion in velocity fields . By now, the mono-scaling behaviour has been challenged to describe the temporal variation of the MSD in biological cells [4, 29, 63, 57], in protein surface water , phospholipids and cholesterols in a lipid bilayer [31, 32, 54], and in viscoelastic flows , see also the reviews [55, 47] and the references therein. These works were the basis for a string of papers [11, 12, 13, 14, 82, 15, 72] in which fractional order derivatives in the fractional diffusion equations (1.3) and (1.4) are replaced by derivatives of distributed-order . Distributed-order time or space fractional diffusion equations are shown to be a versatile tool for the mathematical description of physical processes that become less anomalous in the course of time (accelerating subdiffusion and decelerating superdiffusion) or more anomalous in the course of time (retarding subdiffusion and accelerating superdiffusion).

In  a generalised Cattaneo-type equation was proposed, yielding the dynamic crossover behaviour

xϱ2(t)t2α,t0+,tα,t,(1.7)

where 0 < α < 1. When 0 < α < 0.5, equation (1.7) represents a retarded subdiffusion process, which can be considered a special case of the double-order time fractional diffusion equation of the natural type (β2 = 2β1) , and when 0.5 < α < 1, it represents a transition from the superdiffusive to subdiffusive regime, see also [65, 45, 2]. In the case α = 1 we recover the characteristic property of the classical Cattaneo model, see equation (2.2) below, the crossover from ballistic to normal diffusion. The reader is also referred to generalisations of diffusion, Fokker-Planck, and diffusion-wave equations by introduction of a generalised memory kernel [73, 74, 75]. However, recent numerical simulations have shown different transitions from superdiffusion to subdiffusion which are not, in all cases, covered by equation (1.7), see e.g., quasiperiodic interacting systems  and dewetting-spreading-wetting transitions of self-propelled (run-and-tumble) clustered disks in a substrate with randomly placed pinning sites .

In this work, we present different forms of generalised Cattaneo-type equations, based on single-order and distributed-order time-fractional derivatives, which provide a variety of transitions from superdiffusive to subdiffusive motion. We first outline these forms in Section 2. Conditional solutions for the distributions of the proposed models are derived in terms of Fox H-functions, and the corresponding δth-order moments are discussed in Section 3. In Section 4 we consider the diffusive flux, an important measurable quantity. Moreover, we introduce and study the so-called distribution-like, a quantity similar to the physical distribution, from which we obtain alternative formulations for the diffusion, fractional diffusion, and fractional Cattaneo-type diffusion equations. The subordination of time-fractional Cattaneo-type diffusion to the Gaussian process is also discussed. Finally, we summarise our results, and draw our conclusions in Section 5. We give a brief outline of the mathematical preliminaries in Appendix A and show the methods of deriving our solutions by solving four problems in Appendix B.

## 2 Generalised Cattaneo equations

Several representative models were introduced during the last seven decades in order to address the shortcomings of the Fourier-Fickian law in modelling the necessary finite propagation speed in heat and signal transport. One of the most prominent is the Cattaneo model, also known as Maxwell-Cattaneo or telegrapher’s equation [10, 58], as well as the phonon scattering model , and the parabolic and hyperbolic two-step models [70, 71]. The monograph  contains a comprehensive historical background on these developments, see also .

The Cattaneo equation in dimensionless form reads 

1+tJx,t=xϱx,t,xJx,t=tϱx,t,(2.1)

which encodes the characteristic ballistic-diffusive crossover property of the MSD [65, 46, 87]

xϱ2(t)t2,t0+,t,t,(2.2)

with an exact solution for ϱ(x, t) in terms of modified Bessel functions . The first fractional Cattaneo equation was proposed by Nonnenmacher and Nonnenmacher . The time-fractional Cattaneo equation (TFCE) has two familiar forms [16, 2]. The first form (natural type), uses the Caputo fractional derivative,

GCE-I: 1+0+CDtαJCx,t=xϱx,t,xJCx,t=0+CDtβϱx,t,(2.3)

with ϱ(x, 0+) = δ(x), JC(x, 0+) = 0, x ∈ ℝ, lim|x|→∞ϱ(x, t) = 0, and α, β ∈ (0, 1]. The second form of the time-fractional Cattaneo equation uses the Riemann-Liouville fractional derivative and is given by

GCE-I: 1+0+RLDtαJRLx,t=0+RLDt1βxϱx,t,xJRLx,t=tϱx,t,(2.4)

with ϱ(x, 0+) = δ(x), x ∈ ℝ, lim|x|→∞ϱ(x, t) = 0, and α, β ∈ (0, 1]. The time-fractional Cattaneo equations (2.3) and (2.4) are termed GCE-I with two fractional parameters, where the equivalence between them for the distribution ϱ(x, t) was shown in . Note that we decorate the flux in (1.3) and (2.3) with C, and in (1.4) and (2.4) with RL to point out that these equations are not necessarily equivalent for the flux. Note also that both equations (2.3) and (2.4) reduce to (2.1) in the limiting case α = β = 1, and the generalised version of the Cattaneo-type equation of  is recovered whenever α = β. If we only set α = 1 in equations (2.3) and (2.4), we get

GCE-III: 1+tJRLx,t=0+RLDt1βxϱx,t,xJRLx,t=tϱx,t,(2.5)

with ϱ(x, 0+) = δ(x), JRL(x, 0+) = 0, x ∈ ℝ, lim|x|→∞ϱ(x, t) = 0, and β ∈ (0, 1]. Alternatively,

GCE-III: 1+tJCx,t=xϱx,t,xJCx,t=0+CDtβϱx,t,(2.6)

with ϱ(x, 0+) = δ(x), JC(x, 0+) = 0, x ∈ ℝ, lim|x|→∞ϱ(x, t) = 0, and β ∈ (0, 1], which are known in the literature as GCE-III. Conversely, if we set β = 1 in (2.3) and (2.4), we obtain the GCE-IV

GCE-IV: 1+0+RLDtαJx,t=xϱx,t,xJx,t=tϱx,t,(2.7)

with ϱ(x, 0+) = δ(x), x ∈ ℝ, lim|x|→∞ϱ(x, t) = 0, and α ∈ (0, 1], where 0+RLDtα can be replaced by 0+CDtα and J(x, 0+) = 0 should be added to the initial conditions.

Equations (2.5) to (2.7) can be further generalised using the distri-buted-order derivatives to take the more general forms

1+tJRLx,t=0+RLDtp(1ν)xϱx,t,xJRLx,t=tϱx,t,(2.8)

with ϱ(x, 0+) = δ(x), JRL(x, 0+) = 0, x ∈ ℝ, lim|x|→∞ϱ(x, t) = 0, and sup(p) ⊆ [0, 1],

1+tJCx,t=xϱ(x,t),xJCx,t=0+CDtp(ν)ϱ(x,t),(2.9)

with ϱ(x, 0+) = δ(x), JC(x, 0+) = 0, x ∈ ℝ, lim|x|→∞ϱ(x, t) = 0, and sup(p) ⊆ [0, 1],

1+0+RLDtp(ν)Jx,t=xϱx,t,xJx,t=tϱx,t,(2.10)

with ϱ(x, 0+) = δ(x), x ∈ ℝ, lim|x|→∞ϱ(x, t) = 0, and sup(p) ⊆ [0, 1], where a+CDtp(ν) and a+RLDtp(ν) are the distributed-order time-fractional derivatives in the Caputo and Riemann-Liouville senses respectively, defined through

a+CDtpνft=supppνa+CDtνftdν,a+RLDtp1νf(t)=sup(p)p1νa+RLDtνftdν=supppνa+RLDt1νf(t)dν,(2.11)

and p(ν) is a probability density function, namely, p(ν) ≥ 0∀ν ∈ sup(p), and ∈ tsup(p)p(ν)dν = 1. The non-negativity of the distribution ϱ(x, t) of equations (2.3) to (2.10) was discussed in detail in . Equation (2.10) was proposed in  and analysed in  using the multi-term assumption .

In this work, we are interested in the case

p(ν)=p1δνα+p2δ(νβ),(2.12)

where 0 < α < β ≤ 1 and p1 + p2 = 1. It is salient that when α = β, or p1 = 1, p2 = 0, the distributed-order versions of the time-fractional Cattaneo-type equations (2.8) to (2.10) respectively reduce to the generalised Cattaneo equations of type III and IV, equation (2.5) to (2.7). Let us call equations (2.8) to (2.10) with (2.11) and (2.12) the double-order time-fractional Cattaneo equations of types I, II, and III. The generic form of the analytical solution for the above Cattaneo-type equations in Laplace-Fourier space, (A.5) and (A.6), can be written as

ϱ~^k,s=η(s)/sη(s)+k2,(2.13)

where η(s) is the characteristic function to be defined below, and in Laplace space as

ϱ~(x,s)=η(s)2sexp(|x|η(s)),(2.14)

with δth-order moment

xδs=20xδϱ~x,sdx=η(s)s0xδexpxη(s)dx=Γ(δ+1)sη(s)δ/2,(2.15)

and thus

xδt=Γδ+1L11sηsδ2.(2.16)

## 3 Analytical solution and mean squared displacement

Here we provide the exact solutions of the Cattaneo-type equations featured in the previous section, in addition to their corresponding MSDs.

### 3.1 Generalised Cattaneo equation of type I

In the case of GCE-I with two fractional parameters (TFCE), equations (2.3) and (2.4), the characteristic function η(s) defined in (2.13) is given by

ηs=sβsα+1.(3.1)

Thus, we have

ϱ~^GCEIk,s=sβ1sα+1sβsα+1+k2.(3.2)

Referring to Problem B.1 in appendix B the solution of equation (3.2) can be extracted from (B.1) by setting χ0,1 = 1, and γ = α + β, so that we get

ϱGCEIx,t=n=01nn!tαnϱ1x,t+tαϱ2x,t,
ϱ1x,t=14πtα+βH1,22,0x24tα+β1+αnα+β2,α+β0,1,12+n,1,ϱ2x,t=14πtα+βH1,22,0x24tα+β1+αn+αβ2,α+β0,1,12+n,1.(3.3)

Three subcases of (3.3) can be deduced: the GCE-III (α = 1), the GCE-IV (β = 1), and GCE-I with single fractional parameter (α = β). We refer to the solution of the space-time fractional Cattaneo equation obtained in  which can be considered a special case of (3.3) if the space-fractality is disregarded. Further properties on the space-time Cattaneo equation can be found in  and .

The δth-order moments can be derived by substituting (3.1) into (2.16) and inverting the Laplace transform using (A.19),

xδtGCEI=Γδ+1tα+β2δEα,α+β2δ+1δ2tα,(3.4)

where Eα,βγ(z) is the Prabhakar generalisation of the Mittag-Leffler function (A.16). An equivalent form to (3.4) can be deduced via the Mellin transform of (3.3), see (A.9), and we get

xδtGCEI=20xδϱGCEIx,tdx=2MϱGCEIx,t;xδ+1=2n=0(1)nn!tαn[Mϱ1x,t;x(δ+1)+tαM{ϱ2(x,t);x}(δ+1)]=2δπΓ1+δ2Γ1+δ2tα+β2δn=01+δ2nΓαn+α+β2δ+1tαnn!+tα1+δ2nΓαn+α+β2δ+α+1tαnn!=2δπΓ1+δ2Γ1+δ2tα+β2δ×Eα,α+β2δ+1δ2+1tα+tαEα,α+β2δ+α+1δ2+1tα.(3.5)

Both forms of the MSD, (3.4) and (3.5), should be equivalent. This equivalence of (3.4) and (3.5) can be verified by the aid of Legendre’s duplication rule and the recurrence relation (A.21). Both forms also have the same limiting behaviour, see (A.17),

xδtGCEI2δπΓ1+δ2Γ1+δ2Γ1+α+β2δtα+β2δ,t0+;2δπΓ1+δ2Γ1+δ2Γ1+β2δtβ2δ,t,(3.6)

and when δ = 2, we find the asymptotic behaviour of the MSD,

x2(t)GCEI2tα+βΓ(α+β+1),t0+,2tβΓ(β+1),t.(3.7)

The MSD (3.7) generalises the result of GCE-I with one fractional parameter (1.7). It can also represent a retarded subdiffusion process if 0 < α + β < 1, and a generic crossover from superdiffusion to subdiffusion if 1 < α + β < 2. The case 0 < α + β < 1 is analogous to the finding of Chechkin and colleagues [11, 14, 82, 15, 72]. On the other hand, the case 1 < α +β < 2 can depict many phases in  for certain pinning forces at low disk activity. The probability density (3.3) and the MSD (3.7) are presented in Figures 1a and 1b respectively. In Figure 1a, the subdiffusive transition regime in the short-time limiting included in the Cattaneo-type equation is clear in the blue curve, where 0 < α +β < 1. In the red and green curves, the superdiffusion process is dominant since 1 < α + β < 2. In Figure 1b the blue curve, i.e., the case corresponding to α = 0.5, β = 0.2, the retarded subdiffusion process is salient in the course of time. In the red and green curves of Figure 1b, which consider the choices α = 0.5, β = 0.8 and α = 0.8, β = 0.8, namely 1 < α + β < 2, respectively, the crossover from superdiffusion to subdiffusion in both curves is obvious. Figure 1

(a) Spatial evolution of the probability density ϱGCE-I of GCE-I at short time t = 0.1 for different fractional parameters. (b) Temporal evolution of the MSD of GCE-I. The normal and ballistic diffusion limits are shown in the light grey curves.

### Remark 3.1

In spite of the generality of the limiting behaviour (3.7) covering a wide range of transitions over the subdiffusive regime or from the superdiffusive to the subdiffusive regime, it is still constrained by the fact that the parameter β occurs in both the short-time and the long-time behaviours. In essence, the presence of the parameter β in the short-time behaviour causes an exemplary transition such as

xϱ2(t)t1.8,t0+,t0.3,t,

impossible to be captured by the GCE-I (and thus its special cases; GCE-III and GCE-IV) with two fractional parameters.

### Remark 3.2

It is noteworthy mentioning that GCE-I (3.2) can be written in terms of the generalised diffusion-wave equation with two power-law memory kernels , namely,

0tξtττ2ϱx,τdτ=x2ϱx,t,(3.8)

where ξ(t)=t1βΓ(2β)+t1αβΓ(2αβ),, ϱ(x, 0+) = δ(x), tϱ(x, 0+) = 0, lim|x|→∞ϱ(x, t) = 0, and α, β ∈ (0, 1], and in terms of the generalised diffusion-wave equation with regularised Prabhakar derivative, see (A.4),

0+PDα,1,t1,α+βϱx,t=x2ϱx,t,(3.9)

where ϱ(x, 0+) = δ(x), tϱ(x, 0+) = 0 lim|x|→∞ϱ(x, t) = 0, and provided that 0 < α, β < 1 and 1 < α + β < 2.

### Remark 3.3

The normalisation of (3.3) can be verified through the Mellin transform,

ϱGCEIx,tdx=20ϱGCEIx,tdx=2MϱGCEIx,t;x1=2n=01nn!tαnMϱ1x,t;x1+tαMϱ2x,t;x1=n=01ntαnΓ1+αn+tαn+1Γ1+αn+1=1.

### 3.2 Double-order Cattaneo equation of type I

The double-order time-fractional Cattaneo equation of type I (DTFCE-I) (2.8), (2.11) and (2.12), which is distinguished from the double-order time-fractional diffusion equation of the modified form by the presence of a wave motion in the short-time behaviour, possesses the characteristic function

ηs=s+1p1sα+p2sβ,(3.10)

where 0 < α < β ≤ 1 and p1 + p2 = 1. For sufficiently large values of time, it can be shown that η(s) ∼ (p1sα + p2sβ)−1. The solution of DTFCE-I in the Laplace-Fourier space is thus given by (2.13) and (3.10),

ϱ~^Ik,s=sα+sα1sα+1+sα+k2p1+p2sαβ,(3.11)

which is a special case of Problem B.2, when χ0,1 = 1, and γ = α + 1. Then, one can deduce that

ϱIx,t=n=0tnn!m=0nnmp2p1tβα1mϱ1x,t+tϱ2x,t,

where

ϱ1x,t=14πp1tα+1×H2,32,1x24p1tα+112m,1;12+n+βα1mα2,α+10,1,12+nm,1;12,1,ϱ2x,t=14πp1tα+1×H2,32,1x24p1tα+112m,1;32+n+βα1mα2,α+10,1,12+nm,1;12,1(3.12)

where p1p2. It can be verified that the distribution (3.12) is normalised as follows

ϱIx,tdx=2MϱIx,t;x1=n=0tnn!m=0nnmp2p1tβα1m×Γ1+nmΓmΓ0Γ1+n+βα1m+Γ1+nmΓm×tΓ0Γ2+n+βα1m.

Since 1/Γ(0) = 0, the above series has non-zero value only at m = 0, then

ϱIx,tdx=n=0tn1Γ1+n+tΓ2+n=1.

The δ th- order moments of DTFCE-I can be determined from equations (2.16) and (3.10),

xδsI=Γ(δ+1)s+1δ/2sp1sα+p2sβδ/2=Γ(δ+1)p1δ/2sβ+12δ11+1sδ/2sβα+p2p1δ/2=Γ(δ+1)p1δ/2n=01nδ/2nn!sβ+12δ1nsβα+p2p1δ/2,

where the expansion 1+xr=n=01nrnn!xn,r ∈ ℝ+, |x| < 1, was employed to 1+1sδ/2 for the values Re(s) > 1. Therefore,

xδtI=Γ(δ+1)p1δ/2tα+12δn=0tnδ/2nn!Eβα,α+12δ+n+1δ/2p2p1tβα.(3.13)

Note that the δth-order moments (3.13) represents the exact formula at short and intermediate times time t (since Re(s) > 1).

The reader can analogously deduce a complementary exact formula for long times by assuming Re(s) < 1, and expanding then (1 + s)δ/2 instead of (1 + 1/s)δ/2. The asymptotic behaviour of all the δth-order moments is

xδtIΓ1+δp1δ/2Γ1+α+12δtα+12δ,t0+,Γ1+δp2δ/2Γ1+βδ2tβδ2,t,(3.14)

with MSD limiting behaviours

x2(t)I2p1tα+1Γ(α+2),t0+,2p2tβΓ(β+1),t,(3.15)

for 0 < α < β ≤ 1, which can be generally viewed as a modified transition from superdiffusion to subdiffusion. It is also included in the wide range of crossovers covered by the generic form (3.7), however, the temporal path on which they evolve is distinct. In the case 0.9 < β ≤ 1, it represents processes switching from superdiffusion to normal diffusion, whereas the case 0 < α ≤ 0.1 describes processes crossing over from normal diffusion to subdiffusion. This type of transition can be found in the transport of electrons in an interacting system subject to quasiperiodic potential , specifically, when different interacting strengths U are considered. In Figures 2a and 2b the graphical representations of the distribution (3.12) and MSD of DTFCE-I, 〈x2(t)〉I, are shown. In the course of time, the transition from superdiffusion at small times to subdiffusion at longer times is notable in Figure 2a. The apparent warp in the MSD temporal path of Figure 2b seems almost like the result of changing the interacting strength with fixed quasiperiodic potential amplitude, see Figure 2 of . The curvature of this warp can be increased by decreasing the second constant, e.g., p1 = 0.99 and p1 = 0.01. Figure 2

(a) Spatial distribution of DTFCE-I at different instants of time for fractional parameters α = 0.3, β = 0.8 and constants p1 = p2 = 0.5. (b) MSD of the distribution DTFCE-I versus time for fractional parameters α = 0.3, β = 0.8 and different values of the model constants; p1 = 0.9, p2 = 0.1 (blue curve) and p1 = 0.1, p2 = 0.9 (red curve). The case p1 = p2 = 0.5, lying between blue and red curve, is omitted.

### 3.3 Double-order Cattaneo equation of type II

We consider here the double-order time-fractional Cattaneo equation of type II (DTFCE-II), (2.9), (2.11) and (2.12). It reduces to the double-order time-fractional diffusion equation of the natural form in the relatively long temporal domain which typifies subdiffusive processes that become more anomalous with time progress. The characteristic function η(s) is given as

ηs=s+1p1sα+p2sβ,(3.16)

where 0 < α < β ≤ 1 and p1 + p2 = 1. It appears from (3.16) that η(s) ∼ p1sα + p2sβ for relatively large values of time which coincides the characteristic function of the double-order time-fractional diffusion equation of the natural form. In the long time limit, η(s) ∼ p1sα + p2sβp1sα proviso p1p2. Upon setting η(s) ∼ p1sα in (2.13), one can get 〈x2(t)〉 ∼ tα, and vice versa in the short time limit. This result was analytically shown in [11, 82, 15, 72]. The distribution of DTFCE-II is a direct consequence of using of (3.16) in (2.13),

ϱ~^IIk,s=p1sα+p2sβ+p1sα1+p2sβ1p1sα+1+p2sβ+1+p1sα+p2sβ+k2.(3.17)

A similar case to (3.17) has been processed in Problem B.3. Thus, setting χ0,1 = 1, γ1 = α + 1 and γ2 = β + 1 in (B.3) leads to the solution of DTFCE-II,

ϱIIx,t=n=0tβαnn!m=0nnmp1p2nmtαβ1m×=0mmp1p2tβα××p1p2tβαϱ1x,t+tϱ3x,t+ϱ2x,t+tϱ4x,t,ϱ1x,t=14πtβ+1p2×H1,22,0p2x24tβ+112+m+βαnm++1β2,β+10,1,12+n,1,ϱ2x,t=14πtβ+1p2×H1,22,0p2x24tβ+112+m+βαnm+β2,β+10,1,12+n,1,
ϱ3x,t=14πtβ+1p2×H1,22,0p2x24tβ+132+m+βαnm++1β2,β+10,1,12+n,1,ϱ4x,t=14πtβ+1p2×H1,22,0p2x24tβ+132+m+βαnm+β2,β+10,1,12+n,1,(3.18)

where p1 < p2. The distribution (3.18) is normalised over ℝ. Indeed,

ϱIIx,tdx=n=0tβαnm=0nnmp1p2nmtαβ1m×=0mmp1p2tβα×p1p2tβαΓ1+m+βαnm++1+p1p2tβα+1Γ2+m+βαnm++1+1Γ1+m+βαnm++tΓ2+m+βαnm+=1,

which can be attained by expansion for n = 0, 1, 2, 3. We find that all terms cancel each other except the first term of 1/Γ(1 + m + (βα)(nm + )).

Following similar procedures to the δth-order moments of DTFCE-I, we obtain

xδtII=Γ(δ+1)tβ+12δp2δ/2n=0tnδ/2nn!Eβα,β+12δ+n+1δ/2p1p2tβα,(3.19)

for short- and intermediate- times, with the limiting behaviours

xδtIIΓ1+δp2δ/2Γ1+β+12δtβ+12δ,t0+;Γ1+δp1δ/2Γ1+αδ2tαδ2,t,(3.20)

and MSD asymptotes

x2(t)II2tβ+1p2Γ(β+2),t0+,2tαp1Γ(α+1),t,(3.21)

where 0 < α < β ≤ 1 and p1 + p2 = 1. The limiting behaviour (3.21) proves that the DTFCE-II provides a second modified version of the transitions from superdiffusion to subdiffusion. On the other hand, for 0.9 < β ≤ 1 and 0 < α ≤ 0.9, we have a crossover from ballistic behaviour to subdiffusion, which essentially characterises the DTFCE-II. This crossover could not be performed by the generic behaviour (3.7) nor the specific behaviour (3.15). The distribution (3.18) and the MSD of the DTFCE-II are drawn in figures 3a and 3b, respectively. In accord with the condition (p1 < p2), see (3.18) below, we set p1 = 0.4, p2 = 0.6 in Figure 3a. Also, we choose α = 0.3, β = 0.8 for closer comparisons with Figure 2a. The DTFCE-II curve remarkably records values for the distribution ϱII(x, t) higher than those recorded in figure 2a for DTFCE-I. In the course of time, the transition from superdiffusion at small times to subdiffusion at long times is notable here as well. The presence of p1 and p2, where p1 + p2 = 1, in the MSD covers a significant range in the vicinity of the points tβ+1 and tα, see Figure 3b. Figure 3

(a) Spatial distribution of DTFCE-II at different instants of time for fractional parameters α = 0.3, β = 0.8 and constants p1 = 0.4, p2 = 0.6. (b) MSD 〈x2II of DTFCE-II as a function of time for fractional parameters α = 0.3, β = 0.8 and different values of model constants; p1 = 0.9, p2 = 0.1 (blue curve) and p1 = 0.1, p2 = 0.9 (red curve).

### 3.4 Double-order Cattaneo equation of type III

In this subsection, we study the effect of distributed-order fractional derivatives on the generalised Cattaneo equation of type IV. The resulting model is termed double-order time-fractional Cattaneo equation of type III (DTFCE-III), having the characteristic function

ηs=s1+p1sα+p2sβ,(3.22)

where 0 < α < β ≤ 1 and p1 + p2 = 1. The long time behaviour of (3.22) is distinctly η(s) ∼ s, which in combination with (2.13) conducts as, the most prominent, the Gaussian distribution, i.e. 〈x2(t)〉 ∼ t. Conversely for short times (3.22) behaves as η(s) ∼ s(p1sα + p2sβ). This form resembles the characteristic function of the double-order time-fractional wave equation of natural type [2, 75, 23], which offers decelerating superdiffusive behaviour in the intermediate short-time domain, see Section 3.3 in . Thus, the study of the characteristic function suggests the progress 〈x2(t)〉 ∼ tβ+1t in the course of time.

The characteristic function (3.22) with (2.13) produces the (ks) form

ϱ~^IIIk,s=1+p1sα+p2sβp1sα+1+p2sβ+1+s+k2,(3.23)

which has been considered in a generic form in Problem B.4. Consequently, invoking equation (B.4) after setting χ0,1 = 1, γ1 = α + 1 and γ2 = β + 1, helps in getting the exact (conditional) form

ϱIIIx,t=n=01nn!tβp2nm=0nnmp1tαm×tβp2ϱ1x,t+p1p2tβαϱ2x,t+ϱ3x,t,ϱ1x,t=14πtβ+1/p2×H1,22,0x24tβ+1/p212+βnαm+β2,β+10,1,12+n,1,
ϱ2x,t=14πtβ+1/p2×H1,22,0x24tβ+1/p212+βnαm+β2α,β+10,1,12+n,1,ϱ3x,t=14πtβ+1/p2×H1,22,0x24tβ+1/p212+βnαmβ2,β+10,1,12+n,1,

where p1 ∈ [0, 0.3], p2 ∈ [0.7, 1], p1 + p2 = 1, β ∈ [0.4, 1], 0 < α < β. The normalisation of (3.24) can be checked through the expansion of the series

ϱIIIx,tdx=n=01nm=0nnmp1mp2n+1tβn+1αmΓ1+βn+1αm+p1m+1p2n+1tβn+1αm+1Γ1+βn+1αm+1+p1mp2ntβnαmΓ1+βnαm=1.

The scaling properties for the DTFCE-III can be obtained through its δth-order moments

xδtIII=Γ(δ+1)tβ+12δp2δ/2n=0tβnn!δ/2np2n×Eβα,β+12δ+βn+1n+δ/2p1p2tβα,(3.24)

for short- and intermediate-time values, with short- and long-time behaviours

xδtIIIΓ1+δp2δ/2Γ1+β+12δtβ+12δ,t0+,Γ1+δΓ1+δ2tδ2,t.(3.25)

Hence, the MSD of the DTFCE-III reads

x2(t)III2tβ+1p2Γ(β+2),t0+,2t,t,(3.26)

which confirms the predictions of the characteristic functions (3.22). The domination of the largest exponent in the short-time domain is evident. The behaviour (3.26) represents generally a transition from superdiffusion to strictly normal diffusion. What distinguishes this model from the GCE-IV, in essence, is the presence of another smaller fractional parameter α, which unfortunately disappears from (3.26), and the constant p2 in the short-time limit. This parameter, however, seems to be activated in the short-time domain, which may lead to a retardation in this time scale. In Figure 4a we represent graphically the probability density function ϱIII(x, t) of DTFCE-III at different values of time for the fractional parameters α = 0.3, β = 0.8 and model constants p1 = 0.3, p2 = 0.7, depending on the conditions imposed on the closed-form solution (3.24), and its MSD for fractional parameters α = 0.3, β = 0.8 and different values of model constants in Figure 4b. Due to the resemblance between the MSD of DTFCE-II and DTFCE-III in the short-time behaviour, they record higher values in this time domain compared with the DTFCE-I where the smallest exponent α dominates. Figure 4

(a) Spatial distribution of DTFCE-III at different instants of time for fractional parameters α = 0.3, β = 0.8 and constants p1 = 0.3, p2 = 0.7. (b) MSD 〈x2III of the DTFCE-III as a function of time for fractional parameters α = 0.3, β = 0.8 and different values of model constants; p1 = 0.9, p2 = 0.1 (blue curve) and p1 = 0.1, p2 = 0.9 (red curve).

## 4 Diffusive flux and distribution-like formulation

In this section we provide an alternative formulation for the Green’s function of diffusion (1.1), time-fractional diffusion (1.3) and (1.4), and time-fractional Cattaneo-type diffusion (2.3) and (2.4), by introducing a field variable that behaves like the distribution of the diffusing substance through having the same MSD of ϱ(x, t). The exact formulas for the diffusive flux are also derived.

### 4.1 Normal diffusion

At first, let us start with the normal diffusion equation (1.1) which can be converted using the transform technique to

ϱ~^k,s=1s+k2,J~^k,s=ıks+k2,(4.1)

where ϱ(x, t) and J(x, t) have the same partial differential (diffusion) equation x2fx,t=tfx,t with different initial conditions for the concentration ϱ(x, 0+) = δ(x) and diffusive flux J(x, 0+) = −u1(x), and u1(x) is the unit doublet function, see Section A.4 in Appendix A. Rewriting (4.1) in the form

ϱ~x,s=12sexpxs,J~x,s=x2xexpxs,(4.2)

ensures the non-negativity of ϱ(x, t) and J(x, t). It is clear that ϱ͠(x, λ) is a completely monotone function for λ > 0, see Section A.5 in Appendix A, as a product of two completely monotone functions, and ϱ~x,sdx=1/s which prove that ϱ(x, t) is a probability density function (PDF). On the other hand, evidently, the flux J(x, t) changes its sign according to the sign of x, where sign(x) = |x|/x. Physically, the three-dimensional flux J(x, t), x = (x1, x2, x2), is a vector quantity expressing the rate of movement of diffusing particles through a unit area. So, it is not expected that J(x, t) is a PDF. Furthermore, if we retrieve (4.2) in the space-time domain, we get (see also )

ϱx,t=14πtexpx24t,Jx,t=x2t4πtexpx24t,(4.3)

which shows the Gaussian behaviour for particle distribution with xϱ2(t)t and positive skewness for the diffusive flux in the positive half-real axis. The diffusive flux reverses its sign in the negative half-real axis because of the presence of x. In other words, the overall flux has to vanish, i.e. Jx,tdx=0.

Let us consider the quantity

Qx,t=xJx,t=x22t4πtexpx24t=x22tϱx,t,(4.4)

which takes the dimensional form Q(x,t)=xDJ(x,t)=xxϱ(x,t), where the prime symbols stand for dimensional variables and D is the diffusion coefficient. Hereby, it has the same dimension of ϱ′(x′, t′), i.e. [Q−] = [ϱ−] = mol/cm3. Moreover, its Laplace transform is given as

Q~x,s=12xexpxs.(4.5)

We note that (x, λ) is a completely monotone function for λ > 0, see definitions (A.1) and (A.4), and Q~x,sdx=1/s, thereby, Q(x, t) is a PDF. Interestingly, the MSD of Q(x, t) is given by

xQ2(t)=6t=3xϱ2(t),(4.6)

namely, Q encodes a behaviour similar to the classical distribution ϱ(x, t) although it possesses a different form, and we call it the “distribution-like”.

Using the concept of this distribution-like, we can summarise the main result of this subsection as follows: the classical Fourier-Fick equation (1.1) subject to a Dirac delta initial distribution and unit doublet for the initial flux can be rewritten in the sense of the distribution-like as

Qx,t=xJx,t,xx2Jx,t=tQx,t,(4.7)

subject to the initial conditions Q(x, 0+) = − U1(x) and diffusive flux J(x, 0+) = −u1(x), and boundary condition lim|x|→∞Q(x, t) = 0, where U1(x) is the modified unit doublet function, see equations (A.27) and (A.29).

### 4.2 Time-fractional diffusion

To describe the anomalous diffusive behaviour xϱ2(t)tβ, β > 0, considerable effort has been invested on the mathematical physics side. Owing to the fact that non-Brownian diffusion is non-universal, today a wide variety of anomalous stochastic processes are known  Specifically, fractional diffusion equations were put forward [48, 88, 90, 61], see [40, 43, 59, 28, 27]. Following such ideas, crossovers of the MSD are now frequently described in terms of distributed-order fractional diffusion equations, see the summaries [82, 72, 73]. We here consider single-parameter time-fractional diffusion equations.

Among the single-parameter forms of the time-fractional diffusion equations we pay attention to the two forms (1.3) and (1.4) which result in the same behaviour of the MSD (1.2), 0 < β ≤ 1, however, they produce two different partial integrodifferential equations for the diffusive flux,

x2JRLx,t=0+RLDtβJRLx,t=t0+RLIt1βJRLx,t,(4.8)

which requires the initial condition limt0+0+RLIt1βJRLx,t=u1(x), whereas the other reads

x2JCx,t=0+CDtβJCx,t,(4.9)

that works well with the conventional initial condition J(x, 0+) = − u1(x).

### Remark 4.1

It is worth remarking that the RL form (1.4) can be recovered from the form

0+RLIt1βJRLx,t=xϱx,t,xJRLx,t=tϱx,t,(4.10)

where a+RLDtβ is the left inverse of a+RLItβ,, i.e. a+RLDtβa+RLItβf(t)=ft, with β > 0. Moreover, the form (4.10) produces the well-known Caputo form x2ϱx,t=0+CDtβϱx,t for the distribution and the integrodifferential equation (4.8) for the flux.

Let us now turn our interest to the properties offered by the forms (1.3) and (4.10), where we disregard equation (1.4) because of its equivalence with (59). They possess the transformed distribution

ϱ~^k,s=sβ1sβ+k2,(4.11)

which remarkably can be retrieved from (3.2) if we consider the long time behaviour of the characteristic function (3.1), η(s) = sβ(sα + 1) ∼ sβ, with equation (2.13), and the transformed fluxes

J~^Ck,s=sβ1J~^RLk,s=ιksβ1sβ+k2.(4.12)

By taking the inverse Fourier transform of (4.11) and (4.12) we have

ϱ~(x,s)=sβ2sexp|x|sβ,(4.13)
J~C(x,s)=sβ1J~RL(x,s)=|x|2xsβ1exp|x|sβ.(4.14)

It was established that the distribution ϱ(x, t) is a PDF  with MSD xϱ2(t)tβ, see also  for the subordination technique and  for the Bernstein functions technique. It appears also from (4.14) that both fluxes include the sign function, sign(x) = |x|/x. We introduce the distribution-like for the time-fractional diffusion equation as

Qx,t=0+RLIt1βxJRLx,t=xJCx,t,(4.15)

with its transform

Q~x,s=12sβ1xexpxsβ.(4.16)

One can easily check that (x, λ) ∈ 𝓒𝓜𝓕 for λ > 0, because it is a product of two completely monotone functions as equation (4.16) explains. Also, Q~x,sdx=1/s which proves the normalisation of Q(x, t). Thus, the distribution-like of the time-fractional diffusion equations considered in this work is a probability density function. It can be also shown that the MSD of the distribution-like is

xQ2(t)=6tβΓβ+1=3xϱ2(t),(4.17)

which likens its counterpart in the normal diffusion case, equation (4.6).

Utilising the definition of distribution-like (4.15), we can rewrite expressions (1.3) and (4.10) in the form

Qx,t=0+RLIt1βxJRLx,t,xx2JRLx,t=tQx,t,(4.18)

subject to the initial conditions Q(x, 0+) = − U1(x), the diffusive flux limt0+0+RLIt1βJRLx,t=u1(x), and

Qx,t=xJCx,t,xx2JCx,t=0+CDtβQx,t,(4.19)

xx2JCx,t=0+CDtβQx,t subject to the initial conditions Q(x, 0+) = − U1(x) and the diffusive flux JC(x, 0+) = − u1(x) Both equations are subjected to the boundary condition lim|x|→∞Q(x, t) = 0. In the case β = 1, equations (4.18) and (4.19) reduce to (4.7). We have shifted the exact formulas for Q(x, t), JRL(x, t) and JC(x, t) to the next subsection to avoid reiteration.

### 4.3 Time-fractional finite-speed diffusion

The two forms of the time-fractional Cattaneo equation describing the finite-speed diffusion were introduced in equations (2.3) and (2.4) in the sense of Caputo and Riemann-Liouville fractional derivatives respectively. It was shown that the distribution ϱ(x, t) resulting from (2.3) exactly coincides with the corresponding one based on (2.4), see . Nevertheless, this coincidence is lost when the diffusive flux is considered. In point of fact, in order to mathematically utilise Cattaneo-type equations in studying the itinerant (non-scaling) anomalous behaviour, the initial conditions should be imposed on the rate of the diffusive flux, not the diffusing flux itself. Loosely speaking, equation (2.1) should be subject to t J(x, 0+) = − u1(x), equation (2.3) subject to limt0+0+CDtαJCx,t = − u1(x), and equation (2.4) subject to limt0+0+RLDtα0+RLIt1βJRLx,t = − u1(x) analogously to its counterpart in the time-fractional diffusion equation. These conditions are all leading to zero initial flux or limt0+0+RLIt1βJRLx,t=0. In analogue to remark 4.1, we suggest another form for the time fractional Cattaneo equation.

### Remark 4.2

The following Cattaneo-type equation

1+0+RLDtα0+RLIt1βJRLx,t=xϱx,t,xJRLx,t=tϱx,t,(4.20)

with ϱ(x, 0+) = δ(x), limt0+0+RLIt1βJRLx,t = 0, x ∈ ℝ, as well as lim|x|→∞ϱ(x, t) = 0, and α, β ∈ (0, 1], yields the same distribution and diffusive flux provided by (2.4). The condition limt0+0+RLIt1βJRLx,t=0 guarantees the commutation 0+RLDt1β0+RLDtαf(t)=0+RLDtα0+RLDt1βf(t). Alternatively, it can be written as limt0+0+RLDtα0+RLIt1βJRLx,t=u1(x).

Employing the transform technique to equations (2.3) and (4.20) we obtain

J~^Ck,s=sβ1J~^RLk,s=ιksβ1sα+β+sβ+k2,(4.21)

in the Laplace-Fourier domain, and

J~Ck,s=sβ1J~RLx,s=x2xsβ1expxsβ1+sα(4.22)

in the Laplace domain. Accordingly, we take the distribution-like of the time-fractional Cattaneo equation in the form

Qx,t=1+0+RLDtα0+RLIt1βxJRLx,t=1+0+CDtαxJCx,t,(4.23)

whose Laplace transform is

Q~x,s=12sβ11+sαxexpxsβ1+sα.(4.24)

We now claim that Q(x, t) given by (4.23) is a probability density function if α + β ≤ 1, where α, β ∈ (0, 1]. To justify this claim, we first prove the normalisation of Q(x, t) by integrating (4.24) over ℝ, Q~x,sdx=1/s, namely, Q(x, t) is normalised. Secondly, we study the non-negativity of Q(x, t) by rewriting (4.24) as

Q~x,λ=12f1λf2λf3λ;x,λ>0,f1λ=λβ1,f2λ=λα+β1,f3λ;x=xexpxλβ121+λα12.(4.25)

Obviously, f1(λ) is completely monotone for β ∈ (0, 1] whilst f2(λ) is completely monotone for α + β ∈ (0, 1], refer to definitions (A.1) and (A.3). On the other hand, we have that λβ and 1 + λα are complete Bernstein functions α, β ∈ (0, 1] as definition (A.4) reports, note that 1 ∈ 𝓢𝓕 ∩ 𝓒𝓑𝓕. Therefore, λβ121+λα12CBF, and f3(λ;|x|) ∈ 𝓒𝓜𝓕. This proves that (x, λ) is completely monotone as a product of three complete monotone functions. Hence Q(x, t) ≥ 0 if α + β ≤ 1 and α, β ∈ (0, 1], which completes the proof of our claim. On the other hand, this result can be reinforced by verifying that the process following the distribution-like (4.23) is subordinated to the normal (Gaussian) diffusion provided that α + β ≤ 1 and α, β ∈ (0, 1], using the following criterion , see also . Firstly, equation (4.24) can be rewritten in the form:

Q~x,s=x2πxsβ11+sαsα+β+sβ+k2eιkxdk=x2πxsβ11+sα0esβ1+sα+k2udueιkxdk=x2πx0ek2uG~u,sdueιkxdk=0x22u4πuex24uG~u,sdu,(4.26)

where (u, s), given by

G~u,s=sβ11+sαesβ1+sαu,(4.27)

is the Laplace transform of the function G(u, t) providing the subordination transformation from the time scale t to the time scale u. It remains to prove that G(u, t) is a PDF with respect to u > 0 for any t > 0. Firstly, we note that G(u, t) is normalised with respect to u. Indeed, since 0(u, s)du = 1/s, then 0G(u, t)du = 1. Secondly. The non-negativity of G(u, t) can be proved by verifying that its Laplace transform (4.27) is completely monotone. By rewriting (u, s) as

G~u,λ=f1λf2λg1u,λg2u,λ,λ>0f1λ=λβ1,f2λ=λα+β1,g1u,λ=euλβ,g2u,λ=euλα+β,(4.28)

one concludes that f1(λ) and g1(u, λ) are completely monotone functions for λ > 0 and β ∈ (0, 1], whilst f2(λ) and g2(u, λ) require that λ > 0, α, β ∈ (0, 1] and α + β ≤ 1 to be complete monotone functions, where subsection A.5 of Appendix A has been invoked repeatedly. Therefore, (u, λ), as a product of complete monotone functions, is in turn a completely monotone, which proves the non-negativity of G(u, t). In view of equation (4.4) and the Laplace inversion of (4.26), we can state that the time-fractional Cattaneo equation with two fractional derivatives of different orders α, β ∈ (0, 1] represents a mathematical model for a stochastic process subordinated to Gaussian diffusion provided that the fractional derivative orders satisfy α + β ≤ 1. It is worth mentioning that the same criterion can be implemented for the distribution ϱ(x, t) using the same transformation G(u, t).

The mathematical definition of the distribution-like (4.23) enables us to replace the time-fractional Cattaneo equation in the RL-sense (2.4) and (4.20) by

Qx,t=1+χ0,10+RLDtα0+RLIt1βxJRLx,t,xx2JRLx,t=tQx,t(4.29)

subject to the initial conditions Q(x, 0+) = − U1(x) and diffusive flux limt0+0+RLDtα0+RLIt1βJRLx,t=u1(x), and the time-fractional Cattaneo equation in the C-sense (2.3) by

Qx,t=1+χ0,10+CDtαxJCx,t,xx2JCx,t=0+CDtβQx,t(4.30)

subject to the initial conditions Q(x, 0+) = − U1(x) and diffusive flux limt0+0+CDtαJCx,t=u1(x),χ0,10,1, and both equations are subject to the boundary condition lim|x|→∞Q(x, t) = 0. When χ0,1 = 0 and limα0+0+CDtα=limα0+0+RLDtα=1, equations (4.29) and (4.30) reduce to (4.18) and (4.19), respectively.

Lastly, we complete this section by providing analytical solutions for the diffusive fluxes JRL(x, t) and JC(x, t), and the distribution-like Q(x, t) for TFCE. The diffusive fluxes in RL and C sense are given by

JRLx,t=xφRLx,t,JCx,t=xφCx,t,(4.31)

where φRL(x, t) and φC(x, t) are given by

φ~^RLk,s=1sγ+χ0,1sβ+k2,φ~^Ck,s=sγ1sγ+χ0,1sβ+k2,(4.32)

where χ0,1 ∈ {0, 1}, and γβ such that γ = β if χ0,1 = 0 and γ > β if χ0,1 = 1. By following the same procedures of solving Problem B.1, one derives

φRLx,t=tγ14πtγn=0χ0,1nn!tγβnH1,22,0x24tγγβn+γ2,γ0,1,12+n,1.(4.33)

Then, substituting (4.33) into (4.31) and using (A.13), we arrive at

JRLx,t=tγ1xπtγn=0χ0,1nn!tγβnH1,22,0x24tγγβn+γ2,γ1,1,12+n,1,(4.34)

where the relation Γ(z/2)Γ(1 − z)/Γ (−z) = −2Γ(1 + z/2) has been tacitly utilised. We have two special cases included in equation (4.34): First, when χ0,1 = 0, equation (4.34) presents the flux of time-fractional diffusion equation (1.4) and (4.10),

JRLTFDEx,t=tβ1xπtβH1,22,0x24tββ2,β1,1,12,1,(4.35)

and when β = 1, we get the physical flux of the normal diffusion (1.1),

JNDx,t=1xπtH1,22,0x24t12,11,1,12,1=1xπtH0,11,0x24t1,1=x2t4πtexpx24t,

which coincides with equation (4.3). Here, we have used the Mellin transform with equation (A.15). Secondly, when χ0,1 = 1, then equation (4.34) yields the flux of time-fractional Cattaneo equation (2.4) and (4.20),

JRLTFCE(x,t)=tα+β1xπtα+βn=01nn!tαn×H1,22,0x24tα+βαn+α+β2,α+β1,1,12+n,1.(4.36)

Analogously, we can deduce an exact generic formula for JC(x, t),

JCx,t=1xπtγn=0χ0,1nn!tγβn+1×H1,22,0x24tγ1+γβn+1γ2,γ1,1,12+n,1,(4.37)

which reduces to the flux of TFDE with time-fractality in the continuity equation (1.3) when χ0,1 = 0,

JCTFDEx,t=1xπtβH1,22,0x24tβ1β2,β1,1,12,1.(4.38)

Note that JCTFDEx,t=JRLTFDEx,t=JNDx,t in the case β = 1. Conversely, when χ0,1 = 1, equation (4.37) reduces to the flux of TFCE with time-fractality in the continuity equation (2.3),

JCTFCEx,t=tαxπtα+βn=0tαnn!×H1,22,0x24tα+β1+αn+αβ2,α+β1,1,12+n,1.(4.39)

The diffusive fluxes in the RL sense (without fractality in the continuity equation) (4.36), and in the C sense (with temporal fractality in the continuity equation) (4.39) are drawn in figures 5a and 5b, respectively. Figure 5

(a) Diffusive flux JRLTFCE(without fractality in the continuity equation) at dimensionless time t = 0.1, and different fractional parameters. (b) Diffusive flux JCTFCE(with fractality in the continuity equation) at dimensionless time t = 0.1, and different fractional parameters.

Finally, the distribution-like can be directly derived from the unified nondimensional relation

Qx,t=xxϱx,t(4.40)

together with the universal distribution (B.1), we get, upon using (A.12) and (A.13), the Mellin transform and the known relation of the Gamma function, Γ(z/2)Γ(1 − z)/Γ(−z) = −2Γ(1 + z/2),

Qx,t=n=0χ0,1nn!tγβnQ1x,t+χ0,1tγβQ2x,t,Q1x,t=1πtγH1,22,0x24tγ1+γβnγ2,γ1,1,12+n,1,Q2x,t=1πtγH1,22,0x24tγ1+γβn+1γ2,γ1,1,12+n,1.(4.41)

Again, the last generic formula for the distribution-like contains two special cases of interest: the first case is χ0,1 = 0, that gives the distribution-like for TFDE (4.18) and (4.19)

QTFDEx,t=1πtβH1,22,0x24tβ1β2,β1,1,12,1,(4.42)

which reduces to (4.4) in the limit β = 1. Meanwhile, the second case, χ0,1 = 1, gives the distribution-like for TFCE (4.29) and (4.30)

QTFCEx,t=n=01nn!tαnQ1x,t+tαQ2x,t,Q1x,t=1πtα+βH1,22,0x24tα+β1+αnα+β2,α+β1,1,12+n,1,Q2x,t=1πtα+βH1,22,0x24tα+β1+αn+αβ2,α+β1,1,12+n,1(4.43)

for α + β ≤ 1, see Figures 6a and 6b for the graphical representations of the closed-form (4.43). The subordination of the process governed by TFCE to the Gaussian process is verified whenever the distribution-like Q(x, t) is non-negative. In the green curve of Figure 6a the non-negativity condition has been clearly breached due to the choice 1 < α + β < 2 violating the condition α + β < 1 established above. It comes to our attention that the general form of distribution-like (4.41) can be derived by using the Laplace transform of equations (4.23), (4.34) and (4.37) and utilising the useful relations (2.20) and (2.21) in . As applied in previous sections, the Mellin transform can be employed here to check the normalisation of (4.43) and obtain the δ th-order moment of QTFCE(x, t), one can obtain the most interesting result, Figure 6

(a) Spatial evolution of the distribution-like Q(x, t) of the time-fractional Cattaneo-type equation with two fractional parameters at dimensionless time t = 0.5. (b) Spatial distribution of Q(x, t) of TFCE subordinating to the Gaussian process (non-negative) at different instants of time and fractional parameters α = 0.5, β = 0.4.

xQδt=δ+1xϱδt,(4.44)

which suggests an alternative definition for the distribution of the diffusing substance with alternative governing equations (4.29) and (4.30) for time-fractional finite-speed diffusion, (4.18) and (4.19) for time-fractional diffusion, and (4.7) for the Gaussian diffusion.

### Remark 4.3

The relation between the fluxes, (4.36) and (4.39) and the distribution-like (4.43) can be rewritten in the integral forms

JRLTFCEx,t=1x0tQTFCEx,τtτ2αβEα,α+β1tταdτ,(4.45)
JCTFCEx,t=1x0tQTFCEx,τtτ1αEα,αtταdτ,(4.46)

where Eα,β (z) is the Mittag-Leffler function defined in (A.22). In addition, the fluxes, (4.35) and (4.38) can be given in terms of the distribution-like (4.42) via

JRLTFDEx,t=1x0+RLDt1βQTFDEx,t,(4.47)
JCTFDEx,t=1xQTFDEx,t.(4.48)

## 5 Concluding remarks

In this work we addressed four Cattaneo-type models with time fractality based on Riemann-Liouville and Caputo single-order and distributed-order fractional derivatives. Many processes which transit from ballistic or superdiffusive behaviour in the short-time limit to subdiffusive behaviour, or even without transport 〈x2(t)〉 ∼ 1, in the long-time limit are found to be captured by these Cattaneo-like equations. We summarise the MSD of the proposed models in Table 1.

Table 1

Limiting behaviour of the MSD for the proposed Cattaneo-like models, along with the respective equation numbers.

Model [resp. equations]x2(t)〉x2(t)〉
t → 0+t → ∞
2-parameter TFCE/GCE-I [(2.3), (2.4), (4.20)]tα+βtβ
DTFCE-I [(2.8), (2.11), (2.12)]p1tα+1p2tβ
DTFCE-II [(2.9), (2.11), (2.12)]tβ+1/p2tα/p1
DTFCE-III [(2.10), (2.11), (2.12)]tβ+1/p2t

Conditional closed-form solutions for the proposed models were derived in terms of the Fox H-function. Alternative formulations for the Green’s function of normal, time-fractional, time-fractional finite-speed diffusion were provided. For the generalised Cattaneo equation of type I with two fractional parameters, this modified Green’s function can be viewed as a gauge for the subordination to the Gaussian process, where the non-negativity of the modified Green’s function is broken in the transition from Gaussian to non-Gaussian. It is worth mentioning that we have disregarded the space-fractality which submits a generic crossover from superballistic/subdiffusion to subdiffusion/superballistic. Such space-fractional kinetic models can be found in the double-order space-fractional diffusion equation in the natural/modified forms [82, 15], the focus of a separate future work.

## Appendix A. Mathematical preliminaries

In this appendix we give a short brief for the mathematical preliminaries needed throughout the paper.

### A.1 Fractional derivatives and transform techniques

We start with the left-sided Caputo (C) fractional derivative of order α ∈ (0, 1], defined for any well-behaved function f(t) by [66, 35]

a+CDtαft=a+RLIt1αtft,0<α<1,tft,α=1,(A.1)

a ∈ ℝ, where a+RLItα is the left-sided Riemann-Liouville fractional integral given as

a+RLItαft=1Γαatfτtτ1αdτ,α>0,ft,α=0.(A.2)

The left-sided Riemann-Liouville (RL) fractional derivative of order α ∈ (0, 1] is defined as

a+RLDtαft=ta+RLIt1αft,0<α<1,tft,α=1.(A.3)

The regularised Prabhakar fractional derivative is defined by [67, 22]

a+PDα,λ,tδ,βft=attτn1βEα,nβδλtτατnfτdτ,(A.4)

where n − 1 < β < n, and Eα,βδ(z) is the Prabhakar generalisation of Mittag-Leffler function, see equation (A.16).

The Laplace transform L{f(t);s}=0estf(t)dt=f~(s) any generic function f(t) is given for the RL and C fractional derivatives of f(t) by

L0+RLDtαft;s=sαf~s,L0+CDtαft;s=sαf~ssα1f0+,(A.5)

where α ∈ (0, 1]. The Fourier transform of any generic function g(x) is defined by

Fgx;k=eıkxgxdx=g^k,F1g^k;x=12πeıkxg^kdx=gx.(A.6)

### A.2 Fox H-function

The Fox H-function is defined in terms of the Mellin-Barnes integral 

Hp,qm,nxa1,A1,,(ap,Ap)b1,B1,,(bq,Bq)=12πıΩΘsxsds,(A.7)

where m, n, p, and q are integers satisfying 0 ≤ np, 1 ≤ mq, ai, bj ∈ ℂ, Ai, Bj ∈ ℝ+, i = 1, ⋯, p, j = 1, ⋯, q, and the function Θ(s) is given by

Θs=j=1mΓ(bjBjs)j=1nΓ(1aj+Ajs)j=m+1qΓ(1bj+Bjs)j=n+1pΓ(ajAjs),(A.8)

where Γ(⋅) is the Gamma function. The contour Ω in the right side of equation (1) separates the poles of Γ(bj + Bjs), j = 1, ⋯, m from the poles of Γ(1 – aiAis), i = 1, ⋯, n.

The Mellin transform, defined for any generic function as 𝓜{f(x)}(z) = 0xz1f(x)dx, of the H -function is given by

0xz1Hp,qm,naxap,Apbq,Bqdx=azΘz,(A.9)

where Θ(z) is given by (A.8).

The inverse Fourier transform of the H-function H1,21,1akδ is given by

F1kλH1,21,1akδn,10,1;(β,γ)x=14πaλ+1δH2,32,1xδ2δa1λ+1δ,1;1βλ+1δγ,γ0,δ2,1+nλ+1δ,1;12,δ2,(A.10)

where k is the Fourier variable, a, γ, δ ∈ ℝ+, β, n ∈ ℂ, and λ ∈ ℝ+ ∪ {0}.

### Remark A.1

The proof of relation (A.10) has been omitted for brevity. For more details about the method, the reader may consult [20, 56, 36].

If the poles of j=1mΓ(bjBjs) are simple, the following series expansion holds true

Hp,qm,nxap,Apbq,Bq=h=1mν=01νxbh+νBhν!Bh×j=1,jhmΓbjBjbh+νBhj=1nΓ1aj+Ajbh+νBhj=m+1qΓ1bj+Bjbh+νBhj=n+1pΓajAjbh+νBh.(A.11)

The following properties are used throughout the paper,

Hp,qm,nxδap,Apbq,Bq=1δHp,qm,nxap,Apδbq,Bqδ,δ>0,(A.12)
drdxrxλHp,qm,naxδap,Apbq,Bq=xλrHp+1,q+1m,n+1axδλ,δap,Apbq,Bq,rλ,δ,(A.13)

where a, δ ∈ ℝ+, and r ∈ ℕ ∪ {0},

Hp,qm,nxa1,A1,,(ap,Ap)b1,B1,,bq1,Bq1,a1,A1=Hp1,q1m,n1xa2,A2,,(ap,Ap)b1,B1,,bq1,Bq1,n1,q>m,(A.14)

and

H0,11,0xb,B=1BxbBexpx1B.(A.15)

### A.3 Generalised Mittag-Leffler functions

The Prabhakar generalisation of Mittag-Leffler function (PML) is defined in the series form 

Eα,βγz=n=0γnΓ(αn+β)znn!,α,β,γ,zC,Reα>0,(A.16)

where (γ)n is the ascending Pochhammer symbol defined by (γ)0 = 1, (γ)n = γ(γ + 1) ⋯ (γ + n – 1) = Γ(γ+n)Γ(γ). The PML function, Eα,βγ(–λtα) is a completely monotone function for t ≥ 0, λ is positive constant, 0 < α, β ≤ 1, and 0 < γβ / α, has the asymptotic representation

Eα,βγ(λtα)1ΓβλγtαΓα+β,t0+,λtαγΓβαγ,t,(A.17)

where the short-time behaviour is deduced from the series representation (A.16) and the long-time behaviour can be obtained from the series

Eα,βγz=zγΓ(γ)n=0Γ(n+γ)Γβαn+γznn!,z>1.(A.18)

The most important properties of PML are its Laplace transform

Ltβ1Eα,βγλtα;s=sαγβsα+λγ,(A.19)

where s is the Laplace variable, and its relation to the Fox H-function,

Eα,βγz=1ΓγH1,21,1z(1γ,1)0,1,(1β,α).(A.20)

The following recurrence relation can be validated by definition (A.16)

Eα,βγ+1tα+tαEα,β+αγ+1tα=Eα,βγtα.(A.21)

When γ = 1, the PML reduces to the generalised Mittag-Leffler function with two parameters,

Eα,β1z=Eα,βz=n=0znΓ(αn+β),α,β,zC,Reα>0,(A.22)

and when γ = β = 1, the classical Mittag-Leffler function is recovered,

Eα,1z=Eαz=n=0znΓ(αn+1),α,zC,Reα>0.(A.23)

### A.4 Unit doublet function

The unit doublet function (differentiator)  is defined as the first derivative of the Dirac delta function,

u1x=xδx,(A.24)

with the operational properties

Fu1x;k=ık,fxξu1ξdξ=xfx.(A.25)

It can be generalised to the nth differentiator function

unx=xnδx,(A.26)

where n ∈ ℕ ∪ {0}. When n is negative integers, equation (A.26) defines then the nth integrator function, whereas n = 0 recovers the conventional Dirac delta function (unit impulse function).

Let us introduce the “modified unit doublet” function as a product of the spatial variable x and the unit doublet function,

U1x=xu1x,(A.27)

whose Fourier transform is given by

FU1x;k=1=Fδx;k.(A.28)

Similar to the classical definition of the unit impulse function we consider the following forms for the unit doublet and its modification,

δx=lima0+1aπexpx2a2,u1x=lima0+2xa2aπexpx2a2,U1x=lima0+2x2a2aπexpx2a2.(A.29)

### Definition A.1

The function f : (0, ∞) → [0, ∞) is a completely monotone function if fC and f satisfies the Bernstein–Hausdorff–Widder condition; (–1)nf(n)(t) ≥ 0 for all n ∈ ℕ. The family of all completely monotone functions is denoted by 𝓒𝓜𝓕. The functions λα and exp(–aλ), a, λ > 0, α < 0, are typical examples of 𝓒𝓜𝓕. The product of two complete monotone functions and the linear combination of complete monotone functions are also a complete monotone function. The function φ ∈ 𝓒𝓜𝓕 if, and only if, there exists a function f ≥ 0 such that φ(λ) = 0 exp(–λt) f(t) dt, λ > 0.

### Definition A.2

The function f : (0, ∞) → [0, ∞) is a Bernstein function if fC and f satisfies the condition: (–1)n–1f(n)(t) ≥ 0 for all n ∈ ℕ. The set of all Bernstein functions is denoted by 𝓑𝓕.

### Definition A.3

The function φ : (0, ∞) → [0, ∞) is said to be a Stieltjes function if there exists f ∈ 𝓒𝓜𝓕 such that φ(λ) = 0 exp(–λt) f(t)dt, for λ > 0. Denoting by 𝓢𝓕 to the family of all Stieltjes functions, it is obvious that 𝓢𝓕 ⊂ 𝓒𝓜𝓕 since f ≥ 0. The functions λα, λα–1 and 1u+λ, where α ∈ [0, 1], λ > 0 and u > 0 are typical examples of Stieltjes function. The linear combination of Stieltjes functions are also a Stieltjes function.

### Definition A.4

The Bernstein function φ(λ) is said to be a complete Bernstein function if, and only if, φ(λ)/λ ∈ 𝓢𝓕. Denoting by 𝓒𝓑𝓕 to the set of all complete Bernstein functions, thus 𝓒𝓑𝓕 ⊂ 𝓑𝓕. The functions λα and λ1–α, where α ∈ [0, 1] and λ > 0, are typical examples of complete Bernstein functions. The linear combination of complete Bernstein functions are also a complete Bernstein function. The set of complete Bernstein functions is not in general closed under multiplication, however, if φ, ψ ∈ 𝓒𝓑𝓕, then [φ(λ)]α [ψ(λ)]β ∈ 𝓒𝓑𝓕 provided that α, β ∈ (0, 1) and α + β ≤ 1. If φ(λ) ∈ 𝓒𝓑𝓕, then exp(–(λ)) ∈ 𝓒𝓜𝓕 for a, λ > 0.

## Appendix B. Derivation of the solutions

Here, we consider four different problems of inverting the Laplace-Fourier transform, which are applied in the derivation of our solutions.

### B.1 Problem B.1

Let us consider the inversion of the Laplace-Fourier transform of the following equation

ϱ~^k,s=sγ1+χ0,1sβ1sγ+χ0,1sβ+k2,

where χ0,1 ∈ {0, 1} and γβ such that γ = β if χ0,1 = 0 and γ > β if χ0,1 = 1. Rearranging the above equation, we have

ϱ~^k,s=sγ1+χ0,1sβ1sγ+k21+f0k,s1,

where f0k,s=χ0,1sβsγ+k2. When χ0,1 = 1, we note that f0(k, s) < 1 for γ = α + β, 0 < α, β ≤ 1, Re{s} > 0, and |k| > 0.7. These conditions provide a strong solution for Problem B.1. The long-space behaviours (|k| ≤ 0.7) are not of interest in our investigation. Using the condition that f0(k, s) < 1, then we have

ϱ~^k,s=n=0χ0,1nsβn+γ1+χ0,1sβn+11sγ+k2n+1,

which can be converted to the physical domain by means of (A.10), (A.19) and (A.20), and we obtain

ϱx,t=n=0χ0,1nn!tγβnϱ1x,t+χ0,1tγβϱ2x,t,ϱ1x,t=14πtγH1,22,0x24tγ1+γβnγ2,γ0,1,12+n,1,ϱ2x,t=14πtγH1,22,0x24tγ1+γβn+1γ2,γ0,1,12+n,1.(B.1)

### Remark B.1

When χ0,1 = 0 in (B.1), namely the only non-vanishing term is the first one n=0,n=0χ0,1n=1, and γ = β then we have the classical solution of the time-fractional diffusion equation (TFDE) ,

ϱTFDEx,t=14πtβH1,22,0x24tβ1β2,β0,1,12,1.

### B.2 Problem B.2

If the well-behaved function ϱ(x, t) has the Laplace-Fourier transformed form

ϱ~^k,s=sγ1+χ0,1sα1sγ+χ0,1sα+k2p1+p2sαβ,

where χ0,1 ∈ {0, 1}, γα such that γ = α if χ0,1 = 0 and γ > α if χ0,1 = 1, 0 < α < β ≤ 1 and p1 + p2 = 1, then

ϱx,t=n=0tγαnn!m=0nnmχ0,1nmp2p1tβγm×ϱ1x,t+χ0,1tγαϱ2x,t,ϱ1x,t=14πp1tγ×H2,32,1x24p1tγ12m,1;1+γαn+βγmγ2,γ0,1,12+nm,1;12,1,ϱ2x,t=14πp1tγ×H2,32,1x24p1tγ12m,1;1+γαn+1+βγmγ2,γ0,1,12+nm,1;12,1,(B.2)

provided that p1p2. The proof of the above claim can be accomplished by firstly rewriting ϱ̂͠(k, s) in the following form

ϱ~^k,s=sγ1+χ0,1sα1sγ+p1k21+f1k,s;p1,p21,

where f1k,s;p1,p2=χ0,1sα+p2sαβk2sγ+p1k2. By examining the properties of f1(k, λ; p1, p2) for λ = Re(s) ∈ ℝ+, we have the following cases: If χ0,1 = 0 and γ = α, we have f1k,λ;p1,p2=p2λαβk2λα+p1k2 which satisfies that f1(k, λ; p1, p2) < 1 for λ > 1, |k| ∈ ℝ, provided that p1p2. From these conditions, one can anticipate that the derived solution works well for the small and intermediate values of time, but it diverges for large values of time (for example t = 10). Secondly, if χ0,1 = 1, and γ = α + 1, we have then f1k,λ;p1,p2=λα+p2λαβk2λα+1+p1k2, which also satisfies that f1(k, λ; p1, p2) < 1 for λ > 1, |k| ∈ ℝ, provided that p1p2. Following arguments not dissimilar to those in  and the solution of Problem B.1 lead us to (B.2).

### Remark B.2

If χ0,1 = 0 and γ = α, then χ0,1nm indicates that the only nonvanishing term of the second series of (B.2) is the term m = n, namely,

ϱx,t=14πp1tαn=01nn!p2p1tβαn×H2,32,1x24p1tα12n,1;1+βαnα2,α0,1,12,1;12,1,

where p1p2, which coincides the solution of the double-order time-fractional diffusion equation of the modified form derived in .

### B.3 Problem B.3

If the well-behaved function ϱ(x, t) has the form

ϱ~^k,s=p1sγ11+p2sγ21+χ0,1p1sα1+p2sβ1p1sγ1+p2sγ2+χ0,1p1sα+p2sβ+k2,

in Laplace-Fourier space, where 0 < α < β ≤ 1, p1 + p2 = 1, γ1α and γ2β are given by

γ1=α,if χ0,1=0,α+1,if χ0,1=1,γ2=β,if χ0,1=0,β+1,if χ0,1=1,

then

ϱx,t=n=0tγ2γ1nn!m=0nnmχ0,1mp1p2nm×tγ1βm=0mmp1p2tβα×p1p2tγ2γ1ϱ1x,t+ϱ2x,t+χ0,1tγ2βp1p2tβαϱ3x,t+ϱ4x,t,ϱ1x,t=14πtγ2p2×H1,22,0p2x24tγ21+γ21(n+1)+γ1βm+βαγ22,γ20,1,12+n,1,ϱ2x,t=14πtγ2p2×H1,22,0p2x24tγ21+γ21n+γ1βm+βαγ22,γ20,1,12+n,1,ϱ3x,t=14πtγ2p2×H1,22,0p2x24tγ21+γ21(n+1)+γ1βm+βα+γαγ22,γ20,1,12+n,1,ϱ4x,t=14πtγ2p2×H1,22,0p2x24tγ21+γ21n+γ1βm+βα+γβγ22,γ20,1,12+n,1,(B.3)

provided that p1p2 for χ0,1 = 0, and p1 < p2 for χ0,1 = 1. Here, we used γ21 = γ2γ1, γα = γ1α, and γβ = γ2β.

In order to validate equation (B.3) we first rehash ϱ̂͠(k, s) to take the form

ϱ~^k,s=p1sγ11+p2sγ21+χ0,1p1sα1+p2sβ1p2sγ2+k21+f2k,s;p1,p21,

where

f2k,s;p1,p2=p1sγ1+χ0,1p1sα+p2sβp2sγ2+k2.

Upon investigating the properties of f2(k, λ; p1, p2), λ = Re(s) > 0, we have two special cases of interest: When χ0,1 = 0, we get the characteristic function for DTFDE-II, given by f2k,λ;p1,p2=p1λαp2λβ+k2, with a graphical representation reports that for p1p2, λ > 0, |k| ≥ 0.6 and 0 < α < β ≤ 1, then f2(k, λ; p1, p2) < 1, otherwise, the resulting solution diverges. While, when χ0,1 = 1, we have f2k,λ;p1,p2=p1λα+1+p1λα+p2λβp2λβ+1+k2, which satisfies f2(k, λ; p1, p2) < 1 provided that p1 < p2, λ > 0, |k| ≥ 1 and 0 < α < β ≤ 1. We second apply the same method of proving Problem B.2, we arrive at (B.3).

### Remark B.3

The solution of the double-order time-fractional diffusion equation of the natural form obtained by  can be derived from (B.3) by setting χ0,1 = 0, and thus the only non-vanishing term of the second series is the first term, m = 0, namely,

ϱx,t=n=01nn!p1p2tβαnp1p2tβαϱ1x,t+ϱ2x,t,
ϱ1x,t=14πtβ/p2H1,22,0x24tβ/p21+βαn+1β2,β0,1,12+n,1,
ϱ1x,t=14πtβ/p2H1,22,0x24tβ/p21+βαnβ2,β0,1,12+n,1,

where p1p2.

### B.4 Problem B.4

If the well-behaved function ϱ(x, t) is given by

ϱ~^k,s=χ0,1+p1sγ11+p2sγ21p1sγ1+p2sγ2+χ0,1s+k2,

in Laplace-Fourier space, where p1 + p2 = 1, γ1 and γ2 are given through

γ1=α,if χ0,1=0,α+1,if χ0,1=1,γ2=β,if χ0,1=0,β+1,if χ0,1=1,

and 0 < α < β ≤ 1, then

ϱx,t=n=01nn!tγ21p2nm=0nnmχ0,1nmp1tγ11m×χ0,1tγ21p2ϱ1x,t+p1p2tγ2γ1ϱ2x,t+ϱ3x,t,ϱ1x,t=14πtγ2/p2×H1,22,0x24tγ2/p21+γ21n+1γ11mγ22,γ20,1,12+n,1,ϱ2x,t=14πtγ2/p2×H1,22,0x24tγ2/p21+γ21n+1γ11m+1γ22,γ20,1,12+n,1,ϱ3x,t=14πtγ2/p2×H1,22,0x24tγ2/p21+γ21nγ11mγ22,γ20,1,12+n,1,(B.4)

provided that p1p2 for χ0,1 = 0, and p1 ∈ [0, 0.3], p2 ∈ [0.7, 1], p1 + p2 = 1, and β ∈ [0.4, 1] for χ0,1 = 1.

### Remark B.4

We have written off the solution of Problem B.4. Setting χ0,1 = 0 in (B.4) results in again the solution of double-order time-fractional diffusion equation of the natural form .

## Acknowledgements

RM acknowledges an Alexander von Humboldt Polish Honorary Research Scholarship from the Foundation for Polish Science (Fundacja na rzecz Nauki Polskiej).

## References

 T.M. Atanackovic, S. Pilipovic, On a constitutive equation of heat conduction with fractional derivatives of complex order. Acta Mechanica229, No 3 (2018), 1111–1121.10.1007/s00707-017-1959-4Search in Google Scholar

 E. Awad, On the time-fractional Cattaneo equation of distributed order. Physica A518 (2019), 210–233.10.1016/j.physa.2018.12.005Search in Google Scholar

 E. Awad, On the generalized thermal lagging behavior: Refined aspects. J. Thermal Stresses35, No 4 (2012), 293–325.10.1080/01495739.2012.663682Search in Google Scholar

 E. Barkai, Y. Garini, R. Metzler, Strange kinetics of single molecules in living cells. Phys. Today65, No 8 (2012), 29–35.10.1063/PT.3.1677Search in Google Scholar

 E. Bazhlekova, Subordination in a class of generalized time-fractional diffusion-wave equations. Fract. Calc. Appl. Anal. 21, No 4 (2018), 869–900; 10.1515/fca-2018-0048; https://www.degruyter.com/view/j/fca.2018.21.issue-4/issue-files/fca.2018.21.issue-4.xml.Search in Google Scholar

 J.-P. Bouchaud, A. Georges, Anomalous diffusion in disordered media: statistical mechanisms, models and physical applications. Phys. Rep. 195, No 4-5 (1990), 127–293.10.1016/0370-1573(90)90099-NSearch in Google Scholar

 J.-P. Bouchaud, A. Comtet, A. Georges, and P. Le Doussal, Classical diffusion of a particle in a one-dimensional random force field. Ann. Phys. (N.Y.)201, No 2 (1990), 285–341.10.1016/0003-4916(90)90043-NSearch in Google Scholar

 L. Boyadjiev, Y. Luchko, The neutral-fractional telegraph equation. Math. Modelling Nat. Phenomena12, No 6 (2017), 51–67.10.1051/mmnp/2017064Search in Google Scholar

 M. Caputo, Distributed order differential equations modelling dielectric induction and diffusion. Fract. Calc. Appl. Anal. 4, No 4 (2001), 421–442.Search in Google Scholar

 C. Cattaneo, Sulla conduzione del calore. Atti Sem. Mat. Fis. Univ. Modena3, No 3 (1948), 83–101.10.1007/978-3-642-11051-1_5Search in Google Scholar

 A. Chechkin, R. Gorenflo, I. Sokolov, Retarding subdiffusion and accelerating superdiffusion governed by distributed-order fractional diffusion equations. Phys. Rev. E66, No 4 (2002), Art. 046129.10.1103/PhysRevE.66.046129Search in Google Scholar PubMed

 A.V. Chechkin, J. Klafter, I.M. Sokolov, Fractional Fokker-Planck equation for ultraslow kinetics. Europhys. Lett. 63, No 3 (2003), 326–332.10.1209/epl/i2003-00539-0Search in Google Scholar

 A. Chechkin, V.Y. Gonchar, R. Gorenflo, N. Korabel, I. Sokolov, Generalized fractional diffusion equations for accelerating subdiffusion and truncated Lévy flights. Phys. Rev. E78, No 2 (2008), Art. 021111.10.1103/PhysRevE.78.021111Search in Google Scholar PubMed

 A.V. Chechkin, R. Gorenflo, I.M. Sokolov, V.Y. Gonchar, Distributed order time fractional diffusion equation. Fract. Calc. Appl. Anal. 6, No 3 (2003), 259–280.Search in Google Scholar

 A. Chechkin, I.M. Sokolov, J. Klafter, Natural and modified forms of distributed-order fractional diffusion equations. In: Fractional Dynamics: Recent Advances, Ed. by J. Klafter, S.C. Lim, R. Metzler, World Scientific, Singapore (2011), 107–127.10.1142/9789814340595_0005Search in Google Scholar

 A. Compte, R. Metzler, The generalized Cattaneo equation for the description of anomalous transport processes. J. Phys. A30, No 22 (1997), 7277–7289.10.1088/0305-4470/30/21/006Search in Google Scholar

 A. Dhar, Fractional equation description of an open anomalous heat conduction set-up. J. Stat. Mech. 2019, No 1 (2019), Art. 013205.10.1088/1742-5468/aaf630Search in Google Scholar

 A. Einstein, Über die von der molekularkinetischen Theorie der Wärme geforderte Bewegung von in ruhenden Flüssigkeiten suspendierten Teilchen. Ann. Phys. (Leipzig)322, No 8 (1905), 549–560.10.1002/andp.19053220806Search in Google Scholar

 A. Erdélyi, W. Magnus, F. Oberhettinger, F.G. Tricomi, Tables of Integral Transforms, Vol. 1. McGraw-Hill, New York (1954).Search in Google Scholar

 W.G. Glöckle, T.F. Nonnenmacher, Fox function representation of non-Debye relaxation processes. J. Stat. Phys. 71, No 3-4 (1993), 741–757.10.1007/BF01058445Search in Google Scholar

 A. Godec, A. V. Chechkin, E. Barkai, H. Kantz, and R. Metzler, Localization and universal fluctuations in ultraslow diffusion processes. J. Phys. A47, No 49 (2014), Art. 492002.10.1088/1751-8113/47/49/492002Search in Google Scholar

 R. Garra, R. Gorenflo, F. Polito, Ž. Tomovski, Hilfer–Prabhakar derivatives and some applications. Appl. Math. Comput. 242 (2014), 576–589.10.1016/j.amc.2014.05.129Search in Google Scholar

 R. Gorenflo, Y. Luchko, M. Stojanović, Fundamental solution of a distributed order time-fractional diffusion-wave equation as probability density. Fract. Calc. Appl. Anal. 16, No 2 (2013), 297-316; 10.2478/s13540-013-0019-6; https://www.degruyter.com/view/j/fca.2013.16.issue-2/issue-files/fca.2013.16.issue-2.xml.Search in Google Scholar

 R. Gorenflo, A.A. Kilbas, F. Mainardi, S.V. Rogosin, Mittag-Leffler Functions, Related Topics and Applications. Springer, Berlin (2014).10.1007/978-3-662-43930-2Search in Google Scholar

 R.A. Guyer, J.A. Krumhansl, Solution of the linearized phonon Boltzmann equation. Phys. Rev. 148, No 2 (1966), 766–778.10.1103/PhysRev.148.766Search in Google Scholar

 S. Havlin and G. H. Weiss, A new class of long-tailed pausing time densities. J. Stat. Phys. 58, No 5-6 (1990), 1267–1273.10.1007/BF01026577Search in Google Scholar

 R. Hilfer, Mathematical and physical interpretations of fractional derivatives and integrals. In: Handbook of Fractional Calculus with Applications, Vol. 1, Basic Theory, Ed. by A. Kochubei, Y. Luchko, Berlin, De Gruyter (2019).Search in Google Scholar

 R. Hilfer, Fractional diffusion based on Riemann-Liouville fractional derivatives. J. Phys. Chem. B104, No 16 (2000), 3914–3917.10.1021/jp9936289Search in Google Scholar

 F. Höfling, T. Franosch, Anomalous transport in the crowded world of biological cells. Rep. Prog. Phys. 76, No 4 (2013), Art. 046602.10.1088/0034-4885/76/4/046602Search in Google Scholar

 P. De Jagher, A hyperbolic “diffusion equation” taking a finite collision frequency into account. Physica A101, No 2-3 (1980), 629–633.10.1016/0378-4371(80)90200-9Search in Google Scholar

 J.-H. Jeon, H.M.-S. Monne, M. Javanainen, R. Metzler, Anomalous diffusion of phospholipids and cholesterols in a lipid bilayer and its origins. Phys. Rev. Lett. 109, No 18 (2012), Art. 188103.10.1103/PhysRevLett.109.188103Search in Google Scholar PubMed

 J.-H. Jeon, M. Javanainen, H. Martinez-Seara, R. Metzler, I. Vattulainen, Protein crowding in lipid bilayers gives rise to non-Gaussian anomalous lateral diffusion of phospholipids and proteins. Phys. Rev. X6, No 2 (2016), Art. 021006.10.1103/PhysRevX.6.021006Search in Google Scholar

 D. Jou, G. Lebon, J. Casas-Vázquez, Extended Irreversible Thermodynamics. Springer, Berlin (2014).Search in Google Scholar

 N.G. van Kampen, Stochastic Processes in Physics and Chemistry. North Holland, Amsterdam (1981).Search in Google Scholar

 A.A. Kilbas, H.M. Srivastava, J.J. Trujillo, Theory and Applications of Fractional Differential Equations. Elsevier, Amsterdam (2006).Search in Google Scholar

 T. Langlands, Solution of a modified fractional diffusion equation. Physica A367 (2006), 136–144.10.1016/j.physa.2005.12.012Search in Google Scholar

 P. Langevin, Sur la théorie du mouvement brownien. C. R. Acad. Sci. Paris146, No 10 (1908), 530–533.Search in Google Scholar

 Y.B. Lev, D.M. Kennes, C. Klöckner, D.R. Reichman, C. Karrasch, Transport in quasiperiodic interacting systems: From superdiffusion to subdiffusion. Europhys. Lett. 119, No 3 (2017), Art. 37003.10.1209/0295-5075/119/37003Search in Google Scholar

 Y. Luchko, Initial-boundary-value problems for the generalized multi-term time-fractional diffusion equation. J. Math. Anal. Appl. 374, No 2 (2011), 538–548.10.1016/j.jmaa.2010.08.048Search in Google Scholar

 F. Mainardi, Fractional relaxation-oscillation and fractional diffusion-wave phenomena. Chaos, Solitons & Fractals7, No 9 (1996), 1461–1477.10.1016/0960-0779(95)00125-5Search in Google Scholar

 F. Mainardi, The fundamental solutions for the fractional diffusion-wave equation. Appl. Math. Lett. 9, No 6 (1996), 23–28.10.1016/0893-9659(96)00089-4Search in Google Scholar

 R. Gorenflo, F. Mainardi, Fractional calculus. In: Fractals and Fractional Calculus in Continuum Mechanics, Springer, Berlin (1997), 223-276.10.1007/978-3-7091-2664-6_5Search in Google Scholar

 F. Mainardi, G. Pagnini, R. Gorenflo, Some aspects of fractional diffusion equations of single and distributed order. Applied Mathematics and Computation187, No 1 (2007), 295–305.10.1016/j.amc.2006.08.126Search in Google Scholar

 R.N. Mantegna, H.E. Stanley, Stochastic process with ultraslow convergence to a Gaussian: the truncated Lévy flight. Phys. Rev. Lett. 73, No 22 (1994), 2946–2949.10.1103/PhysRevLett.73.2946Search in Google Scholar

 J. Masoliver, Fractional telegrapher’s equation from fractional persistent random walks. Phys. Rev. E93, No 5 (2016), Art. 052107.10.1103/PhysRevE.93.052107Search in Google Scholar

 J. Masoliver, G.H. Weiss, Finite-velocity diffusion. Euro. J. Phys. 17, No 4 (1996), 190–196.10.1088/0143-0807/17/4/008Search in Google Scholar

 Y. Meroz, I.M. Sokolov, A toolbox for determining subdiffusive mechanisms. Phys. Rep. 573 (2015), 1–29.10.1016/j.physrep.2015.01.002Search in Google Scholar

 R. Metzler, J. Klafter, The random walk’s guide to anomalous diffusion: a fractional dynamics approach. Phys. Rep. 339, No 1 (2000), 1–77.10.1016/S0370-1573(00)00070-3Search in Google Scholar

 R. Metzler, J. Klafter, The restaurant at the end of the random walk: recent developments in the description of anomalous transport by fractional dynamics. J. Phys. A37, No 31 (2004), R161-R208.10.1088/0305-4470/37/31/R01Search in Google Scholar

 R. Metzler, J. Klafter, I.M. Sokolov, Anomalous transport in external fields: Continuous time random walks and fractional diffusion equations extended. Phys. Rev. E58, No 2 (1998), 1621–1633.10.1103/PhysRevE.58.1621Search in Google Scholar

 R. Metzler, Generalized Chapman-Kolmogorov equation: A unifying approach to the description of anomalous transport in external fields. Phys. Rev. E62, No 5 (2000), 6233–6245.10.1103/PhysRevE.62.6233Search in Google Scholar

 R. Metzler, E. Barkai, and J. Klafter, Anomalous diffusion and relaxation close to thermal equilibrium: A fractional Fokker-Planck equation approach. Phys. Rev. Lett. 82, No 18 (1999), 3563–3567.10.1103/PhysRevLett.82.3563Search in Google Scholar

 R. Metzler, E. Barkai, and J. Klafter, Deriving fractional Fokker-Planck equations from a generalised master equation. Europhys. Lett. 46, No 4 (1999), 431–436.10.1209/epl/i1999-00279-7Search in Google Scholar

 R. Metzler, J.-H. Jeon, and A. G. Cherstvy, Non-Brownian diffusion in lipid membranes: experiments and simulations. Biochimica et Biophysica Acta (BBA) - Biomembranes1858, No 10 (2016), 2451–2467.10.1016/j.bbamem.2016.01.022Search in Google Scholar

 R. Metzler, J.-H. Jeon, A.G. Cherstvy, E. Barkai, Anomalous diffusion models and their properties: non-stationarity, non-ergodicity, and ageing at the centenary of single particle tracking. Phys. Chem. Chem. Phys. 16 (2014), 24128–24164.10.1039/C4CP03465ASearch in Google Scholar

 R. Metzler, T.F. Nonnenmacher, Space- and time-fractional diffusion and wave equations, fractional Fokker–Planck equations, and physical motivation. Chem. Phys. 284, No 1-2 (2002), 67-90.10.1016/S0301-0104(02)00537-2Search in Google Scholar

 D. Molina-Garcia, T. Sandev, H. Safdari, G. Pagnini, A. Chechkin, R. Metzler, Crossover from anomalous to normal diffusion: truncated power-law noise correlations and applications to dynamics in lipid bilayers. New J. Phys. 20 (2018), Art. 103027.10.1088/1367-2630/aae4b2Search in Google Scholar

 P.M. Morse, H. Feshbach, Methods of Theoretical Physics. McGraw-Hill, New York (1953).Search in Google Scholar

 F. Mainardi, A. Mura, G. Pagnini, R. Gorenflo, Time-fractional diffusion of distributed order. J. Vibration and Control14, No 9-10 (2008), 1267–1290.10.1177/1077546307087452Search in Google Scholar

 A.M. Mathai, R.K. Saxena, H.J. Haubold, The H-Function: Theory and Applications. Springer, Berlin (2010).10.1007/978-1-4419-0916-9Search in Google Scholar

 R.R. Nigmatullin, The realization of the generalized transfer equation in a medium with fractal geometry. Phys. Status Solidi133, No 1 (1986), 425–430.10.1515/9783112495483-049Search in Google Scholar

 T. Nonnenmacher, D. Nonnenmacher, Towards the formulation of a nonlinear fractional extended irreversible thermodynamics. Acta Phys. Hungarica66, No 1-4 (1989), 145–154.10.1007/BF03155787Search in Google Scholar

 K. Nørregaard, R. Metzler, C. Ritter, K. Berg-Sørensen, L.B. Oddershede, Manipulation and motion of organelles and single molecules in living cells. Chem. Rev. 117, No 5 (2017), 4342-4375.10.1021/acs.chemrev.6b00638Search in Google Scholar

 A.V. Oppenheim, Signals and Systems, MIT OpenCourseWare, Cambridge MA (2011); https://ocw.mit.edu/resources/res-6-007-signals-and-systems-spring-2011.Search in Google Scholar

 E. Orsingher, L. Beghin, Time-fractional telegraph equations and telegraph processes with Brownian time. Probability Theory and Related Fields128, No 1 (2004), 141–160.10.1007/s00440-003-0309-8Search in Google Scholar

 I. Podlubny, Fractional Differential Equations: An Introduction to Fractional Derivatives, Fractional Differential Equations, to Methods of Their Solution and Some of Their Applications. Elsevier, Amsterdam (1998).Search in Google Scholar

 T.R. Prabhakar, A singular integral equation with a generalized Mittag Leffler function in the kernel. Yokohama Math. J. 19, No 1 (1971), 7–15.Search in Google Scholar

 L.F. Richardson, Atmospheric diffusion shown on a distance-neighbour graph. Proc. Roy. Soc. A110, No 756 (1926), 709–737.Search in Google Scholar

 H. Qi, X. Jiang, Solutions of the space-time fractional Cattaneo diffusion equation. Physica A390, No 11 (2011), 1876–1883.10.1016/j.physa.2011.02.010Search in Google Scholar

 T.Q. Qiu, C.L. Tien, Short-pulse laser heating on metals. Intl. J. Heat and Mass Transfer35, No 3 (1992), 719–726.10.1016/0017-9310(92)90131-BSearch in Google Scholar

 T.Q. Qiu, C.L. Tien, Heat transfer mechanisms during short-pulse laser heating of metals. ASME J. Heat Transfer115, No 4 (1993), 835–841.10.1115/1.2911377Search in Google Scholar

 T. Sandev, A.V. Chechkin, N. Korabel, H. Kantz, I.M. Sokolov, R. Metzler, Distributed-order diffusion equations and multifractality: Models and solutions. Phys. Rev. E92, No 4 (2015), Art. 042117.10.1103/PhysRevE.92.042117Search in Google Scholar PubMed

 T. Sandev, A. Chechkin, H. Kantz, R. Metzler, Diffusion and Fokker-Planck-Smoluchowski equations with generalized memory kernel. Fract. Calc. Appl. Anal. 18, No 4 (2015), 1006–1038; 10.1515/fca-2015-0059; https://www.degruyter.com/view/j/fca.2015.18.issue-4/issue-files/fca.2015.18.issue-4.xml.Search in Google Scholar

 T. Sandev, I.M. Sokolov, R. Metzler, A. Chechkin, Beyond monofractional kinetics. Chaos, Solitons & Fractals102 (2017), 210–217.10.1016/j.chaos.2017.05.001Search in Google Scholar

 T. Sandev, Z. Tomovski, J.L. Dubbeldam, A. Chechkin, Generalized diffusion-wave equation with memory kernel. J. Phys. A52, No 1 (2019), Art. 015201.10.1088/1751-8121/aaefa3Search in Google Scholar

 C. Sándor, A. Libál, C. Reichhardt, C. Olson Reichhardt, Dewetting and spreading transitions for active matter on random pinning substrates. J. Chem. Phys. 146, No 20 (2017), Art. 204903.10.1063/1.4983344Search in Google Scholar PubMed PubMed Central

 R.K. Saxena, G. Pagnini, Exact solutions of triple-order time-fractional differential equations for anomalous relaxation and diffusion I: The accelerating case. Physica A390, No 4 (2011), 602–613.10.1016/j.physa.2010.10.012Search in Google Scholar

 R.L. Schilling, R. Song, Z. Vondracek, Bernstein Functions: Theory and Applications. De Gruyter, Berlin (2012).10.1515/9783110269338Search in Google Scholar

 W. Schneider, W. Wyss, Fractional diffusion and wave equations. J. Math. Phys. 30, No 1 (1989), 134–144.10.1063/1.528578Search in Google Scholar

 Y.G. Sinai, The limiting behavior of a one-dimensional random walk in a random medium. Theory of Probability & Applications27, No 2 (1982), 256-268.10.1007/978-1-4419-6205-8_2Search in Google Scholar

 M. von Smoluchowski, Zur kinetischen Theorie der Brownschen Molekularbewegung und der Suspensionen. Ann. Phys. (Leipzig)326, No 14 (1906), 756–780.10.1002/andp.19063261405Search in Google Scholar

 I. Sokolov, A. Chechkin, J. Klafter, Distributed-order fractional kinetics. Acta Phys. Polonica B35, No 47 (2004), 1323–1341.Search in Google Scholar

 I. Sokolov, Thermodynamics and fractional Fokker-Planck equations. Phys. Rev. E63, No 5 (2001), Art. 056111.10.1103/PhysRevE.63.056111Search in Google Scholar PubMed

 P. Tan, Y. Liang, Q. Xu, E. Mamontov, J. Li, X. Xing, L. Hong, Gradual crossover from subdiffusion to normal diffusion: a many-body effect in protein surface water. Phys. Rev. Lett. 120, No 24 (2018), Art. 248101.10.1103/PhysRevLett.120.248101Search in Google Scholar PubMed

 A.M. Tawfik, H. Fichtner, R. Schlickeiser, A. Elhanbaly, Analytical solutions of the space-time fractional telegraph and advection-diffusion equations. Physica A491 (2018), 810–819.10.1016/j.physa.2017.09.105Search in Google Scholar

 D.Y. Tzou, Macro-to Microscale Heat Transfer: The Lagging Behavior. John Wiley & Sons, New York (2014).10.1002/9781118818275Search in Google Scholar

 G.H. Weiss, Some applications of persistent random walks and the telegrapher’s equation. Physica A311, No 3-4 (2002), 381–410.10.1016/S0378-4371(02)00805-1Search in Google Scholar

 W. Wyss, The fractional diffusion equation. J. Math. Phys. 27, No 11 (1986), 2782–2785.10.1063/1.527251Search in Google Scholar

 M.A. Zaks, A. Nepomnyashchy, Subdiffusive and superdiffusive transport in plane steady viscous flows. Proc. Natl. Acad. Sci. USA116, No 37 (2018), 18245–18250.10.1073/pnas.1717225115Search in Google Scholar PubMed PubMed Central

 G.M. Zaslavsky, Hamiltonian Chaos and Fractional Dynamics. Oxford University Press, Oxford (2008).Search in Google Scholar

 V. Želi, D. Zorica, Analytical and numerical treatment of the heat conduction equation obtained via time-fractional distributed-order heat conduction law. Physica A492 (2018), 2316–2335.10.1016/j.physa.2017.11.150Search in Google Scholar 