Poisson-type processes governed by fractional and higher-order recursive differential equations

We consider some fractional extensions of the recursive differential equation governing the Poisson process, by introducing combinations of different fractional time-derivatives. We show that the so-called"Generalized Mittag-Leffler functions"(introduced by Prabhakar [20]) arise as solutions of these equations. The corresponding processes are proved to be renewal, with density of the intearrival times (represented by Mittag-Leffler functions) possessing power, instead of exponential, decay, for t tending to infinite. On the other hand, near the origin the behavior of the law of the interarrival times drastically changes for the parameter fractional parameter varying in the interval (0,1). For integer values of the parameter, these models can be viewed as a higher-order Poisson processes, connected with the standard case by simple and explict relationships.


Introduction
Many well-known differential equations have been extended by introducing fractional-order derivatives with respect to time (for instance, the heat, wave and telegraph equations, as well as the higher-order heat-type equations) or with respect to space (for instance, the equations involving the Riesz fractional operator).
Fractional versions of the Poisson processes have been already presented and studied in the literature: in [9] the so-called fractional master equation was considered.A similar model was treated in [12], where the equation governing the probability distribution of the homogeneous Poisson process was modified, by introducing the Riemann-Liouville fractional derivative.The results are given in analytical form, in terms of infinite series or successive derivatives of Mittag-Leffler functions.We recall the definition of the (two-parameters) Mittag-Leffler function: (see [20], 1.2).
A different definition of Poisson fractional process has been proposed by Wang and Wen [30] and successively studied in [31]- [32]: in analogy to the well-known fractional Brownian motion, the new process is defined as a stochastic integral with respect to the Poisson measure.It displays properties similar to fractional Brownian motion, such as self-similarity (in the wide-sense) and long-range dependence.
Another approach was followed by Repin and Saichev [23]: they start by generalizing, in a fractional sense, the distribution of the interarrival times U j between two Poisson events.This is expressed, in terms of Mittag-Leffler functions, as follows, for ν ∈ (0, 1] : and coincides with the solution to the fractional equation where δ(•) denotes the Dirac delta function and again the fractional derivative is intended in the Riemann-Liouville sense.For ν = 1 formula (1.2) reduces to the well-known density appearing in the case of a homogeneous Poisson process, N (t), t > 0, with intensity λ = 1, i.e. f (t) = e −t .The same approach is followed by Mainardi et al. [14]- [15]- [16], where a deep analysis of the related process is performed: it turns out to be a true renewal process, while it looses the Markovian property.Their first step is the study of the following fractional equation (instead of (1.3)) with initial condition ψ(0 + ) = 1 and with fractional derivative defined in the Caputo sense.The solution ψ(t) = E ν,1 (−t ν ) to (1.4) represents the survival probability of the fractional Poisson process.As a consequence its probability distribution is expressed in terms of derivatives of Mittag-Leffler functions, while the density of the k-th event waiting time is a fractional generalization of the Erlang distribution and coincides with the k-fold convolution of (1.2).The analysis carried out by Beghin and Orsingher [2] starts, as in [11], from the generalization of the equation governing the Poisson process, where the time-derivative is substituted by the fractional derivative (in the Caputo sense) of order ν ∈ (0, 1]: with initial conditions and p −1 (t) = 0.The main result is the expression of the solution as the distribution of a composed process represented by the standard, homogeneous Poisson process N (t), t > 0 with a random time argument T 2ν (t), t > 0 as follows: N ν (t) = N (T 2ν (t)), t > 0.
The process T 2ν (t), t > 0 (independent of N ) possesses a well-known density, which coincides with the solution to a fractional diffusion equation of order 2ν (see (2.15) below).In the particular case where ν = 1/2 this equation coincides with the heat-equation and the process representing time is the reflected Brownian motion.These results are reconsidered here, in the next section, from a different point of view, which is based on the use of the Generalized Mittag-Leffler (GML) function.The latter is defined as , α, β, γ ∈ C, Re(α), Re(β), Re(γ) > 0, (1.6) where (γ) r = γ(γ + 1)...(γ + r − 1) (for r = 1, 2, ..., and γ = 0) is the Pochammer symbol and (γ) 0 = 1.The GML function has been extensively studied by Saxena et al. (in [25]- [26]- [27] and [28]) and applied in connection with some fractional diffusion equations, whose solutions are expressed as infinite sums of (1.6).For some properties of (1.6), see also [29].We note that formula (1.6) reduces to (1.1) for γ = 1.By using the function (1.6) it is possible to write down in a more compact form the solution to (1.5), as well as the density of the waiting-time of the k-th event of the fractional Poisson process.As a consequence some interesting relationships holding between the Mittag-Leffler function (1.1) and the GML function (1.6) are obtained.
Moreover the use of GML functions allows us to derive an explicit expression for the solution of the more complicated recursive differential equation, where two fractional derivatives appear: for ν ∈ (0, 1].As we will see in section 3, even in this case we can define a process governed by (1.7), which turns out to be a renewal.The density of the interarrival times are no-longer expressed by standard Mittag-Leffler functions as in the first case, but the use of GML functions is required and the same is true for the k-th event waiting-time.An interesting relationship between the two models analyzed here can be established by observing that the waiting-time of the k-th event of the process governed by (1.7) coincides in distribution with the waiting time of the (2k)-th event for the first model.This suggests to interpret our second model as a fractional Poisson process of the first type, which jumps upward at even-order events A 2k and the probability of the successive odd-indexed events A 2k+1 is added to that of A 2k .Correspondingly, the distribution of this second process N ν (t), t > 0, can be expressed, in terms of the processes N and T 2ν , as follows: We also study the probability generating functions of the two models, which are themselves solutions to fractional equations; in particular in the second case an interesting link with the fractional telegraph-type equation is explored.
For ν = 1, equation (1.7) takes the following form and the related process can be regarded as a standard Poisson process with Gamma-distributed interarrival times (with parameters λ, 2).This is tantamount to attribute the probability of oddorder values A 2k+1 of a standard Poisson process to the events labelled by 2k.Moreover it should be stressed that, in this special case, the equation satisfied by the probability generating function coincides with that of the damped oscillations.
All the previous results are further generalized to the case n > 2 in the concluding remarks: the structure of the process governed by the equation where ν ∈ (0, 1], is exactly the same as before and all the previous considerations can be easily extended.
2 First-type fractional recursive differential equation

The solution
We begin by considering the following fractional recursive differential equation with p −1 (t) = 0, subject to the initial conditions We apply in (2.1) the definition of fractional derivative in the sense of Caputo, that is, for m ∈ N, We note that, for ν = 1, (2.1) coincides with the equation governing the homogeneous Poisson process with intensity λ > 0, so that our first result generalizes the well-known distribution holding in the standard case, i.e. p k (t) = (λt) k k! e −λt , k ≥ 0, t > 0.
We will obtain the solution to (2.1)-(2.2) in terms of GML functions (defined in (1.6)) and show that it represents a true probability distribution of a process, which we will denote by N ν (t), t > 0 : therefore we will write Proof We take the Laplace transform of equation (2.1) together with the condition (2.2) and consider that, for the Laplace transform of the Caputo derivative, the following expression holds: , where m = ν + 1. Since, in this case, ν ∈ (0, 1] , we get m = 1.Therefore we get the following recursive formula, for k ≥ 1, while, for k = 0, we obtain since equation (2.1) reduces to with initial condition p 0 (0) = 1.By applying (2.7) iteratively we get which can be inverted by using formula (2.5) of [22], i.e.
For any ν ∈ (0, 1], it can be easily seen that it coincides with the result, obtained by a different approach in [2], by noting that The result of Theorem 2.1 shows that the first model proposed by Mainardi et al. [14] as a fractional version of the Poisson process (called renewal process of Mittag-Leffler type) coincides with the solution of equation (2.1) and therefore with (2.5).In the paper cited above the distribution is expressed in terms of successive derivatives of Mittag-Leffler functions as and is obtained by means of the fractional generalization of the Erlang density which represents the distribution of the waiting-time of the k-th event.Clearly, for ν = 1, f k (t) = t k−1 e −t (k−1)! is the Erlang density representing the waiting-time of the k-th event for the standard Poisson process (for λ = 1).Since we can rewrite (2.13) as , it is evident that it coincides with (2.5), for λ = 1.
We derive now an interesting relationship between the GML function in (2.5) and the solution to a fractional-diffusion equation, which is expressed in terms of the Wright function Let us denote by v 2ν = v 2ν (y, t) the solution to the Cauchy problem then it is well-known (see [17], p.142) that the solution of (2.15) can be written as  [18]), the telegraph-type fractional equations (in [17]) and even the higher-order heat-type equations with fractional time-derivative (see [1]).Formula (2.18) was obtained for the first time in [2] and we prove it here in an alternative form, by resorting to the Laplace transform.We compare (2.10) with the Laplace transform of the distribution p k (t), t > 0 of a standard Poisson process with intensity λ > 0, which reads (2.19) Formula (2.10) can be consequently written as which, by inversion, leads to the following convolution: where p k (t), t > 0 represents the distribution of the Poisson process with intensity λ = 1.The inverse transform in (2.21) can be expressed as follows by recalling that where g ν (•; y) is a stable law S ν (µ, β, σ) of order ν, with parameters µ = 0, β = 1 and σ = y λ cos πν We recognize in (2.24) the fractional integral of order ν of the stable law g ν and we apply the result (3.5) of [17], which can be rewritten, in this case, as where v 2ν (y, t) is given in (2.17).We can conclude from (2.5) that Remark 2.3 Formula (2.18) can be useful also in checking that the sum of p ν k (t) (either in the form (2.5) or (2.12)), for k ≥ 0, is equal to one, since it is Moreover, since result (2.25) holds for any t > 0, we can choose t = 1, so that we get, by a change of variable, This shows that the GML function E k+1 ν,νk+1 can be interpreted as the Laplace transform of the function In particular, for ν = 1 2 , since (2.18) reduces to where B λ (t) is a Brownian motion with variance 2λ 2 t (independent of N ), we get (for t = 1) (2.27) The previous relationship can be checked directly, as follows

Properties of the corresponding process
From the previous results we can conclude that the GML function E k+1 ν,νk+1 (−λt ν ), k ≥ 0, suitably normalized by the factor (λt ν ) k , represents a proper probability distribution and we can indicate it as Pr {N ν (t) = k} .
Moreover by (2.18) we can consider the process N ν (t), t > 0 as a time-changed Poisson process.It is well-known (see [8] and [11]) that, for a homogeneous Poisson process N subject to a random time change (by the random function Λ((0, t])), the following equality in distribution holds: where M (t), t > 0 is a Cox process directed by Λ.In our case the random measure Λ((0, t]) possesses distribution v 2ν given in (2.16)-(2.17)and we can conclude that N ν is a Cox process.This conclusion will be confirmed by the analysis of the factorial moments.Moreover, as remarked in [2] and [15], the fractional Poisson process N ν (t), t > 0 represents a renewal process with interarrival times U j distributed according to the following density, for j = 1, 2, ...: with Laplace transform Therefore the density of the waiting time of the k-th event, T k = k j=1 U j , possesses the following Laplace transform (2.31) Its inverse can be obtained by applying again (2.11) for β = ν, γ = νk and ω = −λ and can be expressed, as for the probability distribution, in terms of a GML function as Formula (2.32) coincides with (2.14), for λ = 1.The corresponding distribution function can be obtained in two different ways.The first one is based on (2.32) and yields We can check that (2.33) satisfies the following relationship Pr {T k < t} − Pr {T k+1 < t} = p ν k (t).
(2.34) Indeed from (2.33) we can rewrite (2.34) as The second method of evaluating the distribution function resorts to the probabilities given in the form (2.12): Finally if we rewrite (2.33) as and compare it with (2.35), we extract the following useful combinatorial relationship: Remark 2.4 As pointed out in [5] and [23], the density of the interarrival times (2.29) possess the following asymptotic behavior, for t → ∞: where the well-known expansion of the Mittag-Leffler function has been applied (see the Appendix for a proof of (2.37)).The density (2.36) is characterized by fat tails (with polynomial, instead of exponential, decay) and, as a consequence, the mean waiting time is infinite.For t → 0 the density of the interarrival times displays the following behavior: Remark 2.5 The distribution function (2.33) of the waiting-time of the k-th event coincides, for λ = 1, with the so-called Mittag-Leffler distribution of [13]- [19] (see formula (1) of [13]) when the time argument t takes integer values k ≥ 1.On the contrary the space argument x coincides, in our case, with the time t.Therefore the process X ν (t) with distribution ( 2.39) can be considered as the continuous-time analogue of the process representing the instant of the k-th event of the fractional Poisson process N ν .

Remark 2.6
We observe that also for the waiting-time density (2.32) we can find a link with the solution to the fractional diffusion equation (2.15).This can be shown by rewriting its Laplace transform (2.31) as By using again (2.23) we get Formula (2.40) permits us to conclude that f ν k (t) can be interpreted as a stable law S ν with a random scale parameter possessing an Erlang distribution.

The probability generating function
We consider now the equation governing the probability generating function, defined, for any 0 < u ≤ 1, as (2.41) From (2.1) it is straightforward that it coincides with the solution to the fractional differential equation subject to the initial condition G(u, 0) = 1.As already proved in [2] the Laplace transform of so that the probability generating function can be expressed as (2.44) By considering (2.44) together with the previous results we get the following relationship, valid for the infinite sum of GML functions: (2.45) For u = 1 it shows again that The result (2.45) can be checked by resorting to the Laplace transforms and noting that .
Formula (2.45) suggests a useful general relationship between the infinite sum of GML functions and the standard Mittag-Leffler function: (2.47) For u = 1 it shows again that ∞ k=0 p ν k (t) = 1.By considering the derivatives of the probability generating function (2.44) we can easily derive the factorial moments of N ν (t) which read

.48)
These are particularly useful in checking that N ν (t), t > 0 represents a Cox process with directing measure Λ.Indeed, as pointed out in [11], the factorial moments of a Cox process coincide with the ordinary moments of its directing measure.We show that this holds for N ν , by using the contour integral representation of the inverse of Gamma function, By applying formula (2.47) it is easy to obtain also the moments generating function of the distribution, defined as which is given by The r-th moments of the distribution can be obtained by successively deriving (2.49) as where C j,r are constants.For r = 1 we get , which coincides with the renewal function evaluated in [15], for λ = 1.It is evident also from (2.51) that the mean waiting time (which coincides with lim t→∞ t/m ν (t)) is infinite, since ν < 1.
3 Second-type fractional recursive differential equation

The solution
In this section we generalize the results obtained so far to a fractional recursive differential equation containing two time-fractional derivatives.We show that some properties of the first model of fractional Poisson process are still valid: the solutions represent, for k ≥ 0, a proper probability distribution and the corresponding process is again a renewal process.Moreover the density of the interarrival times display the same asymptotic behavior of the previous model.We consider the following recursive differential equation where ν ∈ (0, 1] , subject to the initial conditions Proof Following the lines of the proof of Theorem 2.1, we take the Laplace transform of equation (3.1) together with the conditions (3.2), thus obtaining the following recursive formula, for k ≥ 1 while, for k = 0, we get Therefore the Laplace transform of the solution reads We can invert (3.6) by using (2.11) with δ = 2k + 2, β = ν and γ = 2kν + 1 or γ = (2k + 1)ν + 1, thus obtaining the following expression We prove now the following general formula holding for a sum of GML functions: which can be proved by rewriting the l.h.s. as follows: x For m = 2k + 2, z = 1, x = λt ν and n = 2k formula (3.8) gives the following identity: ), which coincides with the first term in (3.3).
It remains to check only that the initial conditions in (3.2) hold: the first one is clearly satisfied since it is, for k = 0, = [for t = 0] = 1 and, for k ≥ 1, r!Γ(νr + 2kν + ν + 1) , which vanishes for t = 0.The second condition in (3.2) is immediately verified for k ≥ 1, since it is 9) which for t = 0 vanishes in the interval 1  2 < ν ≤ 1.Then we check that this happens also for k = 0: indeed in this case (3.9) reduces to Therefore it can be interpreted, for k = 0, 1, 2, ..., as the probability distribution Pr N ν (t) = k for a process N ν .Indeed, by (3.10), we get ) Moreover the relationship (3.10) shows that the process governed by the second-type equation can be seen as a first-type fractional process, which jumps upward at even-order events A 2k while the probability of the successive odd-indexed events A 2k+1 is added to that of A 2k .
We can check that expression (3.10) is the solution to equation (3.1), subject to the initial conditions (3.2), by using the form of p ν k appearing in the last line of (2.12) which is more suitable to this aim: The fractional derivatives of (3.13) can be evaluated by applying the definition (2.3), as follows: and, for ν ∈ 0, 1 2 and m = 1, It is easy to check that the last expression is valid also in the case ν ∈ 1 2 , 1 and m = 2.By considering (3.14) and (3.15) together we get where the second step follows by putting r = r − 2 in the first two sums and r = r − 1 in the second ones.
We want to show that (3.16) is equal to Then, by considering together (3.16) and (3.17), we get The sum appearing in (3.18) can be developed by considering that so that, by considering (3.19), we can rewrite (3.18) as follows If we consider now the first two terms of (3.20) we get

The probability generating function
As we did for the first model we evaluate the probability generating function and we show that it coincides with the solution to a fractional equation which arises in the study of the fractional telegraph process (see [17]).

Theorem 3.2 The probability generating function
, coincides with the solution to the following fractional differential equation subject to the initial condition G(u, 0) = 1 and the additional condition G t (u, 0) = 0 for 1/2 < ν < 1.
The explicit expression is given by

Properties of the corresponding process
We can prove that N ν (t), t > 0 represents a renewal process, by showing that, also for this model, the required relationship between p ν k (t) and distribution function of the waiting time T k of the k-th event holds: where or alternatively for the Laplace transform of (3.25) In view of relationship (3.10), we can infer that each interarrival time U j is distributed as the sum of two independent interarrival times U j of the first model and therefore from (2.30) we have that = L f ν 1 (t); s .
By recursively applying (3.26), starting with (3.27), we arrive at By applying again (2.11) for β = ν, γ = 2νk and ω = −λ we invert (3.28) and obtain Moreover from (3.27) or (3.29) we easily get Therefore, in this case, both the waiting-time of the k-th event and the interarrival times possess distributions which are expressed in terms of GML functions.
Remark 3.3 As a check, we derive directly the probability density of the U j from that of U j given in (2.29), without the use of Laplace transform as follows: This confirms that one out of two Poisson events are disregarded in this case (as described in Remark 3.1).

Remark 3.4
We evaluate now the asymptotic behavior of the interarrival-times density (3.30), as follows: By applying (2.36) and (2.37), we finally get If we compare (3.32) with the analogous result (2.36) obtained for the first model, we can conclude that the interarrival-times density displays the same asymptotic behavior, with the power law decay (3.32), so that again the mean waiting time is infinite.
For t → 0, instead of (2.38), we get in this case The behavior near the origin of the density of the interarrival times U j has a different structure for ν < 1/2 (tends to infinity) and ν ∈ (1/2, 1] (vanishes for t → 0 + ).For ν = 1 2 we have instead that Pr U j ∈ dt is constant for t = 0 and equal to λ 2 .
We can conclude that N ν is a renewal process and the corresponding renewal function is given by Formula (3.33) can be obtained either by deriving the probability generating function (3.23) (for u = 1) or by using the well-known relationship between the Laplace transforms of the renewal function and the interarrival-times density: Indeed, in this case, it is which gives (3.33), by applying (2.11) for β = ν, δ = 1 and γ = 2ν + 1.
Remark 3.5 By comparing the second form of (3.33) with (2.51), we can note that the following relationship between the renewal functions of the two models holds: This can be alternatively proved by applying (3.10) as follows: The last term in (3.35), which coincides with the sum of the probabilities of an odd number of events of the first model, can be evaluated as follows: as can be checked by resorting to the Laplace transform: Formula (3.34) confirms that the mean waiting time is infinite (as we have noticed in Remark 3.4): indeed it is where the second limit can be evaluated by the following considerations: (see (5.2) in the Appendix).

The special case ν = 1
We consider now the previous results in the special case ν = 1.Equation (3.1) reduces in this case to the second-order equation: and the corresponding solution (3.3) is given by It is easy to check can show that (3.37) solves (3.36) with initial condition so that we get Therefore, in this case, the random variable T j , representing the instant of the j-th event, is distributed as Gamma(λ, 2j).We derive equation (3.37) in an alternative way, which is similar to the construction of the standard Poisson process, by considering the relationships (3.11) or, equivalently, the property of the interarrival times described in Remark 3.3.We can write that where T j = inf {t > 0 : N (t) = j} which represents the time of the j-th event of N (t). .Theorem 3.3 The probability distribution of the process N (t), t > 0 described by (3.40) and (3.43) solves equation (3.36).Proof Let us consider now the following intervals We evaluate the following probability, by stopping the approximation at the second-order terms: where we used the following well-known approximations valid for the standard Poisson process: Pr (1 Poisson event in ∆t) = λ∆t e −λ∆t λ∆t − λ 2 (∆t) Pr (2 Poisson events in ∆t) = λ 2 (∆t) By ignoring the terms of order greater than (∆t) 2 the probability (3.44) can be rewritten as follows

Since we can write
Pr and, analogously Pr We consider now the following probability, on a single interval of length ∆t We multiply (3.47) by 2(1 − λ∆t) and ignore the terms of order greater than (∆t) 2 , so that we get which can be substituted in the first line of (3.46).Finally we obtain We divide (3.48) by (∆t) 2 , so that we get By letting ∆t → 0 we easily obtain equation (3.36).
Remark 3.6 This special case is particularly interesting because it describes a generalization of the Poisson process which have been used in many articles.Random motions at finite velocities spaced by this particular renewal process have been considered by different authors ([5]-[6]- [10]).In particular in [3] it has been studied a model with uniformly distributed deviations which take place at even-order Poisson events and therefore its interarrival times are distributed as Gamma(λ, 2).

Conclusions
The results of the previous sections can be generalized to the n-th order case, if we consider the following equation where ν ∈ (0, 1) , subject to the initial conditions and p −1 (t) = 0. Following the same steps as for the first two models we get the Laplace transform of (4.1): which can be solved, in view of (2.6) and the initial conditions (4.2), recursively, yielding For n = 2, we obtain from (4.3) formula (3.6).The Laplace transform can be inverted by applying again (2.11) and the solution is given, also in this case, as a sum of GML functions as follows: For n = 1, we get the distribution of the first model (2.5), while, for n = 2, we get the solution of the second-type equation in the form (3.7).In this case the use of (3.8) requires much harder calculations.Nevertheless a relationship similar to (3.10) can be obtained, even for n > 2, by studying the density f ν k (t), t > 0 of the waiting time of the k-th event T k .As already seen in section 3, the following identity must be satisfied by f ν k (t) : so that, by substituting (4.3) in the l.h.s. of (4.5) we get which can be inverted as usual, thus obtaining f ν k (t) = λ nk t nνk−1 E nk ν,nνk (−λt ν ).(4.7) Again the process is a renewal one, since (4.7) coincides with the sum of k independent and identically distributed random variables U j 's (representing the interarrival times) with density given by Pr U j ∈ dt /dt = L −1 λ n (s ν + λ) n ; t = λ n t nν−1 E n ν,nν (−λt ν ) = f nν 1 (t).(4.8) Formula (4.8) shows that each interarrival time of the n-th order case is distributed as the sum of n independent interarrival times of the first model.This suggests that the following relationship between the corresponding probability distributions holds: p ν k (t) = p ν nk (t) + p ν nk+1 (t) + ... + p ν nk+n−1 (t), n > 2. (4.9)

Appendix
We derive the integral form of the Mittag-Leffler function (used in (2.37)), and some generalizations.We start by showing that, for 0 < ν < 1, E ν,1 (−t ν ) = sin (νπ) π e −r 1/ν t dr where X is distributed as a Cauchy with parameters − cos(πν) and sin (νπ) .Formula (5.3) permits us to study the particular case where ν = 1, since we can write (5.3), by means of the characteristic function of a Cauchy r.v., as follows: Analogously we can study the case where ν = 1/2: from the first line of (5.4) and by applying formula 3.896.4,p.480 of [7], we get which coincides with the form given in [17], for x = − √ t.We generalize formula (5.1) to the case of a two-parameters Mittag-Leffler function: for 0 < ν < 1 and 0 < β < ν + 1, we have that