Thomas M Michelitsch

Federico Polito
email: federico.polito@unito.it

Alejandro P Riascos

ON DISCRETE TIME GENERALIZED FRACTIONAL POISSON PROCESSES AND RELATED STOCHASTIC DYNAMICS

Recently the so-called Prabhakar generalization of the fractional Poisson counting process attracted much interest for his flexibility to adapt real world situations. In this renewal process the waiting times between events are IID continuous random variables. In the present paper we analyze discrete-time counterparts: Renewal processes with integer IID interarrival times which converge in well-scaled continuous-time limits to the Prabhakar-generalized fractional Poisson process. These processes exhibit non-Markovian features and long-time memory effects. We recover for special choices of parameters the discrete-time versions of classical cases, such as the fractional Bernoulli process and the standard Bernoulli process as discrete-time approximations of the fractional Poisson and the standard Poisson process, respectively. We derive difference equations of generalized fractional type that govern these discrete time-processes where in well-scaled continuous-time limits known evolution equations of generalized fractional Prabhakar type are recovered. We also develop in Montroll-Weiss fashion the "Prabhakar Discrete-time random walk (DTRW)" as a random walk on a graph time-changed with a discrete-time version of Prabhakar renewal process. We derive the generalized fractional discrete-time Kolmogorov-Feller difference equations governing the resulting stochastic motion. Prabhakar-discrete-time processes open a promising field capturing several aspects in the dynamics of complex systems.

INTRODUCTION

Classically, random occurrence of events in time is modeled as standard Poisson processes, i.e. with exponentially distributed waiting-times and the memoryless Markovian property. Applying this model to random motions in continuous spaces or graphs lead to continuoustime Markov chains. However, one has recognized that many phenomena in anomalous transport and diffusion including the dynamics in certain complex systems exhibit power law distributed waiting-times with non-Markovian long-time memory features that are not compatible with classical exponential patterns [START_REF] Kutner | The continuous time random walk, still trendy: fifty-year history, state of art and outlook[END_REF][START_REF] Gorenflo | The asymptotic universality of the Mittag-Leffler waiting time law in continuous time random walks[END_REF][START_REF] Gorenflo | Mittag-Leffler Waiting Time, Power Laws, Rarefaction, Continuous Time Random Walk, Diffusion Limit[END_REF][START_REF] Metzler | The Random Walk's Guide to Anomalous Diffusion : A Fractional Dynamics Approach[END_REF][START_REF] Metzler | The restaurant at the end of the random walk: recent developments in the description of anomalous transport by fractional dynamics[END_REF].

A powerful approach to tackle these phenomena is obtained by admitting fat-tailed powerlaw waiting-time densities where Mittag-Leffler functions come into play as natural generalizations of the classical exponentials. A prototypical example is the fractional Poisson process (FPP), a counting process with unit-size jumps and IID Mittag-Leffler distributed waiting-times. Fractional diffusion equations governing the fractional Poisson process and many other related properties are already present in the specialized literature -see e.g. [START_REF] Hilfer | Fractional master equations and fractal time random walks[END_REF][START_REF] Repin | Fractional Poisson law[END_REF][START_REF] Laskin | Fractional Poisson process[END_REF][START_REF] Laskin | Some applications of the fractional Poisson probability distribution[END_REF][START_REF] Mainardi | A fractional generalization of the Poisson processes[END_REF][START_REF] Gorenflo | On the Fractional Poisson Process and the Discretized Stable Subordinator[END_REF][START_REF] Zaslavsky | Chaos, fractional kinetics, and anomalous transport[END_REF][START_REF] Beghin | Fractional Poisson processes and related planar random motions[END_REF][START_REF] Meerschaert | The Fractional Poisson Process and the Inverse Stable Subordinator[END_REF][START_REF] Meerschaert | Stochastic Models for Fractional Calculus[END_REF].

Broadly speaking, Markov chains permitting arbitrary waiting times define so called semi-Markov processes [START_REF] Pachon | On Discrete-Time Semi-Markov processes[END_REF]. This area was introduced independently by Lévy [START_REF] Lévy | Processus semi-Markovien[END_REF] Smith [START_REF] Smith | Regenerative Stochastic Processes[END_REF] and Takács [START_REF] Takács | On sojourn time problem[END_REF] and fundamentals of this theory were derived by Pyke [START_REF] Pyke | Markov Renewal Processes with Finitely many states[END_REF] and Feller [START_REF] Feller | On Semi-Markov Processes[END_REF] among others. In these classical models, semi-Markov processes have as special cases continuous-time renewal processes, i.e. the waiting times are IID absolutely continuous random variables. On the other hand, many applications require intrinsic discrete-time scales, and thus semi-Markov processes where the waiting times are discrete integer random variables open an interesting field which merits deeper analysis. Discrete-time renewal processes are relatively little touched in the literature compared to their continuous-time counterparts. A discrete variant of the above-mentioned Mittag-Leffler distribution was derived by Pillai and Jayakumar [START_REF] Pillai | Discrete Mittag-Leffler distributions[END_REF]. An application in terms of a discrete-time random walk (DTRW) diffusive transport model is developed recently [START_REF] Angstmann | A Time-Fractional Generalised Advection Equation from a Stochastic Process[END_REF]. A general approach for discrete-time semi-Markov process and time-fractional difference equations was developed in a recent contribution by Pachon, Polito and Ricciuti [START_REF] Pachon | On Discrete-Time Semi-Markov processes[END_REF]. The aim of the present paper is to develop new pertinent discrete-time counting processes that contain for certain parameter choices classical counterparts such as fractional Bernoulli and standard Bernoulli, as well is in well-scaled continuous-time limits their classical continuous-time counterparts such as fractional Poisson and standard Poisson. A further goal of this paper is to analyze the resulting stochastic dynamics on graphs.

The present paper is organized as follows. As a point of departure, we introduce a class of discrete-time renewal processes which represent approximations of the continuous-time Prabhakar process. The Prabhakar renewal process was first introduced by Cahoy and Polito [START_REF] Cahoy | Renewal processes based on generalized Mittag-Leffler waiting times[END_REF] and the continuous-time random walk (CTRW) model based on this process was developed by Michelitsch and Riascos [START_REF] Michelitsch | Continuous time random walk and diffusion with generalized fractional Poisson process[END_REF]. We describe the Prabhakar renewal process in Section 2. For a thorough review of properties and definitions of Prabhakar-related fractional calculus we refer to the recent review article of Giusti et al [START_REF] Giusti | General Fractional Calculus and Prabhakar Theory[END_REF].

Section 3 is devoted to derive discrete-time versions of the Prabhakar renewal process using a composition of two 'simple' processes. Then, in Section 3.1 we show that under suitable scaling conditions the continuous-time Prabhakar process is recovered.

In Section 3.2 we develop a general procedure to generate discrete-time approximations of the Prabhakar process. We construct these processes in such a way that the waiting time distributions are vanishing at t = 0. The choice of this condition turns out to be crucial to obtain state-probabilities allowing to define proper Cauchy initial value problems in stochastic motions.

As a prototypical example, we analyze in Section 4 the most simple version of discrete-time Prabhakar process with above mentioned good initial conditions. We call this version of Prabhakar discrete-time counting process the 'Prabhakar Discrete-Time Process' (PDTP).

We derive the state-probabilities (probabilities for n arrivals in a given time interval), i.e. the discrete-time counterpart of the Prabhakar-generalized fractional Poisson distribution which was deduced in the references [START_REF] Cahoy | Renewal processes based on generalized Mittag-Leffler waiting times[END_REF][START_REF] Michelitsch | Continuous time random walk and diffusion with generalized fractional Poisson process[END_REF]. The PDTP is defined as a generalization of the 'fractional Bernoulli process' introduced in [START_REF] Pachon | On Discrete-Time Semi-Markov processes[END_REF] which is contained for a certain choice of parameters as well as the standard Bernoulli counting process. We prove these connections by means of explicit formulas. We show explicitly that the discrete-time waiting time and state distributions of a PDTP converge in well-scaled continuous-time limits to their known continuous-time Prabhakar function type counterparts. These results contain for a certain choice of parameters the well-known classical cases of fractional Poisson and standard Poisson distributions, respectively. We show that the well-scaled continuous-time limits yield the state probabilities of Laskin's fractional Poisson [START_REF] Laskin | Fractional Poisson process[END_REF] and standard Poisson distributions, respectively.

In Section 4.2 we derive for the PDTP the discrete-time versions of the generalized fractional Kolmogorov-Feller equations that are solved by the PDTP state-probabilities. These equations constitute discrete-time convolutions of generalized fractional type reflecting long-time memory effects and non-Markovian features (unless in the classical standard Bernoulli with Poisson continuous-time limit case). We show that discrete-time fractional Bernoulli and standard Bernoulli processes are contained for certain choice parameters and that the same is true for their continuous-time limits: They recover the classical Kolmogorov-Feller equations of fractional Poisson and standard Poisson, respectively. Section 4.3 is devoted to the analysis of the expected number of arrivals and their asymptotic features. This part is motivated by the important role of this quantity for a wide class of diffusion problems and stochastic motions in networks and lattices.

As an application we develop in Section 5 in Montroll-Weiss fashion the 'Discrete-Time-Random Walk' (DTRW) on undirected networks and analyze a normal random walk subordinated to the PDTP. We call this walk the 'Prabhakar DTRW'. The developed DTRW approach is a general model to subordinate random walks on graphs to discrete-time counting processes. Although we focus on undirected graphs the DTRW approach can be extended to general walks such as on directed graphs or strictly increasing walks on the integer line. Such an example is briefly outlined in Appendix A.5, namely a strictly increasing walk subordinated to the Sibuya counting process.

Further we derive for the Prabhakar DTRW discrete-time Kolmogorov-Feller generalized fractional difference equations that govern the resulting stochastic motion on undirected graphs and demonstrate by explicit formulas the contained classical cases of fractional Bernoulli and standard Bernoulli and their fractional Poisson and Poisson continuous-time limits, respectively. The applications of this section are motivated by the huge upswing of network science which has become a rapidly growing interdisciplinary field [START_REF] Noh | Random walks on complex networks[END_REF][START_REF] Newman | Networks: An Introduction[END_REF][START_REF] Riascos | Aging in transport processes on networks with stochastic cumulative damage[END_REF][START_REF] Michelitsch | Fractional Dynamics on Networks and Lattices[END_REF][START_REF] Riascos | Fractional dynamics on networks: Emergence of anomalous diffusion and Lévy flights[END_REF][START_REF] Riascos | Random walks on networks with stochastic resetting[END_REF][START_REF] Riascos | Networks and long-range mobility in cities: A study of more than one billion taxi trips in New York City[END_REF] (and see the references therein).

Section 6 is devoted to analyze the influence of the initial condition in the discrete-time waiting time density on DTRW features. We explore resulting effects introducing the feature of uncertainty and randomness in the DTRW initial conditions. All proofs in the paper are accompanied by detailed derivations and supplementary materials in the Appendices.

PRABHAKAR CONTINUOUS-TIME RENEWAL PROCESS

Among several generalizations of the fractional Poisson process which were proposed in the literature, the so called Prabhakar type generalization which we refer to as 'Generalized Fractional Poisson Process (GFPP)' or also 'Prabhakar process' seems to be one of the most pertinent candidates. The GFPP was first introduced by Cahoy and Polito [START_REF] Cahoy | Renewal processes based on generalized Mittag-Leffler waiting times[END_REF] and applied to stochastic motions in networks and lattices by Michelitsch and Riascos [START_REF] Michelitsch | Continuous time random walk and diffusion with generalized fractional Poisson process[END_REF][START_REF] Michelitsch | Generalized fractional Poisson process and related stochastic dynamics[END_REF][START_REF] Michelitsch | Generalized space-time fractional dynamics in networks and lattices Generalized Space-Time Fractional Dynamics in Networks and Lattices[END_REF].

The Prabhakar function which is a three-parameter generalization of the Mittag-Leffler function was introduced in 1971 by Prabhakar [START_REF] Prabhakar | A singular integral equation with a generalized Mittag-Leffler function in the kernel[END_REF] and has attracted recently much attention due to its great flexibility to adapt real-world situations. Meanwhile, the Prabhakar function has been identified as a matter of great interest worthy of thorough investigation. For a comprehensive review of properties and physical applications with generalized fractional calculus emerging from Prabhakar functions we refer to the recent review article by Giusti et al. [START_REF] Giusti | A practical guide to Prabhakar fractional calculus[END_REF] and consult also [START_REF] Garra | Tomovski. Hilfer-Prabhakar derivatives and some applications[END_REF].

The interesting feature of the Prabhakar process is that it contains the fractional Poisson process as well as the Erlang-and standard Poisson processes as special cases. The related Prabhakar-generalized fractional derivative operators may be considered as among the most sophisticated tools to cover certain aspects of complexity in physical systems [START_REF] Santos | Fractional Prabhakar Derivative in Diffusion Equation with Non-Static Stochastic Resetting[END_REF][START_REF] Mainardi | On complete monotonicity of the Prabhakar function and non-Debye relaxation in dielectrics[END_REF]. The continuous-time Prabhakar renewal process is characterized by waiting time density with Laplace transform [START_REF] Cahoy | Renewal processes based on generalized Mittag-Leffler waiting times[END_REF] χα,ν (s) = ξ ν 0 (ξ 0 + s α ) ν , ξ 0 > 0 , 0 < α ≤ 1 , ν > 0 .

(

) 1
Laplace inversion yields the waiting-time PDF of the GFPP [START_REF] Cahoy | Renewal processes based on generalized Mittag-Leffler waiting times[END_REF][START_REF] Michelitsch | Continuous time random walk and diffusion with generalized fractional Poisson process[END_REF] χ α,ν (t) = ξ ν 0 t αν-1

∞ m=0 (ν) m m! (-ξ 0 t α ) m Γ (αm + να) = ξ ν 0 t να-1 E ν α,να (-ξ 0 t α ), t ∈ R + , ( 2
)
which we refer to as Prabhakar-Mittag-Leffler density. The choice of this name is since this expression appears as a generalization of the Mittag-Leffler density (and recovering the Mittag-Leffler density for ν = 1). Expression (2) contains the Prabhakar-Mittag-Leffler function (also referred to as Prabhakar function) [START_REF] Prabhakar | A singular integral equation with a generalized Mittag-Leffler function in the kernel[END_REF] defined by

E c a,b (z) = ∞ m=0 (c) m m! z m Γ (am + b) , ℜ{a} > 0 , ℜ{b} > 0 , c, z ∈ C, (3)
where (c) m indicates the Pochhammer-symbol

(c) m = Γ (c + m) Γ (c) =      1, m = 0, c(c + 1) . . . (c + m -1), m = 1, 2, . . . ( 4
)
Further aspects on generalizations of Mittag-Leffler functions such as the Prabhakar function are outlined in [START_REF] Giusti | General Fractional Calculus and Prabhakar Theory[END_REF][START_REF] Shukla | On a generalization of Mittag-Leffler function and its properties[END_REF], and for an analysis of properties and applications we refer to the references [START_REF] Giusti | A practical guide to Prabhakar fractional calculus[END_REF][START_REF] Garra | Tomovski. Hilfer-Prabhakar derivatives and some applications[END_REF][START_REF] Haubold | Mittag-Leffler functions and their applications[END_REF][START_REF] Beghin | Poisson-type processes governed by fractional and higher-order recursive differential equations[END_REF][START_REF] Mathai | Some properties of Mittag-Leffler functions and matrix variant analogues: A statistical perspective[END_REF]. The GFPP recovers for ν = 1 with 0 < α < 1 the Laskin fractional Poisson process [START_REF] Laskin | Fractional Poisson process[END_REF], for ν > 0 with α = 1 the (generalized) Erlang process, and for α = 1, ν = 1 the standard Poisson process and their related distributions. For details and derivations consult [START_REF] Cahoy | Renewal processes based on generalized Mittag-Leffler waiting times[END_REF][START_REF] Michelitsch | Continuous time random walk and diffusion with generalized fractional Poisson process[END_REF][START_REF] Michelitsch | Generalized fractional Poisson process and related stochastic dynamics[END_REF][START_REF] Michelitsch | Generalized space-time fractional dynamics in networks and lattices Generalized Space-Time Fractional Dynamics in Networks and Lattices[END_REF].

DISCRETE-TIME VARIANTS OF THE GFPP

This section is devoted to the construction of discrete-time variants of the Prabhakar renewal process by means of a composition of two 'simple' discrete-time processes. To this end we evoke first of all the concept of 'discrete-time renewal process' where also term 'discrete-time renewal chain' is used in the literature [START_REF] Pachon | On Discrete-Time Semi-Markov processes[END_REF][START_REF] Barbu | Semi-Markov Chains and Hidden Semi-Markov Models toward Applications[END_REF].

We introduce the strictly increasing random walk X = (X n ) n≥1 , such that

X n = n j=1 Z j , Z j ∈ N, X 0 = 0 , ( 5
)
where the steps are non-zero IID integer random variables Z j = k ∈ N a.s., following each the same distribution P(Z j = k) = w(k). With the choice of w(0) = 0 the walk (5) becomes strictly increasing. A random walk X defined in ( 5) is the natural discrete-time counterpart to a (strictly increasing) subordinator [START_REF] Meerschaert | The Fractional Poisson Process and the Inverse Stable Subordinator[END_REF][START_REF] Pachon | On Discrete-Time Semi-Markov processes[END_REF] (and see the references therein).

In a discrete-time renewal process the random integers X n of ( 5) indicate the times when events occur; we refer them to as 'arrival times' or 'renewal times' where n ∈ N 0 counts the events. We also use the terms 'renewals' and 'arrivals'. The integer IID times Z j ∈ N between the events follow then the same waiting-time distribution P(Z j = k) = w(k). Let us now introduce the generating function of the waiting-time distribution

P(Z = k) = w(k) as Eu Z = w(u) = ∞ k=0 u k w(k) = w(1 )u + w(2 )u 2 + . . . , |u| ≤ 1 , (6)
where w(u)| u=1 = 1 reflects normalization of the w-distribution. Generally generating functions are highly elegant and powerful tools which we will use extensively in the present paper. For some definitions and properties we refer to Appendix A.1.

Consider now two discrete-time renewal processes, I and II, having waiting-time distributions w I , w II , defined in the above general way having both zero initial conditions w I (0) = w II (0) = 0. Then we generate a new discrete-time renewal process resulting from a composition of these two 'elementary' processes. Specifically, its waiting-time distribution W(k)

has generating function such that

W(u) = ∞ n=1 w II (n)(w I (u)) n = wII (w I (u)), (7)
with

(w I (u)) n = E w I u X = ∞ k 1 =1 ∞ k 2 =1 . . . ∞ kn=1 w I (k 1 )w I (k 2 ) . . . w I (k n )u n j=1 k j , X 0 = 0 , (8)
where X = (X n ) n≥1 is the partial sum [START_REF] Metzler | The restaurant at the end of the random walk: recent developments in the description of anomalous transport by fractional dynamics[END_REF] in which the random jumps are w I -distributed and is characterized by the generating function [START_REF] Laskin | Fractional Poisson process[END_REF]. The event counter n is then considered random in [START_REF] Repin | Fractional Poisson law[END_REF] with distribution w II . We observe that W(u = 1) = 1 reflects normalization of the new W-distribution which furthermore fulfills the desired initial condition W(u = 0) = W(t = 0) = 0. This new waiting-time distribution then is characterized by the probabilities

P(Z = t) = W(t) = 1 t! d t du t W(u) u=0 = ∞ n=1 w II (n)(w I ⋆) n (t), t ∈ N 0 , (9)
where (w I ⋆) n (t) stands for convolution power 1 . Now, let us assume that process I is a Sibuya counting process with waiting-times following 'Sibuya(α)' (See Appendix A.4 for definitions and some properties). For the process II we choose the waiting time distribution

w (ν) B (k) k=0 = 0, w (ν) B (k) = p ν q k-1 (-1) k-1 -ν k -1 = (ν) k-1 (k -1)! p ν q k-1 , k ∈ N, ν > 0 , p + q = 1 , ( 10
)
where for ν = 1 (10) yields the geometric waiting-time distribution P(Z = k) = pq k-1 of the Bernoulli process [START_REF] Pachon | On Discrete-Time Semi-Markov processes[END_REF] 2 with w (ν) B (0) = 0. We further employed here the Pochhammer symbol (ν) m defined in [START_REF] Metzler | The Random Walk's Guide to Anomalous Diffusion : A Fractional Dynamics Approach[END_REF]. The waiting-time distribution [START_REF] Mainardi | A fractional generalization of the Poisson processes[END_REF] has generating function

w(ν) B (u) = ∞ k=1 p ν q k-1 u k -ν k -1 = up ν (1 -qu) ν = uξ ν (ξ + 1 -u) ν (11)
where we have put ξ = p q p = ξ ξ+1 , q = 1 ξ+1 and w(ν) B (u) u=1 = 1 reflects normalization. For ν = 1, (11) recovers the generating function of the standard Bernoulli counting process. Now we generate a new process in the above described fashion. The new discrete-time process hence with [START_REF] Repin | Fractional Poisson law[END_REF] has waiting-time generating function

E α,ν u Z = ψ(ν) α (u) = w(ν) B ( wα (u)) = w(ν) B (1 -(1 -u) α ) = (1 -(1 -u) α ) φ(ν) α (u), φ(ν) α (u) = ( p q ) ν p q + (1 -u) α ν , ν > 0, 0 < α ≤ 1. (12)
For ν = 1 and 0 < α < 1, formula [START_REF] Zaslavsky | Chaos, fractional kinetics, and anomalous transport[END_REF] recovers the generating function of a 'discrete-time Mittag-Leffler distribution' (of so-called 'type A') DML A where this process has been named

'fractional Bernoulli process (type A)' in [START_REF] Pachon | On Discrete-Time Semi-Markov processes[END_REF]. The waiting time distribution ϕ ν=1 α (t) defines a discrete-time approximation of the Mittag-Leffler waiting-time distribution [START_REF] Pillai | Discrete Mittag-Leffler distributions[END_REF]. Indeed [START_REF] Zaslavsky | Chaos, fractional kinetics, and anomalous transport[END_REF] is generating function of a discrete-time waiting time distribution which is for ν > 0 and 0 < α < 1 a generalization of discrete-time fractional Bernoulli process (of 'type A') 3 and recovers for ν = 1, α = 1 the generating function of the standard Bernoulli counting process.

Our goal now is to show that the renewal process defined by generating function [START_REF] Zaslavsky | Chaos, fractional kinetics, and anomalous transport[END_REF] is a discrete-time version of the Prabhakar process. To this end we expand [START_REF] Zaslavsky | Chaos, fractional kinetics, and anomalous transport[END_REF] as follows

φ(ν) α (u) = ξ ν (1 -u) -αν (1 + ξ(1 -u) -α ) -ν = ∞ m=0 -ν m ξ m+ν (1 -u) -α(m+ν) , ξ = p q ψ(ν) α (u) = (1 -(1 -u) α ) φ(ν) α (u) = ∞ m=0 (-1) m (ν) m m! ξ m+ν (1 -u) -(m+ν)α -(1 -u) -(m+ν-1)α (13)
1 See Appendix A.1 for details.

2 See Definition 3.1 with Eqs. ( 53), [START_REF] Michelitsch | On Recurrence and Transience of Fractional RandomWalks in Lattices[END_REF] in that paper. 3 See [START_REF] Pachon | On Discrete-Time Semi-Markov processes[END_REF] for details of classification scheme 'type A and B discrete-time processes'. α (t, ξ) versus ξ for different values of t indicated in the colorbar t = 1, 2, . . . , 10 calculated using expression [START_REF] Meerschaert | The Fractional Poisson Process and the Inverse Stable Subordinator[END_REF]. For ξ → 0 they approach a power-law ∼ ξ ν indicated by dashed lines (See also in ( 14)). converging 4 for |(1u) -1 ξ 1 α | < 1 where (ν) m indicates above introduced Pochhammer symbol (4). We notice that in [START_REF] Beghin | Fractional Poisson processes and related planar random motions[END_REF] the powers (1u) -µ can be seen as the generating functions of expected numbers of Sibuya hits (See Appendix A.4). Generating function [START_REF] Beghin | Fractional Poisson processes and related planar random motions[END_REF] yields the probabilities

ϕ (ν) α (t, ξ) = 1 t! d t du t φ(ν) α (u)| u=0 = ∞ m=0 (-1) m (ν) m m! ξ m+ν t! ((m + ν)α) t P α,ν (Z = t) = ψ (ν) α (t, ξ) = 1 t! d t du t ψ(ν) α (u)| u=0 = ∞ m=0 (-1) m (ν) m m! ξ m+ν t! ((m + ν)α) t -((m + ν -1)α) t t ∈ N 0 ( 14
)
where we used that

1 t! d t du t (1 -u) -β | u=0 = Γ(β+t) Γ(β)t! = (β)t
t! . 4 For u = e -hs and ξ(h) = ξ0h α we see that lim h→0 (1e -sh ) -1 hξ

1 α 0 = s -1 ξ 1 α
0 thus the Laplace variable has to fulfill |s| > ξ

1 α 0 .
The discrete-time densities of ( 14) are plotted in Figure [START_REF] Kutner | The continuous time random walk, still trendy: fifty-year history, state of art and outlook[END_REF] for different values of t as functions of the parameter ξ which defines a time scale in the process. The densities behave like a power law similar to ϕ (ν) α (0, ξ) = ξ ν (1+ξ) ν ∼ ξ ν for ξ and t small; the power-law is depicted with dashed lines in Figure 1.

Now we show that both of the distributions in [START_REF] Meerschaert | The Fractional Poisson Process and the Inverse Stable Subordinator[END_REF] are approximations of the Prabhakar-Mittag-Leffler density (2), but only ψ (ν) α (t) per construction fulfills the desired initial condition ψ (ν) α (t) t=0 = 0.

CONTINUOUS-TIME LIMIT

We recommend to consult Appendix A.2 where we outline properties of the shift operator T(..) which we are extensively using to define 'well-scaled' continuous-time limit procedures.

Further we mention that throughout the analysis to follow we utilize as synonymous notations lim x → a ± 0 and lim x → a± for left-and right sided limits, respectively.

Let us introduce the (scaled) 'discrete-time waiting time density' (See (164)-( 168))

χ α,ν (t) h = ψ(ν) α ( T-h )δ h (t) = ∞ k=1 ψ (ν) α (k, ξ 0 h α )δ h (t -kh) = 1 h ψ (ν) α t h , ξ 0 h α , t ∈ hN 0
(15) generalizing the notion of (continuous-time) waiting-time density to the discrete-time cases.

We employ in this paper the notation (..)(t) h for scaled quantities defined on t ∈ hZ 0 and skip subscript h when h = 1

(Appendix A.2). Note that ψ(ν) α ( T-h ) = ∞ k=0 ψ (ν) α (k, ξ)
T-hk is the operator function obtained by replacing u → T-h in the generating function of (13) 5 . In (15) occur the probabilities [START_REF] Meerschaert | The Fractional Poisson Process and the Inverse Stable Subordinator[END_REF] and the discrete-time δdistribution δ h (τ ) (τ ∈ hZ 0 ) is defined in (164) in Appendix A.2. Note that the multiplier h -1 on the right-hand side of (15) comes into play due to the definition (164) of the discrete-time δ-distribution δ h (t) guaranteeing the discrete-time densities indeed have physical dimension sec -1 . We can then write for [START_REF] Meerschaert | Stochastic Models for Fractional Calculus[END_REF] the distributional relations

P(Z = k) = ψ (ν) α (k, ξ = ξ 0 h α ) (k ∈ N 0 ) of
ψ(ν) α ( T-h )δ h (t) = 1 -(1 -e -hDt ) α ξ(h) ν (ξ(h) + (1 -e -hDt ) α ) ν δ h (t), ξ(h) = ξ 0 h α lim h→0 ψ(ν) α ( T-h )δ h (t) = ξ ν 0 (D α t + ξ 0 ) ν δ(t) = ξ ν 0 D -αν t (1 + ξ 0 D -α t ) ν δ(t), t ∈ R (16)
where in the limiting process only fractional integral operators D -β t (β > 0) occur. Indeed the discrete-time density ( 16) can be conceived as a 'generalized fractional integral' (See Appendix A.3, with relations (164) -(166) and consult also [START_REF] Pachon | On Discrete-Time Semi-Markov processes[END_REF][START_REF] Michelitsch | Sur une généralisation de l'opérateur fractionnaire[END_REF]). The scaling of the constant ξ(h) is chosen such that the limit h → 0 of ( 16) exists. Clearly the continuous-time limit in [START_REF] Pachon | On Discrete-Time Semi-Markov processes[END_REF] exists if and only if lim h→0 ξ(h)(1e -hDt ) -α exists and hence ξ(h) = ξ 0 h α is the required scaling where ξ 0 > 0 is an arbitrary positive dimensional constant of physical dimension sec -α and independent of h. We notice that [START_REF] Pachon | On Discrete-Time Semi-Markov processes[END_REF] indeed is a distributional representation of Prabhakar-Mittag-Leffler density (2) which follows in view of Laplace transform of ( 16), namely χα,ν (s) = ξ ν 0 (ξ 0 +s α ) ν coinciding with Laplace transform (1) of the Prabhakar 5 Bear in mind properties of the shift operator such as

( T-h ) a f (t) = T-ah f (t) = f (t -ah) , a ∈ R and T-h = e -hDt (Dt = d dt ).
density. To obtain the limiting density explicitly we introduce the rescaled variable τ k = hk kept finite for h → 0. Hence in [START_REF] Meerschaert | The Fractional Poisson Process and the Inverse Stable Subordinator[END_REF] k = τ k h ∈ N becomes very large for h → 0 thus we can use the asymptotic expression

(β) k k! = Γ(β+k) Γ(k+1)Γ(β) ∼ k β-1 Γ(β)
holding for k large. Then consider first the scaling behavior of the coefficients in [START_REF] Meerschaert | The Fractional Poisson Process and the Inverse Stable Subordinator[END_REF], namely

((m + ν)α) k k! (ξ(h)) m+ν ∼ (ξ 0 h α ) m+ν k α(m+ν)-1 Γ(α(m + ν)) = hξ m+ν 0 τ α(m+ν)-1 k Γ(α(m + ν)) , ((m + ν -1)α) k k! (ξ(h)) m+ν ∼ h 1+α ξ m+ν 0 τ α(m+ν-1)-1 k Γ(α(m + ν -1)) . ( 17
)
It follows from [START_REF] Meerschaert | Stochastic Models for Fractional Calculus[END_REF] (and see also (167)-( 169)) that multiplying (17) by h -1 yields densities which remain finite in the continuous-time limit h → 0. Another important thing here is that the second coefficient tends to zero by a factor h α faster than the first one. Generally we observe that terms of the form ξ µ (1u) -αµ+λ → 0 (λ > 0) giving rise to terms scaling as ∼ h λ → 0 (where τ k = hk and ξ(h) = ξ 0 h α ), namely (See Appendix A.3)

1 h (ξ(h)) µ (αµ -λ) k k! ∼ ξ µ 0 h αµ-1 k αµ-λ-1 Γ (αµ -λ) = h λ τ αµ-λ-1 k Γ (αµ -λ) ∼ h λ τ αµ-λ-1 Γ (αµ -λ) ∼ h λ → 0 (18)
vanishing in the continuous-time limit h → 0. Hence, for the continuous-time limit, only the part ϕ

(ν) α (k) is relevant as 1 -(1 -u) α → 1.
Then consider (15) in the limit h → 0 by using [START_REF] Lévy | Processus semi-Markovien[END_REF] to arrive at

χ α,ν (t) ct = lim h→0 χ α,ν (t) h = lim h→0 ∞ k=1 ψ (ν) α (k)δ h (t -kh) = lim h→0 1 h ψ (ν) α t h , ξ 0 h α (19) with χ α,ν (t) ct = lim h→0 1 h ψ (ν) α k = t h , ξ = ξ 0 h α = ∞ m=0 (-1) m (ν) m m! ξ m+ν 0 t α(m+ν)-1 Γ(α(m + ν)) -h α ξ m+ν 0 t α(m+ν-1)-1 Γ(α(m + ν -1)) = ∞ m=0 (-1) m (ν) m m! ξ m+ν 0 t α(m+ν)-1 Γ(α(m + ν)) , t ∈ R + . ( 20
)
Throughout this paper we commonly use subscript notation (..) ct for continuous-time limit distributions. We can also obtain this result from [START_REF] Takács | On sojourn time problem[END_REF] with h = τ k+1τ k → dτ and τ k = hk → τ and by using the limiting property δ h (tkh) → δ(tτ ) (see again (167)-( 174)).

Hence we can also write [START_REF] Takács | On sojourn time problem[END_REF] in the form

χ α,ν (t) ct = lim h→0 ∞ 0 dτ δ(t -τ ) ∞ m=0 (-1 ) m (ν) m m! ξ m+ν 0 τ α(m+ν)-1 Γ ((m + ν)α) + h α τ α(m+ν-1 )-1 Γ (α(m + ν -1 )) . ( 21
)
The second term tends to zero as h α thus we obtain for h → 0 the density which indeed is the Prabhakar-Mittag-Leffler density [START_REF] Gorenflo | The asymptotic universality of the Mittag-Leffler waiting time law in continuous time random walks[END_REF]. In this way we have shown that the waiting-time probabilities ( 14) are discrete-time approximations of the Prabhakar-Mittag-Leffler density (2), and the underlying discrete-time counting process indeed is a discrete-time version of the GFPP.

χ α,ν (t) ct = ξ ν 0 t αν-1 ∞ m=0 (ν) m m! (-ξ 0 t α ) m Γ ((m + ν)α) = ξ ν 0 t αν-1 E ν α,αν (-ξ 0 t α ), 0 < α ≤ 1 , ν > 0 (22)
In Figure 2 is depicted the behavior of the well-scaled discrete-time density [START_REF] Meerschaert | Stochastic Models for Fractional Calculus[END_REF] versus h for different values of time t. For h → 0 the values converge to the continuous-time Prabhakar density (2) (See [START_REF] Takács | On sojourn time problem[END_REF], [START_REF] Pyke | Markov Renewal Processes with Finitely many states[END_REF]). These plots show as well for fixed h monotonically decreasing values of the well-scaled discrete-time density χ α,ν (t) h for increasing t. This monotonic behavior for the values ν and α used in these plots reflects the complete monotonicity of the continuous-time limit: The Prabhakar density (2) is completely monotonic for αν ≤ 1 with 0 < α ≤ 1 (See [START_REF] Giusti | A practical guide to Prabhakar fractional calculus[END_REF][START_REF] Mainardi | On complete monotonicity of the Prabhakar function and non-Debye relaxation in dielectrics[END_REF] and references therein for a discussion of complete monotonicity of Prabhakar functions). Now we can define the class of Prabhakar discrete-time processes: We call a discrete-time renewal process with waiting times following the distribution P(Z = t) = χ(t) 1 (t ∈ N 0 ) Prabhakar if exists a well-scaled continuous-time limit lim h→0 χ(t) h = χ α,ν (t) ct (in the sense of (167)-( 170)) to the Prabhakar-Mittag-Leffler density [START_REF] Gorenflo | The asymptotic universality of the Mittag-Leffler waiting time law in continuous time random walks[END_REF]. The remaining part of the paper is devoted to the analysis of Prabhakar-discrete time processes and related stochastic motions.

GENERALIZATION

From the above introduced limiting procedures we can infer that further discrete-time generalizations of the Prabhakar process can be obtained by the following class of generating functions

θ(ν) α (u) = f (u) φ(ν) α (u), φ(ν) α (u) = ξ ν (ξ + (1 -u) α ) ν , ν > 0 , α ∈ (0 , 1 ]. (23)
In this expression

f (u) = Eu t = u ∞ k=1 f (k)u k-1 (24)
which can be seen as generating function of any waiting time-distribution, (i.e. with f (u)| u=1 = 1) which fulfills the desired initial condition f (u)| u=0 = f (t) t=0 = 0 (with f (t) ≥ 0 on t ∈ N), we call such a distribution with zero initial condition here simply 'f -distribution'. It is important to notice that f (u) does not contain any scaling parameter ξ, thus in the continuoustime limit f → 1. We will see a little later that the continuous-time limit is uniquely governed by the 'relevant part' φ(ν) α (u) in [START_REF] Angstmann | A Time-Fractional Generalised Advection Equation from a Stochastic Process[END_REF]. We can generate [START_REF] Angstmann | A Time-Fractional Generalised Advection Equation from a Stochastic Process[END_REF] by the following composition procedure 6 leading to Eq. [START_REF] Repin | Fractional Poisson law[END_REF]. Consider the strictly increasing integer time random variable

X n = X 1 , n = 1 X n = X 1 + n j=2 Z j , X 1 , Z j ∈ N, n > 1 ( 25
)
and X 0 = 0. The first step X 1 is a strictly positive random integer following an f -distribution thus Eu X 1 = f (u) whereas the increments Z j for j = 2, . . . n are IID Sibuya(α) with Eu Z = 1 -(1u) α (See Appendix A.4 for details). The integer random variable (25) has generating function

Eu Xn = (w I (u)) (n) = Eu X 1 (Eu Z ) n-1 = f (u) (1 -(1 -u) α ) n-1 . ( 26
)
For the distribution of the events n in [START_REF] Michelitsch | Continuous time random walk and diffusion with generalized fractional Poisson process[END_REF] we choose again [START_REF] Mainardi | A fractional generalization of the Poisson processes[END_REF]. The so defined process has then generating function

Wν α (u) = ∞ n=1 w (δ) B (n)( wI (u)) (n) = ∞ n=1 w (δ) B (n) f (u) (1 -(1 -u) α ) n-1 = f (u) 1 -(1 -u) α ∞ n=1 w (δ) B (n) (1 -(1 -u) α ) n = f (u) 1 -(1 -u) α w(ν) B (1 -(1 -u) α ) = f (u) ( p q ) ν p q + (1 -u) α ν . ( 27
)
In the last line we utilized generating functions [START_REF] Gorenflo | On the Fractional Poisson Process and the Discretized Stable Subordinator[END_REF] and [START_REF] Zaslavsky | Chaos, fractional kinetics, and anomalous transport[END_REF]. Clearly ( 27) is a waiting-time generating function of type [START_REF] Angstmann | A Time-Fractional Generalised Advection Equation from a Stochastic Process[END_REF]. Now it is only a small step to prove that Wν α (u) converges to the Prabhakar density for h → 0. Generating function f (u) (|u| ∈ [0, 1]) can be written as follows

f (u) = ∞ t=1 f (t)(1 -(1 -u)) t = 1 + ∞ k=1 g k (1 -u) k , g k = (-1 ) k ∞ t=1 t k f (t). (28)
We confirm by f (u)| u=1 = ∞ t=1 f (t) = 1 the normalization of the f -distribution 7 . Clearly, in view of the property [START_REF] Smith | Regenerative Stochastic Processes[END_REF] it follows that terms (1u) k produce contributions that tend for h → 0 to zero as h k , namely with u = e -hs we have with [START_REF] Newman | Networks: An Introduction[END_REF]

f (e -hs ) ∼ 1 + ∞ k=1 g k h k s k ∼ 1 + O(h) → 1.
6 See also [START_REF] Pachon | On Discrete-Time Semi-Markov processes[END_REF]. 7 Further, we have

f (u)|u=0 = ∞ k=0 g k = ∞ t=1 f (t)(1 -1) t = 0.

DISCRETE-TIME VERSIONS OF PRABHAKAR-GENERALIZED POISSON DISTRIBUTION

In this section our goal is to analyze a particular important case of Prabhakar discrete-time counting process. We define this process by the strictly increasing random walk

J n = n n=1 Z j , Z j ∈ N, J 0 = 0 , Z j ∈ N (29)
where the Z j are the IID copies of Z (interpreted as waiting time in the related counting process) following a Prabhakar type discrete-time distribution

P(Z = k) = θ (ν) α (t) with
generating function of type [START_REF] Angstmann | A Time-Fractional Generalised Advection Equation from a Stochastic Process[END_REF]. As a proto-example we analyze here the most simple generating function of this type8 namely with f (u) = u, thus

Eu Z = θ(ν) α (u) = ξ ν u (ξ + (1 -u) α ) ν . ( 30
)
For ν = 1 (0 < α < 1) (30) recovers generating function of the fractional Bernoulli counting process (of 'type B') introduced in [16] (Eq. (78) therein). We call the discrete-time counting process with waiting-time generating function [START_REF] Michelitsch | Fractional Dynamics on Networks and Lattices[END_REF] the 'Prabhakar discrete time process' (PDTP). The PDTP stands out by generalizing fractional Bernoulli (type B), and for ν = 1, α = 1 [START_REF] Michelitsch | Fractional Dynamics on Networks and Lattices[END_REF] recovers the generating function θ(1) 1 (u) = pu 1-qu (ξ = p q and p + q = 1) of the standard Bernoulli-process 9 .

The goal is now to derive explicitly the state probabilities of the PDTP and to show that the PDTP converges to the continuous-time Prabhakar renewal process (GFPP) under suitable scaling assumptions. Note also that the PDTP waiting time distribution has the convenient property that it is the distribution {ϕ

α (t) = θ(ν) α ( T-1 )δ t0 = T-1 ξ ν (ξ+(1-T-1 ) α ) ν δ t0 = T-1 ϕ (ν) α (t) = ϕ (ν) α (t -1) = 1 (t -1)! d t-1 du t-1 φ(ν) α (u)| u=0 , φ(ν) α (u) = ξ ν (ξ + (1 -u) α ) ν (31)
where we always utilize causality, i.e. all distributions vanish for negative times. The discrete-time density ϕ (ν) α (t) (t ∈ N 0 ) is evaluated explicitly in Eq. [START_REF] Meerschaert | The Fractional Poisson Process and the Inverse Stable Subordinator[END_REF]. On the other hand relation [START_REF] Riascos | Fractional dynamics on networks: Emergence of anomalous diffusion and Lévy flights[END_REF] is simply reconfirmed by the Leibniz-rule

θ (ν) α (t) = 1 t! d t du t u φ(ν) α (u) u=0 = 1 t! u d t du t φ(ν) α (u) + t d t-1 du t-1 φ(ν) α (u) + 0 u=0 = 1 (t -1)! d t-1 du t-1 φ(ν) α (u) u=0 t ∈ N 0 . ( 32
)
Let us first derive further related distributions such as survival and state probabilities. To this end consider the probability for at least one arrival within [0, t], namely

Ψ (ν) α (t) = t k=1 θ (ν) α (k) = t k=1 ϕ (ν) α (k -1) = 1 t! d t du t u φ(ν) α (u) (1 -u) u=0 = 1 (t -1)! d t-1 du t-1 φ(ν) α (u) (1 -u) u=0 , t ∈ N (33)
with

Ψ (ν) α (t) t=0 = 0 since θ (ν) α (t) t=0 = 0 where the generating function of Ψ (ν) α (t) is Ψ(ν) α (u) = ∞ t=0 u t Ψ (ν) α (t) = u φ(ν) α (u) (1 -u) . ( 34
)
Then the survival probability Φ

(0) α,ν (t) is P(J 1 > t) = Φ (0 ) α,ν (t) = 1 -Ψ ν α (t) = ∞ k=t+1 θ (ν) α (k) = ∞ k=t+1 ϕ (ν) α (k -1 ) (35)
with the generating function

Φ(0) α,ν (u) = ∞ t=0 u t (1 -Ψ ν α (t)) = 1 -u φ(ν) α (u) (1 -u) , |u| < 1 (36)
fulfilling the desired initial condition

Φ(0) α,ν (u) u=0 = Φ (0)
α,ν (t) t=0 = 1 saying that the waiting time J 1 for the first arrival is strictly positive. Then by simple conditioning arguments we obtain the generating function

Φ(n) α,ν (u) of the state probabilities Φ (n)
α,ν (t) (n, t ∈ N 0 ), i.e. the probabilities for n arrivals within [0, t] as 10

Φ(n) α,ν (u) = Φ(0) α,ν (u) u φ(ν) α (u) n = (1 -u φ(ν) α (u)) (1 -u) u n φ(nν) α (u), n ∈ {0, 1, 2, . . .} Φ(n) α,ν (u) = u n φ(nν) α (u) 1 -u - u n+1 φ((n+1)ν) α (u) 1 -u |u| < 1 . (37)
We also mention the normalization of the state probabilities which can be seen by means of the general relation

1 t! d t du t ∞ n=0 Φ(n) α,ν (u) u=0 = ∞ n=0 Φ (n) α,ν (t) = 1 t! d t du t 1 1 -u u=0 = 1 , t ∈ N 0 . ( 38
)
Note that relation [START_REF] Giusti | A practical guide to Prabhakar fractional calculus[END_REF] includes n = 0 where φ0 α (u) = 1 which has a distribution of the form of a discrete-time δ-distribution (See (164) with h = 1)

ϕ (0 ) α (t) = δ 1 (t) = δ t0 ( 39
)
10 Where ( φ(ν)

α (u)) n = ξ nν (ξ+(1-u) α ) nν = φ(nν) α (u). thus θ (0) α (t) = T-1 δ 1 (t) = δ t1
(See also Eq. ( 31)). The 'state-probabilities' are then obtained from

P(J n < t) = 1 t! d t du t Φ(n) α,ν (u) t=0 , t, n ∈ N 0 . ( 40
)
The representation [START_REF] Giusti | A practical guide to Prabhakar fractional calculus[END_REF] is especially convenient for an explicit evaluation of Φ

(n) α,ν (t). Be reminded that the state probabilities are shifted distributions where we account for [START_REF] Santos | Fractional Prabhakar Derivative in Diffusion Equation with Non-Static Stochastic Resetting[END_REF] to arrive at

φ(nν) α ( T-1 ) 1 -T-1 T-n ϕ (0) α (t) = φ(nν) α ( T-1 ) 1 -T-1 δ 1 (t -n) = T-n P (n) α,ν (t) = P (n) α,ν (t -n), P (n) α,ν (k) = P α,νn (k) P α,µ (k) = 1 k! d k du k φ(µ) α (u) 1 -u u=0 , φ(µ) α (u) = ξ µ (ξ + (1 -u) α ) µ . ( 41
)
This result is also obtained from the Leibniz-rule which yields

P α,nν (t -n) = 1 t! d t du t u n φ(nν) α (u) 1 -u u=0 = 1 t! t k=0 t k d k du k u n u=0 d t-k du t-k φ(nν) α (u) 1 -u u=0 = 1 t! t! (t -n)!n! d n du n (u n ) d t-n du t-n φ(nν) α (u) 1 -u u=0 = 1 (t -n)! d t-n du t-n φ(nν) α (u) 1 -u u=0 . ( 42
)
We hence can write for the state-probability distribution

Φ (n) α,ν (t) = P α,νn (t -n) -P α,ν(n+1 ) (t -n -1 ), n, t ∈ N 0 . ( 43
)
To evaluate this expression we account for the expansion with respect to

(1-u) -1 ξ 1 α , namely φ(µ) α (u) 1 -u = ∞ m=0 (-1 ) m (µ) m m! ξ m+µ (1 -u) -α(m+µ)-1 (44)
which converges as [START_REF] Beghin | Fractional Poisson processes and related planar random motions[END_REF] for

|(1 -u) -1 ξ 1 α | < 1. Then we get P α,µ (k) = 1 k! d k du k φ(µ) α (u) 1 -u u=0 = ∞ m=0 (-1 ) m (µ) m m! ξ m+µ (α(m + µ) + 1 ) k k! (45)
where (ρ) m denotes the Pochhammer symbol (4). With relations ( 43) and ( 45) we can write the state-probabilities as 11

Φ (n) α,ν (t) = Φ (n) α,ν (t, ξ) = ∞ m=0 (-1) m m! ξ m+nν × × (nν) m (t -n)! (α(m + nν) + 1) t-n - ξ ν ((n + 1)ν) m (t -n -1)! (1 + α(m + ν(n + 1))) t-n-1 , t ≥ n Φ (n) α,ν (t) = 0, t < n ( 46
)
where n, t = {0, 1, 2, . . . } ∈ N 0 . It follows from the causality of P α,νn (k) in ( 43) that Φ

(n) α,ν (t) = 0 for t < n and hence Φ (n) α,ν (t) t=0 = 0 for n > 0.
It is especially instructive to consider contained special cases in [START_REF] Michelitsch | Sur une généralisation de l'opérateur fractionnaire[END_REF], namely fractional Bernoulli ν = 1 with 0 < α < 1 (subsequent Eq. ( 53)) as well as standard Bernoulli ν = 1 with α = 1 (subsequent Eqs. ( 55), ( 56)). Consider now the survival probability, i.e. n = 0 in [START_REF] Michelitsch | Sur une généralisation de l'opérateur fractionnaire[END_REF], namely

Φ (0) α,ν (t) = 1 - ∞ m=0 (-1) m m! ξ m+ν (ν) m (t -1)! (1 + α{m + ν}) t-1 , t ≥ 1 Φ (0) α,ν (t) t=0 = 1 Φ (0) α,ν (t) = 1 -Ψ (ν) α (t), t ∈ N 0 (47)
where Ψ ν α (t) is the probability of at least one event within [0, t] :

Ψ (ν) α (t) = ∞ m=0 (-1) m m! ξ m+ν (ν) m (t -1)! (1 + α{m + ν}) t-1 , t ≥ 1 Ψ (ν) α (t)| t=0 = 0 (48)
and has generating function [START_REF] Michelitsch | Generalized fractional Poisson process and related stochastic dynamics[END_REF]. Since Ψ ν α (t)| t=0 = 0 we have for initial condition of the survival probability Φ (0) α,ν (t) t=0 = 1. Thus we identify for the state-probabilities (46) the important initial condition

Φ (n) α,ν (t) t=0 = Φ(n) α,ν (u) u=0 = δ n0 , n = {0 , 1 , 2 , . . .} ∈ N 0 . ( 49
)
The initial condition of this form indeed is crucial for many applications of discrete-time renewal processes which come along as Cauchy initial-value problems. By this reason we have constructed generating function [START_REF] Michelitsch | Fractional Dynamics on Networks and Lattices[END_REF] such that it fulfills initial condition θ(ν) α (u) u=0 = 0.

We will come back to this important issue later on in the context of 'discrete-time random walks' (Section 5).

In order to verify that the state probabilities (46) approximate the continuous-time state probabilities of the Prabhakar process, let us consider the continuous-time limiting process 11 We utilize the synonymous notations such as

Φ (n) α,ν (t) and Φ (n)
α,ν (t, ξ), the latter when it is necessary to consider the dependence of parameter ξ (for instance in the continuous-time limit). 52). more closely (Appendices A.2, especially (167)-( 174) for shift-operator properties and general limiting procedures). The continuous-time limit state probabilities are determined by the limiting behavior of the well-scaled state probabilities in the sense of relation ( 172)

Φ (n) α,ν (t) ct = lim h→0 h 1 -T-h T-nh φ(nν) α ( T-h ) -T-(n+1)h φ((n+1)ν) α ( T-h ) δ h (t) = D -1 t ξ nν 0 (ξ 0 + D α t ) nν - ξ (n+1)ν 0 (ξ 0 + D α t ) (n+1)ν δ(t). ( 50
)
Laplace transforming this relation indeed recovers the Laplace transform of the Prabhakar continuous-time state probabilities ( [START_REF] Michelitsch | Continuous time random walk and diffusion with generalized fractional Poisson process[END_REF], Eq. ( 36)). This continuous-time limit is obtained explicitly by performing the well-scaled limit (174) in ( 46) by accounting for the fact that the state probabilities are dimensionless cumulative distributions, namely

Φ (n) α,ν (t) ct = lim h→0 Φ (n) α,ν t h , ξ 0 h α . ( 51
)
By using then the asymptotic relation of the Pochhammer symbol for k = t h large

(µ) k k! = Γ(µ+k) Γ(µ)Γ(k+1) ∼ k µ-1 Γ(µ) we arrive at Φ (n) α,ν (t) ct = (ξ 0 t α ) nν ∞ m=0 (-ξ 0 t α ) m m! (nν) m Γ(αm + ανn + 1) - (ξ 0 t α ) ν ((n + 1)ν) m Γ(αm + αν(n + 1) + 1) = (ξ 0 t α ) nν E nν α,αnν+1 (-ξ 0 t α ) -(ξ 0 t α ) ν E (n+1)ν α,α(n+1)ν+1 (-ξ 0 t α ) , n ∈ N 0 , t ∈ R + (52)
where in this expression appears the Prabhakar function E c a,b (z) [START_REF] Gorenflo | Mittag-Leffler Waiting Time, Power Laws, Rarefaction, Continuous Time Random Walk, Diffusion Limit[END_REF]. Expression (52) indeed coincides with the state probabilities of the continuous-time Prabhakar counting process (Generalized Fractional Poisson process -GFPP) [START_REF] Cahoy | Renewal processes based on generalized Mittag-Leffler waiting times[END_REF][START_REF] Michelitsch | Continuous time random walk and diffusion with generalized fractional Poisson process[END_REF] 12 .

In Figure 3 we draw the GFPP state probabilities [START_REF] Kilbas | Generalized Mittag-Leffler Function and Generalized Fractional Calculus Operators[END_REF] for different states n. The state probabilities exhibit for large t an universal power-law limit which is independent of n (See Eq. ( 61)). The larger n (indicated with brighter colors) the smaller the state probability is for the same time t. This general behavior can be understood with the intuitive picture that states with higher n are less 'occupied' at the same time t.

We show in [START_REF] Michelitsch | Continuous time random walk and diffusion with generalized fractional Poisson process[END_REF] that for ν = 1 relation ( 52) recovers Laskin's fractional Poisson distribution introduced in [START_REF] Laskin | Fractional Poisson process[END_REF] which is also the scaling limit of discrete-time state probabilities [START_REF] Michelitsch | Sur une généralisation de l'opérateur fractionnaire[END_REF] for ν = 1. In order to see explicitly the connection with the fractional Poisson process we obtain for ν = 1 from ( 46) the state probabilities

Φ (n) α,1 (t, ξ) = ξ n (αn + 1) t-n (t -n)! + ∞ m=1 (-1) m ξ m+n (n) m m! (α(n + m) + 1) t-n (t -n)! + (n + 1) m-1 (m -1)! (α(n + m) + 1) t-n-1 (t -n -1)! , t ≥ n Φ (n) α,1 (t, ξ) = 0, t < n ( 53
)
where n, t = {0, 1, 2, . . . } ∈ N 0 . For 0 < α < 1 these are the state-probabilities of the fractional Bernoulli process (type B), and for α = 1 this expression recovers the state probabilities of standard Bernoulli shown subsequently in Eqs. ( 55), (56). The state probabilities [START_REF] Michelitsch | Recurrence of random walks with long-range steps generated by fractional Laplacian matrices on regular networks and simple cubic lattices[END_REF] have with the same limiting procedure the continuous-time limit 13

Φ (n) α,1 (t) ct = lim h→0 Φ (n) α,1 t h , ξ 0 h α = (ξ 0 t α ) n Γ(αn + 1) + ∞ m=1 (-1) m (ξ 0 t α ) m+n Γ(α(m + n) + 1) (n) m m! + (n + 1) m-1 (m -1)! = (ξ 0 t α ) n n! ∞ m=0 (n + m)! m! (-ξ 0 t α ) m Γ(α(m + n) + 1) , t ∈ R + .
(54) We identify [START_REF] Michelitsch | On Recurrence and Transience of Fractional RandomWalks in Lattices[END_REF] indeed with Laskin's fractional Poisson distribution [START_REF] Laskin | Fractional Poisson process[END_REF] which also is recovered for ν = 1 in the expression [START_REF] Kilbas | Generalized Mittag-Leffler Function and Generalized Fractional Calculus Operators[END_REF] for GFPP state probabilities. It follows that the state probabilities (53) of fractional Bernoulli are discrete-time approximations of the fractional Poisson distribution [START_REF] Michelitsch | On Recurrence and Transience of Fractional RandomWalks in Lattices[END_REF]. Now let us consider α = 1, ν = 1, i.e. the case of (standard) Bernoulli more closely. The generating functions of the state probabilities (37) then take the form 12 Eq. (2.13) in [START_REF] Cahoy | Renewal processes based on generalized Mittag-Leffler waiting times[END_REF] and Eq. ( 38) in [START_REF] Michelitsch | Continuous time random walk and diffusion with generalized fractional Poisson process[END_REF]) 13 With

Φ(n) 1 ,1 (u) = (ξ + 1 )(ξu) n (ξ + 1 -u) n+1 , ξ = p q , p + q = 1 , n ∈ N 0 (55)
(n)m m! + (n+1) m-1 (m-1)! = (n+m)! n!m! .
where the state probabilities yield (with p = ξ ξ+1 and q = 1 ξ+1 )

Φ (n) 1,1 (t, ξ) = (ξ + 1)ξ n (t -n)! d t-n du t-n 1 (ξ + 1 -u) n+1 u=0 = (ξ + 1)ξ n (ξ + 1) n+1+t-n n + t -n t -n = t n ξ n (ξ + 1) t = t n p n q t-n , t ≥ n and Φ (n) 1,1 (t, ξ) = 0, t < n (t, n ∈ N 0 ).
(56) We identify (56) with the Binomial distribution (i.e. the state distribution of the Bernoulli process). The continuous-time limit of ( 56) is obtained from p = hξ 0 1+hξ 0 and q = 1 1+hξ 0 thus performing the well-scaled limit (51) yields

Φ (n) 1,1 (t) ct = lim h→0 1 n! t h t h -1 . . . t h -n + 1 h n ξ n 0 (1 + hξ 0 ) -t h t ∈ hN 0 = (ξ 0 t) n n! e -ξ 0 t , t ∈ R + , n ∈ N 0 (57)
which is the Poisson distribution (also recovered in [START_REF] Michelitsch | On Recurrence and Transience of Fractional RandomWalks in Lattices[END_REF] for α = 1) 14 . This reflects the wellknown fact that the standard Bernoulli process converges in a well-scaled continuous-time limit to standard Poisson (see e.g. [START_REF] Pachon | On Discrete-Time Semi-Markov processes[END_REF] and many others).

ASYMPTOTIC FEATURES

In many applications the asymptotic features for large and small observation times are of interest. Clearly the asymptotic behavior for large t in a discrete-time distribution is determined by the leading power in (1u) in the limit u → 1 -0 in the generating function.

The generating function [START_REF] Michelitsch | Fractional Dynamics on Networks and Lattices[END_REF] behaves then as

θ(ν) α (u) = u 1 + 1 ξ (1 -u) α -ν = u ∞ m=0 (-1) m (ν) m m! ξ -m (1 -u) αm θ(ν) α (u) ∼ φ(ν) α (u) ∼ 1 - ν ξ (1 -u) α + O(1 -u) α , (u → 1 -0), α ∈ (0, 1), ν > 0. (58)
The asymptotic behavior of the waiting time density (31) for t large is hence governed by the termν ξ (1u) α . We notice that this term is (up to the scaling multiplier ν ξ ) the same as in the generating function of Sibuya(α) (See Appendix A.4). Thus we get for large t the fat-tailed behavior

ϕ (ν) α (t, ξ) ∼ θ (ν) α (t, ξ) ∼ ∞ m=1 (-1) m (-ν) m m! ξ -m (-αm) t t! ∼ ν ξ (-1) t-1 α t ∼ αν ξ t -α-1 Γ(1 -α) , ϕ (ν) α (t) ct = θ (ν) α (t) ct = lim h→0 1 h ϕ (ν) α t h , ξ 0 h α ∼ αν ξ 0 t -α-1 Γ(1 -α) , α ∈ (0, 1) (59)
where the last line refers to the large-time asymptotic behavior of the continuous-time limit which is of the same type. We notice that this asymptotic distribution is (due to the multiplier ν ξ ) a scaled version of Sibuya(α). It follows that the PDTP converges for large observation times to a scaled version of the Sibuya counting process (Appendix A.5). This is true for all ν > 1 including ν = 1 for the fractional Bernoulli counting process. The asymptotic relation (59) indeed is in accordance with the well-known fat-tailed asymptotic behavior of the Prabhakar density [START_REF] Pillai | Discrete Mittag-Leffler distributions[END_REF] of the GFPP in the continuous-time case (See e.g. [START_REF] Michelitsch | Continuous time random walk and diffusion with generalized fractional Poisson process[END_REF], Eq. ( 35)). It is clear that the same tail (59) occurs in a fractional Bernoulli process (with rescaled parameter ξ ′ = ξ ν ) with discrete-time waiting time density being an approximation of the Mittag-Leffler density. Indeed (59) also is the fat-tail of the Mittag-Leffler density reflecting the asymptotic 'Mittag-Leffler universality' [START_REF] Gorenflo | The asymptotic universality of the Mittag-Leffler waiting time law in continuous time random walks[END_REF].

Let us also consider the asymptotic behavior of the state probabilities. With the same argument we obtain from the generating function [START_REF] Giusti | A practical guide to Prabhakar fractional calculus[END_REF] the leading power in (1u) for u → 1 -0,

namely 15 Φ(n) α,ν (u) = u n (1 -u) 1 + (1 -u) α ξ -nν -u 1 + (1 -u) α ξ -(n+1)ν Φ(n) α,ν (u) ∼ M α,ν (u) ∼ ν ξ (1 -u) α-1 , (u → 1 -0) α ∈ (0 , 1 ). ( 60
)
Here we have accounted for the behavior of the memory generating function M α,ν (u) introduced in (63) for u → 1 -0, with the same leading term as the state probabilities, namely

M α,ν (u) ∼ ν ξ (1-u) α-1 .
We observe in relation (60) that state probability (and memory-) generating functions ∀n ∈ N 0 at u = 1 have a positive weak singularity ∼ (1u) α-1 (α ∈ (0, 1)).

This weak singularity indeed is the origin of the asymptotic large-time power-law decay of the PDTP state probabilities ∀n ∈ N 0 and of the 'memory function' (introduced subsequently in (66)), namely

Φ (n) α,ν (t, ξ) ∼ M α,ν (t) ∼ ν ξ (-1) t α -1 t ∼ ν ξ t -α Γ(1 -α) , Φ (n) α,ν (t) ct = lim h→0 Φ (n) α,ν t h , ξ 0 h α → M α,ν (t) ct ∼ ν ξ 0 t -α Γ(1 -α) (t large) α ∈ (0 , 1 ) (61)
where the last line refers to the continuous-time limit. The power-law (61) is universal and reflects the long-time memory and non-Markovian feature in the PDTP for α ∈ (0, 1). This asymptotic relation is in accordance with the large-time limiting behavior of the continuoustime GFPP state probabilities [START_REF] Michelitsch | Continuous time random walk and diffusion with generalized fractional Poisson process[END_REF]. Physically two 'extreme' regimes are noteworthy: (i) α → 0+: There we have

lim α→0+ t -α Γ(1 -α) = (-1) t α -1 t α=0 = Θ(t) = 1 thus Φ (n)
0+,ν (t, ξ) → ν ξ where extremely long waiting times occur thus the states n 'live' extremely long and the memory of the PDTP becomes infinite. 15 For α = 1 we obtain the finite value

Φ(n) α,ν (u)|u=1 = M α,ν (u)|u=1 = 1 + ν ξ .
In the second 'extreme' regime (ii) we have α → 1 -0 where t -α Γ(1-α) → δ(t) = 0 (t large) indicating lack of memory (for ν = 1) or short-time memory (for ν = 1). This feature also is reflected by the non-singular behavior of the state probability generating functions of standard Bernoulli (55) which yield at u = 1 the finite value

Φ(n) 1,1 (u) u=1 = ξ+1 ξ (∀n ∈ N 0 ).
We observe in these results the following general property of discrete-time counting processes: If the state probability generating functions (and memory generating function) are

weakly singular ∼ (1 -u) λ (λ ∈ (-1, 0)) at u = 1,

GENERALIZED FRACTIONAL DIFFERENCE EQUATIONS GOVERNING THE PDTP STATE

PROBABILITIES

The goal of this part is to derive recursive evolution equations that are solved by the PDTP state probabilities. To this end we utilize the correspondence of generating functions and their operator representations. Whenever we deal with operator-functions of the shift oper-

ator T-h we refer to a renewal chain (29) with rescaled waiting times Z j ∈ hN (h > 0) of the PDTP with generating function [START_REF] Michelitsch | Fractional Dynamics on Networks and Lattices[END_REF]. The state probability generating functions [START_REF] Giusti | A practical guide to Prabhakar fractional calculus[END_REF] fulfill

Φ(n) α,ν (u) φ(ν) α (u) = u Φ(n-1) α,ν (u), n ∈ N (62)
and for n = 0 we have

M α,ν (u) = Φ(0) α,ν (u) φ(ν) α (u) = 1 1 -u 1 + (1 -u) α ξ ν -u . ( 63
)
where we refer (63) to as 'memory generating function'. These simple relations allow us to obtain recursive equations for the state probabilities. We introduce the operator 16 Dν

α (h) = 1 φ(ν) α ( T-h ) = 1 + (1 -T-h ) α ξ ν , Dν α := Dν α (1 ) = 1 φ(ν) α ( T-1 ) (64)
where we skip the argument in Dν α = Dν α (h = 1) when we refer to integer discrete-time processes with h = 1. We can then rewrite (62) and (63) compactly in operator form (See

also Appendix A.2) 17 Dν α • Φ (n) α,ν (t) = Φ (n-1 ) α,ν (t -1 ) + δ n0 M α,ν (t), t, n ∈ N 0 . ( 65
)
(65) can be conceived as generalized fractional Kolmogorov-Feller difference equations governing the time evolution of the PDTP state probabilities. We will determine them subsequently in explicit form including their continuous-time limit representations. For n = 0 (by 16 We utilize as synonymous notations Â-1 = 1 Â and write for unity operator (zero-shift) simply 1 = 1, namely Ta a=0 = 1. 17 Bear in mind the correspondence u ↔ T-h thus 1 u ↔ T+h , see Appendix A.2 and we use in this formulation that Φ (n) α,ν (t) = 0 if at least one of the variables t, n < 0.

accounting for Φ (-1) (t) = 0) we have with (63) the 'memory function' 18

M α,ν (t) = Dν α • Φ (0 ) α,ν (t) = M α,ν 1 δ 1 (t) = 1 1 -T-1 Dν α δ 1 (t) -Θ(t -1 ), t ∈ Z 0 (66) fulfilling initial condition M α,ν (t) t=0 = M ν,α (u) u=0 = (ξ+1) ν ξ ν .
For later use we introduced in (66) the 'memory operator' which has with ξ(h) = ξ 0 h α the well-scaled representation

M α,ν h = h 1 -T-h Dν α (h) -T-h (67) which defines by M α,ν (t) h = M α,ν h δ h (t) (t ∈ hZ 0 , ξ = ξ 0 h α )
the well-scaled memory function maintaining the initial condition of (66). We observe that in the case of Bernoulli α = 1, ν = 1 the memory function becomes 'local', namely M 1,1 h δ h (t) = (1+ξ) ξ hδ h (t) vanishing for t > 0 which reflects the loss of memory in the standard Bernoulli process (See [START_REF] Pachon | On Discrete-Time Semi-Markov processes[END_REF] for general aspects). Now we rewrite above introduced shift-operator functions (64) and (67) in such a way that only existing generalized fractional integrals and derivatives 19 in the continuous-time limit emerge, namely

Dν α (h) = ξ -ν 0 h -να (1 -T-h ) αν 1 + ξ 0 h α (1 -T-h ) -α ν = ξ -ν 0 h -⌈να⌉ (1 -T-h ) ⌈αν⌉ h -(αν-⌈αν⌉) (1 -T-h ) αν-⌈αν⌉ 1 + ξ 0 h α (1 -T-h ) -α ν = ξ -ν 0 h -⌈να⌉ (1 -T-h ) ⌈αν⌉ Bα,ν (h) (68)
where we introduced the ceiling function ⌈µ⌉ which indicates the smallest integer larger or equal to µ (See Appendix A.3). In Appendix A.

Bα,ν (h) = h -(αν-⌈αν⌉) (1 -T-h ) αν-⌈αν⌉ 1 + ξ 0 h α (1 -T-h ) -α ν (69)
which can be conceived as a discrete-time generalized fractional integral operator 20 . In view of (67) we further introduce

K (0 ) αν (h) = h 1 -T-h Dν α (h) = ξ -ν 0 h 1 -⌈να⌉ (1 -T-h ) ⌈αν⌉-1 Bα,ν (h). (70)
All these operators are operator functions of the shift operator T-h = e -hDt and can be considered as discrete-time versions of of generalized fractional derivatives or integrals, respectively. The discrete-time memory function (66) then can be evaluated as (See also 18 In (66) we use

h 1-T-h δ h (t) = ∞ k=0 hδ h (t -kh) = Θ(t) thus h T-h 1-T-h δ h (t) = Θ(t -h), see Appendix A.2. M α,ν (t) = 1 1 -T-1 Dν α δ 1 (t) -T-1 Θ(t) = K (0) α,ν (t) -Θ(t -1) K (0) α,ν (t) = 1 1 -T-1 Dν α δ 1 (t) = ξ -ν 1 -T-1 ⌈αν⌉-1 B α,ν (t), t ∈ Z 0 B α,ν (t) = Bα,ν (1)δ 1 (t) = ∞ m=0 (-ν) m m! (-ξ) m (αm + ⌈αν⌉ -αν) t t! , t ∈ N 0 . (71) Now let us consider Bα,ν (h)Φ(t) = ∞ m=0 (-ν) m m! (-ξ 0 ) m h(1 -T-h ) -1 αm+⌈αν⌉-αν Φ(t) = ∞ m=0 (-ν) m m! (-ξ 0 ) m ∞ k=0 (αm + ⌈αν⌉ -αν) k k! h (αm+⌈αν⌉-αν) Φ(t -kh) t ∈ hZ 0 ( 72
)
where we put ξ = ξ 0 h α in order to have well-scaled operators with existing continuous-time limits. This relation is a generalized fractional difference equation with memory where the whole history {Φ(tkh)} (t ∈ hN 0 ) of the causal function Φ(t) contributes. Now since αm + ⌈αν⌉αν > 0 if αν / ∈ N and αm + ⌈αν⌉αν ≥ 0 if αν ∈ N we can use for the continuous-time limit that

(αm+⌈αν⌉-αν) k k! h (αm+⌈αν⌉-αν) → τ αm+⌈αν⌉-αν-1 Γ(αm+⌈αν⌉-αν) dτ (h → dτ and τ = kh, Appendix A.3) thus D -(αm+⌈αν⌉-αν) t Φ(t) ct = lim h→0 h(1 -T-h ) -1 αm+⌈αν⌉-αν Φ(t) = t 0 Φ(t -τ ) ct τ (αm+⌈αν⌉-αν)-1 Γ(αm + ⌈αν⌉ -αν) dτ (73)
are Riemann-Liouville fractional integrals of orders (αm + ⌈αν⌉αν). Note that for αν ∈ N where ⌈αν⌉αν = 0 the order m = 0 the kernel becomes lim αν→⌈αν⌉ τ (⌈αν⌉-αν)-1 Γ(⌈αν⌉-αν) = δ(τ ) a Dirac's δ-distribution thus (73) recovers Φ(t) ct . The continuous-time limit of (72) then yields

lim h→0 Bα,ν (h)Φ(t) = t 0 B α,ν (τ ) ct Φ(t -τ ) ct dτ = ∞ m=0 (-ν) m m! (-ξ 0 ) m ∞ 0 τ αm+⌈αν⌉-αν-1 Γ(αm + ⌈αν⌉ -αν) Φ(t -τ ) ct dτ (74)
with continuous-time limit kernel

B α,ν (τ ) ct =                  τ ⌈αν⌉-αν-1 ∞ m=0 (-ν) m m! (-ξ 0 τ α ) m Γ(αm + ⌈αν⌉ -αν) , αν / ∈ N δ(τ ) + ∞ m=1 (-ν) m m! (-ξ 0 ) m τ αm-1 Γ(αm) , αν ∈ N τ ∈ R + =          τ ⌈αν⌉-αν-1 E -ν α,⌈αν⌉-αν (-ξ 0 τ α ), αν / ∈ N D τ Θ(τ )E -ν α,1 (-ξ 0 τ α ) , αν ∈ N τ ∈ R + (75) containing the Prabhakar function E c a,b (z) (3). Then it follows that Dν α (h)Φ(t) = ξ -ν 0 h -1 (1 -T-h ) ⌈αν⌉ ∞ m=0 (-ν) m m! (-ξ 0 ) m × ∞ k=0 (αm + ⌈αν⌉ -αν) k k! h (αm+⌈αν⌉-αν) Φ(t -kh), t ∈ hZ 0 ( 76
)
is well-scaled in the sense of limiting equation (174). Thus with (70) we first get for the continuous-time limit

K (0) α,ν (t) ct = lim h→0 K(0) α,ν (h)δ h (t) = ξ -ν 0 D ⌈αν⌉-1 t t 0 B α,ν (τ ) ct δ(t -τ )dτ = ξ -ν 0 D ⌈αν⌉-1 t B α,ν (t) ct t ∈ R + . ( 77
)
Note that, in this limiting procedure derivatives of integer orders come into play with ⌈αν⌉-1 = 0 for 0 < αν ≤ 1 and ⌈αν⌉ -1 ∈ N for αν > 1. We notice that kernel (75) has Laplace transform L(B α,ν (t) ct )(s) = s -⌈αν⌉ (s α + ξ 0 ) ν . The continuous-time limit kernel (75) is for αν / ∈ N by the term m = 0 weakly singular. We then get for the continuous-time limit of (76)

lim h→0 Dν α (h) • Φ(t) = ξ -ν 0 D ⌈αν⌉ t t 0 B α,ν (t -τ ) ct Φ(τ ) ct dτ. (78)
The results (77), (78) are in accordance with the continuous-time expressions derived for the GFPP (Eqs. (66)-(68) in [START_REF] Michelitsch | Continuous time random walk and diffusion with generalized fractional Poisson process[END_REF]) 21 .

Having derived these continuous-time limits, our goal is now to deduce the recursive generalized fractional difference equations governing the state probabilities Φ

(n) α,ν (t) of the PDTP.

To this end let us rewrite Eq. (65) more explicitly. Taking into account the evaluations (71), ( 72) and (76) this equation takes the form (with h = 1)

ξ -ν 1 -T-1 ⌈αν⌉ ∞ m=0 (-ν) m m! (-ξ) m ∞ k=0 (αm + ⌈αν⌉ -αν) k k! Φ (n) α,ν (t -k) = Φ (n-1) α,ν (t -1) + δ n0 ξ -ν 1 -T-1 ⌈αν⌉-1 ∞ m=0 (-ν) m m! (-ξ) m (αm + ⌈αν⌉ -αν) t t! -Θ(t -1)
(79) 21 Where the definition of K where n, t ∈ N 0 . For later use and further analysis of the structure it appears instructive to write this equation by accounting for (71) as

(1 -T-1 ) ∞ k=0 K (0 ) α,ν (k)•Φ (n) α,ν (t -k) = Φ (n-1 ) α,ν (t -1 )+δ n0 K (0 ) α,ν (t) -Θ(t -1 ) , n, t ∈ N 0 . (80)
We observe that the k-series (due to causality and property Φ 65), ( 79) and (80) indeed are equivalent discretetime versions of the generalized fractional Kolmogorov-Feller equations governing the stateprobabilities in a PDTP. These equations are non-local discrete-time convolutions with non-Markovian long-time memory effects where the whole history {Φ

(n) α,µ (t -k) = 0 for t -k < n) has upper limit k = t -n. The Equations (
(n) α,ν (t -k)} comes into play.
We will come back to this issue again later on. With (77) and ( 78) we can directly write the well-scaled continuous-time limit of (80) in the form

ξ -ν 0 D ⌈αν⌉ t t 0 B α,ν (t -τ ) ct Φ (n) α,ν (τ ) ct dτ = Φ (n-1 ) α,ν (t) ct + δ n0 M α,ν (t) ct , t ∈ R + (81)
which is solved by the generalized fractional Poisson distribution functions (GFPP state probabilities) Φ

(n) α,ν (t) ct of ( 52) with the continuous-time limit kernel B α,ν (t) ct (75). This equation contains the continuous-time limit of the memory function

M α,ν (t) ct = lim h→0 M α,ν h δ h (t) = ξ -ν 0 D ⌈αν⌉-1 t B α,ν (t) ct -Θ(t) (82)
where M α,ν h is the well-scaled memory operator of (67).

Continuous-time limit memory function (82) and 'generalized fractional Kolmogorov-Feller equations' (81) indeed are in accordance with the results for the GFPP [START_REF] Michelitsch | Continuous time random walk and diffusion with generalized fractional Poisson process[END_REF].

Let us analyze now more closely the important regime ν = 1 with α ∈ (0, 1].

(i) Fractional Bernoulli counting process (of type B), ν = 1 and α ∈ (0, 1):

We obtain then for the above introduced quantities

K (0) α,1 (t) = 1 ξ (1 -T-1 ) α-1 + (1 -T-1 ) -1 δ 1 (t) = 1 ξ (1 -T-1 ) α-1 δ 1 (t) + Θ(t) t ∈ Z 0 (83) thus (1 -T-1 )K (0 ) α,1 (t) = D 1 α δ 1 (t) = 1 ξ (1 -T-1 ) α + 1 δ t0 . ( 84
)
The memory function (66) then yields

M α,1 (t) = K (0) α,1 (t) -Θ(t -1) = 1 ξ (1 -T-1 ) α-1 + 1 δ 1 (t), 0 < α ≤ 1 M α,1 (t) = 1 ξ 1 t! (1 -α) t + δ 1 (t) = (-1) t ξ α -1 t + δ t0 , t ∈ N 0 M 1,1 (t) = ξ + 1 ξ δ t0 , α = 1, t ∈ N 0 (85)
and the well-scaled form

M α,1 (t) h = h 1 -T-h (D 1 α (h) -T-h )δ h (t) = 1 + 1 ξ (1 -T-h ) α-1 hδ h (t), t ∈ hZ 0 . (86)
For later use we notice the continuous-time limit M α,1 (t

) ct = lim h→0 M α,1 ( t h , ξ 0 h α ) = 1 ξ 0 t -α Γ(1-α)
reflecting long-time power-law memory of fractional Bernoulli process. We used here

(1 -T-h ) α-1 δ h (t) = h -1 (1-α) k k! k= t h (See Appendix A.2).
The discrete-time fractional Kolmogorov-Feller equation then reads (δ

1 (t) = δ t0 ) (1 -T-1 ) α Φ (n) α,1 (t) = ξΦ (n-1) α,1 (t -1) -ξΦ (n) α,1 (t) + δ n0 (-1) t α -1 t + ξδ t0 , α ∈ (0, 1), t ∈ N 0 . ( 87
)
This equation coincides with the fractional difference equation for the fractional Bernoulli counting process (type B) given in [START_REF] Pachon | On Discrete-Time Semi-Markov processes[END_REF] (See Proposition 7, with Eqs. ( 81)-(82) therein). Eq. ( 87) is solved by the state probabilities [START_REF] Michelitsch | Recurrence of random walks with long-range steps generated by fractional Laplacian matrices on regular networks and simple cubic lattices[END_REF]. Indeed for ν = 1 generating function of the state probabilities [START_REF] Giusti | A practical guide to Prabhakar fractional calculus[END_REF] with

φ(1) α (u) = ξ ξ+(1-u) α yields Φ(n) α,1 (u) = (uξ) n [ξ + (1 -u) α-1 ] [ξ + (1 -u) α ] n+1 , n ∈ N 0 ( 88
)
which was also given in [START_REF] Pachon | On Discrete-Time Semi-Markov processes[END_REF]. On the other hand (88) recovers for α = 1 the generating function (55) of standard Binomial distribution. We can also recover (87) from Eq. ( 79) for ν = 1 by accounting for (-1) m is nonzero only for m = 0 and m = 1 and where ⌈α⌉ = 1. The Eq. ( 87) is the discrete-time fractional Kolmogorov-Feller equation of this process where the memory function for t large behaves as M (0)

α,1 (t) ∼ 1 ξ t -α
Γ(1-α) thus fractional Bernoulli has a long-time power-law memory and non-Markovian features. Now let us consider more closely the continuous-time limit of Eq. (87). First let us write Eq. (65) by accounting for the memory operator (67) and with (64) in well-scaled form

(ξ = ξ 0 h α ) to arrive at (1 -T-h ) α Φ (n) α,1 (t) h = ξ 0 h α Φ (n-1) α,1 (t -h) h -ξ 0 h α Φ (n) α,1 (t) h + δ n0 ξ 0 h α + (1 -T -h ) α-1 hδ h (t) t ∈ hZ 0 . (89)
The solution of this equation can be written in well-scaled operator representation (See Eq. ( 50) for ν = 1)

Φ (n) α,1 (t) h = (ξ 0 h α + (1 -T-h ) α-1 ) (ξ 0 h α + (1 -T-h ) α ) n+1 (ξ 0 h α ) n T-hn hδ h (t), t ∈ hZ 0 . ( 90
)
The limit h → 0 of (90) yields for continuous-time limit Φ

(n) α,1 (t) ct = ξ n 0 D α-1 t (ξ 0 +D α t ) n+1 δ(t) (with Laplace transform L{Φ (n) α,1 (t) ct }(s) = ξ n 0 s α-1
(ξ 0 +s α ) n+1 ) which can be identified with Laskin's fractional Poisson distribution [START_REF] Michelitsch | On Recurrence and Transience of Fractional RandomWalks in Lattices[END_REF] where n = 0 recovers the standard Mittag-Leffler survival probability [START_REF] Laskin | Fractional Poisson process[END_REF][START_REF] Michelitsch | Continuous time random walk and diffusion with generalized fractional Poisson process[END_REF] (among many others). The continuous-time limit of Eq. (89) gives

lim h→0 h -α (1 -T-h ) α Φ (n) α,1 (t) h = lim h→0 ξ 0 Φ (n-1 ) α,1 (t -h) h -Φ (n) α,1 (t) h +δ n0 t -α Γ (1 -α) + ξ 0 h δ h (t) . ( 91
)
By using lim h→0 hδ h (t) = 0, however δ h (t) → δ(t) together with

h -α (1 -T -h ) α-1 (hδ h (t)) = h -(α-1) (1 -T -h ) α-1 δ h (t) → D α-1 t δ(t) = t -α Γ(1 -α) ,
we obtain as continuous-time limit of (89) the fractional differential equation

D α t • Φ (n) α,1 (t) ct = +ξ 0 Φ (n-1 ) α,1 (t) ct -ξ 0 Φ (n) α,1 (t) ct + δ n0 t -α Γ (1 -α) , t ∈ R + , n ∈ N 0 . ( 92
)
In this equation occurs the Riemann-Liouville fractional derivative D α t of order α (0 < α < 1)

(See Appendix A.3 for details). Eq. ( 92) coincides with the fractional Kolmogorov-Feller equation given by Laskin [START_REF] Laskin | Fractional Poisson process[END_REF] and is solved by the continuous-time state probabilities of Laskin's fractional Poisson distribution [START_REF] Michelitsch | On Recurrence and Transience of Fractional RandomWalks in Lattices[END_REF].

(ii) Bernoulli counting process, ν = 1, α = 1:

Eq. ( 87) reduces then to

(1 -T-1 )Φ (n) 1 ,1 (t) = ξΦ (n-1 ) 1 ,1 (t -1 ) -ξΦ (n) 1 ,1 (t) + (ξ + 1 )δ n0 δ t0 , α = 1 , t ∈ Z 0 ( 93
)
where the memory function M 1,1 (t) = (ξ+1) ξ δ t0 = 0 is null for t > 0. This observation explains our choice of the name 'memory function' for M α,ν (t). The Bernoulli process is memoryless and Markovian and Eq. ( 93) is solved by the Binomial distribution (56). Eq. ( 93) indeed is the Kolmogorov-Feller difference equation of the corresponding Bernoulli process with ν = 1, α = 1. To see this let us consider (93) for t = 0 and n = 0:

Φ (0 ) 1 ,1 (0 ) -Φ (0 ) 1 ,1 (-1 ) = ξΦ (-1 ) 1 ,1 (-1 ) -ξΦ (0 ) 1 ,1 (0 ) + (ξ + 1 ) (94) by using Φ (-1) 1,1 (0) = Φ (-1)
1,1 (-1) = 0 we recover in this equation the initial condition Φ (0) 1,1 (0) = 1. For t ≥ 1 from Eq. (93) we have then for n = 0

Φ (0 ) 1 ,1 (t) = 1 ξ + 1 Φ (0 ) 1 ,1 (t -1 ), t ≥ 1 . ( 95
)
Now it follows from this recursion and the initial condition Φ (0)

1,1 (0) = 1 that the survival probability in the Bernoulli process is

Φ (0 ) 1 ,1 (t) = 1 (ξ + 1 ) t = q t , t ≥ 0 (96)
also in accordance with (56) for n = 0. Then let us finally consider (93) for n ≥ 1 and t ≥ n which writes

Φ (n) 1 ,1 (t) = 1 ξ + 1 Φ (n) 1 ,1 (t -1 ) + ξΦ (n-1 ) 1 ,1 (t -1 ) , t ≥ n ≥ 1 (97)
which is indeed solved by (56), namely

Φ (n) 1,1 (t) = 1 ξ + 1 t -1 n ξ n (ξ + 1) t-1 + ξ t -1 n -1 ξ n-1 (ξ + 1) t-1 = ξ n (ξ + 1) t t -1 n + t -1 n -1 = t n ξ n (ξ + 1) t = t n p n q t-n , t ≥ n = t n p n q t-n Θ(t -n) (98)
and Φ

(n)

1,1 (t) = 0 for t < n. Hence (98) together with (96) shows that (93) indeed is solved by the state probabilities of the (standard) Bernoulli process (56). We also verify easily the normalization

∞ n=0 Φ (n) 1,1 (t) = t n=0
t n p n q t-n = (p + q) t = 1.

The continuous-time limit of (93) then yields (See also Eq. (89

) for α → 1 -0) lim h→0 (1 -T-h )Φ (n) 1,1 (t) h = ξ 0 h Φ (n-1) 1,1 (t -h) h -ξ 0 h Φ (n) 1,1 (t) h + δ n0 (ξ 0 h + 1)h δ h (t), α = 1, t ∈ hN 0 . ( 99
)
Thus we get with (1 -T-h ) → hD t and δ h (t) → δ(t) the equation

D τ Φ (n) 1 ,1 (t) ct = ξ 0 Φ (n-1 ) α,1 (τ ) ct -ξ 0 Φ (n) α,1 (t) ct + δ n0 δ(t), t ∈ R (100)
which is indeed solved by the Poisson distribution Φ

(n) 1,1 (t) ct = Θ(t) (ξ 0 t) n
n! e -ξ 0 t (See also limiting equation (57)) and with Φ

(0) 1,1 (t) ct = e -ξ 0 t Θ(t) thus D t Φ (0) 1,1 (t) ct = -ξ 0 Φ (0)
1,1 (t) ct + δ(t) occurs in (100) for n = 0. The limiting equation (100) reflects the fact that the Bernoulli process is a discrete-time version of the Poisson process, see also [START_REF] Pachon | On Discrete-Time Semi-Markov processes[END_REF] among many others and consult also [START_REF] Michelitsch | Continuous time random walk and diffusion with generalized fractional Poisson process[END_REF] (Eqs. (79), (80) in that paper.). This limiting equation can be recovered from (92) in the limit α → 1 -0 when we account for (See also Appendix A.3) that lim α→1-0 D α τ → D τ reproduces the standard first-order derivative and memory function

lim α→1-0 M α,1 (t) ct = 1 ξ 0 t -α Γ(1-α) = 1
ξ 0 δ(t) becomes a Dirac distribution (See also Eq. ( 61))

reflecting the memoryless Markovian feature of standard Poisson.

EXPECTED NUMBER OF ARRIVALS IN A PDTP

For many applications especially in problems of diffusive motion, the expected number of arrivals within [0, t] is of utmost importance. For the PDTP this quantity is defined as

n α,ν (t) = ∞ n=0 nΦ (n) α,ν (t), t ∈ N 0 (101)
with the state probabilities Φ

(n) α,ν (t) of Eq. ( 46). We introduce the important generating function

Ḡα,ν (u, v) = Ḡα,ν (u, ξ, v) = Φ(0) α,ν (u) ∞ n=0 (vu φν α (u)) n = 1 -u φν α 1 -u 1 1 -uv φν α (u) = 1 (1 -u) ([1 + 1 ξ (1 -u) α ] ν -u) ([1 + 1 ξ (1 -u) α ] ν -uv) = M α,ν (u) ([1 + 1 ξ (1 -u) α ] ν -uv) , |v| < 1 , |u| ≤ 1 .
(102) We will see subsequently (Section 5) that this generating function is a key-quantity for stochastic motions with PDTP waiting times. We refer this generating function to as 'PDTP generating function' as it contains the complete stochastic information such as state prob- abilities and the expected number of arrivals and also memory generating function (63).

It is now only a small step to directly extract from (102) the asymptotic behavior in the time domain for large t. Consider (102) for u → 1 -0 which yields for α ∈ (0, 1) the same weakly singular behavior as in the PDTP state probability generating functions and memory generating function (See Eqs. ( 60), (61)), namely 22Ḡα,ν (u, v) ∼ M α,ν (u)

1 -v ∼ ν ξ(1 -v) (1 -u) α-1 , u → 1 -0 , |v| < 1 , α ∈ (0 , 1 ). ( 103
)
Accounting for (61) we then get for the PDTP generating function the same type of asymptotic power-law behavior for t large, namely

G α,ν (t, ξ, v) ∼ ν ξ(1 -v) t -α Γ (1 -α) α ∈ (0 , 1 ). ( 104
)
For subsequent use let us also write the scaled time domain representation

G α,ν (t, ξ, v) h = G α,ν t h , ξ, v = h (1 -T-h ) ([1 + 1 ξ (1 -T-h ) α ] ν -T-h ) ([1 + 1 ξ (1 -T-h ) α ] ν -v T-h ) δ h (t), t ∈ hZ 0 (105)
which conserves the initial condition of the process with h = 1, namely G α,ν (t, ξ, v) h t=0 = Ḡα,ν (u, v) u=0 = 1. With the scaling ξ = ξ 0 h α (105) and T-h = e -hDt the well-scaled continuous-time limit writes

G(t, ξ 0 , v) ct = lim h→0 G α,ν t h , ξ 0 h α , v = D -1 t [1 + 1 ξ 0 D α t ) α ] ν -1 [1 + 1 ξ 0 D α t ) α ] ν -v δ(t), t ∈ R (106)
with Laplace transform

L{G(t, ξ 0 , v) ct }(s) = 1 s [1 + 1 ξ 0 s α ] ν -1 [1 + 1 ξ 0 s α ] ν -v . ( 107
)
The expected number of arrivals (101) are then given by

n α,ν (t) = 1 t! ∂ t+1 ∂v∂u t Ḡα,ν (u, v) u=0 ,v=1 = 1 t! ∂ t ∂u t N α,ν (u) u=0 ( 108
)
where

N α,ν (u) = u φν α (u) (1 -u)(1 -u φν α (u)) (109)
is the generating function for the expected number of arrivals in (108). Note that all series are convergent for |v| < 1, |u| ≤ 1 (and for |v| ≤ 1 when |u| < 1) where | φν α (u)| ≤ 1. For the further evaluation it is useful to expand (109) as

N α,ν (u) = 1 (1 -u) ∞ n=1 u n φnν α (u). (110)
From (108) together with relations ( 42), [START_REF] Barbu | Semi-Markov Chains and Hidden Semi-Markov Models toward Applications[END_REF] follows that the terms n > t are zero thus

n α,ν (t) = t n=1 P α,nν (t -n) = ∞ m=0 Q (m) α,ν (t), t, n ∈ N 0 (111) with Q (m) α,ν (t) = t n=1 (-1 ) m m! ξ νn+m (nν) m (t -n)! (α(m + nν) + 1 ) t-n , t, n ∈ N 0 ( 112
)
where P α,nν (k) is the probability P α,nν (tn) of ( 45) for at least n arrivals within [0, t].

Now we derive the continuous-time limit by taking into account that the expected number of arrivals ( 101) is a dimensionless (in the sense of relation (174) a cumulative) distribution.

To this end we account for

lim h→0 ξ νn+m 0 h α(νn+m) ( t h -n)! (α(m + nν) + 1 ) t h -n = (ξ 0 t α ) (m+νn) Γ (αm + ανn + 1 ) , t ∈ R. (113)
The continuous-time limit of (112) then is obtained by virtue of relation (174) leading to

Q (m) α,ν (t) ct = lim h→0 Q t h , ξ 0 h α = ∞ n=1 (-1 ) m m! (nν) m (ξ 0 t α ) (m+νn) Γ (αm + ανn + 1 ) . ( 114
)
Hence the continuous-time limit of (111) yields

n α,ν (t) ct = ∞ n=1 (ξ 0 t α ) νn ∞ m=0 (nν) m m! (-ξ 0 t α ) m Γ(αm + ανn + 1) = ∞ n=1 ξ νn 0 t ανn E νn α,ανn+1 (-ξ 0 t α ) = ∞ n=0 ξ ν(n+1) 0 t αν(n+1) E ν(n+1) α,αν(n+1)+1 (-ξ 0 t α ) t ∈ R + (115)
where E c a,b (z) denotes the Prabhakar function (3). The result (115) coincides with the expression given in [START_REF] Cahoy | Renewal processes based on generalized Mittag-Leffler waiting times[END_REF] (Eq. (2.36) in that paper).

In diffusion problems especially of interest is the asymptotic behavior for large t. With the same argument as in Section (4.1) we infer that this asymptotic behavior is contained in the dominating term for u → 1, i.e. in the leading power of (1u). To cover this part we rewrite generating function (109)

N α,ν (u) = 1 (1 -u) u φν α (1 -u φν α ) ∼ 1 (1 -u) φν α (1 -φν α ) (u → 1 -0) N α,ν (u) ∼ ξ ν (1 -u) -α-1 (u → 1 -0) α ∈ (0 , 1 ]. ( 116
)
where this asymptotic formula also includes α = 1. Hence the large-time behavior of (101) yields n α,ν (t) ∼ ξ ν

(1 + α) t t! ∼ ξ ν t α Γ (1 + α) , (t large), α ∈ (0 , 1 ]. (117)
This expression contains, by accounting for (174), also the continuous-time limit for large observation times n α,ν (t

) ct = lim h→0 n α,ν ( t h , ξ 0 h α ) = ξ 0 ν t α
Γ(1+α) . We note that the average number of arrivals within [0,t] of (117) is a dimensionless function of t (ξ 0 has physical dimension sec -α .). Also we mention that ν for large t enters the power-law (117) only as a scaling parameter where the smaller ν the shorter the waiting times leading to an increase of (117) for a fixed t.

It is further noteworthy that in the fractional interval α ∈ (0, 1) the power-law (117) reflects asymptotic self-similar scaling. In other words: For long observation times the event stream converges to a stochastic fractal on the time line with scaling dimension 0 < α < 1 indicating a disjoint 'dust like' fractal. For α = 1 this behavior turns into linear law indicating compact coverage of the time line by the events. Indeed for α = 1 relation (117) recovers the linear behavior of standard Bernoulli, and of standard Poisson in the continuous-time limit. The asymptotic power-law (117) together with (106), (107) are in accordance with the results obtained for the GFPP [START_REF] Michelitsch | Continuous time random walk and diffusion with generalized fractional Poisson process[END_REF].

Particularly instructive is again the standard Bernoulli case α = 1 and ν = 1. (109) takes then the particularly simple form

N 1 ,1 (u) = ξ (ξ + 1 ) u (1 -u) 2 = ξ ξ + 1 ∞ k=1 ku k . (118) Thus p q = ξ n 1 ,1 (t) = ξ ξ + 1 t = pt, t ∈ N 0 . ( 119
)
We can easily reconfirm this well-known result by employing (101) accounting for the Bernoulli process state probabilities (56). Then we get

n 1,1 (t) = t n=0 t n np n q t-n = p d dp (p + q) t = pt.
By virtue of Eq. ( 174), we can directly derive the continuous-time limit as

n 1 ,1 (t) ct = lim h→0 ξ 0 h (ξ 0 h + 1 ) t h = ξ 0 t, t ∈ R + (120)
in accordance with above considerations. This result indeed is the well-known linear law for the expected number of arrivals of the standard Poisson process (see e.g. [START_REF] Michelitsch | Continuous time random walk and diffusion with generalized fractional Poisson process[END_REF] among many others) 23 .

The results of this subsection are particularly useful in problems of diffusive particle motions with PDTP waiting times between the particle jumps. We devote subsequent section to this subject.

= d dv e -ξ 0 t ∞ n=0 v n (ξ 0 t) n n! v=1 = d dv e ξ 0 t(v-1) v=1 = ξ0t.

PRABHAKAR DISCRETE TIME RANDOM WALK ON UNDIRECTED GRAPHS

The goal of the present section is to analyze the stochastic motion on undirected graphs governed by the PDTP. We consider a walker which performs random steps between connected nodes where the IID waiting times between the steps are drawn from a discrete-time renewal process. If in contrast the sejourn times on the nodes follow a continuous-time renewal process, then the resulting motion is a Montroll-Weiss continuous-time random walk (CTRW) [START_REF] Montroll | Random walks on lattices II[END_REF]. Instead we consider here walks where the waiting times between the steps follow a PDTP. The so defined walk is a discrete-time generalization of the classical Montroll-Weiss CTRW. We call this Montroll-Weiss type of walk 'Discrete-time random walk' (DTRW).

Consider now a random walker on an undirected connected graph of N states (nodes). To characterize the walk on the network we introduce the N × N one-step transition matrix (also referred to as stochastic matrix) H = (H ij ) where the matrix elements 0 ≤ H ij ≤ 1 indicate the conditional probabilities that the walker who is sitting on node i in one step moves to node j. The transition matrix is normalized as N j=1 H ij = 1 (i.e. row-stochastic). From row-stochasticity of H follows row-stochasticity of its matrix powers24 H n (n ∈ N 0 ) with the probability (H n ) ij that the walker in n steps moves from node i to j. The matrix H n is the n-step transition matrix. The one-step transition matrix relates the topological information of the network with the random walk and is defined by [START_REF] Noh | Random walks on complex networks[END_REF][START_REF] Riascos | Fractional dynamics on networks: Emergence of anomalous diffusion and Lévy flights[END_REF][START_REF] Newman | Networks: An Introduction[END_REF][START_REF] Michelitsch | Fractional Dynamics on Networks and Lattices[END_REF]]

H ij = 1 K i A ij = δ ij - L ij K i . ( 121
)
In this relation is introduced the adjacency matrix A ij = A ji (symmetric in an undirected network) with A ij = 1 if the pair of nodes ij is connected by an edge, and A ij = 0 otherwise. Further to force the walker to change the node in each step we have A jj = 0, j = 1, . . . N , i.e. there are no self-connections. The matrix L = (L ij ) is referred to as Laplacian-matrix with L ij = K i δ ij -A ij . Due to this structure we have L ii = K i where K i = N j=1 A ij is the degree of node i indicating the number of neighbor (connected) nodes with node i.

We see in (121) that the inverse degrees play the role of row-normalization factors. In an undirected graph the adjacency matrix and Laplacian matrices are symmetric whereas the transition-matrix (121) in general is not if there is a pair of nodes i, j such that K i = K j .

The spectral structure of Laplacian and transition matrix is analyzed in details in [START_REF] Michelitsch | Fractional Dynamics on Networks and Lattices[END_REF] (and see the references therein). We assume here a connected (ergodic) graph. In such a graph the transition matrix has spectral representation

H = |v 1 v1 | + N m=2 λ m |v m vm | (122)
with unique eigenvalue λ 1 = 1 and |λ m | < 1 (m = 2, . . . , N ) where |v 1 v1 | indicates the stationary (invariant) distribution with lim n→∞ H n = |v 1 v1 |. We ignore here for simplicity cases of so called bipartite graphs where an eigenvalue -1 occurs. For an outline consult [START_REF] Michelitsch | Fractional Dynamics on Networks and Lattices[END_REF] and the references therein. For our convenience we employ in this section Dirac's bra|ket -notation 25 We analyze a Montroll-Weiss DTRW where the IID waiting times between the steps follow a PDTP discrete-time density with generating function [START_REF] Michelitsch | Fractional Dynamics on Networks and Lattices[END_REF]. Let P(t) = (P ij (t)) be the transition matrix of this walk. The element P ij (t) indicates the probability that the walker who is sitting at t = 0 on node i is present at time t on node j. Since the walker is always somewhere on the network this transition matrix as well is row-stochastic, i.e. fulfills 0 ≤ P ij (t) ≤ 1 and is normalized as N j=1 P ij (t) = 1. We assume the initial condition that at t = 0 the walker is sitting at node i which is expressed by P ij (t) t=0 = δ ij . Then by simple conditioning arguments we take into account that the walker can move from i to j within the time interval [0, t] in n = 0, 1, , . . . ∈ N 0 steps where the occurrence of n steps is governed by the PDTP state probabilities P(J n < t) = Φ (n) α,ν (t) (i.e. the probabilities that the walker makes n steps within [0, t]). The transition matrix of a DTRW can be written as a Cox-series [START_REF] Cox | Renewal Theory[END_REF] 26 , namely

P α,ν (t) = ∞ n=0 H n Φ (n) α,ν (t), P ij (t) t=0 = δ ij , t ∈ N 0 (123)
containing here the PDTP state probabilities Φ

(n) α,ν (t) of ( 46). We refer this walk to as Prab- hakar DTRW. Transition matrix (123) is a matrix function of H. It follows from ( 123) that in all equations subsequently derived the matrices have a common base of eigenvectors with H and hence are commuting among each other and with H.

So P α,ν (t) • H = H • P α,ν (t)
is commuting since our initial condition is P(0) = 1 (commuting with H) 27 . Keep in mind that P ij (t) t=0 = δ ij indeed requires the state probabilities to fulfill the initial conditions Φ (n) α,ν (t) t=0 = δ n0 which per construction is fulfilled (See Eq. ( 49)).

For the proofs to follow we utilize the generating function of the transition matrix (123) 28 Pα,ν (u

) = 1 -u φ(ν) α 1 -u ∞ n=0 H n u n φνn α (u) = 1 -u φ(ν) α 1 -u 1 -u φ(ν) α H -1 = M α,ν (u) [1 + 1 ξ (1 -u) α ] ν 1 -uH -1 = Ḡα,ν (u, ξ, H) (124
P α,ν ij (u) = Φ(n) α,ν (u) = 1
1-u . Using (122) with (124) yields for the transition matrix (123) the canonic representation

P α,ν (t) = G α,ν (t, H) = |v 1 v1 | + N m=2 G α,ν (t, λ m )|v m vm |, t ∈ N 0 (125) with time-domain representation G α,ν (t, λ) = G(t, ξ, λ) given in (105). We considered that for λ 1 = 1 we have G α,ν (t, λ) λ=1 = 1, ∀t ∈ N 0 .
26 In a general DTRW the state-probabilities of the respective discrete-time counting process are replacing the PDTP state-probabilities. 27 Clearly P(t) and H do not commute if the initial condition is such that H • P(0) = P(0) • H, i.e. if initial condition P(0) and H do not have the same base of eigenvectors. 28 defined by Pα,ν (u) = ∞ t=0 P α,ν (t)u t It follows from the asymptotic behavior of the state probabilities, memory function (61) and PDTP generating function (asymptotic relations (103) and ( 104)) that the stationary distribution for large t is approached by the power-law29

P α,ν (t) ∼ |v 1 v1 | + ν ξ t -α Γ (1 -α) N m=2 |v m vm | 1 -λ m , α ∈ (0 , 1 ). ( 126
)
This relation shows the non-Markovian long-time memory feature of the Prabhakar DTRW.

Noteworthy here the already mentioned limit α → 0+ where extremely long waiting times between the steps occur with t -α Γ(1-α) ∼ Θ(t) = 1. The walk then becomes 'infinitely slow' thus there is a range for t large where P α,ν (t)

∼ |v 1 v1 | + ν ξ N m=2
|vm vm| 1-λm and thus the walk 'struggles' to take the stationary distribution |v 1 v1 | (which eventually is taken since

t -α
Γ(1-α) → 0+ for infinitesimally small positive α). Now our goal is to derive the evolution equation that governs the PDTP transition matrix. To this end we account for Ḡα,ν (u, λ)

φ(ν) α (u) = uλ Ḡα,ν (u, λ) + M α,ν (u) (127)
with memory generating function M α,ν (u) (63). Rewriting (127) in operator form yields

Dν α G α,ν (t, λ) = λG α,ν (t -1 , λ) + M α,ν (t), t ∈ N 0 (128)
containing the memory function M α,ν (t) defined in (71). We then can write (128) in matricial representation as

Dν α P α,ν (t) -M α,ν (t)1 = H • P α,ν (t -1 ), P α,ν ij (0 ) = δ ij . ( 129
)
This equation is the Kolmogorov-Feller generalized fractional difference equation that governs the Prabhakar DTRW on the network 30 . Since M α,ν (t

) t=0 = D ν α δ 1 (t) t=0 = (ξ+1) ν ξ ν
fulfill the same initial condition, Eq. ( 129) recovers for t = 0 the initial condition of the transition matrix (where due to causality the right-hand side for t = 0 is null).

Eq. ( 129) refers to the general class of equations governing discrete-time semi-Markov chains given in [START_REF] Pachon | On Discrete-Time Semi-Markov processes[END_REF] (See Theorem 3.4.). We can conceive (129) as the discrete-time Cauchy problem which is solved by transition matrix (125) where the complete history {P α,ν ij (tk)} (0 ≤ k ≤ t) of the walk comes into play. This shows the following representation of (129), namely

(1 -T-1 ) ∞ k=0 K (0 ) α,ν (k)P α,ν ij (t -k) -δ ij K (0 ) α,ν (t) -Θ(t -1 ) = N r=1 H ir P α,ν rj (t -1 ), t ∈ N 0
(130) also reflecting Eq. ( 80). This equation contains on the right-hand side the topological information of the graph (See (121)). We observe as a consequence of causality (P α,ν ij (tk) = 0 for k > t) that the upper limit of the k-summation on the left-hand side is k = t. Eqs. (129), (130) are equivalent representations of the discrete-time Kolmogorov-Feller equations governing the stochastic motion of a Prabhakar DTRW. These equations are explicit accounting for relations (71) and (76).

Continuous-time limit

Then by the same scaling arguments as in previous sections and outlined in Appendix A.2, it is not a big deal to establish the continuous time limit of these equations. In view of (80) having continuous-time limit (81) we arrive at

ξ -ν 0 D ⌈αν⌉ t 0 B α,ν (τ ) ct P α,ν ij (t -τ ) ct dτ -δ ij ξ -ν 0 D ⌈αν⌉-1 t B α,ν (t) ct -Θ(t) = N r=1 H ir P α,ν rj (t) ct , ξ -ν 0 D ⌈αν⌉-1 t 0 B α,ν (τ ) ct D t P α,ν ij (t -τ ) ct dτ = N r=1 (H rj -δ ij )P α,ν ir (t) ct t ∈ R + (131)
where δ ij = P α,ν ij (t) ct | t=0 is the initial condition and in the formulation of the last line we used on the right-hand side row-normalization of the transition matrix. The continuous-time limit kernel B α,ν (t) ct was determined in (75) and 131) is in accordance with the 'generalized fractional Kolmogorov-Feller equation' derived for the Prabhakar CTRW on undirected networks [25].

K (0) α,ν (t) ct = ξ -ν 0 D ⌈αν⌉-1 t B α,ν (t) ct in relation (77). Eq. (
Case ν = 1 , 0 < α ≤ 1:

Let us discuss the case ν = 1 with 0 < α < 1 more closely, i.e. the walk subordinated to fractional Bernoulli (type B). We refer this walk to as Fractional Bernoulli Walk (FBW).

(130) then takes with Eqs. ( 83)-( 85) the form of a fractional difference equation:

1 -T-1 α P α,1 ij (t)-δ ij (-1 ) t α -1 t + ξδ t0 = -ξP α,1 ij (t)+ξ N r=1 H ir P α,1 rj (t-1 ), t ∈ N 0 . ( 132
)
This equation for t = 0 yields the initial condition P α,1 ij (t) t=0 = δ ij . The memory term

(-1) t α -1 t ∼ t -α Γ(1 -α) , (t large), α ∈ (0, 1)
reflects the long-time memory and non-Markovianity of the process where this memory for α → 0 + 0 becomes extremely long where t -α Γ(1-α) → Θ(t) = 1; whereas α → 1 -0 represents the memoryless limit with t -α Γ(1-α) → δ(t) = 0. In order to analyze the continuous-time limit we account for the generating function (124) of the transition matrix

Pα,1 (u) = (ξ + (1 -u) α-1 ) [ξ(1 -uH) + 1(1 -u) α ] -1 . ( 133
)
The transition matrix of the FBW has the well-scaled operator representation31

P α,1 (t) h = (ξ + (1 -T-h ) α-1 ) ξ(1 -H T-h ) + 1(1 -T-h ) α -1 hδ h (t), t ∈ hZ 0 . ( 134
)
It follows from representation (134) that P α,1 (t) solves the Cauchy problem

(1 -T-h ) α P α,1 (t) h -1(ξ + (1 -T-h ) α-1 ) hδ h (t) = ξ(H T-h -1)P α,1 (t) h P α,1 (t) h t=0 = 1 t ∈ hZ 0 (135)
where this equation is the scaled version of (132). This equation is also consistent with the fractional difference equations (87), (89) for the fractional Bernoulli state probabilities.

Consider generating function ( 102) with (105) for ν = 1, namely

E α t h , ξ, v = G α,1 t h , v = ∞ n=0 Φ (n) α,1 t h , ξ v n = ξ + (1 -T-h ) α-1 ξ(1 -v T-h ) + (1 -T-h ) α hδ h (t), t ∈ hZ 0
(136) which also is obtained by using Eq. ( 90). The scaled state distribution of the fractional Bernoulli process ( 53) is then obtained from this generating function by

Φ (n) α,1 t h , ξ = 1 n! d n dv n E α t h , ξ, v v=0 . ( 137
)
This is the analogue equation as for the fractional Poisson distribution in the continuoustime limit which is shown a little later. We also show that generating function ( 136) is a discrete-time approximation of the standard Mittag-Leffler function. The transition matrix (134) can then be written in the form of the matrix function

P α,1 (t) h = E α t h , ξ, H , t ∈ hZ 0 . ( 138
)
Let us first consider the continuous-time limit of generating function (136)

E α (t, ξ 0 , v) ct = lim h→0 E α t h , ξ 0 h α , v = lim h→0 ξ 0 h + h 1-α (1 -T-h ) α-1 ξ 0 (1 -v T-h ) + h -α (1 -T-h ) α δ h (t) = D α-1 t ξ 0 (1 -v) + D α t δ(t) = E α (-ξ 0 (1 -v)t α ), t ∈ R + (139)
conserving the initial condition E α (t, ξ 0 , v) ct t=0 = 1 of (136). In this relation the standard

Mittag-Leffler function E α (z) = ∞ m=0 z m Γ (αm + 1 ) (140) comes into play. This result is easily confirmed in view of Laplace transform L{E α (-at α )}(s) = s α-1
a+s α of the Mittag-Leffler function. In this way we have shown that generating function (136) also is a discrete-time approximation of the Mittag-Leffler function (139). Indeed Laskin's fractional Poisson distribution [START_REF] Michelitsch | On Recurrence and Transience of Fractional RandomWalks in Lattices[END_REF] is obtained from the Mittag-Leffler generating function by [START_REF] Laskin | Fractional Poisson process[END_REF][START_REF] Michelitsch | Generalized space-time fractional dynamics in networks and lattices Generalized Space-Time Fractional Dynamics in Networks and Lattices[END_REF] (and many others)

Φ α,1 (t) ct = 1 n! d n dv n E α (-ξ 0 (1 -v)t α ) v=0 (141)
which is also the continuous-time limit of Eq. ( 137). The result (139) allows us to get the continuous-time limit of the transition matrix (138) in the form of the Mittag-Leffler matrix function

P α,1 (t) ct = lim h→0 E α t h , ξ 0 h α , H = E α (-ξ 0 (1 -H)t α ) = |v 1 v1 | + N m=2 |v m vm |E α (-ξ 0 (1 -λ m )t α ), t ∈ R + (142)
conserving initial condition P α,1 (t) ct t = 0 = 1. Accounting for the Mittag-Leffler asymptotic relation (holding for the eigenvalues |λ| < 1)

E α (-ξ 0 (1 -λ)t α ) ∼ 1 (1 -λ)ξ 0 t -α
Γ(1 -α) for large t, we observe (142) agrees with (126) for ν = 1.

The continuous-time limit of Eq. (132) yields

D α t P α,1 ij (t) ct -δ ij t -α Γ(1 -α) = ξ 0 N r=1 (H ir -δ ir )P α,1 rj (t) ct , 1 Γ(1 -α) t 0 (t -τ ) -α D τ P α,1 ij (τ ) ct = - ξ 0 K i N r=1 L ir P α,1 rj (t) ct t ∈ R + , 0 < α < 1 (143)
where the first line contains the Riemann-Liouville fractional derivative D α t of order α whereas in the second line we utilize the Caputo-fractional derivative of order α. The fractional evolution equation ( 143) with initial condition P α,1 ij (t) ct t=0 = δ ij indeed is solved by the the Mittag-Leffler transition matrix (142).

Eq. ( 143) is the fractional Kolmogorov-Feller equation governing fractional diffusion in the network, i.e. a CTRW in the network subordinated to a fractional Poisson process. The fractional differential continuous-limit equation ( 143) refers to the class of equations governing semi-Markov processes related to α-stable subordinators in [START_REF] Pachon | On Discrete-Time Semi-Markov processes[END_REF] (Eqs. ( 14), [START_REF] Meerschaert | Stochastic Models for Fractional Calculus[END_REF] in that paper). Fractional differential equations of this type (mostly for infinite continuous spaces) with Mittag-Leffler solutions occur in a wide range of problems in fractional dynamics and anomalous diffusion (see e.g. [START_REF] Zaslavsky | Chaos, fractional kinetics, and anomalous transport[END_REF][START_REF] Mainardi | A fractional generalization of the Poisson processes[END_REF][START_REF] Michelitsch | Continuous time random walk and diffusion with generalized fractional Poisson process[END_REF][START_REF] Michelitsch | Generalized fractional Poisson process and related stochastic dynamics[END_REF][START_REF] Metzler | The Random Walk's Guide to Anomalous Diffusion : A Fractional Dynamics Approach[END_REF][START_REF] Metzler | The restaurant at the end of the random walk: recent developments in the description of anomalous transport by fractional dynamics[END_REF]).

It remains us to consider the case ν = 1, α = 1 which is a walk subordinated to a standard Bernoulli process. We refer this walk to as Bernoulli Walk (BW). We then get for (130) the difference equation

1 -T-1 P 1 ,1 ij (t) -δ ij δ t0 (ξ + 1 ) = ξ N r=1 (H ir T-1 -δ ir )P 1 ,1 rj (t), t ∈ N 0 ( 144
)
where with T-1 P 1,1 rj (t) = P 1,1 rj (t -1) and the memory term ∼ δ t0 is null for t > 0, i.e. the process is memory-less and Markovian reflecting these properties of standard Bernoulli process. Eq. ( 144) is in accordance with the difference equation given in [START_REF] Pachon | On Discrete-Time Semi-Markov processes[END_REF] holding for Markov chains (See Eq. (65) in that paper). Let us rewrite Eq. ( 144) also in matricial representation

(1 + ξ)P 1 ,1 (t) -δ t0 1(ξ + 1 ) = (1 + ξH) P 1 ,1 (t -1 ), t ∈ N 0 . ( 145
)
By using causality, i.e. P 1,1 (-1) = 0 this equation recovers for t = 0 the initial condition P 1,1 (0) = 1. For t ≥ 1 this equation yields the recursion

P 1 ,1 (t) = 1 1 + ξ (1 + ξH) P 1 ,1 (t -1 ), P 1 ,1 (0 ) = 1. (146)
Thus by iterating this recursion yields for the transition matrix P 1,1 (t) = 1 (1+ξ) t (1 + ξH) t .

On the other hand in view of the Binomial distribution (98) we can obtain this result also by employing the Cox series (123), namely

P 1 ,1 (t) = P 1 ,1 (t, ξ) = 1 (1 + ξ) t t n=0 t n ξ n H n = 1 (1 + ξ) t (1 + ξH) t = (pH + q1) t , t ∈ N 0 (147)
where p = ξ 1+ξ and q = 1 1+ξ . In view of (147) and also with Binomial distribution (98) it follows that p is the probability that the walker makes a step and q that it does not make a step in a time unit 1. For p = 1, q = 0 the walker hence makes (almost surely) t steps up to time t and (147) recovers P 1,1 (t) = H t the definition of the t-step transition matrix. On the other hand for p = 0, q = 1 the walker (almost surely) does not move thus in this case the walker remains on his departure node with P 1,1 (t) = 1.

Now we can directly derive the continuous-time limit of (147) by considering the process on t ∈ hN 0 to arrive at (See Eq. ( 174))

P 1,1 (t) ct = lim h→0 P 1,1 t h , ξ 0 h = lim h→0 1 (1 + ξ 0 h) t h (1 + ξ 0 hH) t h = e -ξ 0 t e ξ 0 Ht = e -ξ 0 t(1-H) , t ∈ R + . ( 148
)
This result is also consistent with the Cox-series generated with the Poisson-distribution state probabilities

P 1,1 (t) ct = Θ(t)e -ξ 0 t ∞ n=0 H n (ξ 0 t) n n! = Θ(t)e -ξ 0 t e ξ 0 Ht = Θ(t)(|v 1 v1 | + N m=2 |v m vm |e ξ 0 (λm-1)t ), t ∈ R (149)
where we added here the Heaviside function to emphasize causality. By accounting for the spectral structure of H, namely λ 1 = 1 and m -1 eigenvalues with λ m -1 < 0 we have lim t→∞ e ξ 0 (λm-1)t → 0 for m = 2, . . . N thus the transition matrix P 1,1 (t) ct → |v 1 v1 | approaches asymptotically the stationary distribution. We also notice that the matrix exponential (148), ( 149) is recovered for α = 1 by the Mittag-Leffler transition matrix (142).

Finally the continuous-time limit of Eq. ( 144) is obtained as

D t P 1,1 ij (t) ct -δ ij δ(t) = -ξ 0 P 1,1 ij (t) ct + ξ 0 N r=1 H ir P 1,1 rj (t) ct = - ξ 0 K i N r=1 L ir P 1,1 rj (t) ct D t P 1,1 (t) ct -1δ(t) = ξ 0 (H -1) P 1,1 (t) ct t ∈ R (150)
with initial condition P 1,1 (t) ct | t=0 = 1. The continuous-time limit (150) is the Kolmogorov- Feller equation governing the transition probabilities in continuous-time Markov chains, i.e. CTRWs on graphs with underlying standard Poisson process. It is straight-forward to see that Cauchy problem (150) is solved by the continuous-time limit exponential transition matrix (148), (149). The Eq. ( 150) can be also be recovered from the fractional case (143) in the limit α → 1 -0. Kolmogorov-Feller equations are widely used as master equations to model Markovian walks on graphs (Markov chains) and normal diffusion in multidimensional infinite spaces [START_REF] Pachon | On Discrete-Time Semi-Markov processes[END_REF][START_REF] Michelitsch | Fractional Dynamics on Networks and Lattices[END_REF] (and see the references therein).

In this section we focused on walks subordinated to PDTPs which we refer to as Prabhakar DTRWs and their special cases such as fractional Bernoulli and standard Bernoulli as well as their classical continuous-time limits. The advantage of the Prabhakar DTRW with the PDTP counting process of three parameters α ∈ (0, 1], ν > 0, ξ > 0 is the great flexibility to adapt real-world stochastic processes.

Another interesting case (though not being Prabhakar) is the DTRW with Sibuya counting process. A detailed analysis of this walk is beyond the scope of the present paper. The essential aspects are analyzed in [START_REF] Pachon | On Discrete-Time Semi-Markov processes[END_REF]. We confine us here to a brief outline in the spirit of a Montroll Weiss Sibuya DTRW model in Appendix A.5.

INFLUENCE OF WAITING TIME INITIAL CONDITIONS ON DTRW FEATURES

This section is devoted to briefly analyze the effect of the initial condition of the discretetime waiting-time density in a DTRW. In this paper we have constructed discrete-time counting processes with discrete-time waiting time densities w(t) on t ∈ {0, 1, 2, . . .} with w(t) t=0 = 0. Here our aim is to consider the effect of a discrete-time counting process which has w(0) = 1ǫ with 0 < ǫ ≤ 1. The waiting time generating function of such a process then is given by

Eu Z = w(u) = 1 -ǫ + u ∞ k=1 u k-1 w(k). (151)
The state probabilities (probabilities for n arrivals within [0, t]) then have the generating function

Φ(n) (u) = 1 -w(u) 1 -u (w(u)) n , n = 0 , 1 , 2 , . . . ( 152
)
Now let us consider the initial condition of the state probabilities

Φ (n) (t) t=0 = Φ(n) (u) u=0 = ǫ(1 -ǫ) n , n = 0 , 1 , 2 , . . . ( 153
)
with normalization

∞ n=0 Φ (n) (t) t=0 = ǫ ǫ = 1.
(i) The first observation is that the initial condition of the survival probability is Φ (0) (t) t=0 = ǫ. The good initial condition Φ (n) (0) = δ n0 is only fulfilled for ǫ = 1 and as a consequence w(0) = 1ǫ = 0.

(ii) The second observation is that for ǫ → 0 we have w(0) = 1-ǫ → 1 and hence Φ (0) (t) t=0 → 0 for the survival probability. Hence at t = 0 at least one event already has arrived (almost surely).

Consider now a Montroll-Weiss DTRW with this discrete-time counting process. Then with (123) the generating function of the transition matrix by accounting for Eq. (152) writes

P(u) = ∞ n=0 Φ(n) (u)H n = 1 -w(u) 1 -u (1 -w(u)H) -1 . ( 154
)
The 'natural' initial condition is P(t = 0) = (δ ij ) when the walker at t = 0 is sitting on a given initial node i. However, this natural initial condition is fulfilled if and only if the state probabilities fulfill the initial condition Φ (n) (t = 0) = δ n0 which is true only for ǫ = 1 (See (153)). Let us now consider the effect of ǫ < 1 on the initial condition of the transition matrix, namely

P(t) t=0 = P(u)| u=0 = ∞ n=0 ǫ(1 -ǫ) n H n = |v 1 v1 | ∞ n=0 ǫ(1 -ǫ) n + ǫ N m=2 ∞ n=0 ((1 -ǫ)λ m ) n |v m vm | (155)
P(0) = |v 1 v1 | + ǫ N m=2 |v m vm | 1 -λ m (1 -ǫ) = |v 1 v1 | + (1 -w(0)) N m=2 |v m vm | 1 -λ m w(0)
.

We observe that for ǫ → 0 (i.e. w(0) → 1 -0) the transition matrix already at t = 0 takes the stationary distribution P(0) → |v 1 v1 | (and not as it is in a 'good' DTRW for t → ∞). Then the departure node of the walker becomes maximally uncertain in the sense of a random initial condition where the walker at t = 0 (almost surely) makes a huge number of steps to reach 'immediately' the stationary equilibrium distribution which remains unchanged

P(t) = |v 1 v1 | ∀t.
In view of this consideration we can formulate an 'uncertainty principle' as follows: The more P(Z = t) = w(t) is 'localized versus t = 0', the more variable is the (random) initial condition of the transition matrix. As a consequence for w(0) > 0 the transition matrix does not solve a Cauchy initial value problem and the departure node is 'uncertain'. On the other hand for ǫ → 1 (i.e. w(0) → 0) the initial condition becomes P(t = 0) = 1 thus the DTRW then is a 'good walk' with a well-defined ('certain') departure node.

We notice that relation (155) is consistent with recent results obtained by a different approach [START_REF] Riascos | Random walks on networks with stochastic resetting[END_REF].

CONCLUSIONS

In this paper we analyzed discrete-time renewal processes and their continuous-time limits. We focused especially on counting processes which are 'Prabhakar'. These processes are discrete-time approximations of the Prabhakar continuous-time renewal process (GFPP). Among the discrete-time Prabhakar processes one process (the PDTP) stands out as it contains for special choices of parameters (namely for ν = 1 and 0 < α < 1) the fractional Bernoulli counting process and (for ν = 1 and α = 1) the standard Bernoulli process. The PDTP and the class of discrete-time 'Prabhakar' counting processes converge in well-scaled continuous-time limits to the continuous-time Prabhakar process.

The PDTP is constructed such that zero waiting times between events are forbidden leading to strictly positive interarrival times. This 'good' initial condition makes the PDTP useful to define a Montroll-Weiss DTRW that solves a Cauchy problem with a well-defined departure node of the walker. We called this walk Prabhakar DTRW. We derived for the Prabhakar DTRW generalized fractional discrete-time Kolmogorov-Feller equations that govern the resulting stochastic dynamics on undirected graphs. The Prabhakar DTRW (unless for α = 1 and ν = 1) is non-Markovian where in the range α ∈ (0, 1) and ν > 0 for long observation times universal (Mittag-Leffler) power-law long-time memory effects emerge with fat-tailed waiting time densities.

We demonstrate explicitly that for certain choices of parameters the generalized fractional difference Kolmogorov-Feller equations of the Prabhakar DTRW turn into their classical counterparts of fractional Bernoulli and standard Bernoulli, and the same is true for their continuous-time limits recovering fractional Poisson and Poisson, respectively.

Although we focused on undirected graphs, a promising field of applications of the PDTP and Prabhakar DTRW arises from biased walks on directed graphs (See [START_REF] Riascos | Non-local biased random walks and fractional transport on directed networks[END_REF] for a recent model of fractional dynamics on directed graphs). Among these cases Prabhakar DTRWs come along as strictly increasing walks with interesting applications such as 'aging in complex systems'. These problems exhibit strictly increasing random quantities ('damagemisrepair accumulation') [START_REF] Riascos | Aging in transport processes on networks with stochastic cumulative damage[END_REF]. Generally the Prabhakar DTRW approach opens a huge potential of new interdisciplinary applications to topical problems as various as the time evolution of pandemic spread, communication in complex networks, dynamics in public transport networks, anomalous relaxation, collapse of financial markets, just to denominate a few examples.

A APPENDICES -SUPPLEMENTARY MATERIALS

A.1 DISCRETE CONVOLUTIONS AND GENERATING FUNCTIONS

In this appendix we review some basic properties of generating functions (discrete Laplace transforms) which are powerful analytical tools. The generating function is defined by

ḡ(u) = Eu Z = ∞ k=0 u k g(k), |u| ≤ 1 ( 156
)
where we assume P(Z = k) = g(k) is a discrete-time waiting time probability density in a discrete-time counting process. Normalization is expressed by ḡ(u)| u=1 = 1, Eu Z stands for expectation value of u Z and Z ∈ {0, 1, 2, ..} = N 0 a.s. indicates a discrete random variable which takes the values Z = k with probability P(Z = k) = g(k). We notice the correspondence of generating functions with Laplace transforms by introducing the density (PDF)

χ(t) = ∞ k=0 δ(t -k)g(k), t ∈ R (157)
defined on continuous-time where δ(τ ) stands for Dirac's δ-distribution. Density (157) has Laplace transform

L{χ(t)}(s) = ∞ 0 - e -st χ(t)dt = Ee -sZ = ḡ(e -s ), ℜs ≥ 0 (158)
which is the generating function (156) for u = e -s and ḡ(e -s ) s=0 = ∞ k=0 g(k) = 1 reflects normalization. The discrete convolution operator 32 of two causal functions g(t), h(t) defined on t ∈ N 0 is defined as

(g ⋆ h)(t) = (h ⋆ g)(t) = t j=0 g(j)h(t -j) = 1 t! d t du t {ḡ(u) h(u)}| u=0 , t ∈ N 0 ( 159
)
where ḡ(u) h(u) = ∞ t=0 u t g(t) ⋆ h(t) is the generating function of the discrete convolution. In the following we denote the nth convolution power with

(g(t)⋆) n = 1 t! d t du t (ḡ(u)) n | u=0 , t, n ∈ N 0 (160)
which includes n = 0 and yields unity (g(t)⋆) 0 = δ t0 (where δ ij denotes the Kronecker symbol, which in the following will also be denoted by the equivalent notation δ i,j ). Relation (160) can be extended to non-integer (especially fractional) convolution powers n ∈ R.

ḡ( T-1 )δ t0 = ∞ k=0 g(k)δ tk = g(t) = 1 t! d t du t ḡ(u)| u=0 , t ∈ N 0 . ( 161
)
In the present paper we extensively make use of the correspondence of generating functions and shift-operator counterparts. Such procedures turn out to be especially useful in the definition and determination of suitably scaled continuous-time limits and governing discrete-time evolution equations.

A.2 SHIFT OPERATORS AND CAUSAL DISTRIBUTIONS

We consider a causal distribution Θ(t)g(t) having all its non-zero values on t ∈ {0, h, 2h, . . .} = hZ 0 (h > 0) and zero values for t < 0. For its definition we make use of the discrete-time Heaviside function defined by

Θ h (t) = Θ(t) =      1, t ∈ {0, h, 2h, . . . } 0, t ∈ {-h, -2h, . . .} t ∈ hZ 0 , h > 0 . ( 162
)
We especially emphasize that with this definition Θ h (0) = 1. In the entire paper for Heaviside functions we write Θ h (t) = Θ(t), i.e. we may skip subscript h. Without loss of generality we may identify the discrete-time Heaviside function with its continuous-time counterpart, i.e. the conventional Heaviside-unit-step function with Θ(t) = 1 for t ≥ 0 (especially Θ(0) = 1) and Θ(t) = 0 for t < 0 defined on t ∈ R. This is necessary when we use T-h = e -hDt leading to the 'distributional representation' e -hDt Θ(t) = Θ(th) which is defined on t ∈ R. This allows to define Laplace transformation of δ h (t) as in (166) which has 'good properties' in the limit for h → 0.

Let δ k,l (we also use notation δ kl ) be the Kronecker-δ defined for any positive and negative integer including zero by

δ i,j = δ i-j,0 =      1, i = j 0, i = j i, j ∈ Z 0 ( 163
)
where we have the translational invariance property δ i+m,j+m = δ i,j (m ∈ Z 0 ). Then we define the 'discrete-time δ-distribution' as follows

δ h (t) = Θ(t) -Θ(t -h) h = 1 -T-h h Θ(t) = 1 h δ t h ,0 =        1 h , t = 0 0, t = 0 t ∈ hZ 0 (164)
where δ t h ,0 in this relation indicates the Kronecker-Symbol (163) where we notice that t h ∈ Z 0 . We observe that for h = 1 we have δ 1 (t) = δ t0 . Formula (164) defines a discrete-time density on t ∈ hZ 0 . However, it makes sense to extend this definition to t ∈ R. With this extended definition δ h (t) = Θ(t)-Θ(t-h) h becomes an integrable distribution in the Gelfand-Shilov sense [START_REF] Gel'fand | Generalized Functions[END_REF] with δ h (t) = 1 h for t ∈ [0, h) and δ h (t) = 0 else. Then we have ∞ -∞ δ h (τ )dτ = h 0 δ h (τ )dτ = 1 and δ h (t) has a well-defined Laplace transform (subsequent relation (166)). Of great importance is the continuous-time limit of the discrete-time δ-distribution

lim h→0 δ h (t) = lim h→0 1 -e -hD t h Θ(t) = D t Θ(t) = δ(t), t ∈ R (165)
recovering Dirac's continuous-time δ-distribution. Further, it appears also useful to briefly consider the Laplace transform of δ h (t), namely

L{δ h (t)}(s) = ∞ 0 - δ h (t)e -st dt = ∞ -∞ e -st (1 -e -hDt ) h Θ(t)dt = ∞ -∞ Θ(t) (1 -e +hDt ) h e -st dt = 1 -e -hs h ∞ 0 e -st dt ℜ{s} > 0 = 1 -e -hs hs (166)
where indeed lim h→0 L{δ h (t)}(s) = L{δ(t)}(s) = 1 recovers the Laplace transform of Dirac's δ-distribution reflecting (165). We applied here partial integration and we used T-h = e -hDt having adjoint operator ( T-h ) † = T+h = e +hDt (with e hDt e -st = e -hs e -st ).

Continuous-time limit of discrete-time densities

Let g(k) be a 'discrete-time density' such as for instance a discrete-time waiting time distribution P(Z = k) = g(k, C) (k ∈ N 0 ) and C stands for an internal parameter (or a set of parameters). We often write g(t) instead of g(t, C) if only the time dependence is of interest. We assume g(t) = Θ(t)g(t) (t ∈ Z 0 ) to be a causal distribution. Then we can define the 'scaled causal discrete-time density' by accounting for (164) as follows

g(t) h = ḡ( T-h )δ h (t) = ∞ k=0 g(k, C) T-hk δ h (t) = ∞ k=0 g(k, C)δ h (t -kh) = Θ(t)h -1 g t h , C , t ∈ hZ 0 = h -1 g t h , C , t ∈ hN 0 (167)
which due to the summation over non-negative k (167) is non-zero only for t ≥ 0 (causality). Here the process takes place on t ∈ hN 0 and the scaled discrete-time density g(t) h = h -1 g( t h , C) (t ∈ hN 0 ) as well as g(t) h=1 = g(t, C) (t ∈ N 0 ) both have physical dimension sec -1 and merit therefore to be called (discrete-time) densities 33 . When h is finite the parameter C does not necessarily need to be scaled with h. However, we will see subsequently that generally we need to introduce a scaling law C = C(h) = C 0 h ρ such that the continuoustime limit h → 0 of (167) exists. The constant C 0 > 0 is independent of h and is a free parameter. C 0 has units sec -ρ (as C = C 0 h ρ is dimensionless). Now let us define the continuous-time limit density of (167) and introduce the variable τ k = hk ∈ hN 0 and rescaled parameter C(h) = C 0 h ρ with scaling dimension ρ such that the limit

g(t, C 0 ) ct = lim h→0 ḡ( T-h )δ h (t) = lim h→0 ∞ k=0 1 h g τ k h , C 0 h ρ δ h (t -hk)h = ∞ 0 g(τ, C 0 ) ct δ(t -τ )dτ = lim h→0 ∞ k=0 g (k, C 0 h ρ ) 1 h δ t h ,k = lim h→0 1 h g t h , C 0 h ρ (168)
exists. In the second line we used δ h (tkh) = 1 h δ t h ,k of (164). We call g(t) ct = g(t, C 0 ) ct (t ∈ R + ) the continuous-time limit density of the discrete-time density g(t, C) (t ∈ N 0 ). We refer 1 h g( t h , C) (t = {0, h, 2h, . . . , } ∈ hN 0 ) to as the 'scaled' discrete-time density and 1 h g( t h , C 0 h ρ ) the 'well-scaled' discrete-time density, the latter with suitably chosen scaling dimension ρ such that the continuous-time limit exists. We further notice that C Limiting relation (168) indeed is useful when we consider processes on discrete times ∈ hN 0 and especially for the transition to the continuous-time limit h → 0 where the asymptotic behavior of g(k) with k = t h large comes into play. To see this we consider the important case of a fat-tailed behavior g(k) = g(k, C) ∼ Ck -α-1 for large k with α ∈ (0, 1) such as for instance in relation (59). Then we need to rescale the internal parameter C(h) = C 0 h ρ such that this limit

g(t, C 0 ) ct = lim h→0 = h -1 g t h , C 0 h ρ lim h→0 h -1 C 0 h ρ t h -α-1 = lim h→0 C 0 h ρ+α t -α-1 , t ∈ R (169)
exists. Clearly this limit exists if and only if ρ + α = 0 and hence ρ = -α thus the required scaling is C(h) = C 0 h -α . Then the continuous-time limit yields

g(t, C 0 ) ct = lim h→0 h -1 g t h , C 0 h -α = C 0 t -α-1 . ( 170
)
We notice that a fat-tailed continuous-time limit of the form (170) originates from a discretetime density having generating function ḡ

(u) = 1-C ′ (1-u) α +0(1-u) α (with C ′ = C Γ(1-α) α > 0)
which corresponds to a rescaled version of Sibuya(α). This underlines the utmost universal importance of Sibuya(α) which is recalled in more details in Appendices A.4-A.5. Furthermore, to establish the connection to relation (59): When we identify the constants

C(h) = αν ξ(h)Γ(1-α) (and C 0 = ν ξ 0 α Γ(1-α)
) then relations (169), (170) recover (59). Well-scaled limiting procedures in the sense of relations (168)-(170) are recurrently performed in the present paper.

Continuous-time limit of 'cumulative discrete-time distributions'

Apart from densities we also deal with 'cumulative discrete-time distributions' which we define as

G(t) = G(t) h=1 = d t du t ḡ(u) 1 -u | u=0 = t k=0 g(t -k)Θ(k) = Θ(t) ⋆ g(t) = 1 1 -T-1 g(t) = t k=0 g(k), t ∈ N 0 (171)
where g(t) = g(t, C) is a causal discrete-time density as defined above in (167) 34

G(t) h = h 1 -T-h ḡ( T-h )δ h (t) = ḡ( T-h ) h 1 -T-h δ h (t) = ḡ( T-h )Θ(t), t ∈ hZ 0 (172)
where 164)) and G(t) h=1 recovers (171). Now in order to consider the continuous-time limit we rewrite (172) as

h 1-T-h δ h (t) = Θ(t) (See (

G(t)

h = h 1 -T-h g(t) h = ∞ k=0 hg(t -kh) h = ∞ k=0 g( t -kh h , C) = t h k=0 g(k, C) = G t h , C t ∈ hN 0 (173) where G(k, C) = G(k, C) 1 (k = t h ∈ N 0 )
is the cumulative distribution (171) and g(t) h = ḡ( T-h )δ h (t) indicates a causal scaled density as in (167). The continuous-time limit then is obtained as

G(t, C 0 ) ct = lim h→0 G(t) h = lim h→0 G t h , C(h) = C 0 h ρ 1 = D -1 t g(t, C 0 ) ct = t 0 g(τ, C 0 ) ct dτ t ∈ R + . ( 174
)
The scaling of the constant(s) C(h) = C 0 h ρ has to be suitably chosen such that the limit (174) exists where the same scaling dimension ρ occurs as for the density in (168). This equation is constructive for the explicit evaluation of continuous-time limits for instance in the state probabilities since it is the limit of the integer time expression G( t h , C 0 h ρ ) 1 ( t h ∈ N large) of Eq. (171) for h → 0.

Further properties of causal discrete-time distributions

Now let us consider discrete-time convolutions of causal distributions G(t) = Θ(t)g(t). By 34 Causality makes in the sum

1 1-T-1 g(t) = ∞ k=0 T-k Θ(t)g(t) = ∞ k=0 Θ(t -k)g(t -k) the upper limit k = t. accounting for the shift operator T-1 such that T-n g(t) = g(t -n), especially we define T-n δ i,j = δ i-n,j = δ i,j+n = δ i,j T+n Tn = ( T-n ) † = ( T-n ) -1 (175)
where with  † we indicate the adjoint (hermitian conjugate) operator to Â. It is well-known that shift operators are unitary and represent a commutative (Abelian) group. From observation (175) we introduce the notation

Ta δ i,j = δ i+a,j = δ i,j-a = δ i,j T-a , δ i,j Tb = δ i,j+b = δ i-b,j = T-b δ i,j . (176)
We introduce the convention that if we write the shift operator T to the left of the Kroneckersymbol δ i,j it acts on the left index (i, ..) and if we write it to the right of the Kroneckersymbol, then it acts on right index (.., j). Let us consider now a chain of discrete-time convolutions

(G 1 ⋆ G 2 ⋆ . . . ⋆ G n )(t) = ḡ1 ( T-1 )ḡ 2 ( T-1 ) . . . ḡn ( T-1 ) δ t,0 = ∞ t 1 =-∞ ∞ t 2 =-∞ . . . ∞ tn=-∞ δ t,t 1 +t 2 +...+tn Θ(t 1 )Θ(t 2 ) . . . Θ(t n )g(t 1 )g(t 2 ) . . . g(t n ), t ∈ Z 0 (177)
where the first line establishes the shift-operator representation. We notice in view of ḡi ( T-1 )ḡ j ( T-1 ) = ḡj ( T-1 )ḡ i ( T-1 ) that discrete-time convolutions (as their continuous-time counterparts) indeed are commuting. The convolution power (160) is a special case of (177). The generating function of (177) can be recovered by

∞ t=0 (G 1 ⋆ G 2 ⋆ . . . ⋆ G n )(t)u t = ∞ t 1 =-∞ ∞ t 2 =-∞ . . . ∞ tn=-∞ Θ(t 1 )Θ(t 2 ) . . . Θ(t n )g(t 1 )g(t 2 ) . . . g(t n ) × ∞ t=0 δ t,t 1 +t 2 +...+tn u t = n j=1 ∞ t j =0 g(t j )u t j = ḡ1 (u)ḡ 2 (u) . . . ḡn (u) (178)
reflecting also shift-operator representation in (177). Then define operator functions A( T-1 ), B( T+1 ) which commute and act on a convolution of the form ( 178)

A( T-1 )B( T+1 )(G 1 ⋆ G 2 ⋆ . . . ⋆ G n )(t) = ∞ t 2 =-∞ . . . ∞ tn=-∞ Θ(t 1 )Θ(t 2 ) . . . Θ(t n )g(t 1 )g(t 2 ) . . . g(t n )A( T-1 )B( T+1 )δ t,t 1 +t 2 +...+tn ∞ t 2 =-∞ . . . ∞ tn=-∞ Θ(t 1 )Θ(t 2 ) . . . Θ(t n )g(t 1 )g(t 2 ) . . . g(t n )δ t,t 1 +t 2 +...+tn A( T+1 )B( T-1 ) = (G 1 ⋆ G 2 ⋆ . . . ⋆ G n )(t)B( T-1 )A( T+1 ) (179)
where it is crucial here that t ∈ Z 0 . All the representations in ( 179

T-n (G 1 ⋆ G 2 ⋆ . . . ⋆ G n )(t) = (G 1 ⋆ G 2 ⋆ . . . ⋆ G n )(t -n) = ∞ t 2 =-∞ . . . ∞ tn=-∞ Θ(t 1 )Θ(t 2 ) . . . Θ(t n )g(t 1 )g(t 2 ) . . . g(t n )δ t-n,t 1 +t 2 +...+tn = Θ(t 1 )Θ(t 2 ) . . . Θ(t n )g(t 1 )g(t 2 ) . . . g(t n )δ t,t 1 +t 2 +...+tn+n = (G 1 ⋆ G 2 ⋆ . . . ⋆ G n )(t) T+n . ( 180
)
Application of any shift-operator functions on discrete-time convolutions of causal distributions is hence well-defined in the sense (179) by reducing the shift operations to shift operations in the argument of the Kronecker-δ. In the distributional sense all operations involving operator functions of shift operators act on convolutions (for instance the discretetime fractional Kolmogorov-Feller equations ( 65), ( 79) and ( 80)) and are well-defined as they can be reduced to shift operations acting on Kronecker-δs defined on Z 0 . It is straightforward to rewrite all these relations for their scaled counterparts defined on hZ 0 involving discrete-time δ-distributions (164) generalizing the Kronecker-δs. We can conceive the causal discrete-time distributions in a wider sense as discrete-time counterparts of Gelfand-Shilov generalized functions and distributions [START_REF] Gel'fand | Generalized Functions[END_REF].

A.3 GENERALIZED FRACTIONAL INTEGRALS AND DERIVATIVES

Here we consider some properties of discrete-time convolution operators and their continuoustime limits in order to establish the connection with fractional calculus. To this end let us introduce the notion of 'generalized discrete-time fractional integrals and derivatives' in the spirit of Ref. [START_REF] Michelitsch | Sur une généralisation de l'opérateur fractionnaire[END_REF]. Let Θ(t)f (t) be a causal discrete-time distribution defined on t ∈ hZ 0 (h > 0). Then we define the well-scaled 'discrete-time integral' as

h 1 -T-h Θ(t)f (t) = h ∞ k=0 T-hk Θ(t)f (t) = t h k=0 hf (t -kh) = t h k=0 hf (kh), t ∈ hZ 0 (181)
which is well-scaled thus its continuous-time limit h → 0 exists for sufficiently 'good', i.e. absolutely integrable distributions f (t) defined on t ∈ R + . The continuous-time limit then yields the standard integral

lim h→0 h (1 -e -hD t ) Θ(t)f (t) = D -1 t Θ(t)f (t) = t 0 f (τ )dτ, t ∈ R + . ( 182
)
Now we introduce 'generalized derivatives and integrals' which include traditional and fractional derivatives and integrals. We generalize (181) by introducing the well-scaled 'discrete-time fractional integral' as

h (1 -T-h ) µ • Θ(t)f (t) = h µ ∞ k=0 (-1) k -µ k f (t -hk)Θ(t -kh) = h µ t h k=0 (-1) k -µ k f (t -hk) µ > 0 , t ∈ hN 0 . ( 183
)
Taking into account that the shift operator can be represented as

T-h = e -hD t , D t = d dt ( 184
)
where lim h→0 h -1 (1 -T-h ) • f (t) denoting the traditional derivative of first order. Then we consider the limit

lim h→0 h µ (1 -e -hD t ) -µ f (t) = D -µ t • f (t) (185)
which defines first of all an operator representation of the integral of order µ > 0. Now we evaluate explicitly the continuous-time limit of (183) to arrive at

D -µ t • f (t) := lim h→0 h µ (1 -e -hDt ) -µ • f (t) = lim h→0 t h k=0 h µ Γ(µ) Γ(µ + k) Γ(k + 1) f (t -kh) = ⌈th -δ ⌉ k=0 h µ Γ(µ) Γ(µ + k) Γ(k + 1) f (t -kh) + th -1 k=1+⌈th -δ ⌉ h µ Γ(µ) Γ(µ + k) Γ(k + 1) f (t -kh) = lim h→0 t h k=1+⌈th -δ ⌉ h (kh) µ-1 Γ(µ) f (t -kh) + O(h (1-δ)µ+2δ ) = lim h→0 t h(1+⌈th -δ ⌉) τ µ-1 Γ(µ) f (t -τ )dτ + O(h (1-δ)µ+2δ ) = t 0 τ µ-1 Γ(µ) f (t -τ )dτ.
δ ∈ (0 , 1 ) (186)

We indeed identify the last line with the Riemann-Liouville fractional integral of order µ (µ > 0) [START_REF] Samko | Fractional Integrals and Derivatives: Theory and Applications[END_REF] where µ ∈ N recovers traditional integer order integration and for µ → 0+ the kernel τ µ-1 Γ(µ) → δ(τ ) thus reproducing the identity D 0+ t f (t) = f (t). For µ > 0 the kernel τ µ-1

Γ(µ)
is on R + a locally integrable function.

In the deduction (186) we employed the ceiling function 35 ⌈..⌉ namely k 0 = ⌈th -δ ⌉ where we have lim h→0 h(1 + ⌈th -δ ⌉) ∼ th 1-δ → 0. We choose δ ∈ (0, 1) such that the first sum S 1 (h) in the second line tends to zero for h → 0 at least as S 1 (h) ∼ h µ(1-δ)+2δ which is shown below.

The crucial point for this choice is that hk 0 (h) ∼ th 1-δ → 0, however k 0 (h) ∼ th -δ → ∞ for h → 0. Hence in the second sum of the second line in all terms k is so can employ the asymptotic relation Γ(k+µ) Γ(µ)Γ(k+1) ∼ k µ-1 Γ(µ) . The scaling behavior of S 1 (h) for h → 0 is then obtained as follows

S 1 (h) = k 0 k=0 h µ Γ(µ) Γ(µ + k) Γ(k + 1) f (t -kh) ∼ f (t -ηk 0 h)h µ k 0 k=0 -µ k (-1) k , k 0 (h) = ⌈th -δ ⌉ = f (t -ηhk 0 )h µ 1 k 0 ! d k 0 du k 0 (1 -u) 1-µ u=0 = f (t -ηhk 0 )h µ (-1) k 0 1 -µ k 0 , η ∈ (0, 1) = f (t -ηhk 0 ) h µ Γ(µ + k 0 -1) Γ(k 0 + 1)Γ(µ -1) = f (t -ηhk 0 )h µ (µ -1) k 0 k 0 ! S 1 (h) ∼ f (t -ηhk 0 ) h µ k µ-2 0 Γ(µ -1) ∼ f (t) t µ-2 Γ(µ -1) h µ h -δ(µ-2) ∼ h µ(1-δ)+2δ → 0, (h → 0) (187)
where in (187) 36 µ = 1. In the last line we used that k 0 = ⌈th -δ ⌉ ∼ th -δ (large) with δ ∈ (0, 1) and µ > 0 and the interim value f (t-ηhk 0 ) ∼ f (t) is approximately constant over the (small) interval [tk 0 h, t] especially for h → 0.

A more sophisticated generalized integral operator which we refer to as 'Mittag-Leffler integral' is of interest [START_REF] Kilbas | Generalized Mittag-Leffler Function and Generalized Fractional Calculus Operators[END_REF]. We can define this convolution operator by the following wellscaled continuous-time limiting procedure

I α • f (t) = lim h→0 ξ(h) ξ(h) + (1 -T-h ) µ
• f (t), µ > 0 , t ∈ hN 0 where we made use of the limiting property lim h→0 L{δ h (t)}(s) = L{δ(t)}(s) = 1 of (166).

We identify (189) with the Laplace transform of the Mittag-Leffler density. Now rewrite (188) as a series of Riemann-Liouville fractional integral operators

I µ • δ(t) = lim h→0 ∞ n=1 (-1) n-1 ξ 0 h µ (1 -e -hDt ) -µ n
• δ h (t) 36 The case µ = 1 of standard integration is considered above (181)-(182).

I µ • δ(t) = ∞ n=1 (-1) n-1 ξ n 0 D -nµ t • δ(t) = ξ 0 ξ 0 + D µ t • δ(t) µ > 0 , t ∈ R + (190)
where the Riemann-Liouville integrals D -nµ t • δ(t) are defined by (186). Thus we obtain the continuous-time limit convolution kernel

I µ • δ(t) = ξ 0 ξ 0 + D µ t • δ(t) = t 0 ∞ n=1 (-1) n-1 ξ n 0 τ nµ-1 Γ(nµ) δ(t -τ )dτ = ξ 0 t µ-1 E µ,µ (-ξ 0 t µ ) = d dt (1 -E µ (ξ 0 t µ )) t ∈ R + (191)
which we identify with the Mittag-Leffler density where here we have µ > 0 and µ = 1 recovers the exponential density I 1 • δ(t) = ξ 0 e -ξ 0 t .

So far in this appendix we generalized well-scaled discrete-time fractional integrals which involve uniquely non-negative powers of the well-scaled discrete-time integral operator

h 1-T-h
. Now we analyze operators involving its negative powers leading to the notion of generalized fractional derivatives. To this end we define the well-scaled 'discrete-time fractional derivative' as37

h 1 -T-h -µ • f (t) = h -µ (1 -T-h ) µ • Θ(t)f (t) = h -µ t h k=0 µ k (-1) k f (t -kh) µ > 0 , t ∈ hN 0 . ( 192
)
Then, recall the ceiling function ⌈µ⌉ indicates the smallest integer larger or equal to µ. We have -1 < µ -⌈µ⌉ ≤ 0 with -1 < µ -⌈µ⌉ < 0 if µ / ∈ N and zero if µ = ⌈µ⌉ is integer. Then

converges by taking into account (186) to the Riemann-Liouville fractional integral of order ⌈µ⌉µ. Hence the continuous-time limit of (192) yields

D µ t • f (t) = lim h→0 h -µ (1 -T-h ) µ • f (t) = D ⌈µ⌉ t t 0 f (t -τ ) τ ⌈µ⌉-µ-1 Γ (⌈µ⌉ -µ) dτ, t ∈ R + (196)
which we identify with the Riemann-Liouville fractional derivative or order µ [START_REF] Samko | Fractional Integrals and Derivatives: Theory and Applications[END_REF] where τ ⌈µ⌉-µ-1 due to -1 < ⌈µ⌉µ -1 < 0 is if µ / ∈ N 0 a weakly singular integrable function. For ⌈µ⌉µ → 0, i.e. in the limit of µ approaching integer we have the distributional relation τ ⌈µ⌉-µ-1 Γ(⌈µ⌉-µ) → δ(τ ) [START_REF] Gel'fand | Generalized Functions[END_REF] thus (196) converges to the integer order derivative D µ t •f (t) → D ⌈µ⌉ t

• f (t).

A.4 SIBUYA WALK

We consider now a strictly increasing random walk with Z j ∈ N are IID steps occurring with the probabilities P(Z j = r) = w α (r) = (-1) r-1 α r = (-1) r-1 α(α -1)..(αr + 1) r! , r ∈ N, w α (r) r=0 = 0, α ∈ (0, 1). (197) This distribution is referred to as Sibuya(α) (See [START_REF] Pachon | On Discrete-Time Semi-Markov processes[END_REF] for a thorough analysis of properties).

The limit α → 1 -0 gives a 'trivial walk' where the walker in each step hops from r to r + 1 (almost surely) with transition probabilities w 1 (r) = δ r1 (r ∈ N 0 ). Despite α = 1 is admissible, per definition Sibuya(α) covers only the fractional interval α ∈ (0, 1).

which is fat-tailed with the same tail as the Mittag-Leffler density.

An interesting quantity in a strictly increasing walk is the expected hitting number τ (r) of a node r. For α = 1 this quantity is trivial and yields τ α=1 (r) = 1 which is the upper bound in a strictly increasing walk. Since (197) allows long-range jumps (with overleaping of nodes) and in view that the walk is strictly increasing so that nodes can be visited either once or never, intuitively we infer that the expected hitting number fulfills 0 < τ α (r) < 1 in the fractional interval α ∈ (0, 1). This will be indeed confirmed subsequently by means of explicit formulas. The expected hitting number τ α (r) indicates the average number of visits of a node r and generally is an important quantity to describe recurrence phenomena [START_REF] Michelitsch | Recurrence of random walks with long-range steps generated by fractional Laplacian matrices on regular networks and simple cubic lattices[END_REF][START_REF] Michelitsch | On Recurrence and Transience of Fractional RandomWalks in Lattices[END_REF][START_REF] Michelitsch | Fractional Dynamics on Networks and Lattices[END_REF]. It is clear that a strictly increasing walk always is transient, since the walker cannot return to a node. Generally the expected hitting number of a node in a strictly increasing walk fulfills 0 < τ (r) ≤ 1. By simple conditioning arguments we have that 38 (201) 38 For an analysis of this quantity for space-fractional Markovian walks on undirected networks, see [START_REF] Michelitsch | Fractional Dynamics on Networks and Lattices[END_REF].

by v → H. Now its is only a small step to write down the Kolmogorov-Feller type equations for the Sibuya DTRW, namely Let us also consider a walk not taking place on an undirected network, but the simplest variant of strictly increasing walk where the walker at each event makes a unit step into the positive integer direction. This walk simply counts the arrivals and has the transition probabilities H ij = δ i+1,j i.e. we have then in (214) (HP S (t)) ij = ∞ k=0 δ i+1,k P kj = P i+1,j .

For this walk Eq. ( 214) then takes the form

(1 -T-1 ) α P i+1 ,j (t) -(-1 ) t α -1 t δ ij = P i+1 ,j (t) -P ij (t)

(215)
where this equation is also obtained in [START_REF] Pachon | On Discrete-Time Semi-Markov processes[END_REF] (Remark 7 in that paper).

Figure 1 :

 1 Figure 1: Discrete-time densities ϕ (ν) α (t, ξ) and ψ (ν)

Figure 2 :

 2 Figure 2: The plots show discrete-time waiting-time density χ α,ν (t) h versus h for different times t. (a) α = 0.5 and ν = 0.5, (b) α = 0.574 and ν = 1.742. The colorbar represents the values of time t. The smaller h the more the GFPP Prabhakar waiting time density (2) is approached (See Eq. (20)).

α

  (t)} just shifted by one time unit into positive timedirection. This shows the (shift)-operator representation θ (ν)

Figure 3 :

 3 Figure 3: State probabilities Φ (n) α,ν (t) ct of the GFPP representing the continuous-time limit of the PDTP state probabilities versus t for different values of n. (a) α = 0.5 and ν = 0.5, (b) α = 0.574 and ν = 1.742. In the colorbar we represent the values n = 1, 2, . . . , 8. The results were obtained numerically using ξ 0 = 1 with Eq. (52).

  α,ν (t)ct there differs by a multiplier ξ ν 0 .

  . Further |v m denote the right-and vm | the left-eigenvectors of the generally non-symmetric transition matrix H. We have the properties vn |v m = δ mn with the N × N unity matrix 1 = N n=1 |v n vn |.

  Now let us see the connection between generating functions and the shift operator T-1 which is defined by T-1 g(t) = g(t-1). Consider again generating function (156) and replace u by the shift operator T-1 and use ( T-1 ) n δ t,0 = T-n δ t0 = δ (t-n),0 = δ tn to arrive at

  scale in the continuous-time density g(t, C 0 ) ct .

  ) are equivalent where we utilized (176) with the adjoint operator function A( T-1 )B( T+1 ) † = B( T-1 )A( T+1 ) (commuting). For instance consider A( T-1 ) = T-n , then we have

  h) = ξ 0 h µ is chosen such that (188) exists. Then (188) has the Laplace transformlim h→0 L{I µ • δ h (t)}(s) = lim h→0 ∞ 0 - e -st ξ 0 h µ ξ 0 h µ + (1e -hDt ) µ δ h (t)dt = lim h→0 ξ 0 ξ 0 + h -µ (1e -hs ) µ L{δ h (t)}(s) = ξ 0 ξ 0 + s µ

( 192 )

 192 can be decomposed as followsh -µ (1 -T-h ) µ = h -⌈µ⌉ (1 -T-h ) ⌈µ⌉ h ⌈µ⌉-µ (1 -T-h ) µ-⌈µ⌉ (193)which is a discrete-time derivative of integer-order ⌈µ⌉ applied to a (generalized fractional discrete-time) integral of order ⌈µ⌉µ. The continuous-time limit of (192) is then obtained by accounting forlim h→0 h -⌈µ⌉ (1 -T-h ) ⌈µ⌉ = d ⌈µ⌉ dt ⌈µ⌉ = D ⌈µ⌉ t (194)converges to the integer-order traditional derivative of order ⌈µ⌉ andlim h→0 h ⌈µ⌉-µ (1 -T-h ) µ-⌈µ⌉ • f (t)

  Sibuya(α) has the generating functionwα (u) = ∞ k=1 (-1 ) k-1 α k u k = 1 -(1u) α , |u| ≤ 1 (198) with w α (r) = 1 r! d rdu r w(u)| u=0 gives the probabilities (197). We notice that for large r using the asymptotic formulaΓ(r+c) Γ(r) ∼ r c we have w α (r) = (-1) r-1 α r = α r Γ(rα) Γ(1α)Γ(r) ∼ αr -α-1 Γ(1α) = -r -α-1 Γ(-α), (r large), α ∈ (0, 1)

  α ⋆) n (r) indicate convolution powers (See Appendix A.1). It follows that (200) for Sibuya(α) has the generating functionG α (u) = ∞ n=0 (w α (u)) n = 1 1wα (u) = (1u) -α , |u| < 1

( 1 -( 1 -

 11 T-1 ) α HP S (t) = (H -1)P S (t) + 1(1 -T-1 ) α-1 δ t0 (213) or T-1 ) α HP S (t) = (H -1)P S (t) + 1(-1 ) t α -1 t . (214)This equation is the Kolmogorov-Feller fractional difference equation solved by the transition matrix P S (t) in a Montroll-Weiss Sibuya DTRW.

  then the discrete-time counting process is non-Markovian and has long-time power-law memory as in (61) with fat-tailed waitingtime density (59). If in contrast the state probability generating functions (and memory generating function) at u = 1 are finite, then the discrete-time counting process either is memoryless with the Markovian property or is non-Markovian and has only a short-time memory with a light-tailed waiting time density.

The analysis to follow can be extended to any f -distribution of the previous section.

For α = 1 and ν > 0, (30) coincides with generating function[START_REF] Gorenflo | On the Fractional Poisson Process and the Discretized Stable Subordinator[END_REF].

and also in[START_REF] Kilbas | Generalized Mittag-Leffler Function and Generalized Fractional Calculus Operators[END_REF] for ν = 1, α = 1, consult also[START_REF] Michelitsch | Continuous time random walk and diffusion with generalized fractional Poisson process[END_REF][START_REF] Michelitsch | Generalized fractional Poisson process and related stochastic dynamics[END_REF]

Involving kernels which are at most weakly singular and hence integrable.

We call an operator a 'generalized fractional integral' if it contains only non-negative powers of the discrete integral operator h(1 -T-h ) -1 , see especially Appendix A.3 and consult also[START_REF] Michelitsch | Sur une généralisation de l'opérateur fractionnaire[END_REF].

In the regime α = 1 we have the finite value Ḡ1,ν(u, v)|u=1 = 1 1-v (1 + ν ξ ).

This result can be reconfirmed by n 1,1(t)ct = e -ξ 0 t ∞ n=0 n (ξ 0 t) n n!

For a proof, see[START_REF] Michelitsch | Fractional Dynamics on Networks and Lattices[END_REF]

In this notation S = |a b| stands for the N × N matrix which has the elements Sij = i|a b|j and i|c = ( c|i ) * where (..) * indicates complex conjugation.

For α → 1 -0 we get t -α Γ(1-α) → δ(t) = 0.

Be reminded that with our initial condition P(0) = 1 the matrices on the right-hand side commute.

We utilize here notation P α,1 (t) h with subscript (. . .) h when we refer to the time scaled walk with t ∈ hZ0.

We use sometimes also the synonymous notation g(t) ⋆ h(t) for (159).

Note that g( t h , C) is dimensionless since t h is dimensionless as t and h have units sec.

For a definition see below(192).

The limit h → 0 of this relation can be identified with the Grünwald-Letnikov derivative, e.g.[START_REF] Meerschaert | Stochastic Models for Fractional Calculus[END_REF].

where wα (u) is the generating function (198) of the Sibuya-steps. The expected hitting number of a node r ∈ N 0 in a Sibuya walk is then obtained as

where (α) r stands for the Pochhammer-symbol [START_REF] Metzler | The Random Walk's Guide to Anomalous Diffusion : A Fractional Dynamics Approach[END_REF]. For α = 1 we have τ 1 (r) = r! r!Γ(1) = 1 ∀r ∈ N 0 which is the above anticipated result when the walker hops in each step from r to r + 1, thus the walker hits each node once (almost surely). Hence τ 1 (r) = 1 is the upper bound for τ α (r) which can be seen accounting for

where each multiplier fulfills 0

The lower bound τ 0+ (r) ∼ δ t0 (i.e. null for r > 0) is approached for α → 0+ where also for r large we then have τ α (r) ∼ r α-1 Γ(α) ∼ δ(t) = 0. The physical interpretation is that the steps then become infinitely long, thus any finite node r > 0 is over-leaped where only the departure node r = 0 is occupied once.

On the other hand for α → 1 -0 we have τ 1-0 (k) = 1 thus each node is hit once (almost surely) since the walker moves in any time step from r to its right neighbor node r + 1 (with transition probabilities w 1 (r) = δ r1 ).

It is instructive to recall briefly the limiting distribution of the rescaled integer random variable σ α (n) → λ n σ α (n) for infinitely many steps n → ∞ where σ(n) = n j=1 Z j is an integer random variable with IID Sibuya steps Z j . Then let us consider the Laplace transform

In this limiting relation we rescale the Sibuya random variable σ(n) = n j=1 Z j → λ n σ(n) such that the limit (204) exists. This requires the scaling factor λ n → 0 and 1-

if and only if we choose λ n = ( ξ n )

1

where ξ is an arbitrary positive constant independent on n. With this choice we get

Equation ( 206) is the Laplace transform of the limiting distribution emerging for many Sibuya steps (n → ∞) of the rescaled random variable ξ

subordinator (see e.g. [START_REF] Pachon | On Discrete-Time Semi-Markov processes[END_REF] and the references therein). 39 i.e. neither null nor infinity.

A.5 SIBUYA-COUNTING PROCESS AND DTRW

In this appendix we review some pertinent results for the 'Sibuya counting process' (i.e. the discrete-time counting process with waiting times following Sibuya(α)) and develop the Sibuya Montroll-Weiss DTRW as a walk subordinated to the Sibuya counting process. The interested reader is referred to [START_REF] Pachon | On Discrete-Time Semi-Markov processes[END_REF] for a profound analysis, and see also the references therein. The generating function of the Sibuya survival probability is then with (198) obtained as

which yields the survival probability

S,α (t) t=0 = 1. We observe that for α → 1 -0 the survival probability Φ S,1 (t) = δ t0 for t > 0 is vanishing since for α = 1 the waiting time is equal to one (almost surely with w 1 (t) = δ t1 ). Now consider the behavior for t large which yields

with the same tail as has the Mittag-Leffler waiting time distribution. Hence for large observation times the Sibuya counting process behaves as a fractional Poisson process and the same is true for the PDTP (See asymptotic relations (61) for the PDTP state probabilities.). We mention the state-probabilities, i.e. probabilities for n arrivals in the discrete time interval [0, t] are given by

with generating function

) n and refer to Ref. [START_REF] Pachon | On Discrete-Time Semi-Markov processes[END_REF] for a thorough analysis.

SIBUYA DTRW

Due to the general importance of Sibuya(α) let us briefly outline the Sibuya Montroll-Weiss DTRW where we subordinate a walk with one step-transition matrix H to a Sibuya counting process. Then we obtain for the generating function of the transition matrix the Cox-series (See (123))

We call the random walk with transition matrix generating function (211) the 'Sibuya DTRW'. Formula (211) is related to the generating function