VMSTA Modern Stochastics: Theory and Applications 2351-60542351-60462351-6046 VTeXMokslininkų g. 2A, 08412 Vilnius, Lithuania VMSTA303 10.15559/26-VMSTA303 Research Article On Matrix Mittag Leffler distribution of PH/PH/1 system https://orcid.org/0000-0000-0000-0000 SinghRatnesh Kumarsingh.106@iitj.ac.in∗ https://orcid.org/0000-0002-1359-1861 VijayVivekvivek@iitj.ac.in Indian Institute of Technology Jodhpur, 342037 Jodhpur, India Corresponding author. 2026 16720260012278202526320262752026 © 2026 The Author(s). Published by VTeX2026 Open access article under the CC BY license.

We study a level-dependent PH/PH/1 queue with an infinite waiting room and derive explicit density formulas for the associated Matrix Mittag-Leffler (MML) random variables. The motivation arises from the use of Phase Type (PH) and Mittag-Leffler type distributions in modeling heavy-tailed and non-Markovian arrival and service processes in queueing systems, where classical exponential assumptions may fail to capture realistic behavior. Our goal is to obtain explicit density expressions for the MML random variables associated with the truncated queueing systems PH/M/ 1 ∼ n and M/PH/ 1 ∼ n, where n denotes the truncation level at which the queueing process is absorbed. We first analyze PH/M/1 queues in which the inter-arrival times are PH-distributed, focusing specifically on exponential and Erlang cases, and derive closed-form MML densities. We then consider M/PH/1 queues with PH-distributed service times and obtain analogous explicit density formulas. Several examples with PH-generators of different orders are presented to illustrate the results.

 Phase type distribution matrix Mittag-Leffler function queueing system cauchy residue principle 60K25 33E12 60J27 30E20 There is no source of funding for this research. Introduction

Matrix Mittag-Leffler (MML) distributions, defined by Albrecher et al. [2, 4], is a wide class of heavy-tailed distributions with some attractive mathematical properties. It is dense in the class of all lifetime distributions, like Phase Type (PH) distribution. This distribution is helpful in simulating a variety of real-world situations, where inter-arrival times of a queueing system have thicker tails. Additionally, the class of MML distributions is a fractional generalisation of a PH distribution (referred to it as a fractional PH distribution [3]). The PH distribution, the Erlang distribution, the exponential distribution, the Mittag-Leffler (ML) distribution, and the fractional Erlang distribution are special cases of this distribution(see, [16, 21, 18, 23, 9, 24, 2, 3]). As introduced by Albrecher et al. [2], PH distribution plays a very important role in MML distribution. A random variable X is said to have a MML distribution if its Laplace transform is given by E ( e − u X ) = π ( u α I − T ) − 1 t and is denoted as X ∼ M M L ( α , π , T ), where ( π , T ) is PH representation and 0 < α ≤ 1. Note that for α = 1, X follows a PH distribution. PH distributions are dense (weakly converging) on the positive real line, allowing them to accurately resemble a positive distribution. However, their light tail may be a challenge in applications that rely largely on tail behaviour [8]. Fitting heavy-tailed distributions with a PH distribution can result in the need for several stages and the resulting model may still fail to correctly represent the tail behaviour. MML distributions possess advantageous characteristics for representing heavy-tailed phenomena, and they can surpass alternative modelling methods in a noteworthy manner. Moreover, this type of distributions is a specific expansion of the Inhomogeneous Phase Type (IPH) class [2].

Due to the heavy tailed nature of MML distributions, they become very useful in modeling events in many real-life applications in queueing systems, especially when inter-arrival times are heavy-tailed. Heavy-tailed distributions appear naturally in queueing models (see [10, 15, 29, 14]). In the past 100 years, researchers and practitioners have given queueing models a lot of attention, resulting in an extensive collection of papers and books (see, [13, 17, 26, 7, 1, 11, 27]). Recently, researchers have examined fractional queue models that involve catastrophes in the model (see [5, 28, 11, 12, 25]). There is a correlated fractional Erlang queue model given in [6]. M/M/1 queue model employs exponentially distributed inter-arrival and service times for the clients, making it suitable for rigorous mathematical analysis and serving as a useful starting point for various scenarios. As a basic model, it does not include certain characteristics, such as memory preservation, which are frequently needed for analysing more intricate systems. There are multiple instances where memory effects are essential. For this reason, more flexible interarrival-time models are needed in queueing theory. PH distributions and, more generally, Mittag–Leffler-type interarrival times are well suited for this purpose, as they can capture heavy-tailed behavior and dependence structures while still retaining analytical tractability. These features make them realistic and powerful tools for modeling complex arrival processes in modern queueing systems.

In this article, we derive explicit forms of the MML density functions for the truncated queueing systems PH/M/ 1 ∼ n and M/PH/ 1 ∼ n, where n denotes the truncation level at which the queueing process is absorbed. We first study PH/M/1 queues with PH-distributed inter-arrival times, focusing on the exponential and Erlang cases, and obtain explicit expressions for the corresponding MML densities. We then consider M/PH/1 queues with PH-distributed service times and derive analogous formulas when the service-time distribution is exponential or Erlang. These results provide tractable representations of the underlying MML distributions and illustrate how PH-structured arrival and service mechanisms naturally give rise to Mittag-Leffler type behavior in queueing systems. Such density representations may be useful in modeling systems exhibiting heavy-tailed inter-arrival or service-time behavior.

The rest of the paper is organized as follows. In Section 2, we define a PH distribution and a matrix Mittag-Leffler function. We derive the density function of MML random variable for the PH / M /1 ∼ n system in Section 3. Section 4 presents the density function of MML random variable for M/PH/1 ∼ n system. The article is concluded with some discussions in Section 5.

 Phase-Type Distributions(PH-distributions) and Matrix Mittag-Leffler function Infinitesimal generator

An infinitesimal generator is a square matrix, of order infinite or finite, such that (i) every diagonal element is non-positive,(ii) every off-diagonal element is non-negative,(iii) each row sum is zero.

Infinitesimal generators arise naturally in continuous-time Markov chains (CTMC). For a CTMC with generator Q, the transition probabilities over time t ≥ 0 are given by the transition semigroup P ( t ) = e t Q . The ( i , j )th entry of P ( t ) represents the probability of being in state j at time t given that the chain started in state i. Thus, the matrix Q characterizes the instantaneous transition rates of the CTMC and determines its entire probabilistic evolution through e t Q . The diagonal element q i i represents the total rate of leaving state i, so the waiting time in state i is exponentially distributed with parameter − q i i . The off-diagonal element q i j for i ≠ j gives the rate of transition from state i to state j, and the ratio q i j − q i i gives the probability of jumping to state j when a transition occurs.

For example, a 5-order infinitesimal generator, corresponding to a continuous Markov chain with 5-states { 1 , 2 , 3 , 4 , 5 }, is given by(see, [19]) Q = − 5 0 1 0 4 1 − 10 0 0 9 0 2 − 5 0 3 0 0 1 − 1 0 0 0 0 0 0 . From the definition of a continuous-time Markov chain, the waiting time in each of the states, { 1 , 2 , 3 , 4 , 5 }, is exponentially distributed with parameters 5 , 10 , 5 , 1 and 0, respectively. State 5 is unique since its waiting time parameter is zero, indicating that there is no state transition after the Markov chain reaches this state. Thereafter, the Markov chain will always be in the state 5, the absorption state.

Now consider a continuous time Markov chain I ( t ) having n + 1 states { 1 , 2 , 3 , … , n + 1 } with its infinitesimal generator Q = T T 0 0 0 , where, T is a PH-generator of order n, that is, T is n × n invertible matrix such that(i) every diagonal element is negative,(ii) every off-diagonal element is non-negative,(iii) each row sum is non-positive.A PH-generator arises in the context of phase-type (PH) distributions, which model the time until absorption in a finite-state CTMC with one absorbing state. A formal definition of Phase-Type (PH) distributions will be given later in this section.

Since each row sum of matrix Q is zero, we have T 0 = − T e, where e is a column vector of order n, having all elements as one. As the total transition rate is zero, one of the states is absorption state, say n + 1.

Let I ( t ) denote the continuous-time Markov chain (CTMC) with state space { 1 , 2 , … , n + 1 } and infinitesimal generator Q, where the initial distribution is ( β , 0 ). Define a random variable X = min { t : I ( t ) = n + 1 , t ≥ 0 } , that is, the absorption time of state n + 1. We say that X is a Phase-Type (PH) random variable with PH representation ( β , T ) (Neuts [22]).

It is not difficult to check that the distribution of X is given by F X ( t ) = P { X ≤ t } = 1 − β e ( T t ) e ≡ 1 − β ( ∑ k = 0 ∞ t k k ! T k ) e , t ≥ 0 , For the Markov chain corresponding to Q, given in (1), state 5 is an absorption state. The absorption time of state 5 has a PH-distribution with m = 4 and T = − 5 0 1 0 1 − 10 0 0 0 2 − 5 0 0 0 1 − 1 , T 0 = 4 9 3 0 .

Matrix Mittag–Leffler(MML) function

Let A ∈ C n × n (set of n × n matrices of complex numbers) and, α and β be two complex numbers such that R e ( α ) > 0. The MML function has the following equivalent forms:([2, 24])

Taylor series form E α , β ( A ) = ∑ k = 0 ∞ A k Γ ( α k + β ) .

Jordan canonical formLet λ 1 , λ 2 , … λ p be distinct eigenvalues of A, then A can be expressed in the Jordan canonical form A = P diag ( J 1 , … , J p ) P − 1 with J k = λ k 1 0 … 0 0 λ k 1 … 0 0 0 λ k … 0 ⋮ ⋮ ⋮ … ⋮ 0 0 0 … λ k ∈ C m k × m k ,

and m 1 + ⋯ + m p = n.

The MML function is then expressed as E α , β ( A ) = P diag ( E α , β ( J 1 ) , … , E α , β ( J k ) ) P − 1 , where, E α , β ( J k ) = E α , β ( λ k ) E α , β ( 1 ) ( λ k ) E α , β ( 2 ) ( λ k ) 2 ! … E α , β ( m k − 1 ) ( λ k ) ( m k − 1 ) ! 0 E α , β ( λ k ) E α , β ( 1 ) ( λ k ) … E α , β ( m k − 2 ) ( λ k ) ( m k − 2 ) ! 0 0 E α , β ( λ k ) … E α , β ( m k − 3 ) ( λ k ) ( m k − 3 ) ! ⋮ ⋮ ⋮ … ⋮ 0 0 0 … E α , β ( λ k ) ,

E α , β ( i ) ( λ k ) is i t h derivative of E α , β ( λ k ) .

Cauchy integral form Consider a simple closed rectifiable curve Γ that encloses the spectrum of A. The MML function is then given by E α , β ( A ) = 1 2 π i ∫ Γ E α , β ( z ) ( z I − A ) − 1 d z .

In this article, we use the Cauchy-integral form of MML function.

 Matrix Mittag–Leffler distribution

As mentioned earlier, the random variable X is said to have a MML distribution if its Laplace transform is given by E ( e − u X ) = π ( u α I − T ) − 1 t , and we write X ∼ M M L ( α , π , T ).

Here, π denotes a stochastic row vector and T is a PH-generator.

The density function of X is given by f ( x ) = x α − 1 π E α , α ( T x α ) t , x > 0 . (see, [2], for more details).

We note that the MML distribution may admit a probabilistic representation of the form X = S α Y, where S α is a positive α-stable random variable and Y is a PH-distributed random variable as introduced in [2, 4, 3]. Such a representation would facilitate simulation of the MML distribution, since S α can be generated using Kanter’s formula and PH distributions can be simulated using standard algorithms. A full derivation of this representation is beyond the scope of this work, but it represents an interesting direction for future research.

 Renewal process (Renewal Process).

A stochastic process { N ( t ) , t ≥ 0 } is called a renewal process if it counts the number of events (renewals) that have occurred by time t, where the interarrival times between consecutive events form an independent and identically distributed (i.i.d.) sequence of non–negative random variables { X n } n ≥ 1 . Let S n = X 1 + X 2 + ⋯ + X n be the time of the n-th renewal. Then the renewal process is defined by N ( t ) = max { n : S n ≤ t } , t ≥ 0 .

Renewal processes are widely used to model systems that reset or renew after random time intervals, such as component replacements in reliability theory or recurrent events in stochastic modeling.

 The PH/M/1 queue

In the PH/M/1 queue customers arrive according to a PH-renewal process with PH- representation ( β , T ) (see, [19] Section 4.1). For each arriving customer, we initiate a CTMC with generator Q given in Section 2.1, and initial distribution ( β , 0 ). The next customer arrival occurs when this CTMC reaches its absorbing state. By concatenating these CTMCs over successive arrivals, we define the phase of the arrival process at time t as I a ( t ). Let q ( t ) denote the length of the queue at time t.

The following result is due to [19].

The stochastic process { ( q ( t ) , I a ( t ) ) , t ≥ 0 } is a continuous time Markov chain with state space { 0 , 1 , 2 , ... } × { 1 , 2 , ... , m a }, and infinitesimal generator Q = T T 0 β μ I − μ I + T T 0 β μ I − μ I + T T 0 β ⋱ ⋱ ⋱ , where, T is PH-generator of order m a .

In this subsection, we consider a truncated version of the PH/M/1 queue.

MML distribution of PH/M/1 <inline-formula id="j_vmsta303_ineq_063"><alternatives><mml:math> <mml:mo stretchy="false">∼</mml:mo> <mml:mi mathvariant="italic">n</mml:mi></mml:math><tex-math><![CDATA[$\sim n$]]></tex-math></alternatives></inline-formula> system

Here, the notation PH/M/1 ∼ n indicates that the underlying PH/M/1 queue is stopped once the queue length reaches level n + 1 that is, the state n + 1 is taken to be absorbing. The process evolves like a standard PH/M/1 queue until this threshold is hit. Let q ( t ) denote the number of customers in the system at time t, and let I a ( t ) be the phase of the arrival process at time t. The pair ( q ( t ) , I a ( t ) ) therefore forms a continuous-time Markov chain on the state space { 0 , 1 , 2 , ... , n } × { 1 , 2 , ... , m a }, with absorption occurring when q ( t ) = n + 1.

The PH-generator for { ( q ( t ) , I a ( t ) ) , t ≥ 0 } is defined as T = 1 T T 0 β 2 μ I − μ I + T T 0 β 3 μ I − μ I + T T 0 β ⋮ ⋱ ⋱ ⋱ n + 1 μ I − μ I + T . The next result proves the existence of a PH-distributed random variable.

Let π 1 = ( 1 , 0 , 0 , … , 0 ) 1 × ( n + 1 ) and, β be a stochastic vector of PH distribution of PH/M/1 ∼ n queue. If π = ( β ⨂ π 1 ), where ⨂ is Kronecker product and, T is the PH-generator of PH/M/1 ∼ n system. Then, there exist a non-negative random variable X that follows PH-distribution with parameter ( π , T ).

The result follows from (2).  □

In Theorems 3–6, we consider the special case where the PH interarrival distribution is exponential with rate λ. Under this assumption, the PH/M/ 1 ∼ n queue reduces to an M/M/ 1 ∼ n queue, where interarrival times are exponentially distributed with rate λ and service times are exponentially distributed with rate μ.

Let X ∼ M M L ( α , π , T ), where ( π , T ) (with PH-generator T having order 2) is PH-representation of PH/M/1 ∼ 1. Then, the probability density function of X is given by f ( x ) = x 2 α − 1 λ 2 ( E α , α ( s 1 ) ( s 1 − s 2 ) + E α , α ( s 2 ) ( s 2 − s 1 ) ) with s 1 , s 2 = ( − 2 λ − μ ± μ 2 + 4 μ λ ) 2 x α , where, λ is inter-arrival rate and μ is service rate.

Note that, for n = 1, T = − λ λ μ − ( μ + λ ) . Clearly, t = − Te = ( 0 , λ ) T .

Thus, ( s I − x α T ) − 1 = 1 s 2 + s ( 2 λ + μ ) x α + λ 2 x 2 α s + μ x α + λ x α x α λ x α μ s + x α λ . Using (3), and the Cauchy integral form, we get f ( x ) = x α − 1 1 2 π i ∫ Γ E α , α ( s ) π ( s I − x α T ) − 1 t d s . This implies, using π = ( 1 , 0 ), f ( x ) = x α − 1 1 2 π i ∫ Γ E α , α ( s ) ( x α λ 2 ) s 2 + s ( 2 λ + μ ) x α + λ 2 x 2 α d s . Using Cauchy’s residue theorem, we get f ( x ) = x 2 α − 1 λ 2 ( E α , α ( s 1 ) ( s 1 − s 2 ) + E α , α ( s 2 ) ( s 2 − s 1 ) ) , where s 1 = ( − 2 λ − μ + μ 2 + 4 μ λ ) 2 x α , s 2 = ( − 2 λ − μ − μ 2 + 4 μ λ ) 2 x α .  □

The following example demonstrates the above result for specific values of parameters.

Consider λ = 2 and μ = 1, that is, T = − 2 2 1 − 3 , and, t = − Te = ( 0 , 2 ) T .

In this case, s 1 = − x α , s 2 = − 4 x α

Thus, we get f ( x ) = 4 x 2 α − 1 ( E α , α ( − x α ) 3 x α + E α , α ( − 4 x α ) − 3 x α ) = 4 3 x α − 1 ( E α , α ( − x α ) − E α , α ( − 4 x α ) ) .

To provide numerical support for the explicit MML density derived in Example 1, we compute and plot the density using the closed-form expression given in that example. This numerical illustration confirms the validity and numerical stability of the derived formula and demonstrates its potential applicability in practical queueing scenarios.

MML density for the PH/M/1 ∼ 1 queue with α = 0.5

Plot of the probability density function f(x) versus x for the PH/M/1∼1 queue with λ=2, μ=1 and α=0.5. The curve decreases monotonically, starting at approximately 1.17 when x=0 and approaching 0 as x increases.

The next result is extension of Theorem (3) for underlying PH-generator of order 3.

Let X ∼ M M L ( α , π , T ), where ( π , T ) (with PH-generator T having order 3) is PH-representation of PH/M/1 ∼ 2. Then, the density function of X is given by f ( x ) = x α − 1 2 π i ∫ Γ E α , α ( s ) ( x 2 α λ 3 ) ( s 3 + s 2 ( 3 λ + 2 μ ) x α + s ( 3 λ 2 + 2 λ μ + μ 2 ) x 2 α d s , + λ 3 x 3 α ) where λ is inter-arrival rate and μ is service rate.

Note that, for n = 2 T = − λ λ 0 μ − ( μ + λ ) λ 0 μ − ( μ + λ ) . Thus, t = − Te = ( 0 , 0 , λ ) T .

Using (3), and the Cauchy integral form, we get f ( x ) = x α − 1 1 2 π i ∫ Γ E α , α ( s ) π ( s I − x α T ) − 1 t d s . This implies, using π = ( 1 , 0 , 0 ) and ( s I − x α T ) − 1 , we get the desired result.  □

The following corollary is a special case when all the roots of s 3 + s 2 ( 3 λ + 2 μ ) x α + s ( 3 λ 2 + 2 λ μ + μ 2 ) x 2 α + λ 3 x 3 α = 0 are distinct.

If all the roots of (4) are distinct then the density function has the following form f ( x ) = x 3 α − 1 λ 3 ∑ k = 0 2 E α , α ( s k ) ( s k − s k + 1 ) ( s k − s k + 2 ) , where s k = − 1 3 ( ( 3 λ + 2 μ ) x α + ξ k C + ( μ 2 + 6 μ λ ) x 2 α ξ k C ) , k ∈ Z 3 , here ξ k is the cube roots of unity and Z 3 is a group under addition mod 3 and C = x α ( 9 μ 2 λ − 2 μ 3 ) ± i 351 μ 4 λ 2 + 108 μ 5 λ + 864 μ 3 λ 3 2 3 .

The result follows by using the Cauchy’s residue theorem.  □

We consider the following example to demonstrate the above result.

Consider λ = 1 and μ = 1, then we have. T = − 1 1 0 1 − 2 1 0 1 − 2 . Thus, t = − Te = ( 0 , 0 , 1 ) T .

Using the above corollary, we get f ( x ) = x 3 α − 1 ∑ k = 0 2 E α , α ( s k ) ( s k − s k + 1 ) ( s k − s k + 2 ) , k ∈ Z 3 where s 0 = − 3.25 x α , s 1 = − 1.55 x α , s 2 = − 0.2 x α . Therefore, f ( x ) = x α − 1 ( E α , α ( − 3.25 x α ) 5.16 − E α , α ( − 1.55 x α ) 2.3 + E α , α ( − 0.2 x α ) 4.14 ) .

The next result is extension of Theorem (4) when the underlying PH-generator has order 4.

Let X ∼ M M L ( α , π , T ), where ( π , T ) (with PH-generator T having order 4) is PH-representation of PH/M/1 ∼ 3. Then, the density function of X is given by f ( x ) = x α − 1 2 π i ∫ Γ E α , α ( s ) ( x 3 α λ 4 ) ( s 4 + s 3 ( 4 λ + 3 μ ) x α + s 2 ( 6 λ 2 + 6 λ μ + 3 μ 2 ) x 2 α d s , + s ( 4 λ 3 + 3 λ 2 μ + 2 λ μ 2 + μ 3 ) x 3 α + λ 4 x 4 α ) where λ is arrival rate and μ is service rate.

Note that, for n = 3, T = − λ λ 0 0 μ − ( μ + λ ) λ 0 0 μ − ( μ + λ ) λ 0 0 μ − ( μ + λ ) . Thus, t = − Te = ( 0 , 0 , 0 , λ ) T .

Using (3), and the Cauchy integral form, we have f ( x ) = x α − 1 1 2 π i ∫ Γ E α , α ( s ) π ( s I − x α T ) − 1 t d s . Using π = ( 1 , 0 , 0 , 0 ) and ( s I − x α T ) − 1 , we get the desired result.  □

The following corollary is a special case when all the roots of s 4 + s 3 ( 4 λ + 3 μ ) x α + s 2 ( 6 λ 2 + 6 λ μ + 3 μ 2 ) x 2 α + s ( 4 λ 3 + 3 λ 2 μ + 2 λ μ 2 + μ 3 ) x 3 α + λ 4 x 4 α = 0 are distinct.

If all the roots of (5) are distinct then the density function has the following form f ( x ) = x 4 α − 1 λ 4 ∑ k = 0 3 E α , α ( s k ) ( s k − s k + 1 ) ( s k − s k + 2 ) ( s k − s k + 3 ) , k ∈ Z 4 , where Z 4 is a group under addition mod 4 and s 0 , s 1 = − ( 4 λ + 3 μ 4 ) x α − S ± 1 2 − 4 S 2 − 2 p + q S s 2 , s 3 = − ( 4 λ + 3 μ 4 ) x α + S ± 1 2 − 4 S 2 − 2 p − q S . Here, p = − 3 μ 8 ( 8 λ + μ ) x 2 α , q = μ 2 8 ( 4 λ − μ ) x 3 α , S = 1 2 − 2 3 p + 1 3 ( Q + ( 21 μ 2 λ 2 + 6 μ 3 λ ) x 4 α Q ) and Q = x 2 α 162 μ 3 λ 3 + 81 μ 4 λ 2 + − 10800 λ 6 μ 6 − 5508 μ 7 λ 5 − 2511 μ 8 λ 4 − 864 μ 9 λ 3 2 3 .

The result follows by using the Cauchy’s residue theorem.  □

Consider λ = 1 and μ = 1, then we have. T = − 1 1 0 0 1 − 2 1 0 0 1 − 2 1 0 0 1 − 2 . Thus, t = − Te = ( 0 , 0 , 0 , 1 ) T .

using the above corollary, we get f ( x ) = x 4 α − 1 ∑ k = 0 3 E α , α ( s k ) ( s k − s k + 1 ) ( s k − s k + 2 ) ( s k − s k + 3 ) , k ∈ Z 4 where, s 0 = − 0.12 x α , s 1 = − x α , s 2 = − 2.35 x α , s 3 = − 3.53 x α . Then f ( x ) = x α − 1 ( E α , α ( − 0.12 x α ) 6.69 − E α , α ( − x α ) 3.01 + E α , α ( − 2.35 x α ) 3.55 − E α , α ( − 3.53 x α ) 10.18 ) .

In Theorem (5), we assume that the PH-generator has order 4, the same can be generalized for PH- generator having order n. For n ≥ 5, obtaining a closed-form expression for the density is generally impossible, since the computation requires solving a polynomial of degree n, and polynomials of degree n ≥ 5 do not admit closed-form roots in radicals (Abel–Ruffini theorem). Nevertheless, studying the polynomial appearing in the denominator of the Cauchy integral representation of the MML function may still reveal useful structural properties and represents an interesting direction for further investigation. We obtain the density for MML-random variable corresponding to n-order PH-generator when μ = 0.

Consider PH/M/1 ∼ ( n − 1 ) with μ = 0. Then, the following random variables follow same distribution with common density function(a) X ∼ M M L ( α , π 1 , T 1 ), and (b) X ∼ M M L ( α , π 2 , T 2 ),where ( π 1 , T 1 ) is the PH-representation, given in Theorem (2), of exponential interarrival time with rate λ and ( π 2 , T 2 ) is the PH-representation of Erlang distribution with shape parameter n and rate λ.

Note that T 1 = T 2 = − λ λ 0 … 0 0 0 − λ λ … 0 0 0 0 − λ … 0 0 ⋮ ⋮ ⋮ … ⋮ ⋮ 0 0 0 … − λ λ 0 0 0 … 0 − λ (see, [2] and Example 3.1 for T 2 )

Now f ( x ) = x α − 1 π E α , α ( T x α ) t = λ n x α n − 1 ( n − 1 ) ! E α , α ( n − 1 ) ( − λ x α ) .  □

So far, we used PH-representation of exponential distribution. We now replace the assumption of exponential by Erlang and obtain the density of MML random variable.

Let X ∼ M M L ( α , π , T ), where ( π , T ) is PH-representation of PH/M/1 ∼ 1, where the interarrival times follow an Erlang-2 distribution with rate λ. Then, the density function of X is given by f ( x ) = x α − 1 2 π i ∫ Γ E α , α ( s ) ( x 3 α λ 4 ) d s ( s 4 + s 3 ( 4 λ + 2 μ ) x α + s 2 ( 6 λ 2 + 6 λ μ + μ 2 ) x 2 α d s . + s ( 4 λ 3 + 4 λ 2 μ + 2 λ μ 2 ) x 3 α + λ 4 x 4 α )

Note that T = T T 0 β μ I − μ I + T . Using the PH-representation of Erlang-2 distribution, we have T = − λ λ 0 0 0 − λ λ 0 μ 0 − ( μ + λ ) λ 0 μ 0 − ( μ + λ ) . Thus, t = − Te = ( 0 , 0 , 0 , λ ) T .

Using (3), and the Cauchy integral form, we get f ( x ) = x α − 1 1 2 π i ∫ Γ E α , α ( s ) π ( s I − x α T ) − 1 t d s . The result follows by using π = ( 1 , 0 , 0 , 0 ) and ( s I − x α T ) − 1 .  □

The following corollary gives an explicit form of density, when the roots of ( s 4 + s 3 ( 4 λ + 2 μ ) x α + s 2 ( 6 λ 2 + 6 λ μ + μ 2 ) x 2 α + s ( 4 λ 3 + 4 λ 2 μ + 2 λ μ 2 ) x 3 α + λ 4 x 4 α ) = 0 are distinct.

If all the roots of (7) are distinct then the density function (6) has the form f ( x ) = x 4 α − 1 λ 4 ∑ k = 0 3 E α , α ( s k ) ( s k − s k + 1 ) ( s k − s k + 2 ) ( s k − s k + 3 ) , k ∈ Z 4 , where, Z 4 is a group under addition mod 4 and s 0 , s 1 = − ( 4 λ + 2 μ 4 ) x α − S ± 1 2 − 4 S 2 − 2 p + q S , s 2 , s 3 = − ( 4 λ + 2 μ 4 ) x α + S ± 1 2 − 4 S 2 − 2 p − q S , where p = − μ 2 2 x 2 α , q = − 2 λ μ 2 x 3 α , S = 1 2 − 2 3 p + 1 3 ( Q + μ 4 x 4 α Q ) and Q = x 2 α 54 μ 2 λ 4 + μ 6 + 2916 λ 8 μ 4 + 108 μ 8 λ 4 3

The result follows by using the Cauchy’s residue theorem.  □

Consider λ = μ = 1, then we have T = − 1 1 0 0 0 − 1 1 0 1 0 − 2 1 0 1 0 − 2 . This implies t = − Te = ( 0 , 0 , 0 , 1 ) T .

Using above corollary and π = ( 1 , 0 , 0 , 0 ), we get f ( x ) = x 4 α − 1 ∑ k = 0 3 E α , α ( s k ) ( s k − s k + 1 ) ( s k − s k + 2 ) ( s k − s k + 3 ) , k ∈ Z 4 where s 0 = − 0.12 x α , s 1 = − 1.47 x α , s 2 , s 3 = ( − 2.21 ± i 0.98 ) x α . Then f ( x ) = x α − 1 ( E α , α ( − 0.12 x α ) 7.19 − E α , α ( − 1.47 x α ) 2.04 + E α , α ( ( − 2.21 + .98 i ) x α ) 5.44 + 1.15 i + E α , α ( ( − 2.21 − .98 i ) x α ) 5.44 − 1.15 i ) .

In Theorem (7), we assume that the number of persons in the queue is 1 and inter-arrival time follows Erlang-2 distribution. The same can be generalized for n persons in the queue and Erlang-n inter-arrival time. However, it is difficult to obtain a closed form expression for the density.

Next we will derive the explicit form of MML distribution for queue in which service time follows a PH distribution.

 The M/PH/1 queue

We assume that customers arrive with arrival process being Poisson with parameter λ and service time has a common PH-distribution with service rate μ, and PH-representation ( γ , S ), with S as a PH-generator of order m s .

Let ( A n ) n ≥ 1 be the sequence of arrival times. We define the sequence of departure times ( D n ) n ≥ 1 together with the service phase process I s ( t ) as follows. We set I s ( t ) = 0 for t < A 1 . For t ≥ A 1 , the process I s ( t ) evolves as a continuous-time Markov chain with generator S and initial distribution ( γ , 0 ), and we define D 1 as its time of absorption. Suppose that D n has been defined. If D n < A n + 1 , then we set I s ( t ) = 0 for t ∈ [ D n , A n + 1 ). At time t = A n + 1 , the process I s ( t ) restarts as a continuous-time Markov chain with generator S and initial distribution ( γ , 0 ), and D n + 1 is defined as its time of absorption. If instead D n ≥ A n + 1 , then I s ( t ) continues evolving as the same Markov chain for t ≥ D n , and D n + 1 is again defined as its time of absorption. In this way, we obtain a process I s ( t ), t ≥ 0, which is equal to 0 when the system is empty and describes the service phase otherwise. Moreover, the service times S n = D n − max { D n − 1 , A n } are PH-distributed with PH representation ( γ , S ).

Consider the stochastic process { ( q ( t ) , I s ( t ) ) , t ≥ 0 }. Because the arrival process is a Poisson process, the inter-arrival time t has exponential distribution with rate λ. Furthermore, because the service phase I s ( t ) is recorded, it is possible to calculate the probability of service time completion. Then, it is easy to observe that { ( q ( t ) , I s ( t ) ) , t ≥ 0 } is a continuous time Markov chain with state space { 0 , 1 , 2 , ... } × { 1 , 2 , ... , m s } ∪ { 0 } having infinitesimal generator Q = − λ λ γ S 0 − λ I + S λ I S 0 γ − λ I + S λ I ⋱ ⋱ ⋱ , where, S 0 = − S e. The process { ( q ( t ) , I s ( t ) ) , t ≥ 0 } is a QBD (Quasi Birth and Death) process, (see, [19, 20]).

MML distribution for M/PH/1 <inline-formula id="j_vmsta303_ineq_169"><alternatives><mml:math> <mml:mo stretchy="false">∼</mml:mo> <mml:mi mathvariant="italic">n</mml:mi></mml:math><tex-math><![CDATA[$\sim n$]]></tex-math></alternatives></inline-formula> system

The notation M/PH/1 ∼ n is defined analogously to the PH/M/1 ∼ n case in Section 3.1 with state n + 1 taken as absorbing.

Note that the PH-generator, for { ( q ( t ) , I s ( t ) ) , t ≥ 0 }, is defined as T = 1 − λ λ γ 2 S 0 − λ I + S λ I 3 S 0 γ − λ I + S λ I ⋮ ⋱ ⋱ ⋱ n + 1 S 0 γ − λ I + S . The next result proves the existence of a PH-distributed random variable.

Consider M/PH/1 ∼ n system with row vector π = ( γ , 0 , 0 , … , 0 ) and PH-generator T. There exists a random variable X which follows PH-distribution with representation ( π , T ).

The result follows from (2).  □

We use Theorem (8) and obtain the distribution of a MML random variable for Erlang underlying PH-distribution. We first obtain the density function of MML random variable, when there is one person in the queue.

Let X ∼ M M L ( α , π , T ), where ( π , T ) is PH-representation of M/PH/1 ∼ 1 system. If the service time follows PH-distribution with PH-representation of an Erlang-2 distribution having service rate μ, that is, γ = ( 1 , 0 ). Then, the density function of X is given by f ( x ) = x α − 1 2 π i ∫ Γ E α , α ( s ) s λ x α + ( λ 2 + 2 μ λ ) x 2 α ( s 3 + s 2 ( 3 λ + 2 μ ) x α + s ( 3 λ 2 + μ 2 + 4 μ λ ) x 2 α d s , + ( λ 3 + 2 λ 2 μ ) x 3 α ) where λ is the arrival rate.

Note that T = − λ λ γ S 0 − λ I + S , or T = − λ λ 0 0 − μ − λ μ μ 0 − μ − λ . Thus, t = − Te = ( 0 , λ , λ ) T . We have ( s I − x α T ) − 1 = 1 ( s 3 + s 2 ( 3 λ + 2 μ ) x α + s ( 3 λ 2 + μ 2 + 4 μ λ ) x 2 α + ( λ 3 + 2 λ 2 μ ) x 3 α ) × A B C D E F G H I where A = s 2 + s ( 2 μ + 2 λ ) x α + ( λ 2 + μ 2 + 2 μ λ ) x 2 α , B = s λ x α + ( λ 2 + μ λ ) x 2 α , C = μ λ x 2 α , D = μ 2 x 2 α , E = s 2 + s ( 2 λ + μ ) x α + ( λ 2 + μ λ ) x 2 α , F = s μ x α + μ λ x 2 α , G = s μ x α + ( μ 2 + μ λ ) x 2 α , H = μ λ x 2 α , and , I = s 2 + s ( 2 λ + μ ) x α + ( λ 2 + μ λ ) x 2 α

Using (3), and the Cauchy integral form, we have f ( x ) = x α − 1 1 2 π i ∫ Γ E α , α ( s ) π ( s I − x α T ) − 1 t d s . The result follows by substituting ( s I − x α T ) − 1 and π = ( 1 , 0 , 0 ) in the above integration.  □

The following corollary is a special case of the above theorem when all the roots of ( s 3 + s 2 ( 3 λ + 2 μ ) x α + s ( 3 λ 2 + μ 2 + 4 μ λ ) x 2 α + ( λ 3 + 2 λ 2 μ ) x 3 α ) = 0 are distinct.

If all the roots of (9) are distinct then the density function of X is given by f ( x ) = x 2 α − 1 λ ∑ k = 0 2 ( s k + ( λ + 2 μ ) x α ) E α , α ( s k ) ( s k − s k + 1 ) ( s k − s k + 2 ) , where s k = − 1 3 ( ( 3 λ + 2 μ ) x α + ξ k C + μ 2 x 2 α ξ k C ) , k ∈ Z 3 . Here, ξ k is the cube roots of unity and Z 3 is a group under addition mod 3, and, C = x α ( − 2 μ 3 − 27 μ 2 λ ) ± 729 μ 4 λ 2 + 108 μ 5 λ 2 3 .

The result follows by using the Cauchy’s residue theorem.  □

The following example demonstrate the above result for specific values of parameters.

Consider λ = 2 and μ = 2. From (8), we have T = − 2 2 0 0 − 4 2 2 0 − 4 . Thus, t = − Te = ( 0 , 2 , 2 ) T .

Using above corollary and π = ( 1 , 0 , 0 ) we get f ( x ) = x α − 1 1 2 π i ∫ γ E α , α ( s ) 2 x α s + 12 x 2 α s 3 + 10 x α s 2 + 32 x 2 α s + 32 x 3 α d s = 2 x α − 1 ( E α , α ( − 2 x α ) − E α , α ( − 4 x α ) − E α , α ′ ( − 4 x α ) ) , where E α , α ′ ( − 4 x α ) is the derivative of E α , α ( − 4 x α ).

The next result is extension of Theorem (9), where the underlying PH-distribution is Erlang-3 distribution.

Let X ∼ M M L ( α , π , T ), where ( π , T ) is PH-representation of M/PH/1 ∼ 1 system. If the service time follows PH-distribution with PH-representation of an Erlang-3 distribution with service rate μ, that is, γ = ( 1 , 0 , 0 ). Then, the density function of X is given by f ( x ) = x α − 1 2 π i ∫ Γ E α , α ( s ) ( λ x α s 2 + s ( 2 λ 2 + 3 λ μ ) x 2 α + ( λ 3 + 3 λ 2 μ + 3 λ μ 2 ) x 3 α ) d s ∑ j = 0 2 ∑ i = 0 4 − j 3 j 4 − j i λ i μ j x ( i + j ) α s 4 − j − i + μ 3 x 3 α s , where λ is the arrival rate.

Note that we have T = − λ λ γ S 0 − λ I + S , that is, T = − λ λ 0 0 0 − ( μ + λ ) μ 0 0 0 − ( μ + λ ) μ μ 0 0 − ( μ + λ ) . Then t = − Te = ( 0 , λ , λ , λ ) T .

Using (3), and the Cauchy integral form, we have f ( x ) = x α − 1 1 2 π i ∫ Γ E α , α ( s ) π ( s I − x α T ) − 1 t d s . The result follows by using π = ( 1 , 0 , 0 , 0 ) and substituting the value of ( s I − x α T ) − 1 for T given in (10).  □

The following corollary is a special case when all the roots of ∑ j = 0 2 ∑ i = 0 4 − j 3 j 4 − j i λ i μ j x ( i + j ) α s 4 − j − i + μ 3 x 3 α s , are distinct.

If all the roots of (11) are distinct then the density function has the form f ( x ) = x 2 α − 1 λ ∑ k = 0 3 ( s k 2 + s k ( 2 λ + 3 μ ) x α + ( λ 2 + 3 λ μ + 3 μ 2 ) x 2 α ) E α , α ( s k ) ( s k − s k + 1 ) ( s k − s k + 2 ) ( s k − s k + 3 ) , k ∈ Z 4 , where Z 4 is a group under addition mod 4, and s 0 , s 1 = − ( 4 λ + 3 μ 4 ) x α − S ± 1 2 − 4 S 2 − 2 p + q S s 2 , s 3 = − ( 4 λ + 3 μ 4 ) x α + S ± 1 2 − 4 S 2 − 2 p − q S , where p = − 3 μ 2 8 x 2 α , q = − μ 3 8 x 3 α , S = 1 2 − 2 3 p + 1 3 ( Q + − 12 μ 3 λ x 4 α Q ) and Q = x 2 α − 27 μ 5 λ + 729 λ 2 μ 10 − 6912 μ 9 λ 3 2 3 .

The result follows by using the Cauchy’s residue theorem.  □

The following example illustrates the above result for specific values of parameters.

Consider λ = 1 and μ = 1. We have, from (10), T = − 1 1 0 0 0 − 2 1 0 0 0 − 2 1 1 0 0 − 2 . Thus t = − Te = ( 0 , 1 , 1 , 1 ) T .

Now f ( x ) = x 2 α − 1 ∑ k = 0 3 ( s k 2 + 5 x α s k + 7 x 2 α ) E α , α ( s k ) ( s k − s k + 1 ) ( s k − s k + 2 ) ( s k − s k + 3 ) , k ∈ Z 4 , Where s 0 − 0.62 x α , s 1 = − 2.82 x α , s 2 , s 3 = ( − 1.78 ± 0.91 i ) x α . Then f ( x ) = x α − 1 ( 4.28 E α , α ( − 0.62 x α ) 4.78126 − 0.85 E α , α ( − 2.82 x α ) 4.2416 + ( 0.44 + 1.31 i ) E α , α ( ( − 1.78 + .91 i ) x α ) 0.2 + 3.7 i + ( 0.44 − 1.31 i ) E α , α ( ( − 1.78 − .91 i ) x α ) 0.2 − 3.7 i )

The PH-distribution used in the above Theorem (10) is Erlang-3, we next generalize the result for Erlang-n distribution.

Let X ∼ M M L ( α , π , T ), where ( π , T ) is PH-representation of M/PH/1 ∼ 1 system. If the service time follows PH-distribution with PH-representation of an Erlang-n distribution with service rate μ, that is, γ = ( 1 , 0 , 0 , … , 0 ). Then the density function of X is given by f ( x ) = x α − 1 2 π i ∫ Γ E α , α ( s ) ∑ j = 0 n − 1 ∑ i = 0 n − ( j + 1 ) n j n − ( j + 1 ) i λ i + 1 μ j x ( i + j + 1 ) α s n − ( j + 1 ) − i d s ∑ j = 0 n − 1 ∑ i = 0 n + 1 − j n j n + 1 − j i λ i μ j x ( i + j ) α s n + 1 − j − i + n n μ n x n α s , where λ is the arrival rate.

Note that, T = − λ λ γ S 0 − λ I + S , that is, T = − λ λ 0 … 0 0 0 − λ − μ μ … 0 0 0 0 − λ − μ … 0 0 ⋮ ⋮ ⋮ … ⋮ ⋮ 0 0 0 … − λ − μ μ μ 0 0 … 0 − λ − μ . Thus, t = − Te = ( 0 , λ , λ , … , λ ) n + 1 × 1 T .

Using (3), and the Cauchy integral form, we have f ( x ) = x α − 1 1 2 π i ∫ Γ E α , α ( s ) π ( s I − x α T ) − 1 t d s . The result now follows by using π = ( 1 , 0 , 0 , … , 0 ) 1 × n + 1 and substituting the value of ( s I − x α T ) − 1 in the above integral.  □

Following corollary is a special case of the above theorem when all the roots of ∑ j = 0 n − 1 ∑ i = 0 n + 1 − j n j n + 1 − j i λ i μ j x ( i + j ) α s n + 1 − j − i + n n μ n x n α s = 0 are distinct.

If all the roots of (12) are distinct, then the density function is given by f ( x ) = x 2 α − 1 λ ∑ k = 0 n ( ∑ j = 0 n − 1 ∑ i = 0 n − ( j + 1 ) n j n − ( j + 1 ) i λ i μ j x ( i + j ) α s k n − ( j + 1 ) − i ) E α , α ( s k ) ( s k − s k + 1 ) ( s k − s k + 2 ) … ( s k − s k + n ) , where k ∈ Z n + 1 , a group under addition mod n + 1 and, s k is a root of (12).

The result follows by using the Cauchy’s residue theorem.  □

 Conclusion

We have considered a random variable arising from a queueing system whose distribution is Phase-Type, and used this PH representation to derive explicit forms of the associated MML distribution. The MML class is heavy-tailed, dense in the class of all lifetime distributions on the positive real line, and mathematically tractable, making it suitable for modeling systems in which inter-arrival or service times exhibit heavy-tailed behavior. In this work, we studied the MML distributions associated with the truncated queueing systems PH/M/ 1 ∼ n and M/PH/ 1 ∼ n, and obtained explicit density formulas when the PH distributions reduce to exponential or Erlang forms. Several examples were provided to illustrate the results. A numerical illustration was also included to demonstrate the stability of the derived expressions. The approach developed here may be extended by considering other subclasses of PH distributions, such as Coxian or generalized Erlang distributions, within the PH/M/ 1 ∼ n and M/PH/ 1 ∼ n frameworks. These extensions may further broaden the applicability of MML distributions in modeling complex queueing and stochastic systems.

 Acknowledgments

The authors would like to thank the anonymous reviewers for their careful reading of the manuscript and for their valuable comments and suggestions, which helped improve the quality of this paper.

 References Al Hanbali, A.: Busy period analysis of the level dependent PH/PH/1/K queue. Queueing Syst. 67(3), 221–249 (2011). https://doi.org/https://doi.org/10.1007/s11134-011-9213-6. MR2800612 Albrecher, H., Bladt, M., Bladt, M.: Matrix Mittag-Leffler distributions and modeling heavy-tailed risks. Extremes 23(3), 425–450 (2020). https://doi.org/https://doi.org/10.1007/s10687-020-00377-0. MR4129559 Albrecher, H., Bladt, M., Bladt, M.: Multivariate fractional phase-type distributions. Fract. Calc. Appl. Anal. 23(5), 1431–1451 (2020). https://doi.org/https://doi.org/10.1515/fca-2020-0071. MR4173831 Albrecher, H., Bladt, M., Bladt, M.: Multivariate matrix Mittag-Leffler distributions. Ann. Inst. Stat. Math. 73(2), 369–394 (2021). https://doi.org/https://doi.org/10.1007/s10463-020-00750-7. MR4233525 Ascione, G., Leonenko, N., Pirozzi, E.: Fractional queues with catastrophes and their transient behaviour. Mathematics 6(9), 159 (2018) MR4092404. https://doi.org/10.1016/j.spa.2019.09.012 Ascione, G., Leonenko, N., Pirozzi, E.: Fractional Erlang queues. Stoch. Process. Appl. 130(6), 3249–3276 (2020). https://doi.org/https://doi.org/10.1016/j.spa.2019.09.012. MR4092404 Asmussen, S.r.: Applied probability and queues. Springer. Stochastic Modelling and Applied Probability (2003). MR1978607 Asmussen, S.r., Albrecher, H.: Ruin probabilities. World Scientific Publishing Co. Pte. Ltd., Hackensack, NJ (2010). https://doi.org/https://doi.org/10.1142/9789814282536. MR2766220 Bean, N.G., Fackrell, M., Taylor, P.: Characterization of matrix-exponential distributions. Stoch. Models 24(3), 339–363 (2008). https://doi.org/https://doi.org/10.1080/15326340802232186. MR2436371 Butt, J., Georgiou, N., Scalas, E.: Queuing models with Mittag-Leffler inter-event times. Fract. Calc. Appl. Anal. 26(4), 1465–1503 (2023). https://doi.org/https://doi.org/10.1007/s13540-023-00161-4. MR4623217 Cahoy, D.O., Polito, F., Phoha, V.: Transient behavior of fractional queues and related processes. Methodol. Comput. Appl. Probab. 17(3), 739–759 (2015). https://doi.org/https://doi.org/10.1007/s11009-013-9391-2. MR3377858 Di Crescenzo, A., Giorno, V., Nobile, A.G., Ricciardi, L.M.: On the M / M / 1 queue with catastrophes and its continuous approximation. Queueing Syst. 43(4), 329–347 (2003). https://doi.org/https://doi.org/10.1023/A:1023261830362. MR1976263 Erlang, A.K.: Solution of some problems in the theory of probabilities of significance in automatic telephone exchanges. Post Off. Electr. Eng. J. 10, 189–197 (1917) Foss, S., Korshunov, D.: On large delays in multi-server queues with heavy tails. Math. Oper. Res. 37(2), 201–218 (2012). https://doi.org/https://doi.org/10.1287/moor.1120.0539. MR2931277 Foss, S., Korshunov, D.: Heavy tails in multi-server queue. Queueing Syst. 52(1), 31–48 (2006). https://doi.org/https://doi.org/10.1007/s11134-006-3613-z. MR2201624 Garrappa, R., Popolizio, M.: Computing the matrix Mittag-Leffler function with applications to fractional calculus. J. Sci. Comput. 77(1), 129–153 (2018). https://doi.org/https://doi.org/10.1007/s10915-018-0699-5. MR3850348 Goyal, T.L., Harris, C.M.: Maximum-likelihood estimates for queues with state-dependent service. Sankhya, Ser. A 34, 65–80 (1972). MR336841 Haubold, H.J., Mathai, A.M., Saxena, R.K.: Mittag-Leffler functions and their applications. J. Appl. Math., 298628–51 (2011). https://doi.org/https://doi.org/10.1155/2011/298628. MR2800586 He, Q.-M.: Fundamentals of matrix-analytic methods. Springer (2014). https://doi.org/https://doi.org/10.1007/978-1-4614-7330-5. MR3112230 Lakatos, L., Szeidl, L., Telek, M.: Introduction to queueing systems with telecommunication applications. Springer (2013). https://doi.org/https://doi.org/10.1007/978-1-4614-5317-8. MR2987305 Matychyn, I., Onyshchenko, V.: Matrix Mittag–Leffler function in fractional systems and its computation. Bull. Pol. Acad. Sci., Tech. Sci., 495–500 (2018) Neuts, M.F.: Computational uses of the method of phases in the theory of queues. Comput. Math. Appl. 1(2), 151–166 (1975). https://doi.org/https://doi.org/10.1016/0898-1221(75)90015-2. MR386055 Pillai, R.N.: On Mittag-Leffler functions and related distributions. Ann. Inst. Stat. Math. 42(1), 157–161 (1990). https://doi.org/https://doi.org/10.1007/BF00050786. MR1054728 Popolizio, M.: On the Matrix Mittag–Leffler function: theoretical properties and numerical computation. Mathematics 7(12), 1140 (2019) Pote, R.B., Kataria, K.K.: On mixed time-changed Erlang queue. Mod. Stoch. Theory Appl., 1–27 (2026) Prabhu, N.U.: Stochastic storage processes. Springer. Queues, insurance risk, dams, and data communication (1998). https://doi.org/https://doi.org/10.1007/978-1-4612-1742-8. MR1492990 Singh, S.K.: Maximum likelihood estimation in single server queues. Sankhya A 85(1), 931–947 (2023). https://doi.org/https://doi.org/10.1007/s13171-022-00283-6. MR4540825 Souza, M.d.O., Rodriguez, P.M.: On a fractional queueing model with catastrophes. Appl. Math. Comput. 410, 126468–14 (2021). https://doi.org/https://doi.org/10.1016/j.amc.2021.126468. MR4281788 Whitt, W.: The impact of a heavy-tailed service-time distribution upon the M / G I / s waiting-time distribution. Queueing Syst. 36(1-3), 71–87 (2000). https://doi.org/https://doi.org/10.1023/A:1019143505968. MR1806961