Non-Local Solvable Birth-Death Processes

In this paper we study strong solutions of some non-local difference-differential equations linked to a class of birth-death processes arising as discrete approximations of Pearson diffusions by means of a spectral decomposition in terms of orthogonal polynomials and eigenfunctions of some non-local derivatives. Moreover, we give a stochastic representation of such solutions in terms of time-changed birth-death processes and study their invariant and their limit distribution. Finally, we describe the correlation structure of the aforementioned time-changed birth-death processes.


Introduction
Birth-death processes constitute a fundamental class of continuous-time Markov chain that are widely used in applications such as, for instance, evolutionary dynamics [37] and queueing theory [24]. In particular, one can achieve a complete characterization of birth-death processes by families of classical orthogonal polynomials of discrete variable. This theory, linked to the solution of the Stieltjes moment problem, has been widely studied by Karlin and McGregor in their seminal papers [20,23]. As birth-death processes are linked to difference-differential equations, fractionalization of such processes can be used to study the solutions of fractional differencedifferential equations. Indeed, with this idea in mind, a fractional version of the Poisson process has been introduced in [8,9] and fractional versions of some birthdeath processes, for instance, in [38,39,40]. In the case of Pearson diffusions, one can use a spectral approach to study strong solutions of fractional backward and forward Kolmogorov equations and, at the same time, define the fractional Pearson diffusions by means of a time-change via an inverse stable subordinator (see, for instance, [29,30,31]). The same approach has been used to study the case of fractional immigration-death processes in [5]. Let us also stress out that this approach can be used to study a fractional M/M/∞ queue or a fractional M/M/1 queue with acceleration of service (for some models of fractional queues, we refer to [4,6,15]). However, one could consider a time-change with a different inverse subordinator. In such case, in place of the fractional derivative in time, one obtains a more general non-local operator. Such kind of operators have been introduced in [25] for the class of complete Bernstein functions and extended in [46] for any Bernstein function. A first step towards the theory of general time-changed Pearson diffusions has been achieved in [18]. In this work, we describe a general theory for non-local solvable birth-death processes in terms of orthogonal polynomials, where such processes are defined by means of a time-change with a general inverse subordinator. In particular, we focus on the strong solutions of general non-local backward and forward Kolmogorov equations associated to such processes and on the stochastic representation of such solutions. In particular, the paper is structured as follows: • In Section 2 we introduce the theory of solvable birth-death processes as discrete approximations of Pearson diffusions and we state the main hypotheses we need on the starting birth-death process; • In Section 3 we give some preliminaries on inverse subordinators and nonlocal time derivatives. In particular we focus on the eigenvalue equation for such derivatives and on some upper bounds for the eigenfunctions. Let us stress out that some properties of such eigenfunctions are expressed in [25,26] in the complete Bernstein case and in [34] in the general case. Moreover, a series expansion in terms of convolutions of potential densities in the special Bernstein case is obtained in [3]; • In Section 4 we focus on the spectral decomposition of the strong solutions of non-local forward and backward Kolmogorov equations in terms of orthogonal polynomials of discrete variable and eigenfunctions of the non-local time derivatives; • In Section 5 we introduce the time-changed birth-death processes and we study the stochastic representation of the aforementioned strong solutions in terms of such processes. In particular we obtain that the time-changed process still admits the same invariant measure that is also the limit measure for any starting distribution; • Finally, in Section 6 we study the correlation structure of the time-changed birth-death processes in terms of the potential measure of the involved subordinator and the eigenfunctions of the non-local time derivatives. In particular, the non-stationarity of the process is evident in the expression of the covariance, thus, to give some information on the memory of the process, we have to refer to a non-stationary extension of the definition of long-range and short-range dependence suggested by the necessary conditions given in [10, Lemma 2.1 and 2.2].

Solvable Birth-Death Processes
Let us fix a filtered space (Ω, F , {F} t∈R + P) and consider a Birth-Death process {N (t), t ≥ 0} on it. Let us denote by E ⊆ Z its state space, that will be finite or at most countable. In particular we can always suppose that E ⊆ N 0 and is a segment, i.e. for any n 1 , n 2 ∈ E and n ∈ N 0 such that n 1 ≤ n ≤ n 2 it holds n ∈ E, with min E = 0. Let us recall that the generator G of a Birth-Death process can be always expressed as where d(x) are the death rates, b(x) are the birth rates (recalling that d, b must be non-negative in E), ∇ ± are the first order forward and bacwkard finite differences defined as and ∆ is the second order central finite difference Here we want to introduce some birth-death version of Pearson diffusions (see, for instance, [30]). To do this, we refer to the theory of birth-death polynomials, whose main papers are [20,23].
• the spectrum of G is purely discrete with non-positive eigenvalues (λ n ) n∈E such that λ 0 = 0 and λ n < 0 for any n ≥ 1; • its eigenfunctions (P n ) n∈E are classical orthogonal polynomials of discrete variable with orthogonality measure m which is the invariant and stationary measure of N (t); • the function m(x) = m({x}) solves the following discrete Pearson equation: is a polynomial of degree at most 2 and b(·) − d(·) is a polynomial of degree at most 1.
Concerning solvable birth-death processes, they arise as lattice approximations of Pearson diffusions. In particular, one has in such case Concerning classical orthogonal polynomials of discrete variable, we mainly refer to [36,43]. Their orthogonality relation is expressed as where δ n,m is Kronecker delta symbol. In particular one obtains that P n ℓ 2 (m) = d n and then we can introduce the normalized polynomials as Q n (x) = Pn(x) dn , such that On the other hand, the function m(n) = 1 d 2 n defines a measure on E. Thus, by proceeding with a Gram-Schmidt orthogonalization procedure on the monomials (1, x, x 2 , ·), we can define a family of orthogonal polynomials P n (x) that satisfies the following orthogonality condition The family of polynomials ( P n ) n∈E are called the dual family of (P n ) n∈E , see [36]. Let us give some examples: • Immigration-death processes (see [2]) are defined by a constant birth rate b and a linear death rate d(x) = dx. In such case the invariant measure is given by the Poisson distribution m(x) = e −ρ ρ x x! for x ∈ N, where ρ = b d . The orthogonal polynomials are Charlier polynomials of parameter ρ, that we denote by P n (x) = C n (x; ρ) and they satisfy the self-duality relation (2.2) P n (x) = P x (n) and the eigenvalues are given by λ n = −bn. In particular the polynomials P n coincide with the family of dual orthogonal polynomials P n . This kind of process arises as a lattice approximation of the Ornstein-Uhlenbeck process. • Let us consider a linear death rate d(x) = dx and a linear birth rate b(x) = (x + β)b, with b, d, β > 0 and b < d. In such case the birth death process admits state space E = N and the generator is given by Defining ρ = b d , we have that the orthogonal polynomials are Meixner polynomials of parameters ρ and β, that we denote by P n (x) = M n (x; ρ, β) and they are orthogonal with respect to the invariant measure , which is a Pascal (or negative binomial) distribution of parameters β and ρ. Also Meixner polynomials satisfy the self-duality relation (2.2) and coincide with their dual polynomials. Finally, let us observe that the eigenvalues are given by λ n = −(d − b)n. This process is called the Meixner process and arises as lattice approximation of the Cox-Ingersoll-Ross process. Meixner processes are discussed for instance in [21].
• Another case is given by a linear death rate d(x) = dx and a linear de- In such case the birth-death process admits finite state space E = {0, . . . , N } and the generator is given by Defining p = b b+d and q = 1 − p to achieve the invariant distribution given by which is a Binomial distribution over E. The orthogonal polynomials are Krawtchouk polynomials of parameters N and p, that we denote by P n (x) = K n (x; N, p) and satisfy the self-duality relation 2.2. The eigenvalues of the generator are given by λ n = −n(b + d). Let us recall that this is actually a time-continuous version of the Ehrenfest urn model (see, for instance, [22]). • Another interesting case is given by a quadratic one. Indeed let us consider for some d > 0 and α, β, N ∈ N Let us observe that b(N ) = 0, thus the state space of the process is given by E = {0, . . . , N }. Moreover, its generator is given by The invariant measure of this birth-death process is an hypergeometric distribution on E given by and the orthogonal polynomials are the Hahn polynomials, that we denote by P n (x) = H n (x; α, β, N ). In this case we do not have self-duality relation, but the family of dual Hahn polynomials, that we denote by P n (x) = R n (x; α, β, N ), is linked to the Hahn polynomials by the relation H n (x; α, β, N ) = R x (n(n + α + β + 1); α, β, N ).
This particular birth-death process is a lattice approximation of the Jacobi process. For such process, we refer directly to [43].
In the examples we have not only the lattice approximations of the light-tailed Pearson diffusions, but also another process, which is the continuous-time Ehrenfest urn process. Hence we can observe that with the definition of solvable birth-death process we do not only cover these lattice approximations, but we also gain some other birth-death processes that are not covered in the theory of light-tailed Pearson diffusions (see [30,31]). Let us observe that for a birth-death process N (t) the forward operator L is defined as . Now let us show what the orthogonal polynomials P n (x) represent for the forward operator.
Lemma 2.1. Let N (t) be a solvable birth-death process with forward operator L, invariant measure m and associated family of orthonormal polynomials Q n (x). Then we have for any n ∈ E and any x ∈ E, L z→x (m(z)Q n (z))(x) = m(x)λ n Q n (x).
Proof. We have, recalling that ∆ = ∇ − ∇ + and that ∇ − is a linear operator, )m(z)Q n (z) + ∇ + y→z (d(y)m(y)Q n (y))(z)](x). Now, by discrete Leibniz rule on ∇ + , we have . Now let us use again Leibniz rule on ∇ − to achieve . Now let us work with ∇ − z→x (d(z + 1)m(z + 1))(x). We have However, by the discrete Pearson equation (2.1) we obtain and then we have concluding the proof.
Thus we have, as a consequence of the discrete Pearson equation, the discrete version of the spectral decomposition for parabolic problems with the generator and the forward operator of a light-tailed Pearson diffusion: Theorem 2.2. Let N (t) be a solvable birth-death process with state space E, generator G, forward operator L, invariant measure m and family of associated orthonormal polynomials (Q n ) n∈E . Then the following assertions hold true: • The transition probability function p(t, x; y) = P(N (t + s) = x|N (s) = y) for x, y ∈ E and t, s ≥ 0 admit the following spectral representation: for any x, y ∈ E and t ≥ 0. • If g ∈ ℓ 2 (m) with g(x) = n∈E g n Q n (x) then the strong solution of the Cauchy problem du dt (t, y) = G u(t, y) t ≥ 0, y ∈ E u(0, y) = g(y) y ∈ E is given by u(t, y) = n∈E g n e λnt Q n (y) = x∈E p(t, x; y)g(x).
In particular p(t, x; y) is the fundamental solution of du dt (t, y) = G u(t, y) and u admits the stochastic interpretation: x ∈ E is given by x; y)f (y).
In particular p(t, x; y) is the fundamental solution of dv dt (t, x) = L v(t, x) and, if f ≥ 0 with f ℓ 1 = 1, v admits the stochastic interpretation: where P f is the probability measure obtained by conditioning with respect to the fact that N (0) admits distribution f .
From this theorem, it is easy to determine the covariance of any solvable birthdeath process in its stationary form. First of all, let us observe that the stationary version of N (t) admits moments of any order. This is obvious if E is finite. To show this when E = N 0 , let us first show the following Proposition. Proposition 2.3. Let N (t) be a solvable birth-death process with invariant measure m and state space E = N 0 . Then there exists a constant ρ < 1 and a state x 0 ∈ E such that for any x ≥ x 0 it holds Proof. First of all, let us observe that since b(x) and d(x) are polynomials, then lim x→+∞ b(x) d(x+1) always exists. Moreover, we can rewrite the discrete Pearson equation as

Thus we have that m is well defined if and only if
It is easy to see that such condition implies lim x→+∞ This could happen only if b(x) and d(x) are polynomials of the same degree and with the same director coefficient. However, if b(x) and d(x) are polynomials of degree at most 1, then λ n = 0 for any n ≥ 1, that is absurd. Thus we have that b(x) and d(x) are polynomials of degree 2. However, since λ n < 0 for any n ≥ 1, it follows that the coefficient director of d(x) must be negative. However, this means that for x big enough it holds d(x) < 0, which is absurd. Thus we conclude that Now let us consider ρ ∈ (l, 1). Then there exists a state for any x ≥ x 0 . Finally, the assertion follows from the previous inequality by induction.
As a direct consequence of the previous proposition we obtain Corollary 2.4. Let N (t) be a solvable birth-death process with invariant distribution m such that N (0) admits m as distribution. Then N (t) admits moments of any order. Now we can focus on the autocovariance of the process N (t).
Corollary 2.5. Let N (t) be a solvable birth-death process with invariant measure m. Then there exists a constant a 1 ∈ R such that, for any t, s ≥ 0 Proof. First of all let us recall that the stationary version of N (t) admits moments of any order, thus in particular also second order moments and then the autocovariance is well-defined. Since we are supposing that N (t) is stationary, we have, for any t ≥ s, Thus let us consider t ≥ 0 and let us evaluate Cov m (N (t), N (0)). To do this, let us rewrite . Since m admits second moment then ι(x) = x is in ℓ 2 (m). Moreover, since deg(ι(x)) = 1, it can be written as a linear combination of Q 0 = 1 and Q 1 . Let then By the previous theorem we have that (recalling that λ 0 = 0 for any solvable birth-death process). Starting from this observation, we have that As we stated before, we can write x = a 0 + a 1 Q 1 (x), thus we have By using the orthonormality relation we have We finally achieve Cov m (N (t), N (0)) = a 2 1 e λ1t concluding the proof.

Classification of solvable birth-death processes.
We can actually improve the result in Proposition 2.3 by obtaining a complete classification of solvable birth-death processes. Indeed we have the following Proposition.
Proposition 2.6. Let N (t) be a solvable birth-death process with state space E. Then one of the following statements holds true: is a Meixner process. In particular, if (P n ) n∈E is the family of orthogonal polynomials associated to N (t), then either E is finite or (P n ) n∈E coincide with its dual family and P n (x) = P x (n) for any x, n ∈ E.
In particular this implies that deg b ≤ deg d. If deg d = 0, then also deg b = 0 and in such case λ n = 0 for any n ∈ N, which is absurd. Thus deg d = 1, 2. However, we have already seen that if deg d = 2, then, since λ n < 0 for any n ≥ 1, the director coefficient of d must be negative and this is a contradiction with the fact that d is non negative on E. Thus we conclude that if E = N 0 , then deg d = 1. Moreover, arguing as before, we know that the director coefficient of d must be positive (since if it is negative then d is negative for big values of x ∈ E) and, being also d(0) = 0, it must hold d(x) = dx for some d > 0. Now let us consider b. Since we want E = N 0 , arguing as we did with d, we need the director coefficient of b to be positive. Hence we have two possibilities: • deg b = 0, thus b > 0 is constant and we get an immigration-death process; for some b, β > 0 and then we get a Meixner process (since also b < d by the condition lim x→+∞ b(x+β) . This concludes the proof.
As we can see from the previous Proposition, we already discussed as examples the unique two cases in which the state space is countably infinite, that are actually the ones for which the proof of the main results (that will follow) are more articulated. Let us also observe some other particular properties concerning the classification of solvable birth-death processes: • In the class of solvable birth-death processes we have lattice approximations of Pearson diffusions of first spectral category. Pearson diffusions are statistically tractable diffusions (see [17]), however, their spectral behaviour can be distinguished in three different classes (see [31]). The first spectral class (the one that contains Pearson diffusions whose generators admit purely discrete spectrum) is composed of the Ornstein-Uhlenbeck, Cox-Ingersoll-Ross and Jacobi diffusions. These diffusions are approximated respectively by immigration-death, Meixner and Hahn birth-death processes; • However, we do not obtain the analogous of Pearson diffusions of the second spectral category: this is due to the fact that to obtain these, we need d(x) to be a polynomial of degree 2 and with positive director coefficient, which goes in contradiction with the request that the spectrum of the generator of a solvable birth-death process is purely discrete and non-positive together with the fact that the state space is infinite. However, if we reduce the state space to a finite segment of N 0 , we can still cover these cases, in which d(x) admits discriminant ∆ ≥ 0 and d(0) = 0, by asking that b(x) admits the first root x 0 ∈ N and E = {0, . . . , x 0 }. If ∆ > 0 we also have to ask that the other root of d(x) is negative. In any case, these lattice schemes do not approximate reciprocal Gamma and Fisher-Snedecor diffusions in the support of their invariant measure, but they still provide, in some cases, an approximation scheme for the backward and forward Kolmogorov equations in a subset of the full support; • Even considering a finite segment in N 0 , we cannot approximate by a solvable birth-death process the Student distribution, which belongs to the third spectral category: this is due to the fact that to achieve this approximation, d(x) must be a polynomial of degree 2, with positive director coefficient and negative discriminant, which is in contradiction with the condition d(0) = 0; • We also get some birth-death processes that are not actual lattice approximation of any Pearson diffusion, such as in the Ehrenfest process case, whose state space is finite but d(x) is a polynomial of degree 1. • In particular, let us notice that for any solvable birth-death processes N (t) such that the polynomial b(x−1)−d(x) is of degree 1, the invariant measure m is in the Ord family (see [19]).
By using the relation m( which is the characterizing equation of the Ord family. If deg(d) = deg(b) = 1, then we can suppose that we still get equation (2.4). In particular we get the form which is the characterizing equation of the Katz family. Finally, last case is deg(d) = 1 and deg(b) = 0, that is to say in which the substitution that leads to equation (2.4) is given by Even in this case we are actually in the Katz family. Finally, let us observe that in such case the state space has to be infinite (since b(x) ≡ b 0 > 0) and then we are considering an immigration-death process (and the invariant measure is a Poisson measure). In particular we cover all the distributions of the Katz family (see [19]).
• It is also interesting to see that we cover the Poisson, Binomial, Negative Binomial and Hypergeometric invariant distribution cases, which are all in the cumulative Ord family (see [1]).

Inverse subordinators and non-local convolution derivatives
Now let us introduce our main object of study. Let us denote by BF the convex cone of Bernstein functions, that is to say Φ ∈ BF if and only if Φ ∈ C ∞ (R + ), Φ(λ) ≥ 0 and for any n ∈ N (−1) n d n Φ dλ n (λ) ≤ 0. In particular it is known that for Φ ∈ BF the following Lévy-Khintchine representation ( [42]) is given where a, b ≥ 0 and ν is a Lévy measure on R + such that The triple (a, b, ν) is called the Lévy triple of Φ. Also the vice versa can be shown, i.e. for any Lévy triple (a, b, ν) such that ν is a Lévy measure satisfying the integral condition (3.2) there exists a unique Bernstein function Φ such that Equation (3.1) holds. We will say that Φ is a driftless Bernstein function if a, b = 0 and ν(0, +∞) = +∞. Actually, the definition of driftless Bernstein function only requires b = 0, but the other two assumptions will be useful in our work. It is also known (see [42]) that for each Bernstein function Φ ∈ BF there exists a unique subordinator σ Φ = {σ Φ (y), y ≥ 0} (i. e. an increasing Lévy process) such that E[e −λσΦ(y) ] = e −yΦ(λ) . In particular we will say that σ Φ is driftless if Φ is driftless. For general notion on subordinators we refer to [11,Chapter 3] and [12]. In particular the hypothesis b = 0 ensure that σ Φ is a pure jump process, a = 0 implies that it is not killed and ν(0, +∞) = +∞ implies that σ Φ is strictly increasing (a. s.) and, for each y > 0, σ Φ (y) is an absolutely continuous random variable, hence it admits a density g Φ (x; y). Let us now fix our driftless Bernstein function Φ and its associated driftless subordinator σ Φ . Now we can define the inverse subordinator E Φ as, for any t > 0 Under our hypotheses, we have that E Φ (t) is absolutely continuous for any t > 0. Let us denote by f Φ (s; t) its density. Let us recall (see [32]) that, denoting by Now let us introduce the non-local convolution derivatives (of Caputo type) associated with Φ. Indeed, for Φ identified by the Lévy triple (0, 0, ν), let us define the Lévy tail ν(t) = ν(t, +∞). Now let us recall the definition of non-local convolution derivative, defined in [25] and [46].
Let us observe that one can define also the regularized version of the non-local convolution derivative as observing that it coincides with the previous definition on absolutely continuous functions. It can be shown, by Laplace transform arguments (see, for instance [26,3]) or by Green functions arguments (see [27]), that the (eigenvalue) Cauchy problem admits a unique solution for any λ > 0 and it is given by e Φ (t; λ) := E[e λEΦ(t) ] (hence, in particular, it is a completely monotone function in λ for fixed t). Let us recall that if Φ(λ) = λ α for α ∈ (0, 1), then ν(t) = t −α Γ(1−α) and d Φ dt Φ coincides with the fractional Caputo derivative of order α. In particular this means that in this case e Φ (t; λ) = E α (λt α ) where E α is the one-parameter Mittag-Leffler function defined, for t ∈ R as Let us recall (see [45]) that The proof is identical to the one of [5,Lemma 4.2]. We want to achieve a similar bound for any inverse subordinator. This is done by means of the following proposition. • We have already referred to the α-stable subordinator, i.e. the one we get when we choose Φ(λ) = λ α for α ∈ (0, 1). In such case, extensive informations on inverse α-stable subordinators are given in [33]. As we stated before, we have in particular e Φ (t; λ) = E α (λt α ), where E α is the one-parameter Mittag-Leffler function (see [13]). As a particular property, let us recall that if we denote by σ α the α-stable subordinator and g α the density of the random variable σ α (1), then the inverse α-stable subordinator E Φ (t) admits density • If we fix a constant θ > 0 and define Φ(λ) = (λ + θ) α − θ α we obtain the tempered α-stable subordinator with tempering parameter θ > 0. Denoting by σ α,θ (t) this subordinator, one can show that the density of the subordinator is given by where g α is the density of the α-stable subordinator σ α (t). An important property to recall is that the introduction of the tempering parameter implies the existence of all the moments of σ α,θ (t) (while this is not true for σ α ). Inverse tempered stable subordinators are studied for instance in [28]. Moreover, it can be shown that the Lévy tail ν is given by is the upper incomplete Gamma function; • For Φ(λ) = log(1 + λ α ) as α ∈ (0, 1) we obtain the geometric α-stable subordinator. From the form of the Bernstein function associated to the geometric α-stable subordinator, one obtains (see [44,Theorem 2.6]) that the density g G,α of the random variable σ G,α (1) (where σ G,α (t) is the geometric α-stable subordinator) satisfies the following asymptotics: Concerning the Lévy tail ν, it cannot be explicitly expressed, but it has been shown in [44, Theorem 2.5] that it satisfies the following asymptotic relation: as t → +∞.
• If in the previous example we consider α = 1 we obtain the Gamma subordinator. In this specific case, one can obtain explicitly the Lévy tail ν as (see, for instance, [44] and references therein) ν(t) = Γ(0; t).

Non-local forward and backward equations
Let us consider N (t) to be a solvable birth-death process with state space E and invariant measure m and let us denote by G and L its backward and forward operators respectively. Let us first focus on the backward equation.

Heuristic derivation of the strong solution. Let us consider a Bernstein function Φ and the Cauchy problem
x ∈ E.
Suppose that g ∈ ℓ 2 (m). Let us consider (Q n ) n∈E the family of orthonormal polynomials associated to N (t). Then we can decompose g as for some coefficients (g n ) n∈E . Let us suppose we want to find a solution u(t, x) by separation of variables. Thus let us suppose u(t, x) = T (t)ϕ(x). If we substitute this relation in the first equation of (4.1) we obtain the following coupled equations Concerning the first equation, let us observe that we need to set ϕ(x) = Q n (x) (up to a multiplicative constant) and λ = λ n for some n ∈ E. Let us also recall that λ n < 0. Concerning the second equation we get and then we have for some n ∈ E u(t, x) = Q n (x) e Φ (t; λ n ). Now, we can also have solutions that are linear combinations of Q n e Φ . Let us suppose that we can consider eventually infinite linear combinations. Then we expect a solution of the form for some coefficients (u n ) n∈E . Finally, let us observe that then, since the components (g n ) n∈E are uniquely determined, then u n = g n for any n ∈ E. Finally, we expect the solution to be of the form u(t, x) = n∈E g n Q n (x) e Φ (t; λ n ).
where Q n are the normalized orthogonal polynomials, absolutely converges for fixed t ≥ 0 and x, y ∈ E.
Proof. Let us first observe that if E is finite, then the summation (4.2) is actually finite. Thus let us consider the case in which E = N 0 . Let us denote by d n = P n ℓ 2 and let us recall that the dual polynomials P n exhibit orthogonality with respect to the measure m(x) = 1 d 2 x . However, since, if E = N 0 , N (t) is either an immigration-death process or a Meixner process, then the dual polynomials P n coincide with the polynomials P n themselves. We have, by using the self-duality relation P n (x) = P x (n), Now let us denote by root(x) the set of all the roots of the polynomial P x (n). These sets are finite with cardinality at most x. Thus we can define n 0 = ⌈max(root(x) ∪ root(y))⌉ + 1.
To show the absolute convergence of the series in p Φ (t, x; y), we only need to show the absolute convergence of +∞ n=n0 m(n) e Φ (t; λ n )P x (n)P y (n).
Let us observe (see [36, Table 2.3]) that the sign of the director coefficient of P x (both in the Charlier than in the Meixner case) depends on the parity of x. In particular we have that for any n ≥ n 0 it holds sign(P x (n)P y (n)) = (−1) x+y . So we get +∞ n=n0 | m(n) e Φ (t; λ n )P x (n)P y (n)| where we have to observe that e Φ (t, λ n ) ≤ 1 since λ n ≤ 0. As before, the series at the right-hand side converges if and only if +∞ n=0 m(n)P x (n)P y (n) converges. However, by the dual orthogonal relation and the self-duality relation of Charlier and Meixner polynomials, we achieve for any x, y ∈ N, concluding the proof.
Before exploiting the strong solution, let us show the total convergence of some auxiliary series of functions.
Lemma 4.2. Let N (t) be a solvable birth-death process with state space E = N 0 , generator G, invariant measure m and family of associated classical orthogonal polynomials (P n ) n≥0 . Let g ∈ ℓ 2 (m) such that g(x) = n≥0 g n Q n (x) for x ∈ E and some constants (g n ) n≥0 , where (Q n ) n≥0 are the normalized orthogonal polynomials. Then (2) For any fixed x ∈ E the sum n≥0 e Φ (t, λ n )g n Q n (x) totally converges for t ∈ [0, +∞); (3) For any fixed x ∈ E and T 1 > 0 the series n≥0 λ n e Φ (t, λ n )g n Q n (x) totally converges for t ∈ [T 1 , +∞).
Proof. Let us show property (1). Since g(x) = n≥0 g n Q n (x) and Q n is an orthonormal basis of ℓ 2 (m), it holds n≥0 g 2 n = g 2 ℓ 2 (m) . In particular it holds, by Cauchy-Schwartz inequality and self-duality relation for Charlier and Meixner polynomials To show property (2), let us just observe that e Φ (t, λ n ) ≤ 1 since λ n ≤ 0 and then where the series on the right-hand side is convergent and independent of t ≥ 0. Concerning property (3), by Equation (3.5) we obtain for some constant K(T 1 ) > 0, since e Φ (t, λ n ) is decreasing, concluding the proof.
Now we are ready to show that for the initial datum g ∈ ℓ 2 (m) our backward problem admits a solution.
Theorem 4.3. Let N (t) be a solvable birth-death process with state space E, generator G, invariant measure m and family of associated classical orthogonal polynomials (P n ) n∈E . Let g ∈ ℓ 2 (m) such that g(x) = n∈E g n Q n (x) for x ∈ E and some constants (g n ) n∈E , where (Q n ) n∈E are the normalized orthogonal polynomials. Then the Cauchy problem y ∈ E admits a unique strong solution of the form (4.4) u(t, y) = n∈E e Φ (t, λ n )g n Q n (y).
In particular sup t≥0 u(t, ·) ℓ 2 (m) ≤ g ℓ 2 (m) . Finally p Φ (t, x; y) is the fundamental solution of (4.3), in the sense that it is the strong solution of (4.3) for g(y) = δ x (y) and for any g ∈ ℓ 2 (m) it holds Before giving the proof of the Theorem, let us state formally what we mean as strong solution.
∂t Φ (t, y) exists for any t > 0 and y ∈ E; • The equations in (4.3) hold pointwise; • u(t; ·) ∈ C([0, +∞); ℓ 2 (m)); • ∂ Φ u ∂t Φ (t; ·) ∈ C((0, +∞); ℓ 2 (m)). Proof of Theorem 4.3. Let us first show that u(t, y) in the form of (4.4) is a strong solution for (4.3). First of all, let us recall that, by definition of e Φ and Q n , it holds . Now let us observe that if E is finite then n E ∈ N and we have that u(t, y) is a strong solution of (4.3) just by linearity of the involved operators. Let us now suppose that E is countably infinite. First of all, let us show that the series (4.4) converges in ℓ 2 (m). To do this, define for N ∈ N u N (t, y) = N n=0 e Φ (t, λ n )g n Q n (y). Now consider N < M in N and observe that, being λ n ≤ 0 and e Φ (t, λ n ) ≤ 1, it holds thus, by Cauchy's criterion, we know that the series converges in ℓ 2 (m). Now let us denote that is increasing and non-negative. Let us then observe that t 0 (u(τ, y) − u(0+, y))ν(t − τ )dτ = t 0 (u(τ, y) − u(0+, y))dI ν (t − τ ).
As next step, we want to differentiate under the series sign. However, we have to show uniform convergence for t in any compact set [T 1 , T 2 ] of the series of the derivatives to use [41,Theorem 7.17]. However, recalling (3.4), one has λ n e Φ (t, λ n )Q n (y)g n that totally converges by statement (3) of Lemma 4.2.
Hence we obtain, differentiating on both sides in (4.5), Now let us recall that G = (b(x) − d(x))∇ + + d(x)∆, hence we have to show that we can exchange ∇ + with the sign of series. This is due to the commutative property of totally convergent series. Indeed we have e Φ (t, λ n )∇ + Q n (y)g n = +∞ n=0 e Φ (t, λ n )∇ + Q n (y)g n where all the passages are justified by the fact that the two series +∞ n=0 e Φ (t, λ n )Q n (y+1)g n , +∞ n=0 e Φ (t, λ n )Q n (y)g n both totally converge by Lemma 4.2. The same holds for ∆. Thus we finally have for any t > 0. Concerning the initial datum, we have u(0, y) = n∈E g n Q n (y) = g(y). Now let us show strong continuity of u in [0, +∞) and of ∂ Φ ∂t Φ u in (0, +∞). These properties are obvious as E is finite, thus let us suppose E = N 0 . Concerning the continuity of u, let us show it in 0 + , since for any other t ∈ (0, +∞) the proof is analogous. We have Now fix ε > 0. Since (g n ) n≥0 belongs to ℓ 2 , there exists n(ε) ≥ 0 such that +∞ n=n(ε) g 2 n ≤ ε. By using the fact that (1 − e Φ (t, λ)) 2 ≤ 1 for any λ < 0, we get (1 − e Φ (t, λ n )) 2 g 2 n + ε.
Sending t → 0 + (using the fact that e Φ (t, λ n ) is continuous in t and that the first summation is finite) and then ε → 0 + we obtain strong continuity of u. Let us discuss the continuity of ∂ Φ ∂t Φ u in (0, +∞). To do this, let us consider t 0 ∈ [t 1 , t 2 ] with t 1 > 0. We have, for any t ∈ [t 1 , t 2 ], arguing as before and using (3.5), Thus, sending t → t 0 (observing that the first sum is finite) and then ε → 0 + we obtain the desired continuity. Uniqueness follows easily from the fact that (Q n ) n≥0 is an orthonormal system in ℓ 2 (m) hence the coefficients are unique. Now let us show the bound of the ℓ 2 (m) norm. We have Now let us consider a function g ∈ ℓ 2 (m). Then we have where we could exchange the order of the series since or E is finite, and then the sums are finite, or E is countably infinite and all the series involved are totally convergent in compact sets containing t.
Finally, let us observe that if for some z ∈ E g(x) = δ z (x), then concluding the proof.
4.3. The forward equation. Now let us apply the same strategy to study the Cauchy problem associated with L.
x ∈ E and some constants (f n ) n∈E , where (Q n ) n∈E are the normalized orthogonal polynomials. Then the Cauchy problem admits a unique strong solution of the form x; y) is the fundamental solution of (4.3), in the sense that it is the strong solution of (4. 3) for f (x) = δ y (x) and for any f /m ∈ ℓ 2 (m) it holds Proof. Let us first observe, by Lemma 2.1 . Thus, the exact same strategy of the proof of Theorem 4.3 leads to formula (4.7) and uniqueness follows as before. Concerning the estimate on the norm, it holds, since m is a probability measure and then m( Moreover, let us observe that Finally, observe that for some fixed z ∈ E, f (y) = δ z (y). Then obviously f /m ∈ ℓ 2 (m) and then we have concluding the proof.
Remark 4.5. One obtains also the following estimate on the norm: Next step is to identify some processes such that p Φ (t, x; y) is its transition probability function (in some sense) and then give some stochastic representation of the strong solutions of the Cauchy problems (4.4) and (4.7).

Non-local solvable birth-death processes
Let us now consider a solvable birth-death process N and the subordinator σ Φ associated to the Bernstein function Φ, with its inverse E Φ . Let us suppose N and E Φ are independent.
Definition 5.1. The non-local solvable birth-death process induced by N and Φ is defined as . Moreover, we define its transition probability function We used the same notation for the transition probability function and the fundamental solution of the Cauchy problems (4.4) and (4.7). Indeed, we can show the following result. If E is finite, we conclude the proof. So let us suppose E = N 0 . Recalling the notation in the proof of Lemma 4.1, let us define n 0 = ⌈max(root(x) ∪ root(y))⌉ + 1. and let us split the summation in two parts. We have Now let us observe that we can exchange the integral sign with the first summation by linearity of the integral, while in the second summation this can be done by Fubini's theorem, being the integrands of fixed sign (exactly the sign is determined by (−1) x+y since Q n (x)Q n (y) = m(n)P x (n)P y (n)). Thus we have Finally, let us recall that, by definition concluding the proof.
By using this result, we can achieve some stochastic representation of the solutions of the Cauchy problems (4.4) and (4.7). Indeed we have the following result.
Proposition 5.2. Let N Φ (t) be the non-local solvable birth-death process associated to N (t) and Φ. Suppose N (t) admits state space E. Then (1) For g ∈ ℓ 2 (m) the function u(t, y) is strong solution of (4.4); (2) For f /m ∈ ℓ 2 (m) such that f ≥ 0 and x∈E f (x) = 1, denoting by P f the probability measure obtained by P conditioning with the fact that N Φ (0) admits distribution f , then the function v(t, x) = P f (N Φ (t) = x) is strong solution of (4.7).
Proof. To show assertion 1, let us observe that obtaining that u(t, y) is the strong solution of (4.3) by Theorem 4.3. Concerning assertion 2, we have obtaining that v(t, x) is the strong solution of (4.6) by Theorem 4.4.
Finally, this proposition allows us to obtain the invariant distribution of N Φ (t) and show that it is also the limit distribution. Corollary 5.3. Let N Φ (t) be the non-local solvable birth-death process associated to N (t) and Φ. Suppose N (t) admits state space E. Then Proof. Let us first show property (1). To do this, let us observe that 1 ∈ ℓ 2 (m) since m is a probability measure. Thus, recall, by Proposition 5.2, that v(t, x) = P m (N Φ (t) = x) is a strong solution of (4.6). Now we need to determine m n such that n∈E m n Q n (x) = 1. Let us recall that Q 0 (x) = 1 while deg(Q n (x)) = n for any n = 1, . . . , n E , thus we have m 0 = 1 and m n = 0 for any n = 1, . . . , n E . Moreover, let us recall that λ 0 = 0, since 1 ∈ Ker(G). Hence we have, by Theorem and N Φ (t) admits m as distribution. Now let us suppose N Φ (0) admits f as distribution with f /m ∈ ℓ 2 (m). By Propo- Let us determine f 0 . We have, by definition of scalar product in ℓ 2 (m), Now let us consider the summation part. First of all, let us recall that m is a probability measure, hence m(x) ≤ 1. We have hence f ∈ ℓ 2 (m). By Lemma 4.2 we know that n∈E n≥1 e Φ (t, λ n )Q n (x)f n totally converges, hence we can take the limit as t → +∞ under the summation sign. Now let us observe that lim t→+∞ E Φ (t) = +∞ almost surely. On the other hand, we have e λnEΦ(t) ≤ 1, thus we can use monotone convergence theorem to achieve Finally, by dominated convergence theorem, we have concluding the proof.

Correlation structure of non-local solvable birth-death processes
Let us consider the potential measure U Φ (t) = E[E Φ (t)] of the subordinator σ Φ (t). As a consequence of Corollary 2.5, we can directly apply [7, Theorem 2] to obtain and expression of the covariance of N Φ (t) in terms of e Φ (t; λ 1 ). In particular, with an easy refinement of the proof, by using the distributional derivative of the potential function U Φ (t), we can get rid of some hypotheses of the aforementioned Theorem.
As we expected, since we are composing a stationary process N (t) (if N (0) admits m as distribution) with a non-stationary one E Φ (t), N Φ (t) is not secondorder stationary. To introduce the notion of memory for our process N Φ (t), we refer to the necessary conditions given in [10, Lemmas 2.1 and 2.2], since, for our processes, it is easier to work with the auto-covariance function. In particular we focus on the long memory with respect to the initial state. Definition 6.1. Given the function γ(n) = Cov m (N Φ (n), N Φ (0)) for n ∈ N: • N Φ (t) is said to exhibit long-range dependence if γ(n) ∼ ℓ(n)n −k where ℓ(n) is a slowly varying function and k ∈ (0, 1); • N Φ (t) is said to exhibit short-range dependence if +∞ n=1 |γ(n)| < +∞. In the specific case of the initial state, we have the following Proposition. Proposition 6.2. Let N (t) be a solvable birth-death process with state space E, invariant measure m and family of associated classical orthogonal polynomials (P n ) n∈E . Let us denote by ι : E → E the identity function and suppose that ι = a 0 + a 1 Q 1 . Then, lim t→+∞ Cov m (N Φ (t), N Φ (0)) = 0. Moreover, if e Φ (t; λ 1 ) ∼ Ct −α as t → +∞ for some α ∈ (0, 1), then N Φ (t) is long-range dependent.
Concerning the asymptotic behaviour of e Φ (t; λ 1 ), one can obtain some information by the behaviour of Φ. Indeed one has the following result. Proposition 6.4. If Φ is regularly varying at 0 + with order α ∈ (0, 1], then, for any fixed λ > 0, e Φ (t; −λ) is regularly varying at ∞ with order −α.
In particular we get the following Corollary.
Let us reconsider the examples given in Section 3 and study the asymptotic behaviour of the covariance.
• In the case Φ(λ) = λ α we can actually obtain an explicit formulation of the autocovariance function for t ≥ s > 0: Cov m (N Φ (t), N Φ (s)) = a 2 1 E α (λ 1 t α ) − where the proof of such formula is identical to the one of [29, Theorem 3.1].