Non-local Solvable Birth–Death Processes

Abstract

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.

1 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 [39] and queueing theory [26]. 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 [22, 25].

As birth–death processes are linked to difference–differential equations, fractionalization of such processes can be used to study the solutions of fractional difference–differential equations. Indeed, with this idea in mind, a fractional version of the Poisson process has been introduced in [10, 11] and fractional versions of some birth–death processes, for instance, in [40,41,42].

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, [31,32,33]). The same approach has been used to study the case of fractional immigration-death processes in [8]. Let us also stress out that this approach can be used to study a fractional $M/M/\infty$ queue or a fractional M/M/1 queue with acceleration of service (for some models of fractional queues, we refer to [6, 7, 17]).

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 [27] for the class of complete Bernstein functions and extended in [49] for any Bernstein function. A first step towards the theory of general time-changed Pearson diffusions has been achieved in [20].

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 Sect. 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 Sect. 3 we give some preliminaries on inverse subordinators and non-local 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 [27, 28] in the complete Bernstein case, in [36] in the special case and in [4] in the general case. Moreover, a series expansion in terms of convolutions of potential densities in the special Bernstein case is obtained in [4];

• In Sect. 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 Sect. 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 Sect. 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 [12, Lemmas 2.1 and 2.2].

2 Solvable Birth–Death Processes

Let us fix a filtered space $(\varOmega , \{\mathcal {F}_t\}_{t \in \mathbb {R}^+}, \mathbb {P})$ and consider a Birth–Death process $\{N(t),t \ge 0\}$ on it. Let us denote by $E \subseteq \mathbb {Z}$ its state space, that will be finite or at most countable. In particular we can always suppose that $E \subseteq \mathbb {N}_0$ and is a segment, i.e. for any $n_1,n_2 \in E$ and $n \in \mathbb {N}_0$ such that $n_1 \le n\le n_2$ it holds $n \in E$, with $\min E=0$. For a Birth–Death process N(t) we can consider the semigroup $T_t:f \in \mathfrak {b} \mapsto T_tf(y)=\mathbb {E}_y[f(N(t))]$ where $\mathfrak {b}$ is a suitable Banach sequence space containing $c_0$ (i.e. the space of infinitesimal sequences) whenever $E=\mathbb {N}_0$, while $\mathfrak {b}=\mathbb {R}^{|E|}$ if E is finite. Since we will consider Birth–Death processes with a stationary measure $\mathbf {m}$, we will always set $\mathfrak {b}=\ell ^2(\mathbf {m})$ (which is equivalent to $\mathbb {R}^{|E|}$ if E is finite). Let us recall that the generator $\mathcal {G}$ of a Birth–Death process can be always expressed as

$$\begin{aligned} \mathcal {G}f(x)=(b(x)-d(x))\nabla ^+ f(x)+d(x)\varDelta f(x), \qquad x \in E, \end{aligned}$$

where d(x) are the death rates, b(x) are the birth rates (recalling that d, b must be non-negative in E), $\nabla ^{\pm }$ are the first-order forward and backward finite differences defined as

$$\begin{aligned} \nabla ^+ f(x)=f(x+1)-f(x)&\nabla ^- f(x)=f(x)-f(x-1) \end{aligned}$$

and $\varDelta$ is the second-order central finite difference

$$\begin{aligned} \varDelta f(x)=f(x+1)-2f(x)+f(x-1)=\nabla ^-\nabla ^+f(x). \end{aligned}$$

Let us stress out that $\mathcal {G}$ is actually a tridiagonal matrix, possibly infinite if $E=\mathbb {N}_0$. Moreover, as $T_t$ is defined on $\ell ^2(\mathbf {m})$, also the generator $\mathcal {G}$ admits $\ell ^2(\mathbf {m})$ as domain.

Here we want to introduce some birth–death version of Pearson diffusions (see, for instance, [32]). To do this, we refer to the theory of birth–death polynomials, whose main papers are [22, 25].

Definition 1

We say the process N(t) is solvable if

• N(t) is irreducible and recurrent;

• the spectrum of $\mathcal {G}$ is purely discrete with non-positive eigenvalues $(\lambda _n)_{n \in E}$ such that $\lambda _0=0$ and $\lambda _n<0$ for any $n \ge 1$;

• its eigenfunctions $(P_n)_{n \in E}$ are classical orthogonal polynomials of discrete variable with orthogonality measure $\mathbf {m}$ which is the invariant and stationary measure of N(t);

• the function $m(x)=\mathbf {m}(\{x\})$ solves the following discrete Pearson equation:

  $$\begin{aligned} \nabla ^+(d(\cdot )m(\cdot ))(x)=(b(x)-d(x))m(x) \end{aligned}$$

  (1)

• $d(\cdot )$ is a polynomial of degree at most 2 and $b(\cdot )-d(\cdot )$ 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

$$\begin{aligned} \lambda _n=n \nabla ^+(b(\cdot )-d(\cdot ))(x)+\frac{1}{2}n(n-1)\varDelta d(x). \end{aligned}$$

Concerning classical orthogonal polynomials of discrete variable, we mainly refer to [38, 46]. Their orthogonality relation is expressed as

$$\begin{aligned} \sum _{x \in E}P_n(x)P_m(x)m(x)=\mathfrak {d}_n^2 \delta _{n,m}, \qquad \ n,m \in E \end{aligned}$$

where $\delta _{n,m}$ is Kronecker delta symbol. In particular one obtains that $\left\| P_n \right\| _{\ell ^2(\mathbf {m})}=\mathfrak {d}_n$ and then we can introduce the normalized polynomials as $Q_n(x)=\frac{P_n(x)}{\mathfrak {d}_n}$, such that

$$\begin{aligned} \sum _{x \in E}Q_n(x)Q_m(x)m(x)=\delta _{n,m}, \qquad \ n,m \in E. \end{aligned}$$

On the other hand, the function $\widetilde{m}(n)=\frac{1}{\mathfrak {d}^2_n}$ defines a measure on E. Thus, by proceeding with a Gram-Schmidt orthogonalization procedure on the monomials $(1,x,x^2,\cdots )$, we can define a family of orthogonal polynomials $\widetilde{P}_n(x)$ that satisfies the following orthogonality condition

$$\begin{aligned} \sum _{x \in E}\widetilde{P}_n(x)\widetilde{P}_m(x)\widetilde{m}(x)=\frac{1}{m(n)} \delta _{n,m}, \qquad \ n,m \in E. \end{aligned}$$

The family of polynomials $(\widetilde{P}_n)_{n \in E}$ is called the dual family of $(P_n)_{n \in E}$, see [38]. 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^{-\rho }\frac{\rho ^x}{x!}$ for $x \in \mathbb {N}$, where $\rho =\frac{b}{d}$. The orthogonal polynomials are Charlier polynomials of parameter $\rho$, that we denote by $P_n(x)=C_n(x;\rho )$ and they satisfy the self-duality relation

  $$\begin{aligned} P_n(x)=P_x(n) \end{aligned}$$

  (2)

  and the eigenvalues are given by $\lambda _n=-bn$. In particular the polynomials $P_n$ coincide with the family of dual orthogonal polynomials $\widetilde{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+\beta )b$, with $b,d,\beta >0$ and $b<d$. In such case the birth death process admits state space $E=\mathbb {N}$ and the generator is given by

  $$\begin{aligned} \mathcal {G}=(\beta b+(b-d)x)\nabla ^++dx\varDelta . \end{aligned}$$

  Defining $\rho =\frac{b}{d}$, we have that the orthogonal polynomials are Meixner polynomials of parameters $\rho$ and $\beta$, that we denote by $P_n(x)=M_n(x;\rho ,\beta )$ and they are orthogonal with respect to the invariant measure

  $$\begin{aligned} m(x)=\frac{(\beta )_x \rho ^x}{x! (1-\rho )^\beta }, \end{aligned}$$

  where $(\beta )_x=\frac{\varGamma (\beta +x)}{\varGamma (\beta )}$, which is a Pascal (or negative binomial) distribution of parameters $\beta$ and $\rho$. Also Meixner polynomials satisfy the self-duality relation (2) and coincide with their dual polynomials. Finally, let us observe that the eigenvalues are given by $\lambda _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 [23].

• Another case is given by a linear death rate $d(x)=dx$ and a linear decreasing birth rate $b(x)=(N-x)b$ with $b,d>0$ and $N \in \mathbb {N}$. In such case the birth–death process admits finite state space $E=\{0,\dots ,N\}$ and the generator is given by

  $$\begin{aligned} \mathcal {G}=(N b-(b+d)x)\nabla ^++dx\varDelta . \end{aligned}$$

  Defining $p=\frac{b}{b+d}$ and $q=1-p$, we achieve the invariant distribution given by

  $$\begin{aligned} m(x)=\left( {\begin{array}{c}N\\ x\end{array}}\right) p^xq^{N-x} \end{aligned}$$

  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). The eigenvalues of the generator are given by $\lambda _n=-n(b+d)$. Let us recall that this is actually a time-continuous version of the Ehrenfest urn model (see, for instance, [24]).

• Another interesting case is given by a quadratic one. Indeed let us consider for some $d>0$ and $\alpha ,\beta ,N \in \mathbb {N}$

  $$\begin{aligned} d(x)=dx(N+\beta +1-x)&b(x)=d[N(\alpha +1)+x(N-1-\alpha )-x^2]. \end{aligned}$$

  Let us observe that $b(N)=0$, thus the state space of the process is given by $E=\{0,\dots ,N\}$. Moreover, its generator is given by

  $$\begin{aligned} \mathcal {G}=d(N(\alpha +1)-(\beta +\alpha +2)x)\nabla ^++dx(N+\beta +1-x)\varDelta \end{aligned}$$

  with eigenvalues

  $$\begin{aligned} \lambda _n=-dn[n+1+\alpha +\beta ]. \end{aligned}$$

  The invariant measure of this birth–death process is an hypergeometric distribution on E given by

  $$\begin{aligned} m(x)=\left( {\begin{array}{c}\alpha +x\\ x\end{array}}\right) \left( {\begin{array}{c}\beta +N-x\\ N-x\end{array}}\right) \end{aligned}$$

  and the orthogonal polynomials are the Hahn polynomials, that we denote by $P_n(x)=H_n(x;\alpha ,\beta ,N)$. In this case we do not have self-duality relation, but the family of dual Hahn polynomials, that we denote by $\widetilde{P}_n(x)=R_n(x;\alpha ,\beta ,N)$, is linked to the Hahn polynomials by the relation

  $$\begin{aligned} H_n(x;\alpha ,\beta ,N)=R_x(n(n+\alpha +\beta +1);\alpha ,\beta ,N). \end{aligned}$$

  This particular birth–death process is a lattice approximation of the Jacobi process. For such process, we refer directly to [46].

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 [32, 33]).

Let us also recall that if we consider the matrix $p(t)=(p(t,x;y))_{x,y \in E}$ where $p(t,x;y)=\mathbb {P}(N(t+s)=x|N(s)=y)$ is the transition probability function, then p(t) solves both the so-called backward and forward equations (see [46])

$$\begin{aligned} p'(t)=\mathcal {G}\cdot p(t)&p'(t)=p(t) \cdot \mathcal {G}. \end{aligned}$$

In particular we can consider the adjoint operator $\mathcal {L}$ of $\mathcal {G}$ (still defined on $\ell ^2(\mathbf {m})$) such that $p'(t)=\mathcal {L}p(t)$. For this reason, we refer to $\mathcal {L}$ as the forward operator. Let us observe that for a birth–death process N(t) the forward operator $\mathcal {L}$ is defined as

$$\begin{aligned} \mathcal {L}f(x)=-\nabla ^-((b(\cdot )-d(\cdot ))f(\cdot ))(x)+\varDelta (d(\cdot )f(\cdot ))(x). \end{aligned}$$

Moreover, let us observe that $\mathcal {L}$ can be considered as a (possibly infinite) tridiagonal matrix such that $p'(t)=\mathcal {L}\cdot p(t)^{\mathbf {t}}$ where with $\mathbf {t}$ we denote the transpose. Now let us show what the orthogonal polynomials $P_n(x)$ represent for the forward operator.

Lemma 1

Let N(t) be a solvable birth–death process with forward operator $\mathcal {L}$, invariant measure $\mathbf {m}$ and associated family of orthonormal polynomials $Q_n(x)$. Then we have for any $n \in E$ and any $x \in E$,

$$\begin{aligned} \mathcal {L}(mQ_n)(x)=m(x)\lambda _nQ_n(x). \end{aligned}$$

Proof

We have, recalling that $\varDelta =\nabla ^-\nabla ^+$ and that $\nabla ^-$ is a linear operator,

$$\begin{aligned} \mathcal {L}(mQ_n)(x)&=-\nabla ^-((b-d)mQ_n)(x)+\varDelta (dmQ_n)(x)\\&=\nabla ^-[-(b-d)mQ_n+\nabla ^+(dmQ_n)](x). \end{aligned}$$

Now, by discrete Leibniz rule on $\nabla ^+$, denoting $\widetilde{d}(x)=d(x+1)$ and $\widetilde{m}=m(x+1)$, we have

$$\begin{aligned} \mathcal {L}(mQ_n)(x)&=\nabla ^-[-(b-d)mQ_n+Q_n\nabla ^+(dm)+\widetilde{d}\widetilde{m}\nabla ^+Q_n](x)\\&=\nabla ^-[Q_n(-(b-d)mQ_n+\nabla ^+(dm))+\widetilde{d}\widetilde{m}\nabla ^+Q_n](x)\\&=\nabla ^-(\widetilde{d}\widetilde{m}\nabla ^+Q_n)(x). \end{aligned}$$

Now let us use again Leibniz rule on $\nabla ^-$ to achieve

$$\begin{aligned} \mathcal {L}(mQ_n)(x)=d(x)m(x)\varDelta Q_n(x)+\nabla ^+Q_n(x)\nabla ^-(\widetilde{d}\widetilde{m})(x). \end{aligned}$$

Now let us work with $\nabla ^-(\widetilde{d}\widetilde{m})(x)$. We have

$$\begin{aligned} \nabla ^-(\widetilde{d}\widetilde{m})(x)=d(x+1)m(x+1)-d(x)m(x). \end{aligned}$$

However, by the discrete Pearson equation (1) we obtain

$$\begin{aligned} d(x+1)m(x+1)-d(x)m(x)=(b(x)-d(x))m(x) \end{aligned}$$

and then we have

$$\begin{aligned} \nabla ^-(\widetilde{d}\widetilde{m})(x)=(b(x)-d(x))m(x). \end{aligned}$$

Hence, we achieve

$$\begin{aligned} \mathcal {L}(mQ_n)(x)&=m(x)[(b(x)-d(x))\nabla ^+Q_n(x)+d(x)\varDelta Q_n(x)]\\&=m(x)\mathcal {G}Q_n(x)=m(x)\lambda _n Q_n(x), \end{aligned}$$

concluding the proof. $\square$

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 (which is a direct consequence of the results in [22] once we notice that in our case the spectral measure coincides with the stationary one):

Theorem 1

Let N(t) be a solvable birth–death process with state space E, generator $\mathcal {G}$, forward operator $\mathcal {L}$, invariant measure $\mathbf {m}$ and family of associated orthonormal polynomials $(Q_n)_{n \in E}$. Then the following assertions hold true:

• The transition probability function $p(t,x;y)=\mathbb {P}(N(t+s)=x|N(s)=y)$ for $x,y \in E$ and $t,s \ge 0$ admit the following spectral representation:

  $$\begin{aligned} p(t,x;y)=m(x)\sum _{n\in E}e^{\lambda _n t}Q_n(x)Q_n(y) \end{aligned}$$

  for any $x,y \in E$ and $t \ge 0$.

• If $g \in \ell ^2(\mathbf {m})$ with $g(x)=\sum _{n \in E}g_nQ_n(x)$ then the strong solution of the Cauchy problem

  $$\begin{aligned} {\left\{ \begin{array}{ll} \frac{\mathrm{d} u}{\mathrm{d} t}(t,y)=\mathcal {G}u(t,y) &{} t \ge 0, y \in E \\ u(0,y)=g(y) &{} y \in E \end{array}\right. } \end{aligned}$$

  is given by

  $$\begin{aligned} u(t,y)=\sum _{n \in E}g_ne^{\lambda _n t}Q_n(y)=\sum _{x \in E}p(t,x;y)g(x). \end{aligned}$$

  In particular p(t, x; y) is the fundamental solution of $\frac{\mathrm{d} u}{\mathrm{d} t}(t,y)=\mathcal {G}u(t,y)$ and u admits the stochastic interpretation:

  $$\begin{aligned} u(t,y)=\mathbb {E}_y[g(N(t))] \end{aligned}$$

  where $\mathbb {E}_y[\cdot ]=\mathbb {E}[\cdot | N(0)=y]$.

• If $f/m \in \ell ^2(\mathbf {m})$ with $\frac{f(x)}{m(x)}=\sum _{n \in E}f_nQ_n(x)$ the strong solution of the Cauchy problem

  $$\begin{aligned} {\left\{ \begin{array}{ll} \frac{\mathrm{d} v}{\mathrm{d} t}(t,x)=\mathcal {L}u(t,x) &{} t \ge 0, x \in E \\ v(0,x)=f(x) &{} x \in E \end{array}\right. } \end{aligned}$$

  is given by

  $$\begin{aligned} v(t,x)=m(x)\sum _{n \in E}f_ne^{\lambda _n t}Q_n(x)=\sum _{y\in E}p(t,x;y)f(y). \end{aligned}$$

  In particular p(t, x; y) is the fundamental solution of $\frac{\mathrm{d} v}{\mathrm{d} t}(t,x)=\mathcal {L}v(t,x)$ and, if $f \ge 0$ with $\left\| f \right\| _{\ell ^1}=1$, v admits the stochastic interpretation:

  $$\begin{aligned} v(t,x)=\mathbb {P}_f(N(t)=x) \end{aligned}$$

  where $\mathbb {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 birth–death 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=\mathbb {N}_0$, let us first show the following Proposition.

Proposition 1

Let N(t) be a solvable birth–death process with invariant measure $\mathbf {m}$ and state space $E=\mathbb {N}_0$. Then there exists a constant $\rho <1$ and a state $x_0 \in E$ such that for any $x \ge x_0$ it holds

$$\begin{aligned} m(x)\le \rho ^{x-x_0}m(x_0) \end{aligned}$$

(3)

Proof

First of all, let us observe that since b(x) and d(x) are polynomials, then $\lim _{x \rightarrow +\infty }\frac{b(x)}{d(x+1)}$ always exists. Moreover, we can rewrite the discrete Pearson equation as

$$\begin{aligned} m(x+1)=\frac{b(x)}{d(x+1)}m(x). \end{aligned}$$

Thus, we have that m is well defined if and only if

$$\begin{aligned} \sum _{x=1}^{+\infty }\prod _{k=0}^{x}\frac{b(k)}{d(k+1)}<+\infty . \end{aligned}$$

It is easy to see that such condition implies $\lim _{x \rightarrow +\infty }\frac{b(x)}{d(x+1)}\le 1$. Let us suppose $\lim _{x \rightarrow +\infty }\frac{b(x)}{d(x+1)}=1$. 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 $\lambda _n=0$ for any $n \ge 1$, that is absurd. Thus, we have that b(x) and d(x) are polynomials of degree 2. However, since $\lambda _n<0$ for any $n \ge 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

$$\begin{aligned} \lim _{x \rightarrow +\infty }\frac{b(x)}{d(x+1)}=l<1. \end{aligned}$$

Now let us consider $\rho \in (l,1)$. Then there exists a state $x_0 \in E$ such that $\frac{b(x)}{d(x+1)}<\rho$ as $x \ge x_0$. Thus, we have

$$\begin{aligned} m(x+1)<\rho m(x) \end{aligned}$$

for any $x \ge x_0$. Finally, the assertion follows from the previous inequality by induction. $\square$

As a direct consequence of the previous proposition we obtain

Corollary 1

Let N(t) be a solvable birth–death process with invariant distribution $\mathbf {m}$ such that N(0) admits $\mathbf {m}$ as distribution. Then N(t) admits moments of any order.

Now we can focus on the autocovariance function of the process N(t).

Corollary 2

Let N(t) be a solvable birth–death process with invariant measure $\mathbf {m}$. Then there exists a constant $a_1 \in \mathbb {R}$ such that, for any $t,s \ge 0$

$$\begin{aligned} Cov_m(N(t),N(s))=a_1^2e^{\lambda _1|t-s|}. \end{aligned}$$

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 \ge s$,

$$\begin{aligned} Cov_m(N(t),N(s))=Cov_m(N(t-s)N(0)) \end{aligned}$$

Thus, let us consider $t \ge 0$ and let us evaluate $Cov_m(N(t),N(0))$. To do this, let us rewrite

$$\begin{aligned} Cov_m(N(t),N(0))= & {} \mathbb {E}_m[N(t)N(0)]-\mathbb {E}_m[N(t)]\mathbb {E}_m[N(0)]\\= & {} \mathbb {E}_m[N(t)N(0)]-\mathbb {E}_m[N(0)]^2. \end{aligned}$$

Now let us first evaluate $\mathbb {E}_x[N(t)]$. Since $\mathbf {m}$ admits second moment then $\iota (x)=x$ is in $\ell ^2(\mathbf {m})$. Moreover, since $\deg (\iota (x))=1$, it can be written as a linear combination of $Q_0=1$ and $Q_1$. Let then

$$\begin{aligned} \iota (x)=a_0+a_1Q_1(x). \end{aligned}$$

By the previous theorem we have that

$$\begin{aligned} \mathbb {E}_x[N(t)]=a_0+a_1e^{\lambda _1 t}Q_1(x) \end{aligned}$$

(recalling that $\lambda _0=0$ for any solvable birth–death process). Starting from this observation, we have that

$$\begin{aligned} \mathbb {E}_m[N(t)N(0)]=\sum _{x \in E}xm(x)\mathbb {E}_x[N(t)]=a_0\sum _{x \in E}xm(x)+a_1e^{\lambda _1 t}\sum _{x \in E}xm(x)Q_1(x). \end{aligned}$$

As we stated before, we can write $x=a_0+a_1Q_1(x)$, thus we have

$$\begin{aligned} \mathbb {E}_m[N(t)N(0)]&=a^2_0\sum _{x \in E}m(x)+a_0a_1\sum _{x \in E}Q_0Q_1(x)m(x)\\&\quad +a_0a_1e^{\lambda _1 t}\sum _{x \in E}Q_0m(x)Q_1(x)+a^2_1e^{\lambda _1 t}\sum _{x \in E}m(x)Q^2_1(x). \end{aligned}$$

By using the orthonormality relation we have

$$\begin{aligned} \mathbb {E}_m[N(t)N(0)]=a^2_0+a^2_1e^{\lambda _1 t}. \end{aligned}$$

Now let us evaluate $\mathbb {E}_m[N(0)]$. We have

$$\begin{aligned} \mathbb {E}_m[N(0)]&=\sum _{x \in E}xm(x)=a_0+a_1\sum _{x \in E}Q_1(x)m(x)\\&=a_0+a_1\sum _{x \in E}Q_0Q_1(x)m(x)=a_0. \end{aligned}$$

We finally achieve

$$\begin{aligned} Cov_m(N(t),N(0))=a_1^2e^{\lambda _1 t} \end{aligned}$$

concluding the proof. $\square$

2.1 Classification of Solvable Birth–Death Processes

We can actually improve the result in Proposition 1 by obtaining a complete classification of solvable birth–death processes. Indeed we have the following Proposition.

Proposition 2

Let N(t) be a solvable birth–death process with state space E. Then one of the following statements holds true:

• E is finite;

• N(t) is an immigration-death process;

• N(t) is a Meixner process.

In particular, if $(P_n)_{n \in E}$ is the family of orthogonal polynomials associated to N(t), then either E is finite or $(P_n)_{n \in E}$ coincide with its dual family and $P_n(x)=P_x(n)$ for any $x,n \in E$.

Proof

First of all, suppose $E=\mathbb {N}_0$. Then we have $\lim _{x \rightarrow +\infty }\frac{b(x)}{d(x+1)}<1$. In particular this implies that $\deg b \le \deg d$. If $\deg d=0$, then also $\deg b=0$ and in such case $\lambda _n=0$ for any $n \in \mathbb {N}$, which is absurd. Thus, $\deg d=1,2$. However, we have already seen that if $\deg d=2$, then, since $\lambda _n<0$ for any $n \ge 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=\mathbb {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 \in 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=\mathbb {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;

• $\deg b=1$, thus $b(x)=b(x+\beta )$ for some $b,\beta >0$ and then we get a Meixner process (since also $b<d$ by the condition $\lim _{x \rightarrow +\infty }\frac{b(x+\beta )}{d(x+1)}=\frac{b}{d}<1$).

This concludes the proof. $\square$

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 [19]), however, their spectral behaviour can be distinguished in three different classes (see [33]). 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 $\mathbb {N}_0$, we can still cover these cases, in which d(x) admits discriminant $\varDelta \ge 0$ and $d(0)=0$, by asking that b(x) admits the first root $x_0 \in \mathbb {N}$ and $E=\{0,\dots ,x_0\}$. If $\varDelta >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 $\mathbb {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 $\mathbf {m}$ is in the Ord family (see [21]). Indeed if $\deg (b)=\deg (d)=2$, we can suppose that (since $d(0)=0$)

  $$\begin{aligned} d(x)=\widetilde{a}x^2+\widetilde{d}_1x, \qquad b(x)=\widetilde{a}x^2+\widetilde{b}_1x+\widetilde{b}_2. \end{aligned}$$

  By using the relation $m(x-1)=\frac{d(x)}{b(x-1)}m(x)$ and setting

  $$\begin{aligned}&k=2\widetilde{a}+\widetilde{d}_1-\widetilde{b}_1,&a=\frac{\widetilde{a}-\widetilde{b}_1+\widetilde{b}_2}{k},\\&b_0=0,&b_1=\frac{\widetilde{d}_1+\widetilde{a}}{k},&b_2=\frac{\widetilde{a}}{k}, \end{aligned}$$

  where $k \not = 0$ since it is the director coefficient of the first degree polynomial $b(x-1)-d(x)$, we get the equation

  $$\begin{aligned} \frac{\nabla ^+m(x-1)}{m(x)}=\frac{a-x}{(a+b_0)+(b_1-1)x+b_2x(x-1)} \end{aligned}$$

  (4)

  which is the characterizing equation of the Ord family.

  If $\deg (d)=\deg (b)=1$, then we can suppose that

  $$\begin{aligned} d(x)=\widetilde{d}_1x, \qquad b(x)=\widetilde{b}_1x+\widetilde{b}_0. \end{aligned}$$

  To have $\deg (b(x-1)-d(x))=1$, we need $\widetilde{d}_1 \not = \widetilde{b}_1$. Thus, setting

  $$\begin{aligned} k=\widetilde{d}_1-\widetilde{b}_1 \not =0,&a=\frac{\widetilde{b}_0-\widetilde{b}_1}{k}, \end{aligned}$$
  $$\begin{aligned} b_0=0,&b_1=\frac{\widetilde{b}_0}{k},&b_2=0, \end{aligned}$$

  we still get Eq. (4). In particular we get the form

  $$\begin{aligned} \frac{\nabla ^+m(x-1)}{m(x)}=\frac{a-x}{a+(b_1-1)x} \end{aligned}$$

  (5)

  which is the characterizing equation of the Katz family. Finally, last case is $\deg (d)=1$ and $\deg (b)=0$, that is to say

  $$\begin{aligned} d(x)=\widetilde{d}_1x\qquad b(x)\equiv \widetilde{b}_0 \end{aligned}$$

  in which the substitution that leads to Eq. (4) is given by

  $$\begin{aligned}&k=\widetilde{d}_1,&a=\frac{\widetilde{b}_0}{k},\\&b_0=0,&b_1=1,&b_2=0. \end{aligned}$$

  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)\equiv \widetilde{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 [21]).

• 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]).

3 Inverse Subordinators and Non-local Convolution Derivatives

Now let us introduce our main object of study. Let us denote by $\mathcal {BF}$ the convex cone of Bernstein functions, that is to say $\varPhi \in \mathcal {BF}$ if and only if $\varPhi \in C^\infty (\mathbb {R}^+)$, $\varPhi (\lambda ) \ge 0$ and for any $n \in \mathbb {N}$

$$\begin{aligned} (-1)^n\frac{\mathrm{d}^{n} \varPhi }{\mathrm{d} \lambda ^{n}}(\lambda )\le 0. \end{aligned}$$

In particular it is known that for $\varPhi \in \mathcal {BF}$ the following Lévy-Khintchine representation ([45]) is given

$$\begin{aligned} \varPhi (\lambda )=a+b\lambda +\int _0^{+\infty }(1-e^{-\lambda t})\nu (\mathrm{d}t) \end{aligned}$$

(6)

where $a,b \ge 0$ and $\nu$ is a Lévy measure on $\mathbb {R}^+$ such that

$$\begin{aligned} \int _0^{+\infty }(1 \wedge t)\nu (\mathrm{d}t)<+\infty . \end{aligned}$$

(7)

The triple $(a,b,\nu )$ is called the Lévy triple of $\varPhi$. Also the vice versa can be shown, i.e. for any Lévy triple $(a,b,\nu )$ such that $\nu$ is a Lévy measure satisfying the integral condition (7) there exists a unique Bernstein function $\varPhi$ such that Eq. (6) holds. In the rest of the paper, we suppose that $a=b=0$ and $\nu (0,+\infty )=+\infty$. It is also known (see [45]) that for each Bernstein function $\varPhi \in \mathcal {BF}$ there exists a unique subordinator $\sigma _\varPhi =\{\sigma _\varPhi (y), y \ge 0\}$ (i. e. an increasing Lévy process) such that

$$\begin{aligned} \mathbb {E}[e^{-\lambda \sigma _\varPhi (y)}]=e^{-y\varPhi (\lambda )}. \end{aligned}$$

For general notion on subordinators we refer to [13, Chapter 3] and [14]. In particular the hypothesis $b=0$ ensure that $\sigma _\varPhi$ is a pure jump process, $a=0$ implies that it is not killed and $\nu (0,+\infty )=+\infty$ implies that $\sigma _\varPhi$ is strictly increasing (a. s.).

Let us now fix our Bernstein function $\varPhi$ and its associated subordinator $\sigma _\varPhi$. Now we can define the inverse subordinator $E_\varPhi$ as, for any $t>0$

$$\begin{aligned} E_\varPhi (t):=\inf \{y \ge 0: \ \sigma _\varPhi (y)>t\}. \end{aligned}$$

Under our hypotheses, we have that $E_\varPhi (t)$ is absolutely continuous for any $t>0$ and its sample paths are almost surely continuous. Let us denote by $f_\varPhi (s;t)$ its density (see [34]). Let us recall (see [34]) that, denoting by $\overline{f}_\varPhi (s;\lambda )$ the Laplace transform of $f_\varPhi (s;t)$ with respect to t,

$$\begin{aligned} \overline{f}_\varPhi (s;\lambda )=\frac{\varPhi (\lambda )}{\lambda }e^{-s\varPhi (\lambda )}. \end{aligned}$$

Now let us introduce the non-local convolution derivatives (of Caputo type) associated with $\varPhi$. Indeed, for $\varPhi$ identified by the Lévy triple $(0,0,\nu )$, let us define the Lévy tail $\overline{\nu }(t)=\nu (t,+\infty )$. Now let us recall the definition of non-local convolution derivative, defined in [27, 49].

Definition 2

Let $f:\mathbb {R}^+ \rightarrow \mathbb {R}$ be an absolutely continuous function. Then we define the non-local convolution derivative induced by $\varPhi$ of f as

$$\begin{aligned} \frac{\mathrm{d}^{\varPhi } }{\mathrm{d} t^{\varPhi }}f(t)=\int _0^t f'(\tau )\overline{\nu }(t-\tau )\mathrm{d}\tau . \end{aligned}$$

(8)

Let us observe that one can define also the regularized version of the non-local convolution derivative as

$$\begin{aligned} \frac{\mathrm{d}^{\varPhi } }{\mathrm{d} t^{\varPhi }}f(t)=\frac{\mathrm{d} }{\mathrm{d} t}\int _0^t (f(\tau )-f(0+))\overline{\nu }(t-\tau )\mathrm{d}\tau \end{aligned}$$

(9)

observing that it coincides with the previous definition on absolutely continuous functions. Here we will always consider the regularized form (9).

It can be shown, by Laplace transform arguments (see, for instance [4, 28]) or by Green functions arguments (see [29]), that the (eigenvalue) Cauchy problem

$$\begin{aligned} {\left\{ \begin{array}{ll} \frac{\mathrm{d}^{\varPhi } }{\mathrm{d} t^{\varPhi }} \mathfrak {e}_\varPhi (t;\lambda )=\lambda \mathfrak {e}_\varPhi (t;\lambda ) &{} t>0\\ \mathfrak {e}_\varPhi (0;\lambda )=1 \end{array}\right. } \end{aligned}$$

admits a unique solution for any $\lambda >0$ and it is given by $\mathfrak {e}_\varPhi (t;\lambda ):=\mathbb {E}[e^{\lambda E_\varPhi (t)}]$ (hence, in particular, it is a completely monotone function in $\lambda$ for fixed t). Let us recall that if $\varPhi (\lambda )=\lambda ^\alpha$ for $\alpha \in (0,1)$, then $\overline{\nu }(t)=\frac{t^{-\alpha }}{\varGamma (1-\alpha )}$ and $\frac{\mathrm{d}^{\varPhi } }{\mathrm{d} t^{\varPhi }}$ coincides with the fractional Caputo derivative of order $\alpha$. In particular this means that in this case $\mathfrak {e}_{\varPhi }(t;\lambda )=E_\alpha (\lambda t^{\alpha })$ where $E_\alpha$ is the one-parameter Mittag-Leffler function defined, for $t \in \mathbb {R}$ as

$$\begin{aligned} E_\alpha (t)=\sum _{k=0}^{+\infty }\frac{t^k}{\varGamma (\alpha k+1)}. \end{aligned}$$

Let us recall (see [48]) that

$$\begin{aligned} E_\alpha (-\lambda t^\alpha )\le \frac{1}{1+\frac{t^\alpha }{\varGamma (1+\alpha )}\lambda } \end{aligned}$$

hence it is not difficult to show the following Proposition.

Proposition 3

For any $\lambda >0$ it holds

$$\begin{aligned} \lambda E_\alpha (-\lambda t^\alpha )\le \frac{\varGamma (1+\alpha )}{t^\alpha }. \end{aligned}$$

The proof is identical to the one of [8, Lemma 4.2]. We want to achieve a similar bound for any inverse subordinator. This is done by means of the following proposition.

Proposition 4

Fix $t>0$. Then there exists a constant K(t) such that

$$\begin{aligned} \lambda \mathfrak {e}_\varPhi (t;-\lambda )\le K(t), \ \forall \lambda \in [0,+\infty ). \end{aligned}$$

(10)

Proof

Let us first recall that $\mathfrak {e}_\varPhi (t;-\lambda )=\mathbb {E}[e^{-\lambda E_\varPhi (t)}]$, thus it is the Laplace transform of $f_\varPhi (s;t)$ with respect to s. In particular it is completely monotone in $\lambda$ and $\mathfrak {e}_\varPhi (t;0)=1$. Now let us recall that $f_\varPhi (0+;t)=\overline{\nu }(t)$ (see, for instance, [49, Theorem 4.1]). On the other hand, by the initial-value theorem (see, for instance, [18, Section 17.8]), we have

$$\begin{aligned} \lim _{\lambda \rightarrow +\infty }\lambda \mathfrak {e}_\varPhi (t;-\lambda )=f_\varPhi (0+;t)=\overline{\nu }(t)<+\infty . \end{aligned}$$

Hence, we can consider the continuous function $\lambda \in [0,+\infty ] \mapsto \lambda \mathfrak {e}_\varPhi (t;-\lambda ) \in \mathbb {R}^+$ and obtain (10) by Weierstrass theorem. $\square$

Let us give some examples of Bernstein functions and associated subordinators.

• We have already referred to the $\alpha$-stable subordinator, i.e. the one we get when we choose $\varPhi (\lambda )=\lambda ^{\alpha }$ for $\alpha \in (0,1)$. In such case, extensive informations on inverse $\alpha$-stable subordinators are given in [35]. As we stated before, we have in particular $\mathfrak {e}_\varPhi (t;\lambda )=E_\alpha (\lambda t^\alpha )$, where $E_\alpha$ is the one-parameter Mittag-Leffler function (see [15]). As a particular property, let us recall that, if we denote by $\sigma _\alpha$ the $\alpha$-stable subordinator and $g_\alpha$ the density of the random variable $\sigma _\alpha (1)$, then the inverse $\alpha$-stable subordinator $E_\varPhi (t)$ admits density

  $$\begin{aligned} f_\alpha (s;t)=\frac{t}{\beta }s^{-1-\frac{1}{\beta }}g_\alpha (ts^{-\frac{1}{\beta }}); \end{aligned}$$

• If we fix a constant $\theta >0$ and define $\varPhi (\lambda )=(\lambda +\theta )^\alpha -\theta ^\alpha$ we obtain the tempered $\alpha$-stable subordinator with tempering parameter $\theta >0$. Denoting by $\sigma _{\alpha ,\theta }(t)$ this subordinator, one can show that the density of the subordinator is given by

  $$\begin{aligned} g_{\alpha ,\theta }(s;t)=e^{-\theta s+t \theta ^\alpha }g_\alpha (s;t) \end{aligned}$$

  where $g_\alpha$ is the density of the $\alpha$-stable subordinator $\sigma _\alpha (t)$. An important property to recall is that the introduction of the tempering parameter implies the existence of all the moments of $\sigma _{\alpha ,\theta }(t)$ (while this is not true for $\sigma _\alpha$). Inverse tempered stable subordinators are studied for instance in [30]. Moreover, it can be shown that the Lévy tail $\overline{\nu }$ is given by

  $$\begin{aligned} \overline{\nu }(t)=\frac{\alpha \theta ^\alpha \varGamma (-\alpha ;t)}{\varGamma (1-\alpha )}, \end{aligned}$$

  where $\varGamma (\alpha ;x)=\int _{x}^{+\infty }t^{\alpha -1}e^{-t}dt$ is the upper incomplete Gamma function;

• For $\varPhi (\lambda )=\log (1+\lambda ^{\alpha })$ as $\alpha \in (0,1)$ we obtain the geometric $\alpha$-stable subordinator. From the form of the Bernstein function associated to the geometric $\alpha$-stable subordinator, one obtains (see [47, Theorem 2.6]) that the density $g_{G,\alpha }$ of the random variable $\sigma _{G,\alpha }(1)$ (where $\sigma _{G,\alpha }(t)$ is the geometric $\alpha$-stable subordinator) satisfies the following asymptotics:

  $$\begin{aligned} g_{G,\alpha }(s)\sim \frac{s^{\alpha -1}}{\varGamma (\alpha )} \qquad \text{ as } s \rightarrow 0^+;\\ g_{G,\alpha }(s)\sim 2\pi \sin \left( \frac{\alpha \pi }{2}\right) \varGamma (1+\alpha )s^{-\alpha -1} \qquad \text{ as } s \rightarrow +\infty . \end{aligned}$$

  Concerning the Lévy tail $\overline{\nu }$, it cannot be explicitly expressed, but it has been shown in [47, Theorem 2.5] that it satisfies the following asymptotic relation:

  $$\begin{aligned} \overline{\nu }(t)\sim \frac{t^{-\alpha }}{\varGamma (1-\alpha )} \qquad \text{ as } t \rightarrow +\infty . \end{aligned}$$

• If in the previous example we consider $\alpha =1$ we obtain the Gamma subordinator. In this specific case, one can obtain explicitly the Lévy tail $\overline{\nu }$ as (see, for instance, [47] and references therein)

  $$\begin{aligned} \overline{\nu }(t)=\varGamma (0;t). \end{aligned}$$

4 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 $\mathbf {m}$ and let us denote by $\mathcal {G}$ and $\mathcal {L}$ its backward and forward operators, respectively. Let us first focus on the backward equation.

4.1 Heuristic Derivation of the Strong Solution

Let us consider a Bernstein function $\varPhi$ and the Cauchy problem

$$\begin{aligned} {\left\{ \begin{array}{ll} \frac{\partial ^{\varPhi } }{\partial t^{\varPhi }}u(t,x)=\mathcal {G}u(t,x) &{} t>0, \ x \in E\\ u(0,x)=g(x) &{} x \in E. \end{array}\right. } \end{aligned}$$

(11)

Suppose that $g \in \ell ^2(\mathbf {m})$. Let us consider $(Q_n)_{n \in E}$ the family of orthonormal polynomials associated to N(t). Then we can decompose g as

$$\begin{aligned} g(x)=\sum _{n \in E}g_nQ_n(x) \end{aligned}$$

for some coefficients $(g_n)_{n \in 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)\varphi (x)$. If we substitute this relation in the first equation of (11) we obtain the following coupled equations

$$\begin{aligned} {\left\{ \begin{array}{ll} \mathcal {G}\varphi (x)=\lambda \varphi (x) &{} x \in E \\ \frac{\mathrm{d}^{\varPhi } }{\mathrm{d} t^{\varPhi }}T(t)=\lambda T(t) &{} t>0. \end{array}\right. } \end{aligned}$$

Concerning the first equation, let us observe that we need to set $\varphi (x)=Q_n(x)$ (up to a multiplicative constant) and $\lambda =\lambda _n$ for some $n \in E$. Let us also recall that $\lambda _n<0$. Concerning the second equation we get

$$\begin{aligned} T(t)=\mathfrak {e}_{\varPhi }(t;\lambda _n). \end{aligned}$$

and then we have for some $n \in E$

$$\begin{aligned} u(t,x)=Q_n(x)\mathfrak {e}_{\varPhi }(t;\lambda _n). \end{aligned}$$

We can also have solutions that are linear combinations of $Q_n\mathfrak {e}_\varPhi$. Let us suppose that we can consider eventually infinite linear combinations. Then we expect a solution of the form

$$\begin{aligned} u(t,x)=\sum _{n \ge 0} u_n Q_n(x)\mathfrak {e}_{\varPhi }(t;\lambda _n) \end{aligned}$$

for some coefficients $(u_n)_{n \in E}$. Finally, let us observe that

$$\begin{aligned} \sum _{n \in E}g_n Q_n(x)=g(x)=u(0,x)=\sum _{n \in E} u_n Q_n(x) \end{aligned}$$

then, since the components $(g_n)_{n \in E}$ are uniquely determined, then $u_n=g_n$ for any $n \in E$.

Finally, we expect the solution to be of the form

$$\begin{aligned} u(t,x)=\sum _{n \in E} g_nQ_n(x)\mathfrak {e}_{\varPhi }(t;\lambda _n). \end{aligned}$$

Now we want to formalize this reasoning.

4.2 The Backward Equation

Before working in the general case, we need to exploit what will be our fundamental solution. To do this, let us show the following Lemma.

Lemma 2

Let N(t) be a solvable birth–death process with state space E, generator $\mathcal {G}$, invariant measure $\mathbf {m}$ and family of associated classical orthogonal polynomials $(P_n)_{n \in E}$. Then the series

$$\begin{aligned} p_\varPhi (t,x;y)=m(x)\sum _{n \in E}\mathfrak {e}_\varPhi (t;\lambda _n)Q_n(x)Q_n(y), \end{aligned}$$

(12)

where $Q_n$ are the normalized orthogonal polynomials, absolutely converges for fixed $t \ge 0$ and $x,y \in E$.

Proof

Let us first observe that if E is finite, then the summation (12) is actually finite. Thus, let us consider the case in which $E=\mathbb {N}_0$. Let us denote by $\mathfrak {d}_n=\left\| P_n \right\| _{\ell ^2}$ and let us recall that the dual polynomials $\widetilde{P}_n$ exhibit orthogonality with respect to the measure $\widetilde{m}(x)=\frac{1}{\mathfrak {d}_x^2}$. However, since, if $E=\mathbb {N}_0$, N(t) is either an immigration-death process or a Meixner process, then the dual polynomials $\widetilde{P}_n$ coincide with the polynomials $P_n$ themselves. We have, by using the self-duality relation $P_n(x)=P_x(n)$,

$$\begin{aligned} p_\varPhi (t,x;y)&=m(x)\sum _{n=0}^{+\infty }\mathfrak {e}_\varPhi (t;\lambda _n)Q_n(x)Q_n(y)\\&=m(x)\sum _{n=0}^{+\infty }\widetilde{m}(n)\mathfrak {e}_\varPhi (t;\lambda _n)P_n(x)P_n(y)\\&=m(x)\sum _{n=0}^{+\infty }\widetilde{m}(n)\mathfrak {e}_\varPhi (t;\lambda _n)P_x(n)P_y(n). \end{aligned}$$

Now let us denote by $\mathrm{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

$$\begin{aligned} n_0=\lceil \max (\mathrm{root}(x) \cup \mathrm{root}(y))\rceil +1. \end{aligned}$$

To show the absolute convergence of the series in $p_\varPhi (t,x;y)$, we only need to show the absolute convergence of

$$\begin{aligned} \sum _{n=n_0}^{+\infty }\widetilde{m}(n)\mathfrak {e}_\varPhi (t;\lambda _n)P_x(n)P_y(n). \end{aligned}$$

Let us observe (see [38, 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 \ge n_0$ it holds $\mathrm {sign}(P_x(n)P_y(n))=(-1)^{x+y}$. So we get

$$\begin{aligned}&\sum _{n=n_0}^{+\infty }\left| \widetilde{m}(n)\mathfrak {e}_\varPhi (t;\lambda _n)P_x(n)P_y(n)\right| \\&\quad \le (-1)^{x+y}\sum _{n=n_0}^{+\infty }\widetilde{m}(n)P_x(n)P_y(n), \end{aligned}$$

where we have to observe that $\mathfrak {e}_\varPhi (t;\lambda _n)\le 1$ since $\lambda _n \le 0$. As before, the series at the right-hand side converges if and only if

$$\begin{aligned} \sum _{n=0}^{+\infty }\widetilde{m}(n)P_x(n)P_y(n) \end{aligned}$$

converges. However, by the dual orthogonal relation and the self-duality relation of Charlier and Meixner polynomials, we achieve

$$\begin{aligned} \sum _{n=0}^{+\infty }\widetilde{m}(n)P_x(n)P_y(n)=\frac{1}{m(x)}\delta _{x,y} \end{aligned}$$

for any $x,y \in \mathbb {N}$, concluding the proof. $\square$

Before exploiting the strong solution, let us show the normal convergence of some auxiliary series of functions.

Lemma 3

Let N(t) be a solvable birth–death process with state space $E=\mathbb {N}_0$, generator $\mathcal {G}$, invariant measure $\mathbf {m}$ and family of associated classical orthogonal polynomials $(P_n)_{n \ge 0}$. Let $g \in \ell ^2(\mathbf {m})$ such that $g(x)=\sum _{n \ge 0}g_nQ_n(x)$ for $x \in E$ and some constants $(g_n)_{n \ge 0}$, where $(Q_n)_{n \ge 0}$ are the normalized orthogonal polynomials. Then

• 1. For any $x \in E$ it holds $\sum _{n \ge 0}|g_nQ_n(x)| \le \frac{ \left\| g \right\| _{\ell ^2(\mathbf {m})}}{\sqrt{m(x)}}$;

• 2. For any fixed $x \in E$ the sum $\sum _{n \ge 0}\mathfrak {e}_\varPhi (t;\lambda _n)g_nQ_n(x)$ normally converges for $t \in [0,+\infty )$;

• 3. For any fixed $x \in E$ and $T_1>0$ the series $\sum _{n \ge 0}\lambda _n\mathfrak {e}_\varPhi (t;\lambda _n)g_nQ_n(x)$ normally converges for $t \in [T_1,+\infty )$.

Proof

Let us show property (1). Since $g(x)=\sum _{n \ge 0}g_nQ_n(x)$ and $Q_n$ is an orthonormal basis of $\ell ^2(\mathbf {m})$, it holds $\sum _{n \ge 0}g_n^2=\left\| g \right\| _{\ell ^2(\mathbf {m})}^2$. In particular it holds, by Cauchy-Schwartz inequality and self-duality relation for Charlier and Meixner polynomials

$$\begin{aligned} \sum _{n\ge 0}|g_nQ_n(x)|&=\sum _{n \ge 0}\sqrt{\widetilde{m}(n)}|g_nP_n(x)|\\&\le \left( \sum _{n\ge 0}\widetilde{m}(n)P^2_n(x)\right) ^{\frac{1}{2}}\left\| g \right\| _{\ell ^2(\mathbf {m})}\\&= \left( \sum _{n\ge 0}\widetilde{m}(n)P^2_x(n)\right) ^{\frac{1}{2}}\left\| g \right\| _{\ell ^2(\mathbf {m})}=\frac{\left\| g \right\| _{\ell ^2(\mathbf {m})}}{\sqrt{m(x)}}. \end{aligned}$$

To show property (2), let us just observe that $\mathfrak {e}_\varPhi (t;\lambda _n)\le 1$ since $\lambda _n\le 0$ and then

$$\begin{aligned} \sum _{n \ge 0}|\mathfrak {e}_\varPhi (t;\lambda _n)g_nQ_n(x)|\le \sum _{n \ge 0}|g_nQ_n(x)| \end{aligned}$$

where the series on the right-hand side is convergent and independent of $t\ge 0$.

Concerning property (3), by Eq. (10) we obtain for some constant $K(T_1)>0$, since $\mathfrak {e}_\varPhi (t;\lambda _n)$ is decreasing,

$$\begin{aligned} \sum _{n \ge 0}|\lambda _n\mathfrak {e}_\varPhi (t;\lambda _n)g_nQ_n(x)|\le \sum _{n \ge 0}|\lambda _n\mathfrak {e}_\varPhi (T_1;\lambda _n)g_nQ_n(x)| \le K(T_1)\sum _{n \ge 0}|g_nQ_n(x)|, \end{aligned}$$

concluding the proof. $\square$

Now we are ready to show that for the initial datum $g \in \ell ^2(\mathbf {m})$ our backward problem admits a solution.

Theorem 2

Let N(t) be a solvable birth–death process with state space E, generator $\mathcal {G}$, invariant measure $\mathbf {m}$ and family of associated classical orthogonal polynomials $(P_n)_{n \in E}$. Let $g \in \ell ^2(\mathbf {m})$ such that $g(x)=\sum _{n \in E}g_nQ_n(x)$ for $x \in E$ and some constants $(g_n)_{n \in E}$, where $(Q_n)_{n \in E}$ are the normalized orthogonal polynomials. Then the Cauchy problem

$$\begin{aligned} {\left\{ \begin{array}{ll} \frac{\partial ^{\varPhi } u}{\partial t^{\varPhi }}(t,y)=\mathcal {G}u(t,y) &{} t>0, \ y \in E \\ u(0,y)=g(y) &{} y \in E \end{array}\right. } \end{aligned}$$

(13)

admits a unique strong solution of the form

$$\begin{aligned} u(t,y)=\sum _{n \in E}\mathfrak {e}_\varPhi (t;\lambda _n)g_nQ_n(y). \end{aligned}$$

(14)

In particular $\sup _{t \ge 0}\left\| u(t,\cdot ) \right\| _{\ell ^2(\mathbf {m})}\le \left\| g \right\| _{\ell ^2(\mathbf {m})}$.

Finally $p_\varPhi (t,x;y)$ is the fundamental solution of (13), in the sense that it is the strong solution of (13) for $g(y)=\delta _x(y)$ and for any $g \in \ell ^2(\mathbf {m})$ it holds

$$\begin{aligned} u(t,y)=\sum _{x \in E}p_\varPhi (t,x;y)g(x). \end{aligned}$$

Before giving the proof of the Theorem, let us state formally what we mean as strong solution.

Definition 3

A function u(t, y) is a strong solution of the Cauchy problem (13) if:

• $\frac{\partial ^{\varPhi } u}{\partial t^{\varPhi }}(t,y)$ exists for any $t>0$ and $y \in E$;

• The convolution integral in $\frac{\partial ^{\varPhi } u}{\partial t^{\varPhi }}(t,y)$ is well-defined as a Bochner integral and the differentiation operator is a strong differentiation in $\ell ^2(\mathbf {m})$;

• The equations in (13) hold pointwise;

• $u(t;\cdot ) \in C([0,+\infty );\ell ^2(\mathbf {m}))$;

• $\frac{\partial ^{\varPhi } u}{\partial t^{\varPhi }}(t;\cdot ) \in C((0,+\infty );\ell ^2(\mathbf {m}))$.

Proof of Theorem 2

Let us first show that u(t, y) in the form of (14) is a strong solution for (13).

First of all, let us recall that, by definition of $\mathfrak {e}_\varPhi$ and $Q_n$, it holds

$$\begin{aligned} \mathcal {G}[\mathfrak {e}_\varPhi (t;\lambda _n)Q_ng_n](x)= & {} \mathfrak {e}_\varPhi (t;\lambda _n)g_n\mathcal {G}Q_n(x)\\= & {} \lambda _n \mathfrak {e}_\varPhi (t;\lambda _n)g_nQ_n(x)=\frac{\partial ^{\varPhi } }{\partial t^{\varPhi }}\mathfrak {e}_\varPhi (t;\lambda _n)g_nQ_n(x). \end{aligned}$$

Now let us observe that if E is finite then $n_E \in \mathbb {N}$ and we have that u(t, y) is a strong solution of (13) 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 (14) converges in $\ell ^2(\mathbf {m})$. To do this, define for $N \in \mathbb {N}$

$$\begin{aligned} u_N(t,y)=\sum _{n=0}^{N}\mathfrak {e}_\varPhi (t;\lambda _n)g_nQ_n(y). \end{aligned}$$

Now consider $N<M$ in $\mathbb {N}$ and observe that, being $\lambda _n\le 0$ and $\mathfrak {e}_\varPhi (t;\lambda _n)\le 1$, it holds

$$\begin{aligned} \left\| u_N(t,\cdot )-u_M(t,\cdot ) \right\| _{\ell ^2(\mathbf {m})}^2\le \sum _{n=N}^{M}g_n^2 \end{aligned}$$

thus, by Cauchy’s criterion, we know that the series converges in $\ell ^2(\mathbf {m})$.

Now let us denote

$$\begin{aligned} I_\nu (t)=\int _0^t\overline{\nu }(\tau )\mathrm{d}\tau \end{aligned}$$

that is increasing and non-negative.

Let us then observe that

$$\begin{aligned} \int _0^t(u(\tau ,y)-u(0+,y))\overline{\nu }(t-\tau )\mathrm{d}\tau =\int _0^t(u(\tau ,y)-u(0+,y))dI_\nu (t-\tau ). \end{aligned}$$

Since we have shown in Lemma 3 that the series defining u(t, y) normally converges for fixed $y \in E$, then we can use [43, Theorem 7.16] to write:

$$\begin{aligned} \int _0^t(u(\tau ,y)-u(0+,y))\overline{\nu }(t-\tau )\mathrm{d}\tau =\sum _{n=0}^{+\infty }\left( \int _0^t(\mathfrak {e}_\varPhi (\tau ;\lambda _n)-1)\overline{\nu }(t-\tau )\mathrm{d}\tau \right) Q_n(y)g_n. \end{aligned}$$

(15)

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 [43, Theorem 7.17]. Recalling (9), one has

$$\begin{aligned} \sum _{n=0}^{+\infty }\frac{\partial }{\partial t}\left( \int _0^t(\mathfrak {e}_\varPhi (\tau ;\lambda _n)-1)\overline{\nu }(t-\tau )\mathrm{d}\tau \right) Q_n(y)g_n&=\sum _{n=0}^{+\infty }\frac{\partial ^{\varPhi } }{\partial t^{\varPhi }}\mathfrak {e}_\varPhi (t;\lambda _n)Q_n(y)g_n\\ {}&=\sum _{n=0}^{+\infty }\lambda _n\mathfrak {e}_\varPhi (t;\lambda _n)Q_n(y)g_n \end{aligned}$$

that normally converges by statement (3) of Lemma 3.

Hence we obtain, differentiating on both sides in (15),

$$\begin{aligned} \frac{\partial ^{\varPhi } }{\partial t^{\varPhi }}u(t,y)=\sum _{n=0}^{+\infty }\frac{\partial ^{\varPhi } }{\partial t^{\varPhi }}\mathfrak {e}_\varPhi (t;\lambda _n)Q_n(y)g_n=\sum _{n=0}^{+\infty }\mathfrak {e}_\varPhi (t;\lambda _n)\mathcal {G}Q_n(y)g_n. \end{aligned}$$

Now let us recall that $\mathcal {G}=(b(x)-d(x))\nabla ^++d(x)\varDelta$, hence we have to show that we can exchange $\nabla ^+$ with the sign of series. This is due to the commutative property of normally convergent series. Indeed we have

$$\begin{aligned} \nabla ^+\sum _{n=0}^{+\infty }&\mathfrak {e}_\varPhi (t;\lambda _n) Q_n(y)g_n=\sum _{n=0}^{+\infty }\mathfrak {e}_\varPhi (t;\lambda _n) Q_n(y+1)g_n-\sum _{n=0}^{+\infty }\mathfrak {e}_\varPhi (t;\lambda _n) Q_n(y)g_n\\&=\lim _{N \rightarrow +\infty }\left( \sum _{n=0}^{N}\mathfrak {e}_\varPhi (t;\lambda _n) Q_n(y+1)g_n-\sum _{n=0}^{N}\mathfrak {e}_\varPhi (t;\lambda _n) Q_n(y)g_n\right) \\&=\lim _{N \rightarrow +\infty }\sum _{n=0}^{N}\mathfrak {e}_\varPhi (t;\lambda _n) \nabla ^+Q_n(y)g_n=\sum _{n=0}^{+\infty }\mathfrak {e}_\varPhi (t;\lambda _n) \nabla ^+Q_n(y)g_n \end{aligned}$$

where all the passages are justified by the fact that the two series $\sum _{n=0}^{+\infty }\mathfrak {e}_\varPhi (t;\lambda _n) Q_n(y+1)g_n$, $\sum _{n=0}^{+\infty }\mathfrak {e}_\varPhi (t;\lambda _n) Q_n(y)g_n$ both normally converge by Lemma 3. The same holds for $\varDelta$. Thus, we finally have

$$\begin{aligned} \frac{\partial ^{\varPhi } }{\partial t^{\varPhi }}u(t,y)=\sum _{n=0}^{+\infty }\mathfrak {e}_\varPhi (t;\lambda _n)\mathcal {G}Q_n(y)g_n=\mathcal {G}u(t,y) \end{aligned}$$

for any $t>0$. Concerning the initial datum, we have

$$\begin{aligned} u(0,y)=\sum _{n\in E}g_nQ_n(y)=g(y). \end{aligned}$$

Now let us show strong continuity of u in $[0,+\infty )$ and of $\frac{\partial ^\varPhi }{\partial t^\varPhi }u$ in $(0,+\infty )$. These properties are obvious as E is finite, thus let us suppose $E=\mathbb {N}_0$. Concerning the continuity of u, let us show it in $0^+$, since for any other $t \in (0,+\infty )$ the proof is analogous. We have

$$\begin{aligned} \left\| u(t,\cdot )-g(\cdot ) \right\| _{\ell ^2(\mathbf {m})}^2=\sum _{n=1}^{+\infty }(1-\mathfrak {e}_\varPhi (t;\lambda _n))^2g_n^2. \end{aligned}$$

Now fix $\varepsilon >0$. Since $(g_n)_{n \ge 0}$ belongs to $\ell ^2$, there exists $n(\varepsilon )\ge 0$ such that $\sum _{n=n(\varepsilon )}^{+\infty }g_n^2\le \varepsilon$. By using the fact that $(1-\mathfrak {e}_\varPhi (t;\lambda ))^2\le 1$ for any $\lambda <0$, we get

$$\begin{aligned} \left\| u(t,\cdot )-g(\cdot ) \right\| _{\ell ^2(\mathbf {m})}^2=\sum _{n=1}^{n(\varepsilon )}(1-\mathfrak {e}_\varPhi (t;\lambda _n))^{2}g_n^2+\varepsilon . \end{aligned}$$

Sending $t \rightarrow 0^+$ (using the fact that $\mathfrak {e}_\varPhi (t;\lambda _n)$ is continuous in t and that the first summation is finite) and then $\varepsilon \rightarrow 0^+$ we obtain strong continuity of u. Let us discuss the continuity of $\frac{\partial ^\varPhi }{\partial t^\varPhi }u$ in $(0,+\infty )$. To do this, let us consider $t_0 \in [t_1,t_2]$ with $t_1>0$. We have, for any $t \in [t_1,t_2]$, arguing as before and using (10),

$$\begin{aligned} \left\| \frac{\partial ^{\varPhi } }{\partial t^{\varPhi }}u(t,\cdot )-\frac{\partial ^{\varPhi } }{\partial t^{\varPhi }}u(t_0,\cdot ) \right\| _{\ell ^2(\mathbf {m})}^2&=\left\| \sum _{n \in E}\lambda _n(\mathfrak {e}_\varPhi (t;\lambda _n)-\mathfrak {e}_\varPhi (t_0;\lambda _n))Q_n(\cdot )g_n \right\| _{\ell ^2(\mathbf {m})}^2\\&=\sum _{n=0}^{+\infty }\lambda _n^2(\mathfrak {e}_\varPhi (t;\lambda _n)-\mathfrak {e}_\varPhi (t_0;\lambda _n))^2g^2_n\\&\le \sum _{n=0}^{n(\varepsilon )}\lambda _n^2(\mathfrak {e}_\varPhi (t;\lambda _n)-\mathfrak {e}_\varPhi (t_0;\lambda _n))^2g^2_n+K(t_1)\varepsilon . \end{aligned}$$

Thus, sending $t \rightarrow t_0$ (observing that the first sum is finite) and then $\varepsilon \rightarrow 0^+$ we obtain the desired continuity. Uniqueness follows easily from the fact that $(Q_n)_{n \ge 0}$ is an orthonormal system in $\ell ^2(\mathbf {m})$ hence the coefficients are unique.

Now let us show the bound of the $\ell ^2(\mathbf {m})$ norm. We have

$$\begin{aligned} \left\| u(t,\cdot ) \right\| _{\ell ^2(\mathbf {m})}^2=\sum _{n=0}^{+\infty }\mathfrak {e}_\varPhi ^2(t;\lambda _n)g_n^2\le \sum _{n=0}^{+\infty }g_n^2=\left\| g \right\| _{\ell ^2(\mathbf {m})}. \end{aligned}$$

We only have to show the second property of Definition 3. Let us first show that the convolution integral involved in the $\varPhi$-derivative $\frac{\partial ^{\varPhi } }{\partial t^{\varPhi }}u(t,\cdot )$ is a Bochner integral. To do this, we can use Bochner’s theorem (see [3, Theorem 1.1.4]). We just have to show that

$$\begin{aligned} \int _0^t \left\| u(\tau ,\cdot )-u(0+,\cdot ) \right\| _{\ell ^2(\mathbf {m})}\overline{\nu }(t-\tau )\mathrm{d}\tau <+\infty ; \end{aligned}$$

to do this, we will work with the square of this quantity. By Jensen’s inequality (see [44, Theorem 3.3]) we have

$$\begin{aligned}&\left( \int _0^t \left\| u(\tau ,\cdot )-u(0+,\cdot ) \right\| _{\ell ^2(\mathbf {m})}\overline{\nu }(t-\tau )\mathrm{d}\tau \right) ^2\\&\quad \le I_\nu (t)\int _0^t \left\| u(\tau ,\cdot )-u(0+,\cdot ) \right\| _{\ell ^2(\mathbf {m})}^2\overline{\nu }(t-\tau )\mathrm{d}\tau \\&\quad =I_\nu (t)\int _0^t \sum _{n=1}^{+\infty }(1-\mathfrak {e}_\varPhi (\tau ;\lambda _n))^2g_n^2\overline{\nu }(t-\tau )\mathrm{d}\tau , \end{aligned}$$

and then, by monotone convergence theorem and the fact that $0 \le 1-\mathfrak {e}_\varPhi (t;\lambda _n)\le 1$,

$$\begin{aligned}&\left( \int _0^t \left\| u(\tau ,\cdot )-u(0+,\cdot ) \right\| _{\ell ^2(\mathbf {m})}\overline{\nu }(t-\tau )\mathrm{d}\tau \right) ^2\\&\quad \le I_\nu (t) \sum _{n=1}^{+\infty }\int _0^t(1-\mathfrak {e}_\varPhi (\tau ;\lambda _n))^2\overline{\nu }(t-\tau )\mathrm{d}\tau g_n^2\\&\quad \le I^2_\nu (t)\left\| g \right\| _{\ell ^2(\mathbf {m})}<+\infty , \end{aligned}$$

concluding the proof of the statement.

Next, let us show that the differential involved in $\frac{\partial ^{\varPhi } }{\partial t^{\varPhi }}u(t,y)$ is a strong differential in $\ell ^2(\mathbf {m})$, that is to say

$$\begin{aligned}&\lim _{h \rightarrow 0}\left\| \frac{1}{h}\left( \int _0^{t+h}(u(t+h-\tau ,\cdot )-u(0+,\cdot ))\bar{\nu }(\tau )\mathrm{d}\tau \right. \right. \nonumber \\&\quad \left. \left. -\int _0^{t}(u(t-\tau ,\cdot )-u(0+,\cdot ))\bar{\nu }(\tau )\mathrm{d}\tau \right) -\sum _{n=0}^{+\infty }\lambda _n\mathfrak {e}_\varPhi (t;\lambda _n)Q_ng_n\right\| {\ell ^2(\mathbf {m})}=0. \end{aligned}$$

(16)

Let us first argue as $h>0$. First of all, let us rewrite the difference of the convolution integrals in the following way

$$\begin{aligned}&\int _0^{t+h}(u(t+h-\tau ,y)-u(0+,y))\bar{\nu }(\tau )\mathrm{d}\tau -\int _0^{t}(u(t-\tau ,y)-u(0+,y))\bar{\nu }(\tau )\mathrm{d}\tau \\&\quad =\int _0^{t}(u(t+h-\tau ,y)-u(t-\tau ,y))\bar{\nu }(\tau )\mathrm{d}\tau \\&\qquad +\int _{t}^{t+h}(u(t+h-\tau ,y)-u(0+,y))\bar{\nu }(\tau )\mathrm{d}\tau , \end{aligned}$$

thus we can use the triangular inequality to achieve

$$\begin{aligned}&\left\| \frac{1}{h}\left( \int _0^{t+h}(u(t+h-\tau ,\cdot )-u(0+,\cdot ))\bar{\nu }(\tau )\mathrm{d}\tau \right. \right. \\&\qquad \left. \left. -\int _0^{t}(u(t-\tau ,\cdot )-u(0+,\cdot ))\bar{\nu }(\tau )\mathrm{d}\tau \right) -\sum _{n=0}^{+\infty }\lambda _n\mathfrak {e}_\varPhi (t;\lambda _n)Q_ng_n\right\| {\ell ^2(\mathbf {m})}\\&\quad \le \left\| \frac{1}{h}\int _0^{t}(u(t+h-\tau ,\cdot )-u(t-\tau ,\cdot ))\bar{\nu }(\tau )\mathrm{d}\tau -\sum _{n=0}^{+\infty }\lambda _n\mathfrak {e}_\varPhi (t;\lambda _n)Q_ng_n \right\| _{\ell ^2(\mathbf {m})}\\&\qquad +\left\| \frac{1}{h}\int _{t}^{t+h}(u(t+h-\tau ,\cdot )-u(0+,\cdot ))\bar{\nu }(\tau )\mathrm{d}\tau \right\| _{\ell ^2(\mathbf {m})}. \end{aligned}$$

Let us first consider the second summand. By Bochner’s theorem and the fact that $\overline{\nu }(t)$ is decreasing we have

$$\begin{aligned}&\left\| \frac{1}{h}\int _{t}^{t+h}(u(t+h-\tau ,\cdot )-u(0+,\cdot ))\bar{\nu }(\tau )\mathrm{d}\tau \right\| _{\ell ^2(\mathbf {m})} \\&\quad \le \overline{\nu }(t)\frac{1}{h}\int _{0}^{h}\left\| u(\tau ,\cdot )-u(0+,\cdot ) \right\| _{\ell ^2(\mathbf {m})}\mathrm{d}\tau . \end{aligned}$$

Since we have already shown strong continuity of $t \in [0,+\infty ) \mapsto u(t,\cdot )\in \ell ^2(\mathbf {m})$, we know that $\tau \mapsto \left\| u(\tau ,\cdot )-u(0+,\cdot ) \right\| _{\ell ^2(\mathbf {m})}$ is a continuous function, hence, by the integral mean theorem, there exists $\theta (h) \in [0,h]$ such that

$$\begin{aligned} \frac{1}{h}\int _{0}^{h}\left\| u(\tau ,\cdot )-u(0+,\cdot ) \right\| _{\ell ^2(\mathbf {m})}\mathrm{d}\tau =\left\| u(\theta (h),\cdot )-u(0+,\cdot ) \right\| _{\ell ^2(\mathbf {m})}. \end{aligned}$$

Taking the limit as $h \rightarrow 0^+$ we conclude that

$$\begin{aligned} \lim _{h \rightarrow 0^+}\frac{1}{h}\int _{0}^{h}\left\| u(\tau ,\cdot )-u(0+,\cdot ) \right\| _{\ell ^2(\mathbf {m})}\mathrm{d}\tau =0, \end{aligned}$$

and then

$$\begin{aligned} \lim _{h \rightarrow 0^+}\left\| \frac{1}{h}\int _{t}^{t+h}(u(t+h-\tau ,\cdot )-u(0+,\cdot ))\bar{\nu }(\tau )\mathrm{d}\tau \right\| _{\ell ^2(\mathbf {m})}=0. \end{aligned}$$

(17)

Now let us consider the first summand. By definition of u(t, y) we have

$$\begin{aligned}&\int _0^{t}(u(t+h-\tau ,y)-u(t-\tau ,y))\bar{\nu }(\tau )\mathrm{d}\tau \\&\quad =\int _0^{t}\sum _{n=1}^{+\infty }(\mathfrak {e}_\varPhi (t+h-\tau ;\lambda _n)-\mathfrak {e}_\varPhi (t-\tau ;\lambda _n))g_nQ_n(y)\overline{\nu }(\tau )\mathrm{d}\tau . \end{aligned}$$

Let us recall that for any $t,s>0$ it holds $|\mathfrak {e}_\varPhi (t;\lambda _n)-\mathfrak {e}_\varPhi (s;\lambda _n)|\le 1$. By Lemma 3, property 1, we have

$$\begin{aligned} \sum _{n=1}^{+\infty }|\mathfrak {e}_\varPhi (t+h-\tau ;\lambda _n)-\mathfrak {e}_\varPhi (t-\tau ;\lambda _n)||g_nQ_n(y)|\le \frac{\left\| g \right\| _{\ell ^2(\mathbf {m})}}{\sqrt{m(y)}}, \end{aligned}$$

thus, being $\overline{\nu }$ integrable in (0, t), we can use dominated convergence theorem to achieve

$$\begin{aligned}&\int _0^{t}(u(t+h-\tau ,y)-u(t-\tau ,y))\bar{\nu }(\tau )\mathrm{d}\tau \\&\quad =\sum _{n=1}^{+\infty }\left( \int _0^{t}(\mathfrak {e}_\varPhi (t+h-\tau ;\lambda _n)-\mathfrak {e}_\varPhi (t-\tau ;\lambda _n))\overline{\nu }(\tau )\mathrm{d}\tau \right) g_nQ_n(y). \end{aligned}$$

Observing that

$$\begin{aligned} \sum _{n=1}^{+\infty }\left| \frac{1}{h}\int _0^{t}(\mathfrak {e}_\varPhi (t+h-\tau ;\lambda _n)-\mathfrak {e}_\varPhi (t-\tau ;\lambda _n))\overline{\nu }(\tau )\mathrm{d}\tau \right| |g_nQ_n(y)|\le \frac{I_\nu (t)\left\| g \right\| _{\ell ^2(\mathbf {m})}}{h\sqrt{m(y)}} \end{aligned}$$

and

$$\begin{aligned} \sum _{n=1}^{+\infty }|\lambda _n \mathfrak {e}_\varPhi (t;\lambda _n)||g_nQ_n(y)|\le \frac{K(t)\left\| g \right\| _{\ell ^2(\mathbf {m})}}{\sqrt{m(y)}} \end{aligned}$$

for a suitable K(t), we have, being all the involved series absolutely convergent,

$$\begin{aligned}&\sum _{n=1}^{+\infty }\left( \frac{1}{h}\int _0^{t}(\mathfrak {e}_\varPhi (t+h-\tau ;\lambda _n)-\mathfrak {e}_\varPhi (t-\tau ;\lambda _n))\overline{\nu }(\tau )\mathrm{d}\tau \right) g_nQ_n(y)\nonumber \\&\qquad -\sum _{n=1}^{+\infty }\lambda _n\mathfrak {e}_\varPhi (t;\lambda _n)g_nQ_n(y)\nonumber \\&\quad =\sum _{n=1}^{+\infty }\left( \frac{1}{h}\int _0^{t}(\mathfrak {e}_\varPhi (t+h-\tau ;\lambda _n)-\mathfrak {e}_\varPhi (t-\tau ;\lambda _n))\overline{\nu }(\tau )\mathrm{d}\tau \right. \nonumber \\&\qquad \left. -\lambda _n\mathfrak {e}_\varPhi (t;\lambda _n)\right) g_nQ_n(y). \end{aligned}$$

(18)

To show that the previous series converges in $\ell ^2(\mathbf {m})$, let us define

$$\begin{aligned} S_N= & {} \sum _{n=1}^{N}\left( \frac{1}{h}\int _0^{t}(\mathfrak {e}_\varPhi (t+h-\tau ;\lambda _n)-\mathfrak {e}_\varPhi (t-\tau ;\lambda _n))\overline{\nu }(\tau )\mathrm{d}\tau \right. \\&\left. -\lambda _n\mathfrak {e}_\varPhi (t;\lambda _n)\right) g_nQ_n(y) \end{aligned}$$

and observe that, for any $M>N$,

$$\begin{aligned} \left\| S_M-S_N \right\| _{\ell ^2(\mathbf {m})}^2&=\sum _{n=N+1}^{M}\Bigg (\frac{1}{h}\int _0^{t}(\mathfrak {e}_\varPhi (t+h-\tau ;\lambda _n)\\&\quad -\mathfrak {e}_\varPhi (t-\tau ;\lambda _n))\overline{\nu }(\tau )\mathrm{d}\tau -\lambda _n\mathfrak {e}_\varPhi (t;\lambda _n)\Bigg )^2g^2_n\\&\le \sum _{n=N+1}^{M}\left( \left( \frac{1}{h}\int _0^{t}(\mathfrak {e}_\varPhi (t+h-\tau ;\lambda _n)\right. \right. \\&\quad \left. \left. -\mathfrak {e}_\varPhi (t-\tau ;\lambda _n))\overline{\nu }(\tau )\mathrm{d}\tau \right) ^2+(\lambda _n\mathfrak {e}_\varPhi (t;\lambda _n))^2\right) g^2_n\\&\le \left( \frac{I^2_\nu (t)}{h^2}+K^2(t)\right) \sum _{n=N+1}^{M}g_n^2 \end{aligned}$$

for a suitable K(t). Being $\sum _{n=0}^{+\infty }g_n^2$ convergent, by Cauchy’s criterion we know that for any $\varepsilon >0$ there exists $N_0>0$ such that, for any $N_0\le N<M$, it holds $\sum _{n=N+1}^{M}g_n^2<\frac{\varepsilon }{\left( \frac{I^2_\nu (t)}{h^2}+K^2(t)\right) }$. In particular, choosing $N_0\le N<M$ we have that $\left\| S_M-S_N \right\| _{\ell ^2(\mathbf {m})}<\varepsilon$, thus $S_N$ converges in $\ell ^2(\mathbf {m})$. In particular, this means that we can use Parseval’s identity on the series in the right-hand side of Eq. (18), hence

$$\begin{aligned}&\left\| \frac{1}{h}\int _0^{t}(u(t+h-\tau ,\cdot )-u(t-\tau ,\cdot ))\bar{\nu }(\tau )\mathrm{d}\tau -\sum _{n=1}^{+\infty }\lambda _n\mathfrak {e}_\varPhi (t;\lambda _n)Q_ng_n \right\| _{\ell ^2(\mathbf {m})}^2\\&=\sum _{n=1}^{+\infty }\left( \frac{1}{h}\int _0^{t}(\mathfrak {e}_\varPhi (t+h-\tau ;\lambda _n)-\mathfrak {e}_\varPhi (t-\tau ;\lambda _n))\overline{\nu }(\tau )\mathrm{d}\tau -\lambda _n\mathfrak {e}_\varPhi (t;\lambda _n)\right) ^2g^2_n\\&=\sum _{n=1}^{+\infty }\left( \frac{1}{h}\left( \int _0^{t+h}(\mathfrak {e}_\varPhi (t+h-\tau ;\lambda _n)-1)\overline{\nu }(\tau )\mathrm{d}\tau \right. \right. \\&\qquad \left. \left. -\int _0^{t}(\mathfrak {e}_\varPhi (t-\tau ;\lambda _n)-1)\overline{\nu }(\tau )\mathrm{d}\tau \right) \right. \\&\qquad \left. -\frac{1}{h}\int _{t}^{t+h}(\mathfrak {e}_\varPhi (t+h-\tau ;\lambda _n)-1)\overline{\nu }(\tau )\mathrm{d}\tau -\lambda _n\mathfrak {e}_\varPhi (t;\lambda _n)\right) ^2g^2_n\\&2\le \sum _{n=1}^{+\infty }\left( \left( \frac{I_n(t+h)-I_n(t)}{h}-\lambda _n\mathfrak {e}_\varPhi (t;\lambda _n)\right) ^2\right. \\&\qquad \left. +\overline{\nu }^2(t)\left( \frac{1}{h}\int _{0}^{h}(\mathfrak {e}_\varPhi (s;\lambda _n)-1)^2\mathrm{d}s\right) \right) g^2_n, \end{aligned}$$

where we set

$$\begin{aligned} I_n(t)=\int _0^{t}(\mathfrak {e}_\varPhi (t-\tau ;\lambda _n)-1)\overline{\nu }(\tau )\mathrm{d}\tau \end{aligned}$$

and we used Jensen’s inequality.

Observe that $I_n(t)$ is a $C^1$ function with $I_n'(t)=\lambda _n\mathfrak {e}_\varPhi (t;\lambda _n)$. By Lagrange’s theorem there exists a function $\theta _n(h)\in [t,t+h]$, for each $n \ge 1$, such that

$$\begin{aligned} \frac{I_n(t+h)-I_n(t)}{h}=I'_n(\theta _n(h))=\lambda _n\mathfrak {e}_\varPhi (\theta _n(h);\lambda _n), \end{aligned}$$

and then we obtain

$$\begin{aligned}&\left\| \frac{1}{h}\int _0^{t+h}(u(t+h-\tau ,\cdot )-u(t-\tau ,\cdot ))\bar{\nu }(\tau )\mathrm{d}\tau \right. \\&\qquad \left. -\sum _{n=0}^{+\infty }\lambda _n\mathfrak {e}_\varPhi (t;\lambda _n)Q_ng_n \right\| _{\ell ^2(\mathbf {m})}^2\\&\quad 2\le \sum _{n=1}^{+\infty }\left( \left( \lambda _n\mathfrak {e}_\varPhi (\theta _n(h);\lambda _n)-\lambda _n\mathfrak {e}_\varPhi (t;\lambda _n)\right) ^2\right. \\&\qquad \left. +\overline{\nu }^2(t)\left( \frac{1}{h}\int _{0}^{h}(\mathfrak {e}_\varPhi (s;\lambda _n)-1)^2\mathrm{d}s\right) \right) g^2_n. \end{aligned}$$

Let us show that the series in the previous equation normally converges for fixed $t>0$. To do this, just observe that

$$\begin{aligned}&\sum _{n=1}^{+\infty }\left( (\lambda _n\mathfrak {e}_\varPhi (\theta _n(h);\lambda _n)-\lambda _n\mathfrak {e}_\varPhi (t;\lambda _n))^2\right. \\&\qquad \left. +\overline{\nu }^2(t)\left( \frac{1}{h}\int _{0}^{h}(\mathfrak {e}_\varPhi (s;\lambda _n)-1)^2\mathrm{d}s\right) \right) g^2_n\\&\quad \le (2K^2(t)+\overline{\nu }^2(t))\left\| g \right\| _{\ell ^2(\mathbf {m})}. \end{aligned}$$

Thus, we can take the limit inside the summation sign to obtain

$$\begin{aligned}&\lim _{h \rightarrow 0^+}\sum _{n=1}^{+\infty }\left( (\lambda _n\mathfrak {e}_\varPhi (\theta _n(h);\lambda _n)-\lambda _n\mathfrak {e}_\varPhi (t;\lambda _n))^2\right. \\&\qquad \left. +\overline{\nu }^2(t)\left( \frac{1}{h}\int _{0}^{h}(\mathfrak {e}_\varPhi (s;\lambda _n)-1)^2\mathrm{d}s\right) \right) g^2_n\\&\quad =\sum _{n=1}^{+\infty }\left( \lambda _n^2\lim _{h \rightarrow 0^+}(\mathfrak {e}_\varPhi (\theta _n(h);\lambda _n)-\mathfrak {e}_\varPhi (t;\lambda _n))^2\right. \\&\qquad \left. +\overline{\nu }^2(t)\lim _{h \rightarrow 0^+}\left( \frac{1}{h}\int _{0}^{h}(\mathfrak {e}_\varPhi (s;\lambda _n)-1)^2\mathrm{d}s\right) \right) g^2_n=0, \end{aligned}$$

where we used the fact that $\mathfrak {e}_\varPhi (t;\lambda _n)$ is continuous. Last relation implies

$$\begin{aligned}&\lim _{h \rightarrow 0^+}\left\| \frac{1}{h}\right. \left. \int _0^{t+h}(u(t+h-\tau ,\cdot )-u(t-\tau ,\cdot ))\bar{\nu }(\tau )\mathrm{d}\tau \right. \nonumber \\&\quad \left. -\sum _{n=0}^{+\infty }\lambda _n\mathfrak {e}_\varPhi (t;\lambda _n)Q_ng_n \right\| _{\ell ^2(\mathbf {m})}^2=0. \end{aligned}$$

(19)

Combining Eqs. (17) and (19) we finally obtain

$$\begin{aligned}&\lim _{h \rightarrow 0^+}\left\| \frac{1}{h}\left( \int _0^{t+h}(u(t+h-\tau ,\cdot )-u(0+,\cdot ))\bar{\nu }(\tau )\mathrm{d}\tau \right. \right. \nonumber \\&\quad \left. \left. -\int _0^{t}(u(t-\tau ,\cdot )-u(0+,\cdot ))\bar{\nu }(\tau )\mathrm{d}\tau \right) -\sum _{n=0}^{+\infty }\lambda _n\mathfrak {e}_\varPhi (t;\lambda _n)Q_ng_n\right\| {\ell ^2(\mathbf {m})}=0. \end{aligned}$$

(20)

Concerning the case $h<0$, after observing that

$$\begin{aligned}&\int _0^{t+h}(u(t+h-\tau ,y)-u(0+,y))\bar{\nu }(\tau )\mathrm{d}\tau -\int _0^{t}(u(t-\tau ,y)-u(0+,y))\bar{\nu }(\tau )\mathrm{d}\tau \\&\quad =\int _0^{t+h}(u(t+h-\tau ,y)-u(t-\tau ,y))\bar{\nu }(\tau )\mathrm{d}\tau \\&\qquad +\int _{t+h}^{t}(u(t-\tau ,y)-u(0+,y))\bar{\nu }(\tau )\mathrm{d}\tau , \end{aligned}$$

one can conclude with the same argument as in the case $h>0$ that

$$\begin{aligned}&\lim _{h \rightarrow 0^-}\left\| \frac{1}{h}\left( \int _0^{t+h}(u(t+h-\tau ,\cdot )-u(0+,\cdot ))\bar{\nu }(\tau )\mathrm{d}\tau \right. \right. \nonumber \\&\quad \left. \left. -\int _0^{t}(u(t-\tau ,\cdot )-u(0+,\cdot ))\bar{\nu }(\tau )\mathrm{d}\tau \right) -\sum _{n=0}^{+\infty }\lambda _n\mathfrak {e}_\varPhi (t;\lambda _n)Q_ng_n\right\| {\ell ^2(\mathbf {m})}=0. \end{aligned}$$

(21)

Combining Eqs. (20) and (21) we finally obtain Eq. (16).

Now let us show that $p_\varPhi (t,x;y)$ is the fundamental solution. Let us consider a function $g \in \ell ^2(\mathbf {m})$. Then we have

$$\begin{aligned} \sum _{x \in E}p_\varPhi (t,x;y)g(x)&=\sum _{x \in E}m(x)\left( \sum _{n \in E}\mathfrak {e}_\varPhi (t;\lambda _n)Q_n(x)Q_n(y)\right) g(x)\\&=\sum _{n\in E}Q_n(y)\mathfrak {e}_\varPhi (t;\lambda _n)\sum _{x \in E}m(x)Q_n(x)g(x)\\&=\sum _{n \in E}Q_n(y)\mathfrak {e}_\varPhi (t;\lambda _n)g_n=u(t,y) \end{aligned}$$

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 normally convergent in compact sets containing t.

Finally, let us observe that if, for some $z \in E$, $g(x) =\delta _z(x)$, then

$$\begin{aligned} u(t,y)=\sum _{x \in E}p_\varPhi (t,x;y)\delta _z(x)=p_\varPhi (t,z;y), \end{aligned}$$

concluding the proof. $\square$

4.3 The Forward Equation

Now let us apply the same strategy to study the Cauchy problem associated with $\mathcal {L}$.

Theorem 3

Let $N(t)$ be a solvable birth–death process with state space $E,\,{ forwardoperator}\mathcal {L},\,{ invariantmeasure}\mathbf {m}{ andfamilyofassociatedclassicalorthogonalpolynomials}(P_n)_{n \in E}.{ Let}f/m \in \ell ^2(\mathbf {m}){ suchthat}f(x)=m(x)\sum _{n \in E}f_nQ_n(x){ for}x \in E{ andsomeconstants}(f_n)_{n \in E},\,{ where}(Q_n)_{n \in E}{}$ are the normalized orthogonal polynomials. Then the Cauchy problem

$$\begin{aligned} {\left\{ \begin{array}{ll} \frac{\partial ^{\varPhi } v}{\partial t^{\varPhi }}(t,x)=\mathcal {L}v(t,x) &{} t>0, \ x \in E \\ v(0,x)=f(x) &{} x \in E \end{array}\right. } \end{aligned}$$

(22)

admits a unique strong solution of the form

$$\begin{aligned} v(t,x)=m(x)\sum _{n \in E}\mathfrak {e}_\varPhi (t;\lambda _n)f_nQ_n(x), \end{aligned}$$

(23)

such that $\sup _{t \ge 0}\left\| v(t,\cdot ) \right\| _{\ell ^2(\mathbf {m})}\le \left\| f/m \right\| _{\ell ^2(\mathbf {m})}.{ Finally}p_\varPhi (t,x;y)$ is the fundamental solution of (13), in the sense that it is the strong solution of (13) for $f(x)=\delta _y(x){ andforany}f/m \in \ell ^2(\mathbf {m})$ it holds

$$\begin{aligned} v(t,x)=\sum _{y \in E}p_\varPhi (t,x;y)f(y). \end{aligned}$$

Proof

Let us first observe, by Lemma 1

$$\begin{aligned}&\mathcal {L}m(x)\mathfrak {e}_\varPhi (t;\lambda _n)f_nQ_n(x)= \mathfrak {e}_\varPhi (t;\lambda _n)f_n\mathcal {L}_{z \rightarrow x}(m(z)Q_n(z))(x)\\&\quad =\lambda _n\mathfrak {e}_\varPhi (t;\lambda _n)f_nm(x)Q_n(x)=\frac{\partial ^{\varPhi } }{\partial t^{\varPhi }}m(x)\mathfrak {e}_\varPhi (t;\lambda _n)f_nQ_n(x). \end{aligned}$$

Thus, the exact same strategy of the proof of Theorem 2 leads to formula (23) and uniqueness follows as before. Concerning the estimate on the norm, it holds, since $\mathbf {m}{ isaprobabilitymeasureandthen}{m(x)\le 1}{}$,

$$\begin{aligned} \left\| v(t,x) \right\| _{\ell ^2(\mathbf {m})}^2&=\sum _{x \in E}m(x)\left( m(x)\sum _{n \in E}\mathfrak {e}_\varPhi (t;\lambda _n)f_nQ_n(x)\right) ^2\\&\le \sum _{x \in E}m(x)\left( \sum _{n\in E}f_n\mathfrak {e}_\varPhi (t;\lambda _n)Q_n(x)\right) ^2\\&=\sum _{n\in E}f^2_n\mathfrak {e}^2_\varPhi (t;\lambda _n)\le \left\| f/m \right\| _{\ell ^2(\mathbf {m})}^2. \end{aligned}$$

Moreover, let us observe that

$$\begin{aligned} \sum _{y \in E}p_\varPhi (t,x;y)f(y)&=m(x)\sum _{y \in E}\left( \sum _{n \in E}\mathfrak {e}_\varPhi (t;\lambda _n)Q_n(x)Q_n(y)\right) \frac{f(y)}{m(y)}m(y)\\&=m(x)\sum _{n \in E}\mathfrak {e}_\varPhi (t;\lambda _n)Q_n(x)\sum _{y \in E}Q_n(y)\frac{f(y)}{m(y)}m(y)\\&=m(x)\sum _{n \in E}\mathfrak {e}_\varPhi (t;\lambda _n)f_nQ_n(x)=v(t,x). \end{aligned}$$

Finally, observe that if, for some fixed $z\in E$, $f(y)=\delta _z(y)$, then obviously $f/m \in \ell ^2(\mathbf {m})$ and we have

$$\begin{aligned} v(t,x)=\sum _{y \in E}p_\varPhi (t,x;y)f(y)=\sum _{y \in E}p_\varPhi (t,x;y)\delta _z(y)=p_\varPhi (t,x;z) \end{aligned}$$

concluding the proof. $\square$

Remark 1

One also obtains the following estimate on the norm:

$$\begin{aligned} \sup _{t \ge 0}\left\| v(t,\cdot )/m(\cdot ) \right\| _{\ell ^2(\mathbf {m})}\le \left\| f/m \right\| _{\ell ^2(\mathbf {m})}. \end{aligned}$$

Next step is to identify some processes such that $p_\varPhi (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 (14) and (23).

5 Non-local Solvable Birth–Death Processes

Let us now consider a solvable birth–death process $N{ andthesubordinator}\sigma _\varPhi { associatedtotheBernsteinfunction}\varPhi ,\,{ withitsinverse}E_\varPhi .{ Letussuppose}N{ and}E_\varPhi$ are independent.

Definition 4

The non-local solvable birth–death process induced by $N{ and}\varPhi$ is defined as

$$\begin{aligned} N_\varPhi (t):=N(E_\varPhi (t)) \end{aligned}$$

(24)

for any $t \ge 0,\,{ where}E_\varPhi (t){ isindependentof}N(t)$. Moreover, we define its transition probability function

$$\begin{aligned} p_\varPhi (t,x;y):=\mathbb {P}(N_\varPhi (t)=x|N_\varPhi (0)=y). \end{aligned}$$

We used the same notation for the transition probability function and the fundamental solution of the Cauchy problems (14) and (23). Indeed, we can show the following result.

Theorem 4

The transition probability function $p_\varPhi (t,x;y)$ coincides with the series defined in Lemma 2.

Proof

Let us first observe that $N_\varPhi (0)=N(0){ almostsurely}({ since}E_\varPhi (0)=0$ almost surely), thus, by a simple conditioning argument, we have

$$\begin{aligned} p_\varPhi (t,x;y)=\int _0^{+\infty }p(s,x;y)f_\varPhi (s,t)\mathrm{d}s. \end{aligned}$$

Now let us recall, by Theorem 1, that

$$\begin{aligned} p(s,x;y)=m(x)\sum _{n \in E}e^{\lambda _n t}Q_n(x)Q_n(y) \end{aligned}$$

thus we have

$$\begin{aligned} p_\varPhi (t,x;y)=m(x)\int _0^{+\infty }\sum _{n \in E}e^{\lambda _n s}Q_n(x)Q_n(y)f_\varPhi (s,t)\mathrm{d}s. \end{aligned}$$

If $E{ isfinite},\,{ weconcludetheproof}.{ Soletussuppose}E=\mathbb {N}_0$. Recalling the notation in the proof of Lemma 2, let us define

$$\begin{aligned} n_0=\lceil \max (\mathrm{root}(x) \cup \mathrm{root}(y))\rceil +1. \end{aligned}$$

and let us split the summation in two parts. We have

$$\begin{aligned} p_\varPhi (t,x;y)&=m(x)\int _0^{+\infty }\sum _{n=0}^{n_0}Q_n(x)Q_n(y)e^{\lambda _n s}f_\varPhi (s,t)\mathrm{d}s\\&\quad +m(x)\int _0^{+\infty }\sum _{n={n_0+1}}^{+\infty }Q_n(x)Q_n(y)e^{\lambda _n s}f_\varPhi (s,t)\mathrm{d}s. \end{aligned}$$

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)=\widetilde{m}(n)P_x(n)P_y(n)$). Thus, we have

$$\begin{aligned} p_\varPhi (t,x;y)=&m(x)\sum _{n=0}^{n_0}Q_n(x)Q_n(y)\int _0^{+\infty }e^{\lambda _n s}f_\varPhi (s,t)\mathrm{d}s\\&+m(x)\sum _{n={n_0+1}}^{+\infty }Q_n(x)Q_n(y)\int _0^{+\infty }e^{\lambda _n s}f_\varPhi (s,t)\mathrm{d}s\\&=m(x)\sum _{n=0}^{+\infty }Q_n(x)Q_n(y)\int _0^{+\infty }e^{\lambda _n s}f_\varPhi (s,t)\mathrm{d}s. \end{aligned}$$

Finally, let us recall that, by definition

$$\begin{aligned} \int _0^{+\infty }e^{\lambda _n s}f_\varPhi (s,t)\mathrm{d}s=\mathbb {E}[e^{\lambda _n E_\varPhi (t)}]=\mathfrak {e}_\varPhi (t;\lambda _n), \end{aligned}$$

concluding the proof. $\square$

By using this result, we can achieve some stochastic representation of the solutions of the Cauchy problems (14) and (23). Indeed we have the following result.

Proposition 5

Let $N_\varPhi (t)$ be the non-local solvable birth–death process associated to $N(t){ and}\varPhi .{ Suppose}N(t){ admitsstatespace}E$. Then

• 1. For $g \in \ell ^2(\mathbf {m}){ thefunction}u(t,y)=\mathbb {E}_y[g(N_\varPhi (t))] ({ where}\mathbb {E}_y[\cdot ]=\mathbb {E}[\cdot |N_\varPhi (0)=y]$) is strong solution of (14);

• 2. For $f/m \in \ell ^2(\mathbf {m}){ suchthat}f \ge 0{ and}\sum _{x \in E}f(x)=1,\,{ denotingby}\mathbb {P}_{f}{} { theprobabilitymeasureobtainedby}\mathbb {P}{ conditioningwiththefactthat}N_\varPhi (0){ admitsdistribution}f,\,{ thenthefunction}v(t,x)=\mathbb {P}_{f}(N_\varPhi (t)=x)$ is strong solution of (23).

Proof

To show assertion $1$, let us observe that

$$\begin{aligned} u(t,y)=\mathbb {E}_y[g(N_\varPhi (t))]=\sum _{x \in E}g(x)p_\varPhi (t,x;y), \end{aligned}$$

obtaining that $u(t,y)$ is the strong solution of (13) by Theorem 2.

Concerning assertion $2$, we have

$$\begin{aligned} v(t,x)=\mathbb {P}_f(N_\varPhi (t)=x)=\sum _{y \in E}f(y)p_\varPhi (t,x;y) \end{aligned}$$

obtaining that $v(t,x)$ is the strong solution of (22) by Theorem 3. $\square$

Finally, this proposition allows us to obtain the invariant distribution of $N_\varPhi (t)$ and show that it is also the limit distribution.

Corollary 3

Let $N_\varPhi (t)$ be the non-local solvable birth–death process associated to $N(t){ and}\varPhi .{ Suppose}N(t){ admitsstatespace}E$. Then

• 1. If $N_\varPhi (0){ admitsdistribution}\mathbf {m},\,{ then}N_\varPhi (t){ admitsdistribution}\mathbf {m}{ forany}t \ge 0$;

• 2. If $N_\varPhi (0){ admitsdistribution}f{ suchthat}f/m \in \ell ^2(\mathbf {m}),\,{ then}\lim _{t \rightarrow +\infty }\mathbb {P}_f(N_\varPhi (t)=x)=m(x)$.

Proof

Let us first show property $(1).{ Todothis},\,{ letusobservethat}1 \in \ell ^2(\mathbf {m}){ since}\mathbf {m}$ is a probability measure. Thus, recall, by Proposition 5, that $v(t,x)=\mathbb {P}_{\mathbf {m}}(N_\varPhi (t)=x)$ is a strong solution of (22). Now we need to determine $m_n{ suchthat}\sum _{n \in E}m_nQ_n(x)=1.{ Letusrecallthat}Q_0(x)=1{ while}\deg (Q_n(x))=n{ forany}n=1,\dots ,n_E,\,{ thuswehave}m_0=1{ and}m_n=0{ forany}n=1,\dots ,n_E.{ Moreover},\,{ letusrecallthat}\lambda _0=0,\,{ since}1 \in Ker(\mathcal {G})$. Hence we have, by Theorem 3,

$$\begin{aligned} v(t,x)=m(x)\sum _{n \in E}\mathfrak {e}_\varPhi (t;\lambda _n)Q_n(x)m_n=m(x) \end{aligned}$$

and $N_\varPhi (t){ admits}\mathbf {m}{ asdistribution}.{ Nowletussuppose}N_\varPhi (0){ admits}f{ asdistributionwith}f/m \in \ell ^2(\mathbf {m})$. By Proposition 5 we have that $v(t,x)=\mathbb {P}_f(N_\varPhi (t)=x)$ is a strong solution of (22) hence it holds

$$\begin{aligned} v(t,x)=m(x)\sum _{n \in E}\mathfrak {e}_\varPhi (t;\lambda _n)Q_n(x)f_n=m(x)f_0+m(x)\sum _{\begin{array}{c} n \in E \\ n\ge 1 \end{array}}\mathfrak {e}_\varPhi (t;\lambda _n)Q_n(x)f_n. \end{aligned}$$

Let us determine $f_0.{ Wehave},\,{ bydefinitionofscalarproductin}\ell ^2(\mathbf {m})$,

$$\begin{aligned} f_0=\sum _{x \in E}m(x)\frac{f(x)}{m(x)}Q_0=\sum _{x \in E}f(x)=1 \end{aligned}$$

hence we have

$$\begin{aligned} v(t,x)=m(x)+m(x)\sum _{\begin{array}{c} n \in E \\ n\ge 1 \end{array}}\mathfrak {e}_\varPhi (t;\lambda _n)Q_n(x)f_n. \end{aligned}$$

Now let us consider the summation part. First of all, let us recall that $\mathbf {m}{ isaprobabilitymeasure},\,{ hence}m(x)\le 1$. We have

$$\begin{aligned} \sum _{x \in E}m(x)f^2(x)\le \sum _{x \in E}m(x)\frac{f^2(x)}{m^2(x)}=\left\| f/m \right\| _{\ell ^2(\mathbf {m})} \end{aligned}$$

hence $f \in \ell ^2(\mathbf {m})$. By Lemma 3 we know that $\sum _{\begin{array}{c} n \in E \\ n\ge 1 \end{array}}\mathfrak {e}_\varPhi (t;\lambda _n)Q_n(x)f_n{ normallyconverges},\,{ hencewecantakethelimitas}t \rightarrow +\infty { underthesummationsign}.{ Nowletusobservethat}\lim _{t \rightarrow +\infty }E_\varPhi (t)=+\infty { almostsurely}.{ Ontheotherhand},\,{ wehave}e^{\lambda _n E_\varPhi (t)}\le 1$, thus we can use monotone convergence theorem to achieve

$$\begin{aligned} \lim _{t \rightarrow +\infty }\mathfrak {e}_\varPhi (t;\lambda _n)=\mathbb {E}[\lim _{t \rightarrow +\infty }e^{\lambda _n E_\varPhi (t)}]=0. \end{aligned}$$

Finally, by dominated convergence theorem, we have

$$\begin{aligned} \lim _{t \rightarrow +\infty }v(t,x)&=m(x)+m(x)\lim _{t \rightarrow +\infty }\sum _{\begin{array}{c} n \in E \\ n\ge 1 \end{array}}\mathfrak {e}_\varPhi (t;\lambda _n)Q_n(x)f_n\\&=m(x)+m(x)\sum _{\begin{array}{c} n \in E \\ n\ge 1 \end{array}}Q_n(x)f_n\lim _{t \rightarrow +\infty }\mathfrak {e}_\varPhi (t;\lambda _n)=m(x), \end{aligned}$$

concluding the proof. $\square$

6 Correlation Structure of Non-local Solvable Birth–Death Processes

Let us consider the potential measure $U_\varPhi (t)=\mathbb {E}[E_\varPhi (t)]{ ofthesubordinator}\sigma _\varPhi (t)$. As a consequence of Corollary 2, we can directly apply [9, Theorem $2]{ toobtainandexpressionofthecovarianceof}N_\varPhi (t){ intermsof}\mathfrak {e}_\varPhi (t;\lambda _1).{ Inparticular},\,{ withaneasyrefinementoftheproof},\,{ byusingthedistributionalderivativeofthepotentialfunction}U_\varPhi (t)$, we can get rid of some hypotheses of the aforementioned Theorem.

Proposition 6

Let $N(t)$ be a solvable birth–death process with state space $E,\,{ invariantmeasure}\mathbf {m}{ andfamilyofassociatedclassicalorthogonalpolynomials}(P_n)_{n \in E}.{ Letusdenoteby}\iota :E \rightarrow E{ theidentityfunctionandsupposethat}\iota =a_0+a_1Q_1.{ Thenitholds},\,{ forany}t \ge s\ge 0$

$$\begin{aligned}&Cov_m(N_\varPhi (t),N_\varPhi (s))\\&\quad =a_1^2\left( -\lambda _1\int _0^{s}\mathfrak {e}_\varPhi (t-\tau ;\lambda _1)dU_\varPhi (\tau )-2+2\mathfrak {e}_\varPhi (s;\lambda _1)+\mathfrak {e}_\varPhi (t;\lambda _1)\right) . \end{aligned}$$

As we expected, since we are composing a stationary process $N(t) ({ if}N(0){ admits}\mathbf {m}$ as distribution) with a non-stationary one $E_\varPhi (t),\,N_\varPhi (t)$ is not second-order stationary. To introduce the notion of memory for our process $N_\varPhi (t)$, we refer to the necessary conditions given in [12, 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 5

Given the function $\gamma (n)=Cov_m(N_\varPhi (n),N_\varPhi (0)){ for}n \in \mathbb {N}$:

• $N_\varPhi (t)$ is said to exhibit long-range dependence if $\gamma (n)\sim \ell (n)n^{-k}{} { where}\ell (n){ isaslowlyvaryingfunctionand}k \in \left( 0,1\right)$;

• $N_\varPhi (t)$ is said to exhibit short-range dependence if $\sum _{n=1}^{+\infty }|\gamma (n)|<+\infty$.

In the specific case of the initial state, we have the following Proposition.

Proposition 7

Let $N(t)$ be a solvable birth–death process with state space $E,\,{ invariantmeasure}\mathbf {m}{ andfamilyofassociatedclassicalorthogonalpolynomials}(P_n)_{n \in E}.{ Letusdenoteby}\iota :E \rightarrow E{ theidentityfunctionandsupposethat}\iota =a_0+a_1Q_1.{ Then},\,\lim _{t \rightarrow +\infty }Cov_m(N_\varPhi (t),N_\varPhi (0))=0.{ Moreover},\,{ if}\mathfrak {e}_\varPhi (t;\lambda _1) \sim Ct^{-\alpha }{} { as}t \rightarrow +\infty { forsome}\alpha \in (0,1),\,{ then}N_\varPhi (t)$ is long-range dependent.

Proof

One has

$$\begin{aligned} Cov_m(N_\varPhi (t),N_\varPhi (0))=a_1^2\mathfrak {e}_\varPhi (t;\lambda _1), \end{aligned}$$

where this identity can be achieved by direct calculations, observing that $Cov_m(N(t),N(0))=a_1^2e^{\lambda _1 t}.{ Thesecondassertioneasilyfollowsfromlastidentity}. \square$

Remark 2

Let us observe that, actually, since $Cov_m(N(t),N(0))=a_1^2e^{\lambda _1 t},\,{ wecanarguebydirectcalculationswithoutusingtheregularityhypotheseson}\mathbb {P}(E_\varPhi (t)\ge s)$.

Concerning the asymptotic behaviour of $\mathfrak {e}_\varPhi (t;\lambda _1),\,{ onecanobtainsomeinformationbythebehaviourof}\varPhi$. Indeed one has the following result.

Proposition 8

If $\varPhi { isregularlyvaryingat}0^+{ withorder}\alpha \in (0,1),\,{ then},\,{ foranyfixed}\lambda >0,\,\mathfrak {e}_\varPhi (t;-\lambda ){ isregularlyvaryingat}\infty { withorder}-\alpha$.

Proof

Let us consider $J(t)=\int _0^t \mathfrak {e}_\varPhi (s;-\lambda )\mathrm{d}s$ and let us observe that the Laplace-Stieltjes transform $\bar{J}(\eta ){ of}J(t)$ is given by

$$\begin{aligned} \bar{J}(\eta )=\frac{\varPhi (\eta )}{\eta (\varPhi (\eta )+\lambda )}. \end{aligned}$$

Since $\lambda { ispositive},\,{ onehasthat}\bar{J}(\eta ){ isregularlyvaryingat}0^+{ ifandonlyif}\frac{\varPhi (\eta )}{\eta }{} { is}({ since}\varPhi (0^+)=0{ bythefactthat}a=0$, see, for instance, [45]). In particular, since $\varPhi { isregularlyvaryingat}0^+{ withorder}\alpha ,\,{ then}\bar{J}(\eta ){ isregularlyvaryingat}0^+{ withorder}\alpha -1$.

By Karamata’s Tauberian theorem (see, for instance, [37] for a compact statement and [16] for the full proof), we know that $J(t){ isregularlyvaryingatinfinitywithorder}1-\alpha .{ Nowletusobservethat}J'(t)=\mathfrak {e}_\varPhi (t;-\lambda )$, that is monotone, hence, by Monotone density theorem (again, see [16, 37]), we have that $\mathfrak {e}_\varPhi (t;-\lambda ){ isregularlyvaryingat}\infty { oforder}-\alpha ,\,{ concludingtheproof}. \square$

Remark 3

Let us observe that the results coincides with [36, Theorem $2.1-(4)],\,{ exceptforthefactthatwearenotasking}\varPhi { tobespecialor}\nu$ to be absolutely continuous and we are not considering any hypothesis on the existence or regularity of the density of the subordinator.

In particular we get the following Corollary.

Corollary 4

Under the hypotheses of Proposition 7, if $\varPhi { isregularlyvaryingat}0^+{ withorder}\alpha \in (0,1),\,{ then}N_\varPhi (t)$ is long-range dependent.

Let us reconsider the examples given in Sect. 3 and study the asymptotic behaviour of the covariance.

• In the case $\varPhi (\lambda )=\lambda ^\alpha { wecanactuallyobtainanexplicitformulationoftheautocovariancefunctionfor}t \ge s >0$:

  $$\begin{aligned}&Cov_m(N_\varPhi (t),N_\varPhi (s))\\&\quad =a_1^2\left( E_\alpha (\lambda _1 t^\alpha )-\frac{\lambda _1 \alpha t^\alpha }{\varGamma (1+\alpha )}\int _0^{\frac{s}{t}}\frac{E_\alpha (\lambda _1t^\alpha (1-z)^\alpha )}{z^{1-\alpha }}dz\right) , \end{aligned}$$

  where the proof of such formula is identical to the one of [31, Theorem $3.1$]. Thus, we can use [31, Remark $3.3$] to obtain directly long-range dependence of the process.

• As $\varPhi (\lambda )=(\lambda +\theta )^\alpha -\theta ^\alpha ,\,{ wehavethat}\varPhi (\lambda ){ isactuallyregularlyvaryingat}0^+{ oforder}1.{ Indeedwehave}\varPhi (\lambda )/\lambda \rightarrow \alpha \theta ^{\alpha +1}{} { as}\lambda \rightarrow 0^+.{ Inparticularthismeansthatthefunction}\overline{J}(\eta )$ defined in Proposition 8 is such that $\overline{J}(0+)=\frac{\alpha \theta ^{\alpha +1}}{\lambda }{}$. By Karamata’s Tauberian theorem, we still have $\lim _{t \rightarrow +\infty }J(t)=\frac{\alpha \theta ^{\alpha +1}}{\lambda }.{ Thismeansinparticularthat}\mathfrak {e}_\varPhi (s;-\lambda ){ isintegrablein}(0,+\infty ){ andthus}\sum _{n=1}^{+\infty }\mathfrak {e}_\varPhi (n;-\lambda )<+\infty .{ Fromthisobservationweobtainthatinsuchcase}N_\varPhi (t)$ is short-range dependent.

• If we consider $\varPhi (\lambda )=\log (1+\lambda ^\alpha ){ for}\alpha \in (0,1),\,{ then}\varPhi (\lambda ){ isregularlyvaryingat}0^+{ oforder}\alpha ,\,{ since}\lim _{\lambda \rightarrow 0^+}\frac{\log (1+\lambda ^\alpha )}{\lambda ^\alpha }=1.{ Thus},\,{ inparticular},\,N_\varPhi$ is long-range dependent by Corollary 4.

• Finally, if $\varPhi (\lambda )=\log (1+\lambda ),\,{ wehave}\lim _{\lambda \rightarrow 0^+}\frac{\log (1+\lambda )}{\lambda }=1{ andthen}\overline{J}(0+)=\frac{1}{\lambda }.{ Thismeansagainthat}\mathfrak {e}_\varPhi (t;-\lambda ){ isintegrablein}(0,+\infty ){ andthen}N_\varPhi$ is short-range dependent.