Semi-Markov models and motion in heterogeneous media

TL;DR

This paper develops mathematical models for continuous time random walks in heterogeneous media, leading to variable order fractional equations that describe anomalous diffusion with state-dependent properties.

Contribution

It introduces integro-differential equations with position-dependent kernels for CTRWs and derives variable order fractional equations for complex media.

Findings

Derived Volterra-type integro-differential equations for CTRWs

Established connection to variable order fractional heat equations

Modeled anomalous diffusion in heterogeneous environments

Abstract

In this paper we study continuous time random walks (CTRWs) such that the holding time in each state has a distribution depending on the state itself. For such processes, we provide integro-differential (backward and forward) equations of Volterra type, exhibiting a position dependent convolution kernel. Particular attention is devoted to the case where the holding times have a power-law decaying density, whose exponent depends on the state itself, which leads to variable order fractional equations. A suitable limit yields a variable order fractional heat equation, which models anomalous diffusions in heterogeneous media.

Equations

$$\frac{\partial }{\partial t}p(x,y,t)=k{\mathcal{D}}_t^{1-\alpha }\frac{{\partial }^2}{\partial x^2}p(x,y,t).$$

$$\frac{d}{dt}{\int }_0^{t}p(x,y,s)k(t-s)ds-k(t)p(x,y,0)=\sum _z{g}_{x,z}p(z,y,t),$$

$${\mathcal{D}}_t^{\alpha (x)}f(x,t)=Gf(x,t)$$

$${\mathcal{D}}_t^{\alpha }u(t):=\frac{1}{\Gamma (1-\alpha )}\frac{d}{dt}{\int }_0^{t}u(s)(t-s{)}^{-\alpha }ds-\frac{t^{-\alpha }u(0)}{\Gamma (1-\alpha )},$$

$${\mathcal{D}}_t^{\alpha (x)}p(x,y,t)=\sum _z{g}_{x,z}p(z,y,t).$$

$$\frac{d}{dt}p(x,y,t)=\sum _z{g}_{z,y}{}^R{\mathcal{D}}_t^{1-\alpha (z)}p(x,z,t)$$

$${}^R{\mathcal{D}}_t^{\beta }u(t):=\frac{1}{\Gamma (1-\beta )}\frac{d}{dt}{\int }_0^{t}u(s)(t-s{)}^{-\beta }ds$$

$$\frac{\partial }{\partial t}p(x,y,t)=\frac{1}{2}\frac{{\partial }^2}{\partial x^2}\left(k(x){\mathcal{D}}_t^{1-\alpha (x)}p(x,y,t)\right).$$

$$\frac{d}{dt}p(x,y,t)=\sum _z{g}_{z,y}\frac{d}{dt}{\int }_0^{t}p(x,z,s)k(t-s,z)ds.$$

$$M(t)=X_n{V}_n\le t<V_{n+1}\text{ where }{V}_0=0V_n=\sum _{k=0}^{n-1}{E}_k$$

$$h_{i,j}=P(X_{n+1}=j|X_n=i),\forall n\in \mathbb{N}i,j\in \mathcal{S},$$

$$P(E_n>t|X_n=i)=e^{-{\lambda }_{i}t}\forall n\in \mathbb{N},t\ge 0.$$

$$p_{i,j}(t)=P(M(t)=j|M(0)=i)$$

$$g_{i,j}={\lambda }_i(h_{i,j}-{\delta }_{i,j})$$

$$p_{i,j}(dt)=\{\begin{array}{ll}{g}_{i,j}dt={\lambda }_i{h}_{i,j}dt,& i\ne j,\\ 1+g_{i,i}dt=1-{\lambda }_{i}dt+{\lambda }_i{h}_{i,i}dt,& i=j.\end{array}$$

$$\zeta :=\underset{n}{\sup }{V}_n=\sum _n{E}_n=\infty ,$$

$$\frac{d}{dt}{p}_{i,j}(t)=\sum _k{g}_{i,k}{p}_{k,j}(t),p_{i,j}(0)={\delta }_{i,j},$$

$$\frac{d}{dt}{p}_{i,j}(t)=\sum _k{p}_{i,k}(t)g_{k,j},p_{i,j}(0)={\delta }_{i,j},$$

$$\frac{d}{dt}P(t)=GP(t)=P(t)G,P(0)=I.$$

$$X(t)=X_n,T_n\le t<T_{n+1},\text{where }{T}_0=0T_n=\sum _{k=0}^{n-1}{J}_k,$$

$$h_{i,j}=P(X_{n+1}=j|X_n=i),\forall n\in \mathbb{N},i,j\in \mathcal{S},$$

$$P(J_n>t|X_n=i)={\overline{F}}_i(t),\forall n\in \mathbb{N},t\ge 0,$$

$$P(J_n>t|X_n=i)=E_{\alpha }(-{\lambda }_i{t}^{\alpha }),\forall n\in \mathbb{N},t\ge 0,$$

$$E_{\alpha }(x):=\sum _{k=0}^{\infty }\frac{x^k}{\Gamma (1+\alpha k)}$$

$$E_{\alpha }(-\lambda t^{\alpha })\sim C_{\lambda }\frac{t^{-\alpha }}{\Gamma (1-\alpha )},C_{\lambda }>0,$$

$$\gamma (s)=s-\sup \left\{w\le s:X(w)\ne X(s)\right\},$$

$$p_{i,j}(t)$$

$$=P(X(t+\tau )=j|X(\tau )=i,\gamma (\tau )=0)$$

$$p_{i,j}(t)=P\{J_i>t\}{\delta }_{i,j}+{\int }_0^t\sum _l{h}_{i,l}{p}_{l,j}(t-s)\text{mathpzc}{f}_i(ds)$$

$$P\left(X(t)=j\mid X(0)=i,\gamma (0)=0\right)=$$

$$+P\left(X(t)=j,J_0\le t\mid X(0)=i,\gamma (0)=0\right).$$

Full text

1. 1 Introduction

2. 2 Semi-Markov models for motion in homogeneous media

3. 3 Semi-Markov models for motion in heterogeneous media

4. 3.1 The time-change by a dependent time process

5. 3.2 Variable order backward equations

6. 4 The forward equations of semi-Markov processes in heterogeneous media

7. 5 Convergence to the variable order fractional diffusion

8. 6 Arbitrary holding times and integro-differential Volterra equations

9. 6.1 Integro-differential Volterra equations with position dependent kernel

1. Introduction

We here consider continuous time random walks (CTRWs) on countable state spaces. It is assumed that every time the walker jumps, the future trajectory becomes independent of its past, namely the next position and the next jump time depend only on the current position; furthermore, in a generic time instant, the future behavior is assumed to be also depending on the time already spent in the current position. Such a process is said to be semi-Markovian. If the waiting times between jumps follow an exponential distribution, then, due to the lack of memory property, the random walk is a Markov process.

It is well known that suitable (Markovian) random walks are good approximations of the Brownian motion. In the last decades it has been noticed that the CTRWs whose waiting times have distribution with a power-law decay, played a central role in statistical physics because they are good approximations of anomalous diffusion processes, where the mean square displacement grows as $\overline{x^{2}}\sim t^{\alpha},\alpha\in(0,1)$, and therefore slower than a standard Brownian motion (for a complete overview on this matter consult [34] and references therein). In these models each site exercises a trapping effect which, in some sense, delayes the time with respect to a corresponding Markov process.

It turns out that these facts can be framed in a nice probabilistic setting: to construct a large class of CTRWs, it is sufficient to replace the deterministic time $t$ of a Markov process by an independent inverse stable subordinator (on this point see, for example, the instructive discussion in [29]). It is well known (see for example [18, page 365] and [33]) that the transition probabilities of the correponding CTRW follows both the fractional backward and forward equations. Such equations are obtained from Kolmogorov backward and forward equations by replacing the time derivative with the fractional one, which introduces a memory effect by a convolution integral with a slowly decaying power-law kernel.

A suitable scaling then leads to anomalous diffusion processes, whose p.d.f. solves the Fokker-Planck equation (see, e.g., [34])

[TABLE]

It has been empirically confirmed (see [34]) that these models are particularly effective in a number of applications, e.g., for modeling diffusion in percolative and porous systems, charge carrier transport in amorphous semiconductors, nuclear magnetic resonance, motion on fractal geometries, dynamics of a bead in a polymeric network, protein conformational dynamics and many others.

We finally stress that many aspects of the theory hold as they are if the distribution of the holding times is arbitrary and not necessarily with a power law decay provided that it satisfies some mild assumptions (see, for example, the discussion in [29, Section 4]). In this case the random time process is given by the inverse of a generic subordinator and the corresponding backward equations have the form of a Volterra integro-differential equation

[TABLE]

where $g_{x,y}:=(G)_{x,y}$ and $G$ is the Markovian generator (see [9; 33; 44; 45] for the general theory and [11; 12] for some particular cases). The fractional case is more familiar in statistical physics because it is widely used in applications.

Up to now, we have only considered the simplest forms of fractional kinetic equations, where the fractional index $\alpha$ is constant. On the other hand, it is clear that further theoretical investigations are required for the description of more complicated (and more realistic) random processes, where the particle moves in an inhomogeneous environment; as we will discuss in the paper it turns out that this leads to equations of multi-fractional type. Equations with time-fractional derivative whose order depends on space have been studied in [38], where the authors considered a CTRW, say $X(t)$, $t\geq 0$, such that the function $f(x,t)=\mathds{E}\{u(X(t))\mid X(0)=x\}$, for a suitable test function $u$, solves the fractional backward equation

[TABLE]

where $\mathcal{D}_{t}^{\alpha(x)}$ denotes the $\alpha$-fractional derivative in the sense of Caputo-Dzerbayshan, i.e., for $\alpha\in(0,1)$,

[TABLE]

for any function $u$ such that the above integral is differentiable. In such a model, the trapping effect is not exercised with the same intensity at all sites. Indeed, when the particle reaches the state $x$, it is trapped for a time interval with density $\psi(t)\sim t^{-1-\alpha(x)}$ before jumping to another point. Thus the time delay is stronger when the particle is located at points with small values of $\alpha$. This leads to the fundamental fact regarding the time-change relation $X(t)=M(L(t))$: the time process $L$ and the Markov process $M$ are not independent. Such a construction is far from trivial, since $L$ is the right continuous inverse of a non-decreasing additive process also called time-inhomogeneous subordinator (for basic information consult [41] and [37]).

In the case of a countable state space $\mathcal{S}$, we here present the derivation of the backward equation

[TABLE]

We further introduce the forward equation

[TABLE]

where ${}^{R}\mathcal{D}_{t}^{1-\alpha(x)}$ denotes the fractional derivative in the sense of Riemann-Liouville, i.e., for $\beta\in(0,1),$

[TABLE]

for any function $u$ such that the above operator is well defined. Further we explain why eq. (1.6) is a true forward equation in the classical sense of Kolmogorov.

Therefore, this paper also creates a further bridge between the theory of semi-Markov processes and models of motion in heterogeneous media (a different theory concerning motions at finite velocity is discussed in [19]): on the one hand, there is the theory of semi-Markov processes, on the other hand, there are recent works concerning the fractional diffusion equation with multifractional index:

[TABLE]

Such equation has been derived in [8], but the related theory is still at an early stage, especially with regard to physical and phenomenological aspects. However, the use of a multi-fractional index $\alpha(x)$ is more realistic in the description of physical phenomena. Indeed, it takes into account the possibility of heterogeneous media, or, more simply, it considers homogeneous media where some impurities are scattered.

Finally, in the same spirit as (1.2) we show that previous models can be generalized by letting the (random) trapping effects having an arbitrary density, subject to some mild assumptions. These models then yield integro-differential equations of Volterra type, with a position dependent kernel of convolution, i.e.,

[TABLE]

The plan of the paper is the following. In sections 2, 3, 4 and 5, we consider the case of power-law holding times, which is the most familiar case in statistical physics. In particular, in section 2 we review (in our notations) some known facts on CTRWs in a homogeneous environment, where the fractional index $\alpha$ is assumed to be constant in space. Section 3 and 4 regard CTRWs in heterogeneous environment, where the fractional index is assumed to be space-dependent. Section 5 deals with the derivation of the multifractional diffusion equation. In section 6 many results are extended to the case where the holding times follow more general distributions.

2. Semi-Markov models for motion in homogeneous media

Before moving to heterogenous media we collect some results from the literature concerning classical models which will be used in the subsequent parts. As we stated in the introduction, the most popular model in statistical physics is related with holding times in each site having a density $\psi(t)\sim Ct^{-\alpha-1}$, $C>0$, $\alpha\in(0,1)$, with a power law decay. So, for example $\psi(t)=-(d/dt)E_{\alpha}(-\lambda t^{\alpha})$ (compare with [42, eq. (26)]), and this is related to fractional processes. Hence we focus the attention on this case to present the results concerning this theory.

First, in order to introduce the notation that we will use hereafter, we recall some basic facts regarding the classical theory of stepped Markov processes. Let us consider a continuous time Markov process $M$ with discrete state space $\mathcal{S}$

[TABLE]

where $X_{n}$ is a homogeneous discrete-time Markov chain on $\mathcal{S}$ with transition probabilities

[TABLE]

and the sojourn times are such that

[TABLE]

Let

[TABLE]

be the transition probabilities. The Markovian generator of $M$ is the matrix with elements

[TABLE]

where $\delta_{i,j}$ denotes the Kronecker symbol. Then the infinitesimal transition probabilities have the form

[TABLE]

It is enough for our models to consider the case in which, a.s.,

[TABLE]

so that the processes here are non explosive and hence we shall not consider what happens to a process after explosion. Under all these assumptions the functions $p_{i,j}(t)$, with $i,j\in S$ solve the Kolmogorov backward equations (e.g. [36, Sec. 2.8])

[TABLE]

as well as the Kolmogorov forward equations

[TABLE]

which can be written in compact matrix notation as

[TABLE]

We now consider a CTRW constructed in the same way of $M$, except for the distribution of the waiting times, which are no longer exponentially distributed. These processes are said to be semi-Markov processes in the sense of Gihman and Skorohod [20, Chapter 3]. Hence let $X(t)$ be

[TABLE]

where $X_{n}$ is a homogeneous discrete time Markov chain on $\mathcal{S}$ with transition probabilities

[TABLE]

and the sojourn times are such that

[TABLE]

where $F_{i}(t)=1-\overline{F}_{i}(t)$ is an arbirtrary c.d.f. We will devote particular attention to the case

[TABLE]

for $\lambda_{i}>0$, where

[TABLE]

is the Mittag-Leffler function. In this case (e.g. [42, eq. (26)])

[TABLE]

and the corresponding equations are fractional. The charactering property of semi-Markov processes is the following: by defining

[TABLE]

the sojourn time of $X$ in the current position, the couple $(X(t),\gamma(t))$ is a (strict) Markov process [20, Chapter 3, Section 3, Lemma 2]. This is to say that, when conditioning on the trajectory up to time $s$, future events depend not only on the current position $X(s)$, as it is for Markov processes, but also on the quantity $\gamma(s)$. Let

[TABLE]

be the transition probabilities (the second equality follows by time-homogeneity). We know from [19, page 20] that the transition probabilities solve the renewal equation

[TABLE]

where here $\mathpzc{f}_{i}(t)$ denotes a smooth density of $F_{i}(t)$. Note that (2.17), which provides a system of integral equations for the transition probabilities (2), comes from a very classical conditioning argument: fixing the time of the first jump $J_{0}$ and using the Markov property of the semi-Markov process at the jump times yields (see [19, page 19] for some details)

[TABLE]

A similar approach on semi-Markov processes, with an interesting discussion on exactly solvable models, can be found in [16].

The process (2.11) is known to have a deep connection to fractional calculus. Indeed, the following result holds.

Proposition 2.1**.**

The transition functions $p_{i,j}(t)$, $i,j\in\mathcal{S}$, defined in (2) solve the following system of backward equations

[TABLE]

as well as the system of “forward” equations

[TABLE]

Proof.

By the convolution Theorem we can compute the Laplace transform in (2.19) and we obtain

[TABLE]

Instead, by applying the Laplace tranform to (2.17) we have

[TABLE]

By setting $g_{i,j}=\lambda_{i}(h_{i,j}-\delta_{i,j})$, it is easy to show that (2.22) reduces to (2.21) and using the uniqueness theorem for Laplace tranforms eq. (2.19) is proved. Now if we apply again the Laplace transform to (2.20) we get that

[TABLE]

and thus the solution of (2.19) and (2.20) coincide. Indeed they can be obtained by solving either the system (2.21) or (2.23) in the variables $\tilde{p}_{i,j}(s)$, which, in compact operator form, reads

[TABLE]

where $I$ is the identity matrix, and this concludes the proof. ∎

Remark 2.2**.**

We observe that (2.19) derives directly by the renewal equation, which has a clear backward meaning. Instead the reason why we call eq. (2.20) “forward equation” is that it is formally obtained by introducing the fractional Caputo derivative in the Kolmogorov forward equation (2.9); the fact that it has a clear probabilistic interpretation has never been proved. A clear probabilistic meaning to (2.20) will be derived later in Section 4 from the general form of the forward equation of semi-Markov processes we will present. Concerning forward equations of semi-Markov processes, see also the discussion in [13].

It is well known that (2.11) can be equivalently constructed by replacing the time $t$ in (2.1) with the right continuous inverse of an independent $\alpha$-stable subordinator. For the sake of clarity we here report a sketched proof of this fact, which essentially follows [29, Theorem 2.2]. Let $H$ and $L$ respectively denote the $\alpha$-stable subordinator and its inverse, i.e.,

[TABLE]

In order to prove that (2.11) is the same process as $M(L(t))$ it is sufficient to prove that $M(L(t))$ has the same Mittag-Leffler intertimes of (2.11). This is clear since to construct a semi-Markov process in the sense of Gihman and Skoroohod, as in Section 2, it is sufficient to have an embedded chain $X_{n}$ and a sequence of independent r.v.’s representing the holding times. Here $M(t)$ and (2.11) have the same embedded chain and thus it remains only to show that they have the same waiting times. Since

[TABLE]

we have

[TABLE]

which is equivalent to (by [29, Lemma 2.1])

[TABLE]

Further, by [1, Lemma 2.3.2] we have that $H(V_{n}-)=H(V_{n})$, a.s., we can rewrite (2.28) as

[TABLE]

Thus the jump times $\tau_{n}$ of $M(L(t))$ are such that $\tau_{n}\stackrel{{\scriptstyle d}}{{=}}H(V_{n})$ and since $H$ has stationary increments, the holding times of $M(L(t))$ become, for any $n$

[TABLE]

where we used that $V_{n+1}-V_{n}$ are exponential r.v.’s $E_{n}$. By a standard conditioning argument we have, under $P\left(\cdot\mid X_{n}=x\right)$,

[TABLE]

Now since [30, eq. (3.4)]

[TABLE]

we have,

[TABLE]

and this shows that the holding times have the same distribution.

Example 2.3** (The fractional Poisson process).**

One of the most popular CTRW with heavy tailed waiting times is the so-called Fractional Poisson process, corresponding to the case where $\lambda_{i}=\lambda$, $h_{i,i+1}=1$ and $X(0)=0$ a.s. It has been studied by a number of authors (see for example [3; 22; 24; 29; 40]).

Its transition probabilities (2) solve the system of fractional Kolmogorov “forward” equations

[TABLE]

as well as the system of Kolmogorov fractional backward equations

[TABLE]

and it is easy to check directly that their common explicit solution in Laplace space is

[TABLE]

However, the equation often reported in the literature (e.g. [4]) is

[TABLE]

which is the “forward” equation corresponding to the special case $i=0$. In [29], the authors proved that the fractional Poisson process can be constructed as a standard Poisson process with the time variable replaced by an inverse stable subordinator.

3. Semi-Markov models for motion in heterogeneous media

We now show how the tools used for modeling homogeneous media can be adapted to include heterogeneity, in the sense that the trapping effect exercised in different sites depends on the current position. To be consistent with the literature introduced in Section 2 we first focus on the case in which the holding time at the position $x$ has a density $\psi(t)\sim t^{-\alpha(x)-1}$. How this can be generalized to different decaying patterns will be showed later. Hence we consider now a CTRW defined exactly as in (2.11), except for the distribution of the waiting times, which here present a state dependent fractional order:

[TABLE]

Use again [42, eq. (26)] to say that, for a constant $C>0$ depending on $\alpha_{i}$ and $\lambda_{i}$, we have, as $t\to\infty$

[TABLE]

and thus this is a model of a motion performed in a medium where the trapping effect has not the same intensity at all sites.

Before moving to the equation, it is usefull to show that also in this situation it is possible to interpret the semi-Markov process $X(t)$ as the time-change of a Markov process. However this is far from trivial and requires some analysis which is carried out in the following section.

3.1. The time-change by a dependent time process

In order to have an interpretation of (3.1) as a time-changed process, we need the notion of multistable subordinator (see for example [35] and [37]). Strictly speaking, a multistable subordinator is a generalization of a stable subordinator, in the sense that the stability index is a function of time $\alpha=\alpha(t)\in(0,1)$. The intensity of jumps is described by a time-dependent Lévy measure

[TABLE]

A multistable subordinator $\sigma(t)$, $t\geq 0$ is an additive process in the sense of [41], i.e., it is right-continuous and has independent but non stationary increments. Hence all the finite dimensional distributions are completely determined by the distribution of the increments which can be obtained from (3.3). Therefore (see [37, Section 2] for details on this point)

[TABLE]

Multistable subordinators are particular cases of a larger class of processes, known as non-homogeneous subordinators, which were introduced in [37].

Definition 3.1**.**

A multistable subordinator $\sigma(t)$, $t\geq 0$, is said to be piecewise stable if there exists a sequence $\alpha_{j}\in(0,1)$ and a sequence $t_{j}\geq 0$ such that the stability index can be written as

[TABLE]

and thus the time-dependent Lévy measure has the form

[TABLE]

Note that to each $\alpha_{j}\in(0,1)$ there corresponds a stable subordinator $H_{\alpha_{j}}(t)$ with index $\alpha_{j}$ in such a way that $\sigma$ is defined as

[TABLE]

The following theorem shows that (3.1) is given by a Markov process time-changed by the inverse of a piecewice stable subordinator. The major novelty consists in the fact that the original process and the random time are not independent as in the classical case. This reflects the fact that the intensity of the trapping effect is not space homogeneous, i.e., the time delay depends on the current position.

Theorem 3.2**.**

Let $M$ be a Markov process defined as in (2.1). Moreover, let $\sigma^{M}(t)$ be a multistable (piecewise stable) subordinator dependent on $M$ whose Lévy measure is given, conditionally on $V_{1}=v_{1},V_{2}=v_{2},\cdots$ and $X_{1}=x_{1},X_{2}=x_{2},\cdots$ by

[TABLE]

Let $L^{M}(t)$ be the right-continuous process

[TABLE]

Then the time-changed process $M(L^{M}(t))$ is the same process as (3.1).

Proof.

The proof is on the line of the discussion at the end Section 2. To prove that (3.1) coincides with $M(L^{M}(t))$ it is sufficient to prove that $M(L^{M}(t))$ has the same Mittag-Leffler intertimes of (3.1) since $M(t)$ and $X(t)$ have the same embedded chain. Let $V_{n}$, $n\geq 1$, be the jump times of $M$. Since

[TABLE]

we have

[TABLE]

Now by [37, Theorem 2.2] we know that $\sigma^{M}(t)$ is strictly increasing and then we can apply [29, Lemma 2.1] to say that (3.10) is equivalent to

[TABLE]

Now use [37, Theorem 2.1] to say that, a.s., $\sigma^{M}(t)=\sigma^{M}(t-)$ and thus to say that (3.12) is equivalent to

[TABLE]

Thus the jump times $\tau_{n}$ of $M(L^{M}(t))$ are such that $\tau_{n}\stackrel{{\scriptstyle d}}{{=}}\sigma^{M}(V_{n})$ and, by (3.7), the holding times are such that, under $P\left(\cdot\mid X_{n}=x\right)$,

[TABLE]

By a standard conditioning argument we have

[TABLE]

This fact together with formula (2.33) concludes the proof. ∎

3.2. Variable order backward equations

We here derive the backward equation for the semi-Markov process (3.1) in this new heterogeneous framework and we show that this equation becomes fractional of order $\alpha(i)$ where $i$ is the state where the transition is started.

Theorem 3.3**.**

The transition functions of (3.1) $p_{i,j}(t)$, $i,j\in\mathcal{S}$, solve the following system of backward equations

[TABLE]

Proof.

We can perform Laplace transform computation similar to that in the proof of Proposition 2.1. By applying the Laplace transform to (3.16) we obtain

[TABLE]

Instead, by applying the Laplace tranform to (2.17) we have

[TABLE]

By setting $g_{i,j}=\lambda_{i}(h_{i,j}-\delta_{i,j})$, it is easy to show that (3.18) can be rewritten as (3.17) and by the uniqueness theorem for Laplace tranform the desired result is immediate. ∎

The explicit form of the transition probabilities is easy obtained in Laplace space. By applying the Laplace transform, the system of fractional equations (3.16) reduces to the system of linear equations (3.17) in the variables $\widetilde{p}_{i,j}(s)$. In compact matrix form, (3.17) can be written as

[TABLE]

where $(\widetilde{P}(s))_{i,j}=\widetilde{p}(s)_{i,j}$, $I$ is the identity matrix, while

[TABLE]

Thus the solution in matrix form is written as

[TABLE]

4. The forward equations of semi-Markov processes in heterogeneous media

In the spirit of what happens in the homogeneous case one can be tempted to look for the forward equation by trying to replace the ordinary time-derivative in (2.9) with a variable-order Caputo derivative $\mathcal{D}_{t}^{\alpha(\cdot)}$, where $(\cdot)$ denotes the final state $j$. However, such an attempt is unsuccessful since it can be shown that the solution in the Laplace space does not coincide with the one of (3.16). We discuss here an example from the literature in which it is showed that this approach fails.

Example 4.1** (The state dependent fractional Poisson process).**

In the pioneering work [15], the authors studied a generalization of the fractional Poisson process in which the waiting times are independent but not identically distributed. For a given sequence $J_{n}$, $n\geq 0$, of independent r.v.’s with distribution

[TABLE]

they defined the state dependent fractional Poisson process as

[TABLE]

where $T_{n}=\sum_{k=0}^{n-1}J_{k}$, $T_{0}=0$. Further they proved that the state probabilities $p_{k}(t):=P(\mathcal{N}(t)=k\mid\mathcal{N}(0)=0)$ are such that

[TABLE]

The authors noticed that apparently $\mathcal{N}(t)$ is not governed by fractional differential equations, since the state probabilities corresponding to (4.3) do not solve the fractional “forward” equation with variable order derivative

[TABLE]

Moreover, the construction of $\mathcal{N}(t)$ as a time-changed process was not clear.

An application of our results to this particular situation also sheds light to such problems. In particular, from Theorem 3.3 it follows that the transition probabilities, which have explicit Laplace transform

[TABLE]

are really related with fractional calculus and indeed they solve the following system of fractional Kolmogorov backward equations

[TABLE]

Moreover, by Theorem 3.2 it follows that to obtain (4.2), a standard Poisson process must be composed with a dependent multistable subordinator. In what follows we determine the structure of the forward equations for semi-Markov processes defined as in (3.1). The state-dependent fractional Poisson process is therefore a particular case, so a fractional forward equation for this process can be really written down and does not coincide with (4.4).

We observe that in this framework the authors of [5] also studied the time-change of a Poisson process by means of the inverse of a multistable subordinator. However in this case the multistable subordinator is assumed independent from the Poisson process. Hence the resulting process is a random walk with independent but non identically distributed sojourn times, whose probability law is a time-inhomogeneous generalization of the Mittag-Leffler distribution, which is not coinciding with (4.1). Hence, the process in [5] is not a model for motions in heterogeneous media, but rather in an environment whose physical conditions change over time.

We now derive the system of forward (or Fokker-Planck) equations governing the process (3.1). Since the random walk has place in a heterogeneous medium, an adequate kinetic description of the process requires variable order fractional operators. The proof of the following theorem is based on the quantities

[TABLE]

where $J_{i}^{+}$ and $J_{i}^{-}$ represent the gain and the loss fluxes for the state $i$. It is intuitive that to deal with $J_{i}^{+}$ and $J_{i}^{-}$ it is convenient to assume that the when the process jumps (so when $t=T_{n}$ for some $n$) it can’t jump in its current position. Hence we will assume in the following theorem that $h_{i,i}=0$. However this is not strictly necessary and by adapting the notations one can generalize and remove the assumption.

Remark 4.2**.**

The Markov property is a consequence of the lack of memory property of the exponential distribution, which roughly states that when the Markov process is at $i$ the probability of having a jump in an infinitesimal interval of time $dt$ is $\lambda_{i}dt$, so the rate $\lambda_{i}$ is constant in time. Of course since the lack of memory is not true for other distributions, we must have here that the rate varies with time, i.e., it is $\lambda_{i}u_{i}(t)dt$. It turns out that in this case the function $u_{i}(t)$ is given by

[TABLE]

This fact will be proved in the general situation in Section 6. We remark that the probability that in an infinitesimal interval there is more than one jump is $o(dt)$. This is because by construction we know that $X(t)$ is the same process as $M(t^{\prime})$ with $t=\sigma^{M}(t^{\prime})$: since $\sigma^{M}(t^{\prime})$ is, a.s., strictly increasing on any finite interval of time and continuous, we have that within an infinitesimal time interval $dt$ also the process $X(t)$ performs at most one jump.

Theorem 4.3**.**

Let $X$ be the process in (3.1). Assume that the $h_{i,i}=0$. Then the transition probabilities

[TABLE]

solve the following system of fractional forward equations

[TABLE]

Proof.

The probability that the process performs more than one jump in an infinitesimal interval is $o(dt)$ is discussed in Remark 4.2. Let $J_{i}^{+}(t)dt$ be the probability of reaching the state $i$ during the time interval $[t,t+dt)$, i.e.,

[TABLE]

and let $J_{i}^{-}(t)dt$ be the probability of leaving the state $i$ during the time interval $[t,t+dt)$, i.e.,

[TABLE]

Then we have under $P:=P\left(\cdot\mid X(0)=l,\gamma(0)=0\right)$

[TABLE]

which, in our notations, reads

[TABLE]

or, equivalently,

[TABLE]

By the total probability law, the ingoing flux can be computed as

[TABLE]

where $h_{k,i}$ is the matrix of the embedded chain. Then we obtain the following balance equation (expressing the conservation of probability mass)

[TABLE]

which can also be written as

[TABLE]

The main goal is to compute the outgoing flux $J_{i}^{-}(t)$. It can be viewed as the sum of two contributions

[TABLE]

where $A_{i}^{1}(t)dt$ is the probability to be initially in the state i and to remain there for a time exactly equal to $t$ and $A_{i}^{2}(t)dt$ is the probability to reach the state i at time $t^{\prime}<t$ and to remain there for a time $t-t^{\prime}$. Thus

[TABLE]

where $p_{l,i}(0)=\delta_{l,i}$ and $f_{i}(t)$ the probability density of the holding time in $i$. We can eliminate $J^{+}_{i}(t)$ since by (4.10) we have $J^{+}_{i}(t)=\frac{d}{dt}p_{l,i}(t)+J^{-}_{i}(t)$. We thus obtain

[TABLE]

which is an integral equation in $J^{-}_{i}(t)$. By applying the Laplace transform to (4.12) we obtain

[TABLE]

which gives

[TABLE]

By assuming that the holding times follow a Mittag-Leffler distribution we have by (2.33) that

[TABLE]

and thus formula (4.13) becomes

[TABLE]

By reminding the definition of Riemann-Liouville derivative (1.7), we have

[TABLE]

and (4.11) reduces to

[TABLE]

and the proof is complete. ∎

Remark 4.4**.**

We remark that the reason why such equation is called forward is that the operator on the right side acts on the “forward” variable $i$ but leaves unchanged the backward variable $l$. Indeed, the equation is derived by conditioning on the event of the last jump (reaching the final state $i$) that may have occurred in a narrow interval near $t$. This last aspect is particularly clear by looking at (4.9).

Remark 4.5**.**

If $X$ is a Markov process, then the sojourn times are exponentially distributed, i.e.

[TABLE]

Then (4.13) reduces to

[TABLE]

namely

[TABLE]

Then the balance equation (4.11) reduces to the forward Kolmogorov equation

[TABLE]

From the physical point of view, the dynamics of the Markovian case and that of the CTRW with Mittag-Leffler waiting times present a wide difference. Indeed, in the Markov case, the outgoing flux $J_{i}^{-}(t)$ from the state $i$ at time $t$ is proportional to the concentration of particles at the present time $t$ (see (4.16)). Instead, in the CTRW with power law waiting times, the outgoing flux depends on the particles concentration at past times, according to a suitable weight kernel (see (4.14)).

Example 4.6** (The state dependent fractional Poisson process (continued)).**

To conclude the discussion of Example 4.1 we remark here that the forward equation for the state-dependent fractional Poisson process can be written down by using Theorem 4.3. We have indeed that the probabilities $p_{l,i}(t):=P\left(\mathcal{N}(t)=i\mid\mathcal{N}(0)=l\right)$ satisfy

[TABLE]

5. Convergence to the variable order fractional diffusion

Suppose that the state space $\mathcal{S}$ is embedded in $\mathbb{R}$. Hence our processes can be viewed as processes on $\mathbb{R}$, whose distribution is supported on $\mathcal{S}$. So in this section we consider a suitable scale limit of a semi-Markov process and we show that the one time distribution converges to the solution of forward heat equation on $\mathbb{R}$ with fractional variable order

[TABLE]

Such equation has been derived for the first time in [8] exactly in the study of anomalous diffusion in heterogeneous media. Hence our method provides a semi-Markov framework to this equation. The homogeneous case is represented by the time-fractional diffusion equation

[TABLE]

which is well known in literature and already has a probabilistic interpretation (see Remark 5.4 below for some details). This equation is related with anomalous diffusion (non-Fickian diffusion), see, for example, [21] for a recent application.

Let us assume that the process defined in (3.1) is a symmetric CTRW with Mittag-Leffler waiting times and with transition probabilities

[TABLE]

and $\lambda_{i}=\lambda$. Since $i,j$ are labels for points on the real line we can safely assume that the distance in $\mathbb{R}$ of two near points of $\mathcal{S}$ is constant and equal to $\epsilon$, i.e., $|i-j|=\epsilon$ for $j=i-1,i+1$ where $|\cdot|$ is the euclidean distance in $\mathbb{R}$. Hence looking at the process in $\mathbb{R}$ the walker performs jumps of size $\epsilon$. Then define

[TABLE]

The forward equation (4.8) reduces to

[TABLE]

By considering now the auxiliary function

[TABLE]

we can rewrite (5.3) as

[TABLE]

By setting $\lambda=1/\epsilon^{2}$ and letting $\epsilon\to 0$, the second derivative $\frac{\partial^{2}}{\partial y^{2}}u(x,y,t)$ arises. We thus obtain

[TABLE]

Note that the same scaling limit of a symmetic CTRW on a $d$-dimensional lattice leads to an analogous equation exibiting the Laplace operator in place of the second order derivative.

Equation (5.1) can be obtained phenomenologically by combining the continuity equation

[TABLE]

with an ad-hoc fractional Fick’s law regarding the flux $q(x,y,t)$:

[TABLE]

The fractional derivative in this expression provides a weighted average of the density gradient over the prior history, provided that the kernel of the average depends on the position $y$.

In terms of probability theory, the picture is completed by the backward heat equation with fractional variable order

[TABLE]

Such equation has been derived in [38], where the authors studied the convergence of the resolvent of semi-Markov evolution operators. Heuristically, eq. $\eqref{H}$ can be obtained from the backward equation (3.16) adapted to the case of a symmetric random walk with Mittag-Leffler waiting times:

[TABLE]

Indeed, passing to a lattice of size $\epsilon$ we have

[TABLE]

By assuming $\lambda=1/\epsilon^{2}$, the limit $\epsilon\to 0$ gives the desired result.

Remark 5.1**.**

In our derivation of the fractional heat equations, the diffusion coefficient is put equal to 1. However, the equations reported in the literature (e.g. [14; 8]) usually exhibit a space-dependent diffusion coefficient. To obtain this from a CTRW scheme, it is sufficient to assume a space dependent intensity $\lambda$ such that, in the limit of small $\epsilon$, it is of order $1/\epsilon^{2}$, that is $\lambda(x)=\widetilde{\lambda}(x)/\epsilon^{2}$. Then, it is natural to define the diffusion coefficient as

[TABLE]

and to repeat the same scaling limit argument in order to obtain both the forward equation

[TABLE]

and the backward equation

[TABLE]

Remark 5.2**.**

The mean square displacement of a subdiffusion in a homogeneous medium grows slower with respect to the Brownian motion, i.e., $\overline{x^{2}}(t)\sim t^{{\alpha}}$, where $\alpha\in(0,1)$. Such a process can be represented as a Brownian motion delayed by an independent inverse stable subordinator. For a study of its long time asymptotic properties, consult [25]. Concerning subdiffusion in heterogeneous media described by our equation, the picture is much more complicated and some unexpected phenomena arise. For example, in [14], the authors find that in the long time limit the CTRW process is localized at the lattice point where $\alpha(x)$ has its minimum, a phenomenon called “anomalous aggregation”. This suggests that the process does not enjoy the same ergodic properties reported in [25] for subdiffusion in homogeneous media.

Remark 5.3**.**

The transition probability $p(x,y,t)$ is the fundamental solution to the partial differential equation (5.9), and thus it is interesting to consider some well-posedness issues. The most recent result in this direction can be found in [17] where the authors consider the following Cauchy problem

[TABLE]

and prove that, under suitable assumptions on the source term $f$ and the initial datum $u_{0}$, there exists a unique weak solution $u(x,t)$ in the sense of [17, Thm 2.3].

Remark 5.4**.**

There are several results on CTRWs limit processes which can be applied in this situation by making some further assumptions [26; 28; 31; 32; 43]. We discuss here an example. Consider the couple process $\left(X_{n},T_{n}\right)$ and introduce a scale parameter $c$, so the process is $\left(X_{n}^{c},T_{n}^{c}\right)$. By making assumptions on the weak convergence of probability measures of the process $\left(X_{[u/c]}^{c},T_{[u/c]}^{c}\right)\to\left(A(t),D(t)\right)$ as $c\to 0$ (e.g. as in [43, Theorem 3.6]) one has that

[TABLE]

where $E(t)$ is the hitting-time of $D$. In our situation the processes are non independent and hence “coupled” in the language CTRW. Of course when $A$ is a Brownian motion and $D(t)$ is an independent $\alpha$-stable subordinator we are in the equivalent homogeneous situation of this section: the one time distribution of $A(E(t))$ solves indeed eq (5.2) in which $\alpha$ is constant [2]. Here we can conjecture that in order to obtain a process governed by (5.1) we must assume that $A$ is still a Brownian motion and that $D(t)$ is a multistable subordinator $\sigma(t)$ obtained as a limit case of the piecewise stable subordinator of Definition (3.1). So we argue that $D(t)$ must be a multistable subordinator (dependent on $A(t)$) whose Lévy measure is the limit of the Lévy measure of a piecewise stable subordinator, i.e., conditionally on a Brownian path $A(t,\omega)$

[TABLE]

Then one can define $E(t)$ as the hitting-time of $\sigma(t)$. This require several further investigations.

6. Arbitrary holding times and integro-differential Volterra equations

The construction of Theorem 3.2 is based on the notion of multistable subordinator. In [37] the authors introduced the more general class of inhomogeneous subordinators, i.e., non decreasing processes with independent and non stationary increments. By using these, it is possible to define a new type of CTRWs, which is constructed in the same way as (3.1), except for the distributions of the waiting times, which are no more Mittag-Leffler.

Indeed, for any $i\in\mathcal{S}$, consider a Lévy measure $\nu(dx,i)$ which defines a homogeneous subordinator $\sigma^{i}$ such that

[TABLE]

where

[TABLE]

is the Laplace exponent of $\sigma^{i}$. Let $L^{i}(t)=\inf\left\{\tau:\sigma^{i}(\tau)>t\right\}$ be the right continuous hitting time of $\sigma^{i}$. For any $i\in\mathcal{S}$ we assume $\nu((0,\infty),i)=\infty$, in such a way that $\sigma^{i}$ is a.s. stricly increasing, $L^{i}$ has a.e. continuous sample paths, and, for any $t>0$, $\sigma^{i}(t)$ and $L^{i}(t)$ are absolutely continuous random variables. We are now ready to define the following CTRW:

[TABLE]

where $T_{n}=\sum_{k=0}^{n-1}J_{k}$, $T_{0}=0$, and

[TABLE]

The generalization of Theorem (3.2) is immediate. Let $M$ be a Markov process defined as in (2.1). Moreover, let $\sigma^{M}(t)$ be an inhomogeneous subordinator dependent on $M$ whose Lévy measure, conditionally on $V_{1}=v_{1},V_{2}=v_{2},\cdots$ and $X_{1}=x_{1},X_{2}=x_{2},\cdots$ is given by

[TABLE]

Let $L^{M}(t)$ be the right continuous inverse of $\sigma^{M}(t)$. Then the time-changed process $M(L^{M}(t))$ is the same process as (6.2). To prove this, the key point is the fact that

[TABLE]

namely, conditionally to $X_{n}=i$, $J_{n}$ has a density $\psi_{i}$ with Laplace transform

[TABLE]

which is coinciding with $\mathds{E}(e^{-s\sigma^{i}(E_{n})}|X_{n}=i)$.

6.1. Integro-differential Volterra equations with position dependent kernel

To obtain a backward equation, we resort again to (2.17) and applying Laplace transform to both sides yields

[TABLE]

which can be rearranged as

[TABLE]

where again $g_{i,j}=\lambda_{i}(h_{i,j}-\delta_{i,j})$. Inverting Laplace transform in (6.8) however does not yield to a time-fractional equation. By using indeed [44, Lemma 2.5 and Proposition 2.7] we get the inverse Laplace transform

[TABLE]

where $\bar{\nu}(t,i):=\nu((t,\infty),i)$, provided that the integral function is differentiable. It is clear that in the situation of Theorem 3.2 one has

[TABLE]

and the operator on the left-hand side of (6.9) becomes a fractional Caputo derivative.

We now also derive a forward equation. Let

[TABLE]

be the number of renewals for the process (6.2) up to time $t$. Of course, conditionally to $X_{1}=x_{1}$, $X_{2}=x_{2},\cdots,$ we have that $N^{*}$ is a birth process (with rates $\lambda_{x_{1}}$, $\lambda_{x_{2}},\cdots$) time changed by the dependent time process $L^{M}$. Our attention focuses on the quantity, conditionally on $\left\{X(t)=i\right\}$,

[TABLE]

which we call renewal density (in the spirit of [10, page 26]) and specifies the mean number of renewals to be expected in a narrow interval near $t$ conditionally on the current position. Since we condition on $X(t)=i$, (6.12) is obviously depending on $i$ and $N^{*}(t+dt)-N^{*}(t)$ behaves like $N(L^{i}(t+dt))-N(L^{i}(t))$ where $N$ is a standard birth process. Thus the limit (6.12) can be computed as

[TABLE]

The function

[TABLE]

on the right-hand side of (6.13) is said to be, in the language of potential theory (e.g. [7]), the potential density of the subordinator $\sigma^{i}$ and is such that (e.g. [6, Section 1.3])

[TABLE]

provided that the derivative in (6.14) exists a.e. Heuristically $u^{i}(t)$ represents the mean of the total amount of time spent by the subordinator $\sigma^{i}$ in the state $dt$.

To obtain a forward equation we can follow the same line of section 4 up to formula (4.13). Then, by using (6.6), the outgoing flux has Laplace transform

[TABLE]

and thus the convolution theorem gives

[TABLE]

Finally, (4.11) reduces to

[TABLE]

which is the forward equation for our process. It is straightforward to prove that in the fractional case the renewal density relatively to the state $k$ reads

[TABLE]

and the operator on the right-hand side of (6.18) reduces to the Riemann Liouville derivative ${}^{R}\mathcal{D}^{1-\alpha_{k}}$. With the above discussion we have proved the following result.

Theorem 6.1**.**

Let $X(t)$ be a process like (6.2) with holding times $\overline{F}_{i}(t)$ given by (6.3). Further assume that $m_{i}(t)$ exists for any $i$. Then the probabilities $p_{i,j}(t)$ satisfy the backward equation (6.9) as well as the forward equation (6.18).