# The probabilistic point of view on the generalized fractional partial differential equations

• Vassili N. Kolokoltsov

## Abstract

This paper aims at unifying and clarifying the recent advances in the analysis of the fractional and generalized fractional Partial Differential Equations of Caputo and Riemann-Liouville type arising essentially from the probabilistic point of view. This point of view leads to the path integral representation for the solutions of these equations, which is seen to be stable with respect to the initial data and key parameters and is directly amenable to numeric calculations (Monte-Carlo simulation). In many cases these solutions can be compactly presented via the wide class of operator-valued analytic functions of the Mittag-Leffler type, which are proved to be expressed as the Laplace transforms of the exit times of monotone Markov processes.

## 1 Introduction

The simplest Cauchy problem for the linear equation of Caputo or Caputo-Dzherbashyan type

Da+βf(x)=λf(x),f(a)=Y,

with β > 0 is known to have the unique solution

f(x)=Eβ(λxβ)Y,

where Eβ is the Mittag-Leffler function. Lots of research in the modern theory of fractional PDEs is devoted to various extensions of this equation, when λ is replaced by an unbounded operator in some Banach space (for instance, a diffusion operator) and Da+β by various versions of generalized derivatives including mixtures of fractional derivatives.

Generalized fractional calculus was initially developed by extending fractional integrals to the integral operators with more general integral kernels and then defining the fractional derivatives as the compositions of these integrals with usual derivatives.

Our alternative approach was suggested in [41] and was motivated by probabilistic interpretation of fractional derivatives. It starts with the definition of the one-sided generalized fractional derivatives as the generators of monotone Markov or sub-Markov processes. The generalized fractional integrals is then defined as the corresponding right inverse operators. The objective of this paper is to overview the development of the probabilistic approach and to explain in detail how the most fundamental examples of generalized fractional operators and related fractional equations fit to the general probabilistic framework. To make the text accessible to the readers with a mild background in probability we stress everywhere when possible the analytic counterparts of the formulas and their interpretations.

The paper is organised as follows. In Section 2 we show how our generalized mixed fractional operators can be equivalently introduced from the three points of views, bringing together the languages and methods of probability theory, operator semigroups and generalized functions (and related pseudo-differential equations) and giving the meaning to the generalized fractional integrals in terms of Dynkin’s martingale, potential operators and fundamental solutions, respectively.

In Section 3 we recall the main generalized fractional operators discussed in the literature and derive their probabilistic representations as the generators or potential operators of appropriate Markov processes or semigroups. For completeness, we explain briefly how the introduction of these operators can be motivated via two analytic approaches: interpolation from integer-valued iterations and the dressing by the operators of multiplication and the change of variables. Thus we show how various different operators scattered through the literature arise in a unified way from simple general ideas.

Section 4 deals with the simplest linear fractional PDEs. Once the generalized fractional derivatives are expressed as the generators of certain Markov processes, probability theory provides a universal tool for obtaining solutions to various boundary-value processes via the so-called Dynkin martingale. This approach leads directly to the uniqueness of the solutions and their integral representations, but often fails to provide the existence of (sufficiently regular) solutions. However, in the case of monotone Markov processes that correspond to the equations with one-sided fractional derivatives the existence problem can be settled in a unified way leading to the introduction of the new class of generalized operator-valued Mittag-Leffler type functions. These functions, as their classical counterpart, can be represented as the Laplace transforms of positive measures, expressed in terms of the transition probabilities of the corresponding stochastic processes, or, in analytic language, the Green functions of the Cauchy problems for the corresponding generalized mixed-fractional derivative operators. These measures turn out to be the distributions of the exit times of the underlying processes. This representation implies both the well-posedness and natural regularity properties of the solutions in various classes of classical and generalized solutions. It is also well suited for numeric schemes, because of its explicit integral form.

Section 5 we touch upon more general classes of fractional linear problems including two-sided problems and the equations with higher order and partial fractional derivatives.

The following notations for function spaces will be used:

For a closed or open subset S of Rd, B(S) and C(S) are the Banach spaces of bounded measurable and continuous functions on S respectively, equipped with the sup-norm, C(S) is the closed subspace of C(S) consisting of functions vanishing at infinity, Cuc(S) is the closed subspace of C(S) consisting of uniformly continuous functions.

Ck(S) is a Banach space of k times continuously differentiable 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 closed S the derivatives on the boundary are understood as the limits of the derivatives defined in their interiors.

Apart from these standard notations, we introduce some more specific ones. For a subset AS let

CkillA(S)={fC(S):f|A=0},CconstA(S)={fC(S):f|Ais a constant}.

Moreover, for AS we shall consider the space C(A) to be the subset of B(S) obtained by setting the values of the functions to zero outside A.

By 1A we shall denote the indicator function of a set A. Specifically, 1a is the indicator of the half-line {ya}.

The letters E and P will be used to denote the expectation and probability with respect to various Markov processes.

## 2 Generalized fractional operators and Markov processes

### 2.1 Preliminaries: Standard fractional derivatives

Let Iaf be the integration operator defined on the set of continuous curves fC([a, b]) as Iaf(x) = axf(t) dt. Integration by parts yields

Ia2f(x)=ax(Iaf)(y)dy=ax(xy)f(y)dy.

Similarly by induction one gets the following formula for the iterated Riemann integral

Ianf(x)=1(n1)!ax(xt)n1f(t)dt,n=1,2,3,....(2.1)

This formula motivates the definition of the (left) fractional or Riemann-Liouville (RL) integral of order β > 0:

Iaβf(x)=Ia+βf(x)=1Γ(β)ax(xt)β1f(t)dt.(2.2)

Noting that the derivation is the inverse operation to usual integration, the definition (2.2) of the fractional integral suggests two notions of fractional derivative, the so-called RL (left) derivatives of order β ∈ (n, n + 1), n a nonnegative integer (where x > a):

Da+βf(x)=dn+1dxn+1Ia+n+1βf(x)=1Γ(n+1β)dn+1dxn+1ax(xt)nβf(t)dt,(2.3)

and the (left) Caputo-Dzherbashyan (CD) derivatives of order β ∈ (n, n + 1):

Da+βf(x)=Ian+1βdn+1dxn+1f(x)=1Γ(n+1β)ax(xt)nβdn+1dtn+1f(t)dt.(2.4)

### Remark 2.1

The CD derivative was initially introduced by Liouville (see [23]) and actively studied by Dzherbashyan, see e.g. [14], and by Caputo (see [65]). For more details on the history and basics of classical (standard) Fractional Calculus (FC) we recommend the book [70].

In this paper we shall be mostly concerned with the case β ∈ (0, 1) and its extensions, and occasionally look at the derivatives of higher order as the compositions of the derivatives of order β ∈ (0, 1). For β ∈ (0, 1) the definitions of the RL and CD derivatives turn to

Da+βf(x)=ddxIa+1βf(x)=1Γ(1β)ddxax(xt)βf(t)dt,x>a,(2.5)

and respectively

Da+βf(x)=Ia1βddxf(x)=1Γ(1β)ax(xt)βddtf(t)dt,x>a.(2.6)

As is seen by direct calculations (see e.g. Appendix in [41]), for smooth enough f,

Da+βf(x)=1Γ(β)0xaf(xz)f(x)z1+βdz+f(x)Γ(1β)(xa)β,(2.7)
Da+βf(x)=1Γ(β)0xaf(xz)f(x)z1+βdz+f(x)f(a)Γ(1β)(xa)β,(2.8)

implying

Da+βf(x)=Da+β[ff(a)](x)=Da+βf(x)f(a)Γ(1β)|xa|β.(2.9)

In particular, it follows that for smooth bounded integrable functions, the RL and CD derivatives coincide for a = −∞, and one defines the fractional derivative in generator form as their common value:

dβdxβf(x)=D+βf(x)=D+βf(x)=D+βf(x)=1Γ(β)0f(xz)f(x)z1+βdz.(2.10)

Another useful rewriting of (2.7) and (2.8) (used in Section 3) is obtained by the change of the variable of integration:

Da+βf(x)=1Γ(β)axf(y)f(x)(xy)1+βdz+f(x)Γ(1β)(xa)β,(2.11)
Da+βf(x)=1Γ(β)axf(y)f(x)(xy)1+βdz+f(x)f(a)Γ(1β)(xa)β.(2.12)

When β ∈ (0, 1) and x < a, the corresponding right derivatives can be introduced by the formulas

Daβf(x)=1Γ(β)0axf(x+z)f(x)z1+βdz+f(x)Γ(1β)(ax)β,(2.13)
Daβf(x)=1Γ(β)0axf(x+z)f(x)z1+βdz+f(x)f(a)Γ(1β)(ax)β,(2.14)

implying

Daβf(x)=Daβ[ff(a)](x)=Daβf(x)f(a)Γ(1β)(ax)β.(2.15)

The right RL and CD derivatives coincide for a = ∞, and one defines the fractional derivative in generator form as their common value:

dβd(x)βf(x)=Dβf(x)=Dβf(x)=Dβf(x)=1Γ(β)0f(x+z)f(x)z1+βdz.(2.16)

### 2.2 Generalized fractional differential operators

The fractional derivative dβf/dxβ, β ∈ (0, 1), from (2.10) was suggested as a substitute to the usual derivative df/dx, which can model some kind of memory by taking into account the past values of f. An obvious extension widely used in the literature (see e.g. [18], [62], [81], [26] and references therein) represent various mixtures of such derivatives, both discrete and continuous,

j=1Najdβjfdxβj,01dβfdxβμ(dβ).(2.17)

To take this idea further, one can observe that dβf/dxβ represents a weighted sum of the increments of f, f(xy) − f(x), from various past values of f to the ‘present value’ at x. From this point of view, the natural class of generalized mixed fractional derivative represent the causal integral operators

D+(ν)=Lν,Lνf(x)=0(f(xy)f(x))ν(dy),(2.18)

with some positive measure ν on {y : y > 0} satisfying the one-sided Lévy condition:

0min(1,y)ν(dy)<,(2.19)

which is just the condition ensuring that Lν is well-defined at least on the set of bounded infinitely smooth functions on {y : y ≥ 0}. The sign − is introduced to comply with the standard notation of the fractional derivatives, so that, for instance,

dβdxβf(x)=D+β=D+(ν)

with ν(y) = −1/[Γ (−β)y1+β] (note that Γ(−β) < 0).

The dual operators to Lν are given by the anticipating integral operators (weighted sums of the increments from the ‘present’ to any point ‘in future’):

D(ν)=Lν,Lνf(x)=0(f(x+y)f(x))ν(dy).(2.20)

As one of the most studied examples of D+(ν), let us mention the so-called tempered fractional derivatives given by the measure

ν(dy)=Ceλyy1β(2.21)

with constants C, λ > 0 (see [6] and [58]). Section 3.4 yields another point of view on this derivative.

Looking at (2.10) with ’probabilistic eyes’ one recognises in the operators −dβf/dxβ, β ∈ (0, 1) the generators of stable Lévy subordinators with the inverted direction (see e.g. [58], Ch.3 and [40], Ch.1 and 8), the operators −dβf/d(−x)β representing the generators of these subordinatots with the standard direction. Thus the probabilistic point of view suggests to take as fully mixed extensions of these operators the general Lévy subordinators, which leads again to the same operators (2.18) and (2.20), the condition (2.19) being the well known condition of the theory of Lévy processes defining the one-sided Lévy measures.

### Remark 2.2

Standard Lévy subordinators are defined as increasing processes, and thus operators (2.18) (which are more natural for fractional derivatives) are marked by primes.

One can weight differently the points in past or future depending on the present position, and one can also add a local part to complete the picture, leading to the operators

D+(ν,γ)=Lν,γl,Lν,γlf(x)=0(f(xy)f(x))ν(x,dy)γ(x)dfdx,(2.22)

with a non-negative b(x) and transition kernel ν(x, .), ∫ min(1, y)ν(x, dy) < ∞, which capture in full the idea of ’weighting the past’ and which can be called the one-sided, namely left or causal, weighted mixed fractional derivatives of order at most one. Symmetrically, one can define the right or anticipatingweighted mixed fractional derivatives of order at most one as

D(ν,γ)=Lν,γr,Lν,γrf(x)=0(f(x+y)f(x))ν(x,dy)+γ(x)dfdx.(2.23)

### Remark 2.3

Notice that Ll and Lr are dual only if ν and b do not depend on x.

From the probabilistic point of view the extension to (2.22) and (2.23) represents the transition from time homogeneous monotonic processes to arbitrary decreasing or increasing Feller processes; operators (2.22) and (2.23) are known to represent the general form of the generators of such processes (Courrège theorem, see e.g. [40]).

To link with the theory of pseudo-differential operators, it is worth noting that the operators D+(ν)andD(ν) are ΨDOs with the symbols −ψν(−p) and −ψν(p), where

ψν(p)=0(eipy1)ν(dy)

is the symbol of the operator Lν.

If ν is finite, the operators D+(ν) are bounded, which is not the case for usual derivatives. Thus the proper extensions of the derivatives represent only the operator D+(ν) arising from infinite measures ν satisfying (2.19). The operators arising from finite ν can be better described as the mixtures of finite differences approximating the derivatives.

The operators D±(ν) represent the extensions of the fractional derivatives D+βandDβ in generator form. Looking for the corresponding extensions of the operators Da±βandDa±β with a finite a we note that, by (2.8) and (2.14), Da+β (resp. Daβ) is obtained from D+β (resp. Dβ) by the restriction of its action to the space Cconst(−∞,a](R) (resp. Cconst[a,∞)(R). Therefore, the analogs of the CD derivatives should be defined as

Da+(ν)f(x)=0xa(f(xy)f(x))ν(dy)xa(f(a)f(x))ν(dy),Da(ν)f(x)=0ax(f(x+y)f(x))ν(dy)xa(f(a)f(x))ν(dy).(2.24)

By (2.9) and (2.13), the operators Da+βorDaβ, the analogs of the Riemann-Liouville derivatives, are obtained by further restricting the actions of D+βandDβ to the spaces Ckill(−∞,a](R) and Ckill[a,∞)(R):

Da+(ν)f(x)=0xa(f(xy)f(x))ν(dy)+xaf(x)ν(dy),Da(ν)f(x)=0ax(f(x+y)f(x))ν(dy)+xaf(x))ν(dy).(2.25)

Similarly, the derivatives with a finite a based on the general one-sided operators (2.23) are obtained by the same reduction of these operators, that is, say for the left derivatives, they are given by

Da+(ν,γ)f(x)=0xa(f(xy)f(x))ν(x,dy)(f(a)f(x))xaν(x,dy)+γ(x)dfdx,Da+(ν,γ)f(x)=0xa(f(xy)f(x))ν(x,dy)+f(x)xaν(x,dy)+γ(x)dfdx.(2.26)

From the probabilistic point of view it is then seen that the operators (2.8) are the generators of the modifications of the stable subordinators obtained by forbidding them (interrupting on an attempt) to cross the boundary x = a with an aR, that is, all jumps aimed to jump over the chosen barrier-point a are forced to land exactly at a. On the other hand, the operators (2.7) are the generators of the modifications of the stable subordinators obtained by killing them on an attempt to cross the boundary x = a. Applying the same procedure to operators (2.23) and (2.22) lead to the operators (2.24) and (2.26), which represent thus the generators of Markov processes interrupted or killed on an attempt to cross the boundary point a.

Finally one can further extend these mixed derivatives by including additional killing mechanisms or by mixing killing and stopping (say, by working with the linear combinations of of CD and RL derivatives). Such extension leads to the following mixed derivative operators:

Da+(ν,γ,S,R)f(x)=0xa(f(xy)f(x))ν(x,dy)(f(a)f(x))S(a,x)+f(x)K(a,x)+γ(x)dfdx,(2.27)

where S(a, x) and K(a, x) are the rates of stopping and killing respectively. Some formulas simplifies (this is especially relevant for the extensions to manifolds) if one counts jumps by their finite points (rather than relative to x) leading to the following modified version of (2.27):

Da+(ν,γ,S,R)f(x)=ax(f(s)f(x))ν~(x,ds)(f(a)f(x))S(a,x)+f(x)K(a,x)+γ(x)dfdx.(2.28)

### Remark 2.4

One can show that these operators Da+(ν,γ,S,R) represent the most general generators (under some mild technical condition) of decreasing sub-Markov Feller processes on [a, ∞).

### 2.3 Generalized fractional integrals: shift invariant case

Let us see what should be the proper analog of the fractional integral exploiting approaches from probability, semigroup theory and the generalized functions/ΨDEs theory. Let us consider first the case of operators (2.24) arising from the Lévy subordinators. We shall talk about the left derivatives for definiteness.

The operator dβ/dxβIβ acts as the identical operator on functions with a compact support. Hence, in the language of the semigroups of operators, the operator Iβ is the potential operator of the strongly continuous semigroup of linear operators in C(R) generated by −dβ/dxβ, that is the limit of the resolvent operator Rλ = (λ + dβ/dxβ)−1, as λ → 0. Notice that Iaβ is just the reduction of Iβ to the space Ckill(−∞,a](R). As is worth noting Iβ is an unbounded operator in C(R), but is bounded when reduced to Ckill(−∞,a]((−∞, b]) for any b > a. Thus we can define the generalized fractional integralIaν as the potential operator of the semigroup generated by Lν reduced to the space Ckill(−∞,a](R).

On the other hand, the usual fractional integral (2.2) solves the boundary value problem Dβg = f with g(x) = 0 for xa. Recall that for a Feller process Xt(x) and a function f from the domain of its generator L the process f(Xt(x)) − 0t Lf (Xs(x)) ds is a martingale, called Dynkin’s martingale (see e.g. [40]). Applying to this martingale Doob’s optional sampling theorem shows that

f(x)=Ef(Xτ(x))E0τg(Xs(x))ds,(2.29)

where g = Lf and τ is the exit time of Xt(x) from a domain (at least if Eτ < ∞). This implies that

Iaβf(x)=E0τxf(xXtβ)dx,(2.30)

where Xtβ is the stable subordinator and τx is the time when Xtβ reaches the point xa (and thus xXtβ reaches a). Therefore, from the probabilistic point of view, the generalized fractional integral, representing the analog of Iaβ for the case of the generalized derivative Da+(ν), is the path integral

Iaνf(x)=E0τxf(xXtν)dt,(2.31)

where Xtν is the Lévy subordinator generated by operator (2.20) and τx the time for this process to reach a.

Finally, as is known (see e.g. [42]), the fundamental solution (vanishing on the negative half-line) to the fractional derivative operator dβ/dxβ is Uβ(x) = x+β1/Γ (β), so that the usual fractional integral Iaβf(x) represents the integral operator with the kernel being the fundamental solution of dβ/dxβ or, in other words, the convolution with this fundamental solution, restricted to the space Ckill(−∞,a](R). Thus, from the point of view of the theory of PDEs and ΨDEs, the generalized fractional integral, representing the analog of Iaβ for the case of the generalized derivative Da+(ν), is the integral operator

Ia(ν)f(x)=0xaf(xz)U(ν)(dz),(2.32)

where U(ν)(dz) is the fundamental solution to the operator Lν (we denoted it U(ν)(dz), as it turns out to be a measure, see below).

These three facets of the generalized fractional integral were given by analogy. Let us see now that they are all well defined and in fact represent the same objects. It is well known that the operators Lν and Lν generate Feller processes on R (called Lévy subordinators) and the corresponding strongly continuous semigroups in spaces C(R) and Cuc(R). The latter space is much less used as the former, but is handy for us as it includes the spaces of functions that are constants on the halflines. Recall that the potential measure is defined as the integral kernel of the potential operator. The potential measure for any Lévy subordinator is known to exist and to equal the vague limit

U(ν)(M)=0G(ν)(t,M)dt

of the measures 0KG(ν)(t, .) dt, K → ∞, such that U(ν)(M) is finite for any compact M whenever ν does not vanish (see e.g. [72] or [44]). More precisely, for any λ > 0,

U(ν)([0,z])eλzϕν(λ),ϕν(λ)=ψν(iλ)=0(1eλy)ν(dy).(2.33)

Here G(ν)(t, dy) is the Green function of the Cauchy problem for the operator Lν. Thus the definition of the Ia(ν) from the point of view of the semigroup theory is correct.

Turning to the probabilistic definition we note that the process obtained from a subordinator by killing on the attempt to cross the boundary is a well defined Feller sub-Markov process, so that formula (2.52) arising from a Dynkin’s martingale is well defined and represents a unique solution to the corresponding boundary problem. From this uniqueness it follows that the definitions of the generalized integral from the semigroup and probabilistic points of view coincide.

Finally, from the definitions of the potential measure it follows that the potential measure represents a fundamental solution vanishing on the negative half-line. Therefore to fix the definition arising from the PDEs theory it is only needed to show the uniqueness of such fundamental solution. This is the content of the following assertion from [43] (see also [44]), which also includes the λ-potential measures defined as

Uλ(ν)(A)=0eλtG(ν)(t,A)dt.(2.34)

These measures are known to be bounded (see e.g. [72] or [44]):

Uλ(ν)=0Uλ(ν)(dx)1ϕν(λ).(2.35)

### Proposition 2.1

Let the measureνon {y : y > 0} satisfy(2.19).

1. For anyλ > 0, theλ-potential measureUλ(ν)represents the unique fundamental solution of the operatorλLν.

2. If the support ofνis not contained in a lattice {αn, nZ}, with someα > 0, the measureU(ν)(dy) represents the unique fundamental solution to the operatorLν, up to an additive constant.

3. Let {αn, nZ} be the minimal lattice (that cannot be further rarified) containing the support ofν, so that for anykZ, k > 1, there existsnZsuch thatαnbelongs to the support ofνandn/kZ. Then any two fundamental solutions to the operatorLνdiffer by a linear combination of the type

G(x)=nZanexp{2πnix/α}(2.36)

with some numbersan. In particular, U(ν)(dy) is again the unique fundamental solution vanishing on the negative half-line.

The following result summarizes the properties of generalized fractional integrals.

### Proposition 2.2

1. Let a measureνon {y : y > 0} satisfy(2.19). For any generalized functiongD′(R) supported on the half-line [a, ∞) with anyaR, and anyλ ≥ 0, the convolutionUλ(ν) * gwith theλ-potential measure(2.34)is a well-defined element ofD′(R), which is also supported on [a, ∞). This convolution represents the unique solution (in the sense of generalized function) of the equation (λLν)f = g, or equivalently

D+(ν)f=λf+g,

supported on [a, ∞).

2. Ifλ > 0 andgC(R) ∩ Ckill(−∞,a](R), then

f(x)=(Uλ(ν)g)(x)=Rλg(x)=g(xy)Uλ(ν)(dy)=0xag(xy)0eλtG(ν)(t,dy)dt(2.37)

belongs to the domain of the operatorLν(considered as the generator of the Feller semigroup onC(R)) and thus represents the classical solution to the equation (λLν)f = g, or equivalently

D+(ν)f=Da+(ν)f=Da+(ν)f=λf+g.(2.38)
3. If reduced to the space Ckill(−∞,a]((−∞, b]) withb > a(this space is invariant underTtand hence under allRλ), the potential operatorR0with the kernelU(ν)becomes bounded and hence

(U(ν)g)(x)=R0g(x)=Ia(ν)g(x)(2.39)

belongs to the domain ofLνand thus represents the classical solution to the equation

Lνf=D+(ν)f=Da+(ν)f=Da+(ν)f=g(2.40)

onCkill{a}([a, b]).

### Proof

1. Since the generalized functions Uλ(ν) and g have support bounded from below, their convolution Uλ(ν) * g is well-defined and solves the equation (λLν)f = g. Uniqueness follows from the uniqueness of the fundamental solution.

2. Since Lν generates a semigroup, which preserves the space Ckill(−∞,a] (R), this space is also invariant under the resolvent Rλ=(λLν)1. Since the image of the resolvent always coincides with the domain of the generator (see e.g. [16]), Rλg belongs to the intersection of the space Ckill(−∞,a](R) and the domain of D+(ν).

3. As was mentioned above, the potential operator

R0g(x)=(U(ν)g)(x)=0xag(xy)U(ν)(dy)

is bounded on Ckill(−∞,a](R) ∩ C((−∞, b]). Hence the required equation can be obtained from (ii) by passing to the limit λ → 0. □

In particular, applying (2.37) to Lν = −dβ/dxβ and comparing with the classical solution of the fractional linear equation via the Mittag-Leffler function yields the integral representation

βzβ1Eβ(λzβ)=0eλtGβ(t,z)dt=Uλβ(z),(2.41)

which is equivalent to Zolotarev’s formula

Eβ(s)=1β0esxx11/βGβ(1,x1/β)dx,(2.42)

thus yielding a proof of this formula from the semigroup theory by-passing subtle analytic manipulations of Zolotarev’s initial derivation (see [77] and [78]).

For general ν the classical interpretation of the solution Rλg(x) is subtle for gC([a, ∞)) not vanishing at a. However, Rλg(x) may well belong to the domain locally, outside the boundary point a. And in fact, the requirement for the solution to belong to the domain outside a boundary point is common for the classical problems of PDEs. The following assertion illustrates this point concretely.

### Proposition 2.3

Under the assumptions of Proposition 2.2 let the potential measureU(ν)(dy) have a continuous density, U(ν)(y), with respect to Lebesgue measure. LetgC1[a, b] and is continued as zero to the left ofa. Then the functionf(x) = R0g(x) is continuously differentiable in (a, b]. Hence it satisfies the equationDa+(ν)f = glocally, at all points from (a, b].

### Proof

From the formula for R0g(x) it follows that

which is well-defined and continuous for xa. The limit from the right of (d/dx)R0g(x) as xa is g(a) U(ν)(0), which may cause a jump when this function crosses the value x = a. □

Apart from generalized solutions arising from duality as considered above, one uses also the notions of generalized solution by approximation. Namely, for a measurable bounded function g(x) on [a, ∞), a continuous curve f(x), ta, is the generalized solution via approximation to the problem D+(ν)f = −λ f + g on C([a, b]), if there exists a sequence of curves gn(.) ∈ Ckill(−∞, a](R) such that gng a.s., as n → ∞, and the corresponding classical (i.e. belonging to the domain) solutions fn(x), given by (2.37) with gn(x) instead of g(x), converge point-wise to f(t), as n → ∞.

The following assertion is a consequence of Proposition 2.2.

### Proposition 2.4

For any measurable bounded functiong(x) on [a, ∞), formula(2.37) (resp. (2.39)) supplies the unique generalized solution by approximation to problem(2.38) (resp. (2.40)) on [a, b] for anyb > a.

### 2.4 Generalized fractional integrals: weighted mixed derivatives

Let us turn to operators (2.22). It is known (see e.g. [40]) that the operator Lν,γl generates a conservative Feller semigroup Tt in C(R) with invariant core C1(R) whenever the following conditions hold:

1. bC1(Rd) and b ≥ 0;

2. ν (x, dy), the gradient of the Lévy kernel with respect to x, exists in the weak sense as a signed measure that depends weakly continuously on x, in the sense that ∫ f(y) ∇ ν(x, dy) is a continuous function for any fC(Rd) with a support separated from zero;

3. supxmin(1,y)ν(x,dy)<,supxmin(1,y)|ν(x,dy)|<,(2.43)

and for any ϵ > 0 there exists a K > 0 such that

supxKν(x,dy)<ϵ,supxK|ν(x,dy)|<ϵ,supx01/Kyν(x,dy)<ϵ.(2.44)

The next result shows the existence of the potential measures describing the integral kernel of the resolvent and potential operators of this Feller semigroup:

Rλg(x)=0eλt(Ttg)(x)dt=0dteλtxP(ν,γ)(t,x,dy)g(y)=xΠ(ν,γ)λ(x,dy)g(y),(2.45)

with λ ≥ 0, where P(ν,γ) denote the transition probabilities of the semigroup Tt and Π(ν,γ)λ(x, dy) the λ-potential measure.

### Proposition 2.5

Let a kernelν(x, dy) and a functionbsatisfy assumptions (i)-(iii) above. Let the inequality

ν(x,dy)ν~(dy)(2.46)

hold with some non-vanishing measureν̃ satisfying(2.19). Then the following holds:

1. For any nonincreasing functionfwe have the comparison principle for semigroups:

TtfT~tf,(2.47)

whereT~tis the semigroup generated by the operatorLν~.

2. The potential operator Π(ν,γ) = Π(ν,γ)0of the semigroupTtis well defined and satisfies the comparison principle for potential operators:

(Π(ν,γ)1a)(x)=Π(ν,γ)(x,[a,x])U(ν~)([a,x])eλzϕν(λ),x>a,(2.48)

where by Π(ν,γ)we denoted both the integral operator and the measure representing its integral kernel (see(2.33)for the last estimate).

3. The same holds for theλ-potential operatorsΠ(ν,γ)λ, λ > 0, x > a:

(Rλ1a)(x)=(Π(ν,γ)λ1a)(x)=Π(ν,γ)λ(x,[a,x])Uλ(ν~)([a,x])1ϕν(λ)(2.49)

(see(2.35)for the last estimate).

4. If additionallyν(x, dy) ≤ ν̄ (dy) with some non-vanishingν̄ (dy) satisfying(2.19), thenTt fT¯tffor a non-increasing functionf.

### Proof

Let us consider the case of b = 0 (by approximating the derivative with finite differences we can reduce the general case to this one).

1. Notice that

(LLν~)f(x)=0(f(xy)f(x))(ν(x,dy)ν~(dy)),

which is positive for any nonincreasing function f. Let us write the difference between the actions of Tt and T~t in the standard form via the difference of the generators (see e.g. [39]):

TtfT~tf=0tTtsf(LLν~)T~sf.

Since T~s preserves the set of nonincreasing functions, T~sf is nonincreasing. Hence, (LLν~)T~sf ≥ 0 and hence Tt fT~tf ≥ 0, yielding (2.47).

2. By changing f to − f it follows from (2.47) that

TtfT~tf(2.50)

for any nondecreasing function f.

Applying (2.50) to the indicator functions of half-lines shows that

P(Xt(x)>c)P(xXt>c),

where Xt is the decreasing process generated by L and Xt is the subordinator defined by the measure ν̃. By the definition of potential operator of the process Xt it equals (whenever it exists)

Π(ν,γ)f(x)=0Ttf(x)dt.

In particular,

(Π(ν,γ)1a)(x)0(T~ν1a)(x)dx=U(ν~)([a,x]),

yielding (2.48).

3. (iii), (iv) are fully analogous to (ii). □

Consequently, under the assumptions of Proposition 2.5, the potential operator Π(ν,γ) is an integral operator, which is well defined for functions with a support bounded below. It is a bounded operator when reduced to the spaces of functions with a fixed compact support. The integral kernel of the potential operator, Π(ν,γ)(x, dy), is such that, for any x, the measure Π(ν,γ)(x, ⋅ .) is supported on the set (−∞, x]. Consequently, from the point of view of the semigroup theory, we can define the generalized mixed fractional (weighted) integralIa(ν,γ) as the potential operator Π(ν,γ) reduced to the space Ckill(−∞,a](R):

Ia(ν,γ)g(x)=axg(y)Π(ν,γ)(x,dy)(2.51)

On the space Ckill(−∞,a](R) the composition Da+(ν,γ)Ia(ν,γ) acts as identity, as it should be.

Since the process generated by L(ν,γ)l and stopped at the attempt to cross the given boundary-point a is a well defined Feller process with a regular boundary point, we can apply the theory of Dynkin’s martingale that implies that if a function g solves the boundary-value problem Dβg = f with g(x) = 0 for xa, then it can be expressed as the path integral yielding the probabilistic definition of the generalized fractional integral:

Ia(ν,γ)g(x)=E0τxg(Xtν(x))dt,(2.52)

where Xtν(x) is the decreasing process generated by operator (2.23) started at x and τx the time for this process to reach a.

As the fundamental solutions are less in use for the operators that are not shift invariant, we shall not give detail on this interpretation here.

The results of the previous section extend now automatically to the case of weighted integrals. For instance, the following holds.

### Proposition 2.6

Let the assumptions of Proposition 2.5 hold.

1. Ifλ > 0 andgC(R) ∩ Ckill(−∞,a](R), then

f(x)=Rλg(x)=xΠ(ν,γ)λ(x,dy)g(y))

belongs to the domain of the operatorLν,γland thus represents the classical solution to the equation (λLν,γl)f = g, or equivalently, to the equation

D+(ν,γ)f=Da+(ν,γ)f=Da+(ν)f=λf+g.(2.53)
2. If reduced to the space Ckill(−∞, a]((−∞, b]) withb > a(this space is invariant underTtand hence under allRλ), the potential operatorR0becomes bounded and hence the integral(2.51)belongs to the domain ofLν,γland thus represents the classical solution to the equation(2.53)withλ = 0.

Similarly, Proposition 2.4 on the generalized solutions has the straightforward extension.

The mixed fractional integrals can be also defined when arising from more general fractional derivatives (2.28). Namely, they represent the potential operators of the semigroups and processes generated by (2.28) and killed at the boundary. Their probabilistic or path integral representation can be written as in (2.52):

Ia+(ν,γ,S,R)f(x)=E0τxf(Xt(x))dt,(2.54)

where Xt is the decreasing sub-Markov process generated by the operator

Da+(ν,γ,S,R)f(x)=ax(f(s)f(x))ν~(x,ds)+f(x)[K(a,x)+S(a,x)]+b(x)dfdx.

Alternatively, using the Feynman-Kac theory, it can be rewritten as

Ia+(ν,γ,S,R)f(x)=E0f(Yt(x))exp{0t[K(a,Ys(x))+S(a,Ys(x))]ds}dt,(2.55)

where Yt is the decreasing Markov process on [a, ∞) generated by the operator

Da+(ν,γ,S,R)f(x)=ax(f(s)f(x))ν~(x,ds)+b(x)dfdx.

### Remark 2.5

It is natural to ask (especially in order to link our mixed fractional integrals with generalized integrals from [4]), which integral operators can be represented as our mixed fractional integrals, that is, as potential operators to the semigroups generated by decreasing processes. This question is essentially answered in Chapter 2, Section 4 of [16]. Namely, in this book the criterion is given for operators to represent potential operators of killed Markov processes (even though this criterion is not very easy to check). In our case just the additional requirement arises that the kernel Π(x, dy) of such operator should have support on the half-line (−∞, x] for any x.

### 2.5 Further extensions of mixed fractional operators

Extensions of the fractional operators and the fractional differential equations introduced above corresponds to monotone stochastic processes. They do not exhaust all interesting problems. Apart from the extensions with β ∈ (1, 2) (which we shall not touch here, see [41]) important additional situations include two-sided problems (very natural in the fractional calculus of variations) and multidimensional problems. Unlike the problems considered above, working with potential operators for these cases is usually much more complicated, because we cannot in general use the potential operator of the free problem (without a boundary) and just reduce it to problems with a one-sided support. On the other hand, the representation via Dynkin’s martingale works usually equally well in all cases. However, it gives only uniqueness of the solution and the integral representation for it, that is, a generalized solution in probabilistic sense. The difficulty lies in the identification of the analytic characteristics of thus obtained generalized solutions and of the conditions when they are classical.

The most general extension of mixed fractional derivatives within the probability theory arises from an arbitrary Markov process by interrupting it on the attempt to cross a boundary. We refer to [41] for detail and only mention here briefly three examples.

1. Two-sided problem mixed with a diffusion (see detail in [21]). The typical case is the problem

(Da+β1+Dbβ2+A)f=g,f(a)=fa,f(b)=fb,

with a given function g and a diffusion operator A. Extending the derivatives Da+β1andDbβ2 to fully mixed derivatives − Lν and Lν leads to the two-sided problem

(L+A)f=g,f(a)=fa,f(b)=fb,,

with

Lf(x)=axbx(f(x+y)f(x))ν(x,dy)+ga(x)dfdx+(f(b)f(x))bxν(x,dy)+(f(a)f(x))axν(x,dy).(2.56)

The generalized solution to this problem (the two-sided analog of the fractional integral) can be given by Dynkin’s martingale (2.29) with τ the exit time from (a, b) of the process Xt(x) generated by L + A, or simply the process generated by

axbx(f(x+y)f(x))ν(x,dy)+ga(x)dfdx+Af(x),

as these two processes coincide before exiting (a, b). Spectral problems for two-sided derivatives will be analyzed in Section 5.1.

2. One-sided multidimensional fractional equations. These are the equations in the orthant

O={x=(x1,,xk):xjaj}

of the type

(Da1+β1++Dak+βk)f(x1,,xk)=g(x1,,xk),

where Dj acts on the variable xj. When f is given on the boundary of the orthant O, the generalized solution is well defined by Dynkin’s martingale (2.29) with τ the exit time from the interior of O.

3. Fully multidimensional case. The analog of RL derivative arising from a process in Rd and a domain DRd is the generator of the process killed on leaving D. For CD version this is more subtle, as we have to specify a point where a jump crosses the boundary. The most natural model assumes that a trajectory of a jump follows the shortest path. Namely, assume that A is a generator of a Feller process Xt(x) in Rd with the generator of type

Af(x)=(γ(x),)f(x)+Rd(f(x+y)f(x))ν(x,dy)(2.57)

with a kernel ν(x, .) on Rd ∖ {0} such that

supxRdmin(1,|y|)ν(x,dy)<,(2.58)

that is, in the terminology of [39], [40], a generator or order at most one.

Let D be an open convex subset of Rd with boundary D and closure D̄. For xRd and a unit vector e let Lx,e = {x + λe, λ ≥ 0} be the ray drawn from x in the direction e. For xD, let

λ(x,y/|y|)=max{R>0:x+Ry/|y|D¯}RD(x,y)=x+y,if|y|λ(x,y/|y|),x+λ(x,y/|y|)y/|y|,if|y|λ(x,y/|y|).(2.59)

The process Xts(x) in D interrupted and stopped on an attempt to cross D can be defined by the condition of stopping at D and the generator

which represents a multidimensional extension of the CD boundary operator. The unique generalized solution for a boundary value problem AD*f = g with a given function g and given values of f on D can be again given via Dynkin’s martingale.

Finally, in the spirit of the most general extension of one-sided mixed derivatives (2.28), the most general multidimensional mixed derivative of order at most one in the open set SRn can be defined as the pseudo-differential operator of order at most one generating (with negative sign) a sub-Markov process in the closure S̄:

Da+(ν,γ,S,R)f(x)=S¯(f(s)f(x))ν~(x,ds)+f(x)K(x)+(b(x)f(x)))(2.61)

with a measure ν such that supx ∫ |sx| ν(x, ds) < ∞.

In all cases above the corresponding mixed fractional integral can be defined as the potential operators of the corresponding sub-Markov processes, with their path integral or probabilistic interpretation given in terms of Dynkin’s martingale.

By the probabilistic representation for various generalized derivatives discussed below we shall understand their representation in forms (2.26), (2.27), (2.28), (2.61), which explicitly reveal the structure of the Markov process they generate.

## 3 Probabilistic representations for the fundamental fractional operators

Apart from the basic RL fractional operators, the next popular fractional operators are possibly the Hadamard and the Erdélyi-Kober operators. We shall now remind their definitions together with their popular extensions and then derive their probabilistic representations. Next we shall describe two general analytic approaches, which can be used to motivate the introduction of these operators and their further extensions. Finally we shall briefly discuss the multiplicative interpolations of derivatives introduced by Hilfer and possible generalizations leading to the random fractional operators.

The Hadamard fractional integral of fractional order q > 0 and a boundary point a > 0 is defined by the formula (where x > a)

HIaqg(x)=1Γ(q)axlnxsq1g(s)sds=1Γ(q)1x/a(lny)q1g(x/y)dyy,(3.1)

the second version being obtained by the change x/s = y, s = x/y, ds = −xdy/y2. The corresponding Hadamard fractional derivative of RL type of order q ∈ (0, 1) is defined as

HDa+qf(x)=xddxHIa1qf(x)=xddx1Γ(1q)1x/a(lny)qf(x/y)dyy,(3.2)

and the Hadamard fractional derivative of CD type of order q ∈ (0, 1) as

HDa+qf(x)=HIa1qxddxf(x)=1Γ(1q)axlnxsqf(s)ds,(3.3)

both again for x > a.

### Proposition 3.1

Iffis continuously differentiable, then

HDa+qf(x)=1Γ(1q)(lnxa)qf(x)
qΓ(1q)0xalnxxyq1(f(xy)f(x))dsxy=1Γ(1q)(lnxa)qf(x)qΓ(1q)axlnxsq1(f(s)f(x))dss,(3.4)
HDa+qf(x)=[HDa+q(f(.)f(a))](x)=1Γ(1q)lnxaq(f(x)f(a))qΓ(1q)axlnxsq1(f(s)f(x))dss.(3.5)

### Proof

We get from (3.2) that

HDa+qf(x)=xa1Γ(1q)lnxaqf(a)ax+xΓ(1q)1x/a(lny)qf(x/y)dyy2=1Γ(1q)lnxaqf(a)+1Γ(1q)axlnxsqf(s)ds.(3.6)

Using f′(s) = (f(s) − f(x))′ and integrating by parts and taking into account that

limsx(f(s)f(x))(ln(x/s))q=0,

yields

HDa+qf(x)=1Γ(1q)lnxaqf(a)1Γ(1q)(lnxa)q(f(a)f(x))qΓ(1q)axlnxsq1(f(s)f(x))dss=1Γ(1q)lnxaqf(x)qΓ(1q)0xa(lnxxy)q1(f(xy)f(x))dyxy.

On the other hand, the r.h.s. of (3.3) equals the second term in (3.6) implying (3.41). □

Since

xaxlnxxyq1dyxy=0alnxsq1dss=1qlnxaq,

one can write equivalently that

HDa+qf(x)=f(x)xaxν(x,y)dy0xa(f(xy)f(x))ν(x,y)dy,ν(x,y)=qΓ(1q)(xy)(lnxxy)q1.(3.7)

This representation shows that the Hadamard derivatives represent particular case of the general operators (2.26), revealing the probabilistic meaning of these operators: HDa+q is seen to generate the process on [0, ∞) obtained from the process generated by the operator

HD0+qf(x)=0x(f(xy)f(x))ν(x,y)dy,(3.8)

by killing the latter process when attempting to cross the boundary x = a.

Notice also that for a = 0 both derivatives are well defined and coincide with operator (3.8) (with inverse sign), which generates a well defined process on [0, ∞) and represents a Hadamard analog of the fractional derivative in generator form (2.10).

Kilbas in [27] introduced the generalized Hadamard (HK) operators:

HKIa,μqg(x)=1Γ(q)ax(s/x)μlnxsq1g(s)sds,x>a>0,(3.9)
HKDa+,μqf(x)=xμ(xddx)xμHIa,μ1qf(x),(3.10)
HKDa+,μqf(x)=HIa,μ1qxμ(xddx)xμf(x).(3.11)

Extending the calculations above yields the following result (we omit similar calculations, more general case is considered below):

### Proposition 3.2

Iffis continuously differentiable, then

HKDa+qf(x)=1Γ(1q)lnxaq(x/a)μf(x)ax(f(s)f(x))νHK(x,s)ds+f(x)μΓ(1q)axlnxsq(x/s)μdss.(3.12)
HKDa+qf(x)=1Γ(1q)lnxaq(x/a)μ(f(x)f(a))ax(f(s)f(x))νHK(x,s)ds+f(x)μΓ(1q)axlnxsq(x/s)μdss(3.13)

with

νHK(x,s)=qsΓ(1q)lnxsq1xsμ.

The distinguishing property of this case as compared with the standard Hadamard case is the presence of the killing term also in the CD version of the derivative.

For a = 0 both derivatives coincide and equal

HKD0+qf(x)=HD0+qf(x)=0x(f(s)f(x))νHK(x,s)ds+f(x)μΓ(1q)0xlnxsq(x/s)μdss=0x(f(s)f(x))νHK(x,s)ds+μqf(x).(3.14)

This operator (with the negative sign) generates a decreasing process on [0, ∞) with the uniform killing rates μq.

This process and its stopped versions arising from a > 0 can be well called the Hadamard-Kilbas (HK) processes. Of course, they belong to the general class of processes generated by (2.28).

### 3.2 Erdélyi-Kober-Kiryakova-Luchko fractional operators

Now we turn to the Erdélyi-Kober (EK) operators. Let us deal directly with the generalized Erdélyi-Kober integrals, see e.g. in the book [70, Sect.18.1]. The corresponding fractional calculus was developed by Kiryakova and Luchko (see e.g. [30], [76], [31], [33], [34]) and we shall refer to these operators as the Erdélyi-Kober-Kiryakova-Luchko (EKKL) operators. The fractional EKKL integrals are defined by the formulas

EKKLIa,βγ,qg(x)=βxβ(γ+q)Γ(q)0x(xβsβ)q1sβγ+β1g(s)ds=Ωβx(γ+q)I0qxγΩβ1g(x),x0,(3.15)

where I0q is the standard RL fractional integral and Ωβ is the operator changing the variable: Ωβf(x) = f(xβ). The choice of the lower bound a = 0 makes these operators quite specific. Namely, by the change of the variable the integral (3.15) rewrites in another useful form:

EKKLIa,βγ,qg(x)=βΓ(q)01(1uβ)q1uβγ+β1g(ux)du.(3.16)

This expression shows that all these integral operators with various β, γ, q commute.

The classical Erdélyi-Kober operators are given by this formula with β = 1 and β = 2, by Kober, Erdélyi and Sneddon, see [70]. For q ∈ (0, 1), the corresponding RL type derivatives can be defined (see [30]) via the three equivalent expressions (simple check shows that these expressions coincide):

EKKLD0+,βγ,qf(x)=x1ββγΓ(1q)ddx0x(xβuβ)quβ(γ+q+1)1f(u)du,(3.17)
EKKLD0+,βγ,qf(x)=xβddx+γ+1βxβ(1+γ)Γ(1q)0x(xβuβ)quβ(γ+q+1)1f(u)du,(3.18)

and

EKKLD0+,βγ,q=Ωβ(xγddxI0+1qxγ+q)Ωβ1.(3.19)

For a Caputo type EK derivative corresponding to (3.17), see [54].

The probabilistic representation of these derivatives is as follows.

### Proposition 3.3

Iffis continuously differentiable, then

EKKLD0+,βγ,q=0x(f(s)f(x))νEKK(x,s)ds+f(x)Γ(γ+q+1)Γ(1+γ),(3.20)

with

νEKKL(x,s)=βxβγΓ(q)(xβsβ)q1sβ(γ+q+1)1.(3.21)

### Proof

Changing the variable of integration in (3.17) to u = xs and then differentiating yields

EKKLD0+,βγ,q=x1ββγΓ(1q)ddxxβ(γ+1)01(1sβ)qsβ(γ+q+1)1f(xs)ds=β(γ+1)Γ(1q)01(1sβ)qsβ(γ+q+1)1f(xs)ds+xΓ(1q)01(1sβ)qsβ(γ+q+1)1sf(xs)ds.

Returning back to the variables s = u/x yields

EKKLD0+,βγ,q=xβ(γ+1)Γ(1q)0x(xβuβ)quβ(γ+δ+1)1[β(γ+1)f(u)+uf(u)]du.

Integrating by parts via f′(u) = (f(u) − f(x))′ yields

EKKLD0+,βγ,q=qβxβ(1+γ)Γ(1q)0x(f(u)f(x))(xβuβ)q1uβ(γ+q+1)uβ1duβxβ(1+γ)Γ(1q)0x(f(u)f(x))(xβuβ)quβ(γ+q+1)1(γ+q+1)du+βxβ(1+γ)Γ(1q)0xf(u)(xβuβ)quβ(γ+q+1)1uβ1(γ+1)du
=qβxβ(1+γ)Γ(1q)0x(f(u)f(x))(xβuβ)q1uβ(γ+q+1)uβ1duqβxβ(1+γ)Γ(1q)0x(f(u)f(x))(xβuβ)quβ(γ+q+1)1du+β(γ+1)f(x)xβ(1+γ)Γ(1q)0x(xβuβ)quβ(γ+q+1)1du.

Performing cancellations and using the integral

0x(xβuβ)quβ(γ+q+1)1du=1βxβ(1+γ)B(1q,γ+q+1)=Γ(1q)Γ(1+q+γ)βΓ(2+γ)xβ(1+γ),

yields (3.20). □

The process generated by (3.20) (with the negative sign) can be called the Erdélyi-Kober-Kiryakova-Luchko process. This is a decreasing process on [0, ∞) with the constant killing rate Γ(γ + q + 1)/Γ(γ).

Of course one can obtain formula for the interrupted version of this process generated by the fractional derivatives EKKLDa+,βγ,qandEKKLDa+,βγ,q with a > 0. We shall do it below in a more general setting. The lower bound a = 0 is specific both for yielding the commuting class of operators and by the fact that the underlying processes have uniform killing rates.

### 3.3 Elementary analytic approaches to generalized fractional operators

The development of the generalized fractional calculus was initially led and motivated by the theory of special functions, namely the class of the so-called G-functions and more general H-functions. This development is well presented in the literature (see e.g. [30], [31], [53], [33], [63], [34] [57], [34], etc.)

We introduce here two elementary approaches both leading rapidly to the wide class of operators that cover the majority of operators discussed in the literature. Then we derive the probabilistic representation for these operators revealing the stochastic processes that govern the solutions of the corresponding PDEs.

1. Method of iterations. Extending the approach leading to the definition (2.1) of the RL integral let us look at a more general integral operator of the Riemann-Stieltjes type

Ia,Gf(x)=axf(u)dG(u),(3.22)

where G(u) is a non decreasing function. Assuming that G is differentiable, G′(u) = g(u), it simplifies to Ia,Gf(x) = axg(u) f(u) du. By direct induction and integration by parts one gets the formula for the iterations:

Ia,Gnf(x)=1(n1)!ax(G(x)G(t))n1g(t)f(t)dt,(3.23)

for differentiable G, and

Ia,Gnf(x)=1(n1)!ax(G(x)G(t))n1f(t)dG(t),(3.24)

in general case. In full analogy with the definition (2.1) formula suggests the definition of the fractional RL-Stieltjes integral by the formula

Ia,Gqf(x)=1Γ(q)ax(G(x)G(t))q1f(t)dG(t)(3.25)

for any q > 0.

If we take g(u) = 1/u and a > 0, formula (3.25) turns to formula (3.1) defining the Hadamard integral.

In [25], Katugampola aimed at obtaining some interpolation between the standard RL integral and the Hadamard integral and for this purpose applied the above scheme of reasoning with the power function g(u) = sρ, ρ > −1, which yields the integral

Ia,ρqf(x)=(ρ+1)1qΓ(q)ax(xρ+1tρ+1)q1tρf(t)dt.(3.26)

Katugampola was seemingly unaware of the Erdélyi-Kober operators’ theory and did not recognise the link with the classical Erdélyi-Kober integral (3.15) with γ = 0:

EKKLIa,β0,qg(x)=xβqβqIa,β+1qf(x).(3.27)

The difference in the so-called (new) Katugampola operator is therefore just in a multiplier that does not seem to be very essential from the first sight. However, as was noted above, with precisely this multiplier the Erdélyi-Kober-Kiryakova-Luchko integrals with various indices become commuting operators, which creates the beautiful link with the theory of H-functions and leads to the multi-indices Mittag-Leffler functions and the multi-indices Erdélyi-Kober integrals (see e.g. [33]).

To extend this commuting class of operators, it was suggested by Kalla [24] (see also references in [33]) to analyse the operators of the type

I0Φg(x)=01Φ(u)g(ux)du,(3.28)

parametrized by functions Φ (can be arbitrary special functions). However, by the direct change of the variables u = eξ, x = ey, these operators turn just to the general shift invariant integral operators

f0f(yξ)K(ξ)dξ,K(ξ)=Φ(eξ)eξ.

Thus by the universal change of variables the whole class of Erdélyi-Kober-Kiryakova-Lucko operators is inserted in the class of such shift invariant operators. However, these concerns only the operators with the lower bound a = 0.

Operators (3.22) were applied in [68], [69] for solving certain ordinary ODEs via the ’fractional tools’.

2. Substitution and dressing. Assuming G is strictly monotone and continuous and rewriting integral (3.22) via the method of substitution,

Ia,Gf(x)=G(a)G(x)f(u(G))dG,(3.29)

with u(G) the inverse function of G, suggests that the fractional RL-Stieltjes integral (3.25) can be obtained directly from the corresponding standard RL integral just by the change of variable. Namely, denoting by Ω the operator acting as ΩGf(x) = f(G(x)) and fixing the function G by the condition G(a) = a, one sees that Ia,G = ΩGIaΩG1 and

Ia,Gq=ΩGIaqΩG1,(3.30)

that is, Ia,G and Ia,Gq (defined by (3.25)) are obtained from the corresponding standard integrals via ’dressing’ with the operator Ω of the change of variable.

Once the idea of dressing is introduced, it is natural to exploit its further by ’dressing’ with respect to other linear transformations of functions, the simplest one being the operators of multiplication Mwf(x) = w(x) f(x) with a function w. Using such dressing in combination with (3.30) leads to the following generalized dressed fractional integral operator

Ia,G,wq=Mw1ΩGIaqΩG1Mw,(3.31)

or concretely,

Ia,G,wqf(x)=1Γ(q)w(x)ax(G(x)G(t))q1f(t)w(t)dG(t).(3.32)

Assuming G to be differentiable, it becomes

Ia,G,wqf(x)=1Γ(q)w(x)ax(G(x)G(t))q1f(t)w(t)G(t)dt.(3.33)

Formula (3.31) suggests the definition of the corresponding derivatives as the dressed versions of the standard fractional derivatives. Namely, working for simplicity with differentiable G and noting that

ΩddxΩ1f(x)=1G(x)f(x),

we obtain, for the case for q ∈ (0, 1), the following RL and CD dressed derivatives:

Notation AD reflects the fact that operators (3.33), (3.34), (3.35) were seemingly first introduced and studied by Agrawal in [3].

### Remark 3.1

The authors of [28] introduced and studied the modified Hadamard integrals, which are seen to be obtained from the usual Hadamard integrals via dressing with the operator Mw with w(x) = x.

As an alternative to (3.31), one can use the dressing by multiplications and change of variables in a different order, and also use different functions in the input and output yielding the operators

Ia,G,w,vq=ΩGMvIaqMwΩG1,(3.36)

the corresponding derivatives being

Da+q;G,w,v=ΩGMvDa+qMwΩG1,(3.37)
Da+q;G,w,v=ΩGMvDa+qMwΩG1.(3.38)

The Erdélyi-Kober-Kiryakova-Luchko fractional operators are obtained via these formulas with all functions G, w, v being power functions.

### 3.4 Probabilistic representations

Of course, differential equations involving dressed operators can be reduced to the usual fractional equations via dressing and redressing. This method was developed in detail in [7] (where it was called the transmutation method) for the case of Erdélyi-Kober operators. However, if one is interested in mixtures of the derivatives of different kind (say, dressed by different operators), one needs to have explicit representation for each mixed component.

### Proposition 3.4

IfGis differentiable increasing function withG(a) = aandfis continuously differentiable, then

+f(x)Γ(1q)1(G(x)G(a))qqaxw(u)w(x)w(x)(G(x)G(u))1+qG(u)du,(3.39)
+f(x)Γ(1q)w(x)w(a)w(x)(G(x)G(a))qqaxw(u)w(x)w(x)(G(x)G(u))1+qG(u)du,(3.40)

with

νA(x,s)=1Γ(q)w(s)w(x)(G(x)G(s))1+qG(s).(3.41)

### Proof

It is most handy to use the first formula in (3.34) and the formula for the standard derivatives (2.11) implying that

=1Γ(q)w(x)aG(x)(wf)(G1(y))(wf)(x)(G(x)y)1+qdy+f(x)Γ(1q)(G(x)a)q.

Changing variable of integration yields

=1Γ(q)ax(w(u)/w(x))f(u)f(x)(G(x)G(u))1+qG(u)du+f(x)Γ(1q)(G(x)a)q,

which can be rewritten as

as yielding (3.39).

Similarly,

ΩGDa+qΩG1f(x)=1Γ(q)ax(f(u)f(x))G(u)du(G(x)G(u))1+q+f(x)f(a)Γ(1q)(G(x)a)q.

Dressing by Mw and rearranging yields (3.40).□

If the function w(x) has a constant sign everywhere and the coefficient at f(x) is positive (for instance, if w(x) is an increasing positive function), the operator ADa+,G,wq generates a decreasing Markov process with killing that can be called Agrawal’s process.

As a hallmark of the dressing with the multiplication operators, the CD derivative also gets a killing term (with the rate given by the last square bracket in (3.40)) and the equation ADa+qf=ADa+q(ff(a)) does not hold anymore. However, the RL and CD derivatives still coincide if f(a) = 0.

Needless to stress that these operators also fit the general framework of the operators (2.28).

Let us turn to the derivatives (3.37), (3.38) representing the extension of the Erdélyi-Kober-Kiryakova-Luchko operators.

### Proposition 3.5

IfGis differentiable increasing function withG(a) = aandfis continuously differentiable, then

Da+q;G,w,v=ax(f(s)f(x))νG,w,v(x,s)ds+f(x)(ΩGMvDa+qw)(x),(3.42)
Da+q;G,w,v=ax(f(s)f(x))νG,w,v(x,s)ds+f(x)(ΩGMvDa+qw)(x)+(f(x)f(a))w(x)v(G(x))Γ(1q)(G(x)a)q,(3.43)

with

νG,w,v(x,s)=v(G(x))w(G(s))G(s)Γ(q)(G(x)G(s))1+q(3.44)

### Proof

By (2.11),

MvDa+qMwΩG1f(x)=v(x)Γ(q)axw(y)f(G1(y))w(x)f(G1(x))(xy)1+q+f(G1(x))w(x)v(x)Γ(1q)(xa)q.

Therefore,

Da+q;G,w,v=v(G(x))Γ(q)aG(x)w(y)f(G1(y))w(G(x))f(x)(G(x)y)1+qdy+f(x)w(G(x))v(G(x))Γ(1q)(G(x)a)q=v(G(x))Γ(q)axw(G(u))f(u)w(G(x))f(x)(G(x)G(u))1+qG(u)du+f(x)w(G(x))v(G(x))Γ(1q)(G(x)a)q,

and consequently,

Da+q;G,w,v=v(G(x))Γ(q)axw(G(u))(f(u)f(x))(G(x)G(u))1+qG(u)du+f(x)v(G(x))Γ(q)axw(G(u))w(G(x))(G(x)G(u))1+qG(u)du+f(x)w(G(x))v(G(x))Γ(1q)(G(x)a)q.

The expression in the square bracket rewrites as

v(G(x))Γ(q)aG(x)w(y)w(G(x))(G(x)y)1+qdy+f(x)w(G(x))v(G(x))Γ(1q)(G(x)a)q

yielding (3.42).

Similarly,

Da+q;G,w,v=v(G(x))Γ(q)aG(x)w(y)f(G1(y))w(G(x))f(x)(G(x)y)1+qdy+(f(x)w(G(x))f(a)w(a))v(G(x))Γ(1q)(G(x)a)q,

yielding (3.43).□

As was noted above, dressing fractional derivatives by the multiplication operators leads to killing even in the CD version. Therefore it is natural to modify the obtained CD-type derivatives by subtracting the killing term. In fact, performing thus procedure with the dressing by the exponential function e–λx leads to the tempered fractional derivatives defined by the Lévy kernel (2.21) (see [6] and [58]).

### 3.5 Multiplicative interpolations of Hilfer

Let us touch upon the multiplicative interpolations between (or multiplicative mixtures of) the CD and RL derivatives introduced by Hilfer (see e.g. [23]), which are defined by the formulas

HiDa+α,β=I(1α)βDIa+(1α)(1β),β(0,1),(3.45)

so that for β = 0 this operator turns to the RL derivative and for β = 1 to the CD one.

From the point of view of probability, these derivatives do not lead to any new processes, as the following statement shows.

### Proposition 3.6

Iffis continuously differentiable andβ < 1, then

HiDa+α,βf(x)=Da+αf(x),(3.46)

that is, for smoothf, the interpolated derivatives turn to the usual RL derivatives.

### Proof

Using the definition and integration by parts we get

Ia+(1α)(1β)f(x)=1Γ[(1α)(1β)]ax(xt)(1α)(1β)1f(t)dt=1Γ[(1α)(1β)]axf(t)ddt(xt)(1α)(1β)(1α)(1β)dt
=1Γ[(1α)(1β)]f(a)(xa)(1α)(1β)(1α)(1β)+(Ia+(1α)(1β)+1f)(x).

Therefore,

DIa+(1α)(1β)f(x)=(xa)(1α)(1β)1Γ[(1α)(1β)]f(a)+(Ia+(1α)(1β)f)(x),

and

I(1α)βDIa+(1α)(1β)f(x)=(Ia+(1α)f)(x)
+f(a)Γ[(1α)β]1Γ[(1α)(1β)]ax(xs)(1α)β1(sa)(1α)(1β)1ds
=Da+αf(x)+f(a)Γ[(1α)β]1Γ[(1α)(1β)](xa)αB[(1α)β,(1α)(1β)]
=Da+αf(x)+f(a)Γ(1α)(xa)α=Da+αf(x).

For completeness, let us note that the working with derivatives (3.45) allows one to distinguish succinctly the scales of the regularity classes of the solutions to the spectral problems:

HiDa+α,βf=Ia+(1α)βDIa+(1α)(1β)f=λf.(3.47)

In fact, applying the operator Ia(1α)(1β) to the both sides of this equation leads to the spectral problem for the CD operator:

Da+αg=λg,g=Ia+(1α)(1β)f,(3.48)

and applying further the usual derivative D leads to to the spectral problem for the RL operator:

Da+αh=λh,h=g=Da+1(1α)(1β)f.(3.49)

The natural boundary value condition g(a) = Y for the CD problem (3.48) leads to the integral boundary condition Ia+(1α)(1β)f(a)=Y for the Hilfer spectral problem for f. In particular, this condition turns to

Ia+(1α)f(a)=Y(3.50)

for the RL spectral problem.

### 3.6 Random fractional integrals

Iterated integrals (3.22) and the corresponding fractional integrals were considered so far only for the case of differentiable G(x). The general approach presented here suggests further extensions. Firstly one can look at non-smooth increasing G, for instance, choosing as G the famous Cantor staircase. Aiming to assess the general features of such integrals, one can look at random fractional integrals arising from choosing G(x) as an increasing process, for instance, a Lévy subordinator. Even more exotic calculus can be introduced attempting to build fractional versions of more general indefinite stochastic integrals. For instance, one can start with the Wiener or Ito integral

IWf(x)=0xf(y)dB(y),

where dB(y) is the Ito stochastic differential of the standard Brownian motion. Simple (but lengthy) calculations show that if f is a smooth function (or even a smooth adapted processes on the Wiener space) vanishing at zero, the iterated integral can be expressed as the following usual (non-stochastic) integral:

(IW)kf(x)=1k!0xf(y)(xy)k/2HkB(x)B(y)xydy,(3.51)

where

Hk(x)=(1)kex2/2dkdxkex2/2=xddxk1,kN,

are the Hermite polynomials. Extending Hk(x) to fractional values of k (by using either fractional derivative in the first equality or the fractional power in the second equality of the last formula) would naturally define the fractional Wiener integral. We shall not explore this topic further here.

### Remark 3.2

In the theory of probability by fractional Wiener integral one understands another object, namely the integral 0xf(y) dBH(y) with BH the fractional Brownian motion.

## 4 Simplest linear equations and generalized Mittag-Leffler functions

### 4.1 Equations with shift invariant mixed derivatives

If a non-negative measure ν on {y : y > 0} satisfy (2.19) and gCconst(–∞, a] (R) ∩ Cuc(R), the resolvent operator Rλ from (2.37) yields the unique solution Rλg of the equation

Da+(ν)f(x)=λf(x)+g(x)(4.1)

in the domain of the generator of the semigroup Tt on Cuc(R). This function equals g(a)/λ to the left of a.

However, Rλg is not the solution we are mostly interested in, as it prescribes the boundary value at a, rather than solves the boundary value problem. By prescribing the boundary value one necessarily takes the solution out of the domain of the generator.

The natural way to state properly the boundary value problem

Da+(ν)f(x)=λf(x)+g(x),f(a)=Y,xa,(4.2)

is by turning it to the problem with the vanishing boundary value, which is a usual trick in the theory of PDEs. Namely, introducing the new unknown function u = fY we see that u must solve the problem

Da+(ν)u(x)=λu(x)λY+g(x),u(a)=0,xa,(4.3)

just with g – λY instead of g. We can thus define the solution to(4.2) to be the function f = u + Y, where u solves (4.3). Let us stress for clarity that in (4.3) the r.h.s. g(x) – λY is considered to be continued as zero to the left of a.

### Remark 4.1

It is also possible to work directly with the resolvents of the CD derivatives (see [20], [21] ), but the approach via equation (4.3) seems to be simpler.

Taking first g = 0 we find the solution to (4.2) to be

f(x)=Y+u(x)=YλY0xa0eλtG(ν)(t,dy)dt=λY0eλtxaG(ν)(t,dy)dt,(4.4)

since

λ00eλtG(ν)(t,dy)dt=1.

Integrating by parts and taking into account that xaG(ν)(0, dy) = 0 for x > a, we get in this case the alternative expression:

f(x)=Y0eλttxaG(ν)(t,dy)dt.

Restoring g we arrive at the following.

### Proposition 4.1

For anygsupported on [a, ∞) the unique solution to problem(4.2)in the sense defined above is given by the formula

f(x)=Y0eλttxaG(ν)(t,dy)dt
+0xag(xy)0eλtG(ν)(t,dy)dt.(4.5)

This solution can be classified as classical (from the domain of the generator) or generalized (in the sense of the generalized functions or by approximation) according to Proposition 2.2 applied to problem(4.3).

In analogy with the derivative dβ/dxβ, one can define the family of generalized Mittag-Leffler functions depending on the positive parameter z as

E(ν),z(λ)=0eλttzG(ν)(t,dy)dt=1λ0eλt0zG(ν)(t,dy)dt=1λ0zUλ(ν)(dy).(4.6)

Since the function xaG(ν)(t, dy) increases with t, its derivative (in the sense of generalized functions) is well-defined as a positive measure (and as a function almost everywhere), and therefore the function E(ν),z(–λ) are well defined and continuous for λ ≥ 0. They are completely monotone function of these λ and are bounded by 1:

|E(ν),z(λ)|0tzG(ν)(t,dy)dt=zG(ν)(t,dy)0=1.(4.7)

Moreover, E(ν),z(0) = 1 and the solution (4.4) to problem (4.3) is expressed as

f(x)=YE(ν),xa(λ)+0xag(xy)Uλ(ν)(dy),(4.8)

where the λ-potential measure is expressed in terms of E(ν),z by the equation

0zUλ(ν)(dy)=(1E(ν),z(λ))/λ.(4.9)

If the measures G(ν)(t, dy) have densities with respect to Lebesgue measure, G(ν)(t, y), then the λ-potential measure also has a density, Uλ(ν)(y), and (4.9) rewrites as

Uλ(ν)(y)=1λE(ν),y(λ)y.(4.10)

However, only for the case of the derivative dβ/dxβ, due to the particular scaling property of Gβ, one has the additional relation E(ν),z(–λ) = Eβ(–λzβ).

In case of Lν = –dβ/dxβ,

1Gβ(t,y)dy=1t1/βGβ(1,t1/βy)dy=t1/βGβ(1,x)dx.

so that

x1Gβ(t,y)dy=1βt11/βGβ(1,t1/β),

and therefore we again arrive at formula (2.42).

Let us now turn to the extension of the linear equations to the Banach-space-valued setting, that is, to the equations

Da+(ν)μ(x)=Aμ(x)+g(x),μ(a)=Y,xa.(4.11)

If μ(a) = Y = 0, this turns to the RL type equation

Da+(ν)μ(x)=Aμ(x)+g(x),μ(a)=0,xa.(4.12)

As above, we shall define the solution to (4.11) as the function μ(x) = Y + u(x), where u(x) solves the problem

Da+(ν)u(x)=Au(x)+AY+g(x),u(a)=0,xa.(4.13)

Notice also that the assumption of etA to be a contraction naturally extends the case A = –λ with λ > 0, as e–λt ≤ 1, and allows one to define the operator-valued generalized Mittag-Leffler functions by the operator-valued integral

E(ν),z(A)=0etAtzG(ν)(t,dy)dt=1+A0etA0zG(ν)(t,dy)dt.(4.14)

### Theorem 4.1

1. Let the measureνon {y : y > 0} satisfy(2.19)and letAbe the generator of the strongly continuous semigroupetAof contractions in the Banach spaceB, with the domain of the generatorDB. Then the 𝓛(B, B)-valued potential measure,

UA(ν)(M)=0etAG(ν)(t,M)dt,(4.15)

of the semigroupTtetAon the subspaceCkill(a)([a, b], B) ofCuc((–∞, b], B) is well-defined as aσ-finite measure on {y : y ≤ 0} such that for anyz > 0, λ > 0,

UA(ν)([0,z])eλz/ϕν(λ).

Therefore, the potential operator (given by convolution withUA(ν)) of the semigroupTtetAonCkill(a)k([a, b], B) is bounded for anyb > a.

2. For any gCkill(a)([a, b], B), theB-valued function

f(x)=0xaUA(ν)(dy)g(xy)=0xa0etAG(ν)(t,dy)dtg(xy)(4.16)

belongs to the domain of the generator of the semigroupTtetAand represents the unique solution to problem(4.12)from the domain. For anygC([a, b], B) (continued as zero to the left ofa), this function represents the unique generalized solution to(4.12), both by approximation and in the sense of generalized functions.

3. For anygC([a, b], B) (continued as zero to the left ofa) andYB, the function

f(x)=Y+0xaUA(ν)(dy)(AY+g(xy))=E(ν),xa(A)Y+0xaUA(ν)(dy)g(xy)(4.17)

represents the unique generalized solution to problem(4.11)or(4.13).

### Proof

(i) The norms of the operator-valued measure UA(ν) are estimated by the measure U(ν), because etA are contractions. (ii) and (iii) follow in the same way as for real-valued A above.□

### 4.2 Equations with weighted mixed derivatives

Let us turn to equations arising from the operator L(ν,γ)l from (2.22) and the corresponding derivative (2.26), assuming everywhere that the assumptions of Proposition 2.5 hold.

The linear problem

Da+(ν,γ)f(x)=λf(x)+g(x),f(a)=Y,xa,(4.18)

will be dealt with in the same way as (4.2). Namely, its solution will be defined as the function f = u + Y, where u solves the equation

Da+(ν,γ)u(x)=λu(x)λY+g(x),u(a)=0,xa.(4.19)

Taking first g = 0 we find (by (2.51) for λ = 0 and by the definition of resolvent (2.45) for λ > 0) that the solution to (4.2) equals

f(x)=Y+u(x)=YλYax0eλtP(ν,γ)(t,x,dy)dt,(4.20)

where P(ν,γ)(t, x, dy) are the transition probabilities for the process generated by L(ν,γ)l. Consequently, we get for x > a that

f(x)=λY0eλtaP(ν,γ)(t,x,dy)dt=Y0eλttaP(ν,γ)(t,x,dy)dt.(4.21)

As in the case of shift invariant mixtures, the expression under the bracket is monotonic in t and hence the derivative is well defined almost everywhere and is positive. Of course, strictly speaking, the last formula holds only in case of absolutely continuous (in t) function aP(ν,γ)(t,x,dy). In the general case, the last formula should be written more precisely as the Stieltjes integral

f(x)=Y0eλtdtaP(ν,γ)(t,x,dy),(4.22)

where dt is the differential with respect to the variable t.

Restoring g, we arrive at the following.

### Proposition 4.2

Let the assumptions of Proposition 2.5 hold. Then for any λ ≥ 0 and anygsupported on [a, ∞) the unique solution to problem(4.18)in the sense defined above is given by the formula

f(x)=Y0eλtdtaP(ν,γ)(t,x,dy)+axg(y)Π(ν,γ)λ(x,dy).(4.23)

We can now define the family of generalized Mittag-Leffler functions as

Ex,a(ν,γ)(λ)=0eλtdtaP(ν,γ)(t,x,dy)=1λΠ(ν,γ)λ(x,[a,x]).(4.24)

In explicitly probabilistic form it can be written also as

Ex,a(ν,γ)(λ)=λYE0eλt1ya(Xxt)dt,(4.25)

where Xxt is the Feller process generated by L(ν,γ)l (a simpler expression will be given below in (4.38)). These functions are completely monotone. Moreover, it follows that the differential of E with respect to a exists as a measure so that

Π(ν,γ)λ(x,dy)=dyEx,y(ν,γ)(λ)

and the solution (4.23) rewrites as

f(x)=YEx,a(ν,γ)(λ)+axg(y)dyEx,y(ν,γ)(λ).(4.26)

Similarly the operator-valued equation

Da+(ν)μ(x)=Aμ(x)+g(x),μ(a)=Y,xa(4.27)

is analyzed for the case of A generating a semigroup of contractions, etA, in a Banach space B.

The operator-valued generalized Mittag-Leffler functions are defined by the operator-valued integral

Ex,a(ν,γ)(A)=0eAtdtaP(ν,γ)(t,x,dy)=1+AΠ(ν,γ)A(x,[a,x]),(4.28)

where Π(ν,γ)A is the operator-valued potential measure,

Π(ν,γ)A(x,dy)=0eAtdtP(ν,γ)(t,x,dy),(4.29)

and the following direct extension of Theorem 4.1 is obtained.

### Theorem 4.2

Let the assumptions of Proposition 2.5 hold and letAbe the generator of the strongly continuous semigroupetAof contractions in the Banach spaceB. Then the 𝓛(B, B)-valued potential measure(4.29)is well defined. Moreover, for anygC([a, b], B) (continued as zero to the left ofa) andYB, the function

f(x)=Ex,a(ν,γ)(λ)Y+ax[dyEx,y(ν,γ)(λ)]g(y),(4.30)

where in the square bracket sits the operator-valued Stiltjes measure, represents the unique generalized solution to problem(4.27).

### 4.3 Power series expansion for generalized Mittag-Leffler functions

We have constructed the solutions to the linear problems (4.12) and (4.11) only for the case of A generating a contraction semigroup (with a direct extension to the case of a uniformly bounded semigroup etA). This restriction was ultimately linked with formula (4.24) for the generalized Mittag-Leffler function, from which it is not seen directly that it can be extended to negative λ. We shall see that this is always possible providing by-passing the power series representation for these functions.

From the definition (4.24) it follows that

Ex,a(ν,γ)(λ)=1λ(Rλ1a)(x)

with the resolvent operator (2.45). But in the space Ckill(–∞,a]((–∞, b]) for any b > a, the potential operator Ia+(ν,γ) is the bounded right inverse to D+(ν,γ), so that the resolvent equals

Rλ=(λ+D+(ν,γ))1=Ia+(ν,γ)(1+λIa+(ν,γ))1=k=0(λ)k(Ia+(ν,γ))k+1,

the series being convergent at least for λIa+(ν,γ)b < 1, where Ia+(ν,γ)b is the norm of the operator Ia+(ν,γ) on the space Ckill(–∞,a]((–∞, b]). Therefore,

Ex,a(ν,γ)(λ)=k=0(λIa+(ν,γ))k1a(x).(4.31)

This yields the series representation for the generalized Mittag-leffler functions. Since the series has a non-vanishing convergence radius (for any finite x), it defines the continuation of the Mittag-leffler functions Ex,a(ν,γ)(–λ) to a neighborhood of zero.

In order for this series to be convergent for all λ ∈ C one has to impose some additional assumptions on the kernel ν, the simplest one being a lower bound of a stable type.

### Proposition 4.3

Under the assumptions of Proposition 2.5 letν (x, dy) have the lower bound of theβ-fractional type

ν(x,dy)(1/Γ(β))Cνy1βdy(4.32)

with someβ ∈ (0, 1) andCν > 0. Then, for anyx > a,

Π(ν,γ)(x,[a,x])CνUβ[0,xa]=Cν(I0β1)(xa)=Cν(xa)β/Γ(β),(4.33)
Π(ν,γ)λ(x,[a,x])CνUλβ[0,xa]Cν|λ|[Eβ(|λ|(xa)β)1],(4.34)
Ex,aν,b(λ)max(1,Cν)Eβ(|λ|(xa)β).(4.35)

### Proof

All inequalities follow from the comparison principle of Proposition 2.5 and equation (4.24).□

Proposition 4.3 implies that the function Ex,aν,b(λ) is well defined as a whole analytic function of λ and series (4.31) converges for all λ. This allows one to get a more or less direct extension of Theorem 4.1 to arbitrary strongly continuous semigroups etA (not necessarily contractions), see detail in [44].

### 4.4 Probabilistic solutions of linear equations and probabilistic representation for generalized Mittag-Leffler functions

Solutions to linear fractional equations constructed above are expressed in terms of the transition probabilities of underlying processes. The derivation was performed analytically, via resolvents. It is instructive to give an alternative, pure probabilistic derivation of these results. this approach leads also to the new representations for the solutions.

The main tool here is the Dynkin martingale. Namely, as is well known in probability theory, if Xt is a Feller process Xtx in Rd generated by the operator L and f is a function from the domain of L, then the process

Mtf,λ=f(Xtx)eλt+0teλs(λL)f(Xs)ds

is a martingale for any λ ≥ 0, called the Dynkin martingale. Let τ be a stopping time with a finite expectation Eτ. Then one can use Doob’ optional sampling theorem to deduce the following Dynkin formula:

f(x)=E[f(Xτx)eλτ+0τeλs(λL)f(Xsx)ds].(4.36)

Suppose now that f is a solution for the problem (4.18) with Y = 0 (chosen for simplicity, as the general case is reduced to this anyway). Recalling that the fractional derivative is the negation of the generator of the corresponding Markov process and choosing the stopping time τx,a to be the time of exit of the process from the interval (a, ∞) we find that

f(x)=E0τx,aeλsg(Xsx)ds.(4.37)

It is of course not surprising that this formula coincides with the expression (4.23) with Y = 0. In fact, recalling that g is considered to be continued as zero to the left of a, we have

E0τx,aeλtg(Xtx)dt=0eλtg(Xtx)dt=0eλtaxP(ν,γ)(t,x,dy)g(y)ds.

Thus (4.37) yields the probabilistic representation for the solution (4.23) with Y = 0.

Next, we get from (4.24) (or from (4.20)) that

Ex,a(ν,γ)(λ)=1λ0eλtaxP(ν,γ)(t,x,dy)ds=1λE0τx,aeλtdt=1λE1eλτx,aλ=Eeλτx,a,

yielding the fundamental probabilistic representation for the generalized Mittag-Leffler function:

Ex,a(ν,γ)(λ)=Eeλτx,a.(4.38)

Thus this function is just the Laplace transform of the exit time of the underlying decreasing Markov process. For instance, for the shift invariant case, the function (4.6) equals

E(ν),x(λ)=Eexp{λτx(ν)},(4.39)

where τx(ν) is the exit time from the interval [0, x) for the Lévy subordinator generated by the operator Lν from (2.20). For the classical Mittag-Leffler function formula (4.39) reduces to the known formula

Eβ(λxβ)=Eexp{λτxβ}(4.40)

(seemingly first derived in [11], see also [60]), where τxβ is the corresponding exit time for the β-stable subordinator.

### Remark 4.2

The remarkable formula (4.40) and its direct application to the Banach space valued problems seem to be poorly appreciated by the ’fractional community’.

The solution to the spectral problem (4.11) with g(x) = 0 can now be written in the following remarkably simple form

μ(x)=Eexp{Aτx(ν)}Y,(4.41)

that explicitly unifies the solutions for all mixed fractional derivatives of order up to and including 1.

As another example let us consider the Cauchy problem

HDa+qf(x)=Af(x),f(a)=Y,(4.42)

in the Banach space B, with the Hadamard derivative (3.41) and A as above. Its solution writes down as

μ(x)=Eexp{Aτx,a}Y,(4.43)

where τx,a is the exit time from the interval (a, x] for the Hadamard process generated by (3.8) and started at x.

### 4.5 Some examples

Let us present some examples.

1. Generalized fractional Schrödinger equation:

Da+(ν,γ)ψ=iHψ,(4.44)

where H is a self-adjoint operator in a Hilbert space 𝓗. Since H generates a unitary group, and hence a contraction semigroup, Theorem 4.2 applies. Similarly one can deal with fractional Schrödinger equation with the complex parameter

Da+(ν,γ)ψ=σHψ,(4.45)

if H is a negative self-adjoint operator. Specific examples of these equations were analyzed recently in [17].

2. Generalized fractional Feller evolution, where the operator A generates a Feller semigroup in C(Rd) and a Feller process, for instance, a diffusion or a stable-like process.

3. Generalized fractional evolutions generated by ΨDOs with spatially homogeneous symbols (or with constant coefficients):

Da+(ν,γ)f=ψ(i)f+g,f|t=a=fa,(4.46)

under various assumptions on symbols ψ(p) ensuring that –ψ(–i∇) generates a semigroup. In this case the solution given by Theorem 4.2 can be constructed explicitly via the Fourier transform, see detail in [44].

## 5 Further linear equations

### 5.1 Two-sided problems

We shall now touch upon the theory of two-sided problems that includes the problem

Da+(ν)f+Db(ν)f=λf+g,f|t=a=fa,f|t=b=fb.(5.1)

In its general form it is the problem of the type

D[a,b](ν)f=Lf=λf+g,f|t=a=fa,f|t=b=fb,(5.2)

with L from (2.56). It turns out that essential simplification for the two-sided problems occurs in case of the Lévy measures with densities. We shall consider only this case (also omitting for simplicity the usual drift) thus choosing

Lf(x)=axbx(f(x+y)f(x))ν(x,y)dy+(f(b)f(x))bxν(x,y)dy+(f(a)f(x))axν(x,y)dy.(5.3)

As follows from the discussions of the previous section, by shifting the unknown function, the problem can be reduced to the problem with the vanishing boundary conditions. To solve the latter problem one has to construct the Feller process and the resolvent generated by the operator

D[a,b](ν)f(x)=Lkillf(x)=axbx(f(x+y)f(x))ν(x,y)dyk(x)f(x),(5.4)

with

k(x)=ka(x)+kb(x)=ax+bxν(x,y)dy,

describing the process killed on the boundary. As was mentioned already, these semigroup and the resolvent cannot be obtained from the corresponding objects for the operator without boundary just by reducing it to the space of functions vanishing outside (a, b).

We shall denote by primes the derivatives with respect to the variable x. For instance, ν′(x, y) = (/∂x)ν(x, y).

The following conditions on ν will be assumed throughout this section:

1. ν(x, y) is a continuous function of two variables on the set (x ∈ [a, b], y ≠ 0) having a continuous derivative ν′(x, y) such that

zν(x,z)ϰzν(x,y)dy(5.5)

with a constant ϰ < 1 and sufficiently small z;

2. ν(x, y) ≤ ν̃(y), ν′(x, y) ≤ ν̃(y) with a function ν̃: ∫ (1 ∧ y)ν̃(y)dy < ∞;

3. k(x) → ∞ as xa or xb.

Note that (C) is just the assumption that the Lévy kernel ν(x, y)dy is unbounded. Assumption (A) is not too restrictive, at least it holds for all standard examples. Say, for α-stable processes (classical fractional derivatives) it holds with ϰ = α.

Let us introduce the special notations for our main Banach spaces: C0 = Ckill{a, b}([a, b]) equipped with the sup-norm ∥.∥, C00 = {fC0 : f′ ∈ C0} equipped with the norm ∥f00 = ∥f∥ + ∥f′∥.

### Lemma 5.1

Under conditions (A)-(C) the operatorK of multiplication byk(x) from(5.15)generates: (i) a semigroup of contractions inC([a, b]), (ii) a strongly continuous semigroup of contractions inC0with the domain

DK={fC0:limxak(x)f(x)=0,limxbk(x)f(x)=0},

(iii) a uniformly bounded semigroup inC0C1([a, b]), (iv) a uniformly bounded strongly continuous semigroup inC00.

### Proof

It is mostly straightforward. Let us prove only (iii) and (iv). We have

[etk(x)f(x)]=etk(x)f(x)tk(x)etk(x)f(x).

The first term is bounded: ∥etk(.)f′(.)∥ ≤ ∥f′∥. The second term is bounded away from the boundaries. Let us estimate it in a neighborhood of x = b (neighborhoods of a are analogous). Near x = b the main (unbounded) part of k(x) is kb(x) = bxν (x, y)dy and thus

tk(x)etk(x)f(x)tν(x,bx)(bx)f(b)etk(x)+t(bx)f(b)etk(x)bxν(x,y)dy.

For bx ≤ 1 we estimate the two terms by (A) and (B) respectively:

tν(x,bx)(bx)|f(b)|etk(x)ϰ|f(b)|tkb(x)etkb(x)ϰ|f(b)|,
t(bx)|f(b)|etk(x)bxν(x,y)dyt|f(b)|0(1y)ν(x,y)dyt|f(b)|0(1y)ν~(y)dy.

Thus both terms are uniformly bounded. Moreover, if f′(b) = 0, then limxb[etk(x)f(x)]′ = 0, so that etK : C0C1([a, b]) → C00 for t > 0. This implies the strong continuity of etK in C00.□

### Theorem 5.1

Assume that (A)-(C). Then the operatorLkillof(5.15)generates a Feller semigroup inC0and a bounded semigroup inC00. The domain of this semigroup inC0contains the space

Dkill={fC00:Lkillf(x)0,asxa,xb}.

### Proof

1. Let us introduce the approximation Lh of Lkill, h ∈ (0, 1), by the formula

Lhf(x)=axbx(f(x+y)f(x))(1χh(y))ν(x,y)dyk(x)f(x),(5.6)

where χh a continuous even function R → [0, 1] such that χh(z) = 1 for z ∈ [–h, h] and χh(z) = 0 for |z| ≥ 2h. Since Lh differs from the operator –K by a bounded operator (in both C([a, b]) and C1([a, b])), we can conclude from Lemma 5.1 and the perturbation theory that Lh generates bounded semigroups Tth in the spaces C([a, b]) and C1([a, b]). Moreover, the perturbation series representation for the semigroup Tth has the form (see e.g. [40])

Tth=etK+m=10s1smtds1dsm×e(tsm)K(Lh+K)e(smsm1)K(Lh+K)es1K.

It follows from this formula and Lemma 5.1 that the spaces C0 and C00 are invariant under Tth. It is also straightforward to see that Tth is strongly continuous in both C0 and C00. Since the operator Lh is conditionally positive, it follows that Tth is a contraction in C0. The domain of the semigroup Tth in C0 is given by those fC0 such that Lhf(x) → 0, as xa and xb. Moreover, (Tthf(x) – f(x))/tLhf(x) uniformly on any closed interval [a′, b′] ⊂ (a, b) and for any fC0C1([a, b]). The domain of the semigroup Tth in C00 is given by those fC00 such that Lhf(x) → 0 and (Lhf)′(x) → 0, as xa and xb.

2. The next key step is to show that the semigroups Tth are bounded in C00 uniformly in h. To this end let us differentiate the equation t = Lhft satisfied by ft = Tthf with f from the C00. Thus for gt = ∂ft/∂x we get the equation

g˙t(x)=Ωhgt(x)=Lhgt(x)+Lhft(x)ft(x)χh(bx)ν(x,bx)+ft(x)χh(ax)ν(x,ax),(5.7)

with

Lhf(x)=axbx(f(x+y)f(x))(1χh(y))ν(x,y)dyft(x)axν(x,y)dy+bxν(x,y)dy.

Writing

f(x)χh(bx)ν(x,bx)+f(x)χh(ax)ν(x,ax)
=xbg(y)dyχh(bx)ν(x,bx)+axg(y)dyχh(ax)ν(x,ax)
=xb(g(y)g(x))dyχh(bx)ν(x,bx)+ax(g(y)g(x))dyχh(ax)ν(x,ax)
+g(x)[χh(bx)(bx)ν(x,bx)+χh(xa)(xa)ν(x,ax)],

we get

Ωhg(x)=L~hg+Lhf,

with

L~hg(x)=axbx(g(x+y)g(x))(1χh(y))ν(x,y)dyk~(x)g(x),(5.8)

where

k~(x)=ka(x)+kb(x)χh(bx)(bx)ν(x,bx)χh(xa)(xa)ν(x,ax).

The key point is now the observation that, by (5.5),

k~(x)(1ϰ)k(x).

Hence we can show by the same perturbation argument as used above for Lh that h generates a strongly continuous semigroup in C0. Moreover, since all terms in the expression for h are conditionally positive operators this semigroup is a group of contractions and thus is bounded uniformly in h. Expressing f via g as above in Lhf, we observe that this operator becomes a uniformly bounded operator in C([a, b]) and hence, using again perturbation argument, we conclude that the semigroups Tht in C00 are uniformly bounded in h.

3. Let us show now that the operators Tth converge strongly in C0, as h → 0. To compare these operators for different h we shall use the following standard (and easy to prove) formula (see e.g. [40])

(Tth1Tth2)f=0tTtsh2(Lh1Lh2)Tsh1ds(5.9)

expressing the difference of the semigroups in terms of the difference of their generators. For arbitrary h1 > h2 and fC00, we have

(Lh1Lh2)ϕ(x)=axbx(f(x+y)f(x))(χh2χh1(y))ν(x,y)dy

and thus,

|(Lh1Lh2)ϕ(x)||y|h1f|y|ν(x,y)dyf|y|h1|y|ν~(y)dy.

Since ∥(Tthf)′∥ is uniformly bounded by Step 2, we have

(Tth1Tth2)f=o(1)tfC1,h10.(5.10)

Therefore the family Tthf converges in C0 to a family Ttf, as h → 0, for any fC00. By the standard density argument this convergence holds also for any fC0 and the limiting family Tt specifies a strongly continuous semigroup in C(Rd).

4. Writing

Ttft=TtTthft+Tthft

and noting that by (5.10) the first term is of order o(1) ∥f′∥, as h → 0, allows one to conclude that, for any fC0C1([a, b]), (Ttf(x) – f(x))/tLkillf(x) uniformly on any closed interval [a′, b′] ⊂ (a, b). It follows that Dkill belongs to the domain of Lkill.□

Having proved Theorem 5.1, we can safely apply the resolvent of the operator Lkill to obtain the solutions for the equations

D[a,b](ν)f=λf+g

with vanishing boundary conditions, which are classical (lie in the domain of the generator of Lkill) for gC0 and generalized otherwise.

By shifting the solutions to (5.2) can be reduced to the boundary problem with vanishing boundary conditions. For instance, writing f = u + f0 with f0 = fa + (xa)(fbfa)/(ba), the spectral problem

D[a,b](ν)f=Lf=λf,f|t=a=fa,f|t=b=fb,(5.11)

is reduced to the problem

D[a,b](ν)u=Lf=λuλf0D[a,b](ν)f0,u(a)=u(b)=0,(5.12)

which is solved by the equation

f(x)=f0(x)+u(x)=f0+Rkillλ(λf0+D[a,b](ν)f0).(5.13)

Using the technique of Dynkin’s martingale (see Section 4.4), this can be expressed in probabilistic (or path integral) terms as

f(x)=f0(x)+E0τxeλs(λf0+D[a,b](ν)f0)(Xsx)ds,(5.14)

where Xsx is the process generated by Lkill and τx is its killing time.

Alternatively, to solve the equations with CD derivatives one can work directly with the process generated by (5.15) by stopping it when reaching the boundary. To justify this approach one needs the following result.

### Theorem 5.2

Assume that (A)-(C). Then the operatorLof(5.15)generates a Feller semigroup inC([a, b]) and also a strongly continuous semigroup in the spaceC0 = {fC1([a, b]) : f′(a) = f′(b) = 0}.

### Proof

It is a consequence of Theorem 5.1. Let us consider the semigroup generated by L in the space C0. Differentiating the equation = Lf we get (after some handy cancellations) for the derivative g = f′ the equation

g˙=(Lf)=axbx(f(x+y)f(x))ν(x,y)dyg(x))×bxν(x,y)dy+axν(x,y)dy+axbx(f(x+y)f(x))ν(x,y)dy+(f(b)f(x))bxν(x,y)dy+(f(a)f(x))axν(x,y)dy.(5.15)

By Theorem 5.1, the first two terms on the r.h.s generates a well-defined Feller semigroup in C0. All other terms (the terms with ν′) are bounded when expressed in terms of g, and hence application of the perturbation theory yields the claimed result.□

### 5.2 Mixed RL-CD-Hilfer boundary-value problems

As was noted (see (3.50)), the usual boundary-value problem for RL derivative imposes a bit artificial integral boundary condition on the unknown function. However, under certain assumptions on the source function g the usual Cauchy problem,

Da+(ν)f(x)=λf(x)+g(x),f(a)=Y,(5.16)

can be solved. To see how it works let us analyse a bit more general case of mixed RL and CD derivatives:

(γD0+(ν)+δD0+(ν))f(x)=λf(x)+g(x),f(0)=Y,(5.17)

where γ, δ are nonnegative constants, g is a locally integrable function vanishing to the left of 0, and the derivative operators are of form (2.24) and (2.25). For γ = 0, problem (5.17) turns to the RL problem (5.16).

Acting as earlier, we make a substitution of the unknown function f(x) = u(x) + Y, which turns (5.17) to the problem

(γ+δ)D0+(ν)u(x)=(γ+δ)D0+(ν)u(x)=λu(x)λY+g(x)δYxν(dy),(5.18)

with u(0) = 0, because RL and CD derivatives coincide for functions vanishing at the boundary and D0+(ν)10(x)=xν(dy). Since g is locally integrable and thus an element of the space of generalized functions D′(R), we can apply Proposition 2.2 to get the unique generalized solution to this problem

u(x)=0[g(xy)δYxyν(dz)λY]Uλ(ν)(dy).

In order to satisfy the boundary requirement classically, that is, to have limx→0u(x) = 0, it is sufficient that

supyg(y)δYyν(dz)<(5.19)

at least for y from some neighborhood of the origin. If g is continuous for x > 0 and satisfies this condition we obtain the classical solution to problem (5.18) and thus to problem (5.17).

Similarly one can solve the boundary-value problem with mixed Hilfer and CD derivatives of the type

(γD0+(ν)+δHiDa+α,β)f(x)=λf(x)+g(x),f(a)=Y,(5.20)

The sufficient condition (5.19) for the classical solvability turns to the condition

supyg(y)δYyαΓ(1α)<,(5.21)

at least for y from some neighborhood of the origin, because HiDa+α,β10(x) = xα/Γ(1 – α).

### 5.3 Higher order and partial derivatives and related equations

Let us mention briefly how the equations with higher order and partial derivatives can be dealt with.

1. Let us consider, for instance, the equation

(Da+(ν))kf(x)=λf(x),f(a)=Y0,Da+(ν))lf(a)=Yl,l=1,k1,(5.22)

where Da+(ν) is given by (2.24). Introducing the vector-valued unknown function F = (f0, ⋯, fk–1) with f0 = f, fl = (Da+(ν))lf, l = 1, ⋯, k – 1, problem (5.22) rewrites in the vector form as

(Da+(ν))kF(x)=AF(x)=010000100001λ0F(x),F(a)=Y=Y0Y1Yk1.(5.23)

By (4.41) the solution to this problem can be written as F(x)=Eexp{Aτx(ν)}Y.

2. Let us consider the partial fractional differential equation

Da2+(ν2),x2Da1+(ν1),x1f(x1,x2)=λf(x1,x2)+g(x1,x2),x1a1,x2a2,(5.24)

with the simplest vanishing boundary conditions f(a1, x2) = f(x1, a2) = 0, where Daj+(νj),xj acts on the variable xj, j = 1, 2. Introducing the vector-valued function F = (f1, f2), f = f1, f2 = Da1+(ν1),x1f, we rewrite equation (5.24) in the matrix form:

DF(x)=Da1+(ν1),x1f1Da1+(ν1),x1f2=AF(x)=01λ0F=f2λf1.(5.25)

The operator D generates a Markov process on the triples (j, x1, x2), or, in other words on the two copies of the orthant {x1a1, x2x2} such that x1 is decreasing according to the generator Da1+(ν1),x1 on one of the copies and x2 is decreasing independently according to the generator Da2+(ν2),x2 on the other copy. The solution can be again expressed either via the resolvent of D or via Dynkin’s martingale.

The literature on fractional calculus is enormous. We shall mention only the sources most closely related to the probabilistic point of view of the present paper. Some historical reviews can be found in e.g. [23], [31], the wealth of physical and economics application, e.g. in [35], [73], [74], [75], numerical methods are dealt with in monographs [9], [50]. One of the main impetus for the physics interest in fractional equations was in fact probabilistic in nature. It was inspired by their appearance as the scaled limits of continuous time random walks, see e.g. [58], [38], [49] and extensive bibliography therein.

Fractional Schrödinger equation is getting popularity in physics community, see e.g. [48], [61], [10], [17] and references therein. For the fractional versions of the wave equations we refer to [66]. Of importance for physics are also fractional kinetic equations with application to statistical mechanics and fractional stochastic PDEs. For these developments we refer to [37], [45] and [51], [55], [79], [80] and references therein.

An insightful collection of references on nonlinear fractional equations and their applications can be found in [64]. Mild forms are studied in detail in [46], [47] and [44], where they are applied to the theory of fractional Hamilton-Jacobi-Bellman equations.

For various approaches to equations in bounded domains we refer to [2], [13], [19], [59] Two-sided and multidimensional problems appear naturally in the fractional calculus of variations, see [56], [5]. The optimization problems of this theory are formulated in terms of the certain class of fractional equations on bounded domains, the so-called fractional Euler-Lagrange equations. Their analysis was initiated seemingly in [12]. Problems with partial fractional derivatives are studied in detail in [1].

A handy tool for dealing with fractional equations is based on the method of duality, which we did not discussed here, see [42].

The method of generalized Mittag-Leffler functions developed above cannot be directly extended to deal with non-autonomous equations of the type Da+(ν)f(x) = A(x)f(x) + g(x) with a family of operators A(x) depending on x. The relevant modification of the theory is developed in [43], [44] and is based on the method of chronological operator-valued Feynman-Kac formulae.

The present paper is mostly based on the ideas suggested by the author in [41] and further developed in [20], [21], [22], [46] and [43].

## Acknowledgements

The author gratefully acknowledge support by the Russian Academic Excellence Project ’5-100’.

## References

[1] S. Abbas, M. Benchohra and G.M. N’Guérékata, Topics in Fractional Differential Equations. Developments in Mathematics, Vol. 27. Springer, New York, 2012.10.1007/978-1-4614-4036-9Search in Google Scholar

[2] O.P. Agrawal, Solution for a fractional diffusion-wave equation defined in a bounded domain. Fractional order calculus and its applications. Nonlinear Dynam. 29 (2002), 145–155.10.1023/A:1016539022492Search in Google Scholar

[3] O. Agrawal, Some generalized fractional calculus operators and their applications in integral equations. Fract. Calc. Appl. Anal. 15, No 4 (2012), 700–711; 10.2478/s13540-012-0047-7; https://www.degruyter.com/view/j/fca.2012.15.issue-4/issue-files/fca.2012.15.issue-4.xml.Search in Google Scholar

[4] O.P. Agrawal, Generalized variational problems and Euler–Lagrange equations. Computers and Math. with Appl. 59 (2010), 1852–1864.10.1016/j.camwa.2009.08.029Search in Google Scholar

[5] R. Almeida, Sh. Pooseh and D.F.M. Torres, Computational Methods in the Fractional Calculus of Variations. Imperial College Press, London, 2015.10.1142/p991Search in Google Scholar

[6] M.S. Alrawashdeh, J.F. Kelly, M.M. Meerschaert and H.-P. Scheffler, Applications of inverse tempered stable subordinators. Computers and Math. with Appl. 73 (2017), 892–905.10.1016/j.camwa.2016.07.026Search in Google Scholar

[7] B. Al-Saqabi and V. Kiryakova, Explicit solutions of fractional integral and differential equations involving Erdélyi-Kober operators. Appl. Math. and Computation95 (1998), 1–13.10.1016/S0096-3003(97)10095-9Search in Google Scholar

[8] T. Atanackovic, D. Dolicanin, S. Pilipovic and B. Stankovic, Cauchy problems for some classes of linear fractional differential equations. Fract. Calc. Appl. Anal. 17, No 4 (2014), 1039–1059; 10.2478/s13540-014-0213-1; https://www.degruyter.com/view/j/fca.2014.17.issue-4/issue-files/fca.2014.17.issue-4.xml.Search in Google Scholar

[9] D. Baleanu, K. Diethelm, E. Scalas and J.J. Trujillo, Fractional Calculus. Models and Numerical Methods, Second Ed. Ser. on Complexity, Nonlinearity and Chaos 5, World Scientific, Hackensack, NJ, 2017.10.1142/10044Search in Google Scholar

[10] S.S. Bayin, Time fractional Schrödinger equation. Fox’s H-functions and the effective potential. J. Math. Phys. 54 (2013), # 012103.10.1063/1.4773100Search in Google Scholar

[11] N.H. Bingham, Limit theorems for occupation times of Markov processes. Z. Wahrsch. verw. Geb. 17 (1971), 1–22.10.1007/BF00538470Search in Google Scholar

[12] L. Bourdin, Existence of a weak solution for fractional Euler-Lagrange equations. J. Math. Anal. Appl. 399, No 1 (2013), 239–251.10.1016/j.jmaa.2012.10.008Search in Google Scholar

[13] Z.-Q. Chen, M.M. Meerschaert and E. Nane, Space time fractional diffusion on bounded domains. J. Math. Anal. Appl. 393 (2012), 479–488.10.1016/j.jmaa.2012.04.032Search in Google Scholar

[14] M.M. Dzherbashian, Integral Transforms and Representations of Functions in the Complex Plane. Nauka, Moscow, 1966 (in Russian).Search in Google Scholar

[15] R. Garra, A. Giusti, F. Mainardi anf G. Pagnini, Fractional relaxation with time-varying coefficient. Fract. Calc. Appl. Anal. 17, No 2 (2014), 424–439; 10.2478/s13540-014-0178-0; https://www.degruyter.com/view/j/fca.2014.17.issue-2/issue-files/fca.2014.17.issue-2.xml.Search in Google Scholar

[16] I.I. Gikhman and A.V. Skorokhod, The Theory of Stochastic Processes, Vol. 2. Moscow, Nauka, 1973 (In Russian).Search in Google Scholar

[17] P. Górka, H. Prado and J. Trujillo, The time fractional Schrödinger equation on Hilbert space. Integral Equations Operator Theory87, No 1 (2017), 1-–14.10.1007/s00020-017-2341-6Search in Google Scholar

[18] R. Gorenflo, Y. Luchko and M. Stojanovic, Fundamental solution of a distributed order time-fractional diffusion-wave equation as probability density. Fract. Calc. Appl. Anal. 16, No 2 (2013), 297–316; 10.2478/s13540-013-0019-6; https://www.degruyter.com/view/j/fca.2013.16.issue-2/issue-files/fca.2013.16.issue-2.xml.Search in Google Scholar

[19] J. Henderson and R. Luca, Boundary Value Problems for Systems of Differential, Difference and Fractional Equations. Positive Solutions. Elsevier, Amsterdam, 2016.10.1186/s13661-016-0569-8Search in Google Scholar

[20] M.E. Hernandez-Hernandez and V.N. Kolokoltsov, On the probabilistic approach to the solution of generalized fractional differential equations of CD and Riemann-Liouville type. J. of Fract. Calc. and Appl. 7, No 1 (2016), 147–175.Search in Google Scholar

[21] M.E. Hernandez-Hernandez and V.N. Kolokoltsov, On the solution of two-sided fractional ordinary differential equations of Caputo type. Fract. Calc. Appl. Anal. 19, No 6 (2016), 1393–1413; 10.1515/fca-2016-0072; https://www.degruyter.com/view/j/fca.2016.19.issue-6/issue-files/fca.2016.19.issue-6.xml.Search in Google Scholar

[22] M.E. Hernandez-Hernandez, V.N. Kolokoltsov and L. Toniazzi, Generalised fractional evolution equations of Caputo type. Chaos Solitons Fractals102 (2017), 184—196.10.1016/j.chaos.2017.05.005Search in Google Scholar

[23] R. Hilfer, Threefold introduction to fractional derivatives. In: Anomalous Transport (Edited by R. Klages, G. Radons, and I.M. Sokolov), Wiley, 2008, 17–74.10.1002/9783527622979.ch2Search in Google Scholar

[24] S.L. Kalla, Operators of fractional integration. In: Lecture Notes in Math. 798 (Proc. Conf. Analytic Functions, Kozubnik 1979), Springer, New York, 1980, 258–280.10.1007/BFb0097270Search in Google Scholar

[25] U.N. Katugampola, New approach to a generalized fractional integral. Appl. Math. and Computation218, No 3 (2011), 860–865.10.1016/j.amc.2011.03.062Search in Google Scholar

[26] M. Kelbert, V. Konakov and S. Menozzi, Weak error for continuous time Markov chains related to fractional in time P(I)DEs. Stochastic Process. Appl. 126, No 4 (2016), 1145-–1183.10.1016/j.spa.2015.10.013Search in Google Scholar

[27] A. Kilbas, Hadamard-type fractional calculus. J. Korean Math. Soc. 38, No 6 (2001), 1191–1204.Search in Google Scholar

[28] A. Kilbas and J.J. Trujillo, Hadamard-type integrals as G-transforms. Integr. Transf. Spec. Funct. 14, No 5 (2003), 413–427.10.1080/1065246031000074443Search in Google Scholar

[29] A.A. Kilbas, H.M. Srivastava and J.J. Trujillo, Theory and Applications of Fractional Differential Equations. Elsevier, Amsterdam, 2006.Search in Google Scholar

[30] V. Kiryakova, Generalized Fractional Calculus and Applications. Pitman Res. Notes in Math. Ser. 301. Longman Scientific, Harlow; Co-publ. with John Wiley and Sons, New York, 1994.Search in Google Scholar

[31] V. Kiryakova, A brief story about the operators of the generalized fractional calculus. Fract. Calc. Appl. Anal. 11, No 2 (2008), 203–220; at http://www.math.bas.bg/complan/fcaa.Search in Google Scholar

[32] V. Kiryakova, From the hyper-Bessel operators of Dimovski to the generalized fractional calculus. Fract. Calc. Appl. Anal. 17, No 4 (2014), 977–1000; 10.2478/s13540-014-0210-4; https://www.degruyter.com/view/j/fca.2014.17.issue-4/issue-files/fca.2014.17.issue-4.xml.Search in Google Scholar

[33] V. Kiryakova, Fractional calculus operators of special functions? The result is well predictable! Chaos, Solitons and Fractals102 (2017), 2–15; 10.1016/j.chaos.2017.03.006.Search in Google Scholar

[34] V. Kiryakova and Yu. Luchko, Multiple Erdélyi–Kober integrals and derivatives as operators of generalized fractional calculus. Ch. 6 in: A. Kochubei, Yu. Luchko (Eds.), Handbook of Fractional Calculus with Applications. Vol. 1: Basic Theory, De Gruyter, Berlin-Boston, 2019, 127–158; 10.1515/9783110571622-006.Search in Google Scholar

[35] R. Klages, G. Radons and I.M. Sokolov, Anomalous Transport: Foundations and Applications. Wiley-VCH, Weinheim, 2008.10.1002/9783527622979Search in Google Scholar

[36] A. Kochubei, General fractional calculus, evolution equations, and renewal processes. Integral Equations Operator Theory71, No 4 (2011), 583–600.10.1007/s00020-011-1918-8Search in Google Scholar

[37] A.N. Kochubei and Y. Kondratiev, Fractional kinetic hierarchies and intermittency. Kinet. Relat. Models10, No 3 (2017), 725–740.10.3934/krm.2017029Search in Google Scholar

[38] V. Kolokoltsov, Generalized Continuous-Time Random Walks (CTRW), Subordination by hitting times and fractional dynamics. Theory of Probability and its Applications53, No 4 (2009), 594–609; arXiv:0706.1928v1[math.PR] 2007.10.1137/S0040585X97983857Search in Google Scholar

[39] V.N. Kolokoltsov, Nonlinear Markov Processes and Kinetic Equations. Cambridge Tracks in Math. 182, Cambridge Univ. Press, 2010.10.1017/CBO9780511760303Search in Google Scholar

[40] V.N. Kolokoltsov, Markov Processes, Semigroups and Generators. DeGruyter Studies in Math. 38, De Gruyter, 2011.10.1515/9783110250114Search in Google Scholar

[41] V.N. Kolokoltsov, On fully mixed and multidimensional extensions of the Caputo and Riemann-Liouville derivatives, related Markov processes and fractional differential equations. Fract. Calc. Appl. Anal. 18, No 4 (2015), 1039–1073; 10.1515/fca-2015-0060; https://www.degruyter.com/view/j/fca.2015.18.issue-4/issue-files/fca.2015.18.issue-4.xml; and http://arxiv.org/abs/1501.03925.Search in Google Scholar

[42] V.N. Kolokoltsov, Stochastic monotonicity and duality of kth order with application to put-call symmetry of powered options. J. of Applied Probability52, No 1 (2015), 82–101; and http://arxiv.org/abs/1405.3894.10.1239/jap/1429282608Search in Google Scholar

[43] V. Kolokoltsov, Chronological operator-valued Feynman-Kac formulae for generalized fractional evolutions. arXiv:1705.08157 (2017); at https://arxiv.org/pdf/1705.08157.pdf.Search in Google Scholar

[44] V. Kolokoltsov, Differential Equation on Measures and Functional Spaces. Birkhäuser, 2019.10.1007/978-3-030-03377-4Search in Google Scholar

[45] V. Kolokoltsov and O.A. Malafeyev, Many Agent Games in Socio-Economic Systems: Corruption, Inspection, Coalition Building, Network Growth, Security. Springer, 2019.10.1007/978-3-030-12371-0Search in Google Scholar

[46] V. Kolokoltsov and M. Veretennikova, Fractional Hamilton Jacobi Bellman equations for scaled limits of controlled Continuous Time Random Walks. Commun. in Appl. and Industr. Math. 6, No 1 (2014), e-484; http://caim.simai.eu/index.php/caim.10.1685/journal.caim.484Search in Google Scholar

[47] V. Kolokoltsov and M. Veretennikova, Well-posedness and regularity of the Cauchy problem for nonlinear fractional in time and space equations. Fract. Differential Calc. 4, No 1 (2014), 1–30; http://files.ele-math.com/articles/fdc-04-01.pdf.10.7153/fdc-04-01Search in Google Scholar

[48] N. Laskin, Fractional Schrödinger equation. Phys. Rev. E66 (2002), # 056108.10.1142/9789813223806_0003Search in Google Scholar

[49] N.N. Leonenko, M.M. Meerschaert and A. Sikorskii, Correlation structure of fractional Pearson diffusions. Comput. Math. Appl. 66, No 5 (2013), 737–745.10.1016/j.camwa.2013.01.009Search in Google Scholar PubMed PubMed Central

[50] Ch. Li and F. Zeng, Numerical Methods for Fractional Calculus. CRC Press, Boca Raton, 2015.10.1201/b18503Search in Google Scholar

[51] W. Liu, M. Röckner and J.L. da Silva, Quasi-linear (stochastic) partial differential equations with time-fractional derivatives. SIAM J. Math. Anal. 50, No 3 (2018), 2588–2607.10.1137/17M1144593Search in Google Scholar

[52] Yu. Luchko, Operational rules for a mixed operator of the Erdélyi-Kober type. Fract. Calc. Appl. Anal. 7, No 3 (2004), 339–364.Search in Google Scholar

[53] Yu. Luchko, Integral transforms of the Mellin convolution type and their generating operators. Integr. Transf. Spec. Funct. 19 (2008), 809–851; 10.1080/10652460802091617.Search in Google Scholar

[54] Yu. Luchko, J.J. Trujillo, Caputo-type modification of the Erdélyi-Kober fractional derivative. Fract. Calc. Appl. Anal. 10, No 3 (2007), 249–267; at http://www.math.bas.bg/complan/fcaa.Search in Google Scholar

[55] G. Lv, J. Duan, H. Gao and J.-L. Wu, On a stochastic nonlocal conservation law in a bounded domain. Bull. Sci. Math. 140, No 6 (2016), 718–746.10.1016/j.bulsci.2016.03.003Search in Google Scholar

[56] A.B. Malinowska, T. Odzijewicz and D.F.M. Torres, Advanced Methods in the Fractional Calculus of Variations. Springer, Heidelberg, 2015.10.1007/978-3-319-14756-7Search in Google Scholar

[57] A.M. Mathai, R.K. Saxena and H.J. Haubold, TheH-Function, Theory and Applications. Springer, 2010.10.1007/978-1-4419-0916-9Search in Google Scholar

[58] M.M. Meerschaert and A. Sikorskii, Stochastic Models for Fractional Calculus. De Gruyter Studies in Mathematics 43, 2012.10.1515/9783110258165Search in Google Scholar

[59] M. Meerschaert, E. Nane and P. Vellaisamy, Fractional Cauchy problems on bounded domains. The Annals of Probability37, No 3 (2009), 979–1007.10.1214/08-AOP426Search in Google Scholar

[60] M.M. Meerschaert and P. Straka, Inverse stable subordinators. Math. Model. Nat. Phenom. 8, No 2 (2013), 1–16.10.1051/mmnp/20138201Search in Google Scholar PubMed PubMed Central

[61] M. Naber, Time fractional Schrödinger equation. J. Math. Phys. 45 (2004), 3339–3352.10.1063/1.1769611Search in Google Scholar

[62] E. Orsinger and B. Toaldo, Space-time fractional equations and the related stable processes at random times. J. Theor. Probab. 30 (2017), 1–26.10.1007/s10959-015-0641-9Search in Google Scholar

[63] J. Paneva-Konovska, From Bessel to Multi-Index Mittag-Leffler Functions. World Scientific, 2017.10.1142/q0026Search in Google Scholar

[64] L. Plociniczak and M. Switala, Existence and uniqueness results for a time-fractional nonlinear diffusion equation. J. Math. Anal. Appl. 462, No 2 (2018), 1425–1434.10.1016/j.jmaa.2018.02.050Search in Google Scholar

[65] I. Podlubny, Fractional Differential Equations. Math. in Science and Engin. Ser. 198, Academic Press, San Diego, 1999.Search in Google Scholar

[66] A.V. Pskhu, Fundamental solution of the diffusive wave equation of fractional order (in Russian). Izvestia RAN, Ser. Math. 73, No 2 (2009), 141–182.Search in Google Scholar

[67] A.V. Pskhu, An initial value problem for a fractional-order linear ordinary differential equation. Mat. Sb. 202, No 4 (2011), 111–122; Engl. transl. in: Sb. Math. 202 (2011), 571–582.Search in Google Scholar

[68] M. Rivero, L. Rodríguez-Germá, J.J. Trujillo and M.P. Velasco, Fractional operators and some special functions. Computers and Math. with Appl. 59 (2010), 1822–1834.10.1016/j.camwa.2009.08.026Search in Google Scholar

[69] M. Rivero, J.J. Trujillo, M.P. Velasco and M. Pilar, On deterministic fractional models. In: New Trends in Nanotechnology and Fractional Calculus Applications, Springer, New York, 2010, 123—150.10.1007/978-90-481-3293-5_10Search in Google Scholar

[70] S. Samko, A. Kilbas, O. Marichev, Fractional Integrals and Derivatives: Theory and Applications, Gordon and Breach, 1993.Search in Google Scholar

[71] E. Scalas, A class of CTRWs: Compound fractional Poisson processes. In: Fractional Dynamics, World Sci. Publ., Hackensack, NJ, 2012, 353–374.10.1142/9789814340595_0015Search in Google Scholar

[72] R.L. Schilling, R. Song and Z. Vondracek, Bernstein Functions. Theory and Applications. Studies in Math. 37, De Gruyter, 2010.Search in Google Scholar

[73] V.E. Tarasov, Fractional Dynamics, Applications of Fractional Calculus to Dynamics of Particles, Fields and Media. Springer, Higher Education Press, 2011.Search in Google Scholar

[74] V.V. Uchaikin, Fractional Derivatives for Physicists and Engineers. Springer, 2012.10.1007/978-3-642-33911-0Search in Google Scholar

[75] B.J. West, Fractional Calculus View of Complexity. Tomorrow’s Science. CRC Press, Boca Raton, 2016.10.1201/b18911Search in Google Scholar

[76] S.B. Yakubovich and Yu.F. Luchko, The Hypergeometric Approach to Integral Transforms and Convolutions. Mathematics and its Applications 287, Kluwer Acad. Publ. Group, Dordrecht, 1994.10.1007/978-94-011-1196-6_21Search in Google Scholar

[77] V.M. Zolotarev, On analytic properties of stable distribution laws (in Russian). Vestnik Leningrad. Univ. 11, No 1 (1956), 49–52.Search in Google Scholar

[78] V.M. Zolotarev, One-Dimensional Stable Distributions. Moscow, Nauka, 1983 (in Russian); Engl. transl. in: Vol. 65 of Transl. of Math. Monographs AMS, Providence, Rhode Island, 1986.10.1090/mmono/065Search in Google Scholar

[79] G. Zou, G. Lv and J.-L. Wu, On the regularity of weak solutions to space-time fractional stochastic heat equations. Statist. Probab. Lett. 139 (2018), 84–89.10.1016/j.spl.2018.04.006Search in Google Scholar

[80] G. Zou, G. Lv and J.-L. Wu, Stochastic Navier-Stokes equations with Caputo derivative driven by fractional noises. J. Math. Anal. Appl. 461, No 1 (2018), 595–609.10.1016/j.jmaa.2018.01.027Search in Google Scholar

[81] S. Umarov and R. Gorenflo, Cauchy and nonlocal multi-point problems for distributed order pseudo-differential equations. Fract. Calc. Appl. Anal. 8, No 1 (2005), 73–86; http://www.math.bas.bg/complan/fcaa.10.4171/ZAA/1250Search in Google Scholar