Spectral Analysis of Fractional Hyperbolic Diffusion Equations with Random Data

Abstract

The paper studies the fundamental solutions to fractional in time hyperbolic diffusion equation or telegraph equations and their properties. Then it derives the exact solutions of the fractional hyperbolic diffusion equation with random data in terms of series expansions of isotropic in space spherical random fields on the unit sphere. Numerical illustration are presented to illustrate the theoretical results.

1 Introduction

Spherical random fields are very useful for modelling some phenomena in areas such as earth sciences (like, for example, in geophysics and climatology [8, 9, 18, 22, 40, 48]) and cosmology (see, for instance, [45]). In fact, the application of statistical methods in cosmology [5] has become increasingly important due to the many experimental data obtained in recent years [1], and spherical random fields are of particular interest regarding the analysis of Cosmic Microwave Background (CMB) radiation [33]. As well-known [12, 54], the CMB is a spatially isotropic radiation field spread throughout the visible universe, originated around 14 billion years ago, and it is the main source of information we have about the evolution of the universe. The CMB radiation can be mathematically modelled as an isotropic, mean-square continuous spherical random field for which there is a spectral representation by means of spherical harmonics. Consequently, the study of models of temporal evolution of spherical random fields, in addition to its innate theoretical interest, is a problem that may have some practical interest in the study of CMB radiation.

One such model was recently provided by [4], where stochastic hyperbolic diffusion equations were used for the spherical random fields. The hyperbolic heat equation is formally identical to the linear telegraph equation introduced by Heaviside in his study of transmission lines, and it was introduced by Cattaneo [7] in order to impose a bounded speed of propagation for the temperature disturbances, in contrast therefore with the classical parabolic heat equation, that has an unbounded propagation speed. The boundedness of the speed propagation is desirable because the large-scale coherent structures that are observed in the CMB are believed to be the remains of acoustic waves in the plasma universe. In [4] the explicit solution of the model was given in terms of series of elementary functions, and therefore it could be useful for various qualitative and numerical studies. For more details and references on the telegraph equations, or hyperbolic diffusion equation, see [4, 16, 28], among others.

On the other hand, deviations from the standard diffusive behaviour are known to occur in many situations [3, 44, 47]. Among the different models of anomalous behaviour, an interesting one is provided by the use of fractional differential equations (FDE) [11]. The calculus of non-integer order, or simply fractional calculus, showed an increasing interest over the last decades, specially for the modelling of phenomena like, for example, processes involving memory effects (see, for example, [52, 53]), anomalous transport [11], problems with dissipation [31], etc. However, these models are not unique in the sense that there are many different definitions of a fractional derivative in the literature [50], and in some cases these definitions are introduced into the equations on an ad hoc basis. Nevertheless, among these different definitions there is one that deserves to be highlighted, which is the so-called Caputo (or Caputo-Djrbashian) derivative [24]. One interesting thing about FDE using Caputo derivative is that they appear in the formalism of continuous time random walk (CTRW) [25, 46] as being associated to a long-tailed probability distribution function (PDF). Caputo fractional derivative is also a particular case of more general non-local operators—see for instance [26, 51].

Let us remember that the normal diffusion can be modelled in terms of a random walk where jumps are taken in equal time intervals. In the CTRW approach, the time interval between sucessive jumps is replaced by a waiting time PDF w(t) and the length of the jumps is given by a PDF $\lambda (x)$. CTRW gives a very general method to discuss anomalous diffusion processes, and as a result many generalizations of classical results have been obtained (see, for example, [2, 19, 36, 38]). For the case of a fractal time random walk, where we have long-tailed waiting time PDF with asymptotic behaviour of the form $w(t) \sim (\tau /t)^{(1+\alpha )}$ (where $0<\alpha <1$ and $\tau$ is a characteristic time), and the usual Gaussian jump length PDF, the so-called master equation of the CTRW approach reduces to an equation of diffusion type with Caputo fractional derivative replacing the usual first order time derivative (see [39, 43] for details).

It is natural at this point to think of looking to Cattaneo’s approach to the hyperbolic diffusion equation under the perspective of the CTRW approach. The key point in Cattaneo’s approach was to modify the constitute equation (Fick’s law) by introducing a term proportional to the first order derivative of the flux. In [10] Compte and Metzler have discussed possibilities for the generalization of the Cattaneo equation, and one of these possibilities is to consider the CTRW scenario of fractal time random walk, which gives an equation which can be written in terms of Caputo fractional derivatives in the time variable. Another CTRW approach to the fractional Cattaneo/telegraph equation is given in [35].

The present paper is a continuation of the line of research of the work [4]. Our main objective is to study the fundamental solutions to fractional hyperbolic diffusion equation in the time variable using the Caputo derivative, and its properties. The exact solutions of the fractional hyperbolic diffusion equation with random data in terms of series expansions of isotropic in space spherical random fields on the unit sphere are derived, and numerical illustration are presented to illustrate the theoretical results. This paper is also an extension of the results of the papers [17, 41], in which the time-fractional telegraph equations and telegraph processes have been considered—see also [13,14,15], among the others.

Basic results and definitions about spherical isotropic random fields and their spectral and covariance representations are given in Sect. 2. Section 3 derives the solution of the fractional hyperbolic diffusion equation on the sphere. Basic results about spherical isotropic random fields which are solutions of random fractional hyperbolic diffusions and their spectral and covariance representations are given in Sect. 4. Section 5 contains the detailed proofs of the main results. Section 6.1 presents numerical illustration of the theoretical results.

2 Isotropic Spherical Random Fields

This section introduces basic notations and background by reviewing some results in the theory of spherical random fields from the monograph [33](see, also [29, 30, 32, 55]).

Consider a sphere in the three-dimensional Euclidean space

$$\begin{aligned} \mathbb {S}^{2}=\left\{ {\mathbf {x}}\in \mathbb {R}^{3}:\Vert {\mathbf {x}}\Vert =1\right\} \subset \mathbb {R}^{3} \nonumber \end{aligned}$$

with the Lebesgue measure (the area element on the sphere)

$$\begin{aligned} \sigma (d\theta ,d\varphi )=\sin \theta d\theta d\varphi ,\quad 0\le \theta \le \pi ,\quad 0\le \varphi \le 2\pi ,\text { or }(\theta ,\varphi )\in \mathbb {S}^{2}. \end{aligned}$$

A spherical random field on a complete probability space $(\Omega ,\mathcal {F },{\mathbf {P}})$, denoted by

$$\begin{aligned} T=\left\{ T(\theta ,\varphi )=T_{\omega }(\theta ,\varphi ):\Omega \times \mathbb {S}^{2}\rightarrow \mathbb {R}^{1},0\le \theta \le \pi ,\quad 0\le \varphi \le 2\pi ,\ \omega \in \Omega \right\} , \end{aligned}$$

that is $\left\{ T(\theta ,\varphi ),(\theta ,\varphi )\in \mathbb {S} ^{2}\right\}$ or ${\widetilde{T}}=\{{\widetilde{T}}({\mathbf {x}})$ ,  ${\mathbf {x}} \in \mathbb {S}^{2}\},$ is a stochastic function defined on the sphere $\mathbb {S}^{2}.$

The field ${\widetilde{T}}({\mathbf {x}})$ is called isotropic (in the weak sense) on the sphere $\mathbb {S}^{2}$ if $\mathrm {E}{\widetilde{T}}({\mathbf {x}} )^{2}<\infty$ and its first and second-order moments are invariant with respect to the group SO(3) of rotations in $\mathbb {R}^{3},$ i.e.

$$\begin{aligned} {\mathbf {E}}{\widetilde{T}}({\mathbf {x}})={\mathbf {E}}{\widetilde{T}}(g{\mathbf {x}} ),\quad {\mathbf {E}}{\widetilde{T}}({\mathbf {x}}){\widetilde{T}}({\mathbf {y}})=\mathbf { E}{\widetilde{T}}(g{\mathbf {x}}){\widetilde{T}}(g{\mathbf {y}}), \end{aligned}$$

for every $g\in SO(3)$ and ${\mathbf {x}},{\mathbf {y}}\in \mathbb {S}^{2}.$ This is equivalent to saying that the mean ${\mathbf {E}}{\widetilde{T}}({\mathbf {x}})= {\mathbf {E}}T(\theta ,\varphi )=c=constant$ (without loss of generality we assume that $c=0),$ and that the covariance function ${\mathbf {E}}{\widetilde{T}} ({\mathbf {x}}){\widetilde{T}}({\mathbf {y}})={\mathbf {E}}T(\theta ,\varphi )T(\theta ^{\prime },\varphi ^{\prime })$ depends only on the angular distance $\Theta =\Theta _{PQ}$ between the points $P=(\theta ,\varphi )$ and $Q=(\theta ^{\prime },\varphi ^{\prime })$ on $\mathbb {S}^{2}.$

We consider a real-valued second-order spherical random field T that is continuous in the mean-square sense. Note that [34] proved that the covariance function of a measurable finite-variance isotropic random field on the sphere is necessarily everywhere continuous.

Under these conditions, the field T can be expanded in the mean-square sense as a Laplace series (see, [55], p. 73, [30], p. 33, or [34], p. 123):

$$\begin{aligned} T(\theta ,\varphi )=\sum _{l=0}^{\infty }\sum _{m=-l}^{l}a_{lm}Y_{lm}(\theta ,\varphi ), \end{aligned}$$

(1)

where $\{Y_{lm}(\theta ,\varphi )\}$ represents the complex spherical harmonics. The spectral representation (1) can be viewed as a Karhunen–Loève expansion, which converges in the Hilbert space $L_{2}(\Omega \times \mathbb {S}^2,\sin \theta d\theta d\varphi )$, that is,

$$\begin{aligned} \lim _{L\rightarrow \infty }{\mathbf {E}}\left( \int \limits _{\mathbb {S}^2}\left( T(\theta ,\varphi )-\sum _{l=0}^{L}\sum _{m=-l}^{l}Y_{lm}(\theta ,\varphi )a_{lm}\right) ^{2}\sin \theta d\theta d\varphi \right) =0. \end{aligned}$$

According to the Peter–Weyl theorem (see [34], p. 69), the expansion (1) also converges in the Hilbert space $L_{2}(\Omega ),$ for every ${\mathbf {x}}\in \mathbb {S}^2$, that is, for each ${\mathbf {x}}\in \mathbb {S}^2,$

$$\begin{aligned} \lim _{L\rightarrow \infty }{\mathbf {E}}\left( {\widetilde{T}}({\mathbf {x}} )-\sum _{l=0}^{L}\sum _{m=-l}^{l}Y_{lm}({\mathbf {x}})a_{lm}\right) ^{2}=0. \end{aligned}$$

Recall that for $-l\le m\le l$ it holds

$$\begin{aligned} {\tilde{Y}}_{lm}({\mathbf {x}})= & {} Y_{lm}(\theta ,\varphi )=d_{lm}\exp (im\varphi )P_{l}^{m}(\cos \theta ),\\ d_{lm}= & {} (-1)^{m}\left[ \frac{(2l+1)(l-m)!}{4\pi (l+m)!}\right] ^{1/2}, \end{aligned}$$

where $P_{l}^{m}(\cdot )$ denotes the associated Legendre functions with the indices l and m,  and $P_{l}(\cdot )$ is the l-th Legendre polynomial, i.e.

$$\begin{aligned} P_{l}^{m}(u)=(-1)^{m}(1-u^{2})^{m/2}\frac{d^{m}}{dx^{m}}P_{l}(u),\quad P_{l}(u)=\frac{1}{2^{l}l!}\frac{d^{l}}{dx^{l}}(u^{2}-1)^{l}. \end{aligned}$$

(2)

The spherical harmonics have the following properties

$$\begin{aligned}&\int _{0}^{\pi }\int _{0}^{2\pi }Y_{lm}^{*}(\theta ,\varphi )Y_{l^{\prime }m^{\prime }}(\theta ,\varphi )\sin \theta d\varphi d\theta =\delta _{l}^{l^{\prime }}\delta _{m}^{m^{\prime }},\nonumber \\&\quad Y_{lm}^{*}(\theta ,\varphi )=(-1)^{m}Y_{l(-m)}(\theta ,\varphi ), \end{aligned}$$

(3)

$$\begin{aligned}&Y_{lm}(\pi -\theta ,\varphi +\pi )=(-1)^{l}Y_{lm}(\theta ,\varphi ),\nonumber \\&\quad Y_{l0}(0,0)=\sqrt{\frac{2l+1}{4\pi }}P_{l}(1)=\sqrt{\frac{2l+1}{4\pi }} ,\quad Y_{lm}(0,0)=Y_{l0}(0,0)\delta _{l}^{m}, \end{aligned}$$

(4)

where $\delta _{l}^{l^{\prime }}$ is the Kronecker delta function and the symbol * means the complex conjugation. The random coefficients $a_{lm}$ in the Laplace series (1) can be obtained through inversion arguments in the form of mean-square stochastic integrals

$$\begin{aligned} a_{lm}=\int _{0}^{\pi }\int _{0}^{2\pi }T(\theta ,\varphi )Y_{lm}^{*}(\theta ,\varphi )\sin \theta d\theta d\varphi . \end{aligned}$$

(5)

As T is real-valued, then, by the property (3), it holds

$$\begin{aligned} a_{lm}=(-1)^{m}a_{l\,(-m)},\quad l\ge 0,\ -l\le m\le l. \end{aligned}$$

(6)

The field is isotropic if and only if

$$\begin{aligned} {\mathbf {E}}a_{lm}a_{l^{\prime }m^{\prime }}^{*}=\delta _{l}^{l^{\prime }}\delta _{m}^{m^{\prime }}C_{l},\quad -l\le m\le l,\quad -l^{\prime }\le m^{\prime }\le l^{\prime }. \end{aligned}$$

Thus, ${\mathbf {E}}|a_{lm}|^{2}=C_{l},$$m=0,\pm 1,...,\pm l.$

From (1) and (5) we deduce that the covariance function of an isotropic random fields has the following spectral representation

$$\begin{aligned} \Gamma (\cos \Theta )={\mathbf {E}}T(\theta ,\varphi )T(\theta ^{\prime },\varphi ^{\prime })=\frac{1}{4\pi }\sum _{l=0}^{\infty }(2l+1)C_{l}P_{l}(\cos \Theta ), \end{aligned}$$

where

$$\begin{aligned} \sum _{l=0}^{\infty }(2l+1)C_{l}<\infty . \end{aligned}$$

(7)

The series $\left\{ C_{0},C_{1},C_{2},\ldots ,C_{l},\ldots \right\}$ is called the angular power spectrum of the isotropic random field $T(\theta ,\varphi ).$

If $T(\theta ,\varphi )$ is an isotropic Gaussian field, then the coefficients $a_{lm},$$m=-l,\ldots ,l,$$l\ge 1,$ are complex-valued independent Gaussian random variables if $m\not =-m^{\prime },$ with

$$\begin{aligned} {\mathbf {E}}a_{lm}=0,\quad {\mathbf {E}}a_{lm}a_{l^{\prime }m^{\prime }}^{*}=\delta _{m}^{m^{\prime }}\delta _{l}^{l^{\prime }}C_{l}, \end{aligned}$$

(8)

if $C_{l}>0$, or degenerate to zero if $C_{l}=0.$

3 Solution for Non-random Point-Source Spherical Fractional Hyperbolic Diffusion Equation

This section derives the fundamental solutions of non-random fractional in time hyperbolic diffusion equations. The obtained results will be used in Sect. 4 to obtain solutions of diffusion equations with random initial conditions.

The fractional in time Caputo or Caputo-Djrbashian derivatives are defined for a nice function $f(\theta ,\varphi ,t),t\ge 0,0\le \theta \le \pi ,\quad 0\le \varphi \le 2\pi$ (see, i.e., [37]) as

$$\begin{aligned}&\mathrm {D}_{t}^{\alpha }f(\theta ,\varphi ,t)=\frac{\partial ^{\alpha }}{ \partial t^{\alpha }}f(\theta ,\varphi ,t)=\frac{1}{\Gamma (m-\alpha )} \int _{0}^{t}\left[ \frac{1}{(t-\tau )^{1+\alpha -m}}\frac{\partial ^{m}}{\partial \tau ^{m}}f(\theta ,\varphi ,\tau )\right] d\tau , \nonumber \\&\quad m-1<\alpha <m, \end{aligned}$$

(9)

and

$$\begin{aligned} \mathrm {D}_{t}^{m}f(\theta ,\varphi ,t)=\frac{\partial ^{m}}{\partial t^{m}} f(\theta ,\varphi ,t),\text { for }\alpha =m, \end{aligned}$$

(10)

where $m=1,2,\ldots ,$ is an integer, $\alpha >0.$

Consider the following fractional in time hyperbolic diffusion equations, also known as the fractional telegraph equations (see [28, 41]) or relativistic fractional diffusion equation ([16]) on sphere $\mathbb {S}^{2}$

$$\begin{aligned} \frac{1}{c^{2}}\frac{\partial ^{2\alpha }p(\theta ,\varphi ,t)}{\partial t^{2\alpha }}+\frac{1}{D}\frac{\partial ^{\alpha }p(\theta ,\varphi ,t)}{ \partial t^{\alpha }}=k^{2}\Delta _{(\theta ,\varphi )}\;p(\theta ,\varphi ,t),0<\alpha \le 1, \end{aligned}$$

(11)

with the initial conditions

$$\begin{aligned} p(\theta ,\varphi ,t)|_{t=0}=\frac{1}{\sin \theta }\delta (\theta )\delta (\varphi ),\qquad \left. \frac{\partial }{\partial t}p(\theta ,\varphi ,t)\right| _{t=0}=0, \end{aligned}$$

(12)

if $\frac{1}{2}<\alpha \le 1,$ and

$$\begin{aligned} p(\theta ,\varphi ,t)|_{t=0}=\frac{1}{\sin \theta }\delta (\theta )\delta (\varphi ), \end{aligned}$$

(13)

if $0<\alpha \le \frac{1}{2},$ where $p(\theta ,\varphi ,t),\theta \in [0,\pi ),\;\varphi \in [0,2\pi ),\;t>0,$ is the real function, $c>0,D>0,k>0$ are constants (see [4]), $\delta (\cdot )$ is the Dirac delta-function, and

$$\begin{aligned} \Delta _{(\theta ,\varphi )}=\frac{1}{\sin \theta }\;\frac{\partial }{ \partial \theta }\left( \sin {\theta }\;\frac{\partial }{\partial \theta } \right) +\frac{1}{\sin ^{2}{\theta }}\;\frac{\partial ^{2}}{\partial \varphi ^{2}} \end{aligned}$$

(14)

is the Laplace–Beltrami operator on the sphere.

It is known (see, i.e., [34], p. 72) that the eigenvalue problem for Laplace operator on the sphere has the following exact solution

$$\begin{aligned} \Delta _{(\theta ,\varphi )}Y_{lm}(\theta ,\varphi )=-l(l+1)Y_{lm}(\theta ,\varphi ),\qquad l=0,1,2,\dots ,\quad m=-l,\dots ,l, \end{aligned}$$

(15)

where $\{Y_{lm}(\theta ,\varphi )\}$ is the system of spherical harmonics. Therefore, it is natural to seek a solution of the problem (11),( 12),(13) in the form of the series

$$\begin{aligned} p(\theta ,\varphi ,t)=\sum _{l=0}^{\infty }\sum _{m=-l}^{l}b_{lm}(t)\;Y_{lm}(\theta ,\varphi ), \end{aligned}$$

(16)

where

$$\begin{aligned} b_{lm}(t)=\int _{\mathbb {S}^{2}}p(\theta ,\varphi ,t)\;Y_{lm}^{*}(\theta ,\varphi )\sin \theta d\varphi d\theta . \end{aligned}$$

(17)

Let

$$\begin{aligned} E_{\alpha ,\beta }(z)=\sum _{k=0}^{\infty }\frac{z^{k}}{\Gamma (\alpha k+\beta )},\ \alpha>0,\beta >0,\ z\in \mathbb {C}. \end{aligned}$$

(18)

be the two-parametrical Mittag-Leffler function ([21]). The case $\beta =1$ corresponds to the usual Mittag-Leffler function $E_{\alpha }(z)=E_{\alpha ,1}(z)$.

Introduce

$$\begin{aligned} \Omega =\left[ \frac{2D}{c^{2}}\sqrt{\frac{c^{4}}{4D^{2}}-k^{2}c^{2}l(l+1)} \right] \mathbf {\ 1}_{\left\{ l\le A\right\} }+\mathrm {i}\left[ \frac{2D}{ c^{2}}\sqrt{k^{2}c^{2}l(l+1)-\frac{c^{4}}{4D^{2}}}\right] {\mathbf {1}} _{\left\{ l>A\right\} }, \nonumber \\ \end{aligned}$$

(19)

where

$$\begin{aligned} A=A(D,c,k)=\frac{\sqrt{D^{2}k^{2}+c^{2}}-Dk}{2Dk}, \end{aligned}$$

(20)

and ${\mathbf {1}}_{\left\{ \cdot \right\} }$ denotes the binary indicator function.

Let

$$\begin{aligned} z_{\pm }=z_{\pm }^{l,\alpha }(t)=-\frac{c^{2}}{D}\frac{t^{\alpha }}{2}(1\pm \Omega ) . \end{aligned}$$

(21)

The proofs of the following two results are given in Sect. 5.

Theorem 1

Consider the functions

$$\begin{aligned} B_{\alpha ,l }^{(1)}(t)= & {} \frac{1}{2}[E_{\alpha }(z_{-}^{l,\alpha }(t))+E_{\alpha }(z_{+}^{l,\alpha }(t))], \end{aligned}$$

(22)

$$\begin{aligned} B_{\alpha ,l }^{(2)}(t)= & {} \frac{1}{2\Omega }[E_{\alpha }(z_{-}^{l,\alpha }(t))-E_{\alpha }(z_{+}^{l,\alpha }(t))], \end{aligned}$$

(23)

where $E_{\alpha }(z)$ denotes the Mittag-Leffler function.

Then, for $\frac{1}{2} < \alpha \le 1$ we have

$$\begin{aligned} \left| B_{\alpha ,l }^{(1)}(t)\right|\le & {} K_{0}^{(1)}\left[ \exp \left[ -K_{1}^{(1)}t(k^{2}c^{2}l(l+1))^{\frac{1}{2\alpha }}\right] +\frac{1}{ 1+t^{\alpha }kc\sqrt{l(l+1)}}\right] , \end{aligned}$$

(24)

$$\begin{aligned} \left| B_{\alpha ,l }^{(2)}(t)\right|\le & {} K_{0}^{(2)}\sqrt{\frac{ c^{2}}{4D^{2}k^{2}l(l+1)}} \bigg [\exp \left[ -K_{1}^{(2)}t(k^{2}c^{2}l(l+1))^{\frac{1}{2\alpha }}\right] \nonumber \\&+\frac{1}{1+t^{\alpha }kc\sqrt{l(l+1)}}\bigg ], \end{aligned}$$

(25)

for some positive constants $K_{0}^{(1)},K_{1}^{(1)},K_{0}^{(2)},K_{1}^{(2)}$ , and for $0 < \alpha \le \frac{1}{2}$ we have

$$\begin{aligned} \left| B_{\alpha ,l }^{(1)}(t)\right|\le & {} K_{2}^{(1)}\frac{1}{ 1+t^{\alpha }kc\sqrt{l(l+1)}}, \end{aligned}$$

(26)

$$\begin{aligned} \left| B_{\alpha ,l }^{(2)}(t)\right|\le & {} K_{2}^{(2)}\frac{1}{ 1+t^{\alpha }kc\sqrt{l(l+1)}} \end{aligned}$$

(27)

for some positive constants $K_{2}^{(1)},K_{2}^{(2)}.$

Theorem 2

The fundamental solution $p(\theta ,\varphi ,t)$ of the fractional hyperbolic diffusion point-source initial-value problems (11), (12), (13) are given by pointwise convergent series

$$\begin{aligned} p(\theta ,\varphi ,t)=\frac{1}{4\pi } \sum _{l=0}^{\infty } (2l+1)P_l(\cos \theta ) F_{l,\alpha }(t) , \end{aligned}$$

(28)

where $0<\alpha \le 1$ and $F_{l,\alpha }(t)$ is defined as

$$\begin{aligned} F_{l,\alpha }(t) = B_{\alpha ,l }^{(1)}(t)+B_{\alpha ,l }^{(2)}(t) , \end{aligned}$$

(29)

where $B_{\alpha ,l }^{(1)}(t)$ and $B_{\alpha ,l }^{(2)}(t)$ are given by eq. (22) and eq. (23), respectively.

Remark 1

The case $\alpha =1$ was considered in [4]. It is well-known that $E_{1}(z)=\exp (z)$, and then

$$\begin{aligned} E_{1}(z_{\pm }^{l,1}(t))=\exp {\left( -\frac{c^{2}t}{2D}(1\pm \Omega )\right) }=\exp {\left( -\frac{c^{2}t}{2D}\right) }\exp {\left( \mp \frac{ c^{2}t}{2D}\Omega \right) } \end{aligned}$$

which gives

$$\begin{aligned} B_{1,l}^{(1)}(t)= & {} \exp {\left( -\frac{c^{2}t}{2D}\right) }\left[ \cosh { (tM_{l})}{\mathbf {1}}_{l\le A}+\cos {(tM_{l}^{(1)})}{\mathbf {1}}_{l>A}\right] , \\ B_{1,l}^{(2)}(t)= & {} \exp {\left( -\frac{c^{2}t}{2D}\right) }\frac{c^{2}}{ 2D}\left[ \frac{\sinh {(tM_{l})}}{M_{l}}{\mathbf {1}}_{l\le A}+\frac{\sin { (tM_{l}^{(1)})}}{M_{l}^{(1)}}{\mathbf {1}}_{l>A}\right] , \end{aligned}$$

where

$$\begin{aligned} M_{l}=\sqrt{\frac{c^{4}}{4D^{2}}-k^{2}c^{2}l(l+1)},\qquad M_{l}^{(1)}=\sqrt{ c^{2}l(l+1)k^{2}-\frac{c^{4}}{4D^{2}}}, \end{aligned}$$

and A is given by 20. Then we obtain the solution given in [4] , that is,

$$\begin{aligned} p(\theta ,\varphi ,t)= & {} \exp {\left( -\frac{c^{2}t}{2D}\right) } \sum _{l=0}^{\infty }\frac{(2l+1)}{4\pi }P_l(\cos \theta ) \bigg \{\bigg [\cosh (tM_{l}) \nonumber \\&+\frac{c^{2}}{2DM_{l}}\sinh (tM_{l})]\bigg ]{\mathbf {1}}_{\left\{ l\le A\right\} }+\bigg [\cos (tM_{l}^{(1)})+\frac{c^{2}}{2DM_{l}^{(1)}}\sin (tM_{l}^{(1)})]\bigg ]{\mathbf {1}}_{\left\{ l>A\right\} }\bigg \},\nonumber \\ \end{aligned}$$

(30)

Note that the formula (30) is a clarification of the formula (21) [4].

Remark 2

The case $\alpha = 1/2$ is another one for which the solution has a nice expression. It is known [21] that $E_{1/2}(z)$ can be written as

$$\begin{aligned} E_{1/2}(z) = {\text{ e }}^{z^2} {\text {erfc}}(-z) , \end{aligned}$$

but the complementary error function does not have properties that allow us to write the combinations of $E_{1/2}(z_\pm ^{l,1/2}(t))$ in the expressions for $B_{1/2,l}^{(1)}(t)$ and $B_{1/2,l}^{(2)}(t)$ in a simple form. However, there is an integral representation for $E_{1/2}(z)$ that gives a better result, that is [41]

$$\begin{aligned} E_{1/2}(z) = \frac{2}{\sqrt{\pi }}\int _0^\infty {\text{ e }}^{-\omega ^2 + 2z\omega }\, d\omega . \end{aligned}$$

Using this integral representation for $E_{1/2}(z)$, we obtain for $B_{1/2,l}^{(1)}(t)$ and $B_{1/2,l}^{(2)}(t)$ that

$$\begin{aligned} B_{1/2,l}^{(1)}(t)= & {} \frac{2}{\sqrt{\pi }} \int _0^\infty {\text{ e }}^{-\omega ^2 - \frac{c^2}{2D}\sqrt{t}\omega } \bigg [ \cosh \left( \sqrt{t}M_l \omega \right) {\mathbf {1}}_{\{l\le A\}} \\&+ \cos \left( \sqrt{t}M_l^{(1)}\omega \right) {\mathbf {1}}_{\{l>A\}} \bigg ] d\omega , \\ B_{1/2,l}^{(2)}(t)= & {} \frac{2}{\sqrt{\pi }} \int _0^\infty {\text{ e }}^{-\omega ^2 - \frac{c^2}{2D}\sqrt{t}\omega } \frac{c^2}{2D} \bigg [ \frac{\sinh \left( \sqrt{t}M_l \omega \right) }{M_l} {\mathbf {1}}_{\{l\le A\}} \\&+ \frac{\sin \left( \sqrt{t}M_l^{(1)}\omega \right) }{M_l^{(1)}} {\mathbf {1}}_{\{l>A\}} \bigg ] d\omega , \end{aligned}$$

and the expression for $p(\theta ,\varphi ,t)$ is

$$\begin{aligned} \begin{aligned}&p(\theta ,\varphi ,t)= \frac{2}{\sqrt{\pi }} \sum _{l=0}^{\infty }\frac{(2l+1)}{4\pi }P_l(\cos \theta ) \int _0^\infty {\text{ e }}^{-\omega ^2 - \frac{c^2}{2D}\sqrt{t}\omega } \\&\quad \bigg \{ \bigg [ \cosh \left( \sqrt{t}M_l \omega \right) + \frac{c^2}{2D M_l} \sinh \left( \sqrt{t}M_l \omega \right) \bigg ]{\mathbf {1}}_{\{l\le A\}} \\&\qquad + \bigg [ \cos \left( \sqrt{t}M_l^{(1)}\omega \right) + \frac{c^2}{2D M_l^{(1)}} \sin \left( \sqrt{t}M_l^{(1)}\omega \right) \bigg ] {\mathbf {1}}_{\{l>A\}} \bigg \} d\omega . \end{aligned} \end{aligned}$$

(31)

Remark 3

It should be noted that Goldstein [20] and Kac [23] proposed the probabilistic interpretation of the solution of the classical one-dimensional telegraph equation as special random walk governing by the Poisson process. This line was continued by many authors; for example Orsingher and Beghin [41] generalised the results for the fractional telegraph equation for some particular values of $\alpha$. In particular for $\alpha = 1/2$, they obtain a nice generalization of Goldstein–Kac formulae replacing in the telegraph stochastic process the time by the reflecting Brownian motion. The multidimensional random motions which govern the hyperbolic diffusion equations is a long-standing problem, which was posed by Mark Kac more than 50 years ago and become a subject of intensive discussions among researchers on whether or not the multidimensional random flights could be described by the telegraph equations similarly to the one-dimensional case. Some exhaustive answers to this question can be found in the papers [27, 28] (see also the references therein). The best of our knowledge the random flight interpretation of the fractional hyperbolic diffusion equation is an interesting open question as well as stochastic solution of the fractional hyperbolic diffusion equations on the compact manifolds such as sphere.

4 Solution for Randomly Forced Fractional Spherical Hyperbolic Diffusion

In this section we use the results of Sect. 3 to derive solutions of the fractional hyperbolic diffusion equations with random initial conditions (or random data).

Consider the following hyperbolic diffusion equation on the sphere

$$\begin{aligned}&\frac{1}{c^{2}}\frac{\partial ^{2\alpha }}{\partial t^{2\alpha }}u(\theta ,\varphi ,t)+\frac{1}{D}\frac{\partial ^{\alpha }}{\partial t^{\alpha }} u(\theta ,\varphi ,t)=k^{2}\Delta _{(\theta ,\varphi )}\;u(\theta ,\varphi ,t),0<\alpha \le 1, \nonumber \\&\quad \theta \in [0,\pi ),\;\varphi \in [0,2\pi ),\;t>0, \end{aligned}$$

(32)

where $\Delta _{(\theta ,\varphi )}$ is the Laplace–Beltrami operator on the sphere given by (14).

Now, the random initial conditions are determined by the Gaussian isotropic random field on the sphere

$$\begin{aligned}&u(\theta ,\varphi ,t)\big |_{t=0}=T(\theta ,\varphi )=\sum _{l=0}^{\infty }\sum _{m=-l}^{l}a_{lm}Y_{lm}(\theta ,\varphi ), \end{aligned}$$

(33)

$$\begin{aligned}&\quad \left. \frac{\partial u(\theta ,\varphi ,t)}{\partial t}\right| _{t=0}=0, \end{aligned}$$

(34)

if $\frac{1}{2}<\alpha \le 1,$ and

$$\begin{aligned} u(\theta ,\varphi ,t)\big |_{t=0}=T(\theta ,\varphi )=\sum _{l=0}^{\infty }\sum _{m=-l}^{l}a_{lm}Y_{lm}(\theta ,\varphi ), \end{aligned}$$

(35)

if $0<\alpha <\frac{1}{2},$where $Y_{l}^{m}(\theta ,\varphi )$ are the spherical harmonics, and $a_{lm},$$m=-l,\dots ,l,$$l\ge 0,$ are complex-valued independent Gaussian random variables satisfying (6) and (8).

The following theorem will be proven in Sect. 5.

Theorem 3

If the angular spectrum $\{C_{l},l=0,1,2\ldots \}$ of the isotropic Gaussian random fields $T(\theta ,\varphi )$ from the initial value problem (32)–(34) satisfies summation (34), then the random solution of $u(\theta ,\varphi ,t)$ of the initial value problem (32)–(34) is given by the convergent (for each $t>0)$ in the Hilbert space $L_{2}(\Omega \times \mathbb {S}^{2},\sin \theta d\theta d\varphi )$ random series

$$\begin{aligned} u(\theta ,\varphi ,t)=\sum _{l=0}^{\infty }\sum _{m=-l}^{l}\xi _{lm}^{(\alpha )}\ Y_{lm}(\theta ,\varphi ), \end{aligned}$$

(36)

where the stochastic processes

$$\begin{aligned} \xi _{lm}^{(\alpha )}= a_{lm}F_{l,\alpha }(t), \end{aligned}$$

where $F_{l,\alpha }(t)$ was defined in eq. (29).

Moreover, the random field (36) is isotropic and Gaussian on $\mathbb {S }^{2},$and its covariance function is given by

$$\begin{aligned} \mathbf {Cov}(u(\theta ,\varphi ,t),u(\theta ^{\prime },\varphi ^{\prime },t^{\prime }))=\frac{1}{4\pi }\sum _{l=0}^{\infty }(2l+1)C_{l}P_{l}(\cos \Theta )F_{l,\alpha }(t)F_{l,\alpha }(t^{\prime }). \end{aligned}$$

Remark 4

The case $\alpha =1$ is related to the random spherical hyperbolic diffusion equation, see [4] . In this case the isotropic spherical random field has the following spectral representation in the Hilbert space $L_{2}(\Omega \times \mathbb {S}^{2},\sin \theta d\theta d\varphi ):$

$$\begin{aligned} u(\theta ,\varphi ,t)=\exp \left( -\frac{c^{2}t}{2D}\right) \sum _{l=0}^{\infty }\sum _{m=-l}^{l}Y_{l0}(\theta ,\varphi )\eta _{lm}(t),\qquad t\ge 0, \end{aligned}$$

(37)

where

$$\begin{aligned} \eta _{lm}(t)= & {} a_{lm}[\cosh \left( tM_{l}\right) +\frac{c^{2}}{2DK_{l}} \;\sinh \left( tM_{l}\right) ]{\mathbf {1}}_{\left\{ l\le A\right\} }+\\&+[\cos \left( tM_{l}^{\prime }\right) +\frac{c^{2}}{2DK_{l}^{\prime }}\sin \left( tM_{l}^{\prime }\right) ]{\mathbf {1}}_{\left\{ l>A\right\} }, \end{aligned}$$

while its covariance function can be written in the form

$$\begin{aligned}&\mathbf {Cov}(u(\theta ,\varphi ,t),u(\theta ^{\prime },\varphi ^{\prime },t^{\prime }))=\exp \left( -\frac{c^{2}}{2D}(t+t^{\prime })\right) \nonumber \\&\quad \times \sum _{l=0}^{\infty }\sum _{m=-l}^{l}Y_{lm}(\theta ,\varphi )Y_{lm}^{*}(\theta ^{\prime },\varphi ^{\prime }){\mathbf {E}}\eta _{lm}(t)\eta _{lm}^{*}(t^{\prime })=\exp \left( -\frac{c^{2}}{2D} (t+t^{\prime })\right) \nonumber \\&\quad \times (4\pi )^{-1}\sum _{l=0}^{\infty }(2l+1)C_{l}P_{l}(\cos \Theta )[A_{l}(t)A_{l}(t^{\prime })+B_{l}(t)B_{l}(t^{\prime })], \end{aligned}$$

(38)

where Ais given by 20, and the series (38) converges for every fixed t and $t^{\prime },$ that is

$$\begin{aligned} \sum _{l=0}^{\infty }(2l+1)C_{l}P_{l}(\cos \Theta )[L_{l}(t)L_{l}(t^{\prime })+B_{l}(t)B_{l}(t^{\prime })]<\infty , \end{aligned}$$

(39)

where

$$\begin{aligned} L_{l}(t)= & {} \left[ \cosh \left( tM_{l}\right) +\frac{c^{2}}{2DK_{l}}\;\sinh \left( tM_{l}\right) \right] {\mathbf {1}}_{\left\{ l\le A\right\} }\\ B_{l}(t)= & {} \left[ \cos \left( tM_{l}^{\prime }\right) +\frac{c^{2}}{2DK_{l}^{\prime } }\sin \left( tM_{l}^{\prime }\right) \right] {\mathbf {1}}_{\left\{ l>A\right\} }. \end{aligned}$$

Note that the formula (68) is a clarification of the formula (25) [4].

Noting that $\left| P_{l}(\cos \Theta )\right| \le 1,$ only a finite number of terms $L_{l}$ is non-zero, and there is a constant K such that $\sup _{t\ge 0}|B(t)|<C,$ we obtain that condition (39) follows from (7). This condition on the angular spectrum $C_{l},l\ge 0,$ guarantees the convergence of the series (68) in the Hilbert space $L_{2}(\Omega \times {\mathbb {S}}^{2},\sin \theta d\theta d\varphi ).$

5 Proofs

5.1 Proof of the Theorem 1

Theorem 1 follows from the asymptotic expansion of Mittag-Leffler functions. It is known [21] that, for $0 \le \alpha \le 2$ and $\pi \alpha /2< \theta < \mathop {\hbox {min}}\{\pi ,\pi \alpha \}$, we have the following estimates:

$$\begin{aligned} |E_\alpha (z)| \le {\left\{ \begin{array}{ll} M_1 \exp {(\mathop {\hbox {Re}}z^{1/\alpha })} + {\displaystyle \frac{M_2}{1+|z|}} , &{} \quad \text {for}\; |z|>0, \; |\mathop {\hbox {Arg}} z| \le \theta , \qquad \qquad \text { (a)} \\ {\displaystyle \frac{M_2}{1+|z|}} , &{} \quad \text {for}\; |z|>0, \; \theta \le |\mathop {\hbox {Arg}} z| \le \pi , \qquad \text {(b)} \end{array}\right. } \end{aligned}$$

with $M_1$ and $M_2$ does not depending on z.

We need an expression for $E_\alpha (z_\pm ^{l,\alpha }(t))$ for $l \rightarrow \infty$ . For l sufficiently large such that $l > A$ we have

$$\begin{aligned} z_{\pm }^{l,\alpha }(t)=-\frac{c^{2}}{D}\frac{t^{\alpha }}{2}\left( 1\pm \mathrm {i}\vartheta \right) , \end{aligned}$$

where

$$\begin{aligned} \vartheta = \sqrt{\frac{4 D^2 k^2 l(l+1)}{c^2}-1} . \end{aligned}$$

Then

$$\begin{aligned} |z_\pm ^{l,\alpha }(t)| = \frac{c^2 t^\alpha }{2D}\sqrt{1+\vartheta ^2} = k c t^\alpha \sqrt{l(l+1)} , \end{aligned}$$

(40)

and

$$\begin{aligned} \mathop {\hbox {Arg}}z_\pm ^{l,\alpha }(t) = \mp (\pi - \arctan \vartheta ) . \end{aligned}$$

(41)

In order to use the above estimates, we need to analyse $|\mathop {\hbox {Arg}} z_\pm ^{l,\alpha }(t)|$. We have

$$\begin{aligned} |\mathop {\hbox {Arg}}z_\pm ^{l,\alpha }(t)| = |\pi - \arctan \vartheta | . \end{aligned}$$

In the limit $l \rightarrow \infty$ we have $\arctan \vartheta = \pi /2 - c/(2 k D l) + \mathcal {O}(l^{-2})$, and then

$$\begin{aligned} |\mathop {\hbox {Arg}}z_\pm ^{l,\alpha }(t)| = \frac{\pi }{2} + \frac{c}{2Dk}\frac{1}{l} + \mathcal {O}(l^{-2}) . \end{aligned}$$

Therefore, for $\frac{1}{2} < \alpha \le 1$ and for l such that

$$\begin{aligned} l > l_0 = \frac{c}{Dk\pi (2\alpha -1)} \end{aligned}$$

we have

$$\begin{aligned} |\mathop {\hbox {Arg}}z_\pm ^{l,\alpha }(t)| \le \theta , \quad \theta = \pi \alpha . \end{aligned}$$

Then, for $\frac{1}{2} < \alpha \le 1$ we use the estimate provided by (a).

On the other hand, for $0 < \alpha \le \frac{1}{2}$, we have

$$\begin{aligned} |\mathop {\hbox {Arg}}z_\pm ^{l,\alpha }(t)| \ge \theta , \quad 0< \theta < \frac{\pi }{2 } . \end{aligned}$$

Then, for $0 < \alpha \le \frac{1}{2}$ we use the estimate provided by (b).

For the case $\frac{1}{2} < \alpha \le 1$ we need the expressions for $|z_\pm ^{l,\alpha }(t)|$ and $\mathop {\hbox {Re}}(z_\pm ^{l,\alpha }(t))^{1/\alpha }$. The expression for $|z_\pm ^{l,\alpha }(t)|$ is already given in eq. 40 and for $\mathop {\hbox {Re}}(z_\pm ^{l,\alpha }(t))^{1/\alpha }$ we obtain for eq. 41 that

$$\begin{aligned} \mathop {\hbox {Re}}(z_\pm ^{l,\alpha }(t))^{1/\alpha } = t (k^2 c^2 l(l+1))^{1/2\alpha } \cos \frac{(\pi -\arctan \vartheta )}{\alpha } . \end{aligned}$$

Then, for $B_{\alpha ,l}^{(1)}(t)$ and $B_{\alpha ,l}^{(2)}(t)$ given in eqs. (22) and (23) we have

$$\begin{aligned} \begin{aligned} |B_{\alpha ,l}^{(1)}(t)|&= \frac{1}{2}|E_\alpha (z_+^{l,\alpha }(t)) + E_\alpha (z_-^{l,\alpha }(t))| \le \frac{1}{2}(|E_\alpha (z_+^{l,\alpha }(t))| + |E_\alpha (z_-^{l,\alpha }(t))|) \\&= \frac{1}{2}(M_1^+ \exp {\left[ t (k^2 c^2 l(l+1))^{1/2\alpha } \cos \frac{ (\pi -\arctan \vartheta )}{\alpha }\right] } + {\displaystyle \frac{M_2^+}{1+k c t^\alpha \sqrt{l(l+1)}}} \\&\quad +M_1^- \exp {\left[ t (k^2 c^2 l(l+1))^{1/2\alpha } \cos \frac{ (\pi -\arctan \vartheta )}{\alpha }\right] } + {\displaystyle \frac{M_2^-}{1+k c t^\alpha \sqrt{l(l+1)}}}) \\&\le K_0^{(1)}\left( \exp {[-K_1^{(1)} t (k^2 c^2 l(l+1))^{1/2\alpha }]} + { \displaystyle \frac{1}{1+k c t^\alpha \sqrt{l(l+1)}}}\right) \end{aligned} \end{aligned}$$

where we denoted $\cos \frac{(\pi -\arctan \vartheta )}{\alpha } = -K_1^{(1)}$ (with $K_1^{(1)} \ge 0$ and $M_{\{1,2\}}^\pm \le 2 K_0^{(1)}$. This gives eq. (24). The proof of eq. (25) for $B_{\alpha ,l}^{(2)}(t)$ is analogous since

$$\begin{aligned} |B_{\alpha ,l}^{(2)}(t)| = \frac{1}{2\Omega }|E_\alpha (z_-^{l,\alpha }(t)) - E_\alpha (z_+^{l,\alpha }(t))| \le \frac{1}{2\Omega }(|E_\alpha (z_-^{l, \alpha }(t))| + |E_\alpha (z_+^{l,\alpha }(t))|) . \end{aligned}$$

The case $0 < \alpha \le \frac{1}{2}$ uses the estimate (b), which contains the second term of the estimate (a), and therefore gives eq. (26) and eq. (27). $\Box$

5.2 Proof of the Theorem 2

Let us look for the solution of eq. (11) satisfying the initial conditions given by eqs. (12) and (13). If we write $p(\theta ,\varphi ,t)$ in the form of a series as in eq. (16) and use eq. ( 15), we obtain that $b_{lm}(t)$ has to satisfy

$$\begin{aligned} \frac{1}{c^{2}}\mathrm {D}_{t}^{2\alpha }b_{lm}(t)+\frac{1}{D}\mathrm {D} _{t}^{\alpha }b_{lm}(t)+k^{2}l(l+1)b_{lm}(t)=0, \end{aligned}$$

(42)

where $\mathrm {D}_{t}^{\alpha }$ denotes the Caputo fractional derivative defined in eq. (9). From eq. (17) and the initial conditions in eqs. (12) and (13), we see that $b_{lm}(t)$ has to satisfy the initial conditions

$$\begin{aligned} b_{l,m}(0)=Y_{lm}^{*}(0,0)=Y_{l0}(0,0)\delta _{l}^{m}=\sqrt{\frac{2l+1}{ 4\pi }}\delta _{l}^{m},\qquad b_{lm}^{\prime }(0)=0 \end{aligned}$$

(43)

for $\frac{1}{2}<\alpha \le 1$, and

$$\begin{aligned} b_{l,m}(0)=Y_{lm}^{*}(0,0)=Y_{l0}(0,0)\delta _{l}^{m}=\sqrt{\frac{2l+1}{ 4\pi }}\delta _{l}^{m} \end{aligned}$$

(44)

for $0<\alpha \le \frac{1}{2}$.

In order to solve eq. (42) we use the Laplace transform. It is known [24] that the Laplace transform $\mathcal {L}$ of a Caputo derivative is given by

$$\begin{aligned} \mathcal {L}[\mathrm {D}_{t}^{2\alpha } b_{lm}(t)](s) = s^{2\alpha } B_{lm}(s) - s^{2\alpha -1} b_{lm}(0) - s^{2\alpha -2}b^\prime _{lm}(0) \end{aligned}$$

(45)

for $\frac{1}{2}< \alpha \le 1$, and

$$\begin{aligned} \mathcal {L}[\mathrm {D}_{t}^{2\alpha } b_{lm}(t)](s) = s^{2\alpha } B_{lm}(s) - s^{2\alpha -1} b_{lm}(0) \end{aligned}$$

(46)

for $0< \alpha \le \frac{1}{2}$, where $B_{lm}(s) = \mathcal {L} [b_{lm}(t)](s)$. However, because of the initial conditions (43) and (44), the cases $0 < \alpha \le \frac{1}{2}$ and $\frac{1}{2} < \alpha \le 1$ can be written in the same form, as in eq. (46). Moreover

$$\begin{aligned} \mathcal {L}[\mathrm {D}_{t}^{\alpha } b_{lm}(t)](s) = s^{\alpha } B_{lm}(s) - s^{\alpha -1} b_{lm}(0) . \end{aligned}$$

(47)

Then, using the Laplace transform in eq. (42) and eq. (46) and eq. (47), and solving the resulting expression for $B_{lm}(s)$, we obtain

$$\begin{aligned} B_{lm}(s) = H(s) b_{lm}(0) , \end{aligned}$$

(48)

where

$$\begin{aligned} H(s) = \frac{s^{2\alpha -1} + c^2 D^{-1} s^{\alpha -1}}{s^{2\alpha } + c^2 D^{-1} s^{\alpha } + c^2 k^2 l(l+1)} . \end{aligned}$$

(49)

In order to calculate $\mathcal {L}^{-1}[H(s)](t)$, we will write H(s) in the form

$$\begin{aligned} H(s)=\frac{1}{s}-\frac{s^{-1}c^{2}k^{2}l(l+1)}{s^{2\alpha }+c^{2}D^{-1}s^{\alpha }+c^{2}k^{2}l(l+1)}, \end{aligned}$$

(50)

and for s such that

$$\begin{aligned} \bigg |\frac{c^{2}k^{2}l(l+1)}{s^{2\alpha }+c^{2}D^{-1}s^{\alpha }}\bigg |<1, \end{aligned}$$

(51)

we have

$$\begin{aligned} H(s)=\frac{1}{s}-c^{2}k^{2}l(l+1)\sum _{n=0}^{\infty }[-c^{2}k^{2}l(l+1)]^{n} \frac{s^{-\alpha (n+1)-1}}{(s^{\alpha }+c^{2}D^{-1})^{n+1}}. \end{aligned}$$

(52)

Let us consider the three-parameter Mittag-Leffler function $E_{a,b}^{c}(z)$ , defined as

$$\begin{aligned} E_{a,b}^{c}(z)=\sum _{n=0}^{\infty }\frac{(c)_{n}}{n!}\frac{z^{n}}{\Gamma (an+b)}, \end{aligned}$$

(53)

where $(c)_{n}=\Gamma (c+n)/\Gamma (c)$ is the Pochhammer symbol. It is known [24] that

$$\begin{aligned} \mathcal {L}[t^{b-1}E_{a,b}^{c}(-\mu t^{a})](s)=\frac{s^{ac-b}}{(s^{\alpha }+\mu )^{c}}. \end{aligned}$$

(54)

Therefore, if we identify $a=\alpha$, $b=2\alpha (n+1)+1$$c=n+1$, and $\mu =c^{2}D^{-1}$, we can express the inverse Laplace transform $\mathcal {L} ^{-1}[H(s)](t)$ in the form

$$\begin{aligned} \mathcal {L}^{-1}[H(s)](t)= & {} 1-c^{2}k^{2}l(l+1)t^{2\alpha }\sum _{n=0}^{\infty }[-c^{2}k^{2}l(l+1)t^{2\alpha }]^{n}\nonumber \\&E_{\alpha ,1+2\alpha (n+1)}^{n+1}(-c^{2}D^{-1}t^{\alpha }). \end{aligned}$$

(55)

We can write the above series in a simple form if we use [49]

$$\begin{aligned} \sum _{n=0}^\infty (-zw)^n E_{a,2an+b}^{n+1}(z+w) = \frac{z E_{a,b}(z) - w E_{a,b}(w)}{z-w} . \end{aligned}$$

(56)

Taking

$$\begin{aligned} zw = c^2 k^2 l(l+1) t^{2\alpha } , \qquad z+w = -c^2 D^{-1} t^\alpha , \end{aligned}$$

(57)

we obtain that

$$\begin{aligned} z = z_+^{l,\alpha }(t) , \qquad w = z_-^{l,\alpha }(t) , \end{aligned}$$

(58)

where

$$\begin{aligned} z_\pm ^{l,\alpha }(t) = -\frac{c^2 t^\alpha }{2D} (1\pm \Omega ) , \qquad \Omega = \sqrt{1-\frac{4k^2 D^2 l(l+1)}{c^2}} . \end{aligned}$$

(59)

Then we have

$$\begin{aligned} \begin{aligned} \mathcal {L}^{-1}[H(s)](t) = 1 - c^2 k^2 l(l+1)t^{2\alpha }&\bigg [ \frac{ (1+\Omega )}{2\Omega } E_{\alpha ,1+2\alpha }(z_+^{l,\alpha }(t)) \\&- \frac{(1-\Omega )}{2\Omega } E_{\alpha ,1+2\alpha }(z_-^{l,\alpha }(t))\bigg ] . \end{aligned} \end{aligned}$$

(60)

Fig. 1

Multiplication factor $F_{l,\alpha }(t)$ for the angular power spectra as a function of l for $t=0.04$, $t=0.08$, $t=0.12$ and $t=0.16$ for $\alpha = 1$ (continuous blue curve), $\alpha = 0.9$ (dashed magenta curve) and $\alpha = 0.8$ (dashed dotted gray curve)

Fig. 2

Multiplication factor $F_{l,\alpha }(t)$ for the angular power spectra for $\alpha = 1$ (top left), $\alpha = 0.95$ (top right), $\alpha = 0.9$ (bottom left) and $\alpha = 0.8$ (bottom right) and different values of l (horizontal axis) and t (vertical axis)

Fig. 3

(Top) Multiplication factor for the angular power spectra for $l=500$ (top left) and $l=1500$ (top right) as a function of t for $\alpha = 1$ (continuous blue curve), $\alpha = 0.95$ (dashed magenta curve) and $\alpha = 0.9$ (dashed dotted gray curve); (Bottom) Multiplication factor for the angular power spectra for $l=500$ (bottom left) and $l=1500$ (bottom right) and different values of t (horizontal axis) and $\alpha$ (vertical axis)

Fig. 4

(Top) Multiplication factor $F_{l,\alpha }(t)$ for the angular power spectra as a function of time t for $l=500$ and $\alpha = 0.7$ (top left) and $\alpha = 0.6$ (top right), and for values $D = 1$ (continuous blue curve) and $D=100$ (dashed magenta curve). (Bottom) Same as above, but with plot ranges restricted to the intervals [0, 0.1] (bottom left) and [0, 0.01] (bottom right). (Color figure online)

Moreover, we know [21] that

$$\begin{aligned} E_{a,b}(z) = z E_{a,a+b}(z) + \frac{1}{\Gamma (b)} , \end{aligned}$$

(61)

from which we can write

$$\begin{aligned} E_{a,1+2a}(z) = \frac{1}{z^2}\left( E_{a,1}(z) - \frac{z}{\Gamma (1+a)} - 1\right) . \end{aligned}$$

(62)

Using this in eq. (60), we obtain, after some simplifications, that

$$\begin{aligned} \mathcal {L}^{-1}[H(s)](t) = \frac{(1+\Omega )}{2\Omega } E_\alpha (z_-^{l, \alpha }(t)) - \frac{(1-\Omega )}{2\Omega } E_\alpha (z_+^{l,\alpha }(t)) , \end{aligned}$$

(63)

where $E_\alpha (z) = E_{\alpha ,1}(z)$ is the usual Mittag-Leffler function. Using the functions $B_{\alpha ,l}^{(1)}$ and $B_{\alpha ,l}^{(2)}$ defined as in eq. (22) and eq. (23), that is,

$$\begin{aligned} B_{\alpha ,l }^{(1)}(t)= & {} \frac{1}{2}[E_{\alpha }(z_{-}^{l,\alpha }(t))+E_{\alpha }(z_{+}^{l,\alpha }(t))], \end{aligned}$$

(64)

$$\begin{aligned} B_{\alpha ,l }^{(2)}(t)= & {} \frac{1}{2\Omega }[E_{\alpha }(z_{-}^{l,\alpha }(t))-E_{\alpha }(z_{+}^{l,\alpha }(t))], \end{aligned}$$

(65)

we obtain that

$$\begin{aligned} b_{lm}(t) = b_{lm}(0) \mathcal {L}^{-1}[H(s)](t) = \sqrt{\frac{2l+1}{4\pi }} \delta _{lm} [B_{\alpha ,l }^{(1)}(t)+ B_{\alpha ,l}^{(2)}(t)] , \end{aligned}$$

(66)

where we have used the initial conditions given eq. (43) and eq. (44). Finally, using this expression for $b_{lm}(t)$ in eq. (16), the definition eq.(29), and the expression for $Y_{l0}(\theta ,\phi )$, we get eq. (28). $\square$

5.3 Proof of the Theorem 3

Let the two functions $f_{1}(\cdot )$ and $f_{2}(\cdot )$ on the sphere $\mathbb {S}^{2}$ belong to the space $L_{2}(\mathbb {S}^{2},\sin \theta d\theta d\varphi )$ and have the Fourier-Laplace coefficients

$$\begin{aligned} a_{lm}^{(i)}=\int _{\mathbb {S}^{2}}f_{i}(\theta ,\varphi )Y_{lm}^{*}(\theta ,\varphi )\sin \theta d\theta d\varphi ,\qquad i=1,2. \end{aligned}$$

Recall (see, i.e., [16]) that their non-commutative spherical convolution is defined as the Laplace series

$$\begin{aligned}{}[\ f_{1}*f_{2}](\theta ,\varphi )=\sum _{l=0}^{\infty }\sum _{m=-l}^{l}a_{lm}^{(*)}\ Y_{lm}(\theta ,\varphi ) \end{aligned}$$

(67)

with the Fourier-Laplace coefficients given by

$$\begin{aligned} a_{lm}^{(*)}=\sqrt{\frac{4\pi }{2l+1}}a_{lm}^{(1)}a_{l0}^{(2)}, \end{aligned}$$

provided that the series (67) converges in the corresponding Hilbert space.

Thus, the random solution $u(\theta ,\varphi ,t)$ of equation (32) with the initial values determined by (33) and (34) can be written as a spherical random field with the following Laplace series representation

$$\begin{aligned} u(\theta ,\varphi ,t)=[\ T*\ p_{t}](\theta ,\varphi )=\sum _{l=0}^{\infty }\sum _{m=-l}^{l}a_{lm}^{(t)}Y_{lm}(\theta ,\varphi ), \end{aligned}$$

(68)

provided that this series is convergent in the Hilbert space $L_{2}(\Omega \times \mathbb {S}^{2},\sin \theta d\theta d\varphi ),$ where $p_{t}=p(\theta ,\varphi ,t)$ is given by Theorem 1, and T is given by (33). The complex Gaussian random variables $a_{lm}^{(t)}$ are given by

$$\begin{aligned} a_{lm}^{(t)}=\sqrt{\frac{4\pi }{2l+1}}a_{lm}a_{l0}^{(p_{t})}, \end{aligned}$$

where $a_{l0}^{(p_{t})}=Y_{l0}^{*}(0,0)F_{l,\alpha }$ (t). This gives the first statement of the theorem.

By using the addition formula for spherical harmonics (see, i.e., [34], p. 66) we obtain the expression for covariance structure.

Noting that $\left| P_{l}(\cos \Theta )\right| \le 1,$ and the condition on the angular spectrum $C_{l},l\ge 0,$ guarantees the convergence of the series (68) in the Hilbert space $L_{2}(\Omega \times {\mathbb {S}}^{2},\sin \theta d\theta d\varphi ).$

Fig. 5

(Top) Multiplication factor $F_{l,\alpha }(t)$ for the angular power spectra as a function of l for $t=0.04$ (top left) and $t=0.08$ (top right) for $\alpha = 0.5$ (continuous blue curve), $\alpha = 0.4$ (dashed magenta curve) and $\alpha = 0.3$ (dashed dotted gray curve). (Bottom) Multiplication factor for the angular power spectra for $l=100$ (bottom left) and $l=1500$ (bottom right) as a function of t for $\alpha = 0.5$ (continuous blue curve), $\alpha = 0.4$ (dashed magenta curve) and $\alpha = 0.3$ (dashed dotted gray curve). (Color figure online)

6 Numerical Illustrations

From eq. (36) we see that the angular spectrum evolution over time is determined by the function $F_{l,\alpha }(t)$, defined by eq. (29 ). Let us call $F_{l,\alpha }(t)$ the multiplication factor. In this section we present some graphs showing the behaviour of this multiplication factor for some different values of its arguments. All plots were made using Mathematica 12.

In order to compare our results with the ones of [4], let us consider the values $c=1$, $D = 1$ and $k = 0.01$. Then the quantity A defined in eq. (20) equals $A = (\sqrt{10001}-1)/2$. In [4] the authors have used the time variable $t^\prime = c^2 t/2D$, but because of the presence of the parameter $\alpha$ in $t^\alpha$ in our expressions, we prefer to use the variable t, so that, for $\alpha = 1$, we have that $t^\prime = t/2$. Moreover, like in [4], we consider $l \le 2500$.

In Fig. 1 we have the plots of the multiplication factor $F_{l,\alpha }(t)$ as a function of l for four different times and for three different values of $\alpha$. The case $\alpha = 1$ corresponds to the case studied in [4]. In Fig. 2, instead of fixing the time variable at some values as in Fig. 1, we consider the contour plots of the multiplication factor $F_{l,\alpha }(t)$ as a function of l and of the time variable t in the interval [0, 1], for four different values of $\alpha$.

In Fig. 3 we have fixed values of l ($l=500$ and $l=1500$) and on the top we have the multiplication factor $F_{l,\alpha }(t)$ as a function of t for three different values of $\alpha$, while on the bottom we have $F_{l,\alpha }(t)$ as a function of t and with a continuous variation of $\alpha$.

It is easy to see from the case $\alpha =1$ that for $l>A$ we have a damped oscillatory behaviour with frequency increasing with the l number. For values $\alpha <1$, we see that we have an increasing attenuation of the amplitude with the decreasing of $\alpha$. The plots in Figs. 2 and  3 suggest that this oscillatory behaviour will turn into a purely diffusive behaviour after some value of $\alpha$, which we expect to be $\alpha = 1/2$. In Fig. 4 we have the plots of the multiplication factor $F_{l,\alpha }(t)$ as a function of time t for the values of $l = 500$ and $\alpha = 0.7$ and $\alpha = 0.6$. The plots displayed in continuous blue curves correspond to the same values used in the above plots, that is, $c=1$, $D = 1$ and $k = 0.01$. A quick visual inspection may suggest that we no longer have an oscillatory behaviour in such cases; however, this false impression is due to the values used for c, D and k. If we increase the value of D, which means that we are decreasing the contribution of the attenuation term $D^{-1} \partial ^\alpha u/\partial t^\alpha$ in eq. (32), we expect that the oscillatory behaviour will be easier to be seen. The plots displayed in dashed magenta curves correspond to the values $c=1$, $D = 100$ and $k=0.01$. Versions with reduced intervals for the plot ranges are shown in the bottom part of the figure in order to clearly visible the small but still present oscillatory behaviour.

The behaviour of the multiplication factor $F_{l,\alpha }(t)$ is expected to transform to a diffusion behaviour for $\alpha \le 1/2$, where we have a fractional deformation of the usual diffusion equation. This can be seen in Fig. 5, particularly in the bottom plots, where we can note the typical anomalous diffusion behaviour of models based on Mittag-Leffler functions [6, 21, 31].