Green's function estimates for time fractional evolution equations

TL;DR

This paper derives comprehensive estimates for Green's functions of time-fractional evolution equations involving Caputo derivatives and various spatial operators, enabling better understanding of their solutions and behaviors.

Contribution

It provides new global and local two-sided estimates for Green's functions of fractional evolution equations with variable coefficients and extends results to stable-like operators.

Findings

Established global two-sided estimates for elliptic operators.

Derived upper bounds for pseudo-differential operators.

Extended estimates to variable coefficient stable-like operators.

Abstract

We look at estimates for the Green's function of time-fractional evolution equations of the form $D^{\nu}_{0+*} u = Lu$, where $D^{\nu}_{0+*}$ is a Caputo-type time-fractional derivative, depending on a L\'evy kernel $\nu$ with variable coefficients, which is comparable to $y^{-1-\beta}$ for $\beta \in (0, 1)$, and $L$ is an operator acting on the spatial variable. First, we obtain global two-sided estimates for the Green's function of $D^{\beta}_0 u = Lu$ in the case that $L$ is a second order elliptic operator in divergence form. Secondly, we obtain global upper bounds for the Green's function of $D^{\beta}_0 u=\Psi(-i\nabla)u$ where $\Psi$ is a pseudo-differential operator with constant coefficients that is homogeneous of order $\alpha$. Thirdly, we obtain local two-sided estimates for the Green's function of $D^{\beta}_0 u = Lu$ where $L$ is a more general non-degenerate second order elliptic operator. Finally we look at the case of stable-like operator, extending the second result from a constant coefficient to variable coefficients. In each case, we also estimate the spatial derivatives of the Green's functions. To obtain these bounds we use a particular form of the Mittag-Leffler functions, which allow us to use directly known estimates for the Green's functions associated with $L$ and $\Psi$, as well as estimates for stable densities. These estimates then allow us to estimate the solutions to a wide class of problems of the form $D^{(\nu, t)}_0 u = Lu$, where $D^{(\nu, t)}$ is a Caputo-type operator with variable coefficients.

Equations

$$D_{\ast }^{\beta }u=\Delta u,$$

$${\psi }_{\alpha }(\xi )=-|\xi {|}^{\alpha }w(\xi /|\xi |),\xi \in {\mathbb{R}}^d,\alpha \in (0,2),w\in C^k({\mathbb{S}}^{d-1}),$$

$$\left\{\begin{array}[]{rll}D^{\beta}_{*}u(t,x)=&Lu(t,x),&\text{in }(0,\infty)\times\mathbb{R}^{d}\\ u(t,x)=&Y(x),&\text{in }\{0\}\times\mathbb{R}^{d},\end{array}\right.$$

$$-D^{(\nu ,t)}u(t,x)=Lu(t,x),$$

$$Lu=\nabla \cdot (A(x)\nabla u),$$

$$G^{(\beta )}(t,x,y)\asymp C\min \left(t^{-\frac{d\beta }{2}}{\Omega }^{1-\frac{d}{2}},t^{-\frac{d\beta }{2}}{\Omega }^{-\frac{d}{2}\left(\frac{1-\beta }{2-\beta }\right)}\exp \{-C_{\beta }{\Omega }^{\frac{1}{2-\beta }}\}\right),$$

$${\psi }_{\alpha }(\xi )=|\xi {|}^{\alpha }{w}_{\mu }(\xi /|\xi |),\xi \in {\mathbb{R}}^d,\alpha \in (0,2),w_{\mu }\in C^k({\mathbb{S}}^{d-1}),$$

$$G_{{\psi }_{\alpha }}^{(\beta )}(t,x-y)\asymp C\min \left(t^{-\frac{d\beta }{\alpha }}{\Omega }^{1-\frac{d}{\alpha }},t^{-\frac{d\beta }{\alpha }}{\Omega }^{-1-\frac{d}{\alpha }}\right)$$

$$Lu=\sum _{i,j=1}^d{a}_{ij}(x){\partial }_{x_i}{\partial }_{x_j}u+\sum _{i=1}^d{b}_i(x){\partial }_{x_j}u+c(x)u.$$

$${\psi }_{\alpha ,x}(x,\xi )=|\xi {|}^{\alpha }{w}_{\mu }(x,\xi /|\xi |),$$

$$D_{0+\ast }^{(\nu )}f(t)=-{\int }_0^t(f(t-s)-f(t))\nu (t,\mathrm{d}s)-{\int }_t^{\infty }(f(0)-f(t))\nu (t,\mathrm{d}s).$$

$$cg(x)\le f(x)\le Cg(x),\forall x\in D,$$

$$\frac{f(x)}{g(x)}\to 1,\text{ as }x\to \infty .$$

$$cg(x)\le f(x)\le Cg(x),x\in (M,\infty ).$$

$$\frac{f(x)}{g(x)}\to 1,\text{ as }x\to 0.$$

$$cg(x)\le f(x)\le cg(x),x\in (0,m).$$

$${\partial }_{t}u-\{a_{ij}(t,x){\partial }_{x_i}{\partial }_{x_j}u+b_i(t,x){\partial }_{x_i}u+c(t,x)u\}=0,u(0,x)=\delta (x-\xi )$$

$$Z(t,x;\tau ,\xi )\asymp (t-\tau {)}^{-d/2}\exp \left\{-c\frac{(x-\xi {)}^2}{t-\tau }\right\}.$$

$${\partial }_{t}u={\partial }_{x_i}(a_{ij}(x){\partial }_{x_j}u).$$

$${\mu }^{-1}|\xi {|}^2\le a_{ij}(x){\xi }_i{\xi }_j\le \mu |\xi {|}^2,\text{ for all }\xi \in {\mathbb{R}}^d.$$

$$c_1{t}^{-d/2}\exp \left\{-c_2\frac{|x-y{|}^2}{t}\right\}\le G(t,x,y)\le C_1{t}^{-d/2}\exp \left\{-C_2\frac{|x-y{|}^2}{t}\right\}.$$

$${\varphi }_{\beta }(y)=\exp \{-\sigma |y{|}^{\beta }{e}^{i\frac{\pi }{2}\gamma \text{ sgn }\gamma }\},y\in \mathbb{R},$$

$$w_{\beta }(x;\gamma ,\sigma )=\frac{1}{2\pi }{\int }_{\mathbb{R}}\exp \{-ixy-\sigma |y{|}^{\beta }{e}^{i\frac{\pi }{2}\gamma \text{ sgn }\gamma }\}\mathrm{d}y.$$

$$w_{\beta }(x)=\frac{1}{\pi }\Re {\int }_0^{\infty }\exp \{-ixy-|y{|}^{\beta }\}\mathrm{d}y,$$

$$w_{\beta}(x)\sim\tilde{c}_{\beta}\left\{\begin{array}[]{lc}x^{-1-\beta},&\text{ as }x\rightarrow\infty,\\ f_{\beta}(x):=x^{-\frac{2-\beta}{2(1-\beta)}}\exp\left\{-c_{\beta}x^{-\frac{\beta}{1-\beta}}\right\},&\text{ as }x\rightarrow 0,\end{array}\right.$$

$$w_{\beta }(x)\asymp c_{\beta }\min \left(x^{-1-\beta },x^{-\frac{2-\beta }{2(1-\beta )}}\exp \left\{-c_{\beta }{x}^{-\frac{\beta }{1-\beta }}\right\}\right),$$

$$w_{\beta}(x)\asymp\tilde{c}_{\beta}\left\{\begin{array}[]{lc}x^{-1-\beta},&\text{ for }x\in(1,\infty),\\ x^{-\frac{2-\beta}{2(1-\beta)}}\exp\left\{-c_{\beta}x^{-\frac{\beta}{1-\beta}}\right\},&\text{ for }x\in(0,1).\end{array}\right.$$

$${\varphi }_{\alpha }(p)=\exp \left\{-|p{|}^{\alpha }{\int }_{{\mathbb{S}}^{d-1}}|(p/|p|,s){|}^{\alpha }\mu (\mathrm{d}s)\right\},$$

$$w_{\mu }(p)={\int }_{{\mathbb{S}}^{d-1}}|(p,s){|}^{\alpha }\mu (\mathrm{d}s),$$

$${\psi }_{\alpha }(p):=\log {\varphi }_{\alpha }(p)=-|p{|}^{\alpha }{w}_{\mu }(p/|p|),p\in {\mathbb{R}}^d.$$

Full text

Keywords: Caputo derivative; Green’s function; Aronson estimates; two-sided estimates; fractional evolution

MSC 2010: Primary: 60J35. Secondary: 60G52; 35J08; 33E12; 44A10; 60J75.

1. 1 Introduction

2. 2 Preliminaries

3. 2.1 Estimates and Stable Processes

4. 2.2 Fractional Derivatives and Their Extensions

5. 3 Global Estimates

6. 3.1 Time-Fractional Diffusion Equation: Divergence Structure

7. 3.2 Time-Fractional Pseudo-Differential Evolution: Constant Coefficients

8. 4 Local Estimates

9. 4.1 Time-Fractional Diffusion Equation: General Non-Degenerate

10. 4.2 Time-Fractional Pseudo-Differential Evolution: Variable Coefficients

11. 5 Generalised Evolution Equations

12. 6 Conclusion

13. 7 Asymptotic Methods

1 Introduction

The area of fractional evolution equations has become hugely popular in recent decades, due to its ability to better model real-world phenomena compared to their non-fractional counterpart which usually model local behaviour. The nature of fractional in-time operators (respectively in space) allow us to model, for example, processes that exhibit some kind of memory (resp. non-local interactions). The processes associated with time-fractional evolution models (in the case of a local spatial operator, they are time-changed diffusion processes), possess some remarkable properties. For some motivation, let us focus for the moment on the particular case of

[TABLE]

where $D^{\beta}_{*}$ is the Captuo fractional derivative in time, $\beta\in(0,1)$ and $\Delta$ is the Laplacian operator (a second order uniformly elliptic operator). This time-fractional diffusion equation is widely used to model anomalous diffusions which exhibit subdiffusive behaviour, which is due to the diffusive particles being trapped. Such fractional time diffusion equations also arise as a scaling limit of random conductance models (random walks in random environments). This point of view is particularly interesting, since the limiting process is a non-Markovian process which arises as the scaling limit of Markovian process, see [barlow2011convergence, meerschaert2012stochastic, leonenko2013correlation] and references therein for a probabilistic account of models related to fractional calculus. Recently in [hairer2018fractional], the authors discussed how a fractional kinetic process (with $\beta=\frac{1}{2}$) emerges as the intermediate time behaviour of perturbed cellular flows. For an extensive account of physical applications see [richard2014fractional], [mainardi2010fractional] or [tarasov2011fractional]. For early applications of continuous time random walks and fractional calulus arising in finance, see the series of articles [scalas2000fractional, mainardi2000fractional, gorenflo2001fractional].

Recently much attention has been given to the Green’s function of fractional differential equations. In [chen2018heat], the authors obtain two-sided estimates for the Green’s functions of fractional evolution equations, under the assumption that the Green’s function of the spatial operator satisfies global (in time) two-sided estimates.

In [grigor2008dichotomy], the authors explore the general structure of two-sided estimates for the transition probabilities associated with local or non-local Dirichlet forms. They show that the bounds for transition probabilities associated with local Dirichlet forms will always be of exponential type, and for those associated with non-local Dirichlet forms the bounds will be of polynomial type. Even more recently in [deng2018exact] they give some exact asymptotic formulas for the Green’s function of fractional evolution equations. The authors in [kelbert2016weak] study error estimates for continuous time random walk (CTRW) approximation of classical fractional evolution equations, for which the heat kernel estimates for $D^{\beta}u=\Delta u$ and $D^{\beta}u=\psi(-i\nabla)u$, where $\psi(-i\nabla)$ generates a symmetric stable process, are obtained as a by-product.

In [eidelman2004cauchy] the authors use the parametrix method (or Levi method) to study the equation $D^{\beta}u(t,x)-Bu(t,x)=f(t,x)$, where the operator $B$ is a uniformly elliptic second order differential operator (which we look at in Theorem 4.1) with bounded continous real-valued coefficients. This is done by looking first at the constant coefficient case then using these estimates to study the variable coefficient case. In the articles [kochubei2019random, kochubei2018random], the authors study the long-time behaviour of the Cesaro mean of the heat kernel of subordinated processes and for this they use a version of a Karamata-Tauberian theorem. The long-time behaviour of solutions to space-time fractional diffusion equations are also considered in [cheng2017asymptotic]. In the articles [kemppainen2017representation, cheng2017asymptotic, kim2016asymptotic], the authors consider the space-time fractional diffusion equation, which involves a Caputo fractional derivative in time and a fractional Laplacian as the spatial operator. This falls under the case of our Theorem 3.3, where we consider as the spatial operator a (non-isotropic) pseudo-differential operator with symbol

[TABLE]

see (2.8). Since this case involves constant coefficients, we could have used the machinery of the Fox H function like they do in [kemppainen2017representation, kim2016asymptotic] (among many others that use them) or Laplace transform arguments. However these approaches would not work for our other results, since the operators involved have variable coefficients.

Diffusion processes in random environments are also closely related objects, and in fact there are many works looking at estimates for the heat kernels of such processes; for example in [cabezas2015sub], the authors obtain sub-gaussian bounds for the transition kernel of a random walk in a random environment.

This article is structured as follows: in Section 2 we begin by recalling definitions that will be used throughout the article. The topic of generalised fractional calculus is briefly covered, which is motived by a probabilistic generalisation of Caputo fractional derivatives. We also recall some important estimates, namely the Aronson estimates for the fundamental solutions of second order parabolic equations and asymptotic estimates for stable densities. In Section 3 (resp. Section 4) we obtain global estimates (resp. local estimates) for the Green’s function of fractional evolution equations of the form

[TABLE]

where $L$ is some differential operator acting on the spatial variable $x$. In the final section we discuss solutions to generalised evolution equations of the form

[TABLE]

where $D^{(\nu,t)}$, acting on the time variable, is the generator of an increasing process, which is comparable to a $\beta$-stable subordinator, and $L$ is a generator of a strongly continuous contraction semigroup. We then conclude by summarising the main results and commenting on some applications.

We summarise briefly some of the results obtained for various spatial operators.

• •

Theorem 3.1: When the spatial operator is given by a second order uniformly elliptic operator in divergence form,

[TABLE]

we obtain the following two-sided estimates for the Green’s function $G^{(\beta)}(t,x,y)$ of (1.1). For $d\geq 3$ and $\beta\in(0,1)$ we find:

[TABLE]

for $(t,x,y)\in(0,\infty)\times\mathbb{R}^{d}\times\mathbb{R}^{d}$ where $\Omega:=|x-y|^{2}t^{-\beta}$.

• •

Theorem 3.3: When the spatial operator is a non-isotropic pseudo-differential operator $\Psi_{\alpha}$, with homogeneous symbol of the form (cf. Equation (2.7))

[TABLE]

we obtain two-sided estimates for the Green’s function $G^{(\beta)}_{\psi_{\alpha}}$ of (1.1). For $d>\alpha>0$ and $\beta\in(0,1)$:

[TABLE]

for $(t,x-y)\in(0,\infty)\times\mathbb{R}^{d}$, where $\Omega:=|x-y|^{\alpha}t^{-\beta}$.

It is well known that the Green’s functions of the above operators $L$ and $\Psi_{\alpha}$ satisfy two-sided estimates for all $(t,x,y)\in(0,\infty)\times\mathbb{R}^{d}\times\mathbb{R}^{d}$, see [aronson1967bounds] or [Stroock1988] for the first case, [kolokoltsov2019differential] or [eidelman2004analytic] for the second. Note that in the case of a lower bound for the Green’s function of $\Psi_{\alpha}$, one also needs to assume that its spectral density $\mu$ is strictly positive, and restrict $\alpha\in(0,2)$.

In Section 4 we consider more general spatial operators, whose Green’s functions satisfy only local (in time) estimates. We obtain local estimates for the Green’s function of (1.1) in the following cases:

• •

Theorem 4.1: A general second order elliptic operator $L$ with variable coefficients of the form

[TABLE]

Unsurprisingly we find that the estimates are the same as (1.2) but only for $(t,x,y)\in(0,T)\times\mathbb{R}^{d}\times\mathbb{R}^{d}$ for some fixed $T>0$.

• •

Theorem 4.6: A non-isotropic pseudo-differential operator $\psi_{\alpha,x}$ with variable coefficients, with homogeneous symbol of the form

[TABLE]

where $w_{\mu}$ is for each fixed $x$ a continuous function on the surface of the sphere $\mathbb{S}^{d-1}$. Note again that the spectral density associated with $\psi_{\alpha}$ (for fixed $x$) must be strictly positive, in order to use two-sided estimates for the Green’s function $G_{\psi_{\alpha},x}$. Again we obtain the same estimates as (1.3) but for $(t,x,y)\in(0,T)\times\mathbb{R}^{d}\times\mathbb{R}^{d}$ for some fixed $T>0$.

For each of the four cases above, we also obtain estimates for the spatial derivatives of the Green’s functions of the fractional evolution equations.

Finally, in the last section we turn our attention to generalised fractional evolutions. This is when one replaces the standard fractional time derivative with a weighted mixture of fractional derivatives,

[TABLE]

The main result is that the solution to such generalised fractional evolutions can be estimated (by means of our Green’s function estimates) by the solutions to classical fractional evolution equations.

2 Preliminaries

2.1 Estimates and Stable Processes

Throughout the article, we will use the notation $f(x)\asymp g(x)$ in $D$, which means that there exists constants $C,c>0$ such that $f$ satisfies the following two-sided estimate,

[TABLE]

for some region $D$. The notation $f(x)\sim g(x)$ for $x\rightarrow\infty$ means that

[TABLE]

Then for each $M>0$ there exists a constant $C>0$ such that

[TABLE]

Similarly, the notation $f(x)\sim g(x)$ for $x\rightarrow 0$ means

[TABLE]

Then for each $m>0$ there exists a $c>0$ such that

[TABLE]

If both $f$ and $g$ on $\mathbb{R}_{+}$ are positive, bounded and satisfy $f(x)\sim g(x)$ for $x\rightarrow\infty$ (resp. $x\rightarrow 0$), then $f(x)\asymp g(x)$ in $(M,\infty)$ for any $M>0$ (resp. in $(0,m)$ for any $m<\infty$). See the excellent De Bruijn book [de1970asymptotic] for more details on asymptotic analysis.

The fundamental solutions $Z(t,x;\tau,\xi)$ of the Cauchy problem for the uniformly parabolic equations

[TABLE]

with bounded and uniformly Hölder continuous coefficients in $x$ defined on $(0,T]\times\mathbb{R}^{d}$ are known, [porper1984two], to satisfy the two-sided estimate

[TABLE]

We will assume that the coefficients do not depend on time, so that the fundamental solution is just a function of $t,x$ and $y$. On the other hand, Aronson [aronson1967bounds] obtained global two-sided estimates for the fundamental solution $G(t,x,y)$ of the divergence equation

[TABLE]

Assume that the coefficients satisfy the uniform ellipticity condition: there exists $\mu\geq 1$ such that

[TABLE]

Further assuming that the coefficients in (2.2) are continuous, then there exists constants $C_{1},C_{2},c_{1}$ and $c_{2}$ such that for $(t,x,y)\in(0,\infty)\times\mathbb{R}^{d}\times\mathbb{R}^{d}$,

[TABLE]

For a discussion on divergence and non-divergence equations, see for example [evans2010partial].

Note that throughout $C$, $\tilde{C}$, $c$ or $\tilde{c}$ denotes some constants. If we wish to stress what they depend on, say $\alpha,\beta$ or $T$, we write $C_{\alpha,\beta,T}$ for example.

Let us recall some basic facts about stable densities; our standard references for these are [zolotarev1986one, kolokoltsov2011markov]. The characteristic function of the general (up to a shift) one-dimensional stable law with index of stability $\beta\in(0,2)$ (but $\beta\neq 1$) is given by

[TABLE]

where the parameter $\gamma\in[-1,1]$ measures the skewness of the distribution and $\sigma>0$ is the scale parameter. The probability density corresponding to the characteristic function $\phi_{\beta}$, which we denote by $w_{\beta}(x;\gamma,\sigma)$, is given by the following Fourier transform:

[TABLE]

We will be using totally positively skewed ($\gamma=1$) normalised ($\sigma=1$) stable densities, which we denote by $w_{\beta}(x)$ and they are given by

[TABLE]

where $\Re(z)$ is the real part of $z\in\mathbb{C}$. We will be using the asymptotic behaviour (as $x\rightarrow 0$ and $x\rightarrow\infty$) of the stable densities $w_{\beta}(x)$, so we state them now, [uchaikin2011chance] (Theorem 5.4.1).

Proposition 2.1**.**

The stable densities $w_{\beta}(x)$ have the following asymptotic behaviour

[TABLE]

where $c_{\beta}=(1-\beta)\beta^{\frac{\beta}{1-\beta}}$.

A key point to note from the above asymptotic behaviour, is that since $w_{\beta}$ is bounded and strictly positive for $x>0$, this gives us a two-sided estimate for $w_{\beta}(x)$ for all $x\in(0,\infty)$. More precisely for $\beta\in(0,1/2)$ we have

[TABLE]

and for $\beta\in[1/2,1)$

[TABLE]

For $\alpha\in(0,2)$, the general symmetric stable density in $\mathbb{R}^{d}$ (up to a shift) has the form

[TABLE]

where the (finite) measure $\mu$ on $\mathbb{S}^{d-1}$ is called the spectral measure, [kolokoltsov2000symmetric]. Let $w$ be a function on $\mathbb{S}^{d-1}$ given by

[TABLE]

so that

[TABLE]

Note that $\psi_{\alpha}$ is the symbol of a pseudo-differntial operator $\Psi_{\alpha}(-i\nabla)$ which we will study later. When $\mu$ is the uniform measure on $\mathbb{S}^{d-1}$ we have $w(p)\equiv 1$, and the operator $\Psi_{\alpha}(-i\nabla)$ is just the fractional Laplacian $(-\Delta)^{\alpha}$ with symbol $\psi_{\alpha}(\xi)=-|\xi|^{\alpha}$.

2.2 Fractional Derivatives and Their Extensions

Fractional derivatives at $a\in\mathbb{R}$ of order $\beta\in(0,1)$ can be viewed probabilistically as the generators of a $\beta$-stable process, interrupted on an attempt to cross a boundary point $a$. This interpretation was extensively explored and extended in [kolokoltsov2015fully], and we recall some of the details here. For a broader background in classical fractional calculus, see any text on fractional calculus [samko1993fractional, diethelm2010analysis, kilbas2006theory, podlubny1998fractional], for example.

The fractional Riemann-Liouville (RL) integral of order $\beta\in(0,1)$ is given by

[TABLE]

and the Caputo fractional derivative of order $\beta\in(0,1)$ is given by

[TABLE]

The fractional derivative in generator from, is written as

[TABLE]

which equals the Caputo derivative for $a=-\infty$. A possible extension for these fractional derivative operators which is widely used in the literature, are various mixtures of these derivatives, for example

[TABLE]

In this article we will look at weighted mixed fractional derivatives, given by

[TABLE]

with some positive kernel $\nu(t,\cdot)$ on $\{s:t>0\}$ satisfying the one-sided Lévy-condition

[TABLE]

The $-$ sign in the definition of $D^{(\nu)}_{+}$ is to comply with the standard notation for fractional derivatives, so that $D^{\beta}_{+}=D^{(\nu)}_{+}$ with $\nu(t,y)=-1/\left[\Gamma(-\beta)y^{1+\beta}\right]$.

Notice that the Caputo derivative $D^{\beta}_{a+*}$ is obtained from $D^{\beta}_{+}$ by the restriction of its action on the space $C^{1}([a,\infty))$ considered as the subpsace of $C(\mathbb{R})$ by extending their values as constants to the left of $a$. Then looking for a generalised Captuo derivative arising from $D^{(\nu)}_{+}$ we define

[TABLE]

These extensions were initially suggested by the second author in [kolokoltsov2015fully], on the basis of the following probabilistic interpretation. It is well known that the fractional derivatives $-d^{\beta}/dx^{\beta}$ for $\beta\in(0,1)$ generate of stable Lévy subordinators with the inverted direction (i.e., decreasing processes). Then the Caputo derivatives $D^{\beta}_{a+*}$ describe the modifications of the stable subordinators obtained by forbidding them to cross the boundary $x=a$ for some $a\in\mathbb{R}$. Applying this modification to the generalised operators $D^{(\nu)}_{+}$ one obtains $D^{(\nu)}_{a+*}$ which are then the generators of a Markov process interrupted on an attempt to cross the boundary point $a$.

The corresponding generalised fractional integrals arising from the generalised fractional derivative $D^{(\nu)}_{+}$ can be defined in a few different ways depending on which point of view one chooses: probability, semigroup theory or generalised functions/$\Psi$DE theory. The objects end up being the same, but we will stick to the semigroup theory prespective in this article, see the upcoming review [kolokoltsov2019mixed] for a full account. In terms of operator semigroups, the operator $I^{\beta}_{-\infty}$ is the potential operator of the semigroup generated by $-D^{\beta}_{+}$. In other words, it is the limit of the resolvent operator

[TABLE]

as $\lambda\rightarrow 0$. Then since $I^{\beta}_{a}$ is the reduction of $I^{\beta}_{-\infty}$ to the space $C_{kill(a)}([a,\infty))$ of functions vanishing to the left of $a$ for any $a\in\mathbb{R}$, we define the generalised fractional integral $I^{\nu}_{a}$ as the potential operator of the semigroup generated by $-D^{(\nu)}_{+}$ reduced to the space $C_{kill(a)}([a,\infty))$. For a background on potential operators and measures, see [schilling2012bernstein] or [van2012potential] (or from a probabilistic point of view, [feller2008introduction]).

Let us denote by $(T^{\nu}_{t})_{t\geq 0}$ the semigroup generated by the operator $-D^{(\nu)}_{+}$. Then for $f\in C_{kill(a)}([a,\infty))\cap C_{\infty}(\mathbb{R})$, the potential operator $U^{(\nu)}$ of the semigroup $T^{\nu}_{t}$ is given by

[TABLE]

where $G_{(\nu)}(r,t,\mathrm{d}s)$ are the transition probabilities of the process generated by $D^{(\nu)}_{+}$. The potential measure is defined as the integral kernel of the potential operator, and by an abuse of notation, we denote this measure by $U^{(\nu)}(t,\mathrm{d}s)$. Thus the generalised fractional integral $I^{(\nu)}_{a}$ is given by

[TABLE]

where the potential measure $U^{(\nu)}(t,\mathrm{d}s)$ is equal to the vague limit

[TABLE]

of the measures $\int_{0}^{K}G_{(\nu)}(r,t,\cdot)\mathrm{d}r,~{}K\rightarrow\infty$ (see [schilling2012bernstein] (p. 63)). Furthermore the $\lambda$-potential measure is defined by

[TABLE]

so that if $\lambda>0$ and $g\in C_{kill(a)}([a,\infty))\cap C_{\infty}(\mathbb{R})$, the convolution $(U^{(\nu)}_{\lambda}\star g)(t)$, which is given by

[TABLE]

is the resolving operator of the semigroup generated by $-D^{(\nu)}_{+}$. That is, $f(x)=(U^{(\nu)}_{\lambda}\star g)(x)$ is the classic solution to the equation

[TABLE]

This also holds for $\lambda=0$, and so the potential operator with kernel $U^{(\nu)}(t,\mathrm{d}y)$,

[TABLE]

represents the classical solution to the equation

[TABLE]

on $C_{kill(a)}([a,b])$.

Example 2.2**.**

For the case $\nu(t,\mathrm{d}y)=-1/[\Gamma(-\beta)y^{1+\beta}]\mathrm{d}y$, (2.9) says that

[TABLE]

where $G_{\beta}(r,s)$ are the transition densities of a $\beta$-stable subordinator, is the solution to the linear fractional equation

[TABLE]

On the other hand, it is well known that the solution to such linear fractional equations are given by

[TABLE]

where $E_{\beta}(z)$ is the Mittag-Leffler function

[TABLE]

Thus we have

[TABLE]

which is equivalent to the Zolotarev-Pollard formula for the Mittag-Leffler function in terms of the transition densities of stable subordinators, which is of great importance to this article,

[TABLE]

For the rest of this article we will take $a=0$ and write $D^{\beta}_{0}:=D^{\beta}_{0+*}$ and $D^{(\nu)}_{0}:=D^{(\nu)}_{0+*}$ to simplify formulas, but keep in mind that the boundary point [math] can be replaced with $a\in\mathbb{R}$ by a shift in the time coordinate.

As noted in [kolokoltsov2017chronological], we can extrapolate from the case

[TABLE]

to the Banach-valued version

[TABLE]

where $L$ is some operator generating a Feller semigroup. One can expect that the solution to this equation can be written in terms of an operator-valued Mittag-Leffler function,

[TABLE]

where $E_{\beta}(s)$ are Mittag-Leffler functions defined by (2.10). However, this series representation does not allow one to define $E_{\beta}(L)$ for an unbounded operator $L$. In both [kolokoltsovwell2014] and [kolokoltsov2017chronological] the authors find that the most convenient way to overcome this difficulty is to use the formula (2.11) for the Mittag-Leffler function. This connection between Mittag-Leffler functions, Laplace transforms and stable densities is due to Zolotarev, [zolotarev1957mellin, zolotarev1961analytic, zolotarev1986one]—although preliminary versions of this formula were also noted almost a decade earlier by Pollard, [pollard1948completely]. Thus formula (2.11) could be called the Pollard-Zolotarev formula. Notice that if an operator $L$ generates a Feller semigroup with transition densities $G(t,x,y)$, then

[TABLE]

With the help of Fubini’s theorem, the solution (2.12) can be written as

[TABLE]

Thus the main aim of this article is to obtain estimates for the Green’s function given by,

[TABLE]

where $G_{L}(z,x,y)$ is the Green’s function associated with the spatial operator $L$, i.e., the fundamental solution of

[TABLE]

3 Global Estimates

We first look at global in time two-sided estimates for $G^{(\beta)}_{L}$ in two special cases. Notice in (2.13) that the integral over the time variable $z$ ranges from [math] to $\infty$, and so in order to perform any estimates on the term $G_{L}(z,x,y)$ one can only use estimates that hold for all $z\in(0,\infty)$. We begin with two such cases, when one has global in time estimates for $G_{L}$. Namely, when $L$ is a second order uniformly elliptic operator in divergence form or when $L$ is a homogeneous pseudo-differential operator (with constant coefficients).

3.1 Time-Fractional Diffusion Equation: Divergence Structure

In this section we consider the time-fractional diffusion equation given by

[TABLE]

where $D^{\beta}_{0}$ is the Caputo fractional derivative acting on the time variable, and the spatial operator is a second order elliptic operator in divergence form. For the conditions on the diffusion coefficient $A$ such that $L$ generates a Feller semigroup, see [Stroock1988]. Recall that the solution of (3.1) is given by

[TABLE]

and the associated Greens function is given by

[TABLE]

where $G(t,x,y)$ is the Greens function associated with the second order elliptic operator in divergence form, (2.1). We have the following two-sided estimates for the Green’s function $G^{(\beta)}$, which are global in time. In the following, we use the notation $\Omega:=|x-y|^{2}t^{-\beta}$.

Theorem 3.1**.**

Assume that the function $A(x)$ is measurable, symmetric and satisfies (2.2) for some $\mu\geq 1$. Then there exists a constant $C$ such that for $(t,x,y)\in(0,\infty)\times\mathbb{R}^{d}\times\mathbb{R}^{d}$, the Green’s function $G^{(\beta)}(t,x,y)$ for the time-fractional diffusion Equation (3.1) satisfies the following two-sided estimates,

• •

For $\Omega\leq 1$,

[TABLE]

• •

For $\Omega\geq 1$,

[TABLE]

Proof.

Let us begin by using the estimate (2.4) for the stable density in (3.2),

[TABLE]

where $f_{\beta}(x)=x^{-\frac{2-\beta}{2(1-\beta)}}\exp\{-c_{\beta}x^{-\frac{\beta}{1-\beta}}\}$. Next we apply Aronsons estimates (2.3) to $G(t^{\beta}z,x,y)$,

[TABLE]

where $\Omega:=|x-y|^{2}t^{-\beta}$. Making a change of variables $z=w^{-1}$ so that $\mathrm{d}z=-w^{-2}\mathrm{d}w$,

[TABLE]

We now estimate $I_{1}$ and $I_{2}$ in two different cases, depending on the behaviour of $\Omega$.

Case 1: $\Omega\leq 1$.** Making a further substitution of $V=\Omega w$ in the integral $I_{1}$ gives us the simpler form of**

[TABLE]

Now if $d=1$, then we have the asymptotic behaviour

[TABLE]

and in particular for $\Omega\leq 1$ there exists a constant $C>0$ such that

[TABLE]

If $d=2$, then we see logarithmic behaviour,

[TABLE]

and in particular for $\Omega\leq 1$ there exists a constant $C>0$ such that

[TABLE]

If $d\geq 3$, then the integral is the so-called upper incomplete gamma function, and has the asymptotic behaviour

[TABLE]

and in particular for $\Omega\leq 1$ there exists a constant $C>0$ such that the two-sided estimate

[TABLE]

holds. Thus we have the following two-sided estimate for $I_{1}$,

[TABLE]

Turning to the integral $I_{2}$,

[TABLE]

due to the fast decay of $f_{\beta}$ in a neighbourhood of [math]. Thus combining the estimates for $I_{1}$ and $I_{2}$ gives (3.3).

Case 2: $\Omega\geq 1$.** In this case we use the Laplace method as described in Appendix 7. Firstly for $I_{1}$, using $g(w)=w^{\frac{d}{2}-1}$, $h(w)=w$ and $b=1$ in (7.1) we have**

[TABLE]

and in particular the estimate

[TABLE]

For the second integral, we use Proposition 7.1 with $N=\frac{d}{2}-1-\frac{1}{2(1-\beta)}$ and $a=\frac{1}{1-\beta}$,

[TABLE]

and again in particular, the two-sided estimate

[TABLE]

Combinging the estimates for $I_{1}$ and $I_{2}$ shows (3.4), and we are done.

∎

If one additionally assumes that the diffusion coefficients $A(x)$ of (2.1) are twice continuously differentiable, the following estimates hold for the spatial derivatives of the fundamental solution of (2.1),

[TABLE]

for $(t,x,y)\in(0,\infty)\times\mathbb{R}^{d}\times\mathbb{R}^{d}$. We next have estimates for the spatial deriviatve of $G^{(\beta)}(t,x,y)$.

Proposition 3.2**.**

Under the same assumptions as Theorem 3.1, assume additionally that $A(x)$ is twice continuously differentiable, then the following estimates for the spatial derivatives of the Green’s function $G^{(\beta)}(t,x,y)$ holds for all $(t,x,y)\in(0,\infty)\times\mathbb{R}^{d}\times\mathbb{R}^{d}$,

• •

For $\Omega\leq 1$,

[TABLE]

• •

For $\Omega\geq 1$,

[TABLE]

Proof.

Recall that

[TABLE]

where $G$ satisfies the global estimate (3.8). Using the triangle inequality after taking the derivative inside the integral,

[TABLE]

where in the above calculations, after using the estimates (3.8) and (2.4), we made the substitution $z=w^{-1}$. Note that the integrals $I_{1}$ and $I_{2}$ differ from those appearing in (3.1) only by replacing $d$ with $d+1$. Thus the only change in the calculations is where the dimension dictates the behaviour of the estimate, namely in the integral $I_{1}$ under the regime $\Omega\leq 1$. In this case, make the substitution $w\Omega=V$,

[TABLE]

For $d=1$, we are in the same situation as (3.6), thus

[TABLE]

and in particular

[TABLE]

Otherwise for $d\geq 2$ we have

[TABLE]

For the integral $I_{2}$, replacing $d$ with $d+1$ in (3.1) does not spoil the estimate, thus

[TABLE]

This shows

[TABLE]

for $\Omega\leq 1$ as required. For $\Omega\geq 1$, the estimates follow again by using the Laplace method. Namely taking $g(w)=w^{\frac{d+1}{2}-1}$, $h(w)=w$ and $b=1$ in (7.1) we have

[TABLE]

Finally using $N=\frac{d+1}{2}-1-\frac{1}{2(1-\beta)}$ and $a=\frac{1}{1-\beta}$ in Proposition 7.1 gives us

[TABLE]

Thus

[TABLE]

as required. ∎

3.2 Time-Fractional Pseudo-Differential Evolution: Constant Coefficients

Next we turn our attention to another class of problems, where the spatial operator is a homogeneous (constant coefficient) pseudo-differential operator. That is, for $\beta\in(0,1)$ and $\alpha>0$,

[TABLE]

where $\Psi_{\alpha}$ is a pseudo-differential operator whose symbol is of the form

[TABLE]

where $w_{\mu}$ is a positive function on $\mathbb{S}^{d-1}$, see (2.8). To this end, we use known properties of the Green’s function $G_{\psi_{\alpha}}(t,x)$ of the evolution

[TABLE]

See for example [kolokoltsov2019differential, Theorem 4.5.1] or [eidelman2004analytic] for the following estimates. Assuming that

• •

The function $w_{\mu}$ belongs to $C^{d+1+[\alpha]}(\mathbb{S}^{d-1})$.

• •

The spectral measure $\mu$ has a density which is strictly positive (see (2.6)).

• •

$\alpha\in(0,2)$,

then the Green’s function $G_{\psi_{\alpha}}(t,x-y)$ of the evolution (3.11) satisifies the following two-sided estimates for $(t,x,y)\in(0,\infty)\times\mathbb{R}^{d}\times\mathbb{R}^{d}$,

[TABLE]

Note that the restriction $\alpha\in(0,2)$ and the positivity of the density of the spectral measure is required for the lower bound of $G_{\psi_{\alpha}}$ - the upper bound is still seen if we drop the strict positivity of the density $\mu$ and take any $\alpha>0$.

If additionally $w$ is $(d+1+[\alpha]+l)$-times continuously differentiable, then $G_{\psi_{\alpha}}(t,x)$ is $l$-times continuously differentiable in $x$ and for $(t,x,y)\in(0,\infty)\times\mathbb{R}^{d}\times\mathbb{R}^{d}$,

[TABLE]

for all $k\leq l$ and $i_{1},\cdots,i_{k}$. As discussed in the introduction, the solution of (3.10) is given by

[TABLE]

where

[TABLE]

Thus the corresponding Greens function $G^{(\beta)}_{\psi_{\alpha}}(t,x,y)$ of (3.10) is given by

[TABLE]

In keeping with the previous section, we denote $\Omega:=|x-y|^{\alpha}t^{-\beta}$. We have the following two-sided estimate for the Green’s function $G^{(\beta)}_{\psi_{\alpha}}(t,x,y)$.

Theorem 3.3**.**

Let $\alpha\in(0,2)$ and $\beta\in(0,1)$. Assume that $w\in C^{(d+1+[\alpha])}(\mathbb{S}^{d-1})$, and that $w\geq w_{0}>0$ for some constant $w_{0}$. Further assume that the spectral measure $\mu$ of the stable operator $\Psi_{\alpha}$ has a strictly positive density. Then there exists a constant $C>0$ such that the Green’s function for the fractional evolution Equation (3.10) satisfies the following two-sided estimates for $(t,x-y)\in(0,\infty)\times\mathbb{R}^{d}$.

• •

For $\Omega\leq 1$,

[TABLE]

• •

For $\Omega\geq 1$,

[TABLE]

Proof.

We begin by splitting up the stable density using (2.4),

[TABLE]

where $f_{\beta}(z)$ is defined in (2.4). Before using the estimates (3.12) for $G_{\psi_{\alpha}}$ (with $t=t^{\beta}z$), note that using the notation $\Omega=|x-y|^{\alpha}t^{-\beta}$, we have

[TABLE]

Thus we have,

[TABLE]

Now we deal with two cases.

Case 1: $\Omega\leq 1$.** Using (3.16), in this case the intergral $I_{1}$ equals**

[TABLE]

Note that for $d=\alpha$, the integral over the interval $(\Omega,1)$ is

[TABLE]

On the other hand, for $d\neq\alpha$ we have

[TABLE]

Thus in this case we have,

[TABLE]

Turning to $I_{2}$, note that the integral does not involve $\Omega$, and is convergent since $f_{\beta}(z^{-\frac{1}{\beta}})$ is bounded and vanishes as $z\rightarrow\infty$. Thus

[TABLE]

Combining the estimates for $I_{1}$ and $I_{2}$ shows (3.14).

Case 2: $\Omega\geq 1$.** In this case, the integral $I_{1}$ is simply**

[TABLE]

For the second integral, we have

[TABLE]

Note that the integral in the first term approaches a convergent integral (for large $\Omega$), while the second can be dealt with by using the Laplace method, see 7.1,

[TABLE]

where we have used (7.1) with $g(x)=x^{-\frac{d}{\alpha}-1+\frac{1}{2(1-\beta)}}$ and $h(x)=x^{\frac{1}{1-\beta}}$. Combining the estimates for $I_{1}$ and $I_{2}$ proves (3.15), which completes the proof. ∎

Proposition 3.4**.**

Under the assumptions of Theorem 3.3, assume additionally that $w\in C^{d+1+[\alpha]+l}(\mathbb{S}^{d-1})$. Then the following estimates hold,

• •

For $\Omega\leq 1$,

[TABLE]

for all $k\leq l$ and $i_{1},\cdots,i_{k}$.

• •

For $\Omega\geq 1$,

[TABLE]

for all $k\leq l$ and $i_{1},\cdots,i_{k}$.

Proof.

Using (2.4) and (3.13),

[TABLE]

where

[TABLE]

and

[TABLE]

Again we are in the situation where these integrals are the same as those found in (3.2) after changing $d\mapsto d+k$. Thus we have for $\Omega\geq 1$,

[TABLE]

and

[TABLE]

Combining the estimates for $I_{1}$ and $I_{2}$ gives us (3.19). For $\Omega\leq 1$ we have

[TABLE]

It only remains to check the estimate for $I_{1}$ when $\Omega\leq 1$,

[TABLE]

Thus we have

[TABLE]

Combining the estimates for $I_{1}$ and $I_{2}$ for $\Omega\leq 1$ gives us (3.18). ∎

4 Local Estimates

In the following two sections we look at two other families of spatial operators which extend the global estimates obtained in the previous sections. Firstly we consider a more general second order elliptic operator (not necessarily in divergence form), then we consider homogeneous pseudo-differential operators with variable coefficients. In both cases we provide local (i.e., small-time) two-sided estimates for the Green’s functions of the associated fractional evolution equations. The key point here is that for these spatial operators, we no longer have global (in time) estimates for the associated Green’s functions. Before going to the new estimates, we describe how one turns local estimates into global estimates. If for some Green’s functions $G_{0}(t,x,y),G_{1}(t,x,y)$, one has the local two-sided estimate for some constant $c>0$

[TABLE]

for some fixed $T>0$, then by taking convolutions and using the Chapman-Kolmogorov equations,

[TABLE]

Repeating this procedure $n$-times,

[TABLE]

By fixing $t$ and setting $\tau=nt$ (so that $\tau\approx n$ for large values of $n$ and $\tau$), we then get

[TABLE]

Applying the same procedure to the lower bound gives us the global two-sided estimate

[TABLE]

for all $(\tau,x,y)\in(0,\infty)\times\mathbb{R}^{d}\times\mathbb{R}^{d}$.

4.1 Time-Fractional Diffusion Equation: General Non-Degenerate

In Section 3.1 we derived global two-sided estimates for the Green’s function of fractional evolution equations involving a fractional derivative in time and a second order elliptic operator in divergence form as the spatial operator. The key point in that case is that Aronsons estimates provides two-sided Gaussian estimates that hold globally for all time $t>0$. In this section we consider the case that the spatial operator is any non-degenerate diffusion operator, which can generally be of the form

[TABLE]

Assuming that $a(x)$ is uniformly elliptic and continuously differentiable, $b(x)$ and $c(x)$ are continuous, and the uniform bound holds:

[TABLE]

Then the Green’s function associated with (4.2), satisfies the following local estimates:

[TABLE]

for $(t,x,y)\in(0,T)\times\mathbb{R}^{d}\times\mathbb{R}^{d}$ for some fixed $T>0$. We also have the following estimates for the spatial derivative of the Green’s function $G$,

[TABLE]

for $(t,x,y)\in(0,T)\times\mathbb{R}^{d}\times\mathbb{R}^{d}$. The main obstacle now is that the estimates for the Green’s function of (4.2) are only for small-time, thus a serious problem seems to arise when trying to insert the local estimate into the Pollard-Zolotarev formula, which involves integrating over all time $z\in(0,\infty)$. However we use the trick described in the previous section to make the local estimates global, in (4.1). To this end, the following two-sided estimate holds for $G(t,x,y)$ for all $(\tau,x,y)\in(0,\infty)\times\mathbb{R}^{d}\times\mathbb{R}^{d}$,

[TABLE]

for some constant $c$. In addition, for all $(\tau,x,y)\in(0,\infty)\times\mathbb{R}^{d}\times\mathbb{R}^{d}$,

[TABLE]

Alternatively we can split the estimates for the spatial derivative up into small-time and large-time - for $\tau\in(0,1)$,

[TABLE]

and for $\tau\in(1,\infty)$,

[TABLE]

Now we proceed to obtain estimates for the Green’s function of the fractional evolution

[TABLE]

where $L$ is defined as above in (4.2). The Green’s function for this fractional evolution equation is given by

[TABLE]

Again let $\Omega=|x-y|^{2}t^{-\beta}$. We have the following local estimates for the Green’s function $G^{(\beta)}$ above.

Theorem 4.1**.**

Assume that $a(\cdot)\in C^{1}(\mathbb{R}^{d})$ is uniformly elliptic and $b(\cdot),c(\cdot)\in C(\mathbb{R}^{d})$. Suppose also that,

[TABLE]

Then for a fixed $T>0$, there exists constants $C_{1},C_{2},C_{3}$ such that for $(t,x,y)\in(0,T]\times\mathbb{R}^{d}\times\mathbb{R}^{d}$ the Green’s function $G^{(\beta)}(t,x,y)$ defined by (4.8) satisfies the following estimates,

• •

For $\Omega\leq 1$,

[TABLE]

• •

For $\Omega\geq 1$,

[TABLE]

where $C_{1},C_{2}$ depends on $T,d,\beta$ and $C_{3}$ depends on $T$ and $\beta$.

Proof.

First splitting up to the stable density,

[TABLE]

Note that on using the estimate (4.3) in $I_{1}$, we have the same integral of the same name appearing in (3.1). Thus for $\Omega\leq 1$,

[TABLE]

In addition for $\Omega\geq 1$,

[TABLE]

Turning our attention to $I_{2}$, let us consider seperately the upper and lower bound.

Upper bound for $I_{2}$

First applying the upper bound from (4.5) to $G$,

[TABLE]

For $\Omega\leq 1$, we have

[TABLE]

for $t<T$ for some fixed $T>0$. Combining this with (4.1) gives (4.9).

For $\Omega\geq 1$, we use again that the decay of $\exp\{-c_{\beta}z^{1/(1-\beta)}\}$ for large $z$ is stronger than the growth of $\exp\{ct^{\beta}z\}$ for large $z$ as long as $t<T$ for some fixed $T>0$. That is,

[TABLE]

where we have made the substitution $w=z^{-1}$ in the last line. Now we apply Proposition 7.1 with $N=\frac{d}{2}-1-\frac{1}{2(1-\beta)}$ and $a=\frac{1}{1-\beta}$,

[TABLE]

**Note the constants in the above estimate depend on $T$. Combining this with (4.15) gives us the required upper bound in (4.10).

Lower bound for $I_{2}$

Using the lower bound from (4.5) in $I_{2}$,**

[TABLE]

Firstly for $\Omega\leq 1$,

[TABLE]

Finally for $\Omega\geq 1$,

[TABLE]

where we have used the fact that $\exp\{-ct^{\beta}z\}\geq\exp\{-ct^{\beta}z^{\frac{1}{1-\beta}}\}$ for $z>1$. After making the substitution $z=w^{-1}$ we apply Proposition 7.1,

[TABLE]

where $C_{1}$ depends on $T,\beta$ and $d$, and $C_{2}$ depends on $T$ and $\beta$. Combining this with (4.15) gives us the lower bound in (4.10), as required. ∎

Next we look at estimating the spatial derivative of the Green’s function $G^{(\beta)}$, firstly for large-time using (4.7) then for small-time using (4.6). As usual, let $\Omega:=|x-y|^{2}t^{-\beta}$. Firstly for large finite time,

Proposition 4.2**.**

Under the same assumptions as Theorem 4.1, suppose further that $a(x)$ is twice continuously differentiable, and $b(x)$, $c(x)$ are continuously differentiable (with all derivatives bounded). Then for a fixed finite $T>1$, the following estimates hold for the spatial derivative of the Green’s function $G^{(\beta)}(t,x,y)$ for $(t,x,y)\in(1,T)\times\mathbb{R}^{d}\times\mathbb{R}^{d}$,

• •

For $\Omega\leq 1$,

[TABLE]

• •

For $\Omega\geq 1$,

[TABLE]

Proof.

We start as usual by first splitting up the integral into small and large $z$, and also use the triangle inequality,

[TABLE]

Note that $t\in(1,T)$ means that $t^{-\beta}\in(T^{-\beta},1)$. Thus for $z\in(1,\infty)$ we have $z\geq t^{-\beta}$. Now we use the local estimate (4.6) for the first integral and (4.7) for the second,

[TABLE]

Note that the integral in $I_{1}$ is the same as (3.1), and thus for $\Omega\leq 1$

[TABLE]

Note however that $t\in(1,T)$, which means that $t^{-\beta}\in(T^{-\beta},1)$. Thus

[TABLE]

For $\Omega\geq 1$,

[TABLE]

As for the integral $I_{2}$, this is the same one which appeared in the previous proof, (4.1), and thus for $\Omega\leq 1$,

[TABLE]

Combining this with (4.19) which gives both (4.17). Finally an application of Proposition 7.1 gives for $\Omega\geq 1$,

[TABLE]

Combining this with the estimate for $I_{1}$, gives the estimate (4.18) for $\Omega\geq 1$, as required. ∎

Next we have the estimates for small-time.

Proposition 4.3**.**

Under the same assumptions as Theorem 4.1, suppose further that $a(x)$ is twice continuously differentiable, and $b(x)$, $c(x)$ are continuously differentiable (with all derivatives bounded). Then the following estimates hold for the spatial derivative of the Green’s function $G^{(\beta)}(t,x,y)$ for $(t,x,y)\in(0,1)\times\mathbb{R}^{d}\times\mathbb{R}^{d}$,

• •

For $\Omega\leq 1$,

[TABLE]

• •

For $1\leq\Omega\leq t^{-\beta\left(\frac{2-\beta}{1-\beta}\right)}$,

[TABLE]

• •

For $\Omega\geq t^{-\beta\left(\frac{2-\beta}{1-\beta}\right)}$,

[TABLE]

Proof.

Splitting the integral up using the stable density then using the estimates (4.6) and (4.7),

[TABLE]

Now we investigate the usual cases.

Case 1: $\Omega\leq 1$****

The integral in $I_{1}$, being the same as the one in (4.19), has the upper bound

[TABLE]

The other two integrals in $I_{2}$ and $I_{3}$ approach convergent integrals for bounded $\Omega$, so

[TABLE]

and

[TABLE]

Thus in this case,

[TABLE]

Case 2: $\Omega\geq 1$****

A direct application of the Laplace method gives

[TABLE]

For the second integral we have,

[TABLE]

where we have used Proposition 7.1 in the last estimate. Finally since $t\in(0,1)$, another application of Proposition 7.1 gives

[TABLE]

Note that

[TABLE]

when $t^{-\frac{\beta}{2}}\Omega^{-\frac{1}{2}\left(\frac{1-\beta}{2-\beta}\right)}\leq 1$. Thus for $\Omega\geq t^{-\beta\left(\frac{2-\beta}{1-\beta}\right)}$,

[TABLE]

while for $1\leq\Omega\leq t^{-\beta\left(\frac{2-\beta}{1-\beta}\right)}$,

[TABLE]

∎

4.2 Time-Fractional Pseudo-Differential Evolution: Variable Coefficients

Finally we derive two-sided estimates for the Green’s function of time-fractional stable-like equations. Stable-like operators are homogeneous pseudo-differential operators with variable coefficients (that depend on the spatial variable $x$, but not time). As noted earlier in (3.12) the fundamental solution $G_{\psi}$ of the evolution equation

[TABLE]

with $\psi_{\alpha}(p)=|p|^{\alpha}w(p/|p|)$, satisfies the following two-sided estimate for all $(t,x-y)\in(0,\infty)\times\mathbb{R}^{d}$,

[TABLE]

When the coefficients of the operator $\psi_{\alpha}$ depends also on the spatial variable, the same kind of estimates hold for small-time. Let $G_{\psi_{\alpha},x}$ denote the fundamental solution to the pseudo-differential evolution equation

[TABLE]

with homogeneous symbol $\psi_{\alpha}(x,p)=|p|^{\alpha}w_{\mu}(x,p/|p|)$, where

[TABLE]

Theorem 4.4**.**

Assume that $w_{\mu}(x,p)$ is a $\gamma$-Hölder continuous function in the variable $x$ taking values in a compact subset of $(0,\infty)$, $\gamma\in(0,1]$. Assume further that $\mu$ has a strictly positive density. Then for some fixed $T>0$, there exists a constant $C>0$ such that for $t\in(0,T)$ and $x,y\in\mathbb{R}^{d}$,

[TABLE]

What this means is that the global in time estimates (4.20) for the Green’s function $G_{\psi_{\alpha}}$, also serve as a small-time estimate for the Green’s function $G_{\psi_{\alpha},x}$. Indeed one would hope that operators with variable coefficients can be approximated by the method of freezing coefficients. So we have the following small-time estimate for $t\in(0,T)$, $x,y\in\mathbb{R}^{d}$

[TABLE]

for some fixed $0<T<\infty$. We also have the following estimates for the spatial derivatives of the $G_{\psi_{\alpha},x}$, see [kolokoltsov2019differential] (Theorem 5.8.3).

Theorem 4.5**.**

Let $\alpha>0$, and denote by $l$ the maximal integer less than $\alpha$. Assume that $\mu\geq\mu_{0}>0$, for some positive number $\mu_{0}$, and $w_{\mu}(x,p)$ is $q$-times differentiable in $x$ and each of these derivatives be $(d+1+(l+q)(\alpha+1))$-times continuously differentiable in $p$ and all bounds uniform in $x,p$. Then for a fixed $T>0$ and any $k\leq l$,

[TABLE]

for $(t,x,y)\in(0,T)\times\mathbb{R}^{d}\times\mathbb{R}^{d}$.

Using the same technique as the previous section to extend these small-time estimates to global estimates, we have the following two-sided estimates for $\tau>0$, $x,y\in\mathbb{R}^{d}$,

[TABLE]

and

[TABLE]

Now consider the following fractional evolution equation,

[TABLE]

with

[TABLE]

where $w_{\mu}$ satisfies the assumptions of Theorem 4.5. The solution of (4.25) is given by

[TABLE]

where $E_{\beta}$ is the Mittag-Leffler function. The Green’s function of Equation (4.25) is then

[TABLE]

Let $\Omega=|x-y|^{\alpha}t^{-\beta}$.

Theorem 4.6**.**

Let $\alpha\in(0,2)$ and $\beta\in(0,1)$. Assume that the function $w_{\mu}$ in (4.26) is $\gamma$-Hölder continuous in the first variable and $k$-times continuously differentiable in the second variable. Assume further that the spectral measure $\mu$ has a strictly positive density. Then for a fixed $T>0$ there exists constants $C$ such that for $(t,x,y)\in(0,T]\times\mathbb{R}^{d}\times\mathbb{R}^{d}$ the following two-sided estimates for (4.27) hold,

• •

For $\Omega\leq 1$,

[TABLE]

• •

For $\Omega\geq 1$,

[TABLE]

where the constants $C$ depend on $d,\alpha,\beta$ and $T$.

Proof.

We start by estimating the stable density with (2.4),

[TABLE]

In the first integral, we use the estimate (4.21), and for the second term we use the global version (4.23) with $\tau=t^{\beta}z$. Starting with the upper bound, we have

[TABLE]

Recall that,

[TABLE]

Then we have

[TABLE]

for the upper bound, and

[TABLE]

for the lower bound. Note that the integral in $I_{1}$ is the as the one appearing in (3.2) and so, for $t<T$,

[TABLE]

for $\Omega\leq 1$, and

[TABLE]

for $\Omega\geq 1$. For the remaining integral $I_{2}$, we have the usual two cases.

Case 1: $\Omega\leq 1$.** In this case we have**

[TABLE]

This integral converges as long as $t<T$, since $\exp\{ct^{\beta}z\}\leq\exp\{Cz^{\frac{1}{1-\beta}}\}$ for sufficiently large $z$. Thus

[TABLE]

On the other hand, we have

[TABLE]

which is strictly positive for $t<T$, thus

[TABLE]

Combining these estimates with those for $I_{1}$ gives the estimates for $G^{(\beta)}_{\psi_{\alpha},x}$ for $\Omega\leq 1$.

Case 2: $\Omega\geq 1$.** In this case we have**

[TABLE]

and

[TABLE]

Firstly we have,

[TABLE]

and

[TABLE]

Next note that $\exp\{t^{\beta}z\}\leq\exp\{T^{\beta}z\}$ and $\exp\{-t^{\beta}z\}\geq\exp\{-T^{\beta}z\}$ for $t<T$. Thus,

[TABLE]

and

[TABLE]

Thus for $\Omega\geq 1$, we have

[TABLE]

as claimed.∎

Next we look at the spatial derivatives, where we consider separately small and large (but finite) time.

Proposition 4.7**.**

Under the same assumptions as Theorem 4.6 and Theorem 4.5, the spatial derivatives of the Green’s function $G^{(\beta)}_{\psi_{\alpha},x}(t,x,y)$ for $(t,x,y)\in(0,1)\times\mathbb{R}^{d}\times\mathbb{R}^{d}$ satisfy,

• •

For $\Omega\leq 1$,

[TABLE]

for all $k\leq l$ and all indicies $i_{i},\cdots,i_{k}$.

• •

For $1\leq\Omega\leq t^{-\beta}$,

[TABLE]

for all $k\leq l$ and all indicies $i_{i},\cdots,i_{k}$.

• •

For $\Omega\geq t^{-\beta}$,

[TABLE]

for all $k\leq l$ and all indicies $i_{i},\cdots,i_{k}$.

Proof.

Splitting up the stable density followed by using the estimates (4.22) and (4.24) we have

[TABLE]

where

[TABLE]

and

[TABLE]

Now note that since $t\in(0,1)$, the integral in $I_{1}$ is the same as that one appearing in (3.20), and thus for $\Omega\leq 1$,

[TABLE]

For $\Omega\geq 1$,

[TABLE]

Turning to $I_{2}$, we need to consider some different cases.

Case 1: $\Omega\leq 1$.** In this case**

[TABLE]

Combining this with the estimate for $I_{1}$ shows (4.30).

Case 2: $1\leq\Omega\leq t^{-\beta}$.** In this case we have**

[TABLE]

Combining this with (4.33) shows (4.31).

Case 3: $t^{-\beta}\leq\Omega$.****

[TABLE]

Finally combining this with (4.33) shows (4.32). ∎

Next, for large (finite) time.

Proposition 4.8**.**

Under the same assumptions as Theorem 4.6 and Theorem 4.5, then for fixed $T>0$, the following estimates hold for the spatial derivatives of the Green’s function $G^{(\beta)}_{\psi_{\alpha},x}(t,x,y)$ for $(t,x,y)\in(1,T)\times\mathbb{R}^{d}\times\mathbb{R}^{d}$,

• •

For $\Omega\leq 1$,

[TABLE]

for all $k\leq l$ and all indicies $i_{i},\cdots,i_{k}$.

• •

For $\Omega\geq 1$,

[TABLE]

for all $k\leq l$ and all indicies $i_{i},\cdots,i_{k}$.

Proof.

As usual we first use 2.5,

[TABLE]

Next, we use the estimate (4.22) for the first term and (4.24) for the second. Note that since $t\in(1,T)$, then

[TABLE]

The integral in $I_{1}$ is the same as that one appearing in (3.20), and thus for $\Omega\leq 1$,

[TABLE]

However in this situation $t\in(1,T)$, so $t$ is away from both [math] and $\infty$. Thus, recalling that $\Omega=|x-y|^{\alpha}t^{-\beta}$,

[TABLE]

For $\Omega\geq 1$,

[TABLE]

Furthermore, the integral $I_{2}$ is the same as the one defined as $I^{up}_{2}$ in (4.2), thus for $\Omega\leq 1$,

[TABLE]

So for $\Omega\leq 1$,

[TABLE]

which thus gives (4.34). Finally for $\Omega\geq 1$, using (4.2),

[TABLE]

thus combining the estimates for $I_{1}$ and $I_{2}$ for $\Omega\geq 1$ gives us (4.35). ∎

5 Generalised Evolution Equations

In this last section, we look at the following generalised evolution,

[TABLE]

where $D^{(\nu)}_{0+*}$ is the Caputo-type operator

[TABLE]

Here $\nu(t,\cdot)$ is a Lévy transition kernel that satisfies $\sup_{t}\int\min(1,r)\nu(t,\mathrm{d}r)<\infty$. The solution to Equation (5.1) is given by

[TABLE]

where $E_{(\nu),t}(A)$ is the operator-valued generalised Mittag-Leffler function which is defined by the operator-valued integral

[TABLE]

where $\Pi^{-A}_{(\nu)}$ is the operator-valued potential measure of the semigroup $T^{(\nu)}_{t}e^{tA}$ generated by $(-D^{(\nu)}_{0+*}-A)$,

[TABLE]

Then we can rewrite this solution to get the Green’s function,

[TABLE]

where $G^{(\nu)}_{A}$ is the Green’s function of the evolution Equation (5.1) given by

[TABLE]

We will use the following comparison principle from [kolokoltsov2019mixed].

Theorem 5.1**.**

Let $\nu$ and $\tilde{\nu}$ be two Lévy measures satisfying

[TABLE]

[TABLE]

and $\nu(t,(0,\infty))=\tilde{\nu}((0,\infty))=\infty$. Then for any non-increasing function $f$ we have the comparison principle for the semigroups:

[TABLE]

where $T^{(\nu)}_{t}$ and $T^{\tilde{\nu}}_{t}$ are the semigroups generated by

[TABLE]

and

[TABLE]

respectively. Moreover, the potential measures of the semigroups $T_{t}^{\nu}$ and $T_{t}^{\tilde{\nu}}$ satisfy the comparison principle,

[TABLE]

A direct application of this comparison principle gives us the main result of this article.

Theorem 5.2**.**

Let $A$ be one of the spatial operators from (3.1), Theorem 3.3, (4.2) or (4.25) along with their relevant assumptions. Let $\nu(t,\mathrm{d}s)$ be a Lévy transition kernel which has upper and lower bounds of $\beta$-fractional type,

[TABLE]

for some $\beta_{1},\beta_{2}\in(0,1)$ and $C_{\nu}>0$. Then

[TABLE]

for a non-increasing function $\phi$, where

[TABLE]

and

[TABLE]

Proof.

This follows from the formula (5.2) and an application of the comparison principle for potential operators. ∎

Thus the estimates obtained in Theorem 3.1, Theorem 3.3, Theorem 4.1 and Theorem 4.6 can be used to estimate solutions of generalised evolutions (5.1).

Remark 1**.**

In order to see why this result is expected, let us give some intuition behind the comparison principle. The assumption that the Lévy kernel is bounded below by the Lévy kernel of a $\beta$-stable subordinator, means that the Lévy subordinator generated by the operator

[TABLE]

has on average jumps that are larger than those of the process generated by

[TABLE]

So on the sample paths level, the jumps of $X^{(\nu)}$ will typically be larger than those of $X^{\beta}$, which means that the inverse process of $X^{(\nu)}$ will typically be constant for longer times than the inverse process of $X^{\beta}$. Thus when we subordinate the spatial process, $Y(t)$, generated by the operator $A$ by the inverse subordinator given by

[TABLE]

and compare its paths to the spatial process subordinated by an inverse stable subordinator $S^{\beta}_{t}$, we will see that $Y(S_{t}^{\nu})$ is dominated by $Y(S_{t}^{\beta})$ in the sense that $Y(S_{t}^{\nu})$ will have longer trapping times.

6 Conclusion

In this article, we have looked at two-sided estimates for the Green’s function of fractional evolution equations of the form

[TABLE]

The solution of such fractional evolution equations can be written with the help of operator-valued Mittag-Leffler functions,

[TABLE]

We have given two-sided estimates for the Green’s function $G^{(\beta)}_{L}(t,x,y)$ (and its spatial derivatives) in several different situations. The situations can be split up into two broad cases: when the Green’s function $G_{L}$ associated with $L$ does or does not have known global in time estimates. In those two cases, we consider generators of diffusion processes in Theorems 3.1 and 4.1; and we consider generators of stable and stable-like processes in Theorems 3.3 and 4.6. Finally, we looked at generalised evolution equations where the operator acting on the time variable is given by a Caputo-type operator

[TABLE]

We concluded that solutions to generalised evolution equations of the form

[TABLE]

where $\nu(t,\mathrm{d}s)$ is a Lévy-type kernel which for fixed $t$ is comparable to the Lévy measure of a $\beta$-stable subordinator, could be estimated using the estimates obtained for $G^{(\beta)}_{L}$. Then whenever one is looking at evolution equations of the form (6.1), or, from the probabilistic point of view, at stochastic processes generated by $-D^{(\nu)}-L$, then under the assumption that $\nu$ is comparable to $\beta$-stable, the estimates shown in this article can be used to gain a lot of information.

Note that in this article we have viewed $G^{(\beta)}_{L}(t,x,y)$ as the Green’s function of the evolution equation

[TABLE]

Probabilistically speaking, $G_{L}^{(\beta)}$ are the transition densities of the process $X_{t}^{L,\beta}$ generated by $-D^{\beta}-L$. The process $X_{t}^{L,\beta}$ is the subordination of the process generated by $L$ by the inverse of the process generated by $D^{\beta}$. In this view one could use the estimates in this article to obtain sample path properties of a subordinated process $X^{L,\beta}_{t}$.

7 Asymptotic Methods

The main idea of the Laplace method for estimating integrals of the form

[TABLE]

is that the major contribution to the asymptotic behaviour comes from a neighbourhood around the point at which the function $h(x)$ in the exponent attains its minimum value. Outside this neighbourhood the contribution is exponentially small, and so when one proves asymptotic formulas using Laplace methods, the integrals are split up into the neighbourhood around which the major contribution occurs (or around each such neighbourhood, if $-h(x)$ is not unimodal) and the regions for which the approximation error is exponentially small. Although we focus on integrals over some interval $(a,\infty)$, the point is that extending the interval only introduces exponentially small errors and so the value of the integral over a larger interval essentially the same. Our standard references for asymptotic methods are [de1970asymptotic], [murray2012asymptotic] or [fedorjuk1987asymptotics]. Consider the integral

[TABLE]

Let us assume that $h$ is a real continuous function, and that it attains its minimum at the boundary point $b$, that $h^{\prime}(b)$ exists and $h^{\prime}(b)>0$. Moreover assume that $h(x)>h(b)$ (for $x>b$) and $h(x)\rightarrow\infty$ as $x\rightarrow\infty$. We will not recount the proof, but we state the asymptotic formula,

[TABLE]

On the other hand, if the function $h$ has a minimum on the interior of the interval $(b,\infty)$, say at the point $\tilde{b}\in(b,\infty)$. Finally, assume that the derivative $h^{\prime}(x)$ exists in some neighbourhood of $x=\tilde{b}$, that $h^{\prime\prime}(\tilde{b})$ exists and that $h^{\prime\prime}(\tilde{b})>0$. Then

[TABLE]

Proposition 7.1**.**

Let $a>0$, $N\in\mathbb{R}$, $c>0$ and $\Omega\geq 1$. Then the following asymptotic formula holds as $\Omega\rightarrow\infty$,

[TABLE]

where $C_{1}(a,N,c)=(ac)^{\frac{2(N+1)-1}{2(a+1)}}\sqrt{\frac{2\pi}{a+1}}$, and $C_{2}(c,a)=(ac)^{\frac{1}{a+1}}[1+a^{-1}]$.

Proof.

Define

[TABLE]

and let $h(w)=-w\Omega-cw^{-a}$. Differentiating $h$ with respect to $w$, one finds the maximum of $h$ at

[TABLE]

Now the trick is to make the substitution $w=w_{*}s$ in the integral $J(\Omega)$, to obtain

[TABLE]

Now we are in a position to apply the asymptotic formula (7.2), with $g(s)=s^{N}$, $h(s)=s+a^{-1}s^{-a}$ and $\lambda=(\Omega^{a}ac)^{\frac{1}{a+1}}$. For this we need some derivatives of $h$,

[TABLE]

[TABLE]

thus $h$ has a minimum at $s=1\in(0,w_{*}^{-1})$. Finally applying (7.2) we have

[TABLE]

as required. ∎

Note that when applying the above asymptotic formula, we will not care so much what the constants are, only what they depend on.