Biased Continuous-Time Random Walks with Mittag-Leffler Jumps

Abstract

We construct admissible circulant Laplacian matrix functions as generators for strictly increasing random walks on the integer line. These Laplacian matrix functions refer to a certain class of Bernstein functions. The approach has connections with biased walks on digraphs. Within this framework, we introduce a space-time generalization of the Poisson process as a strictly increasing walk with discrete Mittag-Leffler jumps time-changed with an independent (continuous-time) fractional Poisson process. We call this process ‘space-time Mittag-Leffler process’. We derive explicit formulae for the state probabilities which solve a Cauchy problem with a Kolmogorov-Feller (forward) difference-differential equation of general fractional type. We analyze a “well-scaled” diffusion limit and obtain a Cauchy problem with a space-time convolution equation involving Mittag-Leffler densities. We deduce in this limit the ‘state density kernel’ solving this Cauchy problem. It turns out that the diffusion limit exhibits connections to Prabhakar general fractional calculus. We also analyze in this way a generalization of the space-time Mittag-Leffler process. The approach of constructing good Laplacian generator functions has a large potential in applications of space-time generalizations of the Poisson process and in the field of continuous-time random walks on digraphs.

1. Introduction

During the last three decades a vast interest in fractional calculus has emerged due to its success to describe so called ‘anomalous phenomena’ such as anomalous transport and diffusion [1,2,3,4,5,6,7], anomalous relaxation in dielectrics [8], creep models [9], and various complex phenomena [10], just to quote a few examples. Many of these models can be traced back to the Montroll-Weiss picture of continuous-time random walk (CTRW) [11]. Whereas classical CTRW models are random walks subordinated to an independent Poisson process, their fractional generalizations lead to fat-tailed waiting-time Mittag-Leffler type densities with non-Markovian long memory behavior governed by evolution equations of fractional or generalized fractional types [1,4,6,12,13]. For fundamental discussions on fractional calculus we invite the reader to consult the references [14,15,16,17]. In the meantime, several generalizations of fractional calculus have been proposed. One of the most pertinent generalizations seems to be the so called Prabhakar generalization [18,19,20]. This generalization involves the Prabhakar function in place of the Mittag-Leffler function and was first introduced by Prabhakar [21]. In the context of continuous-time renewal processes a Prabhakar generalization of the fractional Poisson process was introduced by Cahoy and Polito [22], applied to stochastic motions on undirected graphs by Michelitsch and Riascos [23,24,25], and discrete-time versions were developed in our recent paper [26]. The present paper is mainly concerned with space-time generalizations of the Poisson process which can be traced back to special kinds of biased walks, namely strictly increasing walks on the integer line ${\mathbb{N}}_0$ subordinated to independent renewal processes. Such processes turn out to be governed by Cauchy problems with space-time difference-differential equations of general fractional type (generalized Kolmogorov-Feller or Kolmogorov-forward equation),

$$\begin{array}{ccc}{\displaystyle {\mathcal{D}}_t{p}_n\left(t\right)}& {\displaystyle =-\frac{\xi }{g\left(1\right)}g(1-{\hat{T}}_{-1})p_n\left(t\right)}\hfill & \\ \\ {\displaystyle p_n{\left(t\right)|}_{t=0}}& {\displaystyle ={\delta }_{n0}}\hfill \end{array}\xi >0,t\ge 0,n\in {\mathbb{N}}_0,$$

(1)

for the state probabilities $\left\{p_n\left(t\right)\right\}$ ($n=0,1,2,\dots \in {\mathbb{N}}_0$) of the process. The state probabilities $p_n\left(t\right)=\mathbb{P}(\mathcal{N}\left(t\right)=n)$, $n\in {\mathbb{N}}_0$, $t\in {\mathbb{R}}^{+}$ indicate the probabilities to find a walker (with the indicated initial condition) at time t in state $\mathcal{N}\left(t\right)=n$ in a strictly increasing walk $\mathcal{N}\left(t\right)={\sum }_{j=1}^{N\left(t\right)}{Z}_j\in {\mathbb{N}}_0$ with IID strictly positive integer jumps $Z_j\in \mathbb{N}$, where $N\left(t\right)\in {\mathbb{N}}_0$ denotes the number of arrivals up to time t in a renewal process independent of the jumps $Z_j$. The state probability distribution is normalized on the state space ${\sum }_{n=0}^{\infty }{p}_n\left(t\right)=1,\forall t\ge 0$. On the left-hand side of (1), ${\mathcal{D}}_t$ stands for a general fractional derivative which is related to the mentioned renewal process (see [27] for an outline of the framework of general fractional calculus). ${\hat{T}}_{-1}$ indicates the backward shift operator where shift operators ${\hat{T}}_m$ ($m\in \mathbb{Z}$) act on the state space such that ${\hat{T}}_m{p}_n\left(t\right)=p_{n+m}\left(t\right)$ (we denote with 1 the unity operator ${\hat{T}}_0=1$). The ‘Laplacian operator function’ $g(1-{\hat{T}}_{-1})$ represents the generator of the process. The stochastic process governed by Cauchy problem (1) is a space-time generalization of the Poisson process. Indeed for $g(1-{\hat{T}}_{-1})=1-{\hat{T}}_{-1}$ (with $(1-{\hat{T}}_{-1})p_n\left(t\right)=p_n\left(t\right)-p_{n-1}\left(t\right)$) and ${\mathcal{D}}_t=\frac{d}{dt}$, Equation (1) recovers the classical Kolmogorov-Feller (forward) equation of the homogeneous Poisson process. One of the goals of the present paper is to elaborate a general approach to construct non-trivial generalizations by means of “good Laplacian operator functions” $g(1-{\hat{T}}_{-1})$.

We will show that under suitable “well-scaled” continuous-space limit conditions from the integer-state space to a continuous one (i.e., from $\{0,1,\dots \}={\mathbb{N}}_0$ to $[0,\infty )={\mathbb{R}}^{+}$) the discrete Cauchy problem (1) turns into a Cauchy problem with a diffusion equation of the following general structure:

$$\begin{array}{ccc}{\displaystyle {\mathcal{D}}_t\mathcal{P}(x,t)}& {\displaystyle =-\xi \mathcal{P}(x,t)+\xi {\int }_0^x\mathcal{W}(x-\tau )\mathcal{P}(\tau ,t)\mathrm{d}\tau ,\xi >0,x,t\ge 0}\hfill & \\ \\ {\displaystyle {\mathcal{P}(x,t)|}_{t=0}}& =\delta \left(x\right)\hfill \end{array}$$

(2)

This equation only allows jumps into the positive x-direction, i.e., ‘strictly increasing’ walks. The convolution on the right-hand side describes incoming jumps from $[0,x)$ to state x whereas the term $-\xi \mathcal{P}(x,t)=-\xi \mathcal{P}(x,t){\int }_x^{\infty }\mathcal{W}(\tau -x,t)\mathrm{d}\tau$ accounts for outgoing jumps from x to $(x,\infty )$. Here $\mathcal{P}(x,t)$ stands for the ‘state density’, $\mathcal{W}\left(x\right)$ indicates the transition density kernel and $\delta \left(x\right)$ the Dirac’s $\delta$-distribution.

A further aim of the present paper is to highlight connections of these models with a certain class of biased random walks on directed graphs. This aim is especially motivated by the huge upswing of network science which has become nowadays an immense interdisciplinary field [28,29] last but not least driven by rapidly developing applications in online (social) networks and search engines. There is already a vast amount of specialized literature on various aspects of the subject, see, e.g., [28,30,31,32,33]. For instance, models have been developed to describe Lévy flight dynamics and random walks with long-range jumps on undirected graphs [34,35], random walks with stochastic resetting on graphs [36] and models for long-range mobility in cities [37], among many others.

The present paper is organized as follows. In Section 2 we recall some basic properties of biased walks on directed networks. We focus on ergodic (strongly connected) finite directed graphs and the features of their Laplacian and transition matrices. Having recalled these basic features, we define in Section 3 necessary and sufficient “good Laplacian properties” that a physically admissible ‘good Laplacian matrix’ must fulfill. Starting with a good Laplacian matrix allowing only ‘local’ transitions to next neighbor nodes we use these Laplacian properties to construct non-trivial good Laplacian functions (i.e., those that conserve the “good Laplacian properties”). They will constitute the family of admissible generators $g(1-{\hat{T}}_{-1})$ on the right-hand side of (1). It turns out that the family of good Laplacian functions refers to a certain class of Bernstein functions. For undirected graphs this approach was developed recently [34,38] and has also been extended to a class of biased walks [39].

In Section 4 we recall, within this picture, a class of biased walks on the integer line allowing solely strictly positive integer jumps generated by the fractional Laplacian matrix of a directed line. This leads to the Sibuya walk, which is characterized by a fat-tailed jump distribution, as a proto-typical example of a discrete approximation of the stable subordinator (see, e.g., [13,40,41]).

Section 5 is devoted to explore, by means of this approach, space-time generalizations of the Poisson process leading to Cauchy problems of the general structure (1) involving “good Laplacian matrix functions” of Section 3. Within this framework, we also discuss classical cases such as the space-fractional Poisson process and the space-time fractional Poisson process, introduced by Orsingher and Polito [42].

In the main part of this paper, Section 6, we apply this approach to construct the “space-time Mittag-Leffler process” as a pertinent generalization of the Poisson process governed by an evolution equation of general type (1). This process is a strictly increasing walk on the integer line with (discrete) Mittag-Leffler jumps separated by the Mittag-Leffler waiting times of the fractional Poisson process. We derive explicitly the state probabilities and the ‘well-scaled’ continuous-space limit (diffusion limit) of this process. We obtain a biased (forward) diffusion equation of general space-time fractional type referring to the class (2), which is solved by the continuous-space limit density kernel of the state-probabilities. This kernel is derived in explicit form involving so called Prabhakar kernels. In this way, we show connections with Prabhakar general fractional calculus. We also highlight the equivalence with the Montroll-Weiss CTRW picture.

Finally, as a byproduct, we construct in Section 7 the space-fractional generalization of the space-time Mittag-Leffler process and derive the Cauchy problem governing the state probabilities. The latter as well as its diffusion limit are derived in explicit forms where again Prabhakar kernels and general fractional calculus come into play.

2. Biased Walks on Directed Graphs

In this section we recall some basic notions and properties of biased walks taking place on directed graphs (also referred to as digraphs). We consider first a finite directed graph consisting of N nodes (states) which we denote by $j=1,\dots N$. In a directed graph the edges between nodes have in general direction so that a path may exist allowing the walk $m\to n$ but the return path $n\to m$ does not necessarily exist. If for all pairs of nodes finite paths (and hence return paths) via directed edges exist, then the digraph is said to be strongly connected and equivalently ergodic (see, e.g., [39,43]). In order to define random walks on digraphs we introduce the $N\times N$ Laplacian matrix $\mathbf{L}$ containing this topological information. The Laplacian matrix has the elements [39]

$$L_{ij}=k_i^{\left(\mathrm{out}\right)}{\delta }_{ij}-{\Omega }_{ij}$$

(3)

where we use the synonymous notation ${\delta }_{pq}={\delta }_{p,q}$ for the Kronecker symbol. Matrix (3) is a generalization of the Laplacian matrix for binary undirected networks [28,32,33] to include the possibility of weights in the connections and asymmetry in the flow along the edges. For these particular connections we have ${\Omega }_{ij}\ne {\Omega }_{ji}$ where we assume ${\Omega }_{ij}\ge 0$ is a non-negative matrix. Further, by construction we assume ${\Omega }_{ii}=0$. In (3) is present the so called out-degree defined by [39]

$$k_i^{\left(\mathrm{out}\right)}=\sum _{j=1}^N{\Omega }_{ij}>0$$

(4)

constituting a measure of the number of nodes which can be reached in a single jump from node i. The out-degree is assumed to be strictly positive meaning that we do not allow isolated (disconnected) nodes. Then, we introduce the one-step transition matrix $\mathbf{W}$ [28,29,34]

$$W_{i\to j}={\delta }_{ij}-{\mathcal{L}}_{ij}=\frac{{\Omega }_{ij}}{k_i^{\left(\mathrm{out}\right)}}$$

(5)

denoting the probability of the transition $i\to j$ in one jump. Per construction (5) is row-stochastic, i.e., $0\le W_{i\to j}\le 1$ with ${\sum }_{j=1}^N{W}_{i\to j}=1$ and the transition matrix is such that $W_{i\to i}=0$, i.e., in each jump the walker has to move to a different node. For our convenience and as an equivalent description we introduce here the (non-symmetric) auxiliary Laplacian matrix

$${\mathcal{L}}_{ij}=\frac{L_{ij}}{k_i^{\left(\mathrm{out}\right)}}$$

(6)

having normalized degrees ${\mathcal{L}}_{ii}=1$. We notice that dynamical processes in directed networks have a greater variety than in the undirected case. Since the Laplacian matrices (3) and (6) are not symmetric they have in general complex eigenvalues which are considered more closely in the subsequent Section 3. We can now define a discrete-time Markov chain (Markovian random walk) on this graph governed by the simple master equation (e.g., [29,34])

$$P_{ij}(t+1)=\sum _{r=1}^N{P}_{ir}\left(t\right)W_{r\to j},P_{ij}\left(0\right)={\delta }_{ij},t\in {\mathbb{N}}_0.$$

(7)

We assume here that the walker at $t=0$ is sitting on the departure node i. In the ‘t-step transition matrix’ the element $P_{ij}\left(t\right)$ indicates the probability of the transition $i\to j$ in t jumps. Note that for digraphs $\mathbf{P}\left(t\right)$ is not symmetric and for the initial condition $\mathbf{P}\left(0\right)=\mathbf{1}$ we have

$${\left[\mathbf{P}\left(t\right)\right]}_{ij}={\left[{\mathbf{W}}^t\right]}_{ij},t\in {\mathbb{N}}_0.$$

(8)

An important property of walks on strongly connected digraphs is (aperiodic) ergodicity. A Markov chain (8) is said to be ergodic if it exists $n_0\in \{1,2,\dots \}=\mathbb{N}$ such that

$${\left[\mathbf{P}\left(n\right)\right]}_{ij}>0,\forall i,j=1,\cdots ,Nn\ge n_0,$$

(9)

is strictly positive, i.e., for each pair of nodes $i\to j$ there is at least one connecting path (and return path) not longer than $n_0$ jumps via directed connecting edges. For a discussion of some aspects of ergodicity of biased walks in digraphs we refer to [39] (and see also the references therein).

For later use we consider a Montroll-Weiss CTRW where a biased walk with transition matrix (5) is subordinated to a counting process with state-probabilities $\mathbb{P}(N\left(t\right)=n)={\Phi }^{\left(n\right)}\left(t\right)$. Then the transition matrix (8) with the elements $P_{ij}\left(t\right)$, indicating the probability to find the walker at time t on node j (with the given initial condition $P_{ij}\left(0\right)={\delta }_{ij}$), is generalized by the Cox series [44]

$${\left[\mathbf{P}\left(t\right)\right]}_{ij}=\sum _{n=0}^{\infty }{\Phi }^{\left(n\right)}\left(t\right){\left[{\mathbf{W}}^n\right]}_{ij},t\in {\mathbb{R}}^{+}.$$

(10)

We notice that the discrete-time Markov chain (8) indeed is a special case of (10) when we account for a homogeneous event stream with state probabilities ${\Phi }^{\left(n\right)}\left(t\right)=\Theta (t-n)-\Theta (t-(n+1))$, where $\Theta \left(t\right)$ is the Heaviside step function which is such that $\Theta \left(t\right)=1$ for $t\ge 0$ and $\Theta \left(t\right)=0$ for $t<0$ (in particular $\Theta \left(0\right)=1$). Indeed (10) boils down to ${\left[\mathbf{P}\left(t\right)\right]}_{ij}\left(t\right)={\left[{\mathbf{W}}^n\right]}_{ij}$ and recovers then (8) as ${\Phi }^{\left(n\right)}\left(t\right)=1$ for $t\in [n,n+1)$ and ${\Phi }^{\left(n\right)}\left(t\right)=0$ otherwise.

3. Biased Walks with Long-Range Jumps

One goal of this paper is to introduce new types of biased walks with special emphasis on strictly increasing walks with long-range jumps together with a systematic method for their construction. To this end, we consider a strongly connected (ergodic) digraph with a Laplacian matrix $\mathbf{L}$ defined in (3) characterizing its topology. Now we seek “good” Laplacian matrix functions $\mathbf{L}\to g\left(\mathbf{L}\right)$ defining new topologies such that the following “good Laplacian properties’‘ are retained [34,38,39]:

(i)

${\sum }_{j=1}^N{L}_{ij}=0$ (i.e., the Laplacian matrix has a unique eigenvalue ${\mu }_1=0$ and a right eigenvector with constant (real) components $\langle i|v_1\rangle =1/\sqrt{N}$). This condition ensures that in a strongly connected graph the transition matrix has the unique real eigenvalue ${\lambda }_1=1$ of largest absolute value (reflecting the existence of a unique stationary distribution).

(ii)

$L_{ii}=k_i^{\left(\mathrm{out}\right)}>0$ ($\forall i=1,\dots N$) tells us that each node has neighboring nodes, i.e., there are no isolated nodes.

(iii)

$L_{ij}=-{\Omega }_{ij}\le 0$ for $i\ne j$, i.e., the off-diagonal elements are all non-positive.

Clearly, the conditions (i)–(iii) maintain stochasticity of the transition matrix (5) and hence of (10).

We notice that the auxiliary Laplacian (6) fulfills conditions (i)–(iii) with uniform normalized degree ${\mathcal{L}}_{ii}=1$$\forall i=1,\dots N$. Without loss of generality it is now more convenient to work with the auxiliary Laplacian (6) instead of the Laplacian matrix (3). We will prove subsequently that new non-trivial “good” Laplacian matrices retaining (i)–(iii) are obtained by the class of matrix functions [34,38,39]

$$g\left(\mathcal{L}\right)={\int }_0^{\infty }\left(\mathbf{1}-e^{-\tau \mathcal{L}}\right)\nu \left(\mathrm{d}\tau \right)$$

(11)

containing as argument the auxiliary Laplacian matrix $\mathcal{L}$ (6) and where $\mathbf{1}$ denotes the unity matrix. In the above integral, $\nu \left(\mathrm{d}\tau \right)$ indicates a Lévy measure for which (see e.g., [27])

$${\int }_0^{\infty }min\{1,\tau \}\nu \left(\mathrm{d}\tau \right)<\infty$$

(12)

ensuring convergence of (11). Let us consider for a moment a scalar version $g\left(\mu \right)$ ($\mu \ge 0$) of (11). Then, we can define a special class of Bernstein functions by the Lévy-Khintchine representation

$$g\left(\mu \right)={\int }_0^{\infty }\left(1-e^{-\tau \mu }\right)\nu \left(\mathrm{d}\tau \right),$$

(13)

which fulfills ${g\left(\mu \right)|}_{\mu =0}=0$. This property preserves the zero eigenvalue in the matrix function (11). In subsequent applications we consider Lévy measures $\nu \left(\mathrm{d}\tau \right)=\nu \left(\tau \right)\mathrm{d}\tau$ with non-negative Lévy densities $\nu \left(\tau \right)$. The Laplacian Bernstein function (13) in some part of the literature is also called Laplace exponent [45]. Generally, a non-negative function $f(a,b,\mu )\ge 0$ ($a,b\ge 0$) defined on $\mu \in [0,\infty )$ is said to be a Bernstein function if it fulfills ${(-1)}^{n-1}\frac{d^n}{d{\mu }^n}f(a,b,\mu )\ge 0$ ($n\ge 1$) and hence has the representation [27]

$$f_B(a,b,\mu )=a+b\mu +{\int }_0^{\infty }\left(1-e^{-\tau \mu }\right)\nu \left(\mathrm{d}\tau \right).$$

(14)

We consider Laplacian functions $g\left(\mu \right)$ (13) referring to the class of Bernstein functions with $a=b=0$, i.e., $g\left(\mu \right)=f_B(0,0,\mu )$. The condition $a=0$ is necessary to have ${g\left(\mu \right)|}_{\mu =0}=0$ for the sake of condition (i). The combination $a=0$, $b>0$ gives an admissible Laplacian function fulfilling (i)–(iii). However, it adds a ‘trivial’ contribution of the original Laplacian matrix $\mathcal{L}$. Therefore, we consider here only non-trivial matrix functions (13) with also $b=0$. We notice that the derivative $\frac{d}{d\mu }g\left(\mu \right)$ is a completely monotonic function, i.e.,

$${(-1)}^{n-1}\frac{d^n}{d{\mu }^n}g\left(\mu \right)\ge 0,n\in \mathbb{N}$$

(15)

with $g\left(\mu \right)>0$ for $\mu >0$ and $g\left(0\right)=0$. For theorems and a profound analysis of Bernstein and completely monotonic functions and related subjects, see [46,47]. We now prove that (11) retains the properties (i)–(iii). For a proof in undirected networks with symmetric Laplacian matrices, see [34,38] and for the fractional Laplacian matrix function on digraphs consult [39]. We introduce a constant $\Lambda >1$ such that the matrix $\Lambda \mathbf{1}-\mathcal{L}$ has uniquely non-negative matrix elements, namely $(\Lambda -1){\delta }_{ij}+\frac{{\Omega }_{ij}}{k_i^{\left(\mathrm{out}\right)}}\ge 0$. Then, for ergodic (i.e., strongly connected) digraphs it exists a $n_0>0$ such that all elements of the integer matrix power

$${\left[{(\Lambda \mathbf{1}-\mathcal{L})}^n\right]}_{ij}>0,n\ge n_0,\forall i,j=1,\cdots N$$

(16)

are strictly positive. Condition (16) can be used as an equivalent definition of ergodicity of a digraph. One can further infer that if such a finite $n_0$ exists, then ergodic digraphs must have a finite number N of states (the same holds also for undirected graphs [34]). Then, clearly the matrix exponential

$$e^{\Lambda \mathbf{1}-\mathcal{L}}=\sum _{n=0}^{\infty }\frac{1}{n!}{[\Lambda \mathbf{1}-\mathcal{L}]}^n$$

(17)

is strictly positive (we call a matrix ‘(strictly) positive’ if it has solely positive entries), as there are infinitely many uniquely positive matrices $\frac{1}{n!}{[\Lambda \mathbf{1}-\mathcal{L}]}^n$ contained in this series (namely those for $n\ge n_0$). Since $\mathbf{1}$ is commuting with any matrix $\mathcal{L}$ we have $e^{\Lambda \mathbf{1}-\mathcal{L}}=e^{\Lambda \mathbf{1}}·e^{-\mathcal{L}}$ thus ${\left[e^{\Lambda \mathbf{1}-\mathcal{L}}\right]}_{ij}=e^{\Lambda }{\left[e^{-\mathcal{L}}\right]}_{ij}>0$. Now, since $e^{\Lambda }>0$ the matrix exponential $e^{-\mathcal{L}}$ for an ergodic digraph is indeed strictly positive:

$${\left[e^{-\mathcal{L}}\right]}_{ij}>0,\forall i,j=1,\dots N.$$

(18)

On the other hand, because of condition (i) we have row-stochasticity

$$\sum _{j=1}^N{\left[e^{-\mathcal{L}}\right]}_{ij}=1,\forall i,j=1,\dots N$$

(19)

As a consequence of eigenvalue zero of $\mathcal{L}$ this relation follows from ${\sum }_{j=1}^N{\left[{\mathcal{L}}^n\right]}_{ij}={\delta }_{n0}$. Then since all matrix elements (18) are positive and row-normalized we have that

$$0<{\left[e^{-\mathcal{L}}\right]}_{ij}<1,\forall i,j=1,\dots N.$$

(20)

Hence, it follows that $0<1-{\left[e^{-\mathcal{L}}\right]}_{ii}<1$ thus the Bernstein matrix function $\mathbf{1}-e^{-\tau \mathcal{L}}$ ($\tau >0$) retains the conditions (i)–(iii) of a good Laplacian matrix. By similar considerations for non-ergodic digraphs one can show that the matrix exponential contains zero valued blocks of non-diagonal elements (where zero entries indicate the pairs of nodes where no finite directed paths exist). However, conditions (i)–(iii) in the Bernstein matrix function (11) still are retained [39].

Conditions (i)–(iii) remain also retained by the integration in (11) with a (non-negative) Lévy measure $\nu \left(\mathrm{d}\tau \right)$ if the integral converges. This is definitely the case by virtue of (12) and if and only if the eigenvalues of the auxiliary Laplacian $\mathcal{L}$ have solely non-negative real parts $\Re \left\{{\mu }_m\right\}\ge 0$. This indeed is true as we show by the following brief proof.

Denoting the eigenvalues of the auxiliary Laplacian (6) by ${\mu }_m$ with ${\mu }_1=0$ and the eigenvalues of the one-step transition matrix (5) by ${\lambda }_m$ we have

$${\mu }_m=1-{\lambda }_m,m=1,2,\dots ,N$$

(21)

where ${\mu }_1=1-{\lambda }_1=0$ (thus ${\lambda }_1=1$ which is the unique so called Perron-Frobenius eigenvalue [48]). Then, we have that ${lim}_{n\to \infty }{\mathbf{W}}^n={\mathbf{W}}^{\infty }$ remaining row-stochastic (corresponding to ${\lambda }_1=1$) which contradicts the existence of an exploding eigenvalue ${|\lambda |}^n\to \infty$. We can hence infer that for the complete set of complex eigenvalues it holds

$$|{\lambda }_m|\le 1,\forall m$$

(22)

and with $\Re \left\{{\lambda }_m\right\}\le |{\lambda }_m|\le 1$ we have that

$$\Re \left\{{\mu }_m\right\}=1-\Re \left\{{\lambda }_m\right\}\ge 1-|{\lambda }_m|\ge 0,\forall m.$$

(23)

We notice that in this proof we neither needed ergodicity nor that the digraph is finite, thus it also includes strictly increasing walks on the infinite integer line (see (43)).

This concludes our proof that the Lévy-Khintchine representation (11) indeed also converges for digraphs and is useful to construct good Laplacian matrix functions to define new biased walks with the one-step transition matrix

$$W_{ij}^{\left(g\right)}={\delta }_{ij}-\frac{g_{ij}\left(\mathcal{L}\right)}{g_{ii}\left(\mathcal{L}\right)}.$$

(24)

It follows from the above considerations (see especially (18)–(20) with (11)) for ergodic digraphs that all off-diagonal elements of (24) are strictly positive $W_{ij}^{\left(g\right)}>0$ ($i\ne j$) allowing the walker to reach any node in a single jump introducing new fully connected topologies. In non-ergodic digraphs some off-diagonal elements take values null due to the absence of finite directed paths [39]. We also notice that as $\mathcal{L}$ is not symmetric and in some cases non-diagonisable having Jordan canonic forms [49] where our proof remains valid also for these cases.

4. Circulant Transition Matrices and Strictly Increasing Walks

In this section we consider the class of biased walks on the integer line with IID strictly positive integer increments (‘jumps’) $Z_j>0$. Such a walk is defined by

$$Y_n=\sum _{j=1}^n{Z}_j,Z_j\in \mathbb{N},Y_0=0.$$

(25)

The random walk (25) is a discrete counterpart to a strict subordinator [41]. There is a path $a\to b$ only if $b>a$ but then no return path $b\to a$. Therefore in contrast to the walks on strongly connected (finite) digraphs, the walk (25) is not ergodic. As an example let us recall the Sibuya walk. The Sibuya distribution has generating function [41]

$${\overline{w}}_{\alpha }\left(u\right)=\mathbb{E}{u}^Z=1-{(1-u)}^{\alpha },\alpha \in (0,1),|u|\le 1.$$

(26)

In a Sibuya walk, the jumps $Z_j\in \mathbb{N}$ in (25) follow the Sibuya distribution (also known under the name Sibuya(α)) which is defined for $\alpha \in (0,1)$ as

$$\begin{array}{ccc}{\displaystyle \mathbb{P}(Z_j=k)=w_{\alpha }\left(k\right)=\frac{1}{k!}\frac{d^k}{d{u}^k}{\overline{w}}_{\alpha }\left(u\right){|}_{u=0}}& {\displaystyle ={(-1)}^{k-1}\frac{\alpha (\alpha -1)..(\alpha -k+1)}{k!}}\hfill & \\ \\ {\displaystyle }& {\displaystyle ={(-1)}^{k-1}\left(\genfrac{}{}{0pt}{}{\alpha }{k}\right)=\frac{\alpha }{k}\frac{\Gamma (k-\alpha )}{\Gamma (1-\alpha )\Gamma \left(k\right)},k\in \mathbb{N}}\hfill \end{array}$$

(27)

with $w_{\alpha }\left(k\right){{|}_{k=0}={\overline{w}}_{\alpha }\left(u\right)|}_{u=0}=0$ (crucial for the occurrence of strictly positive jumps $Z_j>0$) and $w_{\alpha }\left(k\right)>0$ for all $k\in \mathbb{N}$. In the limit $\alpha \to 1-$ the ‘trivial’ distribution $w_1\left(k\right)={\delta }_{1,k}$ is recovered (not called Sibuya) where the walker makes unit jumps $Z_j=1$ to its right-sided neighbor node almost surely. For later use we also mention the important property that Sibuya($\alpha$) is fat-tailed, namely for k large by using $\frac{\Gamma (k-\alpha )}{\Gamma \left(k\right)}\sim k^{-\alpha }$ we have the asymptotic behavior

$$w_{\alpha }\left(k\right)\sim \frac{\alpha }{\Gamma (1-\alpha )}{k}^{-\alpha -1}=-\frac{k^{-\alpha -1}}{\Gamma (-\alpha )}>0,(k\to \infty ).$$

(28)

For the following analysis it will be convenient to utilize the connection of generating functions and the shift operators. For a detailed outline we refer to our recent article [26] and see also Appendix A.

Let us introduce the shift operator ${\hat{T}}_a$ ($a\in \mathbb{R}$) which is such that ${\hat{T}}_{a}f\left(x\right)=f(x+a)$ and consider its circulant matrix representation

$${\hat{T}}_{m}f\left(p\right)=f(p+m)=:\sum _{q=-\infty }^{\infty }f\left(q\right){\left[{\hat{T}}_m\right]}_{qp}=\sum _{q=-\infty }^{\infty }f\left(q\right){\delta }_{q,p+m},p,q,m\in \mathbb{Z}$$

(29)

and hence ${\left[{\hat{T}}_m\right]}_{qp}={\delta }_{q,p+m}$. This circulant structure remains true for all matrices of operator functions of shift operators. All matrices $M_{mn}$ (we also write $M_{m,n}$) we are dealing with in the context of increasing walks (25) are characterized by the properties

$$\begin{array}{ccc}{\displaystyle M_{m,n}=M_{m+s,n+s}}& {\displaystyle =M_{0,n-m},n,m,s\in \mathbb{Z}}\hfill & \\ \\ {\displaystyle M_{m,n}=M_{0,n-m}}& {\displaystyle =0,n-m\in {\mathbb{Z}}_{-}}\hfill \end{array}$$

(30)

We call a matrix “upper triangular circulant” if it fulfills (30), i.e., all elements below the main diagonal are strictly null. The Sibuya transition matrix (33) is a proto-typical example for such an upper triangular circulant matrix.

We introduce the generating function of an upper triangular circulant matrix (30) as follows

$$\begin{array}{ccc}{\displaystyle \overline{M}\left(u\right)}& {\displaystyle =\sum _{k=-\infty }^{\infty }{M}_{0,k}{u}^k=\sum _{k=0}^{\infty }{M}_{0,k}{u}^k,}\hfill & {\displaystyle |u|\le 1}\\ \\ {\displaystyle M_{0,n}}& {\displaystyle =\frac{1}{n!}\frac{d^n}{d{u}^n}\hat{M}{\left(u\right)|}_{u=0},}\hfill & n\in {\mathbb{N}}_0\\ \\ M_{0,n}& {\displaystyle =0,}\hfill & n<0\end{array}$$

(31)

thus only non-negative powers $u^k$ ($k\ge 0$) are contained in $\hat{M}\left(u\right)$. Then we observe the important property which holds for upper triangular circulant matrices

$$\begin{array}{ccc}{\displaystyle \frac{1}{n!}\frac{d^n}{u^n}{\left\{{\overline{M}}^1\left(u\right){\overline{M}}^2\left(u\right)\right\}}_{u=0}}& {\displaystyle =\sum _{k=0}^n{M}_{0,n-k}^{\left(2\right)}{M}_{0,k}^{\left(1\right)}=\sum _{k=0}^n{M}_{0,k}^{\left(1\right)}{M}_{k,n}^{\left(2\right)}=\sum _{k=-\infty }^{\infty }{M}_{0,k}^{\left(1\right)}{M}_{k,n}^{\left(2\right)}}\hfill & \\ \\ & {\displaystyle ={\left[{\mathbf{M}}^{\left(1\right)}{\mathbf{M}}^{\left(2\right)}\right]}_{0,n}={\left[{\mathbf{M}}^{\left(2\right)}{\mathbf{M}}^{\left(1\right)}\right]}_{0,n}}\hfill \end{array}$$

(32)

connecting matrix multiplication and discrete (commuting) convolutions (See also [26]).

The operator that emerges when replacing in (26) u with ${\hat{T}}_{-1}$ is hence the upper triangular circulant one-step transition matrix of the Sibuya walk [39], namely

$$\begin{array}{ccc}{\displaystyle {\left[{\overline{w}}_{\alpha }({\hat{T}}_{-1})\right]}_{mn}={\left[{\overline{w}}_{\alpha }({\hat{T}}_{-1})\right]}_{0,n-m}}& =\hfill & {\displaystyle {[1-{(1-{\hat{T}}_{-1})}^{\alpha }]}_{m,n}={\delta }_{mn}-{\left[{\mathcal{L}}^{\alpha }\right]}_{m,n},\alpha \in (0,1]}\\ \\ {\displaystyle }& =\hfill & {\displaystyle \left\{\begin{array}{c}{\displaystyle \sum _{k=1}^{\infty }{(-1)}^{k-1}\left(\genfrac{}{}{0pt}{}{\alpha }{k}\right){\delta }_{m,n-k}={(-1)}^{n-m-1}\left(\genfrac{}{}{0pt}{}{\alpha }{n-m}\right),n-m>0}\hfill \\ \\ {\displaystyle 0,n-m\le 0}\hfill \end{array}\right.}\end{array}$$

(33)

where $m,n\in \mathbb{Z}$. We used the following properties of the (unitary) shift operator with ${({\hat{T}}_{-1})}^k={\hat{T}}_{-k}$ and (29) thus

$${\left[{\hat{T}}_{-k}\right]}_{m,n}={\delta }_{m,n-k}={\delta }_{k,n-m},m,n,k\in \mathbb{Z}$$

(34)

is for $k\ge 0$ upper triangular circulant with entries 1 in the kth upper side-diagonal ($n-m=k$), and 0 elsewhere. Formula (33) contains the upper triangular circulant fractional Laplacian matrix $g\left(\mathcal{L}\right)={\mathcal{L}}^{\alpha }={(1-{\hat{T}}_{-1})}^{\alpha }$ which has the elements

$${\left[{\mathcal{L}}^{\alpha }\right]}_{mn}={\left[{(1-{\hat{T}}_{-1})}^{\alpha }\right]}_{mn}=\left\{\begin{array}{c}{\displaystyle {(-1)}^{n-m}\left(\genfrac{}{}{0pt}{}{\alpha }{n-m}\right),n-m\in {\mathbb{N}}_0}\hfill \\ \\ {\displaystyle 0,n-m<0}\hfill \end{array}\right.\alpha \in (0,1].$$

(35)

For $\alpha =1$ the Laplacian matrix (from now on we refer the auxiliary Laplacian $\mathcal{L}$ to as Laplacian matrix)

$${\left[\mathcal{L}\right]}_{mn}={\left[(1-{\hat{T}}_{-1})\right]}_{m,n}={\delta }_{mn}-{\delta }_{m,n-1},m,n\in \mathbb{Z}$$

(36)

is recovered. It follows from Equation (31) that the generating function of Laplacian (36) yields $\overline{\mathcal{L}}\left(u\right)=1-u$. Thus the generating function of the fractional Laplacian matrix (35) takes the form ${\overline{\mathcal{L}}}^{\alpha }\left(u\right)={(1-u)}^{\alpha }$. The generating function of fractional Sibuya transition matrix (33) per construction recovers the generating function (26) of Sibuya($\alpha$).

Further instructive for subsequent use is to consider the matrix elements of the matrix exponential of the Laplacian (36) which we conveniently obtain from its generating function in the form

$${\left[e^{-\tau \mathcal{L}}\right]}_{0,n}=\frac{1}{n!}\frac{d^n}{d{u}^n}{e}^{-\tau (1-u)}{|}_{u=0}=e^{-\tau }\frac{{\tau }^n}{n!},\tau >0,n\in {\mathbb{N}}_0.$$

(37)

This matrix exponential is upper triangular circulant with a Poisson distribution in the non-vanishing entries. We directly verify in this representation the general properties (18)–(20) of Laplacian matrix exponentials in digraphs. This relation underlines the utmost importance of the Poisson distribution in strictly increasing walks.

These observations suggest that generating functions of good Laplacian matrix functions $g(1-{\hat{T}}_{-1})$ of strictly increasing walks on the integer line such as occurring on right-hand side of (1) are with (11) obtained as

$$g(1-u)={\int }_0^{\infty }(1-e^{-\tau (1-u)})\nu \left(\mathrm{d}\tau \right)=\sum _{k=0}^{\infty }{\left[g(1-{\hat{T}}_{-1})\right]}_{0,k}{u}^k,|u|\le 1$$

(38)

where clearly $g(1-{\hat{T}}_{-1})$ has upper triangular circulant matrix representation. For instance for the Lévy measure ${\nu }_{\alpha }\left(\mathrm{d}\tau \right)=\frac{\alpha }{\Gamma (1-\alpha )}{\tau }^{-1-\alpha }\mathrm{d}\tau$ integral (38) converges within $\alpha \in (0,1)$ and yields generating function ${(1-u)}^{\alpha }$ of the Sibuya fractional Laplacian matrix (35). The generating function of the one-step transition matrix then becomes with (38) and (24)

$$\begin{array}{ccc}{\displaystyle {\overline{W}}^{\left(g\right)}\left(u\right)}& {\displaystyle =1-\frac{1}{g\left(1\right)}g(1-u),|u|\le 1,g\left(1\right)={\left[g\left(\mathcal{L}\right)\right]}_{0,0}}\hfill & \\ \\ & {\displaystyle =\frac{1}{g\left(1\right)}{\int }_0^{\infty }{e}^{-\tau }(e^{\tau u}-1)\nu \left(\mathrm{d}\tau \right)}\hfill \end{array}$$

(39)

with the ‘generalized degree’ ${g\left(1\right)=g(1-u)|}_{u=0}={\int }_0^{\infty }[1-e^{-\tau }]\nu \left(\mathrm{d}\tau \right)>0$. The elements of the transition matrix are then straight-forwardly obtained as

$$\begin{array}{ccc}{\displaystyle W_{0,n}^{\left(g\right)}}& {\displaystyle =\frac{1}{g\left(1\right)}{\int }_0^{\infty }\frac{{\tau }^n}{n!}{e}^{-\tau }\nu \left(\mathrm{d}\tau \right)>0,}\hfill & n\in \mathbb{N}\\ \\ {\displaystyle W_{0,0}^{\left(g\right)}}& =0.\hfill \end{array}$$

(40)

The matrix elements are traced back to Poisson terms $\frac{{\tau }^n}{n!}{e}^{-\tau }$ weighted by the Lévy measure $\nu \left(\mathrm{d}\tau \right)$ resulting in strictly positive $W_{0,n}^{\left(g\right)}>0$ for $n>0$ and $W_{0,n}^{\left(g\right)}=0$ else introducing new topologies with directed edges allowing jumps of any positive integer size $Z_j\in \{1,2,\dots \}=\mathbb{N}$.

For the class of Lévy measures $\nu \left(\mathrm{d}\tau \right)=\nu \left(\tau \right)\mathrm{d}\tau$ with densities $\nu \left(\tau \right)$ which fulfill ${\int }_0^{\infty }\nu \left(\tau \right)\mathrm{d}\tau <\infty$, i.e., having existing Laplace transforms

$$\tilde{\nu }\left(s\right)={\int }_0^{\infty }{e}^{-\tau s}\nu \left(\tau \right)\mathrm{d}\tau <\infty ,\Re \left\{s\right\}\ge 0$$

(41)

we have $g(1-u)=\tilde{\nu }\left(0\right)-\tilde{\nu }(1-u)$ thus the generating function (39) of the transition matrix writes

$${\displaystyle {\overline{W}}^{\left(g\right)}\left(u\right)=\frac{\tilde{\nu }(1-u)-\tilde{\nu }\left(1\right)}{\tilde{\nu }\left(0\right)-\tilde{\nu }\left(1\right)}.}$$

(42)

Due to the restriction of convergence of (41) we notice that relation (39) is more general than (42).

After these considerations let us verify the crucial property (23) when we account for the (in our convention left-) eigenvectors of the unitary shift operator $v_n=C{e}^{i\phi n}$ ($\phi \in (-\pi ,\pi ]$) with $\mathbf{v}·\mathcal{L}={\sum }_{m=-\infty }^{\infty }{v}_m{\left[\mathcal{L}\right]}_{mn}=(1-e^{-i\phi })v_n$. Hence the continuous complex eigenvalue spectrum of the Laplacian matrix (36) is

$$\mu \left(\phi \right)=1-e^{-i\phi },\phi \in (-\pi ,\pi ]$$

(43)

with ${\mu \left(\phi \right)|}_{\phi =0}=0$ and $\Re \left\{\mu \right(\phi \left)\right\}=1-cos\phi >0$ for $\phi \ne 0$ ($\phi \in (-\pi ,\pi ]$) in accordance with (23) (here $\mathbf{W}=\mathbf{1}-\mathcal{L}=\left[{\hat{T}}_{-1}\right]$ is unitary with $\lambda \left(\phi \right)=e^{-i\phi }$ thus for all eigenvalues holds $|\lambda (\phi )|=1$). Hence (11) converges and has the canonical representation

$${\left[g\left(\mathcal{L}\right)\right]}_{n,m}=\frac{1}{2\pi }{\int }_{-\pi }^{\pi }{e}^{i\phi (m-n)}g(1-e^{-i\phi })\mathrm{d}\phi$$

(44)

where the eigenvalues of the Laplacian matrix function write with (11)

$$g(1-e^{-i\phi })={\int }_0^{\infty }[1-e^{-\tau (1-e^{-i\phi })}]\nu \left(\mathrm{d}\tau \right)$$

(45)

which is convergent since $\Re \left\{\mu \right(\phi \left)\right\}=1-cos\phi \ge 0$ ($\phi \in (-\pi ,\pi ]$). Expanding $e^{\tau e^{-i\phi }}={\sum }_{\ell =0}^{\infty }\frac{{\tau }^{\ell }}{\ell !}{e}^{-i\ell \phi }$ leads to

$$\frac{1}{2\pi }{\int }_{-\pi }^{\pi }{e}^{i\tau e^{-i\phi }}{e}^{i(m-n)\phi }\mathrm{d}\phi =\left\{\begin{array}{c}{\displaystyle \frac{{\tau }^{m-n}}{(m-n)!},m\ge n}\hfill \\ \\ {\displaystyle 0,m<n.}\hfill \end{array}\right.$$

(46)

Thus we get for the family of “good Laplacian matrix functions” as generators for strictly increasing walks the upper triangular circulant structure

$${\left[g\left(\mathcal{L}\right)\right]}_{n,m}=\left\{\begin{array}{ccc}{\displaystyle {\int }_0^{\infty }\left[{\delta }_{n,m}-e^{-\tau }\frac{{\tau }^{m-n}}{(m-n)!}\right]\nu \left(\mathrm{d}\tau \right),}& m-n\in {\mathbb{N}}_0\hfill & \\ \\ {\displaystyle 0,}& m-n<0\hfill \end{array}\right.$$

(47)

which clearly fulfills the “good Laplacian properties” (i)–(iii) and is consistent with relation (40). The convergence of these integrals is ensured by the property (12) of the Lévy measure.

Mittag-Leffler Transition Matrix

For later use we derive now in the framework of this approach the one-step transition matrix of a strictly increasing walk on the integer line with jumps following a discrete approximation of the Mittag-Leffler density. Such a distribution was analyzed in [50] and generalizations in recent articles [26,41]. To this end we consider a Lévy measure ${\nu }_{ML,\alpha }\left(\mathrm{d}\tau \right)={\nu }_{ML,\alpha }\left(\tau \right)\mathrm{d}\tau$ with Lévy density ${\nu }_{ML,\alpha }\left(\tau \right)=\lambda {\tau }^{\alpha -1}{E}_{\alpha ,\alpha }(-\lambda {\tau }^{\alpha })$ which itself is a Mittag-Leffler density and has Laplace transform $\frac{\lambda }{\lambda +s^{\alpha }}$ ($\alpha \in (0,1]$, $\lambda >0$). Then we get from the Lévy-Khintchine representation (38) the Laplacian generating function

$$g_{ML,\alpha }(1-u,\lambda )=1-\frac{\lambda }{\lambda +{(1-u)}^{\alpha }}=\frac{{(1-u)}^{\alpha }}{\lambda +{(1-u)}^{\alpha }},\lambda >0,|u|\le 1.$$

(48)

Thus with (39) the generating function of the “Mittag-Leffler transition matrix” yields

$${\overline{W}}_{ML,\alpha }(u,\lambda )=1-\frac{1}{g_{ML,\alpha }(1,\lambda )}{g}_{ML,\alpha }(1-u,\lambda )=\frac{\lambda }{\lambda +{(1-u)}^{\alpha }}(1-{(1-u)}^{\alpha }),\alpha \in (0,1]|u|\le 1.$$

(49)

This is the generating function of a discrete approximation of the Mittag-Leffler density introduced recently [41] in the context of discrete-time renewal processes and named there “Discrete Mittag-Leffler” (${\mathit{DML}}_A$) distribution (of so called “type A”, see Remark 9 in that paper). The elements of the Mittag-Leffler transition matrix are derived explicitly in relation (54).

It can be seen beforehand that (49) generates a discrete approximation of the Mittag-Leffler density by the following consideration. By introducing the scaling $\lambda \left(h\right)=h^{\alpha }{\lambda }_0$ ($h,{\lambda }_0>0$ where ${\lambda }_0$ is a new positive constant independent of h) and with $u=e^{-hs}$ shows that the generating function (49) converges to the Laplace transform of a Mittag-Leffler density, namely

$$\tilde{\mathcal{W}}({\lambda }_0,s)=\underset{h\to 0}{lim}{\overline{W}}_{ML,\alpha }(e^{-hs},{\lambda }_0{h}^{\alpha })=\frac{{\lambda }_0}{{\lambda }_0+s^{\alpha }}$$

(50)

where with $\tilde{f}\left(s\right)$ we denote the (spatial) Laplace transform of $f\left(x\right)$. It is instructive to recall this limit in terms of shift operator representation (See Appendix A and [26] for details)

$$\begin{array}{c}{\displaystyle {\mathcal{W}}_{ML,\alpha }(x,{\lambda }_0)=\underset{h\to 0}{lim}{\overline{W}}_{ML,\alpha }({\hat{T}}_{-h},{\lambda }_0{h}^{\alpha }){\delta }_h\left(x\right)=\underset{h\to 0}{lim}\frac{{\lambda }_0{h}^{\alpha }}{{\lambda }_0{h}^{\alpha }+{(1-{\hat{T}}_{-h})}^{\alpha }}{\delta }_h\left(x\right)}\hfill \\ \\ {\displaystyle =\underset{h\to 0}{lim}\frac{{\lambda }_0}{{\lambda }_0+h^{-\alpha }{(1-{\hat{T}}_{-h})}^{\alpha }}{\delta }_h\left(x\right)=\frac{{\lambda }_0}{{\lambda }_0+D_x^{\alpha }}\delta \left(x\right)=\underset{h\to 0}{lim}\frac{1}{h}{\left[{\overline{W}}_{ML,\alpha }\left({\lambda }_0{h}^{\alpha }\right)\right]}_{0,\frac{x}{h}}}\hfill \end{array}$$

(51)

which has Laplace transform converging to (50) and contains the discrete δ-distribution ${\delta }_h\left(x\right)$ defined in (A3). Relation (51) defines the “well-scaled” continuous-space limit density kernel (‘transition density kernel’) and converges to the Mittag-Leffler density (60). In the second line of (51) appears the (Riemann-Liouville) fractional derivative operator $D_x^{\alpha }$ of order $\alpha$ as ${lim}_{h\to 0}{h}^{-\alpha }{(1-e^{-h{D}_x})}^{\alpha }=D_x^{\alpha }$ (with ${\hat{T}}_{-h}=e^{-h{D}_x}$).

Now we derive the matrix elements of the Mittag-Leffler transition matrix in explicit form. For our convenience and later use we introduce the Pochhammer-symbol

$${\left(c\right)}_m=\frac{\Gamma (c+m)}{\Gamma \left(c\right)}=\left\{\begin{array}{c}1,m=0,\hfill \\ \\ c(c+1)\dots (c+m-1),m=1,2,\dots \hfill \end{array}\right.$$

(52)

where especially ${\left(c\right)}_0=1$ and ${\left(0\right)}_m={\delta }_{m0}$. Then we have with (49) the expansions

$$\begin{array}{cc}{\displaystyle \frac{\lambda }{\lambda +{(1-u)}^{\alpha }}}& {\displaystyle =\left\{\begin{array}{ccc}{\displaystyle \sum _{m=0}^{\infty }{(-1)}^m{\lambda }^{m+1}{(1-u)}^{-(m+1)\alpha },}& {\displaystyle {|1-u|}^{-1}{\lambda }^{\frac{1}{\alpha }}<1}\hfill & \\ \\ {\displaystyle \sum _{m=0}^{\infty }{(-1)}^m{\lambda }^{-m}{(1-u)}^{\alpha m},}& {\displaystyle {|1-u|}^{-1}{\lambda }^{\frac{1}{\alpha }}>1}\hfill \end{array}\right.}\hfill \\ \\ {\displaystyle {\overline{W}}_{ML,\alpha }(u,\lambda )}& {\displaystyle =\left\{\begin{array}{ccc}{\displaystyle \sum _{m=0}^{\infty }{(-1)}^m{\lambda }^{m+1}\left[{(1-u)}^{-(m+1)\alpha }-{(1-u)}^{-m\alpha }\right],}& {\displaystyle {|1-u|}^{-1}{\lambda }^{\frac{1}{\alpha }}<1}\hfill & \\ \\ {\displaystyle \sum _{m=0}^{\infty }{(-1)}^m{\lambda }^{-m}\left[{(1-u)}^{\alpha m}-{(1-u)}^{\alpha (m+1)}\right],}& {\displaystyle {|1-u|}^{-1}{\lambda }^{\frac{1}{\alpha }}>1}\hfill \end{array}\right.}\hfill \end{array}$$

(53)

In order to determine the elements of the Mittag-Leffler transition matrix (since all matrices commute here, we adopt the notation $\frac{\mathbf{A}}{\mathbf{B}}=\mathbf{A}·{\mathbf{B}}^{-1}$) we have to account for the cases of convergence in (53) at $u=0$ to arrive at

$$\begin{array}{ccc}{\displaystyle {\left[W_{ML,\alpha }\left(\lambda \right)\right]}_{0,n}}& {\displaystyle ={\left[\frac{\lambda \mathbf{1}}{\lambda \mathbf{1}+{\mathcal{L}}^{\alpha }}(\mathbf{1}-{\mathcal{L}}^{\alpha })\right]}_{0,n}=\frac{1}{n!}\frac{d^n}{d{u}^u}{\overline{W}}_{ML,\alpha }{(u,\lambda )|}_{u=0},n\in {\mathbb{N}}_0}\hfill & \\ \\ & {\displaystyle =\left\{\begin{array}{ccc}{\displaystyle \frac{\lambda }{n!}\sum _{m=0}^{\infty }{(-1)}^m{\lambda }^m\left[\frac{\Gamma \left(\alpha \right(m+1)+n)}{\Gamma \left(\alpha \right(m+1\left)\right)}-\frac{\Gamma (\alpha m+n)}{\Gamma \left(\alpha m\right)}\right],}& {\displaystyle 0<\lambda <1}\hfill & \\ \\ {\displaystyle \frac{{(-1)}^n}{n!}\sum _{m=0}^{\infty }{(-1)}^m{\lambda }^{-m}\left[\frac{\Gamma (\alpha m+1)}{\Gamma (\alpha m-n+1)}-\frac{\Gamma \left(\alpha \right(m+1)+1)}{\Gamma \left(\alpha \right(m+1)-n+1)}\right],}& {\displaystyle \lambda >1}\hfill \end{array}\right.}\hfill \end{array}$$

(54)

where these series converge absolutely by accounting for the asymptotic behavior of the terms for m large as $\frac{\Gamma \left(\alpha \right(m+1)+n)}{\Gamma \left(\alpha \right(m+1\left)\right)}\sim m^n{\lambda }^m$ for $\lambda <1$, and in the same way $\sim m^n{\lambda }^{-m}$ for $\lambda >1$, respectively. We also have ${\left[W_{ML,\alpha }\left(\lambda \right)\right]}_{0,0}=0$.

It is worthy also to consider the case $\alpha =1$ where (49) takes the form

$$\begin{array}{ccc}{\displaystyle {\overline{W}}_{ML,1}(u,\lambda )}& {\displaystyle =\frac{u\lambda }{\lambda +1-u}=\frac{\lambda }{(\lambda +1)}\frac{u}{(1-\frac{u}{\lambda +1})}}\hfill & \\ \\ & {\displaystyle =\frac{\lambda }{(\lambda +1)}\sum _{k=0}^{\infty }\frac{u^{k+1}}{{(\lambda +1)}^k}=\sum _{k=1}^{\infty }p{q}^{k-1}{u}^k}\hfill & \\ \\ & {\displaystyle p=\frac{\lambda }{\lambda +1},q=\frac{1}{\lambda +1},}\hfill & p+q=1\end{array}$$

(55)

and hence we obtain for $\alpha =1$ the transition matrix

$$\begin{array}{ccc}{\displaystyle {\left[W_{ML,1}\left(\lambda \right)\right]}_{0,n}}& =p{q}^{n-1},n\in \mathbb{N}\hfill & \\ \\ {\displaystyle {\left[W_{ML,1}\left(\lambda \right)\right]}_{0,0}}& =0\hfill \end{array}$$

(56)

which recovers the geometrical distribution (See also [41]). By using the explicit formula (54) for the Mittag-Leffler transition matrix, it is now not a big deal to perform explicitly the “well-scaled” continuous-space limit (51). Accounting for the asymptotic relation $\frac{{\left(\beta \right)}_n}{n!}\sim \frac{n^{\beta -1}}{\Gamma \left(\beta \right)}$ ($n=\frac{x}{h}\to \infty$), and by introducing the scaling $\lambda \left(h\right)={\lambda }_0{h}^{\alpha }\to 0$ which is covered in relation (54) by the case $0<\lambda <1$, we arrive at

$$\begin{array}{ccc}{\displaystyle {\mathcal{W}}_{ML,\alpha }(x,{\lambda }_0)}& {\displaystyle =\underset{h\to 0}{lim}\frac{1}{h}{\left[{\overline{W}}_{ML,\alpha }\left({\lambda }_0{h}^{\alpha }\right)\right]}_{0,\frac{x}{h}},x\in h{\mathbb{N}}_0}\hfill & \\ \\ & {\displaystyle =\underset{h\to 0}{lim}\sum _{m=0}^{\infty }{(-1)}^m{\lambda }_0^{m+1}\frac{1}{\frac{x}{h}!}{h}^{\alpha (m+1)-1}\left[{\left((m+1)\alpha \right)}_{\frac{x}{h}}-{\left(m\alpha \right)}_{\frac{x}{h}}\right].}\hfill \end{array}$$

(57)

We employ here the Pochhammer symbol (52) and we account for

$$\underset{h\to 0}{lim}{h}^{\alpha (m+1)-1}{\left((m+1)\alpha \right)}_{\frac{x}{h}}=\underset{h\to 0}{lim}{h}^{\alpha (m+1)-1}\frac{1}{\Gamma \left[\alpha \right(m+1\left)\right]}{\left(\frac{x}{h}\right)}^{\alpha (m+1)-1}=\frac{x^{\alpha (m+1)-1}}{\Gamma \left[\alpha \right(m+1\left)\right]}$$

(58)

with $x\in {\mathbb{R}}^{+}$. The second term in (57) has a vanishing limit

$$\underset{h\to 0}{lim}{h}^{\alpha (m+1)-1}{\left(m\alpha \right)}_{\frac{x}{h}}=\underset{h\to 0}{lim}{h}^{\alpha (m+1)-1}\frac{1}{\Gamma \left(\alpha m\right)}{\left(\frac{x}{h}\right)}^{\alpha m-1}=\underset{h\to 0}{lim}{h}^{\alpha }\frac{x^{\alpha m-1}}{\Gamma \left(\alpha m\right)}=0.$$

(59)

We hence obtain for the well-scaled continuous-space limit (57) as anticipated the Mittag-Leffler density

$$\begin{array}{ccc}{\displaystyle {\mathcal{W}}_{ML,\alpha }(x,{\lambda }_0)}& {\displaystyle ={\lambda }_0{x}^{\alpha -1}\sum _{m=0}^{\infty }\frac{{(-{\lambda }_0{x}^{\alpha })}^m}{\Gamma (\alpha m+\alpha )}}\hfill & \\ \\ & {\displaystyle ={\lambda }_0{x}^{\alpha -1}{E}_{\alpha ,\alpha }(-{\lambda }_0{x}^{\alpha })=\frac{d}{dx}[1-E_{\alpha }(-{\lambda }_0{x}^{\alpha })],}\hfill & x\in {\mathbb{R}}^{+},\alpha \in (0,1]\end{array}$$

(60)

containing the generalized Mittag-Leffler function $E_{\alpha ,\gamma }\left(z\right)={\sum }_{m=0}^{\infty }\frac{z^m}{\Gamma (\alpha m+\gamma )}$. For $\alpha =1$ (60) recovers the exponential density ${\mathcal{W}}_{ML,1}(x,{\lambda }_0)={\lambda }_0{e}^{-{\lambda }_{0}x}$ which also is the continuous-space limit of (56) obtained with $p\left(h\right)=\frac{{\lambda }_{0}h}{{\lambda }_{0}h+1}$ and $q\left(h\right)=\frac{1}{{\lambda }_{0}h+1}$ thus

$$\begin{array}{ccc}{\displaystyle {\mathcal{W}}_{ML,1}(x,{\lambda }_0)}& {\displaystyle =\underset{h\to 0}{lim}\frac{1}{h}{\left[{\overline{W}}_{ML,1}\left({\lambda }_{0}h\right)\right]}_{0,\frac{x}{h}}}\hfill & \\ \\ & {\displaystyle =\underset{h\to 0}{lim}{\lambda }_0{(1+{\lambda }_{0}h)}^{-\frac{x}{h}}}\hfill & \\ \\ & {\displaystyle ={\lambda }_0{e}^{-{\lambda }_{0}x}=\frac{d}{dx}(1-e^{-{\lambda }_{0}x}),}\hfill & x\in {\mathbb{R}}^{+}.\end{array}$$

(61)

5. Space-Time Generalizations of the Poisson Process

5.1. The Classical Cases

During the last two decades, an increasing interest in generalizations of the Poisson renewal process has emerged. The most natural generalization probably is obtained when the exponential waiting time density is generalized by a Mittag-Leffler density. The resulting generalization is the fractional Poisson process which was introduced and analyzed by several authors [40,51,52,53,54]. Then space-time generalizations of the Poisson process were introduced such as the ‘space-fractional Poisson Process’, the ‘space-time fractional Poisson process’ [42] and further generalizations of the latter [45] were developed within the last decade. Generally these space-time generalizations can be seen as strictly increasing walks on the integer line time-changed with an independent renewal process. Before we introduce in Section 6 such a generalization, let us briefly recall these classical cases.

5.2. Poisson Process

We consider the Cauchy problem

$$\begin{array}{ccc}{\displaystyle \frac{d}{dt}{p}_n\left(t\right)}& {\displaystyle =-\xi p_n\left(t\right)+\xi p_{n-1}\left(t\right)=-\xi (1-{\hat{T}}_{-1})p_n\left(t\right),}\hfill & {\displaystyle p_n\left(0\right)={\delta }_{n0},\xi >0,t\ge 0}\\ \\ {\displaystyle \frac{d}{dt}\mathbf{p}\left(t\right)}& {\displaystyle =-\xi \mathbf{p}\left(t\right)·\mathcal{L}.}\hfill \end{array}$$

(62)

Indeed, the Cauchy problem (62) defines the time-evolution of the state probabilities in the standard Poisson counting process which is the simplest variant of the general class (1) (where $\mathcal{L}$ is the circulant Laplacian matrix (36)). Equation (62) is the Kolmogorov-Feller (also referred to as Kolmogorov-forward) equation which is solved by the state probabilities of the standard Poisson counting process $p_n\left(t\right)=\mathbb{P}(\mathcal{N}\left(t\right)=n)$ ($n\in {\mathbb{N}}_0$, $t\ge 0$) where $\mathcal{N}\left(t\right)\in {\mathbb{N}}_0$ counts the events within the time interval $[0,t]$. We can conceive the Poisson counting process in the Montroll-Weiss CTRW picture [11] as a strictly increasing random walk with one-step transition matrix ${\left[\mathbf{W}\right]}_{m,n}={[{\hat{T}}_{-1}]}_{m,n}={\delta }_{m,n-1}$ subordinated to a Poisson process where at each arrival the walker almost surely makes a jump of increment $Z_j=1$. Clearly the solution of the Cauchy problem (62) is the well-known standard Poisson distribution

$$p_n\left(t\right)={[\mathbf{p}\left(0\right)·e^{-\xi t\mathcal{L}}]}_n=\frac{1}{n!}\frac{d^n}{d{u}^n}{e}^{-\xi t(1-u)}{|}_{u=0}=\frac{{\left(\xi t\right)}^n}{n!}{e}^{-\xi t},n\in {\mathbb{N}}_0,t\in {\mathbb{R}}^{+}.$$

(63)

5.3. Fractional Poisson Process

The most natural time-generalization of the Poisson process is defined by the Cauchy problem

$$\begin{array}{ccc}{\displaystyle \frac{d^{\beta }}{d{t}^{\beta }}{p}_n\left(t\right)}& {\displaystyle =-\xi (1-{\hat{T}}_{-1})p_n\left(t\right),\xi >0,p_n\left(0\right)={\delta }_{n0},\beta \in (0,1]}\hfill & \\ \\ {\displaystyle \frac{d^{\beta }}{d{t}^{\beta }}\mathbf{p}\left(t\right)}& {\displaystyle =-\xi \mathbf{p}\left(t\right)·\mathcal{L}}\hfill \end{array}$$

(64)

which is the well-known fractional Kolmogorov-Feller equation of the fractional Poisson process [52] where $\frac{d^{\beta }}{d{t}^{\beta }}$ denotes the Caputo fractional derivative of order $\beta$ defined as, e.g., [55]

$$\frac{d^{\beta }}{d{t}^{\beta }}p\left(t\right)=\frac{1}{\Gamma (1-\beta )}{\int }_0^t{(t-\tau )}^{-\beta }\frac{d}{d\tau }p\left(\tau \right)\mathrm{d}\tau ,\beta \in (0,1]$$

(65)

and in the limit ${lim}_{\beta \to 1-}\frac{{(t-\tau )}^{-\beta }}{\Gamma (1-\beta )}=\delta (t-\tau )$ (65) recovers the first-order derivative $\frac{d}{dt}p\left(t\right)$; thus, the fractional Poisson process turns into the standard Poisson process. The solution of (64) writes

$$\mathbf{p}\left(t\right)=\mathbf{p}\left(0\right)·E_{\beta }(-\xi t^{\beta }\mathcal{L})=[\dots ,{\left[E_{\beta }(-\xi t^{\beta }\mathcal{L})\right]}_{0,n},\dots ]$$

(66)

where the Mittag-Leffler matrix function comes into play

$$E_{\beta }(-\xi t^{\beta }\mathcal{L})=\sum _{m=0}^{\infty }\frac{{(-\xi t^{\beta })}^m}{\Gamma (\beta m+1)}{\mathcal{L}}^m$$

(67)

with the Mittag-Leffler function $E_{\beta }\left(z\right)={\sum }_{m=0}^{\infty }\frac{z^m}{\Gamma (\beta m+1)}$. Taking into account the generating function ${\sum }_{k=0}^{\infty }{\left[{\mathcal{L}}^m\right]}_{0,k}{u}^k={(1-u)}^m$, we obtain for the generating function of the Mittag-Leffler matrix function (67) the expression ${\sum }_{k=0}^{\infty }{\left[E_{\beta }(-\xi t^{\beta }\mathcal{L})\right]}_{0,k}{u}^k=E_{\beta }(-\xi t^{\beta }(1-u))$. Hence with (31) we get for the components of the state vector (66)

$$\begin{array}{ccc}{\displaystyle p_n^{\beta }\left(t\right)}& {\displaystyle ={\left[E_{\beta }(-\xi t^{\beta }\mathcal{L})\right]}_{0,n}=\frac{1}{n!}\frac{d^n}{d{u}^n}{E}_{\beta }(-(1-u)\xi t^{\beta }){|}_{u=0},n\in {\mathbb{N}}_0,t\in {\mathbb{R}}^{+}}\hfill & \\ \\ {\displaystyle }& {\displaystyle =\frac{{\left(\xi t^{\beta }\right)}^n}{n!}\sum _{m=0}^{\infty }\frac{(m+n)!}{m!}\frac{{(-\xi t^{\beta })}^m}{\Gamma \left(\right(m+n)\beta +1)}.}\hfill \end{array}$$

(68)

We identify (68) indeed with the state probabilities of Laskin’s fractional Poisson distribution [52] which also was derived in different manners by several further authors [51,53,54]. For $\beta =1$ (68) recovers the Poisson distribution (63).

5.4. Space-Time Fractional Poisson Process

Orsingher and Polito have generalized the Cauchy problems (62), (64) by introducing a spatial generalization $\mathcal{L}\to {\mathcal{L}}^{\alpha }$ where the right-hand sides of (62) and (64), respectively, are generalized by the (Sibuya-) fractional Laplacian matrix (35). They named these processes space-fractional Poisson process and space-time fractional Poisson process, respectively [42]. Indeed the space-fractional Poisson process is a Sibuya-walk subordinated to an independent Poisson process, and the space-time fractional Poisson process is a Sibuya walk time-changed with an independent fractional Poisson process. The Cauchy problem defining the space-time fractional Poisson process then writes

$$\begin{array}{ccc}{\displaystyle \frac{d^{\beta }}{d{t}^{\beta }}{p}_n\left(t\right)}& {\displaystyle =-\xi {(1-{\hat{T}}_{-1})}^{\alpha }{p}_n\left(t\right),\alpha ,\beta \in (0,1],\xi >0,t\ge 0}\hfill & \\ \\ {\displaystyle p_n{\left(t\right)|}_{t=0}}& {\displaystyle ={\delta }_{n0}}\hfill & \\ \\ {\displaystyle \frac{d^{\beta }}{d{t}^{\beta }}\mathbf{p}\left(t\right)}& =-\xi \mathbf{p}\left(t\right)·{\mathcal{L}}^{\alpha }\hfill \end{array}$$

(69)

For $\beta =1$ the space-fractional Poisson process is recovered and for $\alpha =1$ the fractional Poisson process. We obtain the state vector solving (69) as

$${\mathbf{p}}^{\alpha ,\beta }\left(t\right)=\mathbf{p}\left(0\right)·E_{\beta }(-\xi t^{\beta }{\mathcal{L}}^{\alpha }),\alpha ,\beta \in (0,1],t\ge 0$$

(70)

having generating function

$${\overline{p}}^{\alpha ,\beta }(u,t)=E_{\beta }(-\xi t^{\beta }{(1-u)}^{\alpha })$$

(71)

where ${\overline{p}}^{\alpha ,\beta }{(u,t)|}_{u=1}=1$ reflects row-stochasticity of the upper triangular circulant transition matrix $E_{\beta }(-\xi t^{\beta }{\mathcal{L}}^{\alpha })$ (normalization of the state probabilities). The state distribution is obtained as

$$\begin{array}{c}{\displaystyle p_n^{\alpha ,\beta }\left(t\right)=\frac{1}{n!}\frac{d^n}{d{u}^n}{E}_{\beta }(-{(1-u)}^{\alpha }\xi t^{\beta }){|}_{u=0},n\in {\mathbb{N}}_0,t\ge 0}\hfill \\ \\ {\displaystyle =\sum _{m=0}^{\infty }\frac{{(-\xi t^{\beta })}^m}{\Gamma (\beta m+1)}\frac{1}{n!}\frac{d^n}{d{u}^n}{(1-u)}^{\alpha m}{|}_{u=0}=\frac{{(-1)}^n}{n!}\sum _{m=0}^{\infty }\frac{\Gamma (\alpha m+1)}{\Gamma (\alpha m+1-n)}\frac{{(-\xi t^{\beta })}^m}{\Gamma (\beta m+1)}}\hfill \end{array}$$

(72)

which is the expression obtained by Orsingher and Polito (Equation (2.28) in [42]) and recovers for $\beta =1$ their expression for the space-fractional Poisson process derived in the same paper [42]. Further generalizations of the space-fractional Poisson are considered in the references [56,57].

5.5. Well-Scaled Diffusion Limits

We consider briefly the diffusion limit of standard Poisson $\alpha ,\beta =1$. To this end we define a well-scaled continuous-space limit in (63) by introducing the scaling $\xi \left(h\right)={\xi }_0{h}^{-1}$ (where ${\xi }_0>0$ is a new dimensional constant independent of h) to arrive at

$${\mathcal{P}}_{1,1}(x,t)=\underset{h\to 0}{lim}{e}^{-{\xi }_{0}t{h}^{-1}(1-e^{-h{D}_x})}{\delta }_h\left(x\right)=e^{-{\xi }_{0}t{D}_x}\delta \left(x\right)={\hat{T}}_{-{\xi }_{0}t}\delta \left(x\right)=\delta (x-{\xi }_{0}t)$$

(73)

which is a moving Dirac $\delta$-distribution propagating with constant velocity ${\xi }_0$ in the positive x-direction. It is immediately checked that (73) solves the continuous-space limit of (62) which is the Cauchy problem

$$\frac{\partial }{\partial t}{\mathcal{P}}_{1,1}(x,t)=-{\xi }_0\frac{\partial }{\partial x}{\mathcal{P}}_{1,1}(x,t),{\mathcal{P}}_{1,1}{(x,t)|}_{t=0}=\delta \left(x\right).$$

(74)

5.6. Diffusion Limit of Space-Time Fractional Poisson

Further consider the space-time fractional Poisson process $\alpha \in (0,1)$ and $\beta \in (0,1]$ with the scaling $\xi \left(h\right)={\xi }_0{h}^{-\alpha }$ (${\xi }_0>0$). Then we can write the continuous-space diffusion equation which emerges from the well-scaled limit of (69). By accounting for the limit ${lim}_{h\to 0}{h}^{-\alpha }{(1-{\hat{T}}_{-h})}^{\alpha }=D_x^{\alpha }$ which takes the Riemann-Liouville fractional derivative of order $\alpha$ (e.g., [55,58,59] and many others) we obtain for the well-scaled limit of (69) the space-time fractional diffusion equation

$$\frac{{\partial }^{\beta }}{\partial t^{\beta }}\mathcal{P}(x,t)=-{\xi }_0{D}_x^{\alpha }{\mathcal{P}(x,t),\mathcal{P}(x,t)|}_{t=0}=\delta \left(x\right).$$

(75)

On the right hand side of (75) occurs the spatial Riemann-Liouville fractional derivative of order $\alpha \in (0,1)$ defined as, e.g., [55,58]

$$D_x^{\alpha }\mathcal{P}(x,t)=\frac{\partial }{\partial x}{\int }_0^x\frac{{(x-\tau )}^{-\alpha }}{\Gamma (1-\alpha )}\mathcal{P}(\tau ,t)\mathrm{d}\tau ,\alpha \in (0,1).$$

(76)

We notice that the fractional Laplacian ${\mathcal{L}}^{\alpha }=\left[{(1-{\hat{T}}_{-1})}^{\alpha }\right]$ on the right-hand side of (69) does not contain its own scaling parameter. Therefore, in order to get an existing limit we have to rescale the constant $\xi \left(h\right)={\xi }_0{h}^{-\alpha }$ (having dimension ${[sec]}^{-\beta }$) where ${\xi }_0$ has units $\left[{sec}^{-\beta }{\mathrm{cm}}^{\alpha }\right]$.

6. Space-Time Mittag-Leffler Process

Having recalled these classical cases, we introduce here a generalization of the Poisson process which is a strictly increasing walk on the integer line with Mittag-Leffler jumps taking place at independent fractional Poisson arrival times. We also highlight the connections with the Montroll-Weiss CTRW picture in more details. We call this process ‘space-time Mittag-Leffler process’ which we define by a Cauchy problem of the general type (1), namely

$$\begin{array}{ccc}{\displaystyle \frac{d^{\beta }}{d{t}^{\beta }}{p}_n\left(t\right)}\hfill & =& {\displaystyle -\xi (\lambda +1)\frac{{(1-{\hat{T}}_{-1})}^{\alpha }}{(\lambda +{(1-{\hat{T}}_{-1})}^{\alpha })}{p}_n\left(t\right),\xi ,\lambda >0,t\ge 0,\alpha ,\beta \in (0,1],n\in {\mathbb{N}}_0}\hfill \\ \\ {\displaystyle p_n{\left(t\right)|}_{t=0}}\hfill & =& {\displaystyle {\delta }_{0,n}}\hfill \\ \\ {\displaystyle \frac{d^{\beta }}{d{t}^{\beta }}\mathbf{p}\left(t\right)}\hfill & =& {\displaystyle -\frac{\xi }{{\left[g_{ML,\alpha }\right]}_{0,0}}\mathbf{p}\left(t\right)·g_{ML,\alpha }\left(\mathcal{L}\right),{[\mathbf{p}\left(t\right)·g_{ML,\alpha }\left(\mathcal{L}\right)]}_n=\sum _{m=0}^n{p}_m\left(t\right){\left[g_{ML,\alpha }\left(\mathcal{L}\right)\right]}_{0,n-m}.}\hfill \end{array}$$

(77)

The right-hand side contains the good Laplacian matrix function $g_{ML,\alpha }\left(\mathcal{L}\right)=\frac{{\mathcal{L}}^{\alpha }}{\lambda +{\mathcal{L}}^{\alpha }}$ of (48) (with generalized degree ${\left[g_{ML,\alpha }\right]}_{0,0}=g_{ML,\alpha }(1-u,\lambda ){|}_{u=0}=\frac{1}{\lambda +1}$) and generates discrete Mittag-Leffler jumps with transition matrix (54). The last line indicates the matrix representation and uses the upper triangular circulant property of $g_{ML,\alpha }\left(\mathcal{L}\right)$ (see (30)) and $\frac{d^{\beta }}{d{t}^{\beta }}$ stands for the Caputo fractional derivative (65). The generating function representation of Cauchy problem (77) then writes

$$\begin{array}{ccc}{\displaystyle \frac{d^{\beta }}{d{t}^{\beta }}\overline{p}(u,t)}\hfill & =& {\displaystyle -\xi \frac{(\lambda +1){(1-u)}^{\alpha }}{\lambda +{(1-u)}^{\alpha }}\overline{p}(u,t),|u|\le 1,}\hfill \\ \\ {\displaystyle \overline{p}{(u,t)|}_{t=0}}\hfill & =& 1.\hfill \end{array}$$

(78)

For the solution of (77) we can write

$${\mathbf{p}}_{\lambda ,\xi }^{\alpha ,\beta }\left(t\right)=\mathbf{p}\left(0\right)·E_{\beta }\left(-\xi t^{\beta }(\lambda +1)\frac{{\mathcal{L}}^{\alpha }}{\lambda \mathbf{1}+{\mathcal{L}}^{\alpha }}\right)$$

(79)

with the generating function of the state-probabilities

$$\begin{array}{ccc}{\displaystyle {\overline{p}}_{\lambda ,\xi }^{\alpha ,\beta }(u,t)}\hfill & =& {\displaystyle E_{\beta }\left(\frac{-\xi t^{\beta }(\lambda +1){(1-u)}^{\alpha }}{\lambda +{(1-u)}^{\alpha }}\right)}\hfill \\ \\ & =& {\displaystyle \sum _{m=0}^{\infty }\frac{{[-\xi t^{\beta }(\lambda +1)]}^m}{\Gamma (\beta m+1)}\frac{{(1-u)}^{\alpha m}}{{[\lambda +{(1-u)}^{\alpha }]}^m}}\hfill \end{array}$$

(80)

involving the standard Mittag-Leffler function $E_{\beta }\left(z\right)$. In the Poisson limit $\beta =1$ due to $E_1\left(z\right)=e^z$ the Mittag-Leffler functions in all relations recover exponentials and the Caputo derivative recovers the standard first order time-derivative. We observe that ${\overline{p}}_{\lambda ,\xi }^{\alpha ,\beta }{(u,t)|}_{u=1}=1$ reflecting normalization of the state probabilities. Expression (80) contains

$$\begin{array}{ccc}{\displaystyle {\overline{\mathcal{E}}}_{\alpha }^{\left(m\right)}(\lambda ,u)}\hfill & =& {\displaystyle {\left[g_{ML,\alpha }(1-u)\right]}^m=\frac{1}{{[1+\lambda {(1-u)}^{-\alpha }]}^m}}\hfill \\ \\ {\displaystyle }\hfill & =& {\displaystyle \sum _{s=0}^{\infty }\frac{{\left(m\right)}_s}{s!}{(-\lambda )}^s{(1-u)}^{-\alpha s},{\lambda }^{\frac{1}{\alpha }}{|1-u|}^{-1}<1}\hfill \end{array}$$

(81)

and plainly

$$\begin{array}{ccc}{\displaystyle {\overline{\mathcal{E}}}_{\alpha }^{\left(m\right)}(\lambda ,u)}\hfill & =& {\displaystyle {\left[g_{ML,\alpha }(1-u)\right]}^m=\frac{{\lambda }^{-m}{(1-u)}^{\alpha m}}{{(1+{\lambda }^{-1}{(1-u)}^{\alpha })}^m}}\hfill \\ \\ & =& {\displaystyle \sum _{s=0}^{\infty }\frac{{\left(m\right)}_s}{s!}{(-1)}^s{\lambda }^{-s-m}{(1-u)}^{\alpha (s+m)},{\lambda }^{\frac{1}{\alpha }}{|1-u|}^{-1}>1}\hfill \end{array}$$

(82)

where ${\left(m\right)}_s=m(m+1)..(m+s-1)$ stands for the Pochhammer-symbol (52). ${\overline{\mathcal{E}}}_{\alpha }^{\left(m\right)}(\lambda ,u)$ is the generating function of a discrete version of a so called Prabhakar kernel [21,26,60]. We will show this connection explicitly in subsequent analysis of well-scaled continuous-space limits. The state-probabilities solving (77) are then with (80) and (81), (82) obtained as

$$\begin{array}{ccc}{\displaystyle p_{\lambda ,\xi ,n}^{\alpha ,\beta }\left(t\right)}\hfill & =& {\displaystyle \frac{1}{n!}\frac{d^n}{d{u}^n}{E}_{\beta }\left(\frac{-(\lambda +1){(1-u)}^{\alpha }\xi t^{\beta }}{\lambda +{(1-u)}^{\alpha }}\right){|}_{u=0},n\in {\mathbb{N}}_0}\hfill \\ \\ & =& {\displaystyle \sum _{m=0}^{\infty }\frac{{[-(\lambda +1)\xi t^{\beta }]}^m}{\Gamma (\beta m+1)}{\mathcal{E}}_{\alpha }^{\left(m\right)}(\lambda ,n).}\hfill \end{array}$$

(83)

Since we take the derivatives at $u=0$ we need to account for that (81) converges at $u=0$ for $0<\lambda <1$, whereas (82) converges at $u=0$ for $\lambda >1$. We hence arrive at

$${\mathcal{E}}_{\alpha }^{\left(m\right)}(\lambda ,n)=\left\{\begin{array}{cc}{\displaystyle \frac{1}{n!}\sum _{s=0}^{\infty }\frac{{(-\lambda )}^s{\left(m\right)}_s}{s!}\frac{\Gamma (\alpha s+n)}{\Gamma \left(\alpha s\right)},}& 0<\lambda <1\hfill \\ \\ {\displaystyle \frac{{(-1)}^n{\lambda }^{-m}}{n!}\sum _{s=0}^{\infty }\frac{{(-1)}^s{\lambda }^{-s}{\left(m\right)}_s}{s!}\frac{\Gamma \left(\alpha \right(s+m)+1)}{\Gamma \left(\alpha \right(s+m)-n+1)},}& \lambda >1\hfill \end{array}\right.n\in {\mathbb{N}}_0.$$

(84)

One can easily verify in view of the asymptotic scaling $\frac{\Gamma (s\alpha +a)}{\Gamma (\alpha s+b)}\sim s^{a-b}$ and $\frac{{\left(m\right)}_s}{s!}\sim s^{m-1}$ for $s\to \infty$ that the series (84) for the two cases converge absolutely. We directly verify the initial condition $p_{\lambda ,\xi ,n}^{\alpha ,\beta }{\left(t\right)|}_{t=0}={\mathcal{E}}_{\alpha }^{\left(0\right)}(\lambda ,n)={\delta }_{0,n}$. For $n=0$ we have ${\mathcal{E}}_{\alpha }^{\left(m\right)}(\lambda ,0)={(1+\lambda )}^{-m}$ thus the ‘survival probability’ (probability that the walker at time t still is in the initial state $n=0$) is Mittag-Leffler, namely

$$p_{\lambda ,\xi ,0}^{\alpha ,\beta }\left(t\right)={\overline{p}}_{\lambda ,\xi }^{\alpha ,\beta }{(u,t)|}_{u=0}=E_{\beta }(-\xi t^{\beta })$$

(85)

which is necessarily the survival probability in the fractional Poisson process. We notice that for $\beta \in (0,1)$ the space-time Mittag-Leffler process is non-Markovian with long memory features, whereas for $\beta =1$ it becomes Markovian due to the memoryless nature of the standard Poisson process, e.g., [7,52].

In Figure 1 the time-dependence of the state probabilities (83) for $n=0,1,2$ is depicted for different values of $\alpha ,\beta \in (0,1]$, respectively. The survival probability $n=0$ is independent of $\alpha$ and given by the Mittag-Leffler survival probability (85) which turns for $\beta =1$ into an exponential (see the plots on the left) with initial conditions $p_{\lambda ,\xi ,0}^{\alpha ,\beta }{\left(t\right)|}_{t=0}=1$ whereas for $n=1,2$ the zero initial conditions can be seen in the plots. For large dimensionless times $\tau =t{\xi }^{\frac{1}{\beta }}\to \infty$ we have for $\beta \in (0,1)$ a power-law decay $p_{\lambda ,\xi ,n}^{\alpha ,\beta }\left(t\right)\sim {\tau }^{-\beta }$ (see asymptotic relation (90)).

6.1. Asymptotic Behavior

It is worthwhile to consider the asymptotic behavior for n large and finite (dimensionless) times ${\xi }^{\frac{1}{\beta }}t$. To this end consider generating function for $u\to 1-0$ for $\alpha \in (0,1)$, namely

$$E_{\beta }\left(\frac{-(\lambda +1)\xi t^{\beta }{(1-u)}^{\alpha }}{\lambda (1+{\lambda }^{-1}{(1-u)}^{\alpha })}\right)=1-\frac{(\lambda +1)\xi t^{\beta }}{\Gamma (\beta +1)\lambda }{(1-u)}^{\alpha }+\mathcal{O}{(1-u)}^{\alpha }$$

(86)

where $\mathcal{O}{(1-u)}^{\alpha }$ indicates higher orders such that ${lim}_{u\to 1-}{(1-u)}^{-\alpha }\mathcal{O}{(1-u)}^{\alpha }=0$. We notice in view of (86) that this asymptotic behavior is for $\alpha \in (0,1)$ of the same fat-tailed type as for the state probabilities (72) of the space-time fractional Poisson process and also as Sibuya($\alpha$), namely

$$\begin{array}{c}{\displaystyle p_{\lambda ,\xi ,n}^{\alpha ,\beta }\left(t\right)\sim -\frac{(\lambda +1)\xi t^{\beta }}{\lambda \Gamma (\beta +1)}\frac{1}{n!}\frac{d^n}{d{u}^n}{(1-u)}^{\alpha }{|}_{u=0}}\hfill \\ \\ {\displaystyle =\frac{(\lambda +1)\xi t^{\beta }}{\Gamma (\beta +1)}{(-1)}^{n-1}\left(\genfrac{}{}{0pt}{}{a}{n}\right)\sim \frac{\alpha (\lambda +1)\xi t^{\beta }{n}^{-\alpha -1}}{\lambda \Gamma (1-\alpha )\Gamma (1+\beta )},n\to \infty }\hfill \end{array}$$

(87)

where $\alpha \in (0,1)$ and $\beta \in (0,1]$. The fat-tailed behavior $p_{\lambda ,\xi ,n}^{\alpha ,\beta }\left(t\right)\sim n^{-\alpha -1}$ reflects the occurrence of long-range forward jumps and is equivalent to the divergence of the first moment, namely $\frac{d}{du}{\overline{p}}_{\lambda ,\xi ,n}^{\alpha ,\beta }{(u,t)|}_{u=1}={\sum }_{n=0}^{\infty }{p}_{\lambda ,\xi ,n}^{\alpha ,\beta }\left(t\right)n\to \infty$ for $\alpha \in (0,1)$. The asymptotic form (87) also contains the ‘well-scaled’ continuous-space limit ($h\to 0$: $\lambda \left(h\right)={\lambda }_0{h}^{\alpha }$ and $x\in h{\mathbb{N}}_0\to {\mathbb{R}}^{+}$), namely

$${\mathcal{P}}_{{\lambda }_0,\xi }^{\alpha ,\beta }(x,t)=\underset{h\to 0}{lim}\frac{1}{h}{p}_{{\lambda }_0{h}^{\alpha },\xi ,\frac{x}{h}}^{\alpha ,\beta }\left(t\right)\sim \frac{\alpha x^{-\alpha -1}\xi t^{\beta }}{{\lambda }_0\Gamma (1-\alpha )\Gamma (1+\beta )},t\ge 0,x\to \infty .$$

(88)

The “well-scaled” continuous-space limiting procedure will be justified in subsequent paragraph in more details. ${\lambda }_0>0$ is a new constant (independent of h) and has units ${\left[\mathrm{cm}\right]}^{-\alpha }$. The constant $\xi >0$ has physical dimension ${[sec]}^{-\beta }$ thus (88) is a spatial density of units ${\left[\mathrm{cm}\right]}^{-1}$.

Then let us also consider the asymptotic behavior for large (dimensionless) time $t{\xi }^{\frac{1}{\beta }}\to \infty$ and finite n. From the asymptotic behavior of the scalar Mittag-Leffler function $E_{\beta }(-a{t}^{\beta })\sim \frac{1}{a}\frac{t^{-\beta }}{\Gamma (1-\beta )}$ with $a=\frac{\xi (\lambda +1)}{1+\lambda {(1-u)}^{-\alpha }}$ follows for the asymptotic behavior of the generating function (80)

$${\overline{p}}_{\lambda ,\xi }^{\alpha ,\beta }(u,t)\sim \frac{(1+\lambda {(1-u)}^{-\alpha })}{(\lambda +1)\xi }\frac{t^{-\beta }}{\Gamma (1-\beta )},t{\xi }^{\frac{1}{\beta }}\to \infty ,\beta \in (0,1),\alpha \in (0,1]$$

(89)

and hence for the state probabilities

$$p_{\lambda ,\xi ,n}^{\alpha ,\beta }\left(t\right)\sim \frac{({\delta }_{n,0}+\lambda \frac{{\left(\alpha \right)}_n}{n!})}{(\lambda +1)}\frac{t^{-\beta }}{\xi \Gamma (1-\beta )},t{\xi }^{\frac{1}{\beta }}\to \infty ,n\in {\mathbb{N}}_0,\beta \in (0,1),\alpha \in (0,1]$$

(90)

where for $n=0$ the pure Mittag-Leffler asymptotic behavior of the survival probability (85) $p_{\lambda ,\xi ,0}^{\alpha ,\beta }\left(t\right)\sim \frac{t^{-\beta }}{\xi \Gamma (1-\beta )}$ is obtained. We also get the asymptotic behavior in view of (89) for the well-scaled continuous-space limit by ($h\to 0$: $\lambda \left(h\right)={\lambda }_0{h}^{\alpha }$ and $x\in h{\mathbb{N}}_0\to {\mathbb{R}}^{+}$) for large dimensionless times and finite continuous state variable $x\in [0,\infty )$

$$\begin{array}{ccc}{\displaystyle {\mathcal{P}}_{{\lambda }_0,\xi }^{\alpha ,\beta }(x,t)}& {\displaystyle \sim \underset{h\to 0}{lim}\frac{1}{h}{p}_{{\lambda }_0{h}^{\alpha },\xi ,\frac{x}{h}}^{\alpha ,\beta }\left(t\right)=\underset{h\to 0}{lim}\frac{\left({\delta }_h\left(x\right)+{\lambda }_0{h}^{\alpha -1}\frac{{\left(\alpha \right)}_{\frac{x}{h}}}{\frac{x}{h}!}\right)}{({\lambda }_0{h}^{\alpha }+1)}\frac{t^{-\beta }}{\xi \Gamma (1-\beta )}}\hfill & \\ \\ {\displaystyle {\mathcal{P}}_{{\lambda }_0,\xi }^{\alpha ,\beta }(x,t)}& {\displaystyle \sim \left(\delta \left(x\right)+{\lambda }_0\frac{x^{\alpha -1}}{\Gamma \left(\alpha \right)}\right)\frac{t^{-\beta }}{\xi \Gamma (1-\beta )},t{\xi }^{\frac{1}{\beta }}\to \infty ,x\in {\mathbb{R}}^{+},\beta \in (0,1)}\hfill \end{array}$$

(91)

where the asymptotic $t^{-\beta }$-power-law decay reflects the non-Markovian long memory feature of the process.

6.2. Well-Scaled Continuous-Space Limit

Having derived the state probabilities (83) solving Cauchy problem (77), we are now interested in the well-scaled continuous-space limit density solving a (forward) diffusion equation which turns out to refer to the general class (2). We can define this diffusion limit by the scaling assumption $\lambda \left(h\right)=h^{\alpha }{\lambda }_0$ where ${\lambda }_0>0$ is independent of h. Here the constant $\xi$ does not need to be rescaled in order to obtain an existing limit. We define the ‘well-scaled’ continuous-space limit state density kernel by

$$\begin{array}{c}{\displaystyle {\mathcal{P}}_{{\lambda }_0,\xi }^{\alpha ,\beta }(x,t)=\underset{h\to 0}{lim}{E}_{\beta }\left(\frac{-\xi ({\lambda }_0{h}^{\alpha }+1)t^{\beta }}{1+{\lambda }_0{h}^{\alpha }{(1-{\hat{T}}_{-h})}^{-\alpha }}\right){\delta }_h\left(x\right)=\underset{h\to 0}{lim}\frac{1}{h}{p}_{{\lambda }_0{h}^{\alpha },\xi ,\frac{x}{h}}^{\alpha ,\beta }\left(t\right)}\hfill \\ \\ {\displaystyle =E_{\beta }\left(\frac{-\xi t^{\beta }}{1+{\lambda }_0{D}_x^{-\alpha }}\right)\delta \left(x\right)}\hfill \end{array}$$

(92)

where we employed the discrete-$\delta$ distribution ${\delta }_h\left(x\right)$ defined in (A3) and its limiting behavior (A4) and in this limiting process $x\in h{\mathbb{N}}_0\to {\mathbb{R}}^{+}$. $D_x^{-\alpha }$ indicates the Riemann-Liouville fractional integral operator of order $\alpha$. By taking into account (83) we get

$${\mathcal{P}}_{{\lambda }_0,\xi }^{\alpha ,\beta }(x,t)=\underset{h\to 0}{lim}\sum _{m=0}^{\infty }\frac{{[-({\lambda }_0{h}^{\alpha }+1)\xi t^{\beta }]}^m}{\Gamma (\beta m+1)}\frac{1}{h}{\mathcal{E}}_{\alpha }^{\left(m\right)}\left({\lambda }_0{h}^{\alpha },\frac{x}{h}\right)$$

(93)

where we use the asymptotic relation $\frac{{\left(\beta \right)}_n}{n!}\sim \frac{n^{\beta -1}}{\Gamma \left(\beta \right)}$. We have to consider in (84) the case $0<\lambda \left(h\right)<1$ as $\lambda \left(h\right)={\lambda }_0{h}^{\alpha }\to 0$ to evaluate the continuous-space limit

$$\begin{array}{ccc}{\displaystyle e_{\alpha ,0}^m(-{\lambda }_0,x)}& {\displaystyle =\underset{h\to 0}{lim}\frac{1}{h}{\mathcal{E}}_{\alpha }^{\left(m\right)}\left({\lambda }_0{h}^{\alpha },\frac{x}{h}\right)}\hfill & \\ \\ & {\displaystyle =\underset{h\to 0}{lim}{\left[1+{\lambda }_0{h}^{\alpha }{(1-{\hat{T}}_{-h})}^{-\alpha }\right]}^{-m}{\delta }_h\left(x\right)}\hfill & \\ \\ & {\displaystyle =\underset{h\to 0}{lim}\sum _{s=0}^{\infty }\frac{{\left(m\right)}_s{(-{\lambda }_0)}^s}{s!}{h}^{\alpha s}\sum _{n=0}^{\infty }\frac{{\left(\alpha s\right)}_n}{n!}{\delta }_h(x-hn)}\hfill & \\ \\ & {\displaystyle =\underset{h\to 0}{lim}{\delta }_h\left(x\right)+\sum _{s=1}^{\infty }\frac{{\left(m\right)}_s{(-{\lambda }_0)}^s}{s!}{h}^{\alpha s-1}\frac{{\left(\alpha s\right)}_{\frac{x}{h}}}{\frac{x}{h}!}}\hfill & \\ \\ {\displaystyle e_{\alpha ,0}^m(-{\lambda }_0,x)}& {\displaystyle =\delta \left(x\right)+\sum _{s=1}^{\infty }\frac{{\left(m\right)}_s{(-{\lambda }_0)}^s{x}^{\alpha s-1}}{s!\Gamma \left(\alpha s\right)},}\hfill & m=0,1,2,..\in {\mathbb{N}}_0,x\in {\mathbb{R}}^{+}\\ \\ {\displaystyle }& {\displaystyle =D_x\left[\Theta \left(x\right)E_{\alpha ,1}^m(-{\lambda }_{0}x)\right]}\hfill \end{array}$$

(94)

and for $m=0$ we have $e_{\alpha ,0}^0(-{\lambda }_0,x)={lim}_{h\to 0}{\delta }_h\left(x\right)=\delta \left(x\right)$. In the last line we have used the property of the Heaviside step function $D_x\Theta \left(x\right)=\delta \left(x\right)$ and it comes into play the so called Prabhakar function which is a generalization of the Mittag-Leffler function introduced by Prabhakar [21] as

$$E_{\alpha ,\beta }^{\gamma }\left(z\right)=\sum _{s=0}^{\infty }\frac{{\left(\gamma \right)}_s}{s!}\frac{z^s}{\Gamma (\alpha s+\beta )},\alpha ,\beta ,\gamma \in \mathbb{C},\Re \left\{\alpha \right\}>0.$$

(95)

We used in the deduction (94) that ${\hat{T}}_{-hn}{\delta }_h\left(x\right)={\delta }_h(x-hn)=\frac{1}{h}{\delta }_{\frac{x}{h},n}$ where $h\to 0:x\in h{\mathbb{N}}_0\to {\mathbb{R}}^{+}$ (see also Appendix A, Equation (A3)). The continuous-space limit expression (94) can be identified with a so called Prabhakar kernel $e_{\alpha ,0}^m(-{\lambda }_0,x)$. The Prabhakar kernel $e_{\alpha ,\beta }^{\gamma }(\zeta ,x)$ was introduced by Giusti [60] (where we use here his definition) as the kernel with Laplace transform ${\tilde{e}}_{\alpha ,\beta }^{\gamma }(\zeta ,x)\left(s\right)=s^{-\beta }{(1-\zeta s^{-\alpha })}^{-\gamma }$. It follows that (84) indeed is a discrete approximation of Prabhakar kernel (94). With (94) and (93) we obtain for the state density kernel

$$\begin{array}{ccc}{\displaystyle {\mathcal{P}}_{{\lambda }_0,\xi }^{\alpha ,\beta }(x,t)}& {\displaystyle =\sum _{m=0}^{\infty }\frac{{(-\xi t^{\beta })}^m}{\Gamma (\beta m+1)}{e}_{\alpha ,0}^m(-{\lambda }_0,x),{\lambda }_0,\xi >0,\alpha ,\beta \in (0,1],x,t\ge 0}\hfill & \\ \\ & {\displaystyle =D_x\left[\Theta \left(x\right)\Theta \left(t\right)\sum _{m=0}^{\infty }\frac{{(-\xi t^{\beta })}^m}{\Gamma (\beta m+1)}{E}_{\alpha ,1}^m(-{\lambda }_{0}x)\right]}\hfill \end{array}$$

(96)

which has units ${\left[\mathrm{cm}\right]}^{-1}$. We directly verify the initial condition ${\mathcal{P}}_{{\lambda }_0,\xi }^{\alpha ,\beta }{(x,t)|}_{t=0}=e_{\alpha ,0}^0(-{\lambda }_0,x)=\delta \left(x\right)$.

The time-dependence of the state density kernel (96) is plotted in Figure 2 for the Poisson limit $\beta =1$ for two different values of $\alpha$. Increasing values of the state variable x are indicated by colors turning from red (small $x>0$) to blue (large x). We observe that for larger $\alpha$ in the lower plot the state density exhibits increasingly oscillating behavior for increasing values of x. We will come back to this issue subsequently when we discuss the emerging continuous-space limit diffusion equation governing the time-evolution of the state density.

It is now only a small step to derive the continuous-space forward diffusion equation of generalized fractional type which is solved by the state density kernel (96). To this end we deduce the convolution kernel of continuous-space limit of the right-hand side in (77), by the well-scaled limit (see also relations (48) and (39))

$$\begin{array}{ccc}{\displaystyle {\mathcal{G}}_{\xi ,{\lambda }_0,\alpha }\left(x\right)}& {\displaystyle =\underset{h\to 0}{lim}\frac{g_{ML,\alpha }(1-{\hat{T}}_{-h},{\lambda }_0{h}^{\alpha })}{g_{ML,\alpha }(1,{\lambda }_0{h}^{\alpha })}{\delta }_h\left(x\right)=\underset{h\to 0}{lim}\frac{({\lambda }_0{h}^{\alpha }+1)}{1+{\lambda }_0{h}^{\alpha }{(1-{\hat{T}}_{-h})}^{-\alpha }}{\delta }_h\left(x\right)}\hfill & \\ & {\displaystyle =\frac{1}{(1+{\lambda }_0{D}_x^{-\alpha })}\delta \left(x\right)=\sum _{m=0}^{\infty }{(-{\lambda }_0)}^m{D}_x^{-m\alpha }\delta \left(x\right).}\hfill \end{array}$$

(97)

Therefore

$$\begin{array}{ccc}{\displaystyle {\mathcal{G}}_{\xi ,{\lambda }_0,\alpha }\left(x\right)}& {\displaystyle =D_x\frac{D_x^{\alpha -1}}{D_x^{\alpha }+{\lambda }_0}\delta \left(x\right)=D_x[\Theta \left(x\right)E_{\alpha }(-{\lambda }_0{x}^{\alpha })]}\hfill & \\ \\ & {\displaystyle =\delta \left(x\right)+\Theta \left(x\right)D_x{E}_{\alpha }(-{\lambda }_0{x}^{\alpha })=\delta \left(x\right)-\Theta \left(x\right){\lambda }_0{x}^{\alpha -1}{E}_{\alpha ,\alpha }(-{\lambda }_0{x}^{\alpha })}\hfill & \\ \\ & {\displaystyle =\delta \left(x\right)-{\mathcal{W}}_{ML,\alpha }(x,{\lambda }_0).}\hfill \end{array}$$

(98)

We call this kernel ‘Laplacian density’. This result also is obtained by taking into account its spatial Laplace transform $\frac{s^{\alpha }}{s^{\alpha }+{\lambda }_0}$ (obtained from the scaling limit of the first line in (98)). We notice that the Laplacian density still maintains the ‘distributional versions’ of the good Laplacian properties (i)–(iii): We have ${\int }_0^{\infty }{\mathcal{G}}_{\xi ,{\lambda }_0,\alpha }\left(x\right)\mathrm{d}x=0$ corresponding to (i), ${\mathcal{G}}_{\xi ,{\lambda }_0,\alpha }\left(x\right)<0$ for $x>0$ (condition (iii)), and ${lim}_{ϵ\to 0+}{\int }_0^{ϵ}{\mathcal{G}}_{\xi ,{\lambda }_0,\alpha }\left(x\right)\mathrm{d}x=1$ corresponds to (ii). In the last line of (98) we account for the Mittag-Leffler transition density kernel ${\mathcal{W}}_{ML,\alpha }(x,{\lambda }_0)={\lambda }_0{x}^{\alpha -1}{E}_{\alpha ,\alpha }(-{\lambda }_0{x}^{\alpha })$ (60) occurring as continuous-space limit of the Mittag-Leffler matrix (54). The continuous-space-time Cauchy problem then writes

$$\begin{array}{ccc}{\displaystyle \frac{{\partial }^{\beta }}{\partial t^{\beta }}\mathcal{P}(x,t)}& {\displaystyle =-\xi {\int }_{-\infty }^{\infty }\mathcal{P}(\tau ,t){\mathcal{G}}_{\xi ,{\lambda }_0,\alpha }(x-\tau )\mathrm{d}\tau ,x,t\ge 0,\alpha ,\beta \in (0,1]}\hfill & \\ \\ & {\displaystyle =-\xi {\int }_0^x\mathcal{P}(\tau ,t){\mathcal{G}}_{\xi ,{\lambda }_0,\alpha }(x-\tau )\mathrm{d}\tau }\hfill & \\ \\ & {\displaystyle =-\xi \mathcal{P}(x,t)+\xi {\int }_0^x{\mathcal{W}}_{ML,\alpha }(x-\tau ,{\lambda }_0)\mathcal{P}(\tau ,t)\mathrm{d}\tau }\hfill & \\ \\ & {\displaystyle =-\xi \frac{\partial }{\partial x}{\int }_0^x{E}_{\alpha }(-{\lambda }_0{(x-\tau )}^{\alpha })\mathcal{P}(\tau ,t)\mathrm{d}\tau }\hfill & \\ \\ & {\displaystyle =-\xi \mathcal{P}(x,t)+\xi {\int }_0^x{\lambda }_0{(x-\tau )}^{\alpha -1}{E}_{\alpha ,\alpha }[-{\lambda }_0{(x-\tau )}^{\alpha }]\mathcal{P}(\tau ,t)\mathrm{d}\tau }\hfill & \\ \\ & {\displaystyle =-\xi \mathcal{P}(x,t)+\xi {\lambda }_0{\mathfrak{E}}_x(\alpha ,\alpha ;1;-{\lambda }_0)\mathcal{P}(x,t)}\hfill \\ \\ {\displaystyle {\mathcal{P}(x,t)|}_{t=0}}& =\delta \left(x\right)\hfill \end{array}$$

(99)

where in the second line it is used that $-D_x{E}_{\alpha }(-{\lambda }_0{x}^{\alpha })={\lambda }_0{x}^{\alpha -1}{E}_{\alpha ,\alpha }(-{\lambda }_0{x}^{\alpha })$ which yields the (weakly singular and hence integrable) Mittag-Leffler density, and where ${\mathfrak{E}}_x$ is the Prabhakar integral acting on the space variable x (see [21]). The integration limits reflect the upper triangular circulant property of the Laplacian matrix function $g_{ML,\alpha }\left(\mathcal{L}\right)$ with ${\mathcal{G}}_{\xi ,{\lambda }_0,\alpha }\left(x\right),{\mathcal{W}}_{ML,\alpha }(x,{\lambda }_0),\mathcal{P}(x,t)=0$ for $x<0$ (See also the last line of (77)). Equation (99) is the well-scaled continuous-space limit of (77) and a generalized Kolmogorov-Feller (forward) diffusion limit equation of general type (2) solved by the state density kernel (96).

Equation (99) has the following physical interpretation: The Mittag-Leffler convolution on the right-hand side is the contribution to $\mathcal{P}(x,t)$ by incoming jumps to x which originate from all states $\tau$ with $0\le \tau <x$. The term $-\xi \mathcal{P}(x,t)=-\xi \mathcal{P}(x,t){\int }_x^{\infty }{\mathcal{W}}_{ML,\alpha }(\tau -x)\mathrm{d}\tau$ describes the “loss” due to outgoing jumps from x by long-range Mittag-Leffler jumps into the infinite half space $\tau >x$.

Coming back to Figure 2: In view of the structure of Equation (99), the emergence of oscillating behavior in the time dependence of the state density, especially visible in the lower plot of this figure, is resulting from the complex interplay of incoming and outgoing jumps to and from a state x where we give a rough qualitative interpretation: For t and x both ‘small’, the incoming jumps at x are ‘mainly’ originating from the initial $\delta$-peak $\mathcal{P}(x,0)=\delta \left(x\right)$ where their accumulation giving rise to the first maximum which decays in time by outgoing jumps and the lack of incoming jumps. This effect becomes the more evanescent the larger x due to dispersion effects caused by jumps of any length drawn from the fat-tailed Mittag-Leffler density. This dispersive behavior is strongly contrasted by the stable dispersion free propagation of the $\delta$-distribution state density in the standard Poisson process (see relation (73)).

It appears worthy also to consider briefly the connection with the Montroll-Weiss CTRW picture [11]. In this picture the Cauchy problem (77) is equivalent to a strictly increasing walk on the integer line with discrete Mittag-Leffler jumps according to the one-step transition matrix (54) subordinated to a fractional Poisson process with Mittag-Leffler waiting time density (having time-Laplace transform ${\tilde{\chi }}_{\beta }\left(s\right)=\frac{\xi }{\xi +s^{\beta }}$). The generating function of the time-Laplace transform of the Cox series (10) then yields the Montroll-Weiss equation

$$\begin{array}{ccc}{\displaystyle {\overline{p}}_{\lambda ,\xi }^{\alpha ,\beta }(u,s)}& {\displaystyle =\frac{1}{s}(1-{\tilde{\chi }}_{\beta }\left(s\right))\sum _{n=0}^{\infty }{\left[{\tilde{\chi }}_{\beta }\left(s\right){\overline{W}}_{ML,\alpha }(u,\lambda )\right]}^n}\hfill & \\ \\ & {\displaystyle =\frac{1}{s}(1-{\tilde{\chi }}_{\beta }\left(s\right))\frac{1}{1-{\overline{W}}_{ML,\alpha }(u,\lambda ){\tilde{\chi }}_{\beta }\left(s\right)}}\hfill & \\ \\ & {\displaystyle =\frac{s^{\beta -1}}{s^{\beta }+\frac{g_{ML,\alpha }(1-u,\lambda )}{g_{ML,\alpha }(1,\lambda )}}=\frac{s^{\beta -1}}{s^{\beta }+\frac{\xi (\lambda +1){(1-u)}^{\alpha }}{\lambda +{(1-u)}^{\alpha }}}}\hfill \end{array}$$

(100)

where $g_{ML,\alpha }(1-u,\lambda )$ is the Laplacian generating function (48) and ${\overline{W}}_{ML,\alpha }(u,\lambda )$ the generating function of the Mittag-Leffler transition matrix (49). We identify (100) indeed with the time-Laplace transform of the Mittag-Leffler generating function (80) of the state-probabilities. It is then straight-forward to see that (100) is equivalent to (78) by rewriting Montroll-Weiss relation (100) in the form

$$s^{\beta }\overline{p}(u,s)-s^{\beta -1}\overline{p}(u,t){|}_{t=0}=-\frac{\xi (\lambda +1){(1-u)}^{\alpha }}{\lambda +{(1-u)}^{\alpha }}\overline{p}(u,s),\overline{p}(u,t){|}_{t=0}=1.$$

(101)

Inverting the time-Laplace transform on the left hand side yields the Caputo-derivative thus the causal time-domain representation of (101) indeed coincides with the generating function representation of the Cauchy problem (78). This concludes our proof of equivalence of the space-time Mittag-Leffler process with a Montroll-Weiss CTRW of a strictly increasing walk with discrete Mittag-Leffler jumps (54) subordinated to a fractional Poisson process.

7. Generalized Space-Time Mittag-Leffler Process

Finally quasi as a byproduct we consider the space-fractional generalization of the space-time Mittag-Leffler process. Let $g_{1,2}\left(\mathcal{L}\right)$ be good Laplacian matrix functions constructed by Lévy measures ${\nu }_{1,2}\left(\mathrm{d}\tau \right)$ with (38). Then a good Laplacian matrix function is obtained also by the chain $g_2\left(g_1\left(\mathcal{L}\right)\right)$. With the choice $g_1\left(\mathcal{L}\right)=g_{ML,\alpha }\left(\mathcal{L}\right)$ and $g_2\left(\mathcal{L}\right)={\mathcal{L}}^{\mu }$ where $\mathcal{L}$ is Laplacian matrix (36), we get (${\nu }_2\left(\mathrm{d}\tau \right)=-\frac{{\tau }^{-1-\mu }}{\Gamma (-\mu )}\mathrm{d}\tau$)

$$\begin{array}{ccc}{\displaystyle {\left[g_{ML,\alpha }\left(\mathcal{L}\right)\right]}^{\mu }}& {\displaystyle =-\frac{1}{\Gamma (-\mu )}{\int }_0^{\infty }(1-e^{-\tau g_{ML,\alpha }\left(\mathcal{L}\right)}){\tau }^{-\mu -1}\mathrm{d}\tau ,\mu \in (0,1)}\hfill & \\ \\ & {\displaystyle =\frac{{\mathcal{L}}^{\alpha \mu }}{{(\lambda \mathbf{1}+{\mathcal{L}}^{\alpha })}^{\mu }}.}\hfill \end{array}$$

(102)

This integral converges for $\mu \in (0,1)$ thus the fractional power ${\left[g_{ML,\alpha }\left(\mathcal{L}\right)\right]}^{\mu }$ in this range is a good Laplacian Bernstein matrix function retaining the Laplacian properties (i)–(iii). The Cauchy problem governing the state probabilities then reads

$$\begin{array}{ccc}{\displaystyle \frac{d^{\beta }}{d{t}^{\beta }}{p}_n\left(t\right)}\hfill & =& {\displaystyle -\xi {(\lambda +1)}^{\mu }\frac{{(1-{\hat{T}}_{-1})}^{\alpha \mu }}{{(\lambda +{(1-{\hat{T}}_{-1})}^{\alpha })}^{\mu }}{p}_n\left(t\right),\xi ,\lambda >0,t\ge 0\alpha ,\mu ,\beta \in (0,1],n\in {\mathbb{N}}_0}\hfill \\ \\ {\displaystyle p_n{\left(t\right)|}_{t=0}}\hfill & =& {\displaystyle {\delta }_{0,n}}\hfill \end{array}$$

(103)

with the generating function of the state probabilities

$${\displaystyle {\overline{p}}_{\lambda ,\xi }^{\alpha ,\mu ,\beta }(u,t)=E_{\beta }\left(\frac{-\xi t^{\beta }{(\lambda +1)}^{\mu }{(1-u)}^{\alpha \mu }}{{(\lambda +{(1-u)}^{\alpha })}^{\mu }}\right)=\sum _{m=0}^{\infty }\frac{{[-\xi t^{\beta }{(\lambda +1)}^{\mu }]}^m}{\Gamma (\beta m+1)}{\left[g_{ML,\alpha }(1-u)\right]}^{m\mu }.}$$

(104)

We hence get for the state-probabilities

$$p_{\lambda ,\xi ,n}^{\alpha ,\mu ,\beta }\left(t\right)=\sum _{m=0}^{\infty }\frac{{[-{(\lambda +1)}^{\mu }\xi t^{\beta }]}^m}{\Gamma (\beta m+1)}{\mathcal{E}}_{\alpha }^{\left(m\mu \right)}(\lambda ,n),n\in {\mathbb{N}}_0$$

(105)

with the cases (See also (84))

$${\mathcal{E}}_{\alpha }^{\left(m\mu \right)}(\lambda ,n)=\left\{\begin{array}{cc}{\displaystyle \frac{1}{n!}\sum _{s=0}^{\infty }\frac{{(-\lambda )}^s{\left(m\mu \right)}_s}{s!}\frac{\Gamma (\alpha s+n)}{\Gamma \left(\alpha s\right)},}& 0<\lambda <1\hfill \\ \\ {\displaystyle \frac{{(-1)}^n{\lambda }^{-m\mu }}{n!}\sum _{s=0}^{\infty }\frac{{(-1)}^s{\lambda }^{-s}{\left(m\mu \right)}_s}{s!}\frac{\Gamma \left(\alpha \right(s+m\mu )+1)}{\Gamma \left(\alpha \right(s+m\mu )-n+1)},}& \lambda >1\hfill \end{array}\right.n\in {\mathbb{N}}_0$$

(106)

where these series converge absolutely in view of the asymptotic behavior of the terms for large s, namely $\sim s^{\mu m+n-1}{\lambda }^s$ ($\lambda <1$) and $\sim s^{\mu m+n-1}{\lambda }^{-s}$ ($\lambda >1$). ${\mathcal{E}}_{\alpha }^{\left(m\mu \right)}(\lambda ,n)$ is a discrete approximation of a Prabhakar kernel. The well-scaled continuous-time limit is now straight-forwardly obtained with $\lambda \left(h\right)={\lambda }_0{h}^{\alpha }\to 0$ in the first line of (106) with $n=\frac{x}{h}$ and yields the state density

$$\begin{array}{cc}{\displaystyle {\mathcal{P}}_{{\lambda }_0,\xi }^{\alpha ,\beta }(x,t)}& {\displaystyle =\underset{h\to 0}{lim}\frac{1}{h}{p}_{{\lambda }_0{h}^{\alpha },\xi ,\frac{x}{h}}^{\alpha ,\mu ,\beta }\left(t\right)}\hfill \\ \\ & {\displaystyle =\sum _{m=0}^{\infty }\frac{{(-\xi t^{\beta })}^m}{\Gamma (\beta m+1)}{e}_{\alpha ,0}^{m\mu }(-{\lambda }_0,x),{\lambda }_0,\xi >0,\alpha ,\mu ,\beta \in (0,1],x,t\ge 0}\hfill & \\ \\ & {\displaystyle =D_x\left[\Theta \left(x\right)\Theta \left(t\right)\sum _{m=0}^{\infty }\frac{{(-\xi t^{\beta })}^m}{\Gamma (\beta m+1)}{E}_{\alpha ,1}^{m\mu }(-{\lambda }_{0}x)\right]}\hfill \end{array}$$

(107)

with the Prabhakar kernel $e_{\alpha ,0}^{m\mu }(-{\lambda }_0,x)$ obtained as (see (94))

$${\displaystyle e_{\alpha ,0}^{m\mu }(-{\lambda }_0,x)=\delta \left(x\right)+\sum _{s=1}^{\infty }\frac{{\left(m\mu \right)}_s{(-{\lambda }_0)}^s{x}^{\alpha s-1}}{s!\Gamma \left(\alpha s\right)},m=0,1,2,..\in {\mathbb{N}}_0,x\in {\mathbb{R}}^{+}.}$$

(108)

Then we further get for the Laplacian density in the same way as in (98) the kernel

$${\mathcal{G}}_{\xi ,{\lambda }_0,\alpha ,\mu }\left(x\right)=e_{\alpha ,0}^{\mu }(-{\lambda }_0,x)=\delta \left(x\right)-{\mathcal{W}}_{\alpha }^{\left(\mu \right)}\left(x\right)=D_x[\Theta \left(x\right)E_{\alpha ,1}^{\mu }(-{\lambda }_0{x}^{\alpha })]$$

(109)

with $E_{\alpha ,1}^{\mu }(-{\lambda }_0{x}^{\alpha })={\sum }_{s=0}^{\infty }\frac{{\left(\mu \right)}_s{(-{\lambda }_0)}^s{x}^{\alpha s}}{s!\Gamma (\alpha s+1)}$. The transition density kernel hence yields

$${\mathcal{W}}_{\alpha }^{\left(\mu \right)}\left(x\right)=-D_x{E}_{\alpha ,1}^{\mu }(-{\lambda }_0{x}^{\alpha })={\lambda }_0{x}^{\alpha -1}\sum _{s=0}^{\infty }\frac{{\left(\mu \right)}_{s+1}}{(s+1)!}\frac{{(-{\lambda }_0{x}^{\alpha })}^s}{\Gamma (\alpha s+\alpha )},x>0$$

(110)

recovering for $\mu =1$ the Mittag-Leffler density. It is easily verified that the transition kernel is normalized ${\int }_0^{\infty }{\mathcal{W}}_{\alpha }^{\left(\mu \right)}\left(x\right)\mathrm{d}x=E_{\alpha ,1}^{\mu }(-{\lambda }_0{x}^{\alpha }){|}_{x=0}=1$. The diffusion-limit Cauchy problem which again is of general type (2) then reads

$${\displaystyle \frac{{\partial }^{\beta }}{\partial t^{\beta }}\mathcal{P}(x,t)=-\xi \mathcal{P}(x,t)+\xi {\int }_0^x{\mathcal{W}}_{\alpha }^{\left(\mu \right)}{(x-\tau )\mathcal{P}(\tau ,t)\mathrm{d}\tau ,\mathcal{P}(x,t)|}_{t=0}=\delta \left(x\right)}$$

(111)

and is solved by the state density kernel (107). For $\mu =1$ all expressions turn into those previously derived for the space-time Mittag-Leffler process.

The fractional generalization of the space-time Mittag-Leffler process shows that the analysis of processes generated by chains of Bernstein matrix functions $g_n(..g_1\left(\mathcal{L}\right))$ which again are of the class of Bernstein functions retaining the good Laplacian properties may open an interesting direction to be further explored.

8. Conclusions

In this paper, we have analyzed space-time generalizations of the Poisson process defined by Cauchy problems with generalized Kolmogorov-Feller difference-differential equations of type (1). These generalizations in the Montroll-Weiss CTRW picture are strictly increasing walks on the integer line time-changed with an independent renewal process. We have shown that this approach can also be applied more generally to biased walks on digraphs and requires the construction of non-trivial ‘good Laplacian matrix functions’ $g\left(\mathcal{L}\right)$. This introduces new topologies of fully connected structures with the small world property if the Laplacian matrix $\mathcal{L}$ is ergodic.

Choosing Lévy densities related to normalized continuous distributions and a Laplacian $\mathcal{L}$ of the trivial strictly increasing walk (36) leads to non-trivial discrete-approximations of these Lévy densities (see (38)–(42)) in the form of upper triangular circulant transition matrices with strictly positive elements above the main diagonal, thus allowing any positive integer jump. These properties are especially useful to construct new space-time generalizations of the Poisson process.

As a pertinent example we have derived in this way a ‘good Laplacian matrix function’ which generates a strictly increasing walk with discrete Mittag-Leffler jumps and introduced the space-time Mittag-Leffler process. For this process, by means of explicit formulae, we derived the state probabilities (Equation (83)) solving the Cauchy problem (77). Further, we developed a well-scaled continuous space limiting procedure and obtained the state density (96) as a limiting expression of the state probabilities where Prabhakar kernels come into play. We derived then the forward diffusion equation of general fractional type (99) which involves the Prabhakar integral and is solved by the state density (96). The continuous-space limit diffusion Equation (99) refers to the general class of generalized Kolmogorov-Feller forward Equation (2). We also showed that the space-time Mittag-Leffler process in the Montroll-Weiss picture is a CTRW with discrete Mittag-Leffler jumps subordinated to an independent fractional Poisson process.

Following this line, we introduced a space-fractional generalization of the space-time Mittag-Leffler process constructed by a chain of two Laplacian Bernstein functions retaining the good Laplacian properties. In this way, we derived the Cauchy problem (103) governing the state probabilities for which we obtained an explicit formula (Equation (105)). We also deduced the well-scaled state density kernel (Equation (107)) which involves Prabhakar kernels and solves the continuous-space Cauchy problem (111).

There is a large potential in the presented approach of constructing space-time generalizations of the Poisson process in terms of strictly increasing walks. Also, new general types of biased walks time-changed with continuous-time or discrete-time counting processes derived in this way may have interesting applications in ‘birth-and-death’ models. On the other hand, construction of new processes involving Prabhakar distributions (see, e.g., [45]) seem to be promising candidates due to their connections to the dynamics of certain complex phenomena. Further applications in the field of ‘aging of complex systems’ which are characterized by strictly increasing random accumulation of damage (‘misrepair’) measures (see [61] for a Markovian model) could open an interesting interdisciplinary field as well.