Time-fractional telegraph equation of distributed order in higher dimensions with Hilfer fractional derivatives

1.   Introduction

Distributed-order fractional calculus (DOFC) is a branch of fractional calculus important for the modeling of complex systems. It generalizes the constant fractional operators by integrating the fractional kernel of these operators over an extended range of orders. The fractional differential operator of distributed order, for orders not great than 2, is given by

where $\frac{d^\alpha}{dt^\alpha}$ stands for a single-order fractional derivative and $b\left({\alpha} \right)$ is a non-negative weight function or generalized function. DOFC takes into account the superposition of orders and is a useful tool for modeling decelerating anomalous diffusion, ultraslow diffusive processes and strong anomaly (see e.g., ^[1,2]). DOFC models systems whose behavior stems from the complex interplay and superposition of nonlocal and memory effects occurring over a multitude of scales ^[3]. Over the last two decades, a significant number of papers appeared focusing on mathematical aspects and real-world applications of fractional partial differential equations with distributed order, see e.g., ^{[2,4,5,6,7,8,9,10,11,12,13,14,15]} for mathematical aspects, ^[16,17] for applications, and the review paper ^[18] about the mathematics of DOFC, including analytical, numerical methods, and the extensive overview of the recent applications of DOFC to fields like transport processes, and control theory. Moreover, DOFC was applied also in the study of composite materials ^[19,20] and viscoelastic materials having spatially varying properties ^[21].

The classical telegraph equation was first derived by Lord Kelvin in the 19th century ^[22]. It is a hyperbolic partial differential equation of the form

This equation was proposed by Cattaneo in 1958 (see ^[23]) to overcome the problem of infinite propagation velocity in heat transmission. Over the years, this equation and its time-fractional versions appeared in the study of several phenomena such as transmission lines for all frequencies ^[24], random walks ^[25], solar particle transport ^[26], oceanic diffusion ^[27], wave propagation ^[28], damped small vibrations, anomalous diffusion and wave-like processes ^{[29,30,31,32]}, scalar part of the Maxwell equations.

The TFTE with time-fractional derivatives of orders $\alpha_1 \in]0, 1]$ and $\alpha_2 \in]1, 2]$ was studied from the analytical, numerical, and probabilistic points of view by several authors. In ^[33], Cascaval et al. discussed the well-posedness of some initial-boundary value problems for the TFTE as well as the asymptotic behaviour of their solutions. In ^[31], the authors studied the neutral case of the TFTE and obtained an explicit Fourier representation of the fundamental solution (FS) and made a probabilistic interpretation of the FS in terms of stable probability density functions. Particular attention was given to the case $\alpha_1 = 1/2$ and $\alpha_2 = 1$ due to its connection of the telegraph process with Brownian motion. Some of these results were generalized by Camargo et al. in ^[34] for general $\alpha_1$ and $\alpha_2$ and studied later by Boyadjiev and Luchko in ^[29]. In ^[35], the authors considered a generalized telegraph equation with time-fractional derivatives in the Hilfer and Hadamard senses and space-fractional derivatives are in the sense of Riesz-Feller. Górska et al., (see ^[36]) considered various types of generalized telegraph equations and determine the conditions under which solutions can be recognized as probability density distributions.

The works ^{[32,37,38,39,40,41]} are examples of works devoted to the study of the TFTE in the multidimensional case with $n$ space variables, where in some cases the second derivative in space is replaced by the Euclidean Laplace operator. In ^[32] the authors solved the multi-dimensional TFTE with multi-term time-fractional derivatives and proved that its fundamental solution is the law of a stable isotropic multi-dimensional process time-changed. Ovidio et al., ^[41] constructed compositions of vector processes whose distribution is related to space-time fractional $n$-dimensional telegraph equations. We refer also the works of Masoliver and his co-workers about the TFTE and its connections with random walks (see ^{[42,43,44,45]}), and the recent survey paper ^[30] where is presented a very complete review of the fractional telegraph process. In ^[38,40] were employed Fourier, Laplace and Mellin transform techniques to obtain the first and second FS. Moreover, the application of the Residue Theorem allows obtaining double series representation for the FS of the TFTE in higher dimensions. Connections of the TFTE with fractional Clifford analysis and Sturm-Liouville theory were presented in ^[39] and ^[37].

In our recent paper ^[46] we studied the time-fractional telegraph equation with generalized distributed order in $\mathbb{R}^n \times \mathbb{R}^+$ finding a representation of the fundamental solution in terms of convolutions involving Fox H-functions. In this work, we extend our analysis to time-fractional telegraph equations of distributed order with Hilfer (or composite) time-fractional derivatives. Hilfer's derivative $_{t}D_{0^+}^{\gamma, \nu}$ was defined by Hilfer as a two-parameter family of fractional derivatives of order $\gamma > 0$ and type $\nu \in [0, 1]$ given by

where $I_{0^+}^\gamma$ denotes the left Riemann-Liouville fractional integral of order $\gamma > 0$ (see Eq (2.1) in Section 2). The Hilfer fractional derivative allows to interpolate smoothly between the Riemann-Liouville and the Caputo fractional derivatives (see ^[16,47,48]). These special cases are obtained when $\nu = 0$ and $\nu = 1,$ respectively. The type-parameter produces more stationary states, provides an extra degree of freedom on the initial condition, and increases the flexibility for the description of complex data. It was first used by Hilfer to describe the dynamics in glass formers over an extremely large-frequency window ^[49]. During the last years fractional differential equations with composite derivatives were studied by several authors, see e.g., ^{[35,50,51,52,53,54,55,56]}. In these works the equations are considered in $\mathbb{R} \times \mathbb{R}^+$, i.e., one single space variable and one time variable. Here, we consider the telegraph equation of distributed order with Hilfer time-fractional derivatives in the higher dimensional case, i.e., $\mathbb{R}^n \times \mathbb{R}^+.$

The paper is organized as follows. In Section 2 we recall some basic facts about fractional derivatives, integral transforms, and special functions, which are necessary for the development of this work. In Section 3 we formulate the problem of generalized distributed order telegraph equation for general density functions. Following the ideas presented in ^[46], we use a combination of Laplace, Fourier and Mellin transforms to obtain a representation of the solution of our equation via convolution integrals involving Fox H-functions. The key points to obtain our main result are the use of the classical Titchmarsh's Theorem to invert the Laplace transform, and the use of the Mellin transform to invert the Fourier transform. Some particular cases of our equation are analyzed by considering specific choices of the parameters of the equation. In Section 4 we compute the expression of the fractional moments of arbitrary order of the first fundamental solution in the Laplace domain. For the particular case of single-order derivatives, we invert the Laplace transform of the second-order moment obtaining an expression in terms of the three-parameter Mittag-Leffler function. For numerical purposes, we study the corresponding asymptotic behavior of the second-order moment in the time domain for the long and short time limit, by using the Tauberian Theorem. In the final part of the paper, we present and analyze some plots of the second-order moment for this particular case. As it will be shown the graphical representations support the analytical conclusions obtained via the Tauberian Theorem.

2.   Preliminaries

Let $a, b \in \mathbb{R}$ with $a < b$ and $\alpha > 0$. The left Riemann-Liouville fractional integral $I_{a^+}^\gamma$ of order $\gamma > 0$ is given by (see ^[57])

The left Hilfer (or composite) fractional derivative $_{t}D_{0^+}^{\gamma, \nu}$ of order $\gamma > 0$ and type $0 \leq \nu \leq 1$ is given by (see ^[16,47,48])

where $m = [\gamma]+1$ and $[\gamma]$ means the integer part of $\gamma$. We observe that in the case when $\nu = 0$ we recover the left Riemann-Liouville fractional derivative and in the case when $\nu = 1$ we have the left Caputo fractional derivative. The previous definitions of fractional integrals and derivatives can be naturally extended to $\mathbb{R}^n$ considering partial fractional integrals and derivatives (see Chapter 5 in ^[58]).

In this work, some integral transforms are used, namely, the Laplace, the Fourier and the Mellin transforms. The Laplace transform of a real-valued function $f\left({t} \right)$ is defined by (see ^[57])

and when it is applied to Eq (2.1) leads to (see ^[50])

where the initial-value terms $\left[{\frac{d^k}{dt^k} \left({{}_{t}I_{0^+}^{\left({1-\nu} \right)\left({m-\gamma} \right)}f} \right)} \right]\left({0^+} \right)$ are evaluated in the limit $t \to 0^+$. Concerning the inverse Laplace transform of functions involving a branch point, we have the following theorem from Titchmarsh (see ^[59]).

Theorem 2.1. Let $\widetilde{f}\left({ \mathbf{s}} \right)$ be an analytic function which has a branch cut on the real negative semiaxis, which has the following properties

for any sector $\left|{\arg\left({ \mathbf{s}} \right)}\right| < \pi-\eta$, where $0 < \eta < \pi$. Then the inverse Laplace transform of $\widetilde{f}\left({ \mathbf{s}} \right)$ is given by

where $\mathrm{Im}\left({\cdot} \right)$ denotes the imaginary part.

The convolution of two integrable functions $f$ and $g$ with support in $[0, +\infty)$ is defined by

and the Convolution Theorem for the Laplace transform is given by

The $n$-dimensional Fourier transform of a real-valued integrable function $f$ in $\mathbb{R}^n$ is defined by (see ^[57])

while the corresponding inverse Fourier transform is given formally by

The convolution operator of two functions in $\mathbb{R}^n$ is defined by the integral

and the Convolution Theorem for the Fourier transform is given by

For the $n$-dimensional Laplace operator $\Delta_x = \sum_{i = 1}^{n} \frac{\partial^2}{\partial x_i^2}$ we have (see formula (1.3.32) in ^[57])

Another integral transform that we use in this work is the Mellin transform. For $f$ locally integrable on $]0, +\infty[$ it is defined by (see ^[57])

and the inverse Mellin transform is given by

The condition for the existence of Eq (2.9) is that $-p < \gamma < -q$ (called the fundamental strip), where $p$, $q$ are the order of $f$ at the origin and $\infty$, respectively. The integration in Eq (2.10) is performed along the imaginary axis and the result does not depend on the choice of $\gamma$ inside the fundamental strip. For more information about this transform and its properties, see e.g., ^[57]. The Mellin convolution between two functions is defined by

and satisfies the Mellin Convolution Theorem (see formula (1.4.40) in ^[57])

The following relation holds (see (1.4.30) in ^[57])

The solution of the time-fractional telegraph equation of distributed order obtained in this work involves the Fox H-function $H_{p, q}^{m, n}$, which is defined, via a Mellin-Barnes type integral, by (see ^[60])

where $a_i, b_j \in \mathbb{C},$ and $\alpha_i, \beta_j \in \mathbb{R}^+,$ for $i = 1, \ldots, p$ and $j = 1, \ldots, q$ and $\mathcal{C}$ is a suitable contour in the complex plane separating the poles of the two factors in the numerator (see ^[60]). The expression of the second-order moment in Section 4.1 is presented in terms of the three parameter Mittag-Leffler function $E_{{\beta_1}, {\beta_2}}^{{\beta_3}}(z)$ (see ^[61]), which is defined, in terms of power series, by

where $\left({\gamma} \right)_k$ is the Pochhammer symbol.

Throughout the paper, we assume that all the involved functions are Laplace and Fourier transformable.

3.   Generalized time-fractional telegraph equation of distributed order

In this work we consider the following generalized time-fractional telegraph equation of distributed order

for given weight functions $b_2\left({\beta, \nu} \right) > 0$ and $b_1\left({\alpha, \mu} \right) > 0$, satisfying

and subject to the following initial and boundary conditions

where $\left({x, \, t} \right) \in \mathbb{R}^n \times \mathbb{R}^+$, $\Delta_x$ is the classical Laplace operator in $\mathbb{R}^n$, the partial time-fractional derivatives of order $\beta \in \, ]1, 2]$ and $\alpha \in \, ]0, 1]$, and types $\mu, \nu \in [0, 1]$ are in the Hilfer sense and given by Eq (2.1), $a \in \mathbb{R}_0^+$, $c \in \mathbb{R} \setminus\{0\}$, $d \in \mathbb{R}$, and $C_1, C_2 \in \mathbb{R}^+$. The positive constants $C_1$ and $C_2$ can be taken as 1 if we assume the normalization condition for the integrals Eq (3.2). Moreover, $q$ belongs to $L_1\left({ \mathbb{R}^n \times I} \right)$, and $f, g_1, g_2 \in L_1\left({ \mathbb{R}^n} \right)$. We look for solutions $u\left({x, t} \right)$ of our problem in the space $C^2\left({ \mathbb{R}^n} \right) \times C^2\left({0, +\infty} \right)$ with possible exception at $x = 0$.

3.1.   Solution in the Fourier-Laplace domain

In order to analytically determine the solution of Eqs (3.1) and (3.2) in the space-time domain we start applying the Fourier and Laplace transforms to Eq (3.1) and solve the equation in the Fourier-Laplace domain. After that, there are two alternative strategies related to the order in carrying out the inversions of the Fourier and Laplace transforms are performed (see ^[2]):

$\rm(S1)$ invert the Fourier transform, yielding $\widetilde{u}\left({x, s} \right)$, and then invert the Laplace transform of the result.

$\rm(S2)$ invert the Laplace transform, yielding $\widehat{u}\left({\kappa, t} \right)$, and then invert the Fourier transform of the result.

In this work, we consider the strategy (S2) where the inversion of the Laplace transform is performed via the classical Titchmarsh's Theorem, and the inversion of the Fourier transform is performed via the Mellin transform.

we start by applying in Eq (3.1) the Laplace transform with respect to the variable $t \in \mathbb{R}^+$ and the $n$-dimensional Fourier transform with respect to the variable $x \in \mathbb{R}^n$. Taking into account relations Eq (2.2) and Eq (2.8), and the initial conditions in Eqs (3.3) and (3.4), we obtain

which is equivalent to

where $\widehat{f}$ and $\widehat{g}_i$ are the Fourier transforms of the functions $f$ and $g_i$, respectively, and

Remark 3.1. When $\mu = \nu = 1$, i.e., the telegraph equation has only Caputo fractional derivatives, we have the following relations between Eqs (3.6) and (3.7)

The previous relations combined with the fact that Eqs (3.3) and (3.4) reduce to

i.e., $f = g_1$ when $\nu = \mu = 1$, allow us to reduce the expression (3.5) to the correspondent one obtained in ^[46].

3.2.   Solution in the space-time domain

In this section we perform the inversion of the Laplace and Fourier transforms in order to obtain our solution in the space-time domain. Let us consider the following auxiliary functions in the Laplace domain

with $p = 0$ or $p = -1.$ To further proceed we make the following additional assumption:

Assumption (H1) holds for the particular cases we consider later on. Applying Theorem 2.1, we have

To evaluate the imaginary parts of the functions $\widehat{\widetilde{u}}_j\left({\kappa, re^{i\pi}} \right)$, $j = 1, 2, 3$, along the ray $\mathbf{s} = re^{i\pi}$, with $r > 0$, we consider the following polar decompositions

After straightforward calculations, we obtain the following expressions

where

Remark 3.2. Taking into account Remark 3.1, when $\mu = \nu = 1$ we have, by straightforward calculations, the following relations

Applying the inverse Laplace transform to Eq (3.5) and taking into account Eqs (3.12), (3.15), (3.16) and (3.17), we obtain

where $\ast_{{t}}$ is given by Eq (2.3) and in the last term me made use of Eq (2.4). For the inversion of the Fourier transform, taking into account Eqs (2.5) and (2.7), we obtain

Using the following formula presented in ^[58] for the inverse Fourier transform of $L_1$-functions

and since we are dealing with radial functions in $\kappa$, Eq (3.20) can be rewritten as

To compute explicitly $\mathbf{I_1}$, $\mathbf{I_2}$, $\mathbf{I_3}$, and $\mathbf{I_4}$ in Eq (3.22) we are going to use the Mellin transform. First, we rewrite these integrals as a Mellin convolution Eq (2.11). In fact, considering the following auxiliary functions

we have for $i = 1, 2, 3, 4$

From the relations Eqs (2.12) and (2.11), we have for $i = 1, 2, 3, 4$

which is equivalent to

Now, we compute the Mellin transforms that appear in (3.24). The Mellin transform of the function $f$ was already calculated in ^[46] (see formula (43)):

To compute the Mellin transform of the function $g_1,$ we take into account Eqs (2.9), (3.15) and (3.18), obtaining

Considering the change of variables $w^2 = z$ in Eq (3.26) we obtain

Integrals $\mathbf{I_5}$ and $\mathbf{I_6}$ were already calculated in ^[46] (see formulas (49) and (50)). Therefore, we have that

and

where

Hence, from Eqs (3.28) and (3.29) we conclude that Eq (3.27) takes the form

Now, we calculate the Mellin transform of $g_2$. Taking into account Eqs (2.9), (3.16)–(3.18) and (3.30) we get

Considering the change of variables $w^2 = z$ in Eq (3.33), we obtain

The two integrals in Eq (3.34) correspond to $\mathbf{I_5}$ and $\mathbf{I_6}$. Hence, from Eqs (3.28) and (3.29) we arrive to

For the calculation of the Mellin transform of $g_3$ we use Eqs (2.9), (3.16)–(3.18) and (3.30) to get

Expression Eq (3.36) is very similar to Eq (3.32) with a difference in the power of $-r$. However, this exponent does not affect any of the performed calculations in obtaining Eq (3.35). Then we get

Finally, we calculate the Mellin transform of $g_4$. Taking into account Eqs (2.9), (3.17), (3.18) and (3.30) we get

Considering the change of variables $w^2 = z$ in Eq (3.38), we obtain

The integral in Eq (3.39) corresponds to the integral $\mathbf{I_6}$. Hence, from Eq (3.29) we arrive to

Now, using the inverse Mellin transform Eq (2.10) applied to (3.24), we obtain the representation of the integrals $\mathbf{I_1}$, $\mathbf{I_2}$, $\mathbf{I_3},$ and $\mathbf{I_4}$ in terms of Mellin-Barnes integrals and, consequently, as Fox H-functions. For the integral $\mathbf{I_1}$, taking into account Eqs (2.10), (3.18), (3.31) and (3.25), we obtain

which is equivalent, by Eq (2.13), to the following expression in terms of Fox H-functions

For the integral $\mathbf{I_2}$, taking into account Eqs(2.10), (3.18), (3.35) and (3.25), we obtain

which is equivalent, by Eq (2.13), to the following expression in terms of Fox H-functions

Due to the similarities between $g_2$ and $g_3$ (and consequently between $\mathbf{I_2}$ and $\mathbf{I_3}$) we have that

Finally, for the integral $\mathbf{I_4}$, taking into account Eqs (2.10), (3.18), (3.40) and (3.25), we obtain

which is equivalent, by (2.13), to the following expression in terms of Fox H-functions

From Eqs (3.41)–(3.44) we conclude that the representation (3.22) of the solution $u\left({x, t} \right)$ of Eqs (3.1) and (3.2) corresponds to the sum of convolution integrals involving Fox H-functions.

In the next subsection we summarize our calculations in the main result of the paper.

3.3.   Main result and corollary

Taking into account Eqs (3.22) and (3.41)–(3.44) we obtain our main result.

Theorem 3.3. The solution of the time-fractional telegraph equation of distributed order Eq (3.1) subject to the conditions Eqs (3.3) and (3.2) and the additional assumption Eq (3.11) is given, in terms of convolution integrals, by

where the functions $G_1$, $G_2$, $G_3$ and $G_4$ are given by

where $\rho_1^\ast,$ $\gamma_1^\ast$, $\rho_2^\ast,$ and $\gamma_2^\ast$, $A$ and $B$, and $\psi$ are given, respectively, by Eq (3.14), (3.18) and (3.30). Moreover, the functions $\mathcal{H}$ and $\mathcal{H}^\ast$ are expressed in terms of the following Fox H-functions

Remark 3.4. If we consider

with $\lambda \in \mathbb{R}^+$ in Eqs (3.1) and (3.2), then the solution $u\left({x, t} \right)$ given by Eq (3.45) corresponds to the eigenfunctions of the generalized time-fractional telegraph equation of distributed order in $\mathbb{R}^n \times \mathbb{R}^+$. Moreover, if additionally $b_2\left({\beta, \nu} \right) = 0$ (resp. $b_1\left({\alpha, \mu} \right) = 0$) we obtain the representation of the generalized eigenfunctions of the time-fractional diffusion (resp. wave) equation of distributed order in $\mathbb{R}^n \times \mathbb{R}^+$.

Putting $a = 0$ in Theorem 3.3, we have the following simplifications

which give the following result.

Corollary 3.5. The solution of the generalized time-fractional wave equation of distributed order in $\mathbb{R}^n \times \mathbb{R}^+$

for a given density function $b_2\left({\beta, \nu} \right)$, subject to the following initial and boundary conditions

is given, in terms of convolution integrals, by

where the functions $G_2$, $G_3$, and $G_4$ are given by

with $\rho$ and $\gamma$, $\rho^\ast$ and $\gamma^\ast$ given by Eqs (3.13) and (3.14), respectively, and the functions $\mathcal{H}$ and $\mathcal{H}^\ast$ are expressed in terms of the following Fox H-functions

Remark 3.6. If we consider $\nu = \mu = 1$ in Theorem 3.3 (i.e., the telegraph equation has only Caputo fractional derivatives) we obtain the main result in ^[46]. For that we need to take into account Remarks 3.1 and 3.2 and to combine the first two integrals in Eq (3.45) into a unique integral. Hence, the solution $u\left({x, t} \right)$ is given by

where for $\nu = \mu = 1$

By Eq (3.18) we have that

and, hence, Eq (3.47) simplifies to

which corresponds to the function $G_1$ presented in the main result of ^[46]. Moreover, in Eq (3.46)

which correspond, respectively, to the functions $G_2$ and $G_3$ that appear in the main result of ^[46]. Therefore, we can claim that there is consistency in our results.

Remark 3.7. The telegraph equations studied in ^[38,40] are particular cases of the equation studied in this paper, for the choices $\nu = \mu = 1,$ $b_2\left({\beta, \nu} \right) = \delta\left({\beta-\beta_1} \right),$ $b_1\left({\alpha, \mu} \right) = \delta\left({\alpha-\alpha_1} \right)$, with $1 < \beta_1 \leq 2$ and $0 < \alpha_1 \leq 1$ $d = 0,$ and $q\left({x, t} \right) = 0.$

The numerical implementation of Eq (3.45) is possible, however, depends substantially on the study of the asymptotic behaviour of $G_1$, $G_2,$ $G_3$ and $G_4$ through the study of the asymptotic behaviour of the associated Fox H-functions. We would like to remark also that Eq (3.45) is a very general solution, but for particular cases of the dimension, of the fractional parameters, and/or of the density functions, it is possible to get simpler expressions.

4.   Fractional moments

In this section we obtain the expression for some fractional moments of the first fundamental solution $\mathbf{G_1}$ of the time-fractional telegraph equation of distributed order Eq (3.1) with $d \, = q\left({x, t} \right) \, = 0,$ and subject to the initial and boundary conditions

Then $\mathbf{G_1}$ is given by $\mathbf{G_1}\left({x, t} \right) = G_1\left({x, t} \right)+G_2\left({x, t} \right)$, where $G_1$ and $G_2$ are given in Theorem 3.3.

It is well known that the Mellin transform (2.9) can be interpreted as the fractional moment of order $s-1$ of the function $f$ (see ^[62]). Therefore, we can calculate the fractional moments $\mathbf{M}_{{n; , \gamma}}^{{\alpha, \mu; \beta, \nu}}$ of arbitrary order $\gamma > 0$ of $\mathbf{\widetilde{G}_1}$, where $\mathbf{\widetilde{G}_1}$ denotes the Laplace transform of $\mathbf{G_1}$. Denoting by $\mathbf{s}$ the variable in the Laplace domain and by $r$ the radial quantity $\left|{x}\right|$, we have, from the definition of the Mellin transform, that

From Eq (3.5), we have that

which is equivalent to

To calculate the inverse Fourier transform we are going to use the Mellin transform, similarly as it was done in Section 3. Taking into account Eq (3.21), we have that

where $\ast_{\mathcal{M}}$ denotes the Mellin convolution given by Eq (2.11) at the point $\frac{1}{r}$ with

and

Denoting by $\mathbf{I_7}$ the integral in Eq (4.3), we have, by relations Eqs (2.12) and (2.11), that

which is equivalent to

From Eq (3.25) we have that

Now, we calculate the Mellin transform of the function $g_5$. Taking into account Eqs (2.9) and (3.8) with $p = 0$, Eqs (3.9) and (3.10) with $p = -1$, we get

The integral in the previous expression was already calculated in ^[46] (see formula (84)). Therefore, we have that

From Eqs (4.5), (4.6), (4.4) and (4.2) we conclude that

By the duplication formula of the Gamma function $\Gamma\left({2z} \right) \, = \frac{2^{2z-1}}{\sqrt{\pi}} \, \Gamma\left({z} \right) \, \Gamma\left({z+\frac{1}{2}} \right)$, we have the following equalities for the Gamma functions that appear in (4.7)

Taking into account Eqs (4.8)–(4.10), expression Eq (4.7) simplifies to

and, consequently, the fractional moments of arbitrary order $\gamma$ in the Laplace domain are given by

If we restrict (4.11) to the time-fractional telegraph equation of distributed order with Caputo fractional derivatives, i.e., if we consider $\mu = \nu = 1$ (which implies that $B_1^\ast\left({ \mathbf{s}} \right) = \mathbf{s}^{-1} B_1\left({ \mathbf{s}} \right)$ and $B_2^\ast\left({ \mathbf{s}} \right) -\frac{d^2}{c^2} = \mathbf{s}^{-2}\left[{B_2\left({ \mathbf{s}} \right)-\frac{d^2}{c^2}} \right]$ by Remark 3.1), then (4.11) becomes equal to

which coincides with the correspondent expression deduced in ^[46]. Let us now analyse expression Eq (4.11) for some special cases:

● When $\gamma = n-2k-3$, with $n > 2k+3$ and $k \in \mathbb{N}_0$, the correspondent moments in the Laplace domain become infinite.

● When $\gamma = 1$ (mean value), we have

which becomes infinite when $n = 4+2k$, with $k \in \mathbb{N}_0$.

● When $\gamma = 2$ (variance), we have

which becomes infinite when $n = 5+2k,$ with $k \in \mathbb{N}_0$.

4.1.   Tauberian analysis for the second-order moment (variance)

When the fundamental solution is a positive function, it is possible to classify the diffusion process by analysing the correspondent second-order moment, also called the mean squared displacement of a particle. This is obtained by comparison with the variance in the normal diffusion process. In this subsection we consider the particular case of double-order distributed fractional derivatives, and we analyze the second-order moment for short and long-time. We separate our analysis between the diffusion and the wave cases. The following Laplace inversion formulas will be needed in the sequel:

● Formula (2.1.1.1) in ^[63]

● Formula (5.1.26) in ^[61]

where $E_{\alpha, \beta}^\gamma$ is the three parameter Mittag-Leffler function given by (2.14).

4.1.1.   Two composite time-fractional derivatives in the diffusion case

Here we consider $b_2\left({\beta, \nu} \right) \, = 0$, which implies that $B_2\left({ \mathbf{s}} \right) = B_2^\ast\left({ \mathbf{s}} \right) \, = 0$. In this case, the second-order moment in the Laplace domain becomes

Further, we assume

with $0 < \alpha_1 < \alpha_2 \leq 1$, $0 \leq \mu_1, \mu_2 \leq 1,$ $\mu_2 < \mu_1 \frac{1-\alpha_1}{1-\alpha_2}$, $k_1, k_2 > 0$, and $k_1+k_2 = 1$. For this $b_1\left({\alpha, \mu} \right)$ we get

Considering Eq (4.18) in Eq (4.16) we get

To invert the Laplace transform of $\widetilde{\mathbf{M}}_{{n; 2}}^{{\left({\alpha_1, \alpha_2} \right), \left({\mu_1, \mu_2} \right)}}$ we first rearrange the expression Eq (4.19):

Taking into account Eq (4.15) with

we get from Eq (4.20) that

Remark 4.1. For $n = 1,$ the expression Eq (4.21) reduces to the expression (28) in ^[51] with suitable identification of the parameters, which indicates consistency in our results.

The graphical representation of $\mathbf{M}_{{n; 2}}^{{\left({\alpha_1, \alpha_2} \right), \left({\mu_1, \mu_2} \right)}}$ as a function of $t$ using Eq (4.21) is not an easy procedure due to the presence of the three parameter Mittag-Leffler function $E_{\alpha, \beta}^\gamma$. The numerical implementation of this special function is possible for some range of the parameter $\alpha$ (see ^[64]), which do not include all the cases studied in this work. Therefore, we will make an asymptotic analysis of (4.19) using the Tauberian analysis. Since $\mathbf{s}^{-\mu_2\left({1-\alpha_2} \right)+\mu_1\left({1-\alpha_1} \right)} \to 0$ and $\mathbf{s}^{\alpha_2-\alpha_1} \to 0$ as $\mathbf{s} \to 0$ then, we get the following asymptotic behaviour of $\widetilde{\mathbf{M}}_{{n; 2}}^{{\left({\alpha_1, \alpha_2} \right), \left({\mu_1, \mu_2} \right)}}\left({s} \right)$ as $\mathbf{s} \to 0:$

Concerning the symbol $\sim$ in the previous and subsequent expressions, we say that $f$ and $g$ are asymptotically equivalent as $w \to \infty$ (resp. as $w \to 0$), i.e., $f \sim g$, if and only if $\lim_{w \to \infty} \frac{f\left({w} \right)}{g\left({w} \right)} = 1$ (resp. $\lim_{w \to 0} \frac{f\left({w} \right)}{g\left({w} \right)} = 1$).

Using (4.14) to invert the Laplace transform in (4.22), we obtain the asymptotic behavior of $\mathbf{M}_{{n; 2}}^{{\left({\alpha_1, \alpha_2} \right), \left({\mu_1, \mu_2} \right)}}$ for $t \to +\infty,$ with $0 \leq \mu_1 \leq 1$ and $n$ and $\alpha_1$ according to the following cases:

To classify the type of diffusion process we need to compare $\mathbf{M}_{{n; 2}}^{{\left({\alpha_1, \alpha_2} \right), \left({\mu_1, \mu_2} \right)}}$ with the moment of the normal diffusion process which is given by

According with the dimension $n$ we have the following cases:

● $n = 1, 2, 3$: $\mathbf{M}_{{n; 2}}^{{\left({\alpha_1, \alpha_2} \right), \left({\mu_1, \mu_2} \right)}}\left({t} \right)/\mathbf{M}_{{n; 2}}^{{\left({1, 1} \right), \left({\mu_1, \mu_2} \right)}}\left({t} \right) \to 0$, as $t \to +\infty,$ for all $0 \leq \mu_1 < 1,$ and all admissible $\alpha_1$ and $n,$ which corresponds to a subdiffusion process in the long time. In the limit case, $\mu_1 = 1$ (Caputo case), for $n = 1, 2$ we still have a subdiffusion process while for $n = 3$ holds $\mathbf{M}_{{3;2}}^{{\left({\alpha_1, \alpha_2} \right), \left({\mu_1, \mu_2} \right)}}\left({t} \right)/\mathbf{M}_{{3;2}}^{{\left({1, 1} \right), \left({\mu_1, \mu_2} \right)}}\left({t} \right) \to k > 0$, as $t \to +\infty,$ thus the process coincides with the normal diffusion in the long time;

● $n = 4$: The classification of the type of diffusion depends on the type $\mu_1.$

▶ For $0\leq \mu_1 < 1/2$ holds $\mathbf{M}_{{4;2}}^{{\left({\alpha_1, \alpha_2} \right), \left({\mu_1, \mu_2} \right)}}\left({t} \right)/\mathbf{M}_{{4;2}}^{{\left({1, 1} \right), \left({\mu_1, \mu_2} \right)}}\left({t} \right) \to 0$, as $t \to +\infty,$ thus corresponding to a subdiffusion process in the long time;

▶ For $\mu_1 = 1/2$ holds $\mathbf{M}_{{4;2}}^{{\left({\alpha_1, \alpha_2} \right), \left({\mu_1, \mu_2} \right)}}\left({t} \right)/\mathbf{M}_{{4;2}}^{{\left({1, 1} \right), \left({\mu_1, \mu_2} \right)}}\left({t} \right) \to k > 0$, as $t \to +\infty,$ thus the process coincides with the normal diffusion in the long time;

▶ For $1/2 < \mu_1 \leq 1$ holds $\mathbf{M}_{{4;2}}^{{\left({\alpha_1, \alpha_2} \right), \left({\mu_1, \mu_2} \right)}}\left({t} \right)/\mathbf{M}_{{4;2}}^{{\left({1, 1} \right), \left({\mu_1, \mu_2} \right)}}\left({t} \right) \to +\infty$, as $t \to +\infty,$ thus corresponding to a superdiffusion process in the long time.

● $n = 6+2k$: The moment $\mathbf{M}_{{n; 2}}^{{\left({\alpha_1, \alpha_2} \right), \left({\mu_1, \mu_2} \right)}}\left({t} \right)$ is positive when $n = 8+4k$ and negative when $n = 6+4k,$ $k \in \mathbb{Z}^+$. This is due to the change of sign of the gamma function in the numerator. In the second case, a probabilistic interpretation is no longer possible.

Finally, we note a different behaviour of the monotonicity of the second-order moment along the dimensions, depending on the values assumed by the several parameters in the expression. When $n = 1, 2, 3, 4$, if we assume that $c > 0$, then the following conclusions about $\mathbf{M}_{{n; 2}}^{{\left({\alpha_1, \alpha_2} \right), \left({\mu_1, \mu_2} \right)}}$ for large values of $t$ can be drawn:

● $n = 1$: $\mathbf{M}_{{n; 2}}^{{\left({\alpha_1, \alpha_2} \right), \left({\mu_1, \mu_2} \right)}}$ is an increasing function when $\frac{1-\mu_1}{2-\mu_1} < \alpha_1 \leq 1$, is constant when $\alpha_1 = \frac{1-\mu_1}{2-\mu_1}$, and is a decreasing function when $0 < \alpha_1 < \frac{1-\mu_1}{2-\mu_1} \leq 0.5$, with $0 \leq \mu_1 \leq 1$.

● $n = 2$: $\mathbf{M}_{{n; 2}}^{{\left({\alpha_1, \alpha_2} \right), \left({\mu_1, \mu_2} \right)}}$ is an increasing function when $\frac{2-2\mu_1}{3-2\mu_1} < \alpha_1 \leq 1$, is constant when $\alpha_1 = \frac{2-2\mu_1}{3-2\mu_1}$, and is a decreasing function when $0 < \alpha_1 < \frac{2-2\mu_1}{3-2\mu_1} \leq \frac{2}{3}$, with $0 \leq \mu_1 \leq 1$.

● $n = 3$: $\mathbf{M}_{{n; 2}}^{{\left({\alpha_1, \alpha_2} \right), \left({\mu_1, \mu_2} \right)}}$ is a decreasing function when $0 < \alpha_1 < 1$ and $0 \leq \mu_1 < 1,$ and is constant when $\alpha_1 = 1$ or $0 \leq \alpha_1 \leq 1$ and $\mu_1 = 1$.

● $n = 4$: $\mathbf{M}_{{n; 2}}^{{\left({\alpha_1, \alpha_2} \right), \left({\mu_1, \mu_2} \right)}}$ is a decreasing function for all $0 < \alpha_1 \leq 1$ and $0 \leq \mu_1 < 1$, and is constant when $\alpha_1 = 0$ and $\mu_1 = 1$.

Now, we study the asymptotic behaviour of $\mathbf{M}_{{n; 2}}^{{\left({\alpha_1, \alpha_2} \right), \left({\mu_1, \mu_2} \right)}}$ for small values of $t.$ From Eq (4.19) we have the following asymptotic behaviour when $\mathbf{s} \to +\infty$:

Using Eq (4.14) to invert the Laplace transform in Eq (4.25), we obtain the asymptotic behavior of $\mathbf{M}_{{n; 2}}^{{\left({\alpha_1, \alpha_2} \right), \left({\mu_1, \mu_2} \right)}}$ for $t \to 0^+$, with $0 \leq \mu_2 \leq 1$ and $n$ and $\alpha_2$ according to the following cases:

Comparing $\mathbf{M}_{{n; 2}}^{{\left({\alpha_1, \alpha_2} \right), \left({\mu_1, \mu_2} \right)}}$ with $\mathbf{M}_{{n; 2}}^{{\left({1, 1} \right), \left({\mu_1, \mu_2} \right)}}$ when $t \to 0^+$, we have the following conclusions:

● $n = 1, 2, 3$: $\mathbf{M}_{{n; 2}}^{{\left({\alpha_1, \alpha_2} \right), \left({\mu_1, \mu_2} \right)}}\left({t} \right)/\mathbf{M}_{{n; 2}}^{{\left({1, 1} \right), \left({\mu_1, \mu_2} \right)}}\left({t} \right) \to +\infty$, as $t \to 0^+$, for all $0 \leq \mu_1 < 1,$ and all admissible $\alpha_1$ and $n,$ which corresponds to a superdiffusion process in the short time. In the limit case, $\mu_1 = 1$ (Caputo case), for $n = 1, 2$ we still have a superdiffusion process while for $n = 3$ holds $\mathbf{M}_{{3;2}}^{{\left({\alpha_1, \alpha_2} \right), \left({\mu_1, \mu_2} \right)}}\left({t} \right)/\mathbf{M}_{{3;2}}^{{\left({1, 1} \right), \left({\mu_1, \mu_2} \right)}}\left({t} \right) \to k > 0$, as $t \to +\infty,$ thus the process coincides with the normal diffusion in the short time;

● $n = 4$: The classification of the type of diffusion depends on the parameter $\mu_1.$

▶ For $0\leq \mu_1 < 1/2$ holds $\mathbf{M}_{{4;2}}^{{\left({\alpha_1, \alpha_2} \right), \left({\mu_1, \mu_2} \right)}}\left({t} \right)/\mathbf{M}_{{4;2}}^{{\left({1, 1} \right), \left({\mu_1, \mu_2} \right)}}\left({t} \right) \to +\infty$, as $t \to 0^+$, thus corresponding to a superdiffusion process in the short time;

▶ For $\mu_1 = 1/2$ holds $\mathbf{M}_{{4;2}}^{{\left({\alpha_1, \alpha_2} \right), \left({\mu_1, \mu_2} \right)}}\left({t} \right)/\mathbf{M}_{{4;2}}^{{\left({1, 1} \right), \left({\mu_1, \mu_2} \right)}}\left({t} \right) \to k > 0$, as $t \to 0^+$, thus the process coincides with the normal diffusion in the short time;

▶ For $1/2 < \mu_1 \leq 1$ holds $\mathbf{M}_{{4;2}}^{{\left({\alpha_1, \alpha_2} \right), \left({\mu_1, \mu_2} \right)}}\left({t} \right)/\mathbf{M}_{{4;2}}^{{\left({1, 1} \right), \left({\mu_1, \mu_2} \right)}}\left({t} \right) \to 0$, as $t \to 0^+$, thus corresponding to a subdiffusion process in the short time.

● $n = 6+2k$: As happen in the long time case the moment is not always positive, and hence also here it is not possible to perform a probabilistic interpretation for all the values of $n$.

4.1.2.   Two composite time-fractional derivatives in the wave case

Here we consider $b_1\left({\alpha, \mu} \right) \, = 0$. Hence, $B_1\left({ \mathbf{s}} \right) = \, B_1^\ast\left({ \mathbf{s}} \right) \, = 0$ and the second-order moment in the Laplace domain becomes

Let us now assume that

with $1 < \beta_1 < \beta_2 \leq 2$, $0 \leq \nu_1, \nu_2 \leq 1$, $\nu_2 < \nu_1 \frac{2-\beta_1}{2-\beta_2}$, $k_1, k_2 > 0$, and $k_1+k_2 = 1$. For this $b_2\left({\beta, \nu} \right)$ we get

Considering Eq (4.28) in Eq (4.27) we get

Inverting the Laplace transform of $\widetilde{\mathbf{M}}_{{n; 2}}^{{\left({\beta_1, \beta_2} \right), \left({\nu_1, \nu_2} \right)}}\left({s} \right)$ following the same steps of the deduction of Eq (4.21) we get

First, we study the asymptotic behaviour of $\mathbf{M}_{{n; 2}}^{{\left({\beta_1, \beta_2} \right), \left({\nu_1, \nu_2} \right)}}$ for large values of $t.$ From Eq (4.19) we have the following asymptotic behaviour when $\mathbf{s} \to 0:$

Using Eq (4.14) to invert the Laplace transform in Eq (4.31), we obtain the asymptotic behavior of $\mathbf{M}_{{n; 2}}^{{\left({\beta_1, \beta_2} \right), \left({\nu_1, \nu_2} \right)}}$ for $t \to + \infty,$ with $\beta_1,$ $\nu_1,$ and $n$ according to the following cases:

To classify the type of diffusion-wave process we need to compare $\mathbf{M}_{{n; 2}}^{{\left({\beta_1, \beta_2} \right), \left({\nu_1, \nu_2} \right)}}$ with the normal diffusion process. The situation here depends on the dimension and we have the following cases:

● $n = 1$: $\mathbf{M}_{{1;2}}^{{\left({\beta_1, \beta_2} \right), \left({\nu_1, \nu_2} \right)}}\left({t} \right)/ \mathbf{M}_{{1;2}}^{{\left({1, 1} \right), \left({\mu_1, \mu_2} \right)}}(t) \to +\infty$, as $t \to +\infty,$ for $1 < \frac{3-2\nu_1}{2-\nu_1} < \beta_1 < 2$ and $0\leq \nu_1\leq 1,$ which corresponds to a superdiffusion process in the long time and $\mathbf{M}_{{1;2}}^{{\left({\beta_1, \beta_2} \right), \left({\nu_1, \nu_2} \right)}}\left({t} \right)/ \mathbf{M}_{{1;2}}^{{\left({1, 1} \right), \left({\mu_1, \mu_2} \right)}} \to 0$, as $t \to +\infty,$ for $1 < \beta_1 < \frac{3-2\nu_1}{2-\nu_1} < 2$ and $0\leq \nu_1\leq 1,$ which corresponds to a subdiffusion process in the long time. In the special case $\beta_1 = \frac{3-2\nu_1}{2-\nu_1}$ the process coincides with the normal diffusion in the long time.

● $n = 2$: $\mathbf{M}_{{n; 2}}^{{\left({\beta_1, \beta_2} \right), \left({\nu_1, \nu_2} \right)}}\left({t} \right)/ \mathbf{M}_{{2;2}}^{{\left({1, 1} \right), \left({\mu_1, \mu_2} \right)}}(t) \to +\infty$, as $t \to +\infty,$ for $1 < \frac{5-4\nu_1}{3-2\nu_1} < \beta_1 < 2$ and $0\leq \nu_1\leq 1,$ which corresponds to a superdiffusion process in the long time and $\mathbf{M}_{{n; 2}}^{{\left({\beta_1, \beta_2} \right), \left({\nu_1, \nu_2} \right)}}\left({t} \right)/ \mathbf{M}_{{2;2}}^{{\left({1, 1} \right), \left({\mu_1, \mu_2} \right)}}(t) \to 0$, as $t \to +\infty,$ for $1 < \beta_1 < \frac{5-4\nu_1}{3-2\nu_1} < 2$ and $0\leq \nu_1\leq 1,$ which corresponds to a subdiffusion process in the long time. In the special case $\beta_1 = \frac{5-4\nu_1}{3-2\nu_1}$ the process coincides with the normal diffusion in the long time.

● $n = 3$ $\mathbf{M}_{{3;2}}^{{\left({\beta_1, \beta_2} \right), \left({\nu_1, \nu_2} \right)}}\left({t} \right)/ \mathbf{M}_{{3;2}}^{{\left({1, 1} \right), \left({\mu_1, \mu_2} \right)}}(t) \to 0$, as $t \to +\infty,$ for all $1 < \beta_1\leq 2$ and $0 \leq \nu_1 < 1$, thus corresponding to a subdiffusion process. For $\nu_1 = 1$ the process coincides with the normal diffusion case.

● $n = 4$: $\mathbf{M}_{{4;2}}^{{\left({\beta_1, \beta_2} \right), \left({\nu_1, \nu_2} \right)}}\left({t} \right)/ \mathbf{M}_{{4;2}}^{{\left({1, 1} \right), \left({\mu_1, \mu_2} \right)}}(t) \to 0$, as $t \to +\infty,$ for all $1 < \beta_1\leq 2$ and $1/2 < \nu_1 \leq 1$, thus corresponding to a subdiffusion process.

Finally, we study the asymptotic behaviour of $\mathbf{M}_{{n; 2}}^{{\left({\beta_1, \beta_2} \right), \left({\nu_1, \nu_2} \right)}}$ for small values of $t$ knowing the asymptotic behaviour of Eq (4.29) when $\mathbf{s} \to +\infty.$ From (4.29), as $\mathbf{s} \to +\infty,$ we have

Using Eq (4.14) to invert the Laplace transform in Eq (4.33), we obtain the asymptotic behavior of $\mathbf{M}_{{n; 2}}^{{\left({\beta_1, \beta_2} \right), \left({\nu_1, \nu_2} \right)}}$ for $t \to 0^+$, with $\beta_2,$ $\nu_2,$ and $n$ according to the following cases:

The analysis of Eq (4.34) is similar to the one performed for Eq (4.32). Regarding the classification of the diffusion-wave process the following conclusions can be taken:

● $n = 1$: $\mathbf{M}_{{1;2}}^{{\left({\beta_1, \beta_2} \right), \left({\nu_1, \nu_2} \right)}}\left({t} \right)/ \mathbf{M}_{{1;2}}^{{\left({1, 1} \right), \left({\mu_1, \mu_2} \right)}}(t) \to 0$, as $t \to 0^+$, for $1 < \frac{3-2\nu_1}{2-\nu_1} < \beta_2 < 2$ and $0\leq \nu_2 \leq 1$, while $\mathbf{M}_{{n; 2}}^{{\left({\beta_1, \beta_2} \right), \left({\nu_1, \nu_2} \right)}}\left({t} \right)/ \mathbf{M}_{{1;2}}^{{\left({1, 1} \right), \left({\mu_1, \mu_2} \right)}}(t) \to +\infty$, as $t \to 0^+$, for $1 < \beta_2 < \frac{3-2\nu_1}{2-\nu_1} < 2$ and $0\leq \nu_2 \leq 1.$ Hence, in the short time, the process is subdiffusive in the first case and is superdiffusive in the second case. In the special case $\beta_1 = \frac{3-2\nu_1}{2-\nu_1}$ the process coincides with the normal diffusion in the short time.

● $n = 2$: $\mathbf{M}_{{2;2}}^{{\left({\beta_1, \beta_2} \right), \left({\nu_1, \nu_2} \right)}}\left({t} \right)/ \mathbf{M}_{{2;2}}^{{\left({1, 1} \right), \left({\mu_1, \mu_2} \right)}}(t) \to 0$, as $t \to 0^+$, for $1 < \frac{5-4\nu_1}{3-2\nu_1} < \beta_2 < 2$ and $0\leq \nu_2 \leq 1,$ which corresponds to a subdiffusion process in the short time, and $\mathbf{M}_{{2;2}}^{{\left({\beta_1, \beta_2} \right), \left({\nu_1, \nu_2} \right)}}\left({t} \right)/ \mathbf{M}_{{2;2}}^{{\left({1, 1} \right), \left({\mu_1, \mu_2} \right)}}(t) \to +\infty$, as $t \to 0^+$, for $1 < \beta_2 < \frac{5-4\nu_1}{3-2\nu_1} < 2$ and $0\leq \nu_2 \leq 1,$ thus corresponding to a superdiffusion process in the short time. In the special case $\beta_1 = \frac{5-4\nu_1}{3-2\nu_1}$ the process coincides with the normal diffusion in the short time.

● $n = 3$: $\mathbf{M}_{{3;2}}^{{\left({\beta_1, \beta_2} \right), \left({\nu_1, \nu_2} \right)}}\left({t} \right)/ \mathbf{M}_{{3;2}}^{{\left({1, 1} \right), \left({\mu_1, \mu_2} \right)}}(t) \to +\infty$, as $t \to 0^+$, for all $1 < \beta_2\leq 2$ and $0 \leq \nu_2 < 1$, thus corresponding to a superdiffusion process. For $\nu_2 = 1$ the process coincides with the normal diffusion case.

● $n = 4$: $\mathbf{M}_{{4;2}}^{{\left({\beta_1, \beta_2} \right), \left({\nu_1, \nu_2} \right)}}\left({t} \right)/ \mathbf{M}_{{4;2}}^{{\left({1, 1} \right), \left({\mu_1, \mu_2} \right)}}(t) \to + \infty$, as $t \to 0^+$, for all $1 < \beta_2\leq 2$ and $1/2 < \nu_2 \leq 1$, thus corresponding to a superdiffusion process in the short time.

Remark 4.2. If we consider the Caputo case in Sections 4.1.1 and 4.1.2, our analysis of the diffusion-wave process for the double-order case improve the results presented in ^[46]. Moreover, considering single order derivatives and the one-dimensional case, it was proved in ^[65] that the fundamental solution corresponds to a probability density function only when the fractional derivatives are in the Caputo sense.

4.1.3.   Graphical representations of the second-order moment for the long time

In this section we present and analyse the plots of the asymptotic behaviour of $\mathbf{M}_{{n; 2}}^{{\left({\alpha_1, \alpha_2} \right), \left({\mu_1, \mu_2} \right)}}\left({t} \right)$, when $t \to +\infty$, for some of the cases studied previously separating the diffusion and the wave cases. The plots were generated using Mathematica software and the commands available in it.

The diffusion case: In the following figures, we show the graphical representation of Eq (4.23) for $n = 1, 2, 3, 4$, $\alpha_1 = 0.25, 0.50, 0.75$, and different values $\mu_1,$ using a logarithmic scale in the axes when needed.

Looking at the plots we see that the classification of the diffusion process in each dimension agrees with the analysis of Eq (4.23) performed previously. The plots show an interpolation between the extreme cases $\mu_1 = 0$ and $\mu_1 = 1,$ which correspond to the Riemann-Liouville (RL) and Caputo cases, respectively. The extreme cases have different behaviour, e.g., the slope of the variance is different and in the dimension, $n = 3$ the variance is constant in the Caputo case. In contrast, in the RL case, the variance decreases for large values of $t.$ Moreover, for $\alpha_1 = 0.25$ we can observe a different behaviour of the diffusion in dimensions $n = 1$ and $n = 2:$ in the RL case the variance decreases for large values of $t$ while in the Caputo case the variance increases.

The wave case: In the following figures, we show the graphical representation of Eq (4.32) for $n = 1, 2, 3, 4$, $\beta_1 = 1.25, 1.50, 1.75$ and different values $\nu_1$.

For each dimension and different values of the fractional parameters, the process classification agrees with the analysis of Eq (4.32). The range of the plots increases with the increase of $\beta_1$ and $\nu_1$. Again it is interesting to observe the different behaviour of the variance for the extreme cases $\nu_1 = 0$ (RL case) and $\nu_1 = 1$ (Caputo case) for the dimension $n = 3.$ Also the slope of the variance is different in the other dimensions.

4.1.4.   Graphical representations of the second-order moment for the short time

In this section we present and analyse the graphical representation of the asymptotic behaviour of $\mathbf{M}_{{n; 2}}^{{\left({\alpha_1, \alpha_2} \right), \left({\mu_1, \mu_2} \right)}}\left({t} \right)$, when $t \to 0^+$, for some of the cases studied previously separating the diffusion and the wave cases. The plots were generated using Mathematica software and the commands available in it.

The diffusion case: In the following figures, we show the graphical representation of (4.26) for $n = 1, 2, 3, 4$, $\alpha_2 = 0.25, 0.50, 0.75$, and different values $\mu_2.$

Looking at the plots we see that the range of the plots decreases with the increase of $\alpha_2$ and $\mu_2$. The type of process is in accordance with the conclusion made in the analysis of Eq (4.26). Again, we can observe a different behaviour of the variance for small values of $t$ in the extreme cases $\nu_2 = 0$ and $\nu_2 = 1,$ corresponding to the RL and Caputo cases, respectively. For $\alpha_1 = 0.25$ and the dimensions $n = 1$ and $n = 2$ the variance decreases in the RL case and increases in the Caputo case, for small values of $t.$ In the case $\nu_2 = 1$ the plots coincide with the correspondent ones obtained in ^[46].

The wave case: In the following figures, we have the graphical representation of Eq (4.34) for $n = 1, 2, 3, 4$, $\beta_2 = 1.25, 1.50, 1.75$, and different values $\nu_2$ and $n.$

Looking at the plots we see that the conclusions are similar to those already exposed. The behaviour of the functions is in accordance with the conclusions made in the analysis of (4.34). When $\nu_2 = 1$ the plots coincide with those presented in ^[46].

5.   Conclusions

The results presented here generalize those obtained in ^[46] by the introduction of the Hilfer derivative, that allows a smooth interpolation between the Riemann-Liouvile and the Caputo fractional derivatives. The solution of the Cauchy problem associated with the telegraph equation was expressed as convolutions with functions that are expressed by Laplace integrals involving Fox H-functions. For particular cases of the equation the solution can be simplified and we showed that we can recover known results presented in the literature, which reveals consistency of our results. The classification of the diffusion-wave process depends, not only on the spatial dimension, but also on the order and type of the derivatives. This is very different from previous works on the literature since in most cases the telegraph equation is studied only for Caputo or Riemann-Liouville fractional derivatives. It would be interesting to consider other types of density functions in addition to those considered in this article, but the calculations would become cumbersome.

Acknowledgments

The work of the authors was supported by Portuguese funds through CIDMA–Center for Research and Development in Mathematics and Applications, and FCT–Fundação para a Ciência e a Tecnologia, within projects UIDB/04106/2020 and UIDP/04106/2020.

N. Vieira was also supported by FCT via the 2018 FCT program of Stimulus of Scientific Employment - Individual Support (Ref: CEECIND/01131/2018).

The authors thank the anonymous reviewers for the careful reading of the manuscript, their suggestions, and additional references that have improved the article.

Conflict of interest

The authors declare there is no conflicts of interest.