Subordination Approach to Space-Time Fractional Diffusion

Abstract

The fundamental solution to the multi-dimensional space-time fractional diffusion equation is studied by applying the subordination principle, which provides a relation to the classical Gaussian function. Integral representations in terms of Mittag-Leffler functions are derived for the fundamental solution and the subordination kernel. The obtained integral representations are used for numerical evaluation of the fundamental solution for different values of the parameters.

1. Introduction

This work is concerned with the n-dimensional space-time fractional diffusion equation

$${\mathbb{D}}_t^{\beta }u(\mathbf{x},t)=-{(-\mathrm{\Delta })}^{\alpha }u(\mathbf{x},t),t>0,\mathbf{x}\in {\mathbb{R}}^n;u(\mathbf{x},0)=v(\mathbf{x});$$

(1)

where $0<\alpha ,\beta \le 1$, ${\mathbb{D}}_t^{\beta }$ is the Caputo time-fractional derivative [1,2]

$${\mathbb{D}}_t^{\beta }f(t)=\frac{1}{\mathsf{\Gamma }(1-\beta )}\frac{d}{dt}{\int }_0^t\frac{f(t)-f(0)}{{(t-\tau )}^{\beta }}d\tau ,t>0,0<\beta <1,$$

(2)

and $-{(-\mathrm{\Delta })}^{\alpha }$, $\alpha \in (0,1)$, is the full-space fractional Laplace operator in ${\mathbb{R}}^n$. Ten equivalent definitions of the fractional Laplacian $-{(-\mathrm{\Delta })}^{\alpha }$ are given in the survey paper [3]. In particular, it can be defined as a pseudo-differential operator, as follows:

$$\mathcal{F}\{-{(-\mathrm{\Delta })}^{\alpha }f;\kappa \}=-{\left|\kappa \right|}^{2\alpha }\mathcal{F}\{f;\kappa \},\kappa \in {\mathbb{R}}^n,$$

where $\mathcal{F}\{f;\kappa \}$ denotes the Fourier transform of a function f at the point $\kappa$. In the one-dimensional case $-{(-\mathrm{\Delta })}^{\alpha }$ is the Riesz space-fractional derivative of order $2\alpha$.

The space-time fractional diffusion Equation (1) has been extensively studied [4,5,6,7,8,9,10,11,12,13,14,15]. The solution $u(\mathbf{x},t)$ of Problem (1) is given in terms of the fundamental solution $G_{\alpha ,\beta ,n}(\mathbf{x},t)$ and the initial function $v(\mathbf{x})$, as follows:

$$u(\mathbf{x},t)={\int }_{{\mathbb{R}}^n}{G}_{\alpha ,\beta ,n}(\mathbf{y},t)v(\mathbf{x}-\mathbf{y})d\mathbf{y},\mathbf{x}\in {\mathbb{R}}^n,t>0.$$

Therefore, the behavior of the solution to Problem (1) is determined by the properties of the fundamental solution. In this paper, we limit our attention to representations of the fundamental solution $G_{\alpha ,\beta ,n}(\mathbf{x},t)$.

In the classical case $\alpha =\beta =1$, Equation (1) reduces to the standard diffusion equation with the fundamental solution $G_{1,1,n}(\mathbf{x},t)$, given by the Gaussian function (see e.g., [16]):

$$G_{1,1,n}(\mathbf{x},t)=\frac{1}{{(4\pi t)}^{n/2}}{e}^{-{\left|\mathbf{x}\right|}^2/4t},\mathbf{x}\in {\mathbb{R}}^n,t>0.$$

(3)

In the fractional-order setting, the following closed-form representations for the fundamental solution are known:

$$\begin{array}{ccc}\hfill G_{\frac{\alpha }{2},\alpha ,1}(x,t)& =& \frac{1}{\pi }\frac{t^{\alpha }{\left|x\right|}^{\alpha -1}sin(\alpha \pi /2)}{t^{2\alpha }+2t^{\alpha }{\left|x\right|}^{\alpha }cos(\alpha \pi /2)+{\left|x\right|}^{2\alpha }},x\in \mathbb{R},t>0,0<\alpha \le 1;\hfill \end{array}$$

(4)

$$\begin{array}{ccc}\hfill G_{\alpha ,\alpha ,2}(\mathbf{x},t)& =& \frac{1}{4\pi t}{\left(\frac{{\left|\mathbf{x}\right|}^2}{4t}\right)}^{\alpha -1}{E}_{\alpha ,\alpha }\left(-{\left(\frac{{\left|\mathbf{x}\right|}^2}{4t}\right)}^{\alpha }\right),\mathbf{x}\in {\mathbb{R}}^2,t>0,0<\alpha \le 1;\hfill \end{array}$$

(5)

$$\begin{array}{ccc}\hfill G_{\frac{1}{2},\frac{1}{2},1}(x,t)& =& \frac{1}{2{\pi }^{3/2}\sqrt{t}}{e}^{x^2/4t}{\mathcal{E}}_1\left(\frac{x^2}{4t}\right),x\in \mathbb{R},t>0;\hfill \end{array}$$

(6)

$$\begin{array}{ccc}\hfill G_{\frac{1}{2},\frac{1}{2},n}(\mathbf{x},t)& =& \frac{\mathsf{\Gamma }\left(\frac{n+1}{2}\right)}{2^n{\pi }^{n/2+1}{t}^{n/2}}U\left(\frac{n+1}{2},\frac{n+1}{2},\frac{{\left|\mathbf{x}\right|}^2}{4t}\right),\mathbf{x}\in {\mathbb{R}}^n,t>0;\hfill \end{array}$$

(7)

where $E_{\alpha ,\alpha }$ denotes the Mittag-Leffler function (see (18)), ${\mathcal{E}}_1$ is the exponential integral [17]

$${\mathcal{E}}_1(r)={\int }_r^{\infty }\frac{e^{-\xi }}{\xi }d\xi ,$$

(8)

and U is the Tricomi’s confluent hypergeometric function [17]

$$U(a,b,r)=\frac{1}{\mathsf{\Gamma }(a)}{\int }_0^{\infty }{\xi }^{a-1}{(1+\xi )}^{b-a-1}{e}^{-r\xi }d\xi ,a>0,r>0.$$

(9)

Representation (4) can be found in a more general setting in [5]. Formula (5) was first derived in the paper [10]. Formula (6) is established in the early work [4]. A derivation of representations (4)–(7) using the subordination relation (see (10)) can be found in [15]. In [13,15], additional closed-form representations for the fundamental solution were derived from (4)–(6) by applying the relations between $G_{\alpha ,\beta ,n+2}(\mathbf{x},t)$ and $G_{\alpha ,\beta ,n}(\mathbf{x},t)$. However, all such simple closed-form expressions for the fundamental solution in terms of known special functions are limited to particular values of the parameters.

Extensive research has been devoted to representations of the fundamental solution in the form of the Mellin-Barnes integral or the Fox H-function, such as in [5,6,11] for the one-dimensional and [12,13,14] for the multi-dimensional space-time fractional diffusion-wave equation. One of the advantages of such representations is that the asymptotic behavior of the fundamental solution can be derived from them, because the asymptotic behavior of the Fox H-function has been well-studied (see e.g., [18] or [19]).

An alternative approach to dealing with the space-time fractional diffusion Equation (1) is based on the subordination formula, which relates the fundamental solution $G_{\alpha ,\beta ,n}(\mathbf{x},t)$ and the Gaussian function $G_{1,1,n}(\mathbf{x},t)$ as follows [14,15]

$$G_{\alpha ,\beta ,n}(\mathbf{x},t)={\int }_0^{\infty }{\psi }_{\alpha ,\beta }(t,\tau )G_{1,1,n}(\mathbf{x},\tau )d\tau ,\mathbf{x}\in {\mathbb{R}}^n,t>0,$$

(10)

where ${\psi }_{\alpha ,\beta }(t,\tau )$ is a unilateral probability density function (pdf) in $\tau$, that is:

$${\psi }_{\alpha ,\beta }(t,\tau )\ge 0,{\int }_0^{\infty }{\psi }_{\alpha ,\beta }(t,\tau )d\tau =1.$$

(11)

The subordination kernel ${\psi }_{\alpha ,\beta }(t,\tau )$ depends on the similarity variable $\tau t^{-\beta /\alpha }$ and admits the representation [14]

$${\psi }_{\alpha ,\beta }(t,\tau )=t^{-\beta /\alpha }{K}_{\alpha ,\beta }(\tau t^{-\beta /\alpha }),$$

(12)

where the function $K_{\alpha ,\beta }(r)$ can be defined as the inverse Laplace transform of the Mittag-Leffler function $E_{\beta }(-{\lambda }^{\alpha })$, that is:

$${\int }_0^{\infty }{e}^{-\lambda r}{K}_{\alpha ,\beta }(r)dr=E_{\beta }(-{\lambda }^{\alpha }).$$

(13)

The Laplace transform pair (13) was first derived in [14] (see Remark 4.4).

It is worth noting that some known basic properties of $G_{\alpha ,\beta ,n}(\mathbf{x},t)$ follow in a straightforward way from the subordination relation (10), taking into account that the subordination kernel is a pdf. In this way, we can prove that for any dimension $n\ge 1$, the fundamental solution $G_{\alpha ,\beta ,n}(\mathbf{x},t)$ is a spatial pdf evolving in time:

$$G_{\alpha ,\beta ,n}(\mathbf{x},t)\ge 0,{\int }_{{\mathbb{R}}^n}{G}_{\alpha ,\beta ,n}(\mathbf{x},t)d\mathbf{x}=1.$$

Therefore, $G_{\alpha ,\beta ,n}(\mathbf{x},t)$, $0<\alpha ,\beta \le 1$, inherits this property of the classical Gaussian kernel $G_{1,1,n}(\mathbf{x},t)$. In a similar way, estimates for the fundamental solution $G_{\alpha ,\beta ,n}(\mathbf{x},t)$ can be derived from known estimates for the Gaussian kernel $G_{1,1,n}(\mathbf{x},t)$. For example, since $\parallel G_{1,1,n}{(·,t)\parallel }_{L^1({\mathbb{R}}^n)}=1$ (see e.g., [16], Remark 3.7.10.), the subordination Formula (10), together with properties (11) imply

$$\parallel G_{\alpha ,\beta ,n}{(·,t)\parallel }_{L^1({\mathbb{R}}^n)}\le {\int }_0^{\infty }{\psi }_{\alpha ,\beta }(t,\tau ){\parallel G_{1,1,n}(·,\tau )\parallel }_{L^1({\mathbb{R}}^n)}d\tau \le {\int }_0^{\infty }{\psi }_{\alpha ,\beta }(t,\tau )d\tau =1.$$

A principle of subordination is closely related to the concept of subordination in stochastic processes [20,21]. It has been extensively studied and employed in the context of fractional order equations. The subordination principle for space-fractional evolution equations has been established in [22] in the setting of abstract Cauchy problems. Subordination formulae for the one-dimensional space-time fractional diffusion equation have been studied in [6,9]. In [14,15], subordination principles for the multi-dimensional space-time fractional diffusion equation are deduced. In the case of time-fractional evolution equations with general time-fractional operators, subordination principles have been studied and employed in [23,24,25,26,27,28].

Based on the subordination principles for space- and time-fractional diffusion equations and the dominated convergence theorem, exact asymptotic expressions for the fundamental solution of the multi-dimensional space-time fractional diffusion equation and more general nonlocal equations have recently been established in [29]. For completeness, we next present the asymptotic expansions for $G_{\alpha ,\beta ,n}(\mathbf{x},t)$, $\alpha ,\beta \in (0,1)$, from [29], Corollary 2.6 (written in our notations and in a slightly more compact form):

If ${\left|\mathbf{x}\right|}^{-2\alpha }{t}^{\beta }\to \infty$, then

$$G_{\alpha ,\beta ,n}(\mathbf{x},t)\sim \left\{\begin{array}{c}{\displaystyle \frac{1}{2\alpha sin\left(\frac{\pi }{2\alpha }\right)\mathsf{\Gamma }(1-\frac{\beta }{2\alpha })}{t}^{-\frac{\beta }{2\alpha }},n=1,\alpha \in (1/2,1),}\hfill \\ {\displaystyle \frac{\beta }{\pi \mathsf{\Gamma }(1-\beta )}{t}^{-\beta }ln\left({t\left|x\right|}^{-1/\beta }\right),n=1,\alpha =1/2,}\hfill \\ {\displaystyle \frac{1}{4^{\alpha }{\pi }^{n/2}}\frac{\mathsf{\Gamma }(n/2-\alpha )}{\mathsf{\Gamma }(\alpha )\mathsf{\Gamma }(1-\beta )}\frac{t^{-\beta }}{{\left|\mathbf{x}\right|}^{n-2\alpha }},n>2\alpha .}\hfill \end{array}\right.$$

(14)

If ${\left|\mathbf{x}\right|}^{-2\alpha }{t}^{\beta }\to 0$, then

$$G_{\alpha ,\beta ,n}(\mathbf{x},t)\sim \frac{4^{\alpha }}{{\pi }^{n/2}}\frac{\alpha \mathsf{\Gamma }(n/2+\alpha )}{\mathsf{\Gamma }(1-\alpha )\mathsf{\Gamma }(\beta +1)}\frac{t^{\beta }}{{\left|\mathbf{x}\right|}^{n+2\alpha }}.$$

(15)

The asymptotic expansions (14) and (15) are in agreement with those obtained for particular ranges of parameter values in, for example, [5,11,15], as well as with the asymptotic behavior of the closed-form solutions (4)–(7), which can be checked by taking into account the asymptotic expansions for the exponential integral ([17], Eqs. 5.1.11 and 5.1.51)

$${\mathcal{E}}_1(r)\sim ln\left(\frac{1}{r}\right),r\to 0;{\mathcal{E}}_1(r)\sim \frac{e^{-r}}{r}\left(1-\frac{1}{r}\right),r\to +\infty ,$$

(16)

and for the Tricomi’s confluent hypergeometric function ([17], Section 13.5)

$$U(a,b,r)\sim \frac{\mathsf{\Gamma }(b-1)}{\mathsf{\Gamma }(a)}{r}^{1-b},r\to 0,b>1;U(a,b,r)\sim r^{-a},r\to +\infty ,$$

(17)

and using some basic properties of the Gamma function.

In the present work, the subordination Formula (10) serves as a starting point for deriving integral representations for the fundamental solution $G_{\alpha ,\beta ,n}(\mathbf{x},t)$. First, an integral representation for the subordination kernel ${\psi }_{\alpha ,\beta }(t,\tau )$ is established in terms of the Mittag-Leffler function of complex argument. Let us note that a study of the function ${\psi }_{\alpha ,\beta }(t,\tau )$ is of interest, since it also plays the role of subordination kernel related to problems with more general spatial operators, such as in [15]. In addition, ${\psi }_{\alpha ,\beta }(t,x)$ coincides with the solution of the one-dimensional space-time fractional diffusion equation with the Riesz-Feller space-fractional derivative of order $\alpha$ and skewness $-\alpha$, as well as the Caputo time-derivative of order $\beta$, studied in [5] (see [15], Remark 3). Next, based on the subordination Formula (10), we derive integral representations for the fundamental solution in terms of Mittag-Leffler functions, which are appropriate for numerical implementation.

The paper is organized as follows. Definitions and basic properties of Mittag-Leffler functions and Bessel functions of the first kind are listed in the next section. In Section 3, an integral representation for the subordination kernel is established. In Section 4, computable integral representations for the fundamental solution are derived for $n=1,2,3$ and used for numerical experiments.

2. Preliminaries

The Mittag-Leffler function is an entire function defined by the series [1,2,30]

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

(18)

For $0<\alpha <2$ and $\beta \in \mathbb{R}$, the following asymptotic expansion for large $\left|z\right|$ holds true in the sector of the complex plane $|arg(-z)|<(1-\alpha /2)\pi$

$$E_{\alpha ,\beta }(z)=-\sum _{k=1}^{N-1}\frac{z^{-k}}{\mathsf{\Gamma }(\beta -\alpha k)}+{O(|z|}^{-N}),|z|\to \infty .$$

(19)

Therefore, taking into account the identity $\mathsf{\Gamma }{(z)}^{-1}=0$ for $z=0,-1,-2,\dots$ we derive from (19) two useful asymptotic expressions for $\left|z\right|\to \infty$ and $|argz|<(1-\alpha /2)\pi$

$$E_{\alpha }(-z)\sim \frac{z^{-1}}{\mathsf{\Gamma }(1-\alpha )};E_{\alpha ,\beta }(-z)\sim -\frac{z^{-2}}{\mathsf{\Gamma }(\beta -2\alpha )},\beta -\alpha =0,-1,-2,\dots$$

(20)

The faster decay for large $\left|z\right|$ of the second function in (20), compared to the first, will be used essentially in this work.

The relations

$$\frac{d}{dz}{E}_{\alpha }(-z^{\alpha })=-z^{\alpha -1}{E}_{\alpha ,\alpha }(-z^{\alpha }),\frac{d}{dz}\left(z^{\alpha -1}{E}_{\alpha ,\alpha }(-z^{\alpha })\right)=z^{\alpha -2}{E}_{\alpha ,\alpha -1}(-z^{\alpha }),$$

(21)

can be derived directly from the definition (18) of the Mittag-Leffler function; (21) and (20) imply that, by differentiation of the Mittag-Leffler function $E_{\alpha }(-z^{\alpha })$, a faster decay for large $\left|z\right|$ can be achieved.

We point out the following representation of the Mittag-Leffler functions as Laplace transforms (see [31]):

$$t^{\beta -1}{E}_{\alpha ,\beta }(-\mu t^{\alpha })=\frac{1}{\pi }{\int }_0^{\infty }{e}^{-rt}\frac{r^{\alpha }sin\beta \pi +\mu sin(\beta -\alpha )\pi }{r^{2\alpha }+2\mu r^{\alpha }cos\alpha \pi +{\mu }^2}{r}^{\alpha -\beta }dr,$$

(22)

where $\mu >0$, $0<\alpha ,\beta \le 1$, excluding the case $\alpha =\beta =1$. Expression (22) is appropriate for numerical computation of the Mittag-Leffler functions. Let us note, however, that (22) is valid only for real values of $\mu$. For computation of the Mittag-Leffler function of complex arguments, another technique should be used (see e.g., [32]).

The Bessel function of the first kind $J_{\nu }(z)$ is defined by the series [17]

$$J_{\nu }(z)=\sum _{k=0}^{\infty }\frac{{(-1)}^k{(z/2)}^{\nu +2k}}{k!\mathsf{\Gamma }(\nu +k+1)}.$$

(23)

The following particular expressions are of interest in the present work:

$$J_{-1/2}(z)=\sqrt{\frac{2}{\pi z}}cosz,J_{1/2}(z)=\sqrt{\frac{2}{\pi z}}sinz,J_0(z)=\frac{1}{\pi }{\int }_0^{\pi }cos(zcos\theta )d\theta .$$

(24)

The asymptotic expansions of the Bessel function $J_{\nu }(r)$ for small and large real arguments are as follows:

$$J_{\nu }(r)\sim \frac{1}{\mathsf{\Gamma }(\nu +1)}{\left(\frac{r}{2}\right)}^{\nu },r\to 0;J_{\nu }(r)\sim \sqrt{\frac{2}{\pi r}}cos(r-\nu \pi /2-\pi /4),r\to +\infty .$$

(25)

For more details on Mittag-Leffler and Bessel functions, we refer to [2,17,30,33,34].

3. An Integral Representation for the Subordination Kernel

Representations of the subordination kernel ${\psi }_{\alpha ,\beta }(t,\tau )$ are useful in view of the integral expression (10) for the fundamental solution. In a limited number of particular cases, the subordination kernel can be expressed in terms of elementary functions [15,16,22]:

$${\psi }_{\frac{1}{2},1}(t,\tau )=\frac{t{e}^{-t^2/4\tau }}{2\sqrt{\pi }{\tau }^{3/2}},{\psi }_{1,\frac{1}{2}}(t,\tau )=\frac{1}{\sqrt{\pi t}}{e}^{-{\tau }^2/4t},$$

(26)

$${\psi }_{\alpha ,\alpha }(t,\tau )=\frac{1}{\pi }\frac{t^{\alpha }{\tau }^{\alpha -1}sin\alpha \pi }{t^{2\alpha }+2t^{\alpha }{\tau }^{\alpha }cos\alpha \pi +{\tau }^{2\alpha }},0<\alpha <1.$$

(27)

However, for arbitrary values of the fractional parameters, explicit expressions are not available and other types of representations are needed.

The following Laplace transform pairs for the subordination kernel ${\psi }_{\alpha ,\beta }(t,\tau )$ can be derived from (12) and (13) (see also [15]):

$${\int }_0^{\infty }{\psi }_{\alpha ,\beta }(t,\tau )e^{-\lambda \tau }d\tau =E_{\beta }(-{\lambda }^{\alpha }{t}^{\beta }),$$

(28)

and

$${\int }_0^{\infty }{\psi }_{\alpha ,\beta }(t,\tau )e^{-st}dt=s^{\beta -1}{\tau }^{\alpha -1}{E}_{\alpha ,\alpha }(-s^{\beta }{\tau }^{\alpha }).$$

(29)

In this section, we deduce an integral representation of the subordination kernel ${\psi }_{\alpha ,\beta }(t,\tau )$ by inversion of the Laplace transform pair (29). We choose (29) instead of (28), because of the faster decay for large arguments of the correponding Mittag-Leffler function (see (20)).

Assume $0<\alpha ,\beta \le 1$ and $\alpha +\beta <2$. Applying the complex Laplace inversion formula to (29) yields:

$${\psi }_{\alpha ,\beta }(t,\tau )=\frac{{\tau }^{\alpha -1}}{2\pi i}{\int }_{c-i\infty }^{c+i\infty }{e}^{st}{s}^{\beta -1}{E}_{\alpha ,\alpha }(-{\tau }^{\alpha }{s}^{\beta })ds,c>0,$$

(30)

where $s^{\beta }=exp(\beta lns)$ means the principal branch of the corresponding multi-valued function defined in the whole complex plane cut along the negative real semi-axis. Since the Mittag-Leffler function is an entire function, $E_{\alpha ,\alpha }(-{\tau }^{\alpha }{s}^{\beta })$ is analytic for $s\in \mathbb{C}\backslash (-\infty ,0]$. Therefore, by the Cauchy’s theorem, the integral in (30) can be replaced by an integral over the composite contour $\mathsf{\Gamma }={\mathsf{\Gamma }}_1^{-}\cup {\mathsf{\Gamma }}_2^{-}\cup {\mathsf{\Gamma }}_3\cup {\mathsf{\Gamma }}_2^{+}\cup {\mathsf{\Gamma }}_1^{+}$, where

$${\mathsf{\Gamma }}_1^{\pm }=\{s=q\pm iR,q\in (0,c)\},{\mathsf{\Gamma }}_2^{\pm }=\{s=r{e}^{\pm i\pi /2},r\in (\rho ,R)\},{\mathsf{\Gamma }}_3=\{s=\rho e^{i\theta },\theta \in (-\pi /2,\pi /2)\},$$

with an appropriate orientation (see Figure 1) and letting $\rho \to 0$, $R\to \infty$.

Since ${(q+iR)}^{\beta }\sim R^{\beta }{e}^{i\beta \pi /2}$ as $R\to \infty$, for the integration on ${\mathsf{\Gamma }}_1^{+}$ as $R\to \infty$ we obtain

$$\left|{\int }_{{\mathsf{\Gamma }}_1^{+}}{e}^{st}{s}^{\beta -1}{E}_{\alpha ,\alpha }(-{\tau }^{\alpha }{s}^{\beta })ds\right|\le C{\int }_0^c{e}^{qt}{R}^{\beta -1}|E_{\alpha ,\alpha }(-{\tau }^{\alpha }{R}^{\beta }{e}^{i\beta \pi /2})|dq\to 0,R\to \infty ,$$

(31)

due to the asymptotic expansion (20) for the Mittag-Leffler function, which is satisfied since $|arg({\tau }^{\alpha }{R}^{\beta }{e}^{i\beta \pi /2})|=\beta \pi /2<(1-\alpha /2)\pi$. The integral on ${\mathsf{\Gamma }}_1^{-}$ is treated in the same way.

Concerning the integral over ${\mathsf{\Gamma }}_3$, we have

$$\left|{\int }_{{\mathsf{\Gamma }}_3}{e}^{st}{s}^{\beta -1}{E}_{\alpha ,\alpha }(-{\tau }^{\alpha }{s}^{\beta })ds\right|\le {\int }_{-\pi /2}^{\pi /2}{e}^{\rho tcos\theta }{\epsilon }^{\beta }\left|E_{\alpha ,\alpha }(-{\tau }^{\alpha }{\rho }^{\beta }{e}^{i\beta \theta })\right|d\theta \to 0,\rho \to 0,$$

(32)

since the Mittag-Leffler function under the integral sign is bounded as $\rho \to 0$. Therefore, (30)–(32) imply that ${\psi }_{\alpha ,\beta }(t,\tau )$ is given by the integral over ${\mathsf{\Gamma }}_2^{+}\cup {\mathsf{\Gamma }}_2^{-}$ along the imaginary axis with $\rho \to 0$ and $R\to \infty$. This implies:

$$\begin{array}{c}{\psi }_{\alpha ,\beta }(t,\tau )=\frac{{\tau }^{\alpha -1}}{2\pi i}{\int }_{-i\infty }^{i\infty }{e}^{st}{s}^{\beta -1}{E}_{\alpha ,\alpha }(-{\tau }^{\alpha }{s}^{\beta })ds\hfill \\ =\frac{{\tau }^{\alpha -1}}{2\pi i}\left({\int }_0^{\infty }exp(rt{e}^{i\pi /2})r^{\beta -1}{e}^{i\beta \pi /2}{E}_{\alpha ,\alpha }(-{\tau }^{\alpha }{r}^{\beta }{e}^{i\beta \pi /2})dr\right.\hfill \\ +\left.{\int }_0^{\infty }exp(rt{e}^{-i\pi /2})r^{\beta -1}{e}^{-i\beta \pi /2}{E}_{\alpha ,\alpha }(-{\tau }^{\alpha }{r}^{\beta }{e}^{-i\beta \pi /2}dr\right).\hfill \end{array}$$

Therefore,

$${\psi }_{\alpha ,\beta }(t,\tau )=\frac{{\tau }^{\alpha -1}}{\pi }{\int }_0^{\infty }{r}^{\beta -1}\Im \left\{e^{i(rt+\beta \pi /2)}{E}_{\alpha ,\alpha }(-{\tau }^{\alpha }{r}^{\beta }{e}^{i\beta \pi /2})\right\}dr.$$

(33)

We observe that the integral in (33) is convergent, since the integrand behaves as $r^{\beta -1}$ for $r\to 0$ and as $r^{-\beta -1}$ for $r\to \infty$ due to the asymptotic Expansion (20) for the Mittag-Leffler function. The representation (33) can also be rewriten in the form

$${\psi }_{\alpha ,\beta }(t,\tau )=\frac{{\tau }^{\alpha -1}}{\pi }{\int }_0^{\infty }{r}^{\beta -1}\left(cos(rt+\beta \pi /2)I_{\alpha ,\beta }(r,\tau )+sin(rt+\beta \pi /2)R_{\alpha ,\beta }(r,\tau )\right)dr,$$

(34)

where

$$I_{\alpha ,\beta }(r,\tau )=\Im \{E_{\alpha ,\alpha }(-{\tau }^{\alpha }{r}^{\beta }{e}^{i\beta \pi /2})\},R_{\alpha ,\beta }(r,\tau )=\Re \{E_{\alpha ,\alpha }(-{\tau }^{\alpha }{r}^{\beta }{e}^{i\beta \pi /2})\}.$$

For the numerical implementation of Formula (34), the real and imaginary parts above can be numerically calculated by employing a method of computation of the Mittag-Leffler function of complex arguments.

In the particular case of $\alpha =1$ (time-fractional diffusion), representation (34) yields the following simpler formula for the subordination kernel

$${\psi }_{1,\beta }(t,\tau )=\frac{1}{\pi }{\int }_0^{\infty }{r}^{\beta -1}sin\left(rt+\beta \pi /2-\tau r^{\beta }sin\beta \pi /2\right)exp(-\tau r^{\beta }cos\beta \pi /2)dr.$$

(35)

Let us recall the relation ${\psi }_{1,\beta }(t,\tau )=t^{-\beta }{M}_{\beta }(\tau t^{-\beta })$, where $M_{\beta }(·)$ denotes the Mainardi function (see [30]). Numerical results based on Formula (35) for the subordination kernel ${\psi }_{1,\beta }(t,\tau )$ are given in Figure 2.

4. Integral Representations for the Fundamental Solution

According to the subordination Relation (10) and the formula for the Gaussian kernel (3), the fundamental solution of Problem (1) admits the representation

$$G_{\alpha ,\beta ,n}(\mathbf{x},t)=\frac{1}{{(4\pi )}^{n/2}}{\int }_0^{\infty }{\psi }_{\alpha ,\beta }(t,\tau ){\tau }^{-n/2}{e}^{-{\left|\mathbf{x}\right|}^2/4\tau }d\tau ,\mathbf{x}\in {\mathbb{R}}^n,t>0.$$

(36)

Subordination Formula (36) yields after the change of variables $\sigma =1/\tau$

$$G_{\alpha ,\beta ,n}(\mathbf{x},t)=\frac{1}{{(4\pi )}^{n/2}}{\int }_0^{\infty }{\psi }_{\alpha ,\beta }(t,{\sigma }^{-1}){\sigma }^{n/2-2}{e}^{-a\sigma }d\sigma ,a={\left|\mathbf{x}\right|}^2/4.$$

(37)

Applying the formula for the Laplace transform ([35], Section 4.1, Eq. (25))

$${\int }_0^{\infty }{\sigma }^{\nu -1}f({\sigma }^{-1})e^{-a\sigma }d\sigma =a^{-\frac{1}{2}\nu }{\int }_0^{\infty }{s}^{\frac{1}{2}\nu }{J}_{\nu }(2\sqrt{as})g(s)ds,Re\nu >-1,$$

where $J_{\nu }(·)$ denotes the Bessel Function (23) and $g(s)=\mathcal{L}\{f;s\}={\int }_0^{\infty }{e}^{-s\sigma }f(\sigma )d\sigma$, we deduce from (37) and (28) the following representation:

$$G_{\alpha ,\beta ,n}(\mathbf{x},t)=\frac{{\left|\mathbf{x}\right|}^{1-\frac{n}{2}}}{{(2\pi )}^{\frac{n}{2}}}{\int }_0^{\infty }{\sigma }^{\frac{n}{2}}{J}_{\frac{n}{2}-1}(|\mathbf{x}|\sigma )E_{\beta }(-{\sigma }^{2\alpha }{t}^{\beta })d\sigma .$$

(38)

The obtained integral representation (38) is not new—see, for example, [12,13,14], where it is deduced by applying a different argument.

Let us first note that for $\beta =1$, the integral in (38) is always convergent and gives the following representation for the fundamental solution to the space-fractional diffusion equation:

$$G_{\alpha ,1,n}(\mathbf{x},t)=\frac{{\left|\mathbf{x}\right|}^{1-\frac{n}{2}}}{{(2\pi )}^{n/2}}{\int }_0^{\infty }{\sigma }^{n/2}{J}_{\frac{n}{2}-1}(|\mathbf{x}|\sigma )exp(-{\sigma }^{2\alpha }t)d\sigma .$$

We observe, however, that if $\beta <1$, the integral in (38) is convergent only for very limited ranges for the values of the other two parameters. Indeed, according to the asymptotic expansions of the Bessel and the Mittag-Leffler functions, (25) and (20), the integral in (38) is convergent only in the following cases: $n=1$ and $\alpha >1/2$ or $n=2$ and $\alpha >3/4$. If $n\ge 3$, the integral is divergent for any $\alpha \in (0,1)$. Our aim here is to derive from (38) convergent integral representations for $n=1,2,3$, which hold for all $\alpha ,\beta \in (0,1)$.

First, let $n=1$. Plugging in (38), the representation for $J_{-\frac{1}{2}}(·)$ from (24) yields:

$$G_{\alpha ,\beta ,1}(x,t)=\frac{1}{\pi }{\int }_0^{\infty }cos(|x|\sigma )E_{\beta }(-{\sigma }^{2\alpha }{t}^{\beta })d\sigma ,$$

(39)

which, according to (20), is convergent at $+\infty$ only if $2\alpha >1$, unless $\beta =1$. However, we can improve the convergence by performing integration by parts in (39). We use the identity

$$\frac{d}{d\sigma }{E}_{\beta }(-{\sigma }^{2\alpha }{t}^{\beta })=-\frac{2\alpha }{\beta }{\sigma }^{2\alpha -1}{t}^{\beta }{E}_{\beta ,\beta }(-{\sigma }^{2\alpha }{t}^{\beta }),$$

(40)

which is derived from (21). In this way, the following integral representation is established:

$$G_{\alpha ,\beta ,1}(x,t)=\frac{2\alpha }{\beta }\frac{t^{\beta }}{\pi \left|x\right|}{\int }_0^{\infty }sin(|x|\sigma ){\sigma }^{2\alpha -1}{E}_{\beta ,\beta }(-{\sigma }^{2\alpha }{t}^{\beta })d\sigma .$$

(41)

The asymptotic Expression (20) indicates that the integral in (41) is convergent for all $0<\alpha ,\beta \le 1$. In the particular case $\alpha =\beta /2$, representation (41) can also be found in [15], Equation 4.13.

Representation (41) was used for the numerical evaluation of the one-dimensional fundamental solution, and the results are given in Figure 3. For numerical computation of the Mittag-Leffler function in (41), the integral representation (22) was used. Figure 3 shows that the numerical results based on Formula (41) are in good agreement with the exact solutions, (4) and (6).

Next, let us consider $n=3$. Plugging in the general Formula (38), the representation for $J_{\frac{1}{2}}(·)$ from (24) yields:

$$G_{\alpha ,\beta ,3}(\mathbf{x},t)=\frac{1}{2{\pi }^2\left|\mathbf{x}\right|}{\int }_0^{\infty }\sigma sin(|\mathbf{x}|\sigma )E_{\beta }(-{\sigma }^{2\alpha }{t}^{\alpha })d\sigma .$$

This integral is divergent for all $0<\alpha ,\beta <1$. Integration by parts gives

$$G_{\alpha ,\beta ,3}(\mathbf{x},t)=\frac{1}{2{\pi }^2{\left|\mathbf{x}\right|}^2}{\int }_0^{\infty }cos(|\mathbf{x}|\sigma )\frac{d}{d\sigma }\left(\sigma E_{\beta }(-{\sigma }^{2\alpha }{t}^{\alpha })\right)d\sigma$$

and, by applying Formula (40), we obtain the following integral expression for the three-dimensional fundamental solution

$$G_{\alpha ,\beta ,3}(\mathbf{x},t)=\frac{1}{2{\pi }^2{\left|\mathbf{x}\right|}^2}{\int }_0^{\infty }cos(|\mathbf{x}|\sigma )F_{\alpha ,\beta }(\sigma ,t)d\sigma ,$$

(42)

where

$$F_{\alpha ,\beta }(\sigma ,t)=E_{\beta }(-{\sigma }^{2\alpha }{t}^{\beta })-\frac{2\alpha }{\beta }{\sigma }^{2\alpha }{t}^{\beta }{E}_{\beta ,\beta }(-{\sigma }^{2\alpha }{t}^{\beta }).$$

(43)

The asymptotic Expansions (20) of the Mittag-Leffler functions imply that the integral in (42) is convergent for $1/2<\alpha <1$ and $0<\beta \le 1$. Again applying integration by parts in (42) yields

$$G_{\alpha ,\beta ,3}(\mathbf{x},t)=\frac{1}{2{\pi }^2{\left|\mathbf{x}\right|}^3}{\int }_0^{\infty }sin(|\mathbf{x}|\sigma )H_{\alpha ,\beta }(\sigma ,t)d\sigma ,$$

(44)

where $H_{\alpha ,\beta }(\sigma ,t)=-\frac{d}{d\sigma }{F}_{\alpha ,\beta }(\sigma ,t)$ and therefore, by (43) and (21),

$$H_{\alpha ,\beta }(\sigma ,t)=\frac{2\alpha }{\beta }{\sigma }^{2\alpha -1}{t}^{\beta }\left(\left(1+\frac{2\alpha }{\beta }\right)E_{\beta ,\beta }(-{\sigma }^{2\alpha }{t}^{\beta })+\frac{2\alpha }{\beta }{E}_{\beta ,\beta -1}(-{\sigma }^{2\alpha }{t}^{\beta })\right).$$

(45)

The asymptotic behavior of the Mittag-Leffler functions (20) implies that the integral in (44) is convergent for all $0<\alpha ,\beta <1$.

In an analogous way, for $n=2$, we deduce from (38) and (24)

$$G_{\alpha ,\beta ,2}(\mathbf{x},t)=-\frac{1}{2{\pi }^2{\left|\mathbf{x}\right|}^2}{\int }_0^{\pi }\frac{1}{{cos}^2\theta }\left(1+{\int }_0^{\infty }cos(|\mathbf{x}|\sigma cos\theta )H_{\alpha ,\beta }(\sigma ,t)d\sigma \right)d\theta ,$$

(46)

where the function $H_{\alpha ,\beta }$ is defined in (45). The integral in (46) is convergent for all $0<\alpha ,\beta <1$.

It is verified numerically that integral representations (41), (46) and (44) for the two- and three-dimensional fundamental solutions are in agreement with the exact Solutions (5) and (7). For the numerical computation of the Mittag-Leffler functions in $H_{\alpha ,\beta }$, the integral representation (22) is used.

Numerical results and comparison of the one- and two-dimensional fundamental solutions are given in Figure 4.

For a discussion on other integral representations for the fundamental solution, we refer to [7,12].

All numerical computations in this work were performed with the help of MATLAB.

5. Concluding Remarks

The subordination principle for space-time fractional diffusion equations is a useful tool for finding integral representations of the fundamental solution. The derived integral representations (41), (46) and (44) for $n=1,2,3$, respectively, are appropriate for numerical implementation. The performed numerical experiments confirm that the analytical findings in this work are in agreement with the known exact solutions.

The technique used in the present work for deriving the integral representation for the subordination kernel does not rely on the scaling property and can be extended to equations with more general nonlocal operators in space, such as those considered in [36], as well as operators with a general memory kernel in time, as in [24,37].