Fractional Queues with Catastrophes and Their Transient Behaviour

Abstract

Starting from the definition of fractional M/M/1 queue given in the reference by Cahoy et al. in 2015 and M/M/1 queue with catastrophes given in the reference by Di Crescenzo et al. in 2003, we define and study a fractional M/M/1 queue with catastrophes. In particular, we focus our attention on the transient behaviour, in which the time-change plays a key role. We first specify the conditions for the global uniqueness of solutions of the corresponding linear fractional differential problem. Then, we provide an alternative expression for the transient distribution of the fractional M/M/1 model, the state probabilities for the fractional queue with catastrophes, the distributions of the busy period for fractional queues without and with catastrophes and, finally, the distribution of the time of the first occurrence of a catastrophe.

1. Introduction

Stochastic models for queueing systems have a wide range of applications in computer systems, sales points, telephone or telematic systems and also in several areas of science including biology, medicine and many others. The well known M/M/1 queueing model [1,2,3,4,5] constitutes the theoretical basis for building many other refined models for service systems.

Due to the Markov nature of its inter-arrival times of the customers and of its service times, the model can be mathematically treated in a simple manner, and, for this reason, it is widely used in many modeling contexts. Nevertheless, in the past few decades, the advent of fractional operators, such as fractional derivatives and integrals (see, for instance, [6] and [7] and references therein), has made it clear that different time scales, themselves random, that preserve memory (therefore not Markovian), allow the construction of more realistic stochastic models.

The introduction of the fractional Caputo derivative into the system of differential-difference equations for an M/M/1-type queue was done in [8], where, for a fractional M/M/1 queue, the state probabilities were determined. In this kind of queue model, the inter-arrival times and service times are characterized by Mittag–Leffler distributions [9]; in this case, the model does not have the property of memory loss that is typical of the exponential distributed times of the classical M/M/1 model. Indeed, a time-changed birth-death process [10,11], by means of an inverse stable subordinator [12], solves the corresponding fractional system of differential-difference equations and fractional Poisson processes [13] characterize the inter-arrival and service times.

The fractional M/M/1 model in [8] is an interesting and powerful model, not only because it is a generalization of the classical one, where the fractional order is set to 1, but also because its range of applications is extremely wide. Its importance can be further augmented by including in the model the occurrence of catastrophes, as it was considered in [14] for the classical M/M/1.

The catastrophe is a particular event that occurs in a random time leading to the instantaneous emptying of the system, or to a momentary inactivation of the system, as, for example, the action of a virus program that can make a computer system inactive [15]; other applications of models with catastrophes can be found in population dynamics and reliability contexts (see [16] and references therein).

Motivated by the mathematical need to enrich the fractional M/M/1 model of [8] with the inclusion of catastrophes, we study in this paper the above model; specifically, we determine the transient distribution, the distribution of the busy period (including that of the fractional M/M/1 queue of [8]) and the probability distribution of the time of the first occurrence of the catastrophe.

For these purposes, we need to guarantee the global uniqueness of the solution of the considered linear fractional Cauchy problem on Banach spaces. After recalling the definitions and known results in Section 2, we address the problem of uniqueness in Section 3. In Section 2, we also provide the transient distribution of the fractional M/M/1 model in an alternative form to that given in [8]. In Section 4, the distribution of the busy period for the fractional M/M/1 queue (without catastrophes) is obtained. Here, the time-changed birth-death process plays a key role to derive the results. In Section 5, we define the fractional queue with catastrophes; we are able to obtain the distribution of the transient state probabilities by following a strategy similar to that in [14]. We also found the distribution of the busy period and of the time of the first occurrence of the catastrophe starting from the empty system. Some special operators and functions used in this paper are specified in the Appendix A and Appendix B.

2. Definition of a Fractional Process Related to M/M/1 Queues

The classical M/M/1 queue process $N(t),t\ge 0$ can be described as continuous time Markov chain whose state space is the set $\{0,1,2,\dots \}$ and the state probabilities

$$p_n(t)=\mathbb{P}(N(t)=n|N(0)=0),n=0,1,2\dots$$

(1)

satisfy the following differential-difference equations:

$$\left\{\begin{array}{cc}{D}_t{p}_n(t)=-(\alpha +\beta )p_n(t)+\alpha p_{n-1}(t)+\beta p_{n+1}(t),\hfill & n\ge 1,\hfill \\ D_t{p}_0(t)=-\alpha p_0(t)+\beta p_1(t),\hfill & \\ p_n(0)={\delta }_{n,0},\hfill & n\ge 0,\hfill \end{array}\right.$$

(2)

where ${\delta }_{n,0}$ is the Kroeneker delta symbol, $D_t=\frac{d}{dt}$ and $\alpha ,\beta >0$ are the entrance and service rates, respectively.

Let $S_{\nu }(t),t\ge 0,\nu \in (0,1)$ be the Lévy $\nu$-stable subordinator with Laplace exponent given by:

$$log\mathbb{E}{e}^{-z{S}_{\nu }(t)}=-t{z}^{\nu },z>0.$$

Consider the inverse $\nu$-stable subordinator

$$L_{\nu }(t)=inf\{u\ge 0:S_{\nu }(u)>t\},t\ge 0.$$

For $0<\nu <1$, the fractional M/M/1 queue process $N^{\nu }(t),t\ge 0$ is defined by a non-Markovian time change $L_{\nu }(t)$ independent of $N(t),t\ge 0$, i.e.,

$$N^{\nu }(t)=N(L_{\nu }(t)),t\ge 0.$$

(3)

This process was defined in [8] and it is non-Markovian with non-stationary and non-independent increments. For $\nu =1$, by definition, $N^1(t)=N(t),t\ge 0$. Then, for a fixed $\nu \in (0,1]$, the state probabilities

$$p_n^{\nu }(t)=\mathbb{P}\{N^{\nu }(t)=n|N^{\nu }(0)=0\},n=0,1,\dots$$

(4)

of the number of customers in the system at time t in the fractional M/M/1 queue are characterized by arrivals and services determined by fractional Poisson processes of order $\nu \in (0,1]$ [13] with parameters $\alpha$ and $\beta$. They are solutions of the following system of differential-difference equations

$$\left\{\begin{array}{cc}{}_0^C{D}_t^{\nu }{p}_n^{\nu }(t)=-(\alpha +\beta )p_n^{\nu }(t)+\alpha p_{n-1}^{\nu }(t)+\beta p_{n+1}^{\nu }(t),\hfill & n\ge 1,\hfill \\ {}_0^C{D}_t^{\nu }{p}_0^{\nu }(t)=-\alpha p_0^{\nu }(t)+\beta p_1^{\nu }(t),\hfill & \\ p_n^{\nu }(0)={\delta }_{n,0},\hfill & n\ge 0,\hfill \end{array}\right.$$

(5)

where ${}_0^C{D}_t^{\nu }$ is the Caputo fractional derivative (see Appendix A).

Using Equation (5) and representation (3), the state probabilities are obtained in [8]:

$$\begin{array}{c}\hfill \begin{array}{cc}\hfill p_n^{\nu }(t)& =\left(1-\frac{\alpha }{\beta }\right){\left(\frac{\alpha }{\beta }\right)}^n\hfill \\ & +{\left(\frac{\alpha }{\beta }\right)}^n\sum _{k=0}^{+\infty }\sum _{m=0}^{n+k}\frac{k-m}{k+m}\left(\genfrac{}{}{0pt}{}{k+m}{k}\right){\alpha }^k\hfill \\ & \times {\beta }^{m-1}{t}^{\nu (k+m)-\nu }{E}_{\nu ,\nu (k+m)-\nu +1}^{k+m}(-(\alpha +\beta )t^{\nu }),\hfill \end{array}\end{array}$$

(6)

as well as its Laplace transform

$$\begin{array}{c}\hfill \begin{array}{cc}\hfill {\pi }_n^{\nu }(z)& ={\int }_0^{+\infty }{e}^{-zt}{p}_n^{\nu }(t)dt=\left(1-\frac{\alpha }{\beta }\right){\left(\frac{\alpha }{\beta }\right)}^n\frac{1}{z}+\hfill \\ & +{\left(\frac{\alpha }{\beta }\right)}^n\sum _{k=0}^{+\infty }\sum _{m=0}^{n+k}\frac{k-m}{k+m}\left(\genfrac{}{}{0pt}{}{k+m}{k}\right){\alpha }^k\hfill \\ & \times {\beta }^{m-1}\frac{z^{\nu -1}}{{(z^{\nu }+\alpha +\beta )}^{k+m}},z>0.\hfill \end{array}\end{array}$$

In Equation (6), the functions $E_{\nu ,\mu }^{\rho }$ are generalized Mittag–Leffler functions (see Appendix B). Note that $p_n^{\nu }(t)\ge 0\forall n\ge 0$ and ${\sum }_{n=0}^{+\infty }{p}_n^{\nu }(t)=1$.

Alternatively, let ${\displaystyle h_{\nu }(t,x)=\frac{d}{dx}\mathbb{P}\{L_{\nu }(t)\le x\},x\ge 0,}$ be the density of $L_{\nu }(t)$; then it is known (see, i.e., [17]) that

$${\int }_0^{+\infty }{e}^{-sx}{h}_{\nu }(t,x)dx=E_{\nu }(-s{t}^{\nu }),s\ge 0,$$

(7)

and (see, i.e., [18], Proposition 4.1)

$$h_{\nu }(t,x)=\frac{1}{\pi }{\int }_0^{+\infty }{u}^{\nu -1}{e}^{-tu-x{u}^{\nu }cos(\nu \pi )}sin\left(\pi \nu -x{u}^{\nu }sin(\pi \nu ))\right)du,x\ge 0.$$

(8)

Using (7) and an analytical expression for $p_n^{\nu }(t)$ given in [19], we can write down an alternative expression for (6) as

$$\begin{array}{c}\hfill \begin{array}{cc}\hfill p_n^{\nu }(t)& ={\left(\frac{p}{q}\right)}^n\sum _{r=n}^{\infty }\frac{{(\alpha +\beta )}^r}{r!}{t}^{r\nu }{E}_{\nu }^{(r)}(-(\alpha +\beta )t^{\nu })\hfill \\ & \times \sum _{r=0}^{\left[\frac{r-n}{2}\right]}\frac{r+1-2k}{r+1}\left(\genfrac{}{}{0pt}{}{r+1}{k}\right)p^k{q}^{r-k},\hfill \end{array}\end{array}$$

(9)

where ${\displaystyle p=\frac{\alpha }{\alpha +\beta }}$, ${\displaystyle q=\frac{\beta }{\alpha +\beta }}$, and $E_{\nu }^{(r)}(-(\alpha +\beta )t^{\nu })$ is the $r$-th derivative of the function $E_{\nu }(z)$ evaluated at $z=-(\alpha +\beta )t^{\nu }$.

Actually, it is easy to see from (7) that

$${\int }_0^{+\infty }{e}^{-sx}{x}^r{h}_{\nu }(t,x)dx=E_{\nu }^{(r)}(-s{t}^{\nu })t^{\nu r};$$

thus, using [19] and (3), we have

$$\begin{array}{cc}\hfill p_n^{\nu }(t)& ={\int }_0^{+\infty }{p}_n(s)h_{\nu }(t,s)ds\hfill \\ & ={\left(\frac{p}{q}\right)}^n\sum _{r=n}^{\infty }\frac{{(\alpha +\beta )}^r}{r!}{\int }_0^{+\infty }{e}^{-(\alpha +\beta )s}{s}^{\nu }{h}_{\nu }(t,s)ds\hfill \\ & \times \sum _{r=0}^{\left[\frac{r-n}{2}\right]}\frac{r+1-2k}{r+1}\left(\genfrac{}{}{0pt}{}{r+1}{k}\right)p^k{q}^{r-k},\hfill \end{array}$$

and formula (9) follows. On the other hand, using (8), we have

$$\begin{array}{cc}\hfill p_n^{\nu }(t)& ={\int }_0^{+\infty }{p}_n(s)h_{\nu }(t,s)ds\hfill \\ & ={\left(\frac{p}{q}\right)}^n\sum _{r=n}^{\infty }\frac{{(\alpha +\beta )}^r}{r!}{\int }_0^{+\infty }{u}^{\nu -1}{e}^{-tu}{F}_{\nu ,r}(u)du\hfill \\ & \times \sum _{r=0}^{\left[\frac{r-n}{2}\right]}\frac{r+1-2k}{r+1}\left(\genfrac{}{}{0pt}{}{r+1}{k}\right)p^k{q}^{r-k},\hfill \end{array}$$

where

$$F_{\nu ,r}(u)=\frac{1}{\pi }{\int }_0^{+\infty }exp\left\{-(\alpha +\beta )x-x{u}^{\nu }cos(\nu \pi )\right\}{x}^{r}sin\left(\pi \nu -x{u}^{\nu }sin(\pi \nu )\right)dx.$$

3. Linear Fractional Cauchy Problems on Banach Spaces

In order to describe the transient probabilities for our queues, we will need some uniqueness results for solutions of linear fractional Cauchy problems defined on Banach spaces. To do that, let us recall the following Theorem (Theorem $3.19$ from [20]):

Theorem 1.

Let $(X,|·|)$ be a Banach space and $J=[0,T]$ for some $T>0$. Consider the ball $B_R=\{x\in X:|x|\le R\}$. Let $\nu \in (0,1)$ and $f:J\times B_R\to X$ and consider the following Cauchy problem:

$$\left\{\begin{array}{c}{}_0^C{D}_t^{\nu }x(t)=f(t,x(t)),\hfill \\ x(0)=x_0,\hfill \end{array}\right.$$

(10)

where ${}_0^C{D}_t^{\nu }$ is the Caputo derivative operator (see Appendix A).

Suppose that:

• $f\in C(J\times B_R,X)$;
• There exists a constant $\overline{M}(R)>0$ such that

  $$|f(t,x(t))|\le \overline{M}(R)$$

  for all $x\in B_R$ and $t\in J$ and such that

  $$R\ge |x_0|+\frac{\overline{M}(R)T^{\nu }}{\mathrm{\Gamma }(\nu +1)};$$
• There exists a constant $\overline{L}>0$ such that $\overline{L}\ge \frac{2\overline{M}(R)}{\mathrm{\Gamma }(\nu +1)}$;
• There exists a constant ${\overline{L}}_0>0$ such that

  $$|f(t,x_1)-f(t,x_2)|\le {\overline{L}}_0|x_1-x_2|$$

  for all $x_1,x_2\in B_R$ and $t\in J$;
• There exist constants ${\nu }_1\in (0,\nu )$ and $\tau >0$ such that

  $${\overline{L}}_A=\frac{{\overline{L}}_0}{\mathrm{\Gamma }(\nu )}\frac{T^{(1+\beta )(1-{\nu }_1)}}{{(1+\beta )}^{1-{\nu }_1}}{\left(\frac{{\nu }_1}{\tau }\right)}^{{\nu }_1}<1,$$

  where $\beta =\frac{\nu -1}{1-{\nu }_1}$.

Then, if $x_0\in B_R$, the problem (10) admits a unique solution $x^{*}\in C_{\nu }(J,B_R)$.

The previous theorem can be easily adapted to the case in which $J=[t_0,T+t_0]$ and the starting point of the derivative is $t_0$. Since we are interested in linear (eventually non-homogeneous) equations, let us show how the previous theorem can be adapted in such a case.

Corollary 1.

Consider the system (10) and suppose $f(t,x)=Ax+\xi$ where $A:X\to X$ is a linear and continuous operator and $\xi \in X$. Then, there exists a $R>|x_0|$ and $T>0$ such that the system admits a unique solution $x^{*}\in C_{\nu }(J,B_R)$.

Proof. 

Observe that, if $|x|\le R$, then

$$|f(x)|\le \parallel A\parallel |x|+|\xi |\le \parallel A\parallel R+|\xi |.$$

Let us choose T such that the conditions of Theorem 1 are verified. To do that, consider $\overline{M}(R)=\parallel A\parallel R+|\xi |$. Fix $R\ge |x_0|$ and define $\tilde{R}=R+\epsilon$ for some $\epsilon >0$. Define then

$$T={\left[\frac{\epsilon \mathrm{\Gamma }(\nu +1)}{\overline{M}(\tilde{R})}\right]}^{\frac{1}{\nu }}$$

and observe that

$$|x_0|+\frac{\overline{M}(\tilde{R})T^{\nu }}{\mathrm{\Gamma }(\nu +1)}=|x_0|+\epsilon \le R+\epsilon =\tilde{R}.$$

Thus, one can fix $\overline{L}=\frac{2\overline{M}(\tilde{R})}{\mathrm{\Gamma }(\nu +1)}$ and ${\overline{L}}_0=\parallel A\parallel$. Moreover, since for fixed ${\nu }_1\in (0,\nu )$ the function $\tau \mapsto {\overline{L}}_A(\tau )$ is decreasing and ${lim}_{\tau \to 0}{\overline{L}}_A(\tau )=0$, then one can easily find a $\tau >0$ such that ${\overline{L}}_A(\tau )<1$. Since we are under the hypotheses of Theorem 1, then we have shown the local existence and uniqueness of a solution $x^{*}\in C_{\nu }(J,B_{\tilde{R}})$. ☐

However, using such corollary, we can only afford local uniqueness. Global uniqueness of the solution of the Cauchy problem (10) can be obtained with the additional hypothesis that such solution is uniformly bounded:

Corollary 2.

Suppose we are under the hypotheses of Corollary 1. If there exists a solution $x^{*}\in C([0,+\infty [,X)$ and a constant $k>0$ such that for any $t\ge 0$ we have $|x^{*}(t)|\le k$, then such solution is unique.

Proof. 

Observe that $|x_0|\le k$ and then fix $\tilde{R}=k+\epsilon$. Define

$$\Delta T={\left[\frac{\epsilon \mathrm{\Gamma }(\nu +1)}{\overline{M}(\tilde{R})}\right]}^{\frac{1}{\nu }}.$$

Fix $T_1=\Delta T$ and observe that, by using Corollary 1, there exists a unique solution in $[0,T_1]$. Since $x^{*}$ is a solution of such problem, we have that $x^{*}$ is unique. Suppose we have defined $T_{n-1}$ such that $x^{*}$ is the unique solution of system (10) in $[0,T_{n-1}]$. Consider the problem

$$\left\{\begin{array}{c}{}_{T_{n-1}}^C{D}_t^{\nu }x(t)=f(x(t)),\hfill \\ x(T_{n-1})=x^{*}(T_{n-1}).\hfill \end{array}\right.$$

(11)

Define then $T_n=T_{n-1}+\Delta T$ and observe that, since $|x^{*}(T_{n-1})|\le k$, by using Corollary 1, there exists a unique solution in $[T_{n-1},T_n]$.

By using a change of variables, it is easy to show that

$${}_{T_{n-1}}^C{D}_t^{\nu }x={}_0^C{D}_{t-T_{n-1}}^{\nu }\tilde{x},$$

where $\tilde{x}:t\mapsto x(t+T_{n-1})$. By using such relation, we have that system (11) is equivalent to

$$\left\{\begin{array}{c}{}_0^C{D}_t^{\nu }\tilde{x}(t)=f(\tilde{x}(t)),\hfill \\ \tilde{x}(0)=x^{*}(T_{n-1}),\hfill \end{array}\right.$$

whose unique solution is $\tilde{x}(t)=x^{*}(t+T_{n-1})$ so that $x(t)=x^{*}(t)$ and $x^{*}(t)$ is the unique solution of system (10) in $[0,T_n]$. Since $T_n\to +\infty$ as $n\to +\infty$, we have global uniqueness of limited solutions. ☐

4. The Fractional M/M/1 Queue

Let us consider again the fractional M/M/1 process $N^{\nu }(t),t\ge 0$ defined by (3) with state probabilites in (6).

Consider the Hilbert space $(l^2(\mathbb{R}),|·{|}_2)$ with the norm ${|x|}_2^2={\sum }_{k=0}^{+\infty }{x}_k^2$ and let $C_{\nu }([0,T],l^2(\mathbb{R}))$ be the space of the $\nu$-Hölder continuous functions from $[0,T]$ to $l^2(\mathbb{R})$. One can rewrite the system (5) in $l^2(\mathbb{R})$ as follows:

$$\left\{\begin{array}{c}{}_0^C{D}_t^{\nu }{\mathit{p}}^{\nu }(t)=A_0{\mathit{p}}^{\nu }(t),\hfill \\ {\mathit{p}}^{\nu }(0)={({\delta }_{n,0})}_{n\ge 0},\hfill \end{array}\right.$$

(12)

where ${\mathit{p}}^{\nu }(t)={(p_n^{\nu }(t))}_{n\ge 0}\in C([0,T],l^2(\mathbb{R}))$ and

$$A_0=\left(\begin{array}{cccccc}-\alpha & \beta & 0& 0& 0& \cdots \\ \alpha & -(\alpha +\beta )& \beta & 0& 0& \cdots \\ 0& \alpha & -(\alpha +\beta )& \beta & 0& \cdots \\ 0& 0& \alpha & -(\alpha +\beta )& \beta & \cdots \\ ⋮& ⋮& ⋮& \ddots & \ddots & \ddots \end{array}\right)$$

is an infinite tridiagonal matrix with $A_0={(a_{i,j})}_{i,j\ge 0}$. Let us show the following:

Lemma 1.

The linear operator $A_0$ is continuous and $\parallel A_0\parallel \le 2(\alpha +\beta )$.

Proof. 

To show that $A_0$ is continuous, let us use Schur’s test (Theorem $5.2$ in [21]). Observe that

$$\begin{array}{cc}\sum _{k=0}^{+\infty }|a_{k,0}|=2\alpha ,\hfill & \hfill \sum _{k=0}^{+\infty }|a_{k,j}|=2(\alpha +\beta )\text{ for }j\ne 0\end{array}$$

so that, in general,

$$\sum _{k=0}^{+\infty }|a_{k,j}|\le 2(\alpha +\beta ).$$

Moreover,

$$\begin{array}{cc}\sum _{k=0}^{+\infty }|a_{0,k}|=\alpha +\beta ,\hfill & \hfill \sum _{k=0}^{+\infty }|a_{j,k}|=2(\alpha +\beta )\text{ for }j\ne 0,\end{array}$$

so that, in general,

$$\sum _{k=0}^{+\infty }|a_{j,k}|\le 2(\alpha +\beta ).$$

By Schur’s test, we have that $A_0$ is a bounded operator on $l^2$ and

$$\parallel A_0\parallel \le 2(\alpha +\beta ).$$

 ☐

Thus, by Corollary 1, we obtain local existence and uniqueness of the solution of system (5). Global uniqueness can be obtained a posteriori, since the solutions of such system are known.

Let us also observe that the distributions of the inter-arrival times are Mittag–Leffler distributions. To do that, consider the system, for fixed $n\ge 0$

$$\left\{\begin{array}{c}{}_0^C{D}_t^{\nu }{b}_n^{\nu }(t)=-\alpha b_n^{\nu }(t),\hfill \\ {}_0^C{D}_t^{\nu }{b}_{n+1}^{\nu }(t)=\alpha b_n^{\nu }(t),\hfill \\ b_n^{\nu }(0)=1,\hfill \\ b_{n+1}^{\nu }(0)=0,\hfill \end{array}\right.$$

which are the state probabilities of a queue with null death rate, fixed birth rate, starting with n customers and with an absorbent state $n+1$. Under such assumptions, $b_{n+1}^{\nu }(t)$ is the probability that a customer arrives before t. Moreover, the normalizing condition becomes

$$b_n^{\nu }(t)+b_{n+1}^{\nu }(t)=1.$$

One can solve the first equation (see Appendix A) to obtain

$$b_n^{\nu }(t)=E_{\nu }(-\alpha t^{\nu }),$$

where $E_{\nu }$ is the one-parameter Mittag–Leffler function (see Appendix B), and then, by using the normalizing condition, we have

$$b_{n+1}^{\nu }(t)=1-E_{\nu }(-\alpha t^{\nu }).$$

In a similar way, let us show that the distributions of the service times are Mittag–Leffler distributions. To show that, consider the system, for fixed $n\ge 0$,

$$\left\{\begin{array}{c}{}_0^C{D}_t^{\nu }{d}_n^{\nu }(t)=\beta d_{n+1}^{\nu }(t),\hfill \\ {}_0^C{D}_t^{\nu }{d}_{n+1}^{\nu }(t)=-\beta d_{n+1}^{\nu }(t),\hfill \\ d_n^{\nu }(t)=0,\hfill \\ d_{n+1}^{\nu }(t)=1,\hfill \end{array}\right.$$

which are the state probabilities of a queue with null birth rate, fixed death rate, starting with $n+1$ customers with an absorbent state n. Under such assumption, $d_n^{\nu }(t)$ is the probability that a customer is served before t. Moreover, the normalizing condition becomes

$$d_n^{\nu }(t)+d_{n+1}^{\nu }(t)=1.$$

One can solve the second equation to obtain

$$d_{n+1}^{\nu }(t)=E_{\nu }(-\beta t^{\nu }),t\ge 0$$

and then, by using the normalizing condition, we have

$$d_n^{\nu }(t)=1-E_{\nu }(-\beta t^{\nu }),t\ge 0.$$

Moreover, since we know that $\forall t\ge 0p_n^{\nu }(t)\ge 0$ and ${\sum }_{n=0}^{\infty }{p}_n^{\nu }(t)=1$, by the continuous inclusion $l^1(\mathbb{R})\subseteq l^2(\mathbb{R})$ (see [22]), ${(p_n(t))}_{n\ge 0}$ is uniformly bounded in $l^2(\mathbb{R})$ and then, by Corollary 2, it is the (global) unique solution of system (5).

Distribution of the Busy Period

We want to determine the probability distribution $K^{\nu }(t)$ of the busy period $K^{\nu }$ of a fractional M/M/1 queue. To do this, we will follow the lines of the proof given in [1] and [4].

Theorem 2.

Let $K^{\nu }$ be the random variable describing the duration of the busy period of a fractional M/M/1 queue $N^{\nu }(t)$ and consider $K^{\nu }(t)=\mathbb{P}(K^{\nu }\le t)$. Then,

$$K^{\nu }(t)=1-\sum _{n=1}^{+\infty }\sum _{m=0}^{+\infty }{C}_{n,m}{t}^{\nu (n+2m-1)}{E}_{\nu ,\nu (n+2m-1)+1}^{n+2m}(-(\alpha +\beta )t^{\nu }),$$

(13)

where

$$C_{n,m}=\left(\genfrac{}{}{0pt}{}{n+2m}{m}\right)\frac{n}{n+2m}{\alpha }^{n+m-1}{\beta }^m.$$

(14)

Proof. 

Let us first define a queue ${\overline{N}}^{\nu }(t)$ such that $\mathbb{P}({\overline{N}}^{\nu }(0)=1)=1$ and ${\overline{N}}^{\nu }(t)$ behaves like $N^{\nu }(t)$ except for the state 0 being an absorbent state. Thus, state probability functions are solution of the following system

$$\left\{\begin{array}{c}{}_0^C{D}_t^{\nu }{\overline{p}}_0^{\nu }=\beta {\overline{p}}_1^{\nu }(t),\hfill \\ {}_0^C{D}_t^{\nu }{\overline{p}}_1^{\nu }=-(\alpha +\beta ){\overline{p}}_1^{\nu }(t)+\beta {\overline{p}}_2^{\nu }(t),\hfill \\ {}_0^C{D}_t^{\nu }{\overline{p}}_n^{\nu }=-(\alpha +\beta ){\overline{p}}_n^{\nu }(t)+\alpha {\overline{p}}_{n-1}^{\nu }(t)+\beta {\overline{p}}_{n+1}^{\nu }(t),\hfill & n\ge 2,\hfill \\ {\overline{p}}_n^{\nu }(0)={\delta }_{n,1},\hfill & n\ge 0.\hfill \end{array}\right.$$

(15)

First, we want to show that, if we consider $L_{\nu }(t),$ the inverse of a $\nu$-stable subordinator that is independent from ${\overline{N}}^1(t)$, then ${\overline{N}}^{\nu }(t)\stackrel{d}{=}{\overline{N}}^1(L_{\nu }(t))$. To do that, consider the probability generating function $G^{\nu }(z,t)$ of ${\overline{N}}^{\nu }(t)$ defined as

$$G^{\nu }(z,t)=\sum _{k=0}^{+\infty }{z}^k{\overline{p}}_k^{\nu }(t).$$

(16)

From system (15), we know that $G^{\nu }(z,t)$ solves the following fractional Cauchy problem:

$$\left\{\begin{array}{c}z{}_0^C{D}_t^{\nu }{G}^{\nu }(z,t)=[\alpha z^2-(\alpha +\beta )z+\beta ][G^{\nu }(z,t)-{\overline{p}}_0^{\nu }(t)],\hfill \\ G^{\nu }(z,0)=z,\hfill \end{array}\right.$$

(17)

which, for $\nu =1$, becomes

$$\left\{\begin{array}{c}z\frac{d}{dt}{G}^1(z,t)=[\alpha z^2-(\alpha +\beta )z+\beta ][G^1(z,t)-{\overline{p}}_0^1(t)],\hfill \\ G^1(z,0)=z.\hfill \end{array}\right.$$

(18)

Taking the Laplace transform in Equation (17) and using Equation (A1), we have

$$z[s^{\nu }{\tilde{G}}^{\nu }(z,s)-z{s}^{\nu -1}]=[\alpha z^2-(\alpha +\beta )z+\beta ][{\tilde{G}}^{\nu }(z,s)-{\overline{\pi }}_0^{\nu }(s)],$$

(19)

where ${\tilde{G}}^{\nu }(z,s)$ and ${\overline{\pi }}_0^{\nu }(s)$ are Laplace transforms of $G^{\nu }(z,t)$ and ${\overline{p}}_0^{\nu }(t)$.

We know that ${\overline{N}}^{\nu }(t)\stackrel{d}{=}{\overline{N}}^1(L_{\nu }(t))$ if and only if

$${\overline{p}}_n^{\nu }(t)=\mathbb{P}({\overline{N}}^{\nu }(t)=n)=\mathbb{P}({\overline{N}}^1(L_{\nu }(t))=n)={\int }_0^{+\infty }{\overline{p}}_n^1(y)\mathbb{P}(L_{\nu }(t)\in dy)$$

(20)

and then if and only if, by Equation (16),

$$G^{\nu }(z,t)={\int }_0^{+\infty }{G}^1(z,y)\mathbb{P}(L_{\nu }(t)\in dy).$$

(21)

Taking the Laplace transform in Equations (20) and (21) for $n=0$ and by using (see, i.e., Equation (10) in [12])

$$ℒ[\mathbb{P}(L_{\nu }(t)\in dy)](s)=s^{\nu -1}{e}^{-ys}dy,$$

(22)

we know we have to show that

$${\overline{\pi }}_0^{\nu }(s)={\int }_0^{+\infty }{\overline{p}}_n^1(y)s^{\nu -1}{e}^{-y{s}^{\nu }}dy$$

(23)

and

$${\tilde{G}}^{\nu }(z,s)={\int }_0^{+\infty }{G}^1(z,y)s^{\nu -1}{e}^{-y{s}^{\nu }}dy.$$

(24)

Since Equation (17) admits a unique solution, then we only need to show that the right-hand sides of Equations (23) and (24) solve Equation (19), that is to say that we have to verify

$$\begin{array}{c}z\left[s^{\nu }{\int }_0^{+\infty }{G}^1(z,y)e^{-y{s}^{\nu }}dy-z\right]\hfill \\ \hfill =[\alpha z^2-(\alpha +\beta )z+\beta ]\left[{\int }_0^{+\infty }{G}^1(z,y)e^{-y{s}^{\nu }}dy-{\int }_0^{+\infty }{\overline{p}}_0^1(y)e^{-y{s}^{\nu }}dy\right].\end{array}$$

(25)

To do that, consider the right-hand side of the previous equation and, recalling that $G^1(z,t)$ is solution of Equation (18):

$${\int }_0^{+\infty }[\alpha z^2-(\alpha +\beta )z+\beta ][G^1(z,y)-{\overline{p}}_0^1(y)]e^{-y{s}^{\nu }}dy={\int }_0^{+\infty }\left(\frac{d}{dy}{G}^1(z,y)\right)e^{-y{s}^{\nu }}dy$$

and then, by integrating by parts, we have Equation (25).

Now remark that ${\overline{p}}_0^{\nu }(t)=B^{\nu }(t)$. Thus, we want to determine ${\overline{p}}_0^{\nu }(t)$. To do that, let us recall, from [1,4] that

$${\overline{p}}_n^1(t)=n{t}^{-1}{\alpha }^{\frac{n}{2}-1}{\beta }^{-\frac{n}{2}}{e}^{-(\alpha +\beta )t}{I}_n(2\sqrt{\alpha \beta }t)\text{ for }n\ge 1$$

from which, explicitly writing $I_n(2\sqrt{\alpha \beta }t)$, we have

$${\overline{p}}_n^1(t)=\sum _{m=0}^{+\infty }\left(\genfrac{}{}{0pt}{}{n+2m}{m}\right)\frac{n}{n+2m}\frac{1}{(n+2m-1)!}{\alpha }^{n+m-1}{\beta }^m{t}^{n+2m-1}{e}^{-(\alpha +\beta )t}\text{ for }n\ge 1.$$

Posing $C_{n,m}=\left(\genfrac{}{}{0pt}{}{n+2m}{m}\right)\frac{n}{n+2m}{\alpha }^{n+m-1}{\beta }^m$, we have

$${\overline{p}}_n^1(t)=\sum _{m=0}^{+\infty }\frac{C_{n,m}}{(n+2m-1)!}{t}^{n+2m-1}{e}^{-(\alpha +\beta )t}\text{ for }n\ge 1$$

and then

$${\overline{p}}_0^1(t)=1-\sum _{n=1}^{+\infty }\sum _{m=0}^{+\infty }\frac{C_{n,m}}{(n+2m-1)!}{t}^{n+2m-1}{e}^{-(\alpha +\beta )t}.$$

(26)

Since ${\overline{N}}^{\nu }(t)={\overline{N}}^1(L_{\nu }(t))$, we have

$${\overline{p}}_0^{\nu }(t)={\int }_0^{+\infty }{\overline{p}}_0^1(y)\mathbb{P}(L_{\nu }(t)\in dy)$$

and then, using Equation (26), we have

$${\overline{p}}_0^{\nu }(t)=1-\sum _{n=1}^{+\infty }\sum _{m=0}^{+\infty }\frac{C_{n,m}}{(n+2m-1)!}{\int }_0^{+\infty }{y}^{n+2m-1}{e}^{-(\alpha +\beta )y}\mathbb{P}(L_{\nu }(t)\in dy).$$

(27)

Taking the Laplace transform in Equation (27), using Equation (22), we have

$${\overline{\pi }}_0^{\nu }(s)=\frac{1}{s}-\sum _{n=1}^{+\infty }\sum _{m=0}^{+\infty }\frac{C_{n,m}}{(n+2m-1)!}{s}^{\nu -1}{\int }_0^{+\infty }{y}^{n+2m-1}{e}^{-(\alpha +\beta +s^{\nu })y}dy$$

and then integrating

$${\overline{\pi }}_0^{\nu }(s)=\frac{1}{s}-\sum _{n=1}^{+\infty }\sum _{m=0}^{+\infty }{C}_{n,m}\frac{s^{\nu -1}}{{(\alpha +\beta +s^{\nu })}^{n+2m}}.$$

Finally, using formula (A2), we have

$${\overline{p}}_0^{\nu }(s)=1-\sum _{n=1}^{+\infty }\sum _{m=0}^{+\infty }{C}_{n,m}{t}^{\nu (n+2m-1)}{E}_{\nu ,\nu (n+2m-1)+1}^{n+2m}(-(\alpha +\beta )t^{\nu })$$

 ☐

Remark 1.

As $\nu \to 1$ we obtain, by using

$$E_{1,n+2m}^{n+2m}(-(\alpha +\beta )t)=\frac{e^{-(\alpha +\beta )t}}{(n+2m-1)!},$$

that ${\overline{p}}_0^{\nu }(t)\to {\overline{p}}_0^1(t)$ and then $K^{\nu }(t)\to K^1(t)$.

5. The Fractional M/M/1 Queue with Catastrophes

Let us consider a classical M/M/1 queue with FIFO discipline and subject to catastrophes whose effect is to instantaneously empty the queue [14] and let $N_{\xi }^1(t)$ be the number of customers in the system at time t with state probabilities

$$p_n^{1,\xi }(t)=\mathbb{P}(N_{\xi }^1(t)=n|N_{\xi }^1(0)=0),n=0,1,\dots$$

Then, the function $p_n^{1,\xi }$ satisfy the following differential-difference equations:

$$\left\{\begin{array}{cc}{D}_t{p}_0^{1,\xi }(t)=-(\alpha +\xi )p_0^{1,\xi }(t)+\beta p_1^{1,\xi }(t)+\xi ,\hfill & \\ D_t{p}_n^{1,\xi }(t)=-(\alpha +\beta +\xi )p_n^{1,\xi }(t)+\alpha p_{n-1}^{1,\xi }(t)+\beta p_{n+1}^{1,\xi }(t),\hfill & n\ge 1,\hfill \\ p_n^{1,\xi }(0)={\delta }_{n,0},\hfill & n\ge 0,\hfill \end{array}\right.$$

(28)

where ${\delta }_{n,0}$ is the Kroeneker delta symbol, $D_t=\frac{d}{dt}$, $\alpha ,\beta >0$ are the entrance and service rates, respectively, and $\xi >0$ is the rate of the catastrophes when the system is not empty.

For $\nu \in (0,1)$, we define the fractional M/M/1 queue process with catastrophes as

$$N_{\xi }^{\nu }(t)=N_{\xi }^1(L_{\nu }(t)),t\ge 0,$$

where $L_{\nu }$ is an inverse $\nu$-stable subordinator that is independent of $N_{\xi }^1(t),t\ge 0$ (see Section 2).

We will show that the state probabilities

$$p_n^{\nu ,\xi }:=\mathbb{P}(N_{\xi }^{\nu }(t)=n|N_{\xi }^{\nu }(0)=0)$$

satisfy the following differential-difference fractional equations:

$$\left\{\begin{array}{cc}{}_0^C{D}_t^{\nu }{p}_0^{\nu ,\xi }(t)=-(\alpha +\xi )p_0^{\nu ,\xi }(t)+\beta p_1^{\nu ,\xi }(t)+\xi ,\hfill & \\ {}_0^C{D}_t^{\nu }{p}_n^{\nu ,\xi }(t)=-(\alpha +\beta +\xi )p_n^{\nu ,\xi }(t)+\alpha p_{n-1}^{\nu ,\xi }(t)+\beta p_{n+1}^{\nu ,\xi }(t),\hfill & n\ge 1,\hfill \\ p_n^{\nu ,\xi }(0)={\delta }_{n,0},\hfill & n\ge 0,\hfill \end{array}\right.$$

(29)

where ${}_0^C{D}_t^{\nu }$ is the Caputo fractional derivative (see Appendix A).

In the classical case, catastrophes occur according to a Poisson process with rate $\xi$ if the system is not empty. In our case, consider for a fixed $n>0$,

$$\left\{\begin{array}{c}{}_0^C{D}_t^{\nu }{c}_0^{\nu }(t)=\xi (1-c_0^{\nu }(t)),\hfill \\ {}_0^C{D}_t^{\nu }{c}_n^{\nu }(t)=-\xi c_n^{\nu }(t),\hfill \\ c_0^{\nu }(0)=0,\hfill \\ c_n^{\nu }(0)=1,\hfill \end{array}\right.$$

which describes the state probabilities of an initially non empty system with null birth and death rate but positive catastrophe rate. In such case, $c_0^{\nu }(t)$ is the probability a catastrophe occurs before time t. Moreover, the normalization property becomes

$$c_0^{\nu }(t)+c_n^{\nu }(t)=1.$$

In such case, we can solve the second equation to obtain

$$c_n^{\nu }(t)=E_{\nu }(-\xi t^{\nu }),t\ge 0.$$

Using the normalization property, we finally obtain

$$c_0^{\nu }(t)=1-E_{\nu }(-\xi t^{\nu }),t\ge 0$$

(30)

and then the distributions of the inter-occurrence of the catastrophes are Mittag–Leffler distributions. We can conclude that, in the fractional case, catastrophes occur according to a fractional Poisson process ([10,11,13]) with rate $\xi$ if the system is not empty. Since the operators ${}_0^C{D}_t^{\nu }$ are Caputo fractional derivatives, we expect the stationary behaviour of the queue to be the same as the classic one. Denoting with $N_{\xi }^1$ the number of customers in the system at the steady state of a classical M/M/1 with catastrophes and defining the state probabilities

$$q_n=\mathbb{P}(N_{\xi }^1=n),n\ge 0,$$

we can use the results obtained in [15] to observe that

$$q_n=\left(1-\frac{1}{z_1}\right){\left(\frac{1}{z_1}\right)}^n,n\ge 0,$$

(31)

where $z_1$ is the solution of

$$\alpha z^2-(\alpha +\beta +\xi )z+\beta =0$$

(32)

such that $z_1>1$. Let us call $z_2$ the other solution of Equation (32) and observe that $0<z_2<1<z_1$. Some properties coming from such equations that will be useful hereafter are

$$\alpha +\beta +\xi =\alpha z_i+\frac{\beta }{z_i}$$

(33)

and

$$\alpha z_i^2=(\alpha +\beta +\xi )z_i-\beta$$

(34)

with $i=1,2$.

5.1. Alternative Representation of the Fractional M/M/1 Queue with Catastrophes

We want to obtain an alternative representation of the fractional M/M/1 queue with catastrophes in a way which is similar to Lemma $2.1$ in [14]. To do that, we firstly need to assure that system (29) admits a unique uniformly bounded solution. To do that, let us write system (29) in the form

$$\left\{\begin{array}{c}{}_0^C{D}_t^{\nu }{\mathit{p}}^{\nu ,\xi }(t)=f({\mathit{p}}^{\nu ,\xi }(t)),\hfill \\ {\mathit{p}}^{\nu ,\xi }(t)={({\delta }_{n,0})}_{n\ge 0},\hfill \end{array}\right.$$

(35)

where ${\mathit{p}}^{\nu ,\xi }(t)={(p_n^{\nu ,\xi }(t))}_{n\ge 0}\in C([0,T],l^2(\mathbb{R}))$, $f(\mathit{x})=A_{\xi }\mathit{x}+\mathit{\xi }$, $\mathit{\xi }=(\xi ,0,\dots ,0,\dots )$ and

$$A_{\xi }=\left(\begin{array}{cccccc}-(\alpha +\xi )& \beta & 0& 0& 0& \cdots \\ \alpha & -(\alpha +\beta +\xi )& \beta & 0& 0& \cdots \\ 0& \alpha & -(\alpha +\beta +\xi )& \beta & 0& \cdots \\ 0& 0& \alpha & -(\alpha +\beta +\xi )& \beta & \cdots \\ ⋮& ⋮& ⋮& \ddots & \ddots & \ddots \end{array}\right)$$

is an infinite tridiagonal matrix with $A_{\xi }={(a_{i,j})}_{i,j\ge 0}$. We need to show the following:

Lemma 2.

The linear operator $A_{\xi }$ is continuous and $\parallel A_{\xi }\parallel \le 2(\alpha +\beta )+\xi$.

Proof. 

To obtain an estimate of the norm of $A_{\xi }$, let us use Schur’s test. Observe that

$$\begin{array}{cc}\hfill \sum _{k=0}^{+\infty }|a_{k,0}|=2\alpha +\xi ,& \hfill \sum _{k=0}^{+\infty }|a_{k,j}|=2\alpha +2\beta +\xi \text{ with }j\ne 0,\end{array}$$

so that, in general,

$$\sum _{k=0}^{+\infty }|a_{k,j}|\le 2\alpha +2\beta +\xi .$$

Moreover,

$$\begin{array}{cc}\hfill \sum _{k=0}^{+\infty }|a_{0,k}|=\alpha +\beta +\xi ,& \hfill \sum _{k=0}^{+\infty }|a_{j,k}|=2\alpha +2\beta +\xi \text{ for }j\ne 0\end{array}$$

so that, in general,

$$\sum _{k=0}^{+\infty }|a_{j,k}|\le 2\alpha +2\beta +\xi .$$

By Schur’s test, we have that $A_{\xi }$ is a bounded operator on $l^2$ and

$$\parallel A_{\xi }\parallel \le 2(\alpha +\beta )+\xi .$$

 ☐

Observe that, if $\xi =0$, the operator $A_0$ is the same of system (12). Let us also observe that by Corollary 1 there locally exists a unique solution. Moreover, if we show that a solution is uniformly bounded, such solution is unique. Now, we are ready to adapt Lemma 2.1 of [14] to the fractional case.

Theorem 3.

Let ${\tilde{N}}^{\nu }(t)$ be the number of customers in a fractional M/M/1 queue with arrival rate $\alpha z_1$ and service rate $\frac{\beta }{z_1}$ such that $\mathbb{P}({\tilde{N}}^{\nu }(0)=0)=1$ and consider N a random variable independent from ${\tilde{N}}^{\nu }(t)$ whose state probabilities $q_n$ are defined in Equation (31). Define

$$M^{\nu }(t):=min\{{\tilde{N}}^{\nu }(t),N\},t\ge 0.$$

Then, the state probabilities of $M^{\nu }(t)$ are the unique solutions of (29).

Moreover, $M^{\nu }(t)\stackrel{d}{=}{N}_{\xi }^{\nu }(t)$, where $\stackrel{d}{=}$ is the equality in distribution, and then $p_n^{\nu ,\xi }(t),n=0,1,\dots$ are the unique solutions of (29).

Proof. 

Define $p_n^{*,\nu }(t)=\mathbb{P}(M^{\nu }(t)=n)$ and ${\tilde{p}}_n^{\nu }(t)=\mathbb{P}({\tilde{N}}^{\nu }(t)=n)$. Since ${\tilde{N}}^{\nu }(t)$ and N are independent, then

$$p_n^{*,\nu }(t)=\mathbb{P}(N=n)\mathbb{P}({\tilde{N}}^{\nu }(t)\ge n)+\mathbb{P}({\tilde{N}}^{\nu }(t)=n)\mathbb{P}(N>n),$$

which, by using the definitions of ${\tilde{p}}_n^{\nu }(t)$ and $q_n$, becomes

$$p_n^{*,\nu }(t)=q_n\sum _{k=n}^{+\infty }{\tilde{p}}_k^{\nu }(t)+\left(\sum _{k=n+1}^{+\infty }{q}_n\right){\tilde{p}}_n^{\nu }(t).$$

(36)

Moreover, by using Equation (31), we have

$$\sum _{k=n+1}^{+\infty }{q}_n=\left(1-\frac{1}{z_1}\right)\sum _{k=n+1}^{+\infty }{\left(\frac{1}{z_1}\right)}^k=\left(1-\frac{1}{z_1}\right){\left(\frac{1}{z_1}\right)}^{n+1}\sum _{k=0}^{+\infty }{\left(\frac{1}{z_1}\right)}^k={\left(\frac{1}{z_1}\right)}^{n+1}$$

(37)

and then, substituting Equation (37) in (36), we obtain

$$p_n^{*,\nu }(t)=q_n\sum _{k=n}^{+\infty }{\tilde{p}}_k^{\nu }(t)+{\left(\frac{1}{z_1}\right)}^{n+1}{\tilde{p}}_n^{\nu }(t).$$

(38)

We want to show that $M^{\nu }(t)=N^{\nu }(t)$. Since, by definition, $p_n^{*,\nu }(t)$ are non-negative and ${\sum }_{n=0}^{+\infty }{p}_n^{*,\nu }(t)=1$, they are uniformly bounded in $l^2(\mathbb{R})$. Thus, we only need to show that ${\mathit{p}}^{*,\nu }(t)={(p_n^{*,\nu }(t))}_{n\ge 0}$ solves system (35).

The initial conditions are easily verified, so we only need to verify the differential relations. Observe that

$$p_0^{*,\nu }(t)=q_0+\frac{1}{z_1}{\tilde{p}}_0^{\nu }(t)$$

and then, applying the Caputo derivative operator, we obtain

$${}_0^C{D}_t^{\nu }{p}_0^{*,\nu }(t)=\frac{1}{z_1}{}_0^C{D}_t^{\nu }{\tilde{p}}_0^{\nu }(t).$$

Since ${\tilde{p}}_0^{\nu }(t)$ is a solution of system (5) with rates $\alpha z_1$ and $\frac{\beta }{z_1}$, we have

$${}_0^C{D}_t^{\nu }{p}_0^{*,\nu }(t)=-\alpha {\tilde{p}}_0^{\nu }(t)+\frac{\beta }{z_1^2}{\tilde{p}}_1^{\nu }(t).$$

Observe also that

$$p_1^{*,\nu }(t)=q_1(1-{\tilde{p}}_0^{\nu }(t))+{\left(\frac{1}{z_1}\right)}^2{\tilde{p}}_1^{\nu }(t)$$

so we have

$$\begin{array}{c}-(\alpha +\xi )p_0^{*,\nu }(t)+\beta p_1^{*,\nu }(t)+\xi \hfill \\ \hfill =-(\alpha +\xi )\left(q_0+\left(\frac{1}{z_1}\right){\tilde{p}}_0^{\nu }(t)\right)+\beta \left[q_1(1-{\tilde{p}}_0^{\nu }(t))+{\left(\frac{1}{z_1}\right)}^2{\tilde{p}}_1^{\nu }(t)\right]+\xi .\end{array}$$

After some calculations, we obtain

$$-(\alpha +\xi )p_0^{*,\nu }(t)+\beta p_1^{*,\nu }(t)+\xi =-(\alpha +\xi )q_0-\frac{\alpha +\xi }{z_1}{\tilde{p}}_0^{\nu }(t)+\beta q_1-\beta q_1{\tilde{p}}_0^{\nu }(t)+\frac{\beta }{z_1^2}{\tilde{p}}_1^{\nu }(t)+\xi .$$

Let us remark that

$$\begin{array}{cc}\hfill q_0=1-\frac{1}{z_1},& \hfill q_1=\left(1-\frac{1}{z_1}\right)\left(\frac{1}{z_1}\right),\end{array}$$

so we have

$$\begin{array}{c}-(\alpha +\xi )p_0^{*,\nu }(t)+\beta p_1^{*,\nu }(t)+\xi \hfill \\ \hfill =\frac{-\alpha z_1^2+(\alpha +\beta +\xi )z_1-\beta }{z_1^2}+\frac{(-(\alpha +\beta +\xi )z_1+\beta ){\tilde{p}}_0^{\nu }(t)}{z_1^2}+\frac{\beta }{z_1^2}{\tilde{p}}_1^{\nu }(t).\end{array}$$

By using Equations (32) and (34), we obtain

$$-(\alpha +\xi )p_0^{*,\nu }(t)+\beta p_1^{*,\nu }(t)+\xi =-\alpha {\tilde{p}}_0^{\nu }(t)+\frac{\beta }{z_1^2}{\tilde{p}}_1^{\nu }(t)={}_0^C{D}_t^{\nu }{p}_0^{*,\nu }(t).$$

Rewrite now Equation (38) in the form

$$p_n^{*,\nu }(t)=q_n\left(1-\sum _{k=0}^{+\infty }{\tilde{p}}_k^{\nu }(t)\right)+{\left(\frac{1}{z_1}\right)}^{n+1}{\tilde{p}}_n^{\nu }(t)$$

(39)

and then apply a Caputo derivative operator to obtain

$${}_0^C{D}_t^{\nu }{p}_n^{*,\nu }(t)=-q_n\sum _{k=1}^{+\infty }{}_0^C{D}_t^{\nu }{\tilde{p}}_k^{\nu }(t)-q_n{}_0^C{D}_t^{\nu }{\tilde{p}}_0^{\nu }(t)+{\left(\frac{1}{z_1}\right)}^{n+1}{}_0^C{D}_t^{\nu }{\tilde{p}}_n^{\nu }(t).$$

Since ${\tilde{p}}_n^{\nu }(t)$ is a solution of system (5) with birth rate $\alpha z_1$ and death rate $\frac{\beta }{z_1}$, then we have

$$\begin{array}{c}\hfill \begin{array}{cc}\hfill {}_0^C{D}_t^{\nu }{p}_n^{*,\nu }(t)& =q_n\left(\alpha z_1+\frac{\beta }{z_1}\right)\sum _{k=1}^{n-1}{\tilde{p}}_k^{\nu }-q_n\alpha z_1\sum _{k=0}^{n-2}{\tilde{p}}_k^{\nu }(t)-\frac{\beta }{z_1}{q}_n\sum _{k=2}^n{\tilde{p}}_k^{\nu }(t)\hfill \\ & +\alpha z_1{q}_n{\tilde{p}}_0(t)-\frac{\beta }{z_1}{q}_n{\tilde{p}}_1^{\nu }-{\left(\frac{1}{z_1}\right)}^{n+1}\left(\alpha z_1+\frac{\beta }{z_1}\right){\tilde{p}}_n^{\nu }(t)\hfill \\ & +\alpha {\left(\frac{1}{z_1}\right)}^n{\tilde{p}}_{n-1}^{\nu }(t)+\beta {\left(\frac{1}{z_1}\right)}^{n+2}{\tilde{p}}_{n+1}^{\nu }(t).\hfill \end{array}\end{array}$$

Remark that, by using Equation (39),

$$\begin{array}{c}\hfill \begin{array}{cc}\hfill -(\alpha +\beta +\xi )p_n^{*,\nu }(t)& +\alpha p_{n-1}^{*,\nu }(t)+\beta p_{n+1}^{*,\nu }(t)=\hfill \\ & -(\alpha +\beta +\xi )\left(q_n\left(1-\sum _{k=0}^{n-1}{\tilde{p}}_k^{\nu }(t)\right)+{\left(\frac{1}{z_1}\right)}^{n+1}{\tilde{p}}_n^{\nu }(t)\right)\hfill \\ & +\alpha \left(q_{n-1}\left(1-\sum _{k=0}^{n-2}{\tilde{p}}_k^{\nu }(t)\right)+{\left(\frac{1}{z_1}\right)}^n{\tilde{p}}_{n-1}^{\nu }(t)\right)\hfill \\ & +\beta \left(q_{n+1}\left(1-\sum _{k=0}^n{\tilde{p}}_k^{\nu }(t)\right)+{\left(\frac{1}{z_1}\right)}^{n+2}{\tilde{p}}_{n+1}^{\nu }(t)\right).\hfill \end{array}\end{array}$$

Then, recalling that by definition $q_{n-1}=z_1{q}_n$ and $q_{n+1}=\frac{q_n}{z_1}$ and doing some calculations, we have

$$\begin{array}{c}\hfill \begin{array}{cc}\hfill -(\alpha +\beta +\xi )p_n^{*,\nu }(t)& +\alpha p_{n-1}^{*,\nu }(t)+\beta p_{n+1}^{*,\nu }(t)=\hfill \\ & (\alpha +\beta +\xi )q_n\sum _{k=1}^{n-1}{\tilde{p}}_k^{\nu }(t)-\alpha z_1{q}_n\sum _{k=0}^{n-2}{\tilde{p}}_k^{\nu }(t)-\frac{\beta }{z_1}{q}_n\sum _{k=2}^n{\tilde{p}}_k^{\nu }(t)\hfill \\ & +\frac{(\alpha +\beta +\xi )z_1-\beta }{z_1}{q}_n{\tilde{p}}_0^{\nu }(t)-\frac{\beta }{z_1}{q}_n{\tilde{p}}_1^{\nu }(t)-(\alpha +\beta +\xi ){\left(\frac{1}{z_1}\right)}^{n+1}{\tilde{p}}_n^{\nu }(t)\hfill \\ & +\alpha {\left(\frac{1}{z_1}\right)}^n{\tilde{p}}_{n-1}^{\nu }(t)+\beta {\left(\frac{1}{z_1}\right)}^{n+2}{\tilde{p}}_{n+1}^{\nu }(t)+\frac{\alpha z_1^2-(\alpha +\beta +\xi )z_1+\beta }{z_1}{q}_n.\hfill \end{array}\end{array}$$

Finally, by using Equations (32), (33) and (34), we have

$$\begin{array}{c}\hfill \begin{array}{cc}\hfill -(\alpha +\beta +\xi )p_n^{*,\nu }(t)& +\alpha p_{n-1}^{*,\nu }(t)+\beta p_{n+1}^{*,\nu }(t)=\hfill \\ & \left(\alpha z_1+\frac{\beta }{z_1}\right)q_n\sum _{k=1}^{n-1}{\tilde{p}}_k^{\nu }(t)-\alpha z_1{q}_n\sum _{k=0}^{n-2}{\tilde{p}}_k^{\nu }(t)-\frac{\beta }{z_1}{q}_n\sum _{k=2}^n{\tilde{p}}_k^{\nu }(t)\hfill \\ & +\alpha z_1{q}_n{\tilde{p}}_0^{\nu }(t)-\frac{\beta }{z_1}{q}_n{\tilde{p}}_1^{\nu }(t)-\left(\alpha z_1+\frac{\beta }{z_1}\right){\left(\frac{1}{z_1}\right)}^{n+1}{\tilde{p}}_n^{\nu }(t)\hfill \\ & +\alpha {\left(\frac{1}{z_1}\right)}^n{\tilde{p}}_{n-1}^{\nu }(t)+\beta {\left(\frac{1}{z_1}\right)}^{n+2}{\tilde{p}}_{n+1}^{\nu }(t)={}_0^C{D}_t^{\nu }{p}_n^{*,\nu }(t).\hfill \end{array}\end{array}$$

We have shown that the state probabilities $p_n^{*,\nu }(t)$ of $M^{\nu }(t)$ are the unique solutions of system (29). Now, we need to show that $M^{\nu }(t)\stackrel{d}{=}{N}_{\xi }^{\nu }(t)$. To do this, consider ${\tilde{N}}^1(t)$ a classical M/M/1 queue with arrival rate $\alpha z_1$ and service rate $\frac{\beta }{z_1}$, N a random variable independent from ${\tilde{N}}^{\nu }(t)$ and ${\tilde{N}}^1(t)$ with probability masses $q_n$ and finally $L_{\nu }(t)$ the inverse of a $\nu$-stable subordinator which is independent from N and ${\tilde{N}}^1(t)$. Define also $M^1(t)=min\{{\tilde{N}}^1(t),N\}$. By Lemma $2.1$ of [14], we know that $M^1(t)\stackrel{d}{=}{N}_{\xi }^1(t)$, so $M^1(L_{\nu }(t))\stackrel{d}{=}{N}_{\xi }^1(L_{\nu }(t))\stackrel{d}{=}{N}_{\xi }^{\nu }(t)$. However, by definition, we know that ${\tilde{N}}^1(L_{\nu }(t))\stackrel{d}{=}{\tilde{N}}^{\nu }(t)$, thus finally

$$M^{\nu }(t)\stackrel{d}{=}{M}^1(L_{\nu }(t))\stackrel{d}{=}{N}_{\xi }^1(L_{\nu }(t))\stackrel{d}{=}{N}_{\xi }^{\nu }(t).$$

 ☐

5.2. State Probabilities for the Fractional M/M/1 with Catastrophes

Since we have defined $N_{\xi }^{\nu }(t)\stackrel{d}{=}{N}_{\xi }^1(L_{\nu }(t))$, where $L_{\nu }(t)$ is the inverse of a $\nu$-stable subordinator, which is independent from $N_{\xi }^1(t)$, we can use such definition and Theorem 3 with the results obtained in [14] to study the state probabilities of $N_{\xi }^{\nu }(t)$. In particular, we refer to the formula

$$\begin{array}{c}\hfill \begin{array}{cc}\hfill p_n^{1,\xi }(t)=q_n& +\sum _{m=0}^{+\infty }\sum _{r=0}^{m+n}\frac{C_{n,m,r}^1}{(m+r-1)!}{t}^{m+r-1}{e}^{-(\alpha +\beta +\xi )t}\hfill \\ & +\sum _{m=0}^{+\infty }\sum _{r=m+n+1}^{+\infty }\frac{C_{n,m,r}^2}{(m+r-1)!}{t}^{m+r-1}{e}^{-(\alpha +\beta +\xi )t},\hfill \end{array}\end{array}$$

(40)

where

$$\begin{array}{c}\hfill \begin{array}{c}{C}_{n,m,r}^1=\frac{z_1-1}{(z_1-z_2)z_1^{n+m+1-r}}\left(\genfrac{}{}{0pt}{}{m+r}{r}\right)\frac{m-r}{m+r}{\beta }^m{\alpha }^{r-1},\hfill \\ C_{n,m,r}^2=\frac{1-z_2}{(z_1-z_2)z_2^{n+m+1-r}}\left(\genfrac{}{}{0pt}{}{m+r}{r}\right)\frac{r-m}{r+m}{\beta }^m{\alpha }^{r-1}.\hfill \end{array}\end{array}$$

(41)

By using such formula, we can show the following:

Theorem 4.

For any $t>0$ and $n=0,1,\dots$, we have

$$\begin{array}{c}\hfill \begin{array}{cc}\hfill p_n^{\nu ,\xi }(t)=q_n& +\sum _{m=0}^{+\infty }\sum _{r=0}^{m+n}{C}_{n,m,r}^1{t}^{\nu (m+r-1)}{E}_{\nu ,\nu (m+r-1)+1}^{m+r}(-(\alpha +\beta +\xi )t^{\nu })\hfill \\ & +\sum _{m=0}^{+\infty }\sum _{r=m+n+1}^{+\infty }{C}_{n,m,r}^2{t}^{\nu (m+r-1)}{E}_{\nu ,\nu (m+r-1)+1}^{m+r}(-(\alpha +\beta +\xi )t^{\nu }),\hfill \end{array}\end{array}$$

(42)

where $C_{n,m,r}^i$ are defined in (41).

Proof. 

From $N_{\xi }^{\nu }(t)\stackrel{d}{=}{N}_{\xi }^1(L_{\nu }(t))$, we have

$$p_n^{\nu ,\xi }(t)={\int }_0^{+\infty }{p}_n^{1,\xi }(y)\mathbb{P}(L_{\nu }(t)\in dy)$$

and then, by using formula (40),

$$\begin{array}{cc}\hfill p_n^{\nu ,\xi }(t)=q_n& +\sum _{m=0}^{+\infty }\sum _{r=0}^{m+n}\frac{C_{n,m,r}^1}{(m+r-1)!}{\int }_0^{+\infty }{y}^{m+r-1}{e}^{-(\alpha +\beta +\xi )y}\mathbb{P}(L_{\nu }(t)\in dy)\hfill \\ & +\sum _{m=0}^{+\infty }\sum _{r=m+n+1}^{+\infty }\frac{C_{n,m,r}^2}{(m+r-1)!}{\int }_0^{+\infty }{y}^{m+r-1}{e}^{-(\alpha +\beta +\xi )y}\mathbb{P}(L_{\nu }(t)\in dy).\hfill \end{array}$$

Taking the Laplace transform and using Equation (22), we obtain

$$\begin{array}{cc}\hfill {\pi }_n^{\nu ,\xi }(s)=q_n& +\sum _{m=0}^{+\infty }\sum _{r=0}^{m+n}\frac{C_{n,m,r}^1}{(m+r-1)!}{s}^{\nu -1}{\int }_0^{+\infty }{y}^{m+r-1}{e}^{-(\alpha +\beta +\xi +s^{\nu })y}dy\hfill \\ & +\sum _{m=0}^{+\infty }\sum _{r=m+n+1}^{+\infty }\frac{C_{n,m,r}^2}{(m+r-1)!}{s}^{\nu -1}{\int }_0^{+\infty }{y}^{m+r-1}{e}^{-(\alpha +\beta +\xi +s^{\nu })y}dy\hfill \end{array}$$

and then, integrating

$$\begin{array}{cc}\hfill {\pi }_n^{\nu ,\xi }(s)=q_n& +\sum _{m=0}^{+\infty }\sum _{r=0}^{m+n}{C}_{n,m,r}^1\frac{s^{\nu -1}}{{(\alpha +\beta +\xi +s^{\nu })}^{m+r}}\hfill \\ & +\sum _{m=0}^{+\infty }\sum _{r=m+n+1}^{+\infty }{C}_{n,m,r}^2\frac{s^{\nu -1}}{{(\alpha +\beta +\xi +s^{\nu })}^{m+r}}.\hfill \end{array}$$

Finally, by using Equation (A2), we obtain

$$\begin{array}{cc}\hfill p_n^{\nu ,\xi }(t)=q_n& +\sum _{m=0}^{+\infty }\sum _{r=0}^{m+n}{C}_{n,m,r}^1{t}^{\nu (m+r-1)}{E}_{\nu ,\nu (m+r-1)+1}^{m+r}(-(\alpha +\beta +\xi )t^{\nu })\hfill \\ & +\sum _{m=0}^{+\infty }\sum _{r=m+n+1}^{+\infty }{C}_{n,m,r}^2{t}^{\nu (m+r-1)}{E}_{\nu ,\nu (m+r-1)+1}^{m+r}(-(\alpha +\beta +\xi )t^{\nu }).\hfill \end{array}$$

 ☐

Remark 2.

From formula (42), we can easily see that ${lim}_{t\to +\infty }{p}_n^{\nu ,\xi }(t)=q_n$ so, as we expected, the steady-state probabilities are the same as the classical ones. For such reason, we can say that the fractional behaviour is influential only in the transient state of the queue.

Remark 3.

As $\nu \to 1$, by using

$$E_{1,m+r}^{m+r}(-(\alpha +\beta +\xi )t)=\frac{e^{-(\alpha +\beta +\xi )t}}{(m+r-1)!},$$

we obtain that ${lim}_{\nu \to 1}{p}_n^{\nu ,\xi }(t)=p_n^{1,\xi }(t)$.

Remark 4.

If $\alpha <\beta$ and $\xi =0$, then $z_1=\frac{\beta }{\alpha }$ and $z_2=1$. For such reason, $q_n=\left(1-\frac{\alpha }{\beta }\right){\left(\frac{\alpha }{\beta }\right)}^n$, $C_{n,m,r}^1={\left(\frac{\alpha }{\beta }\right)}^n\frac{m-r}{m+r}\left(\genfrac{}{}{0pt}{}{m+r}{m}\right){\alpha }^m{\beta }^{r-1}$ and $C_{n,m,r}^2=0$. Then, we have that $p_n^{\nu ,\xi }(t)$ of Equation (42) has the form of Equation (6).

If $\alpha >\beta$ and $\xi \to 0$ then $z_1=1$ and $z_2=\frac{\beta }{\alpha }$. In such case, $q_n=0$, $C_{n,m,r}^1=0$ and $C_{n,m,r}^2={\alpha }^{n+m}{\beta }^{r-n-1}\left(\genfrac{}{}{0pt}{}{m+r}{m}\right)\frac{r-m}{m+r}$. For such case, we have

$$\underset{\xi \to 0}{lim}{p}_n^{\nu ,\xi }(t)={\left(\frac{\alpha }{\beta }\right)}^n\sum _{m=0}^{+\infty }\sum _{r=m+n+1}^{+\infty }{\alpha }^m{\beta }^{r-1}\left(\genfrac{}{}{0pt}{}{m+r}{m}\right)\frac{r-m}{m+r}{t}^{\nu (m+r-1)}{E}_{\nu ,\nu (m+r-1)+1}^{m+r},$$

which is not recognizable as a previously obtained formula. This is due to the fact that the formula

$$\underset{\xi \to 0}{lim}{p}_n^{1,\xi }(t)=\frac{e^{-(\alpha +\beta )t}}{\beta t}{\left(\frac{\alpha }{\beta }\right)}^n\sum _{r=n+1}^{+\infty }r{\left(\frac{\beta }{\alpha }\right)}^{\frac{r}{2}}{I}_r(2\sqrt{\alpha \beta }t)$$

(43)

(which is the one that is obtained from (42) as $\nu =1$ and $\alpha >\beta$, as done in [14]) has no known equivalent in the fractional case. It is also interesting to observe that in [8] another representation of the Laplace transform of $p_n^{\nu }(t)$ is given in formula $2.40$, which is not easily invertible, but has been obtained by using (43) instead of Sharma’s representation of $p_n^1(t)$ ([2])

$$p_n^1(t)=\left(1-\frac{\alpha }{\beta }\right){\left(\frac{\alpha }{\beta }\right)}^n+e^{-(\alpha +\beta )t}{\left(\frac{\alpha }{\beta }\right)}^n\sum _{r=0}^{+\infty }\frac{{(\alpha t)}^r}{r!}\sum _{m=0}^{k+r}(r-m)\frac{{(\beta t)}^{m-1}}{m!}.$$

5.3. Distribution of the Busy Period

Let $B^{\nu }$ denote the duration of the busy period and $B^{\nu }(t)=\mathbb{P}(B^{\nu }\le t)$ be its probability distribution function. Let us observe that, if we pose $N^{\nu }(0)=1$, then the queue empties within t if and only if a catastrophe occurs within t or otherwise the queue empties without catastrophes within t. Let us remark that, if there is no occurrence of catastrophes, the queue behaves as a fractional M/M/1. Let us define $K^{\nu }$ as the duration of a busy period for a fractional M/M/1 queue without catastrophes, ${\mathsf{\Xi }}^{\nu }$ the time of first occurrence of a catastrophe for a non empty queue and $K^{\nu }(t)=\mathbb{P}(K^{\nu }\le t)$ and ${\mathsf{\Xi }}^{\nu }(t)=\mathbb{P}({\mathsf{\Xi }}^{\nu }\le t)$ their probability distribution functions. Thus, we have

$$B^{\nu }(t)={\mathsf{\Xi }}^{\nu }(t)+(1-{\mathsf{\Xi }}^{\nu }(t))K^{\nu }(t).$$

(44)

Remark 5.

If we denote with $b^{\nu }(t)$, ${\xi }^{\nu }(t)$ and $k^{\nu }(t)$ the probability density functions of $B^{\nu }$, ${\mathsf{\Xi }}^{\nu }$ and $K^{\nu }$, we have, by deriving formula (44),

$$b^{\nu }(t)={\xi }^{\nu }(t)(1-K^{\nu }(t))+(1-{\mathsf{\Xi }}^{\nu }(t))k^{\nu }(t),$$

which, for $\nu =1$, is formula (17) of [14].

By using formula (44), we can finally show:

Theorem 5.

Let $B^{\nu }$ be the duration of the busy period of a fractional M/M/1 queue with catastrophes and $B^{\nu }(t)=\mathbb{P}(B^{\nu }\le t)$. Then,

$$B^{\nu }(t)=1-E_{\nu }(-\xi t^{\nu })\sum _{n=1}^{+\infty }\sum _{m=0}^{+\infty }{C}_{n,m}{t}^{\nu (n+2m-1)}{E}_{\nu ,\nu (n+2m-1)+1}^{n+2m}(-(\alpha +\beta )t^{\nu }),$$

(45)

where $C_{n,m}$ is given in (14).

Proof. 

Observe that, by formula (30), we have

$${\mathsf{\Xi }}^{\nu }(t)=c_0^{\nu }(t)=1-E_{\nu }(-\xi t^{\nu })$$

and by formula (13) we also have a closed form of $K^{\nu }(t)$. Thus, by using formula (44), we obtain Equation (45). ☐

5.4. Distribution of the Time of the First Occurrence of a Catastrophe

We already know that if the queue starts from a non-empty state, then the occurrence of the catastrophes is a Mittag–Leffler distribution. However, we are interested in such distribution as the queue starts being empty. To do that, we will need some auxiliary discrete processes.

Theorem 6.

Let ${\mathcal{D}}^{\nu }$ be the time of first occurrence of a catastrophe as $\mathbb{P}(N^{\nu }(0)=0)=1$ and let ${\mathcal{D}}^{\nu }(t)=\mathbb{P}({\mathcal{D}}^{\nu }\le t)$. Then,

$${\mathcal{D}}^{\nu }(t)=1-\sum _{j=1}^{+\infty }\sum _{m=0}^{+\infty }{C}_{m,j}{t}^{\nu (2m+j-1)}{E}_{\nu ,\nu (2m+j-1)+1}^{2m+j}[-(\alpha +\beta +\xi )t^{\nu }],$$

(46)

where

$$C_{m,j}=\frac{j}{2m+j}\frac{{(\beta +\xi )}^j-{\alpha }^j}{\beta +\xi -\alpha }\left(\genfrac{}{}{0pt}{}{2m+j}{m}\right){(\alpha \beta )}^m.$$

Proof. 

Following the lines of [14], let us consider the process ${\overline{N}}^{\nu }(t)$ with state space $\{-1,0,1,2,\dots \}$ such that $\mathbb{P}({\overline{N}}^{\nu }(t)=0)=1$ and posing $r_n(t)=\mathbb{P}({\overline{N}}^{\nu }(t)=n),n\ge -1$ as its state probability, we have

$$\left\{\begin{array}{c}{}_0^C{D}_t^{\nu }{r}_{-1}^{\nu }(t)=\xi [1-r_{-1}^{\nu }(t)-r_0^{\nu }(t)],\hfill \\ {}_0^C{D}_t^{\nu }{r}_0^{\nu }(t)=-\alpha r_0^{\nu }(t)+\beta r_1^{\nu }(t),\hfill \\ {}_0^C{D}_t^{\nu }{r}_n^{\nu }(t)=-(\alpha +\beta +\xi )r_n^{\nu }(t)+\alpha r_{n-1}^{\nu }(t)+\beta r_{n+1}^{\nu }(t),\hfill & n\ge 1,\hfill \\ r_n^{\nu }(0)={\delta }_{n,0},\hfill & n\ge -1.\hfill \end{array}\right.$$

(47)

Let us remark that such process represents our queue until a catastrophe occurs: in such case, instead of emptying the queue, the state of the process becomes $-1$, which is an absorbent state. With such interpretation, we can easily observe that ${\mathcal{D}}^{\nu }(t)=r_{-1}^{\nu }(t)$.

In order to determine $r_n^{\nu }(t)$, we will first show that ${\overline{N}}^{\nu }(t)\stackrel{d}{=}{\overline{N}}^1(L_{\nu }(t))$ where $L_{\nu }(t)$ is the inverse of a $\nu$-stable subordinator which is independent from ${\overline{N}}^1$. To do that, let us consider ${\overline{N}}^{\nu }(t)+1$ instead of ${\overline{N}}^{\nu }(t)$. Let us remark that $\mathbb{P}({\overline{N}}^{\nu }(t)+1=n)=r_{n-1}^{\nu }(t)$. Let $G^{\nu }(z,t)={\sum }_{n=0}^{+\infty }{z}^n{r}_{n-1}^{\nu }(t)$ be the probability generating function of ${\overline{N}}^{\nu }(t)+1$. Multiplying the third sequence of equations in (47) with $z^{n+1}$ and then, summing all these equations, we have

$${}_0^C{D}_t^{\nu }\left(\sum _{n=2}^{+\infty }{z}^n{r}_{n-1}^{\nu }(t)\right)=-(\alpha +\beta +\xi )\sum _{n=2}^{+\infty }{z}^n{r}_{n-1}^{\nu }(t)+\alpha \sum _{n=2}^{+\infty }{z}^n{r}_{n-2}^{\nu }(t)+\beta \sum _{n=2}^{+\infty }{z}^n{r}_n^{\nu }(t).$$

(48)

Now observe that

$$\sum _{n=2}^{+\infty }{z}^n{r}_{n-1}^{\nu }(t)=\sum _{n=0}^{+\infty }{z}^n{r}_{n-1}^{\nu }(t)-r_{-1}^{\nu }(t)-z{r}_0^{\nu }(t)=G^{\nu }(z,t)-r_{-1}^{\nu }(t)-z{r}_0^{\nu }(t);$$

(49)

moreover,

$$\begin{array}{c}\sum _{n=2}^{+\infty }{z}^n{r}_{n-2}^{\nu }(t)=\sum _{n=1}^{+\infty }{z}^{n+1}{r}_{n-1}^{\nu }(t)=z\sum _{n=1}^{+\infty }{z}^n{r}_{n-1}^{\nu }(t)\hfill \\ \hfill =z[G^{\nu }(z,t)-r_{-1}^{\nu }(t)]=z[G^{\nu }(z,t)-r_{-1}^{\nu }(t)-z{r}_0^{\nu }(t)]+z^2{r}_0^{\nu }(t);\end{array}$$

(50)

finally,

$$\begin{array}{c}\sum _{n=2}^{+\infty }{z}^n{r}_n^{\nu }(t)=\sum _{n=3}^{+\infty }{z}^{n-1}{r}_{n-1}^{\nu }(t)=\frac{1}{z}\sum _{n=3}^{+\infty }{z}^n{r}_{n-1}^{\nu }(t)\hfill \\ =\frac{1}{z}[G^{\nu }(z,t)-r_{-1}^{\nu }(t)-z{r}_0^{\nu }(t)-z^2{r}_1^{\nu }(t)]\\ \hfill =\frac{1}{z}[G^{\nu }(z,t)-r_{-1}^{\nu }(t)-z{r}_0^{\nu }(t)]-z{r}_1^{\nu }(t).\end{array}$$

(51)

Using Equations (49), (50) and (51) in Equation (48), we obtain

$$\begin{array}{c}{}_0^C{D}_t^{\nu }[G^{\nu }(z,t)-r_{-1}^{\nu }(t)-z{r}_0^{\nu }(t)]\hfill \\ =\left[\alpha z-(\alpha +\beta +\xi )+\frac{\beta }{z}\right][G^{\nu }(z,t)-r_{-1}^{\nu }(t)-z{r}_0^{\nu }(t)]\\ \hfill +\alpha z^2{r}_0^{\nu }(t)-\beta z{r}_1^{\nu }(t).\end{array}$$

(52)

Finally, by using the first and the second equation of Equation (47) in Equation (52), we obtain

$$\begin{array}{c}{}_0^C{D}_t^{\nu }{G}^{\nu }(z,t)=\left[\alpha z-(\alpha +\beta +\xi )+\frac{\beta }{z}\right][G^{\nu }(z,t)-r_{-1}^{\nu }(t)-z{r}_0^{\nu }(t)]\hfill \\ \hfill +\alpha z(z-1)r_0^{\nu }(t)+\xi [1-r_{-1}^{\nu }(t)-r_0^{\nu }(t)].\end{array}$$

We have obtained that the probability generating function $G^{\nu }(z,t)$ of ${\overline{N}}^{\nu }(t)+1$ solves the Cauchy problem

$$\left\{\begin{array}{c}\begin{array}{c}\hfill z{}_0^C{D}_t^{\nu }{G}^{\nu }(z,t)=\left[\alpha z^2-(\alpha +\beta +\xi )z+\beta \right][G^{\nu }(z,t)-r_{-1}^{\nu }(t)-z{r}_0^{\nu }(t)]\\ \hfill +\alpha z^2(z-1)r_0^{\nu }(t)+\xi z[1-r_{-1}^{\nu }(t)-r_0^{\nu }(t)],\end{array}\hfill \\ G^{\nu }(z,0)=z,\hfill \end{array}\right.$$

(53)

that, for $\nu =1$, becomes

$$\left\{\begin{array}{c}\begin{array}{c}\hfill z\frac{d}{dt}{G}^1(z,t)=\left[\alpha z^2-(\alpha +\beta +\xi )z+\beta \right][G^1(z,t)-r_{-1}^1(t)-z{r}_0^1(t)]\\ \hfill +\alpha z^2(z-1)r_0^1(t)+\xi z[1-r_{-1}^1(t)-r_0^1(t)],\end{array}\hfill \\ G^1(z,0)=z.\hfill \end{array}\right.$$

(54)

Let ${\tilde{G}}^{\nu }(z,s)$, ${\tilde{r}}_0^{\nu }(s)$ and ${\tilde{r}}_{-1}^{\nu }(s)$ be the Laplace transforms of $G^{\nu }(z,t)$, $r_0^{\nu }(t)$ and $r_{-1}^{\nu }(t)$ and let us take the Laplace transform in Equation (53) to obtain

$$\begin{array}{c}\hfill z[s^{\nu }{\tilde{G}}^{\nu }(z,s)-s^{\nu -1}z]=\left[\alpha z^2-(\alpha +\beta +\xi )z+\beta \right][{\tilde{G}}^{\nu }(z,s)-{\tilde{r}}_{-1}^{\nu }(s)-z{\tilde{r}}_0^{\nu }(s)]\\ \hfill +\alpha z^2(z-1){\tilde{r}}_0^{\nu }(s)+\xi z\left[\frac{1}{s}-{\tilde{r}}_{-1}^{\nu }(s)-{\tilde{r}}_0^{\nu }(s)\right].\end{array}$$

(55)

Now, let us remark that ${\overline{N}}^{\nu }(t)+1\stackrel{d}{=}{\overline{N}}^1(L_{\nu }(t))+1$ if and only if for all $n\ge 0$:

$$r_{n-1}^{\nu }(t)={\int }_0^{+\infty }{r}_{n-1}^1(y)\mathbb{P}(L_{\nu }(t)\in dy)$$

(56)

that is to say if and only if

$$G^{\nu }(z,t)={\int }_0^{+\infty }{G}^1(z,y)\mathbb{P}(L_{\nu }(t)\in dy).$$

Taking Laplace transform and using Equation (22), we obtain

$$\begin{array}{c}\hfill \begin{array}{c}{\tilde{G}}^{\nu }(z,s)=s^{\nu -1}{\int }_0^{+\infty }{G}^1(z,y)e^{-y{s}^{\nu }}dy,\hfill \\ {\tilde{r}}_{-1}^{\nu }(s)=s^{\nu -1}{\int }_0^{+\infty }{r}_{-1}^1(y)e^{-y{s}^{\nu }}dy,\hfill \\ {\tilde{r}}_0^{\nu }(s)=s^{\nu -1}{\int }_0^{+\infty }{r}_0^1(y)e^{-y{s}^{\nu }}dy.\hfill \end{array}\end{array}$$

(57)

Thus, by substituting the formulas (57) in (55), we obtain

$$\begin{array}{cc}\hfill z& \left[s^{\nu }{s}^{\nu -1}{\int }_0^{+\infty }{G}^1(z,y)e^{-y{s}^{\nu }}dy-s^{\nu -1}z\right]=[\alpha z^2-(\alpha +\beta +\xi )z+\beta ]\hfill \\ & \times \left[s^{\nu -1}{\int }_0^{+\infty }{G}^1(z,y)e^{-y{s}^{\nu }}dy-s^{\nu -1}{\int }_0^{+\infty }{r}_{-1}^1(y)e^{-y{s}^{\nu }}dy-z{s}^{\nu -1}{\int }_0^{+\infty }{r}_0^1(y)e^{-y{s}^{\nu }}dy\right]\hfill \\ & +\alpha z^2(z-1)s^{\nu -1}{\int }_0^{+\infty }{r}_0^1(y)e^{-y{s}^{\nu }}dy\hfill \\ & +\xi z\left[\frac{1}{s}-s^{\nu -1}{\int }_0^{+\infty }{r}_{-1}^1(y)e^{-y{s}^{\nu }}dy-s^{\nu -1}{\int }_0^{+\infty }{r}_0^1(y)e^{-y{s}^{\nu }}dy\right].\hfill \end{array}$$

Finally, multiplying with $\frac{1}{s^{\nu -1}}$, we have

$$\begin{array}{cc}\hfill z& \left[s^{\nu }{\int }_0^{+\infty }{G}^1(z,y)e^{-y{s}^{\nu }}dy-z\right]=[\alpha z^2-(\alpha +\beta +\xi )z+\beta ]\hfill \\ & \times \left[{\int }_0^{+\infty }{G}^1(z,y)e^{-y{s}^{\nu }}dy-{\int }_0^{+\infty }{r}_{-1}^1(y)e^{-y{s}^{\nu }}dy-z{\int }_0^{+\infty }{r}_0^1(y)e^{-y{s}^{\nu }}dy\right]\hfill \\ & +\alpha z^2(z-1){\int }_0^{+\infty }{r}_0^1(y)e^{-y{s}^{\nu }}dy\hfill \\ & +\xi z\left[\frac{1}{s^{\nu }}-{\int }_0^{+\infty }{r}_{-1}^1(y)e^{-y{s}^{\nu }}dy-{\int }_0^{+\infty }{r}_0^1(y)e^{-y{s}^{\nu }}dy\right].\hfill \end{array}$$

(58)

Now we know that ${\overline{N}}^{\nu }(t)\stackrel{d}{=}{\overline{N}}^1(L_{\nu }(t))$ if and only if Equation (58) is verified. For this reason, we only need to show such equation. To do that, remarking that ${\int }_0^{+\infty }{e}^{-y{s}^{\nu }}dy=\frac{1}{s^{\nu }}$, consider the right-hand side of Equation (58) and observe that

$$\begin{array}{cc}\hfill [\alpha z^2& -(\alpha +\beta +\xi )z+\beta ]\left[{\int }_0^{+\infty }{G}^1(z,y)e^{-y{s}^{\nu }}dy-{\int }_0^{+\infty }{r}_{-1}^1(y)e^{-y{s}^{\nu }}dy-z{\int }_0^{+\infty }{r}_0^1(y)e^{-y{s}^{\nu }}dy\right]\hfill \\ & +\alpha z^2(z-1){\int }_0^{+\infty }{r}_0^1(y)e^{-y{s}^{\nu }}dy+\xi z\left[\frac{1}{s^{\nu }}-{\int }_0^{+\infty }{r}_{-1}^1(y)e^{-y{s}^{\nu }}dy-{\int }_0^{+\infty }{r}_0^1(y)e^{-y{s}^{\nu }}dy\right]\hfill \\ & ={\int }_0^{+\infty }([\alpha z^2-(\alpha +\beta +\xi )z+\beta ][G^1(z,y)-r_{-1}^1(y)-z{r}_0^1(y)]\hfill \\ & +\alpha z^2(z-1)r_0^1(y)+\xi z[1-r_{-1}^1(y)-r_0^1(y)])e^{-y{s}^{\nu }}dy.\hfill \end{array}$$

Thus, by using Equation (54), we have

$$\begin{array}{cc}\hfill [\alpha z^2& -(\alpha +\beta +\xi )z+\beta ]\left[{\int }_0^{+\infty }{G}^1(z,y)e^{-y{s}^{\nu }}dy-{\int }_0^{+\infty }{r}_{-1}^1(y)e^{-y{s}^{\nu }}dy-z{\int }_0^{+\infty }{r}_0^1(y)e^{-y{s}^{\nu }}dy\right]\hfill \\ & +\alpha z^2(z-1){\int }_0^{+\infty }{r}_0^1(y)e^{-y{s}^{\nu }}dy+\xi z\left[\frac{1}{s^{\nu }}-{\int }_0^{+\infty }{r}_{-1}^1(y)e^{-y{s}^{\nu }}dy-{\int }_0^{+\infty }{r}_0^1(y)e^{-y{s}^{\nu }}dy\right]\hfill \\ & =z{\int }_0^{+\infty }\left(\frac{d}{dt}{G}^1(z,y)\right)e^{-y{s}^{\nu }}dy\hfill \\ & =z\left[{\int }_0^{+\infty }{G}^1(z,y)e^{-y{s}^{\nu }}dy-z\right],\hfill \end{array}$$

concluding the proof of our first claim.

From Theorem $3.1$ of [14], we know that

$$r_{-1}^1(t)=1-\sum _{j=1}^{+\infty }\sum _{m=0}^{+\infty }\frac{C_{m,j}}{(2m+j-1)!}{t}^{2m+j-1}{e}^{-(\alpha +\beta +\xi )t}$$

(59)

and, since we know that ${\overline{N}}^{\nu }(t)\stackrel{d}{=}{\overline{N}}^1(L_{\nu }(t))$, we can use (59) in (56) with $n=0$ to obtain:

$$r_{-1}^{\nu }(t)=1-\sum _{j=1}^{+\infty }\sum _{m=0}^{+\infty }\frac{C_{m,j}}{(2m+j-1)!}{\int }_0^{+\infty }{y}^{2m+j-1}{e}^{-(\alpha +\beta +\xi )y}\mathbb{P}(L_{\nu }(t)\in dy).$$

(60)

Taking the Laplace transform in (60) and using formula (22), we obtain

$${\tilde{r}}_{-1}^{\nu }(s)=\frac{1}{s}-\sum _{j=1}^{+\infty }\sum _{m=0}^{+\infty }\frac{C_{m,j}}{(2m+j-1)!}{s}^{\nu -1}{\int }_0^{+\infty }{y}^{2m+j-1}{e}^{-(\alpha +\beta +\xi +s^{\nu })y}dy$$

and then integrate

$${\tilde{r}}_{-1}^{\nu }(s)=\frac{1}{s}-\sum _{j=1}^{+\infty }\sum _{m=0}^{+\infty }{C}_{m,j}\frac{s^{\nu -1}}{{(\alpha +\beta +\xi +s^{\nu })}^{2m+j}}.$$

(61)

Finally, applying the inverse Laplace transform on Equation (61) and using formula (A2), we complete the proof. ☐

6. Conclusions

Our work focused on the transient behaviour of a fractional M/M/1 queue with catastrophes, deriving formulas for the state probabilities, the distribution of the busy period and the distribution of the time of the first occurrence of a catastrophe. This is a non-Markov generalization of the classical M/M/1 queue with catastrophes, obtained through a time-change. The introduction of fractional dynamics in the equations that master the behaviour of the queue led to a sort of transformation of the time scale. Fractional derivatives are global operators, so the state probabilities preserve memory of their past, eventually slowing down the entire dynamics. Indeed, we can see how Mittag–Leffler functions take place where in the classical case we expected to see exponentials. However, such fractional dynamic seems to affect only the transient behaviour, since we have shown in Remark 2 that the limit behaviour is the same.

The main difficulty that is linked with fractional queues (or in general time-changed queues) is the fact that one has to deal with non-local derivative operators, such as the Caputo derivative, losing Markov property. However, fractional dynamics and fractional processes are gaining attention, due to their wide range of applicability, from physics to finance, from computer science to biology. Moreover, time-changed processes have formed a thriving field of application in mathematical finance. Future works will focus on an extension of such results to $E_k/M/1$ and $M/E_k/1$ queues, or even to a generalization of fractional M/M/1 queue to a time-changed M/M/1 queue by using the inverse of any subordinator.