On the fractional Poisson process and the discretized stable subordinator

The fractional Poisson process and the Wright process (as discretization of the stable subordinator) along with their diffusion limits play eminent roles in theory and simulation of fractional diffusion processes. Here we have analyzed these two processes, concretely the corresponding counting number and Erlang processes, the latter being the processes inverse to the former. Furthermore we have obtained the diffusion limits of all these processes by well-scaled refinement of waiting times and jumps


Introduction
Serious studies of the fractional generalization of the Poisson processreplacement of the exponential waiting time distribution by a distribution given via a Mittag-Leffler function with modified argument -have been started around the turn of the millenium, and since then many papers on its various aspects have appeared. There are in the literature many papers on this generalization where the authors have outlined a number of aspects and definitions, see e.g. Repin and Saichev (2000) [44], Jumarie (2001) [22], Wang et al. (2003Wang et al. ( ,2006 [52,53], Laskin (2003Laskin ( ,2009 [28,29],  [32], Uchaikin et al. (2008) [51], Beghin and Orsingher (2009) [3], Cahoy et al. (2010) [5], Meerschaert et al.(2011) [41], Politi et al. (2011) [45], Kochubei (2012) [27], so that it seems impossible to list them all exhaustively. However, in effect this generalization was used already in 1995: Hilfer and Anton [21] (without saying it in our words) showed that the Fractional Kolmogorov-Feller equation (replacement of the first order time derivative by a fractional derivative of order between 0 and 1) requires the underlying random walk to be subordinated to a renewal process with Mittag-Lefler waiing time. Recently, in the context of time change and subordination, the fractional Poisson process has obtained revived interest, see Beghin and Orsingher (2009) [3] and Meerschaert et al.(2011) [41].
Here we will present our formalism for obtaining the essential characteristics of a generic renewal process and apply it to get those of the fractional Poisson counting process and its inverse, the fractional Erlang process. Both of these comprise as limiting cases the corresponding well-known non-fractional processes that are based on exponential waiting time. Then we will analyze an alternative renewal process, that we call the "Wright process", investigated by Mainardi et al (2000), (2005), (2007) [31,33,34], a process arising by discretization of the stable subordinator. In it the so-called M -Wright function plays the essential role. A scaled version of this process has been used by Barkai (2002) [2] for approximating the time-fractional diffusion process directly by a random walk subordinated to it (executing this scaled version in natural time), and he has found rather poor convergence in refinement. In   [19] we have modified the way of using this discretized stable subordinator. By appropriate discretization of the relevant spatial stable process we have then obtained a simulation method equivalent to the solution of a pair of Langevin equations, see Fogedby (1994) [9] and Kleinhans and Friedrich(2007) [26]. For simulation of space-time fractional diffusion one so obtains a sequence of precise snapshots of a true particle trajectory, see for details   [19], and also Gorenflo andMainardi (2011, 2012) [16,17].
However, we should note that already in the Sixties of the past century, Gnedenko and Kovalenko (1968) [11] obtained in disguised form the fractional Poisson process by properly rescaled infinite thinning (rarefaction) of a renewal process with power law waiting time. By "disguised" we mean that they found the Laplace transform of the Mittag-Leffler waiting time density, but being ignorant of the Mittag-Leffler function they only presented this Laplace transform. The same ignorance of the Mittag-Leffler function we again meet in a 1985 paper by Balakrishnan [1], who exhibited the Mittag-Leffler waiting time density in Laplace disguise as essential for approximating time-fractional diffusion for which he used the description in form of a fractional integrodifferential equation. We have shown that the Mittag-Leffler waiting time density in a certain sense is asymptotically universal for power law renewal processes, see Gorenflo and Mainardi (2008) [15], Gorenflo (2010) [12].
The structure of our paper is as follows. In Section 2 we discuss the elements of the general renewal theory and the CTRW concept. In Section 3 we introduce the Poisson process and its fractional generalization then, in Section 4, the socalled Wright process related to the stable subordinator and its discretization. For both processes we consider the corresponding inverse processes, the Erlang processes. In Section 5 we briefly discuss the diffusion limit for all the above processes. Finally, Section 6 is devoted to conclusions. We have devoted the Appendix A to notations and terminology, in particular to the essentials on operators, integral transforms and special functions required to understand our analysis.We have also devoted the Appendix B for providing the reader with an overview of the basic results.
In this paper we disregard the aspects of subordination, thereby we will not touch our method of parametric subordination for which we refer the reader to our papers [16,17,19].

Elements of renewal theory and CTRW
For the reader's convenience let us here present a brief introduction to renewal theory including the basics of continuous random walk (CTRW).
The general renewal process. By a renewal process we mean an infinite sequence 0 = t 0 < t 1 < t 2 < · · · of events separated by i.i.d. (independent and identically distributed) random waiting times T j = t j − t j−1 , whose probability density φ(t) is given as a function or generalized function in the sense of Gel'fand and Shilov [10] (interpretable as a measure) with support on the positive real axis t ≥ 0, non-negative: φ(t) ≥ 0, and normalized: ∞ 0 φ(t) dt = 1, but not having a delta peak at the origin t = 0. The instant t 0 = 0 is not counted as an event. An important global characteristic of a renewal process is its mean waiting time T = ∞ 0 t φ(t) dt. It may be finite or infinite. In any renewal process we can distinguish two processes, namely the counting number process and the process inverse to it, that we call the Erlang process. The instants t 1 , t 2 , t 3 , . . . are often called renewals. In fact renewal theory is relevant in practice of maintenance or required exchange of failed parts, e.g., light bulbs.
The counting number process and its inverse. We are interested in the counting number process x = N = N (t) N (t) := max {n|t n ≤ t} = n for t n ≤ t < t n+1 , n = 0, 1, 2, · · · , (2.1) where in particular N (0) = 0. We ask for the counting number probabilities in n, evolving in t, p n (t) := P[N (t) = n] , n = 0, 1, 2, · · · . (2.2) We denote by p(x, t) the sojourn density for the counting number having the value x. For this process the expectation is p n (t) δ(x − n), see (2.12)] It provides the mean number of events in the half-open interval (0, t], and is called the renewal function, see e.g. [47]. We also will look at the process t = t(N ), the inverse to the process N = N (t), that we call the Erlang process in honour of the Danish telecommunication engineer A.K. Erlang (1878-1929), see Brockmeyer et al. (1948) [4]. It gives the value of time t = t N of the N -th renewal. We ask for the Erlang probability densities q n (t) = q(t, n) , n = 0, 1, 2, . . .
For every n the function q n (t) = q(t, n) is a density in the variable of time having value t in the instant of the n-th event. Clearly, this event occurs after n (original) waiting times have passed, so that q n (t) = φ * n (t) with Laplace transform q n (s) = ( φ(s) n ) .

(2.5)
In other words the function q n (t) = q(t, n) is a probability density in the variable t ≥ 0 evolving in the variable x = n = 0, 1, 2, ....
The continuous random walk. A continuous time random walk (CTRW) is given by an infinite sequence of spatial positions 0 = x 0 , x 1 , x 2 , · · ·, separated by (i.i.d.) random jumps X j = x j − x j−1 , whose probability density function w(x) is given as a non-negative function or generalized function (interpretable as a measure) with support on the real axis −∞ < x < +∞ and normalized: ∞ 0 w(x) dx = 1, this random walk being subordinated to a renewal process so that we have a random process x = x(t) on the real axis with the property x(t) = x n for t n ≤ t < t n+1 , n = 0, 1, 2, · · ·.
We ask for the sojourn probability density u(x, t) of a particle wandering according to the random process x = x(t) being in point x at instant t.
Let us define the following cumulative probabilities related to the probability density function φ(t) For definiteness, we take Φ(t) as right-continuous, Ψ(t) as left-continuous. When the nonnegative random variable represents the lifetime of a technical system, it is common to call Φ(t) := P (T ≤ t) the failure probability and Ψ(t) := P (T > t) the survival probability, because Φ(t) and Ψ(t) are the respective probabilities that the system does or does not fail in (0, t]. These terms, however, are commonly adopted for any renewal process.
In the Fourier-Laplace domain we have and the famous Montroll-Weiss solution formula for a CTRW, see [42,54] u(κ, s) .
In our special situation the jump density has support only on the positive semi-axis x ≥ 0 and thus, by replacing the Fourier transform by the Laplace transform we obtain the Laplace-Laplace solution .
Recalling from Appendix the definition of convolutions, in the physical domain we have for the solution u(x, t) the Cox-Weiss series, see [6,54], This formula has an intuitive meaning: Up to and including instant t, there have occurred 0 jumps, or 1 jump, or 2 jumps, or . . ., and if the last jump has occurred at instant t < t, the wanderer is resting there for a duration t − t .
Renewal process as a special CTRW An essential trick of what follows is that we treat renewal processes as continuous time random walks with waiting time density φ(t) and special jump density w(x) = δ(x − 1) corresponding to the fact that the counting number N (t) increases by 1 at each positive event instant t n . We then have w(κ) = exp (−κ) and get for the counting number process N (t) the sojourn density in the transform domain (s ≥ 0, κ ≥ 0), (2.11) From this formula we can find formulas for the renewal function m(t) and the probabilities p n (t) = P {N (t) = n}. Because N (t) assumes as values only the non-negative integers, the sojourn density p(x, t) vanishes if x is not equal to one of these, but has a delta peak of height p n (t) for x = n (n = 0, 1, 2, 3, · · ·). Hence Inverting (2.11) with respect to κ and s as (2.14) According to the theory of Laplace transform we conclude from Eqs. (2.2) and (2.12) a result naturally expected, and Thus we have found in the Laplace domain the reciprocal pair of relationships , (2.17) saying that the waiting time density and the renewal function mutually determine each other uniquely. The first formula of Eq. (2.17) can also be obtained as the value at κ = 0 of the negative derivative for κ = 0 of the last expression in Eq. (2.11). Eq. (2.17) implies the reciprocal pair of relationships in the physical domain The first of these equations usually is called the renewal equation.
Considering, formally, the counting number process N = N (t) as CTRW (with jumps fixed to unit jumps 1), N running increasingly through the non-negative integers x = 0, 1, 2, ..., happening in natural time t ∈ [0, ∞), we note that in the Erlang process t = t(N ), the roles of N and t are interchanged. The new "waiting time density" now is w(x) = δ(x − 1), the new "jump density" is φ(t).
It is illuminating to consciously perceive the relationships for t ≥ 0, n = 0, 1, 2, . . ., between the counting number probabilities p n (t) and the Erlang densities q n (t). For Eq. (2.5) we have q n (t) = φ * n (t), and then by (2.14) We can also express the q n in another way by the p n . Introducing the All this is true for n = 0 as'well, by the empty sum convention n k=1 T k = 0 for n = 0.

The Poisson process and its fractional generalization
The most popular renewal process is the Poisson process. It is (uniquely) characterized by its mean waiting time 1/λ (equivalently by its intensity λ), which is a given positive number, and by its residual waiting time Ψ(t) = exp (−λt) for t ≥ 0, which corresponds to the waiting time density φ(t) = λ exp (−λt). With λ = 1 we have what we call the standard Poisson process. The general Poisson process arises from the standard one by rescaling the time variable t.
We generalize the standard Poisson process by replacing the exponential function by a function of Mittag-Leffler type. With t ≥ 0 and a parameter β ∈ (0, 1] we take We call this renewal process of Mittag-Leffler type the fractional Poisson process, see e.g. [3,5,18,28,32,41,44,45,49,51], or the Mittag-Leffler renewal process or the Mittag-Leffler waiting time process.
To analyze it we go into the Laplace domain where we have If there is no danger of misunderstanding we will not decorate Ψ and φ with the index β. The special choice β = 1 gives us the standard Poisson process with Ψ 1 (t) = φ 1 (t) = exp (−t).
Whereas the Poisson process has finite mean waiting time (that of its standard version is equal to 1), the fractional Poisson process (0 < β < 1 ) does not have this property. In fact, Let us calculate the renewal function m(t). Inserting φ(s) = 1/(1 + s β ) into Eq. (2.11) and taking w(x) = δ(x − 1) as in Section 2, we find for the sojourn density of the counting function N (t) the expressions and and then . (3.6) Using E β (0) = 1/Γ(1 + β) now yields , 0 < β < 1 . Using general Taylor expansion and, by comparison with Eq. (2.12), the counting number probabilities Observing from Eq. (3.4) and inverting with respect to κ, we finally identify En passant we have proved an often cited special case of an inversion formula by   [46], Eq. (1.80).
For the Poisson process with intensity λ > 0 we have a well-known infinite system of ordinary differential equations (for t ≥ 0), see e.g. Khintchine [23], with initial conditions p n (0) = 0, n = 1, 2, . . ., which sometimes even is used to define the Poisson process. We have an analogous system of fractional differential equations for the fractional Poisson process. In fact, from Eq. (3.13) we have with initial conditions p n (0) = 0, n = 1, 2, . . ., where * D β t denotes the time-fractional derivative of Caputo type of order β, see Appendix A. It is also possible to introduce and define the fractional Poisson process by this difference-differential system. Let us note that by solving the system (3.17), Beghin and Orsingher in [3] introduce what they call the "first form of the fractional Poisson process" , and in [41] Meerschaert et al. show that this process is a renewal process with Mittag-Leffler waiting time density as in (3.1), hence is identical with the fractional Poisson process.
Up to now we have investigated the fractional Poisson counting process N = N (t) and found its probabilities p n (t) in Eq. (3.10). To get the corresponding Erlang probability densities q n (t) = q(t, n), densities in t, evolving in n = 0, 1, 2 . . ., we find by Eq. (2.21) via telescope summation We leave it as an exercise to the readers to show that in Eq. (3.9) interchange of differentiation and summation is allowed.

The stable subordinator and the Wright process
Let us denote by g β (t) the extremal Lévy stable density of order β ∈ (0, 1] and support in t ≥ 0 whose Laplace transform is g β (s) = exp (−s β ), that is The topic of Lévy stable distributions is treated in several books on probability and stochastic processes, see e.g. Feller (1971) [8], Sato (1999) [48]; an overview of the analytical and graphical aspects of the corresponding densities is found in Mainardi et al (2001) [35], where an ad hoc notation is used.
We note that the stable density (4.1) can be expressed in terms of a function of the Wright type. In fact, with the M-Wright function from Appendix A of this paper (see Appendix F of Mainardi's book [30] for more details), we have The renewal process with waiting time density We call this process the Wright renewal process because the corresponding survival function Ψ(t) and the waiting time density φ(t) are expressed in terms of certain Wright functions. So we distinguish it from the so called Mittag-Leffler renewal process, treated in the previous Section as fractional Poisson process. More precisely, recalling the Wright functions from the Appendix A, we have for t ≥ 0, where Θ denotes the unit step Heaviside function.
It is relevant to note the Laplace transform connecting the two transcendental functions M β and E β By the stable subordinator of order β ∈ (0, 1] we mean the stochastic process t = t(x) that has sojourn density in t ≥ 0, evolving in x ≥ 0 provided by the Laplace transform correspondence, This process is monotonically increasing: for this reason it is used in the context of time change and subordination in fractional diffusion processes.
We discretize the process t = t(x) by restricting x to run through the integers n = 0, 1, 2, . . . . The resulting discretized version is a renewal process happening in pseudo-time x ≥ 0 with jumps in pseudo-space t ≥ 0 having Because here the waiting time density is given by a function of Wright type we call this process the Wright renewal process, or simply the Wright process. Immediately we get its Erlang densities (in t ≥ 0, evolving in x = n = 0, 1, 2, . . .) q n (t) = φ * n (t) ÷ e −ns β , (4.9) so that, in view of (4.7) with x = n, q n (t) = f (t, n) = n −1/β g β n −1/β t , (4.10) In the special case β = 1 we have q n (t) = δ(t − n).
We observe that this counting process gives us precise snapshots at x = 0, 1, 2, of the stable subordinator t = t(x).
With the probability distribution function (4.14) In the limiting case β = 1 we have as a function continuous from the right, and we calculate p n (t) = 0 for 0 < t < n , and for t ≥ n + 1 , 1 for n ≤ t < n + 1 . For the renewal function we obtain its Laplace transform from (2.17) so that We do not know an explicit expression for this sum if 0 < β < 1. However, in the limiting case β = 1 we obtain Using (4.17) we investigate the asymptotic behaviour of m(t) for t → ∞. We have for s → 0 m(s) ∼ 1/s 1+β and thus, by Tauber theory, see e.g. Feller (1971) [8], Remember, for the fractional Poisson process, we had found for all t ≥ 0 . Remark A rescaled version of the discretized stable subordinator can be used for producing closely spaced precise snapshots of a true particle trajectory of a space-time fractional diffusion process, see e.g. the recent chapter by Gorenflo and Mainardi (2011) [17] on parametric subordination.

The diffusion limits for the fractional Poisson and the Wright processes
In a CTRW we can, with positive scaling factor h and τ , replace the jumps X by jumps X h = h X, the waiting times T by waiting times T τ = τ T . This leads to the rescaled jump density w h (x) = w(x/h)/h and the rescaled waiting time density φ τ (t) = φ(t/τ )/τ and correspondingly to the transforms w h (κ) = w(hκ), φ τ (s) = φ(τ s).
For the sojourn density u h,τ (x, t), density in x evolving in t, we obtain from (2.9) in the transform domain where, if w(x) has support on x ≥ 0 we can work with the Laplace transform instead of the Fourier transform (replace the by ). If there exists between h and τ a scaling relation R (to be introduced later) under which u(x, t) tends for h → 0, τ → 0 to a meaningful limit v(x, t) = u 0,0 (x, t), then we call the process x = x(t) with this sojourn density a diffusion limit. We find it via v(κ, s) = lim and Fourier-Laplace (or Laplace-Laplace) inversion.
Note: this diffusion limit is a limit in the weak sense (convergence in distribution of the CTRW to the diffusion limit). The mathematical background consists in the application of the Fourier (or Laplace) continuity theorem of probability theory for fixed time t.
We will now find that the fractional Poisson process and the Wright counting process have the same diffusion limit, namely the inverse stable subordinator. The two corresponding Erlang processes have the same diffusion limit, namely the stable subordinator. For t → ∞ the renewal functions have the same asymptotic behaviour, namely m(t) ∼ t β /Γ(1 + β) which, in the case of the fractional Poisson process, is exact for all t ≥ 0.
To prove these statements we need the Laplace transform of the relevant functions φ(t) and w(x). For the fractional Poisson process we have For the Wright process we have In all cases we have, for fixed s and κ φ(τ s) ∼ 1 − (τ s) β as τ → 0 , w(hκ) ∼ 1 − (hκ) as h → 0 , . and straightforwardly we obtain for the sojourn densities in both cases, by use of (5.1) with p in place of u and replaced by Using the scaling relation R h = τ β , we obtain By partial Laplace inversions we get two equivalent representations leading to the density of the inverse stable subordinator where M β and J 1−β t denote respectively the M -Wright function and the Riemann-Liouville fractional integral introduced in Appendix A, and f (t, x) the stable subordinator given by Eq. (4.7).
Remark: In (4.7 and (5.7) the densities of the stable and the inverse stable subordinator are both represented via the M -Wright function.
The diffusion limit for the Erlang process. In the Erlang process the roles of space and time, likewise of jumps and waiting times, are interchanged. In other words we treat x ≥ 0 as a pseudo-time variable and t ≥ 0 as a pseudo-space variable. For the resulting sojourn density q(t, x), we have from interchanging in (5.1) for h → 0 and τ → 0, Again using the scaling relation R in Eq. (5.4) we find q 0,0 (s, κ) = 1 κ + s β , (5.5 ) which is the Laplace-Laplace transform of the density of stable subordinator of Section 4. In fact, by partial Laplace inversion, and it follows that See (4.7) for its explicit representation as a rescaled stable density expressed via a M -Wright function.
We get the same result by continualization of the discretized stable subordinator. Replace in Eqs. (4.9), (4.10) the discrete variable n by the continuous varaiable x.

Conclusions
The fractional Poisson process and the Wright process (as discretization of the stable subordinator) along with their diffusion limits play eminent roles in theory and simulation of fractional diffusion processes. Here we have analyzed these two processes, concretely the corresponding counting number and Erlang processes, the latter being the processes inverse to the former. Furthermore we have obtained the diffusion limits of all these processes by well-scaled refinement of waiting times and jumps.

Fourier and Laplace transforms
By IR (IR + , IR + 0 ) we mean the set of all (positive, non-negative) real numbers, and by C the set of complex numbers. It is known that the Fourier transform is applied to functions defined in L 1 (IR ) whereas the Laplace transform is applied to functions defined in L loc (IR + ). In our cases the arguments of the original function are the space-coordinate x (x ∈ IR or x ∈ IR + 0 ) and the time-coordinate t (t ∈ IR + 0 ). We use the symbol ÷ for the juxtaposition of a function with its Fourier or Laplace transform. A look at the superscript for the Fourier transform, for the Laplace transform reveals their relevant juxtaposition. We use x as argument (associated to real κ) for functions Fourier transformed, and x or t as argument (associated to complex κ or s, respectively) for functions Laplace transformed.
The meaning of the connective * will be clear from the context. For convolution powers we have: where δ denotes the Dirac generalized function.

Mittag-Leffler and Wright functions
The Mittag-Leffler function of parameter α is defined as It is entire of order 1/α. Let us note the trivial cases Without changing the order 1/α the Mittag-Leffler function can be generalized by introducing an additional (arbitrary) parameter β.
The Mittag-Leffler function of parameters α, β is defined as Laplace transforms of the Mittag-Leffler functions function For our purposes we will need, with 0 < ν ≤ 1 and t ≥ 0, the Laplace transform pairs Of high relevance is the algebraic decay of Ψ(t) and φ(t) as t → ∞: is the solution of the fractional relaxation equation with the Riemann-Liouville derivative We refer to [12,13,15] for the relevance of Mittag-Leffler functions in theory of continuous time random walk and space-time fractional diffusion and in asymptotics of power laws.

Laplace transforms of the Wright functions
For the Wright function of the first kind, being entire of exponential type, the Laplace transform can be obtained by transforming the power series term by term: For the Wright function of the second kind, denoting ν = |λ| ∈ (0, 1) we have with µ > 0 for simplicity, we have W −ν,µ (−t) ÷ E ν,µ+ν (−s) , 0 < ν < 1 .
We note the minus sign in the argument in order to ensure the the existence of the Laplace transform thanks to the Wright asymptotic formula valid in a sector about the negative real axis.