Mixed linear fractional boundary value problems

TL;DR

This paper derives two-sided estimates for the Green's function of fractional boundary value problems involving mixed Caputo derivatives and diffusion operators, extending to multiple fractional derivatives.

Contribution

It provides explicit estimates for Green's functions in mixed fractional boundary value problems, including cases with multiple fractional derivatives, advancing analytical understanding.

Findings

Two-sided estimates for Green's functions are established.

The approach extends to problems with multiple fractional derivatives.

Results facilitate analysis of fractional PDEs with boundary conditions.

Abstract

In this article we obtain two-sided estimates for the Greens function of fractional boundary value problems on $\mathbb R_+ \times \mathbb R_+ \times \mathbb R^d$ of the form \[(-{}_{t_1}D^\beta_{0+*} - {}_{t_2}D^\gamma_{0+*})u(t_1, t_2, x) = L_{x}u(t_1, t_2, x),\] with some prescribed boundary functions on the boundaries $\{0\} \times \mathbb R_+ \times \mathbb R^d$ and $\mathbb R_+ \times\{0\}\times \mathbb R^d$. The operators ${}_{t_1}D^\beta$ and ${}_{t_1}D^\gamma$ are Caputo fractional derivatives of order $\beta, \gamma \in (0, 1)$ and $L_{x}$ is the generator of a diffusion semigroup: $L_x= \nabla \cdot(a(x) \nabla)$ for some nice function $a(x)$. The Greens function of such boundary value problems are decomposed into its components along each boundary, giving rise to a natural extension to the case involving $k \geq 2$ number of fractional derivatives on the left hand side.

Equations

$$(-\text{prescript}{t}_1{D}_{0+\ast }^{\beta }-\text{prescript}{t}_2{D}_{0+\ast }^{\gamma })u(t_1,t_2,x)=L_{x}u(t_1,t_2,x),$$

$$-\text{prescript}t{D}_{0+\ast }^{\beta }u(t,x)=L_{x}u(t,x),$$

$$-\text{prescript}{t}_1{D}_{0+\ast }^{\beta }u(t_1,t_2,x)-\text{prescript}{t}_2{D}_{0+\ast }^{\gamma }u(t_1,t_2,x)$$

$$u(0,t_2,x)$$

$$u(t_1,0,x)$$

$$-\sum _{i=1}^k\text{prescript}{t}_i{D}_{0+\ast }^{{\beta }_i}u(t,x)=L_{x}u(t,x),$$

$$-\sum _{i=1}^k\text{prescript}{t}_i{D}_{0+\ast }^{{\nu }_i(t_i,\cdot )}u(t,x)=L_{x}u(t,x),$$

$${\tau }_0^{{\beta }_i}:=\inf \{s>0:X_{t_i}^{{\beta }_i}(s)\le 0\},$$

$$\Psi (t_i,T)=\mathbb{P}[{\tau }_0^{{\beta }_i}<T].$$

$$C_{constA}(S)=\{f\in C(S):{f|}_A\text{ is a constant}\}.$$

$$I_a^{\alpha }f(t)=I_{a+}^{\alpha }f(t)=\frac{1}{\Gamma (\beta )}{\int }_a^t(t-r{)}^{\beta -1}f(r)\mathrm{d}r.$$

$$D_{a+\ast }^{\beta }f(t)=I_a^{1-\beta }\left[\frac{d}{dt}f\right](t)=\frac{1}{\Gamma (1-\beta )}{\int }_a^t(t-r{)}^{-\beta }\left[\frac{d}{dr}f\right](r)\mathrm{d}r,t>a.$$

$$D_{a+\ast }^{\beta }f(t)=\frac{1}{\Gamma (-\beta )}{\int }_0^{t-a}\frac{f(t-r)-f(t)}{r^{1+\beta }}\mathrm{d}r+\frac{f(t)-f(a)}{\Gamma (1-\beta )}(t-a{)}^{\beta }.$$

$$\frac{d^{\beta }}{d{x}^{\beta }}f(x):=D_{-\infty +\ast }^{\beta }f(x)=\frac{1}{\Gamma (-\beta )}{\int }_0^{\infty }\frac{f(x-z)-f(x)}{z^{1+\beta }}\mathrm{d}z.$$

$${\tau }_0^{\alpha }=\inf \{s>0:X_t^{\alpha }(s)\le 0\}.$$

$${\mu }_0^{\alpha }(s)=\frac{\partial }{\partial s}\mathbb{P}\left[{\tau }_0^{\alpha }>s\right]$$

$$=\frac{t}{\alpha }{s}^{-1-\frac{1}{\alpha }}{w}_{\alpha }(t{s}^{-1/\alpha }).$$

$$t^{\alpha }{\mu }_0^{\alpha }(t^{\alpha }s)=\frac{1}{\alpha }{s}^{-1-\frac{1}{\alpha }}{w}_{\alpha }(s^{-\frac{1}{\alpha }}).$$

$$w_{\alpha }(r)\sim C_{\alpha }{r}^{-\frac{2-\alpha }{2(1-\alpha )}}\exp \{-c_{\alpha }{r}^{-\frac{\alpha }{1-\alpha }}\}:=f_{\alpha }(r),r\to 0,$$

$$w_{\alpha }(r)\sim \widetilde{C_{\alpha }}{r}^{-1-\alpha },r\to \infty .$$

$$\exp \{-{\lambda }^{\alpha }\}=\mathbb{E}[\exp \{-\lambda X^{\alpha }(1)\}]={\int }_0^{\infty }{e}^{-\lambda x}{w}_{\alpha }(x)\mathrm{d}x,$$

$$w_{\alpha }(r)\asymp C{1}_{\{r<1\}}{f}_{\alpha }(r)+\tilde{C}{1}_{\{r\ge 1\}}{r}^{-1-\alpha }.$$

$$t^{\alpha}\mu_{0}^{\alpha}(t^{\alpha}s)\sim\left\{\begin{array}[]{lc}c_{\alpha},&s\rightarrow 0,\\ c_{\alpha}s^{-1+\frac{1}{2(1-\alpha)}}\exp\{-cs^{\frac{1}{1-\alpha}}\},&s\rightarrow\infty.\end{array}\right.$$

$$t^{\alpha }{\mu }_0^{\alpha }(t^{\alpha }s)=c_{\alpha }{s}^{-1-\frac{1}{\alpha }}{w}_{\alpha }(s^{-\frac{1}{\alpha }})\sim c_{\alpha }{s}^{-1-\frac{1}{\alpha }}{s}^{1+\frac{1}{\alpha }}=c_{\alpha }.$$

$$t^{\alpha }{\mu }_0^{\alpha }(t^{\alpha }s)=c_{\alpha }{s}^{-1-\frac{1}{\alpha }}{w}_{\alpha }(s^{-\frac{1}{\alpha }})$$

$$=c_{\alpha }{s}^{-1+\frac{1}{2(1-\alpha )}}\exp \left\{-c_{\alpha }{s}^{\frac{1}{1-\alpha }}\right\},$$

$$G^Y(s,x,y)\asymp s^{-\frac{d}{2}}\exp \left\{-c\frac{|x-y{|}^2}{s}\right\}.$$

$$G^{Y,\gamma }(s,r,x,y)=G^Y(s,x,y)s^{-\frac{1}{\gamma }}{w}_{\gamma }(r{s}^{-\frac{1}{\gamma }}).$$

$$G^{Y,\gamma }(s,r,x,y)\asymp C{r}^{-1-\gamma }{s}^{1-\frac{d}{2}}\exp \left\{-c\frac{|x-y{|}^2}{s}\right\}.$$

$$G^{Y,\gamma }(s,r,x,y)\asymp C{r}^{-\frac{2-\gamma }{2(1-\gamma )}}{s}^{\frac{1}{2(1-\gamma )}-\frac{d}{2}}\exp \left\{-c\frac{|x-y{|}^2}{s}-c{s}^{\frac{1}{1-\gamma }}{r}^{-\frac{\gamma }{1-\gamma }}\right\}$$

$${\tau }_0^{\beta ,\gamma }=\min \left({\tau }_0^{\beta },{\tau }_0^{\gamma }\right).$$

$$\begin{array}[]{rcl}(-\prescript{}{t_{1}}{D}^{\beta}_{0+*}-\prescript{}{t_{2}}{D}^{\gamma}_{0+*})f(t_{1},t_{2},x)&=&Lf(t_{1},t_{2},x),\\ f(0,t_{2},x)&=&\phi_{1}(t_{2},x),\\ f(t_{1},0,x)&=&\phi_{2}(t_{1},x).\\ \end{array}$$

Full text

[TABLE]

with some prescribed boundary functions on the boundaries $\{0\}\times\mathbb{R}_{+}\times\mathbb{R}^{d}$ and $\mathbb{R}_{+}\times\{0\}\times\mathbb{R}^{d}$. The operators $\prescript{}{t_{1}}{D}^{\beta}$ and $\prescript{}{t_{1}}{D}^{\gamma}$ are Caputo fractional derivatives of order $\beta,\gamma\in(0,1)$ and $L_{x}$ is the generator of a diffusion semigroup: $L_{x}=\nabla\cdot(a(x)\nabla)$ for some nice function $a(x)$. The Greens function of such boundary value problems are decomposed into its components along each boundary, giving rise to a natural extension to the case involving $k\geq 2$ number of fractional derivatives on the left hand side.

1 Introduction

In the recent article [johnston2019green], we obtained two-sided estimates for the fundamental solution of fractional evolution equations of the form

[TABLE]

where $\prescript{}{t}{D}^{\beta}_{0+*}$ is a Caputo-Dzherbashyan (CD) fractional derivative of order $\beta\in(0,1)$ and $L_{x}$ is either the generator of a diffusion process, or a stable-like process (i.e, either a second order uniformly elliptic operator, or a spatially homogeneous pseudo-differential operator with variable coefficients). In this article we obtain two-sided estimates for the Greens function of the following boundary value problem,

[TABLE]

See Section 3 for details on what is the object that we estimate. In Section 4 we look at higher dimensional version of (1.1) in the sense that we have $k$ fractional derivatives on the left hand side, each acting on a different variable,

[TABLE]

where $(t,x)\in\mathbb{R}^{k}_{+}\times\mathbb{R}^{d}$, with some specified boundary behaviour.

The estimates obtained in this article can be used to study more general CD-type evolution equations (see [kolokoltsov2019differential, Section 8.5]) of the form

[TABLE]

where each $\nu_{i}(t_{i},\cdot)$ is a Lévy-type kernel, under the assumption that each $\nu_{i}(t_{i},\cdot)$ has a density which is comparable to the density of a $\beta_{i}$-stable process. This was done for the case $k=1$ in [johnston2019green], so we do not repeat it here. Another natural (and essentially straightforward) extension of the estimates obtained in this article would be to replace $L_{x}$ with a homogeneous pseudo-differential operator with variable coefficients which generates a stable-like (Feller) process.

Boundary value problems such as (1.1) and (1.2) can be found in many areas of mathematics. A particularly noteworthy application can be found in the mathematics of insurance. Consider $k$ processes $(X_{t_{1}}^{\beta_{1}}(s),\cdots,X_{t_{k}}^{\beta_{k}}(s))$, where each $X^{\beta_{i}}_{t_{i}}(s)$ is generated by $-\prescript{}{t_{i}}{D}^{\beta_{i}}_{0+*}$. If each process corresponds to the wealth of a company, then whenever one of the coordinates hit zero, one of the companies have defaulted. Insurance companies are interested the ruin probability, which is the probability of one of companies defaulting before a finite time horizon $T$. That is, if $\tau_{0}^{\beta_{i}}$ denotes the first time the process $X^{\beta_{i}}_{t_{i}}(s)$ hits zero,

[TABLE]

then the ruin probability is the quantity

[TABLE]

See [djehiche1993ruin, MR3520310, MR3515875, MR3007885, MR3343416, kumar2018fractional] for ruin probabilities of multidimensional risk models, or [Asmussen2010ruin] for a broader treatment of ruin probabilities. Fractional version of compound Poisson processes are also of interest when looking at insurance risk processes, see [leonenko2019limit]. Similar kinds of questions also appear when looking at barrier options under one-dimensional Markov models, see [MR3015232]. It is natural to consider multi-dimensional versions of these, [MR3515877], as investors often deal with basket options. A further natural appearance comes when considering portfolios of credit derivative obligations (CDO), which can be described by a Markov process in $\mathbb{R}^{k}_{+}$. Reaching a boundary of dimension $k-n$ means that $n$ out of $d$ bonds underlying the portfolio of CDOs have defaulted. It is natural in this setting to consider spatially non-homogeneous processes, since the behaviour of the processes should feel the approach to the boundary, which is not the case for Lévy processes. This is then the setting of problem (1.3). The series of articles [scalas2000fractional, mainardi2000fractional, gorenflo2001fractional], give a nice overview of the usage of fractional calculus and jump-diffusion processes in finance.

Another popular model these days is the so called Pearson diffusion, and also the fractional version, which are diffusion processes with polynomial diffusion coefficients, see [leonenko2013correlation]. Fractional models are also finding new footing in theoretical physics, via fractional and non-local Schrödinger operators, see for example the articles [kaleta2019zero, kaleta2018contractivity].

Of more general interest in finance are affine processes which live in $\mathbb{R}^{k}_{+}\times\mathbb{R}^{d}$, see [duffie2003affine]. Our final motivation for considering stable processes on $\mathbb{R}^{2}_{+}$ (i.e, (1.1) without the spatial operator $L_{x}$), is the topic of limit order books. A simplified model would be that one coordinate of $X^{\beta_{1},\beta_{2}}_{t_{1},t_{2}}(s)$ is the volume of trades available at the best buy price while the other is the volume at the best sell price. The event that this process hits the boundary means that there are no more trades offered at that price and thus a price change occurs. These ideas will be developed in a forthcoming paper. See [cont2012order, hambly2018limit] and references therein for related attempts at modelling order books using Brownian motions on the orthant and reflected SPDEs.

2 Preliminaries

2.1 Fractional derivatives and stable processes

We begin by fixing some definitions and notations. For an open or closed convex subset $S$ of $\mathbb{R}^{d}$, $C(S)$ is the Banach space of continuous functions on $S$ equipped with the sup-norm, $C_{\infty}(S)$ is the closed subspace of $C(S)$ consisting of functions vanishing at infinity. $C^{k}(S)$ is a Banach space of $k$ times continuously differential functions with bounded derivatives on $S$ with the norm being the sum of the sup norms of the function itself and all its partial derivatives up to and including order $k$.

For a subset $A\subset S$, let

[TABLE]

The (left) Riemann-Liouville (RL) integral of order $\alpha>0$ is defined as

[TABLE]

Then the fractional Caputo-Dzherbashyan (CD) of order $\beta\in(0,1)$ is defined as

[TABLE]

After some straightforward calculations, the CD derivative can be rewritten for smooth enough $f$ as

[TABLE]

For $a=-\infty$ (and for smooth bounded integrable functions), the operator $D^{\beta}_{-\infty+*}$ is known as the fractional derivative in generator form,

[TABLE]

The operators $-d^{\beta}f/dx^{\beta}$ with $\beta\in(0,1)$ generate stable Lévy subordinators (with inverted direction so that they are decreasing instead of increasing), see [meerschaert2012stochastic, Chapter 3]. Note that by (2.1), the CD derivative $D^{\beta}_{a+*}$ is obtained from $D^{\beta}_{-\infty+*}$ by the restriction of its action to the space $C_{const(-\infty,a]}(\mathbb{R})$. Probabilistically, this means that for $\beta\in(0,1)$ the process $\left(X^{\beta}_{x}(s)\right)_{s\geq 0}$ generated by $-D^{\beta}_{a+*}$ is a decreasing $\beta$-stable processes which is absorbed at $a$ on an attempt to cross it. We denote by $p^{\alpha}_{s}(t,r)$ the transition densities of the process $X^{\alpha}_{t}(s)$, and by $\tau_{0}^{\alpha}$ the time that $X_{t}^{\alpha}$ exits $(0,t]$, that is,

[TABLE]

Let $\mu_{0}^{\alpha}(s)$ be the density of the r.v $\tau_{0}^{\alpha}$.

Remark 1**.**

The density $\mu_{0}^{\alpha}$ exists and is continuous due to classical results from the theory of stable processes.

The density $\mu_{0}^{\alpha}(s)$ is given by,

[TABLE]

see [meerschaert2004limit, Corollary 3.1], where $w_{\alpha}(z)$ is the density of an $\alpha$-stable random variable. We can also write this density conveniently as

[TABLE]

The densities $w_{\alpha}(r)$ are one of the most important tools in studying the Greens function of evolution equations which involve fractional derivatives (of order at most 1). For $\alpha\in(0,1)$ the density $w_{\alpha}(r)$ has the following asymptotic behaviour in a neighbourhood of [math], [zolotarev1999chance, Theorem 5.4.1]

[TABLE]

and in a neighbourhood of $\infty$,

[TABLE]

Remark 2**.**

One can show (2.4) by inverting the Laplace transform

[TABLE]

and applying the method known as the saddle point method.

Due to the positivity of $w_{\alpha}(x)$, we can combine these behaviours so that there exists constants $C,\tilde{C}>0$ such that

[TABLE]

We will also be using the asymptotic behaviour of $\mu_{\alpha}^{0,t}(s)$, which follows from the asymptotic behaviour of $w_{\alpha}(s)$.

Lemma 2.1**.**

For $\alpha\in(0,1)$ the density $\mu_{0}^{\alpha}(s)$ of $\tau_{0}^{\alpha}$ has the following asymptotic behaviour at [math] and $\infty$,

[TABLE]

for some constants $c_{\alpha}>0$.

Proof.

Since $w_{\alpha}(r)\sim r^{-1-\alpha}$ as $r\rightarrow\infty$, then $w_{\alpha}(r^{-\frac{1}{\alpha}})\sim r^{1+\frac{1}{\alpha}}$ as $r\rightarrow 0$. Thus using (2.3), we have for $s\rightarrow 0$,

[TABLE]

Using (2.4), note that $w_{\alpha}(r^{-\frac{1}{\alpha}})\sim f_{\alpha}(r^{-\frac{1}{\alpha}})$ for $r\rightarrow\infty$. Thus for $s\rightarrow\infty$,

[TABLE]

as claimed. ∎

Let $Y_{x}(s)$ be a diffusion process with generator $L=\nabla\cdot(a(x)\nabla)$ for some symmetric measurable function $a$ on $\mathbb{R}^{d}$. The estimates of Aronson, [aronson1967bounds], say that the transition densities $G^{Y}(s,x,y)$ of $Y_{x}(s)$ satisfy the following two-sided Gaussian estimates for all $s>0$,

[TABLE]

Let $X^{\alpha}_{r}(s)$ be the process (independent of $Y_{x}(s)$) generated by $-D^{\alpha}_{0+*}$ (cf. 2.1), which is a decreasing $\beta$-stable process absorbed at [math] on an attempt to cross it. The transition density of the process $(Y_{x}(s),X^{\alpha}_{r}(s))$ is given by

[TABLE]

The following result is obtained by applying Aronsons estimate for $G^{Y}$ and (2.5) for $w_{\gamma}$.

Lemma 2.2**.**

The transition density of $(X^{\gamma}_{r}(s),Y_{x}(s))$ satisfy the following estimates

• •

For $s\leq r^{\gamma}$,

[TABLE]

• •

For $s>r^{\gamma}$,

[TABLE]

2.2 Processes on the orthant

Consider the process living on $\mathbb{R}^{2}_{+}$ defined by $X^{\beta,\gamma}_{t_{1},t_{2}}(s):=(X_{t_{1}}^{\beta}(s),X_{t_{2}}^{\gamma}(s))$, where each coordinate a one-dimensional stable subordinator (with inverted sign) which absorbed at [math], as described in the previous subsection. The process $X^{\beta,\gamma}_{t_{1},t_{2}}(s)$ is generated by $-\prescript{}{t_{1}}{D}^{\beta}_{0+*}~{}-~{}\prescript{}{t_{2}}{D}_{0+*}^{\gamma}$, where $\beta,\gamma\in(0,1)$, and it is started at $(t_{1},t_{2})\in\mathbb{R}_{+}\times\mathbb{R}_{+}$. For clarity, see Figure 1 for a typical sample path of $X^{\beta,\gamma}_{t_{1},t_{2}}(s)$. We assume that the processes $X^{\beta}$ and $X^{\gamma}$ are independent. This independence assumption implies that the first time the process $X^{\beta,\gamma}_{t_{1},t_{2}}$ hits the boundary of $\mathbb{R}_{+}\times\mathbb{R}_{+}$ is given by

[TABLE]

3 Mixed linear evolution

Consider the problem

[TABLE]

Here $L$ is the generator of a Feller process $Y_{x}(s)$ started at $x\in\mathbb{R}^{d}$. For simplicity we take $L=\nabla\cdot(a(x)\nabla)$, where $a(x)$ is a symmetric, uniformly elliptic and measurable function. This means that $A$ generates a non-degenerate diffusion, with transition densities $G^{Y}(s,x,y)$ which satisfy Aronsons two-sided estimates (2.6).

Remark 3**.**

Note that we could also obtain estimates for the Greens function in the case when $L$ is, say, a non-isotropic homogeneous pseudo-differential operator of order $\alpha\in(0,2)$ whose symbol is of the form

[TABLE]

where $w(x,\cdot)$ is some strictly positive function on $\mathbb{S}^{d-1}$. See [eidelman2004analytic, kolokoltsov2000symmetric] for the relevant estimates for $G^{Y}$ in that case.

3.1 Well-posedness of the mixed boundary value problem

Let us briefly discuss the well-posedness of problem (3.1). We only sketch the main steps, but see [kolokoltsov2019differential, Chapter 8], [hernandez2017generalised, Theorem 4.20] or [kolokoltsov2019mixed, Section 4] for a full account of well-posedness for these types of problems.

For even more general operators $A$ generating Feller semigroups (and even generalised versions of Caputo-derivatives), one can obtain both uniqueness and the stochastic representation (3.3) of the solution to (3.1) via the Dynkin formula [dynkin1965markov, Theorem 5.1]. To obtain existence of a classical solution, the main idea is to first transform the problem to an equivalent one involving zero boundary conditions and Riemann-Liouville fractional derivatives (by introducing a new unknown function $u(t_{1},t_{2},x)=f(t_{1},t_{2},x)-\mathbf{1}_{\{t_{2}>0\}}\phi_{1}(t_{2},x)-\mathbf{1}_{\{t_{1}>0\}}\phi_{2}(t_{1},x)$). This equivalent problem is then the following RL-type mixed boundary value problem,

[TABLE]

where

[TABLE]

Notice that here we require $\phi_{1}$ and $\phi_{2}$ to be in the domain of the generators $(-\prescript{}{t_{2}}{D}^{\gamma}_{0+*}-A)$ and $(-\prescript{}{t_{1}}{D}^{\beta}_{0+*}-A)$ respectively. The unique solution in the domain of the generator to (3.2) is then found by applying the potential operator (of the semigroup $T_{s}^{\beta}T_{s}^{\gamma}e^{sA}$ generated by $(-\prescript{}{t_{1}}{D}^{\beta}_{0+}-\prescript{}{t_{2}}{D}^{\gamma}_{0+}-A)$) to the forcing term $g_{\phi}(t_{1},t_{2},x)$. The solution to the Caputo problem (3.1) is then recovered by undoing the shift by $\phi_{1}$ and $\phi_{2}$,

[TABLE]

where $G_{\beta,\gamma}(r,s)$ is the transition density of the process generated by $(-\prescript{}{t_{1}}{D}^{\beta}_{0+*}-\prescript{}{t_{2}}{D}^{\gamma}_{0+*})$. Rearranging and using [kolokoltsov2019mixed, Equation 4.126] we have

[TABLE]

for $L^{\gamma}:=(-\prescript{}{t_{2}}{D}_{0+*}^{\gamma}-A)$ and $L^{\beta}:=(-\prescript{}{t_{1}}{D}_{0+*}^{\beta}-A)$. Here $E^{(\beta,\gamma)}[D]$ are generalised operator-valued Mittag-Leffler functions, which are introduced and extensively studied in the survey [kolokoltsov2019mixed],

[TABLE]

where $\mu_{0}^{\beta}(s)$ is the density of the exit time $\tau_{0}^{\beta}$ (cf. (2.1)).

3.2 Estimates for Greens function

As mentioned in the previous section, an application of the Dynkin formula followed by Doobs optimal stopping theorem gives the following stochastic representation of the solution (whenever it exists) to (3.1),

[TABLE]

A simple conditioning argument (see Appendix A), shows that this solution can be written as

[TABLE]

Inserting (2.1) for $\mu_{0}^{\alpha}$ and $\mu_{0}^{\beta}$,

[TABLE]

where $\mu_{0}^{\beta}(s)$ and $\mu_{0}^{\gamma}(s)$ are the densities of the exit times $\tau_{0}^{\beta}$ and $\tau_{0}^{\gamma}$, and

[TABLE]

and

[TABLE]

Thus, the Greens function associated to (3.1) are the coordinates of the integral kernel of the operator which acts on the boundary functions $\phi_{1}$ and $\phi_{2}$:

[TABLE]

Remark 4**.**

More generally, the function

[TABLE]

solves the boundary value problem

[TABLE]

where $\phi$ is a suitable function on the boundary of $X$.

For this reason, to obtain global two-sided estimates for the full Greens function $G_{full}=(G^{(\beta,\gamma)}_{1},G^{(\beta,\gamma)}_{2})$, it suffices to obtain estimates for $G_{1}^{(\beta,\gamma)}$, since the estimates for $G^{(\beta,\gamma)}_{2}$ will be the same up to exchanging coordinates. For the sake of readability we drop the subscripts from $G^{(\beta,\gamma)}_{1}$ and $t_{1}$ and look only at the function

[TABLE]

Making the substitution $s=t^{\beta}z$, we have

[TABLE]

where $G^{Y,\gamma}$ and $\mu_{0}^{\beta}$ are as in Lemma 2.1 and Lemma 2.2. Let $\Omega:=|x-y|^{2}t^{-\beta}$, $A=r^{\gamma}t^{-\beta}$.

Proposition 3.1**.**

For $(t,r,x,y)\in(0,\infty)\times(0,t_{2})\times\mathbb{R}^{d}\times\mathbb{R}^{d}$ and $t_{2}\in(0,\infty)$, the following estimates hold,

• •

For $\Omega\leq 1$,

[TABLE]

• •

For $\Omega\geq 1$,

[TABLE]

where

[TABLE]

Proof.

We sketch the main ideas of the proof here, see Appendix B for the full details. After applying Lemma 2.1 and Lemma 2.2 in (3.2), we end up with 4 integrals which contribute to the estimate for $G^{(\beta,\gamma)}$. For $\Omega\leq 1$, the main contribution comes from the integral

[TABLE]

After a substitution of $w=\Omega z^{-1}$, we immediately recognise the integral form of the incomplete gamma function, see Appendix C.1,

[TABLE]

Thus we have the two-sided estimate for $I_{1}$ for $\Omega\leq 1$,

[TABLE]

Since the integral $I_{1}$ is the main contributor to the estimate, this proves (3.5). For $\Omega\geq 1$, the main contribution to the estimate comes from the integral

[TABLE]

To estimate this integral, let $\alpha=\min(\beta,\gamma)$ and $\tilde{\alpha}=\max(\beta,\gamma)$. Then as an upper (resp. lower) bound for $I_{4}$ we replace the powers in the exponential term with $\alpha$ (resp. $\tilde{\alpha}$). That is, the upper estimate

[TABLE]

and the lower estimate

[TABLE]

Then an application of Proposition C.1 from the Appendix proves (3.6), and we are done. ∎

4 Extension to higher dimension

Let us outline how to extend the previous sections to the case where we have more than two fractional derivatives. Let $\mathcal{O}$ be the orthant in $\mathbb{R}^{k}$ defined by

[TABLE]

Let $\mathcal{O}_{i,0}$ denote the collection of vectors $t_{i,0}$ from $\mathcal{O}$ whose $i$-th coordinate is zero,

[TABLE]

Define $h_{i,0}(t)$ to be the projection of $\mathcal{O}_{i,0}$ onto the subspace $\mathcal{O}_{i}\subset\mathbb{R}^{k-1}$ by removing the coordinate which is zero, that is, $h_{i,0}(t):\mathcal{O}_{i,0}\mapsto\mathcal{O}_{i}$

[TABLE]

We look at the equations on $\mathcal{O}\times\mathbb{R}^{d}$,

[TABLE]

where each $\phi_{i}$ is a function on $\mathcal{O}_{i}\times\mathbb{R}^{d}$.

Remark 5**.**

In order to have continuity of the solution to the above boundary value problem, we would need to also impose additional boundary conditions in order to ensure that the solution coincides at the points where the boundary meets - i.e, at the origin. Without this additional assumption we only have a generalised solution, which is enough for our purposes.

As before, let $X_{t_{i}}^{\beta_{i}}(s)$ denote the process started at $t_{i}\in\mathbb{R}_{+}$ generated by $D^{\beta_{i}}_{0+*}$ where $\beta_{i}\in(0,1)$, and let $\tau_{0}^{\beta_{i}}$ denote the exit time of this process from $(0,\infty)$,

[TABLE]

Let $X^{\overline{\beta}}_{t}(s)=(X^{\beta_{1}}_{t_{1}}(s),\cdots,X^{\beta_{k}}_{t_{k}})$ be the process on $\mathcal{O}$ generated by

[TABLE]

and due to the independence of each process $X^{\beta_{i}}_{t_{i}}$, the exit time of $X^{\overline{\beta}}_{t}(s)$ from the orthant $\mathcal{O}$ is given by

[TABLE]

For $t\in\mathbb{R}^{k}_{+}$, let $B_{i}(t)$ denote the subset of $\mathcal{O}_{i}$ defined by

[TABLE]

i.e, $B_{i}$ consists of elements of the form

[TABLE]

The solution to (4.1) is given by

[TABLE]

Remark 6**.**

The same kind of conditioning argument works here, once the appropriate notation is adapted, see Appendix A.

Thus the objects we are interested in is

[TABLE]

where $(t_{i},x)\in\mathbb{R}_{+}\times\mathbb{R}^{d}$, and $(r,y)\in\mathcal{O}_{i}\times\mathbb{R}^{d}$. Note that

[TABLE]

and

[TABLE]

thus

[TABLE]

Focusing on the first coordinate, we have

[TABLE]

where $\mu_{0}^{\beta}(s)$ is the density of the exit time $\tau_{0}^{\beta_{1}}$, and $G^{Y,k}(s,r,x,y)$ is the density of the process

$(Y_{x}(s),X_{r_{2}}^{\beta_{2}}(s),X_{r_{3}}^{\beta_{3}}(s),\cdots,X^{\beta_{k}}_{r_{k}}(s))$. We assume that as usual $Y_{x}$ is a diffusion process, so that $p^{Y}$ satisfies Aronsons estimates,

[TABLE]

where the cross terms runs from $n=k-1$ down to $n=1$ in the above above and are the mixtures of long and short tails. Note that we use the convention that $\prod_{i=2}^{1}=1$.

Let $A_{1}=t_{1}^{-\beta_{1}}\prod_{i=2}^{k}r_{i}^{\beta_{i}}$ and $\Omega=|x-y|^{2}t^{-\beta_{1}}_{1}$.

Conjecture 4.1**.**

For $(t_{1},r,x,y)\in(0,\infty)\times\mathcal{O}_{1}\times\mathbb{R}^{d}\times\mathbb{R}^{d}$, we have the following two-sided estimates for the Greens function $G^{(\beta_{1})}$,

• •

For $\Omega\leq 1$,

[TABLE]

where $\Pi_{1}=\prod_{i=2}^{k}t_{1}^{-\frac{\beta_{1}}{\beta_{i}}}A_{1}^{-1-\frac{1}{\beta_{i}}}$.

• •

For $\Omega\geq 1$,

[TABLE]

where $\Pi_{2}=\left(\prod_{i=2}^{k}t_{1}^{-\frac{\beta_{1}}{\beta_{i}}}A_{1}^{-\frac{2-\beta_{i}}{2\beta_{i}(1-\beta_{i})}}\right)$, $\alpha=\min\{\beta_{1},\cdots,\beta_{k}\}$, and the powers $N_{1}$ and $N_{2}$ depend on $k$, $d$ and $\beta_{i}$ for $1\leq i\leq k$.

The calculations used to show the estimates in the case $k=2$ earlier in the article should work the same in this case, since the major contribution to the estimates should be the first term (respectively the last term) in (4.2) for small $\Omega$ (respectively large $\Omega$). We therefore omit the proof to avoid the cumbersome notation.

Appendix A Conditioning argument

Recall that $\mu_{0}^{\alpha}$ is the density of the random variable $\tau_{0}^{\alpha}$ and $p^{\alpha}_{s}(t,r)$ are the transition densities of the monotone process $X^{\alpha}_{t}(s)$ started at $t\in(0,\infty)$.

Proposition A.1**.**

For $\beta,\gamma\in(0,1)$,

[TABLE]

and similarly,

[TABLE]

Proof.

In the LHS of (A.1) condition first on $\{\tau_{0}^{\beta}=s\}$,

[TABLE]

where $p^{Y}(s,x,y)$ are the transition densities of the process $\left(Y_{x}(s)\right)_{s\geq 0}$ started at $x\in\mathbb{R}^{d}$. The proof of (A.2) is similar and is omitted. ∎

Appendix B Proof of Proposition 3.1

Let $A:=r^{\gamma}t^{-\beta}$, $\Omega:=|x-y|^{2}t^{-\beta}$. First we use Lemma 2.1 to estimate the density $\mu_{0}^{\beta}$, then we use Lemma 2.2 to estimate the spatial density

[TABLE]

where

[TABLE]

Now we have 4 regimes to consider, which are

• •

Case 1a): $A\leq 1$ and $\Omega\leq 1$

• •

Case 1b): $A\geq 1$ and $\Omega\leq 1$

• •

Case 2a): $A\leq 1$ and $\Omega\geq 1$

• •

Case 2b): $A\geq 1$ and $\Omega\geq 1$.

By directly comparing the powers of $z,\Omega$ and $A$ in the integrals above, we can reduce our attention to the integrals $I_{1}$ and $I_{4}$. Indeed for $\Omega\leq 1$ we have

[TABLE]

and

[TABLE]

For $\Omega\geq 1$ we have

[TABLE]

and

[TABLE]

Thus we have a preliminary two-sided estimate for $G^{(\beta,\gamma)}(t,r,x,y)$,

[TABLE]

and

[TABLE]

for some constants $C_{1},C_{2},C_{3},C_{4}$.

B.1 Estimates for $I_{1}$

For the first integral, we have for $A\leq 1$,

[TABLE]

Then for $\Omega\rightarrow 0$ and $A\rightarrow 0$,

[TABLE]

For $\Omega\rightarrow\infty$ and $A\rightarrow 0$,

[TABLE]

For $A\geq 1$ we have

[TABLE]

so for $\Omega\rightarrow 0$ and $A\rightarrow\infty$,

[TABLE]

For $\Omega\rightarrow\infty$ and $A\rightarrow\infty$,

[TABLE]

B.2 Estimates for $I_{4}$

For $A\leq 1$,

[TABLE]

Let $\alpha:=\min(\beta,\gamma)$ and $\tilde{\alpha}=\max(\beta,\gamma)$. For bounded $\Omega\leq 1$, we have

[TABLE]

where we have used (C.1). Next we use (C.2) to get for $\Omega\geq 1$,

[TABLE]

where $c=\frac{1}{1-\alpha}$ and $n=-\frac{d}{2}-1+\frac{1}{2(1-\beta)}+\frac{1}{2(1-\gamma)}$. Thus

[TABLE]

Finally for $A\geq 1$, we have

[TABLE]

For bounded $\Omega$, but unbounded $A$ we have

[TABLE]

For unbounded $\Omega$ and $A$, the term $A^{-\frac{1}{1-\gamma}}$ is negligable since $A$ is large, then we apply the usual Laplace approximation Proposition C.1 to get

[TABLE]

for $\Omega\rightarrow\infty$ and $A\geq 1$. For the lower bound of $I_{4}$, simply reverse the role of $\alpha$ and $\tilde{\alpha}$ in each case - otherwise structure of the estimates are the same.

Appendix C Asymptotic behaviour

We describe here an important method used in asymptotic analysis, called the Laplace method. Our main references for asymptotic analysis are [fedorjuk1987asymptotics, de1970asymptotic]. The goal of the Laplace method is to approximate integrals of the form

[TABLE]

as $A\rightarrow\infty$. As a motivating example, let $a=1,b=\infty$, $h(x)=x$ and $g(x)=x^{N}$ for some integer $N>0$111Of course this integral is just the upper-incomplete gamma function, whose asymptotic behaviour is well known.. In this case, one could integrate by parts $N$ times, until the $x^{N}$ term dissapears, and one is left with a final integral

[TABLE]

so that, for sufficiently large $A$,

[TABLE]

Now the main idea is that the major contribution to the asymptotic behaviour of

[TABLE]

comes from a neighbourhood around the point at which the function $h(x)$ attains its miniumum value. Outside this neighbourhood, the contributions to the asymptotic behaviour are exponentially small. Our standard references for asymptotic methods are [de1970asymptotic], [murray2012asymptotic] or [fedorjuk1987asymptotics]. 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 C.1**.**

Let $a>1$, $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)=a^{\frac{2(N+1)-1}{2(a+1)}}c^{\frac{2(N+1)+1}{2(a+1)}}\sqrt{\frac{2\pi}{a+1}}$, and $C_{2}(c,a)=c^{-\frac{1}{a+1}}[a^{\frac{1}{a+1}}+a^{-\frac{a}{a+1}}]$.

Further we have the slight extension to the above.

Proposition C.2**.**

Let $a,b>1$, $n\in\mathbb{R}$, and $c:=\min(a,b)$. Then as $\Omega A^{-1}\rightarrow\infty$,

[TABLE]

C.1 Incomplete Gamma function

Here we give some formulas that we use in the main body of the article. The upper complete Gamma function, is defined by

[TABLE]

Equivalently after a change of variables $y=Aw$,

[TABLE]

We have the following asymptotic behaviour of $\Gamma(s,A)$ for $A\rightarrow 0$,

[TABLE]

Thus, for $A\leq 1$,

[TABLE]

For $A\rightarrow\infty$, we use the Laplace method (C.1) with $h(x)=x$, $b=1$, $g(x)=x^{s-1}$,

[TABLE]